C Denis

E Lebarbier
email: emilie.lebarbier@agroparistech.fr

C L Évy-Leduc

AND O Martin

L Sansonnet
email: laure.sansonnet@agroparistech.fr

A NOVEL REGULARIZED APPROACH FOR FUNCTIONAL DATA CLUSTERING: AN APPLICATION TO MILKING KINETICS IN DAIRY GOATS

Motivated by an application to the clustering of milking kinetics of dairy goats, we propose in this paper a novel approach for functional data clustering. This issue is of growing interest in precision livestock farming that has been largely based on the development of data acquisition automation and on the development of interpretative tools to capitalize on high-throughput raw data and to generate benchmarks for phenotypic traits. The method that we propose in this paper falls in this context. Our methodology relies on a piecewise linear estimation of curves based on a novel regularized changepoint estimation method and on the k-means algorithm applied to a vector of coefficients summarizing the curves. The statistical performance of our method is assessed through numerical experiments and is thoroughly compared with existing ones. Our technique is finally applied to milk emission kinetics data with the aim of a better characterization of inter-animal variability and toward a better understanding of the lactation process.

Introduction

Precision livestock farming is a blooming field grounded in the development of sensors providing high throughput data and thus potentially increasing access to valuable information on biological processes. Therefore, developing methods for data analysis and interpretation has become a challenging issue in animal science. Economic performance of dairy goat farming systems is primarily based on milk production and a large amount of farmers working time is spent milking animals, see [START_REF] Marnet | Machine milking ability in goats: Genetic variability and physiological basis of milk flow rate[END_REF]. Moreover, with the increasing size of goat herds and the rapid growth of the dairy goat industry, more in-depth information on individual milking performance is necessary. In this context, a better understanding of the variability in milk flow kinetics could for instance help refining selection criteria for breeding programs, simplifying milking workload or controlling udder health. Milk emission kinetics recorded during milking of dairy goats are classically described and classified through synthetic parameters such as milking time, maximum and average milk flow rates, time to reach 500 g/min milk flow, see [START_REF] Romero | Relationship of goat milk flow emission variables with milking routine, milking parameters, milking machine characteristics and goat physiology[END_REF]. In this paper, we explore the possibility of considering milk emission kinetics as a whole function, opening new perspectives to study inter-animal variability.

From a statistical point of view, this issue belongs to the general field of functional data analysis, see [START_REF] Ramsay | Functional data analysis[END_REF] for a survey on this subject. In the specific functional data clustering framework, several approaches have been proposed by [START_REF] Abraham | Unsupervised curve clustering using b-splines[END_REF], [START_REF] Jacques | Funclust: a curves clustering method using functional random variables density approximation[END_REF] and [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] among others. For a review on this subject, we refer the reader to Jacques and Preda (2014a) and the references therein. This kind of approaches was extended to deal with multivariate functional data by Jacques and Preda (2014b) who proposed the first model-based clustering algorithm in this multivariate context and more recently by [START_REF] Schmutz | Clustering multivariate functional data in group-specific functional subspaces[END_REF].

To deal with the functional clustering of the milking kinetics of goats, some specific features have to be taken into account, see Figure 1 for some examples of such kinetics. We can see from this figure that these curves are nondecreasing and can be split into two parts, namely an increasing linear part and an almost constant one. Inspired by [START_REF] Abraham | Unsupervised curve clustering using b-splines[END_REF], we propose in this paper a dimension reduction approach based on a continuous piecewise linear function fit to each curve which boils down to a change-point detection issue which will be crucial in our method.

The problem of detecting change-points in the mean of a signal is largely addressed in the literature. In particular, it is now well known that in (penalized-) maximum likelihood frameworks the Dynamic Programming (DP) algorithm [START_REF] Bellman | On the approximation of curves by line segments using dynamic programming[END_REF]; [START_REF] Auger | Algorithms for the optimal identification of segments neighborhoods[END_REF]) and its recent pruned versions [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF]; [START_REF] Rigaill | A pruned dynamic programming algorithm to recover the best segmentations in 1 to Kmax changes[END_REF]; [START_REF] Maidstone | On optimal multiple changepoint algorithms for large data[END_REF] are the only algorithms that retrieve the exact solution very quickly. However, DP can only be used if the contrast to be optimized is additive with respect to the segments, see for example [START_REF] Bai | Computation and analysis of multiple structural change models[END_REF]; [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF]; [START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF]. When detecting changes in the slope with a continuity condition, the segments will unavoidably be linked and therefore the additivity condition is not satisfied. This q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 200 300 0 500 1000

2000

Goat 250053006212047

Time

Volume of milk q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 200 250 300 0 500 1000

2000

Goat 250053006212054

Time

Volume of milk q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 100 200 300 400 0 500 1000 1500

Goat 250053006212064

Time Volume of milk q q q q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 200 250 0 500 1500 2500

Goat 250053006212075

Time Volume of milk q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 200 250 300 0 500 1000 1500 2000

Goat 250053006212084

Time

Volume of milk q q q q q q q q q q q q q q q q q 0 50 100 150 [START_REF] Fearnhead | Detecting changes in slope with an l0 penalty[END_REF] proposed to extend the PELT algorithm [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF] to this problem. Their idea is to include the penalty in the DP algorithm with a pruning strategy. The penalty they proposed is proportional to the number of change-points up to a penalty constant. However, this penalty constant needs to be chosen in advance, which is not easy in practical situations.

In this paper, we first propose a novel change-point estimation in the slope method combining the trend filtering proposed by [START_REF] Tibshirani | Adaptive piecewise polynomial estimation via trend filtering[END_REF] with a (penalized-) maximum likelihood approach which is useful for removing the spurious change-points that may have been proposed by trend filtering. These change-points estimators are then used for devising a new dimension reduction approach: Each curve is summarized by a vector containing the coefficients of its projection onto an order 2 B-spline basis having for knots the obtained change-points and also the change-point locations. Including the change-points both in the features characterizing the curves and in the B-spline knots is the main novelty compared to classical approaches reviewed in Jacques and Preda (2014a).

The paper is organized as follows. The methodology that we propose is decribed in Section 2. The performance of our approach is investigated in Section 3 through numerical experiments. Finally, in Section 4, we apply our method to the data that motivated this study.

Methodology

In this section, we describe our novel functional data clustering approach which consists of two steps which can be summarized as follows:

• First step: Piecewise linear estimation of the curves using a novel change-point estimation method based on the trend filtering approach and B-splines. • Second step: Applying the k-means algorithm to a vector of coefficients summarizing the curves obtained in the first step.

These two steps are further described hereafter.

2.1. First step: Piecewise linear estimation of the curves based on a changepoint estimation method. In the following, we assume that the observations of a given curve Y = (Y 1 , . . . , Y n ) correspond to a noisy function evaluated at the input points x = (x 1 , . . . , x n ). In this step, we aim at estimating each curve by a piecewise linear function using a two-stage approach described below.

2.1.1. First stage: Trend filtering for change-point estimation. We use the trend filtering approach proposed by [START_REF] Tibshirani | Adaptive piecewise polynomial estimation via trend filtering[END_REF] which consists in fitting to the observations Y the vector β = ( β 1 , . . . , β n ) using a regularized method. More precisely, we use

β(λ) = Argmin β∈R n Y -β 2 2 + λ D (2) β 1 ,
where y 2 2 = n i=1 y 2 i , y 1 = n i=1 |y i |, for y = (y 1 , . . . , y n ), λ is a positive constant which has to be tuned and D (2) is the discrete difference operator of order 2 defined by

D (2) =   1 -2 1 0 0 • • • 0 0 1 -2 1 0 • • • 0 . . . . . . . . . . . . . . .   .
The final estimator of β is β( λ) where λ has to be properly chosen. Usually, this parameter is chosen using resampling approaches such as cross-validation or stability selection, see [START_REF] Meinshausen | Stability selection[END_REF]. From β( λ), we define a set of potential changepoint indices as the coordinates where the vector D (2) β( λ) is not equal to zero. However, in change-point estimation frameworks, the performance of such methods may be altered since some change-points may be omitted by subsampling. Moreover, it is well known that such regularization approaches lead to over-segmentation phenomena. Usually, in this case, a DP algorithm is then used on the set of potential change-points obtained with the latter strategy in order to remove the irrelevant ones, see for instance [START_REF] Harchaoui | Catching change-points with lasso[END_REF] and [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF].

We propose following this strategy: In order to avoid the use of a resampling method, we choose a small enough λ in order to obtain a large enough set of potential change-points. More precisely, we set a maximal number of change-points denoted K max and choose λ such that among the λ's leading to K max change-points, λ is the one minimizing Yβ(λ) 2 2 . Let ( n 1 , . . . , n Kmax ) the resulting change-point indices and the associated change-point positions ( t 1 , . . . , t Kmax ) = (x n 1 , . . . , x n Kmax ). For each K in {1, . . . , K max }, we use the DP algorithm to retrieve the K most relevant change-point indices among n 1 , . . . , n Kmax . DP is thus applied to Y n 1 , . . . , Y n Kmax instead of Y 1 , . . . , Y n . Note that a slight modification of the algorithm is considered to make the piecewise linear fit to data continuous. The optimal number of change-points K is then chosen by using the criterion proposed by [START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF].

2.1.2. Second stage: Projection onto the B-spline basis having as knots the obtained changepoints. Each curve will then be summarized by a few coefficients corresponding to the coefficients of its projection onto the B-spline basis (B i,2 ) 1≤i≤ K+2 defined as follows, see (Hastie et al., 2009, p. 206) for a review on the subject. Let t 0 = x 1 and t K+1 = x n . Let us also define the augmented knot sequence τ such that:

τ 1 = τ 2 = t 0 = x 1 , τ j+2 = t j , j = 1, . . . , K, τ K+3 = τ K+4 = t K+1 = x n , namely, (τ 1 , . . . , τ K+4 ) = (x 1 , x 1 , t 1 , • • • , t K , x n , x n ).
The ith B-spline function B i,2 having τ for knot sequence satisfies:

B i,2 (u) = u -τ i τ i+1 -τ i B i,1 (u) + τ i+2 -u τ i+2 -τ i+1 B i+1,1 (u),
where

B i,1 (u) = 1, if τ i ≤ u < τ i+1 0, otherwise with i ∈ {1, . . . , K + 2}.
Thus, each curve is estimated by f defined by:

(1)

f (u) = K+2 i=1 θ i B i,2 (u),
where the θi 's are obtained using a least-square criterion. Hence, the coefficients summarizing each curve is:

(

2) ( θ 1 , . . . , θ K+2 , t 1 , . . . , t K ).
2.2. Second step: Clustering using the k-means algorithm. In order to obtain a clustering of the curves (milking kinetics), we use the k-means algorithm of [START_REF] Hartigan | Algorithm AS 136: A K-means clustering algorithm[END_REF] on the scaled summarized coefficients (2) obtained in the previous step. It has to be noticed that the number of change-points K may change from one curve to the other. Thus, we consider summarized coefficients of length K M corresponding to the largest value of K. For kinetics having a number of change-points smaller than K M , we replace the missing t k and the missing coefficients by 0. Our goal is indeed to propose a strategy which is able to distinguish the curves both thanks to the change-point positions and/or the coefficient values.

The number k of clusters is chosen by using the strategy proposed by [START_REF] Charrad | NbClust: An R Package for Determining the Relevant Number of Clusters in a Data Set[END_REF] which consists in using the majority rule that is taking for k the value chosen by the largest number of criteria among 30 indices such as: CH index, Duda index, Pseudot2 index, C index, Hartigan index, ... Further details on these indices can be found in [START_REF] Charrad | NbClust: An R Package for Determining the Relevant Number of Clusters in a Data Set[END_REF]. Here, we focused on the four following indices: KL index, Hartigan index, SDindex, Ptbiserial index.

Numerical experiments

In this section, we investigate the statistical performance of our procedure. The simulation scheme that we used for this investigation is described in Section 3.1. We also propose in Section 3.2 to benchmark our procedure with existing approaches and to assess our change-point estimation approach in Section 3.3.

3.1. Simulation scheme. In order to be as close as possible to the data coming from our motivating application, we consider two different models for generating the data that we will refer to as Model 1 and Model 2 in the following. For each model, the complete observed data is (Y, Z), where Y is in R n and corresponds to the observations of an underlying function, which we will specify hereafter, at the input points x = (x i ) 1≤i≤n = (10(i -1)) 1≤i≤n with n = 51. Z denotes the label of Y which takes its value in Z = {1, 2, 3, 4}. Moreover, for each z ∈ Z, the associated cluster C z is characterized by a number of change-points K z , a vector of change-points t z , and a vector of parameters θ z ∈ R Kz+1 . Hence, each model is defined by a set of parameters {K z , t z , θ z : z ∈ Z}. The values of the parameters associated to each model are reported in Tables 1 and2. Note that for each model, the clusters are distinguishable by both the change points and the parameters.

For each model, the vector (Y, Z) is simulated according to the following procedure:

(a) The label Z is drawn from a uniform distribution on Z;

(b) We generate tZ = t Z + U, such that U = (U, . . . , U ), where U is a uniformly distributed random variable on {-30, -20, 10, 0, 10, 20, 30}; (c) We generate θZ = θ Z + V, such that V = (V, . . . , V ), where V is a uniformly distributed random variable on [-200, 200]; (d) Then, we consider the sequences ( tZ 0 , . . . , tZ K Z +1 ) = (0, tZ , 500), ( θZ 0 , . . . , θZ K Z +1 ) = (0, θZ ), and define for x ∈ [ tZ j , tZ j+1 ], and j ∈ {0, . . . , K Z }

(3)

f tZ , θZ (x) = ( θZ j+1 -θZ j )
x -tZ

j tZ j+1 -tZ j + θZ j ;
Y i = f tZ , θZ (x i ) + ε i ,
where the ε i 's are i.i.d N (0, σ 2 ) random variables with σ ∈ {1, 5}. Note that the function f defined in (3) can be seen as another way of writing (1).

Figure 2 displays some observations generated using the above simulation scheme for each model and for each σ. We can see from this figure that the clustering problem associated to Model 1 seems to be the most difficult. In Model 1, the clusters are indeed completely mixed whereas in Model 2 Cluster C 4 is well separated from the others. Observe also that the data that is generated has the same behavior as the data coming from our motivating application: They are nondecreasing and piecewise linear constant with a small additive noise, see Figure 1.

3.2. Statistical performance. Following the simulation scheme described in Section 3.1, the performance of our procedure is assessed for each model, each σ and is compared with two different clustering methods: the k-means algorithm applied to the raw data Y and the FunFEM procedure described in [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] and available in the R package FunFEM. The latter method is dedicated to the clustering of functional data and is based on a functional mixture model. All the methods are compared thanks to the Adjusted Rand Index (ARI) defined in [START_REF] Hubert | Comparing partitions[END_REF] which is often used for clustering validation. It is indeed a measure of agreement between two partitions. Note that the number of clusters k in the k-means algorithm is chosen using the same strategy as the one that we considered in our approach. As far as FunFEM is concerned, we used the default parameters.

For each model and for each σ in {1, 5}, we repeat independently 100 times the following steps:

(a) We simulate a sample D N = {(Y 1 , Z 1 ) . . . (Y N , Z N )} of size N = 100 according to the scheme described in Section 3.1; (b) We apply each method to D N ; (c) Based on the obtained clustering, we compute the ARI. The results are displayed in Figure 3 with K max = 10. We can see from this figure that our method outperforms the other ones in all cases except for Model 2 with σ = 5 where the performance of our method is on a par with the one of FunFEM. Note that applying the k-means to a relevant summary measure of Y significantly improves the clustering performance. Moreover, we observe that when σ increases, the performance of our approach is slightly altered since the change-points are more difficult to locate accurately, see Section 3.3.

3.3.

Assessment of our change-point estimation procedure. We provide the following numerical experiments for assessing the change-point estimation stage of our method. We used the parameters associated to Cluster 3 of Model 1, see Table 1. We repeat 100 times (a) We simulate Y according to Equation (4) with σ ∈ {1, 5}; (b) We estimate the change-points according to the procedure described in the first stage of the first step in Section 2. q q q q q q 0.2 0.4 0.6 0.8 kmeans Our method FunFEM ARI q q q q q q 0.2 Some examples of Y for the two values of σ are displayed in Figure 4. We can see from this figure that the change-points located at 300 and 400 are more difficult to detect than the others. It is all the more true when σ = 5.

Figure 5 displays the frequency of the number of times where each position has been estimated as a change-point. We can see that the change-points are all retrieved and that no spurious change-points are provided when σ = 1. In the case where σ = 5, although the

Application

In this section, we apply the methodology described in Section 2 to milking kinetics of dairy goats coming from the experimental herd of the research unit Systemic Modelling Applied to Ruminants (Paris, France). 4.1. Data description. The data set contains 100470 milking kinetics of goats of two different breeds: "Alpine" and "Saanen". All these kinetics are morning milking kinetics and several kinetics are available for each goat. The kinetics can also be separated according to parity which corresponds to the lactation rank i.e. to the number of times a goat has given birth and started a new lactation. In the considered dataset, there are in particular 276 (resp. 191) goats for which we have their milking kinetics for Parity 1 (resp. 2). 4.2. Kinetics clustering. First, note that based on the shapes of the milking kinetics of this data set, the parameter K max defined in the first stage of the first step in Section 2 was set to 2. We obtained three clusters containing 57498, 36757 and 6215 kinetics, respectively. Some examples of kinetics belonging to Clusters 1, 2 and 3 are displayed in Figures 6,7 and 8, respectively. The average of the kinetics estimations obtained within each cluster is displayed in Figure 9. We can observe that the three clusters can be distinguished in terms of quantity of milk production: Cluster 1 has the lowest production, Cluster 3 the highest and Cluster 2 is between them.

Another difference between the three clusters is the number and the positions of changes. The number of changes in the kinetics of Cluster 1 and 2 is mainly one contrary to Cluster 3 where this number is always equal to two. Figure 10 displays the histogram of the changepoint positions for Clusters 1 and 2. We can observe that the change-point having the highest frequency is not located at the same position for these two clusters. Interestingly, our methodology was able to distinguish these two clusters thanks to the change-point position which illustrates the potential of our methodology to extract synthetic traits from raw data.

In practice, such a clustering may be very useful in the precision farming context to refine selection criteria for breeding programs, to simplify milking workload or to control udder health. Thanks to the clustering results, we should be able to define a milking profile for each goat. Moreover, we propose in the next section to characterize dairy goats belonging to a given parity.

Parity characterization.

In order to go further into this analysis, we tried to characterize the parities 1 and 2 in terms of the proportion of kinetics of type 1, 2 or 3 according to the clustering previously obtained. We thus created for each goat belonging to a given parity a vector of proportions corresponding to its belonging frequency to each Cluster 1, 2 or 3. For each parity, the goats are clustered using the k-means algorithm applied to the vectors of proportions. The results are displayed in Figures 11 and12 for Parities 1 and 2, respectively. The number of groups is selected using the method described in Section 2.2: We found 6 (resp. 5) groups for Parity 1 (resp. 2). We can notice that there is one goat which produces a large quantity of milk compared to the others for both parities. In

Goat 250053006212058

RMS residual = 20.7 q q q q q q q q q q q q q q 0 50 100 200 0 500 1500 Goat 250053006212048 1 Volume of milk

Goat 250053006212048

RMS residual = 28.8 q q q q q q q q q q q q q q q q q q q qq qq qqqq

Goat 250053006212064

RMS residual = 42.67 q q q q q q q q qq q q qqq qqqqqqqqq qqqqqqqqqqqqqqqqq 0 20 40 60 80 100

Goat 250053006212083

1 Volume of milk

Goat 250053006212083

RMS residual = 22.8 q q q q q q q q q q q q 0 100 200 300 0 500 1500 Goat 250053006213036

1 Volume of milk

Goat 250053006213036

RMS residual = 37.04 q q q q q q q q qq qqqqq qqq qqqqqqqqqqqqqqqqqqqqq

Goat 250053006213031

RMS residual = 55.93 qq q q q q q q q q q q q q q q q q qqqqqqqqq 0 50 100 150

Goat 250053006213005

1 Volume of milk

Goat 250053006213005

RMS residual = 23.5 q q q q q q q q q q q q q q q q q q q 0 50 100 150 200 0 500 1500 Goat 250053006213008 1 Volume of milk

Goat 250053006213008

RMS residual = 31.2 q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 0 500 1000

Goat 250053006213030

1 Volume of milk

Goat 250053006213030

RMS residual = 22 q q q q q q q q q q q q q q q q q Figure 6. Some examples of milking kinetics belonging to Cluster 1. The data are displayed with 'o', the straight lines correspond to the piecewise linear fit obtained thanks to our method and the vertical line corresponds to the position of the change-point.

Goat 250053006213082

RMS residual = 42.1 q q q q q q q q q q q q q q qqqq q qqqqqq qqq 0 50 150 250 350 0 1000 2000

Goat 250053006212047

2 Volume of milk

Goat 250053006212047

RMS residual = 24 q q q q q q q q q q q q q q q q q qq qqqq q qqqqqqqqqqqq 0 50 100 200 0 500 1500

Goat 250053006213040

2 Volume of milk

Goat 250053006213040

RMS residual = 36.39 q q q q q q q qq q q q q q q q q q qq qq qqqq qqq 0 50 150 250

Goat 250053006212054

2 Volume of milk

Goat 250053006212054

RMS residual = 26.4 q q q q q q q q q q q q q q q q q q qq qqqq q qqqq q 0 50 100 150 0 500 1500 Goat 250053006213039 2 Volume of milk

Goat 250053006213039

RMS residual = 26.9 q q q q q q q q q q q q q q q q q 0 50 100 150 200 250 0 1000 2000

Goat 250053006212075

2 Volume of milk

Goat 250053006212075

RMS residual = 38.92 q q q q q q q q q q q q qqq qqqqqq qqqqq 0 50 150 250

Goat 250053006212084

2 Volume of milk

Goat 250053006212084

RMS residual = 39.65 q q q q q q q q q q q q q q q q q q qq qqqq qqqq qq

Goat 250053006212080

RMS residual = 34.98 q q q q q q q q q q q q q qqqqqqq q qqqqqqqqqqqqqqq

Goat 250053006212091

RMS residual = 31.8 q q q q q q q q q q q q q q q q q q q q q qqq qqqqqqq qqq

Goat 250053006212104

RMS residual = 21 q q q q q q q q q q q q q q q q q q q q q q q q q qq 0 50 150 250 350 0 500 1500 Goat 250053006212103 3 Volume of milk

Goat 250053006212103

RMS residual = 18.6 q q q q q q q q q q q qq q q q q q q q q q q q q q q q qqq qqqq 0 50 150 250 0 1000 2000 3000

Goat 250053006212105

3 Volume of milk

Goat 250053006212105

RMS residual = 26 q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq qq 0 50 100 150

Goat 250053006213006

3 Volume of milk

Goat 250053006213006

RMS residual = 6 q q q q q q q q q q q q q q q q 0 50 100 150 0 500 1000 1500

Goat 250053006213022

3 Volume of milk

Goat 250053006213022

RMS residual = 16.5 q q q q q q q q q q q q q q q q q q q 0 100 200 300 400 0 500 1500 Goat 250053006213073 3

Volume of milk

Goat 250053006213073

RMS residual = 28.8 q q q q q q q q q q q q q q qqqqqqqqq qqq q qq q qq q q q qqqqq qq 0 50 100 150 200

Goat 250053006213066

3

Volume of milk

Goat 250053006213066

RMS residual = 18 q q q q q q q q q q q q q q q q q q q q q q q 0 50 100 150 200 0 500 1000

Goat 250053006213109

3 Volume of milk

Goat 250053006213109

RMS residual = 16.7 qq q q q q q q q q q q q q q q q q q q qq qq 0 20 40 60 80 100 0 500 1500

Goat 250053006213124

3 Volume of milk

Goat 250053006213124

RMS residual = 9 q q q q q q q q q q q Figure 8. Some examples of milking kinetics belonging to Cluster 3. The data are displayed with 'o', the straight lines correspond to the piecewise linear fit obtained thanks to our method and the vertical line corresponds to the position of the change-points. Parity 1, 80% of its milking kinetics belong to Cluster 2 and only 20% to Cluster 1. In Parity 2, 100% of its milking kinetics belong to Cluster 2.

We also observe from Figures 11 and 12 that in both parities, the belonging frequency of the milking kinetics to Cluster 2 is between 50% and 70%. In Parity 2, there is one group (in red) for which the proportion of milking kinetics belonging to Cluster 2 is very high (around 65%) and the proportions of milking kinetics belonging to Cluster 1 and Cluster 3 are very low (around 25% and 13%, respectively). For the other groups the proportions of milking kinetics belonging to Cluster 1 are higher. In Parity 1, the behavior is a little bit different in the sense that the majority of goats have a high proportion of milking kinetics belonging to Cluster 2 (around 65%) and a low proportion of milking kinetics belonging to Cluster 1 (around 25%). Such results may be interesting in the context of precision breeding since they could help to forecast the production of milk at the different parities.

Further analysis should be perfomed in the future to study how evolve the cluster belonging along the lactation course lasting around 150 days in goats. The daily milk yield of a goat for a given parity follows indeed a typical triphasic shape (respectively increasing, plateau and decreasing phase), each daily milk yield being the sum of the total milk produced during each milking (respectively morning and afternoon milking). Being able to link a particular shape at the milking kinetics scale with one at the lactation scale could open perspectives to better characterize individual goats and thus propose options for individual milking management. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q 0.20 0.25 0.30 0.35 0.50 0.60 0.70 0.80

Proportion of kinetics 1 Proportion of kinetics 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q 0.50 0.60 0.70 0.80 0.00 0.05 0.10 0.15 0.20 Proportion of kinetics 2 Proportion of kinetics 3 q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.20 0.25 0.30 0.35 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.5 0.6 0.7 0.8 0.9 1.0 Proportion of kinetics 1 Proportion of kinetics 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 0.6 0.7 0.8 0.9 1.0 0.00 0.05 0.10 0.15 Proportion of kinetics 2 Proportion of kinetics 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 0.00

Figure 1 .

 1 Figure 1. Some examples of milking kinetics of goats.

Figure 2 .

 2 Figure 2. Examples of observations generated from Model 1 (top) and Model 2 (bottom) for σ = 1 (left) and σ = 5 (right). The curves belonging to Cluster 1 (resp. 2, 3, 4) are displayed in red (resp. black, blue and green). The solid lines display the representative curves of each cluster f t z ,θ z and the dashed ones are some examples of the corresponding Y.

Figure 3 .

 3 Figure 3. Boxplots of the ARI for Model 1 (top) and Model 2 (bottom) for σ = 1 (left) and σ = 5 (right).

Figure 4 .Figure 5 .

 45 Figure 4. Examples of Y belonging to Cluster 3 of Model 1 for σ = 1 (left) and σ = 5 (right).

Figure 7 .

 7 Figure7. Some examples of milking kinetics belonging to Cluster 2. The data are displayed with 'o', the straight lines correspond to the piecewise linear fit obtained thanks to our method and the vertical line corresponds to the position of the change-point.

Figure 9 .Figure 10 .

 910 Figure 9. Kinetics average obtained within each of the three clusters.

Figure 11 .

 11 Figure 11. Clustering obtained for goats in Parity 1 (6 clusters) displayed on the plane having for axes the proportion of kinetics belonging to Clusters 1 and 2 (top left), 2 and 3 (top right), 1 and 3 (bottom). The Saanen (resp. Alpine) goats are displayed with '•' (resp. '+').

Figure 12 .

 12 Figure 12. Clustering obtained for goats in Parity 2 (5 clusters) displayed on the plane having for axes the proportion of kinetics belonging to Clusters 1 and 2 (top left), 2 and 3 (top right), 1 and 3 (bottom). The Saanen (resp. Alpine) goats are displayed with '•' (resp. '+').

Table 1 .

 1 Set of parameters for Model 1.

			Model 1
	z K z	t z	θ z
	1 2	(150, 250)	(1600, 1900, 2000)
	2 2	(150, 300)	(1400, 1800, 2200)
	3 4 (100, 200, 300, 400) (300, 1500, 1700, 2000, 2200)
	4 3	(50, 150, 300)	(200, 1300, 1800, 2100)

Table 2 .

 2 Set of parameters for Model 2. Finally, we define Y such that, for i ∈ {1, . . . , n},

		Model 2
	z K z	t z	θ z
	1 2	(150, 250)	(1600, 1900, 2000)
	2 2	(150, 300)	(1400, 1800, 2200)
	3 4 (100, 200, 300, 400) (300, 1500, 1700, 2000, 2200)
	4 3	(150, 250, 300)	(200, 700, 1000, 1600)


Innocent Nomel

Jérôme Lelievre

Laurence Boyer

Olivier Durand-Panteix

Pascal Marchet
email: pascal.marchet@unilim.fr

Synthesis of 0.94 Na0.5Bi0.5TiO3 -0.06 BaTiO3 (NBT-6BT) lead-free piezoelectric powder suitable for aerosol deposition (AD)

Keywords: lead-free piezoelectric material, perovskite, piezoelectric properties

The present work is devoted to the synthesis of 0.94 (Na0.5Bi0.5)TiO3 -0.06 (BaTiO3) powders by solid-state reaction, suitable for Aerosol Deposition method (AD). The powders were synthetized using Na2CO3, Bi2O3, BaCO3 and TiO2. The raw powders were simultaneously milled / mixed by attrition milling in water. Two routes were studied, differing by the order of introduction of the raw materials, for three calcination temperatures. The first route (S1) led to the occurrence of Bi2Ti2O7 pyrochlore secondary phase. Attempts to suppress this phase by annealing were unsuccessful. The second route (S2) gave no secondary phase. The large grain size obtained at 950°C and 1050°C are not suitable for sintering or AD. However, the process carried out at 850°C resulted in powders suitable for both sintering and AD method. The ferroelectric and piezoelectric properties of the ceramic samples obtained using these powders are similar to the best-reported ones. Finally, we successfully used these powders for AD.

Introduction

Currently, piezoelectric ceramic materials are used in a large range of sensors or actuators. We can, for example, mention everyday life devices, such as car parking aid and motor fuel injectors, piezoelectric lighters, piezoelectric microphones and loudspeakers… or industrial and medical ultrasonic applications like non-destructive testing, underwater sonar systems or medical echography… Most of these applications need piezoelectric materials in ceramic form. However, some of them, like piezoelectric sensors or energy-harvesting devices, require thick film shape (1 -100 µm). Only a few elaboration methods are available for this shaping: screen-printing and aerosol deposition method (AD) [START_REF] Akedo | Aerosol Deposition (AD) and Its Applications for Piezoelectric Devices[END_REF]. Introduced around year 2000 by Dr J. Akedo [START_REF] Akedo | Aerosol Deposition (AD) and Its Applications for Piezoelectric Devices[END_REF][START_REF] Akedo | Fabrication of three dimensional micro structure composed of different materials using excimer laser ablation and jet molding[END_REF][START_REF] Akedo | Chapter 10.5 -Aerosol Deposition Method for Room-Temperature Ceramic Coating and Its Applications[END_REF], AD appears now as an interesting emerging process in ceramic technologies, currently co-developed in France by the Center for Transfer of Ceramic Technologies (CTTC) and the Institute of Research on Ceramics (IRCER) since 2010 [START_REF] Piechowiak | Growth of dense Ti3SiC2 MAX phase films elaborated at room temperature by aerosol deposition method[END_REF][START_REF] Henon | Dense and highly textured coatings obtained by aerosol deposition method from Ti3SiC2 powder: Comparison to a dense material sintered by Spark Plasma Sintering[END_REF]. In AD, a dry aerosol of ceramic powder, sent by a carrier gas, is projected onto a substrate using a suitable supersonic nozzle. By fragmentation and rearrangement of the ceramic grains, this process leads directly to a consolidated ceramic film [START_REF] Akedo | Chapter 10.5 -Aerosol Deposition Method for Room-Temperature Ceramic Coating and Its Applications[END_REF]. Unlike bulk ceramics, the main advantage of AD is that it does not require sintering. However, one of the key point of this process is the quality of the ceramic powder. Among the necessary conditions, AD use submicronic or micronic powders, generally between 0.8 µm and 10 µm [START_REF] Akedo | Chapter 10.5 -Aerosol Deposition Method for Room-Temperature Ceramic Coating and Its Applications[END_REF]. Therefore, the elaboration of piezoelectric devices by AD associated to our studies requires a controlled powder.

Nowadays, the most widely used piezoelectric materials are PbZr(1-x)TixO3 (PZT) ceramics [START_REF] Jaffe | Piezoelectric Ceramics[END_REF][START_REF] Waanders | Piezoelectric Ceramics-Properties and Applications, Philips Components[END_REF][START_REF] Moulson | Electroceramics: Materials, Properties and Applications[END_REF][START_REF] Setter | Piezoelectric Materials and Devices[END_REF]. Indeed, PZT are easy to synthesize and to sinter. Moreover, their electromechanical properties can be easily modulated by doping, allowing a large range of properties and applications, from soft to hard piezoelectric materials. However, due to environmental and health concerns, lead-based materials should be avoided in the future [10,[START_REF][END_REF]. Thus, it is necessary to search for lead-free piezoelectric materials suitable for the substitution of PZT, at least for some applications. The most actively investigated lead-free piezoelectric materials currently considered as potential substitutes for PZT are: (i) Aurivillius phases, such as Bi4Ti3O12, (ii) barium titanate BaTiO3 (BT) and its derivatives such as (Ba, Ca)(Ti, Zr)O3 (BCTZ), (iii) alkali-bismuth titanates (Na0.5Bi0.5)TiO3 (NBT) or (K0.5Bi0.5)TiO3 (KBT) and (iv) alkali niobates (K, Na)NbO3 (KNN) and their derivatives [START_REF] Maeder | Lead Free Piezoelectric Materials[END_REF][START_REF] Takenaka | Lead-free piezoelectric ceramics based on perovskite structures[END_REF][START_REF] Zhang | Lead-free piezoelectric ceramics vs. PZT?[END_REF][START_REF] Panda | Review: environmental friendly lead-free piezoelectric materials[END_REF][START_REF] Rödel | Perspective on the Development of Lead-free Piezoceramics[END_REF]. Currently, none of these materials achieves the performance of PZT. However, despite this handicap, some of these lead-free materials have interesting and sufficient electrical properties for applications.

Among them, sodium bismuth titanate is considered as a potential candidate to the substitution of PZT. NBT is easily synthesized, either by the conventional solid-state route [START_REF] Elkechai | Structural and dielectric study of the Na0.5Bi0.5TiO3-PbTiO3 and K0.5Bi0.5TiO3-PbTiO3 systems[END_REF][START_REF] Said | Relaxor behaviour of low lead and lead free ferroelectric ceramics of the Na0.5Bi0.5TiO3-PbTiO3 and Na0.5Bi0.5TiO3-K0.5Bi0.5TiO3 systems[END_REF][START_REF] Hiruma | Thermal depoling process and piezoelectric properties of bismuth sodium titanate ceramics[END_REF], by the hydrothermal / solvothermal route [START_REF] Liu | Hydrothermal synthesis of monosized Bi0.5Na0.5TiO3 spherical particles under low alkaline solution concentration[END_REF][START_REF] Lu | Synthesis of Na0.5Bi0.5TiO3 powders through hydrothermal method[END_REF] or by the sol-gel route [START_REF] Rémondière | Synthesis and crystallization pathway of Na0.5Bi0.5TiO3 thin film obtained by a modified sol-gel route[END_REF][START_REF] Šetinc | Temperature stable dielectric behavior of sol-gel derived compositionally graded SrTiO3 / Na0.5Bi0.5TiO3 / SrTiO3 thin films[END_REF]. However, NBT is difficult to pole due to its large coercive field (Ec=7.5 kV/mm [START_REF] Spreitzer | Sodium deficiency in Na0.5Bi0.5TiO3[END_REF]). Furthermore, NBT presents some shortcomings like narrow sintering temperature range and high conductivity at room temperature, connected to significant dielectric losses [START_REF] Said | Relaxor behaviour of low lead and lead free ferroelectric ceramics of the Na0.5Bi0.5TiO3-PbTiO3 and Na0.5Bi0.5TiO3-K0.5Bi0.5TiO3 systems[END_REF][START_REF] Nagata | Additive effects on electrical properties of (Bi1/2Na1/2)TiO3 ferroelectric ceramics[END_REF][START_REF] Lin | Effects of doping Eu2O3 on the phase transformation and piezoelectric properties of Na0.5Bi0.5TiO3-based ceramics[END_REF]. These result from the vaporization of bismuth and/or sodium during synthesis and sintering and therefore from the creation of oxygen vacancies which fix the domain walls [START_REF] Nagata | Electrical properties and tracer diffusion of oxygen in some Bi-based lead-free piezoelectric ceramics[END_REF]. The reported piezoelectric properties are d33=78 pC/N, d31=15 pC/N, d15=87 pC/N [START_REF] Hiruma | Thermal depoling process and piezoelectric properties of bismuth sodium titanate ceramics[END_REF]. Therefore, the applications of pure NBT ceramics are also limited because of moderate values of piezoelectric charge coefficients, despite the interesting value of coupling coefficient k33=0.47 [START_REF] Hiruma | Thermal depoling process and piezoelectric properties of bismuth sodium titanate ceramics[END_REF].

Fortunately, like other perovskite ferroelectric compounds, the properties of NBT can be easily modified by the formation of solid solutions [START_REF] Said | Relaxor behaviour of low lead and lead free ferroelectric ceramics of the Na0.5Bi0.5TiO3-PbTiO3 and Na0.5Bi0.5TiO3-K0.5Bi0.5TiO3 systems[END_REF][START_REF] Sakata | Ferroelectric and Antiferroelectric Properties of (Na0.5Bi0.5)TiO3-SrTiO3 Solid Solution Ceramics[END_REF][START_REF] Takenaka | Acoustic wave characteristics of lead-free (Bi1/2Na1/2)0.99Ca0.01TiO3 piezoelectric ceramic[END_REF][START_REF] Takenaka | Piezoelectric ceramics of (Bi1/2Na1/2)TiO3-PbTiO3-BaTiO3 system[END_REF][START_REF] Elkechai | Na0.5Bi0.5TiO3-K0.5Bi0.5TiO3 (NBT-KBT) system: A structural and electrical study[END_REF][START_REF] Nagata | Lead-free piezoelectric ceramics of (Bi1/2Na1/2)TiO3-1/2(Bi2O3•Sc2O3) system[END_REF][START_REF] Marchet | Dielectric properties of some low-lead or lead-free perovskite-derived materials: Na0.5Bi0.5TiO3-PbZrO3, Na0.5Bi0.5TiO3-BiScO3 and Na0.5Bi0.5TiO3-BiFeO3 ceramics[END_REF][START_REF] Dorcet | Structural and dielectric studies of the Na0.5Bi0.5TiO3-BiFeO3 system[END_REF]. For the (1-x) (Na0.5Bi0.5)TiO3 -x(BaTiO3) system, hereafter called NBT-xBT, the structural situation is similar to PZT. Indeed, PZT and NBT-xBT are solid solutions between a rhombohedral (PbZrO3 or NBT) and a tetragonal (PbTiO3 or BT) compound. In both case, structural studies evidenced a morphotropic phase boundary (MPB) at the border between these domains, where properties are enhanced. For NBT-xBT, the published results report a MPB composition around x=0.06-0.07 [START_REF] Takenaka | Bi1/2Na1/2)TiO3-BaTiO3 system for lead-free piezoelectric ceramics[END_REF][START_REF] Chu | Electrical properties of Na1/2Bi1/2TiO3-BaTiO3 ceramics[END_REF][START_REF] Gomah-Pettry | Sodium-bismuth titanate based lead-free ferroelectric materials[END_REF][START_REF] Qu | Effect of A-site substitution on crystal component and dielectric properties in Bi0.5Na0.5TiO3 ceramics[END_REF][START_REF] Ranjan | Structure and dielectric properties of (Na0.50Bi0.50)1-xBaxTiO3: 0≤x≤0[END_REF] where the piezoelectric properties are improved: d33=125 pC/N, d31=40 pC/N, d15=194 pC/N and k33=0.55 [START_REF] Takenaka | Bi1/2Na1/2)TiO3-BaTiO3 system for lead-free piezoelectric ceramics[END_REF]. Thus, we retained the NBT-6BT composition for our work on the elaboration of lead-free piezoelectric thick films by AD. Indeed, this material has interesting electrical properties for piezoelectric applications and microelectronics.

The literature reports several methods for synthesizing NBT-6BT powder, in particular the sol-gel route, the hydrothermal route and the solid-state route. Nevertheless, the two first routes are expensive and chemically complicated. Only the solid-state route appears suitable for the synthesis of large amounts of ceramic powders needed for AD. For reasons of cost and simplicity, we retained co-attrition milling of the different raw materials rather than separate milling. Indeed, such process allows a one-step reduction of their grain sizes, with simultaneous homogenization of the mix. For the same reason, we selected water as solvent for attrition because of its low cost and low environmental imprint.

Unfortunately, if many papers devoted to the properties of NBT-6BT are available [START_REF] Takenaka | Bi1/2Na1/2)TiO3-BaTiO3 system for lead-free piezoelectric ceramics[END_REF][START_REF] Chu | Electrical properties of Na1/2Bi1/2TiO3-BaTiO3 ceramics[END_REF][START_REF] Gomah-Pettry | Sodium-bismuth titanate based lead-free ferroelectric materials[END_REF][START_REF] Qu | Effect of A-site substitution on crystal component and dielectric properties in Bi0.5Na0.5TiO3 ceramics[END_REF][START_REF] Ranjan | Structure and dielectric properties of (Na0.50Bi0.50)1-xBaxTiO3: 0≤x≤0[END_REF][START_REF] Chen | Structure and electrical properties of (Na0.5Bi0.5)1-xBaxTiO3 piezoelectric ceramics[END_REF][START_REF] Cordero | Phase transitions and phase diagram of the ferroelectric perovskite (Na0.5Bi0.5)1-xBaxTiO3 by anelastic and dielectric measurements[END_REF][START_REF] Ma | Phase diagram of unpoled lead-free (1-x)(Bi1/2Na1/2)TiO3-xBaTiO3 ceramics[END_REF][START_REF] Ma | Domain structure-dielectric property relationship in lead-free (1-x) (Bi1/2Na1/2)TiO3-x BaTiO3 ceramics[END_REF][START_REF] Jo | Evolving morphotropic phase boundary in lead-free (Bi1/2Na1/2)TiO3-BaTiO3 piezoceramics[END_REF][START_REF] Craciun | Electric-field-induced and spontaneous relaxorferroelectric phase transitions in (Na1/2Bi1/2)1-xBaxTiO3[END_REF][START_REF] Ma | Creation and destruction of morphotropic phase boundaries through electrical poling: A case study of lead-free (Bi1/2Na1/2)TiO3-BaTiO3 piezoelectrics[END_REF][START_REF] Fuentes-Cobas | The 0.96(Bi0.5Na0.5)TiO3 -0.04BaTiO3 crystal structure: A high-Q, high-counting statistics synchrotron diffraction analysis[END_REF], none of them reports a systematic study of the synthesis. Therefore, the aim of the present work is to study the impact of the milling and calcination processes (i) on secondary phases, (ii) on grain size, powder morphology and sintering ability of NBT-6BT powders processed by classical solid-state route and (iii) to test their suitability for the elaboration of films by AD.

Materials and methods

Powder synthesis

Powders of 0.94 Na0.5Bi0.5TiO3 -0.06 BaTiO3 composition were elaborated through solid-state reaction. The selected precursors are Na2CO3 (Aldrich 99.5%), Bi2O3 (Alfa-Aesar 99.9995%), BaCO3 (Alfa-Aesar 99.8%) and TiO2 (Alfa-Aesar 99.8%). The powders were weighed in stoichiometric amount and then homogenized by attrition milling in water, using a polymethacrylate dispersant. The process was carried out at 1000 rpm with 1 mm diameter yttria doped zirconia beads. Two different milling / mixing processes were used, hereafter called S1 and S2, differing in the order of addition of the precursors.

For S1, the idea was to introduce the carbonates first, since Na2CO3 have a coarser particle size than the oxides. Barium carbonate is introduced first and treated for 30 min, followed by Na2CO3, then milled for the same time. After 1 hour, Bi2O3 and TiO2 are added simultaneously to the mix and the whole is further treated for 2 hours. Thus, the order of introduction is BaCO3

→ Na2CO3 → [Bi2O3 + TiO2].
For S2 route, the idea was to introduce the precursors according to their average particle size (D50) in decreasing order. However, since Na2CO3 is soluble in water, we therefore chose to insert chronologically Bi2O3 → BaCO3 → [Na2CO3 + TiO2]. Bi2O3 powder is introduced first and attrition milled for 30 min. Then, BaCO3 is introduced and the milling continues for 30 min. After 1 hour milling, the combination [Na2CO3 + TiO2] is added to the mix. The attrition treatment continues for another 2 hours.

The obtained raw powders were dried overnight at 100°C and then calcined for 10 hours at three different temperatures: 850°C, 950°C and 1050°C (heating rate 5°C/min). These values were selected using results obtained from the literature. The aim was to determine the lowest possible temperature to obtain suitable powders, in order to limit the grain growth, the bismuth evaporation and the energetic cost. After synthesis and characterization, we also examined the effect of annealing on the powders for 1 hour at the same temperature, in order to estimate the influence of such thermal treatment on phase purity and grain growth. Before annealing, the powders were attrition-milled for one hour at 1000 rpm in order to decrease the grain size and the amount of agglomerates.

Characterization techniques

The grain size distribution of the powders was determined using a laser particle size analyzer (Horiba LA 950 V2 Pratica). The measurements were carried out in water. The grain size distribution was extracted from the measurements considering the Fraunhofer approximations. From the literature, the refractive index of NBT-6BT is estimated to be between 2.35 and 2.55. The same device was used to follow the grain size evolution during attrition milling and also after synthesis and annealing of the powders. D10, D50 and D90 are the mean diameters associated respectively to 10%, 50% and 90% of the statistical population.

The morphology of the powders was examined by scanning electron microscopy (SEM, LEO 1530 VP). The structural characterization of the powders and phase identification was performed by X-ray diffraction (XRD), using a θ / 2θ diffractometer (Bruker D8) with Cu Kα1 wavelength (λ = 1.5406 Å). The crystalline phases were determined from XRD patterns recorded for 60 min, for 2θ values ranging from 20° to 80° with a step of 0.02°.

The density of the powders was determined at room temperature using a pycnometer (Accupyc II 1340, Micrometrics). For the best obtained powders, dilatometric measurements were performed up to 1160°C (5 min dwell time), with 5°C/min heating and cooling rate (Setaram, TMA setsys evo 1600).

The ferroelectric cycles of the ceramic samples were obtained at room temperature using the aixACCT device fitted with piezo sample holder unit (PSHU, Aixacct, Germany). The converse d33* piezoelectric coefficient was determined using the slope of the piezoelectric strain cycle. The dielectric properties and resonance curves of the ceramic samples after poling were measured using an impedance meter (Agilent 4294A) between 10 kHz and 5 MHz. The piezoelectric coupling coefficients kp and kt were calculated by the IEEE standard method using the resonance / antiresonance frequencies determined by impedance curves.

Results and discussions

Characterization of the precursors

The particle size measurements and the morphology of raw precursors obtained by SEM are represented in Fig. 1. The results evidence that Na2CO3 present very large grain size (484 µm) with single mode distribution. For BaCO3, the grains present irregular morphology with a mean diameter around 1.6 µm. Bi2O3 powder is associated to a bimodal distribution, showing both faceted large grains (3.3 µm) and irregular small particles (0.3 -0.4 µm). Finally, TiO2 grains are smaller (0.4 µm) than all the other raw materials. The large difference of the measured grain size (D50) confirms the requirement of milling in order to decrease grain size and thus increase reactivity before calcination.

Processing of the NBT-6BT mixture before calcination

Attrition milling and evolution of particle size during attrition-milling

As mentioned hereabove, the synthesis route S1 differs from S2 in the order of addition of the precursors in the attrition system. During the process, samples were drawn from the mix every 15 min, in order to measure the grain size distribution. Each time a precursor is added, a sample was drawn just before the addition (T0 + 30, T1 + 30) and just after (T1, T2). The grain size evolution for S1 and S2 routes (D10, D50 and D90) is presented in Fig. 2 and Fig. 3. For S1, the results evidence an increase of the grain size after the introduction of Na2CO3 (D50 and D90). Since Na2CO3 presents high solubility in water while BaCO3 is only weakly soluble, this grain size increase can be attributed to the precipitation of sodium carbonate, probably onto the surface of barium carbonate. Thus, this phenomenon is detrimental to the chemical homogeneity of the powder. We did not observed a similar phenomenon for S2 after the introduction of [Na2CO3 + TiO2]. Therefore, S2 route will probably induce better chemical homogeneity than S1.

However, for both S1 and S2 routes and after 3 hours of attrition milling, an average grain size (D50) around 1 μm is obtained. The large particles observed before milling, especially for Na2CO3 (D50 around 480 µm) are no more present. At this stage of the powder processing, there is no obvious difference in grain size between the S1 and S2 routes and the D50 is similar. This result evidences that both routes allow a reduction in the grain size toward the desired value, around 1 µm, as expected from attrition milling.

XRD and densities of calcined and annealed powders

The XRD patterns of the powders calcined at 850°C, 950°C and 1050°C for 10h are presented in Fig. 4, together with the patterns of the 1 hour annealed powders. As expected, all XRD patterns reveal a perovskite phase, similar to the one reported for NBT-6BT (PDF4 00-063-0299). We also detected a pyrochlore Bi2Ti2O7 secondary phase in different amounts, depending of the processing parameters. For S1 route, the amount of pyrochlore phase is large after the calcination at 850°C and decrease with increasing synthesis temperature ( The densities of powders resulting from S1 and S2 routes are reported in table 1. The relative values are only indicatives and were calculated considering a theoretical density of 5.95±0.01 g/cm 3 (reference file: PDF4 00-063-0299). Globally, the densities are higher for S1 route (5.86 -5.96) than for S2 route (5.66 -5.77). This reflects the occurrence of pyrochlore secondary phase, since Bi2Ti2O7 presents a higher density of 7.43 g/cm 3 (ISCD #99437) than NBT-6BT. For S1 route, the densities globally increase while increasing calcination temperature (5.87 → 5.94). The annealing globally does not change the density values, apart at 1050°C (5.94 → 5.86). This last result probably comes from bismuth evaporation, since bismuth volatility at high temperature is well known for Bi-containing compounds. For S2, the obtained densities (5.67) are lower than for S1. This result agrees well with the absence of secondary phase observed by XRD (only a small amount of this phase was detected at 850°C/10h). The measurements reveal an almost constant density (5.67) whatever the calcination temperature. Annealing increases the density values, probably indicating a better crystalline order after this further thermal treatment. Contrary to S1, annealing performed at 1050°C improves the density (5.67 → 5.72), probably indicating that bismuth losses are higher for the pyrochlore phase than for perovskite phase.

To sum up, S2 route leads to better results than S1, since no secondary phase appears from 950°C. A short annealing (1h) slightly increases the density, indicating a better crystalline order after annealing. Therefore, these results evidence that S2 route is favorable for obtaining the expected perovskite phase. Study of grain size and morphology will allow determination of the best calcination temperature.

Grain size and grain morphology of the calcined and annealed powders

Figure 5 represents the particle size distribution of the powders after calcination and annealing for both S1 and S2 routes and Figure 6 the corresponding SEM pictures. The analyses reveal a multimodal distribution of grain size for all samples. We remind that before annealing, the powders were attrition milled for one hour at 1000 rpm in order to decrease the amount of agglomerates.

The measurements of S1 powders (Fig. 5.a.) present micrometer-sized particles (1.2 -1.5 µm) together with larger grains (30 -45 µm) and very large ones assimilated to agglomerates (100 -200 µm). The syntheses at the two lowest calcination temperatures, 850°C and 950°C, lead to a high volume fraction of fine particles. The synthesis at 1050°C gives an equivalent volume of fine particles and larger grains. The corresponding SEM pictures (Fig. 6 a,b,c) confirm the particle distribution results. Indeed, the powder synthetized at 850C°/10h presents small grains, around 1 µm, associated to larger agglomerates composed of these grains. For 950°C, the elementary particles have a size with the same order of magnitude (1 µm), but their agglomeration degree increases. For the powders processed at 1050°C, the elementary grains are a little bit larger. Moreover, the grains are not only agglomerated, but a large degree of coalescence is noticed. In addition, we observe agglomerated round-shaped very large grains. This result clearly indicates that sintering has begun for this temperature (1050°C) and to some lower extent for 950°C.

The grain size of S1 annealed powders (Fig. 5.b) still presents a bimodal distribution. For 850°C, two populations are observed (0.45 and 1.5 µm), in agreement with the SEM picture (Fig. 6.d). For 950°C, most of the grain are measured around 1.5 µm, composed of smaller particles partly agglomerated (Fig. 6.e). Logically, for 1050°C (Fig. 6.f), the intermediate milling did not succeed in breaking the large and partly sintered grains observed after calcination.

For S2, the grain-size distribution curves (Fig. 5.c) present particles around 1 µm for calcinations at 850°C and 950°C. For 1050°C, the smallest population is around 13 µm. For all the temperatures, another population larger than 30 µm is also observable. For 850°C (Fig. 6.g), the grain size observed by SEM is around 1 µm, with only low agglomeration degree. For 950°C (Fig. 6.h) and 1050°C (Fig. 6.i), the increasing temperature induces an increase of the coalescence degree of the elementary grains. The bridging between these grains changes and bottleneck grain boundaries appear, as usually reported in the sintering process.

For the S2 sample annealed at 850°C (Fig. 5.d.), the short milling before thermal treatment allowed to reduce the grain size, since a new population is observed around 0.3 µm, while the second one remains around micrometer size (2.6 µm). The population at 77 µm is attributed to agglomerates, according to SEM picture (Fig. 6.j). For 950°C and 1050°C (Fig. 5.d.), the annealing lead to an increase of the grain size with two populations, at 13 µm and ~100 µm for 950°C, and around ~80 µm and 200 µm for 1050°C. These measurements agree well with the SEM pictures (Fig. 6.k and 6.l) which also show that annealing increases the coalescence degree compared to unannealed samples.

To summarize, the results of both grain size measurements and SEM pictures evidence that synthesis at 850°C leads to a grain size around 1µm. The intermediate milling followed by annealing keeps the same grain size, with low degree of agglomeration. As the AD process requires micronic powders, the process with 850°C calcination without annealing, seems to be the most appropriate for powders suitable for the deposition tests. For 950°C and 1050°C, the increase of the temperature lead to larger grain size, with a coalescence degree more and more pronounced with higher temperatures. The SEM pictures evidence bottleneck grain boundaries between rounded grains, showing that sintering has already begun. In addition, the rounded grain shape could suggest the occurrence of a liquid phase. Indeed, the TiO2 -Bi2O3 phase diagram indicates an eutectic point at 835°C for 97 mol% Bi2O3. Thus, it is likely that a transient liquid phase occurs during calcination.

Dilatometric study of the powders and sintering

As only the S2 route gave the NBT-6BT perovskite compound without secondary phase, dilatometric measurements were performed for the S2 powders calcined at 850°C, 950°C and 1050°C (Fig. 7). For these powders, the shrinkage associated to sintering begins around 950°C. Then, its rate increases around 1000°C, giving a maximum rate around 1100°C. However, the final shrinkage is different and depends on the calcination temperature of the powder. The final values are 18.8%, 15.9% and 7.2% respectively for 850°C, 950°C and 1050°C samples. They confirm the previous conclusion obtained by characterizations of powders. Indeed, since shrinkage begin around 950°C (Fig. 7), sintering already occurred during calcination for the 950°C/10h and 1050°C powders. Furthermore, for the 1050°C powder, the calcination led to rounded and large grains (10 -100 µm), with high degree of grain coalescence. Thus, this powder is unsuitable for sintering, as demonstrated by the low value of final shrinkage (7.2%). The coalescence degree of 950°C powders is lower, leading to a higher final shrinkage (15.9%). However, both 950°C and 1050°C powders give lower sintering results than the 850°C ones, as supposed with the previous characterizations. Therefore, only 850°C powders appear to be suitable for sintering, with a final shrinkage value compatible with ceramic densification (18.8%).

Electrical properties of ceramics

In order to evaluate the ferroelectric/piezoelectric properties of ceramics obtained by these powders, ceramic samples were sintered using the S2 850°C/10h powder. Before sintering, the powder was treated by planetary milling (300 rpm, 2x15 min). The powder was then pressed into pellets of 25 mm in diameter and around 1 mm in thickness with a uniaxial pressure of 125 MPa. The sintering was then performed at 1150°C for 10 h. After surface polishing, silver was deposited on the pellets by sputtering. The poling of the ceramic samples was performed in silicon oil at 80°C for 4 kV/mm electric field.

Figure 8 (a) shows a typical ferroelectric / piezoelectric cycle hysteresis loop for the NBT-6BT pellet. The remnant polarization (Pr) and the coercive electric field (Ec) obtained from the P(E) hysteresis loops are 53.5 µC/cm 2 and 2.5 kV/mm, respectively. The Pr value is slightly higher than reported values, Ec≈30-40 µC.cm -2 [START_REF] Wang | Piezoelectric and dielectric properties of CeO2added (Bi0.5Na0.5)0.94Ba0.06TiO3 lead-free ceramics[END_REF][START_REF] Prado | Stabilization of the morphotropic phase boundary in (1 - x)Bi0.5Na0.5TiO3-xBaTiO3 ceramics through two alternative synthesis pathways[END_REF][START_REF] Sundari | Influence of transition metal ions on multiferroic properties of lead-free NBT-BT ceramics[END_REF], while Ec value have the same order of magnitude Pr≈3 kV.mm -1 [START_REF] Wang | Piezoelectric and dielectric properties of CeO2added (Bi0.5Na0.5)0.94Ba0.06TiO3 lead-free ceramics[END_REF][START_REF] Prado | Stabilization of the morphotropic phase boundary in (1 - x)Bi0.5Na0.5TiO3-xBaTiO3 ceramics through two alternative synthesis pathways[END_REF]. Thus, these measured values are characteristics of a wellprepared powder and a well-densified ceramic. However, the polarization loop (Fig. 8 a) appears "rounded" for high fields. The saturation polarization is not fully reached, thus indicating a lossy behavior. This could be due to chemical defects induced by sintering, like oxygen vacancies. Indeed, the Bi-based materials are known as highly sensitive to bismuth losses during calcination and sintering. They induce oxygen vacancies which result in leakage current and pinning of the domain walls, associated to difficulty in poling process [START_REF] Park | Lanthanum-substituted bismuth titanate for use in non-volatile memories[END_REF][START_REF] Noguchi | Defect control for large remanent polarization in bismuth titanate ferroelectrics -doping effect of higher-valent cations[END_REF][START_REF] Noguchi | Impact of defect control on the polarization properties in Bi4Ti3O12 ferroelectric single crystals[END_REF]. Further experiments are in progress in order to improve the sintering and obtain better polarization loops.

Figure 8 (b) represents the S(E) strain -field unipolar response of the NBT-6BT pellet. A maximum strain of 0.11% was obtained for a 2.6 kV.mm -1 field (2500 V). This value makes them similar to the strain of ≈0.1% obtained for PZT ceramic materials [START_REF] Park | Ultrahigh strain and piezoelectric behavior in relaxor based ferroelectric single crystals[END_REF] and for undoped NBT-BT [START_REF] Zhang | Giant strain in lead-free piezoceramics Bi0.5Na0.5TiO3 -BaTiO3 -K0.5Na0.5NbO3 system[END_REF][START_REF] Mcquade | An investigation into group 13 (Al, Ga, In) substituted (Na0.5Bi0.5)TiO3-BaTiO3 (NBT-BT) lead-free piezoelectrics[END_REF]. The d33* converse piezoelectric coefficient, obtained from the piezoelectric cycle, is 358 pm.V -1 , very close to the best-reported values for undoped NBT-6BT ceramics (240-380 pm.V -1 ) [START_REF] Zhang | Giant strain in lead-free piezoceramics Bi0.5Na0.5TiO3 -BaTiO3 -K0.5Na0.5NbO3 system[END_REF][START_REF] Mcquade | An investigation into group 13 (Al, Ga, In) substituted (Na0.5Bi0.5)TiO3-BaTiO3 (NBT-BT) lead-free piezoelectrics[END_REF].

The measurements of the impedance curves around the resonances associated radial mode and thickness mode allowed the calculation of the coupling coefficients kp and kt. The calculated values are kp=0.22 and kt=0.37. These results agree well with the ones previously reported for NBT-6BT ceramic samples: kp=0.29 and kt=0.40 [START_REF] Takenaka | Bi1/2Na1/2)TiO3-BaTiO3 system for lead-free piezoelectric ceramics[END_REF]. These results confirm the high quality of the powder obtain by the S2 route, which appears thus suitable for AD.

Aerosol Deposition tests

Finally, we used the obtained powders for AD tests, using the apparatus reported in a previous paper [START_REF] Piechowiak | Growth of dense Ti3SiC2 MAX phase films elaborated at room temperature by aerosol deposition method[END_REF]. For these purpose, we retained a metallic Kovar substrate (thickness 0.125 mm) for their use in the future elaboration of piezoelectric devices. The XRD pattern of the NBT-6BT film is presented in Fig. 9. As a matter of comparison, the XRD patterns of the powder and of the kovar substrate are also presented, as well as the JCPDS card of NBT-6BT. The main XRD peak of NBT-6BT, observed around 32.5°, indicates that the film deposited by AD has the same perovskite structure as the powder. However, the XRD pattern of the AD film presents low diffracted intensity compared to that of the powder, together with a widening of the diffraction peaks. These phenomena clearly indicate a low grain size for the AD film compared to the powder. There is indeed a reduction in the grain size of the film linked to their distortion during deposition. The SEM top view (Fig. 9 b) evidences a rough surface, with occurrence of craters, one of the characteristic of the AD films. The SEM cross section obtained for a film prepared by Focused Ion Beam (Fig. 9. b) reveals a thickness of 38 µm. Despite the presence of micro porosities, the film presents a high density. These results indicate that we successfully obtained a NBT-6BT thick film by AD using our powder. The study of the properties of such films is under investigation and will be reported latter.

Conclusion

The aim of the present work was the study of the synthesis process of ceramic powders of lead-free piezoelectric material, devoted to the elaboration of piezoelectric thick films by aerosol deposition method (AD). For this purpose, we selected the 0.94 (Na0.5Bi0.5)TiO3 -0.06 (BaTiO3) (NBT-6BT) lead-free composition because of its piezoelectric properties. Since AD requires large batches of powder, we selected the classical solid-state route for the synthesis. Indeed, this method is the easiest and cheapest one in order to obtain large quantities of powder (100-500 g). The selected precursors are Na2CO3, Bi2O3, BaCO3 and TiO2 because of their low cost compared to precursors required by other methods.

Before calcination, stoichiometric amounts of powders were homogenized / milled by attrition milling. For the same reason of simplicity and cost, we retained a one-step milling in water, rather than separate milling of the precursors and a more expensive solvent. We studied two different routes for the milling, which differ by the order of introduction of the precursors. The best results were obtained using the S2 route, which corresponds chronologically to Bi2O3 → BaCO3 → [Na2CO3 + TiO2]. Indeed, this route gave better chemical homogeneity of the precursors mix before calcination. Moreover, it avoids the occurrence of the Bi2Ti2O7 pyrochlore secondary phase (at 950°C and 1050°C) or limits it at 850°C. These powders were calcined at three different temperatures for 10h: 850°C, 950°C and 1050°C. We also studied the influence of a post-calcination milling followed by a short annealing (1h) at the same temperature, in order to reduce the amount of secondary phase. The results evidenced that this elimination cannot be obtained when the chemical homogeneity is poor before calcination (S1 route) whereas this post treatment is efficient for S2 powders. The calcination and annealing at the highest temperature, 950°C/10h and 1050°C/10h, led to large grain size, with high coalescence of the grains. As demonstrated by the dilatometric measurements, this phenomenon comes from the pre-sintering of the NBT-6BT powder during calcination. Indeed, the shrinkage starts around 950°C with a maximum rate around 1100°C. Therefore, the best powders are the ones calcined at 850°C, lower that the starting sintering temperature.

In order to test the ability of the powders for ceramic elaboration and their ability for AD method, ceramic samples were sintered using these powders. The results of the ferroelectric / piezoelectric measurements, d33*=358 pm/V, kp=0.22 and kt=0.37 are similar to the best reported values for NBT-6BT. Because of these electrical properties, and since the grain size of this powder is around 1 µm, this synthesis process appears suitable for its use in aerosol deposition. Finally, we used these powders for AD test and obtained successfully NBT-6BT dense films with a thickness of tens of micrometers. The properties of these films are under investigation and will be reported latter. process Density (g/cm 3 1: Densities of the powders elaborated by S1 and S2 routes after calcination and annealing. The relative densities were calculated with a reference value of 5.95 ± 0.01 g/cm 3 (reference file: PDF4 00-063-0299)

  Fig 4.a). Further annealing at the same temperature for one hour did not succeed in the elimination of this secondary phase (Fig 4.b). This result confirms that S1 route did not succeed in the chemical homogenization of the precursors. For S2 route, we obtained only a minor amount of Bi2Ti2O7 secondary phase after the calcination at 850°C, without secondary phase for the other temperatures (Fig 4.c). A one hour annealing easily eliminates this secondary phase (Fig 4.d).

Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :Figure 8 :Figure 9 :

 123456789 Fig. 1 : Particle size distribution and SEM pictures of a) Bi2O3, b) BaCO3 c) Na2CO3 and d) TiO2 (color online)

Acknowledgments:

This work was supported by Region Nouvelle Aquitaine [project PICORE 2018-1R10112] which funded the postdoctoral position of J. Lelievre