Marie Pierre Krafft

Jean G Riess

Perfluorocubane-a tiny electron guzzler

Perfluorocubane: a tiny electron guzzler

Perfluorination gives cubane the capacity to host an extra electron in its inner structure By Marie Pierre Krafft 1 and Jean Riess 2 1 University of Strasbourg, Institut Charles Sadron (CNRS), 67034 Strasbourg, France. 2 Harangoutte Institute, 68160 Ste-Croix-aux-Mines, France. Email: krafft@unistra.fr High symmetry arouses a sense of perfection and also of perennity. In "Electron in a cube: synthesis and characterization of perfluorocubane as an electron acceptor", Sugiyama et al. [START_REF] Sugiyama | Electron in a cube: synthesis and characterization of perfluorocubane as an electron acceptor[END_REF] offer us an aesthetically appealing small molecular gem, a synthetic feat and, most remarkably, a perfect nest for hosting an added electron, with many inciting perspectives.

Eaton and Cole. succeeded in synthesizing cubane in spite of a priori improbable bonding angles [START_REF] Eaton | The cubane system[END_REF]. Yet the strained, thermodynamically unstable molecule proved kinetically inert. In the molecular world, symmetry tends indeed to impede reactivity, as no weak spot, or symmetry-allowed pathway, is available for intrusion. Replacing all the hydrogens of the utterly strained cubane substrate by much larger fluorine atoms was a further momentous challenge. Perfluorination entails indeed a drastic upheaval of the molecule's structure and electron charge density distribution. Compared to hydrogen, fluorine, with its nine protons and nine electrons, has a much larger size, higher electron density, ionization potential, electron affinity and electronegativity, and lower polarizability. Perfluorocarbons are among the most thermally stable and chemically inert classes of organic molecules ever synthesized. Perfluorination generates new properties or heightens known ones to levels often unattainable otherwise [START_REF] Krafft | Selected physicochemical aspects of poly-and perfluoroalkylated substances relevant to performance, environment and sustainability-Part one[END_REF]. Their outstanding inertness calls, however, for cogent use since their persistence and bioaccumulation characteristics can raise serious environmental issues (4). Regarding perfluorocubane, the predicted perspective of netting an electron in a Platonian carbon polyhedron, thus generating a radical anion [START_REF] Irikura | Sigma stellation: A design strategy for electron boxes[END_REF], offered a tantalizing additional challenge. The Authors selected direct radical fluorination with elemental fluorine, F2, to achieve their aim. Taming the impulsively explosive element is, however, not obvious. Mainly due to steric hindrance, the homonuclear bond in F2 is the weakest in its row in the Periodic Table . If not duly controlled, its reaction with organic material is highly exothermal, provoking profuse C-C bond breaking. Direct fluorination had been pioneered by diluting F2 in an inert gas at low temperature (LaMar process), improved by using a chlorofluorocarbon solvent (Lagow-Exfluor process) [START_REF] Lagow | The discovery of successful direct fluorination synthesis: Three eras of elemental fluorine reaction chemistry[END_REF], and perfected by grafting a heavily fluorinated ester ligand to the substrate to reduce vaporization and enhance solubility in a relatively heavy chlorofluoroether solvent (the "PERFECT" method) [START_REF] Okazoe | Development of the ''PERFECT'' direct fluorination method and its industrial applications[END_REF]. Heptafluorocubane was thus produced and the remaining hydrogen was replaced using standard monofluorination reactions. The precious hexa-and hepta-fluorinated cubanes were also isolated.

Notably, X-ray diffraction showed essentially no C-C bond length change upon perfluorination, likely meaning that the C-C bond-shortening expected from highly electronegative fluorine substitution is counterweighted by the fluorines' increased spatial needs. X-rays also revealed a rare case of directional non-covalent C-bonding in the solid state between a fluorine atom and the center of the face of the nearest cube.

Most importantly, Irikura's prediction that in perfluorocubane the overlap of all eight antibonding orbitals of the C-F bonds (Fig. 1A) would result in the formation of an extended star-shaped low energy antibonding molecular orbital (sigma stellation), and thus in pronounced electron affinity and the possibility to host an electron inside the 8-carbon cubic frame, hence generating a radical anion [START_REF] Irikura | Sigma stellation: A design strategy for electron boxes[END_REF], certainly added considerable incentive to the challenge.

Reduction potentials measured for perfluorocubane by cyclic voltammetry and differential pulse spectrometry confirmed its predicted electron-accepting capacity and the formation of a radical anion [START_REF] Sugiyama | Electron in a cube: synthesis and characterization of perfluorocubane as an electron acceptor[END_REF]. Low-temperature matrix-isolation electron spin resonance spectroscopy, in combination with g-ray radiolysis, determined a rapidly rotating radical anion. Uniquely, the spin density of the sequestered electron is chiefly located inside the cubane frame and displays octahedral (or cubic) symmetry, in contradistinction with the common p-conjugated polyaromatic electron acceptors that store extra electrons on their surface. The new radical anion could, however, not be isolated. About evoking "cages" (aka bird-cages) when portraying molecular structures: such pictorial denominations can be delusive (Fig. 1A). Cages inevitably suggest possibilities for occupancy and incarceration. But electrons can readily clear out, and it takes the 20-carbon dodecahedrane (C20H20) to encage the smallest monoatomic molecule, helium. Besides, the Authors' calculations (Fig. 1A) indicate that part of the caged electron spin distribution is actually located outside the 8-carbon cubic frame, near the fluorine atoms. "Frame" may sound both less engaging and less coercive.

Since its inception almost 60 years ago (2), cubane has inspired an outbreak of investigations, derivatizations, functionalizations, and suggestions for potential uses [START_REF] Chalmers | Validating Eaton's hypothesis: Cubane as a benzene bioisostere[END_REF][START_REF] Wlochal | Cubanes in medicinal chemistry: Synthesis of functionalized building blocks[END_REF][START_REF] Biegasiewicz | Cubane: 50 Years later[END_REF][START_REF] Collin | Cubane electrochemistry: Direct conversion of cubane carboxylic acids to alkoxy cubanes using the Hofer-Moest reaction under flow conditions[END_REF][START_REF] Reekie | Cubanes in medicinal chemistry[END_REF][START_REF] Houston | The cubane paradigm in bioactive molecule discovery: further scope, limitations and the cyclooctatetraene complement[END_REF]. Will perfluorocubane and its derivatives raise similar elation? Perfluorination causes sweeping steric and electronic alterations, entailing profound changes in properties and reactivity, making forecasts hazardous. Only further experimentation and modeling studies can tell.

Future investigations should certainly include the synthesis of additional mono-and 1,4disubstituted functional fluorinated cubane derivatives, developing means of introducing such building blocks into molecules and polymers, and assessing their "electron chemistry". Perfluorocubane-derived compounds may offer a unique opportunity to experimentally assess what differentiates an internally netted, s-star-located electron from an externally located, pconjugated one. Exploration of the rare global carbon-bonding donor capacity of perfluorocubane (revealed in the solid state by the X-ray diffraction studies) towards cyanide and similarly slim anions is desirable. Use of functional perfluorocubane building blocks as linkers, spacers or appendages in polymer chemistry may also be rewarding.

On the applied side, cubane, which can now be produced in kg-size batches, has triggered intensive investigations in medicinal chemistry. Numerous series of derivatives and building blocks with pharmaceutically significant functionalities were synthesized [START_REF] Chalmers | Validating Eaton's hypothesis: Cubane as a benzene bioisostere[END_REF][START_REF] Wlochal | Cubanes in medicinal chemistry: Synthesis of functionalized building blocks[END_REF][START_REF] Biegasiewicz | Cubane: 50 Years later[END_REF][START_REF] Collin | Cubane electrochemistry: Direct conversion of cubane carboxylic acids to alkoxy cubanes using the Hofer-Moest reaction under flow conditions[END_REF][START_REF] Reekie | Cubanes in medicinal chemistry[END_REF][START_REF] Houston | The cubane paradigm in bioactive molecule discovery: further scope, limitations and the cyclooctatetraene complement[END_REF]. Introducing a cubyl in lieu of a phenyl group in pharmaceuticals was attractive because the two groups are considered bioisosteric [START_REF] Chalmers | Validating Eaton's hypothesis: Cubane as a benzene bioisostere[END_REF][START_REF] Houston | The cubane paradigm in bioactive molecule discovery: further scope, limitations and the cyclooctatetraene complement[END_REF][START_REF] Eaton | Cubanes: Starting materials for the chemistry of the 1990s and the new century[END_REF], that is, have similar space requirements (although one is globular and the other rather flat), while the former appears less toxic and slightly more lipophilic. A further incentive for such substitutions was to help rejuvenate "old" drugs. Cubane derivatives helped probing active centers and mechanisms of biological action. In a study on cubyl-analogs of ten clinically used phenyl-containing molecules (anesthetic, analgesic, antidepressant, anticancer, antiparasitic and neotrophic agents) or pesticides, the cubane analog performed as well or better than the benzene analog in seven instances [START_REF] Houston | The cubane paradigm in bioactive molecule discovery: further scope, limitations and the cyclooctatetraene complement[END_REF]. Use as biologically safe cross-linkers in polymeric prodrugs for controlled release was also suggested.

Yet, to our knowledge, no cubane-containing pharmaceutical has so far entered clinical evaluation.

Anyway, those who would consider investigating perfluorocubane along similar lines must realize that the substantially larger perfluorocubane is no longer isosteric to benzene (Fig. 1B) and has significant lipophobic character (in addition to being extremely hydrophobic), forecasting little analogy in behavior. The extra-costs of optimizing and achieving perfluorination may also be prohibitive.

Considering materials science, polymerization of cubane has provided extremely hard, thin and stable nanothreads with exceptional stiffness [START_REF] Huang | Nanoarchitecture through strained molecules: Cubane-derived scaffolds and the smallest carbon nanothreads[END_REF], a feat that, in the case of perfluorocubane, could be hampered by steric hindrance.

But the truly inciting asset of perfluorocubane and its derivatives lies in their distinctive redox activity, electron hosting capacity and, more generally, "electron chemistry". Further exploration of these features can only be encouraged. New practical means, different from aromatics, of grafting a small, space-filling electron-sequestering building block into molecules, would indeed be valuable. Practical, meaning stable and cost-effective, redox-active, functional electron-carrying and transferring molecular materials and polymers may find various applications, including in electronic devices (organic conducting materials, rechargeable batteries), optoelectronics (light-emitting devices), magnetics (magnetic molecular chains and frameworks, NMR standard), magneto-optics, (dopants in organic semiconductors), therapeutics (photothermal materials), and other fields.

Currently, the availability of functional perfluorocubane-derived reagents on an appropriate scale remains an impediment to extensive use. The C8F8 •-radical anion is not stable and bottling cubane radical anions is not yet on the agenda. For most uses, it will be hard to compete, in terms of cost and accessibility, with the over a century old and extensively used polyaromatics-based radicals.

The Authors' findings certainly encourage perfluorination of further polyhedral molecular structures (Fig. 1C), especially those susceptible to achieve sigma stellation, and experimental investigation of their electrochemistry. Valuable candidates include perfluoroadamantane, which is commercially available, perfluorododecahedrane (C20F20) [START_REF] Zhang | Aromatic C20F20 cage and its endohedral complexes X@C20F20[END_REF] for which substantial ability to bind an extra electron has been predicted, star-shaped perfluorobisnoradamantane [START_REF] Adcock | Perfluorinated stellanes: The synthesis of perfluoronoradamantane. Long-range 19 F virtual coupling[END_REF] and perfluoronorbornane among a wealth of possibilities. Calculations on the C20F18(NH)2C20F18 double-cage predict electron transfer between cages and applications in nonlinear optical materials [START_REF] Wang | How does a double-cage single molecule confine an excess electron? Unusual intercage excess electron transfer transition[END_REF].

Meanwhile, this perfluorocubane study contributes brightly to our understanding of the chemistry and potential of electron-accepting perfluorinated polyhedral constructs, and should inspire plenty new explorations. (B) Perfluorocubane (left) and benzene (right) or cubane (center) are neither isosteric nor isolipophilic. (C) Some perfluorinated molecular polyhedrons waiting for their electron chemistry to be explored experimentally, perfluoro-adamantane, -bisnoradamantane,norbornane, -dodecahedrane (the symbol F inside the molecular frame indicates that all carbons are fully fluorinated); all have been synthesized, some are commercially available.

Caption for Figure 1 .

 1 (A) Small polyhedral molecular frames are hardly friendly cages (cartoon, courtesy J.-M. Dupart; they are extremely exiguous (left) and, where electron spin distribution is concerned (right), only partially inclusive, part of the added electron's spin density lying outside the cubane carbon frame (1). Cubane's electron-acceptor capacity derives from the centered 8-branch star-shaped molecular orbital formed by the overlap of the antibonding sigma orbitals of the C-F bonds (center; for clarity, only half of the orbitals are represented).


Yury Grabovsky

Lev Truskinovsky

A nonlinear elasticity problem with no local but many global minimizers unrelated by symmetry

In this paper, dedicated to 85th birthday of R. Fosdick, we touch upon two issues which attracted considerable attention in his own research: nonuniqueness in geometrically linear elasticity and the Clapeyron theorem. To this end we consider a simple model of a solid-solid phase transition with incompatible energy wells which exhibits multiplicity of global minimizers in a hard device, unrelated to either objectivity or material symmetry. At the same time, the same example shows the absence of strong local Lipschitz minimizers which are not global minimizers. Our proof that every strong local minimizer in the corresponding variational problem is also a global one hinges on the new general sufficiency theorem, proved by means of the novel nonlinear generalization of the Clapeyron theorem.

f( )

Introduction

The phenomenon of metastability in elastostatics, manifesting itself through the existence of strong local minimizers which are not global, is usually associated with a Neumann problem (soft device) and is linked to the incompatibility of the energy wells [START_REF] Ball | Incompatible sets of gradients and metastability[END_REF]. In this paper we present an analytically transparent example of an energy with incompatible wells for which one can prove the absence of strong local minimizers which are not global on any domain and for any Dirichlet boundary conditions (hard device).

Moreover, for the same energy one can prove dramatic nonuniqueness of global minimizers. Previously, it has been understood that to ensure uniqueness, the use of Dirichlet boundary conditions and topological, or even geometric simplicity of the domain are essential [START_REF] Post | On homotopy conditions and the existence of multiple equilibria in finite elasticity[END_REF][START_REF] Taheri | Local minimizers and quasiconvexity-the impact of topology[END_REF][START_REF] Emanuele | Non-uniqueness of minimizers for strictly polyconvex functionals[END_REF]. Uniqueness has been established for star-shaped domains, affine displacement boundary conditions, and strictly quasiconvex stored energy functions [START_REF] Knops | Quasiconvexity and uniqueness of equilibrium solutions in nonlinear elasticity[END_REF][START_REF] Taheri | Quasiconvexity and uniqueness of stationary points in the multi-dimensional calculus of variations[END_REF]. Whether uniqueness holds may also depend on the regularity class in which one looks for a minimizer [START_REF] Rj Knops | Uniqueness and complementary energy in nonlinear elastostatics[END_REF], or even on its integrability class [START_REF] Ball | Discontinuous equilibrium solutions and cavitation in nonlinear elasticity[END_REF][START_REF] Sivaloganathan | Singular minimisers in the calculus of variations: a degenerate form of cavitation[END_REF][START_REF] Co Horgan | Cavitation in nonlinearly elastic solids: a review[END_REF]. For mixed boundary-value problems of nonlinear elasticity nonuniqueness is common with the most familiar examples being those associated with buckling, due to the emergence of multiple symmetry-related energy minima [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes sive solutio problematis isoperimetrici latissimo sensu accepti, Additamentum I. De curvis elasticis[END_REF][START_REF] Grabovsky | The flip side of buckling[END_REF]. The possibility of nonuniqueness with Dirichlet boundary conditions was shown for non-homogeneous problems [START_REF] Edelstein | A note on non-uniqueness in linear elasticity theory[END_REF][START_REF] Le | An example of H 1 -unboundedness of solutions to strongly elliptic systems of partial differential equations in a laminated geometry[END_REF][START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF]. Our simple example shows that the multiplicity of global minimizers in a hard device problem can be obtained even in the absence of geometrical complexity of the domain, non-homogeneity, and can be unrelated to either objectivity or material symmetry. It is very similar to examples related to the nonuniqueness of optimal microstructures in composites [START_REF] Avellaneda | Optimal bounds and microgeometries for elastic two-phase composites[END_REF][START_REF] Milton | Variational bounds on the effective moduli of anisotropic composites[END_REF][START_REF] Grabovsky | Bounds and extremal microstructures for two-component composites: A unified treatment based on the translation method[END_REF][START_REF] Lu | Elastic energy minimization and the shape of coherent precipitates[END_REF].

More specifically, we consider a "geometrically linearized" Hadamard material which is a simplification of the fully nonlinear Hadamard material [START_REF] Hadamard | Leçons sur la propagation des ondes et les équations de l'hydrodynamique[END_REF][START_REF] John | Plane elastic waves of finite amplitude. hadamard materials and harmonic materials[END_REF] whose energy density function is of the form

W (F ) = µ|F | 2 + h(det F ), (1.1)
where µ is the measure of rigidity and the non-negative function h(d) is defined on (0, +∞) and has the property that h(d) → ∞, as d → 0 + . If we use the "geometric" approximation det F ≈ 1 + Tr (F -I), which is valid in the limit F → I, such formal asymptotic expansion with respect to a small parameter would also induce physical linearization and would trivialize the problem. To retain physical non-linearity we consider the energy W (H) = g(Tr H) + µ|H| 2 , H = F -I.

(1.2)

We may now perform geometric linearization in the second term as well [START_REF] Khachaturyan | Theory of structural transformation in solids[END_REF][START_REF] Budiansky | Continuum theory of dilatant transformation toughening in ceramics[END_REF][START_REF] Kaganova | Equilibrium between elastically-interacting phases[END_REF][START_REF] Abeyaratne | Dilatationally nonlinear elastic materials-I. Some theory[END_REF][START_REF] Abeyaratne | Dilatationally nonlinear elastic materials-II. an example illustrating stress concentration reduction[END_REF], replacing H with its symmetric part ε = (H + H t )/2 and, to emphasise its isotropic nature, writing the energy density as

W (ε) = f (Tr ε) + µ|dev(ε)| 2 , dev(ε) = ε - 1 3 Tr (ε)I. (1.3)
Mathematically, the analyses of (1.2) and (1.3) are very similar. We will therefore, focus on (1.3), because it preserves at least the linearized version of the frame indifference property. We assume that µ > 0 and the function f has the "double-well" shape. For analytical transparency we use in our explicit constructions the bi-quadratic potential f (θ) = min{κ 0 θ 2 , κ 0 (θ -θ p ) 2 + f 0 }, (1.4) illustrated in Fig. 1. The ensuing model describes a material capable of undergoing a purely dilatational phase transformation between two phases which are both linearly elastic. They have the same moduli and differ only by the transformation strain and the reference (chemical) energy. The relaxation of such energies has been always considered as a benchmark case for all theories of elastic phase transitions [START_REF] Khachaturyan | Theory of structural transformation in solids[END_REF][START_REF] Kohn | The relaxation of a double-well energy[END_REF][START_REF] Allen C Pipkin | Elastic materials with two preferred states[END_REF].

We recall that if the wells are rank-one connected, a quasiconvex energy would have to be convex, due to the rank-one conexity of quasiconvex functions [START_REF] Morrey | Quasi-convexity and the lower semicontinuity of multiple integrals[END_REF]. On the other hand, the incompatibility of the wells can be responsible for quasiconvexity of a non-convex energy density, as we will see in our example, where the energy wells are not rank one connected. In our example, despite the fact that the non-convexity enters (1.3) only through a scalar potential f (θ), the energy density (1.3) is non-quasiconvex for all µ > 0, if f (θ) is given by (1.4), [START_REF] Kohn | The relaxation of a double-well energy[END_REF]. Therefore, even if subjected to homogeneous loading in a hard device, the materials described by (1.4) must undergo non-homogeneous deformations by developing microstructures. In the presence of local minima, such systems can be expected to experience hysteresis [START_REF] Ball | Hysteresis during stress-induced variant rearrangement[END_REF][START_REF] Yongmei M Jin | Macroscopic energy barrier and rate-independent hysteresis in martensitic transformations[END_REF][START_REF] Roytburd | Thermodynamic hysteresis of phase transformation in solids[END_REF][START_REF] Šilhavý | the hysteresis in martensitic transformations[END_REF][START_REF] Zhang | Energy barriers and hysteresis in martensitic phase transformations[END_REF]. However, as we have already mentioned, the local minima in this model are absent. Then, the macroscopic quasistatic mechanical response in the hard device corresponding to global minimization of the total energy at each value of the loading parameter is defined uniquely, and as a result, the deformation path is reversible. The associated stress-strain relation can be obtained from the knowledge of the quasiconvex envelope of the energy (1.3) that can be characterized explicitly in the whole range of parameters and for all double-well shaped potentials f (θ). In this paper we perform the implied relaxation of the energy (1.3) using simple laminates as the optimal microstructure; a different proof based on matching of upper and lower bounds can be found in [START_REF] Grabovsky | When rank-one convexity meets polyconvexity: An algebraic approach to elastic binodal[END_REF].

One of the general results of this paper is the formulation of the necessary and sufficient conditions for global minimizers in a star-shaped domain with affine boundary conditions. This creates a tool to characterize energy minimizing configurations in examples where the energy relaxation, or at least the elastic binodal, is explicitly known. An important technical tool in this analysis is a fundamental nonlinear generalization of the Clapeyron theorem which is presented here for the first time. These results can be viewed as a generalization of optimality conditions for extremal microstructures in composites (e.g. [START_REF] Grabovsky | Bounds and extremal microstructures for two-component composites: A unified treatment based on the translation method[END_REF]).

We show that using our necessary and sufficient conditions one can reduce the problem of finding a global minimizer to the solution of a nonlinear free boundary problem involving a system of linear PDEs in a finite domain. It is remarkable that for the energy density (1.3) the solution of the ensuing problem can be found explicitly, for any double-well scalar potential f (θ), providing us with an explicit example of multiplicity of global minimizers. In particular, we were able to compute explicitly a one parameter family of low surface energy non-affine configurations in the square domain, whose boundary is piecewise smooth. In this way we showed the possibility to always generate in this class of problems an optimal microstructure with low surface area, which is physically advantageous comparing, for instance, to multiscale lamination with infinitely large surface area [START_REF] Kohn | Surface energy and microstructure in coherent phase transitions[END_REF].

The analysis of optimal microstructures suggests that, despite the absence of the total energy hysteresis in this system, the direct and the reverse transformation may follow different transformation paths in the configuration space. Therefore in this problem, in the absence of metastability and "constitutive hysteresis", we encounter a phenomenon of "morphological hysteresis", whereby the forward and reverse transformations occur along different morphological paths, while traversing energetically equivalent configurations.

The paper is organized as follows. In Section 2 we discuss the restrictions on the structure of strong local minimizers. The necessary and sufficient conditions for global minimizers are formulated in Section 3. The issue of attainment and the multiplicity of global minimizers in our model are discussed in Section 4. An explicit construction of a nontrivial global minimizer in a square domain which has expressly low surface area is presented in Section 5. The last Section 6 contains our conclusions. An Appendix A contains a technical discussion of the generic loss of ellipticity for the energy represented by a rank one convex envelope.

Local minimizers

Given that the case of interest involves hard device loading, our goal is to study strong local and global minima of the functional

E 0 (u) = Ω f (∇ • u) + µ e(u) - 1 3 (∇ • u)I 2 dx (2.1)
among all Lipschitz displacement vector fields u ∈ W 1,∞ (Ω; R 3 ) subject to the constraint

u(x) = u 0 (x), x ∈ ∂Ω, (2.2)
where u 0 (x) is a given Lipschitz function on ∂Ω. Restricting attention to Lipschitz minimizers allows us to focus on instabilities caused by failure of quasiconvexity [START_REF] Morrey | Quasi-convexity and the lower semicontinuity of multiple integrals[END_REF][START_REF] Dacorogna | Quasiconvexity and relaxation of nonconvex problems in the calculus of variations[END_REF] and exclude some other instabilities related, for instance, to the mismatch between the integrability of the minimizing sequences and the growth of the energy density at infinity [START_REF] Ball | W 1,p -quasiconvexity and variational problems for multiple integrals[END_REF].

First we recall that one of the necessary conditions of metastability, understood here as the lack of strong local Lipschitz minimizers which are not also global minimizers, is stability with respect to nucleation of a coherent precipitate of the new phase in the interior of the old phase. That means that we must have ∇u(x) ∈ A, for a.e. x ∈ Ω, [START_REF] Taheri | Local minimizers and quasiconvexity-the impact of topology[END_REF]Proposition 4.1], where

A = {H ∈ M : W (H) = QW (H)}, (2.3)
where M denotes the set of all 3 × 3 matrices and QW (H) is the quasiconvexification of W [START_REF] Dacorogna | Quasiconvexity and relaxation of nonconvex problems in the calculus of variations[END_REF].

In this paper we show that if u(x) is a weak solution of the Euler-Lagrange equation for the energy (2.1) with the boundary conditions (2.2) and satisfies ∇u(x) ∈ A for a.e.

x ∈ Ω then, it is necessarily a global minimizer for E 0 (u). However, before we move to the actual analysis, it is appropriate to mention that we are not aware of any examples of proper metastable states in a hard device in domains with trivial topology. Moreover it has been proved for general energies that every strong local minimizer is global in the hard device with affine Dirichlet boundary conditions in star-shaped domains [START_REF] Taheri | Quasiconvexity and uniqueness of stationary points in the multi-dimensional calculus of variations[END_REF]. In our special example, the same result turns out to be true for all Dirichlet boundary conditions in domains with piecewise smooth boundaries and arbitrary topology.

The first step of the analysis is to characterize the set A of admissible displacement gradients for the chosen material model. It can result from the computation of QW (H) which in our case can be done explicitly.

Theorem 2.1. Suppose that W (H) is given by (1.3). Then,

QW (H) = Φ * * (Tr H) + µ 4 |H -H t | 2 -2µJ 2 (H), (2.4)
where

Φ(θ) = f (θ) + 2 3 µθ 2 , (2.5) 2J 2 (H) = (Tr H) 2 -Tr (H 2 ) = 2 i<j h ii h ij h ji h jj (2.6)
is a quadratic null-Lagrangian, and Φ * * denotes the convexification of Φ(θ).

Proof. It will be convenient to use the definition of quasiconvexification based on periodic functions [START_REF] Dacorogna | Quasiconvexity and relaxation of nonconvex problems in the calculus of variations[END_REF]:

QW (H) = inf φ-Q-periodic - Q W (H + ∇φ(x))dx, (2.7)
where

Q = [0, 1] 3 .
The computation of QW is based on the following formula

e(u) - 1 3 (∇ • u)I 2 = 2 3 (∇ • u) 2 + 1 2 |∇ × u| 2 -2J 2 (∇u). (2.8) Let u(x) = Hx + φ(x). Then - Q W (∇u)dx = -2µJ 2 (H) + - Q Φ(∇ • u) + µ 2 |∇ × u| 2 dx. (2.9)
Since Φ * * (θ) ≤ Φ(θ) for all θ ∈ R, we obtain the inequality

- Q W (∇u)dx ≥ -2µJ 2 (H) + - Q Φ * * (∇ • u) + µ 2 |∇ × u| 2 dx.
Now applying the Jensen's inequality we get:

- Q W (∇u)dx ≥ -2µJ 2 (H) + Φ * * (Tr H) + µ 4 |H -H t | 2 . (2.10) Therefore, QW (H) ≥ -2µJ 2 (H) + Φ * * (Tr H) + µ 4 |H -H t | 2 . (2.11)
In order to prove equality, we need to exhibit a periodic function φ such that ∇ × φ = 0 and -

Q Φ(∇ • φ(z) + Tr H)dz = Φ * * (Tr H).
(2.12)

If Φ(θ), given by (2.5) is convex1 , then the original energy W (H) is quasiconvex (polyconvex, actually), and formula (2.4) holds. The case of interest is therefore when Φ(θ) retains the double-well shape, in which case Φ * * (θ) = Φ(θ) if and only if θ ∈ (θ 1 , θ 2 ), as in Fig. 2. For example, for the potential (1.4) we find

θ 1 = f 0 2κ 0 θ p + µθ p 3κ 0 + 2µ , θ 2 = θ p + f 0 2κ 0 θ p - µθ p 3κ 0 + 2µ . (2.13) 1 2 ( ) Figure 2: Definition of θ 1 < θ 2 via the common tangent construction. Thus, if Tr H ∈ (θ 1 , θ 2 ), then φ = 0 attains equality in (2.10). If Tr H ∈ (θ 1 , θ 2 ), then there exists ω ∈ (0, 1) such that Tr H = ωθ 1 + (1 -ω)θ 2 . (2.14)
In order to attain equality in (2.10) we need to split the period cell

Q = [0, 1] 3 into two subsets A 1 , of volume ω and A 2 of volume 1 -ω such that ∇ • u(x) = θ 1 χ A 1 (x) + θ 2 χ A 2 (x) (2.15)
for all x ∈ Q. This is easily achieved with a laminate construction with

A 1 = [0, ω] × [0, 1] 2 and A 2 = [1 -ω, 1] × [0, 1]
2 and function ∇φ being symmetric and piecewise constant:

∇φ(x) = M 1 χ A 1 (x) + M 2 χ A 2 (x).
The symmetric matrices M 1 and M 2 are easily found from the periodicity of φ(x), that is equivalent to -Q ∇φ(x)dx = 0, and (2.15):

M 1 = (1 -ω)[[θ]]e 1 ⊗ e 1 , M 2 = -ω[[θ]]e 1 ⊗ e 1 ,
where

[[θ]] = θ 1 -θ 2 .
These choices guarantee equality in (2.12) and the symmetry of ∇φ

guarantees ∇ × φ = 0. Corollary 2.2. A = {H ∈ M : Tr H ∈ (θ 1 , θ 2 )}. (2.16)
The constructed energy density QW (H) is obviously polyconvex, and coincides with RW (H), the rank-1 convexification of W (H). Indeed, as we see from the proof of Theorem 2.1, the infimum in (2.7) can be achieved by a simple laminate. In view of the simplicity of the optimal microstructures, one can post factum conclude that computing QW (H) explicitly was not really necessary for the characterization of the set A. We could also delineate it by using inequality (2.10) and relying on the knowledge and properties of the jump set introduced in [START_REF] Grabovsky | Roughening instability of broken extremals[END_REF][START_REF] Grabovsky | Normality condition in elasticity[END_REF]. Indeed, inequality (2.10) implies that all H satisfying Φ * * (Tr H) = Φ(Tr H) are admissible, i.e. {H ∈ M : Tr H ∈ (θ 1 , θ 2 )} ⊂ A. We can then show that all H with θ 1 < Tr H < θ 2 fail rank-one convexity by looking for pairs H ± of the displacement gradient values that can occur on the two sides of a smooth jump discontinuity representing a phase boundary. We recall, [START_REF] Grabovsky | Roughening instability of broken extremals[END_REF], that these matrices satisfy the equations

[[H]] = a ⊗ n, [[W H ]]n = 0, [[W H ]] t a = 0, [[W ]] = { {W H } }, [[H]] , (2.17)
where

[[H]] = H + -H -, { {W H } } = 1 2 (W H (H + ) + W H (H -)).
We have shown in [START_REF] Grabovsky | When rank-one convexity meets polyconvexity: An algebraic approach to elastic binodal[END_REF] that the jump set for the energy (1.3) consists of the union of two hyperplanes

J -= {H ∈ M : Tr H = θ 1 }, J + = {H ∈ M : Tr H = θ 2 },
so that if Tr H -= θ 1 , then the set of corresponding H + is a projective plane worth of points

H + = H -+ (θ 2 -θ 1 )n ⊗ n, |n| = 1.
As we have already remarked, inequality (2.10) implies that J ± ⊂ A. In this case the following lemma shows that our results from [START_REF] Grabovsky | Normality condition in elasticity[END_REF] imply that all matrices H t = tH + + (1t)H -would also fail to be rank-one convex for all t ∈ (0, 1 Then all of the matrices H t , t ∈ (0, 1) must fail rank-one convexity in the sense that they must fail to satisfy the inequality

W (H) ≥ W (H t ) + P t , H -H t ∀H : rank(H -H t ) = 1, (2.18)
where

P t = W H (H t ).
Proof. Suppose, with the goal of getting a contradiction, that (2.18) holds. Then, since matrices H -, H t and H + all lie on the same rank-one line, we obtain

W (H -) + t P -, [[H]] ≤ W (H t ) ≤ W (H + ) -(1 -t) P t , [[H]] .
Now, using the jump set equations (2.17), we obtain

(1 -t)( P -, [[H]] -P t , [[H]] ) ≥ 0.
Using the fact that H t -H -= t[[H]] we obtain

P t -P -, H t -H -≤ 0. (2.19)
The variational significance of inequality (2.19), called the phase interchange stability inequality, was understood in [START_REF] Grabovsky | Normality condition in elasticity[END_REF], where, according to [START_REF] Grabovsky | Normality condition in elasticity[END_REF]Lemma 4.1], the pair of rank-one related matrices H -and H t would have to satisfy the jump set equations (2.17), in contradiction with the assumption (ii) of the lemma. Thus, the assumption that (2.18) holds cannot be true, and all matrices H t , t ∈ (0, 1), must fail to be rank-one convex.

Remark 2.4. The energy density QW (H) = RW (H) is non-convex as it has a double-well shape along the multiples of I, specifically, QW (εI) = Φ * * (3ε) -6µε 2 . This is clear from the fact that its bulk modulus is negative in the interval θ ∈ (θ 1 , θ 2 ), where θ 1 and θ 2 are given by (2.13) for the bi-quadratic material model (1.4). Moreover, the value of this modulus is constant and equal to -(4/3)µ, cf. [START_REF] Budiansky | Continuum theory of dilatant transformation toughening in ceramics[END_REF], which corresponds exactly to the threshold for the loss of strong ellipticity of the equilibrium equations (saturation of the Legendre-Hadamard conditions). This is associated with the generic degeneration of the acoustic tensor along rank one envelopes as we explain in our Appendix A. The vanishing of the velocity of the longitudinal waves in our example turns the phase-transforming material into an elastic aether [START_REF] Ericksen | Implications of Hadamard's conditions for elastic stability with respect to uniqueness theorems[END_REF].

We are now ready to show that for a class of energies to which our example belongs all strong local minimizers in a hard device (Dirichlet boundary conditions) must be global. The first observation, leading to this result is that any two energies W and W that have one and the same quasiconvex envelope are equivalent as far as existence of metastable states are concerned. This is because for any metastable state u

Ω W (∇u)dx = Ω QW (∇u)dx.
Thus, any metastable configuration for W is also metastable for QW . Theorem 2.1 shows that our energy has the property

QW (H) = N (H) + C(H), (2.20)
where N (H) is a null-Lagrangian and C(H) is convex and C 1 smooth. Proof. Let v be a Lipschitz competitor that agrees with u on ∂Ω.

Then Ω N (∇v)dx = Ω N (∇u)dx, C(∇v) ≥ C(∇u) + C F (∇u), ∇(v -u) . (2.21)
Since u is an equilibrium we have

0 = ∇ • W F (∇u) = ∇ • QW F (∇u) = ∇ • N F (∇u) + ∇ • C F (∇u) = ∇ • C F (∇u).
The equality of W F (H) = QW F (H) for any H ∈ A follows from the fact that the function W (H) -QW (H) is nonnegative, C 1 smooth and attains its minimum value of 0 at all F ∈ A. Now, the integration by parts in the second term on the right-hand side in the inequality in (2.21) implies Ω C(∇v)dx ≥ Ω C(∇u)dx. It follows that

Ω QW (∇v)dx ≥ Ω QW (∇u)dx.
Remark 2.6. An important observation was made by Sivaloganathan and Spector in [START_REF] Sivaloganathan | On the uniqueness of energy minimizers in finite elasticity[END_REF] that C(H) does not need to be convex for inequality in (2.21) to hold. In fact, as the authors show, it is only required that C * * (∇u(x)) = C(∇u(x)) for all x ∈ Ω, i.e. C(F ) agrees with its convex hull at all values of ∇u.

Global minimizers

Consider next the problem of attainment in the definition of the quasiconvex envelope

QW (F ) = inf y∈F x+W 1,∞ 0 (Ω;R m ) - Ω W (∇y)dx, (3.1)
where Ω ⊂ R d is a star-shaped Lipschitz domain. Below we show that Lipschitz equilibrium configurations for (3.1) must be global minimizers, provided they satisfy the necessary conditions for metastability and are sufficiently regular near ∂Ω.

Recall first the well-known necessary conditions for Lipschitz strong local minima of variational functionals.

The first classical necessary condition is the Euler-Lagrange equation

∇ • P (∇y) = 0, P (F ) = W F (F ), (3.2)
which even weak local minimizers have to satisfy. The second, is the Noether equation The second classical necessary condition is quasiconvexity ∇y ∈ A for a.e. x ∈ Ω, where A is defined in (2.3). It was shown in [START_REF] Grabovsky | Normality condition in elasticity[END_REF] that for configurations, whose only singularities are smooth phase boundaries the combination of (3.2) and quasiconvexity implies stationarity. In general, however, an example in [START_REF] Kristensen | Partial regularity of strong local minimizers in the multidimensional calculus of variations[END_REF] shows that there could be weak local minimizers for strictly quasiconvex energies that are not strong local minimizers; most probably, the configuration in [START_REF] Kristensen | Partial regularity of strong local minimizers in the multidimensional calculus of variations[END_REF] fails stationarity.

∇ • P * (∇y) = 0, P * (F ) = W (F )I d -F t P (F ). (3.3) We recall that if y(x) is of class C 2 then (3.3) is a consequence of (3.2),
It turns out that in the case affine boundary conditions in a star-shaped domain, the two necessary conditions formulated above are sufficient for y(x) to be a global minimizer, provided we assume its regularity near the boundary of the domain. (

Proof. While each ingredient of the proof presented below was in one way or another already present in the uniqueness proof of Knops and Stuart [START_REF] Knops | Quasiconvexity and uniqueness of equilibrium solutions in nonlinear elasticity[END_REF], to express the energy of an equilibrium configuration as a boundary integral, we use a far reaching generalization of a related relation in [START_REF] Knops | Quasiconvexity and uniqueness of equilibrium solutions in nonlinear elasticity[END_REF] which, by itself, can be viewed as a nontrivial nonlinear generalization of the Clapeyron theorem [START_REF] Fosdick | About Clapeyron's theorem in linear elasticity[END_REF].

Lemma 3.3. Suppose y(x) is a Lipschitz stationary equilibrium in a Lipschitz domain Ω. Then E[y] = Ω W (∇y) = 1 d ∂Ω {P n • y + P * n • x}dx, (3.5)
where P n and P * n can be regarded as trace functionals, since they act on Lipschitz functions y(x) and x, respectively.

Proof. We now apply formula (3.5) for the energy of y(x). Using the boundary condition (iv) and the assumption (v) we have at all x ∈ ∂Ω ∇y = F + a ⊗ n, a = ∂y ∂n -F n.

We therefore compute

dE[y] = ∂Ω (W (F + a ⊗ n) -P (F + a ⊗ n)n • a)(n • x)dS
next, we use the local stability condition (iii) and appeal to [START_REF] Grabovsky | Normality condition in elasticity[END_REF]Lemma 4.2]. We quote the part of the lemma we need for the sake of completeness.

Lemma 3.4. Let V (F ) be a rank-one convex function such that V (F ) ≤ W (F ). Let

A V = {F ∈ O : W (F ) = V (F )}, where O is an open subset of M on which W (F ) is of class C 1 . Then for every F ∈ A V , u ∈ R m , and v ∈ R d V (F + u ⊗ v) ≥ W (F ) + P (F )v • u. (3.6)
We now apply Lemma 3.4 by choosing V (F ) to be QW (F ). By assumption (iii) for each x ∈ ∂Ω the field ∇y(

x) = F + a ⊗ n is in A V . Choosing u ⊗ v = -a ⊗ n, inequality (3.6) becomes QW (F ) ≥ W (F + a ⊗ n) -P (F + a ⊗ n)n • a. (3.7)
Finally, we use the assumption of Ω being star-shaped. If we choose the origin at the star point, then the function n(x) • x is always non-negative at all points on ∂Ω. Therefore, inequality (3.7) implies

dE[y] ≤ QW (F ) ∂Ω (n • x)dS = d|Ω|QW (F ).
Since |Ω|QW (F ) is the minimal value of E[y] we conclude that E[y] = QW (F ), and y(x) is the global minimizer.

Multiplicity of global minimizers

We remark that the construction of the periodic laminate in Section 2 provides us with a minimizing sequence in (3.1) in arbitrary domains. In this section we raise the question of attainability of the minimum in (3.1). In this regard, laminates can no longer be used as the boundary conditions (4.2) are inconsistent with ∇u taking exactly two specific values everywhere in Ω. We will therefore attempt to find energy minimizers using the sufficiency conditions from Theorem 3.2. Accordingly, we only need to find solutions u(x) of (3.2) that satisfies (4.2) while respecting the stability condition

∇ • u(x) ∈ (θ 1 , θ 2 ) for a.e. x ∈ Ω. (4.1)
In what follows, in the interest of analytical transparency, we restrict our analysis to specific categories of bounded domains Ω and assume affine boundary conditions

u(x) = H 0 x, x ∈ ∂Ω. (4.2)
If Tr H 0 ∈ [θ 1 , θ 2 ] then, at least for star-shaped domain the theorem of Knops and Stuart [START_REF] Knops | Quasiconvexity and uniqueness of equilibrium solutions in nonlinear elasticity[END_REF] (see also [START_REF] Taheri | Quasiconvexity and uniqueness of stationary points in the multi-dimensional calculus of variations[END_REF]) ensures that the solution u(x) = H 0 x is the unique minimizer. Therefore we only need to consider the case Tr H 0 ∈ (θ 1 , θ 2 ).

Let φ(x) = u(x) -H 0 x. Then (3.2) can be written in terms of φ as follows

µ∆φ + ∇Φ (∇ • u) = µ∇(∇ • φ), x ∈ Ω. (4.3)
Integrating the dot product of φ(x) and ( 4.3) by parts we obtain

Ω µ(|∇ • φ| 2 -|∇φ| 2 ) -Φ (∇ • u)∇ • φ dx = 0.
Now we use the algebraic identity

|∇ • φ| 2 -|∇φ| 2 = 2J 2 (∇φ) -|∇ × φ| 2
and the fact that J 2 (∇φ) given by (2.6) is a null-Lagrangian to obtain

Ω µ|∇ × φ| 2 + Φ (∇ • u)∇ • φ dx = 0. (4.4)
Now let us show that (4.1) implies that

Ω Φ (∇ • u)∇ • φ(x)dx ≥ 0. (4.5)
Indeed, from the geometric interpretation of θ 1 and θ 2 (as the points of common tangency ) we conclude that in view of (4.1) we can write

Φ(∇ • u(x)) = Φ * * (∇ • u(x)) (4.6)
for a.e. x ∈ Ω. For a convex function Φ * * (θ) the tangent line to the graph of that function lies always below the graph which means that

Φ * * (η) ≥ Φ (θ)(η -θ) + Φ(θ)
for any η ∈ R and any θ ∈ (θ 1 , θ 2 ). Substituting η = Tr H 0 and θ = ∇ • u(x) we obtain

Φ (∇ • u(x))∇ • φ(x) ≥ Φ(∇ • u(x)) -Φ * * (Tr H 0 ) = Φ * * (∇ • u(x)) -Φ * * (Tr H 0 )
for a.e. x ∈ Ω. Now (4.5) follows from Jensen's inequality for a convex function Φ * * (θ) and the fact that 1

|Ω| Ω ∇ • u(x)dx = Tr H 0 . (4.7)
Now, inequality (4.5) together with (4.4) imply that ∇ × φ = 0, and in simply connected regions we may use the representation φ(x) = ∇h(x), where h ∈ W 2,2 0 (Ω) is some scalar potential. Substituting u = H 0 x + ∇h into (4.3), we obtain ∇Φ (∇ • u) = 0. Thus, there exists a constant P 0 , such that Φ (∇ • u) = P 0 .

(4.8)

The constraint (4.1) then implies that either ∇ • u = Tr H 0 + ∆h is constant in Ω, which means that u(x) = H 0 x + a 0 , or that P 0 = Φ (θ 1 ) = Φ (θ 2 ), in which case there exists a subset A of Ω such that

∇ • u(x) = θ 1 χ A (x) + θ 2 χ Ω\A (x), (4.9)
where θ 1 < θ 2 are the endpoints of the interval {θ ∈ R : Φ * * (θ) < Φ(θ)} (see Fig. 2).

In terms of the potential h we obtain a free boundary problem

∆h = (θ 1 -θ 2 )(χ A (x) -ω), x ∈ Ω h ∈ W 2,2 0 (Ω), (4.10)
where

ω = |A| |Ω| = Tr H 0 -θ 2 θ 1 -θ 2 (4.11)
is the volume fraction of the phase in which ∇ • u(x) = θ 1 . As formula (4.11) indicates, the volume fraction ω of the precipitate is uniquely determined by the hard device loading H 0 , whose shear component has no effect on the precipitate morphology. We remark that only the divergence of the displacement gradient is constrained by the requirement of optimality. This explains a virtual cornucopia of optimal microstructures, the simplest of which we are going to exhibit. A solution of (4.10) can be easily found in explicit form for some particularly simple domain geometries. For instance, it is easy to construct a radially symmetric solution when Ω is the unit ball, taking as the set A to be either the concentric ball of radius r 0 = ω 1/3 , or its complement, in which case r 0 = (1 -ω) 1/3 . Then (4.10) is solved by h(x) = h ω (|x|), given by

h ω (r) =          1 6 (θ 1 -θ 2 )(1 -ω) r 2 - 3ω(1 -ω 1/3 ) ω 1/3 (1 -ω) , 0 ≤ r ≤ ω 1/3 - 1 6 (θ 1 -θ 2 )ω r 2 + 2 r -3 , ω 1/3 ≤ r ≤ 1, (4.12)
in the former case. In the latter case, h(x) is given by the same formula (4.12) where θ 1 and θ 2 are interchanged and ω is replaced by 1 -ω, i.e. h(x) = -h 1-ω (|x|). It is not clear from the free boundary problem (4.10) if solutions exist in less symmetric domains. In fact, Hashin's "concentric sphere" construction, [START_REF] Hashin | The elastic moduli of heterogeneous materials[END_REF], shows that (4.10) has solutions in any Lipschitz domain Ω. In this construction the set A is a countable union of variously scaled copies of the radially symmetric solution (4.12) filling Ω up to a set of Lebesgue measure zero. Let B(x i , a i ) ⊂ B(x i , R i ), i = 1, 2, . . . are the concentric balls used in Hashin's construction, where

a 3 i R 3 i = ω.
The inner balls B(x i , a i ) in the construction belong to the set A, while the spherical shells

B(x i , R i ) \ B(x i , a i ) belong to Ω \ A. The function h(x) restricted to the ball B(x i , R i ) is given by h(x) = R 2 i h ω |x -x i | R i , x ∈ B(x i , R i ),
where h ω is given by (4.12). The function h(x) defined like this on each of the concentric balls does indeed solve (4.10), since on the boundary of each ball B(x i , R i ) both h(x) and ∇h(x) are zero. The problem with the above "concentric sphere" construction is the same one as with the laminate-based optimal microstructures: they all have an infinite surface area (for the "concentric spheres" see the formal proof in [START_REF] Mergelyan | Uniform approximations to functions of a complex variable[END_REF][START_REF] Wesler | An infinite packing theorem for spheres[END_REF]), which is unacceptable in physical problems of interest. However, we also expect (4.10) to have solutions with finite surface area in all piecewise smooth domains2 . For example, if Ω is an ellipsoid, then A would be a confocal ellipsoid (or its complement) [START_REF] Bergman | Exactly solvable microscopic geometries and rigorous bounds for the complex dielectric constant of a two-component composite material[END_REF][START_REF] Milton | Bounds on complex dielectric constant of a composite material[END_REF][START_REF] Tartar | Estimation fines des coefficients homogénéisés[END_REF][START_REF] Zhikov | Estimates for the homogenized matrix and the homogenized tensor[END_REF][START_REF] Grabovsky | Bounds and extremal microstructures for two-component composites: A unified treatment based on the translation method[END_REF]. In the next section we present a more complex solution with finite surface area in a two dimensional domain with piecewise-smooth boundary.

A finite perimeter global minimizer

For a square, the 2D version of (4.10) can be solved using complex variables. The actual geometry of the phase boundary can be then found numerically.

To obtain the shape of phase boundary we used the conventional complex potentials. Indeed, in view of (4.10) the functions

h + (x, y) = h(x, y) + 1 4 (θ 1 -θ 2 )ω(x 2 + y 2 ), h -(x, y) = h(x, y) - 1 4 (θ 1 -θ 2 )(1 -ω)(x 2 + y 2 )
are harmonic in Ω\A and A respectively. We can then conclude that the functions ∂ h ± /∂x-i∂ h ± /∂y are analytic in the complex variable z = x + iy on their respective domains. Thus, the functions

H + (z) = 2 θ 1 -θ 2 ∂h ∂x -i ∂h ∂y + ωz, z ∈ Ω \ A, H -(z) = 2 θ 1 -θ 2 ∂h ∂x -i ∂h ∂y -(1 -ω)z, z ∈ A
are also analytic. The boundary conditions in (4.10) and the continuity of ∇h across ∂A, representing the kinematic compatibility of the displacement, imply that

H + (z) = ωz, z ∈ ∂Ω, H + (z) -H -(z) = z, z ∈ ∂A. (5.1)
The knowledge of Ω allows one to determine H + . If Ω is a square centered at the origin with diagonal of length 2, then along the bottom side of the square we have

H + (x -i/ √ 2) =
ω(x + i/ √ 2) = ω(z + i √ 2). Therefore H + (z) = ω(z + i √ 2)
, while along the right side of the square we have

H + (1/ √ 2+iy) = ω(1/ √ 2-iy) = ω( √ 2-z). Thus, H + (z) = ω( √ 2-z)
. These contradictory expressions for H + (z) can only be reconciled by a structure with topology indicated in Figure 3, where the set A, in which ∇ • u = θ 1 , is shaded.

The square symmetry of a problem also suggests that we should look for a structure with square symmetry. Thus we just need to find the curve Γ that joins the top and bottom ends of the right side of the square, for example. Knowing H + (z) in all four regions adjacent to the sides of the square gives the boundary values for H -(z) on ∂A. The square symmetry of A together with (5.1) implies that

H -(iz) = -iH -(z), z ∈ ∂A. (5.2)
There is a holomorphic function H -in A with given boundary values if and only if for every n ≥ 0 ∂A H -(z)z n dz = 0.

(5.

3)

The integral in (5.3) can easily be written as an integral along Γ because of (5.2). If n is not a multiple of 4 then the integral in (5.3) will evaluate to zero by virtue of (5.2) alone. If n = 4k, then we require that Γ H -(z)z 4k dz = 0, k ≥ 0.

(5.4)

On Γ we have H -(z) = ω( √ 2 -z) -z. Therefore, (5.4) is equivalent to Γ zz 4k dz = (-1) k iω (4k + 1)(2k + 1)
, k ≥ 0.

(5.5)

If we integrate by parts in (5.5):

Γ zz 4k dz = Γ z d z 4k+1 4k + 1 = - Γ z 4k+1 4k + 1 dz.
and parametrize the curve Γ by

z(t) = a(t) + it, t ∈ [-1/ √ 2, 1/ √ 2]. Then z (t) = z (t) -2i, we obtain Γ zz 4k dz = (-1) k+1 i (4k + 1)(2k + 1) + 2i 4k + 1 1/ √ 2 -1/ √ 2 z(t) 4k+1 dt.
Thus, the equation (5.5) becomes

a(t) = 1 √ 2 + m(t - 1 √ 2 ) + O((t - 1 √ 2 ) 2 ), as t → 1 √ 2 .
(5.7)

Substituting the leading term in (5.6) and computing the integral explicitly we obtain

(-1) k (4k + 2)(m 2 + 1) - (1 -m) 4k+2 m 2 2k+1 (4k + 2)(m 2 + 1) ≈ (-1) k (ω + 1) 2(4k + 2) . (5.8)
If we choose m = 1 -ω 1 + ω (5.9)

then for large k we have equality in (5.6) up to an exponentially small error. Thus, we have determined the slope m with which the curve Γ enters the corner of the square.

In order to find the function a(t) numerically, we approximate it by even polynomials at t = 0 and enforce (5.7) at every approximation. Notice that the error term in (5.8) decays very fast with k. Hence, satisfying (5.6) only for small values of k should give a very good approximate solution. With the knowledge of this, we approximate a(t) by an even polynomial

a n (t) = n j=0 a j t 2j .
In order to satisfy (5.8) we require that

n j=0 a j 2 j = 1 √ 2 , n j=1 ja j 2 j-1 = m √ 2 .
(5.10)

Observe that equation (5.6) for k = 0 is linear. Hence,

n j=0 a j 2 j (2j + 1) = ω + 1 2 √ 2 .
(5.11)

To solve the system (5.6) for k = 0, 1, . . . , K we set n = K + 2 which leaves us with an easily solvable system of K polynomial equations in K unknowns. The configurations of the sets A and Ω \ A for ω = 0.3 (left panel) and ω = 0.8 (right panel) are shown in Fig. 3. When ω is almost zero, we effectively deal with phase nucleation associated with the "direct" (1 → 2) transformation. Phase "1" appears first in the form of two infinitely thin intersecting segments; the subsequent thickening of these segments starts first around their intersection point at the center of the square domain. As ω increases, the segments turn into a 4-cornered star, which eventually takes over the whole body. Instead, at ω almost one, we deal with the nucleation of phase "2" which appears first in the form of a thin "film" along the boundary of the body. In other words, during the "reverse" (1 → 2) transformation a new phase grows away from the boundary and eventually takes over the whole body leaving phase "1" only as a star shaped subdomain centered at the origin.

However, it is important to keep in mind that phases "1" and "2" in these scenaria can be always interchanged. The complete constitutive symmetry between the two phases (the fact that θ 1 < θ 2 is irrelevant) implies that the new phase may first appear either as an infinitely thin coating of the sides of the square (unshaded regions in the right panel in Fig. 3) or as a cross located at the center of the domain and extending its arms to its corners (shaded regions in the left panel in Fig. 3).

Conclusions

In this paper we considered a classical problem of nonlinear elastostatics for a material undergoing phase transition. Mathematically it reduces to a non-convex vectorial problem of the calculus of variations. In this context we presented a simple, yet non-trivial example of an energy density (of a material) for which a large class of local minimizers could be made explicit. In particular, we exhibited a material model, where the absence of metastability coexists with a wild nonuniqueness of global minimizers.

Metastability in elastostatics, understood as the existence of strong local but not global minimizers, is an important problem in nonlinear elasticity because of the ubiquity of hysteresis in martensitic phase transitions. While in Neumann problem the existence of such local minimizers has been linked to the incompatibility of the energy wells, here we showed that in Dirichlet problem the incompatibility of the wells is not sufficient for metastability. More specifically, we presented an example of the energy density of a hyperelastic material with non-rank-one convex, double well energy, for which we could prove the lack of strong local minimizers which are not global on any domain and for any Dirichlet boundary conditions. The analytical transparency of our arguments was due to the choice of the energy density which is geometrically linear and isotropic.

For the same material model we could fully characterize the necessary and sufficient conditions for global energy minima in hard device loading, which allowed us to discover the multiplicity of global minimizers with nonuniqueness unrelated to either objectivity or material symmetry. As an important element of this analysis we used a novel way of expressing the energy of equilibrium configurations as boundary integrals which can be viewed as a nontrivial nonlinear generalization of the classical Clapeyron theorem. First of all, we showed that the quasiconvex envelope in our model can be achieved by simple lamination which can be characterized explicitly. The problem with this construction is that in real physical situations surface energy plays the role of the selection mechanism of the otherwise energetically equivalent configurations [START_REF] Klouček | Computational modeling of the martensitic transformation with surface energy[END_REF][START_REF] Kohn | Surface energy and microstructure in coherent phase transitions[END_REF][START_REF] Dolzmann | The influence of surface energy on stress-free microstructures in shape memory alloys[END_REF], ruling out such well-known optimal microgeometrs as laminates and coated sphere constructions. Therefore, of particular interest in physical applications are global minimizers with finite surface area, and we exhibited their existence for our material model.

An interesting aspect of our highly degenerate problem with elastically indistinguishable phases is a sensitivity of the morphology of the interface between the phases to the global shape of the transforming body. To illustrate this effect, we explicitly computed a one-parameter family of non-affine energy minimizing configurations for the case of a finite domain with piecewise-smooth boundary. Our analysis of this example suggests that despite the absence of metastability and the associated constitutive hysteresis in this model, the direct and the reverse transformation may follow different morphological paths while traversing energetically equivalent configurations.

Acknowledgments. YG was supported by the National Science Foundation under Grant No. DMS-2005538. The work of LT was supported by the French grant ANR-10-IDEX-0001-02 PSL. YG is grateful to Pavel Etingof who pointed out the references for the proof of the infinity of the surface area of Hashin's concentric sphere construction in 1992.

A The degeneracy of acoustic tensors of rank-1 envelopes

Here we prove that the acoustic tensor of a rank-one convex envelope RW of the non rank-one convex energy must have a degenerate direction at all points F 0 where RW (F 0 ) < W (F 0 ). We also show that in particularly simple situations, such as the one discussed in this paper, this property may be even sufficient to compute the whole rank-one convex envelope. We recall that the acoustic tensor of the energy W (F ) at F = F 0 in the direction n is a quadratic form A(n) defined by

A(n)a • a = W F F (F 0 )(a ⊗ n), a ⊗ n .
Theorem A.1. Let F 0 be fixed and suppose RW (F 0 ) < W (F 0 ), where RW denotes the rank-1 convexification of W . Assume further that W and RW are C 2 near F 0 . Let A 0 (n) be the acoustic tensor of RW at F 0 . Then there is a direction n such that det A 0 (n) = 0.

Proof. Assume that there is no such direction n. Since RW is necessarily rank-1 convex we can conclude that for any direction n the matrix A 0 (n) is positive semidefinite. Since the function n → A 0 (n) is continuous there exists a positive number α such that for every unit vector n A 0 (n) ≥ αI.

Since RW ∈ C 2 near F 0 there exists a number δ > 0 such that for every F satisfying |F -F 0 | < δ the following inequalities hold: where φ is a smooth nonnegative function supported on the unit ball in M and such that φ(0) = 1. We can choose so small that the following inequalities hold:

1. φ L ∞ < 1 2 β, 2. √ ∇ F ∇ F φ L ∞ ≤ 1 4 α,
3. supp φ( F -F 0 4 √ ) ⊂ B(F 0 , δ), where B(F 0 , δ) is the ball of radius δ around F 0 in M.

Then T (F ) ≤ W (F ) and the acoustic tensor A (n) of T (F ) satisfies

A (n) ≥ 1 4 αI
in the sense of quadratic forms. Thus T (F ) is rank-1 convex and T (F ) ≤ W (F ) but T (F 0 ) > RW (F 0 ). Contradiction. Thus our assumption is false and at every point F 0 where RW (F 0 ) < W (F 0 ) there is a direction n such that the acoustic tensor A 0 (n) is degenerate.

Remark A.2. It follows that materials that transform by forming microstructures with sharp phase boundaries to accommodate deformations produced by the propagation of sound waves will have a direction with a zero sound speed.

We can apply Theorem A.1 to the energy (1.3). Since our material is isotropic, the degeneration of the acoustic tensor may be either through µ = 0 or through λ + 2µ = 0 with the latter also meaning that the bulk modulus κ = -(4/3)µ. The possibility that µ = 0 is excluded because the tangential shear modulus is the same at every deformation.

Suppose that we have somehow guessed that if

W (H) = f (Tr ε) + µ|dev(ε)| 2 , then RW (H) = F (Tr ε) + µ|dev(ε)| 2
for some function F , yet to be determined. In that case Theorem A.1 will let us determine the function F (θ). It is easy to compute that for any unit vector n and any vector a ∈ R 3

RW HH (H)(n ⊗ a), n ⊗ a = F (Tr ε)(a, n)

Figure 1 :

 1 Figure 1: Double-well nonlinearity in a geometrically linear bi-quadratic Hadamard material.

Theorem 3 . 2 .

 32 Let Ω ⊂ R d be a star-shaped domain with a Lipschitz boundary. Assume that (i) y : Ω → R m is Lipschitz continuous, (ii) y solves (3.2) and (3.3) in the sense of distributions, (iii) ∇y(x) is locally stable for a.e. x ∈ Ω, (iv) y(x) = F x for all x ∈ ∂Ω, (v) There exists > 0, such that y ∈ C 1 (Ω ; R m ), where Ω = {x ∈ Ω : dist(x, ∂Ω) < . Then y(x) is the global minimizer in (3.1), i.e. Ω W (∇y)dx = |Ω|QW (F ).

∂ΩP

  * n • xdS = Ω P * , ∇x dx = Ω Tr (P * )dx = d Ω W (∇y)dx -Ω P , ∇y dx. Equation (3.5) follows from Ω P , ∇y dx = ∂Ω P n • ydS.

Figure 3 :

 3 Figure 3: A coherent precipitate in the square.

  A(n) is the acoustic tensor of RW at F . 2. W (F ) -RW (F ) ≥ 1 2 (W (F 0 ) -RW (F 0 )) = β > 0.Now consider a function T (F ) = RW (F ) + φ( F -

  Theorem 2.5. Let W (H) be the energy density satisfying(2.20). Then any Lipschitz strong local minimizer of Ω W (∇u)dx with prescribed Dirichlet boundary conditions is a global minimizer.

  while there are Lipschitz configurations satisfying (3.2), but not (3.3). For instance, the Maxwell relation, which is the last equation in (2.17), is a consequence of (3.3), but not of (3.2).

	Definition 3.1. A configuration y ∈ W 1,∞ (Ω; R m ) is called a stationary (or an equilibrium)
	state if it satisfies both (3.2) and (3.3).

This is never the case for the potential (1.4).

We can say that this is a conjecture, as we have no formal proof of this statement.

The function y(θ) is affine because the property of F (θ) can be written as y (θ) = 0.

Therefore, we get a formula for the acoustic tensor A(n).

We see that det A(n) = µ 2 (F (Tr ε) + 4µ/3) for all directions n. Thus, for all H for which RW (H) < W (H) we get F (Tr H) = -4µ/3. The continuity of RW H implies that at the boundary points θ 1 and θ 2 of the binodal region we have Φ (θ 1 ) = Φ (θ 2 ). Further, the continuity of W implies that the affine function 3 y(θ) = F (θ)+2µθ 2 /3 would be the equation of the common tangent to the graph of Φ(θ). Thus, we obtain that F (θ) + 2µθ 2 /3 = Φ * * (θ), and we recover the rank-one convex envelope of W (H), which in this case is seen to coincide with its quasiconvexification (2.4).


Trevor J Jones

Etienne Jambon-Puillet

Joel Marthelot

P.-T Brun

Bubble Casting Soft Robotics

fully rationalize the fluid mechanics at play during the 50 bubble injection and the drainage to predict the shape 51 of the actuators' cross-section. We then turn to the elas-52 tic problem to elucidate how the bending motion upon 53 inflation is determined by the shape previously sculpted.

54

We leverage these quantitative results to design actua-55 tors programmed for specific tasks through bubble cast-56 ing such as soft robotic muscles with tunable strength 57 capable of lifting objects (Fig. 1b and Movie S2) and soft 58 fingers with sequential actuation from a single pressure 59 source (Fig. 1c and Movie S3). Additionally, controlled 60 folding in three dimension is demonstrated using curvi-61 linear actuators, either free (Fig. 1d) or attached to thin 62 membranes (Fig. 1e-f). Finally, we show that bending in 63 different directions can be achieved by rotating portions 64 of the mold relative to each other around the gelation 65 point of the polymer [29] (Fig. 1g and Movie S4).

66

Inspired by living organisms, soft robots are developed from intrinsically compliant materials enabling continuous motions mimicking animal and vegetal movement [START_REF] George | Soft robotics[END_REF]. In soft robots, the canonical hinges and bolts are replaced by elastomers assembled into actuators programmed to change shape following the application of stimuli, e.g. pneumatic inflation [START_REF] Polygerinos | Soft 319 robotic glove for combined assistance and at-home rehabilitation[END_REF][START_REF] Robert F Shepherd | Multigait soft robot[END_REF][START_REF] Ramses V Martinez | Robotic tentacles with three-dimensional mobility based on flexible elastomers[END_REF][START_REF] Roche | A bioinspired soft actuated material[END_REF]. The morphing information is typically directly embedded within the shape of these actuators, whose assembly is facilitated by recent advances in rapid prototyping techniques [START_REF] Mosadegh | Pneumatic networks for soft robotics that actuate rapidly[END_REF][START_REF] Hong | High-force soft printable pneumatics for soft robotic applications[END_REF][START_REF] Wehner | An integrated design and fabrication strategy for entirely soft, autonomous robots[END_REF][START_REF] Zhao | Scalable manufacturing of high force wearable soft actuators[END_REF][START_REF] Paek | Microrobotic tentacles with spiral bending capability based on shape-engineered elastomeric microtubes[END_REF][START_REF] Jin | Kirigami-inspired inflatables with programmable shapes[END_REF]. Yet, these manufacturing processes have limitations in scalability, design flexibility and robustness. Here we demonstrate a new all-in-one methodology for the fabrication and the programming of soft machines. Instead of relying on the assembly of individual parts, our approach harnesses interfacial flows in elastomers that progressively cure to robustly produce monolithic pneumatic actuators whose shape can easily be tailored to suit applications ranging from artificial muscles to grippers. We rationalize the fluid mechanics at play in the assembly of our actuators and model their subsequent morphing. We leverage this quantitative knowledge to program these soft machines and produce complex functionalities, e.g. sequential motion obtained from a monotonic stimulus. We expect that the flexibility, robustness and predictive nature of our methodology will accelerate the proliferation of soft robotics by enabling the assembly of complex actuators, e.g. long, tortuous or vascular structures, thereby paving the way towards new functionalities stemming from geometric and material nonlinearities.

Soft robots can achieve complex tasks such as gentle 1 gripping, crawling, or swimming [START_REF] Robert F Shepherd | Multigait soft robot[END_REF][START_REF] Ilievski | Soft 368 robotics for chemists[END_REF][START_REF] Andrew D Marchese | Autonomous soft robotic fish capable of escape 372 maneuvers using fluidic elastomer actuators[END_REF] using low-2 complexity, muscle-like soft actuators that bend, twist, 3 contract, or elongate on demand [START_REF] Polygerinos | Soft 319 robotic glove for combined assistance and at-home rehabilitation[END_REF][START_REF] Roche | A bioinspired soft actuated material[END_REF][START_REF] Yang | Buckling pneumatic linear 378 actuators inspired by muscle[END_REF][START_REF] Johannes Tb Overvelde | Amplifying the re-382 sponse of soft actuators by harnessing snap-through 383 instabilities[END_REF][START_REF] Hawkes | A soft robot that 387 navigates its environment through growth[END_REF]. This niques [START_REF] Zhao | Scalable manufacturing of high force wearable soft actuators[END_REF][START_REF] Paek | Microrobotic tentacles with spiral bending capability based on shape-engineered elastomeric microtubes[END_REF][START_REF] Jin | Kirigami-inspired inflatables with programmable shapes[END_REF] are restricted to simple geometries, while freeform fabrication techniques [START_REF] Hong | High-force soft printable pneumatics for soft robotic applications[END_REF][START_REF] Wehner | An integrated design and fabrication strategy for entirely soft, autonomous robots[END_REF] typically lack scalability and require long printing times. Often, the inflation of these actuators is hard to predict such that trial and error or long simulations are required to tailor the shape of the actuator for specific applications. Additionaly, soft pneumatic robots assembled from these actuators need to perform complex or sequential motion which typically require several actuators with independent fluid sources [START_REF] Vasios | Harnessing viscous flow to simplify the actuation of fluidic soft robots[END_REF].

Here we introduce bubble casting, a simple and versatile fabrication methodology to assemble monolithic actuators programmed using the rules and tools of fluid mechanics. Fig. 1a and Movie S1 illustrate our approach involving elastomers that spontaneously cure when reagents are mixed (see Methods). We first fill a tubular mold by injecting uncured elastomer melt. While the melt is still liquid, we inject air to form an elongated bubble that creates the inner void of the actuator (Fig. 1ai). Gravity then sculpts the actuator by draining the polymer film and allowing the bubble to rise (Fig. 1aii). Eventually, as the melt solidifies this shape is frozen and the actuator can be readily used when de-molded, e.g. as a gripper (Fig. 1aiii). In the following we demonstrate how to control bubble casting to A bubble is injected (i) in a mold previously flooded with polymer melt. The residue of polymer drains (ii) and cure to provide a anisotropic actuator that can be readily used when de-molded, e.g. as a gripper (iii) (scale-bar: 1 cm). (b) Contractile coiling of long actuators yielding: (i) muscle-like contraction and (ii) linear translation (scale-bars: 5cm). (c) The bubble casting dynamics is used to program the deformation of an actuator, i.e. displaying sequential flexion of the four labeled digits following inflation via a simple pressure ramp (scale-bar: 1 cm). Morphing dynamics of (d) a curvilinear actuator (experiments and numerical simulation) submerged in a density matched fluid (scale-bar: 2 cm). (e) (i) Spiral shaped actuator and (ii) cross-like actuator, both attached to a thin membrane (scale-bars: 2 cm). (f) "Fish tail" motion is obtained using bubble-casting in a branched network (scale-bar: 1 cm). (g) Actuator displaying curvature of equal magnitude but opposite signs obtained by rotating the mold at the gelation point of the polymer (scale-bar: 1 cm) (see Movie S4). achieve versatile and programmable actuators. We first The shape of the cross section of our actuators is imparted by a two-fold process (Fig. 1a). First, the injection of the air bubble in the uncured melt results in the deposition of a thin polymer annulus in the mold, then the concomitant gravity-driven drainage and curing of the elastomer sculpt the final shape. We first consider the shape of the cross-section prior to drainage, right after the bubble injection. The flow at the front of the advancing bubble (Fig. 2a) leaves an annulus of constant thickness h i . Fig. 2b shows this thickness in dimensionless form h i /R as a function of the dimensionless bubble speed or capillary number Ca = µU/γ, with R the radius, µ the melt viscosity, U the bubble velocity, and γ the melt surface tension. The data for different polymers and tube radii fall on a master curve showing that the thickness first increases with the velocity until it eventually saturates for Ca 3. ing viscous and capillary forces in the advancing menis-85 cus [START_REF] Patton | The motion of long bubbles in tubes[END_REF][START_REF] Aussillous | Quick depo-454 sition of a fluid on the wall of a tube[END_REF][START_REF] Klaseboer | An extended bretherton model for 458 long taylor bubbles at moderate capillary numbers[END_REF] (see SI S1) yielding:

This master curve can be understood from balanc-

86 h i R = 1.34Ca 2/3 1 + 1.34βCa 2/3 . (1)
We plot Eq. 1 with β = 2.65 along our data in Fig. 2b 87 and find a good agreement. The melt annulus thickness, 88 and hence the void fraction of the actuator, can thus be 89 tuned by controlling the bubble velocity. In particular, 90 using Ca 3 robustly results in h i ∼ 0.3R.

91

The drainage of the annular polymer film left after 92 bubble injection eventually sculpts the cross-section of 93 our actuators shown in Fig. 2a. As evident from Fig. 2c 94 which shows the thickness along the central angle ψ de-95 fined in Fig. 2a, the final shape consists of an upper thin film (ψ π/4) of quasi uniform thickness connected to a thicker region at the bottom.

We first consider the upper thin film drainage driven by gravity g and resisted by viscosity [START_REF] Takagi | Flow and 461 instability of thin films on a cylinder and sphere[END_REF]. To predict the final thickness, we must account for the time varying viscosity of our melt µ(t) as it solidifies. Our rheological measurements (see SI S3.1) are well described by the function µ(t) = µ 0 (1 -t/τ c ) -n which diverges at t = τ c the curing time. The values of the initial viscosity µ 0 , the curing time τ c and the exponent and n ≃ 2 are fitted from the data, while the density of the melt ρ is assumed constant (see Table S1). Accounting that the drainage starts at a time t = τ w after the reagents are mixed and finish at t = τ c , we include the viscosity variations in our drainage model [START_REF] Lee | Fabrication of slen-465 der elastic shells by the coating of curved surfaces[END_REF] (see SI S2.4) and predict the final thickness to be:

h(τ c ) = 3µ 0 R(n + 1)τ n+1 c 2ρgτ c (τ c -τ w ) n+1 . (2)
Note that the final film thickness does not depend on the initial film thickness.

We then consider the lower thick part of the crosssections in Fig. 2a. The shape of this bath results from the competition between gravity and capillarity, yielding the Young-Laplace equation that we solve numerically under the proper constraints (see SI S2.2). Finally, we match the upper film and lower bath to obtain the final thickness profile along the whole cross-section (see SI S2.5). Our model produces the curves overlaid on images of the cross-section in Fig. 2a and the theoretical thickness profile plotted along typical experimental data on Fig. 2c. The agreement is favorable, demonstrating that we have accurately modeled the flow physics of our method.

Note that bubble casting does not require external control so that fluid mechanics alone dictates the actuator shape, thereby making this methodology extremely robust. While the membrane thickness is the order 100 µm, we obtain meter long samples that are virtually uniform (see Fig. 2d). Importantly, we can continuously and predictably tune the membrane average thickness h f through the waiting time τ w since

h f ∼ (τ c / (τ c -τ w )) n+1 2
(see Eq. 2). In turn, the membrane thickness will play a key role in the actuator response to inflation as detailed next.

Upon inflation, the upper membrane of our bubble casted actuators stretches significantly more than its lower part, thereby creating a torque that bends the actuator (Fig. 3a). In absence of other forces, the ac-142 tuator bends on its whole length and adopts a uniform the lower region of the actuator is virtually undeformed while the upper thin membrane is under quasi isotropic stretch (Fig. S5a). This observation allows us to derive a simplified theoretical treatment of the problem (see SI S4.3), which combined with dimensional analysis suggests that the rescaled curvature κR 2 /h f only depends on the rescaled pressure P R/Gh f (see SI S4.4).

In Fig. 3b we show that following this rescaling provides a reasonable collapse of our data. The resulting master curve exhibits a power law in the useful range of curvature (κ/R > 10 -2 ) that we fit and recast to obtain the relationship:

κR ∼ P G R h f α-1 α α , (3)
with α = 3.72 ± 0.07.

Further, we model our actuators as elastic rods [START_REF] Audoly | Elasticity and 468 geometry: from hair curls to the non-linear response 469 of shells[END_REF] with effective natural curvature κ varying according to Eq. 3 when pressure is applied. We solve this centerline based model to predict the deformation of the actuators in the presence of external loads (see SI S5.2), e.g. when blocked from bending by a wall (Fig. 3c). We find that the overall shape of the actuator and the blocking force F (P ) exerted in this setup are well captured by our reduced order model when varying independently the actuator bending stiffness B, shear modulus G, length L, radius R, and thickness h f (see Fig. 3d), thereby validating the approximations introduced in our model. Importantly, this model allows us to program the mechanical response of our actuators via the shape of their cross-section as detailed next. In our framework, solving the inverse design problem, i.e. finding the actuator that will morph to a target curvature, merely requires to invert the scaling law in Eq. 3. Now that we have modeled how our actuators bend and behave as elastic rods, we demonstrate how to leverage this knowledge to design soft robotic muscles inspired by cucumber tendrils [START_REF] Gerbode | How the cu-472 cumber tendril coils and overwinds[END_REF]. These muscles can be used over a wide range of size and strength (Fig. 4a,b) and can be integrated to programmable soft machines (Fig. 4c).

When subject to gravity, the curvature imparted by the actuation pressure competes with the actuator weight and generates a curling motion (see Fig. 4a). These shapes are analogous to the static configurations obtained with naturally curved elastic rods [START_REF] Miller | Shapes of a suspended curly hair[END_REF] and are well captured by our centerline-based theory. In Fig. 4d we show that our robotic muscles can contract up to 1 -ℓ/L ≈ 80% of their initial length before entering in self-contact, exceeding the stroke of typical pneumatic Programming can also be achieved using fluid mechanics alone, i.e. modulating the membrane thickness h f via eq. 2. For instance, using a programmed sequence of injection of the bubble where τ w increases by step across a sample, we readily cast an actuator with four different values of h f , each confined to a "finger" of the actuator (see Fig. 1c). As a result, an increase of the pressure in this monolithic actuator leads to the sequential bending of each finger shown in Fig. 1c and Movie S3. Likewise, we can solve the inverse problem where a simple target shape can be obtained by carefully choosing τ w across the actuator (see Movie S6).

In closing, we note that Bubble casting is a bondingfree fabrication method that relies on the fluid flow instead of an internal template to build a void. This reliance on continuum mechanics as opposed to machined parts allows for defect free fabrication of soft pneumatic actuators over a wide range of sizes (see Fig. 1 and4) and previously impossible aspect ratios. In particular, the thin membranes we achieve optimize the eccentric void topology and outperform most inflatable actuators in terms of bending coefficients (defined as the ratio between the curvature and the applied pressure at a bending deformation of 90 o [START_REF] Gorissen | Elastic inflatable actuators for soft robotic applications[END_REF]). Further, the unconstrained internal surface makes the construction of curved and network actuators that achieve 3D folding shapes (Fig. 1d-f and Movie S4) as simple as building straight actuators. These new capabilities will resonate in the soft matter community, and, in particular could lead to new generation robotic materials suited to move and interact with their environment while maintaining tractable complexity. More generally, our methodology falls under the category of approaches leveraging out-ofequilibrium fluidic processes to revolutionize our ability to build structures, e.g. relying on instabilities [START_REF] Kaufman | Structured spheres 485 generated by an in-fibre fluid instability[END_REF][START_REF] Marthelot | Designing soft materials with interfacial in-489 stabilities in liquid films[END_REF], shear-flow [START_REF] Zhao | Layered nanocomposites by 495 shear-flow-induced alignment of nanosheets[END_REF], Marangoni effects [START_REF] Bryan A Nerger | Marangoni flows drive the alignment of fibrillar cellladen hydrogels[END_REF] and centrifugal forces [START_REF] Cybulski | Concentric liquid reactors for chemical synthesis and separation[END_REF].

Methods

To make our actuators we used vinyl polysiloxane (VPS), Zhermack elite double 8, 16, and 32 silicone elastomers (G = 72.4, 172, 368 kPa respectively). The catalyst and prepolymer base were mixed in a 1:1 weight ratio in a centrifugal mixer for 10 s at 2,000 rpm (clockwise) and 10 s at 2,200 rpm (counterclockwise) at room temperature. The resulting polymer melt (rheology described in SI S3.1) was injected into a mold where we then followed our bubble casting method, as described in the text. For the molds we used cellulose tubes (R = 1.6, 2.4, 3.2, 6.4 mm), glass tubes (R = 0.5 mm), and cast acrylic sheets carved with a cnc milling tool on one side. Schematics for our experiments are shown in Fig. S4.

The bubble velocity, curvature data, and free-hang length was recorded using a camera and image processing in Python. Since the final drainage loses memory of the initial condition (as described in the text), we often injected the bubble early and reset τ w by flipping the mold. This allowed us to maintain constant drainage timing along the length of long actuators. To attach the thin membranes shown in Fig. 1e,f, the two sides of the carved molds (see Fig. 1a) were separated. Some more polymer was then spin-coated on top of the actuator still encased in the mold. We realized the experiment show in Fig. 1g by joining two tubes using a 3D printed tube connector. When the polymer is at the 280 gelation point one tube is carefully rotated relative to 281 the other (see Fig. S8). The cross-section profiles were 282 obtained post curing using an optical microscope and 283 ImageJ. Pressure data was measured using a differen-284 tial pressure sensor (MPX5100dp) connected to a data 285 acquisition board (Arduino). To measure the pulling 286 force, we attached the actuator to an Instron and per-287 formed tensile tests. The blocking force was measured 288 by placing the free end of the actuator on a plate fit-289 ted to a load cell. The curvature of our actuators was 290 measured without external forces (e.g. gravity) by float-291 ing the elastomer on a water bath. Actuator mechanics 292 were measured quasi-statically by slowly injecting air in 293 1 s pulses and waiting 4 s between pulses. All actuators 294 were pierced with syringe needles and inflated by a sy-295 ringe mounted on a syringe pump (Harvard Apparatus). 296

4 unique combination of softness and bioinspired motion 5 make

 45 soft robots appealing for a variety of innovative ap-6 plications where rigid robots would fail[START_REF] Majidi | Soft robotics: a perspec-390 tive-current trends and prospects for the future[END_REF]. This bloom-7ing field is fueled by recent insights in modeling, compu-8 tations and manufacturing that enable the design, the 9 programming and the assembly of various kinds of soft 10 machines. While chemically, thermally, electrically or 11 magnetically activated soft actuators have been demon-12 strated [18-24], silicone bodied pneumatic robots pow-13 ered by pressurized voids have drawn considerable atten-14 tion owing to their simple and rapid actuation [25, 26]. 15 The kinematics of such robots is encoded in the actu-16 ator flesh, i.e. shape or material, such that a change 17 of internal pressure is mechanically converted into spe-18 cific motion [27]. Manufacturing soft pneumatic actu-19 ators, particularly the void, is non-trivial and is usu-20 ally accomplished with sequential molding procedures 21 and removable frameworks that are tailored for specific 22 actuators. Likewise, state-of-the-art film coating tech-23

Figure 1 :

 1 Figure 1: From flow to programmed actuation. (a) Schematic of the bubble Casting fabrication methodology. A bubble is injected (i) in a mold previously flooded with polymer melt. The residue of polymer drains (ii) and cure to provide a anisotropic actuator that can be readily used when de-molded, e.g. as a gripper (iii) (scale-bar: 1 cm). (b) Contractile coiling of long actuators yielding: (i) muscle-like contraction and (ii) linear translation (scale-bars: 5cm). (c) The bubble casting dynamics is used to program the deformation of an actuator, i.e. displaying sequential flexion of the four labeled digits following inflation via a simple pressure ramp (scale-bar: 1 cm). Morphing dynamics of (d) a curvilinear actuator (experiments and numerical simulation) submerged in a density matched fluid (scale-bar: 2 cm). (e) (i) Spiral shaped actuator and (ii) cross-like actuator, both attached to a thin membrane (scale-bars: 2 cm). (f) "Fish tail" motion is obtained using bubble-casting in a branched network (scale-bar: 1 cm). (g) Actuator displaying curvature of equal magnitude but opposite signs obtained by rotating the mold at the gelation point of the polymer (scale-bar: 1 cm) (see Movie S4).

Figure 2 :

 2 Figure 2: (a) Schematic of the Bretherton problem and photographs of the subsequent final cross-section post drainage. VPS-32 actuators with working time τ w = 300s and radius R = {3.2, 2.4, 1.6} mm, the red lines show our prediction (See SI S2). (b) Dimensionless Bretherton film thickness h i /R plotted against the Capillary Number Ca = µ 0 U/γ. The solid line is theory (Eq. 1) with β = 2.65. (c) Film thickness of the samples in (a) plotted against the central angle ψ. The red lines are theoretical predictions (identical to a). (d) Average thickness of the membrane h f over a meter long sample (R=2.4mm). The black line shows our prediction while the green band accounts for parameters uncertainty.

Figure 3 :

 3 Figure 3: (a) Pictures of a VPS-16 actuator at rest and inflated at constant pressure: the actuator bends with a uniform curvature κ (Scale-bar: 1 cm). Inset: actuator curvature κ as a function of applied pressure P as the thickness is independently varied: h f = {130, 165, 273} µm. (b) Rescaled dimensionless curvature κR 2 /h f plotted against the rescaled pressure P R/Gh f . The solid blue lines correspond to the energy minimization theory Eq. S28 for the relevant range of parameters (see SI S4.3). The dashed line is the scaling law Eq. 3 with a prefactor of 0.52. The gray shaded region represents the propagated experimental uncertainty. (c) Pictures of VPS-08 at rest and inflated in a blocking force configuration (Scale-bar: 1 cm). The force F is measured and reported in (d). The black (resp. reddotted) line shows the linear (nonlinear) beam theory prediction (see SI S5.2).

143 curvature κ

  as shown in Fig.3a. The pressure required 144 to bend the actuator is independent of its length but 145 increases with the elastic shear modulus G, outer ra-146 dius R and most importantly with the membrane thick-147 ness h f (Fig.3aand Fig .S6c). Predicting the defor-148 mation of an actuator is amenable to 3D finite element 149 simulations (see SI S4.2), which further confirm that 150

Figure 4 :

 4 Figure 4: (a) Contractile coiling of an actuator formed in a rectilinear mold. Experiments and model are shown side-by-side (see SI S5; scale-bar: 5 cm). (b) Pictures of two VPS-16 robotic muscles of R = 6.4 and 0.5 mm lifting an empty bottle and a paperclip (scale-bars: 4 cm). (c) Pictures of a soft machine gripping and lifting a ball from a narrow cylinder (scale-bar: 5 cm). (d) Contraction ℓ/L (see a) as a function of the applied pressure for actuators with different apex thicknesses h f = {0.14, 0.17, 0.22} mm. (e) Same data shown as a function of the dimensionless pressure. The red curve is the result of Kirchhoff rod simulations (see SI S5). VPS-16, shear modulus G = 172 kPa, outer radius R = 1.6 mm and length L = 320 mm. (f) Force-elongation curves of a typical actuator at different pressure. The solid curves are the result of the equivalent Kirchhoff rod simulations (see SI S5).

Acknowledgments 297 E.J.-P. was partially supported by NSF through the 298 Princeton University Materials Research Science and 299 Engineering Center (Grant DMR-1420541) and TJJ was 300 partially supported by NSF CAREER award (CBET 301 2042930)

Author Contributions 303 T.J.J., J.M., and P.-T.B. conceived the project. T.J.J., 304 E.J.-P., and P.-T.B. conducted the experiments and an-305 alyzed the data. T.J.J. performed the Kirchhoff rod 306 simulations. E.J.-P. performed the finite element sim-307 ulations. All authors wrote the manuscript.

308

Data availability 309

The authors declare that the main data supporting the 310 findings of this study are available within the article 311 and its Supplementary Information files. Extra data are 312 available from the corresponding author upon request. 313


Éric Darrouzet

Laurence Berville

Armel Laot

Marie Legeais

Anne-Geneviève Bagnères

Building behaviour analyzed by X-rays tomography in two subterranean termite species

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Building behaviour analyzed by X-rays tomography in two subterranean termite species.

E. Darrouzet 1 , L. Berville 1 , A. Laot 2 , M. Legeais 2 & A-G Bagnères 1

1 University of Tours, IRBI, Parc de Grandmont, 37200 Tours, France.

2 CHRU of Tours 2 bd Tonnellé 37044 Tours France 2 CHRU of Tours 2 bd Tonnellé, 37044 Tours, France. (eric.darrouzet@univ-tours.fr).

Reticulitermes species are subterranean termites which are found in natural and urban area in Europe. Because of their cryptic life, it is difficult to analyze colony sizes and constructions elaborated in soil and/or pieces of wood, i.e., numbers of chambers, structures elaborated in mud inside the wood, etc. By using X-rays tomography, a non invasive medical imaging technique, we were able to perform an original study on these different parameters in the lab.

We analyzed wood consumption strategies (volumes of wood eaten, kinetic of consumption, how termites colonized the piece of wood)

3D view of chambers in the wood.

and building behaviour by two termite species, Reticulitermes grassei and R. santonensis.

- year.

-In our conditions, the two species present great behavioural differences, as R. santonensis (1) eats more wood than R. grassei and ( 2) is the one who elaborates mud structures inside the wood.

-Our results confirm previous observations than R. santonensis, an invasive species introduced in France 300 years ago, makes more damages than R. grassei.

R. santonensis

workers built mud structures in chambers.

  5 different micro-colonies are analyzed for each species -900 workers/colony in each plastic box (18 x 12 x 7 cm) -a piece of wood (pinus: 12 x 4 x 4 cm) -2 mm of sand in the bottom consumes more wood than R. grassei (Mann-Whitney: U=440261, P<0.001) during 258 days.Different structures built in the woodVolume in cm 3 Number of stThe number of structures built in chambers is correlated to the volume of wood consumed (Pearson correlationrays tomography is a fantastic tool to study insects' architectures and building behaviour in continue.-For the first time, wood consumption and building behaviour of living subterranean termites were observed in situ during one year. r=0.846, P=9.88x10 6 ). During the experiment, only pillars Different mud structures are built inside the chambers in the wood: caps, pillars, floors and walls.