S A Berman

J Dubois

C Chandre

M Perin

T Uzer

Coherent buildup of high harmonic radiation: the classical perspective

We present a classical model for high harmonic generation during the propagation of an intense laser pulse through an atomic gas. Numerical simulations of the model show excellent quantitative agreement with the corresponding quantum model for the blueshift and intensity reduction of the propagating laser pulse over experimentally realistic propagation distances. We observe a significant extension of the high-harmonic cutoff due to propagation effects. Phase-space analysis of our classical model uncovers the mechanism behind this extension.

High harmonic generation (HHG) is the production of coherent high-frequency radiation observed during the ionization of gases by intense laser pulses. The highfrequency part of the spectrum typically consists of a plateau region, where the harmonics are produced with comparable intensity, followed by a cutoff region, which is often harnessed to generate attosecond pulses. To increase the flux of the highest harmonics for more intense attosecond pulses, experiments may be performed under conditions of increased driving laser intensity [1], gas length [2], or gas density [3][4][5]. Under such conditions, the driving laser field undergoes tremendous reshaping while propagating through the gas due to the radiation emitted by the ionizing atoms, leading for example to a blueshift and intensity reduction [6][7][8] throughout propagation. In this case, the high harmonic spectrum measured after propagation crucially depends on which high harmonic frequencies were produced at sufficient intensity all along the gas with just the right phase such that the radiation produced by the many atoms making up the gas adds up coherently, a collective effect referred to as phase-matching [9,10]. Thus, the self-consistent interaction between the ionizing atoms and the laser field plays a decisive role in shaping the high harmonic spectrum [3,4,10].

Ideally, a theoretical or numerical treatment of HHG must bridge the gap between the microscopic response of the atoms to the electromagnetic field and the macroscopic propagation of the field through a gas of billions of atoms. The most rigorous calculation would require the self-consistent solution of Maxwell's equations in three dimensions coupled to time-dependent Schrödinger equations (TDSEs) for the atoms [11,12]. Even today, the computational cost associated with this approach can be prohibitive, precluding the simulation of experimentally relevant sample lengths on the order of mms. Further, solutions of the TDSE provide limited intuition into the electron dynamics behind the single-atom response to the laser field. Alternatively, one can simplify the description of the atomic response, splitting it into a low-frequency part dominated by ionization [7] and a high-frequency part [13,14] comprised of the radiation emitted during repeated encounters between the ionized electrons with their parent ions [START_REF] Kulander | Super-intense laser-atom physics[END_REF][START_REF] Corkum | [END_REF]. The latter may be computed efficiently [9] using a semiclassical approach [17,18] under the assumption that the ionic core potential has a negligible effect on the ionized electron dynamics. This framework allows the simulation of experimental gas lengths, and the semiclassical description of the atomic response in terms of quantum orbits [17,18] facilitates the development of control strategies based on the trajectories of electrons after ionization [1,19,20]. However, these simplifications are inappropriate for the description of the harmonics near and below the atomic ionization threshold I p [21][22][23], which in certain situations can strongly influence the yield of higher harmonics [20,24]. Additionally, they leave out key elements of HHG in elliptically and circularly polarized fields [25,26]. Therefore, a theoretical formulation is needed which simultaneously accounts for the full complexity of the self-consistent atomfield interaction, includes the influence of the core potential on the ionized electrons, and allows for the understanding of the electron dynamics in phase space as the pulse propagates through the gas. Here, we propose a purely classical model which meets these requirements, and we demonstrate its validity and utility by comparing its behavior with a quantum model. In particular, we use it to identify the mechanism behind an intriguing phenomenon -the extension of the cutoff-observed in quantum simulations.

Our model describes the coupled evolution of the timedependent electric field E(τ, z) and the response of the atoms of the gas to the field throughout the laser pulse propagation. The evolution parameter in our model is z, the position along the laser propagation direction. For the atoms located at z, τ is the time relative to the arrival of the laser pulse to their position, i.e. τ = t -z/c. We may specify an arbitrary initial state for the atoms at τ = 0, which is uniform throughout the gas, and we may calculate the response of the atoms to the field in a classical or quantum manner. In what follows, we consider a linearly-polarized laser pulse interacting with single-active-electron (SAE) atoms. As an example, we present simulations of the model for a laser pulse with initial condition E(τ, z = 0) = E 0 cos(ω L τ ) propagating through 1 mm of a gas with density ρ = 5 × 10 17 cm -3 . The atomic initial condition, i.e. at τ = 0, is that of Refs. [27,28] -a fully ionized state with the electron described by a Gaussian wavepacket at rest at a distance of the quiver radius E 0 /ω 2 L from the ion, expressed in atomic units which are used unless stated otherwise. In Fig. 1, we compare the power spectra of the electric field at z = 0.37 mm, using a classical description of the atoms on the one hand and a quantum-mechanical one on the other. The spectra coincide at low frequencies, especially near the laser fundamental ω L , which is the frequency range where the dominant propagation effects on the electric field -the blueshift and intensity reductionare encoded. In the upper inset of Fig. 1, these effects are clearly seen on the time-dependent electric field after propagation through 1 mm of gas: the blueshift is evident from the shift of the field extrema to the left of their initial positions (indicated by dotted lines), while the intensity reduction is seen on the leading edge of the pulse, where the absolute value of the field amplitude relative to the incident amplitude is less than one. These effects are captured equally well by the classical and quantum atomic models: their respective time-dependent electric fields are indistinguishable in the upper inset of Fig. 1 and differ by less than 10 -2 E 0 for all considered times τ and propagation distances z. On the other hand, the classical description does not capture the high-harmonic plateau and cutoff radiation, which is present in the quantum model. There, we also observe a significant extension of the high-harmonic cutoff past the 3.17U p + I p cutoff law, where U p = E 2 0 /4ω 2 L is the ponderomotive energy. This is unexpected, given that the incident laser field is monochromatic and each atom only has a single active electron. Radiation at these anomalously high frequencies only begins to emerge clearly about 0.2 mm into the gas, shown in the lower inset of Fig. 1, indicating that it is truly a propagation effect. In the following, we will show that the purely classical model actually allows us to understand the mechanism of this anomalous high-harmonic radiation, despite its failing to capture the high-harmonic part of the spectrum on a quantitative level.

To begin, we consider the physics behind the classical model. We derived the model from first principles, starting from Maxwell's equations and the Lorentz force law for classical charged particles. The SAE atoms are assumed to have a static ionic core, and the electrons are assumed to be nonrelativistic and moving only in the direction transverse to the laser propagation direction z. Meanwhile, the electric field is assumed to lie in the polarization plane, with no longitudinal component, and its only spatial dependence is assumed to be the propagation coordinate z. Thus, our model neglects three dimensional effects, in particular the focusing of the laser beam and thus the Gouy phase shift. If desired, a z-dependent phase and intensity may be imposed externally to partially account for these effects [29], though we choose not to do this here in order to emphasize the self-consistent interaction between the radiation and the particles. Lastly, we assume that backward-propagating waves may be neglected, i.e. the field propagates solely in the positive z direction. Under these assumptions, the evolution equation for the electric field may be written in a frame moving at the speed of light c with the incident laser pulse as

∂E ∂z = 2πρ c v(τ, z) , (1)
where ρ is the atomic density and v(τ, z) is the ensemble-averaged dipole velocity at time τ of the atoms located at z driven by the field E(τ, z). For simplicity, we consider E to be linearly polarized along the xdirection as it is at z = 0 and take one-dimensional models for the atoms, but the two dimensional generalization is straightforward. For the classical model, v(τ, z) = vf (x, v, τ ; z)dxdv, where f is the probability distribution function to find an electron with position x relative to the ionic core and velocity v. At every z, f satisfies the Liouville equation corresponding to the single-atom Hamiltonian

∂f ∂τ = -v ∂f ∂x + ∂V ∂x + E(τ, z) ∂f ∂v , (2)
H(x, v, τ ; z) = v 2 2 + V (x) + E(τ, z)x. (3)
We use the soft-Coulomb potential V (x) = -(x 2 + a 2 ) -1/2 to describe the electron-ion interaction. For the quantum model, v(τ, z) is obtained from the solution of the TDSE with Hamiltonian (3) at every z [30]. Details on the numerical schemes employed to solve our model equations are provided in the Supplemental Material [START_REF]See Supplemental Material for an illustration of the model geometry, a description of the numerical schemes used to solve the model equations, and more[END_REF].

On a single-atom level, the classical model does not capture the high-harmonic plateau and cutoff because it lacks quantum interference effects [28]. While it may do a fair job for the low order harmonics originating from the nonlinear response of the bounded part of the wavepacket [START_REF] Bandarage | [END_REF]33], its spectrum is not accurate when the harmonics are driven by ionized, recolliding trajectories [27,28,34]. Evidently, this situation is unchanged when the classical single atoms are allowed to interact via the electromagnetic field, as shown in the spectrum of Fig. 1. At the same time, the quantitative agreement between the classical and the quantum models for low frequencies persists during propagation. Because these frequencies are the dominant constituents of the field (see Fig. 1), this suggests that the classical model provides a faithful representation of the true electron dynamics underlying the quantum model all along the propagation. This expectation is indeed borne out by the excellent correspondence between statistics of electron encounters with the core, or recollisions, obtained from the classical model and the dipole radiation spectra obtained from the quantum model, displayed in Fig. 2. According to the semiclassical model [START_REF] Kulander | Super-intense laser-atom physics[END_REF][START_REF] Corkum | [END_REF], an electron which is driven to the ionic core with a kinetic energy κU p may recombine into the atomic ground state with energy -I p and radiate its excess energy as a photon with frequency ω = κU p + I p . This implies that a specific frequency ω is only emitted when an electron trajectory with the appropriate kinetic energy κU p enters the core region [17], a behavior which can be revealed by performing a time-frequency analysis of the dipole acceleration -∂V /∂x(τ, z) from the quantum model [35,36]. In the lower panels of Fig. 2, the dipole acceleration timefrequency maps are displayed, showing which frequencies ω are generated at times τ . On the other hand, the top panels show the probability of an electron recollision with a given kinetic energy κU p at time τ , obtained from the classical model. A strong correspondence between these two figures is expected on the left panels at z = 0, when the classical and quantum atoms are driven by the exact same incident field, E(τ, z = 0) = E 0 cos(ω L τ ). However, after propagating to z = 0.37 mm, the laser fields in each model have been driven by different dipole velocities. Therefore, it is remarkable that the level of agreement between the two fields is so high that the dynamics of the electrons in both the classical and quantum cases continues to be very similar throughout propagation, evidenced by the continuing close correspondence between the two figures in the right panels.

The close correspondence of the two figures does not hold for all τ , however. For example, comparing the upper and lower left panels of Fig. 2, we see that at z = 0, there are recollisions which occur with κ ≈ 2 during the first laser cycle in the classical model without emission of the corresponding high harmonics in the quantum model [27]. The reason is the total depopulation of the ground state [36] at the beginning of the pulse, or more generally speaking the complete lack of an electron wavepacket at least transiently bounded to the ion. Because the highest harmonic frequencies are generated by the quantum interference between a bounded wavepacket and a recolliding wavepacket [28], a complete lack of a bounded wavepacket means very high harmonics cannot be emitted. On the other hand, by comparing the upper and lower right panels of Fig. 2, we see that by z = 0.37 mm, there is high-harmonic recollision-driven radiation emitted during the first laser cycle. This suggests the creation of a bound state earlier in the laser pulse as the propagation proceeds. Here, we exploit the main advantage of the classical model: its ability to confirm and analyze this scenario by visualizing the electron dynamics in phase space.

In Fig. 3, we show snapshots of the classical electron distribution function f (x, v, τ ; z) at particular times τ and propagation positions z. With the initial conditions we have chosen, the electron wavepacket always begins at rest far on the right side of the core, and the field is positive and at a maximum. Therefore, the electron is always initially driven to the left, towards the core. At z = 0, the electron wavepacket is predominantly on the left of the ion by time τ m when the laser field has reversed direction once and is again at an extremum. When the field subsequently goes to zero at time τ 0 , the wavepacket has almost completely vacated the core region. However, as propagation proceeds, the blueshift causes the laser field to reverse direction earlier in the pulse, and this causes the center of the wavepacket at time τ m to be displaced to the right. By z = 0.37 mm, the wavepacket is thus nearly centered over the ion, with the electron velocities distributed about zero. When part of the wavepacket arrives to the core with a low kinetic energy like this, it has a high probability of becoming trapped there [37], and indeed a trapped part of the wavepacket is clearly observed in the subsequent snapshot of the distribution function at τ 0 (lower right panel of Fig. 3). This confirms that a bound state is created earlier in the pulse after propagation through part of the gas. Furthermore, this explains the emergence of recollision-driven high-harmonic radiation for τ < 1 l.c. at longer propagation distances, despite this radiation being absent in this time interval at z = 0. In addition to providing an explanation of the electron trapping, the classical model also allows us to explain the anomalous high harmonic radiation in excess of the 3.17U p + I p cutoff law that we observed in the quantum calculation. At z = 0, the left panels of Fig. 2 indicate that both the recollision kinetic energy and dipole radi-ation cutoffs are in the expected place, with the cutoff frequency being ω c = 160ω L in this case. Thus, for small z, there is no coherent growth of frequencies ω > ω c , as shown in the lower inset of Fig. 1. However, as propagation proceeds, a family of recollisions appears, the first blue arc in the upper right panel of Fig. 2, that returns in the approximate range 0.75 < τ < 1.25 l.c. and has a maximum kinetic energy of about 4.1U p or 190ω L . A trajectory analysis reveals that this new family of recollisions also has a low kinetic energy encounter with the core near τ m , just like the early trapped trajectories. Rather than becoming trapped, however, these trajectories escape the core once more, returning with a range of kinetic energies up to about 4.1U p . This energy agrees well with a calculation of the maximum return kinetic energy of a free electron in the field E(τ, z) ionized at τ ∼ τ m , with the effect of the Coulomb potential treated as a perturbation [38]. Because this calculation predicts a cutoff of 3.17U p + I p for a monochromatic laser field, we conclude that the increase in kinetic energy above the usual high-harmonic cutoff is due to the departure of the field from a monochromatic wave throughout propagation. It is these recollisions' radiation which drives the evolution in energy of the macroscopic electric field in the frequency band of 175ω L to 200ω L starting at z = 0.2 mm, shown in the lower inset of Fig. 1. The circumstances permitting the anomalous high harmonic radiation are maintained over a substantial propagation length, leading first to the coherent growth of the energy in these frequencies, followed by the coherent absorption beginning at about z = 0.5 mm.

In summary, we have presented a purely classical model for HHG during the propagation of intense laser pulses through atomic gases. The model agrees quantitatively with a quantum model for the low frequency components of the laser field in the case of an initially monochromatic pulse propagating through a singly ionized gas, and the phase-space analysis permitted by the classical model explains the extension of the high-harmonic cutoff observed in the quantum simulation. Our model may be useful in the analysis of experimentally-observed excessively high harmonic radiation, though it is rare in gases [39], as well as the analysis and further refinement of trajectory-based semiclassical computational schemes [37,40] and control strategies for HHG [1,20]. Other possible applications of our model include the study of terahertz emission from field-ionized gases, where a classical description of the electron motion is also germane [41], and the study of filamentation, where first-principles descriptions of the atom-field interaction are increasingly sought [42].

We acknowledge François Mauger for helpful discussions. The project leading to this research has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement No 734557. S.A.B.

In this section, we describe the numerical methods we employ to solve our model equations. Equation (1) in the paper is a partial differential equation for the electric field E(τ, z), where z is the evolution parameter. The righthand side of the equation contains the ensemble-averaged dipole velocity v(τ, z) , which represents the response of the atoms to the field at each z. As illustrated in Fig. 1,v(τ,z) must be obtained for a given z, whether classically or quantum-mechanically, and this allows the electric field to be advanced one step ∆z in space to z + ∆z. This space-stepping is performed over the entire length of the gas, from z = 0 to z = 1 mm, using a third-order implicit multistep method.

To obtain v(τ, z) classically, we must first solve Eq. (2) in the paper, which is the Liouville equation

∂f ∂τ = -v ∂f ∂x + ∂V ∂x + E(τ, z) ∂f ∂v . (1)
This may also be written as

∂f ∂τ = -{f, H} (2)
where H(x, v, τ ; z) is Hamiltonian (3) in the paper and {•, •} is the Poisson bracket for canonically conjugate variables (x, v),

{F, G} = ∂F ∂x ∂G ∂v - ∂F ∂v ∂G ∂x . (3)
We solve Eq. (1) using the scheme developed in Ref. [1]: we express f as where δ is the Dirac delta function and the w α are constant weighting factors satisfying α w α = 1. The (x α (τ ; z), v α (τ ; z)) are classical trajectories obtained by numerically solving the equations of motion associated with H(x, v, τ ; z), using a third-order explicit symplectic scheme [2]. Having thus solved Eq. ( 1) for a given z, the ensemble-averaged dipole velocity is obtained as

f (x, v, τ ; z) = N α=1 w α δ(x -x α (τ ; z))δ(v -v α (τ ; z)),
v(τ, z) = vf (x, v, τ ; z)dxdv. ( 4
)
This approach is very similar to the classical trajectory Monte Carlo approach [3][4][5][6]. The main difference is that the initial conditions (x α (0; z), v α (0; z)) are chosen on a uniform grid, improving convergence with increasing N compared to the Monte Carlo case [7], and they are weighted such that the weights w α ∝ f 0 (x α (0; z), v α (0; z)) correspond to a given distribution function f 0 representing the initial state of the atom.

For the simulations presented in the paper, we used N ∼ 4 × 10 6 trajectories.

To obtain v(τ, z) quantum-mechanically, the following time-dependent Schrödinger equation must first be solved for the wavefunction ψ(x, τ ; z):

i ∂ψ ∂τ = - 1 2
∂ 2 ψ ∂x 2 + V (x) + E(τ, z)x ψ.

(5) use a split-operator method to numerically solve Eq. ( 5). The dipole velocity is then obtained from ψ as the integral of the dipole acceleration, v(τ, z) = v(0, z) (7) We take the window function to be W T (τ ) = cos 4 (πτ /T ) for the range -T /2 < τ < T /2 and zero elsewhere. This is done in Fig. 2 in the paper for the dipole acceleration -∂V /∂x(τ, z) computed from the quantum model, with T = 0.35 l.c. To obtain the time-dependent kinetic energy distributions of recolliding electrons from the classical model, we compute R(κ, τ ; z) = f (x, v, τ ; z)Θ(x c -|x|) ×Θ(κ c U p -|v 2 /2 -κU p |) dxdv, (8) where Θ is the Heaviside step function and x c = 5 a.u. and κ c = 0.014 are thresholds.

Because our distribution function f (x, v, τ ; z) is expressed as a sum over Dirac delta functions, it must be coarse-grained in order to make the plots in Fig. 3 of the paper. Thus, what is actually plotted in this figure is f (x g , v g ,τ ; z) = (∆x∆v) -1 f (x, v, τ ; z) ×Θ(∆x/2 -|x -x g |)Θ(∆v/2 -|v -v g |)dxdv (9) for x g and v g on a fine grid of spacing ∆x and ∆v in the x and v directions, respectively.

Figure 1 .

 1 Figure1. Filtered electric field spectra of a 3.5 laser cycle (l.c.) pulse at z = 0.37 mm of a singly ionized 1 mm gas with density ρ = 5 × 10 17 cm -3 . The incident field is E (τ, z = 0) = E0 cos(ωLτ ), with the field amplitude E0 and frequency ωL corresponding to an intensity I = 3.5 × 10 14 W • cm -2 and wavelength λ = 1.2 µm. The softening parameter a 2 = 2 is chosen to correspond to an ionization potential Ip = 0.5 a.u. The electron wavepacket is initialized at the quiver radius. For the blue curve the dipole velocity is computed classically, while for the orange curve it is computed quantummechanically. The dashed black lines correspond to 2 Up + Ip and 3.17 Up + Ip. Upper inset: The time-dependent electric fields from each model at z = 1 mm. The black dotted lines indicate the locations of the extrema of the initial field. Lower inset: Energy (in arbitrary units) in the electric field frequency band between 175ωL and 200ωL, as a function of propagation distance z.

Figure 2 .

 2 Figure 2. Electron dynamics along different points of the propagation. Left panels: z = 0 mm, at the gas entrance. Right panels: z = 0.37 mm. Top panels: Kinetic energy distribution of recolliding as a function time, with the density indicated by a logarithmic color scale, from the model. Bottom panels: maps of the dipole acceleration -∂V /∂x(τ, z) , computed a time window of T = laser cycles, from the quantum model. The intensity of radiation at a given time and frequency is indicated by a logarithmic color scale. In all panels, the red dashed lines indicate the kinetic energies and corresponding frequencies of 2Up + Ip and 3.17Up + Ip.

Figure 3 .

 3 Figure 3. Snapshots of the distribution function f (x, v, τ ; z) in a logarithmic scale. Left panels: f (x, v, τm; z), where τm > 0 is the first field intensity maximum after the start of the pulse. Right panels: f (x, v, τ0; z), where τ0 τm is the first zero of the field following τm. The upper panels are at z = 0, while the lower panels are at z = 0.37 mm.

Figure 1 .

 1 Figure 1. Schematic of the model, Eq. (1) in the paper.

- τ 0 E

 0 (τ ′ , z) + ∂V ∂x |ψ(x, τ ′ ; z)| 2 dx dτ ′ .(6)II. COMPUTATION OF THE OBSERVABLESTo perform the time-frequency analysis of an observable O(τ ; z), we computeÕ(τ, ω; z) = W T (τ ′ -τ )O(τ ′ ; z) exp(-iωτ ′ )dτ ′ 2 .

and T.U. acknowledge funding from the NSF (Grant No. PHY1602823). This material is based upon research supported by the Chateaubriand Fellowship of the Office for Science & Technology of the Embassy of France in the United States.


Olivier Pourret
email: olivier.pourret@lasalle-beauvais.fr.

Raul E Martinez

MODELLING LANTHANIDE SERIES BINDING SITES ON HUMIC ACID

Keywords: Rare earth elements, humic acid, binding, speciation, LPM model

Lanthanide (Ln) binding to humic acid (HA) has been investigated by combining ultrafiltration and ICP-MS techniques. A Langmuir sorption isotherm metal complexation model was used in conjunction with a linear programming method (LPM) to fit experimental data representing various experimental conditions both in HA/Ln ratio (varying between 5 and 20) and in pH range (from 2 to 10) with a ionic strength of 10 -3 mol L -1 . The LPM approach, not requiring prior knowledge of surface complexation parameters, was used to solve the existing discrepancies in LnHA binding constants and site densities. The application of the LPM to experimental data revealed the presence of two discrete metal binding sites at low humic acid concentrations, (5 mg L -1 ), with log metal complexation constants (log K S,j ) of 2.65 ± 0.05 and 7.00 (depending on Ln). The corresponding site densities were 2.71 ± 0.57 x 10 -8 and 0.58 ± 0.32

x 10 -8 moles of Ln 3+ /mg of HA (depending on Ln). Total site densities of 3.28 ± 0.28 x 10 -8 mol mg -1 , 4.99 ± 0.02 x 10 -8 mol mg -1 and 5.01 ± 0.01 x 10 -8 mol mg -1 were obtained by LPM for humic acid, for HA concentration of 5 mg L -1 , 10 mg L -1 and 20 mg L -1 , respectively. These results confirm that lanthanide binding occurs mainly at weak sites (i.e., ca. 80%) and secondly at strong sites (i.e., ca. 20%). The first group of discrete metal binding sites may be attributed to carboxylic groups (known to be the main binding sites of Ln in HA), and the second metal binding group to phenolic moieties. Moreover, this study evidences heterogeneity in the distribution of the binding sites among Ln. Eventually, the LPM approach produced feasible and reasonable results, but it was less sensitive to error and did not require an a priori assumption of the number and concentration of binding sites.

INTRODUCTION

The solution and mineral properties of lanthanides (Ln) make these trace elements excellent probes of low temperature geochemical reactions. Interest in Ln geochemistry comes from their systematic in chemical properties that often leads to fractionation in geochemical systems [START_REF] Henderson | Rare earth element geochemistry[END_REF]. Ln form a coherent geochemical group of trace elements that generally occur in the trivalent oxidation state. The effective ionic radii of Ln systematically decrease when atomic number increases [START_REF] Shannon | [END_REF] producing characteristic regular features of normalised Ln patterns defining the CHArge and RAdius-Controlled process, CHARAC [3]. If a low temperature geochemical system is characterized by CHARAC behaviour, elements of similar charge and radius such as Ln are, should display extremely coherent behaviour. This behaviour disappears when chemical processes are mainly driven by electronic external configuration producing sub-partition (non-CHARAC processes). CHARAC behaviour of Ln should thus generate smooth patterns whereas irregular patterns (excepting redox-related anomalies) may indicate non-CHARAC behaviour.

By studying fractionation trends, it then becomes possible to quantify the underlying fractionation processes in the natural environment (e.g., [3,4]). In addition, quantification of Ln fractionation in natural geochemical systems is essential for the modelling of the immobilization and transport of radioactive elements, as Ln are often used as naturally occurring analogs for the trivalent actinides [5].

In aquatic geochemical systems Ln fractionation occcurs by complexation to organic or inorganic ligands, or adsorption onto aquifer minerals. In natural waters, Ln may be associated with organic colloids which play a key role in complexing Ln elements and facilitate the fractionation of the Ln series (e.g., [6][7][8][9][10]). Dissolved organic matter (DOM), including humic acids (HA) and fulvic acids (FA), are abundant in surface water and groundwater systems [START_REF] Thurman | Organic Geochemistry of Natural Waters[END_REF].

HA are potent sorbents of dissolved metal cations like Ln elements (e.g., [10]). These properties emphasize the importance of HA in regulating the speciation, transport and subsequent fractionation of Ln in pristine and contaminated aquatic environments. The physicochemical quantification of Ln interaction with DOM in aquatic systems is needed to predict Ln fractionation by organic metal complexation in aqueous systems. HA contains an important number of potential metal complexing functional groups including carboxylic sites (e.g., [START_REF] Plaschke | [END_REF]).

Their chemical arrangement and conformation at the molecular level affects the ability of HA to complex Ln elements. For this reason, surface complexation modelling of LnHA acid interactions, using current surface complexation models such as FITEQL [START_REF] Herbelin | FITEQL, A computer program for determination of chemical equilibrium constants from experimental data[END_REF], WHAM and Model V and VI [START_REF] Tipping | [END_REF], or NICA Donnan [15] may provide different discrete metal binding constants for similar LnHA complexation reactions. Few studies have addressed LnHA interactions as a function of pH, ionic strength and different Ln elements. Only a small number of single Ln complexation constants with HA (SmHA, EuHA, TbHA and DyHA) have been evaluated [16][17][18][19][20][21]. Results from recent studies suggest that most of the LnHA complexation constants have to be measured or interpolated, especially for the whole Ln series [22], [23]. In these studies LnHA interactions have been investigated using binding interactions which are based on one-and two-site conditional binding to discretized multi-site (e.g., [24]), and continuous distribution models (e.g., [15]). Various electrostatic models have been developed to account for pH and I effects ion binding (e.g., [24]). Sonke and Salters [25] adopted a monodentate (carboxylic) binding mechanism. The exact binding site nature, multi-dentism, and polyelectrolyte effects of HA are ignored in the monodentate binding concept and therefore implicitly included in the conditional binding constant, which consequently is only valid for a specific pH, I, and temperature. Sonke [26], Pourret et al. [27] and Stern et al. [28] further investigated this feature using Model V and VI. However, such modelling of LnHA complexation only considers specific constants and thus does not implicitly distinguish different binding sites.

The effects of HA on Ln fractionation was previously studied by monitoring the concentration of Ln as a function of increasing concentrations of HA [27]. The resultant sorption data were analyzed to assess multi-site metal surface binding using linear programming regression methods (LPM) [29][30][31]. Most previous studies have investigated HA surface interactions with metals using single-site sorption isotherms that do not adequately describe the heterogeneity of HA surfaces [32]. Despite the previous application of LPM to describe metal binding to organic surfaces [29][30][31], this modelling approach has not yet been applied to study LnHA interactions.

MODELLING OF LnHA COMPLEXATION DATA

The experimental data used in this study has been previously published by Pourret et al. [27]. Ln complexation with HA was investigated using a standard batch equilibration technique.

Ln (ranging from 360 to 286 nmol L -1 depending on Ln) and Aldrich HA (5, 10 and 20 mg L -1 ) were placed together in solution, at an ionic strength of 1 x 10 -3 mol L -1 and at pH values ranging from 2.18 to 10.44. The present work will model Ln complexation to HA data using a surface complexation modelling approach, the LPM method as described previously by Martinez and Ferris [30] and Martinez et al. [31], which does not require prior knowledge of log K or binding site concentrations. In order to describe Ln 3+ complexation to deprotonated HA sites, HA -, a competition reaction was assumed to take place in a 1:1 ratio, as follows:

HHA j + Ln 3+ ⇋ LnHA 2+ j + H + (1)
where HHA j represents HA reactive site. K S,j represents the concentration apparent equilibrium constant for the reaction in Equation 1, conditional on ionic strength. For a j th deprotonated binding site at the i th step of the titration, K S,j can be defined as:

[ ] [ ] [ ] [ ] i meas S i j i meas i j S Ln HHA H LnHA K j , 3 , 2 , .
.

+ + + = (2)
where i = 1, . . ., n titrant additions and j = 1, . . ., m binding sites. In the above expression, K S,j implicitly embodies electrostatic parameters and is a function of experimentally determined proton and metal concentrations, ([H + ] meas,i and [Ln 3+ ] meas,i ) and of the amount of Ln 3+ bound to the j th site at the i th step of the titration, [LnHA 2+ j ] i . The total bound metal at the i th titrant addition, [LnHA 2+ ] T,i , and the total ligand concentration, [HA] T , can be expressed as:

[ ] [ ] [ ] [ ] i meas T i m j j i T Ln Ln LnHA LnHA , 3 3 1 2 , 2 + + = + + - = = ∑ (3) and [ ] [ ] [ ] [ ] [ j j j m j j T HHA HA LnHA HA HA + + = = - + = ∑ 2 1 ] (4)
where [HA - j ] refers to the individual site density for a particular surface functional group type.

The total concentration of bound metal, [LnHA 2+ ] T,i , can be expressed as a sum of complexed metal concentrations for each of the j th surface ligands at the i th step of the titration. However, experimental measurements of total, [Ln 3+ ] T , and free metal concentrations [Ln 3+ ] meas,i , only allow direct determination of [LnHA 2+ ] T,i , as indicated by Equation 3.

The fraction of the total j th ligand concentration, bound by Ln 3+ at the i th step of the titration, α LnHA,ij , can be expressed as a function of the bound metal (i.e., Ln 3+ ) at the i th titrant addition, [LnHA 2+ j ] i and the j th ligand concentration, [HA - j ] as follows:

[ ] [ ] [ ] [ ] [ ] i j i j i j j i j ij LnHA HHA LnHA LnHA HA LnHA + = = + + - + 2 2 2 , α (5)
The protonated j th ligand concentration at the i th step of the titration, [HHA j ] i , can in turn be expressed as a function of [LnHA 2+ j ] i , by rearranging the expression for the equilibrium constant K S,j in Equation 2. The calculated bound metal concentration at the i th titrant addition, [LnHA 2+ ] T,calc,i , can then be determined as a function of measured, ([H + ] meas,i and [Ln 3+ ] meas,i ) and adjustable ([HA - j ]) parameters, as shown below:

] [ ] [ ] [ ] [ ]) [ ( ] [ 1 , , , 3 , , 3 1 , , , 2 -
= + + + = - + ∑ ∑ ⋅ + = = j m j i meas j S i meas j S i meas m j j ij LnHA i calc T HA H K Ln K Ln HA LnHA α (6)
The linear programming approach for multisite metal sorption solves a matrix equation, b = A • x, for x. Here A is an n×m matrix of α LnHA,ij entries as defined in Equations 5 and 6. The b vector is a n×1 vector of calculated bound metal concentrations for each titrant addition, [LnHA 2+ ] T,calc,i , as defined in Equation 6. The m×1 vector x contains the adjustable parameters, [HA - j ], for each of the m binding sites. Numerical difficulties exist in attempting to fit the model in Equation 6because binding constants and site densities are correlated parameters. This problem is solved using a fixed interval grid of log K S,j values and writing the problem in matrix form, as described by Martinez et al. [31]. In addition, the nature of the matrices as described above makes this an ill-posed problem, meaning that more than one error minimum can be found from optimization for x as a solution to the equation b = A• x, unless additional assumptions are made about the nature of the solution [31].

Linear programming regression minimizes the number of binding sites and the absolute error, e = |[LnHA 2+ ] T,calc,i -[LnHA 2+ ] T,i |, using a simplex method [29]. This approach finds one global minimum for the error function, which emphasizes zero as a possible solution and avoids convergence problems such as those found in FITEQL [30,31]. LPM optimizes parameters such as total binding site concentrations. Each site density, [HA -], is assigned a positive value where zero is a possible result. This generates a log K S,j spectrum where discrete metal-binding sites are determined by the number of log K S,j values, which have a corresponding nonzero metal-bindingsite density. When [HA - j ] values are added, their sum should approximate the total available ligand concentration on the sorbent surface, [HA] T , for a maximum experimental pH value.

In our simulations, the binding of the Ln first hydrolysis product to HA was not considered. This choice is supported by the fact that (i) all but two data points (among 28) have pH < 7; yet, it is well established that the proportion of Ln-OH complexes and thus Ln-OH-HA complexes may become important only for water samples having pH > 8 [17,33]; (ii) even for alkaline waters, recent model calculations show that Ln speciation can be reasonably well captured by only considering Ln 3+ complexation with HA [34].

RESULTS AND DISCUSSION

Tables 1 to 3 summarize LPM optimization results of log K S,j, [HA j ] and [HA] T as defined by Equations 2 and 4 respectively for experimental [Ln 3+ ] sorption by three different initial HA concentrations of 5 mg L -1 to 20 mg L -1 , in the pH range of 2-10. Tables 1 and2 indicate the log K S,j and [HA] T values for all Ln elements at the full range of HA concentrations, while Table 3 shows [HA - j ] for 2 sites (i.e., j = 2) at 5 mg L -1 HA. LPM optimization of Equation 6finds an optimal solution set of binding site concentrations, [HA - j ], which are assigned to the corresponding log K S,j on fixed interval grid. This procedure generates a discrete spectrum and by approximating the ideal condition e = |[LnHA 2+ ] T,calc,i -[LnHA 2+ ] T,j | = 0. LPM optimization should generate a unique binding site density corresponding to a single log K S,j on the grid. However, double peaks resulted for particular sites j (data not shown), because the true log K S,j of the sample falls at an intermediate position between two adjacent log K S,j on the grid [31]. Each doublet was converted to a single peak by averaging the two log K S,j values and computing the weighted average of [HA - j ]. The averaged values, along with existing single peaks, in replicate spectra were used to calculate overall average log K S,j and [HA - j ] values (log K S,j(avg) and [HA - j ] (avg) ).

The log K S,j(avg) and binding site densities obtained for the three humic acid concentrations are reported in Tables 1 and2, for the 14 analyzed Ln elements respectively. The LPM fit results for Ln, as a function of pH, are illustrated in Figure 1 pattern showing a middle rare earth elements (MREE) downward concavity (Figure 2). As already proposed by Pourret et al. [27], literature compilation of REE-organic ligand constants as plotted against Gd/Yb ratio (Figure 3) evidenced that new calculated log K S,j values are in the field of natural carboxylic acids and phenolic acids. As the objective of this study was not to document fractionation among Ln series, differences in between Sonke and Salters [25] and Pourret et al. [27] were not discussed (for that see discussions in [27,28]).

As evidenced in Table 1, two log K S,j values are proposed by LPM for HA concentration of 5 mg L -1 . This feature reflects the heterogeneity of HA binding sites which is known to be composed from carboxylic group varying from 3.2 to 4.8 meq g -1 and phenolic ones from 1.4 meq g -1 to 3.4 meq g -1 [35]. This latter characteristic is only evidenced by the first experimental condition as it is the only one that spans up to a pH of 10, and thus only the conditions considering pH greater than 7, can allow development of binding sites at alkaline pH (i.e., phenolic or polycarboxylic ones). Indeed, the contribution of the second site by just doubling the concentration of HA from 5 mg L -1 to 10 mg L -1 , should still be expected at 10 mg L -1 and 20 mg L -1 , since the contribution of the second site may reached up to 30% (i.e., from Er to Lu) at 5 mg L -1 . However, analytical window does not allow this latter feature to express at both HA concentration of 10 mg L -1 and 20 mg L -1 (i.e., pH does not exceed pH values of 7.26 and 6.77, respectively) as phenolic functional groups mostly expressed from pH 8 to 10 [35]. This latter feature is also a verification of the results of Štamberg et al. [36] who previously studied Eu 3+ complexation to Aldrich HA and found that phenolic groups only occur, in Ln complexation, at alkaline pH after neutralization of carboxylic groups. The corresponding total site densities ([HA] T ) are 2.89 x 10 -8 mol g -1 , 5.05 x 10 -8 mol g -1 , and 5.02 x 10 -8 mol g -1 of HA, for 5 mg L -1 , 10 mg L -1 and 20 mg L -1 of HA, respectively for La (Table 2). These latter total site densities only correspond to ca. 1% of the available sites on HA surfaces. Among these sites, as displayed in Table 3 for HA concentration of 5 mg L -1 , two types of sites may be considered. Weaker sites (i.e., sites with a mean log K S,j value of 2.65) are predominant and represent 67 to 94% of the active binding sites whereas stronger sites (i.e., sites with a mean log K S,j value of 7.00) only represent 6 to 33% of the active binding sites. The corresponding site densities were 2.71 ± 0.57

x 10 -8 and 0.58 ± 0.32 x 10 -8 moles of Ln 3+ /mg of HA. As already observed, sites densities decrease with increasing log K S,j values [21].

Even if this study only considered monodentate sites, an interesting feature would be evidenced by the three HA concentrations. The binding of Ln by several separate complexing sites on HA may occur (i.e., carboxylic and phenolic binding sites) resulting in the fact that metal binding will depend upon the total HA concentration. In other words, two monodentate binding sites would result in lower total site densities available to complex Ln, as compared to a unique binding site. This lower total site density observed at the HA concentration of 5 mg L -1 , for all 14

Ln elements, indicates a direct correlation between the amount of Ln sorbed on the external surface microenvironment of HA which would increase with increasing HA concentration. This result confirms the results explained earlier, where the concentration of complexed metal is proportional to the concentration of added HA [START_REF] Tipping | Cation binding by humic substances[END_REF]. Oppositely, for the same conformation (i.e., one single monodentate site at a HA concentration of 10 mg L -1 and 20 mg L -1 ) there is no concentration dependence. Albeit not considered in this study, multidentate binding may also become increasingly important with increasing pH.

As evidenced by Table 3 and Figure 4, site density distribution varies in between Ln.

Sites 1 are more predominant (i.e., represent >80-90% of the binding sites) for the MREE (i.e., from Nd to Dy) whereas for La or Lu they only represent 77% or 67% of the available binding sites. Complementary, "sites 2" in Table 3 are more present for binding La and Lu as regards to MREE (i.e., 23% or 33% compared to <10-20%). This heterogeneity in the distribution of site types may thus reflect the competition among Ln elements which is responsible of the fractionation in between the Ln series. It must be noted that Ce does not follow this general trend which may be ascribed to its redox behaviour. Log K S,j values for the Ln series are maximum for Sm 3+ . Binding site densities display the same feature. It may be the result of the Ln coordination number change in the Ln complexes in solution from 9 to 8, which usually occurs in the Eu 3+ to Gd 3+ range. The water molecules about the inner Ln sphere is 9 from La 3+ to Nd 3+ , and 8 from Tb 3+ to Lu 3+ [START_REF] Brookins | [END_REF]. From a thermodynamical point of view the free energy of the Ln coordination number change in the Eu 3+ to Gd 3+ range as a result of the sudden change of the coordination polyhedrons type. Indeed, Choppin and Peterman [START_REF] Choppin | [END_REF] pointed out that Eu-acetate complexes occur as inner sphere and as acetate represents a model molecule of simple carboxylic sites on complex organic matter such as HA [9], such a coordination number change must be expected in the Eu 3+ to Gd 3+ range. These processes are induced by variations of inter electronic repulsive potentials due to progressive filling of 4f orbital during structural changes that especially involve inner coordination sphere of each Ln caused by equilibrium ligand exchange reaction (e.g., organic phase-aqueous phase equilibrium). Indeed, when processes involve Ln adsorption (or surface complexation) with an inner sphere mechanism, non-CHARAC effect may take place.

The complexation behaviour of Ln does not thus exclusively depend on its ionic charge and radius, but is additionally controlled by its electron configuration and by the type of complexing ligand, since this latter two determine the character of the chemical bonding. Hence, aqueous systems are characterized by non-CHARAC trace element behaviour [3], and electron structure must be considered as an additional parameter. However, such a behaviour need to be further explored and refined.

As already discussed by Pourret et al. [27] and Stern et al. [28] observed differences between experiments presented by Takahashi et al. [40], Yamamoto et al. [41,42], Sonke and Salters [25], and both of these studies (i.e., [27], [28]) strongly suggest the heterogeneity of the complexing sites in HS: high concentration of weak carboxylic sites and low concentrations of strong phenolic sites [35]. The weak sites determine the behaviour of humic complexation at high metal loading, whereas the strong sites determine the complexation strength of humic substances at trace metal concentrations [27], [28], even when their concentration is only in the range of a few percent of the weak sites (i.e., 10 to 22%; [35]). HS/Ln ratios considered in this study are between 5 and 20, whereas in Yamamoto et al. [41,42] experiments they are close to 80, in Sonke and Salters [25] and Stern et al. [28] between 500 and 700, and in Takahashi et al. [40] as high as 100,000. At lower loading, stronger sites are favored whereas at higher loadings, weaker sites are involved in the complexation. Moreover, phenolic sites are predominant at alkaline pH whereas carboxylic are dominant binding sites at acidic pH. Indeed, the pH range in this study spans from 2 to 10 with an analytical window in the 2-6 range, whereas in Yamamoto et al. [41,42] experimental pH varies between 4 and 5.5 and, in Sonke and Salters [25] and Stern et al. [28] it varies between 6 and 10. The present study sheds more light on the process dealing with Ln loading by HA and needs to be extended to other HS by modelling existing data [25], [28] with LPM to generalize this feature. In addition, all results, whether those of our study and Pourret et al. [27] or those of Sonke and Salters [26] and Stern et al. [28], suggest that the multi functionality of the organic matter surface site should be taken into account differently in speciation studies and calculation codes. Indeed, whether Model V and VI, each of the models considered the existence of two groups of surface sites that can form monodentate, bidentate (e.g., as in Model V; [43]) and tridentate (e.g., as in Model VI; [START_REF] Tipping | [END_REF]) complexes. However, their abundance and complexation constants are linked one to another by mathematical expressions and as regards to Ln behaviour (i.e., fractionation among Ln series as regards to binding site affinities), and may be refined considering the results obtained using LPM.

SUMMARY

Linear programming modelling was able to determine differences in Ln speciation as a function of increasing humic acid concentration. At low concentrations of dissolved organic matter two Ln binding sites were found, whereas one site is obtained at higher amounts.

Lanthanide binding occurs mainly at weak sites (i.e., ca. 80%; attributed to carboxylic groups) and secondly at strong sites (i.e., ca. 20%; attributed to phenolic moieties). However, heterogeneity in the distribution of the binding sites among Ln is evidenced, which could explain Ln fractionation patterns. Moreover, this study has been able to quantify LnHA interactions, including the determination of metal binding constants and site concentrations without the need for initial knowledge of these parameters, indicating the usefulness of the linear programming modelling approach to solving surface metal complexation in the presence of complex organic surfaces. Previous models including FITEQL have produced incomparable results emphasizing the use of LPM as a tool for quantification of metal-ligand interactions as previously suggested [29,31].

TABLES AND FIGURES CAPTIONS

Table 1. Log K s,j values.

Table 2. Total site density (in mol of Ln 3+ /mg of HA).

Table 3. Site density for both sites at a HA concentration of 5 mg L -1 (in mol of Ln 3+ /mg of HA).

Figure 1. Binding of (a) La, (b) Eu and (c) Lu to HA for various HA concentrations (i.e., 5 mg L - 1 , 10 mg L -1 and 20 mg L -1 ) as a function of pH. Points correspond to experimental data [22] whereas solid lines correspond to LPM best fits for which the absolute error e is minimal. when necessary; [6,9]). black triangles: amino-carboxylic acids; white squares: iminoacetic acids; black circles: phenolic acids; grey circles: carboxylic acids; white circles: natural carboxylic acids: dark grey circles: humic substances, Aldrich Humic Acid (AHA), Suwannee River Fulvic Acid (SRFA), Suwannee River Humic Acid (SRHA), Leonardite Humic Acid 3).

[HA] mg L -1

Figure 2 .

 2 Figure 2. Log K s,j patterns for the 14 Ln.

Figure 3 .

 3 Figure 3. Literature compilation of REE-organic ligand constants (recalculated at I = 0.1 mol L -1

(

  LHA), Suwannee River Natural Organic Matter (SRNOM) and Summit hill Soil Humic Acid (SSHA) ( a this study; b Yamamoto et al.[42]; c Sonke and Salters[25]; d Stern et al.[28]).

Figure 4 .

 4 Figure 4. Site density distribution (in mol of Ln 3+ /mg of HA) for the 14 Ln (data are from Table

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

  for La3+ . Sample log K S,j values for La 3+ are 2.65 ± 0.21 and 7, and 3.00 ± 0.14 and 3.85 ± 0.07 for[HA]

T = 5 mg L -1 , 10 mg L -1 and 20 mg L -1 , respectively. The main characteristics of the LnHA experiments, reproduced by LPM, are a low proportion of LnHA complexes at relatively low pH and a marked increase of this proportion with increasing pH. The log K S,j values increase with HA concentration of experimental solutions (e.g., for La from 2.65 to 3.85 for HA concentrations of 5 mg L -1 and 20 mg L -1 , respectively). Moreover, log K S,j values increase from La (2.65) to Pr (2.75), then decrease from Eu (2.70) to Lu (2.60) for HA concentrations of 5 mg L -1 . This feature is further marked for a higher HA concentration (i.e., 20 mg L -1 ) where log K S,j values increase from La (3.85) to Eu (4.15) and then decrease from Gd (4.06) to Lu (3.95). It follows a log K S,j

Table 2 .

 2 La 2.89 x 10 -8 5.05 x 10 -8 5.02 x 10 -8 Ce 3.31 x 10 -8 5.00 x 10 -8 5.02 x 10 -8 Pr 3.21 x 10 -8 4.99 x 10 -8 5.01 x 10 -8 Nd 3.65 x 10 -8 4.97 x 10 -8 5.00 x 10 -8 Sm 3.75 x 10 -8 4.98 x 10 -8 5.00 x 10 -8 Eu 3.65 x 10 -8 4.98 x 10 -8 5.01 x 10 -8 Gd 3.50 x 10 -8 4.99 x 10 -8 5.00 x 10 -8 Tb 3.44 x 10 -8 4.98 x 10 -8 5.01 x 10 -8 Dy 3.30 x 10 -8 4.98 x 10 -8 4.98 x 10 -8 Ho 3.20 x 10 -8 5.00 x 10 -8 5.02 x 10 -8 Er 3.04 x 10 -8 4.99 x 10 -8 5.01 x 10 -8 Tm 3.03 x 10 -8 4.98 x 10 -8 5.01 x 10 -8 Yb 3.01 x 10 -8 4.98 x 10 -8 5.00 x 10 -8 Lu 3.00 x 10 -8 4.98 x 10 -8 4.98 x 10 -8

	5	1 0	2 0

Acknowledgements

We thank the Biogeochemical Group from Géosciences Rennes for useful discussions when we were both in Rennes, and the Geologists Rock! Group in Beauvais when finalizing this project.

Two anonymous reviewers are acknowledged for useful comments and Catherine Wilson and Genevieve Green are acknowledged for editorial handling.


Joël Chaskalovic

Franck Assous
email: assous@ariel.ac.il

An optimal first-order Taylor-like formula with a minimized remainder

Keywords: Taylor's theorem, Taylor-like formula, Error estimate, Interpolation error, Approximation error, Finite elements

In this paper, we derive an optimal first-order Taylor-like formula. In a seminal paper [14], we introduced a new first-order Taylor-like formula that yields a reduced remainder compared to the classical Taylor's formula. Here, we relax the assumption of equally spaced points in our formula. Instead, we consider a sequence of unknown points and a sequence of unknown weights. Then, we solve an optimization problem to determine the best distribution of points and weights that ensures that the remainder is as minimal as possible.

Introduction

Even today, improving the accuracy of approximations remains a challenging problem in numerical analysis. In this context, Taylor's formula plays a crucial role in various domains, especially when one considers error estimates in numerical analysis to assess the accuracy of a numerical approximation method (for example, see [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] for finite element methods).

This challenge becomes even more crucial when comparing the relative accuracy between two given numerical methods. All error estimates share a common structure, whether applied to the finite elements method [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF], [START_REF] Ern | Theory and practice of finite elements[END_REF], numerical approximations of ordinary differential equations [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF], or to quadrature formulas used for approximating integrals [START_REF] Chaskalovic | A probabilistic approach for solutions of determinist PDE's as well as their finite element approximations[END_REF].

Let us specify these ideas in this context of numerical integration. Consider, for instance, a composite quadrature rule of order k. For a given interval [a, b], let f be a function in C k+1 ([a, b]). The corresponding error of the composite quadrature rule can be expressed as (refer to, e.g., [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF], [START_REF] Burden | Numerical Analysis[END_REF] or [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF]), for a non-zero integer N :

b a f (x)dx - N i=0 λ i f (x i ) ≤ C k h k+1 .
In this formula, h denotes the size of the N + 1 equally spaced panels [x i , x i+1 ], 0 ≤ i ≤ N , that discretize the interval [a, b], and λ i are N + 1 real numbers. Moreover, C k is an unknown constant, independent of h, but dependent on f and k. This constant is directly linked to the uncertainty associated with the remainder of Taylor's formula.

However, comparing the asymptotic convergence rate is no longer useful when evaluating two composite quadratures rules for a fixed value of h, as is common in various applications. Therefore, we focus our attention to minimize the constants C k by improving the estimation of the remainder involved in Taylor's formula.

In this context, several approaches have been proposed to determine a way to enhance the accuracy of approximation. For example, within the framework of numerical integration, we refer the reader to [START_REF] Barnett | Applications of Ostrowski's version of the Grüss inequality for trapezoid type rules[END_REF], [START_REF] Cerone | Trapezoidal-type rules from an inequalities point of view[END_REF] or [START_REF] Dragomir | Some remarks on the trapezoid rule in numerical integration[END_REF], and references therein. From another point of view, due to the lack of information, heuristic methods were considered, basically based on a probabilistic approach, see for instance [START_REF] Abdulle | A probabilistic finite element method based on random meshes: A posteriori error estimators and Bayesian inverse problems[END_REF], [START_REF] Assous | Indeterminate Constants in Numerical Approximations of PDE's: a Pilot Study Using Data Mining Techniques[END_REF], [START_REF] Hennig | Probabilistic numerics and uncertainty in computations[END_REF], [START_REF] Oates | A modern retrospective on probabilistic numerics[END_REF] or [START_REF] Chaskalovic | A probabilistic approach for solutions of determinist PDE's as well as their finite element approximations[END_REF], [START_REF] Chaskalovic | A new mixed functional-probabilistic approach for finite element accuracy[END_REF] and [START_REF] Chaskalovic | Explicit k-dependence for P k finite elements in W m,p error estimates: application to probabilistic laws for accuracy analysis[END_REF]. This allows to compare different numerical methods, and more precisely finite element, for a given fixed mesh size, [START_REF] Chaskalovic | Numerical validation of probabilistic laws to evaluate finite element error estimates[END_REF].

In this context, we recently developed a first-order Taylor-like formula in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF] and a second-order Taylor-like formula in [START_REF] Chaskalovic | A new second order Taylor-like theorem with an optimized reduced remainder[END_REF]. The goal was to minimize the corresponding remainder by transferring part of the numerical weight of this remainder to the polynomial involved in the Taylor expansion.

In both of these cases, we a priori introduced a linear combination of f (and f in [START_REF] Chaskalovic | A new second order Taylor-like theorem with an optimized reduced remainder[END_REF]) computed at equally spaced points in [a, b], and we determined the corresponding weights in order to minimize the remainder. We demonstrated that the associated upper bound in the error estimate is 2n times smaller than the classical one for the first-order Taylor's theorem, and 3/16n 2 times smaller than the classical second-order Taylors's formula.

In this paper, we relax the assumption of equally spaced points and consider a sequence of unknown points in the interval [a, b], where a given function f needs to be evaluated. Simultaneously, we introduce a sequence of unknown weights to be determined with the goal of minimizing the remainder. Then, we will prove that the remainder of the corresponding first-order expansion is minimized when the points between a and b are equally spaced.

The paper is organized as follows. In Section 2, we present a new first-order Taylor-like formula built on a sequence of given points x k , (k = 1, . . . , n -1), in [a, b], and given weights ω k , (k = 0, . . . , n). In Section 3 we derive the main results of this paper which deal with the optimal choice of points x k and weights ω k that enable us to minimize the corresponding remainder of the first-order Taylor-like formula. Concluding remarks follow.

A first-order Taylor-like theorem

Let us begin by recalling the well-known first order Taylor formula [START_REF] Taylor | Methodus incrementorum directa and inversa[END_REF]. For two reals a and b, a < b, we consider a function

f ∈ C 2 ([a, b]). Then, there exists (m 2 , M 2 ) ∈ R 2 such that, for all x ∈ [a, b], m 2 f (x) M 2
and we have In a previous paper [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF], we derived a first-order Taylor-like formula with the aim of minimizing the classical remainder a,1 (b). More precisely, we proved the following result:

f (b) = f (a) + (b -a)f (a) + (b -a) a,1 (b), (1)
Theorem 2.1 Let f be a real mapping defined on [a, b] which belongs to C 2 ([a, b]), such that: ∀x ∈ [a, b], -∞ < m 2 f (x) M 2 < +∞.
Then, for a given non-zero integer n, we have the following first-order expansion:

f (b) = f (a) + (b -a) f (a) + f (b) 2n + 1 n n-1 k=1 f a + k (b -a) n + (b -a) a,n+1 (b), (2)
where

| a,n+1 (b)| (b -a) 8n (M 2 -m 2 ).
Moreover, for an a priori choice of regularly spaced points a + k (b -a) n , the remainder a,n+1 (b) is minimum.

In the following of the paper, our goal is to relax the assumption of a priori equidistant points to determine the optimal set of points (x k ) k=0,n ∈ [a, b], along with the associated weights (ω k ) k=0,n . This determination will enable us to minimize the remainder a,n+1 (b) in (2).

To derive the main result below, we first introduce the function φ defined by

φ : [0, 1] -→ R t -→ f (a + t(b -a)) , (3)
that satisfies φ(0) = f (a) and φ(1) = f (b). Moreover, we proved in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF] that the remainder a,1 (b) introduced in (1) satisfies the following result:

Proposition 2.2 The remainder a,1 (b) in formula (1) can be expressed as a,1 (b) = 1 0 (1 -t)φ (t)dt . (4)
For a given integer n ∈ N * , we consider the set of points (x k ) k=0,n in the interval [a, b], where x 0 = a and x n = b, on the one hand, and a set of real weights (ω k ) k=0,n , on the other hand.

Then, we define the quantity * n (a, b) by the formula

f (b) = f (a) + (b -a) n k=0 ω k f (x k ) + (b -a) * n (a, b), (5)
where the two sequences (x k ) k=1,n-1 and (ω k ) k=0,n are to be determined to minimize the remainder * n (a, b).

For the upcoming result, we need to introduce some notations. We denote by t k , (k = 0, . . . , n), a sequence of real numbers that allows us to represent the points x k in [a, b] as a barycentric combination of a and b, that is:

x k = a + t k (b -a), (0 ≤ t k ≤ 1). ( 6
)
We also introduce S k as the partial sum of the weights ω j , 0 ≤ j ≤ k,

S k = k j=0 ω j . (7)
Now, we can prove the following result:

Theorem 2.3 Let f be a real mapping defined on [a, b] which belongs to C 2 ([a, b]).
Then, the remainder * n (a, b) defined by ( 5) satisfies:

(b -a) 2 n-1 k=0 m 2 (S k -t k ) 2 -M 2 (S k -t k+1 ) 2 ≤ * n (a, b) ≤ (b -a) 2 n-1 k=0 M 2 (S k -t k ) 2 -m 2 (S k -t k+1 ) 2 . ( 8
)
Proof : Using first the function φ defined in (3), formulas (1) and ( 5) lead to

f (b) -f (a) b -a = φ(0) + a,1 (b) = n k=0 ω k φ(t k ) + * n (a, b),
which can be re-written as:

* n (a, b) = φ(0) + a,1 (b) - n k=0 ω k φ(t k ), = φ(0) + 1 0 (1 -t)φ (t)dt - n k=0 ω k φ(t k ), = φ(1) - 1 0 tφ (t)dt - n k=0 ω k φ(t k ), = φ(1) - 1 0 tφ (t)dt + n k=0 ω k (φ(1) -φ(t k )) - n k=0 ω k φ(1), = 1 - n k=0 ω k φ(1) - 1 0 tφ (t)dt + n k=0 ω k (φ(1) -φ(t k )) . (9)
Let us assume that weights (ω k ) k=0,n fulfill the condition

n k=0 ω k = 1 , (10)
equation ( 9) can be expressed as

* n (a, b) = - 1 0 tφ (t)dt + n k=0 ω k (φ(1) -φ(t k )) , = - 1 0 tφ (t)dt + n k=0 ω k 1 t k φ (t)dt, = - n-1 k=0 t k+1 t k tφ (t)dt + n k=0 ω k 1 t k φ (t)dt. ( 11
)
Considering the second term of [START_REF] Chaskalovic | Numerical validation of probabilistic laws to evaluate finite element error estimates[END_REF], it can be transform as follows:

n k=0 ω k 1 t k φ (t)dt = 1 t 0 ω 0 φ (t)dt + 1 t 1 ω 1 φ (t)dt + • • • + 1 t n-1 ω n-1 φ (t)dt, = t 1 t 0 ω 0 φ (t)dt + t 2 t 1 (ω 0 + ω 1 )φ (t)dt + • • • + 1 t n-1 (ω 0 + • • • + ω n-1 )φ (t)dt. = n-1 k=0 t k+1 t k S k φ (t)dt, (12)
where S k is defined by [START_REF] Cerone | Trapezoidal-type rules from an inequalities point of view[END_REF].

Then, using [START_REF] Chaskalovic | A refined first-order expansion formula in R n : Application to interpolation and finite element error estimates[END_REF] in [START_REF] Chaskalovic | Numerical validation of probabilistic laws to evaluate finite element error estimates[END_REF]

enables to write * n (a, b) as * n (a, b) = n-1 k=0 t k+1 t k (S k -t)φ (t)dt. ( 13
)
In the estimations below, we will use that

∀ t ∈ [0, 1], φ (t) = (b -a)f (a + t(b -a)), and ∀ x ∈ [a, b], m 2 f (x) M 2 .
Next, to derive a double inequality on * n (a, b), we consider the three following cases, depending on the location of S k related to the interval [t k , t k+1 ].

1. If t k ≤ S k ≤ t k+1 , the integral in (13) can be decomposed as follows:

t k+1 t k (S k -t)φ (t)dt = S k t k (S k -t)φ (t)dt + t k+1 S k (S k -t)φ (t)dt. ( 14
)
Now, using that (S k -t) is positive on [t k , S k ], and negative on [S k , t k+1 ], we can write

           (b -a)m 2 S k t k (S k -t)dt ≤ S k t k (S k -t)φ (t)dt ≤ (b -a)M 2 S k t k (S k -t)dt, (b -a)M 2 t k+1 S k (S k -t)dt ≤ t k+1 S k (S k -t)φ (t)dt ≤ (b -a)m 2 t k+1 S k (S k -t)dt.
Summing up these two relations, we obtain first that

t k+1 t k (S k -t)φ (t)dt ≤ (b -a)M 2 S k t k (S k -t)dt + (b -a)m 2 t k+1 S k (S k -t)dt,
and also

t k+1 t k (S k -t)φ (t)dt ≥ (b -a)m 2 S k t k (S k -t)dt + (b -a)M 2 t k+1 S k (S k -t)dt .
By simply computing the integrals involved in these relations, we obtain, in that case, the estimate

(b -a) 2 m 2 (S k -t k ) 2 -M 2 (S k -t k+1 ) 2 ≤ t k+1 t k (S k -t)φ (t)dt ≤ (b -a) 2 M 2 (S k -t k ) 2 -m 2 (S k -t k+1 ) 2 . ( 15
)
2. If S k ≥ t k+1 then S k -t ≥ 0 for all t ∈ [t k , t k+1 ], and we have:

(b -a)m 2 t k+1 t k (S k -t)dt ≤ t k+1 t k (S k -t)φ (t)dt ≤ (b -a)M 2 t k+1 t k (S k -t)dt .
This yields

(b -a)m 2 2 (S k -t k ) 2 -(S k -t k+1 ) 2 ≤ t k+1 t k (S k -t)φ (t)dt ≤ (b -a)M 2 2 (S k -t k ) 2 -(S k -t k+1 ) 2 ,
which also leads to [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF], by simply using that m 2 ≤ M 2 .

3. If S k ≤ t k then S k -t ≤ 0 for all t ∈ [t k , t k+1 ], and in the same way as above, we get

(b -a)M 2 2 (S k -t k ) 2 -(S k -t k+1 ) 2 ≤ t k+1 t k (S k -t)φ (t)dt ≤ (b -a)m 2 2 (S k -t k ) 2 -(S k -t k+1 ) 2 ,
that also gives estimate [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF].

Hence, in all cases, we arrive at the same estimate [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF]. Finally, by summing over [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF] over all values of k from 0 to n-1 in (15), we get inequalities (8) for the remainder * n (a, b), which concludes the proof.

With the aim of minimizing * n (a, b), we introduce the function χ defined by

χ = (b -a)(M 2 -m 2 ) 2 n-1 k=0 k j=0 ω j -t k 2 + k j=0 ω j -t k+1 2 , (16)
which represents to the difference between the right-hand side and the left-hand side in [START_REF] Chaskalovic | A probabilistic approach for solutions of determinist PDE's as well as their finite element approximations[END_REF].

As a consequence, in the next section, we will minimize function χ which depends on 2n -1 variables, namely (t 1 , . . . , t n-1 , ω 0 , . . . , ω n-1 ). This is because the boundary points are known (x 0 = a corresponding to t 0 = 0 and x n = b corresponding to t n = 1), on the one hand, and due to relation [START_REF] Chaskalovic | Explicit k-dependence for P k finite elements in W m,p error estimates: application to probabilistic laws for accuracy analysis[END_REF] which will determine ω n , on the other hand.

3 The optimal first order Taylor-like formula

We begin this section by deriving a lemma that provides a necessary condition for the function χ to have an extremum at the point (t 1 , . . . , t n-1 , ω 0 , . . . , ω n-1 ).

Lemma 3.1 Let (t 1 , . . . , t n-1 , ω 0 , . . . , ω n-1 ) be an extremum of function χ.

Then, we have:

       S k = 1 2 (t k + t k+1
), (k = 0, . . . , n -1),

t k = S k-1 + ω k 2 , (k = 1, . . . , n -1). ( 17
) ( 18
)
Proof : The necessary conditions that guarantee that (t 1 , . . . , t n-1 , ω 0 , . . . , ω n-1 ) is an extremum is written as (see for instance [START_REF] Pshenichnyi | Necessary conditions for an extremum[END_REF]) ∀k = 0, . . . , n -1 : ∂χ ∂ω k = 0 and ∀k = 1, . . . , n -1 :

∂χ ∂t k = 0.
Regarding first the dependence of the function χ on the variables ω k , (k = 0, . . . , n -1), the conditions ∂χ ∂ω k = 0 are expressed as

n-1 m=k 2S m -t m -t m+1 = 0 , or equivalently n-1 m=k S m = 1 2 n-1 m=k (t m + t m+1 ) . ( 19
)
Since this system of equations is triangular, it can be easily solved. Writing two consecutive equations for a given k ∈ {0, . . . , n -2} leads to

S k + S k+1 + • • • + S n-1 = 1 2 (t k + t k+1 ) + (t k+1 + t k+2 ) + • • • + (t n-1 + t n ) , S k+1 + • • • + S n-1 = 1 2 (t k+1 + t k+2 ) + • • • + (t n-1 + t n ) ,
which readily gives, by difference,

S k = 1 2 (t k + t k+1 ) ,
the case k = n -1 corresponding directly to the last equation of the system [START_REF] Oates | A modern retrospective on probabilistic numerics[END_REF]. Now, to study the dependence of the function χ on the variables t k , ∀k = 1, . . . , n -1, we expand formula [START_REF] Dragomir | Some remarks on the trapezoid rule in numerical integration[END_REF], using that t 0 = 0 and t n = 1. We obtain that

2 (b -a)(M 2 -m 2 ) χ = n-1 k=0 ω 0 + • • • + ω k -t k 2 + ω 0 + • • • + ω k -t k+1 2 = (ω 0 -0) 2 + (ω 0 -t 1 ) 2 + (ω 0 + ω 1 -t 1 ) 2 + (ω 0 + ω 1 -t 2 ) 2 + . . . + (ω 0 + • • • + ω k -t k-1 ) 2 + (ω 0 + • • • + ω k -t k ) 2 . + . . . + (ω 0 + • • • + ω n-1 -t n-1 ) 2 + (ω 0 + • • • + ω n-1 -1) 2 .
So, by taking the derivative of the function χ with respect to t k , we obtain, for each k = 1, . . . , n-1:

∂χ ∂t k = 0 ⇔ 2(ω 0 + • • • + ω k-1 -t k ) + 2(ω 0 + • • • + ω k -t k ) = 0 .
This can be expressed as

S k-1 + S k -2t k = 0 that is, using the definition (7) of S k 2t k = 2S k-1 + ω k ,
which corresponds to [START_REF] Hennig | Probabilistic numerics and uncertainty in computations[END_REF].

From Lemma 3.1, we can now state the main result of this paper:

Theorem 3.2 Let f be a real function defined on [a, b] that belongs to C 2 ([a, b]).
Then, the optimal unknown weights (ω m ) m=0,n together with the optimal set of points (x m ) m=1,n-1 determined by the sequence of real numbers (t m ) m=1,n-1 that minimizes the remainder * n (a, b) defined by ( 5) are given by:

ω 0 = ω n = 1 2n and ω k = 1 n , ∀k = 1, . . . , n -1, (20)
t k = k n and x k = a + k n (b -a), ∀k = 1, . . . , n -1. ( 21
)
As a result, the corresponding optimal first-order Taylor-like formula is given by the following expression:

f (b) = f (a) + (b -a) f (a) + f (b) 2n + 1 n n-1 k=1 f a + k (b -a) n + (b -a) * n (a, b), (22)
with

| * n (a, b)| (b -a) 8n (M 2 -m 2 ). (23) Proof :
-We begin by proving that ω k is constant for all values of k belonging to {1, . . . , n -1} From ( 17) and ( 18), we have

2S k = t k + t k+1 = S k-1 + ω k 2 + S k + ω k+1 2 , ∀k = 1, . . . , n -2, that yields S k -S k-1 := ω k = ω k + ω k+1 2 .
Then,

ω k+1 = ω k , ∀k = 1, . . . , n -2, (24)
that corresponds to

ω 1 = • • • = ω n-1 .
-We will now establish the relation between ω 0 and ω k , for k = 1, . . . , n -1.

Firstly, let us write [START_REF] Ern | Theory and practice of finite elements[END_REF] for k = 0 and ( 18) for k = 1. Using that t 0 = 0, we obtain that

2S 0 = t 0 + t 1 = t 1 , t 1 = S 0 + ω 1 2 = ω 0 + ω 1 2 ,
from which we get

ω 1 = 2ω 0 .
This allows us to conclude, using (24), that

ω 1 = ω 2 = • • • = ω n-1 = 2ω 0 . (25)
-Let us compute now the value of ω 0 .

To this end, we write ( 17) and ( 18) for k = n -1. Given that t n = 1 and using (25), this yields

2S n-1 = t n-1 + t n = t n-1 + 1, = S n-2 + ω n-1 2 + 1.
Then, substituting the expressions of S n-1 and S n-2 , we get with (25) that

2ω 0 + 4(n -1)ω 0 = ω 0 + 2(n -2)ω 0 + ω n-1 2 + 1, that leads to 2(2n + 1)ω 0 -2 = ω n-1 = 2ω 0 that is ω 0 = 1 2n . (26)
-It remains now to compute the value of ω n . Using (25) and (26) gives directly that

ω 1 = ω 2 = • • • = ω n-1 = 1 n . (27)
Finally, the value of ω n is obtained from relation [START_REF] Chaskalovic | Explicit k-dependence for P k finite elements in W m,p error estimates: application to probabilistic laws for accuracy analysis[END_REF], that is

ω n = 1 -ω 0 - n-1 k=1 ω k = 1 - 1 2n - n -1 n = 1 2n .
-Let us now consider the values of the t k , for k = 1 . . . , n -1.

Using the expressions (18) of the optimal t k , together with the expressions (26) and ( 27) of the ω k , we obtain that

t k = ω 0 + k-1 j=1 ω j + ω k 2 , (k = 1, . . . , n -1) = 1 2n + k -1 n + 1 2n = k n ,
that yields with [START_REF] Burden | Numerical Analysis[END_REF], the following expressions of the optimal points x k :

x k = a + k n (b -a), (k = 1, . . . , n -1). ( 28
)
-We conclude the proof of this theorem by determining the optimal lower and upper bounds of the remainder * n (a, b) given in [START_REF] Chaskalovic | A probabilistic approach for solutions of determinist PDE's as well as their finite element approximations[END_REF].

1. We first evaluate the quantity

n-1 k=0 (S k -t k ) 2 :
Using the expression (20) of the ω j , we obtain that

S k = k j=0 ω j = ω 0 + k j=1 ω j = 1 2n + k n = 2k + 1 2n , (k = 0, . . . , n -1)
, so that, together with the expression (21) of the t k , we have

n-1 k=0 (S k -t k ) 2 = n-1 k=0 2k + 1 2n - k n 2 = 1 4n . (29)
2. Similarly, we evaluate the quantity

n-1 k=0 (S k -t k+1 ) 2 . The same arguments yield n-1 k=0 (S k -t k+1 ) 2 = n-2 k=0 (S k -t k+1 ) 2 + (S n-1 -t n ) 2 = n-2 k=0 2k + 1 2n - k + 1 n 2 + 2n -1 2n -1 2 = n-2 k=0 -1 2n 2 + 1 4n 2 = 1 4n . (30)
Therefore, combining (29) and (30), we obtain optimal lower bound and upper bounds in [START_REF] Chaskalovic | A probabilistic approach for solutions of determinist PDE's as well as their finite element approximations[END_REF] for the remainder * n (a, b) as follows:

(b -a) 8n (m 2 -M 2 ) ≤ * n (a, b) ≤ (b -a) 8n (M 2 -m 2 ),
which exactly corresponds to the result derived in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF].

Conclusion

In this paper, we have derived a new first-order Taylor-like formula constructed as a linear combination of the first derivatives of a given function, evaluated at specified points x k , (k = 1, . . . , n-1), within an interval [a, b]. These points are weighted by real numbers ω k , (k = 0, . . . , n). Unlike the approach in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF], the positions of these points and their corresponding weights were not fixed a priori. In particular, we relax the assumption of equally spaced points.

Our aim was to determine the optimal positions of these points together with the weights for obtaining the "best formula", in the sense that the corresponding remainder is as small as possible.

To achieve this, we have established an initial result that provides upper and lower bounds for the remainder.

Then, we proved the existence of an optimal set of points and weights that minimize the remainder in the first-order Taylor-like formula. This result corresponds to the formula we derived in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF],

where we explicitely set a priori the values of the points x k , (k = 1, . . . , n -1), as uniformly distributed within the interval [a, b].

So, the consequences derived in [START_REF] Chaskalovic | A new first order expansion formula with a reduced remainder[END_REF], mainly related to applications in error approximations, can be considered as optimal. Mainly, it treats on P 1 -interpolation error estimate, the corrected trapezoidal quadrature formula, and finite element error approximations. For example, using the corrected trapezoidal quadrature formula, we have obtained an upper bound which is two times smaller than the errors obtained by using the classical trapezoidal quadrature formula. It highlights the importance and impact of the new Taylor-like formula ( 22)-( 23) in assessing the accuracy of a given numerical approximation method.

This research can be extended to the case of dimension strictly greater than one. This extension requires a Taylor-like formula we have already derived in [START_REF] Chaskalovic | A refined first-order expansion formula in R n : Application to interpolation and finite element error estimates[END_REF]. Additionally, we could expand this work to a second-order Taylor-like formula, as proposed in [START_REF] Chaskalovic | A new second order Taylor-like theorem with an optimized reduced remainder[END_REF]. Both of these extensions will be explored to assess their impact on error estimates in the context of applications in numerical analysis.

Homages: The authors warmly dedicate this research to pay homage to the memory of Professors André Avez and Gérard Tronel, who largely promote the passion for research and teaching in mathematics of their students.