Huy-Dung Nguyen

Michaël Clément

Boris Mansencal

Pierrick Coupé

Deep Grading based on Collective Artificial Intelligence for AD Diagnosis and Prognosis

Keywords: Deep Grading, Collective Artificial Intelligence, Generalization, Alzheimer's disease classification, Mild Cognitive Impairment

Accurate diagnosis and prognosis of Alzheimer's disease are crucial to develop new therapies and reduce the associated costs. Recently, with the advances of convolutional neural networks, methods have been proposed to automate these two tasks using structural MRI. However, these methods often suffer from lack of interpretability, generalization, and can be limited in terms of performance. In this paper, we propose a novel deep framework designed to overcome these limitations. Our framework consists of two stages. In the first stage, we propose a deep grading model to extract meaningful features. To enhance the robustness of these features against domain shift, we introduce an innovative collective artificial intelligence strategy for training and evaluating steps. In the second stage, we use a graph convolutional neural network to better capture AD signatures. Our experiments based on 2074 subjects show the competitive performance of our deep framework compared to state-of-the-art methods on different datasets for both AD diagnosis and prognosis.

Introduction

The first cognitive symptoms of Alzheimer's disease (AD) appear right after the morphological changes caused by brain atrophy [START_REF] Gordon | Spatial patterns of neuroimaging biomarker change in individuals from families with autosomal dominant Alzheimer's disease: A longitudinal study[END_REF]. Those changes can be identified with the help of structural magnetic resonance imaging (sMRI) [START_REF] Bron | Standardized evaluation of algorithms for computer-aided diagnosis of dementia based on structural MRI: The CADDementia challenge[END_REF]. Recently, with the advances of convolutional neural networks (CNN), methods have been proposed for automatic AD diagnosis using sMRI. Despite encouraging results, current deep learning methods suffer from several limitations. First, deep models lack transparency in their decision-making process [START_REF] Nigri | Explainable deep CNNs for MRI-based diagnosis of Alzheimer's disease[END_REF][START_REF] Zhang | An explainable 3D residual self-attention deep neural network for joint atrophy localization and Alzheimer's Disease diagnosis using structural MRI[END_REF]. Therefore, this limits their use for computer-aided diagnosis tools in clinical practice. Second, for medical applications, the generalization capacity of classification models is essential. However, only a few works have proposed methods robust to domain shift [START_REF] Hosseini-Asl | Alzheimer's disease diagnostics by adaptation of 3D convolutional network[END_REF][START_REF] Wachinger | Domain adaptation for Alzheimer's disease diagnostics[END_REF]. Third, current CNN models proposed for AD diagnosis and prognosis still perform poorly [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF]. Indeed, when properly validated on external datasets, current CNN-based methods perform worse than traditional approaches (i.e., standard linear SVM).

In this paper, to address these three major limitations, we propose a novel interpretable, generalizable and accurate deep framework. An overview of our proposed pipeline is shown in Figure 1. First, we propose a novel Deep Grading (DG) biomarker to improve the interpretability of deep model outputs. Inspired by the patch-based grading frameworks [START_REF] Coupé | Scoring by nonlocal image patch estimator for early detection of Alzheimer's disease[END_REF][START_REF] Hett | Graph of brain structures grading for early detection of Alzheimer's disease[END_REF][START_REF] Tong | A novel grading biomarker for the prediction of conversion from mild cognitive impairment to alzheimer's disease[END_REF], this new biomarker provides a grading map with a score between -1 and 1 at each voxel related to the alteration severity. This interpretable biomarker may help clinicians in their decision and to improve our knowledge on AD progression over the brain. Second, we propose an innovative collective artificial intelligence strategy to improve the generalization across domains and to unseen tasks. As recently shown for segmentation [START_REF] Coupé | AssemblyNet: A large ensemble of CNNs for 3D whole brain MRI segmentation[END_REF][START_REF] Kamraoui | Towards broader generalization of deep learning methods for multiple sclerosis lesion segmentation[END_REF], the use of a large number of networks capable of communicating offers a better capacity for generalization. Based on a large number of CNNs (i.e., 125 U-Nets), we propose a framework using collective artificial intelligence efficient on different datasets and able to provide accurate prognosis while trained for diagnosis task. Finally, we propose to use a graph-based modeling to better capture AD signature using both inter-subject similarity and intrasubject variability. As shown in [START_REF] Hett | Graph of brain structures grading for early detection of Alzheimer's disease[END_REF], such strategy improves performance in AD diagnosis and prognosis.

In this paper, our main contributions are threefold:

-A novel deep grading biomarker providing interpretable grading maps.

-An innovative collective artificial intelligence strategy robust to unseen datasets and unknown tasks. -A new graph convolutional network (GCN) model for classification offering state-of-the-art performance for both AD diagnosis and prognosis.

2 Materials and method

Datasets

The data used in this study, consisting of 2074 subjects, were obtained from multiple cohorts: the Alzheimer's Disease Neuroimaging Initiative (ADNI) [START_REF] Jr | The Alzheimer's disease neuroimaging initiative (ADNI): MRI methods[END_REF], the Open Access Series of Imaging Studies (OASIS) [START_REF] Lamontagne | OASIS-3: Longitudinal Neuroimaging, Clinical, and Cognitive Dataset for Normal Aging and Alzheimer Disease[END_REF], the Australian Imaging, Biomarkers and Lifestyle (AIBL) [START_REF] Ellis | The Australian Imaging, Biomarkers and Lifestyle (AIBL) study of aging: Methodology and baseline characteristics of 1112 individuals recruited for a longitudinal study of Alzheimer's disease[END_REF], the Minimal Interval Resonance Imaging in Alzheimer's Disease (MIRIAD) [START_REF] Malone | MIRIAD-Public release of a multiple time point Alzheimer's MR imaging dataset[END_REF]. We used the baseline T1-weighted MRI available in each of these studies. Each dataset contains AD patients and cognitively normal (CN) subjects. For ADNI1 and AIBL, it also includes mild cognitive impairment (MCI), the early stage of AD composed of abnormal memory dysfunctions. Two groups of MCI are considered: progressive MCI (pMCI) and stable MCI (sMCI). The definition of these two groups is the same as in [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF]. Table 1 summarizes the number of participants for each dataset used in this study. During experiments, AD and CN from ADNI1 are used as training set and the other subjects as testing set.

Preprocessing

All the T1w MRI are preprocessed using the following steps: (1) denoising [START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF],

(2) inhomogeneity correction [START_REF] Tustison | N4ITK: Improved N3 bias correction[END_REF], (3) affine registration into MNI space (181 × 217 × 181 voxels at 1mm × 1mm × 1mm) [START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF], (4) intensity standardization [START_REF] Manjón | Robust MRI brain tissue parameter estimation by multistage outlier rejection[END_REF] and (5) intracranial cavity (ICC) extraction [START_REF] Manjón | NICE: Non-local Intracranial Cavity Extraction[END_REF]. After that preprocessing, we use AssemblyNet [START_REF] Coupé | AssemblyNet: A large ensemble of CNNs for 3D whole brain MRI segmentation[END_REF] to segment 133 brain structures (see Figure 1). The list of structures is the same as in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. In this study, brain structure segmentation is used to determine the structure volume (i.e., normalized volume in % of ICC) and aggregate information in the grading map (see Section 2.3 and Figure 1).

Deep Grading for disease visualization

In AD classification, most of deep learning models only use CNN as binary classification tool. In this study, we propose to use CNN to produce 3D maps indicating where specific anatomical patterns are present and the importance of structural changes caused by AD.

To capture these anatomical alterations, we extend the idea of the patchbased grading (PBG) framework [START_REF] Coupé | Scoring by nonlocal image patch estimator for early detection of Alzheimer's disease[END_REF][START_REF] Hett | Graph of brain structures grading for early detection of Alzheimer's disease[END_REF][START_REF] Tong | A novel grading biomarker for the prediction of conversion from mild cognitive impairment to alzheimer's disease[END_REF]. The PBG framework provides a 3D grading map with a score between -1 and 1 at each voxel related to the alteration severity. Contrary to previous PBG methods based on non-local mean strategy, here we propose a novel DG framework based on 3D U-Nets.

Concretely, each U-Net (similar to [START_REF] Coupé | AssemblyNet: A large ensemble of CNNs for 3D whole brain MRI segmentation[END_REF]) takes a 3D sMRI patch (e.g., 32 × 48 × 32) and outputs a grading map with values in range [-1, 1] for each voxel. Voxels with a higher value are considered closer to AD, while voxels with a lower value are considered closer to CN. For the ground-truth used during training, we assign the value 1 (resp. -1) to all voxels inside a patch extracted from an AD patient (resp. CN subject). All voxels outside of ICC are set to 0.

Once trained, the deep models are used to grade patches. These local outputs are gathered to reconstruct the final grading map (see Section 2.4). Using the structure segmentation, we represent each brain structure grading by its average grading score (see Figure 1). This anatomically driven aggregation allows better and meaningful visualization of the disease progression. In this way, during the classification step (see Section 2.5), each subject is encoded by an n-dimensional vector where n is the number of brain structures.

Collective AI for grading

As recently shown in [START_REF] Bron | Cross-cohort generalizability of deep and conventional machine learning for MRI-based diagnosis and prediction of Alzheimer's disease[END_REF][START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF], current AD classification techniques suffer from a lack of generalization. In this work, we propose an innovative collective artificial intelligence strategy to improve the generalization across domains and to unseen tasks. As recently shown for segmentation [START_REF] Coupé | AssemblyNet: A large ensemble of CNNs for 3D whole brain MRI segmentation[END_REF][START_REF] Kamraoui | Towards broader generalization of deep learning methods for multiple sclerosis lesion segmentation[END_REF], the use of a large number of compact networks capable of communicating offers a better capacity for generalization. There are many advantages to using the collective AI strategy. First, it addresses the problem of GPU memory in 3D since each model processes only a sub-volume of the image. The use of a large number of compact networks is equivalent to a big neural network with more filters. Second, the voting system based on a large number of specialized and diversified models helps the final grading decision to be more robust against domain shift and different tasks.

Concretely, a preprocessed sMRI is decomposed into k × k × k overlapping patches of the same size (e.g., 32 × 48 × 32). During training, for each patch localization in the MNI space, a specialized model is trained. Therefore, in our case (k = 5), we trained m = k × k × k = 125 U-Nets to cover the whole image (see Fig. 1). Moreover, each U-Net is initialized using transfer learning from its nearest neighbor U-Nets in the MNI space, except the first one trained from scratch as proposed in [START_REF] Coupé | AssemblyNet: A large ensemble of CNNs for 3D whole brain MRI segmentation[END_REF]. As adjacent patches have some common patterns, this communication allows grading models to share useful knowledge between them. For each patch, 80% of the training dataset (i.e., ADNI1) is used for training and the remaining 20% for validation. The accuracy obtained on validation set is used to reconstruct the final grading map using a weighted average as follows:

G i = xi∈Pj α j * g ij xi∈Pj α j (1)
where G i is the grading score of the voxel x i in the final grading map, g ij is the grading score of the voxel x i in the local grading patch P j , and α j is the validation accuracy of the patch j. This weighted vote enables to give more weight to the decision of accurate models during the reconstruction.

Graph convolutional neural network for classification

The DG feature provides an inter-subject similarity biomarker which is helpful to detect AD signature. However, the structural alterations leading to cognitive decline could be different between subjects. Indeed, following the idea of [START_REF] Hett | Graph of brain structures grading for early detection of Alzheimer's disease[END_REF], we model the intra-subject variabilities by a graph representation to capture the relationships between several regions related to the disease. We define an undirected graph G = (N, E), where N = {n 1 , . . . , n s } is the set of nodes for the s brain structures and E = s × s is the matrix of edge connections. In our approach, all nodes are connected with each other in a complete graph, where nodes embed brain features (e.g., our proposed DG feature) and potentially other types of external features. Indeed, besides the grading map, the volume of structures obtained from the segmentation could be helpful to distinguish AD patients from CN [START_REF] Hett | Graph of brain structures grading for early detection of Alzheimer's disease[END_REF][START_REF] Tong | A novel grading biomarker for the prediction of conversion from mild cognitive impairment to alzheimer's disease[END_REF]. It is due to the evidence that AD leads to structure atrophy. Age is also an important factor as, within sMRI, patterns in the brain of young AD patients could be similar to elder CN. Indeed, the combination of those features is expected to improve our classification performance. In our method, each node represents a brain structure and embeds a feature vector (DG, V, A) where V and A are respectively the volume of structures and subject's age. Finally, we use the graph convolutional neural network (GCN) [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] as the way to pass messages between nodes and to perform final classification.

Implementation details

First, we downsample the sMRI from 181 × 217 × 181 voxels (at 1mm) to 91 × 109 × 91 voxels to reduce the computational cost, then decompose them into 5 × 5 × 5 overlapping patches of size 32 × 48 × 32 voxels equally spaced along the three axis. For each patch, an U-Net is trained using mean absolute error loss, Adam optimizer with a learning rate of 0.001. The training process is stopped after 20 epochs without improvement in validation loss. We employed several data augmentation and sampling strategies to alleviate the overfitting issue during training. A small perturbation is first created in training samples by randomly translating by t ∈ {-1, 0, 1} voxel in 3 dimensions of the image. We then apply the mixup [START_REF] Zhang | mixup: Beyond Empirical Risk Minimization[END_REF] data augmentation scheme that was shown to improve the generalization capacity of CNN in image classification.

Once the DG feature is obtained, we represent each subject by a graph of 133 nodes. Each node represents a brain structure and embeds DG, volume and age features. Our classifier is composed of 3 layers of GCN with 32 channels, followed by a global mean average pooling layer and a fully connected layer with an output size of 1. The model is trained using the binary cross-entropy loss, Adam optimizer with a learning rate of 0.0003. The training process is stopped after 20 epochs without improvement in validation loss. At inference time, we randomly add noise X ∼ N (0, 0.01) to the node features and compute the average of 3 predictions to get the global decision. Experiments have shown that it helps our GCN to be more stable.

For training and evaluating steps, we use a standard GPU (i.e., NVIDIA TITAN X) with 12Gb of memory.

Experimental results

In this study, the grading models and classifiers are trained using ADNI1 dataset within AD and CN subjects. Then, we assess their generalization capacity in domain shift using AD, CN subjects from ADNI2, AIBL, OASIS, MIRIAD. The generalization capacity in derived tasks is performed using pMCI, sMCI subjects from ADNI1 (same domain) and AIBL (out of domain).

Influence of collective AI strategy. In this part, the DG feature is denoted as DG C (resp. DG I ) when obtained with the collective (resp. individual) AI strategy. The individual AI strategy refers to the use of a single U-Net to learn patterns from all patches of sMRI. We compare the efficiency of DG C and DG I feature when using the same classifier (i.e., SVM or GCN) (see Table 2). These experiments show that using DG C achieves better results in most configurations. When using SVM classifier, we observe a gain of 3.6% (resp. 0.8%) on average in AD/CN (resp. pMCI/sMCI) classification. The efficiency of G C feature is even better with GCN classifier, where a gain of 4.0% (resp. 3.5%) is observed.

Influence of GCN classifier. Besides the DG feature, the intra-subject variabilities are also integrated into our graph representation. Hence, it should be beneficial to use GCN to exploit all this information. In our experiments, GCN outperforms SVM in all the tests using either DG I or DG C feature (see Table 2). Concretely, using DG I feature, we observe a gain of 5.0% (resp. 7.6%) on average for AD/CN (resp. pMCI/sMCI) classification. These improvements are 5.4% and 10.6% when using DG C feature.

Influence of using additional non-image features. Moreover, we analyze the model performance using DG C with the structural volume V and age A as additional node features in our graph representation. By using the combined features, the performance on average in AD/CN and pMCI/sMCI is both improved by 0.3% and 1.4% compared to DG C feature (see Table 2). In the rest of this paper, these results are used to compare with current methods.

Comparison with state-of-the-art methods. Table 3 summarizes the current performance of state-of-the-art methods proposed for AD diagnosis and prognosis classification that have been validated on external datasets. In this comparison we considered five categories of deep methods: patch-based strategy based on a single model (Patch-based CNN [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF]), patch-based strategy based on multiple models (Landmark-based CNN [START_REF] Liu | Landmark-based Deep Multi-Instance Learning for Brain Disease Diagnosis[END_REF], Hierarchical FCN [START_REF] Lian | Hierarchical fully convolutional network for joint atrophy localization and Alzheimer's disease diagnosis using structural MRI[END_REF]), ROI-based strategy based on a single model focused on hippocampus (ROI-based CNN [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF]), subject-based considering the whole image based on a single model (subjectbased CNN [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF], 3D Inception-ResNet-v2 [START_REF] Lu | A practical Alzheimer disease classifier via brain imaging-based deep learning on[END_REF], Efficient 3D [START_REF] Yee | Construction of MRI-based Alzheimer's disease score based on efficient 3D convolutional neural network: Comprehensive validation on 7,902 images from a multi-center dataset[END_REF] and AD 2 A [START_REF] Guan | Attention-guided deep domain adaptation for brain dementia identification with multi-site neuroimaging data[END_REF]) and a classical voxel-based model using a SVM (Voxel-based SVM [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF]).

For AD diagnosis (i.e., AD/CN), all the methods show good balanced accuracy, although some of them failed to generalize on OASIS. In this scenario all the testing subjects. The found structures are known to be specifically and early impacted by AD. These structures are: bilateral hippocampus [START_REF] Frisoni | The clinical use of structural MRI in Alzheimer's disease[END_REF], left amygdala and left inferior lateral ventricle [START_REF] Coupé | Lifespan changes of the human brain in Alzheimer's disease[END_REF], left parahippocampal gyrus [START_REF] Kesslak | Quantification of magnetic resonance scans for hippocampal and parahippocampal atrophy in alzheimer's disease[END_REF], left posterior insula [START_REF] Foundas | Atrophy of the hippocampus, parietal cortex, and insula in alzheimer's disease: a volumetric magnetic resonance imaging study[END_REF], left thalamus proper [START_REF] De Jong | Strongly reduced volumes of putamen and thalamus in Alzheimer's disease: An MRI study[END_REF], left transverse temporal gyrus [START_REF] Liu | Education increases reserve against Alzheimer's disease-Evidence from structural MRI analysis[END_REF], left ventral diencephalon [START_REF] Lebedeva | MRI-based classification models in prediction of mild cognitive impairment and dementia in late-Life depression[END_REF]. While other attention-based deep methods failed to find structures related to AD [START_REF] Bron | Cross-cohort generalizability of deep and conventional machine learning for MRI-based diagnosis and prediction of Alzheimer's disease[END_REF], our DG framework shows high correlation with current physiopathological knowledge on AD [START_REF] Jack | Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade[END_REF].

Conclusion

In this paper, we addressed three major limitations of CNN-based methods by introducing a novel interpretable, generalizable and accurate deep grading framework. First, deep grading offers a meaningful visualization of the disease progression. Second, we proposed a collective artificial intelligence strategy to improve the generalization of our DG strategy. Experimental results showed a gain for both SVM and GCN in all tasks using this strategy. Finally, we proposed to use a graph-based modeling to better capture AD signature using both inter-subject similarity and intra-subject variability. Based on that, our DG method showed state-of-the-art performance in both AD diagnosis and prognosis.

Fig. 1 .

 1 Fig. 1. Overview of our processing pipeline. The MRI image, its segmentation and the deep grading map illustrated are from an AD subject.

Fig. 2 .

 2 Fig. 2. Average grading map per group of subjects.

Table 1 .

 1 Number of participants used in our study. Data used for training is in bold.

	Dataset	CN	AD	sMCI	pMCI
	ADNI1	170	170	129	171
	ADNI2	149	149	-	-
	AIBL	233	47	12	20
	OASIS3	658	97	-	-
	MIRIAD	23	46	-	-

Table 2 .

 2 Validation of the collective AI strategy, GCN classifier, the combination of DG feature with other image and non-image features using GCN classifier. Red: best result, Blue: second best result. The balanced accuracy (BACC) is used to assess the model performance. The results are the average accuracy of 10 repetitions and presented in percentage. All the methods are trained on the AD/CN subjects of the ADNI1 dataset.

	Classifier Features	AD/CN ADNI2 OASIS MIRIAD AIBL ADNI1 AIBL AD/CN p/sMCI pMCI/sMCI Average
	SVM	DG I	83	83	88	79	65	66	83.3	65.5
	SVM	DG C	83	84	91	87	68	64	86.3	66.0
	GCN	DG I	84	88	96	82	68	73	87.5	70.5
	GCN	DG C	87	89	100	88	70	76	91.0	73.0
	GCN DG C , V, A	87	88	98	92	74	74	91.3	74.0

Table 3 .

 3 Comparison of our method with current methods in AD diagnosis and prognosis. Red: best result, Blue: second best result. The balanced accuracy (BACC) is used to assess the model performance. All the methods are trained on the AD/CN subject of the ADNI1 dataset (except[START_REF] Lian | Hierarchical fully convolutional network for joint atrophy localization and Alzheimer's disease diagnosis using structural MRI[END_REF] that is fined-tuned on MCI subjects for sMCI/pMCI task).

Acknowledgments This work benefited from the support of the project Deepvol-Brain of the French National Research Agency (ANR-18-CE45-0013). This study was achieved within the context of the Laboratory of Excellence TRAIL ANR-10-LABX-57 for the BigDataBrain project. Moreover, we thank the Investments for the future Program IdEx Bordeaux (ANR-10-IDEX-03-02), the French Ministry of Education and Research, and the CNRS for DeepMultiBrain project.

Methods

AD/CN pMCI/sMCI ADNI2 OASIS MIRIAD AIBL ADNI1 AIBL Landmark-based CNN [START_REF] Liu | Landmark-based Deep Multi-Instance Learning for Brain Disease Diagnosis[END_REF] 91 -92 ---Hierachical FCN [START_REF] Lian | Hierarchical fully convolutional network for joint atrophy localization and Alzheimer's disease diagnosis using structural MRI[END_REF] 89 (unseen datasets), our method obtained high accuracy for all the datasets. This confirms the generalization capacity of our approach against domain shift. For AD prognosis (i.e., pMCI/sMCI), we observe a significant drop for all the methods. This drop is expected since pMCI/sMCI classification is more challenging and since models are trained on a different task (i.e., AD/CN). For this task, our method is generally robust, especially on AIBL. Moreover, our approach is the only deep learning method that performs competitively with the SVM model [START_REF] Wen | Convolutional neural networks for classification of Alzheimer's disease: Overview and reproducible evaluation[END_REF] on ADNI1, while significantly better on AIBL. In this scenario (unknown task), our method obtains the highest accuracy on average. These results highlight the potential performance of our method on unseen tasks.

Interpretation of collective deep grading. To highlight the interpretability capabilities offered by our DG feature, we first compute the average DG map for each group: AD, pMCI, sMCI and CN (see Figure 2). First, we can note that the average grading maps increase between each stage of the disease. Second, we estimated the top 10 structures with highest absolute value of grading score over


domain: phys spi

text:
Reza Rabani

Mohammad Hassan Saidi

Laurent Joly

Samy Merabia

Ali Rajabpour
email: rajabpour@eng.ikiu.ac.ir

Enhanced local viscosity around colloidal nanoparticles probed by Equilibrium Molecular Dynamics Simulations

Keywords: Green-Kubo, Nanofluids, Heterogeneous systems, Diffusion coefficient, nanolayer, Stokes-Einstein

Nanofluids -dispersions of nanometer-sized particles in a liquid medium-have been proposed for a wide variety of thermal management applications. It is known that a solid-like nanolayer of liquid of typical thickness 0.5-1 nm surrounding the colloidal nanoparticles can act as a thermal bridge between the nanoparticle and the bulk liquid. Yet, its effect on the nanofluid viscosity has not been elucidated so far. In this article, we compute the local viscosity of the nanolayer using equilibrium molecular dynamics based on the Green-Kubo formula. We first assess the validity of the method to predict the viscosity locally. We apply this methodology to the calculation of the local viscosity in the immediate vicinity of a metallic nanoparticle for a wide range of solid-liquid interaction strength, where a nanolayer of thickness 1 nm is observed as a result of the interaction with the nanoparticle. The viscosity of the nanolayer, which is found to be higher than its corresponding bulk value, is directly dependent on the solid-liquid interaction strength. We discuss the origin of this viscosity enhancement and show that the liquid density increment alone cannot explain the values of the viscosity observed. Rather, we suggest that the solid-like structure of the distribution of the liquid atoms in the vicinity of the nanoparticle contributes to the nanolayer viscosity enhancement. Finally, we observe a failure of the Stokes-Einstein relation between viscosity and diffusion close to the wall, depending on the liquid-solid interaction strength, which we rationalize in terms of hydrodynamic slip.

I. INTRODUCTION

Nanofluids are colloidal solutions in which suspended nanometer-sized particles in the liquid medium [1][2][START_REF] Subramanian | Nanofluids Their Eng. Appl[END_REF] are responsible for enhanced thermophysical properties such as density, thermal conductivity, diffusivity, and viscosity compared to the base fluids [START_REF] Taylor | [END_REF][5][6][7][8][9] . The nanoparticles used in nanofluids are typically made of metals [10][11][12] , oxides [13][14][15][16] , carbides [17][START_REF] Ezekwem | Energy Sources, Part A Recover[END_REF][19] , or carbon nanotubes [20][21][22][23] , which are immersed in base fluids such as water [24][25][26] , ethylene glycol [27][28][29] , oil [30][31][32][33] , and refrigerant [34][35][START_REF] Peng | Int. J. Refrig[END_REF] . The main interest behind introducing the nanofluid concept is the increment in the thermal conductivity compared to the conventional fluids [START_REF] Satti | [END_REF][38][39][40][START_REF] Han | Energy Procedia[END_REF][START_REF] Maradiya | [END_REF] . The thermal conductivity of the nanofluids depends on various factors such as Brownian motion [43][44][START_REF] Yang | J. Heat Transfer[END_REF] , thermophoresis [START_REF] Michaelides | [END_REF][47][48] , and clustering [49][50][51] of the solid particles as well as phonon-based nature of heat transfer within solid nanoparticles 52,53 . Besides, liquid molecules form a local structured layer around the nanoparticle, known as nanolayer, which acts as a thermal bridge between the nanoparticle surface and the bulk liquid 54 that can enhance the thermal conductivity of the nanofluids by 10~20% depending on the base fluid type and nanoparticle diameter 55,56 .

The term nanolayer describes the solid-like interfacial layer between the liquid and the nanoparticle, in which liquid molecules form an ordered layer structure similar to shell on the nanoparticle solid surfaces as a result of intermolecular forces existing between liquid and solid atoms 57 . The thickness of the nanolayer vary between 0.3 𝑛𝑚 to 1 𝑛𝑚 (corresponding to 1 to 3 molecular sizes) depending on the interface interaction strength [58][59][60][61][62] . Numerous studies have attempted to model the thermophysical properties of nanofluid by considering the role of this nanolayer. As an example, the role of nanolayer in the enhanced thermal conductivity of nanofluids was investigated thoroughly in a renovated Maxwell model 62 for spherical particles suspensions and a renovated Hamilton-Crosser model 63 for non-spherical particles suspensions.

The effect of nanolayer on the density of nanofluid has been also investigated experimentally and it was shown that in order to reproduce the experimental densities of the nanofluids, the nanolayer should be taken into the account 64 . Molecular dynamics (MD) simulations of silverwater nanofluids were also conducted in order to find the effect of nanolayer on the density of the nanofluid and a new ternary mixture model was proposed for the computation of nanofluids density. It was shown that the nanolayer density is about 1.6 times that of the bulk fluid. Such high density causes a contraction in the fluid which reduces the volume of the base fluid by 3% 60 .

Considering the viscosity of nanofluids, a few models have been proposed during the past decades 65 . For a dilute suspension of small, rigid, spherical particles in the base fluid, a good approximation of the nanofluid viscosity can be obtained by Einstein model 66 . For a moderate particle concentration, Brinkman has modified Einstein's equation and proposed a new model 67 .

Batchelor 68 proposed a new model to account for the effect of Brownian motion on the viscosity of an isotropic suspension of spherical, rigid particles. There are only a few models that include the nanolayer effect in the viscosity of the nanofluid. Considering an elliptical liquid layer shell around the nanoparticles, the Ward model for calculating the viscosity was proposed 69 . The above-mentioned models for the viscosity of the nanofluid were compared with MD simulations 60 . It was shown that the calculated viscosity based on the models which do not include the nanolayer effect is not consistent with MD while the models which account for a nanolayer shell around nanoparticles follow the MD simulations viscosity well. The key role of the nanolayer viscosity highlighted in this work makes it critical to characterize it for a variety of liquid-solid interfaces.

Liquid Viscosity and diffusivity are known to be affected by the fluid heterogeneity at the molecular scale such as density inhomogeneity 70 . The effects of density heterogeneities at the solid/liquid interface on the shear viscosity profile were investigated thoroughly for a sheared confined fluid by the use of non-equilibrium molecular dynamics (NEMD) simulations 71 . It was shown that the local viscosity varied strongly with the distance from the walls for both the dense and dilute fluid, displaying oscillations correlated to the density profile. This work was further developed and a local average density model was developed to describe the distribution of the local shear viscosity in inhomogeneous fluid 72 . It was shown that for a dense inhomogeneous liquid, the proposed model leads to a good description of the viscosity profiles as measured by NEMD. However, some deviations were observed for the low-density inhomogeneous liquid.

Such viscosity variation was shown to be able to significantly reduce the flow rate for a Poiseuille flow as the system dimension became smaller than 3~4 nm for water and 1 nm for nonpolar fluids 73 . The effects of the inhomogeneity length scale on structural and dynamic properties of the fluid were investigated thoroughly by applying sinusoidal density profiles having various wavelengths 74 . It was shown that for a long-wavelength density distribution, bulk-like relationships between thermodynamics, local structure, and diffusivity were observed. However, as it came to small wavelengths, different correlations should be considered between the local static and dynamic quantities, depending on the wavelength.

Due to the importance of the precise viscosity calculation in the nanofluid, the objective of the present research is to characterize the viscosity of the nanolayer around the nanoparticles.

To that aim, we will use equilibrium molecular dynamics (EMD) simulations, where the viscosity of the system is computed with the Green-Kubo (GK) relation. The ability of the GK relation for the calculation of bulk and confined liquid viscosity were investigated by several researchers 75- 83 . The viscosity of confined water inside hydrophilic and hydrophobic nanotubes was computed using the GK method 84 . A great impact of water-wall interactions, confinement size, and density on the viscosity was observed. However, it was shown that for the small nanotube, as the density is increased, the viscosity and the diffusion coefficient are increased together and violates the Stokes-Einstein relation. Zaragoza et al. 85 studied how the planar walls and nanotube confinements affect water viscosity and proposed a confined Stokes-Einstein relation for obtaining the viscosity from diffusivity. It is interesting to note that, while the viscosity computed with the GK method differs from that of the bulk noticeably, the viscosity computed with confined Stokes-Einstein relation is not affected by the confinement. The authors concluded that different methods in calculating the viscosity may provide dramatically different results, which may not easily be related to the standard and experimental definition of viscosity. However, a NEMD simulation of shear-driven liquid argon flow between two graphene walls shows an increment by the factor two for the viscosity of liquid argon in the vicinity of the nanochannel wall compared to the middle of the channel 86 . Recently, Zhou et al. 87 have shown that in the viscosity calculation of the nanoconfined fluid, the wall friction should be decoupled from the fluid viscosity by defining the frictional (near the wall), the transitional, and the viscous (far from the wall) region. This was concluded by a comparison between the fluid viscosities calculated from the GK formula in the entire domain and the viscous region, and the fluid viscosity which was obtained based on the velocity profile of the Hagen-Poiseuille flow. This analogy revealed that only the viscosity in the viscous region coincides with the one deduced from the velocity profile.

The nanolayer definition in our study coincides with the definition of the frictional plus transitional region, whose viscosity was not discussed in previous studies. A closer look at the nanolayer viscosity is the main motivation of our study. Therefore, we should be able to apply the GK method to the separate liquid layers and as the first step of our research, this feasibility is validated. In addition, the diffusion coefficient of the liquid layers is also calculated, which is then connected to the difference between the viscosity of the different layers. The remainder of this article is organized as follows: In Section II, the methodology of the simulation is described in detail. In Section III, the simulation results are presented. In the first part of this section, a validation is presented to assess the accuracy of the applied simulation method. Then the viscosity of the nanolayer of the copper-argon nanofluid is calculated and the results are thoroughly interpreted. Finally, we conclude in Section IV. It is worth mentioning that the rationale behind choosing copper-liquid argon is that this is a model system that can be used as the initial step in the study of new phenomena. Argon, in practice modeled with a Lennard-Jones potential, may be thought of as the simplest liquid one can model, and whose properties can be mapped to any real liquid using the law of corresponding states. Nevertheless, it would be interesting to extend this work by considering models of more complex liquids such as ethylene glycol or water.

II. COMPUTATIONAL DETAILS AND METHODS

In this study, the viscosity of the liquid argon is calculated by EMD method using LAMMPS 88 (Large-Scale Atomic/Molecular Massively Parallel Simulator). The Open Visualization Tool (OVITO) is used to visualize and represent atomic configurations 89 . Lennard-Jones (LJ) potential parametrization is considered for interatomic forces between different types of atoms in this study. The Lennard-Jones (LJ) potential parameters that were used in this study are shown in is used in all simulations, smaller timesteps were also tested and had no significant effect on the results. Considering the temperature of the liquid medium and the walls, a Maxwell-Boltzmann velocity distribution is used to initialize the velocity of the atoms at the beginning of the simulation. In all sets of simulations, two independent Nosé-Hoover thermostats were applied to liquid and wall for 2 𝑛𝑠 to ensure the desired temperature is achieved. After that, the thermostat is removed from the liquid domain while still applied to the walls. Then, the simulation proceeds for 2 𝑛𝑠. This ensured us that no thermostat is applied to the fluid atoms directly and the results are not affected by the selection of thermostat parameters. Additionally, another 2 𝑛𝑠 are performed for the averaging process in which the density, temperature, and viscosity of each computational bin are calculated.

The shear viscosity of a fluid, 𝜇, can be computed using the GK relation which is based on the shear stress 𝜏 𝛼𝛽 autocorrelation function 90,[START_REF] Allen | Computer Simulation of Liquids: Second Edition[END_REF] :

𝜇 = 𝑉 𝑘 𝐵 𝑇 ∫ 〈𝜏 𝛼𝛽 (𝑡)𝜏 𝛼𝛽 (0)〉𝑑𝑡 +∞ 0 , (1)
with 𝛼 ≠ 𝛽,

𝜏 𝛼𝛽 = 1 𝑉 (∑ 𝑃 𝛼 𝑖 𝑃 𝛽 𝑖 𝑚 𝑖 𝑁 𝑖=1 + ∑ ∑ 𝑟 𝑖𝑗 𝛼 𝐹 𝑖𝑗 𝛽 𝑁 𝑗>𝑖 𝑁 𝑖=1 ), (2)
where 𝑉 and 𝑇 are the volume and temperature of the system, respectively, 𝑘 𝐵 the Boltzmann constant, 𝑚 the mass, 𝑃 the momentum, as well as 𝑟 and 𝐹, are the distance and the force between two particles, respectively; 𝑖 and 𝑗 are the atomic indexes and 𝑁 is the number of particles. According to the GK relation, the calculated shear viscosity of equation ( 1) consists of three components, i.e., 𝜇 𝑥𝑦 , 𝜇 𝑥𝑧 and 𝜇 𝑦𝑧 and the shear viscosity of the fluid equals their average.

In order to interpret the result, the diffusion coefficient was also determined using the mean square displacement of the particle 85 :

𝐷 = lim 𝑡→∞ 〈|𝑟(𝑡 0 + 𝑡) -𝑟(𝑡 0 )| 2 〉 2 𝑑𝑖𝑚 × 𝑡 (3)
where r denotes the position vector. The standard formula considering the three-dimensional displacement uses dim=3 while for two-dimensional displacement, dim=2 should be used.

However, since in the current study the MSD is calculated separately in each direction and the diffusion coefficients in perpendicular and parallel directions are computed based on them, dim=1 is considered in equation (3). In order to have a reasonable estimate of the intrinsic diffusion coefficient, the liquid center of mass motion was subtracted before computing the MSD 92 . The validity of the Stokes-Einstein relation was investigated using [93][94][95] :

𝜇 = 𝑘 𝐵 𝑇 3𝜋𝐷𝜎 ℎ (4)
which relates to the viscosity and the diffusivity of the fluid. In this equation, 𝜎 ℎ denotes the hydrodynamic fluid particle diameter. The above-mentioned equations for calculating the viscosity and diffusion coefficient are for the calculation of the entire volume of the liquid.

However, as it comes to layers of liquids, these relations must be used with caution. In particular, since the calculation of the diffusion coefficient is based on the MSD according to equation ( 3), some restrictions must be considered. Actually, the displacement of an atom is calculated based on its reference position which is normally the original position at the onset of the MSD calculations 88 . Consequently, as an atom moves to the adjacent liquid layer, the validity of the MSD calculations become under question. In order to handle such restriction in our study, the MSD calculation time is chosen smaller than the required time that an atom needs to move to its adjacent liquid layer. Considering the diffusion coefficient of bulk liquid argon 96 and liquid layer thickness, the required time for an atom to move to its adjacent liquid layer is on the order of 0.3 𝑛s and in order to be cautious in this study, 0.2 𝑛s is considered as the maximum MSD calculation time.

III. RESULTS

It should be noted that the GK method is derived for a bulk system 104,105 and it was verified that for a system with geometrical confinement, in the direction without confinement where the periodic boundary conditions are applied, calculated thermal conductivity based on the GK equation yields results which agree well with the Non-Equilibrium Molecular Dynamics (NEMD) 106 . Therefore, the first step in this paper is to validate the use of GK viscosity calculations for a system with a non-periodic boundary condition.

A. Validation of the GK methodology for the liquid layer

The simulation domain for the validation is presented in FIG. 1. The system is a cube with sides of 6 𝑛𝑚 filled with 4280 argon atoms at the temperature of 100 𝐾, which corresponds to a density of 1314 𝑘𝑔 𝑚 3 ⁄ . Considering the values of the density and temperature, the viscosity is expected to be 184 µ𝑃𝑎. 𝑠 107 . It is well-known that in EMD calculations, larger system sizes do not significantly affect the calculated shear viscosity [108][109][110][111] . Therefore, we proceed with the above-mentioned dimensions. Periodic boundary conditions are applied along x, y, and z directions. Regarding this configuration, two sets of simulations were conducted. In the first set of simulations, the viscosity of the whole simulation box is considered as a unit configuration, and the components of the viscosity for the whole simulation box are calculated using equations

(1) and (2). In the second set of simulations performed simultaneously, the simulation box is divided into two equal simulation domains considering a virtual plane perpendicular to the x direction. The atoms can freely move between zones from this plane. Since for each instance, the autocorrelation of the shear stress components is used for the calculation of the viscosity, it does not matter if the atoms cross the virtual plane as the calculation is done based on the neighboring atoms. Using such a methodology, the components of the viscosity in each region are calculated separately. The integration time for the viscosity calculation by GK method in this part is 10 𝑝𝑠 according to Fig. 2. The calculated viscosity is averaged over twenty-five different initial independent configurations and the reported viscosity deviation is defined based on them. It should be noted that the error bar of the reported viscosity components is the standard deviation of the tail of the final running 87 .

It is observed that as the GK method is applied on the whole simulation box, the computed viscosity,180 ± 5 µ𝑃𝑎. 𝑠, is close to the expected value while as the GK method is applied on the left or the right box, the calculated viscosity becomes 130 ± 5. A closer look at the boundary conditions might reveal the source of such behavior.

Actually, the periodic boundary condition is applied in 𝑦 and 𝑧 directions, consequently in 𝑦𝑧 plane, which leads to the correct viscosity prediction. On the other hand, the GK method underestimates the 𝜇 𝑥𝑦 and 𝜇 𝑥𝑧 which should be related to the presence of the non-periodic boundary condition, the imaginary boundary, in the middle of the simulation box along 𝑥 direction. It can be concluded that similar to the thermal conductivity 106 , the viscosity calculation using the GK method can be only applied in a direction having periodic boundary conditions. It is also interesting to notice that the boundary condition in the middle of the box is not a wall or any physical confinement. It is just an imaginary boundary that divides the simulation box into two parts. Besides, it is worth mentioning that the convergence of the autocorrelation function of all the three viscosity components does not guarantee the correctness of calculated viscosity in the GK method, and further criteria should also be checked.

B. Copper nanoparticle in liquid argon

In the previous section, we observed that even in the presence of a non-periodic boundary condition, the viscosity component in a plane with periodic boundary conditions in all its directions still gives the local viscosity of the fluid in that region. Considering a nanoparticle and spherical nanolayer around it to calculate the viscosity, such a plane cannot be defined. For the nanoparticles with a diameter larger than 20 𝑛𝑚, a small part of the nanoparticle surface can be modeled as a flat wall approximately as can be seen in FIG. 3. Therefore, the simulation domain in this part is considered as a flat wall consisting of two layers of FCC copper atoms, which is extended 8 𝑛𝑚 in y and z directions, and a 5 𝑛𝑚 liquid domain in the x direction. The temperature of the wall and argon is kept at 100 𝐾. The liquid is divided into 5 bins of 1 nm thickness perpendicular to the wall for the averaging process in the MD method. Therefore, it should be noted that the center of the first, second, third, fourth, and fifth bins is located at 5, 15, 25, 35, and 45 nm from the wall respectively. It should be mentioned that all the presented data for the temperature, density, and viscosity is the value averaged over 20 different initial distributions of the liquid atoms. To find the effect of the interaction strength between nanoparticle surface and the liquid on the density and viscosity of the nanolayer, 𝜀 𝐶𝑢-𝐴𝑟 is considered as 0.25, 0.5, 1.0 and, 1.5 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙. Actually, 1 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙 approximately equals the case with Lorentz-Berthelot mixing rules for interaction strength; √𝜀 𝐶𝑢-𝐶𝑢 × 𝜀 𝐴𝑟-𝐴𝑟 = 1.05 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙. Changing the solid-liquid interaction strength varies the number of liquid atoms in the nanolayer, the first layer adjacent to the wall, which consequently affects the liquid density in the other layers. To be able to compare these solid-liquid interaction strengths together, the number of liquid atoms in the simulation domain has been changed, based on a trial-and-error method, in such a way that the density of the liquid remains approximately 1314 𝑘𝑔 𝑚 3 ⁄ for the third to the fifth layers (which can be considered as the bulk liquid). The calculated number of liquid atoms for each case is shown in Table 4. The integration time for the viscosity calculation based on the autocorrelation of the shear stress components in this section is 50 𝑝𝑠. More details about the viscosity calculation can be found in the supporting information (SI). Meanwhile, the density of the other layers is about 1314 𝑘𝑔 𝑚 3 ⁄ as expected. The increase in the density of the liquid in the first and second layers corresponds to the well-known density layering phenomenon near solid surfaces, 113 which is a direct consequence of interactions with the wall in the near-wall region 97 .

FIG. 4. Local liquid density in the different layers.

As mentioned before, the spatial viscosity distribution in the liquid is calculated using the GK relation as described in the previous section. According to FIG. 5, it is interesting to notice that the viscosity of the first layer is approximately 1.5, 2, 6, and 12 times the bulk viscosity as the interaction strength increased, which shows the great impact of the solid surface on its adjacent liquid viscosity. Interpreting the main reasons behind such anomalous viscosity increment of the first layer is the main subject of this manuscript.

One reason behind such viscosity increment is related to the denser liquid medium of the first layer compared to the other layers. A denser liquid medium increases the collision between liquid atoms which enhances the viscosity accordingly. This can be also clearly observed in the viscosity enhancement of the second layer compared to the third, fourth, and fifth layers as it has higher liquid density according to Fig. 4. It should be noted that the viscosities of the second to the fifth layer are still in the same order of magnitude which implies that a similar mechanism governs the transport phenomenon in these layers. Meanwhile, the viscosity of the first layer is much higher compared to the other layers, which suggests that a different mechanism might govern the transport phenomenon in this layer.

It might be helpful to find the viscosity of the liquid at the corresponding bulk density for the first layer. This information will show us whether the viscosity increment of the nanolayer is only related to the increase in the liquid density or if some other phenomenon is also involved. The viscosities based on the liquid density for the first layer are presented in Table 4. Lemmon and Jacobsen 107 presented these viscosities for the bulk liquid medium based on a combination of the theoretical models and the empirical equations. For density higher than 1481.5 𝑘𝑔/𝑚 [START_REF] Subramanian | Nanofluids Their Eng. Appl[END_REF] (melting density at 𝑇 = 100𝐾), the viscosity is taken from Galliero et al. 112 in which an accurate correlation is presented for the viscosity of the Lennard-Jones fluid. The ratio of this viscosity to the calculated one can be considered as the contribution of the density to the total viscosity. It is observed that the viscosity of the first liquid layer based on its density is 227, 330, 770, and 1227 µ𝑃𝑎. 𝑠 as the solid-liquid interaction strength is increased from 0.25 to 1.5 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙, which shows a notable difference compared to the viscosity of the nanolayer. Therefore, it is concluded that in addition to the liquid density, there must be some other phenomenon that increase the viscosity of the nanolayer. In the liquid nanolayer adjacent to a nanoparticle, the liquid atoms are under the strong attraction of the nanoparticle atoms on the surface. On the other hand, the reported values in Table 4 are for the bulk condition where the liquid atoms can move and interact with each other freely. Therefore, a closer look at the liquid structure and transport properties of the atoms in the nanolayer might reveal the main reason behind such discrepancy between the calculated viscosity of the nanolayer and the expected viscosity based on its density.

Figure 6 displays the distribution of the liquid atoms near the nanoparticle surface. As it can be observed, the liquid atoms in the nanolayer form two ordered liquid atoms layers for the interaction strength of 0.25 kcal/mol while three layers are formed for higher interaction strength. Note that the contact layer always remains liquid: as detailed below, atoms are still diffusing, both within and between the layers. We suggest that the anomalous viscosity enhancement of the nanolayer compared to the other layers can be assigned to such liquid atoms ordering, and we will refer to it as the contribution of the liquid structure to viscosity. This can be computed by subtracting the contribution of the liquid density from total viscosity. Figure 7 compares the contribution of the liquid density and the layered structure in the total viscosity of the nanolayer which shows a greater role played by the liquid structure as the liquid-solid interaction strength increased.

The definition of the nanolayer in our study coincides with the sum of the frictional (distance lower than 1.5σ from the wall) and transitional (distance between 1.5σ and 𝑟 𝐶 𝐿𝐽 from the wall) regions defined by Zhou et al. 87 . It was observed that in the frictional and transitional regions, the liquid atoms are attracted by the wall atoms, inducing additional friction on the liquid atoms compared to the bulk. This was shown to be intensified as the solid-liquid interaction energy increased. The same phenomenon is observed in our study considering Fig. 7 in which the contribution of the liquid structure to the total viscosity is enhanced as the solid-liquid interaction energy increased. It seems that besides the induced wall friction on the liquid atoms in the frictional region which is enhanced as the solid-liquid interaction energy increased, the formation of the ordered liquid structure in the transitional region is another source for the friction increment. This might be related to an increment in the friction between the neighboring liquid atoms as the liquid structure becomes ordered. To go further, the mean square displacement (MSD) and the diffusion coefficient can be considered as the main criteria of the transport characteristics of the liquid atoms. Therefore, a more detailed evaluation of these parameters might help us to have a better understanding of the underlying physics in the nanolayer. The diffusion coefficients are calculated based on the MSD of the liquid atoms using equation (3). According to Fig. 3, the diffusion coefficient perpendicular to the nanoparticle surface is 𝐷 𝑥 and the average of 𝐷 𝑦 and 𝐷 𝑧 is considered as the diffusion coefficient parallel to the surface. Figure 8 shows the variation of the diffusion coefficient parallel and perpendicular to the nanoparticle surface as a function of solid-liquid interaction strength for liquid layers. These values are also shown in Table 6 andTable 7. It is observed in Fig. 8a that as the liquid-solid interaction energy increases from 0.25 to 1.5 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙, the diffusion coefficient parallel to the nanoparticle surface in the third to the fifth layers, is about 3.35 × 10 -9 𝑚 2 𝑠 -1 while a reduction from 3.02 × 10 -9 to 2.76 × 10 -9 𝑚 2 𝑠 -1 and 1.86 × 10 -9 to 0.65 × 10 -9 𝑚 2 𝑠 -1 is observed for the second and first layers respectively. Figure 8b clearly shows that for the same increment of interaction energy, the diffusion coefficient perpendicular to the nanoparticle surface is approximately 3.30 × 10 -9 𝑚 2 𝑠 -1 for the fourth and the fifth layers, while it decreases from 2.92 × 10 -9 to 2.72 × 10 -9 𝑚 2 𝑠 -1 for the third layer. For the second and first layers, the diffusion coefficient decreases more, from 2.16 × 10 -9 to 1.83 × 10 -9 𝑚 2 𝑠 -1 and 1.09 × 10 -9 to 0.5 × 10 -9 𝑚 2 𝑠 -1 respectively. The diffusion coefficient perpendicular direction to the nanoparticle surface is inversely related to the residence time of the liquid atoms in the layers. Therefore, it can be concluded that the liquid atoms remain more time in the first liquid layer compared to the other layers for all interaction strength according to Using equation ( 4 ] × 10 -9 𝑚 2 𝑠 -1 respectively as the solid-liquid interaction strength increases. The effect of liquid-solid slip on the diffusion coefficient close to the wall of the nanoconfined liquid medium was investigated thoroughly 121,122 . Saugey et al. 123 studied the effect of slip length at the solid-liquid interface on the diffusion of the confined liquid numerically, using a hydrodynamic approach based on the Stokes equation. It was shown that for the no-slip boundary condition, the diffusion coefficient in the vicinity of the wall decreased up to 50% compared with the bulk diffusion, while as the slip on the boundary increased, the diffusion coefficient was enhanced by the same amount. Interestingly, the authors have also derived an analytical expression for low and high confinement limits, which were in good agreement with numerical data. Using this analytical approach, we have proposed a similar analytical formula for our simulation domain as:

𝐷 || (𝑥) = 𝐷 𝑏𝑢𝑙𝑘 {1 -( 9 16
) 𝜎 𝑥 + 𝑏 -𝜎 }

where 𝜎 is the diameter of the liquid atoms which is 3.405 Angstrom in our case, 𝑏 is the slip length, and 𝐷 𝑏𝑢𝑙𝑘 is the diffusion coefficient in the bulk region away from the solid surface. It should also mention that according to the simulation domain in Fig. 3, we have also considered the full formula accounting for periodic boundary conditions, and this does not lead to a better agreement. The diffusion coefficients derived based on the MD simulation are shown in Fig. 11a along with the fitted analytical formula for different interaction energy. The bulk diffusion coefficient 𝐷 𝑏𝑢𝑙𝑘 and the slip length 𝑏 in equation ( 5) is considered as the fitting parameters and their corresponding values are obtained for each interaction energy and shown in Fig. 11b. It can be observed that 𝐷 𝑏𝑢𝑙𝑘 is around 3.5 × 10 -9 𝑚 2 𝑠 -1 , independently of the interaction energy.

Meanwhile, Fig. 11b shows that the slip length declines from 2.51 to 0.75 Å as the interaction energy increases which can be attributed to the increment of the induced friction on the frictional region in liquid by the wall 87 .

The difference between 𝐷 𝑏𝑢𝑙𝑘 based on the analytical formula, equation (5), and the calculated 𝐷 || for the fourth and fifth layers based on our numerical simulation presented in Figs.

8a can be considered as a possible manifestation of finite size effect for each layers [START_REF] Simonnin | [END_REF]123 . It should be although noted, however, that for the fourth and fifth layers, the ratio of the distance from the nanoparticle surface to the radius of the argon atom is in the order 25. Therefore, it is expected that no confinement effect exists for the single wall configuration 123 . Considering 3.5 × 10 -9 and 3.3 × 10 -9 𝑚 2 𝑠 -1 as the diffusion coefficients of the 𝐷 𝑏𝑢𝑙𝑘 and 𝐷 || respectively, the finite size effect of each layer should be about 0.2 × 10 -9 𝑚 2 𝑠 -1 . It is difficult to provide a better estimate of possible finite size effects, as analytical formulas in the literature only describe 3D periodic systems 109 or liquids confined by no-slip walls [START_REF] Simonnin | [END_REF] .

IV. SUMMARY AND CONCLUSION

The EMD simulation of liquid argon around a copper nanoparticle reveals that the density increments can be observed at 2 nm from the nanoparticle surface, while it is much more intense in the first 1 nm known as the nanolayer. In this study, we characterize the corresponding viscosity increase using the Green-Kubo approach. It is known that the GK method is traditionally derived for the bulk medium in which all the boundary conditions are periodic. Considering the liquid medium adjacent to the nanoparticle surface, we have shown that the GK method can be applied for such configuration with non-periodic boundary conditions provided that only the viscosity component in the normal plane to the non-periodic direction is considered in the calculation. Then, the viscosity of each layer is calculated and it was shown that for the liquidsolid interaction strength of 0.25, 0.5, 1.0 and, 1.5 kcal/mol, the viscosity is multiplied by approximately 1.5, 2, 6, and 12. We have shown that a fraction of this anomalous viscosity enhancement originates from the increment in the liquid density of the nanolayer. The further increment is suggested to be related to the ordered structure of the distribution of the liquid atoms in the nanolayer. We have shown that the contribution of the liquid structure to the total viscosity of the nanolayer is increased as the solid-liquid interaction strength is enhanced.

Generally, it was inferred that for the low and intermediate interaction energy systems, the enhanced viscosity of the nanolayer can be assigned to the small increase of the local density.

For these systems, the Stokes-Einstein relation -linking viscosity and diffusion-holds approximately in the nanolayer and the effect of the wall friction is small. On the other hand, for the strong interaction energy systems, the increase of the density cannot explain the enhancement of the local viscosity. This enhancement should be rather controlled by the friction with the wall and other neighboring liquid atoms. In this regime, one expects deviations to the Stokes-Einstein relation as it was observed. Finally, we rationalized the deviations from the Stokes-Einstein relation by accounting for the effect of liquid-solid slip-on diffusion close to the wall, with a slip length that decreases with increasing liquid-wall interaction.

FIG. 1 .FIG. 2 .

 12 FIG. 1. (a) Simulation domain filled by liquid argon atoms and (b) schematic sketch of the simulation domain which is divided into two parts.

FIG. 3 .

 3 FIG. 3. Simulation domain considered for viscosity calculation, which is a part of 20 𝑛𝑚 diameter copper nanoparticlesurface that is immersed in liquid argon.

FIG. 5 .

 5 FIG. 5. Normalized viscosity distribution in the liquid argon normalized with 𝜇 𝑏𝑢𝑙𝑘 = 184 µ𝑃𝑎. 𝑠.

FIG. 6 .FIG. 7 .

 67 FIG. 6. Distribution of the liquid atoms near the nanoparticle surface in the steady state for liquid-solid interaction strength of 0.25 (a), 0.5 (b), 1.0 (c), and 1.5 kcal/mol(d).

FIG. 8 .

 8 FIG. 8. Variation of the diffusion coefficient parallel (a) and perpendicular (b) to the nanoparticle surface as a function of solid-liquid interaction strength for liquid layers

  ), the validity of the Stokes-Einstein relation for each layer is also investigated in Fig.9. Considering the nanolayer, it is interesting to notice that the Stokes-Einstein relation is valid for the low and intermediate solid-liquid interaction strength as 𝐷 || × 𝜇 is approximately constant for all liquid layers in Fig.9a. However, it should also be noted that a small reduction in the first layer is observed. The Stokes-Einstein relation implies a hydrodynamic radius which should depend on the local density, but if we take into account this correction as 𝐷 || × 𝜇 × 𝜌 -1/3 , we still have a small deviation for the first layer according to Fig.9b. On the other hand, Figs. 9a and 9b clearly show that for stronger interaction, the violation of the Stokes-Einstein relation becomes noticeable. This violation is mainly related to the different structural and transport properties of the first layer compared to the others as the solid-liquid interaction strength increased which indicates the decoupling of viscosity and diffusion.It is important to notice that the breakdown of the Stokes-Einstein relation has been also observed for a wide range of glassy systems 114 . While Stokes-Einstein breakdown was shown to be related to the fluctuation dominance in the dynamics of the glass formers at lowtemperature 115 , it was also shown that the atomic dynamics is dramatically slowed down[START_REF] Hansen | Theory of Simple Liquids: With Applications to Soft Matter: Fourth Edition[END_REF] for the glassy system. The dynamics slowing down is usually accompanied by the emergence of heterogeneous dynamics, which can manifest in the existence of solid-like clusters characterized by dynamics slower than the average 117 and whose size depends on the temperature118 . such a heterogeneous environment, the Stokes-Einstein relation may be broken as a result of the different spatial averages underlying liquid viscosity and diffusivity 117,119 . The changes in local liquid atomic structure and dynamic of liquid atom motion of the first liquid atom layer in our study are somehow similar to what was observed in the glassy system while their origin is different. In fact for the glassy system, it was shown that as the temperature lowered close to the glass transition, the diffusion coefficient 𝐷 || and 𝑇 𝜇 ⁄ have different functional dependence on the temperature 120 . It should be noted that in our current study violation of the Stokes-Einstein relation is only observed in the first layer and for the strong interaction energy between solid-liquid, 1.0 and 1.5 kcal/mol. Considering Fig. 8a, it is obvious that the diffusion coefficient for 1.0 and 1.5 kcal/mol are even lower than the corresponding values for 0.25 and 0.5 kcal/mol. Therefore, the violation of the Stokes-Einstein relation should be attributed to the anomalous viscosity increment for interaction energy of 1.0 and 1.5 kcal/mol as is shown in Fig. 5. The variation of the viscosity and the diffusivity of the nanolayer against the interaction strength are shown in Fig. 10 in which the solid lines show the exponential behavior. It is observed that the viscosity and the diffusion coefficient of the nanolayer parallel to the nanoparticle surface vary as [

FIG. 9 .FIG. 10 .

 910 FIG. 9. The evolution of the Stokes-Einstein 𝐷 || × 𝜇 (a) and modified Stokes-Einstein with a hydrodynamic radius 𝐷 || × 𝜇 × 𝜌 -1/3 (b) relation as a function of distance from the nanoparticle surface.

FIG. 11 .

 11 FIG. 11. Variation of the diffusion coefficient in different layers and its corresponding fitted analytical formula which is shown by the solid lines (a) and the variation of calculated slip length and bulk diffusion (b) as a function of solid-liquid interaction strength based on the proposed analytical formula

Table 1 ,

 1 in which 𝑟 𝐶 𝐿𝐽 denotes the cut-off distance of the LJ interactions. While a timestep of 1 𝑓𝑠

TABLE 1 .

 1 LJ parameters for non-bonding interactions

	𝑖𝑛𝑡𝑒𝑟𝑎𝑐𝑡𝑖𝑜𝑛	𝜀 (𝑘𝑐𝑎𝑙 𝑚𝑜𝑙 ⁄	)	𝜎 (Å)	𝑟 𝐶 𝐿𝐽 (Å)
	𝐴𝑟 -𝐴𝑟 97-100	0.237		3.405	27
	𝐶𝑢 -𝐶𝑢 101	4.72		2.616	10
	𝐴𝑟 -𝐶𝑢 102,103 0.25, 0.5, 1.0, 1.5	3.010	18

  3 µ𝑃𝑎. 𝑠 and 127 ± 6 µ𝑃𝑎. 𝑠 respectively where the deviation from the desired value is considerable. To find out the source of such deviations, the autocorrelation of the shear stress components, averaged over all cases, is shown in FIG.2. This figure clearly implies that the viscosity component in the 𝑦𝑧 plane should be higher than that of the 𝑥𝑦 and 𝑥𝑧 planes. The viscosity components for the left and right boxes are computed based on these autocorrelations. The calculated viscosity components for 𝑥𝑦,𝑥𝑧, and 𝑦𝑧 planes are 111 ± 5.3 µ𝑃𝑎. 𝑠, 112 ± 4.9 µ𝑃𝑎. 𝑠, and 178 ± 0.9 µ𝑃𝑎. 𝑠 respectively for the left simulation box and 108 ± 6 µ𝑃𝑎. 𝑠, 108 ± 5.7 µ𝑃𝑎. 𝑠, and 175 ± 1.3 µ𝑃𝑎. 𝑠 for the right ones. These values show that while the calculated viscosities in the 𝑥𝑦 and 𝑥𝑧 planes deviate significantly from 184 µ𝑃𝑎. 𝑠, the 𝑦𝑧 component is in good agreement with the expected value.

TABLE 4 .

 4 Calculated density of the nanolayer and its corresponding viscosity extracted from Lemmon and Jacobsen 107 andGalliero et al. 112

	𝜀 𝐶𝑢-𝐴𝑟		# 𝑙𝑖𝑞𝑢𝑖𝑑	𝑑𝑒𝑛𝑠𝑖𝑡𝑦	𝑣𝑖𝑠𝑐𝑜𝑠𝑖𝑡𝑦
	(𝑘𝑐𝑎𝑙 𝑚𝑜𝑙 ⁄	)	𝑎𝑡𝑜𝑚𝑠	(𝑘𝑔 𝑚 3 ⁄ )	(µ𝑃𝑎. 𝑠)
	0.25		6725	1367	227 107
	0.5		6850	1470	330 107
	1.0		6950	1540	770 112
	1.5		7000	1597	1227 112

The density distribution of the liquid in different layers is shown in FIG.

4

. It can be seen that the liquid density in the first layer around the nanoparticle is enhanced notably to 1367, 1470, 1540 and 1597 𝑘𝑔 𝑚

[START_REF] Subramanian | Nanofluids Their Eng. Appl[END_REF]

⁄ as the solid-liquid interaction strength is increased from 0.25 to 1.5 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙. Interestingly, the second layers also show higher densities which are 1331, 1356, 1369 and 1372 𝑘𝑔 𝑚 3 ⁄ as the solid-liquid interaction strength is increased.

  3, 𝜇 𝑦𝑧 is considered as the viscosity of liquid argon in each layer. The autocorrelation function of 𝜏 𝑦𝑧 for the liquid layers is shown in the SI, and the viscosity of each layer is calculated based on them. It is observed that the viscosity of the first liquid layer becomes 275, 442, 1031, and 2284 µ𝑃𝑎. 𝑠 as the solid-liquid interaction strength is increased from 0.25 to 1.5 𝑘𝑐𝑎𝑙/𝑚𝑜𝑙 while for the second layer, the corresponding viscosities are 182, 184, 190, and 190 µ𝑃𝑎. 𝑠. The viscosities of the other layers are around 184 µ𝑃𝑎. 𝑠 107 which is the expected viscosity for 1314 𝑘𝑔 𝑚 3 ⁄ . By normalizing the viscosity of each layer with respect to 184 µ𝑃𝑎. 𝑠, as done in Fig.

TABLE 5 .

 5 Calculated viscosity (µ𝑃𝑎. 𝑠) for the different liquid argon layers as a function of the interaction strength between the wall and liquid atoms

	𝜀 𝐶𝑢-𝐴𝑟 (𝑘𝑐𝑎𝑙 𝑚𝑜𝑙 ⁄	)	First	Second	Third	Fourth	Fifth
	0.25		275±2.1	182±2.7	173±6.6	172±2.1	174±2.0
	0.5		442±1.0	184±1.5	175±1.5	173±3.9	174±2.5
	1.0		1031±2.3	190±0.8	178±1.7	171±2.4	173±1.2
	1.5		2284±2.5	190±2.5	170±5.0	174±2.6	172±1.8

TABLE 6 .

 6 Calculated diffusion coefficient parallel to the nanoparticle surface (10 - 9 m 2 s -1 ) for the different liquid argon layers as a function of the interaction strength between the wall and liquid atoms

	𝜀 𝐶𝑢-𝐴𝑟 (𝑘𝑐𝑎𝑙 𝑚𝑜𝑙 ⁄	)	First	Second	Third	Fourth	Fifth
	0.25		1.87±0.42 3.02±0.15	3.26±0.16	3.28±0.16	3.34±0.13
	0.5		1.04±0.08 2.84±0.11	3.26±0.2	3.34±0.13	3.4±0.19
	1.0		0.75±0.09 2.76±0.14	3.22±0.13	3.3±0.18	3.34±0.18
	1.5			2.76±0.15	3.26±0.13	3.32±0.18	3.4±0.15

TABLE 7 .

 7 Calculated diffusion coefficient perpendicular to the nanoparticle surface (10 -9 m 2 s -1 ) for the different liquid argon layers as a function of the interaction strength between the wall and liquid atoms

	𝜀 𝐶𝑢-𝐴𝑟 (𝑘𝑐𝑎𝑙 𝑚𝑜𝑙 ⁄	)	First	Second	Third	Fourth	Fifth
	0.25		1.09±0.24	2.16±0.08	2.92±0.11	3.22±0.1	3.3±0.14
	0.5		0.75±0.06	1.94±0.1	2.84±0.1	3.22±0.1	3.32±0.19
	1.0		0.59±0.07	1.89±0.07	2.72±0.09	3.14±0.18	3.28±0.16
	1.5		0.5±0.04	1.83±0.09	2.72±0.12	3.2±0.15	3.28±0.19

ACKNOWLEDGMENTS

R.Rabani and M.H.Saidi wish to thank the financial support of the Iran National Science Foundation (INSF). L. Joly is supported by the Institut Universitaire de France.

SUPPORTING INFORMATION

Autocorrelation of the shear stress components in different layers of the liquid and the Evolution of the MSD and the diffusion coefficient of the liquid atoms are presented in this section.

CONFLICT OF INTEREST

There are no conflicts to declare.