{'fraction_non_alphanumeric': 0.031618137363858635, 'fraction_numerical': 0.01294732162702823, 'mean_word_length': 4.31497341996456, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stiffening attractive suspensions with solid inclusions or granularity effect in loaded colloidal gels

The elastic properties of a soft matter material can be greatly altered by the presence of solid inclusions whose microscopic properties, such as their size and interactions, can have a dramatic effect. In order to shed light on these effects we use extensive rheology computer simulations to investigate colloidal gels with solid inclusions of different size. We show that the elastic properties vary in a highly non trivial way as consequence of the interactions between the gel backbone and the inclusions. In particular, we show that the key aspects are the presence of the gel backbone and its mechanical alteration originated by the inclusions. To confirm our observations and their generality, we performed experiments on an emulsion that presents strong analogies with colloidal gels and confirms the trends observed in the simulations.


I. INTRODUCTION
Composite materials have been created since antiquity to combine and enhance the mechanical properties of the constituents. Since Einstein thought about the effects of a low concentration of solid particles on the viscosity of Newtonian fluids, experimental and theoretical works have explored the mechanical properties of a wide range of such materials.
The types of matrix and particle that can be used are enormous and the variability in the properties of the resulting materials equally. For example, the addition of nanoparticles or fibrous objects can be highly efficient in changing the elastic modulus, because of the high surface area of nanoparticles and the capacity of the fibres to create percolated networks even at low volume fractions [1]. In the absence of specific particle-matrix interactions, classical mechanics can be used to estimate the elastic modulus of such composite materials. Those are based on the assumption that the matrix is a continuous linear elastic medium embedding ellipsoidal inclusions [2,3]. The Krieger and Dougherty (KD) equation [4] was found to describe well the elastic modulus of various elastic matrices loaded with monodisperse or bimodal non-Brownian spherical particles [5,6], although the hypothesis of continuous material around the inclusions was not fully justified for some of the studied matrices such as emulsion or bentonite. This could be understood by the fact that the inclusions were large compared to the characteristic length-scale of those matrices. Recent work with liquid foam loaded with solid particles showed that the elastic modulus was correctly predicted by the KD equation for inclusion-to-bubble size ratios close to unity [7]. On the other hand, it has been shown that if the particles are small enough to concentrate between the bubbles and to form a skeleton-like granular packing, stronger strengthening can be obtained [8,9]. When the particles were attractive, and formed a gel within the foam channels, bulk gel elasticity was not sufficient to describe the strengthening as the gel. This is because the gel time (time at which an elasticity can be measured) changes with confinement [10].
Those results revealed that finite size effects are not easy to anticipate for a particle loaded matrix, especially if the matrix can form various different micro-structures. This is the case for colloidal gels for which the individual particles aggregate and cluster to form stress-bearing networks of gels. At low volume fractions (φ g < 0.1) it is understood that the particles form clusters which then percolate [11]. At high volume fractions (φ g > 0.5) the arrested states are spatially homogeneous, the particles are localised between their neighbours and can be described as attractive glasses [12]. Those structures can be described using the cluster size as the characteristic length scale at low φ g , whereas the particle size is the relevant size at high φ g . The transition from one to the other was predicted by Zaccone et al [13], where they also showed how the cluster-cluster interactions control gel elasticity. The role of the number density of the clusters connections was further explored by Whitaker et al [14].
The aim of this work is to discuss how colloidal gels' mechanical properties are modified by the presence of non interacting spherical inclusions of different sizes. In the first part of the manuscript, we will present a simple numerical model in which the colloidal gel is the result of an arrested phase separation [15,16]. Using a finite strain rate algorithm [17,18], we measured the mechanical properties as a function of volume fraction of the inclusions as well as of their size. The picture that emerges suggests that the size ratio between the particle of the gel and the inclusions plays a major role. In fact, when the inclusion's size is comparable to the size of the gel particles, the resulting composite gel is noticeably re- inforced and can sustain much larger stresses. We will show that this effect results from a non trivial coupling between the structure and local mechanical stress that results in an imbalanced distribution of local stresses. In a second part of the paper, we will present an experimental model of a colloidal gel, made of emulsion drops, reinforced with solid polystyrene particles and we show a semi-quantitative agreement with the simulation predictions. The similarity between the results from the experimental and numerical models, despite their intrinsic differences, suggests that the reinforcing effect that we observe in our composite gel could be general and translated in other similar systems.

II. METHODS
In this paragraph we will introduce the numerical approach, as well as the experimental system.

Numerical model
In our simulations, we use a colloidal mixture model made of N spherical particles with two distinguishable species: gel particles and inclusion particles. The particles tagged as g represent the building block of the gel backbone, which is formed by particles with two different diameters, at equal relative concentrations, that have an average diameter of D g . Particles tagged as inc model the spherical solid inclusions of diameter D inc . The interactions are described by the potential where r ij = |r j − r i | denotes the distance between two particles (i and j), D ij is the average diameter of the interacting particles i and j, s(i) is the species of particle i and s (j) is the species of particle j. The constant U ss is chosen accordingly such that u ss (r ss ) = 0. The exponents p and q in Eq. 1 define the range of the attraction, and are set to p = 36 and q = 24. The constant ensures that the minimum of the potential is unitary while the cut-off parameter r ss sets the range of the interaction. For simplicity, we have set = 1.
The cut-off r ss is chosen such that where D min = (3/2) 1/2 D ij is the distance at which the potential in Eq. 1 has its minimum. Then, with Eq. 2, the potential in Eq. 1 describes gel-gel interactions as attractive while all other interactions as completely repulsive. We have chosen the two different diameters of the gel particles as D g1 = 0.909 D 0 and D g2 = 1.091 D 0 . With this, the average diameter of the N g gel particles is D g = (1/2)(D g1 + D g2 ) = D 0 . The unit of lengths are fixed as we have set D 0 = 1. To characterize the size of the inclusions, we use the ratio σ between the diameter of the inclusion particles and the average diameter of the gel particles, σ = D inc / D g . The ratio values σ studied here are 0.8, 1.0, 1.3, 2.0, 3.0, 4.0, 5.0 and 10.0. The total packing fraction of our simulation box is set as φ T = φ g + φ inc , with φ g and φ inc the effective packing fraction of the gel and inclusion particles in the simulation box. To study the effect of different loads of inclusion particles in the final composite structures, we chose the total packing fraction range 0.30 ≤ φ T ≤ 0.60, with 0.10 ≤ φ g ≤ 0.50 and 0.10 ≤ φ inc ≤ 0.40. The total number of particles in our systems is 1.4 × 10 4 ≤ N ≤ 6.6 × 10 4 .
All simulations were performed with the open-source MD simulation package LAMMPS [19]. For each total packing fraction value studied, we start with an initial isotropic configuration, with no overlaps, which we generate randomly at the desired total packing fraction. After this, the interaction described by Eq. 1 is turned on and the system is allowed to evolve at k B T / = 0.01 for a period of at least 4 × 10 3 τ * , where τ * = m(D 0 ) 2 / is the molecular dynamics time unit in our simulations. The time step in the initial sample preparation runs is δt = 5 × 10 −4 τ * . After the initial runs, we apply the following dissipative viscous dynamics to our system [17] m Where m is the particle mass and ξ the friction coefficient. We use m/ξ = 1.0τ * , which corresponds to the overdamped limit [20]. This dynamics is run for at least 1×10 4 τ * , which draws the kinetic energy to a value lower than 10 −10 and leads to an arrested configuration. From this stage onward, all time steps of the simulations are set to δt = 5 × 10 −3 τ * .
In Fig 1, we present a typical final arrested configuration for a gel network (red particles) with solid inclusions (blue particles) at a total packing fraction φ T = 0.60 and equal relative concentrations, φ g = 0.30 and φ inc = 0.30, with a ratio σ = 3 for the inclusions. Here it can be seen that the solid inclusions are scattered in the simulation box, without any apparent preference, and that the gel particles are forming a network. In the following sections, we will expand more on the composite structures explored.

Mechanical tests
The mechanical properties of our composite systems were tested using continuous strain deformation and small amplitude oscillatory rheology (SAOR), as previously reported for colloidal gels [18]. To study samples under continuous strain deformation, we applied a series of incremental strain steps, in simple shear geometry, to each sample analysed [17,18]. This is done by first applying, for each step, the instantaneous affine deformation Γ ∆γ to the simulation box, which represents a simple shear in the xy plane. After a strain step deformation of ∆γ is imposed, the boundary conditions are updated and the simulation box is then relaxed by applying the damped dynamics described in Eq.3. This procedure is repeated for n steps, obtaining a cumulative strain in the system of γ = n∆γ. We have chosen a strain increment of ∆γ = 0.001 and a relaxation time interval ∆t such that the shear rate iṡ γ = ∆γ/∆t = 10 −5 τ −1 0 , with τ 0 = ξD 2 0 / the time it takes a particle subjected to a force /D 0 move a distance of D 0 .
The average state of stress of the system is given by the tensor Σ. In our simulations we compute each component Σ αβ using the standard viral definition [17,18], where α, β stand for the Cartesian components {x, y, z}.
In this analysis, we consider only the pairwise energy contribution, all other terms are assumed negligible as the velocities for all particles are kept small for all shear rates. For the step-strain curves, we compute the shear component Σ xy which corresponds to the same component of the deformation. In the following sections we refer to Σ xy as Σ. For our SAOR analysis, the frequency dependent response is measured by imposing the oscillatory shear strain γ(t) = γ 0 sin(ωt) to the system, with the dynamics given by During each cycle, we compute the shear stress component Σ, as for the step-strain analysis. The first harmonic expansion of the shear stress, for a visco-elastic solid, can be written as: With G (ω, γ 0 ) and G (ω, γ 0 ) the elastic and plastic moduli, respectively. In the present work we will focus our analysis on G , which is refered to as G in the following sections. The modulus is extracted by a simple Fourier transform of Σ and γ.

B. Experiments
We have made weakly elastic emulsions of rapeseed oil (from Brassica Rapa, Sigma-Aldrich) in aqueous solutions of sodium dodecyl sulphate (SDS, Sigma-Aldrich) at 30 g/L. The surfactant solution is freshly prepared by dissolving SDS in deionised water and used within one day. The emulsions are made using the double syringe method [21]: two syringes (60 mL, Codan Medical) are partially filled with oil and SDS solution respectively and connected with a double luer lock. The oil and the aqueous phase are then manually pushed through the connection 20 times, at which point a homogeneous emulsion is formed. A stock emulsion at 0.70 of oil is prepared with an average drop diameter, 2 R g of 5 µm. Emulsions at lower oil volume fractions, between 0.63 and 0.44, are obtained by dilution with the same SDS solution. Solid spheres (polystyrene or poly(methyl methacrylate), (Microbeads) of different sizes (from 6 to 250 µm) are dispersed by carefully mixing them into the diluted emulsions using a spatula. The elastic moduli of the composite materials are measured with a MCR302 rheometer (Anton Paar) by performing oscillatory strain sweeps (frequency 1 Hz, amplitude from 10 −5 to 1) in a Couette geometry with sandblasted surface. All experiments are performed at (20.0 ± 0.5) • C.
The high concentration of surfactant used, ensures the presence of micelles and screening of electrostatic repulsion. Although the exact form of the attraction is not clear, the emulsions exhibit an elastic modulus at packing fractions below random close packing indicating the presence of an attractive interaction [22].
For the experiments with the confocal microscope the oil was dyed with pyrromethene. The inclusion particles were not dyed so appear as black spots. The concentration of emulsion in the confocal microscopy experiment was φ * g = 0.5.

A. Numerical results
We have measured the elastic modulus G at varying gel volume fractions φ g as shown in Figure 2. Our analysis focuses on the limit of small frequency, and all our G values are obtained using an equivalent shear rate ωτ 0 = 1 × 10 −3 and γ 0 = 1 × 10 −4 . The filled black data points in Figure 2 correspond to samples without inclusions, pure gels, and as expected the elastic modulus of the gel increases with φ g . The concentrations studied in this work are high enough so that the percolating structures are not particle fractals, and models based on this assumption are indeed not suitable here. In this volume fraction range, the elasticity of the gel is no longer controlled by the individual particles but by the properties of the clusters [13,23,24]. In the range φ g = 0.2 − 0.5, the data are well described using a 2.8 power law. This power law behaviour is close to that found in experiments by Romer et al [24], where they measured 3.1 at intermediate volume fractions up to around φ g = 0.25.
Above φ g = 0.5, the elastic modulus is higher than predicted by the power-law, probably because we are getting deep into the glassy regime. On the other side, we observe that for our lowest density, i.e. φ g = 0.1, G g is smaller than expected. This could be an indication that, for this specific concentration the fractal nature of the aggregates plays a significant role.
Next, we focus our attention on the mechanical properties of the composite gels that, as discussed before, are made by adding inclusions to the gel (different inclusion sizes σ and volume fractions φ inc ). The corresponding elastic modulus G for a given φ g , as obtained by the computer simulations, is presented as empty symbols in Figure 2. These data points correspond to samples with size ratio values 1.0 ≤ σ ≤ 10.0 and total packing fractions of 0.30 ≤ φ T ≤ 0.6. We first present the general trends observed, later in the paper we will analyse those trends in more detail. In most cases the effect of the inclusions is to increase G. However at low gel volume fractions G is shown to decrease.
To highlight the mechanical effect of the inclusions, we now plot G/G g (φ g ) in Fig 3a as a function of φ inc , the volume fraction of inclusions. The presence of the inclusions leads to a slight stiffening of the material in most cases, but two data series show a stronger increase. These are both at total volume fraction, i.e. gel particles + inclusions, φ T = 0.6 and with the small inclusions: σ = 1 or 2. In such cases the elastic modulus can increase at most by a factor 6.
Comparison with theory leads us to calculate also the ratio G/G * g , where the elastic modulus of the matrix G * g is described by G * g = G g (φ * g ) using φ * g = φ g /(1 − φ inc ), i.e. the volume fraction of the gel embedding the inclusions. This normalisation is shown in Fig 3b. Again, we observe two different behaviours depending on the total volume fraction φ T . For φ T = 0.6, the effect of the inclusions is to stiffen the composite gel with a magnitude that depends strongly on the size ratio σ: the smaller is σ the more pronounced is the variation, reaching a maximum of roughly 2 times the elastic modulus of the pure gel.
As a theoretical reference we will consider the Krieger and Dougherty (KD) equation [4], which was originally proposed to describe the shear viscosity of non-Brownian solid particles in Newtonian fluids [25]. Note that the problem of elasticity of rigid particles embedded with a linear elastic material is formally similar to the problem of the viscosity of rigid particles suspended in a Newtonian liquid. In the KD equation the elastic modulus increases with the volume fraction of solid inclusions φ inc , as: with φ m the close packing volume fraction at which the elastic modulus diverges and where the coefficient x = 2.5 has been conveniently introduced to recover the Einstein equation [26] at low concentration (first order in φ inc ), i.e. G/G * g 1 + 5 2 φ inc . For spherical inclusions characterized by slip condition at their surface, just like the inclusions we have simulated, the equation derived by Taylor [27,28], i.e. G/G * g 1 + φ inc , is more appropriate at low volume fraction. For higher volume fractions, empirical equations were proposed [29] but, for the sake of simplicity, it is sufficient to use Eq. 7 with x = 1 for comparison with our simulations. Such a comparison is presented in Figure 3b, using φ m = 0.59. The elastic modulus is below the KD prediction as a function of φ inc for all the samples except for the σ = 1 at φ T = 0.6.
To discuss further the results, we present G/G * g in Fig. 3c as a function of σ for a fixed φ inc = 0.1 but varying the gel concentration φ g . When the total volume fraction is below 0.6, the elastic modulus displays a weak or moderate dependence on σ, essentially a decrease. On the other side, when the total volume fraction reaches 0.6, the system does display a remarkable sensitivity to the size of the inclusions, with an up to one and a half times larger value of the elastic modulus with respect to the pure gel. We can conclude that the smaller are the inclusions the stronger is the effect on the mechanical properties, a behaviour whose origin can be understood by looking at the microscopic interplay between the gel and the inclusions.
In such a granular (colloidal) system, the elastic modulus depends on the number density of contacts, the stiffness of the contact and the gel structure (texture). By considering a dense homogeneous gel at total volume fraction φ T = 0.6, we can assume that the structure is not significantly modified if gel particles are replaced by inclusions of the same size. In this case, we could write G g (φ T ) ∝ zφ T χ g−g / R g , where z is the coordination number and χ g−g is the contact stiffness between two gel particles [13]. In the composite gel at φ T , volume fraction φ inc of gel particles have been replaced by nonattractive inclusions, and this leads to a decrease of the modulus. With those assumptions, the corresponding modulus writes G ∝ zφT Rg φg φT χ g−g + φinc φT χ g−i , where χ g−i represents the effective contact stiffness between gel particles and inclusion particles. Therefore, the reduced elastic modulus can be expressed as: A single parameter fit to the data, shown by the dashed line in Fig. 4, yields χg−i χg−g = 0.34. This shows how the replacement of g − g interaction by g − i interaction contributes to decrease the shear elastic modulus of those loaded particle gels. As can be seen in Fig. 4, the effect of size ratio σ is to decrease the elastic modulus with respect to Eq. 8. This behavior is expected to be partly related to the well-known increase of the maximum volume fraction of bimodal granular assemblies [30]. Note that here we consider that φ T is constant, so the volume fraction of gel particles between the inclusions φ * g = (φ T − φ inc )/(1 − φ inc ) decreases significantly as a function of φ inc .
More drastic changes can be anticipated in the gel structure due to the introduction of inclusions. As an example, and with contrast to continuous matrices, gel particles can be excluded from the gaps between inclusions. To highlight such structure effect, we performed mechanical tests on samples where we removed the inclusions. In this protocol, the gel assembles in the presence of inclusions, but the latter are then removed before proceeding to the mechanical test. In Fig. 5a, we compare the system with and without the inclusions in a similar fashion to what was done in Fig. 3, but here only for φ T = 0.6 where the size effects appeared to be most prominent. It is shown that the dependence on σ is much weaker than the one observed for systems including the inclusions. At this stage it is interesting to compare these results with the theory for continuous matrix with holes, as given by micromechanics calculation [28,31]: G r /G * g = 1 − 5φ/(3 + 2φ), where φ is the volume fraction of holes (here we have φ = φ inc ). Comparison is shown in Fig 5b for two size ratios (σ = 10 and σ = 1). It appears that results for σ = 10 are compatible with the theory, which can be understood by the fact that φ * g is high enough for the gel particle assembly to be shaped by the precursor inclusions as a thick skeleton (see Fig. 9). In contrast, for σ = 1 such a skeleton does not exist because the precursor inclusions have altered the gel particle network (see Fig. 9) as numerous weakness points. As a result, the second configuration is weaker than the first one although gel volume fractions are the same.
G r appears to be the relevant elastic modulus to be considered when quantifying the effects of the inclusions, so we divided the composite gel modulus G by its modulus with removed inclusions, G r . Moreover, one expects finite size of the gel particles to come into play when the mean gap size e between inclusions is close or smaller than the gel particle size, i.e. 2 R g . According to the parameters defined above, the ratio e/2 R g The resulting values are used to plot G/G r against e/2 R g in Fig. 5c, where a fair collapse is obtained for all the different cases. This indicates that the two contributions, i.e. interaction and structure effects, can be decoupled for this total volume fraction and, once the change in the backbone structure has been taken into account, effects of interactions of the inclusions with the gel particles can be described by the ratio e/2 R g . It appears that although the effective stiffness is weaker for interactions involving inclusions, small inclusions (i.e. e/2 R g < 1) are embedded within the gel particle backbone and provide strengthening by supporting the backbone. On the other hand, larger inclusions (i.e. e/2 R g 1) are not included into the backbone but rather shape it around them, providing strengthening comparable to prediction of the classical theory, i.e. KD.
To pin down the microscopic origin of the mechanical response discussed above, we performed continuous strain deformation simulations to unveil the microscopic stress experienced by each particle in our composite gels. We focus on the case of φ inc = 0.30 and φ g = 0.30 and we simulated, using step-strain deformation, a continuous deformation in the xy plane. The resulting stressstrain curves for different inclusion sizes as well as for the pure gel case are displayed in Fig. 6. We immediately notice a strong difference between the response of the pure gel and the composite gel, however this does not come as a surprise given the fact that for the pure gel there are no inclusions that reinforce the overall structure. A more relevant comparison is for systems with different size of inclusions and, in this case, several different trends are observed. In agreement with the measurement of the elastic modulus discussed above, for small strain the slope of the curve grows when the size of the inclusions decreases. For large inclusions, we notice that, upon straining, stress builds up continuously. For small inclusions, in particular when their size is similar to the size of the gel, the strain-stress curve presents clear signs of yielding. In this case, at γ 0.10 we notice that the curve has leveled off and reached a value that is more than three times larger than the stress observed for larger inclusion sizes, which confirms the differences that we have observed before.
Next, we move to the discussion of the structural properties of the composite gel and, from now on, we will focus only on two representative cases of size ratios, namely σ = 1 and σ = 10. In Fig. 7 (a) and  (b), the inclusion-gel radial distribution functions are presented as a function of the non-dimensional distance (r − D min )/2 R g , where D min is the cut-off distance after which there is no net force between two particles and 2 R g is the average diameter of the gel particles.
A remarkable difference appears between the two cases, despite the fact that the σ = 10 curves present some noise due to finite size. For the large size ratio, the distribution is almost flat while, for smaller size ratio, there is a pronounced peak. This is an important information because only the particles for which (r−D min )/2 R g < 0 will exert a net force on each other.
A more direct way to test this observation is to perform a direct calculation of the local stress acting on the gel backbone. To this aim, we calculated the local stress tensor Σ i for each particle i and its three eigenvalues Σ i 1 ≥ Σ i 2 ≥ Σ i 3 which represent the three invariant principal stresses in the frame of reference of the stress. This stress tensor has been normalized with the volume of the corresponding particle. We focus on Σ i 3 , which indicates the largest compression on each particle, and from now on we will refer to it as Σ 3 . Note that this value is not the same as the shear component Σ used in previous sections. In the calculation of the stress tensor for each particle, all the interactions between particles are taken into account, this includes gel-gel, gel-inclusion and inclusion-inclusion interactions. For different strain values, we have tracked the variation of Σ 3 with respect to the undeformed case γ = 0, i.e. Σ 3 − Σ 3 (0). The histograms for this quantity, for σ = 1 and σ = 10, are presented in Fig. 8a and Fig. 8b respectively for different values of γ. In Fig. 8c we calculate the skewness of the two distributions as a function of the strain. In the case σ = 1, the distribution remains basically symmetric around the minima whereas, for σ = 10, the distribution becomes remarkably skewed toward the negative values.
This trend is clearly related to the structural properties of the composite gel and, in Fig. 9, we present snapshots of the same systems for γ = 0.3. Inclusion particles are coloured in blue while the gel particles are coloured according to their proximity to an inclusion. In particular, if their distance from an inclusion is less than D min (the cut-off of the potential) the particles are coloured in red. It is evident that, for σ = 10, the number of particles in close proximity to an inclusion is much smaller than for the case σ = 1. In the second case almost all particles in the gel are in contact with at least one inclusion. The inclusions reinforce the system and help reduce compression in comparison with the pure system (or the composite mixture with σ = 10). We conclude that the unexpected trend of the mechanical response of our composite gel is the result of a subtle interplay between the structural properties and the mechanical interactions between the inclusions and the gel and for this reason it depends strongly on the size ratio between the two.

B. Experimental results
As explained in the Sec II B, we explored the influence of the matrix grain size using oil-in-water emulsions. Those emulsions are made of polydisperse oil droplets with an average size of around 5 µm. Several common points with our numerical system can be noted as follows: (1) the droplets can be well-approximated by spheres for the oil fractions we studied, and (2) the drops are stabilized using an excess of ionic surfactant, which is at the origin of the attractive interaction between the drops (de- pletion forces). With respect to the second point, note that electrostatics just plays a marginal role in the interaction since it is screened due to the high concentration of surfactant counter-ions.
On the other hand, our experimental system could ap- pear to be very different from the ideal numerical system due to the intrinsic softness of the oil droplets. Actually, we think this difference is an interesting point to investigate. Indeed, the addition of solid spheres (inclusions) in oil-in-water emulsions has been shown to increase the elastic modulus according to the KD model [5,6]. But this was done only in the limit where the size of the inclusions is much bigger than the characteristic size in the matrix (σ 1). Therefore this system is of particular interest to seek for the finite size effects highlighted in the previous section. Besides, note that our results have been obtained out of the size ratio range corresponding to σ 1, for which the inclusion particles are allowed to organize between the particles of the matrix in the shape of a skeleton-like packing, providing strongly amplified strengthening [7][8][9].
A microscope image of an emulsion without inclusions with φ g = 0.5 is shown in Figure 10a. The emulsion drops have been dyed using pyrromethene, which fluoresces in red. In this attractive system, as we are below close packing of the emulsions, the system will present a gel-like structure with its typical heterogenous structure. The size of the pores will depend on the attraction as well as on the oil volume fraction. Once added to the system, the inclusions will sit in these cavities and, as a consequence, will profoundly affect the structure of the gel matrix. An example of the resulting composite material is shown in Figure 10b with φ inc = 0.1 and with σ = 2.5. We measure the elastic modulus of the composite gel at different φ inc for inclusion sizes corresponding to σ = 4 and σ = 20 as calculated with the average emulsion drop size. As it could be expected, the elastic modulus of the composite gel increases with increased φ inc . In the experiments we expect that the boundary condition is a no-slip condition [5,6], and we compare to the KD model with x = 5 2 . The two measurements of G/G * g are shown in Fig. 11a as a function of φ inc . The red line shows the KD prediction with φ m = 0.59 that describes the results well.
Having seen that for σ = 4 the KD model works well, we explore the effect of varying σ. In Fig. 11b, the experimental composite gel elastic moduli are normalised using the KD prediction as a function of inclusion size. For σ larger than 2 or 3, the enhancement of the inclusions to the elastic modulus compares well with the KD model. However, at smaller σ we observe an upturn. As the inclusion size approaches the average drop size we see a strong increase in the elastic modulus. Such a strong strengthening behaviour revealed for σ 1 was not observed by Gorlier et al. [7] for particle loaded foams. This suggests that the nature of interactions between the particles, i.e. here attractive vs repulsive for liquid foams, is a crucial parameter for strengthening. Therefore, and in a more general way, it appears that tuning σ in colloidal gels can be of great interest to optimize their strength.
The σ effect highlighted in Fig. 11b is reminiscent of the simulation results shown in Fig. 5c. In order to compare more closely with our numerical results, we plot in Fig. 11c G/G * g and G/G r as a function of the gap size between inclusion particles, i.e. e/2 R g = σ(( φm φinc ) 1 3 − 1). As can be seen, the agreement with numerical results is remarkable, which gives significant insight into the finite size effects highlighted with the attractive emulsion: small inclusions (i.e. e/2 R g < 1) are allowed to be embedded within the gel particle backbone and to provide optimal strengthening by supporting the backbone, while larger inclusions (i.e. e/2 R g 1) are not included into the backbone and behave like classical solid inclusions, providing strengthening in agreement with continuum mechanics. For practical purposes, the measured finite size effect can be reasonably described by the following approximate function: Although this empirical equation can be useful to evaluate the magnitude of finite size effects for e/2 R g < 2 − 3, we stress that (1) it concerns attractive gels only, and (2) it does not capture properly strengthening in the limit e/2 R g 1, which is better described by the KD equation.

IV. CONCLUSION
We have combined rheology computer simulations with experiments to explore the elastic modulus of composite attractive gels. Our work reveals significant finite size effects in the mechanical response of such systems and gives physical insight into their origin. These effects depend on the total volume fraction (φ T ) in the system (gel + inclusions): At low total volume fraction (φ T < 0.6) the addition of inclusions leads to a weak stiffening of the material, below the predictions of continuum mechanics (KD model). However at φ T = 0.6 the effect of the inclusions becomes strongly size dependent.
With large inclusions σ 1 the elastic modulus of the attractive composite gels is correctly described using the KD prediction. On the contrary, as the inclusion size decreases we observe considerable stiffening of the material. These findings can be explained by the effects of the inclusions, they can change the structure of the gel backbone and the average interactions, as the gel-inclusion interactions are weaker than the gel-gel interactions. If the inclusions are large the gel backbone structure is only weakly modified by the presence of the inclusions, however once they become small they participate actively in the gel structure and hence the mechanical properties.
The degree of stiffening is well accounted for by the distance between the inclusions. Once the distance between the inclusions becomes comparable to the size of the gel particles finite size effects dominate and continuum mechanics approaches do not capture the mechanical properties of the material. This result is confirmed by both the simulations and the experiments. This despite several differences between systems studied in the simulation and experiments, in particular in terms of gel particle polydispersity, interactions and effective modulus of the particles. This suggests that the reported findings concern a broad range of composite materials. It should be stressed however that the attractive nature of interactions in the studied gels seems to be necessary for those finite size effects to be significant.
The results provide guidance for the control over the mechanical properties of colloidal gels using inclusions. However they can also help understand the many practical materials where solid inclusions are found in emulsions or gels, and the size ratios of the particles is no longer large. Such composite materials are present in drilling fluids, food or cosmetic products.