Pauline Martin
email: pauline.martin@inrae.fr

Vincent Ducrocq

Philippe Faverdin

Nicolas N C Friggens

Nicolas C Friggens

Invited review: Disentangling residual feed intake-Insights and approaches to make it more fit for purpose in the modern context

Keywords: feed efficiency, biological processes, individual variation, resilience Invited

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

INTRODUCTION

Improving feed efficiency is of major interest for livestock production, and thus it is not surprising that there has been an increasing number of publications reporting efficiency results. This is especially the case for studies using the residual feed intake (RFI) method. The RFI method is designed to quantify between-animal variation in efficiency, and it is well suited to use in genetic evaluation and breeding schemes because it does not derive efficiency as a ratio between production and intake. However, as has been reviewed elsewhere [START_REF] Berry | Cell biology symposium: Genetics of feed efficiency in dairy and beef cattle[END_REF][START_REF] Kenny | Invited review: Improving feed efficiency of beef cattle -the current state of the art and future challenges[END_REF], there is a confusing array of RFI results, correlations, and heritabilities that are not easy to compare because of differences in the RFI models used, which often makes biological interpretation impossible.

The variability in feed efficiency may hide different biologic realities [START_REF] Cantalapiedra-Hijar | Review: Biological determinants of betweenanimal variation in feed efficiency of growing beef cattle[END_REF]. As shown in Figure 1, animals may vary in their digestive efficiency (i.e., the conversion of feed into nutrients) or in their metabolic efficiency (i.e., the conversion of nutrients into products). It can be envisaged that (in a given feeding environment) the same overall efficiency may be achieved by different combinations of digestive and metabolic efficiencies. Within that feeding environment, these differences may be of no consequence, but may well become important in other environments. For instance, an animal that has a higher metabolic efficiency because it partitions more of its energy to milk and less to maintaining body reserves may not be penalized in abundant nutritional environments, but will probably not cope well in nutritional environments with feed shortages where it needs its body reserves.

A further limitation of most RFI models is that they typically calculate RFI for a fixed time interval, usually a short one (typically a few months), and are not adapted to characterizing long-term RFI (the entire productive life) or time-trends in RFI. These longerterm measures are important in the context of selecting for sustainable efficiency. Indeed, gains in short-term efficiency have been associated with negative conse-quences in the longer term, for instance on fertility and length of productive life [START_REF] Vallimont | Short communication: Feed utilization and its associations with fertility and productive life in 11 commercial Pennsylvania tie-stall herds[END_REF][START_REF] Puillet | Disentangling the relative roles of resource acquisition and allocation on animal feed efficiency: Insights from a dairy cow model[END_REF]. As can be seen in Figure 1, favoring higher production can cause a trade-off with the other life functions that underpin animal robustness and lifetime efficiency. The issue of trade-offs is more fully discussed in [START_REF] Friggens | Review: Deciphering animal robustness. A synthesis to facilitate its use in livestock breeding and management[END_REF].

Accordingly, the aim of this paper is to first examine the RFI method from a biological and theoretical basis to clarify the key questions related to the interpretation of RFI results, before then tackling some of the more pragmatic and statistical considerations. Approaches that overcome the current limitations of RFI methods are identified, and specific focus is given to making the RFI method more adapted to the modern context of precision livestock data and to the need to improve efficiency in a sustainable way. Readers looking for full reviews of the existing RFI literature and its comparison with other efficiency metrics are referred to [START_REF] Berry | Cell biology symposium: Genetics of feed efficiency in dairy and beef cattle[END_REF] or [START_REF] Cantalapiedra-Hijar | Review: Biological determinants of betweenanimal variation in feed efficiency of growing beef cattle[END_REF].

BASIS OF THE RFI EQUATION AND KEY QUESTIONS

The RFI equation is generally constructed from the following base:

Intake = intercept + a × predictor_a + b × predictor_b + c × predictor_c + […] + RFI, [1]
where RFI is the residual of the model. Linear regression is used to determine, for a given group animals on a given diet, the values of the regression coefficients a, b, c... as well as the intercept, but typically with only 1 observation per animal. In theory, the residual (i.e., RFI) captures differences in efficiency. Animals with a negative value of RFI have an observed intake that is less than that predicted for the group and, as such, are deemed to be more efficient because they have eaten less for a given production level (or produced more for a given intake level). Following the same logic, animals with a positive RFI are deemed to be less efficient.

The different predictors (i.e., factors considered in the intake prediction) are energy sinks. These include processes that consume energy such as milk production, growth, and maintenance, or internal energy sources such as use of body reserves. Because it provides a good biological context for highlighting most of the key issues concerning RFI, the rest of this paper focuses on RFI estimation in lactating animals. The RFI equation raises several questions:

• What is actually contained in the residual? As a residual, it is expected to contain both the biological variability in intake not explained by the model and any errors in measurement or in the model structure. The biological part of the residual, which equates to true differences in efficiency, is the part we wish to isolate. Is it possible to distinguish these different components, and if so, how best to do this? Alternatively, can we reduce the error component of RFI? This is explored in the section "Nature of the residual." • What is the biological meaning of the intercept?

From a theoretical point of view, the RFI equation should not have an intercept. This is evident if we consider it on an energy basis because the sum of true energy outputs is always equal to the true energy input. However, in practice there may be reasons to include an intercept in this kind of model, which are discussed in more detail in the section "Nature of the residual" (and "Intake and digestive efficiency") below. • Which predictors for the intake should be included in the equation? In other words, which predictors do we consider as being "required or essential energy sinks" for an animal, and which of these are likely to be the major determinants of feed efficiency differences? When more predictors are included in the model, the model performs better in terms of properly accounting for resource usage. However, when more is explained with the equation, less variability remains in the residuals as well as a greater proportion of the residuals that is simply noise (the error). In addition, many predictors are not easy to measure (e.g., heat production, gestation). How do we include them in the equation? Where do we stop? This issue is detailed in the section "Nature of the predictors." • An issue that emerges when considering how to rationalize the predictors relates to how they change with time. As the animal moves through different physiological stages, the relative importance of different predictors in the RFI estimation will change, and this may impact betweenanimal variation in these changes. For example, animals vary in the extent and duration of body reserves mobilization in early lactation, which has significant consequences on their subsequent reproductive efficiency. Capturing these time-trends requires time-series measures of the predictors, intake, and an expansion of the classical RFI model to deal with these repeated measures. Having such measures offers considerable opportunities for quantifying the relative importance of the different biological components of RFI. Approaches to do this and the benefits of so doing are discussed in the section "Incorporating time-trends in predictors." • Another key question is: how do we consider and incorporate intake data? This is especially important if data from multiple feeds or feeding situations are being used. If we use kilograms of DM, then the coefficients a, b, c… are composite values of both digestive and metabolic efficiency (because they include the conversion from DM to energy equivalents). If net energy (NE) is used, then, in theory, the digestive efficiency disappears from the coefficients. However, as the conversion of kilograms of DM to megajoules of NE assumes a constant digestive efficiency for all, then any interindividual variability for digestibility would "pollute" the residual, and maybe the other coefficients. A reflection on digestibility and its inclusion into the equation is presented in the section "Intake and digestive efficiency."

As several abbreviations are used in the development of the ideas and equations, their meanings are listed in Table 1 for ease of understanding.

NATURE OF THE RESIDUAL

As previously said, the RFI is the residual of the equation (Eq. 1); therefore, it contains not only the feed efficiency variability, but also all of the errors and imprecisions in the data and the RFI model. The residual can be considered to consist of the following: errors in measurement (eim), errors in fitting (eif), unidentified energy sinks (ues), deviations from the assumed average digestibility (ded), individual variability in metabolic efficiency from the identified energy sinks (imf), and a possible nonexplained part (ε i ). We will start by ignoring ded because this is more conveniently treated on the left-hand side of the RFI equation, which contains information about intake. Thus, the equation for the residual is as follows:

Residual = eim + eif + ues (+ ded) + imf + ε i , [2]
where the properties of each of these components of the residual are discussed below.

The eim should be random noise with a population mean of zero, but there may be a nonzero mean if there is a global bias in (one of) the measures. There may be individuals with systematic measurement errors (e.g., a cow who systematically spills more feed than the rest of the group or who systematically steals feed from others). The eif would be due to 1 of 3 issues. The first would be ignoring a nonlinear response to 1 of the predictor energy sinks, which implies variable efficiencies for that process. This is more complex to deal with and should only be considered when the possibilities below have been discounted. The second issue is an unclean calculation of an energy sink that then gives the illusion of a variable efficiency. An example of this, discussed more in the following section, is using the change in body weight (ΔBW) as an energy sink rather than separating it into the underlying changes in body lipid and body protein. Such a separation would then resolve the apparent eif. This example is really due to an ues (changes in body lipid and body protein are not explicitly considered as energy sinks), and it is to be expected that approximately all eif issues can be reconsidered as ues. The third issue is real eif (i.e., efficiencies differ according to performance level). These could result in a nonzero population mean of eif and may include individual differences (e.g., high vs. low producers). Real eif would imply a need for nonlinear RFI models. However, as shown subsequently, the 2 main sources of eif identified in the literature-an effect of intake level on digestibility and changing maintenance requirements as the proportion of metabolically active tissues increases with performance level -can be accommodated without the need for nonlinear models.

The ues can give a nonzero population mean (e.g., when ignoring animal activity levels) and may include individual differences (e.g., an animal with systematically a greater activity than the population mean). Given the technical difficulty of measuring all energy sinks, ues cannot be ignored and may include some "true" components of feed efficiency. For the imf, the individual variability in metabolic efficiency is what is sought to be captured with RFI; it is a "true" component of feed efficiency. The population mean per definition is zero. Given these considerations, it is now possible to look at the residual in terms of these components and specifically identify systematic elements (mean values, indicated by a bar) and animal-specific elements (indicated by subscript i). Thus, the residual is defined as follow:

Residual = + + + + + + + eim eif ues eim eif ues imf i i i i ε i i .
[3]

By definition, the mean value of imf is zero. By including Eq. 3 in Eq. 1, a new intercept is created equal to , eim eif ues + + and the residual term is now centered on zero:

Intake 2 = + × + × + + + + eim eif ues a predictor b predictor eim [...] 1 i i i i i i eif ues imf + + + +ε , [4]
Thus, the intercept adjusts for systematic errors in the RFI data and equation based on a mean value common among all individuals. Its presence is not mandatory (we hope for a zero intercept), but it may help to deal with reality, and, as such, can be useful to provide a meaningful term for these unaccounted systematic effects.

It is interesting to differentiate the components of RFI with a subscript i, and in particular to isolate the "true" individual variation in efficiency (i.e., ded i and imf i ). However, in practice, it is hard to see how to do this (the special case of ded i is dealt with in the "Intake and Digestive Efficiency" section), and the appropriate alternative approach would be considering them all in one. Considering them all in one assumes that all these individual deviations can be considered, in practice, as elements of global differences in efficiency between individuals. This is pragmatic, but may be dangerous if the size of these elements is affected by the local production environment. For example, ues i may include the costs of feeding behavior or activity. These will clearly differ in extensive grazing versus indoors feeding situations, and would probably be considered as beneficial to the efficiency of the grazing animal but not to the housed animal. In the same way, an animal that is systematically good at sorting its diet (resulting in its intake having a higher energy density than that of an animal on the same feed who does not sort) could be considered as being more efficient when the environment provides a heterogeneous feed (i.e., with an opportunity to sort), but would lose this advantage when on homogeneous feeds. Thus, even if no differentiation is done between the components of the residual in RFI calculation, it remains important to keep them in mind together with a detailed description of the environment (and the measurement methods). This applies not just to the interpretation of animal variability in RFI, but also to how farm, year or experiment effects are fitted in these models.

NATURE OF THE PREDICTORS

The predictors of intake in our equation are, by definition, all energy sinks (or supplies in case of body reserve mobilization) that explain changes in intake. If some appear obvious to consider, such as milk production, growth, or body reserves changes, it is not as simple to consider other predictors such as activity and heat diffusion, or other effects such as fertility needs. To deal with this fact, we can propose that energy sinks can be of 2 types: those that are part of production (e.g., milk, meat, fetus) and those that are part of fighting entropy (e.g., maintaining the machinery, protein turnover, immune function). In both cases, the benefits of including additional predictors should be weighed against the methodological difficulties of measuring them and the associated measurement errors that they may bring to the estimation of RFI.

From a biological perspective, it is a useful simplification to postulate that every production is either lipid or protein. This is because the underlying biological efficiencies of protein and lipid production are quite different from each other. Therefore, lipids and proteins can be divided into milk lipid and milk protein (considering that milk lactose is tightly linked to milk protein, and thus can be merged with the protein cost), growth lipid and growth protein, and body lipid reserves and body protein reserves. The second possible simplifying assumption is that the cost of creating 1 g of lipid is the same, regardless of the production it is associated with (e.g., milk, growth). The assumption is also to a large extent the same for protein (we ignore possible minor differences in costs between growth and milk production), but it has been known for a long time that the efficiency of synthesis is much lower for proteins than for lipids (e.g., 0.45 vs. 0.75 in rats; [START_REF] Pullar | The energy cost of fat and protein deposition in the rat[END_REF]and 0.35 vs. 0.60 in cattle;[START_REF] Geay | Energy and protein utilization in growing cattle[END_REF].

Using the assumption above, we accept that a production cost corresponds to the energy of the production itself (39.6 MJ/kg of lipid and 23.8 MJ/kg of protein) plus the work of building and disassembling the product (i.e., the efficiencies of these processes). However, as the work of building and disassembling molecules have different costs (for lipid and protein in ruminants see Friggens et al. 2007b), the equation should differentiate between accretion ⊕ ( ) and catabolism ( ) as follows:

Intake lipids lipids p roteins proteins m = × + × + × + × + ⊕ ⊕ a b c d a aintenance RFI. + [5]
For a product such as milk, a well-established approach to deal with its variation in lipid and protein content is to express milk as fat-and protein-corrected milk (i.e., adjust milk volume to a standard fat and protein content). Alternatively, and in coherence with the above postulate, we could use milk fat yield and milk protein yield in Eq. 5. For growth, where it is much more difficult to precisely measure the composition of growth, it is useful to use BCS in combination with BW to indirectly obtain body protein change and body lipid change. Other methods exist (e.g., calorimetry), but only on limited data. Changes in body protein and especially body lipid will also capture usage of body reserves (by definition, catabolism only occurs for body reserves). The details for doing this are developed elsewhere (e.g., Friggens et al., 2007b), but this boils down to having changes in BW and BCS (ΔBW and ΔBCS), and their interaction in the following equation:

Intake B W B W B CS BCS B W BCS BW B = × + × + × + × + × × + × × ⊕ ⊕ ⊕ a b c d e f " " " " " " " " C CS g MLY h MPY maintenance RFI, + + + + . . [6]
where MLY = milk fat yield, MPY = milk protein yield, and a through h are regression coefficients. Equation 6 raises a problem with respect to its resolution, as 6 different terms need to be estimated from only 2 different measures (BW and BCS). This will cause obvious correlations among the predictors and make resolution difficult, if not impossible. In addition, for any given time interval between 2 consecutive measures, it is not possible to distinguish what was gained and what was lost in that interval, with only the net difference (Δ) between those 2 being accessible. In classical RFI approaches, this will very likely lead to an unwanted ues in the residual, affecting the RFI estimates as it has been shown that there are differences between individual cows in their patterns of body reserve change [START_REF] Banos | Genetic profile of total body energy content of Holstein cows in the first three lactations[END_REF][START_REF] Kessel | Individual variability in physiological adaptation to metabolic stress during early lactation in dairy cows kept under equal conditions[END_REF]. This problem is revisited in the section "Incorporating time-trends in predictors."

With respect to the second component, fighting entropy, the only function routinely considered in the literature is maintenance. Maintenance is usually included through the metabolic weight. But in terms of energy, similar to production, the maintenance costs are different between lipids and proteins; therefore, as a minimum, the maintenance should not only be a function of the metabolic weight but also of the BCS. The equation for adjusting BW to a standard body fat content is

EBW EBW BCS BCS std std = × - - ( ) - - ( ) 1 1 ² ± ² ± , [ 7
]
where EBW is empty body weight (i.e., the animal weight without the digestive content), and std implies a given standard BCS (e.g., 2.75). The assumptions needed to do this (linear relation between BCS and body lipid content) are outlined in Friggens et al. (2007b). In Eq. 7, the coefficients α and β are assumed to be known but depend on the mature size of the breed and on the condition score scale used [START_REF] Zygoyiannis | Estimation of the mature weight of three breeds of Greek sheep using condition scoring corrected for the effect of age[END_REF][START_REF] Roche | Relationships among international body condition scoring systems[END_REF]. Typical values estimated from dairy cattle data are α = 0.5 and β = 0.1 (Friggens et al., 2007b). If the reasonable assumption is made that gutfill is a constant proportion of empty BW, then the equation above is exactly the same for BW. Thus, the predictor for maintenance, BW 0.75 , in the RFI equation can be replaced by BW std 0.75

, thereby eliminating spurious variation in maintenance requirements due to differences in body fatness.

There is also an issue related to maintenance costs associated with "more active than average" tissues or organs, which emphasizes the importance to finding ways to estimate the size of these tissues relative to the whole animal. Maintenance includes the maintenance of the digestive tract, but metabolic BW does not account for the changes in digestive tract size happening throughout the lactation nor the changes in digestive protein secretion. The intake capacity of the cow increases significantly during the first weeks of lactation (i.e., the size of the digestive tract increases), which, as a highly active tissue, implies an increase of maintenance requirements [START_REF] Tulloh | Physical studies of the alimentary tract of grazing cattle[END_REF][START_REF] Smith | Effects of breed, pregnancy, and lactation on weight of organs and tissues in dairy cattle[END_REF]. In practice, this issue is in some way going to be linked to intake levels. This will be dealt with below (section "Intake and Digestive Efficiency").

The question arises: what other aspects of "fighting entropy" could be included (or to what extent they are indexed on maintenance)? From a physiological perspective, it would seem reasonable to include variables such as immunity costs, resilience mechanisms, and activity. However, there is currently almost no quantitative basis for including these specific components, with the possible exception of activity. The extent to which differences in RFI can be explained by accelerometer and body temperature measures can be approached from 2 perspectives, purely in terms of energetics or as proxies for resilience. From the energetics perspective, differences in overall activity will be associated with a direct energy cost of movement, and, likewise, differences in the mean temperature of animals reflect a higher thermodynamic cost. As such, activity and temperature measures could be expected to explain a part of the variance in RFI. From the resilience perspective, highfrequency time-series measures of activity, milk yield, BW, and other indicators such as animal temperature are currently being explored as proxies to estimate the numbers of perturbations an animal experiences during a given time period [START_REF] Adriaens | Productive life span and resilience rank can be predicted from on-farm first-parity sensor time series but not using a common equation across farms[END_REF][START_REF] Poppe | Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding[END_REF]. Assuming that the animal responses to these perturbations have an energetic cost for the animal (i.e., activating the underlying physiological and immunological resilience mechanisms), then it is expected that such proxies (e.g., number of perturbations, cumulated loss, residual standard deviation of production records) may be usefully incorporated into efficiency considerations [START_REF] Berghof | Opportunities to improve resilience in animal breeding programs[END_REF].

Although this offers avenues for including resilience aspects in the efficiency method, it will not cover all aspects as is. For instance, the potential negative effect of increased production on probability of successfully reproducing would not be captured. Intuitively, it does not seem reasonable to directly include these in RFI. However, if RFI is measured over long enough periods, a penalty emerging via those animals that drop out of the herd early should be seen. In other words, resilience comes into the picture via productive lifespan when considering long-term efficiency. With respect to short-term efficiency, especially the case of growing animals reared for meat production, it is harder to see how to incorporate the notion of resilience. The question is important because the aim is to avoid selecting traits that favor short-term efficiency to the detriment of resilience. The issues of resilience and long-term efficiency are discussed in the final section of this paper.

INCORPORATING TIME TREND IN PREDICTORS

Because RFI requires individual data of feed intake, usually expressed as DMI, most of the pre-existing discussion on the time period over which RFI is estimated has focused on how short a time period can be used. These discussions typically pit the negative effect of a shorter period of measure on the accuracy of estimates against the high cost of phenotyping DMI for a longer period. The repeatability of intake measurements is indeed not that high, often ranging from 0.5 to 0.7 in early lactation for weekly measure in dairy [START_REF] Potts | Residual feed intake is repeatable for lactating Holstein dairy cows fed high and low starch diets[END_REF][START_REF] Li | Genetic parameters for dry matter intake in primiparous Holstein, Nordic Red, and Jersey cows in the first half of lactation[END_REF][START_REF] Connor | Defining the optimal period length and stage of growth or lactation to estimate residual feed intake in dairy cows[END_REF]. We can consider from the literature on beef cattle that 5 wk of measurement is a minimum for intake data [START_REF] Archer | Optimum postweaning test for measurement of growth rate, feed intake, and feed efficiency in British breed cattle[END_REF][START_REF] Archer | Duration of performance tests for growth rate, feed intake and feed efficiency in four biological types of beef cattle[END_REF][START_REF] Wang | Test duration for growth, feed intake, and feed efficiency in beef cattle using the GrowSafe System[END_REF], with a consensus that a measurement period of around 70 d is a reasonable length for RFI, corroborated by dairy studies [START_REF] Connor | Defining the optimal period length and stage of growth or lactation to estimate residual feed intake in dairy cows[END_REF]. Furthermore, intake data are not the only sensitive measurements that have an effect on the RFI. Body weight is also a key and tricky measurement. Depending on the feed inside the digestive tract, the live weight may vary up to 40 to 50 kg, enough to hide real changes in body condition, and hence body reserves.

Within the classical RFI model, there are also biological considerations for limiting the measurement period. If the measurement period is sufficiently short to be able to assume that it is not containing a mix of body accretion and mobilization, then the problem of Eq. 6 having too many correlated predictors can be addressed through the following simplification: where BW SIGN is a 0/1 variable, which is 1 when ΔBW is negative, and BCS SIGN is the same but applied to ΔBCS. Clearly, if the intervals over which these deltas are measured are short enough, then in any one interval, there will be no significant mixing of gain and loss. Assuming this was feasible (i.e., that measurement error is adequately dealt with), it allows to simplify further: However, the results suggest that with manual recording of BW and BCS, typically on a weekly or even monthly basis, it is not realistic to expect that the estimates of ΔBW, and especially ΔBCS, are sufficiently accurate for this approach. However, if high-frequency automated BW and BCS recording is done (e.g., at each milking; [START_REF] Fischer | Rear shape in 3 dimensions summarized by principal component analysis is a good predictor of body condition score in Holstein dairy cows[END_REF], then the approach becomes feasible. The logical conclusion of this is to describe these coefficients as function of time (e.g., dependent on the days in milk).

Intake = a′ × ΔBW + b′ × ΔBCS + other predictors + RFI, [ 9
Because precision livestock technologies provide high-frequency repeated measures, it becomes possible to move from classical methods where RFI is a pointmeasure for a given period to methods that explore the differences between animals in their trajectories of RFI through time. For example, [START_REF] Martin | A new method to estimate residual feed intake in dairy cattle using time-series data[END_REF] have recently shown that animal variation in RFI trajectories over the entire lactation can be derived from a multitrait random regression model using repeated measures of milk yield, BW, BCS, and DMI (in this case described by third-order Legendre polynomials). As shown in Figure 2, this enables the estimation of the correlations between the predicted intake and each of its predictors (milk, BW, and BCS) at each time point. Therefore, it is possible to have an indication of the relative importance of the different energy sinks in RFI at different stages of lactation. Because of the smoothing of the time-series within the model, and the continuous daily estimates of RFI, the whole problem of the length of measurement periods is circumvented, as is the problem of mixing body reserve accretion and mobilization. Considering the massive changes expected in early lactation, it could be more efficient to use splines rather than Legendre polynomials. The number of knots and degrees of polynomials between each knot would need to be defined, but there are precedents for this (Friggens et al., 2007a). Models of this type also open the door to a better biological interpretation of RFI. Classical RFI allows identification of animals that are more, or less, efficient (assuming that the errors described in the section "Nature of the residual" are properly dealt with), but cannot shed light on which biological components of the efficiency complex contribute to differences in feed efficiency. The paper of [START_REF] Fischer | Isolating the cow-specific part of residual energy intake in lactating Martin et al[END_REF], building on previous studies [START_REF] Aggrey | Dissection of Koch's residual feed intake: Implications for selection[END_REF]on birds;[START_REF] Savietto | Towards an improved estimation of the biological components of residual feed intake in growing cattle[END_REF]on beef), explored a new method to isolate the cow-specific part of residual energy intake from the residual. It consists of including a random component to the coefficients in the RFI equation, and thus capturing the interindividual variation in efficiency. In theory, models of this type, which exploit the repeated measures structure of precision livestock data, can quantify which biological components are underpinning differences in efficiency. Although the approach is promising, it has shown itself to present some statistical challenges with difficulties to fit the equation because of correlations among the factors. The papers listed above have been able to overcome some of the issues by normalizing the predictor measures, by exploiting different correlation structures, and by using more powerful software. [START_REF] Fischer | Isolating the cow-specific part of residual energy intake in lactating Martin et al[END_REF] used SAS, which is not able to easily take correlations into account; [START_REF] Martin | A new method to estimate residual feed intake in dairy cattle using time-series data[END_REF] used Wombat [START_REF] Meyer | WOMBAT-A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML)[END_REF], which allows the user to better estimate them, and eventually restrict the rank of the matrix. However, these methods are also data hungry, requiring large numbers of animals and large numbers of measures per animal.

INTAKE AND DIGESTIVE EFFICIENCY

The base measure of intake is DMI per day (we assume here that the methods for converting fresh matter intake into DM do not introduce bias in the RFI calculations; [START_REF] Seymour | Invited review: Determination of large-scale individual dry matter intake phenotypes in dairy cattle[END_REF], but of course when efficiency is considered, consideration should be given to the energy density of the diet. If the aim is to pool data sets from different trials or locations-this would certainly be needed in the context of large-scale genetic evaluation of RFI-then the differences in energy density of the feed need to be accounted for. The literature offers ideas on how to deal with this issue [e.g., [START_REF] Tempelman | nongenetic variation and energy sink relationships for residual feed intake across research stations and countries[END_REF] and [START_REF] Lu | Modeling genetic and nongenetic variation of feed efficiency and its partial relationships between component traits as a function of manage-ment and environmental factors[END_REF] use random ration-specific regressions of DMI on milk energy and BW], but to get a comprehensive understanding of the point, let us go into details. Even within a given situation where the feed is the same for all animals, there are important issues for the RFI method. Energy intake can be described in terms of gross energy, digestible energy, metabolic energy, and NE as follows [see [START_REF] Emmans | Effective energy: a concept of energy utilization applied across species[END_REF] for full details]:

Gross energy (GE) intake = intake × heat of combustion of the feed, Digestible energy intake = GE intake -gas -feces, Metabolic energy intake = GE intake -gas -feces -urine, Net energy intake = GE intake -gas -feces -urine -heat.

We can have individual variation in each step of the conversion from one energy level to another, and thus differences in (digestive) efficiency already present (see [START_REF] Phuong | Factors affecting energy and nitrogen efficiency of dairy cows: A meta-analysis[END_REF]. If the feed is described in terms of NE, then it is in the currency that is directly used by the processes that are the energy sinks included as predictors in the RFI equation. This means that the resulting estimated coefficients (a, b, c in Eq. 1) can be interpreted as the metabolic efficiencies for the different energy consuming processes. However, using NE, it is very difficult to isolate any between-animal variation in digestive efficiency. If, on the other hand, the feed measures used in the RFI equation are DMI, then each predictor coefficient becomes a combination of the metabolic efficiency of that process and the overall digestive efficiency. This can make interpretation tricky. It is for these reasons that up until this point in the paper, the left-hand side of the equation has always been called "intake" without any detail.

To deal with this, the effect on the RFI equation of the transformation of the intake into a given energy unit needs to be considered (i.e., the accounting for digestibility and metabolizability). Generally, this value, called energy density (ed), usually comes from literature based on the feed composition or digestibility measures. A priori, it should be assumed that this number is incorrect in 3 ways:

• It can be systematically off for a given feed [i.e., a mean deviation ( ded ) exists from ed]. • It ignores true individual variability in digestive efficiency (i.e., ded i , which is a "true" component of feed efficiency). It should be noted that the individual variability in digestive efficiency may be due to variation in the digestive processes within the digestive tract, but may also come from systematic differences in an animal's ability to select the most digestible parts within the feed offered (i.e., in feed sorting). It is considered here that this variation in feed sorting is one element of the animals feed efficiency. • It ignores the effect of level of intake on digestive efficiency. Let's call this ded FI . When expressing intake as digestible energy intake, this term is, in theory, considered by the conversion into digestible energy. However, assuming ded FI to be zero would be to assume that the digestible energy conversion is perfect across all feeds, which is not done here.

Taking the above into account, the energy density of the feed for a given animal becomes ed ded ded DMI ded

FI i i , + + × + ( )
and thus, when expressing intake in energy terms, the left-hand side of Eq. [1] becomes:

ed ded ded DMI ded DMI FI i . + + × + ( ) × [10]
Given that the aim of RFI is to isolate the individual variation in "true" efficiency (i.e., imf i and ded i ), the rest of this section will first deal with this aspect within a uniform feeding situation, and then will look at possibilities to extend this for multiple feeding situations. In both cases, the aim is that the RFI equation includes the estimation of the energy density of the feed, as well as isolating ded i . Accordingly, DMI can be split into the mean value (indicated with a bar) plus the individual deviations (i.e., DMI

D MI DMI

i = +
) Thus [10], the left-hand side of Eq. 1, becomes:

ed ded ded DMI ded DMI ed ded ded DMI FI i FI + + × + ( ) × + + + × . + ( ) × ded DMI i i [11]
It would be logical to expand the DMI part of ded FI × DMI in the same way. However, let's simplify by assuming that the feeding level effect on ed is sufficiently accounted for by the correction on the mean DMI (this effectively means any individual differences in ded FI are moved into the individual component of ded i , a considerable way is covered toward isolating individual variability in digestive efficiency. For example, a difference between 2 cows in ded i of 0.025 (relative to a typical ed true of 0.66) would give a difference in DMI i of 1 kg of DM/d (for cows eating on average 21 kg of DM/d with identical BW and milk production). To our knowledge, this reformulation of the RFI equation has not been reported in the literature, and thus it remains to be seen if it has practical utility. It should be noted that the above equation includes an animal performance component in the true digestibility (ed true ), the effect of the mean DMI on digestibility ded DMI FJ × (

) . To avoid this, we could redefine ed true as being only ed ded + . This does not change Eq. 15 and the resolution described above; it simply means that the individual variation in digestibility (ded i ) accounts for the variation in DMI in absolute terms rather than as a deviation from the mean DMI.

The second issue deals with the fact that different experiments will use different feeds, and thus the energy density of the feed, ed, and the error associated with its estimation will differ between experiments, as illustrated by the results of [START_REF] Tempelman | nongenetic variation and energy sink relationships for residual feed intake across research stations and countries[END_REF]. Ignoring this could create the situation where the regression coefficients estimated for the pooled data set are markedly different from those estimated within each feed or feeding situation (Simpson's paradox). Clearly, RFI models should be structured to avoid this problem. The idea is to include an element that allows the true value of ed, ed true , to be estimated for each feed or feeding situation within the RFI model. This may be possible, in a relatively straightforward way, by extending the random regression model of [START_REF] Martin | A new method to estimate residual feed intake in dairy cattle using time-series data[END_REF] to include the "experiment" factor among the fixed effects. This relies on the assumption that the average differences in ed can be removed by a linear adjustment, as suggested by the term ed ded DMI

true i + ( ) ×-
in the above equation. However, this assumption should be carefully examined and validated on real data to avoid oversimplification. Another option would be to explore the use of recursive models [START_REF] Gianola | Quantitative genetic models for describing simultaneous and recursive relationships between phenotypes[END_REF] i.e., analyzing one effect after the other iteratively with the aim of having a final convergence).

TOWARD SUSTAINABLE EFFICIENCY: SHORT-VERSUS LONG-TERM EFFICIENCY

So far, the details of the RFI equation have been discussed to contribute to making the RFI method more adapted to the modern context of precision livestock data. In doing this, the technical issues related to the time periods of measurement have been dealt with. Here, the focus is back to the important issue concerning time periods of measurement (i.e., the value and interpretation of short-vs. long-term efficiency measures). The underlying idea is that sustainable efficiency should be the goal of any efficiency improvement strategies. This means considering the longer-term consequences of any efficiency gains. To do this, 2 notions must be considered at the animal level: trade-offs and productive longevity (at the animal level; it gets more complicated at the herd or farm levels; [START_REF] Faverdin | Intégrer les changements d'échelle pour améliorer l'efficience des productions animales et réduire les rejets[END_REF].

As described above, the most efficient animals are those able to produce more at a given level of intake (or those eating less at a given level of production). It seems obvious that animals that put a greater proportion of their acquired resources into product will be more efficient, but this implies a reduced proportion of these resources will go to other life functions, which is known as a trade-off. A trade-off between productive and "fighting entropy" life functions would have possible negative consequences on animal robustness. In dairy cows, a critical function for staying in the herd is reproduction, and it is now well-documented that a negative (phenotypic and genetic) correlation exists between milk production and reproduction [START_REF] Boichard | Genetic analysis of conception rate in French Holstein dairy cattle[END_REF][START_REF] Pryce | Fertility in the high-producing dairy cow[END_REF], and thus intensive selection for milk production has degraded reproductive performance [START_REF] Miglior | A 100-Year Review: Identification and genetic selection of economically important traits in dairy cattle[END_REF]. There has also been a strong reduction in the average productive lifespan over the same period [START_REF] Hare | Survival rates and productive herd life of dairy cattle in the United States[END_REF]. More generally, the majority of evidence suggests that high levels of performance, be it milk, meat, or egg production, are associated with decreased robustness [START_REF] Rauw | Undesirable side effects of selection for high production efficiency in farm animals: A review[END_REF]. [START_REF] Loyau | Heritability of body surface temperature in hens estimated by infrared thermography at normal or hot temperatures and genetic correlations with egg and feather quality[END_REF] reported a negative genetic correlation between egg weight and surface body temperature under heat stress in laying hens (see also review of [START_REF] Moreri | Robustness to chronic heat stress in laying hens: A meta-analysis[END_REF]. In rabbit lines divergently selected for reproductive intensity or reproductive longevity, high reproductive intensity does were less resilient to an environmental perturbation and were less able to use their body reserves [START_REF] Theilgaard | Late reproductive senescence in a rabbit line hyper selected for reproductive longevity, and its association with body reserves[END_REF]. [START_REF] Gilbert | Review: Divergent selection for residual feed intake in the growing pig[END_REF] reported that a pig line selected for low RFI had significantly less body fat and a decreased fasting heat production relative to the high RFI line, although there were no clear differences in the growth between lines when exposed to a sanitary challenge [START_REF] Chatelet | Impact of hygiene of housing conditions on performance and health of two pig genetic lines divergent for residual feed intake[END_REF]. More recently, studies are emerging that seek to characterize resilience to environmental perturbations [START_REF] Ben Abdelkrim | Lactation curve model with explicit representation of perturbations as a phenotyping tool for dairy livestock precision farming[END_REF][START_REF] Adriaens | Productive life span and resilience rank can be predicted from on-farm first-parity sensor time series but not using a common equation across farms[END_REF] and to relate variation in these resilience measures to frequency of health events and productive longevity [START_REF] Berghof | Opportunities to improve resilience in animal breeding programs[END_REF]. [START_REF] Poppe | Exploration of variance, autocorrelation, and skewness of deviations from lactation curves as resilience indicators for breeding[END_REF], building on work of [START_REF] Elgersma | Fluctuations in milk yield are heritable and can be used as a resilience indicator to breed healthy cows[END_REF] and using the residual variance in milk yield from daily milk records as a resilience measure, found favorable genetic correlations between resilience and health, fertility, and longevity. Taken together, this literature suggests that short-term efficiency is unfavorably correlated with resilience, and that resilience contributes positively to a longer productive lifespan, at least in adult producers. From the dairy cattle data of the DEFFILAIT project, we have recently been able to look at the relation between RFI and resilience, measured as the accumulated loss in milk yield during perturbations (based on 237 lactations; unpublished). Because the losses in milk yield were matched by losses in intake during these short-term perturbations (Figure 3; [START_REF] Bareille | Effects of health disorders on feed intake and milk production in dairy cows[END_REF], whole lactation RFI was barely affected. This suggests that the link between resilience and efficiency is mainly through the effect of resilience on productive lifespan.

Returning to consideration of the time-span over which efficiency is measured and going to the extreme of considering lifetime efficiency, it is probable that the more efficient animals may be the resilient ones (i.e., those able to cope well with environmental perturbations and thereby minimize effects on both production and other life functions such as health and reproduction). These longer-lived animals will increase the proportion of their lifespan that is productive, and thereby dilute the efficiency cost of the nonproductive rearing period. Accordingly, as shown in Figure 4, we could hypothesize that the correlation between feed efficiency and resilience changes with the duration of the measurement period. Indeed, results of simulations using a model designed to represent resource allocation between life functions suggest that different resource allocation characteristics underpin short-term efficiency and long-term efficiency [START_REF] Puillet | Disentangling the relative roles of resource acquisition and allocation on animal feed efficiency: Insights from a dairy cow model[END_REF][START_REF] Puillet | Exploring underlying drivers of genotype by environment interactions in feed efficiency traits for dairy cattle with a mechanistic model involving energy acquisition and allocation[END_REF]. However, these results remain to be validated with animal data. Such data would also allow quantification of what measurement duration is needed to avoid a negative correlation between efficiency and resilience. This is important for knowing the period over which we need to measure efficiency (i.e., defining the appropriate phenotypes that would allow selection for sustainable efficiency).

A key question is then, what kind of feed efficiency do we want to select for? Do we want an animal that is highly efficient for a specific period of its life (a few weeks around the peak of lactation for example), or do we want it to have a better accumulated efficiency over its life? Selection for a short-term feed efficiency may select for a switch in the resource allocation that will benefit production but deteriorate resilience and the key life functions that contribute to longevity. On the other hand, if selection for a longer productive lifespan degrades production too much, then the lifetime efficiency of such animals may be not better than that of short-lived high producing animals. In this context, the idea of selecting for sustainable efficiency seems promising. These issues need to be explored not just for adult producers, but also in the context of growing animals reared for meat production. This opens up for consideration of efficiency not just at the animal level, but also at herd or farm level. A final point that should be considered when coming to selection decisions is that in this article the interest was the biological nature of RFI, and thus it was defined at the phenotypic scale. However, it is possible to model RFI to be genetically (and not phenotypically) independent from the energy sinks that are used as predictors [START_REF] Kennedy | Genetic and statistical properties of residual feed intake[END_REF], which ensures the absence of undesired genetic correlations between these traits in selection programs.

SUMMARY

Until very recently, the notion of measuring feed efficiency over a whole lactation, and even more so over a whole lifespan, seemed beyond reach. It is therefore not surprising that data sets with long-term measures of intake, and thus long-term efficiency, are very rare, especially if good resilience measures are also required. However, the burgeoning field of precision livestock technologies, with its promise of automated highfrequency measures, has now developed to the point where not only milk production but also body reserve (weight and condition score) changes [START_REF] Thorup | Robustness and sensitivity of a blueprint for on-farm estimation of dairy cow energy balance[END_REF][START_REF] Fischer | Precision feed restriction improves feed and milk efficiencies and reduces methane emissions of less efficient lactating Holstein cows without impairing their performance[END_REF], and more recently intake measures [START_REF] Lassen | Individual measure of feed intake on in-house commercial dairy cattle using 3D camera system[END_REF], are becoming available as a commercial reality that can be used on sufficient scale to allow phenotyping of large populations of animals for long-term efficiency.

The increasing availability of such time-series data brought into focus several limitations of the traditional RFI methodology, particularly in relation to the time periods of measure and the ability to evaluate the importance of the different biological components of overall feed efficiency. We believe that the technical advances proposed above will contribute to resolving these limitations, and thus to making the RFI method fit for purpose in the precision livestock future. This should allow the relationships between short-and long-term efficiency to be evaluated and to be reframed in terms of sustainable efficiency goals that can be included in future selection and management strategies.

Figure 1 .

 1 Figure 1. Classes of efficiency. "Maintenance" is considered to include all the nonproductive life functions that can contribute positively to longevity.

  Intake = (a + BW SIGN × b) × ΔBW + (c + BCS SIGN × d) × ΔBCS + other predictors + RFI, [8]

  ] with the consequence that the values of a′ and b′ [(that correspond to (a + BWSIGN × b) and (c + BWSIGN × d) at a given time, respectively] differ at different times (e.g., stages of lactation). All potential interactions would be included in the term "other predictors." The possibility of comparing coefficients across different stages of lactation has been explored further (DEFFILAIT project, unpublished; www .deffilait .fr).

Figure 2 .

 2 Figure 2. Correlation between the predicted intake and the animal effects of its predictors [fat-and protein-corrected milk (Cmilk), BW (Weight), and BCS] across the lactation. Figure adapted from Martin et al. (2020).

Figure 3 .

 3 Figure 3. Relationship between the cumulative decreases in DMI during feed intake perturbations and the cumulative losses of milk yield during milk perturbations (average of 250 first days of 237 lactations, DEFFILAIT project data, www .deffilait .fr).

Figure 4 .

 4 Figure 4. Schematic representation of the possible effect of residual feed intake measurement period on the correlation between efficiency and resilience.

Table 1 .

 1 Martin et al.: INVITED REVIEW: DISENTANGLING RESIDUAL FEED INTAKE Abbreviations used in this article

	Abbreviation	Meaning
	BCS	Body condition score
	BW	Body weight
	ed	Energy density (assumed digestive efficiency)
	ded	Deviations from the assumed digestive efficiency
	DMI	Dry matter intake
	eif	Errors in fitting
	eim	Errors in measurement
	FI	Feed intake
	imf	Individual variability in metabolic efficiency from the identified energy sinks
	ues	Unidentified energy sinks
	ΔBW	Changes in body weight (gain or loss)
	ΔBCS	Changes in body condition score (gain or loss)

Journal of Dairy Science Vol. 104 No. 6, 2021

ACKNOWLEDGMENTS

The research leading to these results has received funding from European Union's Horizon 2020 research and innovation programme-GenTORE-under grant agreement No. 727213, and from the national project DEFFILAIT (ANR-15-CE20-0014, France) and from APIS-GENE (Paris, France). The authors have not stated any conflicts of interest.


Riccardo Artoni

Michele Larcher

James T Jenkins

Patrick Richard

Self-Diffusion Scalings in Dense Granular Flows

We report on measurements of self-diffusion coefficients in discrete numerical simulations of steady, homogeneous, collisional shearing flows of nearly identical, frictional, inelastic spheres. We focus on a range of relatively high solid volume fractions that are important in those terrestrial gravitational shearing flows that are dominated by collisional interactions. Diffusion over this range of solid fraction has not been well characterized in previous studies. We first compare the measured values with an empirical scaling based on shear rate previously proposed in the literature, and highlight the presence of anisotropy and the solid fraction dependence. We then compare the numerical measurements with those predicted by the kinetic theory for shearing flows of inelastic spheres and offer an explanation for why the measured and predicted values differ.

Introduction

Collisions between spheres in a dense granular shearing flow induce velocity fluctuations of the grains that drive the diffusion of particles in a fashion that is analogous to the thermal diffusion in a dense gas of elastic molecules, or the diffusion induced by eddies in the turbulent flow of a fluid. Here, we report on measurements of the components of self-diffusion parallel and perpendicular to the flow done in discrete numerical simulations of inelastic spheres in a dense shearing flow. The components are determined by measuring the average squared displacement of spheres as a function of time.

Diffusion in granular shearing flows is important to mixing and segregation and has been studied experimentally in various flow geometries: granular shear cells 1,2 ; vertical channels [START_REF] Hunt | Advances in Micromechanics of Granular Materials[END_REF][START_REF] Natarajan | [END_REF] ; inclined chutes [5][6][7] ; vibrationally excited systems 8 ; free-surface flows 9 ; rotating tumblers [10][11][12][13] ; and rotating tubes 14,15 . Savage & Dai 16 and Thornton, et al. 17 have carried out studies of segregation in discrete numerical simulations; and phenomenological theories exist, such as those described by Gray & Ancey 18 and Fan & Hill 19 , that produce plausible predictions of species' concentrations and mixture velocity for appropriate choices of parameters. Taberlet and Richard 20 studied the spreading of a granular pulse in numerical simulations of bidisperse mixtures in a rotating drum. All of these studies are interpreted in terms of the mechanism of diffusion. The ability of kinetic theory to predict diffusion properties in 2D granular systems sustained by an air table has also been studied extensively in 21 for a wide range of surface fractions .

Direct experimental and numerical measurements of components of the tensor of self-diffusion exist in regimes of flow somewhat different from that considered here. Campbell 22 carried out such measurements in a sheared system of spheres that interacted through frictional, inelastic collisions over a range of restitution coefficients from 0.4 to 1.0 and solid fractions from 0.0001 to 0.5.

a MAST-GPEM, Univ Gustave Eiffel, IFSTTAR, F-44344 Bouguenais, France; E-mail: riccardo.artoni@univ-eiffel.fr b Free University of Bozen-Bolzano, I-39100 Bozen-Bolzano, Italy. c Cornell University, Ithaca, NY 14053, USA.

Macaulay & Rognon 23 investigated the effect of inter-granular cohesive forces on the properties of self-diffusion in dense granular flows. Utter & Behringer 24 performed measurements in experiments on slow, rate-independent shearing of a dense aggregate of disks in a Couette cell. In two-dimensional numerical simulations of a similar system in a periodic cell, Radjai & Roux 25 measured properties of the particle velocity fluctuations and the components of the self-diffusion. Because these systems involve rate-independent interactions, they differ from the system that we consider.

The measurements that we report are similar to Campbell's, in their common range of solid volume fraction; but different, in our focus on the range of volume fractions between 0.49 and 0.6. This is the range of solid fraction important in geophysical flows on Earth. Our interest is in the values of the components of the tensor of self-diffusion, particularly over the range of solid fraction between 0.49, above which long range order might appear in equilibrated system of monosized elastic spheres 26 , and a solid fraction of about φ c = 0.587, at which a collisional flow becomes impossible, when the coefficient of sliding friction is 0.5 27 .

Previous studies have pointed out that self diffusion coefficient scales simply, in dense systems, as D = kd 2 γ, where γ is the shear rate, d the particle diameter, and k approximately a constant of order 0.05 24,28,29 . This scaling is, however, empirical. In addition, in such dense systems, little attention has been given to the tensorial nature of self-diffusivity and to its dependence on solid fraction. Therefore, our first goal is to study the full self-diffusivity tensor over the range of dense solid fractions from 0.49 to 0.587. In addition, to our knowledge, the scaling based on shear rate has not been compared to micromechanical theories, and its general validity may, therefore, be questioned. In complex flows (e.g. flows characterized by shear localization, creep zones, and those influenced by boundaries 30 ) the rheology is known to become nonlocal, and the introduction of velocity fluctuations as an additional variable seems a promising path.

The strength of velocity fluctuations is a classical ingredient of diffusion theories. Dense kinetic theories for the segregation of binary mixtures of inelastic spheres [31][32][33][34][35] predict diffusion coefficients that exhibit explicit dependence on solid fraction, the strength of the particle velocity fluctuations and the particles' size, mass, and collision properties. Consequently, our second goal is to characterize the scaling of self diffusion with respect to granular temperature and compare the measured values of the selfdiffusion coefficient to those predicted by dense kinetic theory.

Discrete numerical simulations

Simulations were performed by means of the open-source molecular dynamics software LAMMPS 36 . A cubic simulation cell was used, with a size of 20x20x20 in particle diameters, in which the solid fraction φ was varied from 0.1 to 0.61. The results presented here are, however, restricted to the dense flow range of φ = 0.49 to 0.586 (number of particles N = 7400 to 8838). The choice of the system size was motivated by the need of having a system larger than the typical correlation length (few particle diameters), but small enough to limit the computational cost. Previous literature 37 suggested that L = 20d may be a good compromise. On our hand, we verified this by running dedicated simulations for a larger system size (L = 40d, not shown here), which gave similar results in terms of self-diffusivity, and therefore support our choice for L. In order to avoid crystallization, a slight polydispersity was introduced: the dimensionless diameter, d, ranged uniformly from 0.9 to 1.1. The mass of the particle with unit diameter was taken as the mass scale m. Then, in dimensionless terms, the dimensionless mass density of the spheres is ρ = 6/π.

For the normal component of contact, a linear spring dashpot model, as employed by Silbert et al. 38 , was used. In this, the force between particles i and j is given by F

i j n = k n δ i j n -m i j e f f γ n δ i j n ,
where k n is the normal stiffness, δ n the particle interpenetration, m i j e f f = m i m j /(m i + m j ) the effective mass of the interaction, and γ n the specific damping coefficient. For the tangential component, an elastic model with stiffness k t , no viscous damping and a frictional threshold was employed. The normal spring stiffness provides an intrinsic time scale for the system; m/k n was used to nondimensionalize the time. The tangential stiffness was set as k t /k n = 2/7, for the periods of normal and tangential oscillations to be equal 39 . The normal damping coefficient, expressed in normalized units through the normalized time, was varied in the range γ n = (0.095 -0.609) k n /m, corresponding to a range of restitution coefficients e n = 0.5 -0.9 for the collision of two unit diameter particles. Note that, due to polydispersity, a slight heterogeneity of the effective restitution coefficient is expected. The coefficient of sliding friction was chosen to be µ = 0.5. In molecular dynamics discrete element simulations, the computational time step ∆t is usually set as a small fraction of the collision time, in order to ensure proper simulation of contact dynamics. The collision time for the linear spring-dashpot model is given by the relation:

t coll = π k n /m e f f -(γ n /2) 2 -1/2
. In our case, the term containing the damping coefficient is quite negligible and therefore t coll ≈ π/ √ 2 m/k n . The time step was therefore taken as ∆t = 2 × 10 -2 m/k n , which corresponds, approximately, to onehundredth of a collision time.

The numerical simulations were performed in simple shear under fully periodic boundary conditions by means of the LAMMPS 'fix deform' scheme which, similar to the method used by Radjai and Roux 25 , applies shear by deformation of the simulation box. Note that this is different from Campbell 22 who used a "sliding blocks", Lees-Edwards 40 scheme for ensuring shear under periodic boundary conditions. Coordinates x, y, and z correspond in the directions of flow, vorticity, and gradient, respectively. Two types of simulations were carried out over the specified range of solid fraction: one with the shear rate constant over the range γ = (10 -6 -10 -3 ) k n /m; the other with the pressure held approximately constant, obtained by decreasing the shear rate while increasing the solid fraction. This was made to check that the results depended on shear rate only trivially, as is expected in the range of shear rates considered. Given that this was verified, we focus on results obtained for a shear rate of γ = 10 -3 k n /m.

The simulation was carried out in two steps: first, the initial state was generated on a lattice partially filled with the desired distribution of particles, and sheared for 10 9 steps to ensure a steady state. A Q6-analysis 41,42 was used to check that the shear applied during 10 9 time steps was enough to remove any trace of the orientational order of the initial state. In the second phase, data analysis was performed, while the system was sheared at the same rate as in the preparation phase, again for 10 9 steps. In this phase, snapshots of positions, velocities and inter-particle forces were extracted every 10 4 steps. The cumulative deformation correponding to each of the two phases, preparation and analysis, was therefore, γ∆t = 200. In contrast, the simulations of Radjai & Roux 25 in two-dimensions had an applied total strain of about two, at an area fraction around 0.8. The relatively large value of cumulative deformation employed here was necessary in order to determine the scaling of self-diffusion with time. The simulation output was treated in two steps: first, particle trajectories were reconstructed from time snapshots, removing the instantaneous mean field; then, the statistics of the velocity probability distribution function and self-diffusion were calculated.

Results

Velocity distribution functions

In Fig. 1a we show the probability distribution function (PDF) of, for example, the instantaneous z-velocity fluctuations for different values of the average solid fraction for e n = 0.7. It is evident that in dense systems the velocity distribution displays an exponential tail, as in the case of rate-independent shearing 25 , which broadens with increasing solid fraction. As the inset of Fig. 1a clarifies, such broadening of the exponential tail is important for φ > 0.55, and quite independent of the restitution coefficient. Next, we describe the probability distribution function of each component. Note that for computing the velocity fluctuation in the x-direction, we subtract the local mean velocity related to the mean shear (δ v x = v xγ(zz 0 ), where z 0 is the center of the cell) 43 . In order to characterize the effect of φ , e n and direction on the velocity PDF, we show two statistical descriptors for each component: the variance δ v departure from a normal distribution, for which it is zero (it is equal to three for a Laplace distribution). In Fig. 1b, the variance of the different components of the velocity fluctuations behaves in a similar way, but the intensity of the fluctuations is different: for relatively low solid fractions, we observe

δ v 2 y < δ v 2 z < δ v 2 x ; while, approaching φ = 0.587, the order is δ v 2 y < δ v 2 x < δ v 2 z .
Nevertheless, the anisotropy of velocity fluctuations appears to decrease with φ . With respect to the restitution coefficient, for strong dissipation (points corresponding to e n ≤ 0.7 in the figure), δ v 2 i /( γd) 2 is monotonically increasing with φ , for the dense flows considered here. For larger values of the restitution coefficient (as in the case e n = 0.9 presented in the figure), the variance first decreases and then increases with φ . Note that for dilute situations, T /( γd) 2 is reported to decrease with φ 22 . So, the value of the solid fraction for which the granular temperature displays a minimum seems to depend on e n . Moreover, for low values of the solid fraction, the variances increase when increasing e n ; while, when approaching φ = 0.587, they appear to become independent of the restitution coefficient.

The excess kurtosis, shown in Fig. 1c, generally increases with solid fraction, which mirrors the broadening of exponential tails with φ observed in Fig. 1a. Deviation from a normal distribution appears to be stronger in the flow direction x. The effect of the restitution coefficient is also stronger for the flow direction: decreasing the coefficient of restitution yields an increase of the excess kurtosis, which seems to be limited to low solid fractions. As was observed for the variance, the anisotropy between the velocity distributions is reduced when increasing φ . Based on the data collected in Fig. 1, we can conclude that the fluctuation velocity vector distribution is non-Maxwellian in both the anisotropy and the exponential tail, and that the importance of these two effects depends on the values of the restitution coefficient and the solid fraction.

Self-diffusion

The components of the diffusion tensor were determined by tracking the movement of the particles relative to their initial position, while taking into account the displacement due to the mean shear flow. We find that the particle self-diffusion, corresponding to correlation of displacements in directions i and j, is proportional to a power of the time:

∆x i ∆x j ∝ ∆t α , (1)
where the exponent is not constant. As Fig. 2 exemplifies for the transverse yy-component, for small cumulative deformations ( γ∆t < 1), an exponent α ≈ 1.8 is found, which corresponds to super-diffusive motions; while for large cumulative deformations ( γ∆t > 1), a simple diffusive behavior is evident, with α ≈ 1. This double scaling is in agreement with previous results for dilute systems by Campbell 22 , and also with those in dilute collisional suspensions 44 . It seems to be the simple consequence, well known in turbulence [45][46][47][48] , of the apparent diffusion associated with a random process with a finite correlation time (which is superdiffusive for short times). For dense granular flows, particle displacements are constrained and frustrated by the mutual hindrance between neighboring particles, and therefore the physical mechanism behind their time evolution is different from turbulence. The present results show, however, that the diffusive behavior is present even in dense systems, and seem to indicate that the super-diffusive behavior observed by Radjai & Roux 25 may be due the small cumulative deformation they employed.

Given that the behavior is diffusive for large cumulative deformations, it is possible to define a self-diffusion tensor, as the limit for large cumulative deformations of the cumulative displacements correlations:

D i j = lim ∆t→∞ ∆x i ∆x j 2∆t . (2)
For dilute granular shear flows (φ < 0.5), Campbell 22 analyzed the components of the self-diffusivity tensor, scaled by γd 2 , and showed that the tensor was anisotropic with a clear hierarchy (D xx > D zz > D yy > D xz ), and that the off-diagonal components other than xz were negligible. In the following, we do not discuss the xy and yz components of the diffusion tensor, because, as in Campbell 22 , these components are negligible. The four other non-zero components are shown in Figs. 3 and4. We first consider them normalized by the shear rate and the squared particle diameter, inspired by the empirical scaling previously proposed 24,28,29 . First, the off diagonal term seems to go to zero in the limit of φ → φ c . Moreover, its magnitude is well below that of the diagonal terms. This appears to be a peculiarity of dense flows, D xz being comparable to D yy and D zz for dilute flows. On the other hand, regarding the dependence on solid fraction, it is evident that the diagonal components, scaled by the shear rate, display a nonmonotonic behavior, the strongest example being given by the streamwise, xx-component, which first and then increases with φ . Then, the diagonal terms display a moderate but evident anisotropy. Similarly to what is obseved in dilute flows 22 , we obtain D xx > D zz > D yy > D xz . Yet the magnitude of the latter components are closer to each other in the dense case.

The anisotropy of the self-diffusion tensor is quantified by the indicator

(D 1 -D 2 ) 2 + (D 2 -D 3 ) 2 + (D 3 -D 1 ) 2 2(D 2 1 + D 3 2 + D 2 3 )
, where the D i are the eigenvalues of the diffusion tensor, which is displayed in the inset of Fig. 3. It is evident that anisotropy decreases with increasing solid fraction, but does not disappear approaching φ c . Finally, the diffusivities, scaled by the shear rate, are independent of the restitution coefficient. Based on these results, we can conclude that the empirical scaling D ≈ 0.05 γd 2 24,28,29 gives the correct order of magnitude for the trace of the diffusivity tensor in the range of dense solid fractions consid-ered. Therefore, we think that the empirical scaling cited above may be employed in approximate analyses. Refined analyses must, however, take into account the nonmonotonic behavior and anisotropy of the D i j .

0

Anisotropy

Fig. 3 The four non-zero components (symbols) of the self-diffusion tensor plotted together for different restitution coefficients (colors, same as in figure 1) and different values of the solid fraction. The tensor components are scaled by the product of the shear rate and squared particle diameter. The inset represents the anisotropy of the tensor as a function of the solid fraction via the measure detailed in the text.

Questions may be raised about the micromechanical origin of the scaling on diffusion on shear rate. A possible micromechanical framework to interpret this empirical result may be found in kinetic theory. As in isotropic turbulence [START_REF] Pope | Turbulent Flows[END_REF] , expressions for the selfdiffusivity in an isotropic dense granular gas can be represented by the formula D = T τ v where τ v is the time of autocorrelation of velocity fluctuations. In the kinetic theory of dense gases of elastic spheres Chapman & Cowling [START_REF] Chapman | The mathematical theory of non-uniform gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases[END_REF] , τ v = d π/T /(16φ g 0 (φ )), where g 0 (φ ) is the radial distribution function at contact. From this the classical scaling is obtained, D = √ πT d/(16φ g 0 (φ )), which was, for example, used by Larcher and Jenkins 34 , 35 for a dense gas of frictional, slightly inelastic spheres. The relevant parameter in such an expression that sets the time scale of diffusion is the strength of velocity fluctuations. Inelasticity in granular gases is known to increase the spatial and temporal span of velocity correlations, therefore increasing self-diffusion. In the isotropic case, this results in a correction to Chapman & Cowling's formula, and the scaling for the self-diffusivity of the kinetic theory of granular gases [START_REF] Brilliantov | Kinetic theory of granular gases[END_REF] is therefore:

D KT = √ π 8(1 + e) d φ g 0 (φ ) √ T . (3)
In order to evaluate this framework and the eventual modifi- Fig. 4 The four non-zero components (symbols, same as in Fig. 3) of the self-diffusion tensor plotted together for different restitution coefficients (colors, same as in Fig. 1) and different values of the solid fraction.

The tensor components are scaled by the square root of the granular temperature. The inset represents the ratio of the numerical D yy to the kinetic theory scaling, as a function of the solid fraction.

cations needed to adapt it to the dense systems considered here, we plot again in Fig. 4 the self-diffusion tensor components, this time scaled by the square root of the granular temperature and by the particle diameter. As in kinetic theory, the ratios D i j /d √ T are decreasing functions of φ . Then, while the off diagonal component goes to zero as φ → φ c , the diagonal components remain finite. Moreover, for low values of the dense solid fraction, the D i j /d √ T depend on the restitution coefficient, which is not the case for high values of the dense solid fraction.

Given that the self-diffusivity tensor displays some anisotropy, we compare isotropic scalings from kinetic theory to the transverse diffusivity component, D yy , which is less affected by the anisotropy of shear. The radial distribution function at contact, g 0 (φ ), is often operationally defined through the equation of state for the pressure. The classical result by Torquato [START_REF] Torquato | [END_REF] for hard elastic frictionless spheres, for solid fractions between freezing and random close packing (φ F < φ < φ RCP , where φ F = 0.49 and

φ RCP ≈ 0.64), is g 0 (φ ) = 5.6916 φ RCP -φ F φ RCP -φ . ( 4
)
We note that in the framework of extended kinetic theory 53,54 , Berzi and Vescovi 27 have discussed constitutive relations for frictional inelastic particles based on a contact radial distribution function that possesses a singularity at a critical solid fraction φ c lower than the random close packing that depends on the friction coefficient. However, in our view, such singular behavior is appropriate for collisional transfers of momentum and energy, but not for those, such as diffusion, that involve transport of mass.

Clearly, our numerical results support the choice of a radial distribution function not divergent at φ c , because the diagonal components remain finite when approaching the critical solid fraction. Therefore, in the following, for estimating the theoretical self-diffusivity, D KT , we combine Eq. 3 with the classical radial distribution function by Torquato [START_REF] Torquato | [END_REF] , Eq. 4.

In the inset of Fig. 4, we plot the ratio between the numerically obtained D yy and the prediction from kinetic theory, D KT , by employing the classical radial distribution function by Torquato [START_REF] Torquato | [END_REF] . It is evident that, for low solid fractions, the kinetic theory prediction is not far from the measured values, and correctly models the effect of the restitution coefficient. However, for denser systems, the deviation increases and the ratio D yy /D KT reaches a value of about three. Note that the deviations from the kinetic theory observed for dense systems have also been reported experimentally and numerically in 2D systems of disks submitted to random fluctuations induced by an air table 21 . In the following subsection, we discuss the possible origin for the differences between kinetic theory and the simulation.

Correlated bulk motion

In order to determine whether part of the deviation from kinetic theory could come from correlated bulk motions, we computed energy spectra of spatial velocity fluctuations. For each time snapshot of the system, we first interpolated the particle velocity fluctuations on a regular grid, then obtained two-point velocity correlations and the energy spectrum through the (spatial) Fourier transform of the interpolated field. As in turbulence 37,[START_REF] Pope | Turbulent Flows[END_REF] , the energy spectrum was spherically averaged with respect to the wavenumber k.

In Fig. 5 we plot the energy spectra, normalized by the value at the smallest wavenumber, k = 2π/L. It is evident that the spectral energy density is not a monotonic function of the wavenumber: particularly for lower dense solid fractions, the spectra display a maximum and then decrease with, ultimately, a power-law cutoff. The energy-containing scale represented by the position of the maximum of the spectrum, slightly depends on solid fraction, as does the integral length scale extracted from the two-point correlation functions, varying between one and two particle diameters. It is evident that velocity correlations associated with several particles exist and contribute energy to the spectra.

In granular shearing flows, there is not the same separation of scales as in molecular gases. However, in dense granular shearing flows, the particles interact over a length scale that is the distance between their edges, which is a fraction of their diameter. At this scale, pre-collisional velocity correlations exist as shown by Mitarai and Nakanishi 55 , which affect the amount of velocity fluctuations. With this in mind, we suggest that a more appropriate measure of the temperature may be the energy of the velocity fluctuations at a scale below the particle diameter, in relation to pre-collisional velocity correlations. Velocity fluctuations evaluated at lengthscales of the order of tens of particle diameters, then, may be analogous to the macroscopic fluid velocity fluctuations that contribute to a turbulent eddy diffusivity. Al-Fig. 5 Power spectra for e n = 0.7 and several solid fractions. The inset represents the integral length scale obtained from two-point correlations, as a function of the solid fraction, for different restitution coefficients (colors, same as in Fig. 1).

though in the present soft particle simulation it was not possible to implement a detailed follow-up of collisions, and therefore it is impossible to access length scales smaller than the diameter of the spheres, we note that the energy of the velocity fluctuations at this scale is an order of magnitude less than that at the largest length scales. Hard sphere simulations of the type carried out by Mitarai and Nakanishi 55 , in a study of pre-collisional velocity correlations, seem to permit access to the velocity fluctuations at smaller length scales.

Conclusions

We have characterized the properties of the velocity distribution function and the components of the self-diffusion tensor in discrete numerical simulations of dense, steady, homogeneous shearing flows of nearly identical, inelastic, frictional spheres. We have discussed the anisotropy of the self-diffusion tensor, and its dependence on solid fraction. Our results provided a test for the empirical scaling D ∼ 0.05d 2 γ in a wide range of dense solid fractions. We found that, although such a scaling gives the correct order for the diagonal components of the self-diffusivity, a nonmonotonic dependence on φ as well as a moderate anisotropy are present, which may be important in refined analyses.

In order to look for micromechanical explanations for the scaling cited above, we compared the values for one of the components of the self-diffusion to that predicted by kinetic theory. When the strength of the velocity fluctuations at all length scales was employed as the granular temperature, kinetic theory was found to under-predict self-diffusion. A spectral analysis of the velocity fluctuations indicated how the strength of the velocity fluctuations varied with their wavelength. This variation led us to suggest that there may be an effective separation of scales between fluctuations at a fraction of diameter, at which the particles interact, and those larger than a particle diameter that may be the analog of the eddies of a turbulent fluid 56,57 . This effective separation of scales might be at the origin of the deviation from kinetic theory.

Our aim is to obtain relations applicable to the evolution in space and time of granular segregation in industrial processes and geophysical flows, in which granular diffusion takes place at high solid fractions. It is well known that in such heterogeneous flows, granular temperature is an important dynamic variable, which can be used to model nonlocal effects. Here we have shown that it is possible to model the magnitude of the self-diffusion tensor by an empirical law based on shear rate 24,28,29 , but also by kinetic theory, provided that a correction for dense systems is introduced. This correction can be written as:

D * KT = ξ (φ ) √ π 8(1 + e) d φ g 0 (φ ) √ T , (5)
where ξ is a correlation factor corresponding to the inset of Fig. 4. It is important to consider such a scaling based on kinetic theory partly because of its micromechanical origin, and partly because at present the effect of nonlocality on self-diffusion is not clear. Further research will deal with the measurement of self-diffusivities in heterogeneous flows in order to provide more evidences concerning the correcting factor ξ and determine the relative validity of the two frameworks discussed above.

2 i 2 - 3 .

 223 2 i and the excess kurtosis δ v 4 i / δ v The average of the variances in the three directions is the usual definition of the granular temperature, T ≡ ∑ i δ v 2 i /3. The excess kurtosis is a measure of the flatness of the distribution, and can be used to characterize the 2 | 1-7 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Fig. 1

 1 Fig. 1 (a) PDF of the instantaneous particle velocities for different values of the average solid fraction for e n = 0.7 (inset: slope of the logarithm of the high energy tail of the PDF as a function of the average solid fraction for three values of the normal restitution coefficient, colors for e n detailed in (b)), (b) Dimensionless variance, and (c) excess kurtosis of the velocity distribution function for the three components of the velocity (symbols, detailed in (c)), and three values of the normal restitution coefficient (colors, detailed in (b)).

2 < ∆y 2 >Fig. 2

 222 Fig.2Average, cumulative squared displacement in the y direction as a function of cumulative deformation, for e n = 0.7, φ = The existence of two different power laws is evident, one with exponent α ≈ 1.8 for low deformations, and one with exponent α ≈ 1 for large deformations. For large deformation the behavior is simple diffusive.

1-7 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

The numerical simulations were carried out at the CCIPL (Centre de Calcul Intensif des Pays de la Loire) under the project "Simulation numérique discrète de la fracture des matériaux granulaires".

Notes and references

Conflicts of interest

There are no conflicts to declare.