Sylvain B Énito

François Molino

Charles-Henri Bruneau

Thierry Colin

Cyprien Gay

Non-linear oscillatory rheological properties of a generic continuum foam model: comparison with experiments and shear-banding predictions

Keywords: numbers: 47, 57, Bc Foams and emulsions 83, 10, Gr Constitutive relations -83, 80, Iz Emulsions and foams in Rheology -83, 50, Ax Steady shear flows, viscometric flow I, SHEAR BANDS AND FOAM RHEOLOGY

The occurence of shear bands in a complex fluid is generally understood as resulting from a structural evolution of the material under shear, which leads (from a theoretical perspective) to a non-monotonic stationnary flow curve related to the coexistence of different states of the material under shear. In this paper we present a scenario for shear-banding in a particular class of complex fluids, namely foams and concentrated emulsions, which differs from other scenarii in two important ways. First, the appearance of shear bands is shown to be possible both without any intrinsic physical evolution of the material (e.g. via a parameter coupled to the flow such as concentration or entanglements) and without any finite critical shear rate below which the flow does not remain stationary and homogeneous. Secondly, the appearance of shear bands depends on the initial conditions, i.e., the preparation of the material. In other words, it is history dependent. This behaviour relies on the tensorial character of the underlying model (2D or 3D) and is triggered by an initially inhomogeneous strain distribution in the material. The shear rate displays a discontinuity at the band boundary, whose amplitude is history dependent and thus depends on the sample preparation.

It may seem paradoxical that a single material, when submitted to a uniform shear stress σ xy , between two parallel plates or two coaxial cylinders, may be observed simultaneously in two distinct states in different regions of the flow. This 'shear bands' observation has nevertheless become common since the early 1990s in a variety of complex fluids: they appear and are stable [1][2][3][4], or sometimes fluctuate [5][6][7][8]. These bands are most of the time parallel with the shearing plates [1], with a different shear rate in each band.

The current understanding of these observations relies in general on two essential ingredients : (i) a structural evolution of the material under shear, and (ii) a stress response that decreases as a function of the shear rate (within a particular range). This decrease is the mechanical signature of the structural evolution of the fluid and is the source of the mechanical instability that triggers the appearance of bands [9].

In polymer melts or entangled polymer solutions [10] and in entangled giant micelle solutions [9], the flow elongates the objects, which alters the apparent viscosity of the material (which must be evaluated after subtracting the effect of wall slip [11]). The fact that this viscosity goes down is principally due to the average orientation of the objects in the shear flow.

In lyotropic lamellar phases, the transition can be associated with the reorganisation of the films into onion-like multilamellar vesicle systems [5,6,12,13], also exhibiting wall slip behavior [5]. In micellar cubic crystals the transition consists in an ordering of the initial polycrystal into a single crystal with specific planes becoming aligned with the plates [14,15]. In the last two cases, no microscopic interpretation of the decrease in effective viscosity occuring during the transition is currently available.

In granular materials, surface flow is a particular case of shear bands: the lower band is in this case blocked (zero shear) while the flowing region is sheared. Again, no complete structural description is available. Nevertheless, it is admitted that through dilatancy, which reflects the necessity for the grains to move a little bit apart in order to move past each other [16], the shear generates a difference in volume fraction between the flowing region and the blocked one. This lower volume fraction tends to facilitate the flow in the flowing region even more as compared to the blocked region, thus stabilizing shearbanding. When it is present, gravity is of course essential: it favours this phenomenon by helping the system segregate into a dense, blocked region (located at the bottom if the particles are denser than the fluid) and a less dense, flowing region. Thus, it allows to determine the concentration profile [17].

In foams and emulsions, the situation is less clear. Shear bands were observed in 2D [18,19]. In some cases, the observed shear-banding could result trivially from shear stress inhomogeneity, due to cylindical Couette geometry (σ(r) ∝ 1/r 2 ) or enhanced (for 2D foams) by the presence of at least one solid boundary [20][21][22][23] whose friction on the foam implies that r 2 σ(r) is not uniform. In some more interesting cases, the shear rate is spatially discontinuous at the boundary between the blocked and the sheared regions [19,24]. This discontinuity may arise from an apparently intrinsic impossibility for the foam to be deformed homogeneously at low shear rates [25], which then implies the presence of shear bands at low shear rates. But at least in some cases, the shear rate at the boundary is not unique for a given system [24] and is thus not intrinsic. Recently, the very existence of such a finite shear rate at the boundary has been seriously questioned [26]. As we shall see, the present work highlights yet another (history dependent) possible origin of the shear rate spatial discontinuity, arising from the tensorial character of the material response. For these materials also, no complete structural description accounts for flow localization in a satisfactory manner. Dilatancy, which corresponds to a local change in water concentration φ, certainly plays an important role by easing the relative motion of bubbles or drops, although it behaves somewhat differently from granular materials depending on the liquid fraction [27][28][29]. The structural disorder is also invoked as a parameter coupled to the flow [30]. In both cases, the local fluidity (ratio of the shear rate and the shear stress) is enhanced.

B. Foam rheology

The specificity of foams as compared to other materials is the following: not only do they flow substantially only above some threshold stress, but they also undergo large elastic deformations at lower stress. Hence, classical models such as visco-elastic fluids (well suited for polymeric fluids) and elasto-plastic solids (well suited for metals) are unsufficient to capture the behaviour of foams and emulsions. In the past few years, much effort has been devoted to address this challenge and describe the richer mechanical behaviour of foams. Several rheological models have thus emerged [31][32][33][34][35]. They all assume that the foam is essentially incompressible. Within this perimeter, some models are purely visco-elastic yet with a non-linear elasticity [31]. As for the models incorporating plasticity, they can be dispatched into two categories: the creep term either depends on the stress and deformation rate [33,35,36] or on the stress only [32,34]. Finally, these models also differ in the tensorial form of elasticity and creep, a feature which is relevant for non strictly 2D systems (or for compressible materials).

Despite this variety of models, most current experiments are not sufficiently stringent to fully test these models and decide which ingredients are indeed relevant to describe the mechanical response of foams. For instance, classical linear rheology experiments, particularly oscillatory measurements, are clearly unable to provide much information about the behaviour under large stress.

Because the constitutive objects of foams are macroscopic and can be observed directly, statistical tools have been elaborated to measure the local deformation and deformation rate. Using these tools, more complex geometries such as flows around obstacles are also used in order to subject the foam to a tensorially broader variety of sollicitations [37]. Yet because these experiments are conducted in a confined geometry, the unknown friction with the walls and the fact that no direct measurement of the total stress is conducted, make it difficult to test stress predictions beyond low velocities.

To complement this, in order to test the full time response of the models even within classical geometries such as those available in a rheometer, a broad range of experiments could be elaborated by choosing many different forms for the time dependence of the applied deformation or stress. As a first step towards this goal, considering the full time-dependent response of a foam subjected to large amplitude oscillatory shear (not only the usually extracted storage and loss moduli), is susceptible to provide more stringent tests for models. Such experiments have been conducted recently: because the flow was observed to remain homogeneous, the measured behaviour can be robustly attributed to the local, yet macroscopic, mechanical response of a 3D foam [38]. The strain-stress (Lissajous) curves display various shapes and show that such data can become available and provide non trivial results. These results will be discussed later in the present work.

C. Scope of the present work

In materials whose plastic threshold corresponds to a small deformation, shear banding requires a fluidizing mechanism such as those mentioned in paragraphe I A above. But for materials whose deformation at plasticity onset is large, like foams, the tensorial nature of the material state, due to stored deformation, is sufficient to obtain shear bands: no extra mechanism is required.

In particular, the model that we suggest does not incorporate such an ingredient as dilatancy. We know that a local plastic event results in an elastic redistribution of stress in the neighbourhood [39,40]. In simple shear geometry, it thus favours flow localization [41][42][43]. In a statistical manner, it then raises locally the material fluidity [44] and generates a non-local material rheology. This non-local character had been observed in concentrated emulsions flowing in microfluidic channels [45]. Let us emphasize that these non-local effects are intrinsically present in our modelling since the underlying elastic propagators [41][42][43][44] result directly from the elastic continuum medium equations that we use.

The paper is organized as follows. We first discuss several types of continuum models for foams (Section II). We then recall (Section III) the construction of a rather generic continuum model for foam or emulsion rheology [34]. It is generic as for the elasticity, the plastic flow rate and the specifically three-dimensional form of the response. We then simulate large amplitude oscillatory shear (Section IV) and conduct a first round of comparison with published data on such experiments at a fixed frequency for various amplitudes [38]. We are not able to reproduce the experimental data in a reasonable way with a single set of parameters, but in the future, fitting similar data over a full range of both frequency and amplitude will be a very effective and stringent method for testing more general models than the present one. Finally, we thoroughly discuss shear banding in our model (Section V). Note that in this work, we restrain ourselves to a strictly mechanical and thermodynamical formulation. The important problems related to the coupling between the rheological behaviour and the structure of the material are not discussed. This coupling is experimentally well documented in various complex fluids systems in which shear bands are de facto associated with structural transitions [9].

Here, we ask a more restricted question: could stationary shear bands in foams and emulsions be accounted for, starting from an inhomogeneous initial stress distribution in the material with otherwise strictly homogeneous mechanical properties? Our main result: shear bands can emerge in a structurally homogeneous material under shear, only due to an inhomogeneous distribution of the initial internal stress in the material. We demonstrate this for a physically very natural form of the elastic and plastic laws.

II. CHOOSING THE TYPE OF CONTINUUM MODEL

As mentioned above, our main point is to explore models without any extra dynamic variables apart from the stored local deformation. In the present Section, we will discuss these models in very general terms, omitting any explicit tensorial features. Here, for simplicity, we assume that the spring and both viscous elements are linear, and that element 3 cannot withstand any stress beyond σ3.

(a)

(b) (d) (c) (e)

A. Models with only one internal variable

The most general arrangement of rheological elements having only one internal deformation variable is represented in Fig. 1a. The deformation of the (not necessarily linear) spring represents the deformation of the local structure. Three creep elements (viscous and/or plastic) can be included, as shown.

At this point, let us recall that foams are viscoelastic under weak stress conditions. Hence, both the spring itself and the combination of spring and creep elements must be able to deform under weak applied stress. As a result, creep elements 1 and 2 must flow under weak stress: we cannot choose them with a stress threshold below which they would not flow at all. In other words, they are purely viscous (although not necessarily linear). By contrast, creep element 3 must have a stress threshold so that the entire system also displays a stress threshold. These considerations are summarized schematically in Fig. 1b.

Although elements 1 and 2 are viscous, they play different roles when some creep motion of element 3 is involved. Let us first illustrate this point by considering an experiment in which we impose a constant deformation rate from t = 0 and reverse the deformation rate as of t = T . For simplicity, we consider a linear spring with modulus G 1 , Newtonian viscosities η 1 and η 2 . Furthermore, we consider that element 3 is simply a solid friction element with threshold σ 3 with no dependence on velocity. Fig. 2 shows the contributions of viscous elements 1 and 2 separately in the response of such a system to a triangle wave deformation. In both cases, the stress jumps up immediately to a finite value at t = 0 due to the viscous elements η 1 and η 2 . The stress then increases at a constant rate as the spring elongates. When the threshold of element 3 is reached, the stress saturates and remains constant. At t = T , when the deformation rate is reversed, the stress jumps down by a finite amount. It then decreases at a constant rate as the spring is relaxed and later stretched in the reverse direction.

There are two differences between the effects of viscous elements 1 and 2. The first difference is that the observed threshold depends on the deformation rate in the case of viscous element 2. But that feature is not essential: one can always decide that the deformation rate of creep element 3 affects its stress (in other words, by considering that it contains not only a solid friction element, but also an extra viscous element in parallel with it).

The second difference between both situations of Fig. 2 is more essential. The jumps in stress at t = 0 and at t = T have equal magnitudes in the case of viscous element 1. By contrast, in the case of element 2, the magnitude of the second jump is twice as large as that of the initial jump (except if T is too short for the system to be able to relax the Voigt element, with timescale η 1 /G 1 , after it has reached the threshold stress). More generally, in such an experiment, one can express each viscosity (at the applied deformation rate) in terms of the magnitudes Spring G ′ 1 elongates immediately (red segment) while spring G1 responds with some delay (green curve) due to viscous element η1. Viscous element η3 gives rise to a constant additional deformation rate (blue segment). (b) Combined Bingham-Burger model: solid friction element σ3 (complemented by viscous element η ′

3 ) is now included so as to provide additional creep above the stress threshold. We believe that some tensorial version of this model could mimic the rheological behaviour of a foam quite adequately. (c) Model studied in the present work. Apart from the additional viscosity η2 introduced in Fig. 1, it represents the combined Bingham-Burger model when the (blue) viscous element has not been deformed yet, either at intermediate time scales when the (green) Voigt element has relaxed (the green and red springs then respond in series) or at short time scales when the Voigt element is still blocked (in which case only the red spring responds).

of the stress jumps at t = 0, when the deformation rate changes from zero to + γ, and at t = T , when it is reversed from + γ to -γ:

η 1 ( γ) = 2|∆σ(0)| -|∆σ(T )| γ ( 1
)
η 2 ( γ) = |∆σ(T )| -|∆σ(0)| γ (2)
When the value of the elastic deformation e can be measured independently, for instance through optical measurements in 2D foams and relevant statistical tools [37], the respective contributions from elements 1 and 2 can be obtained by comparing σ and e (see full lines versus dashed lines in Fig. 2).

B. Weak stress: role of one extra internal variable

Below the flow threshold, the model outlined above behaves like a Voigt element (a spring in parallel with a viscous element). The behaviour of a foam under a weak stress is in fact a little more complex: a linear Burger model (see Fig. 3a) is known to correctly reproduce stepwise creep experiments on liquid foams [46]. Because the Burger model contains two springs, it corresponds to a system with two internal variables (spring elonga-tions e and e ′ ) rather than one. Fig. 3a shows the response of such a model to a stepwise creep experiment. The stress jump generates an immediate elongation of the (red) spring labeled G ′ 1 . The (green) Voigt element then relaxes within a time scale of order η 1 /G 1 . On very long time scales, the (blue) viscous element gives rise to a slow drift with velocity σ/η 3 (depicted by the finite slope of the blue line).

In order to build a model that behaves like the Burger model under weak stresses but which presents a flow threshold like discussed earlier, one could combine the Burger model with the model presented in Fig. 1 for only one internal variable. We would thus obtain the model represented in Fig. 3b (which we already suggested as a generalization of our model [34]).

In the present work, we consider this combined model of Fig. 3b, but we focus on short and intermediate time scales where the highly viscous (blue) element has not moved yet. Because the blue element has not moved, it can be simply omitted. As for the two springs G 1 and G ′ 1 and the viscous element η 1 , their behaviour can be reduced to that of a single spring in two limits. On time scales much shorter than η 1 /G 1 , the (green) Voigt element is blocked due to its viscous part. As a result, the whole system behaves like the red spring G ′ 1 . At intermediate time scales (much longer than η 1 /G 1 but still with a blocked blue viscous element), the Voigt element has relaxed. Hence, both springs are simply in series: they combine into a composite spring. In Fig. 3c, we have represented such a model. The green and red spring represents either the red spring G ′ 1 (on short time scales) or the combination

G 1 G ′ 1 /(G 1 + G ′ 1
) (on intermediate time scales). Meanwhile, the other creep elements provide both the stress threshold (solid friction σ 3 ) and the dependence on deformation rate (viscous element η ′ 3 ). We also include a general viscosity η 2 as in the discussion of Paragraph II A. In this model, the spring and viscous elements must be understood as non-linear, unless stated otherwise. Apart from the viscous element η 2 , the model of Fig. 3c is identical to the model that we constructed earlier and for which we had analysed the local, mean field behaviour [34].

III. CONSTRUCTING THE TENSORIAL MODEL

In the present Section, we will briefly recall how we built the rheological model [34] of Fig. 3c. In particular, it is based on a general nonlinear description of elasticity and plasticity. Indeed, materials such as foams can locally undergo large elastic deformations -located far from the linear regime corresponding to small deformations -before plastic flow occurs [START_REF] Weaire | The Physics of Foams[END_REF][START_REF] Cantat | Les mousses -structure et dynamique[END_REF].

A. General local rheological laws

The relevant framework to describe elastic stresses in a flowing material is the Eulerian one, whether this material possesses elastical properties or not. Indeed, during the flow of a foam or an emulsion, even though elastic stresses exist, any reminiscence of a reference state disappears continuously due to plasticity. The Lagrangian description, which is based on maintaining the correspondance with such an intial state of reference, is formally equivalent, but less adapted conceptually and numerically.

Thus we attach the variables describing the material to a spatial grid (x, y, z), and they correspond to an instantaneous and local description in space.

In this framwork, only two variables are relevant in a strictly mechanical context: the local velocity gradient ∇ v(x, y, z) and the local deformation state stored in the material [34] (green-red spring in Fig. 3c), as described in continuum mechanics by the Finger tensor B(x, y, z) [START_REF] Bertram | Elasticity and Plasticity of large deformations[END_REF]. Note that when the material is at rest, B = I while the stored deformation, depicted schematically in Fig. 3c, vanishes: e = 0.

In this section, we describe our local rheological model. Note that in this local context, the global tensor ∇ v itself has to be considered as an independent local threedimensional tensorial variable, just as B, not as the spatial gradient of a velocity field. Only when we will turn to the description of a spatial system, see section V E, will the vector field v(x, y, z) be introduced. Meanwhile, tensors ∇ v and B will thus be the two variables of our local tensorial model.

The elastic part of the stress, which goes throuth the spring in Fig. 3c depends on the deformation according to the following relation, the most general one compatible with the symmetry constraints in three dimensions [34]:

σ el = a 0 I + a 1 B + a 2 B 2 , (3)
where a 0 , a 1 and a 2 are scalar functions of the invariants of the Finger tensor B.

Turning to plasticity (σ 3 and η ′ 3 in Fig. 3c), we only assume that every event of plastic relaxation is aligned with the stored deformation. The plastic creep D B p should thus be similarly aligned. The most general form compatible with the symmetry constraints is then:

D B p = b0 I + b1 B + b2 B 2 , (4)
where b0 , b1 and b2 are again scalar functions of the invariants of the Finger tensor B.

To complete the model, we gather together in a global viscosity term (which was noted η 2 in Fig. 3c) all the dissipative phenomena which are present even in the absence of any plastic event in the foam. They occur for example at small scales: flows in films squeezed between bubbles or in Plateau borders. We simplify its description in selecting a Newtonian average viscosity η s for these local dissipative phenomena. The list of contributions to the stresses in the material is thus closed. We have:

σ = a 0 I + a 1 B + a 2 B 2 + η s 2 (∇ v + ∇ v T ). (5)
To take into account the incompressible character of foams and emulsions, we add an extra kinematic constraint of strict volume conservation det(B) = 1. Refering to [34] for further details, we take it into account by using only the deviatoric part of the stress:

σ = dev(σ) = σ - I d tr(σ). (6)
The same constraint on plasticity gives the general form [34]:

D B p = B •dev(f (B)) = b 1 B •dev(B)+b 2 B •dev(B 2 ), ( 7
)
where the scalar prefactors b 1 and b 2 are isotropic, and thus depend on the invariants of tensor B.

In what follows, we will use a completely equivalent form of tensor D B p which manifests more clearly that the dissipation is positive (see the discussion in [34]):

D B p = A(B) τ B • G(B) (8)
where A(B) is a scalar isotropic function of B, τ the characteristic time of the dissipative processes; moreover:

G(B) = dev [P(B) • dev(σ el )] tr [P(B) • dev(σ el ) • dev(σ el )] , (9)
with P is a function of the form

P(B) = b(B) B -2 + (1 - b(B)) B 2 [34]
, where b is an isotropic function. In this expression, the total dissipation per unit volume is A(B) and can be chosen as positive.

Eventually one gets the generic local rheological model:

dB dt -∇ v • B -B • ∇ v T = -2 D B p , (10)
D B p = A(B) τ B • G(B), (11)
σ = a 0 I + a 1 B + a 2 B 2 + η s 2 (∇ v + ∇ v T ), (12)
where dB/dt = ∂B/∂t+( v•∇)B is the particulate derivative of the Finger tensor.

B. Complete spatial model

As for any local rheological model, the previous equations must be complemented by field equations which express force balance and mass conservation:

∇ • σ + ρ f = ρ d v dt -∇ p, ( 13
)
∂ρ ∂t + ∇ • (ρ v) = dρ dt + ρ tr 1 2 (∇ v + ∇ v T ) = 0, ( 14
)
where f represents the external forces (per unit mass), and ρ is density. The incompressibility constraint gives here

∇ • v = tr 1 2 (∇ v + ∇ v T ) = 0. ( 15
)
As a result, the density ρ is simply transported by the flow: dρ/dt = 0. In the remaining of this work, we furthermore assume that the density is homogeneous, hence it also remains constant: ∂ρ/∂t = 0.

Last assumption: we restrict ourselves to the Stokes regime, where inertial terms are all negligible in the mass conservation equation. Thus one obtains:

∇ • σ = -∇ p. ( 16
)
The complete system of equations that we need to integrate numerically is thus:

dB dt -∇ v • B -B • ∇ v T = -2 D B p , (17)
D B p = A(B) τ B • G(B), (18)
tr 1 2 (∇ v + ∇ v T ) = 0, (19)
σ = dev(σ) = dev a 1 B + a 2 B 2 + η s 2 (∇ v + ∇ v T ), (20)
∇ • σ = -∇ p, (21)
The initial conditions that must be specified to solve the above system may merely consist in the values of tensor B over the entire sample. Indeed, the value of the velocity and pressure fields can be derived therefrom using equations (21) and (20) which, when combined, are similar to Stokes' equation, using the constraint of equation (19).

C. Selection of a particular form of elasticity and plasticity

Elasticity: Mooney-Rivlin model

We have selected a usual form of incompressible elasticity which has been demonstrated to describe to a good approximation the nonlinear elastic behaviour of foams [START_REF] Höhler | [END_REF]51]: Mooney-Rivlin elasticity. The corresponding elastic energy per unit volume can be written [START_REF] Bertram | Elasticity and Plasticity of large deformations[END_REF]:

ρ E(B) = k 1 2 (I B -3) + k 2 2 (II B -3) (22)
where

I B = tr(B) (23)
II B = 1 2 [tr 2 (B) -tr(B 2 )] = tr(B -1 ) (24)
Going back the coeficients of Eq. ( 3), this corresponds to the following expressions:

a 1 = k 1 + k 2 I B (25) a 2 = -k 2 (26)
Following previous work refs. [START_REF] Höhler | [END_REF]51], we express the values of k 1 and k 2 using an elastic modulus G and an interpolation parameter a as follows:

k 1 = aG (27) k 2 = (1 -a)G. (28)
In the foam modelling litterature, a value a = 1/7 is sometimes recommended [START_REF] Höhler | [END_REF]51]. Keeping in mind our perspective of discussing the conditions for the appearance of shear bands depending on parameter values, in sections IV and beyond, we prefer to keep the parameter a free, although we remain in the framework of the Moonley-Rivlin elasticity.

Plasticity: yield stress fluid

The particular form of plasticity explored in this work is based on a nonlinear threshold-like behaviour. Locally, the plastic reorganisation events only occur in the material when the stored elastic deformation reaches a critical value. We express this transition with a function W y (B) which vanishes linearly at the threshold:

W y (B) = 0, (29)
with, in our case, W y (B) = ρE(B) -K, where ρE is the stored elastic energy per unit volume, and K a constant.

In simple shear from a relaxed state, σ y is the threshold stress: function W y vanishes.

From the point of view of the plastic deformation rate tensor D B p , we have the following expression 8, taking for A(B):

A(B) = (ρE(B) -K) Θ(ρE(B) -K), (30)
where Θ(x) = 1 when x ≥ 0 and Θ(x) = 0 elsewhere. We also set the following form for the polynom:

P(B) = bB -2 + (1 -b)B 2 . ( 31
)
with b between 0 and 1. Our final set of equations is thus:

dB dt -∇ v • B -B • ∇ v T = -2 D B p , (32)
D B p = ρE(B) -K τ Θ(ρE(B) -K) B • G(B), (33)
σ = dev(σ) = dev {(a G + (1 -a) G tr(B)) B -(1 -a)G B 2 +η s (∇ v + ∇ v T )/2 , (34)
∇ • σ = -∇p, (35)
tr 1 2 (∇ v + ∇ v T ) = 0. ( 36
)

Physical parameters and rheological model

In order to be able to highlight physically relevant quantities, we use a non-dimensional form of the above system. The elastic modulus G is taken as the unit of stress, and the weak stress relaxation timescale η s /G as the unit of time, while B is already dimensionless:

σ = σ/G, ( 37
) T = η s t/G, (38)
B = B. ( 39
)
As a result, the various quantities are nondimensionalized as follows:

Ê( B) = ρE(B)/G, (40)
K = K/G, (41) p
= p/G, (42)
Â(B) = A(B)/G, (43)
P(B) = P(B), (44)
Ĝ(B) = G G(B), (45)
∇ v = (η s /G) ∇ v, ( 46
) DB p = (η s /G) D B p . ( 47
)
The system of equations now reads:

d B dT = ∇ v • B + B • ∇ v T -2 DB p , (48) tr
( ∇ v + ∇ v T ) = 0, (49)
∇ • σ = -∇p, (50)
σ = σel + ∇ v + ∇ v T 2 , (51) σel
= (a + (1 -a) tr(B)) dev B -(1 -a) dev B2 , (52)
DB p = Ψ Â( B) B • Ĝ( B), (53)
Â( B) = ( Ê( B) -K) Θ( Ê( B) -K), (54)
Ê( B) = a 2 (I B -3) + 1 -a 2 (II B -3), (55) Ĝ
( B) = dev [P(B) • σel ] tr [P(B) • σel • σel )] , (56)
P( B) = b B-2 + (1 -b) B2 (57)
Ψ = η s G τ . ( 58
)
The new non-dimensional parameter Ψ defined in the last equation above reflects the ratio of the plastic flow rate (proportional to 1/τ ) to the viscoelastic flow rate (proportional to G/η s ) when the other factors have the same order of magnitude. With the present choice for the magnitude of D B p (with  proportional to the distance from the threshold), that occurs when the stored deformation is typically twice the threshold deformation.

D. Simple shear flow

In the remaining of this work, we address specifically the question of shear banding. For this purpose, we consider only simple shear flows. The velocity is oriented along axis x and varies along along axis y. The only non-zero component of the velocity gradient ∇ v is then ∂v x /∂y. The entire system and flow are invariant along x and z. Besides, the force balance given by Eq. ( 35) then implies that σ xy and σ yy are homogeneous at all times. In the axes x, y and z, the non-dimensionalized velocity gradient can thus be written as:

∇ v =   0 Γ(y) 0 0 0 0 0 0 0   ( 59
)
where Γ(y) = (η s /G) γ(y). Let Γ = (η s /G) γ be the macroscopic value of the shear rate at the scale of the entire sample. We now have five non-dimensional parameters: the plastic-to-viscoelastic flow rate ratio Ψ, the threshold K, the Mooney-Rivlin parameter a, the parameter b (which defines the tensorial form of the plasticity G(B)), and the macroscopic shear rate Γ.

For the sake of consistency, let us note that our earlier work [52] discussed parameters

α = 2 Γ (60) We = Γ Ψ (61)
instead of Ψ and Γ.

E. Relation between model parameters and experimentally measurable quantities

In fig. 4, we summarize the physical parameters included in our model.

In experiments, easily accessible dimensional parameters are the viscosity η s and the shear modulus G through linear rheology, as well as the threshold stress σ y . More elaborate setups can yield the value of a. There are indications that a value a = 1/7 is relevant for liquid foams [START_REF] Höhler | [END_REF]51].

our non-dimensional parameters, two can thus be determined easily: a and K. The latter is related to the energy at the threshold K ≈ 1 2 (σ y /G) 2 . As just mentionned Γ(y) = (η s /G) γ(y) is the normalised shear rate. Concerning Ψ = η s /(G τ ) and b, no experiment to our knowledge is able to validate or invalidate the value of the plastic reorganisation time τ at deformations beyond the threshold, or the tensorial form of the plastic flow (here expressed in terms of parameter b). For the time being, we thus consider parameters Ψ and b as free parameters in any comparison of our model with actual data. 27) and (28). The quantities σy and 1/τ constitute a scalar representation of the creep defined by D B p , and parameter b is the equivalent of a for creep, see equation (31). Finally, ηs is a viscosity that is independent of creep.

IV. HOMOGENEOUS FLOW BEHAVIOUR IN LARGE AMPLITUDE OSCILLATORY EXPERIMENTS

A. Method

In the present section, we test the predictions of our model by comparing them to the most stringent available measurements in homogeneous flows, namely the large amplitude oscillatory experiments conducted recently by Rouyer et al. [38]. We have simulated oscillatory shear flow with the local model (no dependence on coordinate y, i.e., homogeneous flow). In other words, in Eq. (59), we choose

Γ(y, t) = Γ(t) = -ωΓ 0 cos(ω t), (62)
which corresponds to the oscillating shear deformation:

Γ(y, t) = Γ(t) = Γ 0 sin(ω t), (63)
B. Typical behaviours Figs. 5 and6 show, in the form of Lissajous curves, the full response of the present model in the (amplitude, frequency)-plane. The normalized stress and strain responses make it easy to discriminate between plastic, elastic, or viscous behaviours of the model. Let us rationalize them in terms of the scalar diagram of Fig. 3c discussed above in Section II B.

A pure elastic behaviour corresponds to an ellipse squeezed into a straight line spanning the diagonal of the diagram. This is obtained at low frequencies and amplitudes. Indeed, deformation rates are then small at all times, hence the viscous elements in diagram 3c play no role. Meanwhile, because deformations remain small, the threshold of the solid friction element is never reached. As a result, the spring alone provides the mechanical response of the system. For the same reason, an elasto-plastic behaviour is expected at low frequency yet large amplitude, since the stress threshold is then reached. A purely elasto-plastic behaviour, as predicted by a scalar model, would correspond to a sharp-cornered parallelogram with two horizontal sides corresponding to the yield stress. The results of our simulation at low frequency and large amplitude differ from this simple picture in the same way as experimental data by Rouyer et al. [38], namely with two main features: (i) the "plastic part" of the Lissajous curve exhibits a slightly negative slope, and (ii) the transition to plasticity is progressive rather than sharp (blunt corners). Feature (i) corresponds to the weakening of the viscous component when the deformation rate decreases along the sinusoidal applied deformation. As for the latter feature, it can result either from viscosity being combined with plasticity (as in the present model [34]) or from a progressive onset of plasticity [36].

As compared to the results by Rouyer et al. [38], our model additionally exhibits an overshoot at very low frequency and large amplitude. Because such a regime is very similar to slow, continuous shear, this response can be understood [35] as a tensorial effect combining the saturation arising from plasticity and the rotation contained in shear (this point is further discussed at the end of Section V B).

This transition between a purely elastic response at low amplitude and an elasto-plastic response at higher amplitude is best illustrated by the top part of Fig. 6. The curves are normalized for clarity, but the actual maximum slope in each curve is essentially identical and is given by the shear modulus G. By contrast, the value of the stress in the most horizontal regions of the curve re- flect both the solid friction element and both viscous elements in diagram 3c. The bottom part of Fig. 6 presents results slightly below and slightly above the plasticity threshold. It shows that weak plasticity causes the deformation cycle to slowly drift towards a limit cycle that differs from the elastic cycle. This very drift, when continued along a longer cycle, is in fact at the origin of the overshoot mentioned above.

If we now turn to higher frequencies, the deformation rate becomes large. As a result, viscous elements now play a role and may even become dominant. Correspondingly, the Lissajous curve tend to become an ellipse whose axes lie along those of the figure. That is particularly clear on diagram 3c at high frequency with a large amplitude, but the trend is very obvious at high frequency and low amplitude, and is also discernable at low frequency and large amplitude.

C. Comparison with experiments

The data obtained by Rouyer et al. [38] correspond to a fixed frequency and different applied strain amplitudes (from 0.055 to 1.2). We have integrated [59] the equations of the present model with the same amplitudes and plotted them together with data. The results are presented in figure 7 (for clarity, stress curves are presented normalized).

Note that in order to obtain a reasonable agreement of our model with the data, we had to artifically choose different values of ω and G for each strain amplitude, while κ and Ψ could be kept constant. This unsatisfying ad hoc parameter adjustment shows the limits of this model in describing the behaviour of the foams studied by Rouyer et al.

V. SHEAR BANDING STUDY

γA γloc 8: Typical form of a stationnary flow curve giving the dependence of the shear stress on the local shear rate. σy is the yield stress as measured under imposed stress, and γc the corresponding shear rate. A macroscopic shear rate γA smaller than γc will not necessarily lead to a homogeneous velocity profile P A , with the expected stress σ A : the flow can separate into a blocked region and a flowing region (profiles P A 2 or P A 3 ). The local shear rate is then faster ( γloc A2 > γA and γloc A3 > γA ), which corresponds to a higher stress (σ A 2 > σ A and σ A 3 > σ A ). Besides, for an average shear rate γB greater than γc, the flow is homogeneous again, which corresponds to the expected stress σ B (greater than σy).

A3 γc γloc A2 γB = γloc B P A 2 P A 3 P B σ B σ y σ A 3 σ A 2 σ A P A σ d FIG.

A. Stationary flow curve and inhomogeneous flow

Let us now turn back to shear banding. For such a peculiar flow to be observed, the same material submitted to the same shear stress σ xy must be simultaneously in two different deformation states. As discussed for many years for various complex fluids [3,53,54], a mathematical condition for this to be possible is the existence of an unstable zone in the local flow curve σ ∞ ( γ) of the material: it must be non-monotous.

In the case of foams, nevertheless, such an unstable portion in the flow curve itself does not exist: how can shear bands with different shear rates coexist?

Foams and emulsions are instances of yield stress fluids, so that there exists a minimal value σ y of the stress σ xy below which no stationary flow occurs. Now when we shear the material, imposing the shear rate, the material has to flow, even for very small γ. The intrinsic flow curve thus possesses an extrapolation in stress when γ → 0. Let's denote it by σ d . Note that σ y and σ d pertain to the local rheology curve σ ∞ ( γ), not to the effective, macroscopic stationary curve as can be measured for example in a rheometer. In this discussion the flow is homogeneous. But the relative values of σ d and σ y , pertaining to the local flow curve, will give us hints about possible conditions for shear banding.

Let us now consider the result of a measurement made on a sample of this material, sheared in a (parallel or very low curvature) Couette cell under imposed shear rate. If σ d ≥ σ y , all parts of the sample will flow, even at low shear rates, since the corresponding stress is necessarily everywhere greater than the yield stress. As mentionned before, since the flow curve has no intrinsic instability for higher γ values, no mechanism is available for shear banding.

The situation is different if σ d < σ y . If we put on the same graph the yield stress σ y and the intrinsic stationary flow curve (figure 8), it is immediately apparent that this configuration allows for the coexistence of zones undergoing shear at rates γ such that 0 < γ ≤ γc , and of blocked zones remaining in the elastic regime at γ = 0. The mechanism is essentially the same as in the classical case of instability in the flow curve (see figure 8). Of course, as soon as γ > γc all regions flow, since γ > γc implies that some regions flow faster than γc . The stress in these regions, as given by the flow curve, has to be above the yield stress σ y . And since the shear stress is the same in the entire material, all regions support a stress greater than σ y and no region can be blocked.

B. Large elastic deformations versus extra dynamic variables

But is the situation where σ d = σ y actually possible? In various complex fluids, the answer is known to be yes. The usual explanation of such a flow curve is to invoque an internal extra variable (of a structural nature in general) which is coupled to the flow. As an example, in a simplified vision, this extra parameter can take one of two values: flowing or non-flowing. Thus the stationnary curve extrapolating to σ d at low γ, and the yield stress value σ y correspond in reality to two different materials, hence σ d and σ y can differ.

But as mentioned in Section I B, foams differ from many other complex fluids in that the deformation that must be reached to trigger plastic flow is large. This feature turns them into an intrinsically tensorial material.

In a stationary situation where shear banding is present, stress conservation implies that the shear stress σ xy is constant along the direction of the velocity gradient, as well as the stress component σ yy . By contrast, the extra components of the stress σ xx (y) and σ zz (y) may vary in an arbitrary manner along the direction of the velocity gradient. Among these, σ xx (y) is present even in a purely 2D system. These extra components will qualitatively play the same role as an extra structural variable in changing the local nature of the material, when viewed as a 1D material (along direction y).

But that only explains how it is possible for shear bands to be present. The reason why the flow curve actually extrapolates below the yield stress at vanishing shear rates (σ d < σ y ) in some tensorial models, thus allowing shear banding, has been shown by Raufaste et al. [35]: as long as the material remains elastic, the local deformation tensor is transported by the shear flow along a path that is not locally aligned with itself: the principal axes of the particulate time-derivative of the deformation do not coincide with those of the deformation. Hence, once plasticity is triggered, it alters the deformation evolution until it progressively reaches the locus where it is aligned with its transport under shear. At least in simple examples of elasticity and plasticity, this migration from the elastic path to the asymptotic locus is the origin of the shear stress overshoot observed during transients [35]. When the plastic flow is triggered rather abruptly, the system is still elastic just before the maximum of this overshoot, and the asymptotic shear stress value can then lie below the last elastic shear stress value. In other words, σ d < σ y .

C. History dependent shear bands

Despite some similarities, the analogy with systems characterised by unstable flow curves has some limitations. In the case of yield stress fluids, there is no unstable range in γ, which would impose phase separation between two phases at different flow rates. Shear bands are possible but not necessary. Also, no lever rule-like criterion can exist to select the relative fraction of the different bands, as have been argued in some fluid systems [3,55].

Rather, the initial distribution of σ xx (y) and σ zz (y) in the material will be of primary importance in the appearance of shear bands even though the material per se remains homogeneous. In other words, it is the material history that will lead to a particular flow profile. We will see that the initial distribution of stress in the material will determine the band structure.

D. 0D flow curve and shear banding criteria

As long as the flow in the material is homogeneous, a local rheological model will be sufficient to describe it. We begin by showing the typical flow curve corresponding to our model (fig. 9). Note that this flow curve is obtained under applied shear rate conditions.

As can be observed, the conditions described in the introduction for the appearance of shear bands are fulfilled: stress σ d is smaller than the static yield stress σ y . In the shear rate range between 0 and γc , the system has the possibility to split the average shear rate γ in different proportions of blocked and flowing bands.

Thus, in the homogeneous case, for any values of the parameters Ψ, K, a, b and Γ, we can use the local rheological model to calculate the static and dynamic thresholds, σ y and σ d , and the critical shear rate γc . Following the line of reasoning developped in the introduction, we can then predict the range of imposed shear rates [0, γc ] inside which shear bands are possible.

The value of σ y can be obtained easily by simulating the system in the elastic regime (D B p = 0) up to the threshold (W y (B) = 0), which corresponds to a state of the system characterized by eigenvalues β y 1 and β y 2 of tensor B, a state for which σ y can be calculated.

The different stationary state values of the shear stress could then obtained independently by continuing the simulation beyond the threshold in the plastic regime for each value of γ, waiting for the stationnary value of the system (dB/dt ≈ 0). The dynamic threshold σ d would then correspond to the limit of σ 12 for small γ. The critical shear rate γc would be obtained when the stress applied to the system in the stationary state would precisely correspond to the plastic threshold: σ stat 12 ( γc ) = σ y . Such a procedure is natural, but requires successive simulations of the system for a large number of γ values.

We have used a more direct approach [34] to obtain σ d and γc (see Appendix). This method relies on the description of the evolution of the system in terms of independent eigenvalues β 1 and β 2 of tensor B (see figure 10).

With the help of this procedure, the three observables which are important for the prediction of shear bands, σ y , σ d , and γc , are obtained directly without the need for simulating separately all the points along the stationary flow curve.

E. Spatial (1D) simulations of stationary flow regimes

As already mentioned, the discussion in the previous paragraph only provides necessary conditions for the appearance of shear bands. In the 1D simulations that will be discussed in the present section, flow inhomogeneities will indeed sometimes emerge in a full spatial simulation of our tensorial model in 1D spatial dimension plus time.

We simulate the full tensorial model in 1D using the equations of paragraph III C 3. The technical details of the numerical scheme can be found in [52]. From a numerical point of view, let us remark in particular that we have checked the the grid used in the discretisation of the equations is fine enough for all simulations presented here.

performing a 1D simulation as a set of partial differential equations, we would obtain the exact same results as in the previous Section.

In other words, since the parameters of the model do not vary in space, shear bands can only appear if initial conditions are, in one way or another, inhomogeneous.

Of course, as mentioned in the Introduction, inhomogeneities could appear in a natural way through an extra state variable coupled to the flow, such as the concentration. This variable could then vary in space and be coupled to the flow. Concerning concentration (a conserved variable), let us mention dilatancy phenomena, imagined for foam [27,28], observed experimentally [56] and interpreted in a geometrical manner [29,57]. Alignement (a non-conserved variable) is another possibility. It has been invoked in the case of wormlike micelle or rigid rod solutions [55].

Here, we focus on inhomogeneous static strain/stress initial conditions, without invoking additional variables, and we will show that they can induce the appearance of persistent inhomogeneities in the flow profile.

The reason for which these initial strain inhomogeneities can induce the appearance of blocked bands can be qualitatively undestood by considering the flow threshold K. Indeed, the stresses generated by the shear combine with the initial stress distribution due to strain inhomogeneities. Depending on its orientation, the initial stress thus precipitates or delays the triggering of the plastic flow.

Initial inhomogeneous strain distribution

First, the existence of stress inhomogeneities stored in the system before it is set into motion is physically well motivated. For instance, introducing a foam sample into an apparatus requires non-homogeneous flows. Inhomogeneous stresses will build up in the sample very likely, except if particular care is taken, such as a slow, in situ drying of an initially wet foam.

We will always assume that the initial state is at rest, that is, that the elastic stresses are at equilibrium in the sample. However, even when this equilibrium is imposed, there exist a large set of possible initial spatial distributions of stresses and strains. For example, if the system is invariant in the xz plane of the shearing walls, some components of the stress must be homogeneous. That is the case for σ xy , σ yy and σ yz . The other stress components, however, can freely vary as a function of y as long as they remain constant in each xz plane. It thus corresponds to a 1D inhomogeneity in the direction of the velocity gradient.

In this paragraph, we examine a very simple case of initial condition, with uniaxial extension along axis x, with B xx = B xx (y) and B yy = B zz = 1/ √ B xx . We used a simple monotonic function:

B xx = 1.1 + ǫ y β (1 -(1 -y) β ). ( 64
)
In practice, in order to prepare a sample in such a state, one must compress the foam in a non-homogeneous manner. Typically, a block of foam with a trapezoidal shape forced to take a rectangular shape will undergo this kind of strain inhomogeneity. In this context, we cannot comment on any relation between these strain inhomogeneities in our continuum model and the local structural disorder existing at the bubble level, as this disorder is averaged out in our continuum description. In particular, there is no clear structural interpretation of the amplitude ǫ of the strain inhomogeneities in the prepared sample.

Characterizing the inhomogeneous flows

A typical sequence of velocity profiles obtained in our numerical simulations displays as follows. The velocity profile is initially homogeneous. It remains homogeneous as long as the entire sample is in the elastic regime. The regions where the initial stress is the highest in the direction of the applied deformation reach the threshold first. The average shear rate being constant, this onset of creep leads both to a higher shear rate in the creeping regions and to a lower one in the others. The high shear rate then induces the saturation of the stress due to creep, and the shear becomes blocked in the region below the threshold. In the stationary regime, a blocked band coexists with a sheared band at the same shear stress.

In the corresponding transient regime, non-trivial phenomena may appear, especially at the boundary of the blocked zone. Transient negative local shear rates are observed due to stored elastic stresses.

Let us now address the characteristics of the stationary velocity profile, again from the behaviour of the local rheological model.

The first feature of interest is that in the flowing regions, the velocity profile is linear, that is, the shear rate is uniform. That can be understood in the following manner. All the regions which, in the stationary state, respond through a non-zero shear rate, correspond to a point located on the stationary flow curve in the β 1 -β 2 diagram of figure 10. Each point of this curve corresponds to a different shear stress. Thus, since each layer of the flow undergoes the same shear stress, they all actually correspond to the same point on the curve and thus respond through the same shear rate.

A second feature results from the fact that in space, while σ xy and σ yy are continuous, σ xx and σ zz can perfectly be discontinuous. That is precisely the case at the boundary between a shear and a blocked region. This is the flow counterpart of the discontinuity in the β 1 -β 2 diagram, between the points below the threshold and the point with a stationary shear rate that corresponds to the flowing region. Actually, the only coupling between the different layers comes from the fact that (i) σ xy and σ yy must be every where the same, and (ii) the integral of γ over the gap thickness is fixed by the imposed wall velocity. As a consequence, in an inhomogeneous flow, the organisation of the blocked and flowing layers is not unique: any permutation of the layers is actually possible. Again, initial conditions decide upon the particular structure adopted by the flow. Two initial conditions corresponding to permutated layers would lead to the same permutation in the stationary flow structure.

F. Parameters affecting the existence of blocked bands

In this section, we want to describe, within the parameter space (Ψ, K, a, b, Γ) the regions inside which shear bands are possible. These domains will be represented through sections in five different planes: (K, Γ), (Ψ, Γ), (a, b) and (Ψ, K). The results are presented in figures 11, 12, 13, 14, and 15.

As will be discussed below (paragraph V F 6), the choice of the initial conditions can have a crucial impact on the existence of shear bands. A complete investigation of the model would therefore require a very thorough exploration not only of the parameters (Ψ, K, a, b, Γ) but also of the shape and amplitude of the initial strain profile. In order to favour the appearance of shear bands without needing to refine the exploration of various strain profiles, we selected a very large amplitude ǫ = 500% for the shape mentioned in Eq. ( 64). This applies to Figs 11-15.

(K, Ψ) plane

In the (K, Ψ) plane, bands are predicted by the local model for large values of Ψ and K. Note that in the figure, zones where no bands can appear are denoted by blue triangles. Small values of Ψ correspond to a situation where the relaxation time ηs G in the absence of plasticity is far smaller than the relaxation time τ corresponding to the plasticity. It is thus a regime dominated by the fluid viscosity, where the creep plays no role. As Ψ increases, the creep becomes dominant, and shear bands can appear for lower values of the threshold (K) (figure 11).

As expected, the shear bands observed in the simulations appear only in regions authorized by the scalar model.

The green dots correspond to values for wich the scalar model allows the presence of shear bands, wich are not observed in the simulations for a specific set of initial conditions. As will be commented further on, the extent of this green zone depends on these initial conditions, demonstrating one of the main points of this work.

(K, Γ) plane

In the (K, Γ) plane, bands should be predicted for small values of Γ (due to the small velocities wich explore regions of the flow curve close to the origin in Γ), and for large values of K. Indeed, in that case the static threshold is large which favours bands since they are possible below this threshold (figure 12).

Concerning the relation between the values predicted for shear bands in the 0D model and the observations in the 1D simulations, the same remarks hold as for the previous paragraph.

(Ψ, Γ) plane

Observations in this plane corroborate the analysis in the two previous planes : bands are allowed (and are observed) for low Γ values and high Ψ values (figure 13).

(K, a) plane

Again bands appear for large values of K. The influence of the a parameter is far more subtle to assess, being related to non-trivial tensorial effects of the elastic (a) and plastic (b taken as 1 -a here) terms.

The same remarks hold concerning the correlation between the 0D model and 1D simulations.

(a, b) plane

Finally, in the (a, b) plane, one is again confronted with 3D effects which are difficult to discuss in intuitive terms (figure 15). The way the elasticity (parameter a) and the plastic deformation rate (parameter b) are coupled in a tensorial way affects the critical rate γc and can be enough to eliminate all possibilities of shear bands.

Dependance on the initial conditions

We have always observed that the regions in which blocked bands actually appeared in the 1D simulations are strictly included, as expected, in the regions authorized by the local rheological model. However, the boundary of these regions do not coincide. In fact, for the same values of the parameters, the extent of the banding zone depends crucially on the initial conditions, while always remaining in the region allowed by the rheological model. In other words, the behaviour of the system is history dependent, a feature realized independently in a recent work on a related tensorial model with plasticity [58].

To illustrate that, we have varied both the form and the amplitude of the spatial modulation of the initial deformation. In Figures 1112131415, the initial profile was given by the non-linear form of Eq. (64) with ǫ = 500%. By contrast, in Figure 16, for the upper graph we chose a simple sigmoidal profile centered around a selected altitude y 0 ,

B xx = 1 + ǫ y β y β 0 + 1 y β 0 + y β , (65)
and for the lower graph we chose a step-like function, both with ǫ = 10%. Comparing figures 12 and 16 shows that the parameter domain where shear bands actually appear can depend in a non-trivial manner not only on the shape but also on the amplitude of the initial strain profile.

Actually, we expect that a thorough exploration of the region authorized for the shear bands could be achieved through a very fine adjustment of the initial condition profile for each set of parameters.

VI. CONCLUSION

The present study on shear bands in liquid foams was conducted on a rheological model [34] whose predictions we here compare to existing rheological measurements under large amplitude oscillations (see Section IV C). The shear bands obtained with the model in parallel Couette geometry display several somewhat unusual features.

1. The shear bands depend on the initial conditions.

More generally, the stationary state is history dependent: only the flowing regions coincide with the simple linear velocity profile obtained in the case of a stationary homogeneous flow.

2. The response of the model in a stationary homogeneous flow is continuous when approaching zero shear rate, with no forbidden region below some finite shear rate.

3. The model does not contain any non-conserved In the upper graph we considered sigmoidal initial conditions, and in the lower graph step-like ones. In both cases, the amplitude of the strain inhomogeneities was reduced to 10% in amplitude, whereas Fig. 12 corresponded to 500% to enhance the effect.

(structural) order parameter.

This study thus shows that shear bands can arise naturally in a fully tensorial rheological model. This departs from most works in the shear banding community which put less emphasis on the tensorial character of the various models. Here, the appearance and persistence of the bands result from the combination of the initial conditions and the difference between the static and the dynamic flow thresholds (in shear geometry), which itself arises from the tensorial character of the model. The shear rate is discontinuous at the boundary between the flowing and blocked regions, but the value of the shear rate near the boundary as well as the band widths depend on the sample history and preparation.

FIG. 1 :

 1 FIG. 1: Schematic, scalar view of rheological models (deformation ε and stress σ) with only one internal degree of freedom, namely the deformation e of the spring. (a) General model: it has one (non-linear) spring and three creeping elements, each of which may have a flow threshold. (b) Foams and emulsions display some (large) deformation before creep is triggered. Hence, creep elements 1 and 2 cannot have a finite threshold. By contrast, element 3 does have a threshold.Here, for simplicity, we assume that the spring and both viscous elements are linear, and that element 3 cannot withstand any stress beyond σ3.

FIG. 2 :

 2 FIG. 2: Response of the model represented in Fig. 1b to a triangle wave deformation. (a) Deformation ε as a function of time. (b) In the limit η2 = 0, stress σ (solid line) and elastic stress G1 e transmitted by the spring (dashed line) as a function of time. When the threshold σ3 is reached, the spring relaxes with the time scale η1/G1 of the Voigt element. (c) Corresponding Lissajous representation of σ and G1 e. The periodic jump in stress (green segment) has the same amplitude as the initial jump (red segment). (d) In the limit η1 = 0, stress σ (solid line) and elastic stress G1 e transmitted by the spring (dashed line) as a function of time. (e) Corresponding Lissajous representation of σ and G1 e. The amplitude of the periodic jump in stress (green segment) is twice that of the initial jump (red segment).

FIG. 3 :

 3 FIG.3:(a) Burger model: schematic diagramme (left) and deformation response ε(t) to a step imposed stress σ(t). Spring G ′ 1 elongates immediately (red segment) while spring G1 responds with some delay (green curve) due to viscous element η1. Viscous element η3 gives rise to a constant additional deformation rate (blue segment). (b) Combined Bingham-Burger model: solid friction element σ3 (complemented by viscous element η ′3 ) is now included so as to provide additional creep above the stress threshold. We believe that some tensorial version of this model could mimic the rheological behaviour of a foam quite adequately. (c) Model studied in the present work. Apart from the additional viscosity η2 introduced in Fig.1, it represents the combined Bingham-Burger model when the (blue) viscous element has not been deformed yet, either at intermediate time scales when the (green) Voigt element has relaxed (the green and red springs then respond in series) or at short time scales when the Voigt element is still blocked (in which case only the red spring responds).

FIG. 4 :

 4 FIG.4: Simplified (scalar) picture of the main rheological parameters. G represents the elastic modulus and a the relative weight of the tensorial components of the elastic deformation, see Eqs. (27) and(28). The quantities σy and 1/τ constitute a scalar representation of the creep defined by D B p , and parameter b is the equivalent of a for creep, see equation(31). Finally, ηs is a viscosity that is independent of creep.

FIG. 5 :

 5 FIG. 5: Large amplitude oscillatory simulations: shape of Γ(t) versus σxy(t) (Lissajous) curves, obtained for a = b = 1/7, Ψ = 0.1 and K = 1. Top: range from ω = 0.01 to ω = 100 and Γ0 = 0.1 to Γ0 = 100. Bottom: zoom on a more restricted range of parameters.

8 FIG. 6 :

 86 FIG.6: Large amplitude oscillatory simulations: shape of σxy(t) versus Γ(t) (Lissajous) curves, obtained for a = b = 1/7, ω = 0.01, Ψ = 0.1 and K = 1. Top: transition from elastic to plastic behaviour as the amplitude is increased. At very large amplitudes, an overshoot is apparent like in continuous shear situations. Bottom: in a slightly plastic situation (amplitude 1.8), it takes several cycles before the system behaves in a periodic manner.

FIG. 7 :

 7 FIG. 7: Comparison of model (curves) with experiments (points). Normalized σ(t) curves for six values of the amplitude Γ0 ranging from 0.055 to 1.2, shifted vertically for clarity. Parameter values are a = 0.14, K = 0.04, b = 0.14, Ψ = 27. For each value of the amplitude Γ we had to select a different value for the frequency ω and for the modulus G. The (Γ, ω, G) values are: (0.055, 0.2, 282), (0.15, 0.4, 246), (0.25, 0.9, 186), (0.43, 0.5, 154), (0.72, 0.5, 146), (1.2, 0.4, 149).

FIG. 9 :

 9 FIG. 9: Top: typical stationary flow curve. Points corresponding to σ d , σy and γc are reported on the curve. Bottom: stress time evolution for different imposed shear rates below or above γc.

β 2 (FIG. 10 :

 210 FIG. 10: Form of the stored elastic deformation in the course of an experiment and in the stationary regime, for three different values of the shear rate γ. The axes are the first two eigenvalues, β1 and β2, of tensor B.

7 FIG. 11 :

 711 FIG. 11: Comparison of 0D and 1D simulations in plane (K, Ψ), using a = b = 1/7 and Γ = 0.05. Blue triangles indicate values for which the 0D model allows only uniform flow. Red squares indicate values for which shear banding was obtained in the 1D simulation for the initial conditions chosen. Green disks indicate additional values for which the 0D model allows shear banding.

7 FIG. 12 :

 712 FIG. 12: Comparison of 0D and 1D simulations in plane ( Γ, K), using a = b = 1/7 and Ψ = 0.1. Blue triangles indicate values for which the 0D model allows only uniform flow. Red squares indicate values for which shear banding was obtained in the 1D simulation for the initial conditions chosen. Green disks indicate additional values for which the 0D model allows shear banding.

K = 1 , 7 FIG. 13 :

 1713 FIG. 13: Comparison of 0D and 1D simulations in plane ( Γ, Ψ), using a = b = 1/7 and K = 1.0. Blue triangles indicate values for which the 0D model allows only uniform flow. Red squares indicate values for which shear banding was obtained in the 1D simulation for the initial conditions chosen. Green disks indicate additional values for which the 0D model allows shear banding.

FIG. 14 :

 14 FIG. 14: Comparison of 0D and 1D simulations in plane (K, a), using b = 1 -a, Γ = 0.05 and Ψ = 0.1. Blue triangles indicate values for which the 0D model allows only uniform flow. Red squares indicate values for which shear banding was obtained in the 1D simulation for the initial conditions chosen. Green disks indicate additional values for which the 0D model allows shear banding.

8 FIG. 15 :

 815 FIG. 15: Comparison of 0D (top) and 1D (bottom) simulations in plane (a, b), using Γ = 0.05 and Ψ = 0.1. Following the indications of Fig. 14 we chose K = 0.3 for 0D simulations and K = 0.8 for 1D simulations. The second diagonals (b = 1-a) in the present diagrams correspond to vertical lines in Fig. 14. Blue triangles indicate values for which the 0D model allows only uniform flow. Red squares indicate values for which shear banding was obtained in the 1D simulation for the initial conditions chosen. Green disks indicate additional values for which the 0D model allows shear banding.

7 FIG. 16 :

 716 FIG.16: Initial conditions dependency in the plane ( Γ, K), using a = b = 1/7 and Ψ = 0.1. In the upper graph we considered sigmoidal initial conditions, and in the lower graph step-like ones. In both cases, the amplitude of the strain inhomogeneities was reduced to 10% in amplitude, whereas Fig.12corresponded to 500% to enhance the effect.

Discussion of the conditions for inhomogeneous flowThe model that we simulate only contains material parameters that are homogeneous in the sample. Hence, if the initial conditions of the flow are also homogeneous, the entire evolution will remain homogeneous. Although

Acknowledgements

We warmly thank Florence Rouyer and her co-authors for transmitting the raw data from the large amplitude oscillatory experiments of Ref. [38] discussed in Section IV C.

Let us start from the point (β y 1 , β y 2 ) and follow the plasticity threshold W y (B) = 0 until the stationarity condition is fulfilled. This condition can be expressed using the following observation: in the stationary regime, there is no plastic flow in the vorticity direction [34]. In other words, the third eigenvalue of tensor G(B) is zero:

with β 3 = 1 β1β2 . We thus directly obtain the dynamic threshold (β d 1 , β d 2 ) of the system. We then follow the same stationarity condition g 3 (β 1 , β 2 ) = 0 until we reach the desired shear stress σ 12 = σ y . We thus directly obtain the stationary state (β cc 1 , β cc 2 ) that corresponds to the critical shear rate γc . In practice, we follow the threshold curve using Ŵy (β 1 , β 2 ) = W y (β 1 , β 2 , β 3 ) = 0 (with β 3 = 1 β1β2 ) by integrating the following differential system:

where the sign of ε Wy = ±1 is chosen in such a way as to follow the curve W y in the desired direction. Similarly, we follow the curve of stationary states, g 3 (β 1 , β 2 ) = G 3 (β 1 , β 2 , β 3 ) = 0 by integrating the following differential system:

where the sign of ε g3 = ±1 is chosen in such a way as to follow the curve g 3 = 0 in the desired direction.


Luca Fabbri
email: luca.fabbri@bo.infn.it

Stefano Vignolo
email: vignolo@diptem.unige.it

Diptem Sez

Metodi

Dirac fields in f (R)-gravity with torsion

Keywords: f (R)-gravity, Dirac field, conservation laws PACS number: 04.50.Kd, 04.62.+v

We study f (R)-gravity with torsion in presence of Dirac massive fields. Using the Bianchi identities, we formulate the conservation laws of the theory and we check the consistency with the matter field equations. Further, we decompose the field equations in torsionless and torsional terms: we show that the nonlinearity of the gravitational Lagrangian reduces to the presence of a scalar field that depends on the spinor field; this additional scalar field gives rise to an effective stress-energy tensor and plays the role of a scale factor modifying the normalization of Dirac fields. Problems for fermions regarding the positivity of energy and the particle-antiparticle duality are discussed.

Introduction

In the last thirty years, several extensions of General Relativity (GR) have been formulated in order to solve some problems left open by Einstein's theory at ultraviolet and infra-red scales; among all the developed extended theories of gravity, one of the simplest approaches is certainly given by the so-called f (R)-theory: it consists in relaxing the hypothesis that the gravitational Lagrangian be a linear function of the Ricci scalar R allowing a dependence on R more general than that of the Einstein-Hilbert action. This approach has acquired great interest in cosmology and astrophysics where f (R)-gravity turned out to be useful in addressing cosmological and astrophysical puzzles such as dark energy and dark matter; for example, it leads to alternative explanations of the accelerated behaviour of the today universe as well as the missing matter at astrophysical scales without postulating the existence of forms of energy and matter unknown so far. Further reasons to investigate extended theories of gravity come from the Field Theory; as it is well-known, in dealing with elementary particles the spin is as fundamental as the momentum. In order to take spin into account, first Cartan introduced torsion in the geometrical background, and then Sciama and Kibble took it to the framework of Einstein gravity following the scheme for which the energy is the source of curvature in the same way in which the spin is the source of torsion; in the resulting Einstein-Cartan-Sciama-Kibble (ECSK) theory, the gravitational fields are given by a metric tensor and a metric compatible connection with torsion, whose dynamics is described by a gravitational Lagrangian linear in the scalar curvature. In this scenario, f (R)-gravity with torsion has been developed as one of the simplest extensions of the ECSK theory, the one relying again on the idea of replacing the Einstein-Hilbert Lagrangian with a non-linear function. An important consequence of this fact is the presence of torsion even without spin as long as the stress-energy trace does not vanish [1,2,3,4]. Nonetheless, the matter field that best exploits the coupling between torsion and spin, and therefore the most general to study, is the spinor field, the simplest of which is the spin- 1 2 spinor field, already studied in the special case of masslessness by Rubilar [5]. In this paper, we consider f (R)-gravity with torsion coupled to the Dirac spinor field in its most general massive case; our approach will face the problem of the consistency of matter field equations with the conservation laws. In fact, the ECSK-like field equations of the theory involve the curvature and the torsion tensors which appear also in the well-known Bianchi identities; the latter are geometrical and they hold independently on the ECSK-like field equations through which they are converted into suitable conservation laws; these conservation laws must be satisfied once matter field equations are eventually given. A point that has to be stressed is that the obtained conservation laws reduce to the ones of the Einstein theory when the spin vanishes (for instance, the particles of a perfect fluid without spin follow the Levi-Civita geodesics, even if the geodesic structure of spacetime differs from the Levi-Civita one). However, in presence of spin all equations are more general than those of the Einstein theory, and for them the check of consistency may not be trivial (for example, for some of the higher-spin supersymmetric models, constraints known as subsidiary conditions need be added in order for such consistency to be obtained); this is a fundamental aspect for the viability of the theory. In particular, in the case of f (R)-gravity coupled to the Dirac field, we verify that the Dirac-like field equations are actually consistent with the aforementioned conservation laws by checking that the former automatically imply the validity of the latter; for this model, as far as our knowledge is concerned, this consistency check that has never been done in literature. After checking the consistency of the field equations, we shall further proceed by decomposing all of them into torsionless and torsional terms, showing that the contributions due to the non-linearity of the function f (R) are equivalent to the presence of a suitable scalar field depending on Dirac field alone. This scalar field is seen to modify the effective stress-energy tensor of the theory and it enters all equations as a scale factor modifying the normalization of Dirac field. We shall see that this fact has consequences for the Dirac field itself: indeed, the possibility to have the two opposite signs for the scale factor that multiplies the terms containing the Dirac field implies that it is possible to have both repulsive and attractive spin-torsion coupling within the Dirac-like field equations and to have both the negative and positive energy, respectively for particles and antiparticles, in the Einstein-like field equations; in particular, this last point is important for the problem of the positivity of the energy of particles and antiparticles, because it solves it in the classical domain. On the same footing, this situation leaves open the path to the possibility of having a particle-antiparticle duality lost in gravitational physics.

Geometrical preliminaries

Throughout the paper we indicate spacetime indices by Latin letters. A metric tensor on the spacetime M is denoted by g ij and a connection by Γ h ij ; metriccompatible connections are those whose covariant derivative applied on the metric tensor vanishes. Any connection Γ h ij defines the associated covariant derivative

∇ ∂ i ∂ j = Γ h ij ∂ ∂x h (2.1)
as well as the corresponding torsion and Riemann curvature tensors expressed as

T h ij = Γ h ij -Γ h ji (2.2a) R h kij = ∂ i Γ h jk -∂ j Γ h ik + Γ h ip Γ p jk -Γ h jp Γ p ik (2.2b)
where contractions

T i = T j ij , R ij = R h ihj and R = R ij g ij
are called respectively the torsion vector, the Ricci tensor and the Ricci scalar curvature. The torsion and curvature tensors are used to express the commutator of covariant derivatives as

[∇ i , ∇ j ]V k = -T h ij ∇ h V k -R a kij V a (2.3)
for any generic vector V k . By considering the commutators of commutators in cyclic permutation one can obtain the so-called Bianchi identities

∇ c T h ij + ∇ i T h jc + ∇ j T h ci -T a ij T h ca -T a jc T h ia -T a ci T h ja - -R h cij -R h ijc -R h jci = 0 (2.4a) ∇ c R p kij + ∇ i R p kjc + ∇ j R p kci -T a ij R p kca -T a jc R p kia -T a ci R p kja = 0 (2.4b)
whose usefulness will be clear in the following. Given a metric tensor g ij , every metric g-compatible connection can be decomposed in the form [6,7]

Γ h ij = Γ h ij -K h ij (2.5)
where

K h ij = 1 2 -T h ij + T h j i -T h ij (2.6)
is called the contorsion tensor and Γ h ij is the Levi-Civita connection induced by the metric g ij ; the contorsion tensor satisfies the antisymmetry property given by K jh i = -K jh i amounting to the metric-compatibility of every connection of the form (2.5). The contorsion has one contraction K ij i = K j , from which we have also the identity K i = -T i . With the contorsion we can decompose the covariant derivative of the full connection into covariant derivative of the Levi-Civita connection ∇ plus contorsional contributions; this gives the following representation of the Ricci curvature

R ij = Rij + ∇j K h hi -∇h K h ji + K p ji K h hp -K p hi K h jp (2.7)
with the Ricci curvature of the Levi-Civita connection Rij and contorsional terms.

In the next sections we shall consider f (R)-gravity coupled with Dirac fields: as it is well known, the most suitable variables to describe fermion fields interacting with gravity are tetrad fields. Tetrad fields possess Lorentz indices denoted by Greek letters as well as spacetime indices. They are defined by e µ = e µ i dx i together with their dual e µ = e i µ ∂ ∂x i where e j µ e µ i = δ j i and e j µ e ν j = δ ν µ and spin-connections 1-forms ω µ ν = ω µ i ν dx i . Metric compatibility conditions are assumed and they are defined by the requirement that the covariant derivatives of tetrads and Minkowskian metric vanish, respectively implying that

Γ h ij = ω µ i ν e h µ e ν j + e h µ ∂ i e µ j (2.8)
and the antisymmetry ω µν i = -ω νµ i of the spin connection. The torsion and curvature tensors induced by the assignment of a spin-connection together with a tetrad field are expressed as

T µ ij = ∂ i e µ j -∂ j e µ i + ω µ i λ e λ j -ω µ j λ e λ i (2.9a) R µν ij = ∂ i ω µν j -∂ j ω µν i + ω µ i λ ω λν j -ω µ j λ ω λν i (2.9b)
and their relationships with the world tensors defined in equations (2.2) are given by the relationships

T h ij := T α ij e h α and R h kij = R µ ij ν e h µ e ν k respectively.

f (R)-gravity with torsion and conservation laws

f (R)-gravity with torsion can be formulated in the metric-affine approach [1] or in the tetrad-affine one [2]. In the first case, the gravitational dynamical fields are represented by the metric g and a metric compatible connection Γ on the spacetime manifold M . In the second case, the gravitational dynamical fields are given by a tetrad field e µ i and a spin-connection ω µν i , still defined on M . Field equations are derived variationally through a Lagrangian function of the kind

L = f (R) -L m (3.1)
where f (R) is a real function of the Ricci curvature scalar R written in terms of the metric and connection, or equivalently tetrad and spin-connection, and L m indicates a suitable matter Lagrangian.

In general the metric-affine formulation is preferred when the connection does not couple with matter, namely when the matter Lagrangian is independent of the connection. In such a circumstance, the field equations of the theory are [1]

f (R)R ij - 1 2 f (R)g ij = Σ ij (3.2a) T h ij = - 1 2f (R) ∂f (R) ∂x p δ p i δ h j -δ p j δ h i (3.2b)
where Σ ij := 1 √ |g| δLm δg ij plays the role of the stress-energy tensor. Equations (3.2a) are by construction symmetric. The symmetrization symbol for the Ricci tensor R ij is omitted because any connection satisfying eqs. (2.5), (2.6) and (3.2b) gives rise to a Ricci tensor that is automatically symmetric [1,3,4]. We also notice that, if f (R) = kR, we have non-vanishing torsion even if the spin density is zero. On the other hand, in case of coupling with spin matter, for instance spinor fields, the tetrad-affine formulation is more suitable. The corresponding field equations have been derived in [2] (see also [5]) and turn out to be

f (R)R λσ µσ e i λ - 1 2 e i µ f (R) = Σ i µ (3.3a) T α ts = 1 f (R) 1 2 ∂f (R) ∂x p + S σ
f (R)R -2f (R) = Σ (3.4)
between the curvature scalar R and the trace Σ of the stress-energy tensor. From (3.4) we see that when the trace Σ is only allowed to be constant (where the constant can assume any value, in particular zero) then also R is forced to be constant; in such a circumstance, the present theory amounts to an Einsteinlike (if S α ts = 0, i.e. ω-independent matter) or an ECSK-like (if S α ts = 0, i.e. ω-dependent matter) theory with cosmological constant [1,2,3,5]. The previous statement holds with the exception of the particular case given, when Σ = 0, by the relationship f (R) = kR 2 , because under these conditions equation (3.4) is a trivial identity imposing no restrictions on the curvature scalar R. Therefore, from now on we shall systematically assume that Σ is not constant as well as that f (R) = kR 2 . Moreover, we shall assume that the relation (3.4) is invertible so that the Ricci curvature is a suitable function of Σ, namely R = F (Σ). With the stated assumptions in mind, defining the world tensors R i j := R λσ µσ e i λ e µ j , Σ i j := Σ i µ e µ j , T h ij := T α ij e h α and S h ij := S α ij e h α , we can rewrite equations (3.3) in their equivalent general covariant form

R ij = 1 f (F (Σ)) 1 4 f (F (Σ))F (Σ) -Σ g ij + Σ ij (3.5a) T h ij = - 1 2f (F (Σ)) ∂f (F (Σ)) ∂x p + S q pq δ p i δ h j -δ p j δ h i + 1 f (F (Σ)) S h ij (3.5b)
giving the Ricci curvature tensor and torsion tensor in terms of the energy and spin densities. In equations (3.5a) one has to distinguish the order of the indices since in general the tensors R ij and is not symmetric if S h ij = 0 while Σ ij is not symmetric because it is obtained by variation with respect to the tetrad fields. It is also important to notice that, as it happens for the purely metric case, any f (R) theory with torsion (with or without spin) is equivalent to a corresponding torsional scalar-tensor theory with Brans-Dicke parameter ω 0 = 0 and with specific potential for the scalar field directly induced by the function f (R) [1,2]. Moreover, recalling the expression (2.6) for the contorsion tensor and making use of equations (3.5b), we obtain the following representations

K h ij = K h ij + Ŝ h ij (3.6)
where one has

Ŝ h ij := 1 2f -S h ij + S h j i -S h ij (3.7a) K h ij := -Tj δ h i + Tp g ph g ij (3.7b) Tj := 1 2f ∂f ∂x j + S q jq (3.7c)
for the decomposed tensors. Inserting equations (3.6) and (3.7) into equations (2.7), we get the following useful decomposition of the Ricci tensor

R ij = Rij + ∇j K h hi + ∇j Ŝ h hi -∇h K h ji -∇h Ŝ h ji + K p ji K h hp + K p ji Ŝ h hp + Ŝ p ji K h hp + Ŝ p ji Ŝ h hp -K p hi K h jp -K p hi Ŝ h jp -Ŝ p hi K h jp -Ŝ p hi Ŝ h jp (3.8)
which we shall employ later. Now, we shall use the Bianchi identities (2.4) to work out the conservation laws for the above outlined theory. To see this point, let us consider the Bianchi identities in their fully contracted form as

∇ c R ci - 1 2 g ci R = + 1 2 T cja R cjai + T ica R ca (3.9a) ∇ c T ijc + g ci T j -g cj T i = -T ija T a -R ij + R ji (3.9b)
together with the field equations (3.5) written in the equivalent form

f R ij - 1 2 g ij f = Σ ij (3.10a) f T tsa -T t g sa + T s g ta = ∇ t f g sa -∇ s f g ta + S tsa . (3.10b)
The divergences of equations (3.10) read respectively as

∇ i f R ij + f ∇ i R ij - 1 2 g ij R = ∇ i Σ ij (3.11a) ∇ a f T tsa -T t g sa + T s g ta + f ∇ a T tsa -T t g sa + T s g ta = = [∇ s , ∇ t ]f + ∇ a S tsa (3.11b)
and by inserting the content of the fully contracted Bianchi identities (3.9) into equations (3.11), we get

∇ i f R ij + 1 2 f T cka R ckaj + f T jca R ca = ∇ i Σ ij (3.12a) ∇ a f T tsa -T t g sa + T s g ta -f T tsa T a + f R st -R ts = = [∇ s , ∇ t ]f + ∇ a S tsa (3.12b)
so that evaluating explicitly the term [∇ s , ∇ t ]f , the expressions (3.12) can be further simplified as

∇ i f R ij + 1 2 f T cka R ckaj + f T jca R ca = ∇ i Σ ij (3.13a) ∇ t f T s -∇ s f T t -f T tsa T a + f R st -R ts = ∇ a S tsa . (3.13b)
At this point, the field equations (3.10) can again be used to calculate: the antisymmetric part of (3.10a)

f R ij -R ji = Σ ij -Σ ji (3.14a)
the contractions of (3.10a) with torsion vector and torsion

f R ij T i - 1 2 fT j = Σ ij T i (3.14b) f R ij T bij - 1 2 fT b = Σ ij T bij (3.14c)
the contractions of (3.10b) with curvature and torsion

1 2 f R tsai T tsa + T t f R ti = -R ti ∇ t f + 1 2 S tsa R tsai (3.14d) f T tsa T a = T s ∇ t f -T t ∇ s f + S tsa T a (3.14e)
and by finally inserting the content of all equations (3.14) into equations (3.13), we obtain the following conservation laws

∇ a Σ ai + T a Σ ai -Σ ca T ica - 1 2 S sta R stai = 0 (3.15a) ∇ h S ijh + T h S ijh + Σ ij -Σ ji = 0 (3.15b)
under which the stress-energy and spin density tensors of the matter fields must undergo once the matter field equations are assigned.

It is worth noticing that in the case of coupling with matter without spin, equations (3.15) reduces to the usual conservation laws holding in general relativity. Indeed, in metric-affine formulation with no spin equations (3.15b) are identities. At the same time, working out equation (3.15a) and taking equations (2.5) and (2.6) into account, we have

0 = ∇i Σ ij -K i ih Σ hj -K j ih Σ ih + T h Σ hj -T j ih Σ ih = ∇i Σ ij . (3.16)
This last result has been already proved in a different way [8,9,10], exploiting the dynamical equivalence of metric-affine f (R)-gravity without spin with certain Brans-Dicke-like theories, for which the standard conservation laws given in equation (3.16) hold [11]. The proof given above is more direct and geometrical.

Next we shall investigate more in detail the coupling with Dirac fields.

Coupling with Dirac matter fields and consistency

Let us consider f (R)-gravity coupled with a Dirac field. The matter Lagrangian function is the Dirac one

L m = i 2 ψγ i D i ψ -D i ψγ i ψ -m ψψ (4.1)
where

D i ψ = ∂ψ ∂x i + ω µν i S µν ψ and D i ψ = ∂ ψ ∂x i -ψω µν i
S µν are the covariant derivatives of the Dirac fields, S µν = 1 8 [γ µ , γ ν ], γ i = γ µ e i µ with γ µ denoting Dirac matrices and where m is the mass of the Dirac field. We use the notation for which

γ µ γ ν γ λ = γ µ η νλ -γ ν η µλ + γ λ η µν + i µνλτ γ 5 γ τ (4.2)
which implicitly defines the γ 5 matrix and which entails the identities

[γ µ , [γ ν , γ λ ]] = 4 γ λ η µν -γ ν η µλ (4.3a) and {γ µ , [γ ν , γ λ ]} = 4i µνλτ γ 5 γ τ (4.3b)
we will employ in the following. By varying (4.1) with respect to ψ, we obtain the matter field equations

iγ h D h ψ + i 2 T h γ h ψ -mψ = 0 (4.4)
which generalize the usual Dirac equations. Due to equations (4.4), it is easily seen that the stress-energy and spin density tensors, derived as in section 3, are respectively given by

Σ ij = i 4 ψγ i D j ψ -D j ψγ i ψ (4.5a) S h ij = i 2 ψ γ h , S ij ψ ≡ - 1 4 η µσ σνλτ ψγ 5 γ τ ψ e h µ e ν i e λ j (4.5b)
where we have used (4.3b) to show the equivalence of the two forms of the spin density, which is then manifestly completely antisymmetric and thus irreducible. Now, it is possible to calculate the divergences of the quantities (4.5) as

D i Σ ij = i 4 D i ψγ i D j ψ + ψγ i D i D j ψ -D i D j ψγ i ψ -D j ψγ i D i ψ (4.6a) D h S ijh = i 2 D h ψγ h S ij ψ + ψγ h S ij D h ψ + D h ψS ij γ h ψ + ψS ij γ h D h ψ (4.6b)
and by adding and removing the terms ψγ i D j D i ψ and D j D i ψγ i ψ in equations (4.6a) as well as the terms ψS ij γ h D h ψ and D h ψγ h S ij ψ in equations (4.6b), we obtain the following expressions

D i Σ ij = i 4 ( ψγ i [D i , D j ]ψ -[D i , D j ] ψγ i ψ + +D i ψγ i D j ψ + ψγ i D j D i ψ -D j D i ψγ i ψ -D j ψγ i D i ψ) ( 4 . 7 a ) D h S ijh = i 2 ( ψ[γ h , S ij ]D h ψ+D h ψ[S ij , γ h ]ψ+2D h ψγ h S ij ψ+2 ψS ij γ h D h ψ). (4.7b)
By employing the matter field equations (4.4), we get

D i Σ ij = i 4 ( ψγ i [D i , D j ]ψ -[D i , D j ] ψγ i ψ + D j ψγ h ψT h -T h ψγ h D j ψ) (4.8a) D h S ijh = i 2 ( ψ[γ h , S ij ]D h ψ + D h ψ[S ij , γ h ]ψ -T h ψ{γ h , S ij }ψ). (4.8b)
Finally, by making use of the commutator of covariant derivatives on spinors and the algebraic identities given by (4.3a), we have

D i Σ ij = -T i Σ ij -T ijk Σ ik + 1 2 S abi R abij (4.9a) D h S ijh = -T h S ijh + (Σ ji -Σ ij ) ( 4 . 9 b )
proving the complete consistency between the matter field equations (4.4) and the corresponding conservation laws (3.15). We now discuss more in detail the Einstein-like equations of the theory, separating the Levi-Civita contributions from the torsional ones. To this end, due to the total antisymmetry of the spin density, we note that the Ricci tensor and the curvature scalar can be expressed as

R ij = Rij -2 ∇j Ti -∇h T h g ij + 2 Ti Tj -2 Th T h g ij -∇h Ŝ h ji -Ŝ p hi Ŝ h jp (4.10a)
and also

R = R -6 ∇i T i -6 Ti T i -Ŝ p hi Ŝi h p (4.10b)
where one has Ti =

1 2f ∂f ∂x i e Ŝ h ij = -1 2f S h ij (compare with equations (3.7) and (3.8)
). Taking the decompositions (4.10) into account, from equations (3.5a) we derive the Einstein-like equations Rij -

1 2 Rg ij = 1 ϕ Σ ij + 1 ϕ 2 - 3 2 ∂ϕ ∂x i ∂ϕ ∂x j + ϕ ∇j ∂ϕ ∂x i + 3 4 ∂ϕ ∂x h ∂ϕ ∂x k g hk g ij -ϕ ∇h ∂ϕ ∂x h g ij -V (ϕ)g ij + ∇h Ŝ h ji + Ŝ p hi Ŝ h jp - 1 2 Ŝ p hq Ŝq h p g ij (4.11)
where the scalar field ϕ := f (F (Σ)) (4.12)

and the effective potential

V (ϕ) := 1 4 ϕF -1 ((f ) -1 (ϕ)) + ϕ 2 (f ) -1 (ϕ) (4.13)
have been introduced [1,2]. Splitting the symmetric part of (4.11) from the antisymmetric one, we get the following two sets of equations Rij -

1 2 Rg ij = 1 ϕ Σ (ij) + 1 ϕ 2 - 3 2 ∂ϕ ∂x i ∂ϕ ∂x j + ϕ ∇j ∂ϕ ∂x i + 3 4 ∂ϕ ∂x h ∂ϕ ∂x k g hk g ij -ϕ ∇h ∂ϕ ∂x h g ij -V (ϕ)g ij + Ŝ p hi Ŝ h jp - 1 2 Ŝ p hq Ŝq h p g ij (4.14a) 1 ϕ Σ [ij] + ∇h Ŝ h ji = 0 (4.14b)
and it is a straightforward matter (left to the reader) to verify that equations (4.14b) are equivalent to equations (4.9b). In other words, the matter field equations (4.4) imply equations (4.14b) and the significant part of the Einstein-like equations reduces to the equations (4.14a) only. The latter can be worked out by giving an explicit representation of the covariant derivative of spinor fields. To this end, we denote by D the spinorial covariant derivative induced by the Levi-Civita connection. Then, we have

D i ψ = Di ψ - 1 4 Kijh + Ŝijh γ h γ j ψ (4.15a) D i ψ = Di ψ + 1 4 ψ Kijh + Ŝijh γ h γ j (4.15b)
which can be used together with of (2.5), (3.6), (3.7) and (4.5b), to obtain the following identities

ψγ i D j ψ -D j ψ γ i ψ = ψγ i Dj ψ -( Dj ψ)γ i ψ - i 8ϕ ( ψγ 5 γ τ ψ)( ψγ 5 γ τ ψ)g ij + i 8ϕ ( ψγ 5 γ i ψ)( ψγ 5 γ j ψ) + i 2ϕ ∂ h ϕ ν σατ e h ν e σ i e α j ( ψγ 5 γ τ ψ) (4.16) Ŝ p hi Ŝ h jp = - 1 32(ϕ) 2 ( ψγ 5 γ τ ψ)( ψγ 5 γ τ ψ)g ij + 1 32(ϕ) 2 ( ψγ 5 γ i ψ)( ψγ 5 γ j ψ) (4.17) - 1 2 Ŝ p hq Ŝq h p g ij = 3 64(ϕ) 2 ( ψγ 5 γ τ ψ)( ψγ 5 γ τ ψ)g ij (4.18)
and finally by inserting (4.16), (4.17) e (4.18) into (4.14) we get the expression for the symmetrized Einstein-like equations Rij -

1 2 Rg ij = 1 ϕ i 4 ψγ (i Dj) ψ -D(j ψ γ i) ψ + 1 ϕ 2 - 3 2 ∂ϕ ∂x i ∂ϕ ∂x j + ϕ ∇j ∂ϕ ∂x i + 3 4 ∂ϕ ∂x h ∂ϕ ∂x k g hk g ij -ϕ ∇h ∂ϕ ∂x h g ij -V (ϕ)g ij + 3 64ϕ 2 ( ψγ 5 γ τ ψ)( ψγ 5 γ τ ψ)g ij . (4.19)
Next, we have that Dirac equations can be decomposed in torsionless and torsional terms; indeed, by using the identity ( ψγ 5 γ τ ψ)γ 5 + ( ψγ τ ψ)I γ τ ψ ≡ 0 together with the Fierz rearrangement, we get the following equivalent form for Dirac equations

iγ h Dh ψ - 1 ϕ 3 16 ( ψψ) + i(i ψγ 5 ψ)γ 5 ψ -mψ = 0 (4.20)
where the coupling between the spin density and the completely antisymmetric irreducible part of torsion generates spinorial self-interactions whose potentials can be recognized to be of the Nambu-Jona-Lasinio type [12,13].

As a final step, by making use of equations (4.20) the symmetrized Einstein-like equations can be written in the equivalent form

Rij = 1 ϕ i 4 ψγ (i Dj) ψ -D(j ψ γ i) ψ - 1 ϕ m 4 ( ψψ)g ij + 1 ϕ 2 - 3 2 ∂ϕ ∂x i ∂ϕ ∂x j + ϕ ∇j ∂ϕ ∂x i + 1 2 ϕ ∇h ∂ϕ ∂x h g ij + V (ϕ)g ij (4.21)
in which all torsional contributions have disappeared. We observe that the nonlinearity of the f (R) function gives rise to the appearance of the scalar field we indicated with ϕ = ϕ( ψψ) in the field equations (4.20) and (4.21). The role of the scalar field ϕ is two-fold: on the one hand, it provides additional terms in the effective stress-energy tensor in the Einstein-like field equations (4.21); on the other hand, it is a scale factor modifying the normalization of the Dirac field contributions in both equations (4.20) and (4.21). Furthermore, we also remark that the scalar ϕ may assume both positive and negative values, but it cannot switch between them due to the fact that it cannot vanish (see for instance equations (4.20) and (4.21)). This has again two implications: it may determine both repulsive and attractive potentials by controlling the sign of the self-interacting term in equations (4.20); it controls the sign of the stress-energy tensor of the Dirac field in equations (4.21). This means that both particles and antiparticles may have positive energies in classical gravity.

To explain this point, consider the circumstance in which f (R) = R -εR 2 , where ε is a fixed positive parameter. In this situation, we have ϕ = 1 + ε m( ψψ), which clearly shows that the sign of ϕ depends on the scalar ψψ built from the spinor field: for particles, we have ψψ 0 so that ϕ 0 and the energy is positive; for antiparticles, we have ψψ 0 so that for ψψ 1 εm it is ϕ 0 and the energy is negative, but for ψψ 1 εm it is ϕ 0 and the energy is positive. The consequence of this fact is that, for extremely high-density states, particles and antiparticles have a similar behaviour both with positive energies, but as the density decreases, their behaviour differentiates, and while particles maintain positive energy, we have that antiparticles would turn out to have negative energy. This mechanism accounts for a spontaneous symmetry breaking that takes place in classical physics. For example, in the instance of cosmological applications, we have that although in the beginning of the universe the particle-antiparticle duality and the positivity of their energies might have been an obvious fact, with the following expansion particles and antiparticles lose their symmetry, and whereas the former are still allowed, the latter are forced to disappear: that the particle-antiparticle duality be broken is not a new feature for torsional gravitational theories [14,15,16]; on the other hand torsional non-linear f (R)-gravity appear to amplify such a mechanism, as we have just discussed here, with no need of any quantization whatsoever. As a conclusive remark, we stress a consequence of having associated particles and antiparticles respectively with positive and negative values of the scalar field ϕ: recalling that the scalar ϕ is allowed to have both positive and negative values but it cannot switch sign, then the model is allowed to describe both particles and antiparticles but it cannot provide a mechanism for transmuting them into each other by means of gravitational effects alone. As the theory describes two sub-theories, one for positive and one for negative ϕ, without the possibility to turn one into the other, then the Dirac matter field gives two kind of solutions, particles and antiparticles, without the possibility to change into one another just in terms of gravitational physics, compatibly with all observations made so far.

Conclusions

In this paper, we have considered f (R)-theories of gravitation with Ricci scalar written in terms of connections having both metric and torsional degrees of freedom. We have studied them coupled to Dirac matter fields: we have seen that for general f (R) functions with R depending on metric and torsion, the corresponding field equations transform the Bianchi geometrical identities into conservation laws; these conservation laws are satisfied for the stress-energy and spin density tensors of the Dirac field, once Dirac matter field equations are used. After the set of field equations has been checked for consistency, we have decomposed all equations in terms of torsionless and torsional contributions. We have seen that within the fermionic field equations the self-interacting term is given by a potential of the Nambu-Jona-Lasinio form scaled in terms of a scalar field ϕ; the Einstein-like field equations are such that the fermionic contributions are also scaled in terms of the same scalar field ϕ, and furthermore an additional stress-energy tensor for the scalar field ϕ itself also appears. In this situation, we have discussed peculiar circumstances given by specific examples of ϕ for which both particles and antiparticles have positive energies, displaying a duality in the beginning of the universe that is eventually lost in further stages of the cosmic evolution. It is remarkable that the positivity of the energy for both particles and antiparticles is ensured and the particle-antiparticle duality is gravitationally spontaneously broken with a classic approach. The fact that no quantization, let alone quantization with anticommutators, is considered at all is important for the foundations of quantum field theory in view of the problems quantum field theory is known to have. In fact, the concern about divergences that are introduced as a byproduct of the quantization of fields can be circumvented by avoiding the quantization of fields in the first place, as long as the problems of the positivity of energies and the particle-antiparticle duality solved by quantizing matter fields find solution in an alternative way. One of these alternative is provided here.

  are the stress-energy and spin density tensors of the matter field. From equation (3.3b) it is seen that, in presence of ω-dependent matter, there are two sources of torsion: the spin density S α ts and the non-linearity of the gravitational Lagrangian. By considering the trace of equations (3.2a) or (3.3a) we obtain the relation

				pσ	(δ p s e α t -δ p t e α s ) + S α ts	(3.3b)
	where Σ i µ := -1 2e	∂Lm ∂e µ i	and S α ts := 1 2e	∂Lm ∂ω µν i	e µ t e ν s e α i