Keywords: adaptive filter, self-mixing, vibration measurement, embedded sensing, optical feedback, laser instrument

In this paper, an adaptive method of cancellation of parasitic vibrations is presented for a Self-Mixing (SM) interferometric laser vibration sensor which has been coupled with a solid state accelerometer (SSA). Previously, this was achieved by using a pre-calibration of phase and gain mismatches over the complete bandwidth of the instrument. Such a precalibration is not only tedious to execute but also hinders a mass production of the instrument as every SSA-SM sensor couple requires customized calibration. On the other hand, the proposed method does not require any pre-calibration as it uses an adaptive filter that self-tunes to match any unknown phase and gain differences between the SSA and the SM sensor. Two different adaptive algorithms, namely Recursive Least Squares (RLS) and Least Mean Squares (LMS) algorithms are tested and a comparison is established on the basis of parameter dependence, convergence time, computational cost and rms error. The proposed algorithms have provided improved results (mean error of 19.1nm and 20.2nm for LMS and RLS respectively) as compared to pre-calibration based results (mean error of 24.7 nm) for a laser wavelength of 785 nm. Simulated and experimental results thus demonstrate the utility of such an approach for embedded vibration sensing corrupted by extraneous parasitic motion.

I. INTRODUCTION

Self-mixing (SM) or optical feedback interferometry technique [START_REF] Donati | Developing self-mixing interferometry for instrumentation and measurements[END_REF][START_REF] Taimre | Laser feedback interferometry: a tutorial on the self-mixing effect for coherent sensing[END_REF][START_REF] Scalise | Selfmixing laser diode velocimetry: application to vibration and velocity measurement[END_REF] is being actively researched for the measurement of displacement [START_REF] Fan | Improving the measurement performance for a self-mixing interferometry-based displacement sensing system[END_REF][START_REF] Magnani | Real-Time Self-Mixing Interferometer for Long Distances[END_REF], distance [START_REF] Norgia | Absolute Distance Measurement with Improved Accuracy Using Laser Diode Self-Mixing Interferometry in a Closed Loop[END_REF], velocity [START_REF] Plantier | Real-time tracking of time-varying velocity using a self-mixing laser diode[END_REF], vibration [START_REF] Bernal | Study of Laser Feedback Phase under Self-Mixing leading to Improved Phase Unwrapping for Vibration Sensing[END_REF], flow [START_REF] Ozdemir | A Comparative Study for the Assessment on Blood Flow Measurement Using Self-Mixing Laser Speckle Interferometer[END_REF], profilometry [START_REF] Wang | Self-mixing microscopic interferometer for the measurement of microprofile[END_REF], range-finding [START_REF] Vogel | A low-cost medium-resolution rangefinder based on the self-mixing effect in a VCSEL[END_REF] and biomedical applications [START_REF] Donati | Self-Mixing Interferometry for Biomedical Signals Sensing[END_REF][START_REF] Milesi | Measurement of Local Chest Wall Displacement by a Custom Self-Mixing Laser Interferometer[END_REF]due to the compact, selfaligned, and low-cost nature of the SM instrument.

These metric instruments have performed well in laboratory conditions. However, more work needs to be done so that these instruments can also be used for industrial and/or embedded applications, where extraneous/parasitic mechanical vibrations can disturb the measurement. In this regard, a solution consists in coupling a Solid-State Accelerometer (SSA) with a SM vibration sensor [START_REF] Zabit | Design and Analysis of an Embedded Accelerometer Coupled Self-Mixing Laser Displacement Sensor[END_REF]. SSA then measures any extraneous/parasitic vibrations affecting the SM laser sensor. It then becomes possible to correct the corrupted measurement of SM sensor, even in real-time [START_REF] Zabit | Real-time accelerometer coupled Self-Mixing laser displacement sensor for embedded applications[END_REF]. Such a scheme thus potentially enables the use of SM metric sensors for embedded/industrial applications.

However, it needs to be highlighted that the SSA-SM sensor required a pre-calibration of its full bandwidth before actual use [START_REF] Zabit | Design and Analysis of an Embedded Accelerometer Coupled Self-Mixing Laser Displacement Sensor[END_REF]. Such a pre-calibration is mandatory in order to match the phase and gain differences between the SSA and the SM sensor. A phase and gain equalization filter was then designed in the light of the pre-calibration [START_REF] Zabit | Design and Analysis of an Embedded Accelerometer Coupled Self-Mixing Laser Displacement Sensor[END_REF][START_REF] Zabit | Real-time accelerometer coupled Self-Mixing laser displacement sensor for embedded applications[END_REF]. Otherwise, a correction of parasitic vibrations becomes ineffective. This pre-calibration, done by mechanically shaking the SSA-SM sensor head over a whole range of frequencies covering its operating bandwidth, can pose four problems. First, it becomes tedious as well as cumbersome as each and every SSA-SM couple must mounted and shaken for its precalibration. This can then hinder a mass production of such an instrument. Second, the mechanical resonances affecting the shaker used for the pre-calibration can falsify the extraction of phase and gain parameters. Third, the phase and gain equalization filter would need to be individually designed for every SSA-SM couple. Consequently, any imprecision either in the extraction of phase and gain parameters or in the design of equalization filter causes added residual error. Fourthly, any subsequent mismatches due to aging or component change would necessitate a re-calibration of the instrument. Therefore, in order to resolve all of these problems, in this paper, a method is proposed that enables cancellation of parasitic vibrations without needing any pre-calibration of the device. This is achieved by incorporating a self-tuning, adaptive filter. Two different adaptive algorithms have been tested, each providing different performance characteristics. The resulting sensor can thus adapt itself to the parasitic vibrations and provide corrected results after achieving convergence in an autonomous manner.

A schematic block diagram of the adaptive SSA-SM sensor is shown in Fig. 1. The SSA-SM sensor was mounted on a mechanical shaker (excited at different frequencies) in order to undergo parasitic vibrations, denoted as D s (t). A piezoelectric transducer (PZT) acted as a vibrating target as well as a reference sensor, denoted as D PZT (t). A photograph of the SSA-SM sensor can be seen in Fig. The paper is organized as follows. A brief introduction to SM interferometry is provided in section II. Then, the signal processing of the adaptive SSA-SM sensor is elaborated in section III. The simulated and experimental results are given in section IV, followed by the discussion and conclusion in Sections V and VI, respectively.

II. SELF-MIXING INTERFEROMETRY

The SM phenomenon happens in a laser when a part of the beam backscattered by a target is fed back into the active laser cavity, thereby causing interference with the emitted beam, which modifies the laser wavelength and output power. The variations in the optical output power P(t) caused by this optical feedback can be written as [START_REF] Donati | Developing self-mixing interferometry for instrumentation and measurements[END_REF]:

 ( ) = [1 + .

cos ( ) ] 

Where P 0 is the emitted optical power under free-running conditions, m is the modulation index and x F (t) is the laser output phase in the presence of feedback, given by:

 ( ) = 2 ( ) ( )/

 

Where D(t) is the target displacement. Under optical feedback, x F (t) is determined by the wellestablished Lang-Kobayashi model [START_REF] Lang | External optical feedback effects on semiconductor injection laser properties[END_REF], given as

 ( ) -( ) -sin[ ( ) + arctan( )] = 0 
where α is the so called Henry's factor [START_REF] Henry | Theory of the linewidth of semiconductor lasers[END_REF], also known as linewidth enhancement factor [START_REF] Yu | Measurement of the linewidth enhancement factor of semiconductor lasers based on the optical feedback self-mixing effect[END_REF], C is the feedback coupling factor [START_REF] Acket | The influence of feedback intensity on longitudinal mode properties and optical noise in index-guided semiconductor lasers[END_REF], also known as Acket's parameter [START_REF] Taimre | On the nature of Acket's characteristic parameter C in semiconductor lasers[END_REF][START_REF] Yu | Optical Feedback Self-Mixing Interferometry With a Large Feedback Factor C: Behavior Studies[END_REF] that determines the SM operating regime [START_REF] Bernal | Classification of laser SMI signal under moderate feedback[END_REF] and x 0 (t) is the laser output phase in the absence of feedback, found by replacing λ F (t) with λ in (2), where λ is the laser diode emission wavelength under free-running conditions.

III. SIGNAL PROCESSING

The signal processing of adaptive SSA-SM sensor (see Fig. 3) can be grouped in three major parts, as explained in the following.

A. Self-mixing Interferometric Signal:

The SM interferometric signal corresponds to the signal P(t). Traditionally, P(t) is easily acquired by using the built-in photodiode (PD) found inside the laser diode (LD) package (see Fig. 1). This SM interferometric signal can then be processed to recover the corresponding displacement or vibration signal. In this paper, we have used the Phase Unwrapping Method (PUM) offering a precision of λ/16 for displacement retrieval [START_REF] Bes | Displacement measurements using a self-mixing laser diode under moderate feedback[END_REF]. Note that for experimental SM signal acquisitions, the PUM automatically takes care of any variations in P(t) caused by possible changes in system parameters, such as injection current, feedback strength, and target external cavity length, through the use of an automatic gain control stage [START_REF] Bes | Displacement measurements using a self-mixing laser diode under moderate feedback[END_REF].

Note that as shown in Fig. 3, the SM signal based measurement is denoted as D ∑ (t) because it represents the sum of true target motion and parasitic motion i.e. D ∑ (t)= D PZT (t)+ D s (t). It may be noted that it is assumed for the rest of this paper that the parasitic vibrations disturbing the SSA-SM sensor are of such an amplitude that the total motion does not influence the laser source and feedback based system parameters (e.g., C and α). It will be seen later in this paper that this assumption holds true for experimental SM signals acquired in the presence of parasitic vibrations.

B. Acceleration Signal:

As already mentioned, the extraneous vibration affecting the sensor (noted as D s (t)) is measured by using a SSA. For this purpose, the acceleration signal a cc (t) is band-pass filtered and double-integrated (as shown in Fig. 3). The band-pass filtering is done so that low-frequency drifts that can cause offsets in the subsequent integration steps as well as high frequency acceleration signal saturation can be dealt with effectively.

C. Adaptive Filter:

An adaptive filter is a self-designing and time-varying system that continuously adapts/tunes its filter coefficients by following an adaptive algorithm with an aim of minimizing a cost/error function. Among other applications, such filters have been used for noise and echo cancellation [START_REF] Paulo | Adaptive Filtering Algorithms and Practical Implementation Third Edition Springer[END_REF].

In our case, the adaptive filter iteratively tunes itself by minimizing the mean square error (MSE) based on e(t) = D ∑ (t) -D s (t) (see Fig. 3). For statistically stationary inputs, the filter is said to have converged (i.e. tuned or designed) when it achieves the minimum MSE. Afterwards, it continues to perform its task (in our case, parasitic vibration cancellation) as long as the inputs retain their statistical nature. However, in case of change in the nature of inputs, the filter starts adapting itself to the new situation, and again achieves convergence to perform its task. For the present proof of concept, a transversal finite impulse response (FIR) adaptive filter structure has been evaluated using two different adaptive algorithms [least mean squares (LMS) and recursive least squares (RLS)] as detailed in the following.

1) Least Mean Square Algorithm (LMS):

The LMS algorithm, known for its simplicity, stability, and robustness [START_REF] Paulo | Adaptive Filtering Algorithms and Practical Implementation Third Edition Springer[END_REF], belongs to the class of adaptive filters that adapts the filter coefficients by minimizing the MSE signal, where the error signal is the difference between the desired signal and the input defined In terms of adaptive filter terminology, the error e(n) between the desired and weighted input signals is given by [START_REF] Paulo | Adaptive Filtering Algorithms and Practical Implementation Third Edition Springer[END_REF] ( ) = ( ) -( )

In our case, the corrupted signal D ∑ (t) becomes desired signal d(n)and the parasitic signal D s (t) becomes input signal X(n).

The filter coefficients in LMS algorithm are updated by:

( + 1) = ( ) + 2 ( ) ( ) (5)
For a given order of the filter N, the rate of convergence and filter's stability are determined by the step size (or convergence factor) denoted as µ, determined by

0 < < ( ) (6)
Where P av is the average power of input signal D s (t). Thus, the LMS filter's stability and convergence rate become a function of average power of signal.

2) Recursive least square algorithm:

The Recursive Least Squares (RLS) adaptive algorithm recursively tunes the filter coefficients in order to minimize a weighted linear least squares cost function with respect to the input signal. The RLS algorithm is known for its excellent performance while working in time varying (non-stationary signals) environment but at the cost of an increased computational complexity and some stability problems [START_REF] Saeed | Advanced Digital Signal Processing and Noise Reduction Fourth Edition[END_REF]. RLS filter coefficients are updated by

( ) = ( -1) + ( ) ( ) (7)
where k(n) is the filter gain given by

( ) = ∅ ( ) ( ) ( )∅ ( ) ( ) ( 8
)
Where Λ is the adaptation factor, and ∅ is the correlation matrix, initially set to to identity matrix δI, and is recursively updated by.

∅ (n) = Λ ∅ ( -1) -Λ ( ) ( )∅ ( -1)(9)
Thus, the major difference between the LMS and RLS algorithms is the use of ∅ in the correction term of the RLS algorithm which helps in achieving better convergence rate by decorrelating the successive inputs [START_REF] Saeed | Advanced Digital Signal Processing and Noise Reduction Fourth Edition[END_REF].

IV. RESULTS

A. Experimental Setup

The experimental set-up deployed for the validation of adaptive SSA-SM sensor has been schematized in Fig. 1. The employed SSA is a SF1500 accelerometer from Colibrys ® (typical noise resolution of 0.3µg/Hz 1/2 and a full-scale range of ±3g). The SM sensor uses a Sanyo ® DL7140 laser diode (λ=785 nm) with an output power of 50 mW. The SSA was glued on the SM sensor laser head so that it could measure as correctly as possible the parasitic movement D s (t) undergone by the SM sensor head.

The SSA-SM sensor head (an approximate size of 3 cm ×6 cm) was mounted on a mechanical shaker (Fig. 2) that was used to generate vibrations disturbing the sensor. A commercial PZT actuator from PhysikInstrumente (P753.2CD) served as target. This device also has a built-in capacitive feedback sensor with 2nm resolution that served as a reference sensor for the PZT target movement D PZT (t). It thus allowed the calculation of error between the corrected signal D c (t) and the reference motion D PZT (t). Simulated and experimental tests, conducted using RLS and LMS algorithms to validate the utility of adaptive filter for the cancellation of parasitic vibrations corrupting the SM sensor output, are detailed in the following.

B. Results of Simulated Signals

In order to test if an adaptive filter can allow parasitic vibration cancellation in a SM vibration sensor, various simulation based cases were evaluated by using LMS and RLS algorithms, as detailed below.

First, it was simulated to test if the adaptive filter can tune itself to provide correction if the target as well as the sensor are both vibrating at two different frequencies. Second, it was simulated to test if the adaptive filter can provide acceptable correction if multiple parasitic harmonic vibrations disturb the SM vibration sensor. Fig. 5 represents case where the parasitic vibrations disturbing the SM sensor are at f s = 61Hz-122Hz-305Hz-427Hz while the target vibration is at f PZT = 81 Hz. As a consequence, the vibration retrieved through the SM signal is heavily corrupted, as seen in Fig. 5 (b). However, the corrected vibration D c (t)again matches very well with D PZT (t)shown in Fig. 5 (d) with an rms error of 8.82nm (LMS) and 12.2nm (RLS). 61Hz-122Hz-305Hz-427Hz d. 41Hz-82Hz-205Hz-287Hz

Various other simulations were also performed (see Table I). All of these gave satisfactory results thereby validating the use of adaptive filter in the absence of pre-calibration. It may be noted that all results based on RLS algorithm, simulated as well as experimental used constant filter order N RLS of 35 and adaptation factor Λ of 1. On the other hand filter order N LMS and convergence factor of LMS algorithm were varied for each case of Table I due to reasons detailed later on in Section V. Using the experimental set-up already detailed, many experimental signals were acquired to validate the principle. The experimental acquisitions are at the same frequencies as already used in the simulated signals. Thus, it would be possible to make a meaningful comparison between the simulated and experimental cases.

Therefore, first, the PZT vibration was set at f PZT = 85 Hz with A p-p = 5 µm while parasitic vibration disturbing the sensor is at f s = 52 Hz with A p-p = 7 µm (measured by SF1500 SSA). The corresponding experimental SM signal is shown in Fig. 6 (a) which has been unwrapped by using PUM to provide D ∑ (t) shown in Fig. 6(b). In spite of the higher vibration amplitudes, the use of adaptive filter has corrected it to provide D c (t)which matches very well with the reference target capacitive feedback sensor vibration D PZT (t) shown in Fig. 6 (d). The error (t) = D PZT (t) -D c (t) is shown in Fig. 6 (e) with an rms value of 13.6 nm(LMS) and 16.9nm (RLS).

Second, an experiment was conducted for the case where the shaker is excited by 229Hz while PZT target was excited by a signal composed of 59Hz-118Hz-295Hz-413Hz (see Fig. 7).In spite of the fact that this case is more complicated than the previous case but still a good correction has been achieved with an rms error of 18.2nm (LMS) and 20.5nm (RLS).

Third, an experiment was conducted to see if multiple parasitic harmonic vibrations disturbing the SM vibration sensor can also be cancelled by using the proposed methods. Fig. 8 thus shows the experimentally acquired signals for f s = 51Hz-102Hz-255Hz-357Hz and f PZT = 41Hz-82Hz-205Hz-287Hz. In spite of multiple mechanical vibrations disturbing the SM sensor, a good recovery of true target vibration has been achieved, as seen in Fig. 8 (e) with an rms error of 28.5 nm(LMS) and 24.2nm (RLS).The frequency spectra of Fig. 9 show corrections of 20.3dB (28.4dB), 24.2dB (21.3dB), 18.3dB (19.1dB), and 20.1dB (18.7dB)at51Hz, 102Hz, 255Hz, and 357Hz respectively using LMS algorithm (RLS algorithm).

Finally, the results of these and all other experimental cases are presented in Table II.

V. DISCUSSION

A. Performance of LMS based Adaptive Sensor

A comparison of Table I and Table II indicates that the error values are expectedly higher for experimental signals as compared to simulated signals due to the addition of different error sources, such as the imprecision of SM displacement retrieval method, the imprecision of the SSA measuring the parasitic movement, the noise of electronic circuits and data acquisition path etc. Similarly, the filter orders of LMS algorithm are generally higher for the experimental signal It can be seen from Table I and Table II that LMS based optimal correction is only achieved when its parameters such as N and are adjusted as function of input signal's power P av as per [START_REF] Norgia | Absolute Distance Measurement with Improved Accuracy Using Laser Diode Self-Mixing Interferometry in a Closed Loop[END_REF].

In order to achieve optimal correction results, note that the filter order for any given signal was kept on increasing so long as final error kept on reducing while maintaining a stable convergence of the algorithms. Thus, the highest filter order leading to the lowest errors are reported. For such a filter order, the final error of the adaptive system reaches the noise floor. Subsequently, any further increase in filter order would not lead to improvement in the final error results. An analysis of tabulated results also brings forth the relationship between filter order, convergence factor, and input signal power for LMS based adaptive FIR filter. It can be seen in Table II that higher filter orders result in smaller convergence factors which cause slower convergence. On the other hand, average power of signal also affects the convergence factor. For example, 143Hz-46Hz case has larger convergence factor as compared to 123Hz-73Hz case even though the filter order is the same. This is due to the difference in their respective average powers. Therefore, it is both filter order and average signal power that determine stable convergence of LMS adaptive filter.

B. Performance of RLS based Adaptive Sensor

In the case of RLS algorithm based SSA-SM sensor, parasitic vibration correction is independent of parameters adjustment. Thus, its filter order N RLS and adaptation factor Λ can be kept constant for all the cases. This makes it suited for

C. Comparison of LMS and RLS based Systems / Sensors

A comparison of rate of convergence of LMS and RLS algorithms is shown in Fig. 10 which uses the experimental data of f s =92Hz and f PZT =65Hz case. When the performance of RLS is compared with that of LMS for the same filter of 80 then it is seen that RLS achieves convergence at a better rate (e.g. 173 iteration for RLS and 1100 iterations for LMS for filter order of 80, leading to an rms error of 14.7nm and 28.9nm respectively) Focusing only on LMS algorithm, the impact of choosing higher filter order (leading to smaller error) on the convergence time (causing slower convergence) of LMS algorithm based adaptive filter is also depicted in Fig. 10. Here, the filter roughly took 3700 iterations to reach convergence resulting in final rms error of 12.5 nm for N = 400.

Thus, the use of RLS algorithm has resulted in achieving convergence after 17.3 ms. For an input sampling rate of 10 000 samples per second, equating to a total computational time per iteration of 0.1 ms in the context of a real-time system, the adaptive SSA-SM system achieves convergence in 17.3 ms for the RLS algorithm for the above-mentioned case. Even in the worst case (i.e., the use of LMS algorithm with NLMS = 400), convergence is achieved after 0.37 s, thereby underlining the performance of the adaptive filter scheme.

Regarding measurement precision, the results of the RLS and LMS algorithms are almost comparable to each other with minor differences (see Tables I andII). Yet, we can claim that the performance of the RLS algorithm is better than that of the LMS algorithm because the filter order of the RLS algorithm is very small compared with that of the LMS algorithm, and it achieves convergence at a faster rate without any external intervention. However, this superior performance of RLS is attained at the expense of a large increase in computational complexity, as detailed below.

The complexity level of the RLS algorithm requires a total of 4N 2 +4N+2 multiplications (division counted as multiplication) while LMS algorithm requires 2N+1 multiplications only, where N is filter's order [START_REF] Ben Allen | Adaptive Array Systems fundamental and its application[END_REF].

The advantage of LMS algorithm is its simple structure requiring only vector operations as seen in [START_REF] Fan | Improving the measurement performance for a self-mixing interferometry-based displacement sensing system[END_REF][START_REF] Magnani | Real-Time Self-Mixing Interferometer for Long Distances[END_REF]. The RLS algorithm structure, on other hand, is much more complex, requiring the calculation and updating of input autocorrelation matrix as seen in [START_REF] Plantier | Real-time tracking of time-varying velocity using a self-mixing laser diode[END_REF][START_REF] Bernal | Study of Laser Feedback Phase under Self-Mixing leading to Improved Phase Unwrapping for Vibration Sensing[END_REF][START_REF] Ozdemir | A Comparative Study for the Assessment on Blood Flow Measurement Using Self-Mixing Laser Speckle Interferometer[END_REF].

D. Comparison with pre-calibration based sensor / system.

Finally, a comparison can be made between the performance of the adaptive filter based instrument and that of the previously proposed pre-calibration based instrument [START_REF] Donati | Self-Mixing Interferometry for Biomedical Signals Sensing[END_REF]. Note that the last column of Table II cites the published results of Zabit et al. [START_REF] Zabit | Design and Analysis of an Embedded Accelerometer Coupled Self-Mixing Laser Displacement Sensor[END_REF] and that the authors have used the same experimental signal acquisitions for the present study. It can be seen that the proposed method enables improved mean rms error of 19.1nm and 20.2nm for LMS and RLS respectively as compared to pre-calibration based results having mean rms error of 24.7nm.

It thus highlights the performance of adaptive filter based instrument which not only removes the problems involved in the design of pre-calibration based instrument but also leads to better correction results. This improved performance of the adaptive filter based instrument can be explained by the absence of errors caused by imprecision either in the extraction of phase and gain parameters or in the design of equalization filter. Furthermore, by its very nature, at any given time, the adaptive filter is able to concentrate its attenuation over a very small set of frequencies whereas the equalization filter provides equally weighted attenuation over the whole operating bandwidth of the instrument. Thus, except for the cases where f s =f PZT which cannot be solved by an adaptive filter, the proposed method not only helps to avoid the difficulties associated with pre-calibration but also provides on average better measurement precision.

Finally, the simple and computationally light LMS based adaptive filter is only suitable for a real-time system implementation without external intervention if additional processing is incorporated in it to adjust filter order and convergence factor as a function of input signal's power. Such an obstacle is completely removed if the RLS algorithm is deployed as it requires no external intervention and works optimally in an autonomous manner.

VI. CONCLUSION

In this paper, it has been demonstrated that the use of adaptive filters can allow us to design a laser vibration sensor that can provide correct measurements even when the sensor itself is disturbed by extraneous parasitic motion. The use of adaptive filter with an accelerometer coupled Self-Mixing laser instrument thus potentially allows sensing (without needing any pre-calibration) in hostile/embedded environment where such extraneous movements which mechanically disturb the sensor cannot be avoided.

For this purpose, LMS algorithm was used with a transversal FIR adaptive filter structure. It was chosen over other advanced algorithms as it requires less computational © IEEE 2017. . Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works DOI: 10.1109/TIM.2016.2626018 8 power and memory while achieving the task of parasitic vibration cancellation in a stable manner. However, the main drawback of LMS algorithm is the difficulty of predefining the value of learning rate or convergence factor µ that ensures the adaptive filter's stability in case of variations in input signal power, as seen in [START_REF] Norgia | Absolute Distance Measurement with Improved Accuracy Using Laser Diode Self-Mixing Interferometry in a Closed Loop[END_REF].

For these reasons, RLS algorithm would be then be preferable choice for a real-time adaptive SSA-SM sensor as its filter parameters are independent of input signal's characteristics. Furthermore it has superior performance as compared to LMS algorithm in terms of faster convergence, the rms error of both algorithms are also comparable, so one can safely say RLS is the better choice but at the cost of additional computational complexity.

The main advantage of the LMS filter is the comparative simplicity of the algorithm. However, for signals with a large eigenvalue spread, the LMS has an unstable and slow convergence rate [START_REF] Saeed | Advanced Digital Signal Processing and Noise Reduction Fourth Edition[END_REF]. In addition to this non-stationary signals with high rate of change with time, the LMS can be an unsuitable adaptation solution, on the other hand, RLS method, with its better convergence rate and less sensitivity to the eigenvalue spread, becomes a more attractive alternative [START_REF] Saeed | Advanced Digital Signal Processing and Noise Reduction Fourth Edition[END_REF].

For the different monotone experimental parasitic vibration cases (tabulated in Table II), the proposed methods has resulted in mean rms error in final displacement correction of better than 15 nm. It is a noteworthy result as it is a precision value similar to what can be obtained with a SM sensor by using PUM in the absence of parasitic vibrations. That is, for monotone parasitic vibrations, the adaptive SSA-SM system has provided almost ideal correction.

Fig. 1 .

 1 Fig. 1. Block diagram of the adaptive Solid-State Accelerometer coupled Self-Mixing (SSA-SM) sensor: photodiode (PD), laser diode (LD), focusing lens (FL), and piezoelectric transducer (PZT).

Fig. 2 .

 2 Fig. 2.Photograph of SSA-SM sensor using the SF1500 accelerometer and a SM displacement sensor based on DL7140 laser diode. The sensor head has been mounted on a mechanical shaker to undergo parasitic vibrations.

Fig. 3 .

 3 Fig. 3.Signal processing (shown in dotted red block) required for the adaptive SSA-SM sensor: PUM (Phase Unwrapping Method [23]).

Fig. 4

 4 Fig. 4 represents such a case where the target vibration is at f PZT = 85 Hz with a peak-to-peak amplitude A p-p = 5 µm while parasitic vibration disturbing the sensor is at f s = 52 Hz again with A p-p = 5 µm. The corresponding SM signal, simulated by

  -295Hz-413Hz c 51Hz-102Hz-255Hz-357Hz b.

Fig. 4 .

 4 Fig. 4. Simulation for fs = 52Hz and fPZT = 85 Hz: (a) simulated SM interferometric signal, (b) corrupted vibration D∑(t)retrieved by PUM, (c) parasitic vibration Ds(t), (d) corrected vibration D c (t)(blue (LMS) , green (RLS)) and reference vibration DPZT(t) (dotted red), and (e) error (t)= DPZT(t)-D c (t) (blue (LMS), dotted green (RLS)) over 20Hz-500Hz.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Simulation for fs = 61Hz-122Hz-305Hz-427Hz and fPZT = 81 Hz: (a) simulated SM interferometric signal, (b) corrupted vibration D∑(t) retrieved by PUM, (c) parasitic vibration Ds(t), (d) corrected vibration D c (t) (blue (LMS) , green (RLS)) and reference vibration DPZT(t) (dotted red), and (e) error (t) = DPZT(t) -D c (t) (blue (LMS), dotted green (RLS)) over 20Hz-500Hz.

Fig. 7 .

 7 Fig.7. Experimental signals for fs = 229Hz and fPZT = 59Hz-118Hz-295Hz-413Hz : (a) experimental SM interferometric signal picked up by DL7140 (b) corrupted vibration D∑(t)retrieved by PUM, (c) parasitic vibration Ds(t) measured by SF1500 SSA, (d) corrected vibration D c (t) (blue (LMS) , green (RLS)) and reference capacitive feedback sensor vibration DPZT(t)(dotted red), and (e) error (t) = DPZT(t) -D c (t) (blue (LMS), dotted green (RLS)) over 20Hz-500Hz.

Fig. 9 .

 9 Fig.9.Experimental signals for fs = 51Hz-102Hz-255Hz-357Hz and fPZT = 41Hz-82Hz-205Hz-287Hz: (a) frequency spectrum of corrupted vibration D∑(t), (b) LMS and (c) RLS algorithm based frequency spectrum of corrected vibration D c (t)indicating an improvement corrections of 20.3dB (28.4dB), 24.2dB (21.3dB), 18.3dB (19.1dB), and 20.1dB (18.7dB) at 51Hz, 102Hz, 255Hz, and 357Hz.

Fig. 8 .

 8 Fig. 8.Experimental signals for fs = 51Hz-102Hz-255Hz-357Hz and fPZT = 41Hz-82Hz-205Hz-287Hz: (a) experimental SM interferometric signal picked up by DL7140, (b) corrupted vibration D∑(t)retrieved by PUM, (c) parasitic vibration Ds(t)measured by SF1500 SSA, (d) corrected vibration D c (t) (blue (LMS) , green (RLS)) and reference capacitive feedback sensor vibration DPZT(dotted red), and (e) error  (t)= DPZT (t)-D c (t) (blue (LMS), dotted green (RLS)) over 20Hz-500Hz.

Fig. 10 .

 10 Fig.10. Evolution of mean square error (MSE) as function of number of iterations needed by the adaptive LMS and RLS algorithms for the experimental case of fs =92Hz and fPZT = 65Hz.

  2.

	Adaptive Cancellation of Parasitic Vibrations
	Affecting a Self-Mixing Interferometric Laser
	Sensor
	ZohaibA. Khan, UsmanZabit, Member, IEEE, Olivier D. Bernal, Member, IEEE, Muhammad O.
	Ullah, Member, IEEE,andThierry Bosch,Senior Member, IEEE

TABLE I .

 I ADAPTIVE FILTER SPECIFICATION AND RMS ERROR RESULTS OF SIMULATED SIGNALS FOR DIFFERENT VIBRATIONS OF TARGET & SHAKER

	Shaker (Hz)	PZT (Hz)	Signal Avg Power Pav	Filter Order N	Conve-gence factor µ	LMS RMS Error (nm)	RLS RMS Error (nm)
	52	85	9.25	300	7e-5	4.41	13.2
	92	65	5.12	350	7e-5	9.74	7.44
	143	46	5.11	170	8e-4	7.97	1.89
	223	131	6.25	200	9e-4	5.61	3.80
	43	105	9.25	280	6e-5	5.88	7.26
	123	73	7.62	360	9e-5	4.51	5.10
	229	59 arb a	7.20	300	7e-5	8.42	3.89
	61 arb b	81	4.12	170	4e-4	8.82	12.2
	51 arb c	41arb d	4.16	700	9e-5	14.1	9.49
					mean	7.71	7.14

OVER (20HZ-500HZ) BANDWIDTH

TABLE II .

 II ADAPTIVE FILTER SPECIFICATION AND RMS ERROR RESULTS FOR EXPERIMENTAL SIGNALS OF DIFFERENT VIBRATIONS OF TARGET & SHAKER OVER (20HZ-500HZ) BANDWIDTH IEEE 2017. . Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works

									7
	practical, real-time experimental measurements scenarios as it
	requires no intervention in terms of its parameters. As a result,
	the SSA-SM sensor using this adaptive filter is then able to
	work	autonomously	while	delivering	comparable
	measurements precision at even better convergence rates.
					Sha-ker (Hz)	PZT (Hz)	Signal Avg Power Pav	Filter order NLMS	Conve-gence factor µ	LMS RMS Error (nm)	RLS RMS Error (nm)	Exemp -lar [14] Error
									(nm)
					52	85	10.37	390	6e-5	13.6	16.9	31.9
					92	65	5.785	400	4e-5	12.5	15.8	24.9
					143	46	5.993	120	3e-4	17.1	17.9	17.0
					223	131	6.384	260	9e-4	14.7	16.7	17.4
					43	105	6.209	380	4e-5	13.7	14.8	22.0
					123	73	7.936	120	2e-4	11.5	16.5	28.1
					229	59 a	3.411	530	7e-5	18.2	20.5	16.9
					61 b	81	7.254	450	7e-5	41.7	38.7	31.8
					41 c	51 d	3.092	500	4e-5	28.5	24.2	32.1
					--	--	--	--	mean	19.1	20.22	24.7
					a .59Hz-118Hz-295Hz-413Hz	c 51Hz-102Hz-255Hz-357Hz
					b. 61Hz-122Hz-305Hz-427Hz	d. 41Hz-82Hz-205Hz-287Hz
	DOI: 10.1109/TIM.2016.2626018

©


Meng Liao

Quy-Dong To
email: quy-dong.to@u-pem.fr

Céline Léonard

Wenlu Yang

Prediction of thermal conductance and friction coefficients at solid-gas interface from statistical learning of collisions

Keywords:

In this paper, we present the construction of statistical models of gas-wall collision based on data issued from Molecular Dynamics (MD) simulations and use them to predict the velocity slip and temperature jump coefficients at the gas-solid interface. The Gaussian Mixture (GM) model, an unsupervised learning technique, is chosen for this purpose. The model shares some similarities with the well-known Cercignani-Lampis model in kinetic theory but it is more robust due to the unlimited number of Gaussian functions used and the ability to deal with correlated data of high dimensions. Applications to real gas-wall systems (Argon-Gold and Helium-Gold) confirm the good performance of the model. The trained GM model predicts physical and statistical properties including accommodation, friction and thermal conductance coefficients in excellent agreement with the MD model.

I. INTRODUCTION

Velocity slip and temperature jump are important interface phenomena which are frequently encountered in micro-nanofluidic systems [START_REF] Karniadakis | Microflows and nanoflows: fundamentals and simulation[END_REF]. The common way to model the imperfect boundary conditions is to use the following phenomenological relations

β 1 v k = -σ kz , β 2 (T -T w ) = -q z , k = x, y, (1)
Here the wall is assumed to be stationary at temperature T w and normal to direction z.

The slip velocity v k is proportional to the (minus) shear stress σ kz and the temperature jump T -T w is proportional to the heat flux q z . The constants, β 1 and β 2 , are respectively the Navier isotropic friction and Kapitza thermal conductance coefficients. For gas flow problems, it is well-known that the interface phenomena become significant when the mean free path λ is comparable to the channel height H, starting from λ/H ≥ 0.01. The origin of the Knudsen layer effect comes from two sources: the finite actual velocity/temperature jump at the boundary and the deviation from the bulk solution within the distance λ from the wall. In this paper, we do not consider the variation of the velocity/temperature in the Knudsen layer and assume that the linear relations Eq. ( 1) are valid at the solid boundary [START_REF] To | [END_REF].

Since both the bulk and interfacial transport properties of gases are governed by collisions at the atomic scale [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction, and diffusion in gases[END_REF], they are expected to be proportional to the collision rate. It is thus natural to non-dimensionalize the boundary coefficients β 1 and β 2 as follows

β 1 = β 1 n π 2k B T m , β 2 = β 2 n mπ 8k 3 B T , k = x, y, (2)
In Eq. [START_REF] To | [END_REF], n is the gas density at the wall, m the atomic mass of the monatomic gas and k B the Boltzmann constant. The other reason to non-dimensionalize the interfacial coefficients is that β 1 and β 2 are connected to tangential momentum accommodation coefficient α 1 and energy accommodation coefficients α 2 [START_REF] Karniadakis | Microflows and nanoflows: fundamentals and simulation[END_REF]

β 1 = α 1 2 -α 1 , β 2 = α 2 2 -α 2 (3)
The origin of the accommodation coefficients comes from the collision models of gas atoms with the solid boundary [START_REF] Maxwell | [END_REF][5][6]. These constants appear in the linear relation between the pre-collision and post-collision momenta and energy of gas atoms with respect to the wall.

However, MD simulations showed that the linearity postulate is not always valid. Errors from the linear regression procedure to obtain α 1 and especially α 2 are observed [7,8]. It is suggested that the more effective way to model slip and jump phenomena is to directly use coefficients β 1 and β 2 to quantify them from realistic collision data whenever possible.

On the other hand, the formulation based on coefficients β1 , β2 can avoid the errors due to the use of (3) and the intermediate accommodation coefficients. It is noted that while Eqs.

(3) are widely used, they are derived under restrictive conditions, e.g the Chapman-Enskog distribution [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction, and diffusion in gases[END_REF], an approximate solution of the Boltzmann equation, and the Maxwell or

Cercignani-Lampis wall model [9,[START_REF] Grad | Thermodynamik der Gase/Thermodynamics of Gases[END_REF]. As a result, if we want to extend or investigate the validity (3) in general situations, more accurate methods (see e.g Refs. [START_REF] Sharipov | [END_REF]12]) should be adopted.

In the framework of probabilistic modeling, there are numerous ways of reconstructing the (collision) behavior of the system based on the available (collision) data, which are realizations via Molecular Dynamics simulations. One of the method, proposed in our previous work is to use the nonparametric (NP) techniques to estimate the probability density functions (PDF) and generate new collision data [8]. To avoid the sparsity of the data in high dimension, the reflection process was assumed to be independent for each velocity components (c x , c x ), (c y , c y ) and (c z , c z ) and the model can be constructed with a good accuracy [START_REF] Hollander | Nonparametric statistical methods[END_REF]. Although such modeling is relatively general, the uncorrelation hypothesis can be strong and oversimplify the true behavior. This motivates us to investigate the use of unsupervised learning techniques and their performance with high dimension data.

Specifically, the present work considers the Gaussian Mixture (GM) model, which is a popular technique and has proved to be successful in many data science applications. The GM are probabilistic models using linear combination of multiple multidimensional gaussians to estimate the probability density of data. Like all collision models (or scattering kernels) in kinetic theory, the construction of GM are based on parameters which can be determined from collision data for each gas-solid couple under consideration. While the existing collision models are limited by a finite number of parameters (e.g accommodation coefficients), the GM models are not. On the other hand, the GM models require a special algorithm, i.e

Expectation Maximization (EM) algorithm to tune the parameters.

The paper is organized as follows. In Sect. II, we present a literature review of statistical models of collisions, the general principles to construct them from MD data and their properties. GM model is then introduced with the EM algorithm to tune the model parameters. Criteria based on both statistical and physical parameters including accommodation coefficients, interfacial coefficients, reflection angles, and Kullback-Leibler divergence will be used to evaluate the performance of the model. Applications to real material systems (Ar and He gases vs Au surface) are presented in Section III, confirming the superiority of the GM model when compared with nonparametric and parametric models of literature.

Finally, some remarks and conclusions are given.

II. DESCRIPTION OF COMPUTATION METHODS

A. Overview of statistical models of collision and their properties

The collision model, also called scattering kernel in kinetic theory, serves as boundary conditions for Boltzmann equations and aims at describing how gas atoms are reflected after collision with the wall at temperature T w . It is defined by the conditional probability distribution, P (output = c|input = c , temperature = T w ) or denoted shortly P (c|c ) of velocity c for a given incoming velocity c . In addition to the usual PDF properties like nonnegativeness and normalization, P (c|c ) is also expected to satisfy the physical reciprocity conditions [START_REF] Cercignani | Theory and application of the Boltzmann equation[END_REF][START_REF] Brull | [END_REF][16]:

P eq (c )P (c|c ) = P eq (-c)P (-c | -c), c ∈ R 2 × R + , c ∈ R 2 × R - (4)
where, P eq (c ) is the distribution associated to the equilibrium gas at the same temperature as the wall T = T w and motionless with respect to the wall (zero mean velocity)

P eq (c ) = 2π θ |c z |M B (c ), M B (c ) = e -c 2 /2θ (2πθ) 3/2 , θ = k B T /m (5)
This condition implies that the incidence-reflection is microscopically reversible when the wall is in equilibrium state.

It can be shown that all popular collision models in literature can be classified as parametric models, in most cases, taking from one to three accommodation coefficients as parameters.

For example, in the Maxwell model [START_REF] Maxwell | [END_REF] where the gas atom can be reflected either specularly or diffusively with percentage 1 -α and α, only one accommodation coefficient is used.

In the anisotropic Cercignani-Lampis (ACL) model [5,6], three different accommodation coefficients associated to tangential momentums along x, y directions and normal kinetic energy along z direction are used. While the Maxwell model exhibits discontinuities in the probability space due to the portion of the mirror reflections (Dirac distribution), the smooth PDF from ACL model agrees better with MD simulations. However, like all parametric models, they are subject to common limitations: the requirement of the existence of a limited number of constant accommodation coefficients and the lack of flexibility in the case of complex collision data. Since the flow solution depends both on the gas-gas and gas-wall collision, those limitations may have consequences on the determination of boundary conditions' coefficients β 1 and β 2 in Navier-Stokes-Fourier equations as mentioned earlier and can cause errors in other simulation methods, e.g direct simulation Monte Carlo (DSMC) [17], MD [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF], moment equations [START_REF] Grad | Thermodynamik der Gase/Thermodynamics of Gases[END_REF][START_REF] Struchtrup | [END_REF], etc. based on the collision model. It suggests that to achieve better accuracy, we need to rely on the whole data range and to not be limited by the number of model parameters. The procedure of reconstructing the scattering kernel from MD simulations of collisions is shown in Fig. 1. The gas atoms are inserted with velocity c at a plane beyond the cut-off distance and interact with the well thermostated solid atoms before escaping by crossing the same plane with velocity c. Each couple of input and output velocities (c, c ) recorded is thus a realization of the joint PDF P (c, c ). In principle, if we can estimate accurately P (c, c ), we can derive accurately the scattering kernel via the equation

P (c|c ) = P (c, c ) P (c ) , (6)
where P (c ) is the distribution of input variable c . While P (c|c ) is a unique function of c and c , the joint PDF P (c, c ) depends on the input distribution P (c ) via (6). In the special case where the input distribution is the equilibrium distribution, P (c ) = P eq (c ), the joint PDF becomes the equilibrium joint PDF P (c, c ) = P eq (c, c ). The reciprocity condition ( 4)

is equivalent to the symmetry condition of P eq (c, c ) in hyperspace

P eq (c, c ) = P eq (-c , -c) (7)
As a result, to model the scattering kernel via P eq (c, c ) and the associated collision data, it is necessary that the above constraints are satisfied.

We remark that in many parametric scattering kernels [5,6,21], the reflection process occurs independently in x, y, and z directions. Under this assumption, both P (c|c ) and P eq (c, c ) can be decomposed as:

P (c|c ) = 3 k=1 P (c k |c k ), P eq (c, c ) = 3 k=1 P eq (c k , c k ), k = x, y, z (8)
A typical example is the well-known Cercignani-Lampis kernel and its associated joint PDF at equilibrium P CL eq (c, c )

P CL (c|c ) = e - [cx-(1-α t )c x ] 2 2α t (2-α t )θw 2πθ w α t (2 -α t ) e - [cy -(1-α t )c y ] 2 2α t (2-α t )θw 2πθ w α t (2 -α t ) c z e -c 2 z +(1-αe)c 2 z 2αeθw α e θ w × I 0 √ 1 -α e c z c z α e θ w P CL eq (c, c ) = e - (c 2 x -2(1-α t )cxc x +c 2 x ) 2α t (2-α t )θw 2πθ w α t (2 -α t ) e - (c 2 y -2(1-α t )cy c y +c 2 y ) 2α t (2-α t )θw 2πθ w α t (2 -α t ) |c z c z |e -(c 2 z +c 2 z ) 2αeθw α e θ 2 w × I 0 √ 1 -α e c z c z α e θ w , (9)
have this property. In particular, they are based on two parameters α t , α e , combination of multivariate Gaussian functions, and I 0 as modified Bessel function of first kind and zeroth order. Additionally, the joint PDF P CL eq (c, c ) also satisfies the symmetry condition (7).

It is interesting to note that a general class of 1D nonparametric collision model can also be constructed based on Eq. ( 8) and collision data. The uncorrelation hypothesis reduces the model from 3D (or 6D in data space c , c) to 1D (or 2D in each of 3 data spaces c x , c x , c y , c y , c z , c z ), which greatly simplifies the kernel construction. Using kernel density estimation (KDE) techniques, the constructed collision model has good performances when determining velocity slip and temperature jump coefficients like for the CH 4 -graphite or Ar-Au systems [8,22]. However, for nearly perfect elastic-reflection system like the He-Au system [22], the independent assumption of Eq. ( 8) can overestimate the reflective energy and can cause errors in the thermal conductance coefficient. In what follows, we shall explore the GM model, a popular and general purpose machine learning technique, which is not subject to above limitations. The model shares some similarities with the CL kernel which is among the best parametric models, by the presence of smooth multivariate Gaussian functions but it is much more powerful and flexible by means of superposition and unlimited number of the parameters used. Furthermore, the algorithms to identify those parameters from data are well developed in literature and in commercial software. Thus the efforts to construct an accurate collision model will be significantly minimized.

B. Preprocessing of training data and Gaussian mixture model (GM)

As mentioned earlier, the data are collected from MD simulations of independent gas-wall collisions at a given temperature T w . The pre-collision velocity is drawn from equilibrium distribution, i.e P (c ) = P eq (c ) and P (c, c ) = P eq (c, c ) and we can fit the GM model with data, the realizations of P eq (c, c ). To achieve the best performance of the method, we propose to precondition the data in the following way -Symmetrization of the data: For each realization (c, c ), a virtual copy (-c , -c) is also added to the existing data. This step will double the size of the training data and ensure the symmetry of P eq (c, c ) as shown in (7).

-Transformation of the random variables c and c : This is because of the support constraints c z > 0 and c z < 0 are not compatible with gaussian functions. Specifically, while c x , c y , c x , c y have Gaussian distributions, c z , c z have Rayleighian distributions. In this case, we propose to adopt a transformation T that changes the distributions of c z , c z to the same type as c x , c y , c x , c y . For example

T (U ) = √ 2θerf -1 1 -2 exp - U 2 2θ , T -1 (U ) = -2θ ln 1 2 - 1 2 erf U √ 2θ , x R =      c x c y T (c z )      , x I =      c x c y T (c z )      , c =      x R1 x R2 T -1 (x R3 )      , c =      x I1 x I2 -T -1 (x I3 )      , (10)
where x R and x I are two vector which are used to record the reflective value or incident value respectively. We note that to map the Gaussian to the Rayleigh distribution, we compute the CDFs of Gaussian distribution and of Rayleigh distribution

F G (T ) = 1 2 1 + erf T √ 2θ , F R (U ) = 1 -exp -U 2 / (2θ) . (11)
and use the relation F G (T ) = F R (U ) to recover the transformation (10) concerning c z , c z .

The GM model is a probabilistic model that assumes all the data points are generated from a mixture of Gaussian distributions [START_REF] Sung | Gaussian mixture regression and classification[END_REF]. More importantly, it is general enough to estimate complex PDF function. Specifically, the GM estimator of P (x R , x I ) can be written as a superposition of M Gaussian functions with weights ϕ 1 , ϕ 2 , .., ϕ M as follows

P (x R , x I ) P GM (x R , x I ) = M m=1 ϕ m p m (x R , x I ), p m (x R , x I ) = N (X|µ m , Σ m ) (12)
Each Gaussian p m (x R , x I ) has its own parameters µ m and Σ m as mean vector and covariance matrix

N (X|µ m , Σ m ) = 1 (2π) d/2 1 |Σ m | 1/2 exp - 1 2 (X -µ m ) T Σ -1 m (X -µ m ) , X =   x R x I   , µ m =   µ mR µ mI   , Σ m =   Σ mRR Σ mRI Σ mIR Σ mII   , ( 13
)
where d is the dimension of variables X. Comparing with the parametric CL model ( 9), the bivariate gaussians also appear in the independent reflection process x and y where the covariance parameters (1 -α t ) are connected to the tangential momentum accommodation coefficients. In the case M = 1, the GM model also shares some strong similarities with Nocilla [START_REF] Collins | [END_REF] model based on a drifted Maxwellian for reflected velocities. However, the application of the GM model on the preprocessed collision data rather than the original data, as mentioned previously has fully eliminated the negative reflective velocity c z that the drifted Maxwellian can not. Furthermore, as seen in (12,[START_REF] Hollander | Nonparametric statistical methods[END_REF], the GM model is much more general. Constituted of multiple multivariate Gaussian functions, it does not rely on independent reflection hypothesis and contains numerous parameters to fit with data.

Those are ϕ m , µ m and Σ m with m = 1, 2, .., M which can be effectively determined using the EM algorithm. (see Appendix 1) [START_REF] Calinon | Robot Programming by Demonstration: A Probabilistic Approach[END_REF][START_REF] Dempster | [END_REF][START_REF] Hastie | The elements of statistical learning[END_REF].

The scattering kernel is the conditional probability function of reflective velocity P (x R |x I )

given incident velocity and this function can also be written as a mixture form [START_REF] Sung | Gaussian mixture regression and classification[END_REF][START_REF] Rasmussen | [END_REF]:

P GM (x R |x I ) = M m=1 φm (x I )p m (x R |x I ). ( 14
)
with the following weights

φm (x I ) = ϕ m p m (x I ) M m=1 ϕ m p m (x I ) . ( 15
)
The marginal and the conditional distributions appearing in equations ( 14) and ( 15) are calculated with the formula

p m (x I ) = p m (x R , x I )dx R = N (x I |µ mI , Σ mII ), (16)
and

p m (x R |x I ) = p m (x R , x I ) p m (x I ) = N (x R |µ mR|I , Σ mR|I ), µ mR|I = µ mR + Σ mRI Σ -1 mII (x I -µ mI ), Σ mR|I = Σ mRR -Σ mRI Σ -1 mII Σ mIR . (17)
In the next section, we will study the performance of the GM model in comparison with the NP and CL kernels.

C. Evaluation of model performance and comparisons

Generally, the criteria to judge the quality of a statistical model is the accuracy of fitting with the given data and more importantly, the ability to generate new data predicting quantities consistent with the realistic system. For evaluation, we propose to compute physical based parameters including accommodation coefficients, the reflective deviation angle and boundary conditions coefficients β1 and β2 at the interface using the GM model and compare with the same parameters issued from NP, CL kernels and the realistic MD simulations.

As mentioned earlier, accommodation coefficients appear in parametric collision models as constants. However, for realistic atomistic systems, these constants may not be properly defined. In literature, there exist several methods to identify those coefficients but they cannot guarantee the uniqueness of the obtained results (see examples in Appendix 2).

Among those works, we consider the method based on correlation analysis [8,[START_REF] Spijker | [END_REF]. In this case, the accommodation coefficients are defined via the correlation parameters of the data In what follows, we assume that those parameters are obtained from the equilibrium distribution P eq (c, c ) which is also the basis of our model. For the sake of simplicity, we continue to call them accommodation coefficients despite knowing that from definition [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF], they are more statistical than physical quantities.

α 1 = 1 - c x c x -c x c x c 2 x -c x 2 , α 2 = 1 - c 2 c 2 -c 2 c 2 c 4 -c 2 2 α n = 1 - |c z |c z -|c z | c z c 2 z -c z 2 , α e = 1 - c 2 z c 2 z -c 2 z c 2 z c 4 z -c 2 z 2 (18)
In addition to the correlation analysis of velocities in association with the accommodation coefficients described above, we also examine in details the statistical relations between the incident γ , reflection γ and deviation η angles as defined in Fig. 2. Other physical parameters to be checked are the interfacial coefficients β1 and β2 . From the kinetic consideration of the terms involved in equations ( 1) and ( 2), β1 and β2 can be computed by the expressions

β1 = π/2θ c x -c x c x /|c z | + c x /|c z | , β2 = 9π/32θ c 2 -c 2 c 2 /|c z | + c 2 /|c z | -3θ w 1/|c z | + 1/|c z | (19)
Because of the non-equilibrium origin of the expressions, those coefficients cannot be determined using the equilibrium distribution for the gas at the same temperature and the same velocity as the wall, i.e P (c ) = P eq (c |θ = θ w , c x = 0), otherwise we shall encounter numerical issues. Instead, for the determination of β1 , the equilibrium distribution with a non zero mean velocity c x = 0 and the temperature of the wall θ = θ w is used as the Another criteria to evaluate the statistical performance of different scattering kernels is based on the relative entropy, the Kullback-Leibler divergence (D KL ) [30]. This quantity is a measure of the differences between the probabilities P and Q. Defined as the expectation of the logarithmic difference, its discrete form is written as:

D KL (P Q) = j P (x j ) log P (x j ) Q(x j ) , (20)
To compute P (x j ) and Q(x i ), we divide the data space into n bins and choose a Parzenwindow size h, then the probability density of data generated by MD simulation or a scattering kernel is estimated. The density in each bin can be written as:

P (x j ) = N j j N j , (21)
where, x j represents the center of bin j, N j is the number of data in the window of length h centered on x j . We note that the above formula can be applied to any distribution of random variables x, for example velocity, kinetic energy, collision angles etc..

III. APPLICATIONS TO REAL GAS-WALL SYSTEMS

A. Construction of scattering kernels from Molecular Dynamics simulation data

Based on the same systems as Refs. [22,31], we shall carry out more detailed analysis of the MD simulation results, in order to construct and evaluate statistical models of collision.

We simulate the process of He or Ar gas atoms impacting smooth/rough Au atomic surface and investigated the gas-solid friction and thermal conductance. The rough surface is generated by deposing randomly atoms on the smooth surface. The He or Ar molecules are inserted one by one in the collision zone, only after the molecule interacts with the wall and goes out of the collision zone, another molecule is inserted in the zone from a random posi-Pair tion with a random velocity. The MD simulations are performed by LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [START_REF] Plimpton | Large-scale atomic/molecular massively parallel simulator[END_REF]. Each collision simulation was run with one gas atom at a time. The simulation parameters will be described briefly as follows as follows:

A[eV] B[eV Å-1 ] C[eV Å] α[ Å-1 ] γ[ Å-2 ] C 6 [eV Å6 ] C 8 [eV Å8 ] β[ Å-
• temperature of the Au surface: 300 K, controlled by the Nosé-Hoover thermostat

• to compute β 1 we apply a average velocity to the incident gas atom: ranging from -0.5 Å/ps to 0.5 Å/ps (with a step of 0.25 Å/ps)

• to compute β 2 we apply a initial temperature to the incident gas: ranging from 100 K to 400 K (with a step of 100 K)

The highly optimized version of multibody EAM potential [START_REF] Sheng | [END_REF] is employed for the description of the interaction between the Au atoms. The Au-Ar and Au-He pairwise potentials, which are recently determined from ab-initio studies [22,31], are implemented in LAMMPS and adopted in this work. The functional form of the pair potentials V g-Au (R) is given by

V g-Au (R) = A + BR + C R e (-αR-γR 2 ) + n=6,8 C n R n f n (βR) . ( 22
)
where R is the distance between the gas atom ( g standing for He or Ar) and a gold atom (Au). The parameters A, B, C, α, and γ are associated with the repulsive part and the remaining parameters correspond to the attractive part of the potential. The function f n is the damping function of Tang and Toennies [34],

f n (x) = 1 -e -x 1 + x + ... + x n n! . ( 23
)
The numerical values of the pairwise potential parameters are given in Tab.I

In this paper, to construct the scattering kernel, we use the velocity data of the MD simulations under the equilibrium state, in which the temperature of gas flow is equal to the temperature of the atomic wall (300 K) and the mean velocity of the flow equals to 0. We can calculate the accommodation coefficients for the CL kernel using Eq. ( 18). For the NP kernel, we followed the model of Ref. [22], the He velocities range from -40 to 40 Å/ps, and the Ar velocities range from -20 to 20 Å/ps. The incident and reflective velocity spaces were discretized into 100 and 1000 intervals, respectively, the Parzen-window size h = 1 Å/ps. The probability that a reflective velocity does not belong to a category is determined by interpolation. In the GM model, the number of mixtures M is a free parameter which needs to be chosen before the learning. To avoid the risk of overfitting or underfitting, there is an optimal value of M for a given MD data size. We tested the number of Gaussian functions from M = 1 to 512 with different initial covariance matrices and mean vectors.

To find the optimal M we calculate the energy accommodation coefficient α 2 using the reflective velocity generated by GM model and Eq. ( 18). In Fig. 3, we found that when M is equal to 64 (2 6 ), the performance of the GM model is the best and much better than CL and NP models. So the GM model with M = 64 has a good accuracy and a reasonable computation time and is used in the present study. We then used each kernel, CL, NP and GM, to generate new incident-reflective velocities (about 10 5 data points) based on velocities under equilibrium or non-equilibrium states. We investigated the characteristics of the scattering kernels according to the following aspects:

• the accommodation coefficients and the momentum/energy distribution

• the incident, reflective, and deflection angles

• the surface friction and the thermal conductance

B. Accommodation coefficients

We decomposed the velocity into the normal and the tangential components and then studied the correlations between the incident-reflective velocities and the kinetic energy.

Using Eq. ( 18), we can calculate TMAC (α 1 ), NMAC (α n ), and EAC (α 2 ) from incidentreflective velocity data. Tab. II presents TMAC, NMAC and EAC of He/Ar gas on the smooth/rough Au surface using the atomic wall and the different statistical models of collision (GM, NP, CL). From this table, we find that TMAC and NMAC obtained by the different kernels are very close to the atomic simulation results. However, significant differences are observed for the coefficient EAC. The result of the GM kernel is consistent with the atomic model, but the NP and CL kernels considerably overestimate EAC, especially for the He-Au rough surface system. We note that, unlike Ar, He is a special gas which has a very small atomic mass and its interaction with the gold surface is very weak (see ref. [22]). These points can explain the nearly perfect elastic collision of an helium atom on a Au wall, even for the rough atomic surface. In other words, while the atomic motion can change its direction after collisions due to the surface corrugation (e.g via the influences of the atomic roughness, the thermal motion, and the energy landscape), the atomic kinetic energy is almost conserved.

Using the density estimator presented in Eq. ( 21), we can estimate the joint probability of the incident-reflective velocities, and of the kinetic energy. The kinetic energy distributions determined by the atomic wall use and by different scattering kernels are plotted in Fig. 4.

The energy clouds of the atomic and GM models are visibly narrower than those of the CL and NP kernels. From the column figures of both atomic and GM surface, we can see that roughness does not significantly change EAC (16.7% for He, 5.9% for Ar). However, this trend cannot be captured by both NP and CL kernels. These phenomenas are further illustrated in Fig. 5. The KL divergence D KL between different statistical models and the atomic model was calculated. We can notice that for tangential and normal momentums, the D KL of the GM and NP kernels is at the same level.

Gas

For the CL kernel, although the D KL for c x and c z are larger than those of the other two kernels, its good performance in determining TMAC and NMAC (Tab. II) indicates that CL kernel is still a good kernel to generate a single velocity component. Meanwhile, Fig. 5b presents anomalies for the kinetic energy. The energy D KL of NP and CL kernels are close to 1 for the He/rough Au case, meaning that the incident-reflective kinetic energy distributions determined by the NP and CL kernels differ greatly from the results of the atomic simulations. Considering the good performance of CL and NP kernels for the momentum reflection, and less for the kinetic energy, it is suggested that the independent assumption of Eq. ( 8) may oversimplify the true behavior especially for nearly perfect elastic systems like He/Au. This deficiency has been properly covered by the GM model. We note that the EAC value difference between the statistical models is due to the correlation parameter As a result, the accuracy of P eq (c, c ) governs the accuracy of the correlation parameter and EAC. While it is not trivial to connect the elastic collision phenomena to the property of the probability function P eq (c, c ) or P (c|c ), we can conclude mathematically that using the decomposition postulate (8) like CL and NP models oversimplifies the real behavior of P eq (c, c ). In other words, due to the exchange between the velocity components in order to conserve the energy, each of post-collision velocity components is expected to depend on all the pre-collision components.

c 2 c 2 = c 2 c 2 P eq (c, c )dcdc

C. Collision angles

To further study the performance of each scattering kernel, the probability densities of deviation angle η and incident-reflective vertical angles γ -γ determined by atomic simulations and scattering kernels were analyzed. The geometric definitions of η, γ , γ are illustrated in Fig. 2.

Angle η can be calculated from velocity components c x , c y , c x , c y and measures the deviation of the gas atom reflection from the incident direction in xy-plane. The incident angle γ and reflective angle γ depend on the incident, c , and reflective, c, velocities, respectively. By calculating the correlation between γ and γ, we can understand the aberration of the reflection in the normal direction.

Figure 6 shows the probability of the reflective direction on the tangential plane obtained from atomic simulation data, GM, NP and CL kernels under equilibrium state. Numerical results show that the rough surface can increase the deviation of the reflective direction with respect the incident direction in xy-plane, especially the percentage of backscattering (π/2 < η < 3π/2). As the result, the overall scattering is more diffusive and the probability contour is shifted towards and envelope the origin. In contrast, the contour of reflection on smooth surfaces is narrower. From analysis of gas atoms, the adhesion between Au wall and Ar gas atoms is stronger than for the He-Au system, so that the reflection direction has a greater probability of deviating from the specular direction. From Fig. 6, we can find that FIG. 6. The probability density of the reflective deviation angle η determined by the atomic simulations, the Gaussian mixture (GM), the nonparametric (NP), and the Cercignani-Lampis (CL) kernels in polar coordinate system. The distance from the curve to the origin, the radius, indicates the probability of associated to angle η, which coincides with the angle coordinate of the system.

both CL and NP kernels overestimate the diffusivity of the reflective direction. Only the probability contour of the GM kernel is identical with the one of the the atomic simulations.

Figure 7 shows the joint probability density of γ and γ. When the deviation angle is small, the data population becomes narrow and shrinks toward the diagonal line. Through these cloud figures, we can also confirm that the roughness of the wall surface and the bonding properties between the gas and the wall affect the reflective direction of the gas atom. To quantitatively compare the differences between the three types of kernels and the atomic simulations, we calculate the KL divergence between the scattering kernels and the atomic wall results. The KL divergence of the reflective angle probability is presented in Fig. 8. The results of the figure confirm the differences observed in Figures 6 and7. The KL divergence of the CL kernel is over 3 times larger than that of the GM kernel, especially for He-Au rough surface system. Basically, the CL kernel incorrectly estimates the distribution of the reflection direction. The results of the NP kernel for the Ar-Au system are better

D. Interfacial friction and thermal conductance

The surface friction and thermal conductance properties can be characterized by the velocity slip coefficient, β1 , and the temperature jump coefficient, β2 , respectively. These two coefficients can be determined by the accommodation coefficients determined previously and relation [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction, and diffusion in gases[END_REF] or by a direct method based on equation ( 19) and the non-equilibrium collision simulation data. This step permits to test the ability of the GM model in terms of generating new data and reproducing these important physical properties for real systems. We note that if the final objective is to construct kernels usable in kinetic simulation methods, this step is crucial and guarantees the reliability of the results.

The results for two coefficients are reported in Tab. III. The coefficient β1 depends on the accommodation momentums and is sensitive to the surface roughness. It means that the friction significantly increases β1 . These effects are observed for both Ar/Au et He/Au systems. From smooth surface to rough surface, the β1 value increases about 34 times for He-Au system and about 2 times for Ar-Au system. In contrast, the contribution of the surface roughness to the increase of the thermal conductance is relative small, about 20%

for He and 9% for Ar.

The numerical tests confirm again the superior performance of the GM model. This is the statistical model that produces the best results when compared with the atomic model. It is surprising to find that the simple analytical relation between β and α (Eq. ( 3)) works relatively well and yields values of the same order as MD and GM models. It is less accurate than the GM model but better than other models in terms of evaluating the thermal coefficient β2 . In Appendix 2, we also carry out theoretical investigation of CL models. Specifically, we compute analytically the slip and jump coefficients β1 and β2 for the CL model with parameters α t , α e using the same approaches as the MD model. Results show that the obtained expressions for β1 and β2 in terms of α t , α e are different but the numerical differences are very small for the whole admissible range (0, 1) of α t and α e . However as discussed previously, the coefficients α 1 and α 2 defined by Eq. ( 18) are more of statistical nature and furthermore α 2 is not associated to any parametric statistical models in literature (except for the Maxwell model). This fact will limit the application of those coefficients in particle and/or statistical based simulation methods like MD or DSMC. In contrary, the GM model is not subject to this limitation. Expectation maximization (EM) algorithm is an iterative method to find the maximum likelihood of parameters in statistical models, which is introduced by Dempster et. al. [START_REF] Dempster | [END_REF].

The parameter set of the GM model is noted as Θ = (θ 1 , θ 2 , ...θ M ) and the m th parameter is θ m = (ϕ m , µ m , Σ m ). Noting samples as X, i.e X i corresponds to the i th training data of total N training data, the likelihood can be written as

L(X|Θ) = N i=1 p(X i |Θ), p(X i |Θ) = M m=1 ϕ m N (X i |µ m , Σ m ) (24)
To determine the parameter Θ, we can maximize the log-likelihood function:

Θ = arg max Θ log L(X|Θ) = arg max Θ N i=1 log p(X i |Θ), (25)
and the sun of log-likelihood function in Eq. ( 25) can be replaced by the expectation of log-likelihood function:

Θ = arg max Θ E (log p(X i |Θ)) . (26)
Then, we can use an iterative algorithm to determine the parameter Θ of the GM model. The EM algorithm can be concluded in four steps: initial state, E-step, M-step and convergence test:

1. Initial state: take the initial value θ

m ) for each Gaussian distribution in the GM model.

2. E-step: use the parameters of the (r -1) th iteration to calculate the mixture weights w im of each sample X i for the (r) th iteration:

µ (r) m = i w im X i i w im Σ (r) m = i w im (X i -µ (r) m )(X i -µ (r) m ) T i w im ϕ (r) m = 1 n n i=1 w im . (28)
4. convergence test: conclude the iterative process if the following conditions are met:

1 -E log p(X i |Θ (r) E (log p(X i |Θ (r-1) )

< ε, (29)
where, ε is a threshold. Repeat the steps 2 to 4 until the expectation of the log-likelihood function converge.

The advantage of the EM algorithm is that it guarantees local convergence, but the drawback is that the algorithm is sensitive to the initial state. In practice, the k-means algorithm [START_REF] Calinon | Robot Programming by Demonstration: A Probabilistic Approach[END_REF] and the iterative pairwise replacement algorithm (IPRA) [START_REF] Sung | Gaussian mixture regression and classification[END_REF] have a good performance in giving an initial state for the GM model. In addition, we can also use multiple random initial states and then compare the distributions of GMs obtained with different initial states.

Comparisons of interfacial coefficients associated to the CL model by different methods

We consider the CL kernel whose analytical expression is given by ( 9) and examine the coefficients α 1 , α 2 , β1 , β2 computed from the correlation analysis and by the non-equilibrium method.

Using the equilibrium distribution P CL eq (c, c ), the moments f (c, c ) in correlation analysis can be calculated by f (c, c ) = f (c, c )P CL eq (c, c )dcdc (30) which yields the explicit expressions

c x = c x = 0, c x c x = (1 -α 1 )θ w , c 2 x = c 2 x = θ w , c 2 z = |c 2 z | = 2θ w , c 2 = c 2 =
under assumption of the validity of the relation βi = α i /(2 -α i ) with i = 1, 2.

Regarding the non-equilibrium method, we impose the input distribution P (c ) at a different temperature or velocity from the wall and use the scattering kernel to compute P (c) P (c) = P CL (c|c )P (c )dc (34) To determine β1 and β2 with [START_REF] Allen | Computer Simulation of Liquids[END_REF], we need to compute the moments

FIG. 1 .

 1 FIG. 1. Procedure of reconstructing the scattering kernel from MD simulations of collisions.

While α 1

 1 and α 2 are the tangential momentum accommodation coefficient (TMAC) and the energy accommodation coefficient (EAC), respectively (see equation (3)), α e is the normal momentum accommodation coefficients (NMAC) and α e is the normal energy accommodation coefficients (NEAC). Those expressions are derived from the linear regression of collision data (c x , c x ), (c 2 , c 2 ), (|c z |, c z ) and (c 2 z , c 2 z ), respectively. Being special composite parameters depending on the moments of P (c, c ), they can be sensitive to the choice of P (c, c ) and P (c ).

FIG. 2 .

 2 FIG.2. A simple scheme of the incident-reflective velocities, c -c, the incident-reflective vertical angles, γ -γ, and the deviation angle η.

  input gas distribution P (c ) = P eq (c |θ = θ w , c x = 0) and for β2 , we used the equilibrium distribution with zero mean velocity c x = 0 and a temperature different from the wall temperature θ = θ w or P (c ) = P eq (c |θ = θ w , c x = 0). We also highlight the importance of transforming the variables c z and c z in the data preconditioning step before training the GM model. This step guarantees the vanishing probability at c z = 0 and c z = 0 like the Rayleigh distribution and the convergence of the rational moments containing 1/c z and 1/c z . Unlike the previous analysis of accommodation coefficients and reflective angle based on the existing collision data at equilibrium, this method requires new realizations, i.e using the new input P (c ) to generate new data P (c , c) with the wall model P (c|c ) (parametric/nonparametric kernels, GM model or atomic wall).

FIG. 3 .

 3 FIG. 3. Number of Gaussian vs. α 2 determined by statistical models and MD

FIG. 4 .

 4 FIG.4. The joint probability density of incident-reflective kinetic energy with unite Å2 /ps 2 (the prefactor m/2 is removed for simplicity) determined by the atomic simulations, the Gaussian mixture (GM), the nonparametric (NP), and the Cercignani-Lampis (CL) kernels. The horizontal axis indicates the incident energy c 2 and the vertical axis indicates the reflective energy c 2 . The red line is the least-square linear regression of kinetic data, its slope is equal to 1-α. When the red line is close to the diagonal dashed line, the reflection is close to be specular (elastic collision). When the red line is close to the horizontal dashed line, the reflection is diffusive.

2 FIG. 5 .

 25 FIG.5. The Kullback Leibler divergence D KL (Atomic||Kernel) of momentum and kinetic energy determined by the atomic simulations, the Gaussian mixture (GM) kernel, the nonparametric (NP) kernel, and the Cercignani-Lampis (CL) kernel.

FIG. 7 .

 7 FIG. 7. The probability density of the incident-reflective vertical angles determined by the atomic simulations, the Gaussian mixture (GM), the nonparametric (NP), and the Cercignani-Lampis (CL) kernels. The horizontal axis indicates the incident angle γ and the vertical axis indicates the reflective angle γ. The incident and reflective angles range from 0 to 90 •

FIG. 8 .

 8 FIG. 8. The Kullback-Leibler divergence D KL (Atomic||Kernel) derived form Fig. 6 and 7.

FIG. 9 .APPENDIX 1 .

 91 FIG. 9. The average relative error of the different schemes, GM, NP, CL, and α, comparing with the atomic simulations.

  -step: use the w im calculated by E-step to determine the new parameters µ

f 2 α

 2 (c ) = f (c )P (c )dc , f (c) = f (c)P (c)dc (35)yielding the final results (after considering θ θ w )β1 = α t 2 -α t , β2 = 3 e + α t (2 -α t ) (1 -α e ) + 2(1 -α t ) 2 + 3 .(36)

TABLE I .

 I The parameters of the He-Au and Ar-Au pairwise interaction potentials[22,31].

TABLE II .

 II The tangential momentum accommodation coefficient (TMAC) α 1 , the normal momentum accommodation coefficient (NMAC) α n , and the energy accommodation coefficient (EAC) α 2 , computed using different statistical models: GM, NP, CL, and the atomic (MD) model. The construction of the statistical models NP, GM and CL from MD data is discussed in sec. IIIA. The accommodation coefficients associated to the statistical models are computed based on generation of new data and Eq. (18). It is noted that the CL model uses tangential momentum accommodation coefficient (TMAC) and normal energy accommodation coefficient (NEAC) computed by MD as parameters. Consequently, numerical TMAC values of CL model are extremely close to MD results.

TABLE III .

 III The friction coefficient ( β1 ) and the thermal conductance ( β2 ) determined by different schemes.

  4θ w , c 2 c 2 = 4θ w 2 α t 2 -2α t -α e + 6 , c 4 = c 4 = 24θ 2 w . (31)Substituting the above results in[START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF], we obtain coefficients α 1 , α 2α 1 = α t , α 2 = 1 2 [α e + α t (2 -α t )] + α t (2 -α t ) 2 + (1 -α e ) + (1 -α t ) 2 ,

				(32)
	and the associated interfacial coefficients, β1 , β2
	β1 =	α t 2 -α t	, β2 =	α e

ACKNOWLEDGMENTS

We thank Dr. Brian Staber for his advice regarding probability calculation.

IV. CONCLUSION

In this paper, we present the application of unsupervised learning techniques in the construction of wall-gas models and the computation of interfacial properties like friction and thermal conductance. The data obtained from the Molecular Dynamics simulations of gas-wall collision is used to train the scattering kernels, i.e statistical collision models that predict the post-collision velocity from a given pre-collision velocity. Due to the common similarities with the well-known Cercignani-Lampis (CL) model in kinetic theory and the good overall performance in general data treatment, the GM model is chosen for this purpose. We also noted that in literature, there also exists different expressions of coefficients β1 , β2 (see Ref. [5,16])

Despite using non equilibrium gas distribution as input distribution, the derivation of these equations is different from the non equilibrium method used in the present work. It is based on the moment balance at the wall and the temperature, velocity and stress and heat flux are parameters of input distribution (Chapman-Enskog) and are not computed from both input and output distributions as in [START_REF] Allen | Computer Simulation of Liquids[END_REF].

A general remark can be made here that while all the methods agree in terms of slip effects (coefficient β1 ), analytical expressions for the interfacial conductance are different (coefficient β2 ) . It is suggested that these discrepancies are due to the fact that the EAC coefficient is undefined for the CL wall kernel [8], except for the diffusive wall α e = α t = 1.

Nevertheless, the discrepancies between the methods are not significant, maximum 6% at the extreme contrast cases α t = 1, α e = 0 or α t = 0, α e = 1 (see Fig. 10). with the same CL models: the equilibrium method (33), the nonequilibirum method 1 in this work (36), the nonequilibirum method 2 in literature (37).