Antoine Haddon
email: ahaddon@dim.uchile.cl

Héctor Ramírez
email: hramirez@dim.uchile.cl

Alain Rapaport
email: alain.rapaport@inrae.fr

Benoit Chachuat

A Note on "Optimal and Sub-Optimal Feedback Controls for Biogas Production"

Keywords: Optimal control, Chemostat model, Singular Arc, Sub-optimality, Infinite horizon Mathematics Subject Classification 49K15, 49N35, 49N90, 93B52

Introduction

In Ref. [START_REF] Haddon | Optimal and sub-optimal feedback controls for biogas production[END_REF], Proposition 3.1 deals with the convergence of the discounted reward (16), the associated value function (17) and optimal trajectories, as the discount factor goes to 0. The proof of the Γ -convergence of the discounted reward is incorrect since, in general, this reward is not monotone with respect to the discount factor δ.

The Correction

Proposition 3.1 can be revised as follows.

Proposition 2.1 For all ξ ∈ D and for all δ > 0, the suprema are attained,

V δ (ξ ) = max ζ(•) J δ (ζ (•)).
If the Γ -limit of J δ (•) exists as δ goes to 0,

J 0 (ζ (•)) := Γ -lim δ→0 J δ (ζ (•)),
then the maxima converge, pointwise in ξ , to the maximum of the limit,

V 0 (ξ ) := lim δ→0 V δ (ξ ) = max ζ(•) J 0 (ζ (•)).
(1)

Furthermore, if ζ δ (•) is an optimal trajectory, i.e. if V δ (ξ ) = J δ (ζ δ (•)), and if ζ δ (•) converges to ζ 0 (•) in S(ξ ), then ζ 0 (•) is an optimal trajectory for (1) and V 0 (ξ ) = J 0 (ζ 0 (•)) = lim δ→0 J δ (ζ δ (•)).
Proof To show that the suprema are attained, we show that the set of all forward trajectories of (3) of [START_REF] Haddon | Optimal and sub-optimal feedback controls for biogas production[END_REF], with initial condition ξ , is compact for the topology on W 1,1 0, ∞; R 2 , e -bt dt given in Definition 3.1 of [START_REF] Haddon | Optimal and sub-optimal feedback controls for biogas production[END_REF].

For each ξ ∈ D we set

F ξ (ζ ) := F P L(ξ ) (ζ ) ,
where P L(ξ ) is the projection on the convex set L(ξ ). Then, F ξ has linear growth, so that we can define

c = sup ζ ∈Dom(F ξ ) ||F ξ (ζ )|| ||ζ || + 1 , where ||F ξ (ζ )|| := sup η∈F ξ (ζ ) ||η||.
Note that F is upper semi-continuous and has compact non-empty convex images (such a map is known as a Marchaud map [START_REF] Aubin | Viability theory[END_REF]). With this, the set S(ξ ) is the set of absolutely continuous solutions of the differential inclusion

ζ (t) ∈ F ξ (ζ (t)), ζ (0) = ξ.
We can therefore use [2, Theorem 3.5.2] to establish that S(ξ ) is compact for the topology of W 1,1 0, ∞; R 2 , e -bt dt for b > c, thereby proving the existence of optimal trajectories in S(ξ ).

123

In addition, this allows us to show that the maxima converge to the maximum of the limit. Indeed, when the rewards Γ -converge, it is sufficient to show that there exists a countably compact set on which the suprema are attained for all δ [START_REF] Dal Maso | An Introduction to Γ -Convergence[END_REF]Theorem 7.4]. The set S(ξ ) is clearly independent of δ and countably compact, since it is compact. Finally, the convergence of optimal trajectories can be shown with [START_REF] Dal Maso | An Introduction to Γ -Convergence[END_REF]Corollary 7.20].

Conclusions

The originally published proof of Proposition 3.1 of [START_REF] Haddon | Optimal and sub-optimal feedback controls for biogas production[END_REF] was incorrect and we have revised here the result to obtain an accurate statement. However, we have not found reasonable assumptions that ensure the existence of the Γ -limit of J δ (•), when δ goes to 0, although it seems to be satisfied in our examples. We thus posit the Γ -convergence as a conjecture, that will be investigated in future research.

Acknowledgements

The authors would like to thank Frédéric Bonnans for pointing out a problem with the originally published proof. The first and second authors were supported by FONDECYT grant 1160567 and 1160204, and by Basal Program CMM-AFB 170001 from CONICYT-Chile. The first author was supported by a doctoral fellowship CONICYT-PFCHA/Doctorado Nacional/2017-21170249. The third author was supported by the LabEx NUMEV incorporated into the I-Site MUSE.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Bruno Leite

Thomas Croguennec

Amira Halabi

Esly Ferreira

Da Costa Junior

Esly Ferreira Da Costa Junior

Comparing different methods for estimating kinetic parameters of whey protein heat-induced denaturation in infant milk formulas

Keywords:

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

Introduction

Modeling, predicting, and controlling the extent of whey protein denaturation is a critical issue to optimize the quality of the dairy products and the performance of the processing unit. Most dairy products are heat processed. Heating takes place immediately after the raw milk is delivered to the dairy companies and/or is a step during the manufacture of dairy products. Aside improving the microbial safety and extending the shelf life of dairy products, heat treatments induce the denaturation of the whey proteins that have positive or negative impacts on the quality of the final products and impair the performance of the processing unit. Heat-induced whey protein denaturation is used to improve yoghurt firmness [START_REF] Sodini | The relative effect of milk base, starter, and process on yogurt texture: a review[END_REF]. However, in the manufacture of high protein yoghurts, the heating of the concentrated milk and consequently the extent of whey protein denaturation has to be reduced to achieve the optimal texture of the final gels. Otherwise, the gel is too strong, too pasty and sticky compared to standard yoghurts [START_REF] Purwanti | New directions towards structure formation and stability of protein-rich foods from globular proteins[END_REF]. Heating cheese milk usually improves cheese yield as more denatured whey proteins and more water are retained in the curd [START_REF] Kelly | Pre-treatment of cheese milk: principles and developments[END_REF]. Whey protein retention in the curd also improves the nutrient value and it modifies the functional properties of the cheeses [START_REF] Hinrichs | Incorporation of whey proteins in cheese[END_REF]. However, if the heating and the extent of whey protein denaturation is too high, the gel time of rennet-induced coagulation is delayed [START_REF] Vasbinder | Impaired rennetability of heated milk; study of enzymatic hydrolysis and gelation kinetics[END_REF], and the gel elasticity decreases. In the context of infant milk formula, the extent of whey proteins denaturation affect the physical stability of the liquid formulation during manufacture [START_REF] Buggy | Stabilising effect of α-lactalbumin on concentrated infant milk formula emulsions heat treated pre-or post-homogenisation[END_REF] and modify the nutritional quality of the proteins [START_REF] Peram | In vitro gastric digestion of heatinduced aggregates of β-lactoglobulin[END_REF]O' Loughlin et al., 2012). Moreover, the denaturation of the whey proteins is mainly responsible for the fouling of the dairy equipment's [START_REF] Burton | Deposit of whole milk in treatment plantsa review and discussion[END_REF]Fickak et al., 2011Fickak et al., 2011;[START_REF] Truong | Role of β-lactoglobulin in the fouling of stainless steel surfaces by heated milk[END_REF].

β-LG is often described as the protein responsible for the abovementioned effects in dairy products. The heat denaturation of β-LG proceeds along various stages and it is considered irreversible above around 75 • C [START_REF] Wit De | Thermal behaviour of bovine b-lactoglobulin at temperatures up to 150 C. a review[END_REF][START_REF] Brodkorb | Heat-induced denaturation, aggregation and gelation of whey proteins[END_REF]. Lactoferrin (LF) is one of the most heat-labile whey proteins with a denaturation temperature around 65 • C [START_REF] Brisson | Heat-induced aggregation of bovine lactoferrin at neutral pH: Effect of iron saturation[END_REF] and it exhibits a wide spectrum of biological properties [START_REF] Baker | A structural framework for understanding the multifunctional character of lactoferrin[END_REF][START_REF] Tomita | Twenty-five years of research on bovine lactoferrin applications[END_REF]. Thus, this study was focused in the denaturation of β-LG and LF.

The relationship between protein denaturation, heating time and temperature is described by a nonlinear mathematical expression. In such case, linear regression techniques must be extended, which introduces considerable complexity [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF]. [START_REF] Fogler | Heat-induced denaturation and aggregation of b-Lactoglobulin: kinetics of formation of hydrophobic and disulphidelinked aggregates[END_REF] suggests rearranging non-linear equations to obtain linear relationships between measured variables and perform a linear regression, then using parameters estimated by the linear regression as initial guesses to perform a nonlinear regression. Other possible approaches are performing a weighted least-squares regression [START_REF] Shalizi | Lecture 24-25: Weighted and Generalized Least Squares[END_REF][START_REF] Dumouchel | Using sample survey weights in multiple regression analyses of stratified samples[END_REF] or a general least-squares regression [START_REF] Mardia | Multivariate Analysis[END_REF].

In some situations, there is no generic formula of one parameter estimation method that overcomes all others. [START_REF] Kohberger | Statistical evaluation of the direct linear plot method for estimation of enzyme kinetic parameters[END_REF] has examined three different methods for evaluating kinetic parameter estimation of an enzymatic reaction described by the Michaelis-Menten equation: a direct linear plot, nonlinear least-squares and weighted least squares. His conclusions reinforce the importance of previous knowledge of the error structure to perform precise estimations, and that the direct linear plot cannot be used as a substitute for the least-squares method when precision is required.

Loveday (2016) has compared estimation of kinetic parameters of β-lactoglobulin (β-LG) denaturation performed in one-step method (estimating the energy of activation, kinetic constant and reaction order simultaneously) to estimation performed in two-step method (estimating kinetic constant for different temperatures in the first step, then estimating the activation energy by the Arrhenius equation). One important issue when estimating kinetic parameters of protein heat-denaturation, is that the kinetic rate constant k n is dependent on the temperature. Thus, to be estimate k n the heating temperature must be stationary. [START_REF] Loveday | β-Lactoglobulin heat denaturation: a critical assessment of kinetic modelling[END_REF] points the uncertainty as to when stationary temperature is reached, because it takes a finite time that is often unknown, and some protein denaturation may occur during temperature equilibration. Although the system quickly reaches stationary temperature, the progress of protein denaturation while temperature increases can't be ignored, so [START_REF] Van Loey | Inactivation of enzymes[END_REF] suggest that to heat a sample then stop heating as soon as stationary temperature is reached and use this time as the time course 'zero point' whose concentration of native proteins is designated C 0 . Another approach is to apply mathematical methods, to convert experimental heating times, in the equivalent heating times that, at the stationary temperature T of each series, would lead to the same denaturation level (Halabi et al., 2020a). When performing estimations redefining the "zero point", as suggested by [START_REF] Van Loey | Inactivation of enzymes[END_REF], there was no evidence of significant difference in the results. Considering the rapid heating-up time of the protein solutions, some authors consider all experimental points in the fit of the kinetic parameters [START_REF] Kehoe | The influence of bovine serum albumin on β-lactoglobulin denaturation, aggregation and gelation[END_REF][START_REF] Croguennec | Heat-induced denaturation/aggregation of b-lactoglobulin A and B: kinetics of the first intermediates formed[END_REF].

The objective of the current study was to evaluate performance of different parameter estimation methods in estimating kinetic parameters of the heat-induced denaturation of individual whey proteins (β-LG and LF) in different infant milk formula. The methods under study were a two-step method, nonlinear least-squares (NLS), one-step linearized, and weighted least-squares (WLS). In case of the latter, applying two different weight functions.

Materials and methods

Preparation of the IMFs

Three different IMFs were used in this study. The IMFs were formulated in agreement with the European regulation (EU, 2016/127) for the bovine-milk-based IMFs regarding the contents of protein (1.3% w/w), lactose (5.7% w/w) and the principal minerals (Calcium: 42, sodium: 22 potassium: 62, iron: 0.42, chloride: 68, inorganic phosphate: 21 mg/100g of IMF). Skimmed milk (SMP), whey protein isolate (WPI, Prolacta®95, Lactalis Ingredients, Bourgbarré, France), LF (Prodiet Lactoferrin®, Ingredia Dairy Experts, Arras, France) and α-lactalbumin were used for IMF formulation. Lactose was first dissolved in Milli-Q water and then protein powders were added. After protein dispersion, the solutions were supplemented by minerals, and stirred during 10 min at room temperature. The IMFs were adjusted at pH 6.8 with 1 M KOH and stored overnight at 4 • C to ensure the complete rehydration and powder equilibration. The only difference between the IMFs is that standard IMF contained β-LG (0.51 ± 0.01%) and no added LF. LF1 and LF2 are IMFs that contained the same amount of added LF (0.16 ± 0.01%) but different amount of β-LG (0.33 ± 0.01% and 0.06%, respectively). α-LA powder was used to vary β-LG content in LF1 and LF2.

Thermal treatment and denaturation kinetics

The IMFs (800 μL) in glass capillary tubes (8 mm inner diameter and 40 mm length, Waters, USA), were heated in a thermally controlled water bath (Fisherbrand™ Isotemp™ water bath, Thermo Fisher Scientific, Newington, USA) set at different temperatures ranging from 67.3 • C to 79.6 • C by intervals of about 2.5 • C. Model IMFs were not heated above 80 • C because the mechanism of β-LG denaturation change leading to a change of bβ-LG activation energy [START_REF] Hillier | Thermal denaturation of α-lactalbumin and β-lactoglobulin in cheese whey: effect of total solids concentration and pH[END_REF][START_REF] Petit | Influence of calcium on β-lactoglobulin denaturation kinetics: implications on unfolding and aggregation reactions[END_REF]. For each temperature, the heating ramp (i.e. the time to reach the target temperature) was controlled with a YC-747UD Data Logger Thermometer inserted in a reference tube. The heating ramp was converted into a heating time at target temperature (Halabi et al., 2020a). At various time intervals, tubes were removed and instantaneously chilled on ice bath to stop whey protein denaturation/aggregation reactions. The unheated IMFs were used as reference. Three independent experiments were performed for the heat treatment at 75 • C, considering the standard deviation equivalents at all heating temperatures investigated.

Determination of the contents of β-LG and LF

Native whey proteins were separated from the caseins and denatured/aggregated whey proteins, which precipitate after IMF pH adjustment at 4.6, by centrifugation at 14,000g for 20 min at 25 • C (Eppendorf Microcentrifuge MicroStar 17R, VWR, Leuven, Belgium). The supernatants of the IMFs were recovered for native β-LG and native LF separation by reverse phase-HPLC (Dionex UltiMate 3000 HPLC System, Thermo Fisher, Dreieich, Germany) using a PLRP-S 300 Å, 8 μm, 150 × 3 mm column (Agilent Technologies, UK). Protein were eluted from the column at a flow rate of 0.2 mL/min using a gradient of acetonitrile obtained by mixing the mobile phase A (0.106% (v/v) of trifluoroacetic acid in Milli-Q water) and the mobile phase B (0.1% (v/v) of trifluoroacetic acid in acetonitrile). Native whey proteins in unheated (C 0 ) and heated samples (C t ) were quantified using a calibration curve established by direct injection of β-LG and LF standards.

Kinetics of denaturation

The kinetics of milk protein denaturation is usually described by equation (1), which integrated for n ∕= 1 and rearranged results in equation (2). [START_REF] Oldfield | Kinetics of Denaturation and Aggregation of Whey Proteins in Skim Milk Heated in an Ultra-high Temperature (UHT) Pilot Plant[END_REF][START_REF] Croguennec | Caseinomacropeptide modifies the heat-induced denaturation-aggregation process of β-lactoglobulin[END_REF]Halabi et al., 2020a).

- dC dt = k n ⋅C n t (1) C t C 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 + (n -1)⋅k n ⋅C (n-1) 0 ⋅t (1-n) √ (2)
In which C t is the native protein concentration in time t, C 0 is the native protein concentration before heating, k n is the reaction rate constant at a defined temperature T n , and n is the reaction order according to the protein concentration. Units might be arbitrary, but, in the current study, concentration is given in (g.L -1 ), time in min, k n in (g (1-n) .L (n-1) .min -1 ), n is dimensionless.

The correlation between the reaction rate constant k n and given temperature T is given by the Arrhenius relationship given by equation (3), as described by [START_REF] Oldfield | Kinetics of Denaturation and Aggregation of Whey Proteins in Skim Milk Heated in an Ultra-high Temperature (UHT) Pilot Plant[END_REF].

k n = k n0 ⋅e -Ea R ⋅ 1 T (3)
In which, Ea corresponds to the activation energy, k n0 to a preexponential constant, R to the ideal gas constant and T to the temperature. Once again, defining units is arbitrary but in the current study Ea is in (J.mol -1 ), k n0 is in (g (1-n) .L (n-1) .min -1 ), R is in (J.mol -1 .K -1 ), and T is in K.

Substituting equation (3) in equation ( 2), results in equation ( 4) that might be rearranged to equation ( 5) by substituting k n0 by k n1 of an arbitrary temperature T 1 from equation (3):

C t C 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ 1 + (n -1)⋅k n0 ⋅e -Ea R ⋅ 1 T ⋅C (n-1) 0 ⋅t (1-n) √ (4) C t C 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 + (n -1)⋅k n1 ⋅e Ea R ⋅ ( 1 T 1 -1 T ) ⋅C (n-1) 0 ⋅t (1-n) √ (5)
2.5. Methods for parameter estimation

Two-step method

The two-step method consists of the estimation of the kinetic parameters in two steps: one for estimating kinetic constants for each temperature and another for estimating the activation energy and global kinetic constant [START_REF] Dannenberg | Reaction kinetics of the denaturation of whey protein in milk[END_REF][START_REF] Petit | Influence of calcium on β-lactoglobulin denaturation kinetics: implications on unfolding and aggregation reactions[END_REF]Halabi et al., 2020a).

Sometimes, the reaction order is fixed arbitrary to 1, 1.5 or 2 in reference to previous studies [START_REF] Dannenberg | Reaction kinetics of the denaturation of whey protein in milk[END_REF]Galani and Owusu Apenten, 1999) to determine the rate constant. In the two-step method used in the current study, the reaction order n for the denaturation of LF and β-LG at each sample was defined plotting (C t /C 0 ) (1-n) versus time for each temperature and retaining the n value in the range 1-2 with the best R 2 .

Using the mean of retained values of n, it was determined the kinetic constant k n for each temperature, by the linear regression slope of (C t / C 0 ) (1-n) versus time, according to equation ( 6).

( C t C 0 ) (1-n) = 1 + (n -1) ⋅ k n ⋅ C (n-1) 0 ⋅t (6)
In the following step, the activation energy Ea and logarithm of the kinetic pre-exponential constant ln(k n0 ) were determined respectively by the slope and intercept of the linear regression for the logarithm of the kinetic constant k n (estimated at each temperature) versus the inverse of the heating temperature, described by equation ( 7)

ln(k n ) = ln(k n0 ) - Ea R ⋅ 1 T n (7)
It is remarkable that a linearization is performed in both steps of this method to correlate to dependent variable to the independent variable. Linearization is a technique described by [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF] that allows to find algebraic solutions to estimate parameters of nonlinear models. It is particularly useful because, since there is an algebraic solution, the method does not demand an initial guess. However, it has the disadvantage of modifying the profile of residuals distribution and so might lead to bad estimations.

Nonlinear least-squares

The nonlinear least-squares (NLS) method was employed to estimate in one step the kinetic parameters activation energy Ea, kinetic constant k n1 (defined at an arbitrary T 1 = 345 K), and reaction order n, for each protein in each IMF, fitting experimental data do equation ( 5).

Comparing to the two-step method we must highlight some differences. The NLS method estimates simultaneously (in a single step) the parameters Ea, n and k n1 , as performed by [START_REF] Oldfield | Kinetics of Denaturation and Aggregation of Whey Proteins in Skim Milk Heated in an Ultra-high Temperature (UHT) Pilot Plant[END_REF]. Moreover, in the NLS method the dependent variable is nonlinearly related to the independent variables, while each step in the two-step method is based on a linear relationship.

According to [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF], a nonlinear regression model can be written as by equation ( 8).

y i = f (Xi, θ) + z (8)
In which y i corresponds to the experimental observation of index i, f to the expectation function of the model, X i is a vector associated to the independent variables of index i, θ is a vector of parameters, and z corresponds to residuals of mean 0 and homogeneous variance σ 2 .

The NLS regression was performed in Python 3.7 using the curve_fit function with least-squares method from scipy library [START_REF] Jones | SciPy: open source scientific tools for Python[END_REF].

As described by [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF], [START_REF] Fogler | Heat-induced denaturation and aggregation of b-Lactoglobulin: kinetics of formation of hydrophobic and disulphidelinked aggregates[END_REF], and by [START_REF] Edgar | Effect of whey protein concentration on the fouling and cleaning of a heat transfer surface[END_REF], a nonlinear least-squares regression consists in, starting from an initial guess, applying optimization algorithms to estimate parameters that minimize the sum of squared residuals from a given model, represented in equation ( 9).

SSe = ∑ no i=1 (y i -ŷi ) 2 (9)
In which, SSe is the sum of squared errors (residuals), i is the index for each observation or estimation, no is the number of experimental points, y i is the experimental observation of index i, and ŷi is the y value calculated by the model for the observation of index i, using estimated parameters.

Another possible approach based on the same assumptions would be minimizing the sum of absolute errors (least absolute deviation), but it is a method of difficult convergence as described by [START_REF] Li | A Maximum Likelihood Approach to Least Absolute Deviation Regression[END_REF].

The parameters errors were calculated according to [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF] from the standard deviation of each parameter, with assumptions of linearity of the model in the neighborhood of the optimal solution, and residuals of mean 0 and constant variance. The standard deviations were obtained as described by [START_REF] Silveira | Modeling of basil essential oil extraction with variable diffusion coefficient[END_REF] and by [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF].

The errors (residuals) variance σ 2 was calculated by equation (10):

σ 2 = SSe no -p (10)
In which, p is the number of estimated parameters and no is the number of observation points (experimental points).

The covariance matrix of the parameters V(θ) (as returned by the curve_fit function), can be calculated by equation (11).

V(θ) = σ 2 ⋅ ( J T ⋅J ) -1 (11)
In which, J is a matrix of shape (no, p) with the partial derivatives of each ŷ of index i, to each parameter θ of index j, according to equation (12).

J ij = ∂ŷ i ∂θ j (12)
In which, i is in range from 1 to no, and j in range from 1 to p. The standard deviation for each parameter s(θ j ) is given by equation (13).

s ( θ j ) = ̅̅̅̅̅̅̅̅̅̅̅̅ V(θ) jj 2 √ (13)
At last, the error associated to each parameter was obtained by the product of s(θ j ) and t-value considering a valid t-distribution.

For a matter of comparison to the parameters estimated by the twostep method, it was calculated the natural logarithm of the preexponential constant k n0 , which was calculated substituting k n by k 1 and T by T 1 at equation (7).

To calculate errors for ln(k n0 ), the same methods applied to Ea, k n1 , and n, were performed. This time, in a model considering Ea, ln(k n0 ), and n, according to equation ( 14).

Ct C 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ 1 + (n -1)⋅e ln(kn0) ⋅e -Ea R⋅T ⋅C 0 (n-1) ⋅t (1-n) √ (14)
As expected, the estimates for Ea and n and their respective confidence intervals remain the same, and the only additional information after the new estimation in the covariance matrix V(θ) is the value associated to ln(k n0 ), used to calculate its confidence interval.

One-step linearized

To distinguish effects of the data linearization from effects of performing estimation in one or two steps, it was developed a one-step method applying nonlinear optimization algorithms to estimate kinetic parameters Ea, k n1 , and n that would best fit the linearized kinetic equation ( 15), in which the dependent variable is linear in respect to time, but not to temperature. In this study this method is referred to as "one-step linearized method" or just "linearized method".

Both the two-step method and the one-step linearized method consider the dependent variable as being linear versus time, while they differ in estimating simultaneously the kinetic parameters Ea, k n0 , and n or estimating in two-steps. Hence, comparing results of the two methods it is possible to assess the difference in performing estimations in one or two steps. Moreover, comparing the one-step linearized method to the NLS method it is possible to assess effects of the linearization with respect to time (both estimations being performed in one step).

( C t C 0 ) (1-n) = 1 + (n -1) ⋅ k n1 ⋅ e Ea R ⋅ ( 1 T 1 -1 T ) ⋅ C (n-1) 0 ⋅t (15)
Although being linear according to time, the linearized method presented in this study was applied to a dataset considering different temperatures, and since the dependent variable has a nonlinear relationship with temperature, there isn't algebraic solution and so it isn't possible to apply standard linear regression methods.

Once again, the solution was obtained using Python 3.7, but this time applying the minimize function using nonlinear default optimization algorithms from scipy.optimize module [START_REF] Jones | SciPy: open source scientific tools for Python[END_REF]. In this method, the objective function was maximizing the R 2 (minimizing -R 2 ), described by equation ( 16), in which y is described by equation ( 17).

R 2 = 1 - ∑ no i=1 (y i -ŷi ) 2 ∑ no i=1 (y i -y i ) 2 (16) y = ( C t C 0 ) (1-n) ( 17
)
It is important to remark that since it is a nonlinear problem, R 2 can't be interpreted as it is in a multivariable linear regression method, where it means variance described by the model over total variance [START_REF] Montgomery | Applied Statistics and Probability for Engineers[END_REF].

The parameters errors were calculated the same way as using the nonlinear least-squares method, and so keeping assumptions of residuals with mean = 0 and constant variance σ 2 .

Weighted least-squares

The weighted least-squares (WLS) is an alternative to unweighted nonlinear least-squares method for reducing a possible bias caused by the sampling scheme [START_REF] Dumouchel | Using sample survey weights in multiple regression analyses of stratified samples[END_REF]. [START_REF] Shalizi | Lecture 24-25: Weighted and Generalized Least Squares[END_REF] suggests weighting residuals according to the inverse of the local variance in order to minimize the parameters variance. Besides, by attributing higher weights for residuals where there is less variance, the model is more likely to match observations where there is less noise and so where observations are more reliable.

The WLS method takes as objective function minimizing the weighted sum of squared errors, described by equation ( 18).

SSe w = ∑ no i=1 w i ⋅(y i -ŷi ) 2
(18) [START_REF] Kohberger | Statistical evaluation of the direct linear plot method for estimation of enzyme kinetic parameters[END_REF] uses the inverse of the predicted values squared as weights. [START_REF] Campos | Performance comparison of different mathematical models in the simulation of a solar desalination by humidification-dehumidification[END_REF] use the inverse of the observation values squared as weights. In this study, weights were proposed by equations ( 19) and (20).

w i = 1 y i ⋅ŷ i ( 19
)
w i = (t i + 0.25⋅t(T)) s(t(T)) ⋅ no no(T)⋅len(Temps) (20)
In which w i is the weight for observation of index i, t i is the time of observation of index i, t(T) is the average time of heating at temperature T, s(t(T)) is the standard deviation of the times of observations at temperature T, no is the total number of observations, no(T) is the number of observations at the given temperature, and len(Temps) is the number of unique temperatures evaluated in the experimental dataset.

The weights of equation ( 19) were chosen to avoid remarkably high relative residuals at low concentrations. The rate of denaturation is lower at lower concentrations, because of the first derivative (equation ( 1)). So, it is expected that the concentration is more stable at lower values and there's lower local variance. It is important to remark that the WLS regression with this weights function is of easier convergence than using the inverse of squared model values, described by [START_REF] Kohberger | Statistical evaluation of the direct linear plot method for estimation of enzyme kinetic parameters[END_REF].

Weights of equation ( 20) were chosen after analyzing the samples distribution to avoid sampling bias of having more observations at initial times than at ending times and of having more observations at some temperatures.

To perform the WLS, experimental data was fitted to equation ( 5), applying least_squares function from scipy.optimize module [START_REF] Jones | SciPy: open source scientific tools for Python[END_REF], to minimize the output of a weighted squared residuals function (equation ( 18)).

The parameters variance matrix was calculated by modifying equation (21) for weighted linear regression [START_REF] Dumouchel | Using sample survey weights in multiple regression analyses of stratified samples[END_REF].

V(β w ) = (X t ⋅W⋅X) -1 ⋅ ( X t ⋅ W 2 ⋅ X ) ⋅(X t ⋅W⋅X) -1 (21)
Substituting X by J in analogy to [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF] to transform from linear to nonlinear regression, it results in equation ( 22).

V(θw) = (J t ⋅W⋅J) -1 ⋅ ( J t ⋅ W 2 ⋅ J ) ⋅(J t ⋅W⋅J) -1 (22)
In which, instead of using β w, it was used θ w as the vector of the estimated parameters. W corresponds to the weight matrix, and J to the matrix of partial derivatives of the model to each parameter.

As to the NLS method, parameters errors were calculated as a product of t-value and the standard deviation calculated for each parameter from the diagonal of the variance matrix.

Methods of evaluation

Statistical

The methods for parameter estimation Ea, n and ln(k n0 ) were evaluated for each protein at each sample by the sum of squared errors from equation ( 9), and by adjusted mean absolute percentage errors (MAPE), given by equation ( 23), described by [START_REF] Armstrong | Long-range Forecasting: from Crystal Ball to Computer[END_REF].

MAPE = 1 no ∑ no i=1 2⋅|y i -ŷi | (y i + ŷi ) (23)

Graphical

Each method was evaluated by graph analysis, plotting C/C 0 versus time for each heating temperature for each protein at each sample. By equation ( 5), it is possible to infer that the shape of the model curves is related to the estimation of n, while the general trend of the curve is related to the kinetic constant, and the accuracy of fitting different temperatures is related to Ea. Underestimations of n lead to model curvatures less accentuated than experimental points while overestimations of n lead to model curvatures more accentuated than experimental points, perhaps intercepting them. Bad estimations of the kinetic constant will probably affect the general trend of the curve. At last, bad estimations of Ea lead the curve not to fit well all temperatures, but also if the Arrhenius equation can't be extended to some temperature interval, it is expected to have a bad fit for estimations at temperatures beyond the limits.

According to [START_REF] Bates | Nonlinear Regression Analysis and its Applications[END_REF], before characterizing the precision of estimates using confidence intervals, one must check the residuals for signs of model inadequacy. [START_REF] Edgar | Effect of whey protein concentration on the fouling and cleaning of a heat transfer surface[END_REF] suggest plotting residuals versus the dependent variable and versus independent variables to check detection of outliers, trends in the residuals, abrupt shifts, and changes in error variance. According to [START_REF] Montgomery | Applied Statistics and Probability for Engineers[END_REF], plotting residuals versus the dependent and independent variables is also a way to check the accuracy of the model, verifying if modeling assumptions are valid, as mean = 0 and constant variance. [START_REF] Shalizi | Lecture 24-25: Weighted and Generalized Least Squares[END_REF] suggests plotting residuals and squared residuals in order to determine an equation to describe weights as the inverse of the local variance.

A graphical analysis was performed in this study by plotting residuals multiplied by the square root of the weight function versus time and versus C/C 0 , for each protein at each sample. Multiplying the residuals by the square root of the weight function represents how the local residuals impact the objective function of each method, and so to what interval each model tends to better adjust.

Results

Applying the methods described above it was possible to estimate the kinetic parameters of the denaturation of β-LG and LF in different IMFs (Table 1). The parameters Ea (J.mol -1 ), k 1 (g (1-n) .L (n-1) .min -1 ), and n, correspond to kinetic parameters from equation ( 5). The arbitrary temperature of k 1 was defined as 345K. The parameter k n0 (g (1-n) .L (n-1) . min -1 ), corresponds to the kinetic pre-exponential term from equation (4).

The procedures for parameter estimation were performed in the temperature range from 67.3 • C to 79.6 • C for LF, while for β-LG the data obtained at 67.3 • C at both samples was not considered by the model because of evidence of poor adjustments at this temperature. Poor adjustments for β-LG at 67.3 • C are consequence of the very slow denaturation of the protein and the possible interference of concomitant reactions such as protein lactosylation through Maillard reactions that prevented the precise quantification of the native protein at longer heating times. To distinguish the two different weight functions used in the WLS method they are referred as WLS1 and WLS2.

For all proteins at all samples, the NLS method produced the lowest values in terms of SSe, as it would be expected once it is the objective function of the method, with mean 0.16. It was followed by the WLS method using the second weight function, that resulted SSe with mean 0.18. This result was significantly better than using the first weight function, that produced SSe with mean 0.27. The linearized method Fig. 2. Experimental points and kinetic model curves for concentration ratio versus heating time using estimations obtained by different methods for β-LG at LF1 sample. performed in one step usually led to the worst estimations in terms of SSe with mean 0.40, except for LF at LF2, in which the two-step method resulted in higher values. The mean SSe considering all methods was 0.27. In terms of adjusted MAPE, results obtained by all methods were quite similar for each protein at each IMF sample. Exceptions were observed for LF at the LF1 sample, in which the NLS method led to significantly higher values, and for LF at the LF2 sample, in which the two-step method led to significantly higher values. Thus, the mean value of adjusted MAPE for the NLS was the highest (15.4%), followed by the two-step method (13.4%). The best results in terms of adjusted MAPE were obtained by WLS1 and WLS2 with mean values of 12.3% for both. The mean adjusted MAPE considering all methods was 13.3%.

Whatever the method, the residuals were the lowest for β-LG at the standard sample, while quite similar for others. It reinforces the reliability of experimental data for β-LG at the standard sample.

Besides analyzing the residuals numerically, it was performed a graphical evaluation. Figs. 1234show the results of different methods in modeling the denaturation of β-LG at the standard sample, β-LG at the LF1 sample, LF at the LF1 sample and LF at the LF2 sample, respectively.

Fig. 1 reinforces the reliability of results obtained in Table 1, once the model curves were precise in fitting experimental points. However, linearized, two-step and WLS1 had slightly worse results at 67.3 • C, in which it is possible to observe a bias of the model curve being above experimental points.

At Figs. 2 and3, it is possible to notice more significant differences between estimation methods, perhaps because of the magnitude of experimental errors. Differences in terms of n observed in Table 1 are reinforced, as the model curves differ in shape, specially obtained using the NLS method, in which there is a remarkable bias of assuming values strictly higher than experimental points for longer holding times (or lower concentrations). When analyzing the shape of the curves using parameters estimated by the linearized method from Figs. 2 and 3 it is possible to notice that the magnitude of residuals is significantly lower for lower concentrations. This can lead to poor fittings at initial holding times, what was specially noticed for lower temperatures. Compared to the two-step method, there was no evidence of improvement, despite what is suggested by [START_REF] Loveday | β-Lactoglobulin heat denaturation: a critical assessment of kinetic modelling[END_REF].

At Fig. 4 it is possible to notice a significantly worse performance of the two-step method compared to others including the linearized onestep method, reinforcing conclusions from analyzing residuals at Table 1. However, this was not a behavior observed at other series.

Comparing the model curves obtained by the two different weight functions WLS1 was more effective in describing the denaturation of β-LG at the standard (Fig. 1) and LF1 (Fig. 2) samples at higher temperatures while WLS2 was more effective in describing it at lower ones. Moreover, the WLS1 was more effective in describing the denaturation of LF at the LF1 sample (Fig. 3) while WLS2 was more effective in describing of LF at LF2 sample (Fig. 4).

Weighted residuals for LF at LF1 sample were of significant relevance to compare methods performances and validate model's assumptions. They can be analyzed by Fig. 5. It is evident that the residuals variance is not constant neither in respect of time nor concentration for any method. It is a remarkable point specially about the NLS method because residuals constant variance is one of the method assumptions. The assumption of residuals with mean 0 also is not valid, once NLS led to residuals not equally distributed between positive and negative values, with bias of negative values at low concentrations (or long holding times) for the NLS method. This bias is reduced by both weight functions, specially the first.

Discussion

Methods performances

By analyzing SSe, adjusted MAPE (Table 1) and by graphical analysis, it is observed that the two-step method resulted in accurate estimations for kinetic parameters. While other methods need accurate initial estimations in order to converge, both steps from the two-step method have analytical solution, so it can be a powerful strategy to perform initial estimations.

Estimations using the two-step method were especially accurate for LF at LF1 sample, while there was lack of accuracy for LF at LF2 sample, as its possible to notice by analyzing Fig. 3 and adjusted MAPE (Table 1), For both linearized and two-step methods the most significant distancing from experimental points was noticed at lower temperatures, especially for LF.

Compared to the two-step method, there was no evidence of consistent improvement performing estimations by the one-step linearized method, despite what is suggested by [START_REF] Loveday | β-Lactoglobulin heat denaturation: a critical assessment of kinetic modelling[END_REF]. However, there was a consistent improvement in terms of residuals when performing the NLS method. This might evidence that the data linearization versus time can have a more negative impact on the accuracy of estimations than performing in two steps instead of one.

Despite the positive effect of minimizing SSe using the NLS method, there was evidence of poor estimations in terms of n and of biased residuals distribution proving the assumptions of residuals' constant variance and constant mean = 0 not to be valid. This might be consequence of sampling and/or of variance bias. However, the WLS method proved to be a good alternative to improve estimations reducing model curves biases observed when using the NLS method. By analyzing residuals (Fig. 5), we can conclude that both weight functions were effective in producing more homeostatic weighted residuals compared to unweighted, what can be a source of reducing bias in parameter estimations.

Comparing the two weight functions in terms of SSe and adjusted MAPE, the second produced more consistent results. Besides, as observed in Table 1, the confidence intervals for parameters estimated using WLS1 were much broader than when estimated using WLS, what evidences that the first weight function is of harder convergence to optimal solution than the second one.

Denaturation of β-LG and LF

Whatever the method and the proteins under study, the reaction order is fractional meaning that the kinetics of denaturation of β-LG and LF in infant milk formulae are complex kinetics. Except for LF in LF1 infant milk formula using NLS method (reaction order of 2.2 ± 0.17), the reaction order ranged between 1.5 and 2. The mechanism of denaturation of the whey proteins is described in two steps, an unfolding step consisting on protein conformational changes leading to reactive molecular species and an aggregation step in which the reactive species participate to several aggregation reactions [START_REF] Roefs | A model for the denaturation and aggregation of β-lactoglobulin[END_REF]. The former reaction is considered of first-order while the latter reactions are considered of second-order. Due to the intricacy of these two steps, unfolding and aggregation, during whey protein denaturation, the reaction order is expected to be in the range 1-2 in agreement with reported values. The reaction orders obtained in this study for β-LG denaturation are slightly higher than the reaction order determined during skimmed milk heating as summarized in Table 1 of the publication of [START_REF] Loveday | β-Lactoglobulin heat denaturation: a critical assessment of kinetic modelling[END_REF]. This difference could be explained by the higher whey protein/casein ratio in infant milk formulae compared to skimmed milk. In the absence of caseins, the reported reaction orders were mainly >1.5 [START_REF] Loveday | β-Lactoglobulin heat denaturation: a critical assessment of kinetic modelling[END_REF]).

Although [START_REF] Liu | Kinetic modelling of the heat stability of bovine lactoferrin in raw whole milk[END_REF] assumed reaction order of n = 1 to predict the kinetic parameters for the denaturation of lactoferrin in bovine milk. The most reliable predictions in this study are close to n = 1.6 for all samples. To our knowledge, no other study investigates the reaction order of the denaturation of LF in infant formulae.

For β-LG, it was observed that the highest denaturation rate was at the LF1 sample, which might evidence that the presence of LF increases the denaturation rate of β-LG. Moreover, in the LF1 sample, the denaturation of β-LG was more sensitive to a variation in temperature as evidenced by the higher Ea. In the temperature range investigated, Ea of β-LG denaturation in skimmed milk is between about 260 and 300 kJ mol -1 [START_REF] Dannenberg | Reaction kinetics of the denaturation of whey protein in milk[END_REF][START_REF] Anema | Effect of milk concentration on the irreversible thermal denaturation and disulfide aggregation of β-lactoglobulin[END_REF]. Ea slightly increases when the total solid of skimmed milk was increased [START_REF] Anema | Effect of milk concentration on the irreversible thermal denaturation and disulfide aggregation of β-lactoglobulin[END_REF]. In the standard IMF, the Ea of β-LG denaturation was around 330-340 kJ mol -1 whereas the β-LG concentration was in the same range in infant milk formula and in skimmed milk. This higher Ea indicates that the denaturation of β-LG is retarded in infant milk formula compared to in skimmed milk. This difference could come from different whey proteins/casein ratio or the difference in the amounts of lactose and minerals between infant milk formulae and skimmed milk. For instance, calcium content was shown to increase the Ea of β-LG denaturation [START_REF] Petit | Influence of calcium on β-lactoglobulin denaturation kinetics: implications on unfolding and aggregation reactions[END_REF]. Ea > 300 kJ mol -1 were determined in some studies dealing with the heating of whey protein solutions or pure β-LG solutions [START_REF] Hoffmann | Aggregation of β-lactoglobulin investigated by in situ light scattering[END_REF][START_REF] Le Bon | Kinetics of aggregation and gelation of globular proteins after heat-induced denaturation[END_REF][START_REF] Tolkach | Reaction kinetic pathway of reversible and irreversible thermal denaturation of β-lactoglobulin[END_REF][START_REF] Petit | Influence of calcium on β-lactoglobulin denaturation kinetics: implications on unfolding and aggregation reactions[END_REF]. The Ea of β-LG denaturation is even higher in the presence of LF (around 400 kJ mol -1 ) whereas β-LG denaturation rate is faster than in the absence of LF. At first sight, it appears contradictory but the presence of LF in infant milk formulae could have different impacts on the IMF components. First, it was observed that in the presence of LF, the casein micelle partially disintegrated (Anema and de Kruifs, 2012;Halabi et al., 2020b) and the presence of free αSand β-caseins is known to interact with partially unfolded β-LG on heating [START_REF] Morgan | Casein proteins as molecular chaperones[END_REF]. Secondly, LF exhibits a lower temperature of denaturation than β-LG [START_REF] Bengoechea | Formation of protein nanoparticles by controlled heat treatment of lactoferrin: factors affecting particle characteristics[END_REF]Fig. 5. Experimental weighted residuals versus C/C 0 and versus time for kinetic model using estimations obtained by NLS (black), WLS1 (magenta), and WLS2 (blue) for LF at LF1 sample. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) [START_REF] Brodkorb | Heat-induced denaturation, aggregation and gelation of whey proteins[END_REF] and a lower Ea than β-LG. Consequently, LF is readily unfolded on heating and it exposed many sites for aggregation reactions with β-LG.

For LF, it was observed that the denaturation rate was slightly lower in the LF2 sample than in the LF1 sample. However, the difference was not as significant as for β-LG and it makes sense if we consider that LF denatured earlier than β-LG.

Conclusions

The assumption of residuals with constant variance in which the nonlinear least-squares method is based proved not to be valid to this dataset, because the variance of residuals was observed to be constant neither versus time nor versus concentration and mean of the residuals was biased not equally distributed between positive and negative values. The linearized method was less effective than NLS and WLS in making good predictions because of uncontrolled distortion of residuals due to the linearization, and, although being performed in one step, it did not result in more accurate estimations when compared to the two-step method. The WLS with the proposed weight functions has been the best method for estimating kinetic parameters, avoiding bias observed in the NLS estimations. Analyzing SSe and confidence intervals of parameters we can infer that the second weight function produced even more accurate results than the first. However, defining the weight function requires previous knowledge of sampling biases and residuals variance, that are particular to each dataset. This can make it harder to implement the WLS without previously performing another regression method. Moreover, to perform the WLS method (in nonlinear problems) it is necessary to provide initial guesses for the parameters, as for any nonlinear regression.

Probably a good approach when performing estimations of kinetic parameters in future studies could be to start by the two-step method, use obtained estimations as initial guesses in the NLS method, analyze residuals distribution, and then define weights and apply WLS.

Declaration of competing interest

The authors whose names are listed immediately below certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Author names: Bruno Leite, Thomas Croguennec, Amira Halabi, Esly Ferreira da Costa Junior.

  (α-LA, Agropur Inc, Appleton, USA) powders, lactose powder (Armor Proteines, Saint-Brice-en-Coglès, France) and Na 3 C 6 H 5 O 7 ⋅2H 2 O (Carlo Erba Reagents, Val-de-Reuil, France), Na 3 PO 4 ⋅12H 2 O (Merck, Darmstadt, Germany), KCl (Panreac, Barcelona, Spain), CaCl 2 ⋅2H 2 O (AnalaR, Leuven, Belgium) and FeSO 4 ⋅7H 2 O (Sigma-Aldrich, St-Louis, USA) salts

Fig. 1 .

 1 Fig. 1. Experimental points and kinetic model curves for concentration ratio versus time using estimations obtained by different methods for β-LG at the standard sample.

Fig. 3 .

 3 Fig. 3. Experimental points and kinetic model curves for concentration ratio versus heating time using estimations obtained by different methods for LF at LF1 sample.

Fig. 4 .

 4 Fig. 4. Experimental points and kinetic model curves for concentration ratio versus heating time using estimations obtained by different methods for LF at LF2 sample.

  Methodology, Software, Formal analysis, Writingoriginal draft, Writing -review & editing, Visualization. Thomas Croguennec: Validation, Writing -original draft, Resources, Project administration. Amira Halabi: Validation, Investigation, Data curation, Resources. Esly Ferreira da Costa Junior: Methodology, Formal analysis, Supervision, Project administration.

Table 1

 1 Kinetic parameters and statistical indicators for estimations performed by different methods.

	Protein	IMF	Method	Ea (x10 3 )	n	k1 (x10 2 )	ln(kn0)	SSe (x10 2 )	MAPE
	β-LG	Std	NLS	326±7	1.81±0.06	0.71±0.05	108±5	4	5.5%
			linearized	346±11	1.61±0.1	0.76±0.09	116±7	11	6.1%
			WLS1	339±15	1.67±0.15	0.77±0.11	113±13	5.9	5.1%
			WLS2	331±8	1.73±0.08	0.75±0.07	110±6	4.4	4.9%
			two-step	337	1.6	0.79	112	8.9	5.8%
		LF1	NLS	390±18	2.03±0.14	2.56±0.29	132±12	11.7	14.6%
			linearized	412±28	1.5±0.21	2.63±0.45	140±19	36.4	14.1%
			WLS1	408±57	1.56±0.43	2.82±0.55	138±49	24.6	12.8%
			WLS2	396±24	1.8±0.19	2.82±0.36	134±16	13.5	12.2%
			two-step	416	1.6	2.79	141	22.5	12.9%
	LF	LF1	NLS	179±19	2.2±0.17	11.78±0.98	60±13	32.6	27.0%
			linearized	189±25	1.53±0.24	7.29±0.78	63±17	91.7	18.3%
			WLS1	180±99	1.6±0.66	8.76±4.63	60±83	59.5	16.3%
			WLS2	177±23	1.9±0.21	10.99±0.99	59±16	36.4	18.4%
			two-step	195	1.65	8.7	65	57.5	16.8%
		LF2	NLS	220±11	1.6±0.12	6.9±0.38	74±8	17	14.3%
			linearized	197±12	1.6±0.13	7.48±0.46	66±8	20.2	13.6%
			WLS1	200±55	1.48±0.7	6.73±3.88	67±68	19.9	14.9%
			WLS2	215±14	1.57±0.17	6.68±0.47	72±9	17.6	13.8%
			two-step	212	1.77	9.2	71	36	18.2%

Acknowledgements

The authors thank Universidade Federal de Minas Gerais and CNPq for granting the scientific initiation fellowship to Bruno Leite. The authors also thank The National Institute for Agricultural and Environmental Research (INRAE) -Science and Technology of Milk and Eggs and Britanny Region.


S Bhattacharyya

M Rejmund

A Navin
email: navin@ganil.fr

E Caurier

F Nowacki

A Poves

D. O' R Chapman

M Donnell

A Gélin

Hodsdon

D O'donnell

M Gelin

A Hodsdon

X Liang

W Mittig

G Mukherjee

F Rejmund

M Rousseau

P Roussel-Chomaz

K-M Spohr

Ch Theisen

Structure of neutron-rich Ar isotopes beyond N=28

Keywords: numbers: 21.60.Cs, 23.20.Lv, 27.40.+z, 25.70.Lm

  L'archive ouverte pluridisciplinaire

The availability of rare isotopic beams has opened up new avenues to study the evolution of nuclear structure as a function of neutron-proton asymmetry in addition to the earlier studied degrees of freedom, namely temperature and angular momentum. The appearance/disappearance of closed shell configurations at extreme values of isospin has attracted a lot of attention. The quest towards the understanding of how and if atomic nuclei near closed shells develop collective excitations/deformation with the available limited singleparticle degrees of freedom is also of current interest. Advances in experimental techniques, especially with inbeam γ -ray spectroscopy, now allow us to probe nuclei having twice the number of neutrons as protons. They provide an opportunity not only to search for "islands of inversion" but also for nuclei that are deformed near a spherical shell closures [1,2]. The deviation from the spherical symmetry (β) and the axial symmetry (γ) are often used to characterize nuclear shapes. Nuclei near "conventional" closed shells at extremes of isospin, thus, could offer a new domain for the study of the role of neutron-proton correlations on collective motion. This could be a fertile region to search for nuclei that violate axial symmetry (first discussed by Davydov and collaborators [3] half a century ago) complementing those made at high angular momentum.

In isotopes (isotones) spanning neutron (proton) numbers around the N (Z) = 20, 28 shell gaps, the couplings and/or the excitations of valence nucleons across a shell play a key role in determining the excitation spectra [4]. The nuclei in this region have been reported to show a wide variety of shapes, including evidence for shape coexistence [5,6] and superdeformation [7]. Shell model calculations have kept pace with experimental developments both to predict and understand the evolution of such structure away from stability. The sdf p effective interaction [8], among others, has been success-fully used to explain regions of deformation like that observed in 42 Si 28 [1,9], the neutron-rich Sulphur isotopes [5,6] and the evolution of the spin-orbit splitting around N = 28 [10]. The present work, in addition to increasing our understanding of cross-shell excitations around the double "magic" 48 Ca, is motivated by the search and understanding of the existence of nuclei having a triaxial ellipsoid shapes in this region [11].

Deep inelastic transfer reactions made at energies near the Coulomb barrier have shown to be an efficient way to study neutron-rich nuclei both with sufficient excitation and angular momentum [12,13]. This technique is complementary to spectroscopy using fast fragmentation beams. Nuclei that have no known spectroscopic information need to be uniquely identified, for example by a spectrometer [1,2,14], before any transitions can be assigned. In the present work a 1.310 GeV 238 U beam, with an intensity of ∼ 2 pnA, from the GANIL cyclotron facility was used on an isotopically enriched 48 Ca target (thickness 1 mg/cm 2 ) to produce neutron-rich nuclei around 48 Ca using deep inelastic transfer reactions. These measurements utilize the advantages of inverse kinematics (higher recoil energies and kinematic focusing of the residues). The large acceptance variable mode spectrometer (VAMOS), placed at 35 • with respect to the beam axis, was used to identify the target-like residues. The focal plane detection system consisted of two position sensitive drift chambers separated by a distance of one meter, one Secondary Electron Detector (SeD), placed between the two drift chambers, a segmented ionization chamber and a 21-element Si-wall [13,14]. The Z identification was obtained from the energy loss measured in the ionization chamber and the residual energy in the Si-wall. The time of flight of the recoils was measured by the SeD, relative to the radio frequency of the cyclotron and was used to obtain the velocity. The measured parameters at the focal plane along with the known mag- netic field were used to reconstruct the magnetic rigidity, the mass (A), mass/charge (A/Q), the scattering angle after the reaction and the path length on an event by event basis. The final selection of a particular nucleus was made from a two dimensional identification matrix of energy loss vs. A/Q for each charge state. Figure 1 shows a typical plot obtained for a single Si detector for a charge selection of Q = 18, where the unique identification of the various nuclei (Z ≥ 18) can be clearly seen.

The Doppler corrected γ-ray spectra, in coincidence with the fragments, were obtained using the EXOGAM array, consisting of 11 segmented clover detectors, placed at a distance of 11.4 cm around the target in a closed-pack configuration with partial Compton suppression. The angle between the γ-ray (from the segment of the relevant clover detector) and the reconstructed velocity vector of the emitting target-like nucleus (obtained from the measurements in the VAMOS spectrometer) were used to obtain the γ-ray energy in the rest frame. A typical energy resolution obtained is about 20 keV at 1 MeV. Figure 2 shows the Doppler corrected spectra for the 47,48 Ar isotopes. It should be pointed out that in such exotic nuclei the use of γ -γ coincidences to generate the level scheme is presently not possible. Further experimental details can be found in Ref. [14]. Gaudefroy et al. [10] reported the excitation energies and spectroscopic factors through the d( 46 Ar, 47 Ar)p reactions for states in 47 Ar. Different states can be expected to be populated in the present work due to the selectivity of the reaction mechanism. The γ-ray spectrum in coincidence with 47 Ar fragments shown in Fig. 2a is reported for the first time. The level scheme shown in Fig. 3a was constructed based on an intensity balance of the γ rays and also from a comparison with the level structure reported in Ref. [10]. The 1200(6) and the 1747 (5) 516(4), 1234(4) and 2190(16) keV in addition to the 1200 and 1747 keV transitions can also be seen in the figure.

The position of the 2190 keV level is tentative and is based on a relatively large decrease in the yield for a higher clover fold selection. The present experiment also reports for the first time the γ-ray spectroscopy of 48 Ar. This to the best of our knowledge represents the spectroscopy of the highest N/Z nucleus obtained using stable beams around the Coulomb barrier. Only a report on the lifetime [15] of the ground state exists in literature. The assignment of the γ-rays was possible due to an unambiguous identification of the coincident 48 Ar fragments in VAMOS. The corresponding Doppler corrected γ ray spectrum is shown in Fig. 2b. Three γ-rays are clearly visible in the spectra, namely 1041(9) keV, 1156 (15) keV and 1729 (13) keV. The corresponding level scheme shown in Fig. 3a was obtained using an intensity balance and using a fold selection of the clover detectors. The relatively large population of non-yrast states in addition to yrast states has been shown in deep inelastic reactions [13]. This was further verified in the Ar isotopes especially in 40 Ar where even the population of the 2 + 3 was observed rather than the population of a 6 + state. This aspect of the reaction mechanism was further used to constrain the proposed level scheme, especially the tentative assignment of the 2 + 2 state at 2197 keV. Fig. 3b shows the systematics for the energies of the first 2 + and 4 + states as a function of the neutron number for even-even isotopes of Ar. Data are taken from Ref. [16]. The energy and proposed assignment of the low-lying states in 48 Ar appear to consistent with the tentative assignments of the known 2 + and 4 + levels of the Ar isotopes.

The theoretical interpretation of the measured levels was carried out within an Interacting Shell Model (ISM) framework using the code NATHAN [17]. The interaction employed includes a recent improvement [18] of the sdf p interaction [8]. The valence space consists of the sd-shell for protons and the pf -shell for neutrons without any restrictions. We also discuss the signatures of triaxiality following the prescription of Davydov and Filippov [3] within this model. The comparison of the experimental level schemes with the results of the shell model calculations is shown in Fig. 3a. There is a very good agreement between experimental and calculated states. The low lying states of 47 Ar can be easily interpreted with the help of the ISM calculations. The ground state In 48 Ar, the data show a strong branch from the 2 + 2 state to the 2 + 1 relative to the direct decay to the 0 + ground state. The intensity limit for the latter was found to be < 30% of that observed for the 2 + 2 → 2 + 1 transition. The inset of Fig. 2b points towards a peak-like structure at around 1085 keV, which could be a candidate for a transition between the 3 + and 2 + 2 states. The ratio of the experimental energies E( 2 + 2 )/E(2 + 1 ), (top panel). In the case of 47 Ar due to the presence of an unpaired neutron, the angular momentum is represented as 2 * Jν -1.

E(4 +

1 )/E(2 + 1 ) are consistent with a signature of triaxiality as given in Ref. [3]. These experimental observations combined with the appearance of a low lying 3 + state in the ISM results (Fig. 3a) strengthen the hypothesis of 48 Ar being triaxial. Further and more importantly, the calculated B(E2 48 Ar is an order of magnitude smaller than the B(E2

)(3 + 1 → 2 + 1 ) in
)(3 + 1 → 2 + 2 ) and the B(E2)(2 + 2 → 0 + 1 )
is also an order of magnitude smaller than the B(E2

)(2 + 2 → 2 + 1 ). In the γ = 30 • case, the B(E2)(2 + 2 → 0 + 1 ) is strictly zero, and in the γ = 0 • (pro- late) or γ = 60 • (oblate) case it is the B(E2)(2 + 2 → 2 + 1
) which is strictly zero. These signatures of collectivity can be also explained by analyzing the properties of the calculated states in terms of (i) the individual angular momentum of neutrons (J ν ) and protons (J π ) coupled to form an observable final state J with |J ν ⊗ J π , J and (ii) the number of the neutron particle-hole (N ph ) excitations across the N = 28 shell gap (equivalent to the number of holes in the ν0f 7 2 orbital). Such an approach provides a greater physical insight into excitations in these nuclei. Fig. 4 shows the calculated shell model probabilities for the ground and excited states in 46-48 Ar. The bottom panel shows them as a function of the decomposition in terms of individual neutron (J ν ) and proton (J π ) angular momenta. The corresponding top panel shows them as function of the number of neutron(s) promoted out of the ν0f 7 2 orbital. As can be seen from the top panel of Fig. 4(a,e) the ground states of the even isotopes of Ar, are described as having only ∼ 50% of a N = 28 closed shell configuration. The corresponding N = 28 closed core configuration for 48 Ca is ∼ 90% [4]. References [4,18,19] discussed some aspects of the excitation in terms of N ph excitations and Ref. [10] indicated that the N = 28 is not fully closed in 46 Ar. The remaining part of the wavefunction of the 0 + states in 46,48 Ar is almost equally shared between one and two neutron excitations (N ph = 1, 2) across the N = 28 gap. The remarkably large contribution of the J ν = J π = 2 excitations (∼ 30%) in the 0 + ground states of 46,48 Ar can be seen from the corresponding bottom panels. Fig. 4(b,f) show that the 2 + 1 states are a mixture of both neutron (J ν = 2) and proton (J π = 2) states. The reason that the N ph = 0 in the 2 + 1 state of 46 Ar is smaller than the corresponding 48 Ar value is because, in order to make a 2 + neutron state in 46 Ar it is necessary to open the 0f 7 2 closed neutron orbit, which is not the case for 48 Ar. Fig. 4(c,d) quantifies the nature of the states in 47 Ar, discussed earlier, namely that the ground state can be represented as a coupling of the 0 + state of 46 Ar to a ν1p 3 2 neutron. The first 7

2

-state can be seen to have a dominant contribution arising from the coupling of 2 + state of 46 Ar and a ν1p 3 2 neutron. The shell model calculations discussed above thus predict that beyond N = 28, the neutron rich isotopes of Argon still have large probability of excitation of one or two neutrons across the N = 28 gap. The implications of the large fraction of N ph excitations seen in Fig. 4 point towards increased correlations which are expected to lead to deformation. The strong presence of N ph excitations from the ν0f 7 2 orbital in 46-48 Ar is a clear manifestation of the strong quadrupole interaction between protons in the sd shell and neutrons in the pf shell that results in energy gains comparable to the N = 28 gap. Some of these N ph excitations are very collective and may produce rapid transitions from spherical to deformed nuclear shapes. The increase of quadrupole collectivity in the heavy Argon and Sulfur isotopes around N = 28 is driven by the sd-shell valence protons outside the Z = 14 sub-shell closure. The relevant orbits, π1s 1 2 and π0d 3 2 , are known to become degenerate at N = 28. This degeneracy maximizes the quadrupole correlation energy of the configurations with open neutron orbits. The Sulfur isotopes with two protons in this pseudo-SU3 space are prolate in this scheme. The ISM calculations for 42 S show a well developed γ-band. However, following the criteria of Ref. [3], its γ is found to be small (5 • -10 • ). Calculations for 44 S show it to be fully axial whereas 46 S develops triaxiality again (β = 0.3, γ = 20 • ). This behavior pertains to variants of Elliott's SU3 coupling scheme; pseudo-SU3 for the protons and quasi-SU3 for the neutrons [20]. The Argon isotopes with two proton-holes in this pseudo-SU3 space are correspondingly oblate, and are found to follow closely the trend of the Sulfur isotopes but with a smaller deformation (β = 0.25). Using the transition probabilities from the ISM, 46 Ar is calculated to be oblate and axial whereas 48 Ar is found to be oblate and triaxial with (β = 0.25, γ = 40 • ). The appearance of triaxial shapes can be traced back to the structure of the Nilsson dia-grams in this region [21]. At N = 26 and N = 30 there are two quasi degenerate Nilsson levels available for the last pair of neutrons, both for oblate and prolate deformations. This can induce the K-mixing characteristic of triaxiality. On the contrary, at N = 28 deformed gaps develop both in the prolate and oblate sides, making 44 S and 46 Ar axial.

In summary, new limits of γ spectroscopy of the most neutron rich Ar isotopes at energies around the Coulomb barrier studied with deep inelastic transfer reactions have been presented. Shell model calculations which are in good agreement with the measured energies predict a large fraction of one and two neutron particle-hole excitations across the N = 28 shell gap, thus predicting a development of collectivity near a closed shell. The effect of these correlations is manifested by the presence of a triaxial intrinsic structure in 48 Ar. The measurements of transition probabilities and the occupancy of the neutron ν0f 7 2 orbital would further confirm the predictions of a non axially symmetric deformation in 48 Ar and also in the less accessible 46 S. This would increase our understanding of shape evolutions far from the valley of stability and lead us to more conclusive signatures for the presence of non axial deformations. Given the advances in beam intensities and associated detector development this appears to be within reach.

We would like to thank P. Van Isacker for many illuminating discussions. This work has been partly supported by the European Commission under contract 506065 (EURONS/INTAG), the IN2P3 (France)-CICyT (Spain) collaboration agreement and the FPA2007-66060 grant, MEC (Spain).

FIG. 2 :

 2 FIG. 2: Doppler corrected γ spectra for the Ar isotopes obtained in coincidence with corresponding residues detected. (a) 47 Ar. The inset shows an expanded view of the high energy region of 47 Ar spectrum. (b) 48 Ar. The inset shows the spectrum when at least 3 clover detectors have fired.

3 2 -and the excited 1 2 - 2 -, 7 2- 2 -state to the 5 2 - 2 -

 32227222 state correspond respectively to the coupling of a 1p3 2 and a 1p 1 2 neutron to the highly correlated 0 + ground state of 46 Ar. On the contrary, the doublet 5 results from the coupling of a 1p3 2 neutron to the 2 + state of 46 Ar. The calculated ratio of the intensities for the decay of 7 (E2) is 0.3 which compares well to the experimental ratio of 0.8(3).

2 -FIG. 4 :

 24 FIG. 4: (color online) Calculated shell model probabilities for 46-48 Ar as a function of the angular momentum of the neutrons (ν) and protons (π) space (bottom panel) and the number of neutron holes in the ν0f 7 2

Energy loss (MeV) 60 80 100 120 140 160 180 200 220 Q=18 48 46 44 42 40 50 52 54 56 58 Ar K Ca Sc Ti V Cr 14 16 18 20 22 Q

	2	2.2	2.4	2.6	2.8	3	3.2
				A/Q
	FIG. 1: (color online) Energy loss measured in the ion cham-
	ber as a function of the derived mass/charge (A/Q) for a
	selected charge state Q = 18. The inset shows the charge
	state distribution, obtained from the derived A and A/Q.