Swantje ; Bargmann
email: swantje.bargmann@uni-dortmund.de

Magnus Ekh

Kenneth Runesson

Bob Svendsen

Modeling of polycrystals with gradient crystal plasticity -A comparison of strategies

Keywords: computational mechanics, grain size, Hall effect computational mechanics, grain size, Hall effect crystal plasticity, gradient hardening, grain modeling, dual mixed finite element, polycrystalline material, constitutive behavior

Modeling of polycrystals with gradient crystal plasticity -A comparison of strategies

Introduction

It is well-known that the macroscopic behavior of a polycrystalline material is influenced by the size and morphology of the grains, the volume fraction of different phases, and the subgrain processes, for example. In order to include micro-effects in a macroscopic model, a complete scale separation is often assumed. With the help of a computational homogenization method the response of a Representative Volume Element (RVE), subjected to suitable boundary conditions that represent the macroscopic deformation (prolongation conditions), is coupled to the macroscopic level. However, conventional models of crystal plasticity, that are often used to model grain behavior, do not contain intrinsic material length-scales and are, therefore, not capable of capturing macroscopic size dependent effects. Thus, models including plastic strain gradients have been introduced in order to overcome this drawback of ordinary crystal plasticity theories. This contribution discusses different approaches to microstructure polycrystal models which include size dependence.

One size effect in metals is known as the Hall-Petch effect, i.e. the influence of the grain size on the macroscopic stress-strain characteristics. It is one of the important unresolved issues in computational material modeling. Hall [START_REF] Hall | The deformation and ageing of mild steel: III discussion of results[END_REF] and Petch [START_REF] Petch | The cleavage strength of polycrystals[END_REF] independently studied different behaviors of steel, but eventually obtained the same result. While Hall [START_REF] Hall | The deformation and ageing of mild steel: III discussion of results[END_REF] focused on factors influencing the mechanical properties of mild steels, Petch [START_REF] Petch | The cleavage strength of polycrystals[END_REF] studied brittle failure of steels. Nevertheless, Petch's results from 1953 are in excellent agreement with the ones of Hall which had been published two years earlier. The fact that the (current) flow stress σ (ǫ) depends on the mean grain size d became known as the Hall-Petch effect:

σ y (ǫ) = σ y,0 (ǫ) + k (ǫ) √ d . (1)
Both parameters, σ y,0 and k, may depend on the strain level ǫ due to hardening of the material.

In the literature, different models have been developed in order to capture size dependence effects in crystals. These include thermodynamic approaches like [7, 13, 14, 17-19, 22, 24, 25, 27] as well as dislocation-field-theory-based approaches like [START_REF] Bayley | A three-dimensional dislocation field crystal plasticity approach applied to miniaturized structures[END_REF] and [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Despite the large amount of models, there only exist few classifications or analyses of the common features and major differences. Examples which discuss similarities between Cosserat, strain gradient crystal plasticity and the statistical theory of dislocations in the case of single slip are the works of Forest and Sedláček [START_REF] Forest | Plastic slip distribution in two-phase laminate microstructures: dislocationbased versus generalized-continuum approaches[END_REF][START_REF] Forest | Some links between Cosserat, strain gradient crystal plasticity and the statistical theory of dislocations[END_REF]. A nonlocal plasticity theory is compared with a discrete dislocation model in [START_REF] Bassani | Plastic flow in a composite: a comparison of nonlocal continuum and discrete dislocation predictions[END_REF] and [START_REF] Bittencourt | A comparison of nonlocal continuum and discrete dislocation plasticity predictions[END_REF]. Numerical results are presented for single slip in a two-dimensional composite material subject to plane strain simple shear. A boundary value problem for a single crystal undergoing small deformations is discussed for a phase field dislocation theory and Gurtin's theory (cf. [START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF]) in [START_REF] Hunter | Direct calculations of material parameters for gradient plasticity[END_REF]. This contribution aims at comparing the three models presented in [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] and a combination of [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] and [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF].

The first model to be investigated in this contribution is the one introduced by Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] and continued by Bayley et al. [START_REF] Bayley | A three-dimensional dislocation field crystal plasticity approach applied to miniaturized structures[END_REF]. As opposed to the two models mentioned below, it does not originate from a thermodynamic setting. 1 Rather, it is based on physical considerations involving dislocation densities. An extended slip law, which also accounts for statistically stored dislocation densities, is introduced as well. Besides the mechanical displacements, the model of Geers and co-workers is based on the fields of geometrically necessary and statistically stored dislocation densities, which are used to affect crystallographic slip resistance in an extended slip law.

The second model is the one introduced in Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], which is an extension of the model for isotropic plasticity proposed by Svedberg and Runesson [START_REF] Svedberg | A thermodynamically consistent theory of gradient-regularized plasticity coupled to damage[END_REF]. In [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], a formulation is proposed which models the size effect by assuming that the hardening stress in the yield function depends on both the associated hardening strain (equal to the accumulated plastic slip) and its Hessian along the slip direction similarly to models for isotropic plasticity originally proposed by Aifantis [START_REF] Aifantis | On the microstructural origin of certain inelastic models[END_REF]. Similar formulations can be found in Gurtin [START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF] and Liebe and Steinmann [START_REF] Liebe | Theory and numerics of a thermodynamically consistent framework for geometrically linear gradient plasticity[END_REF]. The thermodynamic setting of Liebe and Steinmann [START_REF] Liebe | Theory and numerics of a thermodynamically consistent framework for geometrically linear gradient plasticity[END_REF] and Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] is identical. However, in contrast to [START_REF] Liebe | Theory and numerics of a thermodynamically consistent framework for geometrically linear gradient plasticity[END_REF], where the gradient hardening is isotropic, Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] obtain anisotropic hardening by restricting the gradient effect to each slip direction. Note that a similar strategy, also thermodynamically based, is followed by Vrech and Etse [START_REF] Vrech | FE approach for thermodynamically consistent gradient-dependent plasticity[END_REF] for small strains. Gurtin's formulation [START_REF] Gurtin | A gradient theory of small-deformation isotropic plasticity that accounts for the Burgers vector and for dissipation due to plastic spin[END_REF] is motivated from a so-called microforce balance instead of introducing the As a third model, the model of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] is in this paper extended by an idea of Lele and Anand [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF]. The approach of Lele and Anand [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] is formulated for small strains and isotropic materials, with a similar approach to the one of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]. Their idea is to introduce a slip rate with a similar structure to the one of [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]; however, including an additional higher-order dissipative term. Thus, as a third model, we investigate the model of [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] extended by this particular term.

Since we discuss large deformation theories and the field equations are highly coupled, a numerically efficient solution algorithm is needed. A dual mixed finite element algorithm based on [START_REF] Svedberg | An algorithm for gradient-regularized plasticity coupled to damage based on a dual mixed FE-formulation[END_REF] is applied to all three models, rendering a coupled two-field problem. The coupling effect is fully taken into account as the system of equations is solved monolithically. The implicit backward Euler scheme is used for the time integration. The resulting system of equations is then solved using a three-level iteration strategy based on a generalized form of domain decomposition, whereby each grain constitutes a subdomain, as suggested in [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]. Due to the fact that each grain can be considered independently (with prescribed boundary conditions), this algorithm is suitable for parallelization.

The main issue of this paper is the comparison of the three models using this algorithm. For all investigated models, it has been shown in previous publications that they are all capable of handling advanced modeling of the hardening. However, no direct comparison has been possible so far, since the discretization approaches and examples presented differ strongly. Such a comparison is desirable in order to further develop micromechanical models capturing size effects in polycrystals. We therefore choose to investigate all three models in the context of anisotropic polycrystals undergoing large deformations.

The paper is organized as follows. In Section 2, the general underlying mathematical framework is introduced. The basic kinematic and constitutive issues and relations of crystal plasticity relevant for all three models are revised. In a second step, in Section 3 the three different models are presented. Having then all constitutive and governing equations at hand, we outline the dual mixed finite element algorithm in Section 4. Subsequently, a two-dimensional polycrystal is investigated as a numerical example in Section 5. The paper ends with a discussion of the different models.

Basic deformation and stress measures

In this section we reiterate the basic kinematics needed in the following in order to introduce the notation used in this contribution. In large strain plasticity, the main assumption is the classical multiplicative split of the deformation gradient F into an elastic F e and a plastic part F p :

F = F e • F p . ( 2
P := τ • f t = F e • Se • [f p ] t , Se := f e • τ • [f e ] t , (4)
with the notation f x := [F x ] -1 and τ being the Kirchhoff stress. Moreover, the Mandel stress, M e , on the intermediate configuration B, is introduced

M e := [F e ] t • τ • [f e ] t = Ce • Se . (5)
As usual, the model formulation in the context of crystal plasticity is based on the slip-system geometry for each slip-system α: the slip direction s α and the slip plane normal n α which are two orthonormal vectors in the reference configuration B 0 . For the sake of simplicity, it is assumed that the intermediate configuration is isoclinic such that all slip-systems (s α , nα ) are fixed in B and known a priori from the given crystal structure. Therefore, the bars are omitted from now on. Then, following Rice [START_REF] Rice | Inelastic constitutive relations for solids: an internal variable theory and its application to metal plasticity[END_REF], the evolution of the plastic deformation is given by the form

L p = Ḟ p • f p = α γα [s α ⊗ n α ], (6)
with L p denoting the plastic velocity gradient tensor. The resolved shear (Schmid) stress τ α , associated with the slip-system (s α , n α ), is defined as

τ α := s α • M e • n α = M e : [s α ⊗ n α ] . (7)
If γ α is interpreted as the slip-system shear, it can be positive or negative and increase or decrease, in which case γα s α = | γα | sign( γα ) s α = ± | γα | s α . This motivates to consider both s α and -s α as different slip directions. Henceforth, we follow this idea. γ α 0 can then be interpreted as the accumulated slip-system shear and it is always positive and monotonically-increasing, i.e., γα 0 for all α = 1, . . . , n. Of course, from the computational point of view, we have to deal with the expense of doubling the number of slip directions/systems.

Constitutive models

We assume the St. Venant-Kirchhoff law for the elastic behavior is assumed for all three models. Thus the second Piola-Kirchhoff stress Se on the intermediate configuration is given by Se = λ Ī :

Ēe + 2µ Ēe ,

where λ and µ are the Lamé parameters pertinent to small elastic deformations.

Model I

As a first model, we summarize the one developed by Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF], see also [START_REF] Bayley | A three-dimensional dislocation field crystal plasticity approach applied to miniaturized structures[END_REF], which is motivated from physical dislocation densities and evolution equations for these. Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] present a model which quantifies for the densities of geometrically necessary dislocations (GNDs) and statistically stored dislocations (SSDs). Only edge dislocation densities are considered as done in [START_REF] Bayley | A three-dimensional dislocation field crystal plasticity approach applied to miniaturized structures[END_REF] 1 . In addition, they assume that all dislocations of interest lie in the slip planes.

The density of GNDs ρ GND α reads

ρ GND α = ρ GND0 α - 1 b ∇ X γ α • s α , (9)
where b is magnitude of the Burger's vector and ρ GND0 α denotes the initial GND density.

For SSD densities ρ SSD α , Evers et al. apply a generalization of the rate equation proposed by Essmann and Mughrabi [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF]:

ρSSD α = 1 b 1 L α -2y c ρ SSD α | γα |. ( 10
)
The SSD density is initialized by ρ SSD α (t = 0) = ρ SSD0 > 0 and stays positive. The constant y c denotes the critical annihilation length. The average dislocation segment length of statistically stored dislocations L α can be expressed as

L α = K β H αβ |ρ SSD α | + β H αβ |ρ GND α | . ( 11
)
Following the way of Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF], we choose the interaction coefficients H αβ = h 0 ∀α, β with h 0 being the constant immobilization coefficient. The evolution relation for the slip-system shear γ α is given by a power law γα = γ0

|τ eff α | s α m , (12)
relating the slip rates γα to the effective shear stress τ eff α and the slip-system strength s α . The plastic shear rate γ0 and the rate sensitivity parameter m are constant material parameters. The effective shear stress τ eff α is defined as

τ eff α = τ α -τ b : [s α ⊗ n α ] . (13)
where τ b is the global back-stress tensor

τ b = α τ e α [s α ⊗ n α ] + [s α ⊗ n α ] t . ( 14
)
The resulting shear stress τ e α τ e α = develops independently. The circular region around a material point X with radius R accounts for the geometrically necessary dislocations contributing to the shear stress on slip-system α. Moreover, µ is the shear modulus and ν is Poisson's ratio. The slip-system strength s α reads

µbR 2 8 [1 -ν] ∇ X ρ GND α • s α (15)
s α = cµb β a 0 |ρ SSD β | + β a 0 |ρ GND β |, (16)
where c is a constant material parameter. In particular, it is a function of the SSD and GND densities.

Model II

As a second model, we investigate the approach introduced by Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] which is summarized below. In contrast to model I, it is embedded in a thermodynamic framework. Moreover, it is simpler with respect to hardening, but still showed realistic simulation results in [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]. The free energy Ψ0 per unit intermediate volume is introduced as a function of the intermediate Green-Lagrange strain Ēe , the slip-system shears γ α and its gradient ∇ X γ α :

Ψ0 ( Ēe , {γ α }, {∇ X γ α }) = λ 2 tr Ēe 2 + µtr Ēe 2 + 1 2 α H l α γ 2 α + 1 2 α l 2 α ∇ X γ α • [H g α s α ⊗ s α ] • ∇ X γ α . (17)
It is additively decomposed in a hyperelastic contribution following the St. Venant model and a hardening contribution. The hardening contribution to the free energy can be decomposed further into contributions from local and gradient hardening, respectively. The corresponding hardening moduli, H l α and H g α , are chosen as constant and positive semi-definite measures associated with each slip-system α.

This leads to the following hardening stresses κ α which are defined in such a way that they can be derived from the free energy via

κ α := ∂ Ψ0 ∂γ α -Div ∂ Ψ0 ∂(∇ X γ α ) in B 0, grain , α = 1, 2, ..., n slip κ (b) α := N • ∂ Ψ0 ∂(∇ X γ α ) on ∂B 0, grain , α = 1, 2, ..., n slip (18)
see Appendix for further details. Here, N is the outward unit normal to the boundary ∂B 0, grain and the superscript "b" denotes "boundary". The evolution of the plastic slip γ α is assumed to follow a viscoplastic power law Following Gurtin [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF], Lele and Anand [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] propose a thermodynamic-consistent small-deformation strain-gradient theory for isotropic viscoplastic materials similar to the large-deformation model II (which is formulated for anisotropic materials).

γα = 1 t * < τ α -Y α -κ α > C 0 m , (19)
As in the case of model II, the free energy is split additively into an elastic and a plastic contribution and the latter depends on the slip-system shear gradient ∇ X γ α . Moreover, the free energy Ψ0 acts an potential for what they call the energetic microstress ξ en α -in a similar way to the gradient hardening stress κ α of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]. In addition to the approach of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], they introduce a further length scale, the constant dissipative lengthscale l d , corresponding to dissipative effects associated with the slip rate gradient ∇ X γα . Since promising results are presented in [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] for the case of small deformation and isotropy, we extend the flow rule of the large deformation, anisotropic model II in the spirit of Anand, Gurtin and Lele and propose the following expression as the flow rule of model III

γα = 1 t * < τ α -Y α -κ α > C 0 m + l 2 d 0 s α Div d α 0 d q ∇ X γα • s α d α s α . (20)
Hence, the last term in Eq. ( 20) is nonstandard. Interest in such a higher-order dissipative term has arisen only recently, see e.g. the works of [START_REF] Fredriksson | Size-dependent yield strength of thin films[END_REF], [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF], and the very recent [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF]. Here, the scalar 0 s α = s α (t = 0) > 0 represents the initial flow resistance, see also Eq. ( 16). Furthermore, 0 d is a reference rate and d α denotes the scalar non-negative effective flow-rate

d α := | γα | 2 + l 2 d |∇ X γα | 2 . (21)

Model comparison

Here it is appropriate to point out the similarities of the three models. Taking another look at the rate equations of the slip-system shear of model I and II, i.e., Eqs. ( 12) and ( 19), highlights the differences of those two approaches. We rewrite the rate equation of model I by inserting the effective shear stress relations ( 13)-( 15) into the flow rule of model I (Eq. ( 12)) first. Subsequently, the GND density ρ GND α is substituted by relation [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Consequently, we obtain:

γα = γ0   |τ α + µR 2 8[1-ν] β ∇ X [∇ X γ β • s β ] • s β P βα | s α   m = γ0   |τ α + µR 2 8[1-ν] β P βα [s β ⊗ s β ] : [∇ X ⊗ ∇ X ]γ β | s α   m , (22)
where we define the scalar quantity

P βα = [s β ⊗ n β ] + [s β ⊗ n β ] t : [s α ⊗ n α ]
which accounts for cross-gradient hardening. Analogously, inserting the hardening stress relation [START_REF] Levkovitch | On the large-deformation-and continuum-based formulation of models for extended crystal plasticity[END_REF] into the flow rule of model II, i.e., Eq. ( 19), yields Clearly, in both cases, the rate equations contain terms depending on the second spatial derivatives of the slips, i.e., [∇ X ⊗ ∇ X ]γ α . These terms influence the effective resolved stress driving the glide-system dislocation motion. 1Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] use the absolute signs in the numerator whereas Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] employ the Macaulay brackets. Therefore, the rate of the plastic slip can be equal to zero in case of model II as opposed to model I where it is always positive. Thus, in case of model I the slip systems are activated viscously at the stress of s α . On the other hand, model II discretely activates the slip systems at the yield stress Y α .

γα = 1 t * < τ α -Y α -H l α γ α + l 2 α H g α [s α ⊗ s α ] : [∇ X ⊗ ∇ X ]γ α > C 0 m . ( 23
One can easily identify what is called the reference plastic strain rate γ0 in [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] with the inverse of the relaxation time t * in [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF]. Moreover, the material parameter C 0 of model II has the unit of a stress and plays the same role the slip-system strength s α does in the slip-system rate equation of model I. However, C 0 is chosen to be constant and the same for all slip-systems α. On the contrary, the slip-system strength differs for slip-systems and depends on the SSD and GND densities and consequently is an implicit function of the slip-system shear as well, i.e., s α = s α (ρ SSD , ρ GND ; γ). Anticipating the material parameters stated in Table 1 reveals that in the simplest case possible, i.e., no slip and no GND density at initial time t = 0, both are of approximately the same mag-

nitude (1 [MPa] vs. 3.7 [MPa] for C 0 respectively 0 s α = s α (t = 0) = cµb a 0 ρ SSD0 α ).
While model I accounts for local hardening via the fact that the denominator, the slip-system strength s α , depends a.o. on the square root of the SSD densities ρ SSD , Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] include a local hardening contribution in the free energy function Ψ0 and therefore in the hardening stresses κ α -which is part of the numerator. Furthermore, both approaches account for gradient hardening in a similar, but not the same way, via directional gradients of the slip-system shear along the given slip directions via ∇ X γ α • s α .

Lele and Anand's approach [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] is based on the so-called microforce balance representing a viscoplastic yield condition

τ α + Div ξ en α + ξ diss α = f α (γ α , γα ) , (24)
see also Gurtin and Anand [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF]. Here, ξ en α and ξ diss α are referred to as the energetic and the dissipative microstress, respectively, and they expend power over the sliprate gradient ∇ X γα . What is here denoted by a function f is referred to as a scalar microscopic stress in [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] which expends power over the slip-rate γα . In case of model III, the energetic and the dissipative microstress are given by the relations

ξ en α = l 2 α H g α s α [s α • ∇ X γ α ] ξ diss α = l 2 d 0 s α d α 0 d q ∇ X γα • s α d α s α . (25)
Moreover, Clearly, the energetic microstress ξ en α and the dissipative microstress ξ diss α are tangent to the corresponding slip plane α. If plastic loading is obtained in case of model II, i.e., if τ α -Y α -κ α > 0, the microforce balance is obtained with

f α (γ α , γα ) = [t * γα ] 1/m C 0 + Y α + H l α γ α . (26
ξ en α = l 2 α H g α s α [s α • ∇ X γ α ] , ξ diss α = 0, f α (γ α , γα ) = [t * γα ] 1/m C 0 + Y α + H l α γ α . ( 27
)
1

In a similar way, in the case of plastic loading τ eff α > 0, the rate equation ( 12) of model I could be interpreted as a generalized microforce balance with the following choices:

ξ en α = µR 2 8 [1 -ν] β P βα s β [s β • ∇ X γ β ] , ξ diss α = 0, f α (γ α , γα , ∇ X γ α ) = γα γ0 1/m s α . (28)
We choose a more generalized comparison at this point, because the slip strength s α depends on the slip rate γα as well as on the slip gradient ∇ X γ α . Moreover, in the thermodynamic setting of [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF], the energetic microstresses ξ en α are determined in terms of the free energy acting as a potential.

Numerical framework

Governing equations: strong and weak formulation

The mechanical problem is governed by the quasi-static balance of momentum

0 = DivP + b v , (29)
with b v representing the volume forces.

For the slip problem, we apply the following procedure for all three models: In the spirit of Svedberg and Runesson [START_REF] Svedberg | An algorithm for gradient-regularized plasticity coupled to damage based on a dual mixed FE-formulation[END_REF], we introduce an auxiliary scalar field g α for the directional gradient of the slip-system shear γ α along the given slip direction s α ,

g α := ∇ X γ α • s α . ( 30
)
Moreover, following [START_REF] Svedberg | An algorithm for gradient-regularized plasticity coupled to damage based on a dual mixed FE-formulation[END_REF], we approximate the slip-system shear γ α and the scalar g α independently in B 0, grain using a so-called dual mixed space-variational format of the gradient equation g α -∇ X γ α • s α = 0, which is a global equation in space.

1 As Macaulay brackets have been used in the formulation of the flow rule for model II and III, the microforces are determinate because we assume microclamped respectively microfree grain boundary conditions (cf. Eq. ( 31)). In case of an elastically deforming grain with a plastically deforming neighboring grain, the microforces at its grain boundary might be indeterminate -e.g. if continuity conditions are assumed at the grain boundary. In particular, a three-level Newton-Raphson iteration scheme will be applied, see Section 4.2, in which it is possible to choose the FE-approximation of the gradient g α as P q (piecewise polynomial of order q), whereas the FE-approximation of the slip γ α is chosen as P q-1 .

We assume "microclamped" boundary conditions1 , i.e., γ α = 0 on ∂B 0, grain .

Remark: For a quite detailed description of the finite element discretization, the reader is referred to [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], in which the exact FE-approximations, the matrix formulation of the finite element problem, the algorithmic iteration matrices and tool boxes describing the flow of the algorithm are given.

Three-level Newton-Raphson iteration strategy

The basic solution strategy is Newton-Raphson iterations on three levels. Global FE-solutions of the balance of momentum and gradient equations will yield the displacement field u(X) and the gradient field g α (X) at t = t n . There exist different model assumptions on "geometric constraints" on the grain structure, which infer restrictions on the displacement u. It is highly desirable to have a displacement field as unconstrained as possible. Therefore, we only prescribe the displacements at the grain-structure (RVE) boundary. As part of this strategy, we solve for the displacements on the grain boundaries in an outer iteration loop called "grain boundary iteration loop". Then, in the "inner grain iteration loop", updated values of displacements u and gradients g α within each grain are computed. During the "inner grain iteration loop", the displacements on the grain boundaries are kept fixed, given the updated values from the "grain boundary iteration loop". Therefore, the resulting algorithm is suitable for parallelization. Furthermore, a third iteration loop is carried out. The purpose of this "local iteration loop" is to find updated values for the slip γ α (and in case of model I also for the SSD density ρ SSD α ) in each Gauss point, given values of u and g (as provided from the outer loop).

As mentioned above, see Section 4.1, the intergrain conditions are assumed to be micro-clamped, i.e., there is no explicit coupling of slip gradients across the grain boundaries. Consequently, we do not have to solve for the gradients g α during the "grain boundary iteration loop". In the grain iteration loop, the fully discretized system of equations are coupled. Furthermore, in the inner iteration loop of model I, the rate equations ( 12) and ( 10 Although the governing equations ( 29) and (30) are not time-dependent, on the "local iteration loop" level the model contains evolution equations which have to be discretized in time. The state at t = t n-1 is assumed to be known for a given time history of the pertinent loading. Thus, the fully implicit backward Euler scheme is applied to Eqs. ( 6), ( 10), ( 12), [START_REF] Liebe | Theory and numerics of a thermodynamically consistent framework for geometrically linear gradient plasticity[END_REF] and [START_REF] Petch | The cleavage strength of polycrystals[END_REF]. For example, in case of the velocity gradient L p = Ḟ p • f p , i.e., Eq. ( 6), this leads to

I -n-1 F p • n f p - α ∆t γα [s α ⊗ m α ] = 0. (32)
Here, ∆t = [t n -t n-1 ] denotes the current time step and we initialize 0 f p = I and 0 γ α = 0. Furthermore, the slip-system shear rate equations have to discretized in time for all three models. For the sake of comprehensibility, we leave the spatial discretization aside for this moment. Therefore, the semi-discrete formulations of Eqs. ( 12), ( 19) and ( 20) read:

model I 0 = 1 γ0 ∆γ α -∆t |τ eff α | s α m , model II 0 = t * ∆γ α -∆t < τ α -Y α -κ α > C 0 m , (33) model
III 0 = t * ∆γ α -∆t < τ α -Y α -κ α > C 0 m -l 2 d 0 s α Div d α 0 d q ∆g α d α • s α .
In case of the rate equation of model III, we used ∇ X γα • s α = ġα .

5. Numerical example: macro-scale simple shear

Set-up of the model

The grain-structure, occupying the square domain B in 2D with boundary ∂B , as shown in Figure 1, is considered. Dirichlet boundary conditions on ∂B ,D are chosen that correspond to a prescribed macro-scale deformation gradient F m = I + H m (where the subscript "m" indicates a macro-scale quantity). Hence, the boundary conditions for the mechanical problem are

u(X) = H m • X, X ∈ ∂B . ( 34
)
As a consequence, the appropriate spaces of displacements for the sub-scale are defined as

U(H m ) = {u| u = H m • X for X ∈ ∂B }. (35)
The grain-structure consists of 4 grains, see Figure 1 So far, nothing has been said about the number of slip-systems. Henceforth, we limit the numerical study to double slip in order reduce computation time. From the modeling point of view this means α = 4, as explained at the end of Section 2. The directions s α of the slips are randomly distributed, see Figure 1. Values of the material parameters are listed in Table 1.

Plane strain is assumed. In particular, we consider the standard representative example of simple shear on the macro-scale in order to highlight the major characteristics of the different models. Simple shear is defined by F = I + γe 1 ⊗ e 2 , where γ is the macroscopic shear deformation. Computations are carried out in monotonic loading with γmax = 0.15 with the constant loading rate of 0.03 [1/s].

Remark (Material parameters):

The slip rate equation of model I, i.e., Eq. ( 12), is defined in terms of the reference plastic strain rate γ0 and the slip resistance s α . In contrast, the slip rate equations of model II and III, i.e., Eqs. ( 19) and [START_REF] Petch | The cleavage strength of polycrystals[END_REF], are expressed in terms of the relaxation time t * and the constant C 0 . In order to derive comparable results, one has to take care that 1/[t * C m 0 ] ≈ γ0 /s m α (t = 0). In addition, the fraction µR/ [8 [1 -ν]] in Eq. ( 22) (model I) plays a similar role as l 2 α H g α in Eq. ( 23) (model II). As to the choice of the yield stress, a constant value of Y α is chosen in model II, whereas the yield stress of model I is strain-rate dependent and incorporated indirectly via s α .

Discussion of results

First, model I is examined. Figure 2 shows a size effect for fixed material parameters but for increasing side length L of the grain-structure. The plot shows a representation of the plastic strain field in terms of "effective" effective shear strain γ eff def = γ 2 1 + γ 2 2 + γ 2 3 + γ 2 4 . Since the (absolute) size of the boundary layers with large gradient effects should be approximately the same regardless of the actual size of the grains, the boundary layers will appear thinner when the grain size increases (see Figure 2).

The macroscopic stress-strain response in terms of the nominal shear stress component P12 versus the macroscopic shear γ is shown in Figure 3. As expected, stiffer response in the hardening regime can be seen for smaller grain size. This observation is in agreement with the findings of Bayley et al. [START_REF] Bayley | A three-dimensional dislocation field crystal plasticity approach applied to miniaturized structures[END_REF] and Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. The size-dependence of the stress-strain response is due the existence of the slip gradient in the numerator of Eq. ( 22). As already reported by Evers et al. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF], a size dependence for both densities can be found. Figure 4 exemplarily shows the distribution of the GND density ρ GND for two sample lengths. On the left L = 4 µm is depicted and on the right the distribution for a grain-structure with length L = 16 µm is plotted. In agreement with the results of [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF], the GND density ρ GND mainly develops near the grain boundaries. Moreover, we also observe the significant difference between the results for the two lengths. This is due to the fact that the distance between the grain middle and the boundary is much smaller in the case of L = 4 µm -the boundary being the place where the plastic slip is blocked and therefore GNDs develop. On the contrary, the SSD density ρ SSD remains small near the grain boundaries, see Figure 5. Furthermore, the order of magnitude of the GND and the SSD denisties are roughly the same. For the sake of brevity and as we mainly aim at comparing different approaches, we omit figures of this finding for the other examples simulated with model I and instead refer to the discussion in [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF].

One of the differences between model I and II is the use of the absolute signs respectively Macauley brackets. We ran simulations in which we applied Macauley brackets to Eq. ( 12). Due to the fact that load is applied only in one direction, there is no difference for this particular numerical example. However, if loading directions change during the simulation, there should be an impact.

In order to examine the effect of cross-hardening only, simulations are carried out with a constant slip resistance s α = C 0 = 1 [MPa]. In other words, local hardening is omitted. The results are depicted in Figures 6 and7. Due to the gradient contribution, the stress-strain response still depends on the grain-structure side length L. However, the material response is significantly softer than in Figure 3.

In [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF], the rate sensitivity parameter m = 20 is used, which is significantly lower than our choice. Therefore, simulations with a rate sensitivity parameter of m = 5 are run as well in order to investigate its influence, see Figure 8. The higher the rate sensitivy parameter m, the softer responds the material. Now, model II is studied. In all simulations of model II, the rate sensitivity parameter m = 1 is used. Figure 9 presents the distribution of the effective plastic slip γ eff according to Eq. ( 19) with an initial yield stress Y α = 10 3 [MPa] for different grain-structure side lengths. The macroscopic stress-strain response is shown in Figure 10. The results are qualitatively similar to those reported in [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF].

As in the case of model I, the size-dependent stress-strain response is due to the gradient effect. Although both models reach a similar macroscopic stress P12 at 15 % shear, the effective slip-system shear of model I is smaller (cf. Figures 2 and9). For the chosen parameter values less plastic strain develops. Thus, model I accommodates more shear elastically than model II.

In order to investigate the influence of the initial yield stress Y α and the local hardening module H l α (which contributes to the hardening stresses κ α ) in model II, those parameters are neglected one after the other. Figures 11 and12 illustrate the solution with initial yield stress Y α = 0 [MPa] and the local hardening module H l α = 10 4 [MPa] for all slip-systems α = 1, . . . , 4. As expected, the plots of the effective hardening strain γ eff (Figure 11) and the macroscopic stress-strain response (Figure 12) reveal that the initial yield stress Y α = 10 3 [MPa] leads to plastic yielding from the beginning of the deformation. Subsequently, the local hardening module is additionally chosen to be zero, i.e., H l α = 0 [MPa] for all slip-systems α = 1, . . . , 4 and

Y α = 0 [MPa].
It is observed that a material with a vanishing local hardening module H l α still shows a size-dependent stiffness response, see Figures 13 and14. The general stiffness of the material decreases considerably, though. As already mentioned in Section 3.4, model II includes the local hardening contribution in the hardening stress κ α and, thus, in the numerator, whereas model I takes it into account via the slip-system strength s α and, therefore, in the denominator. The comparison of Figures 3 and 7 show the influence of cross-hardening and no cross-hardening. In this special case, the only difference between the two slip laws of model I and II is that cross-hardening is included in the gradient hardening of model I.

In all three cases, the macroscopic stress-strain response is grain-structure size dependent. However, there is no size-dependence in the initial yield limit, see Figures 10,12 and 14. This contradicts the well-known Hall-Petch relation which states that the yield limit is in fact size dependent.

Model II corresponds to what Anand, Lele and Gurtin [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] call the case of gradient energetics respectively energetic hardening. Combining model II with what they refer to as "dissipative strengthening" leads to model III. The results for model III are pointed out in Figures 151617. As stated by Lele and Anand [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF], the divergence term in the flow rule [START_REF] Petch | The cleavage strength of polycrystals[END_REF] leads to a size-dependent initial yield limit, see Figure 17. However, for the parameters chosen in this example of large deformations, Figure 17 shows that the effect unfortunately is not strong as in the example of [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF]. However, in [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] in all the simulations in which the nonstandard higher-order dissipative term was considered, the energetic hardening and standard isotropic hardening were set at zero. In our example, a higher value for the dissipative lengthscale l d would have been desirable in order to receive a considerable side-length dependent yield stress for the case that local hardening is admitted (Figure 16).

As stated in [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF], the model "possesses a mathematically attractive structure, our experience with numerical experiments which use these constitutive equations is that they are too tightly coupled". A completely different numerical approach was applied in this contribution 1, nevertheless and unfortunately, for larger l d the code did not converge.

Conclusion

Three different frameworks for including gradient hardening into crystal plasticity models were studied. Model I is physically based and has already proven to be capable of handling advanced modeling of hardening. As compared to this model, the purpose of model II is to be simple, e.g. with respect to the hardening, but still able to show realistic simulation results. Moreover, it is formulated in a thermodynamic setting -as is model III. In both cases, the inelastic part of the free energy was assumed to include contributions from the gradient of hardening along each slip direction. One can clearly see the different behavior depending on the model assumptions. Model I is capable of mapping more effects than model II, such as cross-hardening. However, it is also much more complex which leads to higher computation time. Model III is a first step towards mapping the Hall-Petch relation, i.e. the squareroot dependence of the yield strength on the grain-size. Unfortunately, numerical convergence problems so far limit the advantage of the mathematically attractive structure. By investigating the meaning and the influence of some of the parameters of the strain gradient crystal plasticity models discussed above, this contribution presents a step towards the enhancement of such theories. Moreover, an efficient discretization scheme suitable for parallelization has been suggested. The arising coupled boundary value problems were discretized with finite elements in space and the backward Euler in time. A dual mixed approach with displacements and gradient hardening variables as degrees of freedom is used. Such an approach is attractive since any crystal plasticity code (for the local problem) can be used by simply modifying the current yield stress due to the given gradient hardening variables in the nodes. The dual mixed finite element algorithm renders a coupled two-field problem which is solved monolithically. In particular, an algorithm suitable for parallelization was presented, where each grain is treated as a subproblem. A fast solution algorithm opens up for a variety of applications, since the resulting finite element problems tend to be rather large. Possible applications are concurrent multiscale modeling, or three dimensional models, for example.

In conclusion, the presented models and the numerical strategy are all suitable for handling, e.g., more advanced modeling of the hardening. Nevertheless, work remains to be done with respect to both, the modeling (e.g. an accurate prediction of the Hall-Petch effect) and the numerical implementation (regarding e.g. the nonstandard higher order dissipative term).

Appendix: Thermodynamic considerations

As mentioned above, there exists a thermodynamic setting for model II and III. First, the ideas of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF] (model II) are presented and a short derivation for the definition of the hardening stresses κ α , i.e., Eq. ( 18), is given. Second, the considerations of Lele and Anand [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF], see also Gurtin and Anand [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF], are summarized.

In rational thermodynamics the second law imposes restrictions on constitutive equations. Whereas there exists a unified opinion about the balance laws, this is not the case for the second law of thermodynamics. Exploitations of the entropy inequality are based on distinctive mathematical procedures with the help of which conditions for a given class of constitutive material behaviors are derived in order to be compatible with the entropy inequality. Probably the best-known principle is that of Clausius and Duhem in the framework established by Coleman and Noll [START_REF] Coleman | The thermodynamics of elastic materials with heat conduction and viscosity[END_REF].

The thermodynamic formulation of field and constitutive relations are ultimately formulated in the context of the total energy balance. Isothermal and quasi-static conditions are considered. In this case, the total energy balance reads

D t B0 Ψ 0 dV + B0 D 0 dV + B0 W int dV = 0 (36)
Here, Ψ 0 is the free energy per unit undeformed volume, D 0 the dissipation power and W int stems from the rate of mechanical work supply. Although the balance equation is stated with respect to integrals over the entire material domain B 0 respectively its entire boundary, it also has to hold for arbitrary bounded and connected parts of the body.

Considerations of Ekh et al. [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF][START_REF] Svedberg | An algorithm for gradient-regularized plasticity coupled to damage based on a dual mixed FE-formulation[END_REF] The free energy Ψ 0 (per unit undeformed volume), in the material configuration B 0 , is proposed as a function of the deformation gradient F and a set of internal variables, which are here taken as F p and the set of slip-system shears {γ α } addition, the free energy Ψ 0 is assumed to depend on the set

{∇ X γ α } nslip α=1
, where ∇ X is the gradient in the reference configuration B 0 . Thus, gradient effects (under isothermal conditions) are included and the scalar variables γ α represent isotropic hardening on the n slip slip-systems. From the principle of (spatial) objectivity, we may conclude that a possible parametrization of the free energy Ψ 0 is given by

Ψ 0 (F , F p , {γ α }, {∇ X γ α }) = Ψ0 ( Ēe , {γ α }, {∇ X γ α }). (37)
Moreover, the the rate of mechanical work supply reads W int def = B0 P : Ḟ dV . For thermodynamical consistency, it is indispensable that the dissipation is strictly nonnegative, corresponding to the second law of thermodynamics. Thus, the fulfilment of the (reduced) dissipation inequality

B0 D 0 dV = B0 M e : L p - α ∂ Ψ0 ∂γ α γα - α ∂ Ψ0 ∂(∇ X γ α ) • ∇ X γα dV ≥ 0 (38)
is an absolutely essential modeling-feature. Upon integrating the last term by parts and using Gauss' theorem, we rewrite Eq. ( 38) as

B0 M e : L p - α κ α γα dV - ∂B0 α κ (b) α γα dA ≥ 0. ( 39
)
This leads to the definition of the auxiliary dissipative hardening stresses as stated in Eq. [START_REF] Levkovitch | On the large-deformation-and continuum-based formulation of models for extended crystal plasticity[END_REF]. It is noted that the hardening stresses κ α may be composed of local as well as gradient contributions, whereas the gradient tractions κ

α on the boundary ∂B 0, grain represent the gradient effect only.

Note that a more restrictive formulation than Eq. ( 38) is obtained, if it is required that the local dissipation contributions are non-negative pointwise, i.e., if D 0 ≥ 0 for all X ∈ B 0 and D (b) 0 ≥ 0 for all X ∈ ∂B 0 .

Considerations of Anand, Lele and Gurtin [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] Anand, Lele and Gurtin [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] apply a special kind of virtual work principle in conjunction with arguments of classical thermodynamics. Key idea is the assumption that each evolution of the body is associated with macroscopic and microscopic forces. The latter depends on a scalar microscopic force f α whose working accompanies slip on slip-system α and a vector microscopic stress ξ α whose working accompanies slip gradients on α. This leads to a different choice of an internal energy

W int def = B0 S e : Ėe + α [f α γα + ξ α • ∇ X γα ] dV = B0 P : Ḟ + α [Div (ξ α ) γα + ξ α • ∇ X γα ] dV, (40)
with S e expanding power over the Green-Lagrange strain Ėe . It is worth noting that, here, the work conjugacy of the scalar microscopic force rate γα , respectively the microscopic stress ξ α and the slip rate gradient ∇ X γα , is assumed first. Then the existence of a free energy function

Ψ 0 (F , F p , {γ α }, {∇ X γ α }) = Ψ0 ( Ēe , {∇ X γ α }). (41)
is established. This is in contrast to the procedure described in Section 6, where first a free energy function depending on the slip-system shear γ α and its gradient ∇ X γ α is assigned which then, as a consequence, leads to the definition of microscopic hardening stresses κ α , i.e., Eq. ( 18). Moreover, the free energy function of Anand, Lele and Gurtin [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF] does not depend on the slip-system shear γ α , only on its gradient. Therefore, the counterpart of dissipation inequality (38) reads

B0 M e : L p + α [Div (ξ α ) γα + ξ α • ∇ X γα ] dV ≥ 0. ( 42
)
The strategy of Anand, Lele and Gurtin leads to the microforce balance, i.e., Eq. ( 24), where an additive decomposition of the microstress ξ α into an energetic and a dissipative component is admitted

ξ α = ξ en α + ξ diss α . (43)
The energetic microscopic stress for slip-system α, i.e., ξ en α , is given by

ξ en α def = ∂ Ψ0 ∂(∇ X γ α ) . (44)
Then, the dissipation takes the form

α f α γα + ξ diss α • ∇ X γα dV ≥ 0. ( 45
)
Inequality (45) imposes restrictions on the constitutive equations for f α and ξ diss α for which a constitutive dependence on the slip rate γα and the slip rate gradient ∇ X γα is suggested, cf. [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF][START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF].

In addition, since the microscopic stress ξ diss α characterizes dissipative microscopic forces associated with the evolution of dislocations on the αth slip plane, and because the motion of such dislocations is tangent to this plane, it sometimes is required that dissipative microscopic stress ξ diss α is tangential, cf. e.g. [START_REF] Gurtin | A theory of strain-gradient plasticity for isotropic, plastically irrotational materials. Part I: Small deformations[END_REF].

= γ 2 1 + γ 2 2 + γ 2 3 + γ 2 4 .
From left to right the solutions for grain-structure side lengths L = 4 µm, L = 8 µm, and L = 16 µm are depicted. The evolution of the slip-system shear is given in Eq. [START_REF] Fredriksson | Size-dependent yield strength of thin films[END_REF]. Computations are carried out in monotonic loading with a maximum macroscopic shear deformation γmax = 0.15. The displacement field inside the grainstructure is unconstrained. The slip resistance sα is computed according to Eq. ( 16), the rate sensitivity parameter reads m = 1. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Macroscopic stress-strain response ( P12 vs. γ) showing the size dependence on the amount of hardening. Under the assumption of increasing slip resistance sα according to Eq. ( 16) and m = 1. [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Slip-system 1 GND densities are plotted for two samples lengths after applying 15% shear. One the left the distribution of the GND density ρ GND for a grain structure with L = 4 µm and on the right with L = 16 µm can be seen. A size-dependence is found for the density distribution. Under the assumption of increasing slip resistance sα according to Eq. ( 16) and m = 1. The sign of the ρ GND is ignored in this plot. The results are in accordance to [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. 16) and m = 1. The results are in accordance to [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF].

F

= γ 2 1 + γ 2 2 + γ 2 3 + γ 2 4 .
From left to right the accumulated plastic strain field for grain-structure side lengths L = 4 µm, L = 8 µm, and L = 16 µm are depicted. The evolution of the slip-system shear is given in Eq. [START_REF] Fredriksson | Size-dependent yield strength of thin films[END_REF]. Computations are carried out in monotonic loading with a maximum macroscopic shear deformation γmax = 0.15. The displacement field inside the grain-structure is unconstrained. The slip resistance sα is assumed to be constant with sα = C 0 = 1 [MPa] and m = 1. α . Again, it can be stated that smaller the grain-structure size the stiffer is the material's response. However, it can clearly be seen that H l α = 0 [MPa] leads to a significant decrease in the stiffness for all lengths L. m = 1. 1. Material parameters. The data is taken from [START_REF] Ekh | Gradient crystal plasticity as part of the computational modeling of polycrystals[END_REF], [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF], [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] and [START_REF] Lele | A small-deformation strain-gradient theory for isotropic viscoplastic materials[END_REF]. Some of the values are adjusted in such a way that comparison between the models is possible, see also the remark on this aspect. During simulations some of the parameters are varied in order to study their influence in the material behavior.

  the free energy.

  ) are coupled. In both cases the coupled system is solved monolithically. http://mc.manuscriptcentral.com/pm-

  , and is discretized by 1990 elements (1384 nodes).

14 Page

 14 versus 12 and 14 underlines that, nevertheless, both models have a similar local hardening response. Figures 3 and 12 include local hardening, whereas Figures 7 and 14 neglect this contribution -leading to a significant decrease of the material's stiffness. Moreover, the cases studied in Figures 6 and 7 respectively 13 and

  f α and the slip http://mc.manuscriptcentral.com/pm-

Figure 1 .Figure 2 .

 12 Figure 1. Left: Discretization and grain geometry used during simulations. The square grain-structure consists of 4 grains and the thick black lines represent the grain boundaries. The grains are discretized with 726 (upper left grain), 620 (upper right), 630 (lower left), and 614 (lower right) elements, respectively. The mesh is refined near the grain boundaries. The side length L of the grain-structure is varied during simulations in order to capture the size effects. Right: Schematic sketch of slip plane directions sα. Double slip is assumed and slip plane directions are randomly distributed.

Figure 3 .

 3 Figure 3. Model I, Evers et al.[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Macroscopic stress-strain response ( P12 vs. γ) showing the size dependence on the amount of hardening. Under the assumption of increasing slip resistance sα according to Eq. (16) and m = 1.

Figure 4 .

 4 Figure 4. Model I, Evers et al.[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Slip-system 1 GND densities are plotted for two samples lengths after applying 15% shear. One the left the distribution of the GND density ρ GND for a grain structure with L = 4 µm and on the right with L = 16 µm can be seen. A size-dependence is found for the density distribution. Under the assumption of increasing slip resistance sα according to Eq. (16) and m = 1. The sign of the ρ GND is ignored in this plot. The results are in accordance to[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF].

Figure 5 .

 5 Figure 5. Model I, Evers et al.[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF]. Slip-system 1 SSD densities are illustrated for two samples lengths after applying 15% shear. One the left the distribution of the SSD density ρ SSD for a grain structure with L = 4 µm and on the right with L = 16 µm can be seen. ρ SSD decreases towards the boundary. Under the assumption of increasing slip resistance sα according to Eq. (16) and m = 1. The results are in accordance to[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF].

Figure 7 .Figure 8 .Figure 9 .+ γ 2 2 + γ 2 3 + γ 2 4 .Figure 10 .Figure 11 .+ γ 2 2 + γ 2 3 + γ 2 4 .Figure 12 .Figure 13 .+ γ 2 2 + γ 2 3 + γ 2 4 .Figure 14 .

 7892241011224121322414 Figure 7. Model I, Evers et al. [9]. Macroscopic stress-strain response ( P12 vs. γ) showing the size dependence on the amount of hardening. m = 1. Constant slip resistance sα = C 0 = 1 [MPa] presumed, i.e., no local contribution to the hardening exists. The grain-structure still shows a size-dependent response.

Figure 15 . 2 1 + γ 2 2 + γ 2 3 + γ 2 4 .LFigure 16 .LFigure 17 .

 1522241617 Figure 15. Model III. Effective hardening strain γ eff = γ 2 1 + γ 2 2 + γ 2 3 + γ 2 4. Left the accumulated plastic strain field according to model III with l d = 0 µm for a grain-structure side length 16 µm is depicted. On the right, the accumulated plastic strain field according with l d = 0.06 µm is shown. Computations are carried out in monotonic loading with a maximum macroscopic shear deformation γmax = 0.15. The evolution of the slip-system shear γα is given in Eq. (20) -with an initial yield stress Yα = 10 3 [MPa] and a local hardening module H l α = 10 4 [MPa] for all slip-systems α = 1, . . . , 4. The influence of the divergence term is very small, compare Figures16 and 17. The displacement field inside the grain-structure is unconstrained and m = 1..

  Next, we recall the stress measures of interest. The 1st Piola-Kirchhoff stress, P , and the 2 nd Piola-Kirchhoff stress, Se on the intermediate configuration B, are defined as:

	where the right Cauchy-Green tensor in the intermediate configuration B is denoted by Ce and Ēe is the Green-Lagrange strain.
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	)

  with < x >:= 1/2[x+|x|] denoting the Macaulay brackets. The material parameter C 0 is constant and the same for all slip-systems and Y α and t * denote the initial yield stress and the relaxation time, respectively.

	3.3. Model III
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

The model introduced in[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] is formulated for edge and screw dislocation densities. However, the numerical example presented in[START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF] is restricted to self-hardening and double slip where only edge dislocations are considered.

Different stress terms and their relations are compared in[START_REF] Forest | Some links between Cosserat, strain gradient crystal plasticity and the statistical theory of dislocations[END_REF] -with the main focus on linking the phenomenological Cosserat theory, the statistical theory of dislocations and strain gradient crystal plasticity.

Both, microclamped and microfree boundary conditions, rule out the possibility of interaction of slip systems across a grain boundary and is therefore quite restrictive. Nevertheless, it is quite popular and has been used by e.g.[7, 13, 14, 17, 19,

[START_REF] Svendsen | On thermodynamic and variational-based formulations of models for inelastic continua with internal lengthscales[END_REF]. There also exist contributions introducing boundary conditions that may account for grain interactions, see e.g.[START_REF] Fredriksson | Size-dependent yield strength of thin films[END_REF]. In this contribution we limit our analysis to the particular case described above, as the influence of different grain boundary conditions is part of future research and a topic of its own.

Lele and Anand implemented the theory by writing a nine-node quadratic user-element subroutine for the commercial finite element package ABAQUS/Standard.

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgement

This work was supported by the Swedish Research Council, which is gratefully acknowledged. This research was carried out while Swantje Bargmann was at Chalmers University, Sweden, as a post-doctoral fellow. The authors are very grateful to the anonymous referees for their helpful comments and suggestions leading to the improvement of this contribution.


Andrey Shishkov

Laurent Véron

Igor Vladimirovich Skrypnik

Propagation of Singularities of Nonlinear Heat Flow in Fissured Media

Introduction

In this paper we investigate the propagation of singularities in a nonlinear parabolic equation with strong absorption when the absorption potential is strongly degenerate following some curve in the (x, t) space. As a very simplified model, we assume that the heat conduction is constant but the absorption of the media depends stronly of the characteristic of the media. More precisely we suppose that the temperature u is governed by the following equation ∂ t u -∆u + h(x, t)u p = 0 in Q T := R N × (0, T ) (

where p > 1 and h ∈ C(Q T ). We suppose that h(x, t) > 0 except when (x, t) belongs to some space-time curve Γ given by Γ = {γ(τ ) := (x(τ ), t(τ

)) : τ ∈ [0, T ]} , (1.2)
where γ ∈ C 0,1 ([0, 1]) with γ(0) = 0 and t(τ ) > 0 for τ ∈ (0, T ]. If there holds

T 0 B R h(x, t)E q (x, t)dxdt < ∞ (1.3)
we first show that, for any k ≥ 0, there exists a unique solution u k of (1.1) such that u k (., 0) = kδ 0 . Furthermore the mapping k → u k is increasing. Because h is continuous and positive outside Γ, we shall show that the set of solutions {u k } remains locally bounded in Q T \ Γ. Therefore u k ↑ u ∞ and u ∞ is a solution of (1.1) in Q T \ Γ. Then either the singularity of the solution issued from the point (0, 0) can propagate along Γ (at least partially), or it remains localized at (0, 0). More precisely we show Theorem A. Assume γ is C 1 and t ′ (τ ) > 0 for any τ ∈ (0, T ]. The following dichotomy of phenomena occurs (i) either u ∞ (x, t) < ∞ for all (x, t) ∈ Q T , (ii) or there exists τ 0 ∈ (0, T ] such that lim sup (x,t)→(y,s)

u ∞ (x, t) = ∞ ∀(y, s) ∈ Γ, 0 ≤ s ≤ τ 0 .

(1.4)

We first prove that the singularity does not propagate from a point where t(τ ) is decreasing. We prove Theorem B. Assume that γ is C 1 and t ′ (τ ) < 0 for all τ ∈ (τ 0 , T ] for some τ 0 > 0. Then u ∞ remains locally bounded in Q T \ γ((τ 0 , T ]).

Due to this fact we shall assume first that the t variable is increasing along Γ, in such a case we can assume that τ is a function of t and, up to a change of parameter, that τ = t and Γ = {(x(t), t)

: t ∈ [0, T ]} , (1.5)
where γ ∈ C 0,1 (0, T ] satisfies γ(0) = 0. In order the singularity to propagate along Γ, h must be very flat near this curve Γ. If we define the parabolic distance between (x, t) and Γ by

d P (x, t); Γ) = inf{d P (x, t); (y, s)) : (y, s) ∈ Γ, s ≤ t} where d P (x, t); (y, s) = |x -y| + √ t -s, for t ≥ s and write h(x, t) = e -ℓ(d P (x,t);Γ)) ,
where ℓ is a positive nonincreasing function. In the next result we shows that propagation of singularity along Γ occurs

Theorem C. Assume t → x(t) belongs to W 2,∞ loc [0, ∞) and lim inf t→0 t 2 ℓ(t) > 0. (1.6) Then lim (x,t)→(y,s)) u ∞ (x, t) = ∞ for all (y, s) ∈ Γ.
The last section is devoted to the case where the curve of degeneracy Γ is a straight line contained in the initial plan. We set

d ∞ ((x, t), Γ) = max{ √ t, x ′ } if dist (x, Γ) = x ′ and we write h(x, t) = e -ℓ(d∞(x,t);Γ))
Up to a rotation, we can suppose that Γ is the x 1 axis and denote

x = (x 1 , x ′ ) the component in R × R N -1 .
Then we prove the following

Theorem D. Assume lim inf t→0 t 2 ℓ(t) > 0. (1.7) Then lim (x,t)→(y 1 ,0)) u ∞ (x, t) = ∞ for all y 1 ∈ R.
Even if this model is a very simplified version of the heat propagation in a fissured absorbing media, it gives interesting insight of the propagation phenomenon which can occur. It is also a starting point for studying other type of propagation of singularities in nonlinear diffusion equations. In a forthcoming article we shall consider the case where the degeneracy line is a surface in R N × [0, ∞) with only one contact point with t = 0 at (0, 0).

Preliminaries and basic estimates

Through out this section we assume that

h ∈ C(R N × [0, ∞)) is nonnegative. Since E(x, t)
is a supersolution for (1.1), the following result holds [10, Theorem 6.12] Proposition 2.1 Assume q > 1 and (1.3) holds. Then for any k > 0 there exists a unique

u = u k ∈ C(Q T \ {(0, 0)}) ∩ L 1 (Q T ), such that h|u| q ∈ L 1 (Q T ), satisfying Q T (-u∂ t ζ -∆ζ + h(x, t)|u| q-1 uζ)dx dt = kζ(0, 0) (2.1) for all ζ ∈ C 2,1 (Q T ) which vanishes at t = T . More generaly, if µ ∈ M(R N ) and ν ∈ M(R N × (0, T )) are two positive bounded measures, the solution v = v µ,ν of ∂ t v -∆v = ν in R N × [0, T ) v(., 0) = µ in R N , (2.2)
is expressed by

v µ,ν (x, t) = 1 (4πt) N/2 R N e -4|x-y| 2 /4t dµ(y) + t 0 1 (4π(t -s)) N/2 R N
e -4|x-y| 2 /4(t-s) dν(y)ds.

(2.3) Actually, by direct adaptation to the parabolic case of [START_REF] Véron | Singularities of Second Order Quasilinear Equations[END_REF]Theorem 4.2], combined with [START_REF] Véron | Singularities of Second Order Quasilinear Equations[END_REF]Theorem 6.12], one can proves that, for any bounded Radon measures µ on R N ×(0, ∞) and ν on R N which satisfy

Q T v q µs,νs h(x, t)dxdt < ∞ (2.4)
where µ s and ν s are the singular parts (with respect to respective Lebesgue measures) µ and ν respectively, there exists a unique weak solution u to problem

∂ t u -∆u + h(t, x)|u| q-1 u = ν in R N × [0, T ) u(., 0) = µ in R N , (2.5)
The next result is an adaptation of Brezis-Friedman a priori estimate [START_REF] Brezis | Nonlinear parabolic equations involving measures as initial coinditions[END_REF].

Proposition 2.2 Let Q = Q r,a t 0 ,t 1 := B r (a) × (t 0 , t 1 ) for some a ∈ R N , t 1 > t 0 ≥ 0 and r > 0 and assume β = min{h(x, t) : (x, t) ∈ Q} > 0. Then any solution of (1.1)

in Q satisfies |u(x, t)| ≤ C β 1/(q-1) 1 t -t 0 + 1 (r -|x -a|) 2 1/(q-1) ∀(x, t) ∈ Q, (2.6)
for some C = C(N, q) > 0.

Proof. The maximal solution of the parabolic equation y ′ + βy q = 0 on 0, ∞) is expressed by

y M (t) = 1 β(q -1)t 1/(q-1)
.

(2.7)

By Keller-Osserman estimate, the maximal solution v M of -∆v + β|u| q-1 u = 0 in B r , satisfies

v M (x) ≤ C 1 β(r -|x|) 2 1/(q-1) (2.8) in B r . Since y M (t -t 0 ) + v M (x -a
) is a supersolution of (1.1) in Q which blows up on the parabolic boundary, an easy approxition argument (just replacing r by {r n } ↑ r and t 0 by {t n } ↓ t 0 ) leads us to (2.6).

Proposition 2.3 Let 0 < r 0 < r 1 , 0 < t 0 < t 1 and Θ := Θ r 0 ,r 1 ,a t 0 ,t 1 := Q r 1 ,a t 0 ,t 1 \ Q r 0 ,
a t 0 ,t 1 for some a ∈ R N and assume β = min{h(x, t) : (x, t) ∈ Θ} > 0. Then any solution of (1.1 )

in Q r 1 ,a t 0 ,t 1 such that |u(x, t 0 )| ≤ µ for x ∈ B r 1 (a) satisfies |u(x, t 1 )| ≤ µ + C (β(r 1 -r 0 ) 2 ) 1/(q-1)
∀x ∈ B r 0 (a), (2.9)

for some C = C(N, q) > 0.

Proof. Let b ∈ R N such that |b -a| = r 0 +r 1 2 and v M the maximal solution of -∆v

+ β|u| q-1 u = 0 in B r 1 -r 0 2 (b). Then v M (x) ≤ C β r 1 -r 0 2 -|x -b| 2/(q-1) ∀x ∈ B r 1 -r 0 2 (b). Since v M + µ is a super solution of (1.1 ) in Q r 1 -r 0 2
,b t 0 ,t 1 which dominates u at t = t 0 and for

|x -b| = r 1 -r 0 2 , it follows that u(x, t) ≤ v M (x, t) + µ in Q r 1 -r 0 2 ,b t 0 ,t 1 . In particular, if x = b, we get u(b, t) ≤ µ + C β r 1 -r 0 2 2/(q-1)
This estimate is valid for any b ∈ R N with |b -a| = r 0 +r 1 2 . Since u is a subsolution of the heat equation in Q r 0 ,a t 0 ,t 1 , (2.9) follows by the maximum principle.

The previous estimates are based upon constructions of supersolutions in cylinders. In the next result we construct estimates in tubular neighborhood of Γ. If t(τ ) is increasing, we can take τ = t and

Γ = γ([0, T ]) := {(x(t), t) : t ∈ [0, T ]} (2.10)
Proposition 2.4 Assume Γ is C 1 and parametrized by t as in (2.10), and for ǫ > 0 denote

T Γ T 0 s := {(x, t) : 0 ≤ t ≤ T, |x -x(t)| < s}. For 1 ≥ r 1 > r 0 > 0 we set η = min h(x, t) : (x, t) ∈ T Γ T 0 r 1 \ T Γ T 0 r 0 . Then there exists a constant C depending on c := max{|x ′ (t)| : 0 ≤ t ≤ T } such that if u(x, 0) ≤ m + C 2 (r 1 -r 0 ) 2 q-1 ∀x ∈ B r 1 +r 0 2 (x(0)), (2.11) then u(x, t) ≤ m q-1 1 + η(q -1)m q-1 t 1 q-1 + C 2 (r 1 -r 0 ) 2 q-1 ∀t ∈ [0, T ] and x ∈ B r 1 +r 0 2 (x(t)) (2.12) Proof. Consider the change of space variable x = y + x(t) and u(x, t) = v(y, t). Then v satisfies ∂ t v -∆v + x ′ (t), ∇v + h(y + x(t), t)v q = 0 in Q r 0 ,r 1 0,t(T ) .
(2.13)

Thus ∂ t v -∆v -c|∇v| + ηv q ≤ 0 in Q r 0 ,r 1 0,t(T ) .
In the ball B r 1 -r 0 2 (z) where |z| = r 1 +r 0 2

It is standard easy to construct a radial function

ψ satisfying -∆ψ -c|∇ψ| + ηψ q ≥ 0 in B r 1 -r 0 2 (z) (2.14)
under the form

ψ(y) = C ρ 2 q-1 (ρ 2 -|y -z| 2 ) 2 q-1 where ρ = r 1 -r 0 2 and C = C(N, q, c, η). Therefore ψ(z) = C 1 ρ 2 q-1 = C 2 (r 1 -r 0 ) 2 q-1 The solution φ m = φ of φ ′ + ηφ q = 0 on (0, T ) φ(0) = m > 0 (2.15)
is expressed by

φ m (t) = m q-1 1 + η(q -1)m q-1 t 1 q-1 .
Consequently, the function φ m (t) + ψ(y) is a supersolution for (2.13) in Q r 0 ,r 1 0,t(T ) which domintates v at t = 0 and on the lateral boundary. Therefore it is larger than v and in particular v(z, t) ≤ φ m (t) + ψ(z).

(2.16)

Consequently

u(x, t) ≤ φ m (t) + ψ(z) ∀t ∈ (0, T ) and |x -x(t)| = r 1 + r 0 2 .
(2.17)

We derive (2.12) by the maximum principle.

Geometric obstruction to propagation

We assume that h vanishes on a continuous curve

Γ ⊂ R N × [0, T ) defined by parametri- sation Γ = {γ(t) := (x(τ ), t(τ ) : τ ∈ [0, T ]} (3.1)
issued from (0, 0) (i.e. x(0), t(0) = (0, 0) with t(τ ) > 0 if τ ∈ (0, T ). and τ → γ(τ ) is Lipschitz with no self intersection, which means that τ → γ(τ ) is one to one.

Proof of Theorem A. Since t ′ is continuous and increasing, we apply Proposition 2.4 with 0 replaced by τ 0 , we set t 0 = t(τ 0 ) and write Γ under the form (2.10). If we assume that lim sup

(x,t)→(x(t 0 ),t 0 ) u(x, t) < ∞,
there exists r 1 > 0 and µ > 0 such that u(x, t 0 ) ≤ µ if x ∈ B r 1 (x(t 0 )). Then, for 0 < r 0 < r 1 there exists m > 0 such that (2.10) is verified. Thus (2.13) holds. This implies that the blow-up set of u along Γ is empty if t ≥ t 0 .

Proof of Theorem B. The proof is based upon the same ideas than in Theorem A above except that we only study the part of Γ between τ 0 and T , where t ′ (τ ) < 0. Let t 0 = t(τ 0 ) and t * = t(T ). We parametrized Γ by t between t * and t 0 , thus

γ([τ 0 , T ]) = {(t, x(t)) : t * ≤ t ≤ t 0 }.
For s > 0 and t * < t ′ ≤ t 0 , we set T

Γ t ′ t * s := {(x, t) : t * ≤ t ≤ t 0 , |x-x(t)| < s}, then, for t * < t 1 < t 0 , there exists r 1 > 0 and τ 1 ∈ (τ 0 , T ) such that t 1 = t(τ 1 ) and Γ ∩ T Γ t 1 t * r 1 = γ([τ 1 , T ]). If t * = 0, then u(x, t * ) = 0 for |x -x(t * )| ≤ r 1 . If t * > 0, then h(x, t) > 0 in the cylinder Q r 1 ,x(t * ) 0,t * , thus there exists β > 0 such that inf h(x, t) : (x, t) ∈ Q r 1 ,x(t * ) 0, t * 2 = β. Up to replacing r 1 by some r ′ 1 > r 1 such that Γ ∩ Q r 1 ,x(t * ) 0, t * all x ∈ B r 1 (x(t * )) from Proposition 2.2. In both cases u ∞ (x, t * ) is bounded in B r 1 (x(t * )).
Replacing (0, T ) by (t * , t 1 )) in Proposition 2.4, it follows that u ∞ remains bounded in

T Γ t 1 t * r 1 +r 0 2
for some 0 < r 0 < r 1 . Since t 1 < t 0 is arbitrary and u ∞ is locally bounded in Q T \ Γ, the proof follows.

In the next case the monotonicity of τ → t(τ ) is replaced by a box-assumption.

Proposition 3.1 Assume that γ is continuous and there exists a ∈ R N , r 0 > 0 and

τ 0 ∈ (0, T ) such that t(T ) ≤ t(τ 0 ) and γ([τ 0 , T ]) ⊂ B r 0 (a) × [t(T ), t(τ 0 )]. Then u ∞ is bounded in B r 0 (a) × [t(T ), t(τ 0 )].
Proof. There exist r ′ 0 < r 0 and β > 0 such that γ([τ 0 , T ]) ⊂ B r ′ 0 (a) × [t(T ), t(τ 0 )] and min h(x, t) : (x, t) ∈ Θ r ′ 0 ,r 0 t(T ),t(τ 0 ) = β. Therefore the conclusion follows from Proposition 2.3.

In the next case we show that non-propagation of singularities may occur even if τ → t(τ ) is increasing after some τ 0 , provided there is a local maximum in (0, τ 0 ). We put

Γ τ 0 = γ([0, τ 0 ]) Theorem 3.2 Let γ be C 1 and γ ′ (τ ) = 0 on [0, T ]. Assume there exist τ 0 > 0, a ∈ R N and r > 0 such that t(τ ) ≤ t(τ 0 ) on [0, τ 0 ], τ → t(τ ) is decreasing on (τ 0 , τ 0 + δ) for some δ > 0, γ((τ 0 , T ]) ⊂ Q 0,∞ r,a , Γ τ 0 ∩ Q 0,∞ r,a = {γ(τ 0 )} and for any τ ∈ (τ 0 , T ], |x(τ ) -a| < r. Then u ∞ is locally bounded in Q \ Γ τ 0 , where Γ τ 0 := γ([0, τ 0 ]).
Proof. Since t(τ ) is decreasing on (τ 0 , τ 0 + δ), for any τ 1 ∈ (τ 0 , τ 0 + δ), the set of τ ∈ (τ 1 , T ] such that t(σ) < t(τ 1 ) for all σ ∈ (τ 1 , τ ) is not empty. Its upper bound τ * 1 is less or equal to T and t(τ * 1 ) = min{t(τ 1 ), t(T )}. If lim τ 1 →τ 0 t(τ * 1 ) = t(T ), then t(τ 0 ) = t(T ), the box-assumption holds and the conclusion follows from Proposition 3.1. If t(T ) > t(τ 0 ), then t(τ 1 ) = t(τ * 1 ) < t(T ) for any τ 1 ∈ (τ 0 , τ 0 + δ). Since lim τ 1 →τ 0 t(τ 1 ) = t(τ 0 ) and γ is continuous with γ((τ 0 , T ]) ⊂ Q 0,∞ r,a , there some fixed constants λ > 0 and ρ > 0 such that x(τ ) ∈ B ρ (x(τ 1 )), for any τ ∈ [τ 1 , T ] verifying t(τ ) ≤ t(τ 1 ) + λ. This means that the part of Γ starting from γ(τ 1 ) for which t(τ ) belongs to (t(τ 1 ), t(τ 1 ) + λ) remains in B ρ (x(τ 1 )). Moreover we can assume that t(τ 1 )+λ > t(τ 0 ). By restricting ρ, we can assume that B ρ (x(τ 1 )) ⊂ B r ′ (a). Since u k (x, t(τ 1 ) is uniformly bounded when x ∈ B ρ (x(τ 1 )), it follows from Proposition 2.3 that for any ρ ′ ∈ (0, ρ), u k remains uniformly bounded in

Q t(τ 1 ),t(τ 1 )+λ ρ ′ ,x(τ 1 )
. Moreover, for any compact subset K ∈ B c ρ ′ (x τ 1 ), there holds

u k (x, t(τ 1 ) + λ) ≤ C (λβ) 1/(q-1)
, where β = min{h(x, t) : (x, t) ∈ K × [t(τ 1 ), t(τ 1 ) + λ]}. Iterating this construction, we can construct a finite number of cylinders Q t(τ 1 )+jλ,t(τ 1 )+(j+1)λ ρ ′ ,a j containing Γ and in which u k remains uniformly bounded. Since local uniform boundedness holds also outside such cylinders, the proof follows.

Remark. In full generality we conjecture that if γ is C 1 with γ ′ (τ ) = 0 and there exists τ 0 ∈ (0, T ) such that t(τ ) admits a local strict maximum on the right at τ 0 ∈ (0, T ), then u ∞ is locally bounded Q T \ Γ τ 0 .

Propagation of singularities in the space

In this section, we assume that the degeneracy curve Γ is parametrized by the variable t ∈ [0, T ] and defined by (2.10) with |x(t)| > 0 if t > 0. We denote by

d P [(x, t), (y, s)] := |x -y| + √ t -s if t ≥ s
the parabolic distance and we assume that h(x, t) depends on d P [(x, t), Γ] under the following form h(x, t) = e -ℓ(d P [(x,t),Γ]) (4.1)

where ℓ ∈ C([0, ∞)) is positive, nonincreasing and lim r→0 ℓ(r) = ∞. For ǫ > 0, we recall

T Γ T 0 ǫ
denotes the ǫ -spherical tubular neighborhood of Γ between t = 0 and t = T defined by T

Γ T 0 ǫ := {(x, t) ∈ Q T : |x -x(t)| < ǫ}. (4.2)
The basis of T

Γ T 0 ǫ , in R N ×{0}, is the ball B ǫ . Since d P [(x, t), Γ] ≤ |x-x(t)|, d P [(x, t), Γ] ≤ ǫ in T Γ T 0 ǫ and ℓ(ǫ) ≤ ℓ(d P [(x, t), Γ]). Then, u ∞ is bounded from below in Γ ǫ by the solution v ǫ of          ∂ t v ǫ -∆v ǫ + e -ℓ(ǫ) v p ǫ = 0 in T Γ T 0 ǫ v ǫ = 0 in ∂ P T Γ T 0 ǫ v ǫ = ∞δ 0 in B ǫ , (4.3)
where

∂ P T Γ T 0 ǫ := {(x, t) : |x -x(t)| = ǫ}
is the lateral parabolic boundary of T Γ T 0 ǫ . The formal aspect comes from the fact that the existence of v ǫ has to be proved.

For the sake of simplicity we assume that 1 < p < 1 + 2 N , the case p ≥ 1 + 2 N needing a simple adaptation of the type which is developed in the proof of Theorem D-Case 2.

We consider the following change of variable x = y + x(t) and v ǫ (x, t) = ṽǫ (y, t). Then ṽǫ satisfies in

Q Bǫ τ := B ǫ × (0, τ )        ∂ t ṽǫ -∆ṽ ǫ + x ′ (t)|∇ṽ ǫ + e -ℓ(ǫ) ṽp ǫ = 0 in Q Bǫ τ ṽǫ = 0 in ∂ P Q Bǫ τ ṽǫ = ∞δ 0 in B ǫ . (4.4)
In particular u(x(α), α) ≥ v ǫ (x(α), α) = ṽǫ (0, α). We set ω ǫ (x, t) = ǫ 2/(p-1) e -ℓ(ǫ)/(p-1) ṽǫ (ǫx, ǫ 2 t), and x ǫ (t) = ǫ -1 x(ǫ 2 t). Then ω ǫ satisfies

       ∂ t ω ǫ -∆ω ǫ + x ′ ǫ (t)|∇ω ǫ + ω p ǫ = 0 in Q B 1 ǫ -2 τ ω ǫ = 0 in ∂ P Q B 1 ǫ -2 τ ω ǫ = ∞δ 0 in B 1 , (4.5)
and for any α > 0

u ∞ (x(α), α) ≥ ǫ -2/(p-1) e ℓ(ǫ)/(p-1) ω ǫ (0, ǫ -2 α). (4.6)
Therefore, the problem is reduced to showing that

lim ǫ→0 ǫ -2/(p-1) e ℓ(ǫ)/(p-1) ω ǫ (0, ǫ -2 α) = ∞. (4.7)
We associate the following parabolic equation [START_REF] Friedman | Partial differential equations of parabolic type[END_REF], [START_REF] Ladyshenskaya | Linear and quasilinear equations of parabolic type[END_REF]).

       ∂ t v -∆v + β(t)|∇v + v p = 0 in Q B 1 ∞ v = 0 in ∂ P Q B 1 ∞ v(y, 0) = v 0 in B 1 . (4.8) Existence of solution is classical when β ∈ L ∞ loc ([0, ∞)) (see
Proposition 4.1 Let β ∈ L ∞ loc ([0, ∞)).
Then the following estimates holds v(., t) L 2 ≤ e -λ 0 t v(., 0) L 2 (4.9)

where λ 0 is the first eigenvalue of -∆ in B 1 , and, for some

C = C(N ) > 0, v(., t) L ∞ ≤ C min{t -N/4 , e -λ 0 t } v(., 0) L 2 . (4.10)
Proof. Since for any continuous function r → g(r) = G ′ (r), we have

B 1 β(t)|∇v g(v)dx = B 1 β(t)|∇G(v) dx = - B 1 G(v) divβ(t)dx = 0, we have 2 -1 d dt B 1 v 2 dx + B 1 |∇v| 2 dx ≤ 0 =⇒ d dt B 1 v 2 dx ≤ -2λ 0 B 1 v 2 dx.
Thus (4.9) follows. Furthermore, for any q > 1,

1 q + 1 d dt B 1 |v| q+1 dx + 4q (q + 1) 2 B 1 |∇v (q+1)/2 |dx ≤ 0,
it follows that t → v(., t) L q+1 is decaying. Therefore, using Gagliardo-Nirenberg estimate,

d dt v q+1 L q+1 + 4qC q + 1 v q+1 L (q+1)N/(N-2) ≤ 0.
We assume N ≥ 3, the cases N = 1, 2 needing a simple modification. Thus, for any

t > s > 0, 4qC(t -s) q + 1 v(., t) q+1 L (q+1)N/(N-2) ≤ v(., s) q+1 L q+1 .
By a standard use of Moser iterative scheeme, we derive

v(., t) L ∞ ≤ Ct -N/4 v(., 0) L 2 . (4.11)
Consequently, for any 0 < s < t, we have

v(., t) L ∞ ≤ C(t -s) -N/4 v(., s) L 2 ≤ C(t -s) -N/4 e -λ 0 s v(., 0) L 2 .
Noting that, if t > N/4λ 0 , min{(ts) -N/4 e -λ 0 s 0 < s < t} = e N/4 (4λ 0 /N ) N/4 e -λ 0 t , we derive (4.10). and the proof follows.

In the sequel we denote by λ β the first eigenvalue of L β , thus

λ β = λ 0 + |β| 2 4 . If ψ β is a corresponding positive eigenfunction, then ψ β (x) = e 1 2 β|x ψ 0 (x)
where ψ 0 is a positive first eigenfunction of -∆ in H 1,0 (Ω) Proposition 4.3 Under the assumptions of Proposition 4.2 and for p > 1, we denote by v be the solution of

       ∂ t v -∆v + β|∇v + v p = 0 in Q B 1 ∞ v = 0 in ∂ P Q B 1 ∞ v(y, 0) = v 0 in B 1 .
(4.14)

where v 0 ∈ L 2 (B 1 ) is nonnegative. Then there exists some c = c(v 0 ) > 0 such that lim t→∞ e λ β t v(., t) = cψ (4.15)
uniformly in B 1 .

Proof. We write v(x, t) = e 1 2 β|x w(x, t), thus (4.14) turns into

∂ t w -∆w + |β| 2 4 w + e p-1 2
β|x w p = 0. (4.16)

The proof of [3, Th 3.1] applies easily and the result follows.

Proposition 4.4 Assume β ∈ W 1,∞ loc ([0, ∞)). For τ > 1, we set sup 1≤t≤τ |β(t)| = β τ and sup 1≤t≤τ |β ′ (t)| = δ τ . If v is the solution of (4.8) where v 0 ∈ L 2 (B 1 ) is nonnegative, we denote σ(τ ) = sup e p-1 2 |β(t)| w p-1 (x, t) : (x, t) ∈ B 1 × [1, τ ] . Then if v(., 1) satisfies c 1 ψ 0 ≤ v(., 1), there holds v(x, t) ≥ c 1 e -(λ 0 + β 2 τ 4 + δτ 2 +στ )(t-1) ψ 0 (x) ∀(x, t) ∈ B 1 × [1, τ ] (4.17) Proof. Let w(x, t) = e -1 2 β(t)|x v(x, t). Then ∂ t w -∆w + |β| 2 4 + 1 2 β ′ |x + e p-1 2
β(t)|x w p-1 w = 0. (4.18) Since |x| ≤ 1 there holds in B 1 × [1, τ ] |β(t)| 2 4 + 1 2 β ′ (t)|x + e p-1 2
β(t)|x w p-1 ≤ β 2 τ 4 + δ τ 2 + σ τ . (4.19)
Therefore

∂ t w -∆w + β 2 τ 4 + δ τ 2 + σ τ w ≥ 0 in B 1 × [1, τ ]. (4.20) Since (x, t) → e -(λ 0 + β 2 τ 4 + δτ 2
+στ )(t-1) ψ 0 (x) satisfies the equation associated to (4.20), (4.17) follows.

Proof of Theorem C. Step 1: Initialization of the blow-up. With our previous notations,

β(t) = β ǫ (t) = x ′ ǫ (t) = ǫx ′ (ǫ 2 t) and β ′ ǫ (t) = ǫ 3 x ′′ (ǫ 2 t). Since x ′ ǫ is locally bounded, it follows
by Hopf lemma that there exists c 1 > 0 such that ω ǫ (., 1) ≥ c 1 ψ 0 . We take τ = ǫ -2 α where α > 0 is fixed. Then

β τ = sup{|ǫx ′ (ǫ 2 t)| : 1 ≤ t ≤ ǫ -2 α} = ǫ sup{|x ′ (t)| : ǫ 2 ≤ t ≤ α}, δ τ = sup{|ǫ 3 x ′′ (ǫ 2 t)| : 1 ≤ t ≤ ǫ -2 α} = ǫ 3 sup{|x ′′ (t)| : ǫ 2 ≤ t ≤ α}.
Therefore ǫ

-2 p-1 e ℓ(ǫ) p-1 ω ǫ (0, ǫ -2 α) ≥ c ′ 1 e A(ǫ) ψ 0 (0) (4.21)
where

A(ǫ) := - 2 p -1 ln ǫ + ℓ(ǫ) p -1 -λ 0 + β 2 τ 4 + δ τ 2 + σ τ α ǫ 2 (4.22)
Since lim inf ǫ→0 ǫ 2 ℓ(ǫ) > 0 there exists α 0 > 0 such that for any 0 < α < α 0 , there holds lim ǫ→0 A(ǫ) = ∞. This implies

u ∞ (x(α), α) = ∞ ∀ 0 < α < α 0 . (4.23)
Step 2: Propagation. In order to prove that the blow-up propagates along Γ we have replace t = 0 by t = α < α 0 . We claim that

Bσ(x(α)) u ∞ (x, α)dx = ∞ ∀σ > 0. (4.24)
Actually, it is sufficient to prove the result with σ = ǫ and with u ∞ replaced by v ǫ . Then

Bǫ(x(α) v ǫ (x, α)dx = Bǫ ṽǫ (x, α)dx = ǫ N -2 p-1 e ℓ(ǫ) p-1 B 1 ω ǫ (x, ǫ -2 α)dx (4.25)
Using (4.17), (4.21) we have

ǫ N -2 p-1 e ℓ(ǫ) p-1 B 1 ω ǫ (x, ǫ -2 α)dx ≥ c ′ 1 e A ′ (ǫ) B 1 ψ 0 (x)dx
where

A ′ (ǫ) := (N - 2 p -1 ) ln ǫ + ℓ(ǫ) p -1 -λ 0 + β 2 τ 4 + δ τ 2 + σ τ α ǫ 2 (4.26)
Thus lim ǫ→0 A ′ (ǫ) = ∞ for any α < α 0 . This implies the claim.

Step 3: End of the proof. For k > 0, we denote by u α,kδ the solution of

∂ t u -∆u + h(x, t)u p = 0 in R N × (α, ∞) u(., α) = kδ x(α) in R N . (4.27) We claim that u ∞ (x, t)) ≥ u α,kδ (x, t) ∀(x, t) ∈ R N × [α, ∞). (4.28)
We fix k > 0, then for any σ > 0, there exists m = m(σ) > 0 such that

Bσ(x(α)) min{m, u ∞ (x, α)}dx = k.
Furthermore lim σ→0 m(σ) = ∞. Let u α,k,σ be the solution of

∂ t u -∆u + h(x, t)u p = 0 in R N × (α, ∞) u(., α) = min{m, u ∞ (., α)}χ Bσ (x(α)) in R N . (4.29)
By the maximum principle u ∞ ≥ u α,k,σ in R N × (α, ∞). But min{m, u ∞ (., α)}χ Bσ (x(α)) converges to kδ x(α) in the weak sense of measure when σ → 0. By stability, since we have assumed p < 1 + 2 N , u α,k,σ → u α,kδ locally uniformly in R N × (α, ∞) (see [START_REF] Marcus | Initial trace of positive solutions to semilinear parabolic inequalities[END_REF] e.g.). Therefore (4.28) follows. Since k is arbitrary, it follows

u ∞ (x, t)) ≥ u α,∞δ (x, t) ∀(x, t) ∈ R N × [α, ∞). (4.30)
By

Step 1, u α,∞δ (x(α + β), α + β) is infinite for 0 < β < α 0 . This implies that (4.23) holds for any 0 < α < 2α 0 . Iterating this process we conclude that u ∞ blows-up on whole Γ and (4.24) holds for any α ∈ (0, T ].

Propagation of singularities in the initial plane

In this section, we consider

∂ t u -∆u + h(x, t)u p = 0 in Q ∞ := R N × (0, ∞) (5.1)
where h ∈ C(Q ∞ ). We set x = (x 1 , ..., x N ) = (x 1 , x ′ ) and we suppose that h(x, t) > 0 except when (x, t) belongs to some straight line Γ that we can assume to be the x 1 axis in the plane t = 0. We set d ∞ ((x, t), Γ) = max{ √ t, |x ′ |} and write h under the form h(x, t) = e -ℓ(d∞(x,t);Γ) (5.2)

where ℓ : (0, ∞) → (0, ∞) is continuous and nonincreasing with limit ∞ at 0.

Proof of Theorem D. Case 1: 1 < p < 1 + 2 N . For ǫ > 0, we consider the "tunnel" with axis Γ defined by

T ǫ := {(x, t) : x 1 ∈ R, (x ′ , t) ∈ B ′ ǫ × (0, ǫ 2 )} where B ′
ǫ is the ball in R N -1 with center 0 and radius ǫ. Since ℓ is decreasing, there holds

∂ t u -∆u + e -ℓ(ǫ) u p ≥ 0 in T ǫ . Thus u ∞δ 0 ≥ v ∞δ 0 in T ǫ ,
where v ∞δ 0 = lim k→∞ v kδ 0 and v kδ is the solution of

   ∂ t v -∆v + e -ℓ(ǫ) v p = 0 in T ǫ v = kδ 0 in R × B ′ ǫ v = 0 in R × ∂B ′ ǫ .
(5.3)

We put v ∞δ (x, t) = ǫ -2/(p-1) e ℓ(ǫ)/(q-1) w(x/ǫ, t/ǫ 2 )

Then w = w ǫ satisfies    ∂ t w -∆w + w p = 0 in R × B ′ 1 × (0, 1) w = ∞δ 0 in R × B ′ 1 w = 0 in R × ∂B ′ 1 .
(5.4)

We denote by W the solution of

   ∂ τ W -∆W + W = 0 in R × B ′ 1 × (0, 1) W (ξ 1 , ξ ′ , 0) = ψ(ξ 1 )φ(ξ ′ ) in R × B ′ 1 W = 0 in R × ∂B ′ 1 × (0, 1). (5.5)
where φ is the first eigenfunction of -∆ ξ ′ in W 1,2 0 (B ′ 1 ) with maximum 1 and corresponding eigenvalue λ and ψ(ξ

1 ) = cos(ξ 1 )χ [-π 2 , π 2 ] (ξ 1 ). Then W (ξ 1 , ξ ′ , τ ) = e -(λ+1)τ φ(ξ ′ ) √ 4πτ π/2 -π/2 e -|ξ 1 -ζ| 2 /4τ ψ(ζ)dζ.
Since 0 ≤ W ≤ 1, W is a subsolution for the equation satisfied by w and there exists a > 0 and c > 0 such that

w ǫ (ξ, τ + a) ≥ cW (ξ, τ ) ∀(x, t) ∈ R × B ′ 1 × (0, 1).
Returning to v ∞ δ 0 , we derive v ∞δ 0 (x, t + aǫ 2 ) ≥ cǫ -2/(p-1) e ℓ(ǫ)/(q-1) W (x/ǫ, t/ǫ 2 )

∀(x 1 , x ′ , t) ∈ R × B ′ ǫ × (0, ǫ 2 ], (5.6)
which implies, with t = ǫ 2 and x ′ = 0, v ∞ (x 1 , 0, (a + 1)ǫ 2 ) ≥ cǫ -2/(p-1) e ℓ(ǫ)/(q-1) e -λ-1 φ(0)

√ 4π π/2 -π/2 e -|x 1 /ǫ-ζ| 2 /4 ψ(ζ)dζ. (5.7) But π/2 -π/2 e -|x 1 /ǫ-ζ| 2 /4 ψ(ζ)dζ ≥ e -x 2 1 /2ǫ 2 π/2 -π/2 e -|ζ| 2 /2 ψ(ζ)dζ (5.8)
If we fix in particular |x 1 | ≤ δ where

|x 1 | < δ = 2ǫ 2 ℓ(ǫ) q -1 , we derive lim ǫ→0 v ∞ (x 1 , 0, (a + 1)ǫ 2 ) = ∞.
(5.9) Furthermore this limit is uniform for x 1 ∈ [-δ ′ , δ ′ ], where δ ′ < δ. Furthermore the interval [-δ ′ , δ ′ ] does not shrink to {0} when ǫ → 0, since it is assumed that lim inf ǫ→0 ǫ 2 ℓ(ǫ) > 0.

(5.10)

Replacing [-δ ′ , δ ′ ] by [δ ′ , 3δ ′ ] and [-3δ ′ , -δ ′ ] and iterating, we conclude that lim t→0 u ∞δ 0 (x 1 , 0, t) ≥ lim t→0 v ∞ (x 1 , 0, t) = ∞.

(5.11)

From this, it is easy to obtain that u ∞δ 0 (x, t) = u ∞δ 0 (0, x ′ , t) is independent of x 1 and coincide with U (x ′ , t) where U is the solution of (5.13) where l(s) = ℓ(s)γ ln s and γ > N (p -1) -2. Then γ > 0, (1.3) is satisfied and for any k > 0 there exists a unique solution v = v kδ to

     ∂ t U -∆U + e -ℓ(max{ √ t,|x ′ |}) U p = 0 in (0, ∞) × R N -1 U = ∞δ 0 in R N -
     ∂ t v -∆v + e -l(ǫ) max{ √ t, |x ′ |} γ v p = 0 in T ǫ v = kδ 0 in R × B ′ ǫ v = 0 in R × ∂B ′ ǫ .
(5.14) Furthermore u ∞ (x, t) ≥ v ∞δ in T ǫ . We set v ∞δ (x, t) = ǫ -(2+γ)/(p-1) e ℓ(ǫ)/(q-1) w(x/ǫ, t/ǫ 2 ), and w = w ǫ satisfies

   ∂ τ w -∆w + (max{ √ τ , |ξ ′ |}) γ w p = 0 in R × B ′ 1 × (0, 1) w = ∞δ 0 in R × B ′ 1 w = 0 in R × ∂B ′ 1 × (0, 1).
(5.15)

We denote by W the solution of

   ∂ τ W -∆W + (max{ √ τ , |ξ ′ |}) γ W = 0 in R × B ′ 1 × (0, 1) W (ξ 1 , ξ ′ , 0) = ψ(ξ 1 )φ(ξ ′ ) in R × B ′ 1 W = 0 in R × ∂B ′ 1 × (0, 1).
(5.16)

where ψ and φ are as in the first case. This equation admits a separable solution W (ξ, τ ) = W 1 (ξ 1 , τ )W ′ (ξ ′ , τ ) where

   ∂ τ W ′ -∆ ξ ′ W ′ + (max{ √ τ , |ξ ′ |}) γ W ′ = 0 in B ′ 1 × (0, 1) W ′ (ξ ′ , 0) = φ(ξ ′ ) in B ′ 1 W ′ = 0 in ∂B ′ 1 × (0, 1),
(5.17 The exact expression of W ′ is not simple but since (max{ √ τ , |ξ ′ |}) γ ≤ 1 in B ′ 1 × (0, 1), there holds

W (ξ 1 , ξ ′ , τ ) ≥ e -(λ+1)τ φ(ξ ′ ) √ 4πτ π/2 -π/2
e -|ξ 1 -ζ| 2 /4τ ψ(ζ)dζ.

From this point, and since (5.10) holds with ℓ replaced by l, the proof is the same as in Case 1.

Proposition 4 . 2 4 : 1 2

 4241 Let Ω be a bounded open domain in R N , β ∈ R N and L β the operator v → -∆v + β|∇v . Then the spectrum of L β in H 1,0 (Ω) is the given by σ(L β ) = λ + |β| 2 λ ∈ σ(L 0 ) . (4.12) Proof. Put v(x) = e β|x w(x). Then ∇v(x) = e

) and ∂ τ W 1 -∂ x 1 x 1 W 1 + 2 -π/ 2 e

 11122 W 1 = 0 in R × (0, 1) W 1 (ξ 1 , 0) = φ(ξ 1 ) in R (5.18) Thus W (ξ 1 , ξ ′ , τ ) = W ′ (ξ ′ , τ ) √ 4πτ π/-|ξ 1 -ζ| 2 /4τ ψ(ζ)dζ.

= ∅, we obtain that u ∞ (x, t * ) ≤ µ for