A Gerbi

C González

R Buzio

N Manca

D Marrè

L D Bell

D G Trabada

S Di Matteo

P L De Andres

F Flores

Accurate ab initio determination of Ballistic Electron Emission Spectroscopy: application to Au/Ge

Keywords: numbers: 73.23.Ad, 73.30+y, 73.40.-c, 68.37.-d Condensed Matter and Materials Physics. Research Areas: Transport, nano, interfaces. Physical Systems: Au, Ge, metal-semiconductor. Properties: Schottky barrier. Techniques: Ballistic Electron Microscopy

Ab initio non-equilibrium Keldysh formalism based on an N-order renormalisation technique is used to compute I(V ) Ballistic Electron Emission Microscopy characteristics at the Au/Ge(001) interface. Such a formalism quantitatively reproduces precise experimental measurements under ultra-high vacuum and low-temperature conditions. At T = 0 K the ballistic current follows the law (V -V SB ) 2.1 , V SB being the Schottky barrier. At T > 0 K, temperature effects become significant near the onset and must be taken into account to identify an accurate value for V SB from a best-fit procedure. We find two values for V SB : 0.67 and 0.75 eV, which we associate with two different atomic registries at the interface.

I. INTRODUCTION

Scanning Tunneling Microscopy (STM) is a technique that has revolutionised the characterisation of nanostructures on surfaces due to its atomic resolution and its ability to manipulate atoms and molecules on a surface. 1 A significant extension of STM is Ballistic Electron Emission Microscopy (BEEM), that was introduced by Bell and Kaiser to study buried objects like the metal-semiconductor interface (MSI). 2 In particular, BEEM is an invaluable tool with potential to characterize with nanometric resolution the all-important technological system formed by the metal-semiconductor interface, as it has been experimentally demonstrated 3,4 and, theoretically rationalized 5 . Indeed, the MSI is the basis for rectifying devices widely found in electronics; the knowledge of their characteristics with nanoscopic resolution is of primary importance as their typical sizes diminish towards the nanoscale. Therefore, it is clear the importance of the determination of V SB , which can only be achieved with precision by a suitable combination of carefully controlled experiments and accurate theoretical simulations.

However, a significant drawback to faithfully extract information from BEEM has been the lack of a comprehensive first-principles theory able to describe with internal consistency the whole tip-metal-semiconductor system in a BEEM setup. 6 Here, we introduce an ab initio formalism based on non-equilibrium Keldysh's Greens functions, that consistently includes multiple-scattering and coherence along all the steps in the BEEM experiment. Furthermore, since it only involves N-order calculations, it allows the study of large enough threedimensional atomistic models that can be realistically compared to experiments. 7 Therefore, we combine our own ultra-high vacuum (UHV), low-temperature BEEM experiments for the Au(001)/Ge(001) interface with the ab initio theory presented below to extract accurate values for the Schottky barrier.

Good agreement between theory and experiments requires two slightly different values for V SB : 0.67 ± 0.01 and 0.75 ± 0.01 eV. We attribute such a double barrier to the existence of different atomic registries at the interface (details are given in the Appendix A). Our calculations for alignment of levels across the interface for two Density Functional Theory (DFT) optimized models predict a difference in their respective V SB of 0.10 ± 0.01 eV, which is comparable to the one obtained from our best-fit to experiments. The proposed procedure yields a remarkable accuracy in the determination of V SB (≈ ±k B T eV) due to the ab initio superlayer of Ge(100). We have considered the BB configuration, with a tip-sample distance of 5 Å, V SB = 0.649 V, η = 0.005 V, and 9216 k-points in the IBZ. The dotted blue line is the best-fit using Eq. ( 2) with µ = 2.1 and α = 6.5 pA V -µ . Lower inset shows the same plot in log scale.

Upper inset: bands alignment between tip-metal-semiconductor (notice the strong inversion).

II. EXPERIMENTAL: AU/GE

Because its higher electron and hole mobility relative to silicon, Ge has received considerable attention for high-speed metal-oxide-semiconductor devices (e.g. MOSFETs). [8][9][10][11][12] However, a notable disadvantage that has so far limited performance in Ge MOSFET nchannel devices is the strong Fermi Level Pinning at the charge neutrality level (CNL), [13][14][15] creating a high Schottky barrier and preventing the formation of low resistance contacts. 9,16 Curiously enough, no BEEM data have been reported so far in literature for the metal/Ge interface, although some works have been published on buried Ge dots, and Si 1-x Ge x strained interfaces 17,18 .

In our experiments, a Ge(100) wafer (n-type, Sb-doped, 3.97-4.46 Ωcm, MTI Corporation) was cut into pieces and cleaned to remove the native oxide and to obtain a hydrogenterminated Ge(100) surface. 12 The cleaned Ge pieces were loaded within few minutes into the UHV deposition chamber for the Au contact fabrication, obtained by thermal evaporation through a shadow mask (nominal thickness and area of the metal contact are ≈ 20 nm and 2.1 mm 2 , respectively). The Ohmic back contact was fabricated by depositing a thick Al film by pulsed laser deposition from a high purity target. 19 Interpolation of room-temperature current-voltage characteristics using conventional thermionic emission theory provided an effective Schottky barrier height V SB = 0.62 eV and an ideality factor n = 1.04, in line with previous reports. 20 BEEM measurements have been performed under UHV with a base pressure of 3 × 10 -10 mbar using a modified commercial STM equipped with an additional low-noise variable-gain current amplifier. [21][22][23] To improve the signal-to-noise ratio, data was taken under dark, with tunneling current I T in the range 0.5 to 5 nA, and at T = 80 K. For the acquisition of each BEEM spectrum, the tip voltage V was ramped under feedback control, keeping the tunneling current I T constant. A typical spectrum consisted of 100-3600 averaged curves, and each one was acquired in ≈ 1 min to 10 s. Individual BEEM spectra were acquired from a grid on an area of lateral size at least 400 × 400 nm 2 , with spacing between two neighbor spectra of at least 10 nm (to reduce local damage induced by the electron beam). For a typical Au thickness of about 10 nm and propagation perpendicular to a {001} plane we estimate that the BEEM electronic spot covers an area of about 10×10 nm 2 at the interface; such a region is not uniformly illuminated, it displays lines with a typical width of ≈ 1 nm due to diffraction-like elastic scattering from electronic bands, e.g. see Fig. 3 in 24 or Fig. 9 in 6 . Although the lines are clearly thin enough to yield nanometric resolution at the interface while scanning a step, a grain boundary or even a point defect, the existence of symmetry related lines separated by about 10 nm add contributions from SBHs of regions apart by the same distance. Therefore even a single I(V) curve should show the signature of both domains for an interface made of two domains separated by ≈ 10 nm with different SBHs.

Such multi-contribution from several domains is reinforced when I(V)s taken in points on the metal surface separated by more than 10 nm are averaged together to increase the signalto-noise ratio, as it has been made in this work. Noise current fluctuations in individual raw spectra amounted to ≈ 12 fA rms. Such a low noise level enables an accurate comparison between theory and experiment. The many different spectra acquired on different locations of the gold electrode indicated a high spatial homogeneity for the Schottky barrier height distribution (i.e., good reproducibility). Typical spectra showed a threshold around 0.66 eV. However, interpolation with simplified BEEM models for a single Schottky barrier was not satisfying in some respects (see below). In what follows, we consider the BEEM-current for electrons only, because Ge is n-doped and has an inversion layer (the Fermi level is very close to the valence band edge), making the transport by holes difficult to analyze without including the long semiconductor depletion region in the theory.

III. THEORY

To obtain an accurate value for V SB it is necessary to resort to first-principles theory to get the best possible fit to experimental data, avoiding unnecessary and unwarranted assumptions. Previous theoretical approaches have addressed the BEEM current as a sequence of steps for (i) tunnelling, (ii) propagation through the metal base, (iii) transmission at the metal-semiconductor interface, and (iv) propagation in the semiconductor collector. 5,25 Each of the steps has been treated to different sophistication levels, sometimes mixing quantummechanics in some places with semi-classical approaches in others.

Here we introduce a coherent approach based in Eq. (1) (below) that improves over previous calculations by joining all processes in a single self-consistent step, therefore increasing our ability to determine an accurate theoretical I(V) characteristics that could be compared with experimental data to extract the much sought value of V SB . Indeed, an approach breaking down separately the processes contributing to BEEM could not determine the Schottky barrier with better accuracy than the model (i) used to describe the transmission through the interface and (ii) the currents coming in and out from that region. In particular, to reproduce our experimental data we build a 3D atomistic model consisting of a Au(100) slab made of 50 superlayers (cf. Appendix B for an explanation of the term superlayer) That region, which includes the metal-semiconductor interface, is grown on a semi-infinite Ge(100) surface (collector), which is computed using the same scheme but keeping bulk-like values for the interaction between layers. Therefore, we incorporate a detailed ab initio description of the atomic positions and the charge redistribution at the interface region, and the currents in that region are computed by bundling everything together in Eq. ( 1). Such a specific description of the system allows us to validate a particular power-law for the given interface under consideration. Figure 1-(a) shows a scheme for our BEEM model and details of the different Green's functions involved in the calculation are given in Appendix A.

consisting
For electrons propagating with energies between E and E + dE at T = 0K (assuming no exchange of energy with the medium), the ballistic current injected from the STM tip into the n + 1 superlayer (i.e., p semiconductor superlayers below the MSI; a place where the system can be safely matched to the semi-infinite semiconductor acting as collector) is

given by, 7 m+1,m +1 are retarded and advanced Green's functions of the semiconductor computed using a decimation technique (cf. Fig 1). 26 On the other hand,

∂I ∂E = 4e IBZ d 2 k (l,l )<n+1 Tr D R l,1 Γ 1,1 D A 1,l Γ l ,l (1)
D R l,1 = δ l,l -G R l,l t S l ,l +1 g R (S) l +1,l +1 t S l +1,l -1 G R l ,1 D A 1,l = G A 1,m δ m ,l -t S m ,m +1 g A (S) m +1,l +1 t S l +1,l G A l ,l -1 G R,A l,l are
Γ l ,l = t l ,n+1 ρ (S) n+1,n+1 t n+1,l Γ 1,1 = t 1,0 ρ (T ) 0,0 t 0,1
are injection rate matrices from the metal to the semi-infinite semiconductor region (S), and from the tip (T) to the metal. We notice that injection of tunneling electrons can be treated using similar techniques. Therefore, we take advantage of introducing the Green's functions G that describe accurately the metal base and a few semiconductor superlayers (p = 4) that correspond to the geometrical configuration we have independently determined using ab-initio plane-waves DFT (cf. Appendix A). Finally, the resulting Green's functions are coupled to the semi-infinite semiconductor via the density of states ρ (S) n+1,n+1 in Eq. ( 1) that is calculated using the decimation technique. The sheer advantage of using N-order techniques is that the computational effort scales linearly with the size of the system, allowing us to set up models with a realistic enough size so they can be safely compared with experiments (i.e., including hundreds of layers or more).

We notice that even if the ballistic current above is actually a parameter-free expression, it is usual to add an small imaginary part to the energy (E → E + iη) to: (1) smooth out peaks, substituting the k integration over the 2D interface Brillouine zone (IBZ) by a summation on a set of special k-points, and (2) ensure convergence of outgoing waves.

These calculations have been programmed so they can use as input either ab initio bulk bands described by tight-binding parameters, [START_REF] Papaconstantopoulos | Handbook of the Band Structure of Elemental Solids[END_REF] or electronic structure bands calculated with a local combination atomic orbitals (LCAO) DFT-code. [START_REF] Lewis | [END_REF] We have checked both and we have found only minute differences when the agreement between both approaches is good.

The former approach is probably simpler to use, while the latter would be necessary when dealing with materials not listed in the literature or if a complete description of the atomic configuration and charge distributions is required. For our case of study all the results in this work have been obtained relying on Ge-bands in [START_REF] Papaconstantopoulos | Handbook of the Band Structure of Elemental Solids[END_REF] while the Au-bands as well as the interface are parametrized using the LCAO-DFT data.

Since the tip-metal distance in a BEEM experiment is large and the tip geometry is sharp, the use of a single apex with an s-orbital (ρ T 0,0 ) is justified. On the other hand, several atoms on the periodic surface can contribute to the injection, possibly bringing about coherence effects which are fully accounted for by this formalism. In fact, because k is conserved at the ideal MSI, the initial energy and momentum distribution at the surface influences the BEEM current and permits to study with submolecular resolution organic molecules deposited at the surface. 29,30 Finally, we remark the instrumental advantage of dealing with an N-order method to calculate the Green's functions, which becomes crucial to obtain the BEEM current on a realistic MSI. 31 The atomic configuration at the experimental interface is not currently available from structural techniques. Therefore, we have optimized plausible models using DFT. 32 Fig. 2 shows the result of using Eq. ( 1) to compute the ballistic current for the BB Au(001)/Ge(001) interface (black continuous line). Such a plot immediately suggests that the available phase space, ∂I ∂E ∝ (E -V SB ) µ-1 , plays a significant role (µ takes a value that only depends on details of the system like the band structure of the semiconductor or the atomic arrangements at the interface). Therefore, the ballistic current in a small interval of energies around V SB can be easily obtained from the differential conductance (blue dashed curve in Fig. 2):

I B (V -V SB ) ∝ V V SB (E -V SB ) µ-1 dE ∝ (V -V SB ) µ ; V SB ≤ V (2)
Such a free-electron model has been proposed in the literature with µ = 2.0, 2 or µ = 2.5, 33 depending on different hypothesis on the transmission coefficient at the MSI, cf. Appendix C.

It allows a most simple procedure to get a value for V SB by a best-fit to experimental data, with an error bar of ≈ ±0.05 eV which is mainly originated in the different values available for µ. 34 Therefore, to improve the accuracy on the determination for V SB it is desirable to know a better value for µ corresponding to a given metal-semiconductor combination;

we obtain its specific value for Au(001)/Ge(001) by comparing our ab initio Eq. ( 1) to the simplified model of Eq. ( 2). We find (Fig. 2) that µ = 2.1 ± 0.1 yields a best-fit (the incertitude is associated to different tip-sample distances). It is interesting to notice that this implies for our particular Au/Ge system a quantum-mechanical transmission at the interface quite constant, with an smoother energy dependence than the square root associated with a square barrier, cf. Appendix C.

To determine V SB with an accuracy of ≈ ±0.01 eV it is also necessary to care at least about two non-ballistic contributions: inelastic losses and thermal effects. By limiting the energy interval to | V = V -V SB |< 0.2 eV around the onset secondary electrons are mostly removed since their probability to be reinjected is low. Therefore, they can be accounted for by an approximately constant factor describing attenuation, σ = e -l λ , where l represents a typical length of the system and λ is the mean free path related to the optical potential due to inelastic interactions, λ(E) = v(E) 2η(E) . [START_REF] Pines | Elementary Excitations In Solids[END_REF][START_REF] Pendry | Low Energy Electron Diffraction[END_REF][START_REF] Reuter | [END_REF][38] We notice that the energy dependence of the optical potential is considerably smoothed if | V |< 0.2 eV, which allows us to take a constant value for η. Under these conditions, the limiting value η = 0 is simply recovered by multiplying the current in Eq. ( 1) by the inverse of σ. Since we use in the IBZ around 10000 special points to cover a band of approximately 1 eV, a typical value for η would be 1 100 eV. In practice, we use η = 0.005 eV, that corresponds to an attenuation length of λ ≈ 1500 Å and is much larger than the metallic base (l ≈ 200 Å), which is another reason to use a constant value within ∆V < 0.2 eV from the onset. Finally, η = 0.005 eV ensures both good convergence with our choice for special points and describes well the damping of the ballistic current in this small energy range.

Thermal effects, on the other hand, are most easily incorporated by taking into account the Fermi-Dirac distribution, f (E, T ). Then, Eq. ( 2) can be re-written:

I B (V -V SB , T ) = α ∞ 0 (E -V SB ) µ-1 [f (E -V ) -f (E)] dE = (3) = -α (k B T ) µ Γ[µ] Li µ -e V -V SB k B T (4)
where α = I EXP (0, T ) is a free parameter to adjust the absolute value of the model to the experiment at a single energy value (V = V SB has been used here), Γ is the Legendre's Gamma function, and Li µ (z) is the polylogarithm function. [START_REF] Weisstein | Polylogarithm[END_REF] While this function is amenable enough to be used around V = V -V SB = 0 to obtain best-fits like the ones given in Fig. 3, an even simpler expression based on a power expansion is provided in Appendix C. We remark that thermal effects cannot be neglected for | V | 2k B T , a region significant for the sake of a better determination of V SB . In addition, it is interesting to notice the exponentially decaying contribution of electrons excited by thermal effects over the Schottky barrier in the sub-threshold region (V < 0). [START_REF] Qin | [END_REF] Naturally, this contribution can only be used above the background noise level, which in our experiments appears as a fluctuating value around 10 -6 < I B I T < 10 -7 .

IV. RESULTS AND DISCUSSION

A simple procedure to get a value for the Schottky barrier would be to locate the separation between the two regimes in Eq. (C2) of Appendix C, e.g. as signaled by a change from a constant value in the derivative of d log I dV (sub-threshold exponential law) to a decreasing value (above-threshold power law). In the inset of the lower panel of Fig. 3 we plot that value for the experimental data (black line) to be compared with the theoretical result from Eq. (3) -red line-. We conclude that the sub-threshold regime in our experiments displays a noise level that makes difficult to identify the onset from the break between both regimes. Therefore, we concentrate on extracting the Schottky barrier by utilizing mostly the above-threshold region (Eq. (C2) of Appendix C). In the upper panel of Fig. 3 we show a fit to experimental values (magenta stars) using Eq. ( 3) with µ = 2.1 (black lines).

The slope of I(V ) near V ≈ 0 fixes T ≈ 110 K. The significance of these thermal effects can be seen by comparing with the value T = 0 K (blue line). Our determination of the effective temperature predicts a slightly larger value than the nominal experimental one; for comparison, the green line shows that T = 80 K yields a slightly poorer fit, indicating that we cannot determine the effective temperature operating on the ballistic electron beam with an incertitude better than a few tens of K. However, one expects the ballistic current to induce slightly larger local temperatures within the electron beams due to inelastic effects.

The log-log plot in the lower panel of Fig. 3 makes clear the necessity to include a second onset near V ≈ 0.07 eV to explain the experimental data,

I EXP B = I B (V -V SB 1 , T ) + I B (V -V SB 2 , T
). Therefore, our best fit includes two Schottky barriers located at 0.67 and 0.75 eV, in the relative proportion of 1:5 (dotted and dashed black lines respectively in the upper panel of Fig. 3). The first Schottky barrier helps to fit well the experimental data in the interval [0, 0.05] eV (upper inset in Fig. 3), including the initial slope expected from the thermal contribution. The second Schottky barrier provides an excellent fit to the experimental data beyond that region, up to ≈ 0.2 eV from the onset.

According to our DFT calculations, the reason for the two Schottky barriers can be physically assigned to the presence of two different atomic configurations coexisting at the interface (BB and TH). Our calculations find a difference between the electrostatic dipole in these two geometries that amounts to ≈ 0.1 eV, which is in good agreement with the value obtained by fitting the experimental data to Eq. (3). 14

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ 0.
I B (V',T)/I T ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ 0.01 0.02 0.05 0.1 0.2 V'(V) 10 T 10 5 10 R 10 Q I B (V',T)/I T ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★

V. CONCLUSIONS

In conclusion, we have developed an N-order parameter-free non-equilibrium Keldysh

Green's functions formalism that accurately reproduces the BEEM current normalized to the injected tunnelling intensity. Such a theoretical methodology can be applied to deduce a power-law approximation for the Au(001)/Ge(001) interface, As a method capable of finding reliable geometry configurations we utilize DFT Therefore, we explore the two possible high-symmetry models: either the two Au atoms in the surface unit cell take positions over both bridge positions in the Ge surface unit cell(BB), or one is located atop and the other one over the hollow (TH), cf. Fig. 1-(b).

I B ∝ (V -V SB ) 2.
Notice that these two models are simply related by a two-dimensional registry shift of ( 1 2 , 0) (fractional coordinates), but imply different relaxations and corrugations of layers. The BB configuration is preferred over the TH one by -0.64 eV, which is small enough to allow the coexistence of patches of both structures. In the BB geometry the closest Au-Ge is 2.75 Å, while in the TH configuration is 2.43 Å. The two Au closest layers to Ge display corrugations of 0.3 and 0.01 Å in the BB configuration; for the TH configuration corrugations are below 0.02 Å for all layers. These calculations have been performed using the CASTEP program with norm-conserving pseudopotentials (*OP.recpot), the Perdew-Burke-Ernzerhof GGA-XC 44 , a 610 eV energy cutoff, and a 6 × 6 × 1 Monkhorst-Pack grid 45 . Threshold values for convergence were 10 -7 eV for the total energy at each self-consistency cycle, 10 -6 eV for the total energy after each geometry optimization, a maximum force of 0.01 eV/ Å, and a maximum displacement of any atom in the optimization cycle of 5 × 10 -4 Å.

Once the interface has been accurately relaxed with the plane wave code, we have exported the atomic reconstruction to the LCAO code in order to obtain the tight-binding parameters required by our formalism. We improve the electronic distribution in the first four layers of each material by adding new layers at both sides in the unit cell. The resulting slab is formed by an interface of 42 atoms: 24 Au atoms in 12 layers, and 16 Ge atoms (one per layer) plus two final H atoms saturating the Ge dangling bonds at the end of the slab. We use an spd basis for both Au and Ge that provide a correct band structure in the bulk comparing to the results shown in previous very accurate calculations [START_REF] Papaconstantopoulos | Handbook of the Band Structure of Elemental Solids[END_REF] equation reduces to Eq. ( 1), where two kind of Green's functions are needed: g

(S)
n+1,n+1 to describe propagation through the semiinfinite semiconductor acting as collector, plus G l,m to describe the propagation through the metal base plus a few semiconductor superlayers to better describe the metal-semiconductor interface. The former is obtained from an efficient decimation technique well described in the literature. 26 The latter are computed iteratively solving Dyson's equation as explained in what follows. 7 The Green's functions, g l,l , for a system composed of l superlayers are obtained from the Green's functions for a system composed of l -1 superlayers and the interaction with a single non-interacting superlayer, t l-1,l :      g l,l = δ l,l -g 0 l,l t l ,l-1 g 0 l-1,l -1 t l -1,l -1 g 0 l ,l g l,1 = g l,l t l ,l-1 g 0

l-1,1 (B2)
This equation is solved iteratively starting from the diagonal resolvent for a single isolated superlayer, g 0 l,l = E -H 0 l,l -1 , to construct successively the pairs (g 21 , g 22 ), (g 31 , g 33 ), etc, to finally get (g l,1 , g l,l ) = (G l,1 , G l,l ).

FIG. 1 .FIG. 2 .

 12 FIG. 1. (color online) (a) Schematic BEEM setup: carriers are injected from a STM tip (green, layer n = 0) into a metal slab (yellow, from superlayer 1 to superlayer m) adsorbed on a thin semiconductor slab representing the interface region (blue, from superlayer m + 1 to n, p = n -m superlayers) that is connected to the semi-infinite semiconductor acting as collector (blue, from superlayer n + 1 to ∞). Green's functions G l,m are computed adding one superlayer at a time recursively, while g (S) n+1,n+1 are computed via a decimation technique (cf. Appendix B). (b) Two √ 2 × √ 2R 45 • optimized interfaces ( = 0.016). BB: Au on both bridges of Ge(100). TH: Au on atop and hollow sites.

× √ 2 R 45 •

 245 of bilayers of Au(001) planes with AB stacking and √ 2 2D unit cell (i.e. 200 atoms in the slab unit cell) in close contact to 4 superlayers made of Ge(100) planes with ABCD stacking sequence and 1 × 1 2D unit cell (i.e. 16 atoms in the slab unit cell).

  Green's functions taking care of propagation and multiple-scattering through the m superlayers of the metal base plus p superlayers of the semiconductor. These functions are computed adding one superlayer at a time using an N-order recursive technique based on Dyson's equation as explained in Appendix B. g R,A(S)

  We construct candidate interfaces from 1 × 1 surfaces: to compensate for Ge and Au 2D different lattice constants, we rotate one against other by 45 • to minimize interfacial strain. The resulting √ 2 × √ 2 R 45 • Au(100)/Ge(100) displays a = 1.9% mismatch between both lattices across the interface; the remaining stress can be reasonably accommodated by de-fects, in particular by different stackings. DFT optimization yields two competing local minima with an enthalpy difference of less than 0.1 eV per atom in the interface, supporting the appearance of both patches in the microscopic interface in the experiment. In the first configuration (frontal view in the upper part of Fig. 1-(b), both Au atoms in the 2D unit cell are located over the two bridges (named BB, the preferred one). In the second one (labelled TH), one Au atom is located atop, while the other is on hollow, cf. the lower part in Fig. 1-(b). The red squares indicate the size of the unit cell in the DFT simulation. These geometries have been fully relaxed to analyze the electronic structure of the MSI (details in Appendix A).

  05 0.10 0.15 0.20 0.25 0.30 V

FIG. 3 .

 3 FIG. 3. (color online) Upper panel. Stars (magenta): experimental values. Black continuous line is the best fit to experimental values using Eq. (3) (µ = 2.1), assuming T = 110 K and two Schottky barriers at 0.67 and 0.75 V (dotted and dashed black lines respectively). Inset: close-up near the region V ≈ 0 V (log plot). Bottom panel. Log-log plot showing the appearance of the second barrier. The green curve is T = 80 K, and the blue one is T = 0 K. Inset: Derivative of experimental Log(I) (black), compared to the model based in Eq. (3) -red-; notice the exp vs power-law behavior and the fluctuations for V < 0.

1

 1 , that can be flexibly used to determine an accurate value for the Schottky barrier at the MSI.For V ≈ V SB temperature effects become significant, and they should not be ignored. A detailed comparison with experimental data for Au/Ge shows the existence of two different Schottky barriers at 0.67 and 0.75 eV that have been fitted with approximately a ±0.01 accuracy. These values agree within the estimated error bars with the reported gap of Ge of 0.74 eV at low temperature,41 and can be associated with the existence of different registries at the interface. The procedure takes into account the effect of secondary electrons on the attenuation of currents due to inelastic e-e interactions. We have chosen to work in the interval | V -V SB |< 0.2 eV because for larger biases, secondary electrons start to play a significant role, and because we can neglect other effects, like electron back-scattering in the semiconductor associated with the electron-phonon interaction.Acknowledgments. This work has been funded by MAT2017-85089-C2-1-R, MAT2017-88258-R, Maria de Maeztu Programme (MINECO, Spain), MIUR Progetto Premiale 2012 EOS: organic electronics for advanced research instrumentation (Italy), and ERC-2013-SYG-610236 (EU). A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contact with the National Aeronautics and Space Administration.

  . The corresponding cutoff radii for Au (Ge) are: r c [s] = 4.5 (5.7) a.u., r c [p] = 4.9 (5.7) a.u. and r c [d] = 4.3 (5.7) a.u. and r c [H -s] = 3.8 a.u. for hydrogen. The first surface Brillouin Zone is sampled with 9216 k-points (the same used in the BEEM simulation), allowing an accurate descritption of the band gap at the bulk.

FIG. 4 .

 4 FIG. 4. (color online) 3D view for Ge(001)/ √ 2 × √ 2R45 • -Au(001) models. Left panel: Top-Hollow (tH). Right panel: Bridge-Bridge (BB). To increase visibility only two Au-layers and five Ge-layers are shown.

Forg 22 = ( 1 -g 0 22 t 21 g 0 11 t 12 ) -1 g 0 22 g 21 = g 22 t 21 g 0 11 ( 47 -FIG. 5 .

 22112222111475 FIG. 5. (color online) Eq. (4) -blue-vs the power-expansion Eq. (C2) -red dashed-for T = 110 K. Notice the change from the sub-threshold exponential regime to the above-threshold potential regime. Black dots are experimental values.

Appendix B: Green's functions transport formalism

The differential current of carriers with energy between E and E + d E injected from an STM tip into a semi-infinite semiconductor acting as a collector can be written as:

where G +-are Keldysh's Greens functions 46 and t i,j are hopping matrices describing the interaction between nearby superlayers, with j (≥ n + 1) located in the semi-infinite semiconductor acting as collector and i (≤ n) is located in the region of the metal base plus a few superlayers of semiconductor used to describe the details of the interface, cf. Fig. 1. This

Appendix C: Ballistic Model

Eq. ( 2) can be interpreted in physical terms for ballistic electrons by integrating the flux of carriers perpendicular to the interface,

where µ takes a value that depends on the model for the transmission coefficient, T , e.g., 2,33

Therefore, Eq. (C1) shows how the available phase space for injection of carriers (including a transmission coefficient) determines the purely ballistic I(V) characteristics near the onset, as written in Eq. ( 2).

Finally, we notice that Eq. ( 4) can be expanded around V = V -V SB ≈ 0 for µ > 2 to yield an approximate expression including thermal effects that only involve powers of V , the exponential function, and the Gamma (Γ) and Riemann's Zeta (ζ) functions. In the simplest case we have for V < 0,

Where the arbitrary constant 2 µ -2 2 µ ζ(µ) has been included so both approximating functions, Eqs. (C2) and (C3), coincide with Eq. ( 4), e.g. at V = 0 (cf. Fig 5).

On the other hand, for V > 0 we find the following expansion in powers of V :

Fig. 5 shows a comparison between the exact value from Eq. ( 4) and the approximations in Eqs. (C2) and (C3) for V SB = 0.67 eV and T = 110 K. We remark that for T = 0 K, in the region V ≈ 0 + the BEEM current acquires a quasi-linear dependence on the voltage that cannot be neglected if an accurate determination of V SB is sought. * andrea.gerbi@spin.cnr.it, cesar.gonzalez@uam.es


Guilherme De

Oliveira Reis

Thomas Gibaud

Brice Saint-Michel

Sébastien Manneville

Mathieu Leocmach

Laurent Vaysse

Frédéric Bonfils

Christian Sanchez

Paul Menut
email: paul.menut@supagro.fr

Irreversible hardening of a colloidal gel under shear: the smart response of natural rubber latex gels

Keywords: Isotropic particles, strain stiffening/hardening, mechano-responsive gels, re-structuring

Natural rubber is obtained by processing natural rubber latex, a liquid colloidal suspension that rapidly gels after exudation from the tree. We prepared such gels by acidification, in a large range of particle volume fractions, and investigated their rheological properties. We show that natural rubber latex gels exhibit a unique behavior of irreversible strain hardening: when subjected to a large enough strain, the elastic modulus increases irreversibly. Hardening proceeds over a large range of deformations in such a way that the material maintains an elastic modulus close to, or slightly higher than the imposed shear stress. Local displacements inside the gel are investigated by ultrasound imaging coupled to oscillatory rheometry, together with a Fourier decomposition of the oscillatory response of the material during hardening. Our observations suggest that hardening is associated with irreversible local rearrangements of the fractal structure, which occur homogeneously throughout the sample.

INTRODUCTION

Natural Rubber Latex (NRL) is naturally produced by the Hevea brasiliensis tree, in laticiferous cells located in the tree bark. This colloidal suspension is released by wounded trees after tapping and rapidly solidifies to form a protective layer against further attacks. The composition of the exudate is complex [1], but the main physical properties are attributed to rubber particles accounting for more than 90% of the dry mass. Their diameters range from a few hundred nanometers to a few micrometers. Rubber particles mostly consist of a poly(cis-1,4-isoprene) rubber core surrounded by a membrane of phospholipids [2,3] and proteins [4]. Their surface is negatively charged at native pH (pH ≈ 7), which ensures net repulsion and therefore latex stability within the tree. After exudation, the pH rapidly decreases under bacterial acidification [5]; the particle surface charge density drops and the suspension is consequently destabilized. Particles aggregate in fractal-type clusters, which eventually form a gel [6].

To produce Natural Rubber (NR), NRL gels are compressed and dried so that the residual water and its soluble nonisoprene content are extracted. During this process, rubber particles lose their individuality and the material is described as a structured network of cis-poly-isoprene chains loosely cross-linked with non-isoprene compounds (mainly lipids and proteins) [7,8]. The cross-link density of these polymeric materials is further increased during a process called vulcanization. Vulcanized NRs can sustain up to 800% deformation before rupture, and their stress-strain relationship is highly nonlinear. Above a critical strain, whose value depends on the preparation process, but is generally above 300%, vulcanized NR exhibits strain stiffening: its apparent elastic modulus increases rapidly with strain. This behavior, also known as strain hardening, is attributed to strain-induced crystallization [9], [10]. It is highly sensitive to the supramolecular organization in the natural product and is called "green strength" [11]. NR exhibits such unique thermo-mechanical properties that more than 40% of all the rubber used worldwide is, even today, of natural origin, despite fundamental drawbacks, such as its variability. Understanding the structure-properties relationship in such materials is still a challenge, and requires an understanding of how the structure is formed at the different stages of NRL processing into NR.

We focus here on the mechanical properties of NRL gels, which are soft materials in which rubber particles form a percolating network in water. As NRL suspensions are very variable due to agronomic parameters (clone, tapping system), or edaphoclimatic factors (mainly the season) [12], we work with a commercially available NRL suspension of pre-vulcanized rubber particles, which is stabilized in ammonia and far less variable. In spite of a slightly different composition of the continuous phase, such commercial suspensions behave similarly to the fresh exudate: gelation is observed under similar acidification conditions and the elastic modulus of the gels formed from commercial or native NRL are roughly the same. We induce continuous and controlled acidification of our NRL suspensions in the Couette geometry of a rheometer, in order to form gels that can be mechanically characterized in situ, with no stress history. In previous work focusing on gel formation [6], we showed that such gels display strain hardening behavior, which is examined in detail in this article.

In the following, we first describe the main features of the hardening phenomenon under oscillatory shear. We demonstrate that this behavior is irreversible and occurs over a large range of deformations. We then show that NRL gels display a fractal scaling of their elastic modulus at rest. The spring constant between particles determined from this analysis is com-patible with known values of rubber elasticity. In the third section, we show that the hardening amplitude decreases with the particle volume fraction. We further investigate the nonlinear behavior through a Fourier analysis of the strain response and show that hardening is associated with the emergence of odd harmonics, whose strength grows with the imposed stress. These observations are further discussed in the last section, where we propose a simplified picture of the structural reorganization that might lie behind this unique behavior.

MATERIALS AND METHODS

Materials

Natural rubber latex was obtained from Dalbe, France (ref. 4770002). The volume-based diameter 2a of the prevulcanized particles is 1 µm, as measured previously by static light scattering [6]. To remove ammonia, suspensions were first diluted with an equivalent amount of Tris-HCl buffer (pH 8.5, ionic strength = 7 mM). TRIS (purity ≥99%, CAS 77-86-1) and HCl (purity 37%, ≤1 ppm free chlorine) were obtained from Merck, Germany.

Then, they were dialyzed for 7 days against the same buffer using a dialysis tube (SpectraPor, cut-off 12-14kDa, diameter 29 mm), before dilution with the buffer to reach the desired volume fraction φ v . Suspensions with φ v > 0.3 were obtained by direct dilution of the commercial suspension with the buffer, without dialysis, which does not impact the final properties as will be shown in this paper. All the prepared suspensions were stored at 5 • C and used within three weeks.

Gelation

Glucono-δ-lactone (GDL, ref G2164, purity ≥99%, Sigma-Aldrich) was used to reduce the pH of the suspension uniformly and gradually. For φ v < 0.1, a concentration of 1 %wt GDL set as a standard condition was used to reach a final pH close to 4 in less than 1 h. For φ v > 0.1, the amount of GDL was increased to reproduce the standard acidification kinetics [6]. GDL was dispersed in the suspension by stirring for 1 min, before loading into the Couette cell of the rheometer.

Rheological characterization

Rheological properties were determined using three stresscontrolled rheometers: an AR2000ex (TA Instruments), equipped with a stainless-steel smooth Couette concentric cylinder geometry (inner radius 14 mm, outer radius 15 mm and height 42 mm) and a home-made stainless-steel serrated Couette geometry (same dimensions with saw tooth serrations 1 mm in depth and 1 mm base on each side) and an ARG2 (TA Instruments) equipped with a smooth PMMA (Poly(Methyl MethAcrylate)) Couette geometry (inner radius 23 mm, outer radius 25 mm and height 60 mm). Thanks to their feedback loop, these instruments can be used either in strain-or stresscontrolled mode.

Gelation was first monitored for 4 hours thanks to smallamplitude oscillatory shear at a constant frequency of 1 Hz and a constant strain amplitude of 0.5%. The results showed no measurable impact of the oscillation amplitude of 0.5% on the gelation process and on the final gel properties, particularly on the elastic modulus at rest G 0 . The gel was then characterized with a frequency sweep from 0.01 to 10 Hz at a constant strain of 0.5%, which was followed by a stress sweep at a constant frequency of 1 Hz, starting from 0.02 Pa up to the fracture of the gel.

Local rheology

We used local oscillatory rheology from echography (LORE) to reconstruct time-resolved local displacements during oscillatory stress experiments [13], [14]. In LORE, rheological measurements are synchronized with high-frequency ultrasonic echography in order to map the local strain response of the NRL gel subject to an oscillatory stress. LORE is performed with the ARG2 (TA Instruments) rheometer equipped with the PMMA Couette geometry of gap e = 2 mm. The gels are seeded with almost density matched polystyrene particles (Dynoseeds TS-20 of diameter 20.00(1) µm) to provide ultrasound contrast. Their concentration of 1% wt yields sufficient ultrasound intensity yet is small enough to prevent multiple scattering and to ensure that the polystyrene particles do not affect the mechanical properties of the sample.

In brief, our high-frequency ultrasonic imaging device relies on a linear array of piezoelectric transducers with a total active length of 32 mm [14,15]. The transducer array is immersed in a water tank surrounding the Couette cell. Short ultrasonic plane pulses with a central frequency of 15 MHz propagate across the gap. These pulses get scattered by the seeding particles within the gel and the backscattered pressure signal is recorded by the transducer array, leading to an "ultrasonic speckle" signal with 128 measurement lines in the vorticity direction z and 640 sampling points in the radial direction r. The spatial resolution along the z-direction is 250 µm and 100 µm in the r-direction with a sampling of 600 images per oscillation period [14]. For each oscillatory stress experiment, a 2D-ultrasonic scan with 20 sequences of five oscillation periods each is recorded. The analysis of the ultrasonic data gives access to the displacement maps ∆ within the Couette cell as a function of time t and of the position (r, z) across the gap, where r is the radial distance to the stator and z the position along the vertical direction. We checked that such displacement maps are always invariant along the z-direction. Hence, we focus thereafter on their average over z, ∆(r, t), from which we compute the local strain γ [14]. Under an imposed oscillatory stress, the strain response of the sample may be nonlinear in contrast to the sinusoidal stress input. This results in the presence of harmonics in the Fourier series decom-position of the strain γ(t) = k γ k cos 2πk f t + ϕ k , where γ k and ϕ k are the amplitude and phase of the k th harmonics.

RESULTS AND DISCUSSION

NRL gels exhibit irreversible hardening under shear

When submitted to some mechanical strain, Natural Rubber Latex gels harden. This is seen through oscillatory shear rheology where the material is submitted to oscillations of increasingly large amplitudes in either strain-controlled or stress-controlled experiments. When the amplitude of the imposed oscillations is small, the elastic modulus G remains independent of the oscillation amplitude. This corresponds to the linear viscoelastic regime. However, for large enough strain amplitudes, typically γ > 0.2, G steeply increases until material rupture as illustrated for φ v = 0.04 in Figure 1. This behavior clearly differs from the strain hardening of dense rubber materials, a phenomenon observed at much higher deformations (γ > 3) and associated with a transient, and reversible, local crystallization of polymer chains [9,10]. In the present hardening regime, which is observed on NRL gels in both strain-and stress-controlled modes, the strain-stress relationship is no longer linear: above the onset of hardening, for σ > σ 0 and γ > γ 0 , one has γ ∼ σ α with α ≈ 0.3 (inset of Fig. 1). Correspondingly, the elastic modulus grows as G = σ/γ ∼ γ 1/α-1 ∼ γ 2.3 at large strains as seen in Fig. 1. For the stress-controlled experiments, this corresponds to G ∼ σ 1-α ∼ σ 0.7 . Further decomposing of the strain response in a Fourier series shows that strain hardening is in line with the emergence of odd harmonics (insets of Fig. 1). Lastly, for γ > 1, the gel fails.

We conducted LORE experiments to map local displacements ∆(r, t) in the sample and thereby investigate whether hardening is localized or homogeneously distributed throughout the material. For each imposed stress amplitude, local deformation is recorded during five successive oscillations at f = 1 Hz. As shown in Fig. 1b (left) for a stress amplitude σ=1.8 Pa, corresponding to a strain amplitude γ = 0.23 at the beginning of the hardening regime, the strain response is distorted. Thus the strain response is clearly nonlinear and contains significant higher-order harmonics. Yet the local strain amplitude ∆ varies linearly across the gap (Fig. 1b, right): the bulk material undergoes homogeneous deformation. There is no sign of slip at the walls nor any significant curvature in ∆(r) that could demonstrate some degree of strain localization, e.g. through localized strain stiffening or strain thinning, over length scales larger than the LORE spatial resolution of about 100 µm. This suggests that the structural reorganization responsible for hardening is localized at scales below 100 µm and is homogeneously distributed over the entire material. In addition, once a stress oscillation is imposed, we do not observe any temporal change in the strain response over the first five oscillations. Hardening therefore occurs within the first oscillation in the bulk material.

Strain hardening has been reported for colloidal gels obtained by fractal aggregation of nanometer-sized polystyrene beads [16] or globular proteins [17,18]. Hardening has also been observed in biopolymer-based hydrogels [19][20][21]. NRL gels however exhibit two very specific properties that we now discuss in more details.

First, the hardening reported here is irreversible: if at a given strain γ max during material hardening under an upward strain sweep, one decides to decrease the strain amplitude, the material remains "harder" and the elastic modulus G remains roughly constant whatever the strain amplitude γ < γ max , as illustrated in Fig. 2. Increasing the strain amplitude again and above γ max is associated with further increase of the elastic modulus, as long as the material is not irreversibly fractured. To our knowledge, such striking irreversible hardening has never been reported before: it suggests that in NRL gels, strain hardening is not associated with a transient increase in σ/γ, but instead with a permanent increase in the elastic modulus, together with an extension of the material linear regime up to γ max , the maximum strain amplitude reached in the material previous history. Hardening in NRL gels should therefore be associated with some irreversible modifications of the microstructure at the local scale. Note also that during the last cycles in Fig. 2, a slight dip is observed in G' when the material is deformed in a range of strains where it has been previously hardened. This reproducible behavior was observed identically during a further gradual increase or decrease in strain.

Second, hardening is observed in deformation ranges that are much larger than previously reported for colloidal gels. This allows us to investigate large increases in the stress. In particular, when switching to stress-imposed experiments, we observe in Fig. 3 that during hardening, the gels elastic modulus remains of the same order of magnitude as the imposed stress amplitude. Thus the material adapts to the stress in order to sustain it as long as fracture does not occur. This provides an elegant way of tuning the elastic modulus of NRL gels at will, simply by applying an oscillatory stress with an amplitude of the order of the target modulus. Of course, this behavior has some limitations: the maximum stress and therefore the elastic modulus that a gel can reach depend on its volume fraction, as seen in Fig. 3. The maximum stress at rupture is reached for φ v ≈ 0.162, and decreases for φ v > 0.32. Below, we first characterize NRL gels behavior at rest, before investigating the hardening behavior in more details.

Fractal scaling of the elastic modulus at rest

We investigated NRL gels for volume fractions ranging from 0.011 to 0.53. As the minimum volume fraction for gel formation is φ v ≈ 0.01, and the volume fraction at random closed packing is φ rcp ≈ 0.715 [6], this roughly covers the full range of volume fractions accessible for gel formation. The elastic modulus at rest spans more than five decades, from 7 × 10 -2 Pa to 2.3 × 10 4 Pa, and scales with the particle volume fraction φ v as G 0 ∼ φ 3.3 v (Figure 4). We checked that measurements conducted in geometries with different materials and different roughnesses lead to the exact same results, confirming the absence of wall slip. The power-law behavior of G 0 with φ v is reminiscent of fractal gels obtained by colloidal aggregation in the diffusion-limited cluster aggregation regime [22].

In fractal colloidal gels, the elastic modulus is given by

G 0 = (k 0 /a) φ (1+2 +d b )/(3-d f )
v [16,[23][24][25]. In this equation, k 0 is the spring constant between a pair of interacting particles, a is the particle radius, ε expresses the anisotropy of the backbone (from ε = 0 for a straight chain to ε = 1 for a purely isotropic chain), and d b and d f are the fractal dimensions of the backbone and the cluster structure respectively. If the in-FIG. 4. Elastic modulus at rest G 0 (filled symbols) and maximum elastic modulus before rupture G max (empty symbols) as a function of the rubber particles volume fraction φ v . The symbols and colors correspond to different Couette geometries: smooth stainless-steel (blues circles), serrated stainless steel (red diamonds) and smooth PMMA (green triangles). The red vertical dashed line shows the volume fraction at random close packing, estimated at 0.715 based on sample polydispersity (see Ref. [6]). The dashed lines represent the best power-law fits of the data, and the vertical gray line shows the limit of hardening. Inset: change in the hardening amplitude, G max /G 0 , with φ v .

teraction potential between particles is isotropic, aggregation leads to isotropic structures, hence ε = 1. In this case, the network response to shear stress is controlled by bond-bending interactions between the primary particles that constitute the gel, and we obtain (3 [24], it follows that 1.59 ≤ d f ≤ 1.79. This is characteristic of a diffusion-limited cluster aggregation regime, for which d f ≈ 1.7 -1.8 [26][27][28], while the reaction-limited cluster aggregation regime is associated with d f ≈ 2.0 -2.1 [28,29].

+ d b )/(3 -d f ) = 3.3. As 1 ≤ d b ≤ 5/3
From the values of the elastic modulus, it is also possible to estimate the spring constant between particles. Previous works on colloidal gels generally assumes that bending mechanics is dominated by contact interactions localized at the interface between two particles, which strongly depend on the physico-chemical environment [30]. For rubber particles, the surface is covered by proteins and polar lipids, such as phospholipids and glycolipids [31,32]. They play a key role in latex stability and in the specific mechanical properties of dense rubber. We might therefore expect mechanical properties that strongly depend on pH or temperature, which strongly affect protein-protein or protein-phospholipid interactions.

However, we surprisingly observed previously [6] that NRL gel mechanical properties are not affected by variations in the physico-chemical environment, for a pH ranging from 3.5 to 5.5, and temperatures ranging from 5 to 40 • C. This unique behavior might originate from the nature of the rubber particles. In NRL, the core of the particles, composed of poly(cis) isoprene, is in a rubbery state at room temperature, well above the polymer glass transition that occurs around -73 • C [33]. This is in sharp contrast with most of the colloidal gels in-vestigated so far that are made from "hard" colloidal particles whose core is in a glassy or crystalline state. Consequently, in NRL gels, polymer chains from two different particles might locally interpenetrate at their contact point, so that the path from one center of mass to the other consists of a rubber continuum. In such a case, the mechanical properties associated with the overall inter-particle bending rigidity are not controlled by surface interactions, as usually assumed in colloidal gels, but instead by the rubber elastic modulus itself. If this holds true, the spring constant k 0 between two particle centers of mass should be equal to the dense rubber elastic modulus, G rubber , times the distance between two centers of mass, here equal to two particle radii. Taking G rubber ≈0.372+/-0.006 MPa [34] for the rubber elastic modulus of polyisoprene chains, a value which can vary depending on rubber properties, we find k 0 ≈ 0.37 Pa m.

From our experimental results, we determine the prefactor of the power law fitted in Fig. 4 as k 0 /a = 269 kPa. With a = 0.5 µm, this provides an estimate for k 0 of 0.54 Pa m, in reasonable agreement with the above estimate from rubber elastic modulus. This strongly suggests that the surface properties of the primary particles, covered in this case by proteins and phospholipids, do not play a major role in the linear mechanical properties of NRL gels, unlike what is usually observed in colloidal gels.

It is noteworthy that the scaling of the elastic modulus with φ v holds for the full range of accessible volume fractions, while the fractal description should fail at high φ v . Indeed, the average cluster size at the gel point, R c , follows

R C = aφ -1/(3-d f )
v [16]. For volume fractions above 0.14 (for the case d b = 1.59) or 0.28 (for the case d b = 1.79), the cluster radius drops below four particle radii, for which the fractal description is irrelevant. Still, as shown by others studies focused on internal dynamics of colloidal gels, the scaling for G 0 remains valid, and there is no sharp mechanical signature of the transition from the dilute to the concentrated (φ v > 0.1) regime [35].

Hardening amplitude is volume fraction dependent

Let us now turn to the dependence of hardening on the volume fraction of rubber particles. As shown in the inset of Fig. 4, the relative amplitude of hardening, measured as the ratio between the elastic modulus before rupture G max and the modulus at rest G 0 , decreases sharply with φ v . For the lowest accessible volume fraction φ v = 0.01, where the strain hardening amplitude is maximum, G max is about 200 times greater than G 0 (see also Fig. 3). By contrast, for φ v ≥ 0.25, no strain hardening is observed. This limit is in the volume fraction range previously determined for which the fractal description does not hold anymore, suggesting that the hardening phenomenon originates from the fractal nature of the particulate stress-bearing network.

The absolute value of G max also depends on the particle volume fraction, with G max φ 1.64 v . The value of 1.64 is not compatible with a behavior dominated by bending interactions, for which ε = 1 and the exponent of φ v is given by (3

+ d b )/(3 -d f )
. Considering the uncertainty surrounding the measurement of the exponent (about +/-0.05), a precise quantification of ε, d b and d f in hardened materials is difficult. Nevertheless, this value implies that during hardening, ε decreases, i.e. that structural reorganizations occurring during hardening lead to increasing backbone anisotropy. If we assume that before rupture, the mechanical response is dominated by bond-stretching, i.e. ε ≈ 0 and d b ≈ 1 [24,36], we obtain d f ≈ 1.7 -1.8, which constitutes the upper limit for d f , suggesting that the cluster fractal dimension is not significantly modified by hardening.

Characterization of the hardening effect in the nonlinear regime

To characterize the nonlinear features of the hardening phenomenon, we now focus on the Fourier components of the strain response as a function of the imposed stress. As shown previously in Fig. 1, the strain response is distorted during hardening. In Fig. 5a, the amplitude γ k of the harmonics of the strain response is plotted as a function of the amplitude σ of the sinusoidal stress input.

As already described above in Fig. 1, the low-stress regime is characterized by a linear scaling of the strain fundamental amplitude γ 1 with the stress amplitude, γ 1 = σ/G 0 , where G 0 is the elastic modulus of the gel at rest (Fig. 5a). Above the onset of hardening, i.e. for σ > σ 0 and γ > γ 0 , one has γ ∼ σ α with α ≈ 0.3 consistently with the scaling observed in Fig. 1. In the nonlinear regime, hardening is also characterized by the presence of odd harmonics, the amplitude of which γ k increases with σ. At high stress amplitudes, the harmonic content converges toward γ k /γ 1 ∼ k -β with β close to 2 (Fig. 5b).

Not surprisingly, the particle volume fraction affects the properties of the strain response (Fig. 6a). In particular, the onset stress for hardening σ 0 increases with φ v (Fig. 6b), as already observed in Fig. 3. Figure 6b shows that σ 0 ∼ φ 3.3 v i.e. the onset stress follows the same scaling as G 0 . In other words, σ 0 is linearly related to G 0 . In contrast, the strain amplitude at the onset of hardening γ 0 remains constant, γ 0 ≈ 0.2, for all but the smallest volume fraction (φ v = 0.011) for which the gels withstood a higher deformation before hardening (Fig. 6c). We conclude that above φ v = 0.02, the parameter that controls the hardening is the strain amplitude.

Deep into the hardening regime, the harmonic content behaves as γ k /γ 1 ∼ k -β for all the volume fractions (inset of Fig. 6a). We observe a slight increase in the exponent β from 1.9 to 2.1 with φ v in Fig. 6d. Unfortunately, non-linear measurements at φ v > 0.05 were out of reach in the PMMA Couette cell since the torque exceeded the limit of the rheometer. Thus, it is not clear whether the exponent β, and therefore the shape of the strain response, changes significantly with the particle volume fraction. Still, φ v strongly affects the hardening exponent α which increases with φ v as α ∼ φ 0.42 v . This

γ k /γ 1 ∼ k -2 .
means that the gel hardening declines as φ v increases. This is also seen in Fig. 3 where the slope of G as a function of σ decreases with φ v . These observations are in stark contrast with the reversible stress stiffening that characterizes collagen tissues [21], for which the elastic modulus scales linearly with the stress, implying that α ≈ 0. In the present NRL gels, a very low value of α, which also implies that the strain amplitude at rupture is very close to γ 0 , is only observed for the weakest NRL gel at φ v = 0.01. For all the other gels, hardening is observed along a significant range of deformations, typically from γ 0 ≈ 0.2 to γ max ≈ 0.6 -2.

Discussion

In the following, we use the general framework of fractal gels to interpret our data and draw a simplified picture of the structural rearrangements that could lie behind strain hardening. This picture is consistent with our rheological observations but remains hypothetical in the absence of direct structural characterization of NRL gels during hardening. It should therefore be considered as a first step prior to complete system characterization. As already mentioned above, colloidal gels can be described as a dense packing of fractal clusters produced in situ by the fractal aggregation of colloidal particles (Fig. 7a). They have been extensively studied in the litera-ture [37,38]. The present results show that the spectacular irreversible hardening behavior of NRL gels depends on two control parameters: the particle volume fraction φ v and the strain amplitude γ.

We propose to describe the changes in the gel structure as a function of strain in four steps, as illustrated in Fig. 7b: 1. At rest, for γ = 0, the gel is described as a continuous network of rubber particles. The characteristic size is the cluster size, above which the system appears homogeneous. Clusters are characterized by a fractal dimension, d f , which is independent of φ v , and by their radius, R c , which decreases with the volume fraction [39]. A fraction of the particles constituting the cluster is involved in the stress-bearing backbone, which ensures the stability of the self-supported gel. This backbone is characterized by a fractal dimension, d b , with d b < d f . Another fraction of the particles is included in dangling strands that do not support stresses. At very high φ v , clusters are made of a very small number of particles. In such dense materials, it is likely that there are no dangling strands. However as φ v decreases, R c increases, and the internal density of the clusters decreases. Consequently, the proportion of particles involved in dangling strands might be maximum in the gels prepared at the lowest volume fraction, and almost negligible in concentrated gels, for which the fractal theory does not hold anymore.

2. Under small strain, for γ < 0.2, the response to mechan- ical shear stress or strain is that of a fractal gel. The elastic modulus scales with φ 3.3 v , a relationship that suggests a cluster fractal dimension below 1.8, characteristic of gels formed in the diffusion-limited cluster aggregation regime. The mechanical response is governed by the resistance to bending of the interparticle bonds that are included in the cluster backbone, i.e. the particles that are involved in the stress-bearing network. In this regime, strain and stress are linearly related.

3. In the hardening regime, for γ > 0.2, the material opposes an increasing resistance to deformation. Hardening is however not observed for high volume fractions, which suggests that it is associated with rearrangements of the fractal structure rather than with local modifications of the particles properties themselves. Our hypothesis is that under shear, local displacements gather dangling chains on the backbone. As contacts between rubber particles lead to irreversible aggregation [6], the backbone reinforcement is also irreversible. As the proportion of dangling strands should decrease with φ v , so does the maximum relative amplitude of the hardening phenomenon (inset of Fig. 4). This scenario is compatible with the results obtained from numerical simulations of colloidal gels by Colombo et al. [40], which predict that before yielding, strain can induce the formation of new bonds responsible for material stiffening. For even larger deformations, dangling strands could also interact with a neighboring backbone, creating a new connection between two adjacent backbones. These new bridges could be under tension when the overall material is at rest, but relaxed when the material is returned to the strain under which they were formed. This would explain the slight but very reproducible dip in G observed during the last cycles in Fig. 2 for strains corresponding to previous hardening, a behavior that was also reported in recent simulations of soft gels [41].

4. Lastly, when all the dangling strands are included in the backbone network, stress-bearing backbones behave as stretched chains, and further deformation leads to gel fracture.

To summarize, the above mechanism is quite unique and may be linked to two specificities of rubber particles. First, the irreversible nature of their interaction is such that when a dangling strand is captured by a neighboring backbone, it remains "stuck" to the backbone. Second, the core of NRL particles is in a rubbery state so that particles can locally in-FIG. 7. Structural changes in NRL gels during hardening as suggested from their rheological characterization. (a) NRL gels are formed by space-filling clusters, whose size R c decreases with the particle volume fraction φ v , but whose fractal dimension d f is assumed to remain constant with φ v . (b) In the linear regime, material deformation is locally associated with bending between particles. In the hardening regime, dangling strands are gradually incorporated into the backbone. At very large deformations, the backbone is stretched until it breaks. (c) The interparticle potential can be represented by a deep attractive well followed by a long tail: the harmonic regime is probed at small deformations while large deformations probe the anharmonic part of the potential.

terpenetrate and interparticle bonds can sustain a very large strain before breaking. We can therefore probe large deformations and access the anharmonic part of the interparticle potential, leading to the emergence of harmonics in the rheological response (Fig. 7c). These two specificities are directly linked to the response of the stress-bearing part of the network under stress, and not to the fractal flocks interactions. They are therefore not specific to the general framework that we choose to describe our system, and would remain valid even if the local structure were better described by the so-called bondbending percolation model [42] or by that of a heterogeneous glass [43].

CONCLUSIONS

In this study, we investigated the rheological properties of Natural Rubber Latex gels. We discovered a striking hardening behavior under oscillatory shear. Although strain-and stress-controlled experiments are equivalent, we found that strain controls hardening as it appears above a characteristic strain γ 0 ≈ 0.2 independent of the volume fraction of rubber particles. We showed that this strain-hardening displays two unique features: (i) it is fast and irreversible and (ii) it occurs on an extended range of strains, from γ 0 ≈ 0.2 up to γ ≈ 0.6 -2. The relative amplitude of hardening decreases with the particle volume fraction, and ultrasonic imaging un-der oscillatory stress showed it to be homogeneously distributed in the gel on scales larger than a few tens of micrometers. Compared to previous reports of strain hardening in colloidal gels [16,18], the irreversible nature of this behavior in NRL gels enables precise tuning of the material's mechanical properties. We proposed a simplified picture based on a fractal description of the gel to explain these observations. Our hypothesis is that hardening is associated with the gradual incorporation of dangling strands into the stress-bearing backbone. In NRL gels, the persistence of dangling strands at rest (before deformation) could be associated with the ability of particles to partially fuse after a prolonged contact. This would explain the irreversibility of their aggregation, but would also result in oriented, non-mobile attractive interactions. Indeed, the Peclet number associated with bending a single bond can be written as kT/ G rubber a 3 ≈ kT/ k 0 a 2 ≈ 10 -7 in the present case. It results that a dangling strand, once formed, should be extremely rigid and only weakly subjected to Brownian forces.

Strain hardening/stiffening plays a key role in the mechanical behavior and properties of living cells [44], but it has been rarely reported in colloidal gels composed of isotropic particles. As NRL gels are composed of large particles that can be easily observed by confocal microscopy, and because their structural reorganization under strain is irreversible, a precise microstructural characterization of these materials before and after hardening is now within reach. Thus, NRL gels appear as a model tunable system for thoroughly investigating the local features of strain-stiffening in future work. Numerical simulations should also help to confirm the microscopic scenario for hardening proposed in the present work. This work was carried out with support from CNPq, National Council for Scientific and Technological Development -Brazil.

FIG. 1 .FIG. 2 .

 12 FIG. 1. (a) Evidence for hardening in NRL gels: change in the elastic modulus, G , under imposed oscillatory strain (open symbols) and under imposed oscillatory stress (filled symbols) for a NRL gel at φ v = 0.04. The continuous line indicates the value G 0 = 5 Pa of the elastic modulus at rest and the dashed line is the best power-law fit G ∼ γ 1/α-1 with α=0.32 in the hardening regime, before rupture. Their intersection defines the onset of hardening (γ 0 ,G 0 ). The relative amplitude of the harmonics γ k /γ 1 obtained from the Fourier decomposition of the strain response is shown for two stress-imposed measurements in the linear and in the hardening regimes. Inset: stress-strain relationship for the same data set. (b) local rheology from echography (LORE) in Couette geometry for a stress amplitude σ = 1.8 Pa (corresponding to a strain amplitude γ = 0.23). Left panel: Local displacements ∆(r, t) (linear color scale from brown at the stator to yellow at the rotor) and strain γ(t) recorded by the rheometer (black line) in response to oscillatory stress σ(t) (blue dashed line) as a function of the normalized time ft. All quantities are normalized. Right panel: the amplitude ∆ 1 of the fundamental Fourier mode of ∆(r, t) as a function of r/e where r is the distance from the stator and e is the gap width. The dashed line shows the theoretical profile for a linear homogeneous strain field in Couette geometry.

FIG. 5 .

 5 FIG. 5. Strain-stress relationship for φ v = 0.04. (a) Change in strain harmonics γ k as a function of the applied stress amplitude σ. The dotted line is the best linear fit γ 1 = σ/G 0 in the linear regime with G 0 = 5 Pa. The dashed line is the best power-law fit in the nonlinear regime γ 1 ∼ σ α with α = 0.32. (b) Change in the relative intensity of odd harmonics γ k as a function of their order k. The dashed line is γ k /γ 1 ∼ k -2 .

FIG. 6 .

 6 FIG. 6. (a) Normalized strain-stress relationships for different particle volume fractions. The black dotted line is σ/σ 0 = γ/γ 0 where σ 0 and γ 0 are the stress and strain amplitudes at the onset of hardening. Colored lines are power-law fits γ ∼ σ α in the hardening regime. Inset: Change in the relative intensity of the odd harmonics γ k as a function of their order k. Color codes for the volume fraction. The Fourier series are computed deep into the nonlinear regime for σ ≈ 100σ 0 . The dashed line is γ k /γ 1 ∼ k -β with β = 2. (b) σ 0 as function of φ v . (c) γ 0 as function of φ v . (d) β as function of φ v . (e) α as function of φ v .