C Montella
email: claude.montella@lepmi.grenoble-inp.fr

Voigt circuit representation model for electrochemical impedances under finite-length diffusion conditions

Keywords: Modeling, Impedance, EIS, Finite-length diffusion, Finite-length Warburg, Voigt circuit

The Voigt circuit, which is the serial association of a number of RC parallel elements, should provide an adequate description of the finite-length diffusion impedance Z d , as conjectured and numerically experimented in previous electrochemical works. In this article, it is proved that the Voigt circuit, with infinitely many RC elements, is an exact representation model for Z d . Next, a generalization is proposed to establish the appropriate Voigt circuit for the electrode impedance, taking into account the electron transfer reaction kinetics and the double-layer charging process at the electrode/electrolyte interface. Theoretical derivation of the Voigt circuit parameters in terms of the kinetic and mass transport parameters of the electrochemical system makes it possible to discuss the dependence of the largest time constant of the circuit on the rate constant of the electrochemical reaction, the diffusion parameters of redox species, the double-layer capacitance at the interface, and the steady-state potential imposed for EIS calculation or measurement. The theoretical predictions in this article are used in an attempt to analyze some experimental results from the electrochemical literature.

Introduction

The Voigt circuit, also called Foster 1 realization of RC circuits in circuit synthesis theory [START_REF] Temes | Introduction to Circuit Synthesis and Design[END_REF], results from n serially connected RC elements, i.e. n R k C k ( ) elements where R k and C k denote the k-th resistance and capacitance, respectively, and the symbol '/' means parallel association. The Voigt circuit is often employed as a representation model for Electrochemical Impedance Spectroscopy (EIS) data [START_REF] Orazem | Electrochemical Impedance Spectroscopy[END_REF]. A series resistance R 0 is generally added to the circuit 1 to take into account the Ohmic resistances (electrolyte and connection), thus leading to the impedance:

Z V ω ( ) = R 0 + R k 1+ iωτ k k=1 n ∑ (1)
where i = -1 ( ) 1 2 , ω = 2π f is the angular frequency, f is the frequency, and τ k = R k C k is the k-th time constant of the Voigt circuit. The condition (τ 1 > τ 2 > L > τ n-1 > τ n ) is assumed in the above equation to order its successive terms.

In their pioneering work focusing on the so-called measurement model, Agarwal et al. [START_REF] Agarwal | Measurement Models for Electrochemical Impedance Spectroscopy, I. Demonstration and Applicability[END_REF] postulated that the Voigt circuit is a good representation model for a wide variety of electrochemical impedance spectra. As an example, starting from synthetic impedance data generated from the equivalent circuit proposed for the hydrogen evolution reaction on LaNi 5 , which includes a Finite Length Warburg (FLW) impedance, these authors demonstrated numerically that such impedance data are well fitted by the impedance in Eq. [START_REF] Temes | Introduction to Circuit Synthesis and Design[END_REF]. They concluded that, by fitting experimental EIS data to the Voigt circuit impedance at an increasing number (n) of RC elements, the fit residuals should give access to the measurement noise of the EIS data. This was the starting point for using the Voigt circuit as the relevant measurement model to study the error structure of experimental impedance data [START_REF] Agarwal | Measurement Models for Electrochemical Impedance Spectroscopy, I. Demonstration and Applicability[END_REF][START_REF] Agarwal | Application of Measurement Models to Impedance Spectroscopy: II. Determination of the Stochastic Contribution to the Error Structure[END_REF][START_REF] Agarwal | Application of Measurement Models to Impedance Spectroscopy: III. Evaluation of Consistency with the Kramers-Kronig Relations[END_REF]. Orazem and Tribollet [2] developed this approach by taking the example of EIS measurements for the ferro-ferricyanide redox couple/Pt system. 1 The condition ( R 0 = 0 ) used below in this article does not affect the theoretical derivations in the following Sections.

Agarwal et al. [START_REF] Agarwal | Measurement Models for Electrochemical Impedance Spectroscopy, I. Demonstration and Applicability[END_REF] also noted that the Voigt circuit impedance satisfies the Kramers-Kronig transformation. They concluded that, if experimental EIS data are well fitted by the Voigt circuit impedance, these data should also satisfy the Kramers-Kronig relationships, which is an alternative test for validation of EIS measurements.

Boukamp [START_REF] Boukamp | A linear Kronig-Kramers Transform Test for Immitance Data Validation[END_REF] proposed a linear fit procedure in the parameters R k , so-called linear KK test, assuming that the time constants τ k of the Voigt circuit are distributed logarithmically according to the inverse of the angular frequencies of the input signal. Hence, this author showed that a Voigt circuit including n RC elements could exactly represent any linear circuit at n experimental impedance points although no discussion of the values of the circuit parameters was provided. As pointed out by a referee of this article, "Boukamp wrote a KK test program that carries out exact approximation in all the frequency points of the data file (to check KK transforms). The program KKtest.exe might be found on the Boukamp Web page and it is also included in the free Extra Material to the EIS book [START_REF] Lasia | Electrochemical impedance spectroscopy and its applications[END_REF]. The above program deals with both diffusion and kinetic-diffusion impedances".

More recently, Schönleber and Ivers-Tiffée [START_REF] Schönleber | Approximability of impedance spectra by RC elements and implications for impedance analysis[END_REF] questioned about the approximability of impedance spectra by serially connected RC elements, and its implications for impedance analysis.

They demonstrated that "the impedance of every possible non-oscillating electrochemical process can be approximated by a simple serial connection of generalized RC elements". They also acknowledged that RQ-, Gerischer-or Warburg-type elements could also be approximated by generalized RC elements, because of the possibility of obtaining approximated integer order models for fractional order systems [START_REF] Vinagre | Some Approximations of Fractional Order Operators used in Control Theory and Applications[END_REF].

The experimental approach by Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF] is quite different. Using a rotating disk electrode (RDE) device, these authors evaluated the low-frequency (LF) time constant of the ferriferrocyanide redox couple/Pt-RDE system using the Voigt circuit and factorized modeling to fit their EIS measurements, with application to impedance extrapolation to zero frequency and determination of the LF capacitance of the electrochemical system. For this purpose, they focused on the best-fit estimates of Voigt circuit parameters, and, more especially, on the best-fit estimate of the first (largest) time constant τ1 of this circuit. They observed that, by increasing the number of RC elements to n = 6, τ1 seemingly converges to a constant value, unlike the other time constant esti- mates for k = 2, 3, ... They also noted that τ1 is nearly constant under different steady-state polari- zation conditions, i.e. under equilibrium as well as non-equilibrium conditions, although the associated resistance and capacitance are potential-dependent components of the Voigt circuit. Some questions then arise. Why does the estimate of the first time constant of the best-fit Voigt circuit converge to a constant value when the number of RC elements increases to at least n = 6, unlike the other time constant estimates? How do the components R 1 , C 1 and the time con- stant τ 1 connect to the mass transport process near the electrode? In other words, is it possible to predict theoretically the numerical values of R 1 , C 1 and τ 1 from the electrochemical system pa- rameters? An additional question focuses on the invariance (or not) of the time constant τ 1 upon the steady-state value of potential imposed to the electrode prior to EIS measurements.

In order to provide some answers to the above questions, we need a model for the diffusionconvection impedance of RDE. Although the FLW model is a crude approximation of the mass transport process near RDE, it provides the simplest way to explain, at least qualitatively, the experimental observations by Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF]. This point deserves to be mathematically clarified, first in the simple case where the electrode impedance is equal to the FLW impedance, and next in the more general case where the FLW impedance is an electrochemical component of equivalent circuits involving other components like the electron transfer resistance and the double-layer capacitance. It is the aim of the present work, taking the example of a one-step electrochemical reaction.

Finite-length Warburg impedance

The electrochemical reaction dealt with in this work is the E reaction [START_REF] Bard | Electrochemical Methods. Fundamentals and Applications[END_REF] involving two soluble species, O and R, in the presence of a large excess of supporting electrolyte:

O + z e ↔ R (2)
The E reaction takes place on the surface of a uniformly accessible electrode, so that edge effects can be disregarded and one-dimensional diffusion conditions are satisfied in the electrolytic solution. The FLW impedance Z d ω ( ) results from the assumption of a constant concentration of redox species at the distance δ from the electrode surface, under static as well as dynamic condi- tions. This impedance is the sum of the two contributions from species O and R:

Z d ω ( ) = Z dO ω ( ) + Z dR ω ( ) (3)
For the sake of simplicity, we assume below that the diffusion coefficient is the same for both redox species, D O = D R = D , so that the impedances Z dO , Z dR and therefore Z d have the same diffusion time constant, τ dO = τ dR = τ d = δ 2 D, where δ is the diffusion layer thickness.

The FLW impedance has the well-known formulation [START_REF] Llopis | Study of the Impedance of a Platinum Electrode Acting as Anode[END_REF][START_REF] Drossbach | Elektrochemische untersuchungen an kohleelektroden-I: Die überspannung des wasserstoffs[END_REF][START_REF] Schuhmann | Sur l'impédance de diffusion en basse fréquence[END_REF]:

Z d ω ( ) = R d tanh iωτ d ( ) 1 2     iωτ d ( ) 1 2 (4)
and the associated diffusion resistance R d is defined as the low-frequency real limit:

R d = lim ω→0 Z d ω ( ) (5)
We will first show in Section 3 that the Voigt circuit impedance with infinitely many terms, Eq. (1) with R 0 = 0 and n = ∞ , is an exact representation model for the FLW impedance, Eq. (4), in the ideal case of Nernstian redox systems investigated without double-layer charging complication. Next, we will investigate in Section 4 the effect of the double-layer charging process on the Voigt circuit parameters. Finally, we will focus in Section 5 on the more general case of slow electron-transfer reaction kinetics for non-Nernstian redox systems.

Nernstian redox systems investigated without double-layer charging effect

Further simplification for EIS investigation of the E reaction results from the assumption of very fast electron transfer kinetics at the electrode surface, so-called E r reaction [11] involving a Nernstian redox system. The Faradaic impedance for the E r reaction is equal to the diffusion im- pedance, in the presence of a large excess of supporting electrolyte. Moreover, in this Section, we disregard the Ohmic potential drop in the electrolytic solution, as well as the double-layer charging process at the electrolyte/electrode interface, so that the electrode impedance is equal to the Faradaic impedance, and, therefore, to the diffusion impedance, Z ω

( ) = Z d ω ( ).
Starting from the infinite product expansions of hyperbolic sine and cosine functions [15],

the following formulation applies to the function tanh z ( ) z, with complex argument 2 z:

tanh z ( ) z = 1+ z 2 kπ ( ) 2 1+ 4z 2 2k -1 ( )π     2 k=1 ∞ ∏ (6)
as used in our previous work [START_REF] Berthier | Hopf bifurcation and sign of the transfer resistance[END_REF]. Partial-fraction expansion of the above function follows directly from logarithmic differentiation of the infinite product expansion of hyperbolic cosine:

tanh z ( ) z = 8 1 2k -1 ( ) 2 π 2 + 4z 2 k=1 ∞ ∑ (7)
The argument iωτ d ( ) 1 2 can be substituted for z in the above equation to derive the equiva- lent formulation of the FLW impedance, Eq. ( 4), as:

Z d ω ( ) = 8R d 1 2k -1 ( ) 2 π 2 + 4iωτ d k=1 ∞ ∑ (8)
2 The usual notation for the complex argument, z ≠ ± 2k -1 ( )iπ 2 , should not be confused with the stoichiometric number of electron for the E reaction in Eq. [START_REF] Orazem | Electrochemical Impedance Spectroscopy[END_REF].

which is the same as the Voigt circuit impedance in Eq. ( 1) where R 0 = 0 , n is infinite, and the parameters of the k-th RC element have the relevant expressions:

R k = 8 2k -1 ( ) 2 π 2 R d (9)
and:

τ k = 4 2k -1 ( ) 2 π 2 τ d (10)
In addition, because of τ k = R k C k , the Voigt circuit capacitances have the same formulation irrespective of the index k:

C k = τ d 2R d (11)
It is worth noticing that the time constant formulation in Eq. ( 10) is formally the same as that derived recently by Boukamp [START_REF] Boukamp | Derivation of a Distribution of Relaxation Times for the (fractal) Finite Length Warburg[END_REF] in his work focusing on the distribution of relaxation times for the FLW impedance, although different notation is used in the two articles. In addition, the resistances R k and the time constants τ k decrease rapidly as the index k rises in Eqs. ( 9) and (10), according to the sequence 1, 1 9, 1 25, 1 49, 1 81, … { } for both resistance and time constant ratios, R k R 1 and τ k τ 1 , at k = 1, 2, 3, 4, 5, …, so that a Voigt circuit composed of the serial association of a few RC elements should be enough to represent the FLW impedance in practice.

Effect of double-layer charging for Nernstian redox systems 4.1. Appropriate Voigt circuit

Disregarding the Ohmic potential drop in the electrolytic solution, the electrode impedance

for the E r reaction, with double-layer charging effect, reads:

Z ω ( ) = 1 iωC dl + 1 Z d ω ( ) (12)
where C dl is the double-layer capacitance and the diffusion impedance Z d ω ( ) is given by Eq. ( 4).

Introducing the time constant ratio:

r τ = R d C dl τ d ( 13
)
and using the extension of the Heaviside method [START_REF] Angot | Compléments de Mathématiques[END_REF] to deal with irrational impedances, as presented in the Appendix, the electrode impedance, Eq. ( 12), admits the partial-fraction expansion:

Z ω ( ) = 2R d 1 1+ r τ + r τ 2 q k 2 ( ) q k 2 + iωτ d ( ) k=1 ∞ ∑ ( 14
)
where q k is the k-th positive root of the generating equation:

cot q q = r τ (15)
The electrode impedance in Eq. ( 14) is the same as the Voigt circuit impedance in Eq. ( 1)

where R 0 = 0 , n is infinite, and the parameters of the k-th RC element have the relevant expres- sions:

R k = 2R d q k 2 1+ r τ + r τ 2 q k 2 ( ) (16)
and:

τ k = τ d q k 2 (17)
In addition, because of τ k = R k C k , the k-th capacitance of the Voigt circuit can be derived as:

C k = τ d 2R d 1+ r τ + r τ 2 q k 2 ( ) (18)
which depends both on the index k and the steady-state potential imposed to the electrode through the diffusion resistance R d , and, therefore, the time constant ratio r τ in Eq. (13).

The previous results in Section 3 are recovered setting r τ = 0 in Eqs. ( 14)- (18). In contrast, at large enough r τ values, Eq. ( 15) predicts that its first (smallest) positive root tends towards 1 r τ , while the other positive roots satisfy q k = k -1 ( )π for k = 2, 3, … Using Eq. ( 17), the time constant ratio τ 1 τ d approaches r τ , and therefore τ 1 is close to R d C dl , while the other time con-

stants satisfy τ k = τ d k -1 ( ) 2 π 2     at k >1.
In the same way, a glance at Eq. ( 16) shows that the resistance R 1 approaches R d , while the other resistances of Voigt circuit tend towards zero. Hence, at large enough r τ values, the appropriate Voigt circuit reduces to its first element whose imped-

ance is R d 1+ iω R d C dl ( ) .
In a more general way, it can be concluded that the number of RC ele- ments with non-vanishing impedances should decrease in the Voigt equivalent circuit for the electrode impedance at increasing r τ values.

Effect of the steady-state current/potential

The FLW resistance has been derived for Nernstian redox systems in our textbook [START_REF] Montella | Exercices de cinétique électrochimique, II. Méthode d'impédance[END_REF].

Assuming that D O = D R = D , as indicated above, one obtains the simplified formulation:

R d = 2 1+ coshξ ( ) z 2 f N FAm c O * + c R * ( ) (19)
Here z is the stoichiometric number of electron in Eq. ( 2), F is Faraday's constant,

f N = F RT ( ) is the Nernst constant, R is the ideal gas constant, T is the absolute temperature, ξ = zf N E -E°'
( ) is the dimensionless potential, E is the steady-state potential imposed to the elec- trode, E°' is the formal potential of redox couple, A is the electrode surface area, c O * and c R * are the initial/bulk concentrations of species O and R, respectively, and m = D δ is the mass transport constant for redox species under steady-state conditions.

The FLW resistance for Nernstian redox systems depends on the steady-state potential through the dimensionless potential ξ in Eq. (19). Exponential increase of R d is predicted at large ξ values, i.e. when the steady-state current approaches the mass-transport-limited current plateau in the anodic or cathodic direction. All parameters of the Voigt circuit, R k , C k andτ k , for k = 1, 2, ..., depend on the steady-state potential through Eqs. ( 13) and ( 15)- [START_REF] Montella | Exercices de cinétique électrochimique, II. Méthode d'impédance[END_REF]. For shortening reasons, the discussion below is limited to the first time constant of the Voigt circuit, which is its largest time constant evaluated from Eq. ( 17) as τ 1 = τ d q 1 2 , where q 1 is the first (smallest) positive root of Eq. ( 15).

The equality of the diffusion coefficients of redox species has been tacitly assumed in this article. If the following additional condition is imposed on the initial/bulk concentrations,

c O * = c R * = c * ,
for the sake of simplification, the time constant ratio derived from Eqs. ( 13) and (19) can be reformulated as a function of the dimensionless current ψ as:

r τ = 2C dl z 2 f N FAδ c * 1-ψ 2 ( ) (20)
where ψ = I I lim is the ratio of the steady-state current I to the mass-transport-controlled current limit in the anodic direction ( I d,a = I lim ), as well as minus the current limit in the cathodic direction

( I d,c = -I lim ).
For illustration purpose, let us assume that EIS measurement are carried out at 25 °C for a one-electron transfer reaction (z = 1) with very fast electron-transfer reaction kinetics. Using the double-layer capacitance per unit of electrode surface area ( C dl A = 3.5 × 10 -5 F cm -2 ), the initial concentration of electroactive redox species ( c * = 2 × 10 -5 mol cm -3 ) and the diffusion layer thick- ness (δ = 1.6 × 10 -3 cm ), we get r τ = 5.825 × 10 -4 1-ψ 2 ( ) from Eq. ( 20). The minimum r τ value is observed at equilibrium ( ξ = 0 and ψ = 0 ). At the opposite, r τ becomes infinite when ψ → ±1, in the direction of oxidation (plus sign) or reduction (minus sign). The time constant ratio, τ 1 τ d , can be evaluated as a function of ψ solving numerically Eqs. ( 15), ( 17) and (20). The diffusion time constant, τ d = δ 2 D , only depends on the diffusion parameters D and δ , which are invariant with the electrode potential. Hence, the graphical plot in Fig. 1 predicts that the first time constant of the Voigt circuit is nearly constant along the steady-state current-potential curve, except when the current is very close to the diffusion current plateau in the anodic or cathodic direction (ψ → ±1).

In addition, because of the very small r τ value observed at equilibrium, i.e. r τ = 5.825 × 10 -4 , the nearly constant value predicted for τ 1 along the steady-state current-potential curve is close to its asymptotic value derived in Section 3 in the absence of double-layer charging effect, which is

τ 1 = 4τ d π 2 ≈ 0.4053 τ d .
Please insert Fig. 1 at this place

Effect of slow electron-transfer reaction kinetics

Appropriate Voigt circuit

The E r reaction is the hypothetical version of the E reaction pertaining to infinitely fast electron-transfer kinetics at the electrode surface. Such a condition does not hold for real systems.

As a consequence, some additional complication for the Voigt circuit representation model should result from the assumption of slow electron-transfer reaction kinetics. The electrode impedance then includes the electron transfer resistance R t , which is a potential-dependent impedance compo- nent. Disregarding the Ohmic potential drop in the electrolytic solution, the electrode impedance for the E reaction reads:

Z ω ( ) = 1 iωC dl + 1 R t + Z d ω ( )     (21)
where the diffusion impedance is still given by Eq. ( 4). We introduce the resistance ratio:

r R = R t R d ( 22
)
which depends on the electron transfer reaction kinetics, the mass transport conditions and the steady-state potential imposed to the electrode. Using the same simplifications as in Section 4, i.e.

D O = D R = D and c O * = c R * = c * ,
the resistance ratio is maximum at the equilibrium potential of the electrode. This ratio tends towards zero as the current approaches the diffusion current plateau in the anodic or cathodic direction (ψ → ±1). Derivation of the Voigt circuit representation model for the electrode impedance in Eq. ( 21) follows the same calculation procedure as in Section 4.1 and in the Appendix, although the theoretical formulations are more intricate than the previous ones.

The successive positive roots q k of the following equation, which is the generalization of Eq. ( 15) to the case of non-Nernstian redox systems ( r R ≠ 0 ), generate the poles p k = -q k 2 τ d :

cot q q 1+ r R q cot q ( ) = r τ (23)
while the generating equation for the impedance zeros, z k =q z,k 2 τ d , is:

q z cot q z = -1 r R (24)
irrespective of r τ . The electrode impedance in Eq. ( 21) admits the following infinite product expan- sion where R p denotes the polarization resistance, R p = R t + R d :

Z ω ( ) = R p 1-iω z k 1-iω p k k=1 ∞ ∏ (25)
In addition, the locations of poles and zeros in Fig. 2 have the same features as those in 21) and ( 25) can be realized as a RC network.

Please insert Fig. 2 at this place

In fine, the Voigt circuit, with n = ∞ , is an exact representation model for the electrochemi- cal impedance Z ω ( ) in Eq. ( 21). Its k-th time constant is given by Eq. ( 17) using the k-th positive root q k of Eq. ( 23). The k-th resistance can be derived as:

R k = 2R d 1+ r R q k cot q k ( ) q k 2 1+ cot q k q k + cot q k ( ) 2 -r R r τ q k 2 1-cot q k q k + cot q k ( ) 2             (26)
used together with Eq. ( 23). Finally, because of τ k = R k C k , the k-th capacitance of the Voigt circuit is given by:

C k = τ d 2R d ⋅ 1+ cot q k q k + cot q k ( ) 2 -r R r τ q k 2 1- cot q k q k + cot q k ( ) 2       1+ r R q k cot q k ( 27
)
also used with Eq. ( 23). The Voigt circuit parameters formulated from Eqs. ( 17), ( 23), ( 26) and [START_REF] Curtain | Transfer Functions of Distributed Parameter Systems: A Tutorial[END_REF] depend on the index k, the mass transport parameters, the electron transfer reaction kinetics, the initial concentration of redox species, the double-layer capacitance, and the steady-state potential imposed to the electrode, which are involved in the time constant ratio r τ , Eq. ( 13), and the re- sistance ratio r R , Eq. ( 22). The simplified formulations of R k , C k and τ k derived in Section 3 are recovered setting r R = 0 and r τ = 0, while the intermediate formulations derived in Section 4.1 result from the condition r R = 0 .

The approximating (truncated) version of Eq. ( 25), limited to its first n terms, is:

Z ap,1 n,ω ( ) = R p 1-iω z k 1-iω p k k=1 n ∏ (28)
It can be compared to the following approximation resulting from the Voigt circuit representation model:

Z ap,2 n,ω ( ) = R p -R k k=1 n ∑ + R k 1-iω p k k=1 n ∑ (29)
The rates of convergence of the two approximations towards the electrode impedance Z ω ( ) are compared in Fig. 3. Clearly, the approximation based on the Voigt circuit is the most efficient because of its higher convergence rate.

Several asymptotic limits are of interest for this time constant. First, τ 1 = 4τ d π 2 has been predicted in Section 3 for Nernstian redox systems ( r R = 0 ) investigated without double-layer charging effect ( r τ = 0). Next, τ 1 = R d C dl has been predicted in Section 4.2 at large r τ values and negligible r R values. A third possible expression is τ 1 = R t C dl when the resistance ratio r R is very large. Another limiting expression is τ

1 = R t + R d ( )C dl .
The validity domains for those formulations can be visualized by plotting the zone diagram in Fig. 4 with log-log coordinates. In each zone of Fig. 4, the relative deviation of the time constant τ 1 from its asymptotic formulation is less than 5 %. It is worth noticing that the validity domain for the limiting formulation (τ 1 = 4τ d π 2 ) on bottom side of the zone diagram covers the case of ideal Nernstian systems ( r R = 0 ) investigated in Section 3, as well as the case of non-Nernstian systems ( r R > 0 ) investigated in this Section, pro- vided that the time constant ratio,

r τ 1+ r R ( ) = R t + R d ( )C dl τ d , is low enough.

Please insert Fig. 4 at this place

Given numerical values for the dimensionless kinetic parameter, Λ = k°m , and the dimen- sionless capacitive parameter, Κ = C dl z 2 f N FAδ c * ( ) , the equilibrium condition ( ξ = 0 , subscript 'eq') for τ 1 is represented in the zone diagram by the dot with coordinates log r R,eq = -log 2Λ ( ) and log r τ ,eq = log 2Κ

( ), on bottom side of Fig. 4. If the steady-state potential E is increased or de- creased with respect to its equilibrium value, which corresponds to increasing ξ values in Eqs. ( 30) and (31), the characteristic point, with coordinates log r R and log r τ , moves away from its equilibrium position in the zone diagram. The trace left by the moving point is the curve, so-called trajectory below, whose equation can be derived by elimination of the dimensionless potential ξ between Eqs. (30) and (31) as:

log r τ = log 2Κ ( ) + log 1+ 1 r R 1+ 4Λ 2 r R (32)
It is worth noticing that the asymptotic moving direction at very small r R values is log 2Κ

( ) -log r R , irrespective of the Λ value.

Three examples of trajectories are plotted in the zone diagram of Fig. 4, for Κ = 10 -3 and Λ = 10 -3 , 2 and 8.5 , which corresponds, respectively, to very slow, intermediate and fast electron- transfer reaction kinetics. The intersections of those trajectories with the limiting zones of the zone diagram make it possible to predict the possible sequences of τ 1 formulations/values that can be observed when the steady-state potential moves away from its equilibrium value. Keeping in mind that the electron transfer resistance R t and the diffusion resistance R d are potential-dependent components of the electrode impedance, the first time constant of the Voigt equivalent circuit is nearly constant only when it is close to 4τ d π 2 in the bottom zone of Fig. 4. Another information, derived from the observation of Fig. 4, is that the asymptotic formulations/values of the first time constant, R t + R d ( )C dl and R d C dl , although mathematically predictable, are unlikely to be observed in practice.

Attempt analysis of previous experimental results from the literature

Here we focus on the recent article by Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF], which is concerned with the Voigt circuit representation model of the electrode impedance for the ferro-ferricyanide redox couple/Pt-RDE system. As indicated in the introduction Section of this article, the above authors investigated this electrochemical system using the Voigt circuit and factorized modeling to fit their EIS data. They focused more especially on the best-fit estimate τ1 of the first (largest) time con- stant of the Voigt circuit 3 . They observed that, by increasing the number of RC elements to n = 6, τ1 seemingly converges to a constant value, unlike the other estimates of time constants for k = 2, 3, … They also noted that the converged value of τ1 is nearly constant under different steady-state polarization conditions, i.e. under equilibrium as well as non-equilibrium conditions.

The finite-length diffusion model is often employed as an approximation of the diffusionconvection process near RDE, with the diffusion layer thickness δ being replaced by the character- istic length of diffusion-convection derived by Levich [START_REF] Levich | Physicochemical Hydrodynamics[END_REF]. Thus, the theoretical predictions in this work can be applied, at least as a first approach, to the experimental observations by Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF].

Let us consider the Nyquist graph for the electrode impedance Z ω ( ) in Fig. 5(a), which has been plotted from Eqs. ( 4) and ( 21 Please insert Fig. 5 at this place 4 "Noise-free impedance data" means that no random noise has been added before fitting the impedance. Of course, the computational noise of the synthetic impedance data is present. The impact of additional random noise on the complex nonlinear least squares fit could be evaluated numerically without major difficulty, but it is not the aim of this work.

The use of noise-free EIS data makes it possible to study the convergence of the best-fit estimates of the Voigt circuit parameters to their asymptotic values at large enough numbers of Voigt elements (see Fig. 6).

Good visual agreement is observed between the synthetic EIS data and the Voigt circuit impedance at n ≥ 4 , in agreement with Lagonotte et al.'s observations [10]. A Voigt circuit with 6 serially connected RC elements has been used to fit the impedance data in Fig. 5. At increasing values of n, the best-fit estimate τ1 converges to its asymptotic value τ 1 predicted from the set of Eqs. ( 17), ( 23), ( 30) and (31). The rate of convergence is illustrated in Fig. 6 by plotting the decimal logarithm of the absolute value of relative deviation of τ1 from τ 1 as a function of n. Rapid conver- gence is observed. Similar convergence is noted for R1 to R 1 (not shown in the Figure ). In contrast, the convergence rate is slower for the parameter estimates of the second RC element, as illustrated for τ 2 in the same Figure . At small values of n, e.g. n = 6 in Ref. [10], the estimate of the first time constant approaches its asymptotic value with less than 1 ‰ relative deviation, unlike the estimate of the second time constant because of its slower convergence rate. Indeed, the relative deviation of τ 2 from τ 2 is close to 5 % at n = 6. This qualitatively agrees with the experimental observations by Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF]. However, the experimental rate of convergence of τ1 to τ 1 in Ref. [10] is slower than that predicted in Fig. 6 of this article. We should be aware, here, of the weakness of the FLW model to accurately describe the diffusion-convection process near RDE. This point deserves to be further analyzed in a forthcoming article [START_REF] Montella | [END_REF], starting from a more efficient model of the diffusion-convection impedance for RDE.

Please insert Fig. 6 at this place

Next, we look at the dependence of τ1 on the steady-state potential imposed for EIS compu- tation/measurement. The left-side trajectory (blue curve) in the zone diagram of Fig. 4 has been plotted using the same parameter values as in Fig. 5. The time constant τ 1 evaluated along this trajectory is predicted to be nearly constant (with less than 5 % deviation) in the interval ψ ≤ 0.93.

The same condition applies to the converged value of its best-fit estimate τ1 . This should explain the nearly constant value of τ1 obtained by Lagonotte et al. at the three steady-state potentials imposed in Fig. 2 of Ref. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF].

Conclusion

First, the Voigt circuit, with infinitely many RC elements, is an exact representation model Exact representation, at any frequency, of electrochemical impedances involving a finitelength diffusion component requires infinitely many RC elements in the Voigt equivalent circuit.

Fairly accurate approximation of such impedances by the Voigt circuit is less demanding and requires only a few RC elements, as illustrated by the complex nonlinear least squares fits of synthetic EIS data performed in this work. The minimum number ( n min ) of RC elements required for ap- proximation of the electrode impedance for the E reaction on RDE at any frequency depends on the approximation level demanded, and, in addition, on the resistance ratio, r R = R t R d , and the time constant ratio, r τ = R d C dl τ d . As a general rule, the number n min of RC elements should decrease at increasing values of r R and/or r τ . At the limit, the condition ( n min = 1) should apply in the upper

Figures

Fig. 1: Dependence of the time constant ratio (τ 1 τ d ) on the normalized steady-state current (ψ = I I lim ). The time constant ratio is calculated from Eqs. ( 15), ( 17) and ( 20), using the parame- ter values given in the text. Same parameter values as in Fig. 5. , calculated from Eqs. ( 4) and ( 12) at r τ = 0.01.

Fig. A. 1

 1 Fig. A.1 of the Appendix. Hence, the electrode impedance in Eq. (21) and (25) can be realized as a

4 (

 4 ) using the parameter values given in the captions. The frequency range explored satisfies -2 ≤ log f Hz() ≤ ) at a spacing of 9 impedance calculations per decade.The impedance graph mimics, at least roughly, the Nyquist graph plotted by Lagonotte et al. at the equilibrium potential of the Pt electrode (see Fig.3of Ref.[START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF]).The EIS data can be fitted using the impedance model in Eq.[START_REF] Temes | Introduction to Circuit Synthesis and Design[END_REF]. For this purpose, we use noise-free4 unweighted 5 impedance data together with the Levenberg-Marquardt algorithm implemented in the 'NonlinearModelFit' built-in function of Mathematica [22]. The fit residuals are shown in Fig. 5(b).

  for the electrode impedance pertaining to the E reaction investigated on a uniformly accessible electrode under finite-length diffusion conditions for redox species in the electrolyte. This fully validates the conjecture in the seminal work by Agarwal et al.[START_REF] Agarwal | Measurement Models for Electrochemical Impedance Spectroscopy, I. Demonstration and Applicability[END_REF]. The generalized Heaviside method presented in the Appendix provides a general way to derive the analytical expressions of the Voigt circuit parameters (resistances, capacitances and related time constants) as a function of the kinetic rate constant of the electrochemical reaction, the mass transport parameters, the initial concentrations of redox species, the double-layer capacitance at the interface, and the steady-state potential imposed to the electrode. Application to Nernstian redox systems investigated without or with double-layer charging effect, as well as to non-Nernstian redox systems involving slow electron-transfer reaction kinetics, has been presented in this article as illustration examples.

Fig. 2 :

 2 Fig. 2: (a) Real-p impedance and (b) the first three poles (crosses) and three zeros (circles) of the dimensionless impedance, Z p( ) R d = pτ d r τ + 1 r R + pτ d ( ) -1 2 tanh pτ d ( ) 1 2         { } -1, calculated from

Fig. 3 :Fig. 4 :

 34 Fig. 3: Convergence test for the two truncated formulations of the electrode impedance either based on (a) product expansion of the impedance in Eq. (28), or (b) partial fraction expansion of the impedance in Eq. (29). Parameter values: r R = 0.2 and r τ = 0.01. Nyquist graph of the normalized electrode impedance (red curve). Nyquist graphs of the normalized impedance approximations (blue curves) with (a) n = 1, 3, 9, 18 from right to left, and (b) n = 1, 3, 7, 9 from right to left. For interpretation of the references to colors in the figure captions, the reader is referred to the Web version of this article.

Fig. 5 :

 5 Fig. 5: (a) Nyquist plot for the impedance data (dots) calculated from Eqs. (4) and (21), and bestfitting impedance graph (solid line) obtained from a Voigt circuit with 6 RC elements. Parametrization of the impedance graph with log f Hz ( ) . (b) Residual errors (dots) of the complex nonlinear least squares fit for the real parts (data N° 1 to 64) and the imaginary parts (data N° 65 to 128) of the electrode impedance at increasing frequencies from left to right. Parameter values: R t = 12.5 Ω, R d = 210 Ω , δ = 1.8 × 10 -3 cm , D = 6.5 × 10 -6 cm 2 s -1 , and C dl = 4.0 × 10 -6 F .

Fig. 6 :

 6 Fig. 6: Decimal logarithm of the absolute value of relative deviation of the best-fit estimates τ1 (filled circles) and τ 2 (open circles) from the time constants τ 1 and τ 2 , respectively, plotted as a function of the number of RC elements used in the Voigt circuit. Best linear fits (dashed lines).

Fig. A. 1 :

 1 Fig. A.1: (a) Real-p impedance and (b) the first three poles (crosses) and three zeros (circles) of the normalized electrode impedance, Z p( ) R d = pτ d r τ + pτ d ( ) 1 2 coth pτ d ( ) 1 2     { } -1

In practice, Lagonotte et al.[START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF] evaluated the estimates of the poles ( pk for k = 1, 2, ..., n ) of the Voigt circuit impedance, using an appropriate fitting procedure. The time constant estimate τ k is equal to minus the reciprocal of pk .

Proportional and modulus-based weighting strategies yield nearly the same best-fit parameter estimates as unitary weighting for the leading RC elements of the Voigt circuit.

Please insert Fig. 3 at this place

Effect of the steady-state current/potential

We focus again on the first (largest) time constant of the Voigt circuit, τ 1 = τ d q 1 2 , where q 1 is the first (smallest) positive root of Eq. ( 23), which depends both on the time constant ratio r τ in Eq. ( 13) and the resistance ratio r R in Eq. ( 22). Usual impedance calculation rules [20] provide the theoretical formulations of the resistances R t and R d for the E reaction, and, therefore, the theoreti- cal formulations of the dimensionless parameters r R and r τ . For the sake of simplicity, we assume below the same diffusion coefficient for both redox species, D O = D R = D , the same initial concen- tration of redox species, c O * = c R * = c * , and the same symmetry factor for the electron transfer reac- tion in the direction of oxidation (subscript 'o') and reduction (subscript 'r'), α o = α r = 1 2 . The resistance ratio, r R = R t R d , takes on the following simplified formulation where k° is the standard rate constant of the electrochemical reaction, m = D δ is the mass transport constant of redox species under steady-state conditions, and ξ = zf N E -E°'

( ) is the dimensionless potential:

Using the same simplifications as above, the time constant ratio, r τ = R d C dl τ d , can be derived as:

where sech(.) denotes the hyperbolic secant that is the reciprocal of hyperbolic cosine. The first time constant of Voigt circuit can be evaluated from Eqs. ( 17), ( 23), (30) and (31) as a function of the standard rate constant ( k°) of the electrochemical reaction, the mass transport parameters (D and δ ), the initial concentration ( c * ) of redox species, the double-layer capacitance ( C dl ) and the steady-state potential (E) imposed to the electrode through the dimensionless potential ξ .

zones of the zone diagram of Fig. 4 where the Nyquist graph for the electrode impedance approaches a semi-circle graph.

Next, we have been focusing on the first (largest) time constant τ 1 of the Voigt circuit, and on its dependence on the steady-state potential imposed to the electrode for EIS calculation or measurement. Nearly invariance of τ 1 with respect to the electrode potential has been predicted for Nernstian redox systems investigated without double-layer charging effect, as well as for non-Nernstian redox systems when the time constant ratio, R t + R d ( )C dl τ d , takes on a low enough value. This applies in practice to fast or intermediate electron-transfer reaction kinetics when EIS measurements are carried out not too far from the equilibrium potential, because of the linear proportionality of τ 1 with the diffusion-convection time constant in this case, i.e. τ 1 = 4τ d π 2 . Clear- ly, such conditions were satisfied in the experimental work of Lagonotte et al. [START_REF] Lagonotte | Extrapolation to zero frequency with first time constant of ferri-ferrocyanide/Pt interface using Foster equivalent circuit and factorized Modeling[END_REF], which explains the nearly constant value of τ1 observed by these authors under equilibrium as well as nonequilibrium conditions. Other asymptotic formulations are possible for this time constant depending on the electrochemical system considered and the experimental conditions, as summarized in the zone diagram of Fig. 4. Looking at this diagram, all possible sequences of asymptotic formulations/values of the time constant τ 1 , and therefore of the converged value of its best-fit estimate τ1 , can be predicted when the electrode potential moves away from its equilibrium value.

The theoretical derivations and conclusions presented in this article apply to finite-length diffusion of redox species in the electrolytic solution. In addition, they provide a first (rough) approach to the Voigt circuit representation model for the impedance of uniformly accessible rotating disk electrodes. More elaborate models have been proposed in the electrochemical literature [START_REF] Orazem | Electrochemical Impedance Spectroscopy[END_REF][START_REF] Tribollet | The modulated Flow at a Rotating Disk Electrode[END_REF][START_REF] Michel | Diffusion-convection impedance using an efficient analytical approximation of the mass transfer function for a rotating disk[END_REF][START_REF] Diard | Re-examination of the diffusion-convection impedance for a uniformly accessible rotating disk. Computation and accuracy[END_REF] to accurately describe the diffusion-convection impedance for RDE. The impact of such models on the appropriate Voigt circuit representation model for the electrode impedance will be analyzed in a forthcoming article [START_REF] Montella | [END_REF].

Appendix

Using the Heaviside expansion method [START_REF] Angot | Compléments de Mathématiques[END_REF], strictly proper rational functions of the Laplace variable p:

whose only singularities are simple real poles pk, for k = 1, 2, ..., n, can be expanded into partial fractions as:

where the residue of the k-th pole is given by:

Let us consider the class of electrochemical models whose impedance Z p ( ) can be substituted for H p ( ) in Eqs. (A.1) and (A.2). We get the Voigt circuit impedance in Eq. ( 1), with R 0 = 0 :

where the k-th resistance, capacitance and time constant have the following expressions:

Many irrational transfer functions for spatially distributed systems admit a partial-fraction expansion like that in Eqs. (A.1)-(A.5), with infinitely many poles ( n → ∞) [START_REF] Curtain | Transfer Functions of Distributed Parameter Systems: A Tutorial[END_REF]. The electrochem- ical impedances investigated in this article belong to this class. The electrode impedance derived from Eqs. ( 4) and ( 12) takes on the form of Eq. (A.1) where f p ( ) = R d and g p ( ) is the function:

The impedance Z p ( ) has an infinite number of simple real poles p k , which are the roots of the equation:

Those poles should have negative values so that the condition (τ k > 0) applies to the RC elements of the Voigt circuit. Hence, we can write p k τ d ( ) 1 2 = iq k , with real q k , and, therefore:

Using the identity coth iq k ( ) = -icot q k , Eq. (A.7) provides the generating equation for q k :

and the impedance poles follow from Eq. (A.8) for k = 1, 2, …, while the impedance zeros are equal to -k 2 π 2 τ d , irrespective of the r τ value. The locations of poles and zeros in Fig. A.1(b) shows that (i) all poles and zeros are simple and lie on the negative real axis of the complex (p) plane, (ii) the poles and zeros are interlaced, and (iii) the first location on the negative real axis (i.e.

the smallest absolute value) is that of a pole.

Please insert Fig. A.1 at this place

Calculus of residues, from Eq. (A.3), requires the derivative of the function g p ( ). One obtains from Eq. (A.6):

where csch z ( ) = 1 sinh z ( ) denotes the hyperbolic cosecant. Using again p k τ d ( ) 1 2 = iq k , and coth iq k ( ) = -icot q k , together with Eq. (A.9) and the identity csch iq k ( ) = -icsc q k = -i sin q k , the residue of the k-th pole can be derived from Eq. (A.3) as:

From Eq. (A.11), (iv) the residues at the poles of the electrode impedance are real and positive. The properties (i)-(iv) above are those of RC network impedances [1, p. 73]. The relevant expressions for the k-th parameters of the Voigt circuit reported in Eqs. ( 16)-( 18) of the article result from Eqs. (A.5), (A.8) and (A.11). Moreover, setting r τ = 0 in the absence of double-layer charging effect, we get q k = 2k -1 ( )π 2 from Eq. (A.9), p k = -2k -1 ( ) 2 π 2 4τ d ( ) from Eq. (A.8), and res p k ( ) = 2 R d τ from Eq. (A.11), so that the partial-fraction expansion of the finite-length diffusion impedance in Eq. ( 8) is recovered using Eqs. (A.4) and (A.5).

Declaration of competing interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Claude Montella: Conceptualization, Methodology, Formal analysis, Visualization, Writing original draft.


Gérard Gouesbet
email: gouesbet@coria.fr

Leonardo André

Axicon terms associated with gradient optical forces in generalized Lorenz-Mie theory

Keywords: generalized Lorenz-Mie theory, Bessel beams, optical forces, scattering forces, gradient forces, axicon forces. 1 Introduction

In a recent paper devoted to the study of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in generalized Lorenz-Mie theory, new kind of optical forces are discovered. They are di¤erent than the well established scattering and gradient forces, and called axicon forces. In the present paper, within the same framework, we demonstrate the existence of a new kind of axicon terms, associated with gradient forces, which are zero both when the axicon angle is zero and when an on-axis con…guration is considered.

manipulation. In [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF], beside classical scattering and gradient forces, a new kind of optical forces associated with scattering forces, named axicon (scattering) forces, has been discovered. Concerning gradient forces, they have been found to be zero in agreement with the fact the longitudinal gradients of (i) the electric intensity (ii) the energy density and (iii) the longitudinal component of the Poynting vector are zero for Bessel beams. In the present paper, we proceed to a re…ned investigation of the gradient forces identi…ed in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] and demonstrate that they are actually accompanied by axicon terms which are zero if the axicon angle is zero and/or if the con…guration is an on-axis con…guration.

Before proceeding further, it is important to insist on the meaning of new expressions introduced in the present paper and in the previous paper [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. The new terminology introduces "axicon forces" and "axicon terms". Axicon forces are forces which are deduced from the expressions for pressure radiation crosssections (or forces) but which have no counterpart in expressions dealing with scattering expressions, e.g. such as expressed by the Poynting vector, which may be called axicon scattering forces as examined in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. Conversely, axicon terms are extra-terms which do not arise from the expressions for pressure radiation cross-sections (or forces) but which arise from gradient expressions. Axicon terms associated with gradient expressions are examined in the present paper. . The paper is organized as follows. Section 2 presents a background already presented in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] but which is repeated here for the convenience of the reader, complemented by recalling a few results taken from [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] but useful for the present paper. Section 3 deals with the longitudinal gradient of the electric intensity. Section 4 deals with the longitudinal gradient of the time averaged energy density. Section 5 deals with the longitudinal gradient of the longitudinal component of the Poynting vector. In these sections 3, 4, and 5, it is demonstrated that some axicon terms must be associated with the gradient force. Section 6 provides an extra-discussion of some issues enlightening the results obtained in the present paper and in the previous paper [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. Section 7 is a conclusion.

Physical and mathematical background.

The GLMT stricto sensu describes the interaction between an arbitrary (structured) shaped beam and a homogeneous spherical particle de…ned by its diameter and its complex refractive index, e.g. [START_REF] Gouesbet | Light scattering from a sphere arbitrarily located in a Gaussian beam, using a Bromwich formulation[END_REF], [START_REF] Gouesbet | Combustion measurements[END_REF], [START_REF] Gouesbet | Generalized Lorenz-Mie theories, second edition[END_REF], and references therein. Also, see reviews [START_REF] Gouesbet | Latest achievements in generalized Lorenz-Mie theories: A commented reference database[END_REF], [START_REF] Gouesbet | T-matrix methods for electromagnetic structured beams: A commented reference database for the period 2014-2018[END_REF], in particular for other GLMTs devoted to the interaction between laser beams and other kinds of scatterers. In the present paper, we consider a Rayleigh particle located at the origin O P of a Cartesian coordinate system O P xyz illuminated by an o¤-axis Bessel beam propagating along the z-direction.

Optical forces.

In the GLMT framework, we express the optical force components F i by using (unnormalized) pressure radiation cross-section components C pr;i (i = x; y; z). In the present paper, we deal only with the longitudinal cross-section (along z direction) reading as:

C pr;z = E 0 H 0 2 2 1 X n=1 n X m= n f 1 (n + 1) 2 (n + 1 + jmj)! (n jmj)! (1)
Re

[(a n + a n+1 2a n a n+1 )g m n;T M g m n+1;T M +(b n + b n+1 2b n b n+1 )g m n;T E g m n+1;T E ] +m 2n + 1 n 2 (n + 1) 2 (n + jmj)! (n jmj)! Re[i(2a n b n a n b n )g m n;T M g m n;T E )]g
with notations being the ones of [START_REF] Gouesbet | Light scattering from a sphere arbitrarily located in a Gaussian beam, using a Bromwich formulation[END_REF]. However, note that, strictly speaking, C pr;z of Eq.1 would be a pressure radiation cross-section, i.e. would be homogeneous to an area, if we used the normalization condition E 0 H 0 =2 =1 (which was used in original works, e.g. Eq.(3.106) in [START_REF] Gouesbet | Generalized Lorenz-Mie theories, second edition[END_REF]) . Here, we did not use this condition for the sake of consistency with later expressions concerning the Poynting vector. When using the normalization condition, the forces are related to the cross-sections according to F i = C pr;i =c ( [START_REF] Van De Hulst | Light scattering by small particles[END_REF], p.14). The time-dependence of the wave is exp(+i!t) which is the usual choice in GLMT. Also, note that the denominator 1=(2 ) was missprinted to 1=(4 ) in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF].

Furthermore, speci…cally, a n and b n are the usual Mie coe¢ cients of the usual Lorenz-Mie theory, g m n;T M and g m ;;T E , with TM and TE standing for "Transverse Magnetic" and "Transverse Electric" respectively, are the beam shape coe¢cients (BSCs) encoding the description of the beam, either in terms of scalar potentials, e.g. [START_REF] Gouesbet | Generalized Lorenz-Mie theories, second edition[END_REF] or in terms of vector spherical wave functions [START_REF] Gouesbet | T-matrix formulation and generalized Lorenz-Mie theories in spherical coordinates[END_REF], the star denotes a complex conjugation and is the wavelength. In the present paper, we shall consider only the longitudinal force of Eq.1 because it will lead to simpler computations than for the two other force components along x and y, and, more important, because we expect a natural separation between gradient and scattering forces along the z-direction.

2.2

Poynting vector.

For the same reasons, we shall only consider the z component S z of the Poynting vector de…ned as:

S z = 1 2 Re(E x H y E y H x ) (2)
which has already been evaluated to [START_REF] Gouesbet | Partial wave expansions and properties of axisymmetric light beams[END_REF], [START_REF] Gouesbet | Poynting theorem in terms of beam shape coe¢ cients and applications to axisymmetric, dark and non-dark, vortex and nonvortex beams[END_REF]:

S z = E 0 H 0 2r 2 Re 1 X n=1 +n X m= n 1 X p=1 +p X q= p ic pw n c pw p e i(m q)' (3)
(sin S mq np + cos C mq np )
in which (r; ; ') are spherical coordinates attached to the Cartesian coordinates (x; y; z), c pw n are prefactors appearing in the Bromwich formulation of the Lorenz-Mie theory [START_REF] Gouesbet | Sur la généralisation de la théorie de Lorenz-Mie[END_REF], reading as:

c pw n = 1 ik ( i) n 2n + 1 n(n + 1) (4) and :
S

+g m n;T M g q p;T E 0 n 0 p (m jmj n jqj p + q jqj p jmj n ) g m n;T E g q p;T M p n (m jmj n jqj p + q jqj p jmj n ) +g m n;T E g q p;T E n 0 p (mq jmj n jqj p + jmj n jqj p ) (6)
in which k is the wavenumber, P = jmj n (cos ) are generalized Legendre functions, and n = n (kr) are Ricatti-Bessel functions, while a prime denotes a derivative with respect to the argument (and a double prime a second derivative). In Eqs.5-6, arguments are omitted for convenience. The omission of arguments will be recurrent in the sequel whenever we …nd it convenient.

Mathematical background.

Associated Legendre functions are de…ned according to Hobson's notation [START_REF] Robin | Fonctions sphériques de Legendre et fonctions sphéroidales[END_REF]:

P jmj n (cos ) = ( 1) jmj (sin ) jmj d jmj P n (cos ) (d cos ) jmj (7)
in which P n (cos ) are the usual Legendre polynomials. Generalized Legendre functions may then be evaluated according to:

m n (cos ) = dP m n (cos ) d (8) m n (cos ) = P m n (cos ) sin (9)
Ricatti-Bessel functions n (kr) may be expressed in terms of spherical Bessel functions j n (kr) according to:

n (x) = xj n (x) (10)
Furthermore, we have [START_REF] Arfken | Mathematical methods for physicists[END_REF]:

j 1 (x) = sin x x 2 cos x x (11)
2.4 Bessel beams.

Bessel beams have been introduced by Durnin and co-workers [START_REF] Durnin | Di¤raction-free beams[END_REF], [START_REF] Durnin | Exact solutions for nondi¤racting beams. I. The scalar theory[END_REF]. They possess the appealing property of being self-healing and non di¤racting, and, more important in the context of the present paper, they furthermore possess a propagation invariance property, namely the intensity of the beam is constant along the direction of propagation. There exist an in…nity of kinds of Bessel beams, in particular depending on the value given to an arbitrary function g( 0 ), in which 0 is the axicon angle (or half-cone angle) of the beam, with di¤erent linear and circular polarizations [START_REF] Lock | Angular spectrum and localized model of Davis-type beam[END_REF], [START_REF] Wang | General description of circularly symmetric Bessel beams of arbitrary order[END_REF], [START_REF] Wang | General description of transverse mode Bessel beams and construction of basis Bessel …elds[END_REF], [START_REF] Wang | Multipole expansion of circularly Bessel beams of arbitrary order for scattering calculations[END_REF]. All of them generically exhibit the following structure:

K i = k i (x; y)e ikzz , i = x; y; z, K = E or H (12)
in which k z = k cos 0 is the longitudinal wavenumber. Among all of them, we shall pay a particular attention to circularly symmetric Bessel beams of arbitrary order whose BSCs in an o¤-axis con…guration read as [START_REF] Wang | Multipole expansion of circularly Bessel beams of arbitrary order for scattering calculations[END_REF]:

g m n;T M = g( 0 )e ikzz0 [i l m+1 e i(l m+1) 0 A m n + i l m 1 e i(l m 1) 0 B m n ] ( 13
)
g m n;T E = ig( 0 )e ikzz0 [i l m+1 e i(l m+1) 0 A m n i l m 1 e i(l m 1) 0 B m n ] (14)
in which:

A m n = ( 1) (m jmj)=2 (n m)! (n + jmj)! J l m+1 ( 0 )[ m n (cos 0 ) + m m n (cos 0 )] (15)
B m n = ( 1) (m jmj)=2 (n m)! (n + jmj)! J l m 1 ( 0 )[ m n (cos 0 ) m m n (cos 0 )] (16)
in which l is the order of the Bessel beam under consideration, 0 = k t 0 , k t = k sin 0 is the transverse wavenumber, 0 = (x 2 0 +y 2 0 ) 1=2 , 0 = tan 1 (y 0 =x 0 ), (x 0 ; y 0 ; z 0 ) denotes the location of the beam origin with respect to the origin of the particle coordinate system at which the scatterer is located, and therefore speci…es the o¤-axis location of the beam, and J k (:) is the k-order Bessel function of the …rst kind. When g( 0 ) = (1+ cos 0 )=4, the beam reduces to a Davis circularly symmetric Bessel beam as discussed in [START_REF] Lock | Angular spectrum and localized model of Davis-type beam[END_REF], [START_REF] Mishra | A vector wave analysis of a Bessel beam[END_REF]. When g( 0 ) = 1=2, the beam reduces to another kind of Bessel beam discussed in [START_REF] Cizmar | Sub-micron particle organization by self-imaging of non-di¤racting beams[END_REF], [START_REF] Taylor | Multipole expansion of Bessel and Gaussian beams for Mie scattering calculations[END_REF], [START_REF] Chen | Analytical partial wave expansion of vector Bessel beam and its application to optical binding[END_REF].

2.5

Summary of some previous results.

In this subsection, we recall a few results obtained in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] related to gradient forces. In the Rayleigh regime which was considered in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] and which is still considered in the present paper, only the (n = 1)-partial waves are to be considered and we only have to retain, in the expansions, the Mie coe¢ cient a 1 which simpli…es to :

a 1 = 2i 3 m 2 1 m 2 + 2 3 + O(i 5 ) + 4 9 ( m 2 1 m 2 + 2 ) 2 6 (17) b 1 = O(i 5 ) (18)
in which m is the refractive index (here taken to be real) with respect to the surrounding medium and the size parameter d= . The other coe¢ cients a n and b n (n > 1) involves still higher powers of : Real parts are then proportional to 6 and imaginary parts are proportional to 3 while higher powers are discarded and we therefore retain only:

Im(a 1 ) = 2 3 m 2 1 m 2 + 2 3 (19) Re(a 1 ) = 4 9 ( m 2 1 m 2 + 2 ) 2 6 (20)
The scattering force represented by C s pr;z and the gradient force represented by C g pr;z are proportional to 6 and 3 respectively. It is then demonstrated that, by evaluating C g pr;z with the proviso that the Rayleigh particle is only sensitive to the (n = 1)-partial waves, we have:

C g pr;z = 0 (21)
The nullity of this term identi…ed as a gradient term is consistent with the fact that, due to Eq.12, the longitudinal gradients of the electric intensity E:E , of the time-averaged energy density (" 0 E:E + 0 H:H ) and of the z component of the Poynting vector are zero. However, the Rayleigh particle does not perceive all the partial waves involved in Eq.12 but only the (n = 1)partial waves. This implies that a re…nement of the analysis is required in which the electric intensity, the time-averaged energy density and the z component of the Poynting vector would be evaluated by only considering the contributions of the (n = 1)-partial waves. The aim of the present paper is to proceed to this re…ned analysis. Quantities restricted to (n = 1) partial waves are decorated with a tilde.

Electric intensity.

When restricted to the (n = 1)-partial waves, the electric components in spherical coordinates read as (e.g. [START_REF] Gouesbet | Generalized Lorenz-Mie theories, second edition[END_REF], pp. [START_REF] Gouesbet | Finite series expressions to evaluate the beam shape coe¢ cients of a Laguerre-Gauss beam focused by a lens in an on-axis con…guration[END_REF][56]:

f E r = kE 0 c pw 1 ( 00 1 + 1 ) +1 X m= 1 g m 1;T M P jmj 1 e im' (22)
f E = E 0 r c pw 1 +1 X m= 1 [g m 1;T M 0 1 jmj 1 + mg m 1;T E 1 jmj 1 ]e im' (23)
f E ' = iE 0 r c pw 1 +1 X m= 1 [mg m 1;T M 0 1 jmj 1 + g m 1;T E 1 jmj 1 ]e im' (24)
From these equations, we may evaluate:

Êr E r = 9 4 jE 0 j 2 k 2 k 2 ( 00 1 + 1 ) 2 [e 2i' g 1 1;T M g 1 1;T M P 1 1 P 1 1 ( 25
)
+e i' (g 1 1;T M g 0 1;T M + g 0 1;T M g 1 1;T M )P 1 1 P 0 1 +e 0' (g 1 1;T M g 1 1;T M P 1 1 P 1 1 + g 0 1;T M g 0 1;T M P 0 1 P 0 1 + g 1 1;T M g 1 1;T M P 1 1 P 1 1 ] +e i' (g 0 1;T M g 1 1;T M + g 1 1;T M g 0 1;T M )P 1 1 P 0 1 +e 2i' g 1 1;T M g 1 1;T M P 1 1 P 1 1 ]
in which we have used:

c pw 1 c pw 1 = 9 4k 2 (26)
and:

Ê E = 9 4 jE 0 j 2 k 2 f 0 1 0 1 r 2 [e 2i' g 1 1;T M g 1 1;T M 1 1 1 1 ( 27
)
+e i' (g 1 1;T M g 0 1;T M + g 0 1;T M g 1 1;T M ) 1 1 0 1 +e 0' (g 1 1;T M g 1 1;T M 1 1 1 1 + g 0 1;T M g 0 1;T M 0 1 0 1 + g 1 1;T M g 1 1;T M 1 1 1 1 ) +e i' (g 1 1;T M g 0 1;T M + g 0 1;T M g 1 1;T M ) 1 1 0 1 +e 2i' g 1 1;T M g 1 1;T M 1 1 1 1 ] + 0 1 1 r 2 [e 2i' (g 1 1;T M g 1 1;T E g 1 1;T E g 1 1;T M ) 1 1 1 1 +e i' (g 0 1;T M g 1 1;T E g 1 1;T E g 0 1;T M ) 1 1 0 1 +e 0' (g 1 1;T E g 1 1;T M g 1 1;T M g 1 1;T E + g 1 1;T M g 1 1;T E g 1 1;T E g 1 1;T M ) 1 1 1 1 +e i' (g 1 1;T E g 0 1;T M g 0 1;T M g 1 1;T E ) 1 1 0 1 +e 2i' (g 1 1;T E g 1 1;T M g 1 1;T M g 1 1;T E ) 1 1 1 1 1 1 r 2 [e 2i' g 1 1;T E g 1 1;T E e 0' (g 1 1;T E g 1 1;T E + g 1 1;T E g 1 1;T E ) + e 2i' g 1 1;T E g 1 1;T E ] 1 1 1 1 Ê' E ' = 9 4 jE 0 j 2 k 2 f 0 1 0 1 r 2 [ e 2i' g 1 1;T M g 1 1;T M 1 1 1 1 (28) +e 0' (g 1 1;T M g 1 1;T M + g 1 1;T M g 1 1;T M ) 1 1 1 1 e 2i' g 1 1;T M g 1 1;T M 1 1 1 1 ] + 0 1 1 r 2 [e 2i' (g 1 1;T E g 1 1;T M g 1 1;T M g 1 1;T E ) 1 1 1 1 +e i' (g 0 1;T E g 1 1;T M g 1 1;T M g 0 1;T E ) 1 1 0 1 +e 0' (g 1 1;T E g 1 1;T M g 1 1;T M g 1 1;T E + g 1 1;T E g 1 1;T M g 1 1;T E g 1 1;T M ) 1 1 1 1 +e i' (g 1 1;T M g 0 1;T E g 0 1;T E g 1 1;T M ) 1 1 0 1 +e 2i' (g 1 1;T M g 1 1;T E g 1 1;T E g 1 1;T M ) 1 1 1 1 + 1 1 r 2 [e 2i' g 1 1;T E g 1 1;T E 1 1 1 1 + e i' (g 1 1;T E g 0 1;T E + g 0 1;T E g 1 1;T E ) 1 1 0 1 +e 0' (g 1 1;T E g 1 1;T E 1 1 1 1 + g 0 1;T E g 0 1;T E 0 1 0 1 + g 1 1;T E g 1 1;T E 1 1 1 1 ) +e i' (g 0 1;T E g 1 1;T E + g 1 1;T E g 0 1;T E ) 1 1 0 1 ) + e 2i' g 1 1;T E g 1 1;T E 1 1 1 1 ]
We may then evaluate the corresponding electric intensity according to:

f I E = Êr E r + Ê E + Ê' E ' (29)
which is found to read as:

f I E = +2 X k= 2 f I k E ( 30
)
Specifying the values of associated Legendre functions and of generalized Legendre functions according to:

P 0 0 = P 0 0 (cos ) = 1 (31)
P 0 1 = P 0 1 (cos ) = 1 1 (cos ) = cos (32)
P 1 1 = P 1 1 (cos ) = 0 1 (cos ) = sin (33)
1 1 = 1 1 (cos ) = 1 (34)
we obtain:

g I 2 E = f I 2 E = 9 4 jE 0 j 2 k 2 e 2i' sin 2 f[k 2 ( 00 1 + 1 ) 2 0 1 0 1 r 2 ]g 1 1;T M g 1 1;T M + 1 1 r 2 g 1 1;T E g 1 1;T E g (35) g I 1 E = f I 1 E = 9 4 jE 0 j 2 k 2 e i' sin fcos [ 0 1 0 1 r 2 k 2 ( 00 1 + 1 ) 2 ](g 1 1;T M g 0 1;T M + g 0 1;T M g 1 1;T M ) + 0 1 1 r 2 (g 0 1;T M g 1 1;T E g 1 1;T E g 0 1;T M g 1 1;T M g 0 1;T E + g 0 1;T E g 1 1;T M ) (36)
+ cos 1 1 r 2 (g 1 1;T E g 0 1;T E + g 0 1;T E g 1 1;T E )g f I 0 E = 9 4 jE 0 j 2 k 2 e 0' f[k 2 sin 2 ( 00 1 + 1 ) 2 + cos 2 0 1 0 1 r 2 ](g 1 1;T M g 1 1;T M + g 1 1;T M g 1 1;T M ) +[k 2 cos 2 ( 00 1 + 1 ) 2 + sin 2 0 1 0 1 r 2 ]g 0 1;T M g 0 1;T M ( 37
)
2 cos 0 1 1 r 2 (g 1 1;T M g 1 1;T E g 1 1;T M g 1 1;T E + g 1 1;T E g 1 1;T M g 1 1;T E g 1 1;T M ) (cos 2 + 1) 1 1 r 2 (g 1 1;T E g 1 1;T E + g 1 1;T E g 1 1;T E ) + sin 2 1 1 r 2 g 0 1;T E g 0 1;T E + 0 1 0 1 r 2 (g 1 1;T M g 1 1;T M + g 1 1;T M g 1 1;T M )g
We are now interested in the calculation of the derivative of the (restricted) electric intensity with respect to z, taken at z = 0, i.e. at the place where the particle is located. This requires in a …rst step to express cos , sin and the radial coordinate r in terms of Cartesian coordinates, a trivial task indeed. Furthermore, in evaluating these derivatives, and other derivatives in subsequent subsections, derivatives with respect to ' may be omitted. In a few words, this comes from the fact that z dependence is only encoded in the spherical coordinates r and , but not in ': The suspicious reader may have a direct check by using the equation @f (r; ) exp( ik')=@z = exp( ik')@f (r; )=@z + f (r; )@ exp( ik')=@z, leaving the derivatives evaluated below unchanged.

Then, we consider g I 2 E and f I 2 E , and evaluate:

f @ @z sin 2 [k 2 ( 00 1 + 1 ) 2 0 1 0 1 r 2 ]g z=0 = 0 (38)
( @ @z sin 2 1 1 r 2 ) z=0 = 0 (39)
in which we expressed the terms dependent on r and in terms of x; y; z and recalled that, using Eqs.10 and 11:

1 (kr) = kr[ sin(kr) (kr) 2 cos(kr) kr ] (40)
Hence:

( @ g I 2 E @z ) z=0 = ( @ f I 2 E @z ) z=0 = 0 (41)
For g I 1 E and g I 1 E , we have:

1 = f @ @z sin cos [ 0 1 0 1 r 2 k 2 ( 00 1 + 1 ) 2 ]g z=0 6 = 0 (42)
( @ @z sin 0 1 1 r 2 ) z=0 = 0 (43) 2 = ( @ @z sin cos 1 1 r 2 ) z=0 6 = 0 (44)
so that we obtain:

( @ g I 1 E @z ) z=0 = [( @ f I 1 E @z ) z=0 ] = 9 4 jE 0 j 2 k 2 e i' ( 1 + 2 ) (45)
in which (not to be confused with the size parameter) and are given by:

= g 1 1;T M g 0 1;T M + g 0 1;T M g 1 1;T M (46) = g 1 1;T E g 0 1;T E + g 0 1;T E g 1 1;T E (47)
For (@ f I 0 E =@z) z=0 , all derivatives with respect to z, taken at z = 0, of the di¤erent terms involved in Eq.37 are found to be 0, excepted one of them, namely:

3 = ( @ @z cos 0 1 1 r 2 ) z=0 6 = 0 (48)
leading to:

( @ f I 0 E @z ) z=0 = 9 2 jE 0 j 2 k 2 3 ( 49
)
in which:

B 1 1 = J l ( 1 1 + 1 1 ) = 1 2 J l (C + 1) (58)
B 0 1 = J l 1 0 1 = J l 1 S (59)
B 1 1 = 1 2 J l 2 ( 1 1 1 1 ) = 1 2 J l 2 (C 1) (60)
in which C = cos 0 , S = sin 0 and in which we have used

1 1 = C=2, 0 1 = S, 1 1 = C, 1 1 = 1=2 and 1 1 = 1.
Furthermore the argument 0 of the Bessel functions is omitted. We then obtain:

= iS 2 jg( 0 j 2 fe i 0 [(C 1)(J l+1 J l+2 + J l 2 J l 1 ) + (C + 1)(J l 1 J l + J l J l+1 )] e i 0 (C + 1)(J l J l+1 + J l 1 J l ) (61)
e 3i 0 (C 1)(J l 1 J l+2 + J l 2 J l+1 )g = iS 2 jg( 0 j 2 fe i 0 [(C 1)(J l+1 J l+2 + J l 2 J l 1 ) + (C + 1)(J l 1 J l + J l J l+1 )] +e i 0 (C + 1)(J l J l+1 + J l 1 J l ) ( 62
)
+e 3i 0 (C 1)(J l 1 J l+2 + J l 2 J l+1 )g = jg( 0 )j 2 (C 1)(C + 1)(J l 2 J l J l J l+2 ) sin 2 0 (63)
Therefore, although (@I E =@z) z=0 = 0, it happens that (@ e I E =@z) z=0 6 = 0 and contains extra axicon terms. However these extra-terms are zero whevener (i) the axicon angle is zero, i.e. C = 1 and S = 0 (ii) the con…guration is an on-axis con…guration. The fact that the extra-terms are zero for an on-axis con-…guration results from the mathematical property that J k (0) = 0 excepted for k = 0, and from the physical property that 0 = 0 for an on-axis con…guration.

Energy density.

To deal with the energy density "0 E:E + 0 H:H , we have to evaluate the magnetic intensity f I H = Ĥ:H . Instead of Eqs.22-24, we now have to deal with (see again [START_REF] Gouesbet | Generalized Lorenz-Mie theories, second edition[END_REF], pp.55-56):

f H r = kH 0 c pw 1 ( 00 1 + 1 ) +1 X m= 1 g m 1;T E P jmj 1 e im' (64)
f H = H 0 r c pw 1 +1 X m= 1 [g m 1;T E 0 1 jmj 1 mg m 1;T M 1 jmj 1 ]e im' ( 65
) f H ' = iH 0 r c pw 1 +1 X m= 1 [mg m 1;T E 0 1 jmj 1 g m 1;T M 1 jmj 1 ]e im' (66)
We then observe that the magnetic intensity is deduced directly from the electric intensity by carrying out the following changes:

E 0 ! H 0 , g m 1;T M ! g m 1;T E , g m 1;T E ! g m 1;T M ( 67
)
so that we directly obtain:

( @ g I 2 H @z ) z=0 = ( @ f I 2 H @z ) z=0 = 0 (68) ( @ g I 1 H @z ) z=0 = [( @ f I 1 H @z ) z=0 ] = 9 4 jH 0 j 2 k 2 e i' ( 1 + 2 ) (69) ( @ f I 0 H @z ) z=0 = 9 2 jH 0 j 2 k 2 3 (70)
Therefore, although (@I H =@z) z=0 = 0, it happens that (@ e I H =@z) z=0 6 = 0 and contains extra axicon terms. However, once again, these extra-terms are zero whevener (i) the axicon angle is zero, i.e. C = 1 and S = 0 (ii) the con…guration is an on-axis con…guration. Because this is the same conclusion that for the electric intensity, it immediately extends to the time averaged energy density "0 E:E + 0 H:H . More speci…cally, adding " 0 E:E and ^ 0 H:H , and using H 0 =E 0 = p " 0 = 0 , we may establish that, after simplifying and rearranging:

[ @ @z ( "0 E:E + 0 H:H )] z=0 = 2 Re( e i' ) 4 ( 71
)
in which:

4 = 1 + 2 = f @ @z sin cos [ ( 1 ) 2 + ( 0 1 ) 2 r 2 k 2 ( 00 1 + 1 ) 2 ]g z=0 6 = 0 (72)
and:

= + = iS jg( 0 )j 2 e i 0 [(C 1)(J l+1 J l+2 +J l 2 J l 1 )+(C+1)((J l 1 J l )+J l J l+1 )] (73)
to which the above comments still apply.

5 Poynting vector.

Finally, we have to evaluate and discuss (@ f S z =@z) z=0 . We start from Eq.( 3) in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] and specify it to the Rayleigh situation when only the partial waves with n = 1 are to be retained, that is to say we now consider:

( @ f S z @z ) z=0 = 9E 0 H 0 8k 2 Re +1 X m= 1 +1 X q= 1 ie i(m q)' [ @ @z 1 r 2 (sin S mq 11 + cos C mq 11 )] z=0 (74)
in which we have used Eq.26. The coe¢ cients S mq 11 and C mq 11 are given by Eqs.5-6. Each of them contains four terms. Therefore we have 9x2x4=72 derivatives to evaluate with the fortunate fact, however, that there are many redundant terms. Once these derivatives are done, the only nonzero derivatives are found to be:

[ @ @z 1 r sin 1 ( 00 1 + 1 )P 1 1 1 1 ] z=0 = [ @ @z 1 r sin 2 cos 1 ( 00 1 + 1 )] z=0 (75) [ @ @z 1 r sin 0 1 ( 00 1 + 1 )P 1 1 0 1 ] z=0 = [ @ @z 1 r sin cos 0 1 ( 00 1 + 1 )] z=0 (76) [ @ @z 1 r sin 0 1 ( 00 1 + 1 )P 0 1 1 1 ] z=0 = [ @ @z 1 r sin cos 0 1 ( 00 1 + 1 )] z=0 (77) [ @ @z 1 r 2 cos 1 0 1 )( 1 1 1 1 + 1 1 1 1 )] z=0 = [ @ @z 1 r 2 cos (cos 2 + 1) 1 0 1 ] z=0 (78) [ @ @z 1 r 2 cos 1 0 1 )( 1 1 1 1 1 1 1 1 )] z=0 = [ @ @z 1 r 2 cos (cos 2 1) 1 0 1 ] z=0 (79) [ @ @z 1 r 2 cos 0 1 0 1 1 1 0 1 ] z=0 = [ @ @z 1 r 2 sin cos 0 1 0 1 ] z=0 (80)
[ @ @z 1 r 2 cos 1 1 1 1 0 1 ] z=0 = [ @ @z 1 r 2 sin cos 1 1 ] z=0 (81)
[ @ @z 1 r 2 cos 1 0 1 0 1 0 1 ] z=0 = [ @ @z 1 r 2 sin 2 cos 1 0 1 ] z=0 (82)
These quantities appear in S mq 11 with mq = 1 1, 11; 11, 1 1, 10, 10, 0 1, 01, 00 and in C mq 11 with the same values of mq. Inserting these nonzero values in Eq.74 and rearranging we obtain, after a fairly lengthy and tedious computation:

( @ f S z @z ) z=0 = 9E 0 H 0 8k 2 Re if(cos 'G 1 + i sin 'G 2 ) (83)
[ @ @z 1 r 2 sin cos (kr

0 1 ( 00 1 + 1 ) 0 1 0 1 )] z=0 +(cos 'G 1 i sin 'G 2 ) [ @ @z 1 r 2 sin cos 1 1 ] z=0
in which we used the fact that Re i(z + z ) = 0, and in which:

G 1 = g 0 1;T E (g 1 1;T M g 1 1;T M ) + g 0 1;T M (g 1 1;T E g 1 1;T E ) (84)
G 2 = g 0 1;T E (g 1 1;T M + g 1 1;T M ) + g 0 1;T M (g 1 1;T E + g 1 1;T E ) (85)
G 1 and G 2 are now evaluated using the now usual manner, leading to:

G 1 = jg( 0 )j 2 (G 1 1 e i 0 G 1 1 e i 0 G 3 1 e 3i 0 G 3 1 e 3i 0 ) (86)
in which:

G 1 1 = S[ C 1 2 (J l+1 J l+2 + J l 2 J l 1 ) (C + 1)J l 1 J l ] (87)
G 1 1 = S[ C 1 2 (J l+1 J l+2 + J l 2 J l 1 ) (C + 1)J l J l+1 ] (88)
G 3 1 = G 3 1 = S(C 1) 2 (J l 1 J l+2 J l 2 J l+1 ) (89)
Therefore, G 1 = 0 if (i) the axicon angle is zero (S = 0, C = 1) and (ii) if the o¤-axis con…guration becomes an on-axis con…guration. As for G 2 , we …nd:

G 2 = jg( 0 )j 2 (G 1 2 e i 0 G 1 2 e i 0 G 3 2 e 3i 0 + G 3 2 e 3i 0 ) (90)
in which:

G 1 2 = G 1 1 , G 1 2 = G 1 1 , G 3 2 = G 3 1 , G 3 2 = G 3 1 (91)
which is commented similarly as G 1 . Therefore, as a whole, although (@S z =@z) z=0 is equal to zero, it happens that (@ f S z =@z) z=0 is in general not equal to 0, although it is zero (i) if the axicon angle is zero (ii) for an on-axis con…guration.

6 Discussions related to axicon forces and axicon terms.

The contents of this section arises from reviewer comments and from mail exchanges with Pr. James Lock from Cleveland University, USA.

6.1

Relationship with the standard theory of scattering and optical forces.

The previous paper [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] introduces new axicon forces which have no counterpart in the usual standard theory of optical forces expressed in terms of scattering and of gradient forces. We begin with axicon forces discussed in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. We may then observe that the newly discussed axicon forces indeed arise from higher-order multipole contributions. Such high-order contributions arise from the structure of Eq.( 1) in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] which expresses the longitudinal radiation pressure cross-section C pr;z . BSCs then occur in a coupled way in terms reading as g m n;T M g m n+1;T M , g m n;T E g m n+1;T E and g m n;T M g m n;T E (m = n; :::; +n). In the Rayleigh limit of the GLMT, we only retain (n = 1)-partial waves so that we have to deal with BSCs coupling terms of the form g m 1;T M g m 2;T M , g m 1;T E g m 2;T E and g m 1;T M g m 1;T E (m = 1; 0; +1). The standard scattering force C s11 pr;z is expressed by Eq.(67) of [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] containing couplings of the form g 1 1;T M g 1 1;T E and g 1 1;T M g 1 1;T E which therefore do not involve any higher-order multipole contributions. Conversely, the new axicon forces C s12 pr;z and C s0 pr;z are expressed by Eqs.( 49) and ( 50) of [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF], which are expressed in terms of quantities denoted G 12 and G 0 given by Eqs.( 45) and ( 46) of [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. We then observe that G 12 involves couplings of the form g 1 1;T M g 1 2;T M and g 1 1;T M g 1 2;T M while G 0 involves g 0 1;T M g 0 2;T M , that is to say involves coupling with higher-order multipole contributions, namely between (n = 1) and (n = 2) multipoles. In other words, the Rayleigh limit of the GLMT does not exactly identify with the Rayleigh dipole theory of optical forces.

6.2

Relationship with the localization principle of van de Hulst.

A plane wave of in…nite extent may be thought of as being made up of separate rays of light each independently pursuing its own path. According to the van de Hulst principle of localization, e.g. p.208 of [START_REF] Van De Hulst | Light scattering by small particles[END_REF], a partial wave of order n in the electromagnetic description of the plane wave corresponds to a ray passing at a distance from the origin equal to (n + 1=2)( =2 ). Here, the partial wave of order n refers to a term in the classical Lorenz-Mie theory which contains Bessel functions and spherical harmonics, with the integer n ranging from 1 to 1. This van de Hulst principle of localization has been successfully extended from the case of plane waves to the case of structured beams and has been at the origin of the localized approximation to the evaluation of BSCs, e.g. [START_REF] Gréhan | Scattering of laser beams by Mie scatter centers: numerical results using a localized approximation[END_REF], [START_REF] Maheu | Generalized Lorenz-Mie theory: …rst exact values and comparisons with the localized approximation[END_REF], see [START_REF] Gouesbet | Generalized Lorenz-Mie theories and description of electromagnetic arbitrary shaped beams: localized approximations and localized beam models, a review[END_REF] for a review.

Let us …rst consider an on-axis location of the particle or, equivalently, in the present context, a coordinate system centered on the beam, and let us call z on its axis which is as well the direction of propagation of the beam (and which is parallel to the z-axis considered in the o¤-axis con…guration studied in the present paper). The GLMT expands the on-axis wave in terms of partial waves characterized by two subscripts n on (the order of the partial wave) and m on (the azimuthal mode). The van de Hulst principle then associates the partial wave n on with a ray passing at a distance (n on + 1=2)( =2 ) from the axis z on . Note however that the value of this distance is questionable because its derivation, using the method of stationary phase -not available from van de Hulst's book -requires that n >> 1 (James Lock, private communication), but this feature does not change the gist of the argument discussed in this subsection. In the Rayleigh limit, only the n on = 1 partial wave interacts with the particle because only the Mie coe¢ cient a 1 is to be retained. This corresponds to a ray passing at a distance equal to d on = 3 =(4 )

=4 from the z on -axis. The o¤-axis coordinate system is the one used in the present paper, with Cartesian coordinates (x; y; z) and partial waves characterized by the order n of f and the azimuthal mode m of f (more simply denoted as n and m in the present paper). The o¤-axis description of the beam is more complicated than the on-axis description of the beam, but such a complication is required in order to express the light scattered by the particle in terms of spherical coordinates (r; ; ') centered on the particle which is located at the origin of the Cartesian coordinates of the o¤-axis beam con…guration. We now consider the van de Hulst principle of localization in the o¤-axis coordinate system. Similarly as above, in the Rayleigh regime, only the n of f = 1 partial wave interacts with the particle, corresponding to a ray passing at a distance d of f = =4 from the z of f -axis which may be at a distance d on=of f which is large with respect to the z on -axis.

The fact that only n on = 1 and n of f = 1 are to be retained in the on-axis and o¤-axis con…gurations respectively does not depend on the respective values of BSCs associated with these particular partial waves, but results from the fact that only the Mie scattering coe¢ cient a 1 is to be retained in the calculation. But the fact that optical forces generated by rays passing at d on from z on in the on-axis con…guration or by rays passing at d on=of f from the same axis in the o¤-axis con…guration is linked to the speci…c values taken by the BSCs in these di¤erent con…gurations.

Furthermore, in both cases, a large enough particle would interact with all partial waves, whether it is in an on-axis or in an o¤-axis con…guration, and therefore the gradient force would be zero insofar as there is no longitudinal gradient of E:E associated in the whole beam, in agreement with the standard theory of optical forces expressed in terms of scattering and gradient forces. Conversely, a Rayleigh particle interacting only with the n = 1 partial waves does not interact with the whole beam but only with a part of the whole beam. The axicon terms revealed in the present paper then indeed show that the longitudinal gradient of ] E:E is not zero, from which we could expect a longitudinal gradient force which would be non zero as well. The fact that, nevertheless we have C g pr;z = 0, see Eq.( 42) in [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF], that is to say that the longitudinal gradient force is still zero may then be surprising. This is to be viewed as a coincidence or as a signi…cant result requiring an explanation still to be revealed. In particular, the evaluation of C g pr;z accounts for higher-order multipoles with n = 2 which does not occur in the gradient of ] E:E and which are not included in the use of the van de Hulst principle. Whatever the future of this story, our calculation exhibits a case when there is a zero longitudinal gradient force although the longitudinal gradients "felt" by the particle are not zero, but are indeed expressed by what we have called axicon terms.

6.3

Further works.

Further works are required to complement what is already known and understood. The next step should be a study of transverse optical forces similar to the previous work [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF] and to the present work, both of them devoted to longitudinal optical forces [START_REF] Gouesbet | Axicon optical forces and other kinds of transverse optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. Next, the on-axis con…guration of Bessel beams should be studied (i) to con…rm in a simpli…ed framework that axicon forces and terms indeed cancel in such a con…guration and (ii) to provide a comparison between the cases of dark and non dark beams which exhibit structural di¤erent features [START_REF] Gouesbet | Poynting theorem in terms of beam shape coe¢ cients and applications to axisymmetric, dark and non-dark, vortex and nonvortex beams[END_REF], [START_REF] Gouesbet | A darkness theorem for the beam shape coe¢cients and its relationship to higher-order non vortex Bessel beams[END_REF]. We may anticipate that speci…c departures from the standard theory could occur for dark beams, although the anticipation could be deceived. Also, the case of o¤-axis Gaussian beams which can be analytically handled using a localized beam model [START_REF] Gouesbet | Rigorous justi…cation of the localized approximation to the beam shape coe¢ cients in generalized Lorenz-Mie theory. II. O¤-axis beams[END_REF] should be investigated in order to check whether the newly axicon forces and terms are speci…c of axicon beams, or whether they occur in more general situations. Our current anticipation, which could again be deceived, is that indeed the newly studied terms occur in other situations than the one of axicon beams.

Conclusion.

Since the work of Ashkin, see [START_REF] Ashkin | Optical trapping and manipulations of neutral particles using lasers: A reprint volume with commentaries[END_REF] for a compilation, we are used to think of optical forces in terms of gradient and scattering forces. However, in a recent paper studying longitudinal optical forces in the case of o¤-axis Bessel beams in the Rayleigh regime, a new unexpected kind of optical forces, associated with the scattering force procedure, named axicon forces have been uncovered. In the present paper, we demonstrated the existence of axicon terms associated with gradient forces. These axicon terms are zero if the axicon angle is zero and/or if we are facing an on-axis con…guration. It is expected that such axicon forces and terms would occur whenever axicon angles are involved in the beam description. This would encompass Bessel beams whatever the polarization, all kinds of beams obtained from Bessel beam superpositions, e.g. frozen waves, Mathieu and Lommel beams [START_REF] Zamboni-Rached | Stationary optical wave …elds with arbitrary longitudinal shape by superposing equal frequency Bessel beams: Frozen waves[END_REF], [START_REF] Ambrosio | Analytical approach of ordinary frozen waves for optical trapping and micromanipulation[END_REF], [START_REF] Ambrosio | Optical forces experienced by arbitrary-sized spherical scatterers from superpositions of equal-frequency Bessel beams[END_REF], [START_REF] Ambrosio | Circularly symmetric frozen waves: Vector approach for light scattering calculations[END_REF], [START_REF] Ambrosio | Discrete vector frozen waves in generalized Lorenz-Mie theory: linear, azimuthal and radial polarization[END_REF], [START_REF] Ambrosio | Assessing the validity of the localized approximation for discrete superposition of Bessel beams[END_REF], [START_REF] Ambrosio | Zeroth-order continuous vector frozen waves for light scattering: exact multipole expansion in the generalized Lorenz-Mie theory[END_REF], [START_REF] Cha…q | Paraxial approximation of Mathieu beams through an apertured ABCD optical system[END_REF], [START_REF] Cha…q | Radiation pressure cross section exerted on homogeneous dielectric spherical particle by zeroth order Mathieu beams[END_REF][49], [START_REF] Cha…q | Scattering of Lommel beams by homogeneous spherical particle in generalized Lorenz-Mie theory[END_REF], or Laguerre-Gauss beams focused by a lens [START_REF] Van De Nes | On the conservation of fundamental optical quantities in non-paraxial imaging systems[END_REF], [START_REF] Van De Nes | Rigorous analysis of spheres in Gauss-Laguerre beams[END_REF], [START_REF] Török | The use of Gauss-Laguerre vector beams in STED microscopy[END_REF], [START_REF] Ambrosio | On localized approximations for Laguerre-Gauss beams focused by a lens[END_REF], [START_REF] Gouesbet | Finite series expressions to evaluate the beam shape coe¢ cients of a Laguerre-Gauss beam focused by a lens in an on-axis con…guration[END_REF].

Acknowledgments.

We acknowledge the support of "National Council for Scienti…c and Technological Development (CNPq) (426990/2018-8, 307898/2018-0)." We also heartly acknowledge useful reviewers who cleverly point out to a few relevant issues, and Pr. James Lock from Cleveland University, USA for his inquisitive mail exchanges.

We may then evaluate , and using Eqs.( 15)-( 18) of [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF]. In terms of A m n and B m n , we obtain:

Recalling that the A m n 's and B m n 's are real numbers, Eq.53 simpli…es to:

Next, we express the A m n 's and B m n 's using expressions ( 17)-( 18) of [START_REF] Gouesbet | Gradient, scattering and other kinds of longitudinal optical forces exerted by o¤-axis Bessel beams in the Rayleigh regime in the framework of generalized Lorenz-Mie theory[END_REF], in terms of generalized Legendre functions m n and m n , and, expressing these Legendre functions themselves, we obtain:

A 0 1 = J l+1 0 1 = J l+1 S (56)