Catherine Larzul
email: catherine.larzul@jouy.inra.fr

PIG GENETICS, INSIGHT IN MINIPIGS

Keywords: Genetic, minipig, selection

The sequence of the pig has now been published. Associated with genomic tools such as high density genotyping, (re)sequencing, comparative genetic, it offers new opportunity to understand the genetic basis of complex traits. The challenge is now to get adequate phenotypes. For minipigs, focus is made on traits relevant for medicine, in order to use minipigs as a biomedical model. Because of the number of research fields, specific lines had to be developed and within the last decades, several minipig lines were made available. For example, to specifically study development and regression of cutaneous melanoma, the MeLiM line was built up to largely express this particular trait. It was then possible to study the genetic determinism of this trait and acquire knowledge relevant for human disease, and in a lesser extent also relevant for pig production. Currently, one of our aims is to develop a minipig line adequate for studies on human nutrition related diseases, such as obesity or atherosclerosis. This approach first requires the characterization of minipigs to understand the genetic determinism of their behavior but also some physiological parameters related to human nutrition.

INTRODUCTION

Swine breeding has a long history. The process of domestication, which begun circa 7000 b.c. includes selective breeding for specific characteristics, and may be considered a first step in genetic improvement. The process of domestication includes selective breeding for specific characteristics, and may be considered a first step in genetic improvement. Pig improvement has first been the result of empirical methods employed by individual farmers. Then selection became a technique scientifically based, progressing with the knowledge of the biology involved, and particularly with the development of genetics. Through centuries and geographical areas, the variety of breeding goals has conducted to a large number of pig breeds. Crossbreeding has also been used to create synthetic lines. Those breeds have been mainly developed for food subsistence. In the mid 20 th century, pig has been viewed as an alternative to other animals for human oriented studies [START_REF] Bustad | Use of pigs in biomedical research[END_REF]. Most of the strains derived from a small number of reproducers taken from local populations of small size (e.g. Ossabaw, Lan Yu, Yucatan) and raised in a controlled environment, others were successively crossed with Asian populations and/or European populations for increasing diversity, introducing new coat colours or desirable characteristics. Some of the strains were then selected, generally for small size but also for other specific traits. Minipig strains were created under laboratory conditions in the USA, Russia, Japan, Taïwan, Germany, Czech Republic, Australia, sometimes derived from one another. In 1996, Thikonov enumerated 20 breeds of miniature pigs used for biomedical and biological studies. Some strains were constituted for general research purposes (e.g. Hormel, Göttingen) and other for very specific studies (e.g. Yucatan strains, MeliM). Recently, [START_REF] Köhn | History and development of miniature, micro-and minipigs[END_REF] described the constitution and major characteristics of nine micro-or minipig breeds including Clawn, Göttingen, Hanford, Troll, Sinclair, Westran and 3 Yucatan breeds. He also offered a schematic overview of minipig breeds used in biomedical research.

Advances in molecular biology as well as in statistics and computational facilities offer a new challenge for pig breeding and genetics in the next few years, which might increase the interest of minipig for biomedical research.

CONSTITUTION OF SOME MINIPIG STRAINS

A description of the first small size pig strains developed for biomedical and pharmaceutical research purpose is given by [START_REF] Bustad | Miniature Swine: development, management, and utilization[END_REF]. The very first project was initiated in 1949 at the Hormel Institute of the University of Minnesota, the strain between refered as Minnesota Miniatures or Hormel Miniatures. As detailed by [START_REF] Dettmers | Minnesota's Miniature pigs[END_REF], it was started the crossbreeding of four American breeds; first 4 Guinea hogs females were mated to one wild boar from Catalina Island, then two and four years later, Piney Woods breeders were introduced and in 1957, the fourth breed stock was constituted of Ras-n-lansa (Guam pigs), introduced to reduce size. In total, only 4 males and 11 females were used from original breeds. The Hormel minipig strain was further on selected for small size (weight at 140 days of age) (Dettmers et al., 1965). In order to develop a white color strain, a Yorkshire boar was mated to Minnesota sows in 1963. As underlined by Simianer and Khön (2010), two founder breeds of the Minnesota Minipig are island-based strains, making use of the evolutionary phenomenon of so-called island dwarfism, reflecting the observation that isolated populations on small islands tend to develop dwarf variants in many mammalian species. Sinclair pigs derived directly from Hormel pigs.

Pigs from Florida swamps were bred by Pitman-Moore pharmaceutical company, which gave the name of the strain. Breeders were distributed to other laboratories to found new minipig strains (Hanford, VitaVet, Labco). The Hanford miniature swine was obtained by mating Pitman-Moore breeders with Palouse swine and a Louisiana swamp hog and later with introduction of Mexican pigs (Yucatan pigs) in their breeding stock.

The Berlin Miniature Pig breed (also known as Mini-Lewe) has been developed 20 years ago at the experimental Station of the Humboldt-University from a cross of Vietnamese Pot Belly Pigs, Saddle Back Pigs and German Landrace followed by consistent selection for small body weight [START_REF] Hardge | Maternal influences on litter size and growth in reciprocal crossed Miniature Pigs and Durocs[END_REF].

A detailed description of one the minipig breeds mostly used for pharmaceutical and biomedical research, the Göttingen minipig, was recently provided by [START_REF] Simianer | Genetic management of the Göttingen Minipig population[END_REF]. The development of the Göttingen minipig was initiated in the 1960's at the Friedland Experiment Station of the University of Göttingen, Germany. It was inspired by the existence of similar miniature breeds, like the Minnesota Minipig. The initial breeding stock was constituted of four male and two female Minnesota Minipigs and three male and eight female, and later four additional female potbelly Vietnamese pigs from German zoos. Göttingen Minipigs were selected for low body weight on the basis of the 154-day weight, then on weaning weight. In order to obtain a white color strain, German Landrace pigs were introduced by artificial insemination during the years 1965 to 1969. Separate lines of coloured and white minipigs were established. In 1969, the proportional representation of the three original breeds in the white line was 60% Vietnamese potbelly pigs, 33% Minnesota Minipigs and 7% German Landrace (reported by [START_REF] Simianer | Genetic management of the Göttingen Minipig population[END_REF]. From 1992, efforts were conducted in the Göttingen minipig population to improve population uniformity, health status, recorded breed characteristics, and to exclusively produce white animals. The current selection objective is mainly based on body weight and litter size. Additional traits are taken into account such as the number of teats, the hair coat, black hair or spots, abnormalities and different body shape parameters like saddle back, position of legs and jaws. As temperament is considered of major relevance for minipigs used as laboratory animals, further studies were conducted showing that behaviour scoring is of moderate heritability and that temperament can be included in selection objectives (Köhn et al., 2009).

Pure Yucatan strains were also developed from original 25 Yucatan pigs (also mentioned as hairless Mexican pigs) imported from Yucatan Peninsula of Mexico. [START_REF] Panepinto | The Yucatan miniature pig as a laboratory animal[END_REF] constituted an original breeding stock at Colorado State University from these imported pigs. Another miniature pure breed regained some interest within the last decade: the Ossabaw breed [START_REF] Sturek | Ossabaw Island Miniature Swine: Cardiometabolic Syndrome Assessment[END_REF]. Pigs originated from the Ossabaw Island (United States). Research on the obesity-prone Ossabaw miniature swine clearly indicates that these animals develop metabolism syndrom and cardiovascular disease when fed a high-calorie atherogenic diet, and might prove to be superior to Yucatan pig [START_REF] Neeb | Metabolic Syndrome and Coronary Artery Disease in Ossabaw Compared with Yucatan Swine[END_REF].

Selection for other criterion than small size (or coat color) was applied in some minipig populations. [START_REF] Phillips | Genetic selection for diabetogenic traits in Yucatan miniature swine[END_REF] developed two divergent Yucatan lines selected for low or high glucose-tolerance. These lines were developed in order to provide a model for diabete studies. [START_REF] Tikhonov | The Siberian miniature pig, its development , genetics and use in biomedical research[END_REF] established a minipig line (Minisibs), from Vietnamese, Swedish Landrace and wild boar, with chromosomal polymorphism (2n=38, 37, 36). MeliM breed, Munich miniature (Troll) and one Sinclair line were dedicated to melanoma research.

THE MELIM PIG BREED AND MELANOMA PREDISPOSITION

An example of a pig breed developed for a specific biomedical research, cutaneous melanomas, is the Melanoblastoma-bearing Libechov Minipig (MeLiM) breed [START_REF] Horak | Hereditary melanoblastoma in miniature pigs and its successful therapy by devitalization technique[END_REF]. This breed was constituted at the Libechov Institute of Animal Physiology and Genetics in the Czech Republic because of the similarities between human and porcine skin [START_REF] Zhang | Comparison of integrins in human skin, pig skin, and perfused skin: an in vitro skin toxicology model[END_REF]. Initially, an inheritable form of cutaneous melanoma has been described in the Hormel miniature pig strain. A strain of Hormel miniature pigs, known as Sinclair pigs were founders of a herd that was bred to have a high incidence of spontaneous melanomas. A similar pig model of melanoma was established in Germany Munich Miniature Swine. The history of MeliM breed was described by [START_REF] Hruban | Inheritance of malignant melanoma in the MeLiM strain of miniature pigs[END_REF]. The experimental herd of laboratory pigs was founded by importation of 5 animals of the Hormel strain from the USA in 1967. These animals were crossbred with Landrace, Large White, Cornwall, Vietnamese pigs and Göttingen pigs. A few black piglets with melanoma which occurred in this genetically heterogeneous population and the MeLiM strain was then established using selective breeding for melanoma. In 2004, the contribution of each ancestral breed was: Hormel strain 34-53%, Vietnamese pig 10-18%, Landrace 3-7%, Cornwall 1-6%, Large White 0-6% and Göttingen breed 19-34%, depending on the family.

For further studies, MeliM breeders were imported from the Czech Republic to the CEA/INRA Laboratory of Radiobiology and Genome Study in Jouy-en-Josas (France). With such a strain, it was possible to further describe tumour development, but as melanoma incidence was demonstrated as heritable, it was possible to explore the genetic determinism of melanoma susceptibility and unravel genes involved in incidence of spontaneous cutaneous melanoma. Studies has been conducted both for gene localisation with linkage and association analyses and gene expression in melanoma. Consequently, for genetic analysis, affected boars were crossed with healthy sows of the Duroc strain to create an F1 and a backcross generation [START_REF] Bourneuf | Genetic and functional evaluation of MITF as a candidate gene for cutaneous melanoma predisposition in pigs[END_REF] These crossbred animals were first analysed with 3 microsatellite markers linked to CDKN2A gene, a candidate gene for human melanoma [START_REF] Piepkorn | Melanoma genetics: an update with focus on the CDKN2A(p16)/ARF tumor suppressors[END_REF]. Allelic association, linkage analysis and haplotype analysis of these data led to exclusion of the CDKN2A gene as a candidate for melanoma susceptibility (Le Chalony et al., 2003). A second analysis with 137 microsatellites revealed 4 chromosomal regions potentially involved in melanoma susceptibility on chromosomes (SSC) 1, 2, 6 (near MC1R), 7 and 8. The 3 regions on SSC 1, 6 and 7 have counterparts on human chromosomes harboring melanoma candidate loci [START_REF] Geffrotin | Identification of five chromosomal regions involved in predisposition to melanoma by genome-wide scan in the MeLiM swine model[END_REF]. A third QTL analysis was performed on additional backcross individual with a higher number of markers ( 153). The results revealed QTL on SSC 1, 2, 13, 15 and 17. A focus on MC1R gene showed that one allele, which determines black coat colour in pigs, predisposes significantly to melanoma. A refined QTL analysis of the most significant QTL confirmed a linkage between SSC13 and melanoma susceptibility [START_REF] Bourneuf | Genetic and functional evaluation of MITF as a candidate gene for cutaneous melanoma predisposition in pigs[END_REF]. In order to investigate a potential association of MITF gene, an identified candidate gene localised on SSC13, with melanoma development. A single-marker association analysis with 23 SNPs conducted to exclude MITF gene as a causal gene. New developments to localize QTL focused on structural variation of genomes such as copy number variants (CNV) or insertion-deletions. A higher density genotyping, provided by the Illumina 64k pig Beadchip, was performed on the backcross pedigree to detect potential CNV events which might lead to new QTL localisations and candidate genes for melanoma susceptibility [START_REF] Corominas | Detection of CNVs throughout the genome of a porcine melanoma model. XXIst International Pigment Cell Conference (IPCC)[END_REF]. Gene expression in melanoma also allowed to identified putative candidate genes.

PIG GENETIC TOOLS

The pig genome sequence has been released [START_REF] Groenen | Analyses of pig genomes provide insight into porcine demography and evolution[END_REF]. This is a major step for pig genetic studies in general, and will prove to be an important resource in use of pig model for biomedical research. As an example, [START_REF] Groenen | Analyses of pig genomes provide insight into porcine demography and evolution[END_REF] studied predicted porcine protein sequences and compared them with their human orthologues. They observed 112 positions where the porcine protein has the same amino acid that is implicated in a human disease and most of these changes in humans have been shown to increase risk in multifactorial traits, including obesity and diabetes. Also the pig genome sequence was mainly based on a Duroc animal, a Göttingen pig has also been sequenced (Vamathevan et al., 2013). From comparison between genome sequences, one of their conclusions was that the differences identified between minipig when compared to the Duroc pig, highlight the fundamental importance of having gene information from the exact breed used in biomedical research. The consequence is that further genetic knowledge should be seeked out in each particular minipig strain, either in a very general way by resequencing, or in a more focus manner on gene/region of interest. Additionally, they provided evidence, at least for drug development, that some target genes for cardiovascular drugs differed from human leading to potential functional divergence. Thus binding assays or dose-ranging investigations of drug target relationships may need to be carefully interpreted.

The constitution and use of the 60k single-nucleotide polymorphism (SNP) chip (Ramos et al., 2009) improved our possibility to localize region of interest (for an update of detected QTL on various traits, see pigqtldb http://www.animalgenome.org/cgi-bin/QTLdb/SS/index), but also to better understand pig evolution through domestication and selection. A recent example is given by Burgo-Paz et al. (2013) on American pig populations, including Ossabaw and Yucatan pigs. They observed that Ossabaw and Yucatan pigs were clearly differentiated, despite an assumed shared Iberian ancestry. Yucatan was the closest breed to Spanish Iberian, whereas Ossabaw clustered among other American village pigs. An analysis of selection signatures revealed there could have been a common response among American populations involving a modification of their pattern of development and, perhaps, also by how they respond to external stimuli.

Another field of interest is the study of gut microbiota (the metagenome). Recent evidence suggests that the gut microbiota is an important contributing factor to obesity and related metabolic disorders. A study was conducted to characterise the intestinal microbiota in Göttingen minipigs and the Ossabaw minipigs [START_REF] Pedersen | Characterisation of Gut Microbiota in Ossabaw and Göttingen Minipigs as Models of Obesity and Metabolic Syndrome[END_REF]. Obesity-related changes in the composition of the gut microbiota were found in lean versus obese Göttingen and Ossabaw minipigs, buth the changes differed from one strain to the other. In both pig models diet seems to be the defining factor that shapes the gut microbiota as observed by changes in different bacteria divisions between lean and obese minipigs. They showed that both Ossabaw and Göttingen minipigs may be useful animal models in studying gut microbiota and provide a better understanding of obesity and metabolic syndrome and the relationship to the gut microbiota. It provided new insight in exploring this relationship in other minipig strains.

CONSTITUTION OF A NEW MINIPG BREED

As it has been explained, there are already numerous pig strains already used for biomedical research, some of them locally used, some of them largely used throughout the world for different research fields. For research involving physiological or genetic studies, the adequacy of the model, i.e. the pig strain, has to be carefully evaluated considering that the characteristics of the animals might less relevant in one strain than in another one. The constitution of new minipig breed at INRA was driven by the aim to produce a strain adequate for human nutrition studies oriented towards diabetes and obesity. The first step is to characterized different breeds (Vietnamese, Pitman-Moore and Yucatan) in order to estimate the (dis)advantages of each breed. Yucatan appeared to be a good basis, especially regarding behavior. Crossbreeding of different breeds at the beginning to benefit from the advantages of each foundation lines, to increase the variability of the traits and to decrease inbreeding, at least at the constitution of the line. The second step is to define a crossbreeding and/or selection strategy. Up to now, except selection conducted by [START_REF] Phillips | Genetic selection for diabetogenic traits in Yucatan miniature swine[END_REF] in Yucatan, no specific selection objective has been established in order to improve the adequacy of a pig line to human nutrition studies. Such a selection objective would have to take into account relevant metabolic parameters. This step has yet to be evaluated (which traits, which genetic parameters). Then the selection has to conducted minimizing the increase of inbreeding (using optimal contribution or a at least a family-based selection). Additional crossbreeding could be taken into consideration for introgressing interesting QTL or genes, using high density genotyping first for introgression then for controlling the initial genetic background.


G Machado

G Chagnon

D Favier

Theory and identification of a constitutive model of induced anisotropy by the Mullins effect

Keywords: Mullins effect, stress-softening, strain-induced anisotropy, constitutive equation

Rubber-like materials present a stress softening phenomenon after a first loading known as the Mullins effect. Some recent experimental data on filled silicone rubber is presented in literature, using uniaxial and biaxial tests to precondition samples thus induce some primary stress softening. A generic modeling based on the polymer network decomposition into an isotropic hyperelastic one, and a stress-softening evolution one, is proposed taking into account the contribution of many spatial directions. A new stress softening criterion tensor is built by means of a tensor that measures the repartition of energy in space. A general form of the stress softening function associated to a spatial direction is written by the way of two variables: one, the maximal eigenvalue of the energy tensor; the other, the energy in the considered direction. Finally, a particular form of constitutive equation is proposed. The model is fitted and compared to experimental data. The capacities of such modeling are finally discussed.

Introduction

Rubber-like materials present a stress softening after a first loading cycle, known as the Mullins effect [START_REF] Muhr | Effect of stretching on the properties of rubber[END_REF]. Different definitions have been given to the Mullins effect, in this paper the Mullins effect is considered as the difference between the first and second loadings. Moreover, different studies have highlighted that this phenomenon induces anisotropy, since the stress softening is strongly dependent on the second load direction.

In a first approach, many isotropic models were proposed in the literature to describe stress softening. First, physical models taking into account the evolution of the chain network were proposed. [START_REF] Govindjee | A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins's effect[END_REF] proposed a model based on the macromolecular network evolution by decomposition into a hyperelastic network and an evolving network. [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] considered that the macromolecular network can be represented by the eight chains model [START_REF] Arruda | A three dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF]. The model containing chain lengths and chain densities evolving with the maximal principal stretch. In another way, double network theory was developed [START_REF] Green | A new approach for the theory of relaxing polymeric media[END_REF] considering that the rubber-like material can be decomposed into a hard and a soft phase; the hard phase is transformed into soft phase with the stress softening. Different equations were proposed [START_REF] Beatty | A theory of stress-softening in incompressible isotropic materials[END_REF]Zuñiga and Beatty, 2002). At the same time, the damage theory was often used to describe the stress softening [START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects[END_REF][START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden type large strain elastic materials[END_REF][START_REF] Chagnon | On the relevance of Continuum Damage Mechanics as applied to the Mullins effect in elastomers[END_REF]. In another way, [START_REF] Li | A constitutive model dealing with damage due to cavity growth and the Mullins effect in rubber-like materials under triaxial loading[END_REF] associated the Mullins effect to the growth of cavities in the material and a compressible model was proposed. In a last point of view, [START_REF] Ogden | A pseudo-elastic model for the Mullins effect in filled rubber[END_REF] and [START_REF] Dorfmann | A pseudo-elastic model for loading, partial unloading and reloading of particle-reinforced rubbers[END_REF] proposed models based on pseudo-elasticity. All these models fit experimental data more or less accurately in one loading direction, i.e., without changing loading direction between first and second loadings. For a more exhaustive review about these isotropic models, the reader can refer to [START_REF] Diani | A review on the Mullins effect[END_REF].

To improve the modeling and to fit anisotropic stress softening, new approaches were developed taking into account the difference of stress softening in each strain direction.

At first, [START_REF] Göktepe | A micro-macro approach to rubber-like materials. Part III: The microsphere model of anisotropic Mullins-type damage[END_REF] generalized the approach proposed by [START_REF] Govindjee | A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins's effect[END_REF] taking into account a spatial repartition of the chains. In the same way, Diani et al. (2006a) proposed a generalization of the [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] model by means of chains oriented into 42 or more directions in space. Using a phenomenological damage function, this model can describe different stress softening in different directions, with permanent deformation after unloading. [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF] proposed a similar approach by taking into account the existence of different chains with different lengths in each direction. They integrate the density of probability in each direction, by taking into account the network evolution at each step. [START_REF] Shariff | An anisotropic model of the Mullins effect[END_REF] proposed an anisotropic damage model that describes transverse anisotropy of Mullins effect, taking into account different damages in the three principal strain directions using a second-order damage tensor. In the same way, [START_REF] Itskov | A thermodynamically consistent phenomenological model of the anisotropic Mullins effect[END_REF] proposed three damage evolution functions for the three principal strain directions. These functions are formulated in terms of material parameters that partly depend on the maximal principal stretch. Recently, [START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF] proposed a phenomenological model, based on the theory of pseudo-elasticity, which includes scalar variables in the strain energy function to account for stress softening and changes in material symmetry.

Most of the anisotropic models mentioned above are proposed by analyzing successive tensile tests performed along different directions. In spite of that, Machado et al. (2012a) has recently performed other original tests based on preconditioning with uniaxial tension and biaxial tension tests. Based on Machado et al. (2012a) experimental results using silicone rubber, this paper proposes a new approach for modeling the induced anisotropy by the Mullins effect. In Section 2, the global framework of the Mullins effect modeling is presented. In Section 3, a new approach is proposed to write constitutive equations by introducing a tensor that describes the strain energy repartition in the space directions.

The conditions to be verified by the equations are detailed. In Section 4, a first constitutive equation is proposed. It is fitted and compared to experimental data. Finally, Section 5 contains some concluding remarks and outlines some future perspectives.

Macromolecular approach to model Mullins effect

Filled silicone behavior

In the last few years, different tests highlighting the stress softening anisotropy have been presented in the literature for different rubber-like materials, see for example [START_REF] Muhr | Experimental determination of model for liquid silicone rubber: Hyperelasticity and Mullins effect[END_REF][START_REF] Besdo | Modelling inelastic stress-strain phenomena and a scheme for efficient experimental characterization[END_REF][START_REF] Hanson | Stress softening experiments in silica-filled polydimethylsiloxane provide insight into a mechanism for the Mullins effect[END_REF][START_REF] Diani | Observation and modeling of the anisotropic viscohyperelastic behavior of a rubberlike material[END_REF][START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF]. In this paper, attention is focused on the largest and most diverse database concerning Mullins effect anisotropy of a rubber-like material to the best of our knowledge. These data concern the results for the RTV3428 filled silicone rubber [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF](Machado et al., , 2012a)). First classical experimental tests, i.e., cyclic experiments with an increasing deformation after each cycle, were realized during tensile, pure shear and equibiaxial tensile tests.

The data are reported in [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF]. Second, stress softening anisotropy is presented in Machado et al. (2012a) induced by two distinguished preconditioning methods. The first one (noted as TT in the following), consists in a first loading in tension and a second loading also in tension, in four different directions. The second preconditioning method (noted as BT), consists in a first biaxial extension loading with constant principal strain direction. It is followed by a second loading in tension along the two principal strain directions of the first biaxial loading. The originality of these data is that loading states are very different between the first and second loads.

These new experimental results question the existing anisotropic constitutive equations and the main reasons are detailed here. The first reason is, Diani et al. (2006a) and [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF] models present an important permanent deformation that is related to the stress softening. But here, the material exhibits an important stress softening without permanent deformation. The other reason, for a second tensile loading orthogonal to the first loading, the models of [START_REF] Shariff | An anisotropic model of the Mullins effect[END_REF] and [START_REF] Itskov | A thermodynamically consistent phenomenological model of the anisotropic Mullins effect[END_REF] present a stiffer behavior than the virgin material, which is not the case of the filled silicone rubber. Last, all these models are based on a set of material directions and Mullins effect is controlled in each direction only by the maximum stretch reached during the deformation history along the considered direction. Recently, [START_REF] Merckel | A simple framework for full-network hyperelasticity and anisotropic damage[END_REF] analyzed the damage spatial repartition and proposed a softening criterion [START_REF] Merckel | A Mullins softening criterion for general loading conditions[END_REF] that is still the maximum stretch in each direction. Therefore, as pointed out in Machado et al. (2012a), a maximal deformation criterion that depends only on the considered direction is not enough to describe the stress softening for an arbitrary second load direction. This means that, if the maximal principal direction remains the same during the first and second load cycles, the strain energy can be a measure to quantify the Mullins effect in this direction. In the other directions, a coupling effect exists between different directions and it influences the stress softening. Under these circumstances, a new way to handle Mullins effect should be proposed at the sight of Machado et al. (2012a) experimental data.

Two networks theory

The results presented using silicone rubber-like materials highlight that unfilled silicone rubbers do not present stress softening [START_REF] Rey | Influence of the temperature on the mechanical behaviour of filled and unfilled silicone rubbers[END_REF] whereas filled silicone rubbers [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF] present an important one. For this silicone rubber material, it can be argued that the Mullins effect is principally due to the presence of filler in the material, which is not the case for every rubber-like material. In the light of these findings, a model based on [START_REF] Govindjee | A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins's effect[END_REF] theory is proposed. The main feature retained from their theory is the additive split of the strain energy into two contributions, motivated by the micromechanical structure represented in the Fig. 1. The hypothesis is based on distinguishing different macromolecular chains in the network, those that are linked to filler and those that are only linked to other macromolecular chains. It is assumed that only chains that are linked to fillers are concerned with Mullins effect. In [START_REF] Govindjee | Transition from micro-mechanics to computationally efficient phenomenology: Carbon black-filled rubbers incorporating Mullins's effect[END_REF] authors modified the initial approach into a phenomenological isotropic frame, introducing a normalized stress function that governs the damage level. Later, [START_REF] Göktepe | A micro-macro approach to rubber-like materials. Part III: The microsphere model of anisotropic Mullins-type damage[END_REF] conceptually extend the isotropic theory of [START_REF] Govindjee | Transition from micro-mechanics to computationally efficient phenomenology: Carbon black-filled rubbers incorporating Mullins's effect[END_REF] to the anisotropic case where damage history is described by one scalar for each material direction. The same additive split from [START_REF] Govindjee | A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins's effect[END_REF] is considered here, that means that the strain energy density (per unit of undeformed volume) over a representative elementary volume (REV) is decomposed into two parts

W = W cc + W cf (1)
whereW cf and W cc denote the energy densities of chains linked to filler network and chains linked to other chains network, respectively. On one hand, as it is considered that chains linked to other chains do not undergo Mullins effect, W cc is therefore represented by an isotropic hyperelastic energy density. On the other hand, W cf represents the anisotropy of stress softening induced by Mullins effect contributions in different directions of the REV.

Choice of the governing parameters of the Mullins effect

Analysis of literature experimental data

The different conclusions identified by Machado et al. (2012a) are analyzed and the consequences for the modeling are here in detail. During TT tests (uniaxial tension preconditioning followed by second tensile tests) it is shown that whatever is the second loading orientation, the strain-stress curves come back on the first loading curve at the same point corresponding to the maximum tensile stretch encountered during the first tensile loading.

Nevertheless the amount of stress softening depends on the angle between the first and second loadings. This means that the return on the tensile virgin curve is controlled by the maximal tensile stretch. However, the amount of stress softening depends on the relative orientation of first and second loadings.

The BT tests (biaxial tension followed by second tensile test) use preconditioning circular bulge test. Displacement and strain fields were obtained using three-dimensional image correlation measurements. In the preconditioning step, the circular bulge test specimen underwent biaxial loadings with different biaxiality ratios along a meridian. At the top of the bulge specimen, an equibiaxial loading is generated whereas a planar tension state is generated near the grips. Between these two points different biaxial states are generated [START_REF] Machado | Membrane curvatures and stress-strain full fields of axisymmetric bulge tests from 3D-DIC measurements. Theory and validation on virtual and experimental results[END_REF]. For different points along the bulge meridian, pairs of specimens were cut along circumferential and meridional directions and they were tested in tension. For each pair, the two different second stress-strain tensile curves, in the same way of TT tests, come back at the same point on the virgin tensile loading curve but with a different amount of stress softening according to the direction (circumferential or meridional). The conclusions are thus the same as TT tests but with a biaxial loading as the preconditioning test.

These results encourage to consider that the strain energy in the maximal principal strain direction is the governing parameter for the come-back on the first loading curve whatever is the second loading. Moreover, the stress softening amount in the other directions are linked to this parameter but it is attenuated if the direction of the maximal principal strain is not the same between first and second loadings. A measure for these quantities should thus be introduced.

New measure definition to quantify the stress softening

The strain energy is often used to describe the Mullins effect (see references in [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF]; [START_REF] Diani | A review on the Mullins effect[END_REF]) but its use is limited to the isotropic approach, since energy is a scalar global measure of the deformation state. Therefore, to compare the strain energy in different directions can be the clue. It is thus proposed to introduce a measure of the strain energy given by the contribution of each material direction.

Figure 2 illustrates the kinematics of an infinitesimal cone element extracted from the initial spherical representative volume element centered in P of radius dR 0 . In the initial configuration the slant height, surface area and volume are dR 0 , dS 0 and dV 0 = 1 3 dR 0 dS 0 respectively; the unity vector a 0 defines the material direction in the undeformed REV.

Considering the point Q, lying within an infinitesimal neighborhood of P , defined by the vector dx 0 = dR 0 a 0 . Under the deformation, this vector is mapped into dx = F dx 0 ,

where F is the deformation gradient. Thus, one obtains the following relation

dx = F dx 0 = dR 0 Fa 0 (2)
In the deformed configuration points p and q are referenced by the position vectors x and x + dx respectively; and the vector normal to the deformed surface dS given by the Nanson's relation

n = det(F) dS 0 dS F -T a 0 . ( 3
)
The velocity field within the infinitesimal neighborhood of x, with respect the reference frame R is given by

v q (x + dx) /R = v p (x, t) /R + L (x, t) /R dx (4) with L (x, t) /R = W (x, t) /R + D (x, t)
where L, W and D are the velocity gradient, spin and rate of deformation tensors. In any motion, the velocity field is locally decomposed as a sum of a rigid velocity v p (x, t) /R + W (x, t) /R dx and a straining velocity D (x, t) dx.

Considering Figure 3, the power P /R = t • (v q ) /R is expended by a force t = σn dS acting at point q, where σ is the Cauchy stress tensor. The interest is the expended power associated only with deformations. Then, it is possible to write the strain energy increment dM during a time increment dt, excluding rigid velocity, by the scalar product

dM = σn dS • Ddx dt (5)
where Ddx is the straining velocity field associated exclusively to the rate of deformation tensor D. Replacing Eq. (3) into Eq. ( 5), one obtains

dM = 3dV 0 det(F) [ F -1 σDF ] : (a 0 ⊗ a 0 ) dt (6)
Finally, the strain energy contribution in the a 0 direction is written, per unity of undeformed volume dV 0 , as

M(a 0 ) = 3 [∫ t 0 det(F) F -1 σDF dt ] : (a 0 ⊗ a 0 ) = 3 Ma 0 • a 0 (7)
This permits to introduce a tensor M in Eq. ( 7), defined in the reference configura-tion. It is decomposed into a symmetric M sym and skew-symmetric M skew part. Note that the product of a symmetric tensor and a skew-symmetric tensor has zero trace, i.e., M skew : (a 0 ⊗ a 0 ) = 0. Thus, in a general formulation, M sym describes the contribution of each material direction in the total strain energy. As symmetric tensor M sym possesses three real eigenvalues (M I > M II > M III ), where the maximal principal value M I is connected with an eigenvector determining the direction of maximum strain energy. Different noticeable parameters can be defined. The maximum strain energy for each conical elementary volume in direction a 0 along the history is

M max (a 0 ) = max τ ≤t M(a 0 , τ ) (8)
At the current time t, the maximum instantaneous strain energy in any direction is defined as

I(t) = M I (t) (9)
And last, the maximum strain energy in any direction in the history is defined as

G = max τ ≤t M I (t).
(10)

Construction of the evolution equation

An evolution function F is introduced along each direction, it describes the evolution of the network in the considered direction. The global strain energy is then rewritten as

W = W cc + ∫ V REV 0 F(a 0 )W cf (a 0 ) dV 0 . ( 11
)
where V REV 0 is the undeformed REV volume. At the sight of the silicone filled rubber experimental data, the function F(a 0 ) can be written as a function of the characteristic energy measures previously introduced

F(a 0 ) = F(M(a 0 ), M max (a 0 ), I, G) (12)
The main difference with the models from the literature is that the function F does not only depend on what happens in the considered direction a 0 but also on the global strain energy in the material, i.e., I and G. Then, different forms can be proposed.

In this paper, the concept for evolution function [START_REF] Beatty | A theory of stress-softening in incompressible isotropic materials[END_REF]Zuñiga, 2005;Zuñiga and Rodríguez, 2010) is used. This concept describes the Mullins effect by comparing the current deformation state and the maximum one, that means that the evolution function depends on the difference between the maximum and current deformation state. Thus, the value of introduced function remains one during the first load and it decreases when the current state differs from the maximum state. Thus, during a first loading I = G the function F(a 0 ) should not evolve if the material is stretched in a given direction for the first time, then

F(M(a 0 ), M max (a 0 ), I = G, G) = 1 (13)
During a second loading curve the function evolves, as the difference between the current and the maximum strain increases, in the interval given by

F(M(a 0 ), M max (a 0 ), I, G) ∈ [0, 1] (14)
This approach leads to a dependence in (G-I) and (M max (a 0 )-M(a 0 )) of the constitutive equation. Moreover, the amount of stress softening is directly linked to the orientation of the loading, the ratio of what happens in each direction compared to the maximum deformation should be taken into account. In this way, a general form is proposed as follows

F = 1 -F 1 (G -I) F 2 (M max (a 0 ) -M(a 0 )) F 3 ( M max (a 0 ) G ) (15)
where F 1 , F 2 and F 3 are functions to be determined. This multiplicative decomposition, also used in [START_REF] Rebouah | Anisotropic mullins stress softening of a deformed silicone holey plate[END_REF], is principally phenomenological, since stress softening is treated as a multiplicative function of the strain energy, see Eq. 11. The conditions evoked in Eq. ( 13) lead to

F 1 (G -I) = 0 if G = I (16) F 2 (M max (a 0 ) -M(a 0 )) = 0 if M max (a 0 )= M(a 0 ). ( 17
)
Now, different constitutive equation forms can be proposed for each function F 1 , F 2 and F 3 in Eq. ( 15).

The anisotropic constitutive equation

Hyperelastic constitutive equation

The advantage of such formulation is that the first loading curve is independent of the evolution function on the contrary of damage mechanics [START_REF] Lemaitre | Mechanics of solid materials[END_REF].

Thus, the choice of the hyperelastic energy only depends on the first loading curves. In a first approach, it is proposed to use the classical [START_REF] Mooney | A theory of large elastic deformation[END_REF] constitutive equation to represent the isotropic energy density, then

W cc = C 1 (I 1 -3) + C 2 (I 2 -3) (18)
For the anisotropic part of the constitutive equation the material could be represented by an infinite number of directions, introducing a probability density. But in this study, it is preferred to use a distribution of n-direction a (i) 0 oriented in any direction throughout the three-dimensional space instead of an integral formulation (Eq. 11). [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF] proposed different orientation schemes, that define the set of vectors a (i) 0 with different weight ω (i) for each direction to obtain a material as close as possible to an isotropic material when all the chains have the same mechanical behavior. W cf is then written as

W cf = n ∑ i=1 ω (i) F (i) W (i) cf (a (i) 0 ) ( 19
)
where n is the number of considered directions and

W (i) cf (a (i) 0
) is the hyperelastic strain energy of the chain in the initial direction a (i) 0 . The classical centrally symmetric n = 2×21 scheme was chosen to represent the material directions. The vector and weight of each direction can be found in [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF]. All the other direction distribution schemes could also be used. A comparative study of recently proposed integration schemes in application to a full network model of rubber can be found in [START_REF] Ehret | Numerical integration on the sphere and its effect on the material symmetry of constitutive equations -A comparative study[END_REF].

The non-Gaussian theory is classically used to capture the anisotropy. Diani et al. (2006a) and [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF] use the Langevin chain representation for

W (i) cf
energy. The great advantage of this choice is that it brings physical understanding to the modeling and it presents two main consequences. The first is that the zero-stress state is only ensured by the compensation of all the directions contribution as ∂W(λ (i) )/∂λ (i) ̸ = 0 if λ (i) = 1. Hence, this formulation could hardly be used for an initially non-isotropic material. The second one is that it allows to capture an important permanent deformation of the material after a loading cycle. However, in filled silicone experiments, it was shown that the permanent deformation is quite negligible. To this purpose, the classical hyperelastic anisotropic approach using the strain invariant

I (i) 4 = a (i)T 0 C a (i) 0 is used,
where C = F T F is the right Cauchy-Green strain tensor. The function should verify the following conditions

W (i) cf (I (i) 4 ) = 0 if I (i) 4 = 1 (20) ∂W (i) cf (I (i) 4 ) ∂I (i) 4 = 0 if I (i) 4 = 1. ( 21
)
In a first approach, an ordinary constitutive equation is used, considering that the chains are only stretched by tensile stresses. Otherwise, it is considered that compressive stretches lead to buckling. Thus, one may write

W (i) cf = K (i) 2 ( I (i) 4 -1 ) 2 if I (i) 4 ≥ 1 else 0. (22)
This formulation can be adapted to non-initially isotropic materials by choosing different functions for

W (i)
cf . As the filled silicone rubber is initially isotropic, every

W (i)
cf is initially the same in all directions, i.e., ∀i, j K (i) = K (j) .

stress softening constitutive equation

In part 3.3, a multiplicative decomposition was postulated in Eq. ( 15). The use of simple power functions, for F 1 , F 2 and F 3 , is proposed to represent the stress softening, given by

F (i) = 1 -η √ G -I G √ M (i) max -M (i) G ( M (i) max G ) 2 (23)
where η is a material parameter. The functions F 1 and F 2 are normalized according to the maximum strain energy G to ensure a normalized evolution function for each second loading curve. It is important to note that, even if the objective is to describe Mullins effect anisotropy, the constitutive equation for stress softening only depends on one parameter η.

All the other parameters describe the hyperelastic first loading. The evolution functions have the same form in all directions, but this approach could be extended to non-isotropic stress softening function by defining different values for the parameter η in the different directions.

There remains to verify that the presented model is in agreement with the requirements of thermodynamics (see e.g. [START_REF] Coleman | Thermodynamics with internal state variables[END_REF]. If only isothermal processes is considered, the Clausius-Duhem inequality must be satisfied by the conditions

- ∂W ∂M (i) max Ṁ(i) max 0 (24) - ∂W ∂G Ġ 0 (25)
where Ṁ(i) max 0 and Ġ 0. By means of straightforward manipulations of Eqs. ( 11), ( 24)

and ( 25) one can easily establish the above relations in terms of the evolution function i) . It is also important to show that

F (
∂F (i) ∂M (i) max 0, ∀i (26)
∂F (i) ∂G 0, ∀i. (27)
Considering the form of Eq. ( 23), the explicit form for Eq. ( 26) is given by

-η √ G -I G 1 G   1 2 ( M (i) max -M G ) -1 2 ( M (i) max G ) 2 + 2 ( M (i) max G ) √ M (i) max -M G   0, ∀i.(28)
First, an elementary study of the Eq.( 28) shows that all fractions terms are positive.

Second, when the stress softening evolves, i.e., G increasing, the maximum instantaneous strain energy I is equal to the maximum G. Thus, the function remains equal zero and the condition of Eq. ( 26) is automatically satisfied. In this way, the choice of F respects the conditions of Eqs.( 26) and ( 27) and consequently the Clausius-Duhem inequality is satisfied.

Comparison of the modeling with experimental data

The model is fitted on all the experimental data presented in [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF]Machado et al. ( , 2012a)), i.e., tests where the principal stretch directions remain unchanged during first and second load or tests where the principal stretch directions are not necessarily the same during first and second loads. First, the parameters of the hyperelastic constitutive equations are fitted on the different first loading curves. Different parameters can be chosen according to the repartition of the strain energy between W cc and W cf .

Considering a second tensile loading immediately before the sample rupture, the most stress softening level is obtained and this state corresponds to the strain energy of chains that were not affected by the Mullins effect. Thus, to ensure a good balance between W cc and W cf the portion W cc , i.e., the Mooney model, is fitted on the beginning of the second tensile loading curve at the higher deformation achieved before rupture.

Next, the part of W cf is fitted to complete the stress amount of the first loading curves.

The fitted parameters are presented in Table 1. (i) 0.20 MPa η 1.0

The last parameter that describes the stress softening is fitted on the second loading curves for all the tests, the value η = 1.0 is obtained. The condition in Eq. ( 14) must be satisfied, and as explained the function F (i) cannot be negative. If its softening is too large, i.e., F (i) < 0, the value F (i) = 0 is imposed. That means that in the considered direction a great number of chain-filler links were broken. In the second load, for the same direction, the suspended chains are no longer acting enough to impose a force on the macromolecular network, i.e., they do not contribute to the network entropic energy any more and their energy is thus lost [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF]. This assumption is consistent with the two networks theory and justified for relative short chains. Note that for longer molecular chains bonded at different places to fillers this assumption can be relaxed, for example, to take into account permanent set.

The simulations of the cyclic uniaxial tensile, pure shear and equibiaxial tensile tests are presented in Fig. 4. Concerning the first load, it appears that the model describes adequately uniaxial and pure shear tests whereas equibiaxial tests are underestimated.

This phenomenon is due to the hyperelastic equation and not to stress softening equation.

As pointed out by [START_REF] Marckmann | Comparison of hyperelastic models for rubber-like materials[END_REF] and [START_REF] Boyce | Constitutive models of rubber elasticity: A review[END_REF], there are very few hyperelastic constitutive models able to simultaneously simulate the both multi-dimensional data with a unique set of material parameters. Concerning the cyclic behavior, the form of the stress softening for all tests is quite well described. For uniaxial tensile and pure shear curves, the model has a slight tendency to underestimate stress softening and it is even more pronounced for equibiaxial tensile test. To finish, the model is used to simulate BT tests, the first load being a biaxial test and the second load being a tensile test. The comparison of the second loading curves is presented in Fig. 6. It appears that stress softening is moderately overestimated by the model, however the come-back on the first loading curve is perfectly described. All these simulations emphasize that the use of this new elongation energy measure is a good point to describe the come back of the second loading curves on the virgin one.

The amount of stress softening is well described for cyclic loading experiments and for TT tests, where the principal stretch directions were not the same during first and second loads. Nevertheless, the stress softening during BT test is overestimated. As observed in Fig. 4(c), stress was underestimated for the equibiaxial state, but this is due to the underestimation of hyperelastic strain energy obtained at the first load.

It can be noticed that the model describes correctly all the experimental tests with a simple constitutive equation that only depends on one parameter. Evidently, the results can be improved by proposing more complex constitutive equations, that consequently would lead to a significant increase in the number of parameters. Nevertheless, the presented results allowed to demonstrate the efficiency of this new approach.

Conclusion

This paper presents an original approach to model the stress induced anisotropy by the Mullins effect, by the definition of a tensor to measure the repartition of the strain energy in space. The comparison of the strain energy in different directions with the maximal principal strain energy permits to create a new formulation for stress softening modeling.

In this approach, the constitutive equation is written in function of the variation of strain energy in each direction and the variation of strain energy in the maximal principal strain direction. This new approach captures the principal characteristics of the Mullins effect underlined in literature. This new way of describing Mullins effect anisotropy can be a good starting point to elaborate new constitutive equations.

In this paper, a simple constitutive equation to describe the stress softening evolution was proposed. It clearly appears that the results are quite encouraging for a model that can describe many different types experimental tests, with very different strain histories, and the models presents only one material parameter. Of course, the agreement with the experimental data can be improved by using more sophisticated constitutive equation forms.

Acknowledge

We would like to thank the French ANR for supporting this work through the project RAAMO ("Robot Anguille Autonome pour Milieux Opaques")

Figure 1 :

 1 Figure 1: Representation of the silicone organization with macromolecular chains and filler particles

Figure 2 :Figure 3 :

 23 Figure 2: Kinematics of an infinitesimal cone element from the spherical REV.

Figure 4 :

 4 Figure 4: Comparison of the model (solid lines) with experimental data from Machado et al. (2010) (dotted lines) for: (a) cyclic uniaxial tensile test, (b) cyclic pure shear test and (c) cyclic equibiaxial test.

Figure 5 :

 5 Figure 5: Comparison of the model (solid lines) with TT uniaxial prestretching experimental data (dotted lines). (a) simulation of the model for different orientations of the second load. Details of the experimental (dotted lines) and modeled (solid lines) first and second load curves with an angle between stretch direction of: (b) 90 • , (c) 45 • , (d) 30 • , (e) 15 • and (f) 0 • .

Figure 6 :

 6 Figure6: Comparison of the model (solid lines) with BT biaxial prestretching experimental data (dotted lines). Curve a: simulation of the model for the second tensile load after an equibiaxial test; curve g: simulation of the model for the second load after a biaxial test of biaxiality ratio µ = 0.7; curve h: simulation of the model for the second load after a biaxial test of biaxiality ratio µ = 0.5

Table 1 :

 1 Values of the constitutive equation parameters

	Parameter	Value
	C 1	0.05 MPa
	C 2	0.03 MPa
	∀i K


Kostyantyn Slyusarenko

Doru Constantin

Patrick Davidson

A two-dimensional nematic phase of magnetic nanorods a)

Keywords: numbers: 82.70.Dd, 61.30.-v, 61.05.cf, 78.20.Fm

We report a hybrid mesophase consisting of magnetic nanorods confined between the non-ionic surfactant bilayers of a lamellar phase. The magnetic field-induced ordering of the nanorods was measured experimentally and modeled by a two-dimensional Onsager theory including the third virial coefficient. The nanorods are strongly confined in layers, with no orientational coupling from one layer to the next. At high volume concentration they exhibit spontaneous in-plane orientational ordering and form a stack of independent twodimensional nematic systems. This isotropic-nematic transition is first-order.

I. INTRODUCTION

Ordered phases in reduced dimensions are fascinating systems, with physical properties that can be qualitatively different from those of their three-dimensional counterparts. One of the most (intuitively) straightforward examples would be the two-dimensional nematic phase, with quadrupolar orientational order but no positional order. This deceivingly simple system has been the topic of extensive theoretical 1,2 and numerical [START_REF] Cosentino Lagomarsino | Isotropic-nematic transition of long, thin, hard spherocylinders confined in a quasi-two-dimensional planar geometry[END_REF][START_REF] Martínez-Ratón | Orientational ordering in hard rectangles: the role of the three body correlations[END_REF][START_REF] Vink | Liquid crystals in two dimensions: First-order phase transitions and nonuniversal critical behavior[END_REF] work, with a particular emphasis on the first-or second-order character of the isotropic-nematic transition.

There are however surprisingly few experimental realisations of such phases. The first qualitative observation involved rigid phospholipid tubules at the air-water interface [START_REF] Fisch | Observation of quasi twodimensional nematic order in a system of rigid rods[END_REF] . Similar results were obtained for Langmuir-Blodgett films of nanorods [START_REF] Kim | Langmuir-Blodgett nanorod assembly[END_REF][START_REF] Li | Semiconductor nanorod liquid crystals and their assembly on a substrate[END_REF] , but the particles could only be observed after deposition on a substrate. Another strategy relies on inserting DNA molecules within stacks of neutral [START_REF] Pott | DNA intercalation in neutral multilamellar membranes[END_REF][START_REF] Pott | DNA intercalation in neutral multilamellar membranes: Experiments and theory[END_REF] or charged [START_REF] Bouxsein | Two-dimensional packing of short DNA with nonpairing overhangs in cationic liposome-DNA complexes: From Onsager nematics to columnar nematics with finite-length columns[END_REF] lipid bilayers. These disordered mixtures were studied by X-ray scattering and the nematic phase was identified indirectly, by analyzing the shape of the DNA interaction peak. It is not clear whether these are purely two-dimensional phases or whether there is some interlayer coupling. We emphasize that the nematic order parameter was not measured in any of the systems above.

In this paper, we present a system where rigid magnetic nanorods are inserted into a soft lamellar matrix of nonionic surfactant. The resulting phase has three key advantages:

• It is easily aligned by thermal treatment, so that in the small-angle x-ray scattering (SAXS) images we can clearly discriminate between the direction of the smectic director and that in the plane of the layers. We can thus explore the anisotropy of the phase.

a) Originally published in J. Chem. Phys. 140 (10), 104904 (2014); http://dx.doi.org/10.1063/1.4867790 b) Electronic mail: doru.constantin@u-psud.fr

• The nanorods have magnetic properties, hence their orientation can be controlled using an applied magnetic field.

• When the phase is aligned in homeotropic anchoring, we can perform optical birefringence measurements (in particular under magnetic field), which yield a completely independent estimate of the orientational order of the particles.

The formulation of the phase [START_REF] Béneut | Magnetic nanorods confined in a lamellar lyotropic phase[END_REF] and its structural study [START_REF] Constantin | Lyotropic lamellar phase doped with a nematic phase of magnetic nanorods[END_REF] have already been presented. Here, we are concerned with the orientational order of the particles under magnetic field, that we describe using an Onsager-type model (up to the third virial coefficient), accounting for the electrostatic effects and for an applied magnetic field.

We conclusively demonstrate a first-order phase transition between two-dimensional isotropic and nematic phases, settling a long-standing debate in the theoretical and numerical literature [START_REF] Cosentino Lagomarsino | Isotropic-nematic transition of long, thin, hard spherocylinders confined in a quasi-two-dimensional planar geometry[END_REF][START_REF] Martínez-Ratón | Orientational ordering in hard rectangles: the role of the three body correlations[END_REF][START_REF] Vink | Liquid crystals in two dimensions: First-order phase transitions and nonuniversal critical behavior[END_REF][START_REF] Van Enter | First-order transitions for n-vector models in two and more dimensions: Rigorous proof[END_REF][START_REF] Almarza | Phase behavior of the confined Lebwohl-Lasher model[END_REF] .

II. MODEL AND ANALYSIS

The nanorods can be seen as parallelepipeds, with dimensions L × W × H and volume V = LW H. They bear a permanent magnetic moment µ along their long axis and exhibit negative susceptibility anisotropy ∆χ = χ -χ ⊥ < 0. For simplicity, we will describe them in the following as cylinders of length L and diameter D = 4W H/π.

The particle orientation is quantified by the distribution f (Ω), such that a fraction f (Ω)dΩ of particles have their long axes e within the solid angle element dΩ around direction Ω.

Our main result concerns the orientational distribution of the particles under confinement, but we start the discussion with their behavior in aqueous solution, a configuration that we will use as reference.

In isotropic solution at a volume concentration φ and in the presence of a magnetic field (taken along the polar axis, B = B ẑ), f (Ω) can be described by an extended Onsager model [START_REF] Vroege | Phase transitions in lyotropic colloidal and polymer liquid crystals[END_REF][START_REF] Lemaire | Physical properties of aqueous suspensions of goethite (α-FeOOH) nanorods; Part I: In the isotropic phase[END_REF] :

f (Ω) = 1 A exp KB cos θ + JB 2 P 2 (cos θ) - 32 π φ φ * |sin γ(Ω, Ω ′ )| f (Ω ′ )dΩ ′ (1)
with K = µ kB T and J = ∆χV 3µ0kB T , where P 2 is the Legendre polynomial of the second order, γ is the angle between two rods oriented along directions Ω = (θ, ϕ) and Ω ′ = (θ ′ , ϕ ′ ), while

φ * = 4D 2 L (D + X)(L + X) 2 (2)
is the spinodal concentration (above which the isotropic state is absolutely unstable), X is an effective length describing the electrostatic interaction between rods 16 (X = (ln A ′ + γ + ln 2 -0.5)/κ, with A ′ the amplitude of the electrostatic interaction, γ ≈ 0.577 Euler's constant, and κ -1 the Debye screening length), and A is a normalization constant imposing f (Ω)dΩ = 1. Equation ( 1) is implicit in f (Ω) and must be solved numerically. The resulting order parameter S is used to fit the experimental data.

When the particles are confined in the lamellar phase (with interlayer spacing d), we consider that within the same layer they have the same effective interaction as in solution so we model them using the 2D analogue of Equation ( 1), but this time the polar axis is along the director ẑ n and the magnetic field is applied in the plane of the layers B⊥n :

f (ϕ) = 1 A exp KB cos ϕ + JB 2 P 2 (cos ϕ) -η 2π 0 k 2 (ϕ, ϕ ′ )f (ϕ ′ )dϕ ′ -η 2 2π 0 2π 0 k 3 (ϕ, ϕ ′ , ϕ ′′ )f (ϕ ′ )f (ϕ ′′ )dϕ ′ dϕ ′′ + η M k B T 2π 0 cos 2 (ϕ -ϕ ′ )f (ϕ ′ )dϕ ′ , (3)
where the in-plane angle ϕ is measured with respect to the orientation of the magnetic field and η,

η = 4φd(L + X)(D + X) πD 2 L ,
is the effective surface fraction. By using D, we implicitly consider that the particles can rotate freely along their long axis, as in the isotropic solution. If they were to adopt a particular configuration, for instance with the largest face (L × W ) in contact with the surfactant bilayer, the values of X and η would be slightly changed. However, we have no experimental evidence for or against this assertion and, for simplicity's sake, we conserve the "cylindrical" model used above.

In (3) the third and fourth terms on the right-hand side describe the second and third virial coefficients of the Onsager interaction. The coefficient k 2 (ϕ, ϕ ′ ) has an analytical form,

k 2 (ϕ, ϕ ′ ) = λ + 1 λ | sin(ϕ -ϕ ′ )| + 2[1 + | cos(ϕ -ϕ ′ )|],
where λ = L + X D + X is the effective aspect ratio, and

the coefficient k 3 (ϕ, ϕ ′ , ϕ ′′ ) is calculated numerically 4 .
Note that Onsager's argument for neglecting higher-order virial terms does not apply in 2D [START_REF] Martínez-Ratón | Orientational ordering in hard rectangles: the role of the three body correlations[END_REF][START_REF] Straley | Liquid crystals in two dimensions[END_REF] , especially for rods of moderate aspect ratio (such as those used in our experiments). We therefore kept the third term in the virial expansion.

As presented so far, the 2D model considers that particles in the same plane interact by the effective hard-core model measured in solution, while particles in different layers do not interact at all. This latter feature is indispensable for having a true 2D phase, so it should be thoroughly checked. That is why we added in (3) a possible orientational coupling 19 between the nanorods in neighboring layers, described by the last term (with M the coupling amplitude). Below, we discuss possible physical origins for such a coupling, estimate its strength and show that it is negligible.

A. X-ray scattering

The X-ray scattering form factor of an individual particle is:

F (q, β, L, D) = F 0 sinc qL cos β 2 J 1 (qD/2 sin β) qD/2 sin β (4)
where β is the angle between the long axis of the rod e and the scattering vector q: cos β = cos θ sin ψ + sin θ cos ψ cos ϕ (see also Figure 1).

For a monodisperse population, the scattering signal is given by 20 : (5) while in the general (polydisperse) case Equation ( 5) must be averaged over a radius distribution g(D) (the polydispersity in the length L does not contribute over the available q-range.) In ( 5) we approximate the orientational distribution of the rods by the analytical Maier-Saupe form:

I(q, ψ) = 2 2π 0 dϕ π/2 0 f MS (θ)F 2 (q, β) sin(θ)dθ
f MS (Ω) = exp(m cos 2 θ) Z ( 6
)
where Z is a normalization constant and fitted the data with only one parameter m. The order parameter S can be expressed analytically as a function of m.

B. Optical birefringence

The birefringence of the samples can generally be expressed as:

∆n = ∆n sat S(B)φ (7)
where the field-dependent order parameter S(B) is the uniaxial one S for the case of water suspensions and the in-plane component 2 3 P of the biaxial order when the particles are confined within the lamellar matrix. ∆n sat is the specific birefringence, depending on the geometry of the particles and on their dielectric permittivity tensor.

III. RESULTS AND DISCUSSION

From electron microscopy and X-ray diffraction measurements, the average particle dimensions are:

L = 315 nm × W = 38 nm × H = 18 nm , (8)
and the relative polydispersity is of the order of 0.3 [START_REF] Constantin | Lyotropic lamellar phase doped with a nematic phase of magnetic nanorods[END_REF] . They correspond to cylinders with a diameter D = 4W H/π = 29.5 nm.

A. In solution

Under a magnetic field B, aqueous solutions of nanoparticles become uniaxial, with a director n B and an order parameter S that can be positive or negative, depending on the field amplitude [START_REF] Lemaire | Physical properties of aqueous suspensions of goethite (α-FeOOH) nanorods; Part I: In the isotropic phase[END_REF] .

From the SAXS images we determined the order parameter S(B) as a function of the field (Figure 2) by fitting the data to Equation (5) and using for the orientation distribution the Maier-Saupe form (6). The polydispersity g(D) is estimated from the scattering spectrum of a dilute suspension in the absence of the field. The order parameter is then computed over the distribution f (θ) determined from the complete model (1) as a function of B and φ, yielding estimates for the material constants K and J.

The spinodal concentration φ * = 7.1% was estimated from the independent measurement of the phase diagram of the aqueous suspension in zero field: φ * ≈ 0.95φ 3D N , where φ 3D N is the volume concentration of the nematic phase at coexistence (in aqueous solution). From (2) we then obtain X = 72 nm.

We also measured the birefringence ∆n(B) of the solutions (Figure 3) and fitted it to Equation (7), which involves the additional constant ∆n sat .

Finally, we obtain the material constants:

K = 8.2 ± 0.5 T -1 (9a) J = -13.8 ± 1 T -2 (9b) ∆n sat = 0.8 (9c) φ * = 7.1% ⇒ X = 72 nm (9d)
that we use alongside the particle dimensions (8) to study the hybrid system (particles inserted in the lamellar matrix).

B. In the lamellar phase

The SAXS measurements can be performed in either homeotropic or planar anchoring, giving access to both S and P , but the optical birefringence can only be determined in the homeotropic configuration (sensitive to P ). 7) and ( 1) with parameter values (9) In the planar configuration, with B = 0, for φ = 0.5%, 1.5%, and 3.5% we obtained S = -0.44, -0.46, and -0.47, respectively, very close to the perfect confinement limit S = -0.5. We conclude that, even in the absence of the field, the particles are practically contained within the plane of the layers (and the optical axis is along n).

If we now apply a field B⊥n the system becomes biaxial, with a second order parameter P .

In the following we treat the system as purely twodimensional, with S = -0.5. The particle population is then completely described by the orientation distribution f (ϕ), determined from the model (3) as a function of B, M , and φ. The material constants are those measured in aqueous solution (9), the periodicity of the lamellar phase d = 45 nm, [START_REF] Constantin | Lyotropic lamellar phase doped with a nematic phase of magnetic nanorods[END_REF] and M is the sole fitting parameter.

Isotropic two-dimensional phase

We describe the birefringence ∆n(B) using Equation (7) (with S(B) replaced by (2/3)P (B)) coupled with (3). The simultaneous fit of the three curves yields M = (-3 ± 3) k B T see Figure 4, and the resulting order parameter is compared to that measured by SAXS in Figure 5. The M values for the individual fits at 1.5 and 3.5 % are shown in Figure 7. At 0.5 % the uncertainty on M is very large, so we excluded this point from the analysis.

Nematic two-dimensional phase

The most concentrated system, with φ = 8 vol % in the lamellar phase is a biaxial nematic, with a spontaneous second order parameter P . It can thus be seen as a stack of two-dimensional nematic layers. To measure P we applied a small magnetic field of 30 mT in the plane of the layers to orient all nematic domains along the field, and measured the dependence of the x-ray scattering intensity I on ψ. Fitting I(ψ) by the Leadbetter method 21 , which takes into account the interaction between nanorods (and is more appropriate at high concentration than the model ( 5)), we obtained P = 0.51. Both in aqueous solution and in the lamellar phase the isotropic-nematic transition is first-order, as the two phases can coexist over a certain concentration range. By preparing several samples within this domain and estimating the fraction occupied by each phase we determine the nematic concentration at coexistence φ 3D N (in the aqueous solution) and φ N (for the two-dimensional nematic). In our system, these two distinct parameters happen to have the same numerical value, 7.5±0.5 vol %. In the lamellar phase, the isotropic-nematic coexistence range is approximately φ = 4.5 -7.5 vol %, corresponding to η = 0.38 -0.62 [START_REF] Constantin | Lyotropic lamellar phase doped with a nematic phase of magnetic nanorods[END_REF] . A sample at coexistence, with φ = 6.6 vol %, i.e. η = 0.54, is shown in Figure 6. The interface is not sharply defined, mainly due to the presence of the lamellar defects which create arbitrarily shaped domains of the two phases. The two-dimensional model (3) does indeed predict a first-order isotropic-nematic transition, with a φ N that corresponds to the experimental one (φ exp N = 7.5 ± 0.5 vol %) for M = (0.5 ± 0.5) k B T . For our most concentrated sample, with φ = 8 vol %, the same model reproduces the order parameter P = 0.51 for M = (0.07 ± 0.01) k B T . These two values correspond to the high-concentration points in Figure 7.

I/Lα N/Lα

Coupling

The coupling coefficient M is not significantly different from zero: at low concentration (in the two-dimensional isotropic phase) M = (-3±3) k B T from the birefringence data and at high concentration (in the two-dimensional nematic phase) M = (0.5 ± 0.5) k B T from the transition concentration and M = (0.07 ± 0.01) k B T from P (φ = 0.08). We conclude that, within the experimental precision, the orientation of particles in different layers is uncoupled, see Figure 7. Is this result in agreement with the theoretical estimations? The main difference between the particles in lamellar phase and those in isotropic solution is that the former are constrained to lie in parallel planes, and thus can more easily become parallel. Magnetic and electrostatic interaction might then be stronger than in three dimensions. We estimate their amplitude in the twodimensional nematic phase, at φ = 8 vol %:

The magnetic interaction between the permanent dipole moments of the particles W dd favors an antiparallel orientation and its magnitude is:

W dd k B T ≈ µ 0 µ 2 4πr 3 k B T ≈ µ 0 K 2 k B T φ 4πLDH = 6 • 10 -6 ,
where r is an average interparticle distance. The electrostatic twist energy 22 per molecule W tw (which favors a perpendicular orientation) is given by 23 :

W tw k B T = ηh tw 2 k tw (ϕ, ϕ ′ )f (ϕ)f (ϕ ′ )dϕdϕ ′ , (10)
where h tw = (κX) -1 . The coefficient k tw (ϕ, ϕ ′ ) is:

k tw (ϕ, ϕ ′ ) = λ + 1 λ [-| sin(ϕ -ϕ ′ )| ln | sin(ϕ -ϕ ′ )| -α| sin(ϕ -ϕ ′ )|], with α = - | sin(ϕ -ϕ ′ )| ln | sin(ϕ -ϕ ′ )|dϕdϕ ′ | sin(ϕ -ϕ ′ )|dϕdϕ ′ ≈ 0.3065.
The main uncertainty is in evaluating the Debye length: κ -1 = 2-14 nm, where the lower value was measured in a different goethite batch [START_REF] Lemaire | Physical properties of aqueous suspensions of goethite (α-FeOOH) nanorods; Part I: In the isotropic phase[END_REF] and the higher one corresponds to an ionic strength of 1 mM (the minimum possible at pH = 3). We then have W tw ≈ (4 -25) • 10 -3 k B T . Both effects are therefore negligible compared to the thermal energy, in agreement with our experimental findings.

IV. CONCLUSION

We formulated stable anisotropic materials where goethite nanorods are confined between the bilayers of a soft lamellar phase. The nanorods form a nematic system with director along the director of the lamellar phase and a negative order parameter S≤ -0.45 (to be compared with S = -0.5 for perfect confinement).

At low volume concentration φ < 7.5 vol %, the system is uniaxial in the absence of a magnetic field and forms a two-dimensional isotropic phase. Under a small magnetic field (0 < B < 350 mT) parallel to the layers, the distribution of the nanorods acquires a biaxial character (the long axis is preferentially oriented along the field), with a low order parameter P (P ≤ 0.05 for φ = 3.5 vol %).

At higher concentration φ ≥ 7.5vol % the biaxial order is spontaneous, yielding a two-dimensional nematic phase in the plane perpendicular to n, with an order parameter P ≃ 0.51.

We modelled the magnetic field-induced ordering of the nanorods P (B) by a two-dimensional Onsager theory and obtained good agreement with the experiment. In the simulation we took into account the magnetic properties of the nanorods and the effect of the electrostatic interaction and expanded the free energy to the third virial coefficient. We also included a possible orientational coupling M between nanorods in neighboring layers.

We conclude that M is not significantly different from zero (and, at any rate, much smaller than k B T ) and that the phase is a stack of isolated layers. The particles can be modelled as rigid rectangles with aspect ratio 3.8; at low concentration they exhibit a twodimensional isotropic phase which (on increasing the concentration) undergoes a first-order phase transition to a two-dimensional nematic phase with an order parameter P ∼ 0.5. This result is in very good agreement with numerical simulations [START_REF] Martínez-Ratón | Orientational ordering in hard rectangles: the role of the three body correlations[END_REF]24 .

In our analysis, we chose a particular effective shape (rigid rectangles) for the particles, mainly due to its tractability (very important for including the third virial coefficient).

Other choices can lead to a different phase diagram, in particular to a second-order phase transition 25 , in contrast with the experimental findings. We also ignored the effect of polydispersity, which could have profound consequences 26 . Hopefully, our experimental results will motivate further theoretical and numerical research in the field of two-dimensional ordered phases. 19 In principle, one could also have a positional coupling between layers. We do not include this type of term, because here we concentrate on the orientational order and because an extensive X-ray scattering investigation showed no characteristic sign (peak with qz = 0 and qr = 0) of such coupling. 20 B. J. Lemaire

ORDER PARAMETER

The order parameter tensor Q of the phase describes the distribution of the particle direction e: Q αβ = 1 2 (3e α e β -δ αβ ) (where the average • is taken over the distribution f (Ω)) and can be written in the principal axis frame as (Palffy-Muhoray and Hoatson 27 , Eq. 7):

Q αβ =   -1 2 (S -P ) 0 0 0 -1 2 (S + P ) 0 0 0 S   (11)
In spherical coordinates, with the polar axis along ẑ, the values of P and S are related to the distribution function f (θ, ϕ) via: S = 3 2 f (θ, ϕ) sin θ cos 2 θ dθdϕ -1 2

P = 3 2 f (θ, ϕ) sin 3 θ cos(2ϕ) dθdϕ (12)
For a uniaxial system the distribution f depends only on the polar angle, f (Ω) = f (θ), P = 0 and S is defined via the simplified formula: S = 3 2 f (θ) sin θ cos 2 θ dθ -1 2 . In the "complete confinement" case (relevant for the lamellar system), f (Ω) = δ(θ -π/2)f (ϕ), S = -1 2 , and P reduces to: P = 3 2 f (ϕ) cos(2ϕ) dϕ.

FIG. 1 .

 1 FIG.1. Experimental geometry. The scattering object is a rod with orientation e given by the polar angle θ (measured with respect to the nematic director n ẑ) and the azimuthal angle ϕ. The incident beam is along the ŷ axis. The detector records in the (x, z) plane (normal to the incident beam) the scattered intensity I(q, ψ) given by Equation(5)

FIG. 2 .

 2 FIG.2. Order parameter S(B) (determined from the SAXS images) for the aqueous suspensions of nanorods at various concentrations φ (symbols). The curves are fits derived from Equation (1) with parameter values(9)

FIG. 3 .

 3 FIG.3. Birefringence ∆n(B) of the aqueous suspensions of nanorods at various concentrations φ (symbols). The curves are fits derived from Equations (7) and (1) with parameter values(9)

FIG. 4 .FIG. 5 .

 45 FIG.4. Birefringence ∆n(B) of the suspensions of nanorods in lamellar phase at various concentrations φ (symbols). The curves are fits derived from Equations (7) (with S(B) replaced by (2/3)P (B)) and (3) with the material parameters in (9), d = 45 nm, and M = -3 kBT .

FIG. 6 .

 6 FIG. 6. Sample at coexistence (φ = 6.6 vol %, η = 0.54). The dashed line (with a greyed uncertain area) delimits the 2D isotropic phase I/Lα (upper left) from the 2D nematic phase N/Lα (lower right). The width of the capillary is 1 mm.

FIG. 7 .

 7 FIG. 7.Best estimate of M for various concentrations φ. For φ = 1.5 and 3.5 % we used the birefringence data (Figure4). The two high-concentration points are obtained from the transition concentration and from P (φ = 0.08) (see text).

  , P. Panine, J. C. P. Gabriel, and P. Davidson, "The measurement by SAXS of the nematic order parameter of laponite gels," Europhys. Lett., 59, 55 (2002). 21 A. J. Leadbetter and E. K. Norris, "Distribution functions in three liquid crystals from X-ray diffraction measurements," Molecular Physics, 38, 669 (1979). 22 This is simply the repulsion between identical and uniformly charged rods placed on top of each other in parallel planes. 23 A. Stroobants, H. N. W. Lekkerkerker, and T. Odijk, "Effect of electrostatic interaction on the liquid crystal phase transition in solutions of rodlike polyelectrolytes," Macromolecules, 19, 2232 (1986). 24 J. A. Cuesta and D. Frenkel, "Monte Carlo simulation of twodimensional hard ellipses," Phys. Rev. A, 42, 2126 (1990). 25 Y. Martínez-Ratón, E. Velasco, and L. Mederos, "Effect of particle geometry on phase transitions in two-dimensional liquid crystals," The Journal of Chemical Physics, 122, 064903 (2005). 26 Y. Martínez-Ratón, "Phase behaviour of mixtures of hard ellipses: a scaled particle density functional study," Liquid Crystals, 38, 697 (2011). 27 P. Palffy-Muhoray and G. L. Hoatson, "Mean-field theory of binary mixtures of nematic liquid crystals consisting of biaxial molecules," Phys. Rev. A, 44, 5052 (1991), in our notation symbols S and Q are swapped with respect to this reference.

ACKNOWLEDGMENTS

We acknowledge support from the Triangle de la Physique (project 2011-083T). The ESRF is acknowledged for the provision of beamtime (experiment SC-2393, ID02 beamline). The authors thank C. Chanéac for the goethite suspension, P. Boesecke, M. Impéror, A. Poulos, and B. Pansu for assistance with the SAXS experiments, J. Andrieu for performing some of the optical birefringence measurements and L. Navailles for stimulating discussions.