Amandine Regis
email: amandine.regis@univ-amu.fr

Santiago Arroyave-Tobon

-Jean-Marc Linares

Emmanuel Mermoz

Physic-based vs data-based digital twins for bush bearing wear diagnostic

Keywords:

Industry is increasingly turning to predictive maintenance by using digital twins (DTs) to follow and predict evolution of mechanical system. This article presents, compares and discusses two DTs to diagnose wear of bush bearings under dynamic loads. The first DT is driven by a model based on data analysis using statistical process control (SPC). The second DT is based on physical laws: Boussinesq's and Archard's equations. Both DTs are fed by data recorded on a test bench instrumented with sensors of temperature, acceleration and displacement. Rules for fault detection were identified, explained and applied to the two DTs implemented. The two implemented DTs detected abnormal wear behaviours. The data-based DT using the SPC is easier to implement and it detects change in wear behaviour earlier. In contrast, the physic-based DT has the advantage of being predictive, so it can be used when only the operating conditions of the system are known. This work is a contribution for new wear diagnostic tools.

Introduction

Wear is a critical phenomenon in dynamic systems of helicopters. Currently, wear monitoring is done by preventive maintenance at regular intervals. However, this type of maintenance is expensive [START_REF] Kraft | Engine Fleet-Management: The Use of Digital Twins From a MRO Perspective[END_REF]. Nowadays, industries are increasingly turning to predictive maintenance by using, for example, digital twins (DTs).

A DT is a virtual replica of a physical system and both exchange data and information throughout the whole life cycle of the physical system [START_REF] Grieves | Digital Twin: Manufacturing Excellence through Virtual Factory Replication[END_REF]. The DT perceives its environment through sensors [START_REF] Rosen | About The Importance of Autonomy and Digital Twins for the Future of Manufacturing[END_REF] and the data flow between the virtual replica and the physical system can be either unidirectional or bidirectional [START_REF] Kritzinger | Digital Twin in manufacturing: A categorical literature review and classification[END_REF].

In mechanical engineering, DTs have the advantage reducing maintenance costs because they allow to mirror systems where direct inspection is not easy [START_REF] Kraft | Engine Fleet-Management: The Use of Digital Twins From a MRO Perspective[END_REF]. In consequence, it allows increasing safety by detecting problems earlier [START_REF] Tao | Digital twin driven prognostics and health management for complex equipment[END_REF][START_REF] Wang | Digital Twin for rotating machinery fault diagnosis in smart manufacturing[END_REF]. They can be used to facilitate the lifecycle management [START_REF] Tao | Digital twin driven prognostics and health management for complex equipment[END_REF][START_REF] Tuegel | Reengineering Aircraft Structural Life Prediction Using a Digital Twin[END_REF], to support manufacturing processes [START_REF] Söderberg | Toward a Digital Twin for real-time geometry assurance in individualized production[END_REF][START_REF] Tao | Digital Twin Shop-Floor: A New Shop-Floor Paradigm Towards Smart Manufacturing[END_REF][START_REF] Zhang | A multi-scale modeling method for digital twin shop-floor[END_REF] and to monitor or predict mechanical damages [START_REF] Karve | Digital twin approach for damage-tolerant mission planning under uncertainty[END_REF][START_REF] Yeratapally | A digital twin feasibility study (Part I): Non-deterministic predictions of fatigue life in aluminum alloy 7075-T651 using a microstructurebased multi-scale model[END_REF][START_REF] De C.-E. Kerman | Towards a Circular Rotating Blade Wear Assessment Digital Twin for Manufacturing Lines[END_REF][START_REF] Desai | Lifetime Prediction Using a Tribology-Aware, Deep Learning-Based Digital Twin of Ball Bearing-Like Tribosystems in Oil and Gas[END_REF][START_REF] Rajesh | Digital Twin of an Automotive Brake Pad for Predictive Maintenance[END_REF][START_REF] Qin | Data-Model Combined Driven Digital Twin of Life-Cycle Rolling Bearing[END_REF].

DTs are supported by models based on data analysis or on physical laws.

Data-based models are created from the analysis of data recorded by the studied system sensors. The analysis of recorded data can be done through different approaches. Neural network and machine learning have been used to analyse the degradation performance of rolling bearing and to predict the remaining useful life of rolling bearing [START_REF] Desai | Lifetime Prediction Using a Tribology-Aware, Deep Learning-Based Digital Twin of Ball Bearing-Like Tribosystems in Oil and Gas[END_REF][START_REF] Qin | Data-Model Combined Driven Digital Twin of Life-Cycle Rolling Bearing[END_REF]. In [START_REF] De C.-E. Kerman | Towards a Circular Rotating Blade Wear Assessment Digital Twin for Manufacturing Lines[END_REF], the authors predicted the wear of a circular blade from historic data collected with a SCADA system. In [START_REF] Chen | Anomaly detection of the tapered roller bearings with statistical data-driven approaches[END_REF], a probabilistic Gaussian mixture model is used to develop a specific strategy for bearing fault detection.

Physic-based models use physical laws to describe the behaviour of the system studied. One way to implement physic-based models is by means of the finite element method (FEM). In [START_REF] Rajesh | Digital Twin of an Automotive Brake Pad for Predictive Maintenance[END_REF], the authors used FEM in a DT to simulate automotive brake pad wear. FEM with mesh updating was used to predict fretting wear of spline couplings [START_REF] Ratsimba | Measurement, analysis and prediction of fretting wear damage in a representative aeroengine spline coupling[END_REF] and to predict wear on classic bush bearings [START_REF] Schmidt | Transient wear simulation based on three-dimensional finite element analysis for a dry running tilted shaftbushing bearing[END_REF] and bio-inspired bush bearings [START_REF] Sysaykeo | Wear Behavior of a Bio-inspired Bearing for off-center Loads[END_REF]. Specific numerical models are another way to develop physic-based model. Specific numerical models have been developed to predict wear of bearings under dynamic loads [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF], revolute joint in 4-bar mechanisms [START_REF] Flores | Modeling and simulation of wear in revolute clearance joints in multibody systems[END_REF][START_REF] Lai | Computational prediction and experimental validation of revolute joint clearance wear in the low-velocity planar mechanism[END_REF], gears [START_REF] Akbarzadeh | Prediction of Steady State Adhesive Wear in Spur Gears Using the EHL Load Sharing Concept[END_REF] and crossed steel wires [START_REF] Chen | Bearing Condition Monitoring Using Multiple Sensors and Integrated Data Fusion Techniques[END_REF]. Specific numerical models have the advantage of being less time consuming than FEM [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF][START_REF] Flores | Modeling and simulation of wear in revolute clearance joints in multibody systems[END_REF]. However, in the case of wear monitoring, FEM is well adapted to complex contact geometries.

In this work, two DTs (a data-based and a physicbased) are developed, implemented and compared. This with aim of monitoring wear of a bush bearing system mounted on a test bench. The first DT relies on a data-based model using the statistical process control (SPC). The SPC allows to continuously monitor the evolution of wear rate of the bearings. The second DT uses the specific numerical model developed by [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF]. Rules of fault detection were defined to monitor a variable with linear evolution such as wear. These rules were applied on both DTs implemented.

The objective of this article is to monitor wear of a bush bearing system by means of DTs. Two DTs are implemented and compared in terms of implementation time, execution time and detection capability. This with the aim of identifying the advantages of each type of DT in the case of wear monitoring.

Material and Methods

The global architecture of both bush bearing DTs implemented in this work is presented in Fig. 1. This architecture can be divided into four components.

The first component (green box in in Fig. 1) is the physical system. In this work, it is a test bench of a bush bearing instrumented with temperature, acceleration and displacement sensors to monitor wear effects.

The second component (black box in in Fig. 1) corresponds to the transmission of data collected on the physical system to the DTs model. Data collected by the optical displacement sensor was used to identify the operating configuration (rotational speed and force loading) and to monitor wear effects over time. The other sensors were employed to control the stability of the temperature and vibrations during experiments.

The third component (blue box in in Fig. 1) is the DT to follow wear of the bush bearings. Two models were implemented in this work: one model based on data analysis using Statistical Process Control (SPC) and one model based on physic laws using a specific numerical model [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF].

The last component (red box in in Fig. 1) apply the rules of fault detection.

In this section, each component of the architecture of the DTs is detailed.

Physical system 2.1.1. Experimental test bench

The test bench of the bush bearing system is illustrated in Fig. 2. It is composed of a steel shaft OTMT 81001300 (item n°6 in Fig. 2) guided by two bronze bearings SINT A50 (item n°7 in Fig. 2). The material properties of the shaft and the bearings are summarized in Table 1. For the tests considered in this work, no lubrication was used, only dry contact to consider the worst operating case. The initial clearance on diameter between the shaft and the bearings was 0.144 mm. An unbalanced mass (50 or 75 g) fixed at the end of the shaft (item n°1 in Fig. 2) was used to create a dynamic load. A double universal joint was used (item n°4 in Fig. 2) to connect the shaft to the motor (item n°5 in Fig. 2) and reduce parasitic dynamic loads due to misalignment between the shaft and the electric motor shaft. The frame (item n°3 in Fig. 2) of the test bench was made of aluminium alloy to diffuse the temperature generated by the friction. Before each trial, the test bench was heated to the operating temperature (38-40°C) by using an external heat source. The rotational speed was constant and equal to 1000 rpm. The two unbalanced masses were alternatively used during tests. Table 2 presents the details about the operating conditions (OC) through time. The test bench was instrumented with several sensors. Two temperature sensors (IEC type K silicone rubber patch thermocouple) were fixed on the frame (item n°8 in Fig. 2). Two accelerometers (Model 603C01 -PCB Piezotronics) were used to monitor the stability of the experiment (item n°2 in Fig. 2): one was positioned to detect vertical vibrations and the other for horizontal vibrations. A confocal chromatic displacement sensor (Micro-Epsilon IFS2405-10) was mounted perpendicular to the shaft axis to measure its displacement (run-out) (item n°9 in Fig. 2). During the experiments, data was collected every 10 minutes for 2 seconds with a NI data acquisition card. The acquisition frequencies and the monitored parameters are presented in Table 3. 3. Acquisition parameters of the used sensors.

In addition, the wear depth of the shaft and the bearings was measured after disassembly. To allow the disassembly of the bearings, the frame was designed in two parts and assembled using 4 screws (ISO screw M6). The screws were tightened each time they were assembled using a torque wrench in order to guarantee an identical tightening torque at each reassembly. The tightening torque was set to 5 Nm. Fig. 3 illustrates the worn shaft and bearings. For the shaft, the mean wear depth was obtained by measuring profiles of 50 mm long on the wear zones (blue line in Fig. 3a) with a micromesure station (STIL Marposs, Micromesure 2). For the bearings, the inner surface was acquired with a coordinate-measuring machine (Messwelk, MM 1004 E Trimesures). One circular path at each end of the bearings was measured (blue lines in Fig. 3b). Each circular path was composed of 36 points uniformly distributed. The wear depth of the bearing was assumed to be the diameter variation between the initial and the final circular paths.

Identification of operating configurations

In this work, it is assumed that the loading to which the system is subjected is not known during the operation of the system.

In order to launch the digital twin, the operating conditions of the monitored system are automatically identified from experimental data. The general workflow of this procedure is represented in Fig. 4a. The features of the shaft displacement signal (Fig. 4b) and its Fast Fourier Transformation (FFT) (Fig. 4c) were used to do so. The loading was deduced from the amplitude of the first peak of the FFT of the signal. The signal processing toolbox of Matlab, and more precisely, the classification learner tool was used. A feed forward network with a single internal layer was trained with sets of data where the unbalanced mass was known.

Digital twins

In this section, the implementation of the DTs to diagnose wear of bush bearings under dynamic loads is described. DT models can be based on data or on physic laws. This article presents one example of each model. The data-based model was developed using SPC of the data recorded by the sensors. The physicbased model corresponds to a specific numerical model that predict the wear evolution of bearing under dynamic loading. All results presented in this article are for the stable wear stage.

Controlled variables

The aim of the DTs developed in this work is to control wear through time. The controlled variables are the wear depth and the wear rate. The wear depth is deduced from the peak-to-peak value of the shaft displacement signal (Fig. 5). From the displacement sensor measure d1 and by geometric relations, the displacement d0 at the end of the bearing can be calculated. This distance d0 corresponds to the sum of the wear depths of the two bearings and the clearance on the diameter. The wear depth equation is defined in Fig. 5. The wear rate is defined as the wear depth per time unit (µm/h).

Data-based model

SPC was used to control the mean and the standard deviation of the wear rate to study its trend and variation. Both, the mean and the standard deviation of the wear rate, were deduced from the measured displacements. As this work is focused on the stable wear stage, the wear depth is supposed to evolve linearly through time. A linear regression of the wear depth through the time was done to determine the mean wear rate which corresponds to the slope of the regression. This regression was carried out using both the least squares and the Moore-Penrose inverse methods [START_REF] Courrieu | Fast Computation of Moore-Penrose Inverse Matrices[END_REF][START_REF] Linares | Fatigue life optimization for 17-4Ph steel produced by selective laser melting[END_REF]. This allowed to determine the mean and the standard deviation of the wear rate.

The data controlled were supposed to follow normal distribution. The upper control limit (UCL) and the lower control limit (LCL) for the mean wear rate, respectively UCLµ and LCLµ, and the upper control limit of the standard deviation of the wear rate, UCLσ, are calculated as follows [START_REF]Shewhart control charts[END_REF]:

𝑈𝐶𝐿 𝜇 = 𝜇 + 𝑡 𝜎 √𝑛 (1)
𝐿𝐶𝐿 𝜇 = 𝜇 - 𝑡 𝜎 √𝑛 ( 2
)
𝑈𝐶𝐿 𝜎 = 𝐵 4 𝜎 (3)
where n is the size of the sample, µ the mean wear rate target (see Table 4 for numerical values), σ the standard deviation target of the wear rate (see Table 4 for numerical values), t the coverage factor determined by the inverse of Student's law with a risk of 0.25% and B4 the coefficient for the calculation of the control limit from [START_REF]Shewhart control charts[END_REF]. The mean target µ and the standard deviation target σ were obtained from calibration tests.

Physic-based model

The numerical model to predict wear of bush bearings under off-centred dynamic loads is summarized hereafter. The algorithm behind this model is presented in Fig. 6. For more details readers can refer to [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF]. Several input parameters are necessary to run the model. The time is initialized to t0. The rotational speed of the shaft ω, the unbalanced mass m and the wear coefficients K are required. The contact pressure distribution field P(x,z) and the wear distribution field h(x,z) are supposed to be sinusoidal [START_REF] Colbert | Edges, clearances, and wear: Little things that make big differences in bushing friction[END_REF] in the plane with normal x (Fig. 7a) and linear in the plane with normal z (Fig. 7b). The boundary of the fields on the shaft corresponds to the area and is assumed to be parabolic (Fig. 7c). The force balance between the bearings and the shaft is made following these assumptions: bodies are supposed rigid and the effects of the clearance and friction are neglected. The contact reaction forces are obtained by solving the force balance in the fixed frame considering the weight of shaft, the weight of the unbalanced mass and the centrifugal force generated by the unbalanced mass fixed at the end of the rotating shaft.

Then, a recursive approach is used to calculate the cumulative wear though time. In each iteration the contact pressure and area, the wear volume, the wear depth and the cumulative wear are computed.

To calculate the contact pressure and area, the contacts are supposed frictionless and under normal pressure only. The following Boussinesq's equation is used:

𝛿 𝑚𝑎𝑥 = 1-𝜈 2𝜋𝐺 𝑒𝑞 ∬ 𝑃(𝑥,𝑧) √𝑥 2 +𝑧² 𝑑𝑥 𝑑𝑧 ( 4
)
where δmax is the maximum elastic strain, ν the Poisson's ratio and Geq equivalent shear modulus. The contact pressure and area are obtained by solving by dichotomy the system of equations. This numerical resolution is used because of the non-linearities of the equations.

The wear volume ΔV is calculated using the Archard's law:

𝛥𝑉 = 𝐾 𝑅 𝛥𝑑 ( 5
)
where K is the wear coefficient, R the magnitude of normal reaction force between the shaft and the bearings and Δd is the sliding distance.

The wear depth for one revolution is computed from the wear distribution field h(x,z) (represented in red in Fig. 7):

ℎ(𝑥, 𝑧) = ℎ 𝑚𝑎𝑥 (1 - 𝑥 𝐿 𝑐 ) cos ( 𝜋 𝑧 2𝑟 𝑒𝑞 𝜑(𝑥) ) (6)
𝜑(𝑥) = ± sin -1 ( √ 1 - 𝑥 𝐿 𝑐 sin(𝜑 0 )) ( 7
)
where hmax is the maximum wear depth, Lc the contact length, req the equivalent radius, φ(x) the half aperture angle of the contact surface and φ0 the maximum half aperture angle. By double integrating Eq. ( 6) and knowing the wear volume, the maximum wear depth hmax can be computed.

To obtain the cumulative wear depth corresponding to each time step ti, the wear depth for one revolution is multiplied by the number of revolutions done during the time step.

Finally, the geometry is updated to take into account the evolution of the contact area and the evolution of the clearance between the shaft and the bearings that increases with wear.

For the stable wear stage, the wear coefficient was supposed constant as done in other works [START_REF] Schmidt | Transient wear simulation based on three-dimensional finite element analysis for a dry running tilted shaftbushing bearing[END_REF][START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF]. Calibration tests were carried out to determine the wear coefficient according to the operating conditions used in this work. The mean value of the wear coefficients of the calibration tests was then used. The obtained value was equal to 9.955 E-11 mm²/N. By using the calibrated model, the wear depth can be predicted for a given operating condition.

The variability of the input parameters (unbalanced mass, rotational speed, position of the mass, wear coefficient) was taken into account in uncertainty analysis. This was done using the Monte Carlo's method and allowed to derive the confidence intervals of the predicted wear depth. The assumed variations of the input parameters used for Monte Carlo's method are the same as those in [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF].

Digital twin fault diagnosis

The aim of the implemented DTs is to detect an unexpected behaviour of the wear in stable stage. An unexpected behaviour is detected when control limits are exceeded. In this work, the fault diagnosis monitors the wear depth and its rate of change. For each DT, the scenarios of detection were identified and explained.

Fault diagnosis rules for data-based DT

For the DT using the data-based model, the wear depth, the mean wear rate and the standard deviation of the wear rate are followed. The wear depth is systematically compared against a fixed threshold. The mean and the standard deviation of the wear rate are compared against the UCL and the LCL defined by the SPC. The DT detects a change of behaviour when at least one of these three parameters exceeds the control limits. In this work, a fault scenario is defined as the excess of one or more control limits.

In order to identify all possible scenarios, a crosstabulation matrix was build (Table 5). For example, if the wear depth threshold and the control limits of the mean wear rate are exceeded, it corresponds to the scenario D2. A detailed description of these scenarios and the risk of incident that they occur are presented in Table 6. The level of risk that an incident appears is proposed from the authors' appreciation.

D7

All

Unlikely situation *** (3/3) Table 6. Details of scenarios of detection for the data-based DT. The risk of incident is noted from 1 star (low risk) to 3 stars (high risk).

Fault diagnosis rules for physic-based DT

For the DT using the physic-based model, the wear depth and the mean wear rate are monitored. The wear depth is compared against a threshold and to confidence interval. The mean wear rate is compared against to confidence interval too. The confidence intervals were calculated from the numerical model. The DT detects a change of behaviour when at least one of the control limits is exceeded. Possible scenarios are listed in Table 7. A detailed description of these scenarios and the risk of incident that they occur are presented in Table 8. 8. Details of scenarios of detection for the physic-based DT. The risk of incident is noted from 1 star (low risk) to 3 stars (high risk).

Fault diagnosis test cases

The rules for fault detection described previously were tested in both DT under two situations: normal and abnormal conditions of wear. The aim is to highlight the capacity of the two DTs to detect a change in behaviour that could affect the integrity of the monitored mechanical system.

The normal conditions of wear represent the evolution of the system under the operating conditions described in Section 2.1.1. An unexpected change of behaviour, was designed to obtain the abnormal conditions of wear. In this work, this situation was simulated by altering the wear rate.

In this article, the wear depth threshold was fixed to 10 µm for both models. In abnormal conditions of wear, the wear rate was accelerated by adding 6 µm/h after 6h20 of use to the measured displacement. This acceleration represents an increase loadings of about 3 times the current loads or a deterioration of the contact conditions (i.e. presence of a third body, increase of the wear coefficient).

Results

Data-based DT

Fig. 8 shows the results of the wear evolution monitoring and the fault diagnosis with the DT using the SPC.

The evolution of the wear under normal conditions is presented in Fig. 8a-c. The mean (Fig. 8b) and the standard deviation (Fig. 8c) of the wear rate stay inside the control limits. The threshold of the wear depth (red line in Fig. 8a) is exceeded after 10h of use. This situation correspond to scenario D1 explained in Section 2.3.1.

The evolution of the wear under abnormal conditions is presented in Fig. 8d-f. The change of behaviour is clearly visible after 6h20 of use in contrast with the normal condition of wear. The control chart of the mean wear rate (Fig. 8e) is the first to detect the fault at about 6h30. This situation corresponds to scenario D4a (Section 2.3.1). Then the threshold of the wear depth (red line in Fig. 8d) is also exceeded at about 7h20. Therefore, at this moment, scenario D4a evolves to scenario D2 (Section 2.3.1). On the control chart of the standard deviation of the wear rate (Fig. 8f), a change of behaviour is visible but not sufficient to exceed the UCL.

Physic-based DT

The results of the wear evolution monitoring and the fault diagnosis with the numerical model are presented in Fig. 9.

In normal conditions of wear (Fig. 9a-b), all recorded data stay inside the confidence intervals. The threshold of the wear depth (red line in Fig. 9a) is exceeded after 10h of use. This situation corresponds to scenario P1 explained in Section 2.3.2.

The evolution of the wear under abnormal conditions is presented in Fig. 9c-d. The mean of the wear rate (Fig. 9d) exceeds the control limits at 7h. This situation corresponds to scenario P6 (Section 2.3.2). Then the wear depth (Fig. 9c) exceeds the confidence interval and the threshold simultaneously at about 7h20. At this moment, scenario P6 evolves to scenario P7 (Section 2.3.2).

Discussion

Comparison of detection capacity

Both DTs detect wear in normal and abnormal conditions. The monitoring of the wear rate permits to detect faults earlier than with only the monitoring of the wear depth.

The data-based DT allows to detect faults earlier because the size of the control limits is smaller (except at the very beginning) than those of the numerical model. In order to illustrate that, an overlapping of both control limits for the trial OC 1 (Table 2) is presented in Fig. 10. It can be noticed, that at 5h, e2 is 7 times larger than e1.

The SPC technique is well adapted to the implementation of DT if the target value of the monitored variable is known and if its evolution is linear. If the monitored variable is not linear (i.e. varying speed), the use of higher order statistical tools are necessary. Moreover, to use the SPC a large amount of data is necessary compared to the numerical model.

In this work, the monitored variable can be estimated experimentally, which is not the general case. There exist cases in which only the operating conditions are known. In those cases, the data-based DT cannot be used; in contrast, the physic-based can be, thanks to its predictive capabilities.

Comparison with other DTs

In the literature, several DTs of wear based on data analysis or on physical laws were reported. Some of these works are summarized in Table 9. They are compared regarding the cost of implementation and calibration, the computation time, the capability to be re-use and the capability to detect and predict wear. For the comparison, two data-based DTs and two physicbased DTs were analysed. The data-based DTs analysed rely on machine learning [START_REF] Desai | Lifetime Prediction Using a Tribology-Aware, Deep Learning-Based Digital Twin of Ball Bearing-Like Tribosystems in Oil and Gas[END_REF] and on SPC (as presented in this work). The physic-based DTs analysed rely on FEM [START_REF] Mukras | Numerical integration schemes and parallel computation for wear prediction using finite element method[END_REF] and on numerical model (as presented in this work).

From this comparison, it can be concluded that the data-based DTs are easier to implement and execute.

Physic-based DTs are able to detect and predict wear behaviour. If an analytical model is too hard to implement, then a finite elements model is an alternative.

Limitations

The limitations linked to the choices and assumptions to implement the DTs models are discussed in this section.

Only dry contacts were considered in this work to represent the worst operating conditions for a bush bearing system. Indeed, on helicopters, the loss of lubrication leads to run in dry contact conditions that alter the mechanical systems.

In the case studied in this work, the run-out due to clearance is very small regarding the length of the bearings mounted on the frame so the effect of clearance were neglected. About the effects of friction not taken into account, another work [START_REF] Mattei | Wear Simulation of Metal-on-Metal Hip Replacements With Frictional Contact[END_REF] shows that friction influences only the location of wear and not the pressure distribution.

To instrument the test bench, the choice was made to use a displacement sensor and a FFT analysis to collect data and nourish the DTs model. Another possibility was to use speed and force sensors to directly have the two input parameters (rotational speed and unbalanced mass) of the DT models implemented.

Conclusion

In this work, two DTs for monitoring wear were developed and compared. The aim of this work was to demonstrate the proof of concept of DT implemented diagnostic tool for bush bearing wear. The studied physical system is a bush bearing mounted on a test bench. This test bench was instrumented with temperature, acceleration and displacement sensors to monitor wear effects. Both DTs are fed by data recorded on a test bench to monitor the wear depth, the mean wear rate and the standard deviation of the wear rate. To detect wear faults, different scenarios of detection were identified, explained and applied to the two DTs.

The first DT is a data-based model using SPC and the second DT is a physic-based model relying on a numerical model. The two implemented DTs detect abnormal wear but not in the same delay; the databased allows detects abnormal wear faster. The databased DT implemented in this work is easier to develop but it is limited to monitoring variable with linear behaviour. If the monitored variable is not linear, the use of higher order statistical tools are necessary (i.e. machine learning). The physic-based DT can predict wear evolution and detect abnormal behaviour by means of confidence intervals. The physic-based DT needs more time to be developed and it is adapted only to contact surfaces defined by analytical geometries.

In future work, the models could be completed to be more realistic by limiting hypotheses taken to implement the model, to predict more variables (temperature, shaft displacement, etc.), to work with more wear cases (running-in, varying of speed, etc.) or even to use more complex input data (vibration signals with harmonics). This work is a contribution for new wear diagnostic tools.

Data-based models

Physic-based models

Machine learning SPC (from this work) FEM Numerical model (from this work)

Studied case

Remaining useful lifetime of ball bearingtype contacts [START_REF] Desai | Lifetime Prediction Using a Tribology-Aware, Deep Learning-Based Digital Twin of Ball Bearing-Like Tribosystems in Oil and Gas[END_REF] Wear of bearings under dynamic loading

Wear on an oscillatory pin joint [START_REF] Mukras | Numerical integration schemes and parallel computation for wear prediction using finite element method[END_REF] Wear of bearings under dynamic loading [START_REF] Regis | Numerical model to predict wear of dynamically loaded plain bearings[END_REF]

Costs of implementation and calibration

Training time of about 3 min [START_REF] Desai | Lifetime Prediction Using a Tribology-Aware, Deep Learning-Based Digital Twin of Ball Bearing-Like Tribosystems in Oil and Gas[END_REF] Only need target values of the monitored variables

Well adapted for complex geometries, need to calibrate wear coefficient on experimental data [START_REF] Mukras | Numerical integration schemes and parallel computation for wear prediction using finite element method[END_REF] Adapted only for analytical geometries, need to calibrate wear coefficient on experimental data [

Fig. 1 .

 1 Fig. 1. Global architecture of DTs implemented.

Fig. 2 .

 2 Fig. 2. Experimental test bench: (a) test bench photo; (b) test bench. (1) unbalanced mass; (2) accelerometers; (3) frame; (4) double universal joint; (5) motor; (6) shaft; (7) bearings; (8) thermal sensors; (9) displacement sensor.

  Mean of the peak-to-peak of the recorded signal Table

Fig. 3 .

 3 Fig. 3. Picture of (a) worn shaft and (b) worn bearing with the measurement paths.

Fig. 4 .

 4 Fig. 4. Detection of the operating configurations using neural network of Matlab: (a) flow chart used to identify operating configurations; (b) shaft displacement signal; (c) Fast Fourier Transformation (FFT) of the shaft displacement signal.

Fig. 5 .

 5 Fig. 5. Wear depth from displacement sensor measurement.

Fig. 6 .

 6 Fig. 6. Pseudo algorithm of the wear numerical model.

Fig. 7 .

 7 Fig. 7. Details of the wear distribution and the contact pressure field between the shaft and the bearings: (a) Plane with normal x; (b) Plane with normal z; (b) Plane with normal y.

Table 7 .

 7 Scenarios of detection for the physic-based DT.

Fig. 8 .

 8 Fig. 8. Detection of wear using the data-based DT: left column is the detection in normal use for (a) the wear depth, (b) the mean wear rate and (c) the standard deviation of the wear rate; right column is the detection in abnormal use for (d) the wear depth, (e) the mean wear rate and (f) the standard deviation (SD) of the wear rate. On each graph, the operating conditions (OC) are marked by vertical dotted grey lines.

Fig. 9 .

 9 Fig. 9. Detection of wear using the physic-based DT: left column is the detection in normal conditions for (a) the wear depth and (b) the mean wear rate; right column is the detection in abnormal conditions for (c) the wear depth and (d) the mean wear rate. On each graph, the operating conditions (OC) are marked by vertical dotted grey lines.

Fig. 10 .

 10 Fig. 10. Overlapping of control limits of both DTs implemented.

Table 1 .

 1 Material properties of the shaft and the bearings.

		Shaft		Bearing
	Material	Steel 115Cr2	Bronze CuSn10
	Young's modulus	206 GPa	65 GPa
	Roughness (Ra)	0.3 -0.6	0.3
	Hardness	54 -56 HRC	25 HB
		OC 1	OC 2	OC 3	OC 4
	Start time	0h	5h30	8h20	10h
	End time	5h30	8h20	10h	13h
	Unbalanced mass	50 g	75 g	50 g	75 g

Table 2 .

 2 Operating conditions (OC) of the experiments.

Table 4 .

 4 Target values for the mean and the standard deviation of the wear rate for unbalanced mass of 50 g and 75 g.

		µ (µm/h)	σ (µm/h)
	m = 50 g	0.933	0.320
	m = 75 g	1.412	0.389

Table 5 .

 5 Scenarios of detection for the data-based DT.

		Wear depth threshold	Limits of the mean wear rate	Limits of the the wear rate standard deviation of	All
	Wear depth threshold	D1	D2	D3	NA
	Limits of
	the mean	D2	D4a, D4b	D5a, D5b	NA
	wear rate
	Limits of
	the standard	D3	D5a, D5b	D6	NA
	deviation of

Table 9 .

 9 Comparison of wear DT models.

	21]

Acknowledgements

This work was supported by Airbus Helicopters, Aix-Marseille University and ANRT (Association Nationale pour la Recherche et la Technologie) by the grant n° 2019/1242. The experimental devices were founded by: European Community, French Ministry of Research and Education and Aix-Marseille Conurbation Community.


Moriya Interaction M Debbichi
email: mourad.debbichi@fsm.rnu.tn

H Saïd

H Garbouj

S El Hog

Van An Dinh

A new ternary Pentagonal Monolayer based on Bi with large intrinsic Dzyaloshinskii

We systematically investigate the structural, electronic, and magnetic properties of a new pentagonal CoBiS monolayer using first-principles and Monte Carlo simulations. We find that Penta-CoBiS is stable mechanically, dynamically, and thermally and is an antiferromagnetic semiconductor with an indirect band gap of 0.5 eV with HSE functional. In addition, the band-gap increased by applying in-plane biaxial strain. We further show that this monolayer has an in-plane easy axis and possesses large intrinsic Dzyaloshinskii-Moriya interaction (DMI) because of the broken inversion symmetry, and strong spin-orbit coupling (SOC) originated from the Bi atoms. Moreover, the Néel temperature is also predicted using Monte Carlo simulations. An out-of-plane magnetic field B is then applied to compensate the in-plane anisotropy. It is found that for B = 1.72 T the spins are fully polarized to the out-of-plane direction. Our results demonstrate that Penta-CoBiS monolayer may find numerous applications in flexible spintronics and electronics.

I. INTRODUCTION

Two-dimensional (2D) materials with trigonal (2H-type hexagonal) prismatic symmetry have received a lot of attention due to their outstanding physicochemical properties and their varieties of applications from sensing to data storage [START_REF] Debbichi | Tuning the magnetic and electronic properties of monolayer chromium tritelluride through strain engineering[END_REF][START_REF] Chae | A new family of two-dimensional crystals: Open-framework T3X(T = C, Si, Ge, Sn; X = O, S, Se, Te) compounds with tetrahedral bonding[END_REF][START_REF] Liu | Screening magnetic two-dimensional atomic crystals with nontrivial electronic topology[END_REF][START_REF] Wang | Large magnetic anisotropy and its strain modulation in two-dimensional intrinsic ferromagnetic monolayer RuO2 and OsO2[END_REF]. Despite those efforts, the database size of the 2D materials is still limited compared to that of the three-dimensional.

The prediction of new 2D prototypes for a number of novel materials is a strategy to tune their properties. In order to discover novel 2D materials, it has been demonstrated that it is important to explore structural property relationships by changing the building blocks of a 2D material from hexagons to pentagons [START_REF] Houlong | From pentagonal geometries to two-dimensional materials[END_REF][START_REF] Lan | Penta-PdX2(X = S, Se, Te) monolayers: promising anisotropic thermoelectric materials[END_REF]. This new allotrope, composed of pentagonal building blocks was proposed from first-principles calculations and confirmed experimentally by exfoliation method from bulk crystal [START_REF] Lan | Penta-PdX2(X = S, Se, Te) monolayers: promising anisotropic thermoelectric materials[END_REF][START_REF] Zhang | Penta-graphene: A new carbon allotrope[END_REF][START_REF] Oyedele | Pdse2: Pentagonal two-dimensional layers with high air stability for electronics[END_REF].

Since the prediction of penta-graphene in 2015 [START_REF] Zhang | Penta-graphene: A new carbon allotrope[END_REF] , several new 2D nanosheets constructed exclusively from pentagonal units have gained a lot of attention. Yang et al. predicted by first-principles calculations that penta-PdX 2 (X = S, Se, Te) monolayer is a potential anisotropic thermoelectric material making it a candidate for field-effect transistors applications [START_REF] Lan | Penta-PdX2(X = S, Se, Te) monolayers: promising anisotropic thermoelectric materials[END_REF][START_REF] Oyedele | Pdse2: Pentagonal two-dimensional layers with high air stability for electronics[END_REF]. N-type few-layers penta-GeP 2 can be useful in photocatalyzed water splitting and in CO 2 splitting to Co [START_REF] Shojaei | Electronic structure and photocatalytic band offset of few-layer GeP2[END_REF]. SnX 2 (X = S, Se, Te) is identified as a promising candidate for the application of low-power-consuming electronic devices [START_REF] Ma | Room temperature quantum spin hall states in two-dimensional crystals composed of pentagonal rings and their quantum wells[END_REF]. Penta-graphene was proposed from the first-principles calculation as a promising anode material for Li/Na-ion battery [START_REF] Xiao | Penta-graphene: A promising anode material as the Li/Na-ion battery with both extremely high theoretical capacity and fast charge/discharge rate[END_REF]. Other new 2D monolayer pentagonal structures with the MX 2 formula were also predicted and studied by density functional theory (DFT) calculations such as SiX 2 (X = S, N) [START_REF] Liu | Disparate strain dependent thermal conductivity of two-dimensional penta-structures[END_REF],CoS 2 [START_REF] Liu | Single-layer antiferromagnetic semiconductor CoS2 with pentagonal structure[END_REF], MX 2 (M = Ni, Pd and Pt; X = P and As) [START_REF] Wei | The thermoelectric properties of monolayer MAs2 (M = Ni, Pd and Pt) from first-principles calculations[END_REF][START_REF] Qian | Penta-MX2 (M = Ni, Pd and Pt; X = PandAs) monolayers: direct band-gap semiconductors with high carrier mobility[END_REF], MSe 2 (M =Cd, Zn and X =S, Se) [START_REF] Yang | Two unexplored two-dimensional MSe2(M = Cd, Zn) structures as the photocatalysts of water splitting and the enhancement of their performances by strain[END_REF]17].

Single-layer pentagonal with the formulas MX 2 has an inversion center of symmetry, but there is no such symmetry in single-layer MXY (X =Y), opening up possibilities for novel properties. Recently it has been demonstrated experimentally and theoretically that many noncentrosymmetric 2D monolayer materials exhibit considerable piezoelectricity making them very promising for sensor and energy conversion system applications [START_REF] Ahammed | Ultrahigh out-of-plane piezoelectricity meets giant rashba effect in 2d janus monolayers and bilayers of group iv transition-metal trichalcogenides[END_REF]. It has been also predicted that 2D magnet without the inversion symmetry are good candidates for the formation of magnetic skyrmionic and chiral domain wall, making them promising as energyefficient logic and memory devices [START_REF] Liang | Electrically switchable rashba-type dzyaloshinskii-moriya interaction and skyrmion in two-dimensional magnetoelectric multiferroics[END_REF][START_REF] Zhang | Electronic structure, magnetic anisotropy and dzyaloshinskii-moriya interaction in janus Cr2I3X3(X = Br, Cl) bilayers[END_REF]. Their use in skyrmion-based spintronics has been facilitated by the presence of the Dzyaloshinskii-Moriya interaction (DMI), which is responsible for the creation of exotic chiral magnetic states [START_REF] Vedmedenko | Interlayer dzyaloshinskii-moriya interactions[END_REF]. This behavior has been discovered in a variety of 2D materials systems such as, MnXY (X,Y= S, Se, Te, X = Y)monolayers [START_REF] Liang | Very large dzyaloshinskii-moriya interaction in two-dimensional janus manganese dichalcogenides and its application to realize skyrmion states[END_REF][START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF] , Cr 2 I 3 X 3 (X = Br, Cl) bilayers [START_REF] Zhang | Electronic structure, magnetic anisotropy and dzyaloshinskii-moriya interaction in janus Cr2I3X3(X = Br, Cl) bilayers[END_REF] and 2H-VSeX (X = S, Te) monolayers [START_REF] Qi | Tunable valley polarization, magnetic anisotropy and dzyaloshinskii-moriya interaction in two-dimensional intrinsic ferromagnetic janus 2h-vsex (x = s, te) monolayers[END_REF].

Most of the work focused on 2D pentagonal structures have so far been restricted to the binary system. However, by substituting one chalcogen atom (X) in penta-MX 2 monolayer by an isovalent atom it is plausible to stimulate novel physical phenomena, and some new properties may appear (DMI, valley polarization, Intrinsic ferromagnetism and spontaneous vertical electric dipole moment ), which have practical applications in multifunctional spintronic devices.

In this paper, by combing first-principles calculations and Monte Carlo simulation we have explored the stability, electronic structure, and magnetic properties of Penta-monolayer CoBiS.

II. COMPUTATIONAL DETAILS AND METHODS

Most of the calculations were performed on the basis of the density functional theory as implemented in QUANTUM-ESPRESSO ab-initio simulation package [START_REF] Giannozzi | Quantum espresso: a modular andopen-source software project for quantumsimulations of materials[END_REF] with projected augmented wave method (PAW). For the exchange-correlation part of the calculations, the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] form is employed with a kinetic cutoff energy of 60 Ryd. The first Brillouin zone (BZ) integration was sampled by 12×12× 1 Monkhorst-Pack k-point grid for unit cell. For the localization of the d electrons of cobalt atom, the Hubbard-model-based DFT+U corrective functional was employed as introduced by Cococcioni et al. [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u[END_REF] with the on-site Coulomb U = 3.32 eV [START_REF] Liu | Single-layer antiferromagnetic semiconductor CoS2 with pentagonal structure[END_REF]. Furthermore, the hybrid Heyd-Scuseria-Ernzerhof (HSE) functional is used to achieve more accurate results of the electronic properties. The monolayer was constructed by inserting a vacuum space of 18 Å along the z direction to avoid the interactions between two neighboring images. For the structure optimisation, we have used a tolerance for total-energy of 10 -7 Ryd and a total force in the unit cell less than 10 -4 Ryd/ Å. Car-parrinello molecular dynamics (CPMD) simulation as implemented in the QUANTUM-ESPRESSO package [START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum espresso[END_REF] was carried out by using the Nose-Hoover thermostat for a total time of 6 ps with a step of 1 fs to control the temperature. Phonon dispersion spectra are calculated by using Density-Functional Perturbation Theory (DFPT) [START_REF] Sohier | Density functional perturbation theory for gated twodimensional heterostructures: Theoretical developments and application to flexural phonons in graphene[END_REF] with an ultrasoft pseudopotentials (USPPs) and a 4×4× 1 q-point grid. The strain coefficients ε and the lattice parameter are related via the equation:

ε = [(a -a 0 )/a 0 ] × 100. ( 1
)
where a 0 and a are lattice constants of the unstrained and strained cell, respectively. The values of ε considered in our calculations are ranging from -4% to 4% with incremental step of 0.5%. The negative value of ε represents compression, while the positive one is the tension.

III. RESULTS AND DISCUSSION

A. Crystal structure and stability

The starting values of the lattice parameters and atomic positions used in our calculations are those of the pentagonal monolayer CoS 2 presented in the reference [START_REF] Liu | Single-layer antiferromagnetic semiconductor CoS2 with pentagonal structure[END_REF]. To form CoBiS, the inversion symmetry of CoS 2 is broken by partially replacing the S atoms connecting the Co atoms with Bi atoms. The blue, yellow and purple spheres represent Co, S, and Bi atoms, respectively.

Structural stability of of CoBiS is examined by it's formation energy expressed as:

E f = E(CoBiS) -N Bi E(Bi) -N S E(S) -N Co E(Co) N Bi + N S + N Co (2)
where E(CoBiS) is the total energy of the CoBiS monolayer while E(Bi), E(Co) and E(S) (N Bi , N Co and N S )are respectively, the chemical potential (total number) of Bi, Co and S atom in their standard bulk structures. We found that CoBiS have the E f of -133 meV/atom. Also, the cohesive energy (E coh ) is calculated by considering the total energy of the CoBiS and total energies of the isolated constituents atoms. The obtained value is -2.85 eV/atom.

To check the mechanical stability, we calculated the 2D elastic constants. For 2D system there are only four independent constants, noted in Voigt notation [START_REF] Zhao | Penta-BCN: A new ternary pentagonal monolayer with intrinsic piezoelectricity[END_REF] [START_REF] Ying | Screening magnetic two-dimensional atomic crystals with nontrivial electronic topology[END_REF], our calculated elastic constants show that CoBiS monolayer is mechanically stable. Furthermore, the 2D Young's moduli along the x and y directions (E x and E y ) are also evaluated according to the formulas: E x = (C [START_REF] Alyörük | Promising piezoelectric performance of single layer transition-metal dichalcogenides and dioxides[END_REF][START_REF] Deniz C ¸akır | Mechanical and thermal properties of h-MX2 (M = Cr[END_REF][START_REF] Wang | The hydrogen-induced structural stability and promising electronic properties of molybdenum and tungsten dinitride nanosheets: a first-principles study[END_REF], indicating that CoBiS monolayer is a flexible system.

The thermal stability of the CoBiS monolayer is investigated by the Car-parrinello molecular dynamics (CPMD) simulation using a 2×2×1 supercell. As shown in Fig. 2 By considering the stability of this monolayer (thermodynamic, mechanical, thermal, and dynamic), we can be affirmed that the synthesis of the penta-CoBiS monolayer will be realistic.

To explore the preferred magnetic ground state of the CoBiS monolayer, we considered three configurations, ferromagnetic (FM) and the antiferromagnetic (AFM) spin coupling, and nonmagnetic (NM). Through the total energy calculation, we found that the NM state was the most unstable by an energy difference of approximately 657.61 meV compared with the antiferromagnetic state. While, the calculated energy difference (∆E = E F M -E AF M ) demonstrated that the AFM state was more stable with an energy difference of 53. [START_REF] Yang | Significant dzyaloshinskii-moriya interaction at graphene-ferromagnet interfaces due to the rashba effect[END_REF] We next study the biaxial strain effects on magnetic ground state of CoBiS monolayer.

Indeed and for any 2D materials, tuning the magnetic and electronic properties with the application of external strain is very important because it undergoes structural modification during synthesis while growing on a substrate. Fig. 3 shows the energy differences of single-layer pentagonal CoBiS with the FM and AFM orderings under the biaxial strains. For the compressive strain, we find that the energy difference decreases with the increasing strains while keeping a positive sign, showing that the AFM ground state. On the other hand, during the tensile strain, the energy difference decreases with the applied strain until it becomes negative from ε=2.75% indicating the appearance of the magnetic transition from AFM to FM.

B. Electronic properties

In Fig. 4 is shown the band structure along the high symmetries points (Γ-X-M-Γ-Y-M) and density of states (DOS) of CoBiS in it's stable magnetic configuration. As clearly shown the system is an indirect semiconductor with a narrow band-gap of 0.4 eV for the GGA+U.

The valence band maximum (VBM) is located at the X point, whereas the conduction band minimum (CBM) is along Γ -X. The density partial (PDOS) analysis shows that the VBM of CoBiS monolayer is mainly contributed by 3d state of Co and 6p state of Bi, while the CBM is mainly from the mixed states 3d of Co and 3p of S and 6p of Bi. To obtain more reliable band structure, the hybrid functional HSE+U is employed and the corresponding band structure is also given in Fig. 4, showing a band gap of 0.55 eV.

When the spin-orbit coupling (SOC) is included there are no remarkable changes for the electron properties near the Fermi level as shown in Fig. S1, the pronounced change is observed especially at high symmetry points.

In addition to the band structure, other electronic properties such as the Work func-tion (Φ), Ionization Potential (I.P.), and electron affinity (χ) are also computed using the following relations:

Φ = E vac -E F er

I.P. = E vac -E V BM χ = E vac -E CBM ; (3)
where E vac , E F er , E V BM and E CBM refer to the vacuum potential, Fermi energy level, valence band maximum and the conduction band minimum, respectively. The calculated parameters from the electrostatic potential are Φ=5.02 eV, I.P.=2.38 eV and χ = 2.05 eV. The calculated work function is in good agreement with that of MoX(S,Se) [START_REF] Rawat | A comprehensive study on carrier mobility and artificial photosynthetic properties in group vib transition metal dichalcogenide monolayers[END_REF], while the other electronic parameters such as ionization potential and electronic affinity for the monolayers have not been systematically studied. The change of band-gap with applied biaxial strain ranging from -4%-2.75% is described in Fig. 5. As shown, the band-gap energy decreases monotonically with increasing compression strain, and she increases with tensile strain, while remaining the semiconducting behavior with a narrow band-gap. It is interesting to note that the biaxial strain tends to shift the positions of the bottom of the conduction band from Γ -X path to Γ point, as depicted in Fig. S2. As a consequence, the indirect gap transition Γ -X has been found for ε = 2.5%.

C. Magnetic properties

As indicated in the last section the CoBiS monolayer is antiferromagnetic with a magnetic moments mainly from the Co atoms (1.72 and 1.96 µ B /Co with GGA+U and HSE+U, respectively). However the contribution of S and Bi atoms is very small. Fig. 6 shows the evolution of the absolute magnetic moments of the Co atoms as functions of applied strains. The origin of this enhancement was explored by the Bader charge calculations (Fig. 6). We see that magnetic moment correlate strongly with the charge transfer from the Co atoms to the S atoms. Magnetic anisotropy energy (MAE) is a required parameter for magnetic information storage. This parameter is obtained by computing the energy difference between two perpendicular magnetic phases of the Co atoms in the unit cell. MAE is calculated by using the force theorem as implemented in the QUANTUM ESPRESSO package [START_REF] Li | Out-versus in-plane magnetic anisotropy of free Fe and Co nanocrystals: Tight-binding and first-principles studies[END_REF]. This method consists of the following steps. Firstly, a self-consistent (SCF) calculation with scalar-relativistic pseudo-potential ; (without SOC)is performed to obtain the charge density and the spin moment, then a non-SCF non-collinear calculations is done by including the SOC with a full relativistic pseudo-potential and much denser k-mesh of 20 × 20 × 1.

Then MAE is evaluated according to the following equation:

M AE = E [uvw] -E min (4)
with E min the energy of the easy axis magnetization and E [uvw] is the energy for a magnetization along the [uvw] crystallographic direction. Our calculations show that b-axis is the easy axis from the antiferromagnetic configuration, showing an in-plane magnetization feature.

The calculated MAE is found to be 0.554 meV/Co atom. This MAE value is comparatively larger than those of some 2D system such as, CrX 3 ,(X= Br, Cl, Te and F)(0.0315-0.185 meV per Cr atom) [START_REF] Debbichi | Tuning the magnetic and electronic properties of monolayer chromium tritelluride through strain engineering[END_REF][START_REF] Bhattacharyya | Density functional theory study of defect induced ferromagnetism and half-metallicity in CaI2 based monolayer for spintronics applications[END_REF], CrX, (X=P,As) (0.217-0.404 meV per Cr atom) [START_REF] Ma | Intrinsic ferromagnetism with high temperature, strong anisotropy and controllable magnetization in the CrX(X = P, As) monolayer[END_REF] and in the same order of magnitude as CrI 3 (0.685 meV per Cr atom) [START_REF] Zhang | Robust intrinsic ferromagnetism and half semiconductivity in stable two-dimensional single-layer chromium trihalides[END_REF]. Indeed, compared to other 2D systems with 4d or 5d electrons which possesses a strong spin-orbit coupling (SOC) interaction, this MAE value remains a low value [START_REF] Wang | Large magnetic anisotropy and its strain modulation in two-dimensional intrinsic ferromagnetic monolayer RuO2 and OsO2[END_REF][START_REF] Bhattacharyya | Density functional theory study of defect induced ferromagnetism and half-metallicity in CaI2 based monolayer for spintronics applications[END_REF][START_REF] Hang | Prediction of high-temperature chern insulator with half-metallic edge states in asymmetry-functionalized stanene[END_REF].

The low-lying excited states of a magnetic material can be generally described by the Heisenberg Hamiltonian [START_REF] Debbichi | Tuning the magnetic and electronic properties of monolayer chromium tritelluride through strain engineering[END_REF]:

H = H 0 - <i,j> J ij S i S j ( 5
)
where H 0 is the nonmagnetic part of the Hamiltonian, S i and S j are the vector spin operators on sites i and j, and J ij is the spin coupling constant between i and S j . Correspondingly, a positive (negative) J ij indicates that the coupling is ferromagnetic (antiferromagnetic). For the spin state in a 2× 2× 1 supercell as shown in Fig. 7, and according to this model, the exchange coupling parameters can verify the following relations:

E (a) = E 0 -(8J Co-S-Co + 8J Co-Bi-Co ) (6)
E (b) = E 0 -(-8J Co-S-Co -8J Co-Bi-Co ) (7)
E (c) = E 0 + 4J Co-Bi-Co (8)
By combining the three equations, we obtained:

J Co-Bi-Co = 2E (c) -E (b) -E (a) 8
and

J Co-S-Co = E (b) -E (a) 16 -J Co-Bi-Co
The obtained exchange coupling constant of J Co-S-Co and J Co-Bi-Co are -40.79 meV et 27.43 meV, respectively. These values indicate that the Co-S-Co interaction is antiferromagnetic while that of Co-Bi-Co is ferromagnetic. Consequently, we can state that the magnetic order (AFM) of the CoBiS monolayer is mainly due to the Co-S-Co superexchange. This result follows the empirical Goodenough-Kanamori-Anderson rule [START_REF] Zhang | Robust intrinsic ferromagnetism and half semiconductivity in stable two-dimensional single-layer chromium trihalides[END_REF][START_REF] Goodenough | Theory of the role of covalence in the perovskite-type manganites[END_REF][START_REF] Kanamori | Superexchange interaction and symmetry properties of electron orbitals[END_REF]. According to this law, for a superexchange interaction angle between two magnetic ions greater than 90°, AFM coupling is favored, whereas for superexchange interaction close to 90°, FM coupling is preferred. The bond-angle, α 1 of Co-S-Co, as shown in TableI, is higher than 90°(∼ 130°), which means Co-d (dz 2 ) and S-p (pz) orbitals have a significant overlap integral (see Fig. S3). The critical transition temperature of a material is an important parameter to defines the possibility of its use in a spintronic device. In the mean field approximation and for an antiferromagnetic configuration, the critical transition temperature T N can be expressed as [START_REF] Liu | Screening magnetic two-dimensional atomic crystals with nontrivial electronic topology[END_REF]:

E F M -E AF M = N 3 2 k B T N (9)
where N is the number of magnetic atoms in the unit cell and k B is the Boltzmann constant.

Therefore, the calculated Néel temperature T N value is 206.59 K. This value is four times higher than that of CoS 2 [START_REF] Liu | Single-layer antiferromagnetic semiconductor CoS2 with pentagonal structure[END_REF].

The theories of the mean field neglect correlation and fluctuation, thus overestimating the T N [START_REF] Bhattacharyya | Density functional theory study of defect induced ferromagnetism and half-metallicity in CaI2 based monolayer for spintronics applications[END_REF][START_REF] Bhattacharyya | Ferromagnetism in magnesium chloride monolayer with an unusually large spin-up gap[END_REF]. For more realistic results, Monte Carlo (MC) simulations based on the Heisenberg Hamiltonian. In this step the Hamiltonian, which includes only the exchange interaction is adopted [START_REF] Hog | Stability and phase transition of skyrmion crystals generated by dzyaloshinskii-moriya interaction[END_REF]. A supercell of 70×70 with periodic boundary conditions was created for the MC simulation, for more detail see supplementary material. In Fig. 8 is plotted the heat capacity C v versus temperature (K) of the CoBiS system. From the plot, we noted that there is a transition at T N =124 K, and 204 K from antiferromagnetic to paramagnetic phase for the Heisenberg and XY model, respectively. Dzyaloshinskii-Moriya interaction (DMI) may appear in CoBiS due to the broken of inversion symmetry. To evaluate DMI we have used the chirality-dependent total energy difference approach as implemented in VASP [START_REF] Qi | Tunable valley polarization, magnetic anisotropy and dzyaloshinskii-moriya interaction in two-dimensional intrinsic ferromagnetic janus 2h-vsex (x = s, te) monolayers[END_REF] by taking into account the spin-orbit coupling (SOC). The DMI coefficient D which represents the antisymmetric exchange effects can be expressed as: D=

√ 3d rt F
, where t F and r represent the thickness of the magnetic atomic layer and the distance between the nearest neighbors of magnetic element, respectively. d is the norm of the DMI vector d ij , given by:

d ij = d (z × u ij ) + d ⊥ z.
Here u ij is the unit vector from site i to site j and z is the out-of-plane unit vector. d and d ⊥ are the in-plane and out-of-plane components of d ij , respectively. d ⊥ is negligible [START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF][START_REF] Qi | Tunable valley polarization, magnetic anisotropy and dzyaloshinskii-moriya interaction in two-dimensional intrinsic ferromagnetic janus 2h-vsex (x = s, te) monolayers[END_REF] and we can note d = d = (E CW -E ACW )/n , where E CW and E ACW are the energies of clockwise and anticlockwise spin configurations(Fig. 9), respectively, and n is a number related to the atomic and magnetic chirality of the system. Using this method and for the nearest neighbors of Co atoms, we find the magnitude of the d parameter to be 4.53 meV. This value is comparable to those of MnSTe and MnSeTe and higher to VSeTe and MnSSe monolayers [START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF]. To elucidate the origin of the large DMI, we have extracted the energy of SOC from opposite chiralities in the monolayer, ∆E SO =2.03 meV. So we can state that the large DMI in CoBiS is associated with a large variation ∆E SO . This phenomenon has been identified

for the FM/HM heterostructures and some other 2D materials such as MnXY(X,Y= S, Se, Te, X = Y ), (CrXTe, X = S, Se) [START_REF] Liang | Very large dzyaloshinskii-moriya interaction in two-dimensional janus manganese dichalcogenides and its application to realize skyrmion states[END_REF][START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF][START_REF] Yang | Anatomy of dzyaloshinskii-moriya interaction at Co/Pt interfaces[END_REF]. Indeed, the heavy atom Bi (5d) functions as active spin-orbit sites to induce the spin-orbit dispersion needed for the DMI [START_REF] Liang | Very large dzyaloshinskii-moriya interaction in two-dimensional janus manganese dichalcogenides and its application to realize skyrmion states[END_REF][START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF]. The micromagnetic DMI coefficient D, which means the average DMI in the CoBiS monolayer is 5.49 mJ/m 2 , which is higher to that of Cr 2 I 3 X 3 (X = Br, Cl) bilayers [START_REF] Zhang | Electronic structure, magnetic anisotropy and dzyaloshinskii-moriya interaction in janus Cr2I3X3(X = Br, Cl) bilayers[END_REF] and in the same order of magnitude to those of MnSeTe monolayer [START_REF] Yuan | Intrinsic skyrmions in monolayer janus magnets[END_REF]and graphene/Co interface [START_REF] Yang | Significant dzyaloshinskii-moriya interaction at graphene-ferromagnet interfaces due to the rashba effect[END_REF].

CW ACW

FIG. 9: The spin configurations used to calculate the DMI parameter.

To inspect the effect the magnetocrystalline anisotropy, DMI and the magnetic field, we have adapted a Monte Carlo simulation based on this Hamiltonian [START_REF] Hog | Theory of magnons in spin systems with dzyaloshinskii-moriya interaction[END_REF]:

H = i,j -J i,j S i . S j -D i,j δr( S i × S j ) - i k(S y i ) 2 -B i S z i (10)
where δr = (r i -r j )

|r i -r j | u is a unitary vector. Here, J i,j is the Heisenberg exchange interaction between S i and S j , D is DMI coefficient, k is the anisotropic exchange coefficient, and B represents the uniform external magnetic field.

The DMI favors a noncollinear alignment among the spins, while the anisotropy determines the desired direction for the spin. The competition between these two-parameter and the exchange interaction determines the characteristics of the noncollinear states and tends to minimize the total energy (FigS.4).

The ground state together with heat capacity for B= 0 T are shown in Fig. 10, we obtain a transverse conical spin-spiral structure with an antiferromagnetic coupling. The same phenomena is also observed in Mn/W(110) [START_REF] Yoshida | Conical spin-spiral state in an ultrathin film driven by higher-order spin interactions[END_REF]49]. The extracted transition temperature is about 136 K. In fact, T N increased by about 10%, indicating that DMI have a limited impacts on the temperature dependence of magnetization.

IV. CONCLUSION

On the basis of first-principles DFT calculations, we have systematically studied the stability, electronic and magnetic properties of a new ternary pentagonal monolayer CoBiS.

We have found that, penta-CoBiS is an antiferromagnetic and stable (thermodynamically, mechanically, dynamically, and thermally) 2D crystal. Our electronic structure calculations, on the basis of HSE hybrid functional, show that CoBiS is semiconductor with a narrow indirect band-gap of 0.5 eV. The variation of band gap is found to be very sensitive to strain and can be tuned by the applied biaxial strain. Furthermore, our results also suggest that this monolayer exhibit significant intrinsic DMI and in-plane magnetic anisotropy. By MC simulations based on the Heisenberg model we find that CoBiS monolayer has a T N of 136 K. For a magnetic field parallel to the z-axis less than 1.72 T, the magnetic moments are tilted, forming a conical magnetic state. However, all the spin becomes completely aligned along the z-axis for B = 1.72 T. Ou results indicate that CoBiS can be considered as a

Fig. 1

 1 shows the top and side views of the relaxed 2×2×1 supercell of CoBiS monolayer. It has a monoclinic lattice symmetry and possesses the space group of C2(space group N o .5). Each unit cell contains 6 atoms (two from each element). Each Co atom is surrounded by tow S atoms and two Bi atoms, while each Bi(S) atom is three-fold coordinated with two Co atoms and one S(Bi) atom, forming an intriguing pentagonal ring network known as the Cairo pentagonal tiling. The structural parameters of the optimized structures are summarized in TableI. By comparing the structural parameters to that of CoS 2[START_REF] Liu | Single-layer antiferromagnetic semiconductor CoS2 with pentagonal structure[END_REF], we can state that the Bi atoms have a greater impact on the Co atoms than the S atom.

FIG. 1 :

 1 FIG. 1: (a) Top and (b) side views of the optimized geometry of penta-CoBiS monolayer.

  (a), after atomic relaxation there is no appreciable energy change throughout the simulation and the structure recovers it's initial configuration at 300 K. To examine the dynamic stability we have calculated phonon dispersion as shown in Fig.2(b). The absence of imaginary phonon modes confirms the dynamical stability of the CoBiS monolayer.

FIG. 2 :

 2 FIG. 2: (a) The energy evolution of CoBiS monolayer during CPMD of 5600 ps at 300 K. The insets are the initial and final atomic structure. (b)phonon spectrum calculated by DFPT.

FIG. 3 :

 3 FIG. 3:The evolution of the energy difference between FM and AFM configurations with the applied strain.

FIG. 4 :

 4 FIG.4: Band structure and density of states.

FIG. 5 :

 5 FIG.5: Variation of band gap energy in CoBiS with application of biaxial strain (-4% to 2.5%).

FIG. 6 :

 6 FIG. 6: The magnetic moment of the Co atom and Bader charge transfer from the Co atoms the S atoms as a function of strain.

FIG. 7 :

 7 FIG. 7: Different spin configurations used for the calculation of the exchange parameters. The arrows indicate the spin orientation.

FIG. 8 :

 8 FIG. 8: (Top) Heat capacity C v versus Temperature using the Heisenberg and XY models of CoBiS monolayer; (Bottom) Ground State spin configuration for the Heisenberg model(b), and XY model (c).

FIG. 10 :

 10 FIG. 10: (a)Monte Carlo simulations on the heat capacity C v using the Heisenberg model and (b) the spin ground state for B = 0 T.

FIG. 11 :

 11 FIG. 11: Ground state spin configurations as a function of the external magnetic field for CoBiS monolayer. (a) B= 0.86 T, and for (b)and (c) B=1.72 T.

  2 11 -C 2 12 )/C 11 and E y = (C 2 22 -C 2 21 )/C 22 , their corresponding values are 69.88 and 47.25 N/m, respectively. These values are much smaller than those of H-MoS 2 , H-MoN 2 H 2 , T d -MoN 2 H 2 , graphene and single layer h-BN

  meV.

	Configuration	NM	FM	AFM
	a( Å)	5.925	5.847	5.987
	b( Å)	5.950	5.909	5.867
	∆E(meV)	657.61	53.46	0.0
	m(µ B )	0.0	1.82	1.72
	Eg(eV)	0.36 0.54 a ,0.66 b 0.40
	d 1 ( Å)	2.09	2.15	2.13
	d 2 ( Å)	2.42	2.54	2.54
	α 1 (deg)	136.86 132.02 130.66
	α 2 (deg)	138.54 122.68 125.93
	(a) Spin-up, (b) Spin-down

TABLE I

 I

: Calculated in-plane lattice parameters (a and b ), relative energy ∆E, localized magnetic moments for each Co ion, energy gap Eg(eV), bond lengths d 1 =Co-S and d 2 = Co-Bi and bond angles α 1 :( Co -S -Co) and α 2 :( Co -S -Co) for spin and non-spin polarized configurations.


Nolan Mestres

Romain Vergne

Joëlle Thollot

Arthur Novat

Arthur Novat Controllable

Controllable Lighting Model for Designing Digital Panorama Maps in the Style of Novat

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Controllable

Multi-scale Adjusted Shadows

We compute a set of light directions to compute multi-scale shadows. We control their length to solve masking while aligning shadows with the relief.

To be continued...

Alpe

  Lighting Model for Designing Digital Panorama Maps in the Style of Novat Nolan Mestres 1 , Romain Vergne 1 , Joëlle Thollot 1 , Arthur Novat 2 | 1 Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP | 2 Atelier Novat Analytical Shading Multi-scale Local Light Alignment Whatever the shading model (diffuse, aspect based, global illumination,...) having a single light direction results in masking effects and contrast issues. For each detail scale, we move the light at each point of the surface to maximize contrast between both sides of features. References A Stylistic Study of the Winter Panorama Maps of Pierre Novat Nolan Mestres Cartographic Perspectives 2022, 10.14714/CP100.1753 Local Light Alignment for Multi-Scale Shape Depiction Nolan Mestres, Romain Vergne, Camille Noûs, Jöelle Thollot Eurographics 2021, 10.1111/cgf.142656