Cyrine Ernandes

Hung-Ju Lin

Michel Mortier

Patrick Gredin

Mathieu Mivelle

Lionel Aigouy

Exploring the magnetic and electric side of light through plasmonic nanocavities

Light-matter interactions are often considered to be mediated by the electric component of light only, neglecting the magnetic contribution. However, the electromagnetic energy density is equally distributed between both parts of the optical fields. Within this scope, we experimentally demonstrate here, in excellent agreement with numerical simulations, that plasmonic nanostructures can selectively manipulate and tune the magnetic versus electric emission of luminescent nanocrystals. In particular, we show selective enhancement or decay of magnetic and electric emission from trivalent europium-doped nanoparticles in the vicinity of plasmonic nanocavities, designed to efficiently couple to either the electric or magnetic emission of the quantum emitter. Specifically, by precisely controlling the spatial position of the emitter with respect to our plasmonic nanostructures, by means of a near-field optical microscope, we record local distributions of both magnetic and electric radiative local densities of states (LDOS) with nanoscale precision. The distribution of the radiative LDOS reveals the modification of both the magnetic and electric optical quantum environments induced by the presence of the metallic nanocavities. This manipulation and enhancement of magnetic lightmatter interaction by means of plasmonic nanostructures opens up new possibilities for the research fields of opto-electronics, chiral optics, nonlinear&nano-optics, spintronics and metamaterials, amongst others.

KEYWORDS: Nano-photonics, electric and magnetic dipole transitions, luminescence enhancement, plasmonic cavities, electric and magnetic LDOS, near-field optics.

TEXT. The constant advances in nanophotonics have recently led to the development of optical nanostructures allowing to manipulate the emission of fluorescent emitters at the nanoscale. [1][2][3][4] However, the underlying interactions between quantum emitters and photonic nanostructures have often been assumed to be only mediated by the optical electric field, discarding the magnetic component of light. 5 Indeed, at optical frequencies, magnetic dipole (MD) transitions are at least two orders of magnitude weaker than electric dipole (ED) transitions. [START_REF] Cowan | The theory of atomic structure and spectra[END_REF][START_REF] Giessen | [END_REF][8] For this reason, most of the past studies involving single emitters (molecules [9][10][11] or quantum dots 12,13 ) have solely been focused on the modification of the excitation or emission properties of ED transitions. [14][15][16][17][18] Recently, it was shown that strong MD transitions could as well be found in a certain class of metallic materials such as rare earth ions. [19][20][21][22][23] Researchers have thus naturally explored the manipulation of the excitation and emission of the MD transitions from these ions. [23][24][25][26][27][28] In particular, in trivalent europium (Eu 3+ ), it was demonstrated that by using a simple metallic surface, spontaneous emission of a MD transition could be manipulated in the same way as an ED transition in the visible range. 21,29 In particular, we previously demonstrated that by changing the distance between a nanoparticle containing electric and magnetic dipolar transitions and a gold strip of micrometer size, we could engineer their respective electric and magnetic emission, 29 in the same way it was predicted and demonstrated by K. H Drexhage [START_REF] Drexhage | IV interaction of light with monomolecular dye layers[END_REF] and Chance et al. 14 This approach Nevertheless, in order to go further in the engineering of the interactions between magnetic light and matter, it is necessary to design nanostructures allowing to boost these interactions. [START_REF] Baranov | [END_REF][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] Along those lines, we recently demonstrated for the first time how dielectric nanoantennas were able to enhance the magnetic emission of quantum emitters 34 and researchers have demonstrated, through assemble measurements, that metallic nanostructured films could promote the transmission of magnetic emission. 26 But so far, only theoretical works have studied and demonstrated the ability of single plasmonic nanostructures [35][36][37][38] to manipulate the spontaneous emission of MD transitions by manipulating the magnetic LDOS at the nanoscale.

In this article, by coupling a nanocrystal that holds strong intrinsic magnetic and electric dipole transitions to plasmonic nanostructures, we experimentally demonstrate the manipulation of magnetic and electric light emission, in excellent agreement with numerical simulations. To this end, we consider the near field interaction between plasmonic cavities and a single Eu 3+doped nanoparticle. While our previous study 29 demonstrated manipulation of electric and magnetic emission in a longitudinal plane of a gold mirror, we showed in here that the luminescence associated to magnetic and electric transitions of the emitters can be modified and manipulated in the same way by changing the lateral position of the particle inside nanostructured plasmonic cavities. Moreover, by scanning the particle over the nanostructures, we recover from the experimental data the spatial distribution of the electric and magnetic relative radiative local densities of state (EMLDOS) at the surface of the cavities.

In order to find the appropriate parameters for which a plasmonic cavity can efficiently modify the electric and magnetic transitions of Eu 3+ nanocrystals, we first carried out Finite Difference Time Domain (FDTD) numerical simulations to explore the evolution of the electric and magnetic branching ratios for different types of cavities (see supporting information). To simplify the problem, we focused on linear, one-dimensional structures, formed by parallel grooves in a thin metal film. We show in Figure 1a the evolution of the electric and magnetic branching ratios ((relative weight of each transition (see supporting information)) for isotropic electric and magnetic dipoles centered on a cavity of variable length. In here, the relative experimental strength of each transition is not yet considered. The highest variations of the ED and MD transitions occur for cavity lengths of 230 nm, 570 nm and 870 nm, which correspond to the situation for which the nanoemitters are efficiently coupled to it. Based on these numerical simulations, we designed and nanofabricated (see supporting information) two cavities with lateral lengths of 570 nm and 870 nm, corresponding to the maximum coupling of the ED and MD transitions in their centers, respectively. Figure 1b depicts a scanning electron microscope (SEM) image of both plasmonic cavities (see supporting information for nanofabrication details). variable length (see supporting information). In here, the value of 230 nm will not be considered in this study. Indeed, the size of the europium-doped nanocrystals (see supporting information) is only slightly smaller, so only weak spatial variations could be expected for these structures due to average effects. (b) SEM image and schematic of the nanofabricated sample. Each cavity length is indicated above the image.

To study the interaction between the plasmonic cavities and the europium emitter, we used the experimental setup shown in Figure 2a. It consists of a sharp atomic force microscope (AFM) tip with a single Eu 3+ doped nanocrystal attached at its extremity (inset in Figure 2a). The tip is positioned in close proximity to the plasmonic nanostructures and both the tip and the sample are illuminated by a continuous 532 nm laser focused at oblique incidence (see supporting information).

The photoluminescence emission spectrum of the Eu 3+ nanocrystal far away from the surface and its associated band diagram are shown in Figure 2b andc. We can distinguish the three fundamental emission lines associated to the different electric transitions (ET) or magnetic transition (MT) in the nanocrystal: 20 Using the setup described in Figure 2a, we then scanned the cavities shown in Figure 1b with a Eu 3+ -doped nanoparticle and measured its emission spectrum at every position above the structures. Before quantitative analyzing the variations of each photoluminescence lines separately, we show in Figure 3a and c the evolution of the distribution of the total nanocrystal photoluminescence on the surface of the cavities. The excitation being made through an electric dipolar transition (λ = 532 nm), 24 these distributions are strongly related to the intensity of the electric near field located at the surface of the cavity and exciting the particle at this wavelength. This is corroborated by numerical simulations displayed in Figure 3b andd andobtained by exciting the cavities at λ = 532 nm by an incident plane wave launched from the upper side of the cavities, with an angle of 70 degrees with respect to the normal of the sample and with a TM polarization. The variations of the collected signal along the scanned area are in fairly good agreement with the simulations, and the maxima and minima of photoluminescence correspond well to the ones of the total theoretical electric field intensity. However, there are some differences, which can be attributed to the presence of the tip and the particle which were both not taken into account in the simulations. The observed optical oscillations indicate that the particle is placed in close enough proximity to the surface since it senses the stationary surface plasmon waves located inside the cavities. In Figure 3e and 3f, we show the normalized emission spectra of the nanocrystal at two different positions of the cavities indicated by colored dots in Figure 3a and3c. These spectra exhibit two peaks associated to the magnetic (labeled by an M) and electric (labeled by an E) dipole transitions of the Eu 3+ , respectively. As it can be seen, the spectra do not overlap perfectly according to the particle position inside the cavities. This result indicates that the relative importance of the magnetic and electric dipolar transitions (MT and ET) of the europium ions can manipulate only by changing the particle position inside the nanostructures. To evaluate the relative contribution of the electric and magnetic transitions to the total luminescence, we then calculated the branching ratios (β ",$ ) associated to these transitions.

Experimentally, the relative weight of the MT (β $ ) and ET (β " ) along the plasmonic structures are calculated as: 25,29 β ",$ p = ( ),* + ( -.-(+) (1) where p is the position of the particle along the cavity, I ",$ is the luminescence intensity of the electric (I " ) or magnetic (I $ ) transition and I 121 is the total luminescence intensity. We show in Figure 4a and c the evolution of the experimental electric and magnetic branching ratios for both cavities. These curves were obtained by averaging 20 scans along the plasmonic nanostructures (x-direction). The corresponding theoretical calculations of the branching ratios are displayed in Figure 4b and d (see supporting information). We can observe an excellent agreement between the numerical simulations and the experimental data: the spatial oscillations, the maxima and the minima of β " and β $ are all located, for each case, at the same positions inside the cavities. Meanwhile, in accordance with the result displayed in Figure 1, the small cavity exhibits a maximum value of β " in its center whereas the large cavity shows an opposite behavior, namely, a maximum value of β $ in the center (see supporting information for more examples). The small differences in amplitude between the experiments and the simulations are attributed to the presence of the tip inside the cavities, which might slightly change their behavior. Another example shown in Figure S3, made with another tip and particle but in the same experimental conditions, exhibits as well a variability in the amplitude with respect to Figure 4. This indicates that the probe is indeed most likely the variable factor causing this fluctuation. These results clearly demonstrate that the electric and magnetic emission of trivalent europium ions can be manipulated, in a controlled manner, by means of plasmonic nanostructures. The different emission behaviors of the electric and magnetic transitions, given by the branching ratios, are due to the non-spatial overlap of the electric and magnetic radiative LDOS at the cavities surface. To illustrate this relationship, we calculated the relative contribution of the electric and magnetic components of the radiative LDOS using the values of the experimental branching ratio, as we previously described in: 29 ρ ",$ = β ",$ . f $," β ",$ . f $," + β $," . f ",$ (2) where ρ ",$ are the electric (ρ " ) or magnetic (ρ $ ) relative radiative LDOS and f ",$ are the oscillator strengths for the electric (f " = 0.67) or magnetic (f $ = 0.33) transitions respectively (calculated in far field (Figure 2c)).

The spatial distributions of the radiative magnetic local density of state (further referred as MLDOS) and the radiative electric one (further referred as ELDOS) in the cavities plane are displayed in Figure 5 a-d. As expected, the ELDOS and MLDOS do not overlap in the cavities plane, and exhibit several maxima and minima situated at different positions. In Figure 5e and 5f, we show the evolution of the experimental ELDOS and MLDOS by averaging 20 scans along the cavities. Although these curves seem similar to those already shown in Figure 4a and4c, they represent the electric and magnetic radiative LDOS inside the cavities and are independent of the oscillator strengths of the transitions and therefore from the emitter considered. Thus, they are directly related to the electric and magnetic optical quantum environment of the nanocavities.

These results highlight the ability of such nanostructures to influence and tune the emission of rare-earth ions, by maximizing or minimizing either the electric or magnetic LDOS at different locations in their near-field. In summary, we demonstrated experimentally, in excellent agreement with numerical simulations, that the emission of a single Eu 3+ doped nanocrystal that contains magnetic and electric transitions, can be perfectly manipulated and controlled by plasmonic cavities.

According to the position of the trivalent europium nanoparticle inside the nanostructures, we showed that both ED and MD transitions emission can be manipulated with an increase, or alternatively a decrease, of their respective near-field luminescence. Furthermore, we retrieved the spatial distribution of radiative electric and magnetic parts of the local density of state at the surface of the cavities, which is a fundamental quantum characteristic of these structures. This study shows that the magnetic part of light can be engineered in the same way as its electric counterpart using metallic nanostructures and we are confident that this work will pave the way for the design of new plasmonic geometries enabling the boost of magnetic light and matter interactions, allowing applications in fields as diverse as spintronics, photochemistry, 50 sensing, 51 molecular chirality 52 amongst others.

Figure 1 .

 1 Figure 1. (a) Theoretical branching ratios associated with magnetic (green) and electric (red) isotropic dipoles, placed in the center of a cavity made of a 40 thick layer of aluminum and of

CavityFigure 2 .

 2 Figure 2. (a) Schematic of the near field experimental setup. An AFM tip holding a single Eu 3+ doped nanoparticle (150-200 nm diameter, see Figure S1 in the supporting information) is scanning a plasmonic nanocavity at a constant distance of about 20 nm. Both the particle and the cavity are excited by a 532 nm laser, and for each position of the emitter with respect to the nanostructure, a spectrometer is recording the emission spectrum of the nanosource. Inset: SEM image of a single nanocrystal doped with Eu 3+ and attached at the apex of a tungsten AFM tip. (b) Partial band diagram of the Eu 3+ ions. (c) Emission spectrum of the Eu 3+ doped nanoparticle recorded in far field.

Figure 3 .

 3 Figure 3. Two-dimensional distribution of the total collected photoluminescence (a,c) and total simulated electric intensity (λ = 532 nm) (b,d) in the cavities plane, for the cavity lengths of 570 nm (a,b) and 870 nm (c,d). (e,f) Normalized emission spectra of MT (1) and ET (2) measured at two different positions inside each cavity lengths, (e) 570 nm and (f) 870 nm, in each case along the same transversal dashed line and indicated by blue and pink dots in (a,c).

Figure 4 .

 4 Figure 4. Magnetic and electric branching ratios. (a,c) Experimental branching ratios associated with magnetic (green) and electric (red) transitions of the Eu 3+ according to the nanoparticle positions inside the 570 nm (a) and 870 nm (c) long cavities. The red and green dots represent the experimental data, the full line is a polynomial fit and the error bars represent the standard of the mean. (b,d) Numerical counterparts of (a,c) respectively.

Figure 5 .

 5 Figure 5. Radiative local density of state. 2D experimental maps of the radiatives (a,b) MLDOS and (c,d) ELDOS in the plane above the 570 nm (a,c) and 870 nm (b,d) long cavities. (e,f) Experimental magnetic (green) and electric (red) radiative LDOS for each cavity lengths, (e) 570 nm and (f) 870 nm. Dots are the experimental data, full lines represent a polynomial fit and the error bars are the standard error of the mean.

F 2 from 600-625nm and ET(3) 5 D 0 7 F 4 from 685-705nm. It is important to mention that, for the present work we will ignore the contribution of the last ET (3). Since it is well separated from the higher energy ones, it was found to be weakly influenced by the environment, as observed in previous studies.29 And because we chose to focus mainly on the spectrally close magnetic (MT(1)) and electric (ET(2)) dipolar transitions, in this study, all other decay channels from 5 D 0 to the ground state, including ET(3), are considered as losses. Using this spectrum recorded in far field, we then calculated the effective values of the oscillator strengths of MT(1) and ET (2), found to be respectively 0.33 and 0.67. These values will be then considered later on in the simulations, in order to quantitatively compare our experimental and theoretical results.

ACKNOWLEDGMENTS

The authors thank Maria F. Garcia-Parajo for fruitful discussions and her help on the nanostructuration of the sample. M.M acknowledges support from the French "Investissements d'Avenir" program (Labex MATISSE), the DIM Nano-K program from "Région Ile de France", the CNRS emergence program and the ANR tremplin under reference ANR-17-ERC3-0006-01.

The authors declare no competing financial interest.


Patrick Maillé

Bruno Tuffin

Analysis of Sponsored Data Practices in the Case of Competing Wireless Service Providers

With wireless sponsored data, a third party, content or service provider, can pay for some of your data traffic so that it is not counted in your plan's monthly cap. This type of behavior is currently under scrutiny, with telecommunication regulators wondering if it could be applied to prevent competitors from entering the market, and what the impact on all telecommunication actors can be. To answer those questions, we design and analyze in this paper a model where a Content Provider (CP) can choose the proportion of data to sponsor and a level of advertisement to get a return on investment, with several Internet Service Providers (ISPs) in competition. We distinguish three scenarios: no sponsoring, the same sponsoring to all users, and a different sponsoring depending on the ISP you have subscribed to. This last possibility may particularly be considered an infringement of the network neutrality principle. We see that sponsoring can be beneficial to users and ISPs depending on the chosen advertisement level. We also discuss the impact of zero-rating where an ISP offers free data to a CP to attract more customers, and vertical integration where a CP and an ISP are the same company.

Introduction

With the improving capacity of smartphones, wireless data consumption is exponentially increasing. According to previsions 1 , mobile data consumption will be sevenfold larger by 2021 than by 2017. Wireless communication subscription offers are Patrick Maillé IMT Atlantique / IRISA, Rennes, France, e-mail: patrick.maille@imt.fr Bruno Tuffin Inria, Rennes, France e-mail: bruno.tuffin@inria.fr often made of unlimited telephony but caps on data over which a volume-based fee is applied. Users may therefore be forced to monitor and even limit their data usage.

Content providers (CPs) on the other hand are economically dependent on the amount of data consumed by users, typically via displayed advertisements: the more volume consumed, the more advertisement and then revenue. For this reason, there might be an incentive for content providers to sponsor data, that is, to at least partially pay the Internet Service Provider (ISP) for the user's consumption of its own services. This way, users may expand their usage, CPs may get more revenue due to advertisement, and ISPs are paid for the traffic increase. Examples in the USA are: Netflix or Binge-On with T-Mobile, DIRECTV and U-verse Data Free TV with AT&T, etc. ISPs are even offering services such that CPs can, if they wish, apply sponsored data: examples are Verizon Wireless, AT&T Mobility and T-Mobile US, or Orange with DataMI in France, among others.

But sponsored data bring concerns from user associations and small content providers. It is claimed to give an unfair advantage to some content/service providers and to eventually prevent some actors from entering the market, due to lower visibility and high entrance costs to the sponsored data system if they want to get the same service as incumbent and big providers already present. The culmination of the sponsored data principle is the so-called zero-rating where ISPs (freely) remove some content providers from their data caps, hoping more customers will subscribe due to potentially unlimited usage. Those providers could then hardly be challenged. The French ISP SFR did it with YouTube in its offer RED a few years ago. Similarly, Facebook, Google and Wikipedia have built special programs in developing countries, with the claimed goal to increase connectivity and Internet access. For those reasons, sponsored data and zero-rating are currently under investigation by regulators to determine if rules should be imposed to prevent or limit their use. In Chile for example, the national telecom regulator has stated that it was violating net neutrality laws and that it should be forbidden. Net neutrality means that all packets/flows are treated the same, independently of their origin, destination, and type of service [START_REF] Maillé | Internet governance and economics of network neutrality[END_REF][START_REF] Maillé | Toward a net neutrality debate that conforms to the 2010s[END_REF].

The purpose of this paper is to design and analyze a mathematical model of sponsored data to understand if this type of behavior should be banned or encouraged. To our knowledge, very few similar works exist, that we now summarize and from which we differentiate ourselves. The notable works are first [START_REF] Joe-Wong | Sponsoring mobile data: An economic analysis of the impact on users and content providers[END_REF], where the authors consider a model with a discrete set of users, a single ISP, and several (complement) CPs able to choose a different level of sponsoring per user. They show that sponsoring can benefit more to users than to CPs. Our model presents similarities with that one but strongly differs on many aspects: we use a different model for CPs' advertisement level and determination of the best strategy; our model also represents the negative externality of advertisements on users, something often forgotten, and a different model of volume consumption highlighting the correlation between willingness-to-pay for connectivity and consumed volume of data; finally and more importantly, we model competition between ISPs to better link the sponsored data practice to the network neutrality debate. We can also mention [START_REF] Somogyi | The economics of zero-rating and net neutrality[END_REF], dealing with several substitutable CPs in competition, or [START_REF] Xiong | Competition and Cooperation Analysis for Data Sponsored Market: A Network Effects Model[END_REF][START_REF] Zhang | Sponsored data plan: A two-class service model in wireless data networks[END_REF] including network externalities, but again a single ISP, or [START_REF] Pang | When data sponsoring meets edge caching: A gametheoretic analysis[END_REF] which nicely combines sponsored data and caching strategies, but not considering any ISP in the equation.

Our contributions are the following: i) the use of a model different from the literature to represent user behavior. This model separates the willingness-to-pay for connectivity and for CP content. Remark that using a new model is important to better understand the assumptions that may change the results, or see the robustness of conclusions. ii) We focus on CPs' strategical decisions and their impact on all other actors (ISPs and users). The decision variables are the amount (proportion) of data to be sponsored and the amount of advertisement CPs will set; our study allows us to see if sponsoring data also means more advertisement for a return on investment for CPs. iii) We investigate the best sponsoring and advertisement strategies in the case of competition between ISPs; three situations are compared: (a) no sponsored data (b) the same level of sponsoring at all ISPs (c) the possibility to sponsor differently according to the ISP; these situations also correspond to three different levels of neutrality in the net neutrality debate. iv) We investigate and compare with all other situations the case of zero-rating, when the CP content is not counted in the user data cap; ISPs could indeed be in favor of such a strategy if it brings additional subscriptions. v) Finally we look at the case when the CP is managed by one ISP, that is the so-called vertical integration; does it lead to an unfair advantage if data are sponsored with respect to a non-sponsored situation? This corresponds to practical situations within an increasingly vertically-integrated ecosystem.

The remainder of the paper is organized as follows. Section 2 presents the model, the decision variables and their hierarchy. Section 3 discusses the last level, how users distributes themselves among ISPs depending on prices and strategies of CPs. Then Section 4 describes the three sponsoring strategies and compares their impact on all actors in the case of independent CPs and ISPs. The particular situation of zero-rating, with the CP content not included in one ISP data plan but included in the other is investigated in Section 5. Section 6 presents the case of a verticallyintegrated CP-ISP. Finally, from all the described scenarios, we conclude in Section 7 on the need to regulate sponsored data.

Model

We consider three types of actors: users, CPs and ISPs. We describe here the model with full generality, and will later restrict the number of actors to focus on specificand to our knowledge un-investigated-aspects of the sponsored data debate.

Consider M CPs indexed by j and N ISPs indexed by i.

Users

We index users by θ ∈ R + , which characterizes their type. We consider a continuum of users of total mass 1, and denote by F the cumulative distribution function of θ .

Each user θ has several variables of decision: its choice of ISP, that we will denote by i(θ ), but also the volume of CP j data (without advertisements) he will consume if using ISP i, which we will denote by v i, j (θ ). The dependence on i comes from the fact that prices for data may differ between ISPs, and the consumed volume too as a consequence.

To simplify, let θ be the willingness-to-pay for connectivity and p i the subscription price at ISP i, which would be giving a "utility for connectivity" θp i . But we additionally weigh the difference using an ISP reputation parameter a i ≥ 0, meaning that users value more the gain obtained with incumbent providers than with newcomers, resulting in a "weighted utility for connectivity"

a i (θ -p i ).
In addition, User θ gains satisfaction from using each CP j. The cost per unit of volume he pays in his data plan is denoted by c i, j , which depends on j because that usage can be sponsored by CP j as we will see in the next subsection. Moreover there is a valuation for each unit of volume of CP j. Let r θ , j be the marginal valuation, that we will consider to be

r θ , j (x) = [θ -(α j s 2 j )x] +
for the x-th unit of useful volume, i.e., without advertising, where [y] + := max(y, 0), α j is a fixed parameter, and s j ≥ 1 corresponds to the relative increase of volume due to advertisement displayed by CP j, expressed as the total downloaded volume divided by the volume of "useful" data, excluding ads. The larger this value of s j , the smaller the service valuation for users. We consider a square value of s j to assure the reasonable assumption that at one point, for the CPs, the loss due to users unpleased by advertisements exceeds the gain from those ads (see later ( 2) and ( 3) with s j tending to infinity); any other choice can be used without difficulty though. Then the willingness-to-pay r θ , j (x) for consuming a volume x of CP j data over a month is

r θ , j (x) =    θ x -α j s 2 j x 2 2 if x ≤ θ α j s 2 j θ 2 2α j s 2 j otherwise.
This form is commonly adopted in [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF][START_REF] Maillé | Users facing volume-based and flat-rate-based charging schemes at the same time[END_REF] and presents the advantage of positively correlating willingness-to-pay for data and willingness-to-pay for connectivity. We have added in the model the negative correlation due to advertisement (some kind of "pollution" which can advert users). Note that a positive effect (only) from advertisement was considered in [START_REF] Joe-Wong | Sponsoring mobile data: An economic analysis of the impact on users and content providers[END_REF] due to relevant ads which can be clicked, but ads seems to us rather negatively perceived.

Our assumptions give a utility for user θ at ISP i

U i (θ ) := a i (θ -p i ) + M ∑ j=1 r θ , j (v i, j (θ )) -c i, j v i, j (θ )s j , (1)
using the fact that the total downloaded volume is actually v i, j (θ )s j from the definition of s j . Users indeed also download advertisements, which are not differentiated from "real" content by ISPs. The additive expression indicates that CPs are assumed to be independent in terms of content.

Our modeling, unlike those of the literature ( [START_REF] Joe-Wong | Sponsoring mobile data: An economic analysis of the impact on users and content providers[END_REF] for example), chooses to represent user preferences for ISPs, in terms of connectivity and a correlation between willingness-to-pay for connectivity and consumed volume of data; it separates what comes from ISPs from what comes from CPs.

If subscribing to ISP i, the volume v i, j (θ ) User θ chooses is the one maximizing r θ , j (v i, j (θ ))c i, j s j v i, j (θ ), and can easily be computed [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF] to be

v i, j (θ ) = θ -c i, j s j α j s 2 j + , (2)
leading to r θ , j (v i, j (θ )) -c i, j s j v i, j (θ ) = (θ -c i, j s j ) 2 2(α j s 2 j )
1l {θ >c i, j s j } . Summarizing, User θ chooses ISP i(θ ) maximizing his overall utility:

i(θ ) = argmax i U i (θ ) = argmax i a i (θ -p i ) + M ∑ j=1 (θ -c i, j s j ) 2 2(α j s 2 j ) 1l {θ >c i, j s j }
if the max is non-negative, otherwise i(θ ) = 0, meaning no subscription at all.

CPs

CP j is assumed to have a utility (a revenue) linearly increasing with the volume of displayed advertisement, with CP-dependent linear parameter β j . The advertisement volume is the total volume s j v i(θ ), j (θ ) minus the "real" data volume v i(θ ), j (θ ), hence (s j -1)v i(θ ), j (θ ). As a consequence, the gain for CP j is

θ β j (s j -1)v i(θ ), j (θ )dF(θ ).
Each CP j can decide to sponsor a fraction γ i, j of the data usage cost of ISP i's users. This could be an incentive to consume more CP j content and therefore generate more revenue from advertisement. Of course, there is a trade-off, which is one purpose of this paper, between sponsoring cost and generated volume of consumed data. We will investigate three different sponsoring policies:

1. No sponsored data: γ i, j = 0 ∀i, j; 2. The same data sponsoring level for all ISPs: γ i, j = γ j ∀i, j;

3. A possible differentiation between ISPs, with γ i, j = γ i , j for i = i .

In terms of the network neutrality debate, those three scenarios correspond to three different levels of neutrality. Option 1 corresponds to full neutrality. Option 2 seems neutral with respect to ISPs, even if non-neutral between applications because different applications can have different levels of sponsoring (note on the other hand that advocates of sponsored data also claim that packets are treated the same within the network, but it is not the issue we want to address here). Option 3 is the fully non-neutral one.

Let q i be the unit price ISP i is (officially) charging users for data. After data sponsoring, the unit cost seen by users when consuming data of CP j is then

c i, j = q i (1 -γ i, j ).
If sponsoring data, CP j has to pay to each ISP for the proportion of volume it has chosen to sponsor, hence a total cost: θ γ i(θ ), j q i(θ ) s j v i(θ ), j (θ )dF(θ ), where the proportion of volume paid, is the proportion γ j of the total volume s j v i(θ ), j (θ ), (which includes ads).

Compiling all the elements, the revenue of CP j is

G j = θ (β j (s j -1) -γ j q i(θ ) s j )v i(θ ), j (θ )dF(θ ). (3)
CP j has several decision variables: the sponsoring levels γ i, j ∀i and the advertisement volume increase level s j ≥ 1. Note that we could also consider advertising levels that would depend on the user ISP (i.e., s i, j instead of s j ): we rather impose an equal advertising level, since our focus is on the sponsoring strategies, and CPs may not know the ISP of the user upon a request for content.

ISPs

Each ISP i tries to maximize its revenue. Revenue comes from subscriptions and consumed data:

R i = θ p i + ∑ j q i s j v i, j (θ ) 1l {i(θ )=i} dF(θ ).

What we analyze

Our purpose in this paper is to analyze the impact of the decisions of CPs on all actors, in the context of competing ISPs. In particular, we would like to see the impact of a "neutrally" sponsoring CP (that is, γ i, j = γ i , j for any ISPs i, i ) with respect to sponsoring differently depending on ISPs and with respect to no sponsoring.

Three scenarios are analyzed in the next sections:

1. In the case of independent ISPs and CPs, what are the best advertising and sponsoring levels of CPs and the impact on all actors? We can compare the outputs with optimal solutions for the three sponsoring possibilities. 2. Does zero-rating, where an ISP leaves a CP out of its data plan to try to attract customers, harm users and the other ISPs? Should it be forbidden? 3. Does vertical integration of a CP and an ISP have a positive or negative impact on the other ISPs and on users? This, with a dominant or a non-dominant ISP.

Since CPs are complements (due to the additive form of users' utility), we will consider without much loss of generality in the next sections the case of a single CP. This in order to focus on our purpose, that is the impact of advertisement and sponsoring strategies on competition between ISPs, and on users. Hence we can simplify the writing by removing the CP-relative indices in the notations defined before. We will also limit ourselves to 2 ISPs, labeled 1 and 2.

User subscription decisions

User decisions depend on prices set by ISPs, but also on decisions of the CP. Even if the CP "plays" first, its optimal strategy should anticipate users' choices. For this reason, we discuss here the repartition of users depending on fixed CP decisions. The optimal CP decisions are analyzed in the next sections, anticipating this subsequent choice (this implicitly assumes that the CP knows the distribution F).

In other works using the same type of models for user preferences [START_REF] Cho | Pricing for mobile data services considering service evolution and change of user heterogeneity[END_REF][START_REF] Maillé | Users facing volume-based and flat-rate-based charging schemes at the same time[END_REF], we had depending on the value of θ an interval over which users do not subscribe, then an interval (possibly empty) over which one ISP is chosen, and then another one for the other ISP (that is, users with small valuations do no subscribe, those with "intermediate" valuation go with one ISP and those with the largest go with the other). Here things can be more complicated, with interestingly more intricate intervals. To illustrate this, consider the values a

1 = 1.3, a 2 = 1, s = 1, α = 1, c 1 = 2, c 2 = 1, p 1 = 1,
and p 2 = 0.97, i.e., ISP 1 has a better reputation but higher subscription and usage prices. Figure 1 displays the difference max(U 1 (θ ), 0)max(U 2 (θ ), 0) in terms of θ . We consider the max with 0 to show only the situations where a user θ is willing to subscribe to an ISP.

We can see (when the curve gives zero) that users with small valuations θ prefer not to subscribe to any ISP; then there is an interval where ISP 2 is preferred, then an interval where it is ISP 1, and then again ISP 2. An interpretation is as follows: users with small valuation (but large enough to subscribe) are first interested in connectivity with a smaller connectivity price at ISP 2. At one point, due to its larger reputation, connectivity attractiveness of ISP 1 becomes larger, but it is counterbalanced by the data cost of ISP 2 (the parabola).

In general, in terms of the combination of prices, we can express how many such successive intervals above the "no subscription" possibility can be found, between

(θ )] + -[U 2 (θ )] + Fig. 1 [U 1 (θ )] + -[U 2 (θ )
] + in terms of θ . In the left zone (small θ values), users do not subscribe to any provider.

one and three. We do not give such details here since it is a list of cases which is not very instructive. The situation becomes even more intricate if we increase the number of ISPs (and CPs). Solving it numerically in the next section does not pose any problem.

Analysis with independent CP and ISPs

In all our numerical investigations, θ follows an exponential distribution with rate 1 truncated at 2, i.e., F(θ ) = (1e -θ )/(1e -2 ) for θ ∈ [0, 2]. Unless specified otherwise (i.e., when a parameter explicitly varies), the parameter values we consider throughout all the paper are summarized in Table 1. ISP 1 has the position a 1 a 2 p 1 p 2 q 1 q 2 β α 1.5 1 1 0.6 0.6 0.4 1 1

Table 1 Parameter values for the numerical investigations of an incumbent: with a better reputation than its competitor, it can charge higher subscription and per-usage prices. We discuss in this section the case where the ISP and CP are independent. We display in Figure 2 the CP revenue as s varies for fixed values of the other parameters, including fixed values of the sponsoring coefficient γ i on ISP i = 1, 2, but using the subsequent user subscription decision. It can be checked that there is an optimal level of advertising s * around 2.7 above which increasing the displayed ads diminishes revenue.

In Figure 3 we display the CP revenue-maximizing sponsoring levels γ i as s varies, and the corresponding CP Revenue. We can see that the optimal γ i is in-

s γ γ 1 = γ 2 γ 1 γ 2 2 4 0 2 4 •10 -2 s CP revenue G γ 1 = γ 2 = 0 Opt. γ 1 = γ 2 Opt. (γ 1 , γ 2 )
Fig. 3 CP revenue-maximizing sponsoring coefficient γ (left), and corresponding CP revenue (right) in each scenario.

creasing with the level of advertisement s: the more ads you put, the more you can sponsor data because you earn more, and the more you have to do to compensate the negative externality for users. Typically also, applying the same level of sponsoring (γ 1 = γ 2 ) leads to a sponsoring level in between the differentiated levels of sponsoring (γ 1 = γ 2 ). Similarly to Figure 2 where sponsoring was fixed, in the present case of optimal sponsoring for each s, there is an optimal advertising level for each sponsoring scenario in Figure 3 (right). As expected, the CP revenue with differentiated sponsoring is larger than the one with identical sponsoring, itself larger than the one with no sponsoring, because the optimization is on a larger set of sponsoring parameters. When the advertisement level is low, there is no significative difference between the three options; it is significant when s is large. The optimal advertising levels s * are 1.34 with no sponsoring, 1.34 too with identical sponsoring, and 1.44 with differentiated sponsoring, a larger value. "No sponsoring" or "identical sponsoring" yield a similar optimal revenue, while it is significantly better for differentiated sponsoring which then seems the appropriate option for the CP.

γ 1 = γ 2 = 0 γ 1 = γ 2 Opt. (γ 1 , γ 2 )
Fig. 4 ISP revenues (left) and Consumer Surplus (right) with CP-revenue-maximizing γ values.

Figure 4 displays the impact of sponsoring on the revenues of ISPs and on consumer surplus. In terms of ISP revenue first, we can see that ISP 1, the one with the highest reputation, makes more money than ISP 2. For very low values of s the three policies give the same revenues. ISP 2 always prefers sponsoring to no sponsoring, while for ISP 1 an identical sponsoring is always the best option, something not so obvious at first sight. As a second choice for ISP 1, depending on the advertisement level, preference varies between no sponsoring and a differentiated one. Sponsoring is also always preferred by users, since giving a higher consumer surplus. On the other hand, the consumer surplus decreases as the advertisement level increases: the increased sponsoring level does not compensate enough the negative externality of ads. Interestingly, there is no significant difference for users between the two sponsoring options, and no dominance either.

From the figures, sponsoring seems a relevant option for all actors, with some differences in preferences between the two sponsoring options depending on the actors.

Table 2 compares the outputs at the optimal advertising levels selected by the CP. Note that in the case of identical sponsoring policy, the optimal advertising level is

s * γ i G R 1 R 2 CS
No sponsoring 1.34 (0,0) 0.0378 0.4124 0.0931 0.3813 Id. sponsoring 1.34 (0,0) 0.0378 0.4124 0.0931 0.3813 Dif. sponsoring 1.45 (0,0.41) 0.0411 0.3317 0.1559 0.3770 Table 2 Comparison of the three sponsoring scenarios with CP-revenue-maximizing advertisement and sponsoring levels, for the parameter values of Table 1.

such that there is no sponsoring at all. Hence no sponsoring and identical sponsoring lead to the same output with our set of parameters. Differentiated sponsoring leads to slightly more ads. With this optimal choice, only customers of ISP 2 get sponsored data, 41% of data being sponsored. It leads to a larger CP revenue than the two other options. The revenue of ISP 1 decreases because attracting less customers, to the benefit of ISP 2. Consumer surplus is then slightly decreased. In conclusion, with this set of parameters and optimal advertising, identical sponsoring is of no help since leading to no sponsoring at all, and differentiated sponsoring just slightly changes consumer surplus but really modifies ISP revenues to the advantage of the one with the lowest reputation; this could be beneficial for competition with incumbent providers having usually a better reputation than newcomers.

Zero-rating

Zero-rating consists in leaving a CP out of the data plan, making it free of access for users. The potential interest for an ISP is to attract more customers, and therefore compensate the volume-based revenue loss by a subscription revenue increase. We study here whether this kind of strategy makes sense and what is the impact on all actors. Implementing zero-rating translates into c i, j = 0. In the case of a single CP, it is equivalent to considering q i = 0.

Figures 5 to 8 alternatively display the outputs in the cases where ISP 1 or ISP 2 implement zero-rating. Looking at the two scenarios helps to see if zero-rating is more relevant and has more impact depending on the market position of the ISP. Indeed, by keeping all values of parameters the same as in the previous section, ISP 1 is more established than ISP 2.

We see in Figure 5 (left) that in the case of zero-rating for ISP 1, the optimal sponsoring strategy for ISP 2 does not impact much the CP revenue, whatever the level of advertising. With respect to no zero-rating in Figure 3, the CP gain is substantial. Similar results can be observed when it is ISP 2 which applies zero-rating: substantial gain for the CP, but at a smaller advertising level than when it is applied by ISP 1, while at high advertising levels the difference between no sponsoring or sponsoring for ISP 1 is significant. At the optimal advertising level, there is no dif-ference though. In Figure 6, we can see that, again, the sponsoring level increases with the level of advertisement. But it can be checked that for smaller advertising levels, there is no sponsoring (which is understandable due to the low(er) revenue of the CP based on visits) and that sponsoring (for ISP 2 users) starts at a smaller advertising level for zero-rating with ISP 1 than in the opposite situation. In terms of ISP revenue in Figure 7, a zero-rating at ISP i induces for small advertising levels no revenue at the other ISP, which can be seen as against competition and could lead to action from a regulator, especially given that the optimal advertising levels are actually quite low from Figures 5 to 8. Sponsoring tends to reduce the gap in revenue with respect to no sponsoring when zero-rating is applied at ISP 1 only, while it is the opposite when applied at ISP 2.

Interestingly in Figure 8, consumer surplus does not change much with the policy if zero-rating is applied at ISP 1, the largest ISP. It changes a bit more when it is Advertisement level s Consumer Surplus γ 2 = 0, zero-rating for ISP 1 γ 2 > 0, zero-rating for ISP 1 γ 1 = 0, zero-rating for ISP 2 γ 1 > 0, zero-rating for ISP 2

Fig. 8 Consumer Surplus with zero-rating for one ISP.

applied at ISP 2, with large levels of advertising only. In both cases anyway, note that sponsoring is beneficial to users, but at the advertising level chosen optimally by the CP, there is no real difference (because the selected value is γ 1 = 0 in case of zero-rating at ISP 2).

Vertical Integration

We know look at the case of vertical integration: when the CP and one of the ISPs are the same firm. This situation becomes increasingly frequent in practice, with ISPs offering video-on-demand services, among other. In that case the optimization of parameters (sponsoring and advertising) is based on the combined revenues of the CP and the ISP, by just adding them. Figure 9 (left) displays the revenue of the integrated CP-ISP in terms of s, in the cases where it is either ISP 1 or ISP 2 which is integrated, and with fixed sponsoring levels. ISP 1 being more established, it is not surprising that integrating it yields a larger revenue. But the optimal advertising level is larger in the case of ISP 2 with those parameters. Figure 9 (right) shows the revenue of the integrated entity (CP-ISP 1, or CP-ISP 2) with the corresponding optimal γ. Recall that in the case of non-integration, the optimal advertisement level was 1.35 for no or identical sponsoring, and 1.45 for differentiated sponsoring. We still have the dominance such that differentiated sponsoring is better than identical one, itself better than no sponsoring. It can be seen that for an integrated CP-ISP 1, identical sponsoring leads to the most ads (a level around 1.6, and 1.0-i.e., no ads-for no sponsoring and 1.4 for a differentiated one), while it is differentiated sponsoring in the case of an integrated CP-ISP 2 (a level around 1.65, and 1.05 both for no sponsoring and for an identical one). So the optimal advertising level maybe less or more for the same sponsoring policy if the CP is integrated or not. But if choosing the optimal policy, advertising is larger if vertical integration is applied.

γ 1 = γ 2 = 0 γ 1 = γ 2 Opt. (γ 1 , γ 2 )
Fig. 9 CP+ISP revenue in the two integration scenarios, when γ 1 = 0.7, γ 2 = 0.4 (left), and with optimized values of (γ 1 , γ 2 ) (right). Figure 10 gives the consumer surplus in terms of s for the two integration cases, and compares it with the optimal non-integrated strategy obtained in Section 4. It allows to see if integration has a positive or negative impact on users, something of interest for the regulator. We observe that vertical integration is good for users, especially if sponsoring data is allowed. An identical sponsoring is the best option with this regard. The reason is that it induces much more data sponsoring resulting in lower user cost.

Finally, Figure 11 displays the revenue of the non-integrated ISP for all strategies, in order to illustrate whether integration brings down competition. We also again plot the revenues in cases of no integration for comparison sake. An integrated ISP 1 reduces drastically ISP 2 revenue, which was expected, but at a lesser extent for the identical sponsoring strategy because of less freedom. Integrating the incumbent ISP is typically against competition. Surprisingly, when it is ISP 2 who is integrated and with identical sponsoring, ISP 1 can benefit from the integration of the competitor. Note it is not the case for the other sponsoring strategies.

Table 3 compares the outputs at the optimal advertising levels selected by the CP in the two integrated scenarios. They can be compared with the non-integrated case of Table 2. Vertical integration, except partially for the last line of Table 3, 3 Output at optimal s when the CP is integrated with an ISP. always significantly benefits to users and competitors, something not so intuitive. We see in Table 3 that vertical integration leads to more sponsoring than in Table 2, so that consumer surplus is increased. With an integrated CP-ISP 2, the increased satisfaction in the first two lines is due to the smaller optimal advertising level; for the last line, the advertising level is larger but only partially counter-balanced by the increased sponsoring.

CP integrated with ISP 1 CP integrated with ISP 2 s * (γ 1 , γ 2 ) G + R 1 R 2 CS s * (γ 1 , γ 2 ) G + R 2 R

Conclusions

The purpose of this paper was to study the impact of sponsored data, differentiated or not, on all telecommunication actors, in the presence of competing ISPs. Our main results are the following:

• Sponsoring can be beneficial to users and ISPs and the impact depends on the chosen advertisement level.

• Zero-rating at one ISP can lead to no revenue at the other ISP, which can be seen as harming competition. • Vertical integration can be beneficial for all actors in this context, if the CP is integrated to the incumbent ISP (the one with the largest reputation).

Conclusions of course depend on the arbitrarily chosen set of parameters, but the paper provides insight on the potential effects of sponsored data, and a model to perform such an analysis. We plan to extend the promising first results of this work in one main direction: We aim at inserting another round of decisions, on ISP prices (for subscription and for data usage). Our purpose here was to see the impact of having several ISPs on the CP decisions leading to important conclusions. Inserting a pricing game will complete the analysis but is computationally demanding. Other interesting extensions include considering a larger set of ISPs to show that the obtained results are without loss of generality in this respect, as well as considering CPs which can be substitutes.

Fig. 2

 2 Fig.2CP revenue, with γ 1 = 0.7, γ 2 = 0.4.

γ 1 =

 1 γ 2 = 0 γ 1 = γ 2 opt. (γ 1 , γ 2 )

Fig. 5

 5 Fig.5CP revenue with zero-rating for ISP 1 (left) or for ISP 2 (right).

γFig. 6

 6 Fig.6Optimal γ for an ISP, with zero-rating for the other ISP.

Fig. 7

 7 Fig.7ISP revenues with zero-rating for ISP 1 (left) and for ISP 2 (right).

1 = 1 =

 11 Fig.10Consumer surplus with integrated CP-ISP 1 (left) or integrated CP-ISP 2 (right), for β A = 1.00, β B = 1.50, r 1 = 1.50, r 2 = 1.00, p 1 = 1.00, p 2 = 0.60, q 1 = 0.60, q 2 = 0.40, α A = 1.00, α B = 1.00, s B = 5.00.

γ 1 =Fig. 11

 111 Fig. 11 Revenue of the non-integrated ISP if ISP 1 (left) or ISP 2 (right) is integrated.

  . 1.4 (0.783,0.995) 0.6697 0.1925 0.4882 1.65 (0.243,0.997) 0.4644 0.2021 0.3978 Table

see http://www.businessinsider.fr/us/mobile-data-will-skyrocket-700-by-2021-2017-2/ among others 1


Pauline Simonnin

Virginie Marry

Benoit Noetinger

Carlos Nieto-Draghi

Benjamin Rotenberg
email: benjamin.rotenberg@sorbonne-universite.fr

Mineral-and Ion-Specific Effects at Clay-Water Interfaces: Structure, Diffusion and Hydrodynamics

We use molecular dynamics to investigate how the structure, diffusion and hydrodynamic properties of clay interfaces with aqueous solution depend on the nature of the clay, the nature of the counterions and the salt concentration in the solution. Specifically, we study water-filled nanopores between uncharged (pyrophyllite) and charged (montmorillonite and beidellite, with susbtitutions located in the octahedral and tetrahedral layers, respectively) clays, with sodium or cesium as counterions, in the absence and in the presence of added salt. We discuss how the balance between solvation and attraction of the cations to the surface results in various distributions between innerand outer-sphere complexes and how this influences the dynamics of water near the surface, as well as the hydrodynamic flow in the presence of an external force. In the latter case, the discussion based on mapping the molecular velocity profiles to a continuous description (parabolic Poiseuille flow) shows that the larger effects come from the presence/absence of charge in the mineral, as well as the localization of substitutions within the clay layer. The salt concentration and the nature of the counterions have a comparatively less important impact far from the surface -even though some differences are observed in its close vicinity, which are not properly captured by the continuous description.

Introduction

The transport of fluids through nanopores plays a role in a many environmental, biological and technological processes, ranging from clean water production or the harvesting of blue energy using membranes 1 to oil and natural gas recovery, 2 carbon sequestration 3 or the geological disposal of radioactive waste. 4 Several fundamental questions arise for the understanding and prediction of the transport properties in these contexts, such as the validity of continuum hydrodynamic theories when the confining length scale becomes comparable to the molecular size, as well as the dominant role of interfacial effects.

Molecular simulations played a key role in addressing the first aspect with the finding that continuous hydrodynamics holds down to remarkably small length scales, 5 even though it breaks down beyond 1-2 nm in particular due to the layering of the fluid at the interface, or for pore sizes of the order of that of the confined fluid molecules. 6 They also provide molecular scale information on the local viscosity of the confined fluid [7][8][9][10] as well as the relevant boundary conditions (stick or slip) which should be applied for a given interface, depending on the nature of both the solid wall and the confined fluid, with the quantification of the slip length [11][12][13] and its correlation with other properties such as the contact angle. 14 Due to their importance in the above-mentioned geological applications, the dynamics of fluids in nanopores of clay and related natural minerals [15][16][17][18][19][20] have attracted a lot of attention from the experimental and modelling point of view. We will not discuss here the interlayer pores, where only 1 or 2 fluid layers are present (see e.g. Refs. 21-26 for recent molecular simulation studies), but rather focus on larger pores (a few nanometers), where pressuredriven flow may provide a contribution in addition to mere diffusion. Molecular dynamics have already provided detailed information on the structure and dynamics of the interface between such charged layered aluminosilicates and aqueous solutions. Previous studies focussed e.g. on the structure and energetics of ion sorption onto external surfaces, 27,28 on the relevance of continuous double-and triple-layer models to describe the ion distribution, 29 on the structure and energetics of ions in the Stern layer 30 or on the diffusion properties on the molecular (to investigate the solvent and ion dynamics within the first adsorbed layers of fluids [31][32][33][34][35] ) and larger scales (from the local interlayer and interparticle dynamics to that through a porous medium consisting of clay particles 36 ).

Using non-equilibrium molecular dynamics simulations, we have also investigated pressuredriven flow in sodium montorillonite clay nanopores. 37 Due to the negative charge of the surface and the presence of counterions, electro-osmotic flows can also be induced by applied electric field, and the corresponding transport properties have also been investigated with molecular dynamics simulations. [38][39][40][41][42] It was found in particular that the slip length measured for both types of flow are similar. Despite the small value (smaller than the size of a water molecule), neglecting such slip may result in sizeable errors in the prediction of the overall flow in nanopores.

Clay minerals further provide an ideal situation to investigate the influence of the nature of the surface and of the fluid on the hydrodynamic properties of experimentally relevant interfaces. Indeed, they can differ by the presence or absence (and density) of substitutions in the mineral framework, leading to various charge densities, as well as by the localization of these substitutions within the framework (closer to or further from the surface in contact with the fluid). Finally, the nature of the counterion compensating the charge of the mineral may also be varied, resulting in ion-specific properties. Following previous work e.g. on ion sorption 28 or clay swelling, 43 here we investigate the effect of the nature of the clay surface, that of the counterion and the presence of salt on the properties of clay mineral -aqueous interfaces. We focus in particular on hydrodynamic properties and on the possibility to describe the flow in the presence of an external perturbation with continuum hydrodynamics, using appropriate boundary conditions derived from molecular simulations. This requires understanding first the structure of the interfacial fluid. Systems and methods are presented in the next section. We then present equilibrium density and water diffusion profiles for the various systems, and finally discuss the hydrodynamic properties.

Methods

Systems

We consider three different clay minerals: pyrophyllite, montmorillonite and beidellite. All three belong to the family of 2:1 phyllosilicates (layered aluminosilicates), where an octahedral aluminum oxide sheet is sandwiched between two tetrahedral silicon oxide sheets (see Figure 1). While pyrophillite with unit formula Al 4 Si 8 O 20 (OH) 4 is electrically neutral, montmorillonite and beidellite are negatively charged due to the presence of substitutions in the octahedral and tetrahedral sheets, respectively. This charge is then compensated by counterions in the interlayer. Specifically, we simulate montmorillonite and beidellite clays with identical charge densities with unit formula X 0.75 (Al , respectively, with X=Na + or Cs + the compensating counterion. This allows us to investigate the combined effects of the clay charge (or absence thereof), its localization within the mineral layer and the nature of the counterion.

Computational Details

For each system, the simulation box (see Figure 1) contains one clay layer of 4x8 unit cells.

The layer is divided in two halves and periodic boundary conditions in all directions are used.

A single layer may seem insufficient to describe the interface between a clay particle (with several layers and interlayers) and a solution. Indeed, it is known that for low interlayer content the properties of bulk clay particles may depend on the number of layers used, even though the influence is already small in the bihydrated states (two water layers). 44 We have checked for one of the systems that no differences were observed when considering: (a) a single clay layer, without interlayer (b) two layers and an interlayer in the bihydrated state and (c) three layers and two interlayers in the bihydrated state. This suggests that the first setup, which was then chosen for its lower computational cost, is sufficient to mimic the interface between a fluid and a clay particle in the bihydrated state. Such a state if of course relevant in real materials only for swelling clays such as montmorillonite and beidellite. Pyrophyllite doesn't swell in practice, but we keep the same setup for a consistent comparison between the three clays.

The negative charge of the mineral is compensated by 24 Na + or Cs + counterions. In the following, we will refer to the simulated systems as 0 M (resp. 1 M) if they contain no additional ions (resp. 40 NaCl or CsCl ion pairs). The box size is 41.44 × 35.88 × 58.4 Å3 for all systems, with a common length in the direction perpendicular to the surfaces determined from a prior 500 ps simulation for Na-montmorillonite in the N P T ensemble with P = 1 bar and T = 298 K using Nosé-Hoover barostat and thermostat with time constants of 2 ps and 1 ps, respectively. The salt-free systems all contain 2436 water molecules, while their number for the salt-containing systems was adjusted to recover a concentration of 1 M far from the surface (2408 and 2332 molecules with NaCl and CsCl, respectively). This choice of imposing the bulk density in the central region, rather than the water and ions chemical potentials, is motivated by the objective of a simple hydrodynamic description of a fluid characterized by its bulk properties (most importantly for hydrodynamics, its density and viscosity), together with an interface-specific boundary condition, we have preferred to impose. We note that this results in slightly different pressures for a given pore size (to facilitate the comparison between systems). Interactions between atoms are described by a classical, non-polarizable force field including Coulomb and Lennard-Jones (LJ) interactions. The partial charges and LJ parameters are taken from the CLAYFF force field for the clay atoms, 15 from the SPC/E water model 45 and from Ref 46 for the ions, together with the Lorentz-Berthelot mixing rules. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] The water molecules are treated as rigid via the SHAKE algorithm. [START_REF] Ryckaert | [END_REF] Long-range electrostatic interactions are computed using the Particle-Particle-Particle-Mesh method, [START_REF] Hockney | Computer simulation using particles[END_REF] while a cut-off radius of 12 Å is used for the computation of LJ interactions. Molecular dynamics simulations are performed using a time step of 1 fs. A Nose-Hoover 50 thermostat with a time constant of 1 ps is applied on the clay atoms only. Energy exchanges with the fluid are then sufficient to thermalize the latter for both equilibrium and non-equilibrium simulations.

For each simulation box, the system is first equilibrated for 1 ns. The density profiles are then computed from production runs of 5 ns. The final configuration is then used as initial configuration for two sets of simulations. On the one hand, four independent trajectories of 4 ns each (following 2.5 ns of equilibration after generating different sets of random velocities), used to compute diffusion coefficients in various layers (see below) using the method described in Ref. 31. The reported uncertainties correspond to the standard deviation between these 4 independent trajectories. On the other hand, non-equilibrium simulations in which a force is applied to all atoms in the fluid (water oxygen and hydrogen, ions) are performed to sample the steady-state velocity profiles for 6 ns, after a transient regime of 5 ns. This common approach is not strictly equivalent to an applied pressure gradient, which corresponds to a uniform volumic force, but it has the advantage of simplicity as it does not require prior knowledge of the density profiles. In such non-equilibrium simulations, a single Al atom in each clay layer is constrained at its initial position to prevent the drift of the whole solid due to the liquid-solid friction.

Results and discussion

3.1 Equilibrium properties

Structure of the interface

We first investigate the structure of the various interfaces. Fig. 2 shows the density profiles of water oxygen and hydrogen for the three different clays with Na + counterion (similar results are obtained with Cs + ), as well as the density profiles of the cations in the presence of 1 M of added salt for both the Na + and Cs + cases. As reported previously on similar surfaces (see e.g. Ref. 31-33), water is organized with 2-3 well defined layers at the surface, beyond which it behaves as a bulk-like fluid. The combined oxygen and hydrogen positions indicate that the orientation of water in the first adsorbed layer depends on the charge of the surface.

While in pyrophyllite hydrogen bonds are formed almost exclusively among water molecules, as on a hydrophobic surface, 51 a significant fraction of OH bonds are oriented toward the surface on beidellite and montmorillonite. Na + . 53 These observations are consistent with the distribution of these cations in bihydrated interlayers of these clays. 43 We finally note an asymmetry between both sides of the ion density profiles in the case of beidellite, which results from the asymmetry in the distribution of the substitutions within the tetrahedral sheet: Even though the number of substitutions and hence that of counterions is the same on both sides of the layer (12 each), they are more clustered on one side and more evenly distributed in the other (see Figure 3). The repulsion between adsorbed counterions then modifies their distribution in the vicinity of the attractive charge defects within the solid. A difference between NaCl and CsCl salts is also observed in the case of the uncharged pyrophyllite. Specifically, the strongly hydrated Na + cations are somewhat depleted from the first adsorbed water layers. In contrast, an excess of Cs + is seen near the surface. Even though Cs + remains solvated, this is reminiscent of the behaviour observed at the liquidvapour interface of water. 54 This further supports the hydrophobic behaviour of the neutral pyrophyllite surface, which was already confirmed in molecular simulation in particular by the computation of contact angles. 55

Al 3+ Td Na +

Water diffusion

The nature of the ion and of the clay mineral also has an influence on the dynamics of the The results are presented in Figure 4 for both cations in the absence of added salt (note that the system is the same for pyrophyllite). In the bulk region, we obtain the same diffusion coefficient for all surfaces, as expected. The common value of 2.9 10 -9 m 2 s -1 is larger than the experimental value for bulk water (2.3 10 -9 m 2 s -1 ). Two factors may explain such a discrepancy. On the one hand, the SPC/E water model is known to overestimate the bulk diffusion coefficient (the extrapolated value for an infinite cubic box size is close to 3.0 10 -9 m 2 s -1 , see Ref.

∆D (H > L) = k B T η 3 40
H L 2 - 3 ln(1 + √ 2) 4πL (1)
which was shown to be in very good agreement with molecular dynamics simulation results for a Lennard-Jones fluid between walls of Lennard-Jones particles. In the present case, this expression leads to corrections ranging from -1.1 × 10 -10 to -1.7 × 10 -10 m 2 s -1 for all systems. These corrections are smaller than most variations of the diffusion coefficients presented below to discuss the effects of the nature of the surface and counterions (see the results of Figure 4).

In all cases, D decreases as water approaches the surface (the only exception, discussed below, is the first water layer on Na-montmorillonite, which diffuses slightly faster than the second layer). For both counterions, the presence of substitutions in the mineral results in smaller water diffusion coefficients near the surface, possibly due to the change in orientation of molecules due to the surface electric field, 39 but also to the associated presence of counterions in the adsorbed water layers. Indeed, the residence of water in the solvation shell of cations tends to slow down their overall diffusion. Within the first layer, D is reduced more dramatically for both counterions on the surface of beidellite than that of montmorillonite. This can be explained by the localization of the substitutions closer to the fluid in the former case, resulting in a stronger interaction of the cations with the surface. Overall, water diffuses faster with Cs + as a counterion on both beidellite and montmorillonite. This subtle combination of surface-and ion-specific effects can be explained by the larger density of Na + in that layer compared to Cs + as discussed in the previous section (the density profiles of Figure 2 correspond to the case with added salt, but the results are similar without salt). The virtual absence of inner-sphere complexes on Na-montmorillonite then explains the non-monotonic behaviour of D in the vicinity of that surface. Similar observations can be made in the case of added salt (results not shown). One observes however a global decrease of D in that case. In the central regions, the computed diffusion coefficients are 2.5 and 2.7 10 -9 m 2 s -1 for Na + and Cs + , respectively. A larger decrease in the former case is consistent with previous simulations 58 and with experiments [START_REF] Mills | Self-diffusion in Electrolyte Solutions[END_REF] for bulk electrolyte solutions.

The reported results for beidellite are averages over both surfaces. As explained above, the ion density profiles are in fact asymmetric due to the different distribution of the substitutions within the tetrahedral layers (see Figure 3). The difference in cation distributions also results in an asymmetry in the diffusion coefficients in the first (and to a lesser extent second) adsorbed layer. For example, in the Na + case without added salt we obtain (0.76 ± 0.02) and (0.70 ± 0.01) 10 -9 m 2 s -1 in the first adsorbed layer on the left and right surfaces, respectively. The asymmetry is less pronounced with added salt, with (0.72 ± 0.01) vs (0.68 ± 0.02) 10 -9 m 2 s -1 .

To conclude this first part on the equilibrium properties, the present results indicate that the charge of the clay surface and its localization within the mineral layer has no effect on water diffusion away from the first few molecular layers, but it does have a significant impact within the latter. In addition, the nature of the cation has an influence on the diffusion of water within the first adsorbed layer. These features can be explained by the mineral-and ion-specific structure of the interface.

3.2 Hydrodynamics: slip, stick and stagnation

Continuous hydrodynamics

We now turn to the hydrodynamic properties of these interfaces. Following previous studies, we apply a body force on water molecules and ions and measure the steady-state velocity profile. At the continuum level, the balance between viscous and other forces is then described by Stokes equation η∆v + f V = 0, where η is the fluid viscosity, v is the velocity and f V the volumic force acting on the fluid. In the case of the Poiseuille flow (uniform

f V = f V e x
in a slit pore), the profile v = v(z)e x is parabolic. Defining z = 0 as the mid-plane of the pore and ± L hyd 2

the position of the hydrodynamic interfaces, v(z) reads:

v(z) = v L hyd 2 - f V 2η z 2 - L hyd 2 2 (2)
Therefore the viscosity η can be computed from the curvature of the profile and the velocity at the interface v(L hyd /2) can be used to determine appropriate boundary conditions. When the velocity of the fluid does not vanish at the interface (slip boundary condition), a slip length b can be defined as:

5 1 b = - 1 v(L hyd /2) dv dz z=L hyd /2 . ( 3
)
It can therfore be determined from the molecular simulation velocity profile when L hyd is known. Several choices are possible. In particular, Bocquet and Barrat provided a Green-Kubo expression for the location of the "hydrodynamic position". 11 Following previous work, 37,39 we consider here a more phenomenolgical approach, based on the partitioning of the pore space between a region of homogeneous fluid and a region without fluid. Therefore we identify in this slip case the position of the hydrodynamic interfaces ±L hyd /2 with the positions ±z GDS of the Gibbs Dividing Surfaces (GDS):

±z GDS 0 [ρ bulk -ρ(z)]dz = +∞ ±z GDS ρ(z)dz (4)
where ρ bulk and ρ(z) are respectively the fluid density in the bulk, i.e. in the central region of the pore, and that at z. Thus z GDS can be calculated from the density profile of the fluid ρ(z) obtained by simulation (see Figure 5). In contrast, when the fluid sticks to the surfaces, a layer of thickness d is immobile. The positions of the shear planes can be determined from the MD velocity profiles as v(L hyd /2) = 0. In this stick case they are located further from the surface than the GDS. The distance between the shear plane and the GDS then defines a stagnation length d as:

d = z GDS - L hyd 2 (5)
The location of the GDS and the hydrodynamic width for all systems will be reported below, together with the corresponding slip or stagnation lengths.

As can be seen on Figure 5, the parabolic fit in the central region is a very good approximation, even if a common value of the viscosity is used for similar systems (see next section). It should be noted, however, that deviations from the parabolic profile may be observed near the wall where the fluid is layered. The resulting estimate of the flux near the wall by the parabolic profile is therefore inaccurate, but this contribution to the overall flux is small -even more so that the pore is larger. 37 It is more important to accurately describe the profile further from the surface by introducing an appropriate boundary condition.

Viscosity of the confined fluid

A force of 0.75 cal.mol -1 Å-1 per atom is applied in the x direction. This value was chosen in the range where the response of the system is linear. The steady-state velocity profiles for the three clays and both cations, are reported in Figure 6 without and with added salt in panels 6a and 6b, respectively. The nature of the clay, the nature of the counterion and the Figure 5: The velocity profile obtained by non-equilibrium molecular dynamics (solid blue line) here for montmorillonite with Na + as counterion without added salt, is fitted to a parabola in the central region (green dashed-dotted line) and the viscosity is determined from the curvature. The figure also shows the fit to a parabola enforcing a value of the viscosity common for similar systems (see text, red dashed line), which is then used to determine the boundary conditions. The Gibbs Dividing Surfaces (GDS, vertical black dashed line, defined by Eq. 4), determined from the density profile (grey dotted line) is used to localize the solid/fluid interface.

salt concentration all contribute to the hydrodynamic response. Nevertheless, in all cases the velocity profile is parabolic in the central region, so that the viscosity can be computed from the curvature (in practice, the fit is done in a central slab of width 3 nm). The corresponding values are summarized in Table 1. The reported error bars on the viscosity were obtained by comparing the estimates for the first and second half of the production runs.

Despite some variations, the results for a given family of systems are similar. In the saltfree case, regardless of the nature of counterion, the values are consistent with our previous work on pyrophyllite and Na-montmorillonite. 37,39 In the following analysis of the velocity profile, we therefore use the average value in all salt-free systems which is in perfect agreement with the bulk value (η = 0.68 cP) for the SPC/E water model at room temperature. 56 Similarly, the viscosities for a given salt are consistent for all clays and in the following we use a common value of 0.89 cP (resp. 0.78 cP) for all systems with 1 M of NaCl (resp Figure 6: Steady-state velocity profiles under an applied force parallel to the surface applied each fluid atom (the same force is applied in all cases), in the absence of added salt (a) and with 1 M of salt (b). Dashed and solid lines correspond to Na + and Cs + counterions, respectively, while colors refer to the clay type: green for pyrophyllite, blue for montmorillonite and red for beidellite. Note that in (a) no cations are present for pyrophyllite, so that symbols are used instead of lines.

CsCl

). An increase in the viscosity with NaCl concentration is consistent with experimental results. 60 In contrast, the less pronounced increase in the CsCl case points to a limitation of the force field, since the experimental viscosity decreases in the presence of CsCl salt. 61 Nevertheless, the smaller viscosity with Cs + compared to Na + is in the same order as in the experiments. We finally note that the choice of common values to model all velocity profiles of systems without salt, with 1 M of NaCl or CsCl, even though it does not change the conclusions of the discussion below, may have a quantitative influence on the resulting slip and stagnation length and therefore contributes to the reported uncertainties.

Hydrodynamic boundary conditions

From the velocity profiles, and under the above assumptions on where to apply the boundary condition and on the value of the viscosity, we determine in each case the slip length b or the stagnation length d. The results are summarized in Tables 2 and3 for Na + and Cs + , respectively, together with the corresponding location of the GDS and hydrodynamic width.

Reported uncertainties are based on comparing estimates for two half trajectories as well as Montmorillonite with the CLAYFF force field, we had determined a slip boundary condition with a slip length of b = 1.4 ± 0.3 Å. Here we find instead a small stagnation length d = 0.5 ± 0.3 Å. While such a difference is not huge, it is statistically significant and may be due to the fact that in the former study the clay layers were treated as rigid, while in the present work we account for their flexibility and their finite temperature. It is therefore not surprising that the resulting energy dissipation at the wall results in less slippage at the interface. However the effect of wall flexibility on the flow is not straightforward, as recent work on carbon nanotubes suggests e.g. that it could reduce the friction at the interface, 62 and coupling of the fluid dynamics with phonons in such tubes may significantly enhance the diffusion of water. 63,64 We also note that the difference between neutral and charged surfaces is much larger than between the two charged clays. Nevertheless, the small stagnation length in the case of beidellite (compared to b ≈ 0 or d ≈ 0 for montmorillonite) reflects the stronger interaction of cations with the substitutions which are closer to the fluid in that case, as explained above.

The effect of surface charge on hydrodynamic slippage was analyzed quantitatively by Effect of salt concentration and ion specificity In the absence of added salt, the results are very similar for both counterions. There is little effect of adding NaCl salt, except for a slight decrease in the slip length for pyrophyllite. Adding CsCl slightly increases the stagnation length with montmorillonite and beidellite and slightly decreases the slip length with pyrophyllite. Overall, the effects of added salt and of the nature of the counterions are less pronounced than that of the surface charge and its localization.

As a final remark, we note that the above discussion of hydrodynamic boundary conditions is based on a parabolic fit of the flow profile in the central region of the pore. While this is the most relevant from a practical point of view in order to introduce mineral-and ionspecificity into hydrodynamic descriptions for larger scales, it should be reminded that the molecular velocity profile deviates from this parabolic fit in the close vicinity of the surface. This is obvious in Figure 5, where the velocity in the first adsorbed layer is larger than the profile extrapolated from the bulk. In this region, the shape of the molecular velocity profile may also be different from what would be inferred from the bulk viscosity. While introducing a local viscosity to account for such variations is not sufficient to obtain a quantitative description, 37 an illustration of such an effect can be found in Figure 6a with the case of montmorillonite. In this salt-free case, the velocity increases more slowly close to the surface with Na + as a counterion than with Cs + , even though the fluid far from the surface is the same in both cases (hence identical viscosities). Such a difference probably arises because fully hydrated Na + near the surface (outer-sphere complexes) have a stronger effect on the interfacial fluid dynamics than Cs + ions which form inner-sphere complexes.

Conclusion

The structure, diffusion and hydrodynamic properties of clay interfaces with aqueous solution depends on the nature of the clay, the nature of the counterions and the salt concentration in the solution. We used molecular dynamics to investigate water-filled nanopores between uncharged (pyrophyllite) and charged (montmorillonite and beidellite, with susbtitutions located in the octahedral and tetrahedral layers, respectively) clays, with sodium or cesium as counterions, in the absence and in the presence of added salt. We discussed how the balance between solvation and attraction of the cations to the surface results in various distributions between inner-and outer-sphere complexes and how this influences the dynamics of water near the surface, as well as the hydrodynamic flow in the presence of an external force.

In the latter case, the discussion based on a mapping of the molecular velocity profiles to a continuous description (parabolic Poiseuille flow) shows that the larger effects come from the presence/absence of charge in the mineral, as well as the localization of substitutions within the clay layer. The salt concentration and the nature of the counterions have a comparatively less important impact far from the surface -even though some differences are observed in its close vicinity, which are not properly captured by the continuous description.

In the future, this study could benefit from recent developments of improved force fields which allow to better capture ion-and mineral-specificity in clays and related aluminosilicate minerals, [66][67][68] starting with the structure and diffusion under equlibrium conditions. Such models may also improve the description of the structural relaxation of the mineral itself, induced by the truncation of the mineral periodicity at the interface with the liquid, demonstrated e.g. in X-ray reflectivity experiments on mica. 69 This work could also be extended using force fields which better describe concentrated electrolyte solutions, e.g. via rescaled charges 70 or introducing polarization effects, 71 as well as by considering other ions, including divalent ones, and of course other simple or more complex including binary fluids. [72][73][74][75] As a first step, one could for example examine mixtures of sodium and cesium counterions.

Finally, in the broader perspective of transport on the macroscopic scale, the results of the present work could be introduced in a bottom-up approach combining molecular simulation results with transport models on larger scales. [76][77][78]

Figure 1 :

 1 Figure 1: Top Pyrophyllite, montmorillonite and beidellite clay layers differ by the presence or absence of substitutions in the octahedral (Al 3+ →Mg 2+ , montmorillonite) and tetrahedral (Al 3+ →Mg 2+ , beidellite) sheets. Bottom Simulation box (see text for details).

Figure 2 :

 2 Figure2: Distribution of (a) water oxygen, (b) water hydrogen, (c) Na + and (d) Cs + as a function of the position z across the pore, for a 1 M salt concentration and the three clays: pyrophyllite (green), montmorillonite (blue) and beidellite (red). The water densities are similar for both counterions and only the case of Na + is shown.

Figure 3 :

 3 Figure3: The asymmetric ionic density profile for beidellite (here for Na + counterions with 1 M NaCl) arises from the asymmetry in the substitutions within the tetrahedral layer close to the surface. The bottom panels indicate the localization of the Al susbtitutions (green) within the surface Si and O atoms (grey). The positions of Na + ions (red) within the first adsorbed layers, for 400 configurations, is also reported. Clustering of substitutions (left side) results in a decrease in the ionic density near the surface due to ion-ion repulsion.

  adsorbed fluid. Specifically, we discuss here the component D of the diffusion tensor of water along the clay surfaces. Since it depends on the distance of the water molecules from the surface, we first define relevant regions in which D is computed, based on the density profiles.31 These regions, illustrated by vertical lines in Figure4, are defined according to the position z of the water oxygen atoms: within 4 Å of the plane containing the surface oxygen atoms, between 4 and 7 Å, between 7 and 10 Å; the central region has a width of 30.8 Å. D is then calculated in each layer as explained in Ref.31.

Figure 4 :

 4 Figure 4: Water diffusion coefficient along the surface, D , in the regions defined by the dashed vertical lines, for the three different clays (green for pyrophyllite, blue for montmorillonite and red for beidellite) and the two cations (Na + left, Cs + right), in the absence of added salt. The corresponding total density profiles are also shown in the bottom panel to clarify the definition of the various regions, in the case of montmorillonite. Error bars are smaller than the symbols.

  56). On the other hand, the diffusion coefficient computed by molecular dynamics depends on the size and the shape of the simulation box, due to hydrodynamic finite-size effects. We derived in Ref.57 an analytical expression for the correction ∆D to the average diffusion coefficient along planar walls separated by a distance H, with stick boundary conditions, due to the finite length L of the simulation box in the directions along the walls. Even though in the present case the box lengths in the x and y directions are not equal, this is approximately the case. When H > L (elongated systems), the correction reads:

Joly et al. 65

 65 More precisely, they predicted and confirmed by MD simulations for a Lennard-Jones (LJ) fluid with LJ walls that the slip length varies b= b LJ /[1+(1/α)(Σσ 2 /e) 2 (l B /σ)b LJ ]with b LJ the slip length for a neutral surface, α ≈ 1 a constant, Σ the surface charge density, σ the LJ diameter of the wall atoms and l B the Bjerrum length. This simple expression correctly predicts the decrease in b from neutral pyrophyllite to charged montmorillonite and beidellite, even though it is not sufficient to quantatively describe the subtle difference between the latter two.

  3.25 Mg 0.75 )Si 8 O 20 (OH) 4 and X 0.75 Al 4 (Si 7.25 Al 0.75 )O 20 (OH)

  4

Table 1 :

 1 Fluid viscosities for the three clays with Na + or Cs + as counterion and two salt concentrations.

	Clay	Salt	Viscosity (cP)
	Pyro	-	0.62 ± 0.03
	Pyro	NaCl 1M	0.95 ± 0.02
	Pyro	CsCl 1M	0.74 ± 0.02
	Mont-Na	-	0.70 ± 0.02
	Mont-Cs	-	0.71 ± 0.02
	Mont-Na NaCl 1M	0.84 ± 0.02
	Mont-Cs CsCl 1M	0.79 ± 0.06
	Beid-Na	-	0.74 ± 0.02
	Beid-Cs	-	0.68 ± 0.03
	Beid-Na NaCl 1M	0.88 ± 0.02
	Beid-Cs CsCl 1M	0.80 ± 0.04
	the choice of viscosity (see previous section).

Table 2 :

 2 Position of the Gibbs dividing surface z GDS , Slip length b or stagnation length d of the fluid in Å, for the three Na-clays and two NaCl concentrations. The boundary condition is determined from the velocity profile using a viscosity of 0.68 cP for without added salt and 0.89 cP with 1 M NaCl.

	Clay	NaCl salt z GDS ( Å) L hyd /2 ( Å)	b ( Å)	d ( Å)
	Pyro	0 M	24.2 ± 0.2 24.2 ± 0.2 5.9 ± 0.4	-
	Pyro	1 M	24.1 ± 0.2 24.1 ± 0.2 5.6 ± 0.3	-
	Mont-Na	0 M	24.8 ± 0.2 24.3 ± 0.5	-	0.5 ± 0.3
	Mont-Na	1 M	24.9 ± 0.2 24.9 ± 0.2 0.2 ± 0.5	-
	Beid-Na	0 M	24.8 ± 0.2 23.6 ± 0.5	-	1.2 ± 0.3
	Beid-Na	1 M	24.6 ± 0.2 23.7 ± 0.6	-	0.9 ± 0.5
	Effect of surface charge and its localization With both cations, with and without added
	salt, we observe a clear pattern of decreasing slip from the neutral pyrophyllite to the charged
	montmorillonite (with octahedral substitutions) and beidellite (with tetrahedral substitu-
	tions). More precisely, a slip boundary condition with a slip length of approximately two
	water layers applies to pyrophyllite; a stick boundary condition applies to montmorillonite
	with b ≈ 0 or d ≈ 0 (except with 1 M CsCl); a small stagnation layer of approximately
	one half of water diameter should be introduced for beidellite. In a previous study of Na-

Table 3 :

 3 Slip length b or stagnation length d of the fluid in Å, for the three Cs-clays and two CsCl concentrations. The boundary condition is determined from the velocity profile using a viscosity of 0.68 cP for without added salt and 0.78 cP with 1 M CsCl.

	Clay	CsCl salt z GDS ( Å) L hyd /2 ( Å)	b ( Å)	d ( Å)
	Pyro	0 M	24.3 ± 0.2 24.3 ± 0.2 5.9 ± 0.4	-
	Pyro	1 M	24.1 ± 0.2 24.1 ± 0.2 5.2 ± 0.4	-
	Mont-Cs	0 M	25.3 ± 0.2 25.0 ± 0.5	-	0.3 ± 0.3
	Mont-Cs	1 M	25.3 ± 0.2 24.2 ± 0.6	-	1.1 ± 0.4
	Beid-Cs	0 M	25.5 ± 0.2 24.2 ± 0.6	-	1.3 ± 0.4
	Beid-Cs	1 M	25.5 ± 0.2 23.8 ± 0.6	-	1.7 ± 0.4

Acknowledgement

The authors acknowledge financial support from IFPEN and access to computing resources on Curie (TGCC, French National HPC) via the GENCI project x2016087684 and on MeSU (UPMC).

Pyrophyllite