Emeric Bouin
email: bouin@ceremade.dauphine.fr

Christopher Henderson
email: christopher.henderson@ens-lyon.fr

Super-linear spreading in local bistable cane toads equations

Keywords: Structured populations, reaction-diffusion equations, front acceleration AMS Class. No: 35Q92, 45K05, 35C07

In this paper, we study the influence of an Allee effect on the spreading rate in a local reaction-diffusion-mutation equation modeling the invasion of cane toads in Australia. We are, in particular, concerned with the case when the diffusivity can take unbounded values. We show that the acceleration feature that arises in this model with a Fisher-KPP, or monostable, non-linearity still occurs when this non-linearity is instead bistable, despite the fact that this kills the small populations. This is in stark contrast to the work of Alfaro, Gui-Huan, and Mellet-Roquejoffre-Sire in related models, where the change to a bistable non-linearity prevents acceleration.

1 Introduction

The biological model

The invasion of cane toads in Australia has interesting features different from the standard spreading observed in most other species. The experimental data [START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF][START_REF] Shine | An evolutionary process that assembles phenotypes through space rather than through time[END_REF] show that the invasion speed has steadily increased during the eighty years since the toads were introduced in Australia. In addition, the younger individuals at the edge of the invasion front have a significantly different morphology compared to other populations -their legs tend to be on average much longer than away from the front. This is just one example of a non-uniform space-trait distribution -see, for instance, a study on the expansion of bush crickets in Britain [START_REF] Thomas | Ecological and evolutionary processes at expanding range margins[END_REF]. Several works have addressed the front invasions in ecology, where the trait is related to the dispersal ability [START_REF] Arnold | Existence of nontrivial steady states for populations structured with respect to space and a continuous trait[END_REF][START_REF] Champagnat | Invasion and adaptive evolution for individual-based spatially structured populations[END_REF]. It has been observed that selection of more mobile individuals can occur, even if they have no advantage in their reproductive rate, due to the spatial sorting [START_REF] Thomas | Changes in dispersal during species range expansions[END_REF][START_REF] Kokko | From individual dispersal to species ranges: Perspectives for a changing world[END_REF][START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF][START_REF] Ronce | How does it feel to be like a rolling stone? Ten questions about dispersal evolution[END_REF].

In this paper, we focus on a model for this phenomenon (proposed in [START_REF] Bénichou | Front acceleration by dynamic selection in fisher population waves[END_REF], based on the classical Fisher-KPP equation [START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kolmogorov | Étude de l'équation de la chaleurde matière et son application à un problème biologique[END_REF] ). The population density is structured by a spatial variable, x ∈ R, and a motility variable θ ∈ Θ def = [θ, ∞), with a fixed θ > 0. This population undergoes diffusion in the trait variable θ, with a constant diffusion coefficient α > 0, representing mutation, and in the spatial variable, with the diffusion coefficient θ, representing the effect of the trait on the spatial spreading rates of the species.

In the model proposed in [START_REF] Bénichou | Front acceleration by dynamic selection in fisher population waves[END_REF], the reaction is of non-local, monostable type. In the case where Θ is bounded, the problem has been shown to be well-posed [START_REF] Turanova | On a model of a population with variable motility[END_REF] and exhibit travelling-wave solutions [START_REF] Bouin | Travelling waves for the cane toads equation with bounded traits[END_REF]. When Θ is unbounded, sharp results have been proven for the local analogue while weaker results have been obtained for the non-local equation [START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF]. In particular, the Cauchy problem exhibits propagation of of the order t 3/2 .

Motivation

With a monostable non-linearity, acceleration, i.e. super-linear spreading of level sets, has been shown in various settings in mathematical biology [7, 11,13, 14,16, 18,20-22, 26,32]. We mention, in particular, the existence of acceleration for Fisher-KPP with fat-tailed initial data and for Fisher-KPP with the fractional Laplacian. In these models, acceleration is related to the notion of pulled fronts. In particular, we see that small populations far from the origin grow exponentially, causing the acceleration. Recently, in these settings, it has been shown that when the non-linearity is ignition or bistable type, there is no acceleration [START_REF] Alfaro | Slowing allee effect vs. accelerating heavy tails in monostable reaction diffusion equations[END_REF][START_REF] Gui | Traveling wave solutions to some reaction diffusion equations with fractional Laplacians[END_REF][START_REF] Mellet | Existence and asymptotics of fronts in non local combustion models[END_REF]. Indeed, the small populations that drive the acceleration are killed when the non-linearity is bistable. In other words, acceleration is a tail phenomenon, so when a bistable non-linearity kills the tails, there is no acceleration.

In this article, we obtain results in stark contrast with the results discussed above. Indeed, we demonstrate propagation of the order t 3/2 in the local cane toads equation with a bistable nonlinearity. This shows that, for the cane toads equation, acceleration is a bulk phenomenon, i.e. that it is not driven by small populations far from the origin.

Main results

Fix a bistable non-linearity f : [0, 1] → [0, ∞] that is Lipschitz continuous function such that there exists α ∈ (0, 1 2 ) with

f (u) ≥ u(u -α)(1 -u).
We assume also that f (0) = f (1) = 0. We are interested in the long-time asymptotics of solutions to the Cauchy problem

         u t = θu xx + u θθ + f (u), (t, x, θ) ∈ R + × R × Θ, u θ (t, x, θ) = 0, (t, x) ∈ R + × R, u(0, x, θ) = u 0 (x, θ) (x, θ) ∈ R × Θ. (1.1)
The initial data function u 0 is assumed to satisfy

u 0 ≥ 1 R -×[θ,(1+λ)θ] ,
where λ > 0.

Our interest is in understanding where the "front" of u is. In other words, we will fix a level set and understand its propagation. Using the same equation with a monostable non-linearity as a super-solution to (1.1), the results from [START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] show that no level set can move faster than O(t 3/2 ). Our main result is to show that a lower bound of the same order holds as well.

Theorem 1.1. Fix any m ∈ (0, 1). There exists λ 0 ∈ R * + , depending only on α and θ, and γ ∈ R * + , depending only on α, such that if λ ≥ λ 0 then lim inf

t→∞ max{x : ∃θ ∈ Θ, u(t, x, θ) = m} t 3/2 ≥ γ.
Before we continue, we note that we expect that this theorem holds in much greater generality. However, our interest is in showing that the acceleration in the cane toads equation is a also bulk phenomenon, as opposed to a tail phenomenon. Hence, we seek to simply show that acceleration occurs despite the bistable non-linearity, and we do not attempt to find the most general setting or the sharpest bounds.

To prove this theorem, we will follow a similar strategy as in a previous paper by the authors and Ryzhik [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] and that also appears in [START_REF] Berestycki | Asymptotic spreading in heterogeneous diffusive excitable media[END_REF][START_REF] Roquejoffre | Lecture notes in a Toulouse school on KPP and probability[END_REF]. In rough words, we will slide a suitable "bump" along suitable trajectories in the phase plane R × Θ, making sure that the ball remains below the solution of the original Cauchy problem (1.1). Once the trajectories are well chosen, this will imply the acceleration phenomena claimed in theorem 1.1. However, the "bump" sub-solution is significantly more complicated to create in this setting than in the monostable one. Indeed, in [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF], the sub-solution is created almost entirely with the linearized (around zero) problem. In our setting, of course, the linearized problem decays to zero on any traveling ball. We describe how to overcome this difficulty below. This is the main technical difficulty in the present article.

Strategy of the proof: constructing the sub-solution

We will now end this introduction by describing and presenting more precisely the major objects that are needed to prove theorem 1.1, namely the suitable trajectory and "bump". To begin, we fix a large T > 0 and any level set height m ∈ (0, 1). We then define a trajectory in the space R × Θ by t → (X T (t), Θ T (t)), for any t ∈ [0, T ]. The coordinate functions (X T , Θ T ) will be defined later on (see section 2), and are crucial for our analysis. In order to slide a bump along this trajectory, we define the moving (growing) ellipse

E t,Λ = (x, θ) ∈ R × Θ : (x -X T (t)) 2 Θ T (t) + (θ -Θ T (t)) 2 ≤ Λ 2 (1.2)
and the moving (growing) annulus

A t,Λ = (x, θ) ∈ R × Θ : Λ 2 ≤ (x -X T (t)) 2 Θ T (t) + (θ -Θ T (t)) 2 ≤ 4Λ 2 , (1.3)
where Λ is a positive constant to be chosen later that encodes the sizes of this two objects. For now, our only assumption on Λ is that Λ ≤ λθ/8. Our goal is to build a sub-solution to (1.1) on the bigger ellipse A t,Λ ∪ E t,Λ , what we called a "bump" above. Since it is the main issue of the paper, let us now explain how we build a sub-solution on A t,Λ ∪ E t,Λ . We shall patch together a solution w

+ ≥ α on E t,Λ of    w + t = θw + xx + w + θθ + f r (w + ), (x, θ) ∈ E t,Λ , w + ≡ α, (x, θ) ∈ ∂E t,Λ , (1.4)
and a positive solution

w -≤ α on A t,Λ of          w - t = θw - xx + w - θθ + f r (w -), (x, θ) ∈ A t,Λ , w -≡ α, (x, θ) ∈ ∂E t,Λ , w -≡ 0, (x, θ) ∈ ∂A t,Λ \ ∂E t,Λ . (1.5)
Notice that in these two previous definitions, there is a small but important discrepancy between the nonlinearities used and in (1.1). The new non-linearity f r is defined as follows for any r ∈ (2α, 1]:

f r (u) def =    u(u -α)(1 -u), for u ≤ α, c r u(u -α)(r -u), for u ≥ α. (1.6)
where we define c r = (1 -α)(rα) -1 . It is easily verified that f r (u) ≤ f (u) for all u, and it is clear that f r is Lipschitz continuous. For technical reasons, we are required to take r to be slightly less than 1 in the sequel. The bolded green line denotes the set A t , the red dotted one the set E t . Now that we have defined w ± , we obtain a sub-solution w of (1.1) on R + × R × Θ defined by

α w + 0 w - A t,Λ A t,Λ E t,Λ
w(t, x, θ) = w + 1 E t,Λ + w -1 A t,Λ . (1.7)
By construction, w is Lipschitz continuous. However, it need not be a C 1 function along ∂A t,Λ \∂E t,Λ and ∂E t,Λ . In order to be a sub-solution, we must check the convexity of w along both boundaries. The positivity of w -ensures that the convexity of w is correct at the boundary ∂A t,Λ \ ∂E t,Λ . However, to make sure that it is indeed a sub-solution at the boundary ∂E t,Λ we need to check properly that the normal derivatives along this boundary are well ordered:

|∂ n w + | ≥ |∂ n w -| on ∂E t,Λ , (1.8)
where ∂ n is the (outward) normal derivative to the boundary of E t,Λ . This is the main technical issue at hand in the proof.

Numerics and comments

Numerics

Let us enlighten theorem 1.1 by showing in fig. 1.2 some numerical insights of the acceleration behavior. One can compare these to [START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF][START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF] that are numerics related to the study of the acceleration in the nonlocal cane toads equation.

A nonlocal bistable model

We conclude the introduction by discussing a nonlocal version of (1.1):

                   n t = θn xx + n θθ + n(ρ -α)(1 -ρ), (t, x, θ) ∈ R + × R × Θ, ρ(t, x) := ˆ∞ θ n(t, x, θ) dθ, (t, x) ∈ R + × R, n θ (t, x, θ) = 0, (t, x) ∈ R + × R, n(0, x, θ) = n 0 (x, θ) (x, θ) ∈ R × Θ.
(1.9)

It is easy to see that one many bound the propagation rate above by using the monostable model with growth rate r = sup ρ∈R + ((ρα)(1 -ρ)). This implies that the propagation can be no faster than O(t 3/2 ). However, obtaining a lower bound is significantly more complicated. We are led to the following conjecture, which we are unable to prove at this time.

Conjecture 1.2. The model (1.9) exhibits acceleration. That is, the level sets of ρ move superlinearly in time.

We provide some numerics in fig. 1.3 to support this conjecture. However, the proof of such a result is far beyond the scope of this paper. Indeed, it requires to study the dynamics of the zone where ρ is greater that α, which is not at the edge of the front at all, so that an argument by contradiction as in [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] is impossible. Moreover, to our knowledge, it is not possible to derive a probabilistic framework such as in [START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF] (and the references therein) to be able to study the nonlocal model (1.9). As new tools will be needed to study this model, we thus leave this conjecture for further investigation.

Outline of the article

In section 2, we give the main outline of the proof of theorem 1.1 assuming the main proposition 2.1. In section 3, we prove this proposition by finishing the construction of the subsolution (1.7), and reducing it to showing a proposition 3.2 that gives (1.8). This main lemma involves understanding a related steady state problem, which we analyze in sections 4 and 5. One can track the accelerated behavior, for example on the space axis. The invasion at the back in the θ-direction is expected to be linear in time. This pattern is very similar as for the monostable cane toads equation, see also [START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF][START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF]. 2 The acceleration feature -proof of theorem 1.1

In this section, we prove theorem 1.1. For this purpose, we first define carefully the trajectory that has been introduced in the previous section. Second, we take for granted the existence of a sub-solution of the form (1.7) (this will the object of section 3) and we show how this implies the acceleration feature stated in theorem 1.1.

We now define the trajectory t

∈ [0, T ] → (X T (t), Θ T (t)) ∈ R × Θ.
For this purpose, we concatenate two trajectories t → (X

T (t), Θ (1)
T (t)) and t → (X (2) T (t), Θ (2)
T (t)):

X T =    X (1) T , t ∈ 0, T 2 , X (2)
T , t ∈ T 2 , T , Θ T =    Θ (1)
T , t ∈ [0, T 2 ], Θ (2)
T , t ∈ [ T 2 , T ].
These latter trajectories are defined by

∀t ∈ 0, T 2 , X (1)
T , Θ

T (t) = - λ 4 θ, ct + 1 + 3λ 4 θ , (1)
and

∀t ∈ T 2 , T , X (2)
T , Θ

T (t) =   c t -T 2 T cT 2 + 1 + 3λ 4 θ - λ 4 θ, cT 2 + 1 + 3λ 4 θ   . (2)
In related works [START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] about the standard local and non-local cane toads equations, some trajectories of the same kind have been introduced.

Before continuing, we point out that our choice of Λ ≤ λθ/8 gives us that

E 0,Λ ∪ A 0,Λ ⊂ (-∞, 0] × [θ, (1 + λ)θ].
This implies that w will initially sit below u, allowing us to apply the maximum principle. With this in mind, we state in (a very minimalist fashion) a proposition that will imply our main theorem 1.1.

Proposition 2.1. Fix m ∈ (0, 1). There exist positive constants λ 0 , Λ 0 , T 0 , and c 0 , such that if Λ ≥ Λ 0 , T ≥ T 0 , and c ≤ c 0 , then there exists a sub-solution w of the form (1.7) satisfying (1.8)

such that ∀(x, θ) ∈ E T, Λ 2 , w + (T, x, θ) ≥ m. (2.1)
In addition, all constants λ 0 , Λ 0 , c 0 , and T 0 are bounded for m bounded away from 1, and λ 0 depends on Λ.

x θ

(X T , Θ T ) T 2 (X T , Θ T ) (0) (X T , Θ T ) (T ) (1 + λ)θ θ Figure 2
.1: Plot of the trajectory in the phase plane (x, θ). We emphasize that the first part of the trajectory corresponds to a movement towards large traits only in the θ-direction, whereas the second part of the trajectory corresponds to a movement towards large space positions only in the x-direction, at an accelerated spreading rate. The red dotted line is the support of the initial condition. The red bold line is the initial support of the subsolution w : it sits inside the support of the initial condition.

We will prove this proposition in section 3. Let us now conclude the proof of theorem 1.1. We simply write, as a consequence of the comparison principle, the following:

lim inf T →∞ max {x : ∃θ ∈ Θ, u(T, x, θ) ≥ m} T 3/2 ≥ lim inf T →∞ max {x : ∃θ ∈ Θ, w + (T, x, θ) ≥ m} T 3/2 ≥ lim inf T →∞ X (2)
T (T )

T 3/2 = lim inf T →∞     c 2 T 2 T 3/2 cT 2 + 1 + 3λ 4 θ - 1 + 3λ 4 θ T 3/2     = c 2 .
After calling the latter constant γ, this finishes the proof.

Remark 2.2. To find optimal constants, one could write a general framework for the trajectories such as for the cane toads equation in [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] and try to optimize afterwards. However, in this paper, the suboptimal trajectories are not the only reason for sub-optimality. One would have to quantify all the constants in various lemmas in this paper, that are not strongly related to the definition of the trajectories themselves. Since this is not our main interest, we opt for the simplest possible presentation.

3 Building the sub-solution w along trajectories

In this section, we will prove proposition 2.1, by working on the two components of the trajectories separately. Since we want to build a sub-solution along a trajectory (X T , Θ T ), we work in the moving coordinates introduced in [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF])

y = x -X T √ Θ T , η = θ -Θ T .
The moving ellipse and annulus E t,Λ and A t,Λ are respectively transformed into the following stationary circle and annulus E Λ and A Λ :

E Λ = (y, η) : y 2 + η 2 ≤ Λ 2 , A Λ = (y, η) : Λ 2 ≤ y 2 + η 2 ≤ 4Λ 2 .
The Dirichlet boundary condition that we enforce on ∂E t,Λ pushes us to make the following change of unknown functions:

w + (t, x, θ) = α + v + (t, y, η) , w -(t, x, θ) = α -v -(t, y, η) .
We see directly that v + and v -solve respectively, for t ∈ [0, T ],

   v + t -c 1 (y, t)v + y -c 2 (t)v + η = d(η, t)v + yy + v + ηη + f r (α + v + ), (y, η) ∈ E Λ , v + (t, y, η) = 0, (y, η) ∈ ∂E Λ . (3.1)
         v - t -c 1 (y, t)v - y -c 2 (t)v - η = d(η, t)v - yy + v - ηη -f r (α -v -), (y, η) ∈ A Λ , v -(t, y, η) = 0, (y, η) ∈ ∂E Λ , v -(t, y, η) = α, (y, η) ∈ ∂A Λ \∂E Λ . (3.2)
where we have defined

c 1 (y, t) := ẊT √ Θ T + y ΘT 2Θ T , c 2 (t) := ΘT , d(η, t) := 1 + η Θ T .
The maximum principle assures us that 0 ≤ v -≤ α holds everywhere in time and space as long as the initial conditions of v -satisfy the same bound. Similarly, The maximum principle assures us that 0 ≤ v + ≤ rα holds everywhere in time and space as long as the initial conditions of v + satisfy the same bound.

We now need to fix the initial conditions. It is easy to see that the non-constant coefficients in the equation vanish as λ tends to infinity. Hence, it behooves us to look at the steady state solutions to (3.1) and (3.2) where the dependence on λ disappears.

We define for i = 1, 2, corresponding to each part of the trajectory,

c ∞ (i) := c ∞ 1,(i) , c ∞ 2,(i) .
For clarity, we emphasize that c ∞ (1) corresponds to the limit of the functions c 1 and c 2 on the interval [0, T /2] and c ∞

(2) corresponds to the limit of the functions c 1 and c 2 on the interval [T /2, T ].

For any Λ ∈ R + , we thus look for solutions ϕ ±

(i) of            -∆ϕ - (i) -c ∞ (i) • ∇ϕ - (i) = -f r (i) α -ϕ - (i) , (y, η) ∈ A Λ , ϕ - (i) (y, η) = 0, (y, η) ∈ ∂E Λ , ϕ - (i) (y, η) = α, (y, η) ∈ ∂A Λ \∂E Λ . (3.3) and    -∆ϕ + (i) -c ∞ (i) • ∇ϕ + (i) = f r (i) α + ϕ + (i) , (y, η) ∈ E Λ , ϕ + (i) (y, η) = 0, (y, η) ∈ ∂E Λ , (3.4)
where we choose r (1) = 1 and we leave r (2) ∈ (max{m, 2α}, 1) to be determined. We now state the lemma regarding these steady solutions that we require.

Lemma 3.1. There exists positive constants c 0 and Λ 0 such that if c ∞ (i) ≤ c 0 and Λ > Λ 0 , then there exist non-trivial solutions ϕ - (i) and ϕ + (i) to (3.3) and (3.4), respectively. Moreover, if ϕ + (i) is any steady solution to (3.4), then for any > 0, there exist Λ and L ≤ Λ /2, both depending only on such that if

Λ > Λ then ϕ + (i) (y, η) ≥ r (i) -α - (3.5) whenever dist((y, η), ∂E Λ ) > L .
We postpone the proof of this lemma to section 4. We now complete the definition of the function w by setting the initial conditions

v - (1) (0, •) = ϕ - (1) , v + (1) (0, •) = ϕ + (1) , and v - (2) (T /2, •) = ϕ - (2) , v + (2) (0, •) = ϕ + (2) , With v ±
(i) defined, we may now state a proposition guaranteeing that w will satisfy the condition (1.8) ensuring that it is a sub-solution. This proposition will be proved at the end of this section. Proposition 3.2. There exists positive constants Λ 0 , T 0 , and c 0 , such that if Λ ≥ Λ 0 , T ≥ T 0 , and c ≤ c 0 , then v + and v -, defined above, satisfy

|∂ n v + (i) | ≥ |∂ n v - (i) | (3.6) on E Λ .
With the initial conditions we have fixed for v ± , it is clear that as Λ tends to infinity, v ± should tend to ϕ ± uniformly in the C 1 topology. The following lemma ensures that this convergence holds and is uniform in T . Lemma 3.3. Fix > 0. There exists Λ and λ such that if Λ ≥ Λ and λ ≥ λ then

v ± (i) (t, •) -ϕ ± (i) (•) C 1 ≤ .
Moreover, this convergence holds uniformly in t.

We postpone the proof of this Lemma until the dedicated section 5. We now show how the combination of lemma 3.1 and lemma 3.3 finally implies proposition 2.1. After changing variables it is equivalent to proving the same result for v + . This is exactly now that we are going to travel along the trajectory (X T , Θ T ) and play with parameters.

We start with the first part of the trajectory. Fix a positive constant < (1 -m)/3. The discussion above implies that we may find Λ 0 , T 0 and δ 0 such that if Λ ≥ Λ 0 , T ≥ T 0 , and c ≤ c 0 , then we have that v + gives a subsolution to u such that

u T 2 , • ≥ α + v + T 2 , • ≥ ϕ + (1) (•) -ε + α ≥ 1 -2ε (3.7)
on the set E Λ/2 . This concludes our sliding over the first part of the trajectory. First, by increasing Λ 0 if necessary, we may find r (2) := r ≥ + m in the interval (2α, 1 -2ε). By the maximum principle any solution to (1.4) will be less than or equal to r.

With the choice Λ = Λ/2, we find that, by increasing T 0 and Λ 0 , if necessary, if T ≥ T 0 and Λ ≥ Λ 0 , then w is a sub-solution to u such that

u(T, •) ≥ w + (2) (T, •) = v + (2) (T, •) + α ≥ ϕ + (2) (•) -ε + α ≥ m.
on the set E Λ 2

. We point out that (3.7) is crucial to guarantee that w is a sub-solution to u on the interval [T /2, T ].

We finish this section by reducing proposition 3.2 to the following lemma.

Lemma 3.4. For Λ sufficiently large,

|∂ n ϕ + (i) (Λv)| > |∂ n ϕ - (i) (Λv)| (3.8) holds uniformly in v ∈ S 1 (R 2 ).
We postpone the proof of this lemma to section 4 below.

Proof of proposition 3.2. Fix any unit vector v ∈ S 1 (R 2 ) and look at any point Λv on the boundary of E Λ . Then (3.6) follows from lemma 3.3 and lemma 3.4.

Lemmas about the steady solutions

In this section, we prove lemma 3.1 and lemma 3.4, which give the desired behavior for the steady states ϕ ± (i) .

Proof of lemma 3.1. The existence of ϕ ± is an easy result of a fixed point argument using the principle Dirichlet eigenfunction as a sub-solution to obtain non-triviality. As such, we omit it here. We now show the lower bound for ϕ + . We first show a preliminary lower bound for ϕ + that we will bootstrap to the lower bound (2.1). Let L be a constant yet to be determined. Let (y L , η L ) be the location of the maximum of

e -1 2 c ∞ •(y,η) cos πy 2L cos πη 2L .
Now fix any point (y 0 , η 0 ) to be determined and define

ψ(y, η) = e -1 2 c ∞ •(y+y L -y 0 ,η+η L -η 0 ) cos π 2L (y + y L -y 0 ) cos π 2L (η + η L -η 0 ) , (4.1)
Notice that, by construction, the maximum of ψ occurs at (y 0 , η 0 ). Define

ψ = γ ψ(y 0 , η 0 ) ψ, (4.2)
where γ is a small parameter to be chosen later. It is easy to check that

-∆ψ -c ∞ • ∇ψ -f r (α + ψ) ≤ c ∞ 2 4 + π 2 2L 2 - c r α(r -α) 2 ψ ≤ 0, on the set C = {(y, η) | |y -(y 0 -y L )| = L or |η -(η 0 -η L )| = L} .
The first inequality follows by a simple computation and by choosing γ sufficiently small, depending only on r and f , such that c r α(rα)

2 u ≤ f r (α + u) for all u ∈ [0, γ].
The last inequality follows by decreasing c 0 and increasing L, if necessary. Hence, ψ is sub-solution to ϕ + which satisfies Dirichlet boundary conditions on the square C. This implies that as long as C ⊂ E Λ , we have ψ ≤ ϕ + on C. We point out that a weak condition for

C ⊂ E Λ is that dist((y 0 , η 0 ), ∂E Λ ) ≥ 2L + |y L | + |η L |. (4.3)
This, in turn, implies that ϕ + (y 0 , η 0 ) ≥ γ, and this may be valid as soon as (y 0 , η 0 ) satisfy (4.3). Now we claim that ϕ + tends to rα locally uniformly in the set E Λ-L as Λ and L tend to ∞. To prove this, we argue by contradiction. If this is not true, then we may find a positive > 0, sequences L k ≤ Λ k /2 tending to infinity, and a sequence of points (y

k , η k ) ∈ E Λ k -L k such that ϕ + (y k , η k ) ≤ r -α -.
Let ϕ + k (y, η) = ϕ + (y+y k , η+η k ). Elliptic regularity ensures that ϕ + k converges locally uniformly in C 2 to a function ϕ + ∞ . By assumption, we have that ϕ + ∞ (0, 0) ≤ rα -. By our work above, we have that ϕ + ∞ ≥ γ everywhere. By the maximum principle, we have that ϕ + ∞ ≤ rα everywhere. Finally, we have that ϕ + ∞ satisfies (3.4). We may assume that (0, 0) is a positive global local minimum for ϕ ∞ since, if not, we may simply re-center the equation again, taking limits if necessary. It is apparent, however, that a solution to (3.4) can only have a minimum at the zeros of f r (α + •), which are 0 and rα. Thus, we arrive at a contradiction. Hence, it must be that ϕ + tends uniformly to rα, finishing the proof of the claim.

We are now in a position to prove lemma 3.4, showing the ordering of the normal derivatives for the steady states ϕ ± (i) .

Proof of lemma 3.4. We assume the contrary. Suppose there is a sequence of vk and Λ k , with Λ k tending to infinity, such that

∂ n ϕ + (Λ k vk ) ≤ ∂ n ϕ -(Λ k vk ) + o k→+∞ (1), (4.4)
By the compactness of the unit circle, vk converges (up to a subsequence) to some unit vector v∞ . Define (y k , η k ) = Λ k vk . We center the solutions on it. In other words, we define

z - k (y, η) = ϕ -(y k + y, η k + η), z + k (y, η) = ϕ + (y k -y, η k -η). The new function z + k satisfies    -∆z + k + c ∞ • ∇z + k = f r (α + z + k ), (y -y k , η -η k ) ∈ E Λ k , z + k (y, η) = 0, (y -y k , η -η k ) ∈ ∂E Λ k , (4.5)
We point out that the condition

(y -y k , η -η k ) ∈ E Λ k ⇐⇒ (y, η) 2 -2(y, η) • (Λ k vk ) + Λ 2 k ≤ Λ 2 k , ⇐⇒ (y, η) 2 ≤ 2Λ k (y, η) • vk . The new function z - k satisfies          -∆z - k -c ∞ • ∇z - k = -f r (α -z - k ), (y + y k , η + η k ) ∈ A Λ k , z - k (y, η) = 0, (y + y k , η + η k ) ∈ ∂E Λ k , z - k (y, η) = α, (y + y k , η + η k ) ∈ ∂A Λ k \∂E Λ k . (4.6)
We point out that the condition

(y + y k , η + η k ) ∈ A Λ k ⇐⇒ Λ 2 k ≤ (y, η) 2 + 2(y, η) • (Λ k vk ) + Λ 2 k ≤ 4Λ 2 k , ⇐⇒ 0 ≤ (y, η) 2 + 2Λ k (y, η) • vk ≤ 3Λ 2 k ,
By elliptic regularity theory, we may take local uniform limits to obtain z -and z + as the limits of z - k and z + k , respectively, which solve

   -∆z + + c ∞ • ∇z + = f r (α + z + ), (y, η) ∈ P, z + (y, η) = 0, (y, η) ∈ ∂P, (4.7)
and

   -∆z --c ∞ • ∇z -= -f r (α -z -), (y, η) ∈ P, z -(y, η) = 0, (y, η) ∈ ∂P, (4.8)
where we have defined the half-plane

P def = (y , η ) : (y , η ) • v∞ > 0 .
Due to lemma 3.1, we conclude the following about z ± . In addition, z ± (y, η) is uniformly positive away from this boundary. Furthermore, z + (y, η) converges to rα and z -(y, η) converges to α as the distance between (y, η) and ∂P tends to infinity.

Finally, we point out that (4.4) implies that ∂ n z + (0, 0) ≤ ∂ n z -(0, 0), where we have removed the absolute values since the sign of both quantities is clearly positive.

This contradicts the following lemma, which we prove in the following section as it requires the results from lemma 5.1.

Lemma 4.1. There exists

δ > 0 such that if c ∞ ≤ δ then ∂ n z + (0, 0) > ∂ n z -(0, 0).
Since we have reached a contradiction, it must be that (3.8) holds for all Λ sufficiently large. This finishes the proof of lemma 3.4.

5 Convergence to the steady states: Proof of lemma 3.3

Before we prove lemma 3.3, we prove a few useful lemmas. The first of these is the uniqueness of solutions to the related steady half-plane problem. We require this to prove lemma 4.1, stated above.

We also note that lemma 5.1 with c ∞ = 0 is well-known. We include the proof of the case with constant c ∞ , which is virtually unchanged, for completeness.

Lemma 5.1. Let (y b , η b ) ∈ R 2 and v ∈ S 1 (R 2 ). Define P = {(y, η) : 0 ≤ (y + y b , η + η b ) • v}.
For any choice of c ∞ , there is at most one solution to

   -∆z + + c ∞ • ∇z + = f r (α + z + ), (y, η) ∈ P, z + (y, η) = 0, (y, η) ∈ ∂P, (5.1)
which satisfies Dirichlet boundary conditions on ∂P and which tends uniformly to rα as dist((y, η), ∂P ) tends to infinity. In addition, for any choice of c ∞ , there is at most one solution to

   -∆z -+ c ∞ • ∇z -= -f r (α -z -), (y, η) ∈ P, z -(y, η) = 0, (y, η) ∈ ∂P, (5.2)
which satisfies Dirichlet boundary conditions on ∂P and which tends uniformly to α as dist((y, η), ∂P ) tends to infinity.

Proof of lemma 5.1. We show the argument for z + but the argument for z -is identical. By rotation and translation, we may assume that v = (1, 0) and that (y b , η b ) = (0, 0). We will use a sliding method [START_REF] Berestycki | On the method of moving planes and the sliding method[END_REF]. We proceed by contradiction and assume that z 1 and z 2 solve (5.1) and satisfy z 2 (y 0 , η 0 ) < z 1 (y 0 , η 0 ) for some point (y 0 , η 0 ). Fix u 0 ∈ (0, rα) large enough such that f (u) := u(α + u)(rαu)

Lastly, we need one more ingredient in order to prove lemma 3.3. This is a uniform lower bound to v ± to ensure that when taking limits in the arguments, we do not lose positivity of v ± . Lemma 5.3. There exists Λ 0 and λ 0 such that if Λ ≥ Λ 0 and λ ≥ λ 0 , then v + and v -are uniformly positive for any subset K of E and A, respectively. Moreover, the lower bound may be chosen in such a way to depend only on the distance between the subset K and the boundary ∂E and ∂A, respectively.

Proof of lemma 5.3. This follows easily by using the time dependent eigenvalues of Lemmas 5.1 and 5.2 from [START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF] along with the strategy from lemma 3.1. As such we omit the proof.

We now finally tackle the convergence result of lemma 3.3.

Proof of lemma 3.3. Before we begin, we point out that, by parabolic regularity, it is enough to show that v ± converges to ϕ ± uniformly in L ∞ . In addition, we prove the claim for v + and ϕ + , but the proof is identical for v -and ϕ -.

Fix > 0. Choose Λ large enough that we may apply lemma 5.2 to get that all non-trivial solutions to (3.4) are within /4 in the L ∞ norm. Now we proceed by contradiction. Namely, we assume that there exists y k , η k , λ k , T k and t k ∈ [0, T k ] such that λ k and T k tend to infinity and such that

|v + (t k , y k , η k ) -ϕ + (y k , η k )| ≥ . Define v k (t, •) = v + (t + t k , •).
By the compactness of E, we have that (y k , η k ) tends, along a sub-sequence if necessary, to a point (y, η). In addition, lemma 5.3 ensures that v k is uniformly positive on any compact subset of E.

There are two cases: either t k is bounded or not. If it is bounded, then the choice of initial conditions for v + along with parabolic regularity, see e.g. [START_REF] Krylov | Lectures on elliptic and parabolic equations in Hölder spaces[END_REF], implies that v k converges to ϕ + . Hence it must be that t k tends to infinity.

In this case, it is easy to check that v + k falls within 3 /4 of a steady state of (3.4). Indeed, consider the limit v + ∞ (t, y, η). First, notice that this must solve the parabolic analogue of (3.4). Second, notice that, for any t 0 > 0, w 0 (y, η) ≤ v ∞ (-t 0 , y, η) for a positive lower bound w 0 , independent of t 0 , given by lemma 5.3. Using the linearized equation for v ∞ it is clear that we may find w 0 ≤ w 0 which is a super-solution to the related elliptic problem (3.4). Letting w be the solution to the parabolic analogue of (3.4) starting from w 0 we have that w(t 0 , y, η) ≤ v ∞ (0, y, η) by the maximum principle. By our choice of initial conditions, it must be that w τ > 0 for all times and for all y and η. By the boundedness of w, it follows that w τ tends uniformly to zero and hence, to a solution of eq. (3.4). We may similarly find an upper bound for v ∞ (0, y, η) in terms of a non-trivial solution of eq. (3.4). Hence, choosing t 0 large enough and applyinglemma 5.2 with the choice = /8, we have that v ∞ is within /4 of a steady state.

Figure 1 . 1 :

 11 Figure 1.1: Schematic plot of the flying saucer like sub-solution w that is slid along the trajectory. The bolded green line denotes the set A t , the red dotted one the set E t .

Figure 1 . 2 :

 12 Figure 1.2: Numerical simulations of the Cauchy problem of equation (1.1) at a fixed time, in the phase space R × Θ. From top left to bottom right: the initial data, t = 30, t = 60, t = 90, t = 120.One can track the accelerated behavior, for example on the space axis. The invasion at the back in the θ-direction is expected to be linear in time. This pattern is very similar as for the monostable cane toads equation, see also[START_REF] Berestycki | Existence of self-accelerating fronts for a non-local reaction-diffusion equations[END_REF][START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF][START_REF] Bouin | Super-linear spreading in local and non-local cane toads equations[END_REF].

Figure 1 . 3 :

 13 Figure 1.3: Numerical simulations of the Cauchy problem of equation (1.9) at a fixed time, in the phase space R × Θ. From top left to bottom right: the initial data, t = 30, t = 60, t = 90, t = 120.Here we see the distinctive "up-and-over" behavior that leads to acceleration for the non-local, monostable model. One can track the accelerated behavior, for example on the space axis.

Acknowledgments

EB is very grateful to Cambridge University for its hospitality during the second semester of the academic year 2015-2016. EB and CH acknowledge the support of the ERC Grant MATKIT (ERC-Applying lemma 5.2 again implies that v + and ϕ + are within /2, using again our choice of Λ above. However, this contradicts our assumption that |v + (0, y, η) -ϕ + (y, η)| ≥ , finishing the proof.

is monotonically decreasing for u ≥ u 0 . Since z 1 and z 2 tend uniformly to rα as they tends to infinity, we may find L such that if y > L, then z i (τ, y, η) ≥ u 0 for i = 1, 2.

We define z τ 2 (y, η) := z 2 (y + τ, η). Since z 2 converges locally uniformly to rα as y tends to infinity, we may find τ 0 sufficiently large that z τ 0 2 > z 1 on the set P L = {(y, η) ∈ P : 0 < y < L} .

We claim that z τ 0 2 ≥ z 1 on P . If not, let ψ = z τ 0 2 -z 1 and let G = {ψ < 0}. It is easy to see that G ⊂ P \ P L , by construction. Also, it is easy to check that ψ satisfies

on the set G. Since ψ is zero on ∂G and since f ≤ 0 by our assumptions on L and u 0 , the maximum principle ensures that ψ ≡ 0 on G, which is a contradiction. Hence, we have that z τ 0 2 ≥ z 1 on P . Let τ = inf{τ : z τ 2 ≥ z 1 on P }, and notice that τ ≤ τ 0 . By continuity, it follows that z τ 2 ≥ z 1 . We claim that τ = 0, and we proceed by supposing the opposite, that is τ > 0. There are two cases: either inf P L z τ 2 -z 1 > 0, or inf P L z τ 2 -z 1 = 0. We handle the former case first. In this case, elliptic regularity ensures that z 1 , z 2 are uniformly Lipschitz continuous, and that we may then decrease τ to a positive constant τ < τ and preserve the fact that z τ 2 ≥ z 1 on P L . Arguing as above, this implies that z τ 2 ≥ z 1 on P . This is contradicts the definition of τ . Now we handle the latter case. Assume that inf P L z τ 1 -z 2 = 0. Then there exists a sequence of points (y k , η k ) ∈ P L such that z τ 1 (y k , η k ) -z 2 (y k , η k ) tends to zero. Defining

By elliptic regularity, we may, taking a subsequence if necessary, find ψ i which are the local limits of ψ i,k as k tends to infinity. In addition, y k converges to some y since (y k , η k ) ∈ P L implies that y k is bounded. Hence we have that ψ i satisfy (5.1) on P L and satisfy that

By the maximum principle, this implies that ψ 1 ≡ ψ 2 on P L . On the other hand, since τ > 0, we have that, ψ 2 (y, η) is positive for y = -y while ψ 1 (y, η) is zero for y = -y. This is a contradiction. Hence, it must hold that τ = 0. By the definition of τ , this implies that z 2 = z τ =0 2 ≥ z 1 . This contradicts our original assumption that z 2 (y 0 , η 0 ) < z 1 (y 0 , η 0 ). This finishes our proof.

With lemma 5.1 in hand, we may now conclude the proof lemma 4.1.

Proof of lemma 4.1. Suppose that this is not true. Then there exist a sequence δ k tending to 0 as k goes to infinity and a sequence c ∞ k such that

and functions z + k and z - k solving (4.7) and (4.8) respectively which satisfy

By elliptic regularity theory, we may take local uniform limits to obtain that z + k and z - k converge to z + ∞ and z - ∞ , respectively, which are the unique solutions to (4.7) and (4.8) with c ∞ = 0. Rotating if necessary, we may now assume that v = (1, 0). Using lemma 5.1, it is easy to see that z + ∞ and z - ∞ depend only on the y-variable. It is easily checked that

Hence, we may use a sliding technique [START_REF] Berestycki | On the method of moving planes and the sliding method[END_REF] (see also lemma 5.1 below) to show that z + ∞ > z - ∞ . Finally, the Hopf Lemma implies that ∂ n z + ∞ > ∂ n z - ∞ , which contradicts (5.3). This finishes the proof.

The following lemma addresses the issue of uniqueness for the steady problem (3.4) on the bounded ellipse E Λ . While we expect uniqueness to hold in this setting, we are unable to prove it because f r need not satisfy the condition that f r (α + u)/u that is usually used for uniqueness results of this type. As such, we settle for a type of asymptotic uniqueness, stated below. Lemma 5.2. For any > 0, there exists Λ such that if Λ ≥ Λ any two positive solutions to (3.4) must be -close in the L ∞ norm.

Proof of lemma 5.2. Suppose the opposite, that there is a positive constant , a sequence of Λ k tending to infinity, and two sequences of functions ϕ 1,k and ϕ 2,k which both solve (3.4) such that ϕ 1,k -ϕ 2,k ∞ ≥ . Then we may find a sequence of points (y k , η k ) such that ϕ 1,k (y k , η k ) ≥ ϕ 2,k (y k , η k ) + . We may select (y k , η k ) to be a local maximum of ϕ

By elliptic regularity, it follows that φi,k converge, along a subsequence, if necessary, locally uniformly to ϕ i . There are two cases to consider here: either (y k , η k ) remains a finite distance from the boundary of B Λ k or it does not. The latter cannot happen due to lemma 3.1. Hence it follows that ϕ i solve (3.4) on a half-space P = {(y, η) : (y + y ∞ , η + η ∞ ) • v > 0} for some unit vector v and some translation (y ∞ , η ∞ ), with Dirichlet boundary conditions on ∂P . From lemma 5.1, we know that (3.4) has a unique positive solution on P . From this, it follows that ϕ 1 = ϕ 2 . This contradicts the fact that ϕ 1 (0, 0) ≥ ϕ 2 (0, 0) + . This contradiction finishes the proof of the claim.


Ivan Giorgio
email: ivan.giorgio@uniroma1.it

Alessandro Della Corte
email: alessandro.dellacorte.memocs@gmail.com

Francesco Dell'isola
email: francesco.dellisola@uniroma1.it

David J Steigmann
email: dsteigmann@berkeley.edu

Buckling modes in pantographic lattices

Keywords: Pantographic structures, Buckling, Second gradient models, Geodesic bending, Non-linear elasticity Structures pantographiques, Flambage, milieux continus du second gradient, flexion géodésique, Élasticité non

We study buckling patterns in pantographic sheets, regarded as two-dimensional continua consisting of lattices of continuously distributed fibers. The fibers are modelled as beams endowed with elastic resistance to stretching, shearing, bending and twist. Included in the theory is a non-standard elasticity due to geodesic bending of the fibers relative to the lattice surface.

Introduction

We have recently developed a two-dimensional continuum theory of elastic surfaces [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] to describe the deformations of pantographic lattices. The lattices are composed of intersecting fibers, or rods, that form two orthogonal curvilinear families in a reference plane. Every material point of the continuum is a point of intersection of two fibers -one from each family. In a general three-dimensional deformation of the sheet, the convected fibers of the pantograph flex to conform to the evolving surface in three-space. This is modelled by assigning strain energy to the normal curvatures and twists of the fibers. In contrast to conventional plate theories, the present model also accounts for the geodesic bending of the fibers. This mode of deformation occurs when, for example, the lattice is bent while deforming in a fixed plane. In the continuum theory, geodesic bending manifests itself as a strain-gradient e↵ect. The lattice substructure thus confers strain-gradient elasticity to the continuum [START_REF] Toupin | Theories of elasticity with couple stress[END_REF][START_REF] Mindlin | E↵ects of couple-stresses in linear elasticity[END_REF][START_REF] Koiter | Couple-stresses in the theory of elasticity[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics, part 2: microstructure[END_REF][START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF][START_REF] Grillo | Growth, mass transfer, and remodeling in fiber-reinforced, multi-constituent materials[END_REF][START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity and damage[END_REF][START_REF] Fosdick | A generalized continuum theory with internal corner and surface contact interactions[END_REF][START_REF] Camar-Eddine | Determination of the closure of the set of elasticity functionals[END_REF]. Strain gradient models, and in general higher gradient ones, are today a very topical research subject, and many relevant results are available in both theoretical [START_REF] Dell'isola | Cauchy Tetrahedron Argument Applied to Higher Contact Interactions[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Placidi | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coe cients[END_REF] and applicative [START_REF] Dell'isola | Elastic pantographic 2D lattices: a numerical analysis on the static response and wave propagation[END_REF][START_REF] Madeo | Wave propagation in relaxed micromorphic continua: modeling metamaterials with frequency band-gaps[END_REF][START_REF] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF][START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic second-gradient continuum model[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF][START_REF] Yang | Higher-Order Continuum Theory Applied to Fracture Simulation of Nanoscale Intergranular Glassy Film[END_REF][START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF][START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF] directions. These models naturally arise as homogenized limits of microstructured plates (see e.g. [START_REF] Altenbach | On the linear theory of micropolar plates[END_REF][START_REF] Altenbach | Acceleration waves and ellipticity in thermoelastic micropolar media[END_REF][START_REF] Altenbach | Chapter in Generalized Continua from the Theory to Engineering Applications 541 of the series[END_REF][START_REF] Harrison | Modelling the forming mechanics of engineering fabrics using a mutually constrained pantographic beam and membrane mesh[END_REF][START_REF] Eremeyev | On free oscillations of an elastic solids with ordered arrays of nano-sized objects[END_REF][START_REF] Boubaker | Discrete models of woven structures. Macroscopic approach[END_REF][START_REF] Reis | Equivalent mechanical properties of auxetic lattices from discrete homogenization[END_REF]).

A simpler version of the theory was used in [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF] to model deformations of a lattice of initially straight rods in which the shear strain segregates into coexistent phases separated by narrow transition layers in which geodesic bending predominates. The present study extends and generalizes the results of [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF], by using the framework outlined in [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]. Specifically, we consider herein fibers which are not straight lines in the reference configuration (sinusoids, and logarithmic spirals). The interest in these geometries lies in their peculiar mechanical e↵ects, induced by arch-like response to external stresses and other e↵ects already mentioned in the literature (see below, Section 3). Moreover, we performed here a numerical study that exhibits interesting and heretofore unknown two-and three-dimensional buckling patterns in presence of both rectilinear and curvilinear fibers. In Section 2, we recall the basic kinematics of the model, following [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]. In Section 4, we specify the basic constitutive structure in terms of a strain-energy function, measured per unit area of the initial plane. This is applied to initially curved lattices conforming to two types of coordinate systems: logarithmic spirals and sinusoidal coordinates. The relevant details are outlined in Section 5.

Kinematics of pantographic lattices

We adopt the model developed in [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF][START_REF]dell'Isola F and Steigmann DJ, A two-dimensional gradient-elasticity theory for woven fabrics[END_REF], in which the pantographic sheet is regarded as an elastic surface endowed with suitable kinematic descriptors and an associated constitutive response in the form of an areal strain-energy density. The latter is assumed to depend on the stretches of the fibers and on their curvatures and twists. Further, as tangential stretch gradients appear in the constitutive equations for one-dimensional models of thin fibers that account for finite-thickness e↵ects [START_REF] Coleman | Necking and drawing in polymeric fibers under tension[END_REF][START_REF] Coleman | On the Rheology of Cold Drawing: I. Elastic Materials[END_REF], we also include a constitutive sensitivity to the tangential derivatives of the stretches. In addition, we incorporate a nonstandard constitutive sensitivity to the cross derivatives of the fiber stretches; these record the variation in the stretch of a fiber family as the (initially) orthogonal family is traversed. This e↵ect is discussed more fully in [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF]. To make the paper reasonably self contained, we devote the remainder of the present section to a brief summary of the theory proposed in [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF].

The deformation of the surface is described by a map r(x): ⌦ 7 ! E 3 , where ⌦ ⇢ E 2 is a connected plane region in the 2D Euclidean space E 2 and E 3 is Euclidean 3-space. Let L, M 2 E 2 be orthogonal families of unit vectors defining the fiber trajectories prior to deformation. The fiber stretches { , µ} and fiber trajectories {l, m} induced by the deformation are l = (rr)L and µm = (rr)M,

where r(•) is the (two-dimensional) gradient with respect to x. Here l and m are unit vectors spanning the deformed tangent plane at the material point x 2 ⌦. We use them to define the fiber shear angle by sin = l • m.

(

Because {L, M} is an orthonormal basis at any point of ⌦, we may use [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] to derive the representation

rr = l ⌦ L + µm ⌦ M (3)
and thus obtain the Cauchy-Green deformation tensor, C = (rr) t rr, in the form

C = 2 L ⌦ L + µ 2 M ⌦ M + µ sin (L ⌦ M + M ⌦ L). (4)
Further, Eqs. [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] provide

Jn = (rr)L ⇥ (rr)M, ( 5
)
where n is the unit normal to the deformed surface, and

J = µ |cos | (6)
is the local areal dilation induced by the deformation. Accordingly,

l ⇥ m = |cos | n. (7)
To describe the fiber curvatures and twist we require the second gradient rrr of the deformation. This may be represented in the form [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] rrr

= r |↵ ⌦ e ↵ ⌦ e , (8)
where

r |↵ = S ↵ r , + b ↵ n. (9)
Here (•) |↵ is the second covariant derivative with respect to an arbitrary convected curvilinear coordinate system on the reference plane ⌦, the b ↵ are the coe cients of the second fundamental form (the covariant curvature tensor) of the deformed surface, and the S ↵ are terms involving the strain and strain gradient. The rather complicated explicit expressions are [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]

S ↵ r , = ( L ↵ ) | l + ( 2 ⌘ l L ↵ L + µ l L ↵ M )p +(µM ↵ ) | m + (µ 2 ⌘ m M ↵ M + µ m M ↵ L )q, ( 10
)
and

b ↵ = K L L ↵ L + K M M ↵ M + T (L ↵ M + M ↵ L ), (11)
where ⌘ l and ⌘ m are the geodesic curvatures of the deformed fibers; l and m are the so-called Tchebychev curvatures;

p = n ⇥ l and q = n ⇥ m (12)
define the orthogonal trajectories of the fibers on the deformed surface; and,

K L = 2  l , K M = µ 2  m and T = µ⌧, (13)
where  l and  m are the normal curvatures of the deformed fibers and ⌧ measures the twist of the deformed surface. These are non-zero if the deformation is such as to generate a curvature of the surface in three-space. Accordingly, they describe those parts of the fiber curvatures that can be attributed to surface flexure, whereas the geodesic curvatures represent the components of fiber curvatures in the tangent planes of the deformed surface. The latter arise from the surface strain and the strain gradient; the explicit expressions are [START_REF] Steigmann | Equilibrium of elastic nets[END_REF] J⌘

l = div[(µ sin )L M] and J⌘ m = div[µL ( sin )M], ( 14
)
where div is the divergence on ⌦, whereas the Tchebychev curvatures are given by [START_REF] Steigmann | Equilibrium of elastic nets[END_REF] J l = J⌘ m + M • r(sin ) and J m = J⌘ l µL • r(sin ).

Later, we will have occasion to exploit the fiber decomposition [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]

rrr = (L ↵ L r |↵ ) ⌦ L ⌦ L + (M ↵ M r |↵ ) ⌦ M ⌦ M + (L ↵ M r |↵ ) ⌦ (L ⌦ M + M ⌦ L), (16)
which follows from [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity and damage[END_REF], with

e ↵ = L ↵ L + M ↵ M, (17)
and from the symmetry of r |↵ with respect to interchange of the subscripts.

In the reference configuration, i.e., in the absence of deformation, we have rr = 1, the two-dimensional identity on ⌦. It follows that rrr vanishes in this configuration and hence that S ↵ and b ↵ also vanish. Moreover, ⌘ l and ⌘ m reduce to

⌘ L = divM and ⌘ M = divL, (18)
respectively. These are the geodesic curvatures of the fibers on the reference plane. They account for the initial curvatures of the fibers, which are assumed to be initially untwisted. The Tchebychev curvatures l and m reduce, respectively, to ⌘ M and ⌘ L in the reference configuration. Because the S ↵ and b ↵ are Galilean invariant [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] and vanish in the reference configuration, they furnish natural strain measures in the present theory. The vanishing of b ↵ in the reference configuration follows trivially from the fact that the latter is a region on a plane. We have proved that the S ↵ also vanish there, but nevertheless pause to verify this directly -as a check on (10)-before proceeding.

In the reference configuration (10) reduces to

S ↵ e = (L ↵| ⌘ M M ↵ M ⌘ L M ↵ L )L + (M ↵| + ⌘ L L ↵ L + ⌘ M L ↵ M )M. (19)
Using

1 = L ⌦ L + M ⌦ M we can write rL = (rL)L ⌦ L + (rL)M ⌦ M, (20)
and then conclude that divL = tr(rL

) = L • (rL)L + M • (rL)M. ( 21
)
The first term on the right is (rL) t L • L, in which (rL) t L = 0 because L(x) is a unit-vector field. It then follows from [START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF] 2 that ⌘ M = M • (rL)M. In the same way, we use [START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF]

2 to prove that ⌘ L = L • (rM)L = M • (rL)L,
where the second equality is a consequence of (rL) t M = (rM) t L, which in turn follows by computing the gradient of L • M(= 0). These results yield

(rL)L = [L • (rL)L]L + [M • (rL)L]M = ⌘ L M (22)
and

(rL)M = [L • (rL)M]L + [M • (rL)M]M = ⌘ M M, ( 23
) yielding rL = ⌘ L M ⌦ L + ⌘ M M ⌦ M, (24)
and implying that the first parenthesis in [START_REF] Yang | Higher-Order Continuum Theory Applied to Fracture Simulation of Nanoscale Intergranular Glassy Film[END_REF] vanishes identically. The second parenthetical term is disposed of similarly, and our claim is proved. The present model di↵ers from that presented in [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF] in that {L, M} are allowed here to be functions on ⌦.

Lattice arrangement

To specify the distribution of curvilinear lattice fibers on the reference plane, we identify any particular fiber, with unit-tangent field M, say, as the locus of points on which a scalar field '(x) assumes a constant value. The family of all such fibers is then obtain by varying the value of the constant. The unit-tangent fields of the fibers are given by

L(x) = r'/ |r'| and M(x) = k ⇥ L(x), (25)
where k is the unit normal to the reference plane (i.e., the value of the surface-normal field n in the reference configuration).

Technological applications for di↵erent geometries of the fibers

In the opinion of the authors, it has been proven by countless examples that byproducts of what looked like pure 'amusement' paid o↵ extremely well in applications. More often than not, pursuing the aim of Jean Dieudonné, i.e. studying a subject only pour l'honneur de l'esprit humain, has been the best policy in order to concretely improve our lives. This mere consideration could serve as a justification for studying a certain scientific subject, as long as its theoretical foundation is clear and sound enough. However, in our specific case, direct reference to existing or conceivable application can be mentioned concerning particular geometries for the fibers.

In recent years, spiral woven fabrics are increasingly employed and examples of these artifacts can be found in many fields concerning e.g. biomedical, aerospace and automotive applications, including spinning components like flywheels or turbine bliscs [START_REF] Clarke | Net Shape Woven Fabrics-2D and 3D[END_REF]; the reason behind the interest for spiral geometry lies in the possibility of combining circumferential fiber orientation for strength and radial fiber to manage interlaminar stresses [START_REF] Clarke | Net Shape Woven Fabrics-2D and 3D[END_REF]. However, theoretical coverage of this type of geometry for the fibers is still not very developed. This seems indeed one of the cases in which technological potential is being exploited even before a fully comprehensive understanding is reached.

Sinusoidal fibers furnish continuous-curvature approximations to the curved lattices adopted in the design of flexible electronics for health-monitoring purposes [START_REF] Kim | Epidermal electronics[END_REF][START_REF] Dell'isola | Synthesis of complex structures. Designing micro-structure to deliver targeted macro-scale response[END_REF]. Such devices are applied to skin tissue and are conceived to stick perfectly on it, possibly exhibiting complex large-deformation behavior, including out-of-plane buckling. The aim of the following sections is to provide some insight in the possible behavior of the two relatively complex geometries before mentioned.

Logarithmic spirals

These woven fabrics are associated with the orthogonal curvilinear coordinates

' = log(r/r 0 ) k✓ and = k log(r/r 0 ) + ✓, (26)
where r 0 and k are fixed parameters and {r, ✓} are plane polar coordinates. We compute

r' = 1 r (e r ke ✓ ), (27)
where {e r , e ✓ } is the standard orthonormal basis associated with the polar coordinates. Thus,

L(x) = 1 p 1 + k 2 (e r ke ✓ ) and M(x) = 1 p 1 + k 2 (e ✓ + ke r ). ( 28
)
This is expressed directly in terms of x(= x ↵ i ↵ ) by substituting

r = |x| , e r = r 1 x and e ✓ = k ⇥ e r . (29)

Oscillatory lattice

The fiber net which we consider here is described by

x 2 = A sin(!x 1 ) + ', (30)
in which the amplitude and frequency parameters A and ! are taken to be fixed. This immediately determines '(x), and hence L(x) and M(x) via Eq. ( 25). We omit the details of the trivial calculation.

Strain-energy function

All constitutive information is contained in an areal strain-energy density W that depends on the first and second gradients of the deformation. To non-dimensionalize the variables appearing in this function it is necessary to introduce a local length scale. This could be the thickness of the actual sheet, the characteristic spacing of the internal points of intersection of rods of the actual pantographic lattice, or the widths of the constituent fibers. If any of these is used as the length scale, and if the couple stresses and bending/twisting moments vanish when the lattice is undeformed, then in practice the non-dimensionalized second gradient of the deformation is su ciently small as to justify a quadratic dependence in the strain-energy function. A simple strain-energy function of this type, incorporating the curvilinear orthotropic symmetry associated with the initial fiber geometry, is [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF]

W = w(✏ L , ✏ M , J) + 1 2 A L L ↵ L r |↵ 2 + 1 2 A M M ↵ M r |↵ 2 + 1 2 A L ↵ M r |↵ 2 , (31)
where the coe cients A L,M, are constants; this form is assumed for the sake of definiteness and tractability. Other forms are, of course, possible. We take A L,M, and to be strictly positive and observe, from the fiber decomposition [START_REF] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF], that the part of the energy depending on the second gradient rrr is then non-negative, vanishing if and only if all r |↵ vanish simultaneously. It is thus a convex function of rrr, ensuring the existence of energy-minimizing deformations in conservative boundary-value problems [START_REF] Ball | Null Lagrangians, weak continuity, and variational problems of arbitrary order[END_REF]. Further, the existence of a solution to the weak form of the equilibrium equations for similar problems in second-gradient elasticity has been proved in [START_REF] Healey | Injective weak solutions in second-gradient nonlinear elasticity[END_REF], albeit under hypotheses slightly more stringent than those satisfied here.

The constitutive sensitivities to geodesic and normal curvatures may be understood in terms of the mechanics of the lattice substructure. Thus, fibers o↵er resistance to bending in three-space, which may be resolved into geodesic bending and surface flexure. Fiber twist is determined by the twist of the surface because of the presumed connectivity of the lattice; that is, both families of intersecting fibers are assumed to pivot about the normal to the deformed surface at their points of intersection. This constraint implies that fiber twist is controlled entirely by surface geometry and is therefore not an independent kinematic variable [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF].

For the strain-dependent function w, we adopt

w(✏ L , ✏ M , J) = 1 2 (E L ✏ 2 L + E M ✏ 2 M ) G(ln J + 1 J), (32)
where

✏ L = 1 2 ( 2 1), ✏ M = 1 2 (µ 2 1) (33)
are the extensional fiber strains and E L,M and G are positive constants. This energy does not include a term proportional to ✏ L ✏ M and therefore does not exhibit a Poisson e↵ect with respect to the fiber axes.

Poisson e↵ects are typically negligible in pantographic lattices. Finally, the term involving J penalizes fiber collapse (J ! 0) by imposing unbounded growth of the attendant energy, whereas the remaining terms are appropriate for small-to-moderate fiber strains. Our solution procedure, discussed in Section 5, is a finite-element scheme based on the weak form of the equilibrium equations derived from the proposed strain-energy function. This is not presented here, but reference may be made to [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF], where the weak and strong formulations are given explicitly. Recent applications of the finite element method and related methods to second-gradient elasticity are discussed in [START_REF] Turco | A three-dimensional B-spline boundary element[END_REF][START_REF] Fischer | On the C1 continuous discretization of non-linear gradient elasticity: a comparison of NEM and FEM based on Bernstein-Bézier patches[END_REF][START_REF] Fischer | Isogeometric analysis of 2D gradient elasticity[END_REF][START_REF] Greco | An isogeometric implicit G1 mixed finite element for Kirchho↵ space rods[END_REF][START_REF] Greco | An implicit G1 multi patch B-spline interpolation for Kirchho↵-Love space rod[END_REF][START_REF] Rudrarajua | Three-dimensional isogeometric solutions to general boundary value problems of Toupin's gradient elasticity theory at finite strains[END_REF][START_REF] Cazzani | Isogeometric analysis: a powerful numerical tool for the elastic analysis of historical masonry arches[END_REF][START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF][START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF].

We observe that in Eq. ( 31) the contravariant components L ↵ and M ↵ , and the covariant derivatives r |↵ , pertain to any convected coordinate system. Moreover, the symmetry r |↵ = r | ↵ implies that the ordering in the final term of Eq. ( 31) is irrelevant; i.e., L ↵ M r |↵ = M ↵ L r |↵ . It would be natural to identify the coordinates with the fiber parameters ' and of Section 3, but here, for the sake of tractability, we exploit the underlying coordinate invariance and instead choose them to be Cartesian on ⌦. In this case the distinction between co-and contra-variance disappears, covariant derivatives reduce to partial derivatives, and Eq. ( 31) becomes

W = w(✏ L , ✏ M , J) + 1 2 A L |L ↵ L r ,↵ | 2 + 1 2 A M |M ↵ M r ,↵ | 2 + 1 2 A |L ↵ M r ,↵ | 2 , (34)
with

J = |L ↵ M r ,↵ ⇥ r , | , ✏ L = E ↵ L ↵ L and ✏ M = E ↵ M ↵ M , (35)
where

E ↵ = 1 2 (r ,↵ • r , ↵ ), (36)
with ↵ the Kronecker delta, is the strain. Here, r ,↵ = r i,↵ i i , with i 3 = i 1 ⇥ i 2 ,; {i k } a fixed orthonormal basis associated with the Cartesian coordinates, and r ,↵ = r i,↵ i i . In the examples considered, the Cartesian components L ↵ (x) and M ↵ (x) are easily read o↵ from the formulas provided in Section 3.

Examples

Numerical simulations here shown are based on Finite Element Analysis and performed with COMSOL Multiphysics R in the framework of a weak formulation. Because of the nature of the energy (31), we employed triangular Argyris elements having C 1 continuity [START_REF] Giorgio | Isola F and Steigmann DJ, Pattern formation in the three-dimensional deformations of fibered sheets[END_REF]. The number of mesh elements is ⇡ 10 4 , which considering the 21 degrees of freedom of each element, provide a total of ⇡10 5 degrees of freedom.

In the performed simulations, we adopt the representation [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]:

L ↵ L r ,↵ = g l + K L n, M ↵ M r ,↵ = g m + K M n, L ↵ M r ,↵ = + T n (37)
where g l,m and account for geodesic bending in the tangent plane of the deformed surfaces and n is the (varying) unit normal to the surface. The variables K L,M and T are the normal curvatures of the embedded fibers and the twist induced by the bending of the surface in three-space. With this notation, in Eq. ( 31), the part of the energy depending on the second gradient rrr can be decomposed into the three following contributions which allow to emphasize the role of each term: i) the strain-gradient energy,

W SGE = 1 2 A L |g l | 2 + 1 2 A M |g m | 2 + 1 2 A | | 2 ; (38)
ii) the normal bending energy,

W NBE = 1 2 A L K 2 L + 1 2 A M K 2 M ; (39)
and iii) the twisting energy,

W T E = 1 2 A T 2 . ( 40
)
Because the normal derivative of position is not assigned, implicitly we are assigning zero double force on edges where position is assigned. Elsewhere we assign zero double force and zero traction unless otherwise specified. In the considered examples, a non-dimensional form of the problem is achieved by normalizing the strain energy with respect to a reference sti↵ness E 0 , i.e.:

W = W E 0 (41)
As a consequence of this normalization, the material parameters associated with Eq. ( 34) become:

ẼL = E L E 0 , ẼM = E M E 0 , G = G E 0 , ÃL = A L E 0 , ÃM = A M E 0 , Ã = A E 0 (42)
The lengths are normalized with respect to a reference length L 0 , whose value will be specified case by case. In what follows, quantities in non-dimensional form are denoted with a superimposed tilde.

Herein we study the buckling mode shapes of the considered pantographic fibered structures. These two-dimensional continua are highly sensitive to geometrical and mechanical imperfections. Specifically, defects that may be encountered are: some di↵erences in the directions of fiber pattern, delaminations or initial transverse deformations, variations in thickness, inclusions, eccentricity of loading and so forth. For numerical purposes, rather than perturbing the geometry, it is often easier to obtain the perturbation using an additional load or a further constraint on the displacement. This latter approach is the one followed in this paper. The heuristic method employed here for fibered pantographic sheets is su cient for the immediate goal of determining some equilibrium shapes, but of course a rigorous analysis is needed in order to successfully predict buckling modes (see e.g. [START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF][START_REF] Luongo | On the e↵ect of twist angle on nonlinear galloping of suspended cables[END_REF][START_REF] Rizzi | The e↵ects of warping on the postbuckling behaviour of thin-walled structures[END_REF][START_REF] Gabriele | On the imperfection sensitivity of thin-walled frames[END_REF][START_REF] Rizzi | Initial postbuckling behavior of thin-walled frames under mode interaction[END_REF] for specific analytical tools concerning buckling and postbuckling).

Geodesic and out-of plane buckling of a lattice with initially straight fibers

In this section, we consider a rectangle plane region ⌦ whose sides are in the ratio 1:3, with the longest of them being parallel to the bisector of the first quadrant of the reference plane. The two orthogonal families of fibers are arranged so as to form an angle of 45 degrees with the sides of this region; thus, L, M are constant unit vectors and coincide with the basis elements e 1 , e 2 associated with a Cartesian system. The shorter edge is used as the length scale for non-dimensionalization.

A first simulation is performed by imposing a shear displacement (ũ 0 = 0.1 p 2) on one of the long sides and by fixing the opposite one, while the short sides are left free. The non-dimensional material parameters employed in this case are: ẼL = ẼM = 100; G = 0.03; Ã = 0.1; ÃL = ÃM = 0.001. In this case to take imperfections into account, we impose on the displaced long side the additional constraint on the derivatives of the displacement lying on the plane of the pantographic sheet and, in particular we impose: ũ↵, n = " with ↵, = {1, 2}, " being a positive quantity equal to 2.22 ⇥ 10 16 . In Fig. 1 is shown an example of buckling occurred in the plane of the two-dimensional continuum considered, that we call "geodesic buckling" because it is related to the strain-gradient energy of Eq. ( 38) and therefore to the geodesic curvatures which appear in this term of the energy. In particular, on the left side of the figure the fiber pattern of the 2D continuous sheet is shown, while on the right side is plotted the strain-gradient energy [START_REF] Healey | Injective weak solutions in second-gradient nonlinear elasticity[END_REF] which, in this case, represents the only term in the second gradient energy, the contributions related to the normal bending and twisting being zero.

The curves which are visible in the deformed shape of the body are simply selected mono-dimensional sets of material points which, in the reference configuration, are parallel to the field lines relative to the fields L and M. Figure 2 shows an example of out-of-plane buckling. In Fig. 2a the equilibrium shape is exhibited (colors indicate the out-of-plane component of displacement). In this case, in addition to a shear displacement (ũ 0 = 0.05 p 2) we have also imposed a displacement of compression (ṽ 0 = 0.007 p 2) on the same long side and the constraint: ũ↵, n = 0, with ↵, = {1, 2}; the other sides are subject to the same constraints of the case previously considered. The non-dimensional material parameters used in this case are: ẼL = ẼM = 100; G = 0.02; Ã = 0.1; ÃL = ÃM = 0.01. To simulate imperfections, we set the condition for the derivatives of the displacement out-of-plane as: ũ3, n = 1 ⇥ 10 4 with = {1, 2}.

In this simulation, we can see how the strain-gradient energy (Fig. 2b) is localized mainly near the corners of the rectangular fibered sheet, while the normal bending energy (Fig. 2c) and the twisting energy (Fig. 2d) have greater relevance.

In the next simulations, we consider two examples with di↵erent material parameters to investigate how they a↵ect the solution. The simulations are performed by imposing a compressive displacement (ũ 0 = 0.07 p 2) on one of the long sides, by fixing the opposite one and leaving the short sides free. In addition, on the moving long side we assign the constraint: ũ↵, n = 0, with ↵, = {1, 2}. The non-dimensional material parameters employed in these cases are:

Case I. ẼL = ẼM = 100; G = 0.02; Ã = 0.1; ÃL = ÃM = 0.1;

Case II. ẼL = ẼM = 100; G = 0.02; à = 0.1; ÃL = ÃM = 0.01. To take imperfections into account, we impose on the moved long side the additional constrain on the derivatives of the displacement lying on the plane of the pantographic sheet and, in particular, we impose: ũ↵, n = 1 ⇥ 10 4 with ↵, = {1, 2} Figures 3 and4 show the out-of-plane buckling modes in the cases examined. Specifically, these figures give details of the particular fiber pattern assumed in the current configuration (Figs. 3a and4a), and of the terms of the second gradient energy above mentioned (i.e., Eqs. ( 38), ( 39) and ( 40)). Also in these simulations, we can see how the strain-gradient energy (Figs. 3b and4b) is localized mostly near the corners of the rectangular fibered sheet and is of the same level in both cases; in the second case, when ÃL = ÃM = 0.01 are smaller than à = 0.1 the twisting energy attains a maximum in the center very high if compared with the first case (Figs. 3d and4d); on the contrary, the normal bending energy achieves its maximum level in the first case (Figs. 3c and4c). As a result, the two buckling modes obtained are enough di↵erent given that the relative rigidities of the terms of the second gradient energy promote a type of deformation rather than the other one.

Figure 5 shows the plot of convergence of the total energy stored versus mesh refinement, i.e. the inverse of the mesh size h in the two cases before examined. Of course, the convergence analysis here shown is not to be intended in the sense of a convergence toward a known limit solution, as there is no chance of getting a close form solution as a benchmark even with the simplest possible boundary conditions and constitutive assumptions. However, the analysis just shown strongly supports the claim that the set of numerical solutions one has decreasing the mesh size is a subset of a Cauchy sequence, i.e. the di↵erence between any two solutions can be made arbitrarily small if the largest mesh size between the two is small enough.

Buckling of logarithmic spirals

In this section, we consider an annular two-dimensional sheet with fibers forming logarithmic spirals (Fig. 6). This fibered sheet is deformed by fixing it at the outer boundary and imposing a displacement at the inner boundary corresponding to a counterclockwise rotation of 30 degrees and a radial opening of amplitude ũ0 = 0.06. The radius of the inner circle is 0.05 and that of the outer boundary is 0.2. The logarithmic spirals in the initial configuration (see Fig. 6-left) are characterized by the parameter k = tan(⇡/4) and r0 = 0.1 (see Eq. ( 26)). The non-dimensional material parameters utilized in this case Figure 7 shows the deformation shape of the annular fibered sheet, and in particular the pattern of fibers (Fig. 7a) -colors indicate the out-of-plane component of the displacement. Figure 7 shows also the three terms of the second gradient energy which in this case are significant.

Oscillatory lattice

In this case, we consider a rectangle plane region ⌦ whose sides are in the ratio 0.75:1. The longer edge is used as the length scale for non-dimensionalization. The parameters characterizing the oscillatory fibers which appear in Eq.( 30) are assumed to be à = A/L 0 = 0.25, ! = !L 0 = 4⇡.

In a first case, we examine the standard bias extension test in which one of the longer side is fixed and on the other edge a displacement ũ0 = 0.2 is imposed in the direction orthogonal to the same side. On the remaining edges, we assign zero double force and zero traction. 8a), the equilibrium shape of the sample after the deformation and the new disposition of the net (Fig. 8b), the second gradient energy (Fig. 8c) and the measure of the shear strain relative to the initial fiber axes (Fig. 8d). The more easily detectable feature in the plot of the shear strain (Fig. 8d) is the presence of three distinct bands in which the shear strain presents a distribution almost uniform.

In the next examples, we impose a shear displacement (ṽ 0 = 0.2) and a displacement of compression (ũ 0 = 0.2) on one of the longer side being fixed the opposite one. We study two cases of buckling by assuming the following material parameters: Case a. ẼL = ẼM = 100; G = 0.2; Ã = 0.1; ÃL = ÃM = 0.01; Case b. ẼL = ẼM = 100; G = 0.2; Ã = 0.1; ÃL = ÃM = 0.001. To simulate imperfections, we impose on the moved long side the constraint: With this kind of application at hand, it is very interesting to investigate how changes in the constitutive parameters can a↵ect the out-of-plane displacement of sinusoidal lattices. Indeed, epidermal electronics strongly requires a close to perfect adherence of the lattice to the skin surface, which means that a set of constitutive assumptions that robustly resist to buckling is to be preferred over a generic one. The simulations shown in Fig. 9 and 10 compare the case with ÃM = 0.01 and the one with ÃM = 0.001. As one can see, in the latter case the average out-of-plain displacement is one order of magnitude smaller (0.0467 vs 0.00339), which means that a smaller second gradient sti↵ness ÃM of the fibers may ensure a better adherence of the lattice to the skin surface. Specifically, these figures give details of the particular fiber pattern assumed in the current configuration -colors indicate the out-of-plane component of displacement-(Figs. 9a and10a), and of the terms of the second gradient energy already considered, i.e. Eqs. [START_REF] Healey | Injective weak solutions in second-gradient nonlinear elasticity[END_REF], [START_REF] Turco | A three-dimensional B-spline boundary element[END_REF] and [START_REF] Fischer | On the C1 continuous discretization of non-linear gradient elasticity: a comparison of NEM and FEM based on Bernstein-Bézier patches[END_REF]. Analogous considerations apply in these last cases as for the straight fibers regarding the second gradient terms of the energy considered.

Conclusion

The numerical simulations performed on the basis of the theoretical model briefly described in Sections 2 and 4 provided some examples of mechanical instability. Specifically, out-of-plane buckling was observed with straight and curved fibres geometries. Moreover a more exotic in-plane buckling was also obtained in Fig. 1 concerning the simple case of straight orthogonal fibers. This buckling is linked to the fact that our continuum model is able to store energy in connection with the geodesic bending of the fibers, which in turn is modeled by means of a strain gradient theory. The amount of energy stored in this way, when compared to the standard out-of-plane flexural and torsional energies, can be non-negligible, as demonstrated in simulations shown in Figs. 2,7, 9 and 10.

Figure 1 .

 1 Figure 1. Geodesic buckling for the case of plane shear: fiber pattern (left) and the strain-gradient energy (right).

Figure 2 .

 2 Figure 2. Case of initially straight fibers with compression (ṽ 0 = 0.007 p 2) and shear (ũ 0 = 0.05 p 2) displacement imposed.

Figure 3 .

 3 Figure 3. Case I ( ÃL = ÃM = 0.1): compressive displacement imposed (ũ 0 = 0.07 p 2). are: ẼL = ẼM = 500; G = 0.1; Ã = 0.1; ÃL = ÃM = 0.01. Further, to simulate imperfections, we impose: ũ3, x |x| 1 = 1 ⇥ 10 3 with = {1, 2} on the outer boundary; ũ3, x |x| 1 = 1 ⇥ 10 3 with = {1, 2} on the inner boundary.

Figure 4 .

 4 Figure 4. Case II ( ÃL = ÃM = 0.01): compressive displacement imposed (ũ 0 = 0.07 p 2).

Figure 5 .

 5 Figure 5. Convergence plots: Case I (left); Case II (right).

Figure 6 .Figure 7 .Figure 8 . 2 Figure 8

 67828 Figure 6. Buckling example with logarithmic spirals arrangement of fibers in the case of radial displacement (ũ 0 = 0.06) and rotation (30 degrees) of the inner circle: Fiber pattern in the reference configuration (left); the shear strain relative to the initial fiber directions (right).

Figure 9 .

 9 Figure 9. Case a: Compressive (ũ 0 = 0.2) and shear (ṽ 0 = 0.2) displacement imposed ( ÃM = 0.01).

Figure 10 .

 10 Figure 10. Case b: compressive (ũ 0 = 0.2) and shear (ṽ 0 = 0.2) displacement imposed ( ÃM = 0.001).

Acknowledgements

The authors would like to thank Tomasz Lekszycki for his helpful comments and advice. The authors also thank the reviewer for having drawn their attention on the importance of explicitly showing connections between the object of their investigations and possible real-world applications.