Andi Cuko

Antoni Maciá

Monica Calatayud
email: calatayu@lct.jussieu.fr

Stefan T Bromley
email: s.bromley@ub.edu

Antoni Macià Escatllar

Properties of hydrated TiO2 and SiO2 nanoclusters: dependence on size, temperature and water vapour pressure

Introduction

Nanoparticles of titania (TiO2) and silica (SiO2) are two of the most massively produced nanomaterials in the world. 1 Nano-TiO2 and nano-SiO2 are both widely used as fillers to improve the properties of materials, whereby the former typically imparts enhanced (photo)electroactive properties and the latter increased structural resistance and visible light transmittance (e.g. proton conducting electrolytes, 2 photo-resistive rubbers, 3 pollutant-degrading durable mortars 4 ). This synergy has also led to a number of mixed-oxide titanosilicate nanomaterials with applications in (photo)catalysis, 5 water remediation, 6 and self-cleaning. 7 In all the above applications, the nano-TiO2 and nano-SiO2 components typically, and often deliberately, come into contact with watercontaining environments. Understanding and controlling such interactions is thus essential for optimising the performance of these systems. The aqueous chemistry of nano-SiO2 and nano-TiO2 is also important for understanding the bottom-up solution syntheses of nanoporous silicates 8 and TiO2-based nanomaterials 9 and the possible eco-toxicity of these species. 10 Many studies have specifically focussed on the interaction of water with nano-TiO2 due to the potential these systems have for fuel production through photo-dissociative water splitting. 11 A number of recent theoretical studies have used top-down bulk cut nanoparticles containing 10s -1000s of atoms in order to investigate the electronic properties of nano-TiO2 photocatalysts 12,13,14,15 and their interactions with water. 16,17,18,19,20,21 Although such calculations can provide important new insights, in most cases the nanoparticles used are still relatively small idealised model systems compared to experimental systems where larger and more complex nanoparticles are typically used. Detailed spectroscopic experiments probing the interaction of water with small oxide (i.e. ~10-30 atoms) systems have tended to focus on hydroxylated silica (e.g. 29 Si nuclear magnetic resonance spectroscopy 22,23,24 and mass spectroscopy 25 ). As such there are a number of theoretical studies proposing corresponding structures for small (SiO2)M(H2O)N species with M=2-4 26,27,28 and their oligomerisation reactions. 27,29,30 Recently, cluster beam experiments by Weichman et al. 31 reported infrared (IR) photodissociation spectra for gas phase [(TiO2)M(D2O)N] -for M=2-4 and N=1-3. In their work, direct comparison of the experimental IR spectra with theoretically calculated IR spectra from low energy cluster isomers was used to assign structures to the experimental clusters. The assignments utilised cluster structures from a previously reported global optimisation search. 32 The use of global optimisation provides a means to find stable structures of small nanoclusters which are generally unrelated to simple cuts from the respective bulk crystal. The global optimized structures show a huge range of variability with respect to size and chemical composition, yet, for small sizes they are more stable than crystal cuts of the same size. 33,34 Herein, following our previous theoretical work on anhydrous SiO2, 35,36,37 TiO2 33 , and TiO2-SiO2 38 nanoclusters and hydroxylated SiO2 39,40,41 nanoclusters, we use global optimisation to provide detailed bottom-up insights into the structures and stabilities of (TiO2)M(H2O)N nanoclusters with sizes from M = 4 -16 and N/M ratios of ≥ 0.5 (i.e. up to 72 atoms). We compare the obtained structures with the corresponding global minima (SiO2)M(H2O)N nanoclusters. Further, we calculate how the free energies of both (TiO2)M(H2O)N and (SiO2)M(H2O)N nanoclusters depend on size, degree of hydroxylation, temperature and water vapour pressure. In this way we provide an unprecedented detailed comparison between nano-TiO2 and nano-SiO2 in terms of how they interact with water, while bridging the gap between very small molecular species and bulk-like nanoparticles for these important hydrated nano-oxide systems.

Methodology

In this study we apply the Monte Carlo Basin Hopping (MCBH) algorithm 42 for exploring the potential energy surface (PES) of hydrated (TiO2)M(H2O)N nanoclusters. This global optimisation methodology consists of repeated steps whereby in each step random structural distortions are applied to a cluster followed by optimisation of the cluster structure to a local energy minimum on the PES. The Metropolis criterion 43 is applied after each step to decide whether or not to accept the new minimum energy structure. This is a probabilistic criterion which depends on the simulation temperature and the energy difference between the present cluster structure and the previously accepted one. An accepted cluster is used as initial structure in the following step. In each step, we use interatomic potentials (IPs) for the structural optimisation. The IPs we use are empirical functions that provide a computationally efficient means to evaluate the relationship between structure and energy. Although the IP-MCBH approach is a very effective means to quickly assess many low energy minima on the PES, IPs generally provide only a reasonable estimate of energetic stability of any cluster structure. In order to refine the accuracy of our predictions we selected a set of low energy structures from each global optimisation run which were then fully optimised using density functional theory (DFT) based calculations (see details below). We have successfully applied this method in several previous studies, often finding new improved global minima candidates with respect to other reported cluster isomers. 35,36,37 ,38,39,40

Global optimization of (TiO2)M(H2O)N systems

The choice of the IP plays a decisive role in determining the accuracy of predicting the structure of low energy clusters. For titania-water systems, although there are a number of empirical IPs available in literature, 44,45,46,47 such IPs are parameterized for bulk titania and were found to be inappropriate for finding low energy nanoclusters which have non-bulk-like structures and in which the coordination of Ti cations ranges from 4 to 6. Therefore, in order to have an IP suitable for global optimisation of (TiO2)M(H2O)N nanoclusters we constructed a new IP. Due to the structural similarities between globally optimized nanoclusters of bare silica and titania, 38 we based our IP on a well-established silica IP for hydroxylated silica systems (FFSiOH). 48 In particular, we have previously shown that FFSiOH is very effective in accurately predicting low energy structures and stabilities of hydroxylated silica nanoclusters. 49,50 In order to provide relevant data to fit our new IP we first selected global minima candidates of anhydrous (TiO2)M from our previous work 33 with sizes between M = 4 -14. All the clusters were then manually hydrated by adding molecular water in a dissociative way for low degrees of hydration and via molecularly adsorbed water for highly hydrated clusters. For each cluster we then ran DFT-based ab initio molecular dynamics (AIMD) for 3 ps at 300 K employing the PBE functional 51 with a light/tier-1 basis of numerical atom-centred basis functions employing the FHI-AIMS code. 52 In this way we produced a large and varied training set of structures for fitting the parameters of our IP. By adapting the original FFSiOH parameters to reproduce the radial distribution functions of the AIMD runs we produced a new parameterisation which we refer to as FFTiOH. Using FFTiOH in the MCBH runs, we found that it was often difficult to optimise randomly distorted geometries when large core-shell separations occurred. To overcome this issue, we generated a simpler version of FFTiOH where core-shell polarizable ions were removed and the O-H interaction was fitted to a Buckingham potential to mimic as much as possible the original Morse potential. We We used the MCBH runs to globally optimise (TiO2)M(H2O)N clusters with M = 4, 8, 12 and 16 and for each M respectively N = 1 -4, 1 -5, 1 -6 and 1 -8. A typical MCBH run employed either FFTiOH* or FFTiOH*-mod for 20000-100000 steps depending on the system size. Our MCBH algorithm uses the Atomic Simulation Environment (ASE) 53 for manipulating structures and the GULP 54 code for structural optimisations. During the MCBH procedure the temperature was dynamically adapted to maintain the ratio of accepted steps close to 65%. The MCBH procedure was repeated a number of times for each distinct cluster (TiO2)M(H2O)N stoichiometry using both FFTiOH* and FFTiOH*-mod to better ensure a thorough search of the respective PES. Our MCBH searches were also further extended by data mining from hydrated silica nanoclusters from our previous works. 39,40,41 From each MCBH run a selection of isomers were optimised using progressively more accurate methods. Firstly, the 200-500 lowest energy cluster isomers optimised using FFTiOH* and FFTiOH*-mod during each MCBH run were first post-optimised optimized with FFTiOH. Afterwards, between 50 and 150 (depending on the system size) of the energetically lowest lying structures were fully optimized with DFT based calculations using the PBE functional and a tier1/light numeric atomcentred orbital basis set. Then, all cluster isomers within 0.5 eV of the energetically lowest PBEoptimised structure were further optimized using the PBE0 55 functional and a tier1/tight basis. These calculations were also used to extract the difference in energy between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) in order to estimate the vertical electronic energy gap magnitudes. All DFT based calculations were performed with FHI-AIMS package. 52 We note that the employed basis set approximately provide results of a similar or higher quality to those obtained with valence triple-zeta plus polarization Gaussian type orbitals. 56 Generally, the higher the hydration degree and/or the system size the higher the number of structures were optimized using DFT because of the correspondingly higher density of distinct cluster structures with respect to energy.

Energetics of hydration

As discussed above, both nano-TiO2 and nano-SiO2 are prone to incorporate water which we can formally write as

( ) + ( ) ⇄ ( ) ( ) (eq 2.1)
Where is the cation (i.e. Si, Ti), M and N are the number of formal oxide units and the number of in-cluster incorporated water molecules, respectively. The first measure of stability of water incorporation can be assessed by calculating the energy change for each hydration step, which is known to be system size dependent. Using the lowest energy (XO2)M(H2O)N structures for each degree of hydration (RN/M = N/M) for any single system size (M) we can calculate the energy change upon hydration with relative to the respective anhydrous bare (XO2)M cluster via

= ( ) ( ) -( ) - (eq 2.2)
where ( ) ( ) and ( ) are the energies of hydrated and anhydrous nanoclusters respectively, and is the energy of the water molecule. provides the energetic strength of the interactions between water and silica/titania for each nanocluster size, M, considered. A negative indicates thermodynamic favorability of cluster hydration and vice versa.

For a more accurate assessment of the thermodynamics of hydration, free energies must be considered. Following a similar approach to that previously applied to the oxidation of small oxide clusters, 57 using standard statistical thermodynamics we can calculate the free energies of hydroxylation of our titania and silica nanoclusters entirely from the relevant partition function. The partition function for N non-interacting indistinguishable particles is defined as:

( , , ) = [ ( , )] ! (eq 2.3)
where #( , ) is the molecular partition function and N is the number of particles in the system.

#( , ) = # $%&'( ( , ) # %)$ ( )# *+, ( )# ( -( )# ./.0 ( ) # '10 (eq 2.4)

From the partition function, we can calculate the free energy as:

;( ) %)$ = -3 4 56 E8A @ 2A3 4 ℎ D G H + 1 2 3 4 567J K J J L 8 ;( ) *+, = M ℎN + 2 + + M 3 4 56 O1 -PQ9 @- ℎN + 3 4 DR + ;( ) ( -= 3 4 56(S) (eq 2.5)
where B is the system total mass, J K, ,L are the system inertial moments along the x, y and z axes, N + are the harmonic vibrational frequencies and S is the symmetry number. Frequency rotational contributions are calculated within the harmonic and rigid rotor approximations respectively. All systems are treated as ideal gases. Symmetry contributions should formally be considered in order to avoid over-counting the number of possible states. However, for the hydrated systems the presence of hydroxyl groups usually breaks any symmetry and even for the highest symmetric anhydrous systems effect of symmetry has an almost negligible effect on the cluster free energies.

Considering nanoparticles to be in equilibrium with water vapour, the actual degree of hydration of a system will depend on the experimental temperature and water vapour pressure conditions. In a gas-phase hydration reaction, the free energy of the process is described as:

2 T ( , U) = 2 ( ) ( ) ] ( ) -2 ( ) ( ) -V ( , U) (eq 2.6)
where V ( , U) is the water chemical potential in equilibrium with the clusters dependent on temperature T and water vapour partial pressure P. Both clusters free energies and the water chemical potential can be computed from ab initio data by means of statistical thermodynamics as shown in equation 2.5. The water vapour partial pressure dependency arises from -3 4 56( ) term which can be written, according to the ideal gas law, as -3 4 56(3 4 /U). The same term, in the cluster case can be omitted because the reference volume cancels out since we consider differences in free energy. All quantities were calculated using DFT based calculations using the PBE0 functional and using tier1/tight numerical basis set. Harmonic frequencies were calculated by the finite difference method using a displacement of 0.0025 Å for each degree of freedom.

Results and Discussion

In the following, we provide a comparison of structures and thermodynamic stability of known and newly found globally optimized hydrated titania and hydrated silica nanoclusters. found to be very slightly more stable (0.07 eV in total energy) with the method employed herein.

Hydrated nanocluster structures

Both crystalline bulk SiO2 and TiO2 have very low solubility in water but both can dissociatively adsorb water on their surfaces. Amorphous phases and nanoclusters absorb dissociated water more readily due to the presence of reactive surface defects such as terminal non-bridging oxygens (NBOs)

or strained rings of different sizes. 58,59,60,61 Hydration of bare (SiO2)M nanoparticles proceeds by saturating terminal oxygen defective sites first, as they are generally the most reactively unstable defects. 62 Depending on system size, one (for the large particle sizes M=12 and 16) or two (for the small particle sizes M=4, M=8) water molecules are enough to fully saturate such surface defects on the silica nanoclusters. Further incorporation of water molecules leads to hydroxylation of silicon atoms until each silicon atom has one hydroxyl group forming cage-like nanoclusters. Beyond this level of hydration geminal silanol groups begin to be formed (i.e. two hydroxyl groups attached to a single silicon centre) which tends to open up the cage-like structure. The limit of this behaviour is the complete dissolution of each silicon atom to generate silicic acid (Si(OH)4) monomers. We note that this tendency is size dependent whereby silica nanoclusters with larger sizes tend to form denser Hydration of titania nanoclusters proceeds in a similar way to silica nanoclusters at low hydration degree and small nanocluster size. As for silica nanoclusters, depending on the system size, one or two incorporated water molecules are enough to hydroxylate the terminal oxygen defects on anhydrous titania nanoclusters. With further hydration, titania also generates cage-like structures, in some cases even identical to silica (e.g. for M=4, N=2 and M=8, N=4). For relatively higher degrees of hydration titania and silica nanocluster systems are no longer so similar. Unlike Si in silica, the Ti cations in titania readily assume a higher coordination than four (i.e. Ti can be 4-, 5and 6-coordinated) and thus titania can more easily retain its -Ti-O-bonding network with increasing hydration. Indeed, highly hydroxylated titania systems, especially those with M=16, become more compact with 6-coordinated Ti atoms in their centres. With respect to silica, titania nanoclusters can also more easily adsorb molecular water without breaking their bonded -Ti-Onetworks. This capacity is linked to the character of d orbitals in titanium centres which can participate in dative bonds with lone pairs of water molecules.

Despite of these differences some hydrated titania and silica nanocluster structures are found to be very similar. The structures of the (XO2)4(H2O)2, (XO2)4 (H2O)3, (XO2)8(H2O)2, (XO2)8(H2O)4, species ,for example, are essentially the same for both systems. Larger sized nanoclusters with different degrees of hydration also are found to have a similar structures, such as (e.g. (TiO2)12(H2O)4

and (SiO2)12(H2O)5, (TiO2)16(H2O)2 and (SiO2)16(H2O) ).

Hydration free energies

The hydration free energy (∆Ghyd) indicates the thermodynamic tendency of an anhydrous nanocluster to be hydrated at specific temperature and pressure conditions (eq. 2.6). In Fig 5 we plot the ∆Ghyd of silica (solid lines) and titania (dashed lines) with respect to incorporated water ratio (RN/M) at several temperatures (0, 300, 600, 900 K). All hydration energies are normalized by number of oxide units (i.e. M). In the 0 K case each term of eq. 2.6 is approximated by only DFT energies and to distinguish this quantity from free energies we indicate it as ∆Ehyd. In all other cases, ∆Ghyd is calculated under the assumption that the system is in equilibrium with water vapour at a partial pressure of 1 Pa. at different temperatures while keeping the water partial pressure at 1 Pa which is sufficiently low to consider water as an ideal gas even at ambient conditions. Values of total energy, hydration energy at 0K and hydration free energy at 300 K for all systems reported in these graphs are reported in Table S1 of the Supporting Information.

As we can see from Fig. 4, the hydration reaction is a more exothermic process on silica nanoclusters than on the correspondingly sized titania nanoclusters. As the system size (M) increases, the magnitude of the hydration energies decreases. At 0K, for both titania and silica systems irrespectively of the size, the hydration energy decreases indicating a higher stability of hydrated systems with respect to anhydrous ones. The energetic gain through hydration for both systems can be roughly divided into two regions. At first, the slopes are rather steep, and the energy drop per incorporated water molecule is rather constant, although of different magnitude, in both cases. For silica we find a decrease of -0.4 eV for each 0.1% increase of water content, and a decrease of -0.2 eV per 0.1% increase of water content for TiO2. Thus in all cases hydration of silica nanoclusters is more energetically favourable than for respectively sized titania nanoclusters. At a certain degree of

M=4 M=8 M=16 M=12 a) b) c) d)
hydroxylation, depending on the size and system type, the magnitude of the energy change due to hydration begins to reduce and the hydration energies appear to saturate towards a constant value.

In this saturation region, the difference in hydration energy between similar sized titania and silica systems stabilizes to a constant maximum value. This difference varies from 1 eV for system sizes of M = 4 to 0.6 eV for sizes M = 16. When including thermal contributions to the hydration free energy, the magnitude of the energy changes due to incremental hydration are reduced compared to the 0K case. Although, as in the 0K case, the ∆Ghyd data initially decreases with increasing water incorporation, when reaching hydration level marking the beginning of the 0K saturation regime, the energies pass through minimum point beyond where ∆Ghyd starts to increase again. Both these effects are enhanced with increasing temperature. Indeed, at 900K for most of titania systems the ∆Ghyd is nearly zero or positive showing that hydration is thermodynamically unfavourable. This is caused mostly by the vibrational contributions to the entropy and enthalpy of highly hydrated systems. The minimum in ∆Ghyd tends to move towards smaller values of RN/M as the system sizes increase. This indicates that as the system becomes larger, its tendency to be hydrated decreases, or in other words, it tends to be more hydrophobic. Moreover, it is interesting to notice that the silica and titania ∆Ghyd minima for small sizes (i.e. M=4, 8) corresponds to the same water incorporation ratio, RN/M, but, for larger sizes (i.e. M=12, M=16), RN/M is shifted towards smaller values for titania nanoclusters relative to silica. This indicates that, already at the nanoscale, titania tends to be more hydrophobic than silica.

To best of our knowledge, the work of Dixon et al. 32 is the only other study attempting to determine global minima of small hydrated titania nanoparticles. Here the authors reported several low lying structures for (TiO2)M(H2O)N systems with M = 1 -4 and N = 1 -2 and few with M = 8 and N = 3, 4, 7 and 8. In this section we compare their reported structures for M = 4, and some of the predictions and trends they reported for their systems, with our findings. For systems with M = 8

and N = 3 and 4 we obtain the same set of global minima structures, hence we will limit our comparison to systems with M = 4. For systems with M = 4 and N = 1 and 2 we obtain the same lowest energy structures as those reported from Dixon et al. which further enhance the likelihood that such systems are global minima. For M = 4 and N = 3 and 4 we find isomers not reported in their work whose energetic stability is found to be essentially degenerate with their systems (i.e.

differences of only 0.0076 eV for N=3 and in 0.029 eV for N=4 in total energy in favour of our structures). In Fig. 5 we report both structures found by our search as compared with those found in ref [32]. In our M = 4 and N = 3 nanoclusters, the water molecule is incorporated into the nanocluster in a dissociative way, opposed to structures reported in ref 32 where the water adsorbs in a molecular way via a Ti-O(H2) Lewis acid-base complex. Therefore, since the water dissociation is energetically competitive with the molecular adsorption, the hydrolysis limit (i.e. maximum number of water molecules that can be dissociated upon hydration of titania nanocluster) is likely to be in equilibrium between two and three dissociated water molecules rather than simply two as predicted in ref [32]. One factor that can affect the hydrolysis limit of titania nanoclusters is the presence of "unsaturated" Ti centres. Such species are either under coordinated Ti atoms with less than 4-fold coordination, and are mainly present in small nanoclusters (i.e. M ≤ 4) and/or highly distorted tetracoordinated Ti centres. In M = 4, N = 3 systems, the hydrolysis limit appears to be determined through a competition between factors that favour the fully dissociated form such as the cage strain and hydroxyls hydrogen bond and factors favouring the molecular adsorption of water such as the Ti-O(H2) Lewis acid-base interaction. We note that other factors affecting the hydrolysis limit such as the hydroxylation of bridging oxygens, flexibility of Ti centres to adopt different coordinations (i.e. 4-, 5-and 6-coordinated) etc., appear as the system size increases. The contribution of such factors is not trivial to evaluate, however, (TiO2)8(H2O)5 isomers found from our search, are small enough to qualitatively highlight their relevance. In Fig. 6 we show three different (TiO2)8(H2O)5 isomers for which in the most energetically stable case (i.e. left-side of Fig. 6), the hydration occurs via fully dissociated water molecules. In this case the hydration reaction induces hydroxylation of bridging oxygens and a small degree of structural rearrangement in which three Ti centres increase their coordination from 4 to 5-fold. This structure is energetically favoured over nanoclusters which adsorb water molecularly (Fig. 6 middle isomer) or in which the Ti centres are all 4 coordinates and the cage-like cluster is opened (Fig. 6 right-hand isomer).

We further notice that the hydrolysis limit, in addition to being system size dependent, is also temperature and water vapour pressure dependent. Indeed, as the temperature increases, gasphase titania nanoclusters with the inclusion of N waters displaying one molecularly adsorbed H2O (i.e. (TiO2)4(H2O)4 and (TiO2)16(H2O)8) will present progressively more positive and therefore endothermic ∆Ghyd with respect the nanocluster with N-1 waters, as shown from hydration free energies in Fig. 4.

Aggregation of hydrated nanoclusters

Using the data in the previous section we can compare the energetics of nanocluster aggregation for silica and titania, and how it depends on temperature. For simplicity, we consider only the aggregation reactions of ( ) ( ) nanoclusters to form a twice larger nanoclusters with the same degree of hydration, i.e. reactions of the type: above each plot. The energy of the initial aggregating species is set to zero in all cases for ease of comparison of the energy profiles.

2( ) ( ) ⇄ ( ) ( ) = X, YX (
In Fig. 7 we plot the energies of aggregation of (XO2)4(H2O)N species with N=1 (left) and 2 (right)

where in each case the (free) energy of the corresponding global minimum candidate nanocluster is used. The corresponding energy profiles thus track the energies of nanoclusters with a RN/M value of 0.25 and 0.5 respectively. For 25% hydration the aggregation reactions at 0K for both hydrated titania and silica systems are all highly exothermic (between 5.2-8.6 eV). Comparatively, however, the energies for hydrated silica aggregation are always significantly more exothermic (by between 2.1-3.5 eV) than the corresponding hydrated titania aggregation. The energy drop upon aggregation of 25% hydrated silica nanoclusters is found to be monontonically decreasing with increasing nanocluster size, with a slightly higher exothermicity for the aggregation of two (SiO2)8(H2O)2 species as compared to that of two (SiO2)4(H2O) species. However, for the 25% hydrated titania nanoclusters, the aggregation of two (TiO2)8(H2O)2 is between 0.6-1.45 eV less exothermic than for the aggregation of two (TiO2)4(H2O) species. With increasing temperature all energy profiles are still always exothermic but are shifted upwards in energy. These reaction energy profiles can be qualitatively rationalised by inspection of the comparative tendencies of silica and titania nanoclusters in the free energy of hydration plots in Fig. 4, and noting that at 25% hydration we are always in the highly exothermic pre-saturation regime.

In the case of nanoclusters with 50% hydration all aggregation reaction energy profiles are qualitatively similar in that the aggregation of two (XO2)4(H2O)2 species is always more energetically favourable than the aggregation of two (XO2) the energy decrease upon the initial aggregation of two (XO2)4(H2O)2 nanoclusters at 0K is quite moderate and very similar in magnitude for both silica (2.3 eV) and titania (2.7 eV). Interestingly, in this case both silica and titania nanoclusters exhibit the same global minimum isomer (except for the orientation of hydroxyl groups) for both reactant and product.

For small systems such as those with M=4 that aggregate to form a larger nanocluster with M=8, the exothermicity of the reaction mainly arises by lowering the nanocluster ring and cage tension present in the reactant.

For the aggregation of slightly larger hydrated nanoparticles such as those having M=8, for poorly hydrated systems (i.e. (XO2)8(H2O)1 and (XO2)8(H2O)2) there is a strong aggregation energy. In such systems the products of their aggregation do not display unsaturated dangling oxygens contrary to their reactant. Such behaviour is particularly enhanced in the aggregation of silica nanoclusters given the destabilizing effect of dangling oxygens.

Electronic energy gaps

(SiO 2 ) M (H 2 O) N (TiO 2 ) M (H 2 O) N a) b)
In Fig. 8 we show the vertical electronic energy gaps of silica (Fig. 8a) and titania (Fig. 8b) systems with respect to RN/M. We note that our discussion is only related to trends and not on absolute quantitative values of the electronic gaps for which DFT calculations are often unreliable. Generally, silica nanoclusters show a much higher variability in their gap magnitudes with respect to the degree of hydration than titania. Indeed, as example case, in (SiO2)12 systems, the magnitude of the gap ranges from 4 eV for the anhydrous case to around 9 eV at 33% hydration, while for the corresponding titania case the gap varies from a minimum of 5.3 eV at 8% hydration to a maximum of 6.2 eV at 33% hydration.

In silica nanoclusters the electronic energy gap is extremely sensitive to the type of oxygen defective sites. For instance, for bare nanoparticles with M = 4, 8 and 16 and N = 0, which have defective oxygen atoms in silanone species, we find that the energy gaps range between 6.5 and 7

eV, while for (SiO2)12 which has NBO defective sites the gap is only 4 eV. Generally, with increasing size these defects will be saturated with hydroxyl groups tending to heal the defects. Thus, as the hydration increases, the gap generally increases until a maximum value around 9 eV is reached at around 35% and 40% of hydration depending on the system size (50% for the M = 4 system). Further increase of hydration does not affect the gap magnitude with the only exception of (SiO2)12(H2O)6 which shows a gap decrease of 3 eV at the highest hydration considered, probably due to the presence of geminal hydroxyl groups.

In titania systems there are only NBO defective sites, however there is much greater variety of titanium cation coordination (i.e. 4-, 5-and 6-coordinated Ti centres) and therefore the change in the gap magnitude with increasing hydration is less easy to rationalise than in the silica case. Overall, the titania gap magnitude evolution shows some similarities with silica systems: (i) there is a maximum value of the band gap as hydration increases and (ii) after the maximum the band gap remains rather constant (i.e. M=4 and 8). However, in some cases there is a slight decline (i.e. M=12 and 16) in the band gap.

Generally, the larger the gap magnitude the more electronically stable the system should be., Adsorption is known to help stabilise a large band-gap. 63 In both silica and titania nanoclusters we find that hydration contributes to maximise the gap magnitude and to stabilize the electronic state. In silica nanoclusters this effect is particularly pronounced showing that silica is more readily stabilized by hydration than titania. Conversely, our results also highlight the remarkable stability of the electronic state of titania nanoclusters with respect to hydration.

Ab initio thermodynamics of hydration

According to equation 2.6 (∆G(T, P)) it is possible to define the hydration reaction of clusters as function of both temperature (T) and partial pressure of water vapour (P). Since the most thermodynamically stable nanocluster at specific conditions of T and P is the one that minimizes the free energy of hydration, we can plot hydration phase diagrams where each region represents the most thermodynamically stable nanoclusters with specific amount of incorporated water within a certain T and P range. All necessary energy contributions can be derived from our DFT calculations within the framework of classical statistical thermodynamics through the canonical ensemble (NVT) partition function formulation as described above. Due to the assumed applicability of Boltzmann statistics and the harmonicity of molecular vibrations, we consider a T range between 100 and 1000 K for where the validity of the approach should be acceptable, in line with other studies. 64 Moreover, in our case we consider water as an ideal gas for a wide range of temperatures, although this assumption depends on the water partial pressure. Water vapour is in equilibrium with its condensed state at ~3kPa at 300K. In order to consider water as an ideal gas, however, we should consider a much lower partial pressure. In the phase diagrams reported in Fig. 9, the water partial pressure that we consider does not go above the 1x10 5 Pa (e.g. ~1 atm). At this pressure we can consider water as an ideal gas from ~450-500 K and higher. The saturation water pressure at 273K (0 C) is about 0.6 kPa. In order to treat water as ideal gas in these conditions we should consider either temperatures higher than 400K or pressures lower than 1Pa. Although considering water as an ideal gas at relatively low temperatures and high pressures is conceptually wrong, these conditions represent a small region of the phase diagrams reported in Fig. 9 and therefore do not affect most of the P versus T range considered.

(N=4, 100%) (N=4, 50%) (N=3, 38%) (N=3, 38%) (N=2, 25%) (N=0, 0%) (N=4, 50%) c) d) (N=0, 0%) (N=5, 63%) (N=4, 33%) (N=4, 33%) (N=0, 0%) (N=5, 42%) e) f) (N=6, 50%) (N=5, 42%) (N=2, 17%) (N=4, 25%) (N=5, 31%) (N=0, g) h) (N=8, 50%) (N=6, 38%) (N=7, 44%) (SiO 2 ) M (H 2 O) N (TiO 2 ) M (H 2 O) N M=4 M=8 M=12 M=16
Pa (from ~1x10 -15 to ~1 atm). Blue and black dashed lines indicate the equilibrium vapour pressure with respectively ice and liquid water respectively. The colours of the P-T phase diagram regions are determined by hydration degree intervals in all the graphs. The cooler the colour (towards dark blue) the lower the hydration degree, the warmer the colour (towards dark red) the higher the hydration degree. From top to bottom graphs are arranged from that smallest sized nanoclusters (XO2)4(H2O)N to the largest sized nanoclusters (XO2)16(H2O)N, where X=Si, Ti.

In Fig. 9 we plot P versus T phase diagrams for both titania and silica systems in thermodynamic equilibrium with water vapour. As discussed above, the pressure and temperature ranges have been chosen in order that water may be considered as an ideal gas in most of the conditions represented in the plots. Each region of each plot in fig. 9 represents the most thermodynamically stable hydrated cluster at the particular conditions of T and PH2O according to equation 2.6. The region are coloured consistently for each hydration degree in each plot (e.g. the 50% hydrated structures have same colour plots). The phase diagrams show that, silica and titania, which can have very different properties at macroscopic scale (amorphous silica can be hydrophilic and hydrophobic, while TiO2 is always hydrophobic), have quite similar hydration behaviour at the nanoscale. This similarity is particularly pronounced at the smallest size considered, M = 4, and becomes less pronounced as the system size increases. For the systems with M = 4, most of the phase diagrams (see Fig. 9a and9b) are dominated by the same tetrahedral isomer for both hydrated titania and silica (i.e. M=4; N=2 with 50% hydration degree.). However, the most stable structure at ambient temperature conditions for silica corresponds to the 75% hydrated system. At relatively higher PH2O and lower T silica becomes 100% hydrated while titania barely reaches such a high degree of hydration. as temperature increases. In the M = 8 phase diagrams we note that: (i) systems with N=1 do not appear in the P-T range considered in both titania and silica, and (ii) at very low temperatures and moderate to high pressures titania has a larger degree of hydration than silica. This second feature is not exclusive for M=8, but rather appears in all phase diagrams. A possible explanation is that TiO2 is able to adsorb water molecules more easily than SiO2 (SiO2 needs formation of hydroxyls, TiO2 protonates bridging oxygens). It must be noted though that the such very low temperature regions tend to fall outside that within which water vapour can be strictly considered to be an ideal gas.

For M = 12, the silica phase diagram (Fig 9e) is mainly populated by two structures as in the case of M = 8 but with a the relatively lower hydration degrees of 42% and 33%, compared to the previous hydration degrees of 50% and 38% at M = 8. Less hydrated structures are found at higher temperatures (> ~700K). For titania systems of the same size (Fig 9f) the phase diagram is slightly more complex than silica and has more stable phases. Despite that, only 33% and 0% (i.e. N=4 and N=0 respectively) hydrated phases seems to play a dominant role. We note that (SiO2)12(H2O)5 and (TiO2)12(H2O)4 which occupy most of the area of their respective phase diagram and are close to the region representing ambient conditions contain the cube motif (i.e. as exemplified by the M=8, N=4 structures) incorporated in the cluster. This is an indicator of the stability of such a motif, not only as known for silica systems, 22,23,25,40 but for titania systems as well. This suggests that such a structural motif could thus perhaps be used as a building block to synthesize more complex materials such as titanosilicates in a controlled way. For M=16, silica systems (Fig 9g ) show one dominant phase at 31% hydration degree (i.e. N = 5). The maximum hydration degree reaches 50% for narrow regions in the phase diagram at low temperature as in the case of M =12. As the temperature rises, the hydration decreases to 13%. This trend is similar to the case of titania systems (Fig 9h) but is shifted towards lower hydration degrees (i.e. from 25% to 0% hydration). Only the low temperature M = 16 titania phases (50% and 38% hydration) show similar stabilities with silica phases, albeit with a larger area of the respective phase diagram for the most hydroxylated case for titania.

Overall, the P-T phase diagrams for both systems highlight that the the tendency to be hydrated decreases with size. As the system size increases, the dominant colours in the phase diagrams gradually blueshift indicating a tendency towards lower degrees of hydration for both silica and titania. However, the shift towards less hydrated systems with size increase is less pronounced for silica than for titania systems indicating that silica nanoclusters are more hydrophilic. Indeed, at high temperature and low water vapour partial pressure silica becomes completely anhydrous at much higher temperatures than titania. Indeed, regions of the P-T diagrams with completely anhydrous silica are hardly present in the T and P range considered. We also note that at ambient temperature conditions, even with the presence of trace humidity, clusters tend to hydrate. Also, for all conditions of temperature and water vapour pressure considered in this work in which the hydration of nanoclusters is thermodynamically favoured, the adsorption of water molecules occurs exclusively in a dissociative way. In Fig. 10 we plot the calculated infrared harmonic vibrational frequencies of all global optimized titania (black line) and silica (blue lines) nanoclusters. In all spectra the region of hydroxyl O-H stretching (νOH) around 3900 cm -1 is easily discernible. This peak is persistent with system size and hydration degree and is consistent in frequency in both titania and silica. As the hydration increases there are other stretching modes between 3000 and 3900 cm -1 due to involvement of OH groups in hydrogen bonds (νOH--O). As hydration increases, for a fixed system size M, the νOH intensity increases with the increasing number of OHs. As the system becomes larger (increasing M) the OH peaks weaken because the spectra are dominated by the more numerous cation-O vibrational modes.

Vibrational frequencies

Differences between silica and titania cannot be easily appreciated by inspection of the hydroxyl stretching region but rather in the cation-O stretching/bending one. Indeed, between 600 to 1400 cm -1 in the spectra there are intense IR peaks mostly dominated by complicated coupling of cation-O stretching and bending modes involving several centres (δν-O-X-O). In silica systems this set of peaks is highly blue shifted (i.e. mostly with vibrational modes higher than 1000 cm -1 ) with respect to titania ones (i.e. mostly with vibrational modes lower than 1000 cm -1 ). Moreover, these δν-O-X-O become slightly blue shifted in the case of silica systems as the hydration degree increases. O distance distributions used to justify the choice cut-off distance adopted in this work. In all cases Si atoms tends to display at most a 4-coordinated anion environment rapidly reaching the bulk value as the number of included water molecules increases. In the case of TiO2, M = 4 and M = 8 systems Ti cations tend to display similar coordination environments and after a small increase with the addition of the first one or two water molecules their ACN does not change. This behaviour is mainly due to the large structural strains possessed by these small systems, which is mitigated by the nanoclusters opening up after a dissociative water adsorption rather than allowing for increased Ti coordination. Larger systems, such as M=12 and 16, are found to be more structurally flexible allowing for higher coordination of Ti cation centres without opening up their structures. Here the CAN of Ti centres tends towards the bulk value as the number of included water molecules increase.

Average Coordination Number (ACN)

However, as shown from Fig. 4 and Fig. 8, the hydration reaction free energy for highly hydrated systems (close to RN/M~50%) tends to be unfavourable, therefore such species will not form under realistic equilibrium conditions.

Conclusions

In this work we have investigated the hydration of titania and silica nanoclusters and their associated relative structural end energetic properties of the by means of accurate ab initio modelling. All nanocluster structures were generated by employing an IP-based global optimization search which yielded several new global minima candidate isomers for (TiO2)M(H2O)N with M=4, 8, 12, 16, each for a range of N values. In order to improve the efficiency reliability of these searches a newly parameterised IP for hydrated titania nanosystems was also developed. For the nanocluster size regime considered the hydration energy at 0K is found to be exothermic over a wide range of hydration both for nano-SiO2 and nano-TiO2 with characteristic dependencies on nanocluster size and degree of hydration. We have also calculated hydration free energies which are found to thus showing how the 0K exothemicity of hydration decreases with increasing temperature. In particular, some of the highly hydrated titania systems are found to be thermodynamically unstable as the temperature increases. By ab initio atomistic thermodynamics we also provide phase diagrams showing the stability dependence on both temperature and water vapour pressure for hydrated silica and titania nanoclusters. Our study shows that under ambient conditions, silica and titania nanoclusters would be hydrated even with trace amounts of water vapour. In addition to the thermodynamic properties we also calculate how the structural and electronic properties of all these considered nanoclusters depend on the degree of hydration and nanocluster size. Overall, our study provides a wide range of detailed new insights into the interaction of water with titania and silica

  refer to this simplified version as FFTiOH*. To overcome bias of using a single IP in the MCBH global optimisation procedure we also slightly modified the degree of O-O repulsion in FFTiOH* to generate FFTiOH*mod. The use of FFTiOH*-mod generates slightly more open cluster structures allowing the MCBH searches to explore PES regions not accessible with only FFSTiOH*. The full set of parameters for FFTiOH, FFTiOH* and FFTiOH*-mod are provided in the ESI.

2 (

 2 ) = -3 4 567#( , )8 + 9 ≃ ;( , ) $%&'( + ;( ) %)$ + ;( ) *+, + ;( ) ( -+ <= ;( , ) $%&'(

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Global optimized (TiO2)M(H2O)N and (SiO2)M(H2O)N clusters for M = 4 and 8. Numbers in the cluster labels show M, the number of XO2 units, and N, the number of water molecules incorporated in the cluster. (SiO2)M(H2O)N clusters are obtained from refs 37,39,40 except for the (SiO2)8(H2O)5 nanocluster which was found in this study.

Figure 4 .

 4 Figure 4. Hydration energies free energies of silica (dashed lines) and titania (solid lines) nanoclusters for sizes for: a) M=4, b) M=8, c) M=12, d) M=16. Hydration free energies are reported

Figure 5 .

 5 Figure 5. Titania nanocluster structures reported from Dixon et al. in ref[32] (right side) as compared with geometries found in this study (left side). In parenthesis we report our calculated relative electronic energies in eV.

Figure 6 .

 6 Figure 6. (TiO2)8(H2O)5 isomers found from our global optimization search. Number in parentheses show our calculated relative total energies in eV.

Figure 7 .

 7 Figure 7. Reaction energies for the aggregation of (XO2)4(H2O) nanoclusters (left) and of (XO2)4(H2O)2 nanoclusters (right) at temperatures of 0K, 300K and 600K. The corresponding reactions are shown

Figure 8 .

 8 Figure 8. Electronic energy gaps of silica (a) and titania (b) nanocluster systems with respect to the incorporated water ratio, RN/M.

Figure 9 . 5 (

 95 Figure 9. P-T phase diagrams for silica (left side) and titania nanoparticles (right side). The temperature range goes from 100 to 1000 K and P (H2O partial pressure) range from 1x10 -10 to 1x10 5

  For M = 8, the silica phase diagram (Fig 9c) is dominated by two structures (i.e. for N=3, 4) while titania phase diagram (Fig 9d) shows four prominent structures (i.e. for N=0, 2, 3, 4). The M = 8, N = 4 isomer has the same structure again for both systems but for silica the P-T stability region is wider. At ambient conditions (i.e. T=300 K and 1k Pa water vapour partial pressure -the equilibrium vapour partial pressure with its condensed state at 300K is 4.2 kPa) this last structure is thermodynamically favourable for both systems. At temperatures below 900K silica is stablised by losing one water molecule to form the M = 8, N = 3 structure. Titania for M = 8 shows a relatively higher variety of hydration degrees at high temperature, with structures going from N = 3 to N = 2 and finally to N = 0

Figure 10 .

 10 Figure 10. Simulated, infrared harmonic vibrational frequencies of silica (blue lines) and titania (black lines) systems for different numbers of incorporated water molecules (hydration degrees in parenthesis). a), b), c), d) represents IR spectra of system sizes M = 4, 8, 12, 16 respectively. Each plot is scaled with respect to the most intense peak. Intensities are in arbitrary units.

Figure 11 .

 11 Figure 11. Average Coordination Number (ACN) as a function of N incorporated water molecules for Si and Ti cations. Respective bulk coordination numbers (for quartz and rutile respectively) are also indicated by dashed lines. The CAN number is defined by taking all nearest neighbours within a 2.1 Å cut-off in the case of SiO2 nanoclusters and within a 2.5 Å cut-off in the case of TiO2 nanoclusters.

  TiO 2 ) 16 (H 2 O) N non-cage-like cluster form cages and are more resistant to dissolution.39 For these larger sizes and at at high degrees of hydration, silica nanoclusters start to absorb water rather in a molecular fashion rather than form germinal silanols (e.g. for (SiO2)24 see ref.39). Atomic coordinates of all clusters reported in Figs 1, 2 and 3 are reported in the Supporting Information.

	1: M=16, N=0	2: M=16, N=1	3: M=16, N=2	4: M=16, N=3	5: M=16, N=4

6: M=16, N=5 7: M=16, N=6 8: M=16, N=7 9: M=16, N=8 (SiO 2 ) 16 (H 2 O) N (

nanoclusters and how their behaviour compare with one another. We hope that our comprehensive and systematic work will inspire experimental studies, especially from the cluster beam community, in order to further our understanding of nanoscale hydration of these important systems.

Acknowledgements

This research was supported by the Spanish MINECO / AEI-FEDER grant: CTQ2015-64618-R and, in part, by Generalitat de Catalunya grants: 2014SGR13 and XRQTC. We acknowledge the Red Española de Supercomputación for provision of computing time. AC acknowledges PhD funding from the ITN-EJD in Theoretical Chemistry and Computational Modelling with EU H2020 Marie Curie Sklodowska-Curie grant agreement No.642294.

Conflict of interests

There are no conflicts of interest to declare

Supporting information

In Supporting Information, we report Ti-O and Si-O distance distributions over the whole set of reported global minima candidate nanoclusters. For these nanoclusters we also include, atomic coordinates, total energy, hydration energy at 0K and hydration free energy at 300K. Atomic coordinates will be also uploaded to the open access Nomad repository (https://nomadrepository.eu/).


F Decremps

J Pellicer-Porres

A M Saitta

J.-C Chervin

A Polian

High pressure Raman spectroscopy study of Wurtzite ZnO

Keywords: numbers: : 78.30.-j, 64.70.Dv, 91.60.GF

The high pressure behavior of optical phonons in wurtzite zinc oxide (w-ZnO) has been studied using room temperature Raman spectroscopy and ab-initio calculations based on a plane wave pseudopotential method within the density functional theory. The pressure dependence of the zonecenter phonons (E2, A1 and E1) was measured for the wurtzite structure up to the hexagonal→cubic transition near 9 GPa. Above this pressure no active mode was observed. The only negative Grüneisen parameter is that of the E low 2 mode. E1(LO) and (TO) frequencies increase with increasing pressure. The corresponding perpendicular tensor component of the Born's transverse dynamic charge e * T is experimentally found to increase under compression like e * T (P ) = 2.02 + 6.4 • 10 -3 .P whereas calculations give e * T (P ) = 2.09-2.5•10 -3 .P (in units of the elementary charge e, P in GPa). In both cases, the pressure variation is small, indicating a weak dependence of the bond ionicity with pressure. The pressure dependence of the optical mode energies is also compared with the prediction of a model that treats the wurtzite-to-rocksalt transition as an homogeneous shear strain. There is no evidence of anomaly in the E2 and A1 modes behavior before the phase transition.

Zinc oxide belongs to the wide band-gap semiconductor family with large ionic characters of chemical bonds 2 . The ionic size or relative electronegativity has been used to explain high pressure structure changes in A N B 8-N semiconductors. First principles calculations have shown that the zinc-blende (or wurtzite)→rocksalt→ β-Sn transition sequence is a common feature for most semiconductors. However, recent experimental identification of new intermediate phases (like cinnabar in CdTe for example) and the systematic absence of some phases (rocksalt (rs) for covalent compounds and β-Sn for all except the most covalent) invalidate the traditional transition sequence and, consequently, question the standard theoretical approach. In a recent Letter, Ozoliņš and Zunger 3 suggested that phase transitions are accompanied by phonon softening, whose instability has to be taken into account in the calculation to correctly predict the phase diagram. The latter outcome can also be discussed in terms of the transition mechanism, where the presence of negative Grüneisen parameters of phonon modes may play an important role. Recently, a new intermediate structure has been proposed in the wurtzite to rocksalt transformation path of GaN, along which the optical A 1 and E high 2 modes are expected to be affected by the new bond formations 4 . Wurtzite ZnO, which transforms under pressure into the rs phase at 9 GPa, may provide a good example of this trend.

Wurtzite ZnO belongs to the C 4 6v (P6 3 mc) space group. The primitive cell includes two formula units, with all atoms occupying 2b sites of symmetry C 3v . At the Γ point of the Brillouin zone, group theory predicts the existence of the following optics modes:

Γ opt =A 1 +2B 1 +E 1 +2E 2 .
Frequency from the B low 1 and B high 1 silent modes has been calculated at 260 and 540 cm -1 respectively. A 1 , E 1 and E 2 modes are Raman active. In addition A 1 and E 1 are infrared active, and therefore split into longitudinal and transverse components (LO and TO). The mode assignment at ambient conditions is well established in the literature 5 . To the best of our knowledge, no calculation and only two experimental attempts have been made to study the w-ZnO phonon frequency shift under pressure 6,[START_REF] Minomura | High Pressure in Science and Technology[END_REF] . In the first reference, the authors report the evolution of two phonon frequencies, E high 2 and E low 2 , over a relative small pressure range (0-1 GPa). In the second one, the pressure dependence of four Raman modes (E high 2 , E low 2 , A 1 (TO) and E 1 (TO)) are given up to ∼7.5 GPa. However no experimental details is described, neither on the sample nor on the apparatus (high pressure cell, pressure transmitting medium or pressure measurement). This lack of information is detrimental to any discussion, especially in our case where, on compression at room temperature, nonhydrostatic stresses are known to affect considerably the phonon response of the crystal. Here, we report the results on the high pressure behaviour of optical phonons in w-ZnO from both Raman experiments and ab-initio calculations.

For the high pressure Raman experiments, two samples were cut as platelets from a large, colorless and transparent single-crystalline cube of w-ZnO purchased from SPC Goodwill. The sample was tailored to about 30x20 µm 2 surface and 20 µm thick. One crystal was polished with plane parallel surfaces normal to the hexagonal c-axis, and the other one with surfaces normal to the a-axis. The high pressure cell was a membrane diamond anvil cell (DAC) [START_REF] Letoullec | [END_REF] . A stainless steel gasket was preindented to 50 µm and a 200 µm hole was drilled in the center by spark-erosion. Neon gas loaded at high pressure (0.15 GPa) was used as a pressure transmitting medium because : (i) it is hydrostatic up to 16 GPa 9 , (ii) it is chemically inert, (iii) it has no luminescence and no Raman activity. Pressure was systematically measured be-fore and after each Raman spectra using the fluorescence emission of a ruby sphere 10 placed close to the sample into the gasket hole. The accuracy was better than 0.1 GPa at the maximum pressure reached. Raman spectra were excited with the 514.5 nm line of an Ar + laser focused down to 5 µm with a power level of about 500 mW at the entrance of the DAC. The scattered light was analyzed with a Dilor XY triple spectrometer and a liquidnitrogen-cooled CCD multichannel detector. All spectra were recorded in the backscattering geometry with unpolarized light. In this geometry, E 2 , A 1 (TO) and E 1 (TO) are allowed with incident radiation perpendicular to the c-axis, and E 2 , A 1 (LO) with light parallel to the c-axis. The Raman frequencies were determined from a computer fit of the peaks with a Lorentzian profile. The accuracy was better than 1 cm -1 .

Phonon frequencies, Born effective charges, and the dielectric constant ∞ have been calculated from firstprinciples in the framework of density-functional theory and local-density approximation, by using the PWSCF code 11 . We adopted the density-functional pertubation theory (DFPT) approach 12 and a standard local-density approximation plane-wave/pseudopotential scheme.

Norm-conserving pseudopotentials of the Troullier-Martins form 13 have been used, with the inclusion of 3d electrons of Zn in the valence shell. Tests carried out with zinc 3d electrons in the core and nonlinear core-correction gave an unsatisfactory accuracy on phonon frequencies. We used a cutoff of 70 Ry and a 8×8×4 mesh of special k-points for Brillouin zone integration. At each pressure, w-ZnO has been relaxed with respect to the c a ratio and the internal parameter u, and phonon frequencies have been calculated at these optimized structural parameters. Because of the polar character of ZnO, the dynamical matrix displays nonanalytic behavior in the limit q→0, arising from the long-range character of the Coulomb forces. The calculations of LO modes is nonetheless straightforward, since the LO-TO splitting only depends on the phonon frequencies, the dielectric constant and the Born effective charges, which can all be directly calculated in the DFPT framework.

Before performing the high pressure experiment we examined the optical zone-center phonons of w-ZnO at ambient pressure with light parallel and perpendicular to the c-axis. We employed the large single crystal from which the micrometric samples used in the high pressure experiments were extracted. The E 2 , A 1 (TO) and E 1 (TO) modes are clearly seen in the spectra. On the contrary, the A 1 (LO) mode is hardly observable and not seen at high pressure (crystal in the DAC). The frequency assignment, as well as a comparison with ab-initio results and previous works, is presented in Table I. The Raman spectra in both configurations also show additional features at 331 and 552 cm -1 . In reference 5 they are tentatively assigned to second order structures by comparison with flat regions in the dispersion curves.

Under pressure, with incidence parallel to the c-axis, only E 2 single phonon mode was apparent when focusing at the center of the sample. In the same configuration, A 1 (TO), E 1 (TO) and a new mode at 580 cm -1 appear with the laser light focused on the sample border. The new mode is much more intense than the A 1 (LO) found in ambient conditions. Its frequency is also slightly higher.

The E 1 (LO) mode is forbidden in the backscattering geometry. However, it has been observed in ZnO in the forbidden X(ZZ)X geometry 5 . We have thus assigned the mode appearing at 580 cm -1 to the E 1 (LO) mode.

Our assignment should be considered as tentative since the new mode is observed only when we focused on the sample border. In this situation, the incidence direction with respect to the crystal axes is not defined and mixed modes may be obtained. In that geometry, typical Raman spectra of ZnO at various pressures are shown in Fig. 1. The pressure dependence of the phonon frequencies for the two high-pressure runs are shown in Fig. 2. No difference between the two runs is observed. Around 8.7 GPa, the disappearance of all Raman peaks reflects the onset of the wurtzite-rocksalt structural transformation in excellent agreement with the previous high-pressure studies. E low 2 is the only mode which exhibit a negative pressure dependence. The frequency shift of E high 2 , A 1 (TO), E 1 (TO) and E 1 (LO) modes increase with pressure and no anomaly up to the transition pressure is detected. The frequency dependence of E low 2 and E high 2 are in perfect agreement with the low pressure values obtained from Mitra et al 6 . On the other hand, while our frequency data seem to agree with those of Minomura 7 (except A 1 (TO)), the discrepancy in the corresponding Grüneisen parameters is important (except for E 1 (TO), see Table II). It is difficult to account for that disagreement since no calculation details are given in Ref. [START_REF] Minomura | High Pressure in Science and Technology[END_REF] . We suspect the author to have used a different bulk modulus. Here, the mode Grüneisen parameters γ i are defined as: γ i = -dlnνi dlnV P =0 = B0 νi dνi dP . The isothermal Birch-Murnaghan equation of state was used to determine B 0 =170 GPa and B'=4 from the X-ray diffraction data V(P) 15 .

The experimental LO-TO splitting of the E 1 phonon mode increases with a small pressure dependence (+2% at 5 GPa, see Fig. 2), in good agreement with theoretical results where no variation of the E 1 or A 1 LO-TO splitting is obtained in the same pressure range. It is noteworthy that this splitting increases (nearly constant from calculations) with pressure, an uncommon behavior for A N B 8-N semiconductors except SiC 16 , AlN and GaN 17 . The change of the Born's transverse dynamic charge e * T under compression can be determined from the measured frequencies using (in S.I. units):

(e * T ) 2 = 4π 2 V µε 0 ε ∞ (ν 2 LO -ν 2 T O
), where V is the volume per formula, µ is the reduced mass, ε 0 is the vacuum permittivity and ε ∞ is the high frequency (optical) dielectric constant [START_REF] Born | Dynamical theory of crystal lattices[END_REF] . At ambient pressure, and with ε ∞ (P=0)=3.95•ε 0 , the experimental effective charge values are 2.02 and 2.4 (in units of the elementary charge e) for the E 1 and A 1 modes respectively. The evolution of the high frequency dielectric constant ε ∞ under high pressure has not been measured. To our knowledge only the dependence of the ordinary (n o ) and extraordinary (n e ) refractive indices at 0.5893 µm has been measured [START_REF] Vedam | [END_REF] . In order to estimate the pressure behavior of e * T from our measurements, we suppose that the pressure dependence of ε ∞ and n 2 o are equal. With this assumption, we obtain (in S.I. units): ∂ε∞ ∂P P =0 = 2n 0 ∂n0 ∂P P =0 = -0.014 • ε 0 GPa -1 . The behavior under pressure of the normalized perpendicular tensor component e * T /e * T (0) is plotted in Fig. 3. The effective charge in w-ZnO increases under compression with ∂e * T /∂P = 6.4 • 10 -3 (in units of the elementary charge e per GPa). Theoretically, ε ∞ and the effective charge have been found to decrease linearly with pressure: ∂ε∞ ∂P P =0 = -0.025 • ε 0 GPa -1 and ∂e * T /∂P = -2.5 • 10 -3 . However, if we constrain the calculations to use the experimental pressure dependence of ε ∞ , we obtain ∂e * T /∂P 0. These results may indicate a small variation of the w-ZnO bond ionicity when compressed between 0 and 4 GPa. In comparison, the ∂e * T /∂P experimental value for SiC 16 is 5.3 • 10 -3 , 0.15 • 10 -3 for AlN 17 and -1.3 • 10 -3 for GaN 17 . Above 4 GPa, Fig. 3 shows a saturation of e * T . However, the refractive indices (and thus the dielectric constant) have only been measured between 0 and ∼1 GPa [START_REF] Vedam | [END_REF] , which means that our e * T (P ) values above that pressure range have been obtained from an extrapolation of ε ∞ (P ). We now turn the discussion to the simple homogeneous orthorhombic shear strain path proposed to picture the w to rs mechanism of transition 4 . This model predicts that under pressure, a relaxation of the structural parameter u from ∼ 3 8 up to 1 2 is expected to occur in the wurtzite structure (i.e. below the phase transition) which gives rise to a new intermediate phase isomorphic to h-MgO. In addition to the expected (and experimentally observed) negative mode Grüneisen parameter of E 2 low at Γ, the previous model predicts the optic A 1 and E high 2 modes at Γ to be affected by the w →h-MgO transition. However, in the pressure range probed during the present work (from atmospheric up to the transition pressure, i.e. 9 GPa), no instability of these modes is observed. Our results show a discrepancy with the A 1 and E 1 anomaly expected from calculations and agree with the conclusions derived from long-range order probe (x-ray diffraction experiments) by Desgreniers 20 , where the evaluation of the normalized integrated intensity of the (002) and (101) diffraction lines indicated no deviation of the u parameter from its ideal value.

In summary, the LO-TO splitting of the E 1 phonon mode has been experimentally and theoretically found to be weakly pressure dependent. This unusual behavior is related to a small variation of the chemical bonds ionicity of w-ZnO with pressure. Moreover, the present observation of optical phonon under pressure does not show any softening of the optic A 1 and E high 2 modes as theoretically expected 4 . All Grüneisen parameters show positive values except that of E low 2 , which is a quite common feature in II-VI compounds. However, from first-principles calculation on w-ZnO, we predict (in good agreement with Ref. 4 ) u parameter to shift up to 0.5 at a volume of 24.08 Å3 (i.e. ∼21 GPa). Experimentally, the w-to-rs firstorder transition limits the pressure range of w structure existence and, consequently, does not allow the secondorder process observation (via optical phonons instabilities for example). However, if we refer to the phase diagram of ZnO 21 , it may be possible to constrain the structure to be still tetrahedrally coordinated at low temperature in order to make the wurtzite → NaCl transition to occur at higher pressure, rendering the second-order process directly observable for the first time in semiconductors.

FIG. 1 :

 1 FIG. 1: Raman spectra of single-crystalline ZnO under different hydrostatic pressures in the backscattering geometry. The spectra are shifted upward for clarity.

FIG. 2 :

 2 FIG. 2: Top: Pressure dependence of the observed optical phonons. Open (full) symbols: propagation of light perpendicular (parallel) to c. Bottom: (νLO -νT O ) E1 phonon mode splitting versus pressure. Solid lines are linear least-square fits to the experimental points.

FIG. 3 :

 3 FIG. 3: Experimental pressure and volume dependence of the perpendicular component e * T under compression normalized to its ambient pressure value e * T (0) = 2.02. For comparison, dotted lines represent the experimental values of e * T for SiC 16 , AlN 17 and GaN 17 .

TABLE I :

 I Experimental and theoretical Raman-active w-ZnO Γ-point phonon frequencies (in cm -1 ) at ambient conditions. Mode this work (exp.) this work (theo.) Ref.14 Ref.5

	E low 2	99	92	101 101
	E high 2	439	449	444 437
	A1(TO)	382	397	380 380
	E1(TO)	414	426	413 408
	A1(LO)	574	559	579 574
	E1(LO)	580	577	591 584

TABLE II :

 II

	E low 2	-1.6	-1.67	-1.75
	E high 2	2.0	1.84	1.62
	A1(TO)	2.1	1.70	1.70
	E1(TO)	1.8	1.80	1.76
	E1(LO)	1.4	1.30	-

ZnO mode Grüneisen parameters γi of the zonecenter phonons symmetry this work (exp.) this work (theo.) Ref.

[START_REF] Minomura | High Pressure in Science and Technology[END_REF]

Acknowledgments

We wish to acknowledge Bernard Couzinet for the neon DAC loading. One of us (A.M.S.) acknowledges the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) for computer time allocation (project 11387-CP9). This research has been partly supported by a Marie Curie fellowship of the European Community number HPMF CT2000-00764.


Rémi Federicci

Stéphane Holé

Vincent Démery

Brigitte Leridon

Memory effects in the ion conductor Rb 2 Ti 2 O 5

Recent studies on Rb 2 Ti 2 O 5 crystals have demonstrated remarkable electrical properties. This material exhibits colossal electrical polarization between 200 K and 330 K. In the present work, we report on the observation of memory effects in Rb 2 Ti 2 O 5 due to charge accumulation and we discuss the genuine memristive character of this material. An analytical model is proposed for the system, which takes into account the ionic diffusion and ionic migration and is in good agreement with the observed volatile memristive properties of the material.

I. INTRODUCTION

Memristive systems are ubiquitous in nature, from the simple tungsten filament 1 to the complex synaptic junction in neuronal systems 2 . Predicted by Chua in 1971 3 , and claimed to be discovered at the nanoscale in 2008 by Strukov and coworkers 4 , the memristor is the electronic component whose state equation is expressed by a linear dependence between the charge q and the flux Φ = V (t)dt. The consequence is that V = R(q)I where R (the memristance) depends on the total charge q received by the device. Therefore the history of the currents and voltages experienced by the material impacts its current internal state, giving rise to memory effects 5,[START_REF] Chua | Proceedings of the IEEE[END_REF] . In the recent years, the research about memristors and more generally about memristive, memcapacitive and meminductive systems-where respectively C or L are a function of charge or flux-has become a topic of huge interest especially in regards of neuromorphic implementation [START_REF] Wang | [END_REF][8][9][10][11] and data storage (ReRAM) 5,12,13 .

Recent studies have shown remarkable electrical transport properties in Rb 2 Ti 2 O 5 (RTO) crystals or ceramics 14 . An ionic conductivity up to 10 -3 S cm -1 between 200 K and 330 K was found, making the RTO a superionic conductor in this temperature range. Ionic transport in this highly electronically-insulating material creates accumulation of ions at the edges of the system, giving rise to an unscreened macro-dipole with a polarization larger than 0.1 C cm -2 and a colossal equivalent relative permittivity of about 10 8 at low frequency 14 . The possibility of a conventional ferroelectric phase associated to a structural phase transition was completely discarded by Raman spectroscopy and XRD measurements. Indeed, RTO structure was found to remain consistent with the centrosymmetric C2/m space group between 10 K and 450 K 15 .

In this article, we demonstrate intrinsic memory effects in bulk RTO at room temperature. We also discuss the intrinsic memristor character of the material and the volatile nature of the memory effects. In addition, a simple analytical model based on the electrolytic nature of the material is presented. The I-V curves calculated from the model are in very good agreement with the experimental measurements, and reproduce the observed non-transverse memristive behavior, further confirming the ionic mobility origin of these effects.

II. SAMPLE PREPARATION AND EXPERIMENTAL SET-UP

The RTO crystal synthesis was performed by melting TiO 2 and MNO 3 powders in a Pt crucible 14 and characterization was performed as described in previous work 14 . The crystals, which size is typically 1 mm × 300 µm × 300 µm, were then extracted from the synthesis batch and silver wires were glued with carbon paste at both ends of the elongated platelet-type crystals. The entire system was then annealed at 400 K under room atmosphere to polymerize the carbon paste. Once connected, in order to activate the electrical transport properties 14 , the crystals were annealed under vacuum at 400 K for two hours.

Due to the very high impedance of our samples, we used a two-contact geometry. Voltage biased I-V curves were acquired using a Lecroy LT224 oscilloscope for measuring the voltage drop across a resistance of 200 kΩ in series with the sample, and a Stanford SI1200 Lock-In was used as voltage source. For triangular voltage measurements and current relaxations a Keithley 6517b electrometer was used to measure the current while applying a voltage. All the transport measurements presented in the present article were performed under room pressure and room temperature conditions.

III. ELECTRICAL MEASUREMENTS

III.1. Memristive behavior

Voltage biased I-V curves were acquired for different frequencies ranging from 1 mHz to 500 mHz. Typical results are displayed in Fig. 1. The results show non linear non transverse -the two branches do not cross at the origin -hysteretical loops that are pinched in zero. The area of the lobes increases with frequency from 1 mHz to 10 mHz, and decreases above 10 mHz. When the frequency of the applied voltage reaches 500 mHz, the I-V curve shows an ohmic behavior with a single-valued conductivity of about 10 -4 S cm -1 .

For a system to be qualified as an ideal memristor, three conditions must be met 16 . The first is that the I-V curve displays a pinched (either transverse or non-transverse) hysteretical loop, i.e., I = 0 when V = 0. The second condition is that the lobe area enclosed by the I-V curve decreases with increasing frequency. And the third condition is that at infinite frequency, the I-V curve should be single-valued. The first criterium is met since to a first approximation the loops are pinched in zero. We discuss later some experimental cases when the loops are not exactly pinched in zero and when the I-V lobes are not exactly restricted to the first and third quadrants. The second criterion is not exactly fulfilled since the area of the butterfly lobes monotonously A) under room conditions for frequencies ranging from 1 mHz to 500 mHz. The sample is a single crystal which dimensions are 1000 µm × 300 µm × 300 µm. This crystal demonstrates genuine memristive behaviour, with an almost ideal "memristor" character according to the criteria defined by Adhikari et al. 16 (see text). decreases with frequency only at high frequencies. The third criterion for ideal memristor is fulfilled since an ohmic regime is observed at high frequency.

These findings demonstrate that a properly annealed crystal of RTO connected to conducting wires is a genuine memristive device, close to an ideal memristor, at low frequencies and at room temperature.

III.2. Memory-volatile character

As a matter of fact, as can be seen from the arrows drawn in Fig. 1, the I-V curves are of a non-transverse memristor, which is, for instance, different from all known examples of resistance switching memory devices based on ion conductors 12 . Consequently, it was observed in some of the samples, either single crystals or polycrystals, that the loop was slightly opened (not fully pinched)-which violates the first criterion. In this case when the bias voltage goes from positive to zero, a negative current is created at zero voltage, releasing charge into the outer circuit. We refer to this current as the "relaxation current". The necessary conditions to produce this cycle opening and the relaxation current (which varies from sample to sample) are not fully understood, but seem to be related to the annealing conditions, the presence of defects in the sample and possibly the nature of the electrodes. In this case, the device can be modeled by adding a memcapacitive component to the memristive one, and the system then partially behaves like a supercapacitor.

Fig. 2 displays a set of eight I-V curves measured on the polycrystalline sample B where the relaxation currents can be observed. The frequency of the applied field is about 1.2 mHz. The perfect coincidence between these eight curves demonstrates the reproducibility of the measurements. The inset in Fig. 2 shows a close-up around the origin in order to outline the relaxation current.

In order to better evidence the memristive properties of RTO, a triangular voltage signal was used to study the evolution of the internal state of the material as a function of time. The measured current (red) and the applied voltage (blue) are displayed as a function of time in Fig. 3. A serie of 5 negative cycles between -5 V and 0 V was performed, followed by a series of 5 positive cycles between 0 V and 5 V. These alternatively negative and positive cycles show how the conducting properties of RTO are impacted by the voltage history.

First, during the series of negative cycles, the peak current amplitude decreases regularly in absolute value between the first and the last cycle and the peak current for the last cycle is observed to be almost 10 times smaller than for the first one. Once the series of positive cycles starts, the conductivity of the material is suddenly increased, being almost five times higher than during the previous last negative cycle. During this series of positive cycles, the conductivity is also observed to decrease in a similar way as for the series of negative cycles.

After the set of alternatively negative and positive cycles, the applied voltage was reduced to zero during a time ∆t for a relaxation to occur. After this elapsed time, a series of 5 positive cycles was performed again to measure the internal conducting state of the RTO. The first cycle of the series exhibits the highest conductivity which is used to estimate the internal conducting state. The process was performed with different elapsed times and the current intensity of the first positive cycle was observed to increase with increasing elapsed time. This measurement shows that the conductivity of the RTO decreases with the increasing flux and that a slow leaky process exists when the applied voltage is zero which tends to restore the conductivity, attesting to the volatile nature of the memory effects at room temperature in this material.

III.3. Origin of the memory effects

Memory effects in devices comporting ion conductors are largely documented. Both unipolar and bipolar I-V characteristics are reported, depending on the microscopic nature of the underlying mechanims 12 . Memory effects are also largely studied in perovskite-based devices 13 .

A very commonly reported mechanism is the creation of microscopic conduction channels through the application of a high voltage, which provokes a local breakdown of the insulator 17,18 . After application of a threshold voltage, these channels are destroyed because of local Joule heating. This gives a very typical I-V curve which depends on the compliance current, with an either bipolar or unipolar switching type. In any case the I-V curve is always transverse, i.e. the two branches cross around the origin. This is not what is observed in our system.

A second possible mechanism relies on redox processes associated with cation electromigration. In this case, the process involves the oxidation of a metallic electrode such as Ag, the migration of the formed cations such as, e.g., silver ions Ag + through the ion conductor and the deposition at an inert electrode whose role is simply to provide electrons. This process, which creates silver dendrites, is reversible and leads to switching behavior, which was first reported in 1976 by Hirose 19 . In our case, the nature of the wires (Ag or Au) was varied without affecting the results. Gold contacts were evaporated on the cystals and used to perform the same measurements, leading to identical I-V curves so the role of Ag + ions can be excluded.

A third reported mechanism is electrochemical migration of oxygen ions. Such mechanism was observed for example in STO thin films 20 . It is very tempting indeed to consider such a mechanism in the present system, since we know that ionic motion is present, and the preparation process (annealing under oxygen-poor atmosphere and application of large voltage) reported for STO is very similar to the one applied on RTO. However, again, the I-V curves reported in STO systems [20][21][22] are transverse in contradiction to the present observation. The direction of variation of the cycle in the first quadrant is opposite to the one reported in Fig. 1.

As a matter of fact a notable difference between the two systems is present: in the above-reported mechanism, the mobile oxygen ions experience charge transfer through redox reaction at the anode interface. This means that the following reaction takes place,

O 2-→ 2e -+ 1 2 O 2 , (1)
and a production of oxygen gas is observed at the anode. In this case, ions are not accumulated and there is no building up of charge at the interface. This is in contradiction to what is observed in the present system where ions remain accumulated, a giant polarization (of the order of 0.1 C cm -2 ) is created 14 , and where no particular gas production is observed at the anode. Moreover, the phenomenon is also observed under helium atmosphere, which discards a possible role played by oxygen gas production and absorption. In order to explain the observed I-V curve, we therefore developped a simple model taking into account the ion mobility, ion accumulation, and supposing a total absence of charge transfer at the interface. This means that we assumed that redox reaction are completely absent at the interface where the ions accumulate and that for some reason remaining to be elucidated if the mobile ions are oxygen ions O 2-, which is very likely, they cannot escape the material and transform into oxygen gas.

IV. IONIC ELECTROMIGRATION MODEL

IV.1. Model

In order to better understand the microscopic mechanisms at play in this material, we compare the experimental results with a simple model of ionic migration. First, we give the general relations satisfied by this kind of system, and then make additional hypothesis to arrive at a system of equations that can be integrated numerically.

We assume that all the quantities only depend on the coordinate x ∈ [0, L], where L is the length of the sample, and that all the fields are oriented along this direction. The Poisson's law thus reads

∂ x E(x, t) = ρ(x, t) ǫ , (2)
where E is the electric field and ρ the charge density. The evolution of the charge density is given by the current j:

∂ t ρ(x, t) = -∂ x j(x, t). (3)
The electric field at the boundary x = L is related to the surface charge density σ(t) of the right electrode through

E(L, t) = - σ(t) ǫ . ( 4
)
The potential difference between the electrodes can be written as

V (t) = - L 0 E(x, t)dx = Lσ(t) ǫ + 1 ǫ L 0 xρ(x, t)dx, (5)
where we have integrated by parts and used Poisson's law. This equation shows how imposing V (t) sets the charge on the electrodes, given the distribution of charges. Moreover assuming that there are no charges entering or leaving the sample at the electrodes,

j(0, t) = j(L, t) = 0, (6)
we can derive the previous relation and get for the current

I(t) = A dσ dt (t) = Aǫ L dV dt (t) - A L L 0 j(x, t)dx, ( 7
)
where A is the area of the electrodes, or equivalently the area of the section of the sample. The current thus consists in two terms (Eq. ( 7)): the capacitive current I cap (t) and the ionic current I ion (t).

We now make further assumptions to describe ionic migration in the sample. We assume that there is a single species of mobile ions, carrying a charge q, and we denotes its number density n(x, t). We further assume that the material is neutral when the mobile ions are uniformly distributed in the material, meaning that there is a uniform background of non-mobile charges, with number density n. The charge density thus reads:

ρ(x, t) = q[n(x, t) -n]. (8)
We make the hypothesis that the mobile ions move according to the Nernst-Planck equation:

∂ t n(x, t) = -∂ x j n (x, t), (9) j n (x, t) = -D∂ x n(x, t) + µE(x, t)n(x, t), ( 10
)
where µ is the mobility of the mobile ions, and D is their thermal diffusion coefficient; they are related through the Einstein relation D = k B T q -1 µ. The ionic flux is related to the current through j(x, t) = qj n (x, t).

Finally, in order to introduce the caracteristic amplitude V 0 and frequency ω of the driving, we write it as

V (t) = V 0 f (ωt), (11)
where f (u) is a 2π-periodic function with absolute maximum 1.

IV.2. Dimensionless parameters and equations

The model contains many parameters; in order to understand its behavior, we introduce the following dimensionless quantities

x = x L , (12) t
= µqn ǫ t, (13)
ñ = n n , (14)
Ẽ = ǫ qnL E, (15) σ
= σ qnL , (16)
Ĩ = ǫI µq 2 n2 LA . (17)
These quantities allow to rewrite the equations as

∂ t ñ(x, t) = -∂ xj n (x, t), (18) jn
(x, t) = Ẽ(x, t)ñ(x, t) -T ∂ x ñ(x, t), (19)
∂ x Ẽ(x, t) = ñ(x, t) -1, (20)
Ẽ(0, t) = Ẽ(1, t) = -σ( t), (21)
Ṽ (t) = σ( t) + 1 0 x ñ(x, t) -1 dx, (22)
Ĩ( t) = d Ṽ d t ( t) - 1 0 jn (x, t)dx, (23)
Ṽ ( t) = Ṽ0 f (ω t). (24)
These equations depend only on three dimensionless parameters, which determine the behavior of the system:

Ṽ0 = ǫV 0 qnL 2 , (25)
T = ǫk B T q 2 nL 2 , (26) ω
= ǫω µqn . (27)
Since the parameters n, µ and ǫ are not known, the dimensionless parameters cannot be computed directly. However, the ratio T / Ṽ0 can be estimated for the experiment corresponding to Fig. 1 at room temperature (300 K), assuming that the mobile ions are oxygen ions carrying two elementary charges, and a 10 V voltage drop:

T Ṽ0 = k B T qV 0 ≃ 10 -3 . (28)
Furthermore, the ohmic behavior observed at high frequency in Fig. 1 can be exploited. The order of magnitude of the capacitive current in Eq. ( 7) is

I cap ∼ AǫωV 0 L , (29)
and the ionic current may be estimated by assuming that the charge displacement is small, so that the effect on the electric field is weak, hence

I ion ∼ AµqnV 0 L . (30)
In the experiment corresponding to Fig. 1, a triangular signal is imposed, so that the capacitive current, which is proportional to the derivative of the imposed voltage, should be constant. This constant is not observed, meaning that the capacitive current is much smaller than the ionic current, I cap ≪ I ion . Hence, the measured current is the ionic current, which allows to estimate the ionic conductivity

µqn = LI ion AV 0 ≃ 10 -4 S cm -1 . ( 31
)
We can check the consistency of our relation by evaluating the ratio of the currents at the highest frequency, which corresponds to the dimensionless frequency:

ω = I cap I ion = ǫω qµn ≃ 5 10 -10 ǫ r , (32)
where ǫ r is the relative permittivity of the medium. Note that this permittivity is evaluated without the effect of the mobile charges, and should thus be of order one, so that the parameter ω should be very small.

IV.3. Numerical integration

We now turn to the numerical integration of the dimensionless equations (18-24). We define the value of the density on a lattice and use a finite differences scheme that conserves the number of ions. The ionic current Ĩion is plotted as a function of the voltage drop Ṽ for a triangular signal with Ṽ0 = 100, T = 0.1 and different frequencies in Fig. 4. We note the close resemblance with Figs. 1 and2. The simulations reproduce several features of the experimental findings, such as the non monotonous variation of the cycle area with frequency and the existence of relaxation currents.

In order to get insight onto the physical mechanism responsible for the current drop during the cycles, we plot the charge distribution at different points of the cycle for ω = 50 in Fig. 5. We observe that during a cycle, the mobile charges move as a whole from the left electrode to the right one, and spread slightly before hitting the right electrode. The current drop indeed occurs when the charges meet the right boundary of the system.

IV.4. Discussion of the model

Because of the diffusion force coming from the ionic concentration gradient, the ions motion is getting slower and slower. This regime is nothing but the appearance of the Warburg diffusion 23,24 , caused by the ionic accumulation and already observed in the low-frequency behavior of the dielectric constant 14 . Before the mobile ions begin to feel the edges of the sample they are in a free motion regime with a conductivity of 10 -4 S cm -1 related to the ionic conductivity 14 .

However, important discrepancies remain between the model and the experimental results. First, the value of the dimensionless parameter ω used in the model, 10 ≤ ω ≤ 300, is incompatible with the experimental estimation for a reasonable value of the relative permittivity of the medium, Eq. (32). Second, we plot only the ionic current in the model, and discard the capacitive current. Yet, for the values of the dimensionless parameters used, the capacitive current is much larger than the ionic current. A possible explanation could be a bad choice of dimensionless parameters, but we could not recover the experimental behavior for other values of the parameters. A key parameter is Ṽ0 , which compares the voltage drop to the mobile charges available in the system. The current drop at low frequency is observed only when Ṽ0 ≫ 1, when the voltage drop is able to drive all the mobile ions to a boundary of the system. We conclude that this model contains the basic ingredients to explain the experimental behavior, but misses some important physical mechanism. Identifying such mechanism remains an open question and will certainly provide further insight onto the nature of ionic migration in such system.

The fact that the system is found to not strictly obey the first criterion for memristors 16 in the sense that the I-V curves are not strictly pinched in zero, is clearly captured by the model (Fig. 4). Indeed, two mechanisms homogenize the ionic distribution in the absence of electric field: the electrostatic interaction between the ions and the uniform background charge, and the thermal diffusion of the individual ions; in the model, the ratio of these two effects is given by T . Due to these effects, the memory of the material has a volatile character, making RTO material a good candidate for the applications involving volatile memory components.

Further work has to be carried out in order to elucidate the conditions of creation of non-volatile (or less volatile) memory effects, that are also observed. This involves a large number of parameters such as annealing conditions, pressure and temperature, nature of the electrodes, quality of the surface contact and will be the object of future studies. For instance, one may consider using porous electrodes with ionic adsorption or intercalation processes. In this case non volatile memory could be stabilized, opening the way to another range of applications.

V. CONCLUSION

In this work we have demonstrated experimentally and confirmed by simulations that the important ionic electromigration dynamics prevailing in bulk Rb 2 Ti 2 O 5 at room temperature leads to non-transverse memristive properties. Because of the finite dimensions on which the ionic migration occurs, two important behaviors exist: one regime where the ions are freely migrating and another one where the ions are accumulating at the edges of the system, which yields a hysteretic current drop, giving rise to the memristive properties. Electrostatic interactions and thermal diffusion both contribute to the volatility of the memory effects. A microscopic model based on ionic migration and diffusion of one-type charged species is able to reproduce the observed nontransverse I-V curves but points to a missing physical phenomenon that still remains to be identified.

BIBLIOGRAPHY

Figure 1 .

 1 Figure1. I-V curves for Rb2Ti2O5 single crystal (sample A) under room conditions for frequencies ranging from 1 mHz to 500 mHz. The sample is a single crystal which dimensions are 1000 µm × 300 µm × 300 µm. This crystal demonstrates genuine memristive behaviour, with an almost ideal "memristor" character according to the criteria defined by Adhikari et al.16 (see text).

Figure 2 .

 2 Figure 2. I-V curves acquired at 1.2 mHz on Rb2Ti2O5 (sample B) under room atmosphere. The sample is a polycrystal whose dimensions are 1500 µm × 700 µm × 800 µm. There are eight I-V cycles that appear superimposed demonstrating the excellent reproducibility of the measurements.

Figure 3 .

 3 Figure 3. Evolution of the current measured on a Rb2Ti2O5 sample at room conditions as function of applied voltage. The blue curve is for the triangular applied voltage and the red curve is the measured current. The crystal dimensions are 3500 µm ×400 µm × 1500 µm.

Figure 4 .

 4 Figure 4. Ionic current versus voltage drop in the Nernst-Planck-Poisson model Eqs. (18-24) with a triangular driving for Ṽ0 = 100, T = 0.1 and different values of ω given in the legend. The circles indicate the positions where the charge distribution is plotted in Fig. 5 for ω = 50.

Figure 5 .

 5 Figure 5. Mobile charge distribution in the sample for Ṽ0 = 100, T = 0.1 and ω = 50 along the cycle. The values of Ṽ are indicated in the legend, and correspond to the circles in Fig. 4.