Guillaume Péron
email: guillaume.peron@univ-lyon1.fr

Benjamin Vollot

Differences in phenotypic variance between old and young congeneric species on a small island Running title: Phenotypic variance in island endemics

Keywords: individual heterogeneity, biometry, island syndrome, taxon cycle, functional trait

Aim: Insular environments theoretically promote within-species phenotypic variance, in part but not only because insular species enjoy a lower number of trophic interactions than their mainland counterparts. However, environmental factors apparently confound empirical comparisons between insular and mainland systems. To address this issue, we studied sympatric, congeneric species from the same small island, but that stemmed from separate colonization events. We predicted that the time since the species had colonized the island would correlate positively with the within-species variance in morphometrics. Location: Lifu island, New Caledonia, South Pacific Taxa: Zosterops inornatus, Z. minutus, Z. lateralis (Passeriformes) Methods: We measured 7 morphological traits and used principal component analysis to compute axes of variation. Results: The within-species variance in bill size and in feet/tail ratio corrected for body size increased with the estimated time since the species had been on the island. A similar but nonsignificant trend existed for body size and bill width. Main conclusions: The morphometrics of the three species of Zosterops birds from Lifu supported the hypothesis that insular environments promote within-species phenotypic variance, and cast a new light on the role of competition in the erosion of phenotypic variance.

Introduction

Instead of consistently converging towards a locally optimal phenotype [START_REF] Aubret | Genetic Assimilation and the Postcolonization Erosion of Phenotypic Plasticity in Island Tiger Snakes[END_REF][START_REF] Waddington | Genetic Assimilation of an Acquired Character[END_REF], most populations exhibit significant within-species phenotypic variation, which has been the topic of a rich literature [START_REF] Barabás | The evolution of trait variance creates a tension between species diversity and functional diversity[END_REF][START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF][START_REF] Clegg | 4000 Years of phenotypic change in an island bird: Heterogeneity of selection over three microevolutionary timescales[END_REF][START_REF] Ibarra-Isassi | Community-wide trait adaptation, but not plasticity, explains ant community structure in extreme environments[END_REF][START_REF] Violle | The return of the variance: Intraspecific variability in community ecology[END_REF]. In this paper, we first conducted a review of the reasons why island/mainland comparisons can provide insights into the maintenance of within-species phenotypic variation, but also highlighted the issue of environmental confounding factors in empirical island/mainland comparisons. We then described an approach to circumvent that bottleneck by comparing the species inside single islands.

Six reasons why insular populations should be more phenotypically diverse than mainland ones

Decreased trophic degree

In a species-rich system, the more variable a species is, the more it overlaps in niche with its competitors [START_REF] Fridley | Community and ecosystem effects of intraspecific genetic diversity in grassland microcosms of varying species diversity[END_REF] and the more diverse are its potential predators [START_REF] Schoener | Models of Optimal Size for Solitary Predators[END_REF] and pathogens [START_REF] Ganz | Benefits of host genetic diversity for resistance to infection depend on parasite diversity[END_REF]. By contrast, in species-poor island systems, the trophic degree does not increase much with the phenotypic variance [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. In addition, the more variable the phenotypes are, the more unique individuals are relative to conspecifics, and the less they compete with conspecifics [START_REF] Polis | Age structure component of niche width and intraspecific resource partitioning-Can age-groups function as ecological species?[END_REF][START_REF] Shine | Ecological Causes for the Evolution of Sexual Dimorphism: A Review of the Evidence[END_REF]. In other words, within-species variance is also a way to escape intra-specific density-dependence [START_REF] Tinker | Food limitation leads to behavioral diversification and dietary specialization in sea otters[END_REF], which is more intense in island than mainland systems [START_REF] Macarthur | Density Compensation in Island Faunas[END_REF][START_REF] Warren | Islands as model systems in ecology and evolution: Prospects fifty years after MacArthur-Wilson[END_REF]. These two mechanisms are traditionally invoked to explain the greater variability in island than mainland systems [START_REF] Macarthur | Density Compensation in Island Faunas[END_REF][START_REF] Soule | The "Niche-Variation" Hypothesis: A Test and Alternatives[END_REF]. More recently, this theory was updated with contributions from the niche variation hypothesis that species with broad niches are made of variable individuals with narrow niches [START_REF] Bolnick | Comparative support for the niche variation hypothesis that more generalized populations also are more heterogeneous[END_REF], and with the corollary that increased phenotypic variation feeds back into increased resistance to invasion by new competitors [START_REF] Barabás | The evolution of trait variance creates a tension between species diversity and functional diversity[END_REF]).

Jensen's inequality

Considering a mathematical function f of a trait x, if the function is concave, then its value at the population average 𝑥̅ is larger than the population average of the function value: 𝑓(𝑥̅ ) > 𝑓(𝑥) ̅̅̅̅̅̅ [START_REF] Jensen | Sur les fonctions convexes et les inégalités entre les valeurs moyennes[END_REF], and more importantly, the difference between 𝑓(𝑥̅ ) and 𝑓(𝑥) ̅̅̅̅̅̅ increases with the variance in x. Applying this rule to the strength of species interactions in trophic networks, one reaches the theory that the odds of any two species coexisting increases with their within-species phenotypic variance [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. Since island systems are extinction-filtered [START_REF] Warren | Islands as model systems in ecology and evolution: Prospects fifty years after MacArthur-Wilson[END_REF], the species that survived on islands would then be more variable than those that went extinct.

Portfolio effect

Negative correlations between the demographic performance of different phenotypes in a population, such as those caused by frequency-dependence and temporal fluctuations in the optimal phenotype [START_REF] Clegg | 4000 Years of phenotypic change in an island bird: Heterogeneity of selection over three microevolutionary timescales[END_REF][START_REF] Schindler | Population diversity and the portfolio effect in an exploited species[END_REF] can help different genotypes coexist in a population, but also shelter variable species from extinction [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. This is especially true in spatially-constrained systems, in which individuals cannot escape adverse conditions by moving away. Therefore, under the portfolio effect, the species that survived on islands should be more variable than those that went extinct.

Additive genetic variance

A high heritability makes it possible for species to evolve rapidly when environmental conditions change, including when new competitors enter the stage, thereby stabilizing community equilibria [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. By contrast, if heritability was low, ecological dynamics would dominate and competitive exclusion would ensue irrespective of within-species variation [START_REF] Warren | Islands as model systems in ecology and evolution: Prospects fifty years after MacArthur-Wilson[END_REF]. In island systems, the species that survived the colonization of their home island by new competitors should then be the ones with the most additive genetic variance.

Demographic stochasticity

In small populations, some phenotypes may contribute disproportionately to the next generation by chance only, thereby helping maintain the phenotypic variance (trait sampling; [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. This demographic stochasticity plays a larger role in small than large populations [START_REF] Gauthier | Partitioning prediction uncertainty in climate-dependent population models[END_REF][START_REF] Lande | Risks of population extinction from demographic and environmental stochasticity and random catastrophes[END_REF]. Therefore, it would contribute more phenotypic variation in island than mainland populations.

Island syndrome

In island systems, the low inter-specific density-dependence and high intra-specific densitydependence favor the evolution of "island syndromes", which involve long lives and low fecundity rates [START_REF] Covas | Evolution of reproductive life histories in island birds worldwide[END_REF][START_REF] Terborgh | The 'island syndrome' is an alternative state[END_REF]. As part of this life history strategy, traits that are associated with longevity should be buffered against environmental and genetic variation [START_REF] Péron | Evidence of reduced individual heterogeneity in adult survival of long-lived species[END_REF] and therefore exhibit low phenotypic variance. By contrast, traits that are mostly associated with fecundity should remain individually variable and may actually increase in variance compared to shorter-lived species, as a result of life history tradeoffs that preserve longevity. This would contribute to phenotypic variance in islands, at least for fecundity-associated traits.

Island/mainland comparisons can be confounded

The average environment may differ between island and mainland systems, as well as across islands.

Island systems may also exhibit less spatial and temporal heterogeneity in environmental conditions than mainland ones. These differences could inflate the variance in mainland populations, interfere with the above mechanisms, and explain why empirical island/mainland comparisons have so far yielded many negative results (Table 1). To control the influence of as many confounding factors as possible, we compared sympatric, congeneric species from the same small island. We especially looked for a set of species that exhibited a lot of spatial overlap even inside their home island. At first glance, the fact that all the species would be insular might seem counter-productive. Yet, in a given island system, there can be "young" and "old" species (sensu [START_REF] Fjeldså | Geographical patterns of old and young species in African forest biota: The significance of specific montane areas as evolutionary centres[END_REF], i.e., species that colonized the island at different points in time and have experienced the island environment for different amounts of time [START_REF] Ricklefs | Stage of Taxon Cycle, Habitat Distribution, and Population Density in the Avifauna of the West Indies[END_REF]. In some cases, this makes it possible to replace island/mainland comparisons with old/new comparisons, as described in the present study.

Material and methods

The Lifu Zosterops and their morphometrics

The Zosterops genus or white-eyes, known as the "great speciator" [START_REF] Diamond | Continental and Insular Speciation in Pacific Land Birds[END_REF][START_REF] Moyle | Explosive Pleistocene diversification and hemispheric expansion of a "great speciator[END_REF], has colonized the Old World in multiple waves, yielding complex patterns of secondary contacts that are just being deciphered (e.g., [START_REF] Cowles | Rapid, complete reproductive isolation in two closely related Zosterops White-eye bird species despite broadly overlapping ranges[END_REF][START_REF] Wickramasinghe | Non-sister Sri Lankan white-eyes (genus Zosterops) are a result of independent colonizations[END_REF].

Lifu island (1200km²) is part of the Loyalty archipelago in New Caledonia (South Pacific) that was historically free of terrestrial mammals. It lies about 100km away from the main island of New Caledonia. It hosts three endemic species or subspecies of white-eyes that are found nowhere else:

Z. inornatus, Z. minutus, and Z. lateralis melanops. These three species are not each other closest relatives [START_REF] Cornetti | The genome of the "great speciator" provides insights into bird diversification[END_REF] (Fig. 1). They did not diverge in situ and instead stem from three independent colonization events. Cornetti et al. estimated the colonization dates at 1.5Ma, 1.15Ma, and 0.5Ma respectively. This rare case of three-way insular sympatry after secondary contact raises the question of how the species coexisted and what were the consequences on their phenotypes. (2015). Lifu species are indicated in bold font with a thick tree branch. The kilometric distances correspond to the oversea distances between Lifu and the ranges of the species. Based on the genetic divergence with Zosterops albogularis, the first species to colonize Lifu was Z. inornatus. Z. minutus colonized next based on the divergence with Z. xanthochroa. Z. lateralis melanops was clearly the last species to colonize given its shallow divergence with the subspecies lateralis from mainland Australia and paraphyly with other subspecies in the region (not shown here) [START_REF] Cornetti | The genome of the "great speciator" provides insights into bird diversification[END_REF]. The cross symbols indicate recently extinct species. The inset represents the capture locations on Lifu (black for forest locations, white for edge locations).

We focused on the morphometrics in this study. In passerines, measurements of the limbs, tail, and bill are often good indicators of the ecological niche, including adaptation to various habitats and foraging strategies [START_REF] Marchetti | Correlates of wing morphology with foraging behaviour and migration distance in the genus Phylloscopus[END_REF][START_REF] Miles | Concordance of Ecomorphological Relationships in Three Assemblages of Passerine Birds[END_REF][START_REF] Pigot | Functional traits reveal the expansion and packing of ecological niche space underlying an elevational diversity gradient in passerine birds[END_REF]. Zosterops are mixed nectarivore-frugivore-insectivores, and when feeding on insects they exhibit a broad range of maneuvers including gleaning from branches and foliage, deaf-leaf probing, and aerial sallying [START_REF] Winkler | White-eyes, Yuhinas, and Allies (Zosteropidae)[END_REF]. In such generalist species, each individual could theoretically specialize on a few foraging maneuvers, or on a subset of the species' diet, depending on its individual morphology [START_REF] Bolnick | Comparative support for the niche variation hypothesis that more generalized populations also are more heterogeneous[END_REF][START_REF] Soule | The "Niche-Variation" Hypothesis: A Test and Alternatives[END_REF].

Due to the earlier colonization date of Z. inornatus, the aforementioned mechanisms that promote within-species phenotypic variation operated for a longer period on Z. inornatus than the other two species. In addition, under the taxon cycling hypothesis [START_REF] Ricklefs | Stage of Taxon Cycle, Habitat Distribution, and Population Density in the Avifauna of the West Indies[END_REF], Z. inornatus corresponds to a late-stage species (single island endemic), but it evolved to exploit most naturallyoccurring landcovers on its home island, and it apparently escaped interspecific competition by specializing on the lower strata of the vegetation and being the largest member of the genus [START_REF] Winkler | White-eyes, Yuhinas, and Allies (Zosteropidae)[END_REF]. For these reasons, we expected Z. inornatus to be the most phenotypically variable of our three study species.

By contrast, the newest species, Z. lateralis, would correspond to an early-stage species, i.e., a species that specializes on early-successional, non-forest landcovers of many islands rather than occupying all the landcovers of a single island (cf. results). Due to the need to fit in an already packed community module, Z. lateralis likely experienced more constraints on phenotypic variance on Lifu than anywhere else in its large geographic range. The mechanisms that promote within-species phenotypic variation thereby operated not only for the least amount of time on Z. lateralis, but also under the most constraints associated with inter-specific competition. For these reasons, we expected Z. lateralis to be the least variable of the three species.

In a nutshell, based on the combination of the different colonization times and the different ecological niches, we predicted that Z. inornatus should be more variable than Z. minutus which should be more variable than Z. lateralis.

Field work

All the study individuals were free-flying birds that we captured over a period of 14 days in November 2022, using mist-nets and playbacks of their vocalizations. We worked at 3 human-made forest edge locations (around the coastal town of Wé) and 2 relatively untouched forest locations (resp. near the settlement of Wedjumel and north of Wé) (Fig. 1). Within each location, we operated 3 to 8 mist-netting stations. The opening hours of the stations were not standardized in an attempt to maximize the number of first captures as opposed to the number of recaptures.

The same experienced person (B.V.) measured individual morphology using the same devices. We measured: tarsus length, tarsus width, wing length, tail length, bill length (the distance between the distal extremity of the nostril and the bill tip), bill height (at the base of the bill, along the vertical plane), and bill width (at the base of the bill, along the horizontal plane). We also determined the age (juvenile or adult) using skull ossification in 95 individuals. Owing to the limited sexual dimorphism, we did not document the sex of the study individuals, but we note that our prediction applies to all the within-species variance, including that caused by sexual dimorphism if any, and that both sexes incubate and display intraspecific aggression in this genus [START_REF] Winkler | White-eyes, Yuhinas, and Allies (Zosteropidae)[END_REF], limiting the risk that one of the sexes was over-represented in the sample.

Statistical analysis

We performed a principal component analysis (PCA) separately for the three bill measurements (bill length, bill height, and bill width) and for the four body measurements (wing length, tail length, tarsus length, and tarsus width). We obtained the three bill measurements for 149 individuals (39 Z.

inornatus, 65 Z. minutus, and 45 Z. lateralis) and the four body measurements for 144 individuals (38,63,42). The PCA for bill measurements yielded two independent axes of variation. The PCA for body size yielded three axes, of which we selected the first two that captured most of the variance (cf.

results). We used the R-package FactoMineR [START_REF] Husson | FactoMineR: Multivariate Exploratory Data Analysis and Data Mining with R[END_REF].

We did not log-transform the measurements before analysis. When a dataset contains outliers, logtransformation offset their disproportionate contribution to empirical moments and improves the fit of the data to the normal distribution. Since there was no outlier in our dataset however, logtransformation would hardly have improved the fit to the normal distribution, but would have reduced the variance in some of the traits, eventually influencing the conclusions (not shown) (see also [START_REF] Feng | Log-transformation and its implications for data analysis[END_REF].

We applied the Fligner-Killeen nonparametric test of homogeneity of variances [START_REF] Conover | A Comparative Study of Tests for Homogeneity of Variances, with Applications to the Outer Continental Shelf Bidding Data[END_REF], on each of the 4 axes of variation (two for bill and two for body measurements) to determine whether the phenotypic variance differed between species. We also tested each of the 12 pairwise differences individually, using variance F-tests with the usual statistical threshold (P = 0.05) and the Bonferroni-adjusted threshold (P = 0.004, likely to be too conservative here).

Finally, we investigated whether the prediction was better matched by the absolute values of the traits, or by their values relative to the individuals' body sizes. To do that, we divided each trait by the tarsus length (TL). A high trait/TL ratio indicates an individual that exhibits a larger trait than expected from its body size. To also rescale relative to the species' average, we computed the coefficient of variation of the trait/TL ratio: 𝑅𝐶𝑉 𝑠 (𝑡𝑟𝑎𝑖𝑡) = 𝑉𝑎𝑟 𝑠 (𝑡𝑟𝑎𝑖𝑡 𝑇𝐿) ⁄ 𝑀𝑒𝑎𝑛 𝑠 (𝑡𝑟𝑎𝑖𝑡 𝑇𝐿) ⁄

. We tested whether the RCV decreased from Z. inornatus to Z. minutus to Z. lateralis using [START_REF] Feltz | An asymptotic test for the equality of coefficients of variation from k populations[END_REF] multipopulation asymptotic test of equality between coefficients of variation (R package cvequality:

Marwick & Krishnamoorthy, 2019).

Comparison with published data

We were aware of another set of morphometric data, that was collected in Oct-Nov 2011 and available in the form of data summaries (PCA results) (Cornetti et al., 2015, Fig. 2a). There were some technical differences: Cornetti et al. included the body mass in the PCA, measured the bill at different cross sections than we did, and did not measure the tarsus width. In addition, several different field workers contributed to that dataset. Observer bias and differences between measuring devices can contribute significant amounts of artefactual variance in the morphometrics of small passerines. For these reasons, we did not pool the two datasets together, but instead analyzed them separately using the same procedure.

Results

The proportion of Z. lateralis in the captures increased from 0-4% at forest locations to 36-46% at edge locations, while the ratio between Z. inornatus and Z. minutus decreased from 6:5 to 1:5. These data indicated that Z. lateralis specialized on edge locations. There was no evidence of a difference in morphology between edge and forest locations in either species (all t-tests P > 0.25). The proportion of juveniles was 61%, 67% and 41% juveniles in Z. inornatus, Z. minutus and Z. lateralis respectively.

These data suggested a possible difference in breeding phenology between Z. lateralis and the others. There was no evidence of a difference in morphology between adults and (fledged) juveniles in either species (all ANOVA P > 0.20).

Four axes of variation among morphological features

In the PCA of the bill measurements, the first axis of variation, which we interpreted as an index of individual bill size, accounted for 95% of the variance in three measurements. The second axis accounting for the remaining 5% captured variation in bill width not explained by bill size.

In the PCA of the body measurements, the first axis of variation, which we interpreted as an index of individual body size, accounted for 92% of the variance in four measurements. The second axis accounted for 6% of the variance and captured a negative correlation between the remaining variation in feet size and tail length.

Comparing the within-species variance across species

The data contained very strong evidence that the variance in bill size decreased from Z. inornatus to Z. minutus to Z. lateralis (Fig. 2 and Table 2: "Bill 1"; Fligner-Killeen test: ²(2) = 9.76, P = 0.008), and moderate evidence regarding the feet/tail ratio (Fig. 2 and Table 2: "Body 2"; ²(2) = 6.84, P = 0.03).

These results supported our prediction. Bill width (Fig. 2 and Table 2: "Bill 2") and overall body size (Fig. 2 and Table 2: "Body 1") on the contrary provided little support for the prediction (Fligner-Killeen tests: both P > 0.2).

Table 2: Pairwise F-tests of differences between variances in the four axes of variation and between the 3 pairs of species. Statistically significant differences (P-value < 0.05) are highlighted in bold font, and those that remain significant after Bonferroni adjustment (P-value < 0.004) are underlined. (a) Variance in the bill measurements. (b) Variance in the body measurements. In both cases, the upper triangle is for the first axis of the PCA and the lower triangle is for the second axis of the PCA. Thus the conclusions from the two datasets were similar: some aspects of the morphology supported the prediction and the aspects that did not support the prediction did not contradict it either.

After the rescaling procedure (Eq. 1), the trends hold (Fig. 3) but the statistical support disappeared

Discussion

After reviewing the mechanisms that lead to the prediction that insular populations should be more phenotypically diverse than mainland ones, and after highlighting that empirical island/mainland comparisons seem to be confounded by environmental factors, we set out to Lifu island. The Lifu Zosterops birds provided an opportunity to test the prediction while controlling the influence of many potentially confounding environmental factors. We predicted that Zosterops inornatus should be more morphologically variable than Z. minutus which should be more variable than Z. lateralis, at least for those traits that indicate the individual ecological niche or niche breadth. The morphometrics of ~140 individuals supported that prediction in two out of four axes of variation, with nonsignificant trends in the other two, prompting us to call for continued research into the morphological variance inside insular guilds rather than between insular and mainland systems.

Ecological function of the four axes of variation in morphology

Bill PCA1 corresponded to the bill size, which across insectivorous passerine species correlates positively with the frequency of probing maneuvers (for nectar, fruit pulp, in the leaf litter, etc.) and negatively with the frequency of gleaning from branches and foliage [START_REF] Marchetti | Correlates of wing morphology with foraging behaviour and migration distance in the genus Phylloscopus[END_REF][START_REF] Miles | Concordance of Ecomorphological Relationships in Three Assemblages of Passerine Birds[END_REF][START_REF] Pigot | Functional traits reveal the expansion and packing of ecological niche space underlying an elevational diversity gradient in passerine birds[END_REF]. Thus, we might speculate that the variance we observed in bill size may correspond to varying levels of individual specialization into probing vs. gleaning.

Body PCA2 corresponded to a negative correlation between feet size and tail length after correcting for body size. This covariation between traits might have to do with the birds' agility when foraging for fruit and insects [START_REF] Marchetti | Correlates of wing morphology with foraging behaviour and migration distance in the genus Phylloscopus[END_REF] or when escaping avian predators [START_REF] Adriaensen | Stabilizing selection on blue tit fledgling mass in the presence of sparrowhawks[END_REF]. Predators of forest-dwelling passerines are indeed relatively numerous on Lifu, (Accipiter fasciatus, Circus approximans). If the feet/tail ratio provided anti-predator benefits, we would have expected the most typically insular species and therefore the species with the most buffering against lifespan-reducing threats, Z. inornatus, to exhibit the least variance the feet/tail ratio. This is not what we observed. Therefore, we might speculate that the feet/tail ratio mostly helped with foraging efficiency.

Body PCA1 corresponded to body size, which is highly diverged in our three species. Indeed, Z. inornatus and Z. minutus are respectively the largest and smallest members of a 105+ species-strong genus [START_REF] Winkler | White-eyes, Yuhinas, and Allies (Zosteropidae)[END_REF], which suggests that these two species underwent an episode of character displacement. Since traits that are under strong directional selection are more buffered against environmental and genetic variation than traits under weak selection [START_REF] Waddington | Genetic Assimilation of an Acquired Character[END_REF], it is maybe not surprising that the variation in body size was more constrained than other aspects of the morphology. Our data suggested that individual variation in body size within species did not reflect individual variation in ecological niche.

Bill PCA2 corresponded to bill width corrected for bill size, which, among passerines and near passerines, often indicates the diet. Very thin bills can indicate a nectarivore diet whereby the bird probe flowers for nectar, whereas very broad bills usually indicate a sallying insectivore [START_REF] Pigot | Functional traits reveal the expansion and packing of ecological niche space underlying an elevational diversity gradient in passerine birds[END_REF]. The Lifu Zosterops do exhibit these two foraging maneuvers, but infrequently only (pers. obs.).

The infrequent use of foraging maneuvers that require a particularly thin or broad bill could have contributed to the relative lack of variance in bill width.

A caveat regarding the focus on morphology as opposed to directly quantifying the ecological niche

In this study we assumed but did not demonstrate that the individual morphology conditioned the individual behavior (e.g., the preferred foraging maneuver, preferred foraging strata, or the preferred diet) and in turn the individual ecological niche. However, individuals may conceivably buffer the effect of their traits on their fitness. For example individuals may buffer the effect of their morphology on their survival rate by adopting risk-minimizing behaviors [START_REF] Péron | Evidence of reduced individual heterogeneity in adult survival of long-lived species[END_REF].

Alternatively, individual behaviors and individual morphologies may be uncorrelated, in which case behavioral traits if ever documented would offer the opportunity for an independent test of the prediction we made in this study. In any case, we note that the use of morphological traits as functional traits is traditional and well substantiated in avian ecology (e.g., [START_REF] Luck | Improving the application of vertebrate traitbased frameworks to the study of ecosystem services[END_REF][START_REF] Pigot | Functional traits reveal the expansion and packing of ecological niche space underlying an elevational diversity gradient in passerine birds[END_REF].

Character displacement in Z. lateralis

Our data indicated that Z. lateralis melanops was an edge specialist on Lifu, whereas Z. inornatus and Z. minutus occurred in both edge and interior forest but were more frequent in the interior, especially the former species. Z. inornatus and Z. minutus were thus single-island endemics that used all the naturally-occurring landcovers on their home island, which corresponded to a late-stage species under the taxon cycling hypothesis [START_REF] Ricklefs | The concept of the taxon cycle in biogeography[END_REF][START_REF] Ricklefs | Stage of Taxon Cycle, Habitat Distribution, and Population Density in the Avifauna of the West Indies[END_REF]. By contrast, we would interpret Z. lateralis as an early-stage species, because it is widespread and expanding (e.g., reaching New Zealand in historical time), but on Lifu it was specialized on earlysuccessional landcovers. Another perspective on the same pattern is that Z. lateralis was relegated to early-successional landcovers by competition with the older species. In any case, the result that Z. lateralis was the least morphologically diverse helped explain its ability to coexist with two older congeneric species on a single small island, in line with broader theories of character displacement and coexistence [START_REF] Taper | Models of Character Displacement and the Theoretical Robustness of Taxon Cycles[END_REF].

To further support this view that interspecific competition played a role in the erosion of phenotypic variation in Z. lateralis, and that this role was additive to the effect of the short time that Z. lateralis spent on Lifu, we turned to the long term study of Z. lateralis chlorocephalus on Heron island (off eastern Australia) [START_REF] Clegg | 4000 Years of phenotypic change in an island bird: Heterogeneity of selection over three microevolutionary timescales[END_REF]. The Heron island population has been exposed to the island environment for a shorter period of time than the Lifu population (4000 years: [START_REF] Clegg | 4000 Years of phenotypic change in an island bird: Heterogeneity of selection over three microevolutionary timescales[END_REF]. But on the other hand, Z. lateralis is the only representative of the Zosterops genus on Heron island and does not have to compete with two congeneric species as is the case on Lifu. If only the time since the populations had colonized their respective islands played a role, then the Heron island population should have been unambiguously less variable than the Lifu one. Although differences in methodology prevented direct comparisons between our results and Clegg et al.'s, that prediction did not seem to be met, especially by the bill and tarsus lengths. One of the reasons for the rapid emergence of large within-species variation on Heron island was reported to be the displacement of the optimal phenotype from one year to the next, owing to environmental stochasticity [START_REF] Clegg | 4000 Years of phenotypic change in an island bird: Heterogeneity of selection over three microevolutionary timescales[END_REF] (termed the portfolio effect in our introduction, following [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF]. On Lifu, we can speculate that the competition with conspecifics stabilized the optimal phenotype of Z. lateralis.

Comparing old and young insular taxa instead of mainland and insular taxa

Island/mainland comparisons have so far yielded many negative results regarding the predictions that insular taxa should be more phenotypically diverse than mainland taxa (Table 1). Our study helped contextualize some of these previous results by highlighting not only the role of environmental confounding factors, but also the consequences of the colonization sequence. We demonstrated that different species from the same small island may exhibit different amounts of phenotypic variance. Therefore, island/mainland comparison studies should strive to include all the species from insular guilds, and to factor in the colonization sequences, which might be achieved by treating the different stages of the taxon cycles separately in meta-analyses. In extreme cases, the island systems may actually have operated as a source for the mainland systems [START_REF] Moysi | Genomic data reveal contrasting patterns of divergence among island and mainland birds of the Eastern Mediterranean[END_REF][START_REF] Vinciguerra | Island life accelerates geographical radiation in the white-eyes (Zosteropidae)[END_REF], rather than the opposite which is the default working hypothesis. Such a situation would lead to completely opposite predictions regarding phenotypic variance. Lastly, we speculated a role for the shape of the selection gradient. In our study, body size did not verify our prediction and seemed likely to be, or at least have been, under strong directional selection. By contrast, the feet/tail ratio and bill size verified our prediction and seemed likely to be under stabilizing selection

Conclusion

This study advanced our understanding of the role of the insular environment in the maintenance of phenotypic variance, and specifically the lack of interspecific competition. We did that both theoretically through our review of the mechanisms and empirically in a non-model system that we explored for the first time. Our results based on a moderate sample size were encouraging but also illustrated the challenges of working with empirical variances, and the fact that our working hypothesis rested on the assumption that the morphological traits indicated the individual ecological niche, making it critical to select appropriate traits. We recommend publishing the raw phenotypic variances in addition to the population means in all studies of insular populations, and also envision insights from longer-term or more in-depth studies into the behavioral, demographic, and genetic basis of phenotypic variance.

Fig. 1 .

 1 Fig. 1. Simplified cladogram redrawn from[START_REF] Cornetti | The genome of the "great speciator" provides insights into bird diversification[END_REF]. Lifu species are indicated in bold font with a thick tree branch. The kilometric distances correspond to the oversea distances between Lifu and the ranges of the species. Based on the genetic divergence with Zosterops albogularis, the first species to colonize Lifu was Z. inornatus. Z. minutus colonized next based on the divergence with Z. xanthochroa. Z. lateralis melanops was clearly the last species to colonize given its shallow divergence with the subspecies lateralis from mainland Australia and paraphyly with other subspecies in the region (not shown here)[START_REF] Cornetti | The genome of the "great speciator" provides insights into bird diversification[END_REF]. The cross symbols indicate recently extinct species. The inset represents the capture locations on Lifu (black for forest locations, white for edge locations).

Fig. 2 .

 2 Fig. 2. (a) Individual ordinations in the PCA of the bill measurements. (b) Same for body measurements. (c) Within-species variance in the resulting four axes of variation.

  al. data contained moderate support for the prediction in their second axis of variation (Fig.S1: "PCA2"; ²(2) = 3.77, P = 0.05) and none in their first axis of variation (P = 0.53).

(

  Feltz-Miller CV tests: P > 0.25) except for tail length (P = 0.02; Fig. 3b: thin line). In other words, our prediction was better matched by the absolute values of the traits than by their values relative to individual body sizes.

Fig. 3 :

 3 Fig. 3: Relative coefficient of variation (RCV) for 5 morphometric measurements. The ratio between tarsus width and tarsus length was almost invariant in all three species (RCV=0; not shown).

Table 1 : Overview of previous studies that compared the within-species phenotypic variance between island and mainland systems

 1

	Reference	Location	Taxa	Trait	Main result
	(Kavanagh & Burns, 2015) (Madeira, 2018) African continent/ New Zealand /outlying islands São Tomé (Meiri et al., 2005) Worldwide (Novosolov et al., 2018) Worldwide	dioecious plants Organ sizes passerines Body size and plumage carnivores Cranial features Mainland No predictable difference in the degree of sexual dimorphism No predictable difference in the degree of sexual dimorphism species are more variable on average lizards Niche breadth Niche breadths are similar on average

provided insightful comments that greatly improved a previous version. This research received no external funding.

the Loyalty Islands gave permit for the captures. Two reviewers

Data availability:

The source data will be deposited onto Dryad or similar upon acceptance Biosketch: The Lifu Zosterops project aims at leveraging the configuration of the avian community of Lifu island to test hypotheses about mutation accumulation, Allee effects, and species coexistence.

Author contributions: G.P. conceived the ideas; B.V. and G.P. collected the data; G.P. analyzed the data; and G.P. led the writing.

Supplementary material

Fig. S1 and Table S1: reanalysis of [START_REF] Cornetti | The genome of the "great speciator" provides insights into bird diversification[END_REF] morphometric datasets using the same procedure as the present study.


Josephine Faddoul
email: josephine.faddoul@insa-rennes.fr

Pierre Rahme

Dominique Guines

Lionel Leotoing

Keywords: Glass fiber reinforced thermoplastics, Strain rate, Temperature, Mechanical characterization, Heat assisted forming

Thermoplastic composite materials are being increasingly used in many domains, such as automotive, marine and aeronautics. This growing interest is due to the relatively good mechanical properties and the recycle-ability of these materials. In various heat assisted forming technologies, thermoplastic composites may be subjected to severe conditions such as high temperatures, complex strain paths and potentially significant strain rates. This is expected to highly influence the overall behavior of these materials and make the characterization of their mechanical behavior more complicated. In this paper, the mechanical behavior of discontinuous glass fiber reinforced polypropylene composite is investigated at temperatures going from room temperature (RT) to 120 • C and strain rates ranging from quasi-static conditions to 10 s -1 by means of uniaxial tensile tests. The effect of material orientation is also investigated at RT and quasi-static strain rate of 0.001 s -1 pointing out a clear anisotropic behavior. The tensile properties including Young's modulus E, ultimate tensile strength σ ult and strain at failure ε f are determined based on the Standard ISO-527 reporting an obvious sensitivity of the material to temperature and strain rate. Based on literature review, the phenomenological isotropic constitutive model of G'Sell and Jonas, originally designed for unreinforced semicrystalline polymers, is selected to describe the mechanical behavior of the studied material. A modification of the original constitutive equation is proposed for a better representation of the overall behavior of the material at different temperature and strain rate conditions. Finally, the validity of this phenomenological law calibrated by means of the uniaxial tests, is evaluated for an equi-biaxial loading state close to that encountered by the material during heat assisted sheet forming process. The significant disagreement highlighted between numerical and experimental results proved that the uniaxial tensile tests are not sufficient to characterize the behavior of the material when a complex loading condition is imposed.

Introduction

Composite materials are now widely used in many fields, such as aeronautics, marine and automotive. This is mainly due to the relatively high strength/weight ratio that these materials can offer [START_REF] Tanasa | Fiber-reinforced polymer composites as structural materials for aeronautics[END_REF] . There exist several types of composites in which the structure can be composed with different combinations of reinforcement/matrix. Constituents should be wisely selected in order to fulfill the requirements of the intended application as well as the product lifetime, complexity of shape, costs... [START_REF] Verma | Others Coir fibre reinforcement and application in polymer composites[END_REF] More particularly, thermoplastic composites exhibit good mechanical properties as well as recycle-ability. This makes them an eco-friendly alternative to other plastics and metals [START_REF] Stewart | Thermoplastic composites-recyclable and fast to process[END_REF] .

Although the use of thermoplastic composites can be advantageous, their sensitivity to temperature and strain rates complicates the study of their mechanical behavior [START_REF] Schoßig | Mechanical behavior of glass-fiber reinforced thermoplastic materials under high strain rates[END_REF] . Such complex conditions can be encountered by the material during heat assisted forming such as hot stamping [START_REF] Greco | Thermoforming of thermoplastics matrix composites[END_REF][START_REF] Zhao | A Study on Hot Stamping Formability of Continuous Glass Fiber Reinforced Thermoplastic Composites[END_REF] and more particularly for heat assisted incremental sheet forming, already tested for metal sheets but relatively new for thermoplastics. The latter is drawing considerable interests and significant progress has been made recently [START_REF] Zhu | Incremental sheet forming of thermoplastics: a review[END_REF] . The sensitivity of the matrix to severe conditions of temperature and strain rates is expected to be a significant part of the overall mechanical behavior of a fiber reinforced thermoplastic polymer [START_REF] Nasraoui | Influence of strain rate, temperature and adiabatic heating on the mechanical behaviour of poly-methyl-methacrylate: experimental and modelling analyses[END_REF][START_REF] Atmani | Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plugassisted thermoforming[END_REF] . Atmani et al. [START_REF] Atmani | Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plugassisted thermoforming[END_REF] performed an experimental and numerical study to describe the behavior of an unreinforced thermoplastic polymer (high impact polystyrene HIPS) during plug-assisted thermoforming. In the experimental part of their study, they investigated the elastic properties by performing uniaxial tensile tests at different temperatures, from 80 • C to 140 • C, and at three different strain rates: 0.0167s -1 , 0.033s -1 , and 0.167s -1 . They noticed that the elastic modulus and the yield strength increased by 67% and 50%, respectively, as the strain rate increased by a factor of 10. Moreover, a temperature increment of 20 • C at the lowest strain rate, caused a decrease in the values of Young's modulus and yield strength by 67% and 54%, respectively. Furthermore, Wang et al. [START_REF] Wang | Effects of temperature and strain rate on the tensile behavior of short fiber reinforced polyamide-6[END_REF] , studied the effect of temperature (ranging from 21 • C to 100 • C) and strain rate (0.05 min -1 , 0.5 min -1 , 5 min -1 ) on the mechanical behavior of short glass fiber reinforced polyamide-6. They concluded that increasing the strain rate by a factor of 100 at a temperature of 21 • C, induces an increase of 16% and 18% for both Young's modulus and tensile strength, respectively. Whereas, an elevation of temperature from 21 • C to 100 • C causes a decrease of 75% in Young's modulus and 43% in the tensile strength. Duan et al. [START_REF] Duan | Experimental and numerical investigations of strain rate effects on mechanical properties of LGFRP composite[END_REF] , investigated the effect of uniaxial strain rate on the mechanical properties of glass fiber reinforced polypropylene at RT. They showed that increasing the strain rate from 0.001s -1 to 50s -1 increases the ultimate tensile strength by 84% while having a negligible effect on Young's modulus. Similarly, Cui et al. [START_REF] Cui | The effects of strain rates on mechanical properties and failure behavior of long glass fiber reinforced thermoplastic composites[END_REF] , studied the tensile behavior of long glass fiber reinforced polypropylene at room temperature at a strain rate ranging between 0.001s -1 and 400s -1 . Results showed that the ultimate strength, fracture strain, and Young's modulus increase by increasing the strain rate. The effect of strain rate on the mechanical behavior was also studied for the case of hybrid composites, i.e. reinforced with natural and synthetic fibers. Basrani et al. [START_REF] Basrani | An Experimental Study of Static and Dynamic Behavior of the Hybrid Composite[END_REF] investigated the influence of strain rate on a composite reinforced with four different natural fibers in addition to top and bottom layer of glass fibers. It was found that for strain rates ranging from 0.00022 s -1 to 0.0888 s -1 , the tensile strength of the composite increased by 8% with the increase of the deformation rate.

On the other hand, several studies have been carried out on polymers as well as their corresponding composites to model their mechanical behavior specifically for thermoforming applications. Mooney-Rivlin's model [START_REF] Kumar | Hyperelastic Mooney-Rivlin model: determination and physical interpretation of material constants[END_REF] is an example of hyperelastic models that can predict the stress flow of polymeric materials. Gong et al. [START_REF] Gong | An anisotropic hyperelastic constitutive model for thermoplastic woven composite prepregs[END_REF] developed an anisotropic hyperelastic constitutive model to predict the large deformation behavior of woven prepregs during thermo-forming process. Pham et al. [START_REF] Pham | Modeling and simulation of stretch blow molding of polyethylene terephthalate[END_REF] predicted the strain rate and temperature dependent stretching behavior of a thermoplastic polymer (polyethylene therephtalate PET) using a visco-hyperelastic model based on Mooney-Rivlin's model. They performed biaxial stretching tests using the "Bruckner stretcher" machine. The model parameters were identified using non-linear curve fitting of biaxial results in MATLAB, based on the Levenberg-Marquardt algorithm. Additionally, Guzman et al. [START_REF] Guzman-Maldonado | Thermomechanical analysis, modelling and simulation of the forming of pre-impregnated thermoplastics composites[END_REF] developed a non-linear viscohyperelastic model based on the generalisation of the Maxwell's rheological model. The objective is to model the forming of pre-impregnated thermoplastic composites. Also, Peng et al. [START_REF] Peng | A phenomenological thermal-mechanical viscoelastic constitutive modeling for polypropylene wood composites[END_REF] proposed a non-linear thermo-viscoelastic constitutive model, used in thermoforming simulation of a car interior part, to describe the strain rate and temperature dependent behavior of a polypropylene wood composite. The Cowper-Symonds viscoplastic model, originally designed for metals, was applied in some studies to model the behavior of fiber reinforced polymers [START_REF] Wang | Mechanical characteristics, constitutive models and fracture behaviors of short basalt fiber reinforced thermoplastic composites under varying strain rates[END_REF][START_REF] Chen | The compressive and tensile behavior of a 0/90 C fiber woven composite at high strain rates[END_REF] . Another popular viscoplastic model, G'Sell and Jonas, was mainly designed for modeling the mechanical behavior of solid semi-crystalline polymers where the authors studied the plastic behavior of poly(vinyl chloride) and high density polyethylene [START_REF] G'sell | Determination of the plastic behaviour of solid polymers at constant true strain rate[END_REF] . Consequently, many studies were inspired by this model to describe the behavior of different polymers such as as PMMA, HDPE, PP, HIPS [START_REF] Atmani | Identification of a thermo-elasto-viscoplastic behavior law for the simulation of thermoforming of high impact polystyrene[END_REF][START_REF] Sala | A numerical and experimental approach to optimise sheet stamping technologies: polymers thermoforming[END_REF][START_REF] Abbès | Experimental and numerical study of the aging effects of sorption conditions on the mechanical behaviour of polypropylene bottles under columnar crush conditions[END_REF][START_REF] Abbès | Modelling of the degradation of mechanical properties of highdensity polyethylene based-packaging exposed to amyl acetate solution[END_REF] . The numerical work of Atmani et al. [START_REF] Atmani | Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plugassisted thermoforming[END_REF] consisted in modeling the thermoforming of a HIPS yogurt container by an elasto-viscoplastic model based on G'Sell and Jonas. Morover, the latter gave acceptable results describing the mechanical behavior of fiber reinforced polymers. In the work of Schossig et al. [START_REF] Schoßig | Mechanical behavior of glass-fiber reinforced thermoplastic materials under high strain rates[END_REF] , the strain rate dependent behavior of glass fiber reinforced polypropylene and polybutene-1 was experimentally investigated and modeled using the phenomenological G'Sell and Jonas law. Concerning the mechanical characterization techniques, conventional uniaxial tests are widely applied to identify the mechanical behavior of fiber reinforced thermoset polymers such as epoxy. Less frequently, biaxial testing methods have also been employed with different types of specimen including cruciform shaped samples [START_REF] Kobeissi | Strength characterization of glass/epoxy plain weave composite under different biaxial loading ratios[END_REF][START_REF] Bhatnagar | Development of a biaxial tensile test fixture for reinforced thermoplastic composites[END_REF][START_REF] Kobeissi | Caractérisation mécanique des matériaux composites à partir d'un essai de traction biaxiale intégrant l'effet des chemins de déformation[END_REF][START_REF] Chen | A review of multiaxial/biaxial loading tests for composite materials[END_REF][START_REF] Hartmann | Basic studies in biaxial tensile tests[END_REF] . In this paper, the mechanical behavior of glass fiber reinforced polypropylene (40% GF/PP) is investigated. The combined effect of strain rate and temperature is considered. Normally, in the majority of the published studies related to the mechanical characterization fiber reinforced polymers, the effect of temperature is monitored at one strain rate (mainly quasi-static) and the strain sensitivity analysis is performed at room temperature. Therefore, three strain rates and three temperatures are considered in this study. For each temperature, uniaxial tensile tests are conducted for three strain rates. The anisotropic behavior is examined as well. Additional tests are performed using uniaxial specimens extracted at different orientations to examine the anisotropy of the material. A phenomenological constitutive model, inspired by G'Sell and Jonas model [START_REF] G'sell | Determination of the plastic behaviour of solid polymers at constant true strain rate[END_REF][START_REF] G'sell | Effect of stress triaxiality on neck propagation during the tensile stretching of solid polymers[END_REF] , originally designed for unreinforced semi-crystalline polymers is proposed. A modification is applied to the original equation for a better description of the behavior experimentally observed. The proposed law is then numerically validated by finite element simulations of the uniaxial tensile test on the commercial software Abaqus. In the last part of this paper, the model identified from uniaxial tests is used to simulate an equi-biaxial tensile test where the state of loading is supposed to be close to that encountered during the forming process of the material. The results presented in this article lay the groundwork for a more advanced characterization technique based on biaxial tensile tests on cruciform shaped specimens. The key purpose behind this investigation is to model the complex behavior of glass fiber reinforced polypropylene during heat assisted incremental sheet forming.

Material and experimental testing

Material

The composite material used in this work is a 40% glass fiber reinforced polypropylene thermoplastic (40%GF/PP). Injection molded rectangular plates, of dimension 510x300x2 mm 3 , are supplied by SABIC company (Stamax, 40YM240). Specimens for uniaxial tests are extracted by means of water-jet cutting machine in the 0 • , 90 • and 45 • directions, where the 0 • direction represents Prepared using sagej.cls the length of the plate. Figure 1 shows the dimensions of the uniaxial tensile specimen made of 40%GF/PP.

Uniaxial Tensile tests

Uniaxial tensile tests are performed on a hydraulic machine [START_REF] Liang | Thermo-viscoplastic behavior of AA6061 under dynamic biaxial loadings[END_REF] equipped with a heating system. In this work, load velocities of 0.5, 50 and 500 mm/s are tested, corresponding to initial strain rate values of 0.001, 1 and 10 s The experimental forces are measured by a load sensor mounted on the upper fixed crosshead. Digital Image Correlation (DIC) is applied for strain measurements. This approach is nowadays extensively applied to obtain displacement and strain fields from images of a loaded part. It is applicable not only for metals, but also for polymers such as polymethyl methacrylate (PMMA) [START_REF] Alavi | Experimental stress determination of blunt notches under combinations of modes I and II loading[END_REF][START_REF] Alavi | On the applicability of digital image correlation method in extracting the higher order terms in stress field around blunt notches[END_REF] as well as for fiber reinforced polymers e.g., Basalt fiber reinforced polypropylene and glass fiber reinforced polypropylene [START_REF] Duan | Experimental and numerical investigations of strain rate effects on mechanical properties of LGFRP composite[END_REF][START_REF] Wang | Mechanical characteristics, constitutive models and fracture behaviors of short basalt fiber reinforced thermoplastic composites under varying strain rates[END_REF] . To perform DIC in this work, a speckle pattern (gray scale) is applied to the surface of the object (a white layer is first applied and then sprayed with random black dots). A highspeed/resolution camera (PHOTRON FASTCAM NOVA S9) is fixed above the specimen along with a proper lighting system ensuring that no surface heating is caused by the light. GOM Correlate software was used for the correlation procedure. Relying on A DIC guideline [START_REF] Bigger | A Good Practices Guide for Digital Image Correlation[END_REF] , a facet size of 21 pixels and a point distance of 7 pixels (0.12mm/pixel) are applied. A virtual extensometer is employed to calculate the length change in the region of interest, i.e. calculating the average deformation of the specimen's gauge length as illustrated in Figure 2. The use of the virtual extensometer to calculate the average deformation makes the results less sensitive to the correlation parameters. The tensile properties including Young's modulus E, ultimate tensile strength σ ult and strain at failure ε f are determined based on the Standard ISO-527 for plastics [START_REF]Organization, I. ISO 527-1: 2019-Plastics-Determination of Tensile Properties-Part 1: General Principles[END_REF] . In the following tests, the effect of specimen orientation, strain rate and temperature are studied. It's important to state that the repeatability of the results is verified with three tests for each condition. An example is illustrated in Figure 3 showing the repeatability of the stress strain curves resulting from three tests performed at the same experimental conditions at RT and quasi-static strain rate (0.001 s -1 ). The error bars represent the standard deviation calculated from the three tests. Good repeatability is achieved for all tested conditions.

Effect of specimen orientation

In order to investigate the effect of the material orientation and examine the anisotropic behavior of the material, uniaxial tensile tests are conducted at RT and quasi-static strain rate of 0.001 s -1 , using specimens extracted from three different orientations: i) along the length of the plate (0 • orientation), ii) along the width (90 • orientation) and iii) with an angle of 45 • . Figure 4 gives the stress-strain curves for the three directions showing a clear difference in the material behavior. The lowest value of the stiffness is attributed to the 90 • direction. The stiffness at 45 • is represented by the curve located between the 0 • and 90 • direction curves. Moreover, the ultimate tensile strength decreases with the increase of the specimen orientation angle. The highest value of 91 MPa is reached for the 0 • direction. On the other hand, the elongation at failure shows Prepared using sagej.cls a slight increase when the angle increases; it reaches a value of 2.3% for the 90 • direction. This behavior is thought to be related to the manufacturing of the structure and the induced fibers orientation/distribution, given that the plates are rectangular and the injection port is located in the middle. The discontinuous fibers are more likely to be oriented along the longitudinal direction at the expense of the other. In this case, it is along the length of the plate (0 • direction). Table 1 lists the values of obtained tensile properties for the three orientations at room temperature.

Effect of strain rate

To study the effect of strain rate on the mechanical behavior, test specimens extracted from the 0 • direction are considered. The influence of the strain rate is clearly shown in Figure 5 where stress-strain curves are plotted for three different temperatures (20 • C, Figure 5a; 70 • C, Figure 5b; 120 • C, Figure 5c). Note that the stress-strain behavior of the material exhibits a limited linear behavior at the beginning followed by a non-linear response.

For each temperature, three different strain rates are tested: 0.001, 1, and 10 s -1 . As a result, a positive strain rate effect is observed. This is remarkable for the three applied temperatures. Looking at Figure 5a, it is clear that when considering 0.001s -1 as a reference strain rate, σ ult increases by 16% and 43% when the reference strain rate is multiplied by a factor of 10 3 and 10 4 , respectively. The resultant maximum stress is then 130 MPa while the stress corresponding to a strain rate of 0.001s -1 is 91 MPa. In the considered strain rate range, the strain at failure and the stiffness vary slightly. Indeed, compared to the lowest strain rate, a small increase of 7% and 14% is noticed for the elongation at failure ε f and Young's modulus E, respectively, at the highest strain rate. For temperatures T = 70 • C and T = 120 • C, Figure 5b and Figure 5c respectively, the strength and the stiffness also increase with the strain rate. However, despite the clear Prepared using sagej.cls difference between the curves corresponding to strain rates of 0.001 and 1 s -1 , the curves obtained at 1 and 10 s -1 are relatively close. This shows that the sensitivity of the tested reinforced thermoplastic material to the strain rate becomes low at high temperatures. Note that for a strain rate of 0.001s -1 , elongation at failure increases from 2% to 4.06% and 4.85% for T = 20 • C, T = 70 • C and T = 120 • C, respectively.

Effect of temperature

The thermoplastic materials can be molded as often as desired by cooling and reheating, as long as the material is not overheated. Consequently, it is important to study the effect of temperature on the mechanical behavior of the tested material. The uniaxial tests already performed in the previous section are presented in this section as a function of the temperature. For every strain rate, three stress-strain curves corresponding to the three tested temperatures (20 • C, 70 • C and 120 • C) are plotted in Figure 6.

Figure 6a, Figure 6b and Figure 6c show the stress-strain curves corresponding to strain rates of 0.001, 1 and 10 s -1 , respectively. According to these curves, the effect of temperature is noteworthy. For a quasi-static uniaxial test ( ε = 0.001s -1 ), the ultimate strength decreases from 91 to 25 MPa while the elongation at failure increases from 1.94% to 4.85%. Moreover, when increasing the temperature, the stiffness decreases, for the three tested strain rates.

Figure 7 shows the combined effect of temperature and strain rate on the tensile properties of the thermoplastic material for the 0 • direction.

Table 2 summarizes the different values of the tensile properties as a function of the velocity for the three tested temperatures such as the ultimate stress, the elongation at failure and Young's modulus which is defined as the slope of stress-strain curves in the deformation range 0.05%-0.25% according to ISO 527 [START_REF]Organization, I. ISO 527-1: 2019-Plastics-Determination of Tensile Properties-Part 1: General Principles[END_REF] .

Table 2. Effect of strain rate on the tensile properties at three different temperatures.

ε(s -1) E(MPa) σ ult (MPa) ε f (%) T = 20

Identification of thermomechanical model

Isotropic model

The experimental investigations performed in this work underlined the significant effect of strain rate and temperature on the mechanical behavior of discontinuous glass fiber reinforced polypropylene. In this section, the

σ(ε, ε) = K • exp hε 2 • (1 -exp(-W ε)) • εm (1)
Where:

• σ is the equivalent stress in MPa,

• ε is the equivalent strain,

• ε is the equivalent strain rate,

• m is the strain rate sensitivity, • K is a scaling factor in MPa,

• h and W are material parameters.

Behavior at low strain levels is represented by the term (1 -exp(-W ε)). The term exp hε 2 represents the stiffening of the material at large deformations [START_REF] Atmani | Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plugassisted thermoforming[END_REF] which is neglected (h = 0) in this work due to the low level of deformation at failure observed from the experimental results. Additionally, the effect of the temperature can be incorporated in the model equation by performing a temperature sensitivity analysis for the material parameters and expressing them as a function of temperature. Therefore, the constitutive equation considered in this work becomes (equation 2):

σ(ε, ε, T ) = K(T ) • (1 -exp(-W (T ) • ε)) • εm(T ) (2)
As proposed in many studies [START_REF] Atmani | Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plugassisted thermoforming[END_REF][START_REF] Atmani | Identification of a thermo-elasto-viscoplastic behavior law for the simulation of thermoforming of high impact polystyrene[END_REF][START_REF] Sala | A numerical and experimental approach to optimise sheet stamping technologies: polymers thermoforming[END_REF][START_REF] Abbès | Experimental and numerical study of the aging effects of sorption conditions on the mechanical behaviour of polypropylene bottles under columnar crush conditions[END_REF][START_REF] Abbès | Modelling of the degradation of mechanical properties of highdensity polyethylene based-packaging exposed to amyl acetate solution[END_REF] , the equation of the model can be modified to clearly identify the linear and nonlinear behaviors. Consequently, in this work, the constitutive equation of the model is modified as stated in equation 3:

   σ = E( ε, T ) • ε for σ < σ y σ = σ y ( ε, T ) +K(T ) • (1 -exp (-W (T ).ε p )) εm(T ) for σ > σ y
(3) The equivalent deformation and equivalent plastic deformation are represented by ε and ε p , respectively. The latter is determined from the experimental stress-strain curves based on equation 4:

ε p = ε - σ E (4)
Accordingly, the yield strength σ y is identified as the stress level at which plastic deformation initiates, indicating that ε p exceeds zero. The strain rate sensitivity m is determined in the non-linear region (equation 5) as the slope of the plot ln(σ -σ y ( ε, T )) vs. ln ε, at a given temperature T and a plastic strain level ε p , as stated in equation ( 6). The value of m is assumed to be constant for all strain levels; an average of m values is then taken at five different strain levels, as suggested by [START_REF] G'sell | Determination of the plastic behaviour of solid polymers at constant true strain rate[END_REF] :

σ -σ y ( ε, T ) = K(T ) (1 -exp (-W (T ).ε p )) εm (5) m = ∂ ln (σ -σ y ( ε, T )) ∂ ln ε ε p,T (6)
The parameters K and W are analytically determined by parameter regression using a nonlinear optimization function in MATLAB. The cost function to be minimized is defined by the error δ as follows,

δ = min n i=1 (σ model -σ experimental ) 2 (7)
Where n is the total number of data points, σ model is the predicted stress, and σ experimental is the stress calculated from the force measured experimentally.

Prepared using sagej.cls

Accepted manuscript

To generalize equation 3, Young's modulus E and yield strength σ y are approximated by the following relations,

E(T, ε) = E 0 (T ) εm (8) σ y (T, ε) = σ 0 (T ) εm (9)
The parameters of G'Sell and Jonas model identified for each temperature are summarized in Table 3.

Anisotropy modeling

It is clear from the stress-strain curves at 0 • , 45 • and 90 • that the anisotropy in this material cannot be neglected. As suggested by some authors [START_REF] Van Erp | Prediction of yield and long-term failure of oriented polypropylene: Kinetics and anisotropy[END_REF][START_REF] Abida | Anisotropic visco-elastoplastic modeling of quasiunidirectional flax fiber reinforced epoxy behavior: an investigation on low-velocity impact response[END_REF][START_REF] Amiri-Rad | An anisotropic viscoelastic-viscoplastic model for short-fiber composites[END_REF][START_REF] Amiri-Rad | Constitutive modeling of injectionmolded short-fiber composites: Characterization and model application[END_REF][START_REF] Joki | Coupling of plasticity and damage in glass fibre reinforced polymer composites[END_REF][START_REF] Pastukhov | Plasticity-controlled failure of fibre-reinforced thermoplastics[END_REF] , the equivalent stress ψ (σ) can be described by the conventional Hill 48 yield criterion [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] (equation 10):

ψ (σ) = F (σ 22 -σ 33 ) 2 + G (σ 33 -σ 11 ) 2 +H (σ 11 -σ 22 ) 2 + 2(Lσ 2 23 + M σ 2 + N σ 2 12 ) (10)
In the case of plane stress the equation reduces to:

ψ (σ) = (G + H)σ 11 2 + (H + F )σ 22 2 -2Hσ 11 σ 22 + 2N σ 12 2 (11)
The parameters F , G, H, and N are determined from the yield stress ratio R(θ) at a given orientation θ with respect to the 0 • direction (equation 12):

R(θ) = 1 F • sin 2 θ + G • cos 2 θ + H+ (2N -F -G -4H) sin 2 θ cos 2 θ = σ y (θ) σ y (0 • ) (12)
The plastic potentials in ABAQUS are defined as:

R 11 = 1 G + H , R 22 = 1 F + H , R 33 = 1 G + F , R 12 = 3 2N , R 13 = R 23 = 1 (13)
Table 4 shows the Hill's parameters representing the anisotropy of the studied material at RT and quasi-static strain rate. The same parameters are applied for all the tested conditions.

Results and Discussions

Numerical validation of the identified model

By applying the same uniaxial experimental conditions (homogeneous temperature and velocity), the calibrated model is numerically validated by performing finite element simulations of the uniaxial tensile tests on the commercial FE software ABAQUS. The model is implemented by FORTRAN code as a user subroutine UHARD. Only the gauge region of the specimen is modeled. Due to the homogeneous strain field, the model is not sensitive to the element size, number and type. Figure 8 shows the stressstrain variations obtained experimentally for 0 • orientation when compared with the simulation results. As shown in the figures, the experimental and numerical results are in good agreement. The mean relative error between experimental and simulation stress-strain curves is indicated in the figures.

On the other hand, the anisotropy is numerically validated by performing finite element simulation using Hill's yield criterion while integrating the parameters determined in the previous section. The three directions (0 • , 45 • and 90 • ) are then numerically simulated and compared to the uniaxial experimental results, in Figure 9. It shows that Hill's anisotropic criterion associated to the identified numerical model agrees well with the uniaxial experimental results.

Numerical validation of the identified model using biaxial results

As already discussed, for the future work, the identified phenomenological model will be used to simulate the forming of a plate made of 40% GF/PP, through the heat assisted incremental sheet forming process. The in-plane biaxial tensile test on cruciform specimens is a promising characterization tool in which the material can be loaded from shear to expansion, according to strain states similar to those encountered in forming processes. Therefore, it is worth checking the ability of the model identified from the uniaxial tensile tests to successfully model biaxial stress states such as those encountered by the material during the envisaged forming process. In the following section, a biaxial testing setup used to perform an in-plane equibiaxial tensile test on a 40%GF/PP thermoplastic cruciform specimen at room temperature and quasi-static condition is briefly presented. Then, experimental results are compared with numerical ones using the model identified from the uniaxial characterization.

Experimental biaxial tensile tests

In-plane biaxial tensile testing on flat cruciform specimens is a promising technique to evaluate the performance of sheet metal/composite components under various loading conditions. This loading configuration is representative of the stress state that the components experience in real-world applications [START_REF] Hartmann | Basic studies in biaxial tensile tests[END_REF][START_REF] Chen | Development and Experimental Verification of a High-Temperature and In-Plane Biaxial Testing Apparatus[END_REF] . In the frame of the biaxial tests, Liu et al. [START_REF] Liu | Identification of sheet metal hardening for large strains with an inplane biaxial tensile test and a dedicated cross specimen[END_REF] proposed a cruciform specimen to identify the hardening behavior of metallic sheets subjected to large deformations. As a first approach, this optimized shape of the cruciform specimen is considered in this study. Figure 10a shows the used cruciform specimen with its dimensions in Figure 10b. The specimens are machined from the same plates used for the uniaxial tests. From a 2mm initial thickness, a thinning Prepared using sagej.cls of the central zone is made (up to a final thickness of 0.625 mm) to ensure failure in this zone.

A servo-hydraulic biaxial tensile machine [START_REF] Kobeissi | Strength characterization of glass/epoxy plain weave composite under different biaxial loading ratios[END_REF] , with four independent actuators, is used to conduct the biaxial tensile tests. The force is measured on each axis of the machine using gauge sensors. To measure the strain, the same digital image correlation technique (DIC), presented in the uniaxial tests, is used. The high speed camera is placed above the central zone of the specimen. The equi-biaxial tensile test is performed at room temperature and at a tensile speed of 0.1mm/s imposed on each specimen arm. This loading induces an average equivalent strain rate in the central zone of approximately 0.003s -1 which is close to the strain rate range obtained from the uniaxial tensile test realized at RT Prepared using sagej.cls

Accepted manuscript

and for quasi-static condition. Figure 11 shows the force and the strain curves resulting from the biaxial tensile test. As shown in the figure, maximum forces of Fx = 1400 N and Fy = 1100 N are obtained at failure where x-axis is along the length of the plate i.e, 0 • and y-axis is the width of plate i.e, 90 • , respectively. The difference between the forces is due to the anisotropy of the tested material. Also, the obtained maximum strains at failure are ε f,x = 1% and ε f,y = 2%, corresponding for x-axis and y-axis, respectively. The equivalent Von-Mises strain ε M ises is found to have a maximum value of 3.1% at failure. This value is relatively higher than the deformation at failure attained during the uniaxial tensile test (strain of 2%) for the same experimental conditions.

Simulation of biaxial tensile tests

The validity of the phenomenological model identified from the uniaxial characterization is evaluated by the FE simulation of the experimentally performed in-plane equibiaxial tensile test. The material parameters outlined in Table 3 along with Hill's anisotropic potentials summarized in Table 4 are used to simulate the biaxial tensile test. As a first step, the test conducted at RT and quasi-static velocity is simulated. Due to the symmetry, only a quarter of the cruciform specimen is considered. Therefore, on each arm, half of the experimentally obtained forces are applied i.e. Fx/2 along x-axis and Fy/2 along y-axis. More details concerning the FE simulation of the biaxial test simulations are given in the work of Liu et al. [START_REF] Liu | Identification of sheet metal hardening for large strains with an inplane biaxial tensile test and a dedicated cross specimen[END_REF][START_REF] Liu | Identification of strain rate dependent hardening sensitivity of metallic sheets under in-plane biaxial loading[END_REF] . A comparison between the experimental and numerical equivalent Mises strain curves is shown in Figure 12. It is clearly obvious that simulating the biaxial test using the parameters identified from the uniaxial tests leads to poor agreement between experimental and numerical results. This discrepancy shows that the behavior of the material under biaxial loadings must be integrated in the calibration stage of the material model. Therefore, an inverse analysis based on the in-plane biaxial test is required for the re-identification of the model's parameters.

Conclusion

The mechanical behavior of glass fiber reinforced polypropylene is investigated, focusing on the combined effect of strain rate/temperature on the material response. A series of uniaxial tensile tests are performed at room temperature up to 120 • C, for strain rates ranging from quasi-static conditions to 10 s -1 . As suspected, tensile properties are very sensitive to temperature and strain rate. A phenomenological constitutive behavior law inspired by the G'Sell and Jonas model is identified. The original equation is modified according to the experimental observations of the uniaxial characterization. The model parameters are calibrated from the uniaxial results. Moreover, in order to take into account the anisotropy, the effect of the material orientation is examined by performing uniaxial tests at 0 • , 45 • and 90 • with respect to the length of the plate. The influence of the orientation is detectable on the stress-strain evolutions and tensile properties. Consequently, the results are utilized to calibrate the parameters of Hill's criterion. The model is then implemented in Abaqus to simulate the uniaxial tensile tests and good agreement is shown between experimental and numerical results. The validity of the identified model should be inspected when the material is subjected to a biaxial stress state as is the case during heat assisted forming. Therefore, as an introduction to future tasks, an equi-biaxial tensile test on a dedicated cruciform specimen is simulated at room temperature and quasi-static strain rate conditions using the identified uniaxial model. Numerical results are compared to the experimental ones showing that the predicted model does not fit with the actual biaxial behavior. This proves the importance of integrating the biaxial characterization in the model identification procedure. In future work, the material parameters of G'Sell and Jonas model will be recalibrated based on both uniaxial and biaxial tensile tests for temperatures and strain rates corresponding to the conditions experienced by the material during its forming stage.

Prepared using sagej.cls

Figure 1 .

 1 Figure 1. Dimensions of uniaxial tensile specimen (in mm).

  -1 , respectively. Experiments are carried out at three different temperatures: 20 • C, 70 • C and 120 • C. The room temperature RT (20 • C) is a conventional temperature of reference to test. Also, polypropylene is known to have a maximum continuous service temperature between 100 • C and 130 • C. For this reason, a temperature of 120 • C is selected in this range of maximum continuous operating temperature. An intermediate temperature of 70 • C is also tested. The temperature and strain rate ranges for the characterization are then directly linked to the conditions of heat forming processes.

Figure 2 .

 2 Figure 2. Strain measurement using virtual extensometer in DIC Correlate for the uniaxial tensile test at RT and 0.001 s -1 .

Figure 3 .

 3 Figure 3. Repeatability test conducted at RT and 0.001s -1 .

Figure 4 .

 4 Figure 4. Uniaxial stress-strain curves at room temperature and quasi-static strain rate for three different orientations.

( a )Figure 5 .

 a5 Figure 5. Effect of testing velocity on the behavior for three different temperatures.

Figure 6 .Figure 7 .

 67 Figure 6. Effect of temperature on the behavior for three different strain rates.

(Figure 8 .

 8 Figure 8. Comparison between experimental and numerical simulation of stress-strain curves at three different temperatures for 0 • orientation. The mean relative error between the experimental and simulation curves is also indicated.

Figure 9 .

 9 Figure 9. Comparison between experimental and simulation results for the three orientations 0 • , 45 • and 90 • taking into account the anisotropy (Hill's criterion) of the material along with the mean square error percentage.

( a )

 a Biaxial tensile test cruciform specimen designed by Liu et al.[START_REF] Liu | Identification of sheet metal hardening for large strains with an inplane biaxial tensile test and a dedicated cross specimen[END_REF][START_REF] Liu | Identification of strain rate dependent hardening sensitivity of metallic sheets under in-plane biaxial loading[END_REF] .(b) Dimensions of the cruciform specimen (in mm).

Figure 10 .

 10 Figure 10. Cruciform specimen used for biaxial tensile testing

Figure 11 .

 11 Figure 11. Force and strain curves for biaxial tensile test at room temperature and for a tensile test velocity of 0.1mm/s imposed on each specimen arm.

Figure 12 .

 12 Figure 12. Comparing the experimental equivalent Mises strain with the numerical equivalent Mises strain simulated based on the material parameters identified from the uniaxial characterization.

Table 1 .

 1 Tensile properties for three different material orientations at RT.

	ε f (%)

θ( • ) E (MPa) σ ult (MPa)

Table 3 .

 3 Material parameters determined for the three different temperatures.

	T ( • C) m(-)	K(M P a) W (-) σ 0 (M P a) E 0 (M P a)
	20	0.0347 80	419.5	34.96	8621
	70	0.074	51.3	290.7	19.9	5387
	120	0.07	27.3	212.9	13.2	3827

Table 4 .

 4 Parameters of Hill's plasticity for anisotropy.

	Direction Yield stress Yield stress Hill	Plastic
	θ( • )	σ y (M P a)	ratio R(θ)	coefficients potentials
	0 •	29.37	1	F = 2.7	R 11 = 1
	90 •	16.48	0.56	G = 0.53	R 22 = 0.56
	45 •	20.78	0.707	H = 0.47	R 33 = 0.6
				N = 2.38	R 12 = 0.79
					R 13 = R 23 = 1

Journal of Composite Materials XX(X)

Prepared using sagej.cls

Acknowledgments

This project has been funded by "Region Bretagne" of France. It was also supported by the Lebanese University. The authors also acknowledge "SABIC" company for providing the tested material.

Conflict of Interest

None declared.

Accepted manuscript


Étienne Payet

Binary Non-Termination in Term Rewriting and Logic Programming

Keywords: 2012 ACM Subject Classification Theory of computation → Constraint and logic programming, Theory of computation → Rewrite systems, Theory of computation → Program analysis Non-Termination, Term Rewriting, Logic Programming 1

We present a new syntactic criterion for the automatic detection of non-termination in an abstract setting that encompasses a simplified form of term rewriting and logic programming.

Introduction

This paper is concerned with non-termination in structures where one rewrites elements using indexed binary relations. Such structures can be formalised by abstract reduction systems (ARSs) [START_REF] Baader | Term Rewriting and All That[END_REF], i.e., couples (A, ⇒ I ) where A is a set and ⇒ I (the rewrite relation) is the union of binary relations on A, indexed by a set I, i.e., ⇒ I = {⇒ ι | ι ∈ I}. Non-termination in these structures can be formalised as the existence of an infinite rewrite sequence a 0 ⇒ ι0 a 1 ⇒ ι1 • • • . Term rewrite systems (TRSs) and logic programs (LPs) are concrete instances of ARSs: A is the set of finite terms and I indicates what rule (= a couple of finite terms) is applied at what position. A crucial difference is that the rewrite relation of TRSs relies on instantiation while that of LPs relies on narrowing, i.e., on unification. In this paper, we present a new syntactic criterion for the automatic detection of non-termination in an abstract setting that encompasses a simplified form of term rewriting and logic programming. Namely, we suppose that the rewriting always takes place at the root position of terms (see Def. 4 below). There exist program transformation techniques that make it possible to place oneself in such a context, e.g., the overlap closure [START_REF] Guttag | On proving uniform termination and restricted termination of rewriting systems[END_REF] in term rewriting or the binary unfoldings [START_REF] Codish | A semantic basis for the termination analysis of logic programs[END_REF][START_REF] Gabbrielli | Goal independency and call patterns in the analysis of logic programs[END_REF] in logic programming preserve the non-termination property of the original program.

Preliminaries

We let N denote the set of non-negative integers.

Binary Relations

If ⇒ and → are binary relations on a set A, then ⇒ • → denotes their composition.

We let ⇒ 0 be the identity relation and, for all n ∈ N,

⇒ n+1 = (⇒ n • ⇒). Moreover, ⇒ * = {⇒ n | n ≥ 0}
is the reflexive and transitive closure of ⇒. We formalise nontermination as the existence of an infinite sequence of connected elements:

▶ Definition 1. Let ⇒ be a binary relation on a set A. A ⇒-chain is a (possibly infinite) sequence a 0 , a 1 , . . . of elements of A such that a n ⇒ a n+1 for all n ∈ N. We simply write it as a 0 ⇒ a 1 ⇒ • • • .

Terms

We use the same definitions and notations as [START_REF] Baader | Term Rewriting and All That[END_REF] for terms. From now on, we fix a signature Σ (the function symbols) together with an infinite countable set X of variables, with Σ ∩ X = ∅.

We let f, g, s be function symbols of positive arity and 0 be a constant symbol. The set of all terms built from Σ and X is denoted by T (Σ, X). A context is a term with at least one "hole", represented by □, in it. For all terms or contexts t, we let Var(t) denote the set of variables occurring in t and, for all contexts c, we let c[t] denote the term or context obtained from c by replacing all the occurrences of □ by t. For all contexts c, we let c 0 = □ and, for all n ∈ N, c n+1 = c[c n ]. Terms are generally denoted by a, s, t, u, v, variables by x, y and contexts by c, possibly with subscripts and quotes. The set S(Σ, X) of all substitutions consists of the functions θ from X to T (Σ, X) such that Dom(θ) = {x ∈ X | θ(x) ̸ = x} is finite. A substitution θ is usually written as {x → θ(x) | x ∈ Dom(θ)} and its application to a term s as sθ. A renaming is a substitution that is a bijection on X. The composition of substitutions σ and θ is denoted as σθ. We say that σ is more general than θ if θ = ση for some substitution η. We let θ 0 = ∅ (the identity substitution) and, for all n ∈ N, θ n+1 = θ n θ.

A term s is an instance of a term t if s = tθ for some θ ∈ S(Σ, X). On the other hand, s unifies with t if sθ = tθ for some θ ∈ S(Σ, X); then, θ is called a unifier of s and t and mgu(s, t) denotes the most general unifier (mgu) of s and t.

Term Rewriting and Logic Programming

We refer to [START_REF] Baader | Term Rewriting and All That[END_REF] (resp. [START_REF] Apt | From Logic Programming to Prolog[END_REF]) for the basics of term rewriting (resp. logic programming). 2 , every element (u, v) of which is called a rule, where u (resp. v) is the left-hand side (resp. right-hand side). For each program P , we let P denote the set of all finite, non-empty, sequences of elements of P .

▶ Definition 2. A program is a subset of T (Σ, X)
In this paper, we only consider ARSs (A, ⇒ I ) such that A = T (Σ, X) and I is a program. Hence the following simplified definition. ▶ Definition 3. An abstract reduction system (ARS) is a union of binary relations on T (Σ, X) indexed by a program, i.e., it has the form ⇒ P = {⇒ r ⊆ T (Σ, X) 2 | r ∈ P } for some program P . For each ARS ⇒ P and each ω = (r 1 , . . . , r n ) in P , we let

⇒ ω = (⇒ r1 • • • ••⇒ rn ).
The next definition introduces term rewrite systems and logic programs as concrete instances of ARSs. For all terms s and rules (u, v) and (u ′ , v ′ ), we write (u, v) ≪ s (u ′ , v ′ ) to denote that (u, v) is a variant of (u ′ , v ′ ) variable disjoint with s, i.e., for some renaming γ,

we have u = u ′ γ, v = v ′ γ and Var(u) ∩ Var(s) = Var(v) ∩ Var(s) = ∅.
▶ Definition 4. For each program P , we let → P = {→ r | r ∈ P } and ⇝ P = {⇝ r | r ∈ P } where, for all r ∈ P ,

→ r = uθ, vθ ∈ T (Σ, X) 2 (u, v) = r, θ ∈ S(Σ, X) (Term Rewriting) ⇝ r = s, vθ ∈ T (Σ, X) 2 (u, v) ≪ s r, θ = mgu(s, u) (Logic Programming)
We say that → P (resp. ⇝ P ) is a term rewrite system (resp. a logic program).

▶ Example 5. Let r = f(x), s(x) = (u, v). Then, f 2 (x) → r s(f(x)) because f 2 (x) = uθ and s(f(x)) = vθ for θ = {x → f(x)}. Let r ′ = f(g(x, 0)), f(x) and s = f(g(x, x)). The rule (u ′ , v ′ ) = (f(g(x ′ , 0)), f(x ′ )) is a variant of r ′ variable disjoint with s. Let θ ′ = {x → 0, x ′ → 0}. Then, θ ′ = mgu (s, u ′ ) and we have s ⇝ r ′ v ′ θ ′ , i.e., f(g(x, x)) ⇝ r ′ f(0).
In term rewriting and in logic programming (modulo a condition), the left-hand side of a rule can be rewritten to the corresponding instance of the right-hand side.

▶ Lemma 6. Let r = (u, v) be a rule and θ be a substitution. We have uθ → r vθ and, if Var(v) ⊆ Var(u), uθ ⇝ r vθ.

Binary Non-Termination

We are interested in binary chains, i.e., infinite chains that consist of the repetition of two sequences of rules. There are ARSs that admit such chains but no infinite chain consisting of the repetition of a single sequence (see, e.g., → P in Ex. 8 and Ex. 9 below). More precisely: ▶ Definition 7. Let ⇒ P be an ARS and ω 1 , ω 2 ∈ P . A (ω 1 , ω 2 , ⇒ P )-chain is an infinite

(⇒ * ω1 • ⇒ ω2 )-chain.
▶ Example 8. Let ⇒ P ∈ {→ P , ⇝ P } where P is the program that consists of the rules

r 1 = f(x, s(y)), f(s 2 (x), y) r 2 = f(x, 0), f(s(0), x)
(see [START_REF] Zantema | Non-looping rewriting[END_REF] and TRS_Standard/Zantema_15/ex12.xml in [START_REF]Termination Problems Data Base[END_REF]). We have the (r 1 , r 2 , ⇒ P )-chain:

f(s(0), 0) 0 ⇒ r1 f(s(0), 0) ⇒ r2 f(s(0), s(0)) 1 ⇒ r1 f(s 3 (0), 0) ⇒ r2 f(s(0), s 3 (0)) 3 ⇒ r1 • • • ▶ Example 9.
Let ⇒ P ∈ {→ P , ⇝ P } where P is the program that consists of the rules

r 1 = f(x, s(y)), f(s(x), y) r 2 = f(x, 0), f(x, s(x))
(see [START_REF] Zantema | Non-looping rewriting[END_REF] and TRS_Standard/Zantema_15/ex14.xml in [START_REF]Termination Problems Data Base[END_REF]). We have the (r 1 , r 2 , ⇒ P )-chain:

f(0, s(0)) 1 ⇒ r1 f(s(0), 0) ⇒ r2 f(s(0), s 2 (0)) 2 ⇒ r1 f(s 3 (0), 0) ⇒ r2 f(s 3 (0), s 4 (0)) 4 ⇒ r1 • • •
Now, we present a criterion for the detection of binary chains. It is tailored to deal with specific sequences ω 1 and ω 2 that each consist of a single rule of a particular form. Intuitively, the rule r 1 of ω 1 and the rule r 2 of ω 2 are mutually recursive; in r 1 , a context c is removed from the left-hand side to the right-hand side while, in r 2 , c is added again. Ex. 8 and Ex. 9 are concrete instances, with c = s(□). This is formalised as follows.

▶ Definition 10. A recurrent pair for a program P is a pair (r 1 , r 2 ) ∈ P 2 such that

r 1 = f(x, c[y]), f(c n1 [x], y) and r 2 = f(x, s), f(c n2 [t], c n3 [x]) x ̸ = y Var(c) = Var(s) = ∅ t ∈ {x, s}
▶ Example 11. In Ex. 8, we have (n 1 , n 2 , n 3 ) = (2, 1, 0), c = s(□) and s = t = 0. In Ex. 9, we have (n 1 , n 2 , n 3 ) = (1, 0, 1), c = s(□), s = 0 and t = x.

We show that the existence of a recurrent pair leads to that of a binary chain (see Prop. 20), provided that property (1) below is satisfied. The rest of this section is parametric in an ARS ⇒ P and a recurrent pair (r 1 , r 2 ) for P as in Def. 10, with r 1 = (u 1 , v 1 ) and r 2 = (u 2 , v 2 ). We suppose that we have

∀θ ∈ S(Σ, X) (u 1 θ ⇒ r1 v 1 θ) ∧ (u 2 θ ⇒ r2 v 2 θ) (1)
As Var(v 1 ) ⊆ Var(u 1 ) and Var(v 2 ) ⊆ Var(u 2 ), by Lem. 6 both → p and ⇝ P satisfy (1).

For the sake of readability, we introduce the following notation.

▶ Definition 12. For all m, n ∈ N, we let f(m, n) denote the term f(c m [s], c n [s]).
Then, we have the following two lemmas. Lem. 13 states that r 1 allows one to iteratively move a tower of c's from the second to the first argument of f. Conversely, Lem. 14 states that r 2 allows one to copy a tower of c's from the first to the second argument of f in just one step.

▶ Lemma 13. For all m, n ∈ N, f(m, n) ⇒ n r1 f(n 1 × n + m, 0).
Proof. We proceed by induction on n.

(Base: n = 0) Here, ⇒ n r1 is the identity. Hence, for all m ∈ N, we have

f(m, n) ⇒ n r1 f(m, n), where f(m, n) = f(n 1 × n + m, 0). (Induction) Suppose that for some n ∈ N we have f(m, n) ⇒ n r1 f(n 1 × n + m, 0) for all m ∈ N. Let m ∈ N. Then, f(m, n + 1) = f(c m [s], c n+1 [s]) = u 1 {x → c m [s], y → c n [s]}. Therefore, by (1), we have f(m, n + 1) ⇒ r1 v 1 {x → c m [s], y → c n [s]} where v 1 {x → c m [s], y → c n [s]} = f(c n1+m [s], c n [s]) = f(n 1 + m, n). But, by induction hypothesis, we have f(n 1 + m, n) ⇒ n r1 f(n 1 × n + (n 1 + m), 0), i.e., f(n 1 + m, n) ⇒ n r1 f(n 1 × (n + 1) + m, 0). Finally, f(m, n + 1) ⇒ n+1 r1 f(n 1 × (n + 1) + m, 0). ◀ ▶ Lemma 14. For all m ∈ N, f(m, 0) ⇒ r2 f(m ′ + n 2 , m + n 3 ) where m ′ = 0 if t = s and m ′ = m if t = x. Proof. Let m ∈ N. We have f(m, 0) = f(c m [s], s) = u 2 {x → c m [s]}. Hence, by (1), we have f(m
, 0) ⇒ r2 v 2 {x → c m [s]}. If t = s then v 2 {x → c m [s]} = f(c n2 [s], c m+n3 [s]) = f(n 2 , m + n 3 ). If t = x then v 2 {x → c m [s]} = f(c m+n2 [s], c m+n3 [s]) = f(m + n 2 , m + n 3 ). ◀
We consider the following polynomials in the indeterminate i ∈ N. We define them in a mutually recursive way, which reflects the mutually recursive nature of r 1 and r 2 and hence facilitates the proof of the existence of a (r 1 , r 2 , ⇒ P )-chain (Prop. 20 below).

▶ Definition 15. We let

Π 0 (i) = n 2 and Π ′ 0 (i) = n 3 Π n+1 (i) = ∆ n (i) + n 2 and Π ′ n+1 (i) = ∆ ′ n (i) + n 3 for all n ∈ N where, for all n ∈ N, ∆ n (i) = 0 if t = s and ∆ n (i) = ∆ ′ n (i) if t = x ∆ ′ n (i) = iΠ ′ n (i) + Π n (i).
▶ Example 16. In Ex. 9, we have t = x and (n 1 , n 2 , n 3 ) = (1, 0, 1). Hence:

Π 0 (i) = n 2 = 0 and Π ′ 0 (i) = n 3 = 1 Π 1 (i) = ∆ 0 (i) + n 2 = ∆ ′ 0 (i) = iΠ ′ 0 (i) + Π 0 (i) = i Π ′ 1 (i) = ∆ ′ 0 (i) + n 3 = i + 1 Π 2 (i) = ∆ 1 (i) + n 2 = ∆ ′ 1 (i) = iΠ ′ 1 (i) + Π 1 (i) = i 2 + i + i = i 2 + 2i Π ′ 2 (i) = ∆ ′ 1 (i) + n 3 = i 2 + 2i + 1
The next lemma provides a simpler form of Π and Π ′ for the case t = s (the case t = x is more intricate).

▶ Lemma 17. If t = s then, for all n ∈ N, Π n (i) = n 2 and Π ′ n (i) = n 3 i n + n-1 k=0 (n 2 + n 3 )i k .
Proof. Suppose that t = s. Then, for all n ∈ N, ∆ n (i) = 0, so Π n+1 (i) = n 2 . As Π 0 (i) = n 2 also, for all n ∈ N we have Π n (i) = n 2 . Now, we prove that Π ′ n (i) = n 3 i n + n-1 k=0 (n 2 + n 3 )i k . We proceed by induction on n.

(Base: n = 0) We have Π ′ n (i) = n 3 = n 3 i n + n-1 k=0 (n 2 + n 3 )i k .
(Induction) Suppose that the property holds for some n ∈ N. We have Π

′ n+1 (i) = ∆ ′ n (i) + n 3 = iΠ ′ n (i) + Π n (i) + n 3 . But, as t = s, Π n (i) = n 2 and, by induction hypothesis, Π ′ n (i) = n 3 i n + n-1 k=0 (n 2 + n 3 )i k . So, Π ′ n+1 (i) = i(n 3 i n + n-1 k=0 (n 2 + n 3 )i k ) + n 2 + n 3 = n 3 i n+1 + n k=0 (n 2 + n 3 )i k . ◀ ▶ Example
(i) = 1 and Π ′ n (i) = n-1 k=0 i k for all n ∈ N.
Using Π and Π ′ , we define the set of terms A:

▶ Definition 19. We let A = {a n = f(Π n (n 1 ), Π ′ n (n 1 )) | n ∈ N}.
Now we prove the existence of the (r 1 , r 2 , ⇒ P )-chain

a 0 ( Π ′ 0 (n1) ⇒ r1 • ⇒ r2 ) a 1 ( Π ′ 1 (n1) ⇒ r1 • ⇒ r2 ) a 2 ( Π ′ 2 (n1) ⇒ r1 • ⇒ r2 ) • • • ▶ Proposition 20. For all n ∈ N, we have a n (⇒ Π ′ n (n1) r1 • ⇒ r2 ) a n+1 . Proof. Let n ∈ N. We have a n = f(Π n (n 1 ), Π ′ n (n 1
)). By Lem. 13 and Lem. 14,

a n Π ′ n (n1) ⇒ r1 f n 1 × Π ′ n (n 1 ) + Π n (n 1 ) ∆ ′ n (n1) , 0 ⇒ r2 f m, ∆ ′ n (n 1 ) + n 3 Π ′ n+1 (n1) where m = n 2 = Π n+1 (n 1 ) if t = s and m = ∆ ′ n (n 1 ) + n 2 = Π n+1 (n 1 ) if t = x. Hence, a n (⇒ Π ′ n (n1) r1 • ⇒ r2 ) a n+1 . ◀ ▶ Example

Future Work and Implementation

We plan to investigate how our work relates to the forms of non-termination detected by the approaches of [START_REF] Emmes | Proving non-looping non-termination automatically[END_REF][START_REF] Geser | Non-looping string rewriting[END_REF][START_REF] Wang | On non-looping term rewriting[END_REF]. We have no clear idea for the moment.

Our tool NTI (Non-Termination Inference) [START_REF]NTI (Non-Termination Inference[END_REF] is designed to automatically prove the existence of infinite chains in TRSs and in LPs. It first transforms the original program P into a program P ′ : for TRSs, it uses the dependency pairs combined with a variant of the overlap closure [START_REF] Payet | Guided unfoldings for finding loops in standard term rewriting[END_REF] and, for LPs, it uses the binary unfolding [START_REF] Codish | A semantic basis for the termination analysis of logic programs[END_REF][START_REF] Gabbrielli | Goal independency and call patterns in the analysis of logic programs[END_REF]. By [START_REF] Arts | Termination of term rewriting using dependency pairs[END_REF][START_REF] Codish | A semantic basis for the termination analysis of logic programs[END_REF][START_REF] Guttag | On proving uniform termination and restricted termination of rewriting systems[END_REF], non-termination of P ′ implies that of P . Then, it detects recurrent pairs (Def. 10), hence binary chains (Prop. 20), in P ′ .

▶

  21. In Ex. 8, we have Π n (i) = 1 and Π ′ n (i) = n-1 k=0 i k for all n ∈ N (see Ex. 18). We also have n 1 = 2 and the (r 1 , r 2 , ⇒ P )-chain: Example 22. In Ex. 9, we have Π 0 (n1 ) = 0, Π ′ 0 (n 1 ) = 1, Π 1 (n 1 ) = 1, Π ′ 1 (n 1 ) = 2, Π 2 (n 1 ) = 3, Π ′ 2 (i) = 4, . . . (see Ex. 16). We have the (r 1 , r 2 , ⇒ P )-chain:

  18. In Ex. 8, we have t = s and (n 1 , n 2 , n 3 ) = (2, 1, 0). Hence, by Lem. 17, we have Π n