Judy Najnudel

Rémy Müller

Thomas Hélie

David Roze

IDENTIFICATION OF NONLINEAR CIRCUITS AS PORT-HAMILTONIAN SYSTEMS

This paper addresses identification of nonlinear circuits for power-balanced virtual analog modeling and simulation. The proposed method combines a port-Hamiltonian system formulation with kernel-based methods to retrieve model laws from measurements. This combination allows for the estimated model to retain physical properties that are crucial for the accuracy of simulations, while representing a variety of nonlinear behaviors. As an illustration, the method is used to identify a nonlinear passive peaking EQ.

INTRODUCTION

Virtual analog modeling is an active research field, in particular within the audio community. Indeed, vintage analog audio effects are still appreciated among musicians, but original devices are now rare and delicate to maintain. A virtual replica then becomes a compelling alternative.

Historically, modeling methods have been roughly classified into two categories, white box and black box. White box modeling relies heavily on physics, and requires extensive knowledge on the circuit, from components datasheets to circuit schematics. Such modeling include Wave Digital Filters [START_REF] Fettweis | Wave digital filters: Theory and practice[END_REF][START_REF] Werner | Modeling circuits with arbitrary topologies and active linear multiports using wave digital filters[END_REF][START_REF] Bernardini | Wave digital modeling of nonlinear 3-terminal devices for virtual analog applications[END_REF] and State-Space representations [START_REF] Cohen | Real-time simulation of a guitar power amplifier[END_REF][START_REF] Macak | Simulation of Fender type guitar preamp using approximation and state space model[END_REF][START_REF] Holters | A generalized method for the derivation of non-linear state-space models from circuit DAFx[END_REF]. Black box modeling on the other hand is more concerned with reproducing a global input-output behavior, without necessarily capturing local phenomena taking place inside the system. It is usually more adaptable, but less physically interpretable. Neural Networks [START_REF] Damskägg | Real-time modeling of audio distortion circuits with deep learning[END_REF][START_REF] Martínez Ramírez | Deep learning for black-box modeling of audio effects[END_REF] are popular black box models. Volterra-based methods can be either white box oriented [START_REF] Hélie | Sound synthesis of a nonlinear string using Volterra series[END_REF][START_REF] Hélie | Volterra series and state transformation for real-time simulations of audio circuits including saturations: Application to the Moog ladder filter[END_REF], or black box oriented [START_REF] Boyd | Measuring Volterra kernels[END_REF].

A crucial counterpart of modeling is system identification, that is, retrieving model laws and parameters from measurements. Obviously, one designs an identification method with a specific model -white box or black box-in mind (see e.g. [START_REF] Orcioni | Identification of Volterra models of tube audio devices using multiplevariance method[END_REF][START_REF] Bouvier | Phasebased order separation for Volterra series identification[END_REF] for Volterra, or [START_REF] Sondhi | Lattice Wave Digital Filter based IIR system identification with reduced coefficients[END_REF] for WDF). In recent years nonetheless, hybrid modeling, or grey box, has gained considerable momentum for identification. Indeed, it often successfully combines desirable features from both paradigms. Such in-between modeling associates State-Space representations and polynomial models [START_REF] Paduart | Identification of nonlinear systems using polynomial nonlinear state space models[END_REF], State-Space representations and Neural Networks [START_REF] Parker | Modelling of nonlinear state-space systems using a deep neural network[END_REF], digital filters and Neural Networks [START_REF] Nercessian | Lightweight and interpretable neural modeling of an audio distortion effect using hyperconditioned differentiable biquads[END_REF] to list a few.

In this paper, we consider an identification method relying on the State-Space representation known as port-Hamiltonian systems (PHS) [START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF] on one hand, and kernel-based methods [START_REF] Park | Universal approximation using radial-basis-function networks[END_REF][START_REF] Vapnik | Support vector method for function approximation, regression estimation, and signal processing[END_REF] on the other hand. PHS are built as interconnected components with physical constitutive laws, so that the power balance of the system is structurally satisfied. This energy-based formulation can be associated with numerical methods that preserve the power balance and passivity in the discrete-time domain, for both linear and nonlinear systems. The PHS approach has proved relevant to simulate audio systems [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach[END_REF]. Although extensive work is concerned with linear PHS identification [START_REF] Medianu | Identification for port-controlled Hamiltonian systems[END_REF][START_REF] Benner | Identification of port-Hamiltonian systems from frequency response data[END_REF], nonlinear PHS identification is still little explored (see e.g. [START_REF] Cherifi | An overview on recent machine learning techniques for port-Hamiltonian systems[END_REF] for an overview). To address nonlinear behavior, we rely on reproducing kernels. Reproducing kernels have long proven to be a valuable and flexible tool for function approximation and solving of differential equations [START_REF] Schaback | Kernel techniques: From machine learning to meshless methods[END_REF]. As such, they already have succesfully been implemented for audio circuit modeling [START_REF] Daniel | Modeling nonlinear circuits with linearized dynamical models via kernel regression[END_REF]. Reproducing kernels are privileged here over deep learning approaches such as in [START_REF] Lutter | Deep lagrangian networks: Using physics as model prior for deep learning[END_REF], because the number of parameters to infer is much smaller, and because they are generally more interpretable. The specificity of our method resides in choosing a kernel and tailoring the regression so that key physical properties of the system, such as power balance and passivity, are retained.

This paper is organized as follows. In section 2, we give a brief overview of PHS. In section 3, we propose a parametrization of the PHS interconnection matrix. In section 4, we present an energy modeling based on reproducing kernels. Section 5 describes an optimization procedure to retrieve PHS parameters from measurements. In section 6, the complete method is tested on a virtual passive peaking EQ. Finally, some conclusions and perspectives are given in section 7.

PORT-HAMILTONIAN FORMULATION AND WORKING ASSUMPTIONS

Port-Hamiltonian formulation

The identification method described in this paper relies on port-Hamiltonian systems [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF], under a differential-algebraic formulation (PHS). This formulation allows the representation of a dynamical system as a network of:

1. storage components of state x and energy H (x), storing the power P stored = ∇H (x) ⊺ ẋ;

2. dissipative components of flow w and effort z (w), such as the dissipated power P diss = z (w) ⊺ w is non-negative for all w;

3. connection ports conveying the outgoing power Pext = u ⊺ y where u are inputs and y are outputs.

The flows f and efforts e of all components are coupled through a skew-symmetric interconnection matrix S = -S ⊺ (possibly de-

  ẋ w y   f = S   ∇H (x) z (w) u   e . (1)
Here, flows can either be currents (e.g. for capacitors) or voltages (e.g. for inductors), and vice versa for efforts. Such systems satisfy the power balance P stored + P diss + Pext = 0. Indeed, P stored + P diss + Pext = e ⊺ f = e ⊺ Se is zero since e ⊺ Se = (e ⊺ Se) ⊺ = -(e ⊺ Se) due to the skew-symmetry of S.

Under some additional assumptions (not detailed here), this formulation can be reduced to the following PH-ODE representation [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF]:

ẋ y = (J -R) ∇H (x) u , (2)
where J = -J ⊺ is skew-symmetric, and R = R ⊺ is positive semi-definite. Both matrices are possibly dependent on x. In particular, formulation (2) is possible when the dissipative law z is linear.

Throughout this paper, we adopt the passive sign convention for all components, including ports and external sources. This means that the current is defined positive when entering the component through the positive voltage terminal [START_REF] Bigelow | Power and energy in electric circuits[END_REF].

Example

As an illustration, consider the linear RLC circuit in series (Fig. 1a). The resistor obey Ohm's law, with flow iR (current through the resistor) and effort z(iR) = R iR = vR (voltage across the resistor). The inductor has state ΦL (magnetic flux linkage through the inductor), energy HL(ΦL) = Φ 2 L /2 L and effort H ′ L (ΦL) = ΦL/L = iL (current through the inductor). The capacitor has state qC (electric charge in the capacitor), energy HC(qC) = q 2 C /2 C and effort H ′ C (qC) = qC/C = vC (voltage across the capacitor). Kirchhoff's laws in receptor convention yield the PHS and PH-ODE representations in Fig. 1b and Fig. 1c, respectively.

Working assumptions

In the following, we limit ourselves to circuits verifying that:

1. The dissipation law z is linear.

2. The coefficients of J -R are constant.

3. Separability: the energy law is a separable function of the state, i.e. it takes the form

H (x) = Nx k=1 H k (x k ) , Nx = dim (x).

Smoothness:

H k is at least C 2 ∀k ∈ {1, ..., Nx}.

Convexity: the energy law is convex, i.e. H

(2)

k (x k ) ≥ 0 ∀k, x k , where H (2)
k denote the second derivative of H k (a formal definition can be found in Appendix A).

Assumptions 1 to 3 cover electronic circuits constituted of one-ports such as (possibly nonlinear, see [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF]) inductors, (possibly nonlinear, see [START_REF] Biolek | Real-world capacitor as a memcapacitive element[END_REF]) capacitors, linear conductors and linear resistors, which admit a PH-ODE formulation such as Eq. ( 2). This deliberately excludes nonlinear dissipative components such as diodes, transistors and vacuum tubes, which will be the object of future work. Assumptions 4 and 5 are stricter than necessary: actually, a sufficient condition to ensure passivity of storage components is for the energy to have an inferior bound [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF]. Nevertheless, enforcing convexity results in a desirable asymptotic behavior, and most energy laws are convex anyway (a notable exception concerns meta-stable ferromagnetic cores [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF]).

In the next sections, we assume that we can measure x, u and y, and look for an estimation of the reduced matrix J -R and an approximant of H, verifying both Eq. (2) and our working assumptions.

INTERCONNECTION MATRIX MODELING

Decomposition of J

Since J = -J ⊺ is skew-symmetric, it can be written as:

J = J (θ J ) := N J k=1 θ J k J k , θ J = θ J 1 , . . . , θ J N J ∈ R N J , (3)
where {J k } is the canonical base of skew-symmetric matrices, and

NJ = dim (J ) dim (J ) -1 / 2.
For the example of section 2.2, we have:

{J k } =      0 1 0 -1 0 0 0 0 0   ,   0 0 1 0 0 0 -1 0 0   ,   0 0 0 0 0 1 0 -1 0      , θ J = [1, 0, 1].

Decomposition of R

Since R is positive semi-definite, it admits a Cholesky decomposition [START_REF] Nicholas | Analysis of the Cholesky decomposition of a semi-definite matrix[END_REF][START_REF] Benoit | Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues. Application de la méthode à la résolution d'un système défini d'équations linéaires (procédé du Commandant Cholesky)[END_REF] and can be written as:

R = R(θ R ) = T (θ R )T (θ R ) ⊺ , T (θ R ) := N R k=1 θ R k T k , θ R = θ R 1 , . . . , θ R N R ∈ R N R , (4)
where {T k } is the canonical base of lower triangular matrices, NR = dim (R) dim (R) + 1 / 2, and diagonal coefficients are non-negative. For convenience, we choose the first dim (R) coefficients to be the diagonal coefficients.

For the example of section 2.2, we have:

{T k } =      1 0 0 0 0 0 0 0 0   ,   0 0 0 0 1 0 0 0 0   ,   0 0 0 0 0 0 0 0 1   ,   0 0 0 1 0 0 0 0 0   ,   0 0 0 0 0 0 1 0 0   ,   0 0 0 0 0 0 0 1 0      θ R = 0, √ R, 0, 0, 0, 0 .

ENERGY LAW MODELING

Optimal approximant

We consider the reproducing kernel theory to build an approximant of H. We refer to [START_REF] Schaback | Native Hilbert spaces for radial basis functions I[END_REF][START_REF] Schaback | Kernel techniques: From machine learning to meshless methods[END_REF] for a complete theoretical overview of Reproducing Kernel Hilbert Spaces (RKHS). In this section, we only recall practical results. Consider the classic scattered data interpolation problem, which is finding a function f verifying:

f (xi) = yi ∀i ∈ {1, ..., N }, (5)
for some given data {(x1, y1), ..., (xN , yN )}, [x1, ..., xN ] ∈ Ω N . Assuming that f belongs to a RKHS K, the optimal approximant of f in K is the function f defined as:

f : x → N j=1 λjK (x, xj) , (6)
with K : Ω × Ω → R the positive definite kernel inducing K. In Eq. ( 6), the vector λ is the solution of the linear system Kλ = y, where Ki,j = K (xi, xj). For instance, the kernel defined as:

K (xi, xj) = a π sinc a(xi -xj) ,
where 0 < a < +∞ and sinc(u) = sin(u)/u, is the reproducing kernel of the space of bandlimited continuous functions of bandwith in [-a, a].

In practice, the RKHS K of interest for a given problem is often unknown. Therefore what motivates the choice of kernel boils down to the properties one wishes to attach to the approximant, such as locality 1 , smoothness, interpolation behavior, sensitivity to noise, etc. However, there is an inherent tradeoff between these 1 that is, how much influence a data point has over its neighborhood. properties. Indeed, if locality is desirable (typically, because f contains high frequencies), choosing a local kernel is necessary. Nonetheless good reconstruction with a highly local kernel necessitates a high number of samples N . But a high number of samples can in turn lead to overfitting, especially with noisy data. Therefore, N is usually chosen low enough so that the approximant generalizes well and small changes in the given data do not impact the reconstruction too much. Concomitantly, the kernel support is adjusted through a scaling parameter ρ > 0, so that the approximant is reasonably smooth for the chosen N .

Once N is fixed, several strategies are available in order to find optimal parameters (the N samples and the scaling parameter ρ) for a given kernel, see e.g. [START_REF] Nelles | Nonlinear system identification[END_REF][START_REF] Gregory | On choosing "optimal" shape parameters for rbf approximation[END_REF][START_REF] Mongillo | Choosing basis functions and shape parameters for radial basis function methods[END_REF]. In this work, as the input space dimension is only 1 due to the separability of the energy law, we adopt a grid-based approach: the samples are distributed evenly along the measured data range, and the scaling parameter is adjusted so that some test function is approximated with an arbitrary degree of accuracy. The advantage of this approach resides in its simplicity. Should the energy law not be separable though, a more sophisticated strategy would be preferable. Figure 2 shows approximations from noisy data (SNR = 14 dB with a normal distribution) for different values of N , and K defined as:

K (xi, xj) =
(1 -r) 4 (1 + 4r) r = ∥xi -xj∥ < 1, 0 otherwise.

Choice of kernel

The only assumption we have on the energy H is the convexity assumption. Therefore, any kernel reproducing the convexity of H should be relevant. A simple way of enforcing the convexity of H is to construct a positive approximant of ∇ 2 H. An approximant of H can then be obtained by integration. To ensure the positivity of each H

k , we look for approximants of the form:

H (2) k (x k ) ≈ N H i=1 θ H k,i K (x k , x k,i ) , with θ H k,i ≥ 0 ∀k, i, (7)
where K is positive definite and continuous. Additionally, we choose K radial, of the form K (xi, xj) = ϕ ∥xi -xj∥ . Indeed, a radial kernel is local by construction. We also choose K compactly supported, so that interpolation matrices are sparse and computation is efficient. A possible choice of kernel verifying these properties is one of the Wendland functions [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF][START_REF] Schaback | Kernel techniques: From machine learning to meshless methods[END_REF] defined as:

ϕ(ρ; r) =      1 ρ 1 - r ρ 2 r = ∥xi -xj∥ < ρ, 0 otherwise, (8)
where the scaling parameter ρ is different for each component.

Energy law model

To simplify notations in the following, we denote

ϕ k,i : x → ϕ ρ k ;∥x -x k,i ∥ . (9)
Finally, the energy law approximant we look for has the form:

H θ H ; x := Nx k=1 N H i=1 θ H k,i ϕ (-2) k,i (x k ) , θ H k,i ≥ 0 ∀k, i, (10)
where ϕ

(-2)
k,i denotes the second antiderivative of ϕ k,i (formal definition in Appendix A). The {ϕ k,i } constitute a compactly supported radial basis of ∇ 2 H. The {x k,i } are called the centers of the radial basis [START_REF] Schaback | A practical guide to radial basis functions[END_REF].

PARAMETER ESTIMATION

Objective function

We denote f [j] the measured average flows at sample j:

f [j] = δx[j]fs y[j] , (11)
where δx

[j] = x[j + 1] -x[j]
and fs is the sample rate. Similarly we denote f (Θ)[j] the estimated flows at sample j:

f (Θ)[j] = J (θ J ) -R(θ R ) ∇H θ H , x[j], δx[j] u[j] , (12)
where θ H ⪰ 0, diag T (θ R ) ⪰ 0, and ∇H is the discrete gradient [START_REF] Mclachlan | Discrete gradient methods have an energy conservation law[END_REF] defined component-wise as:

∇H k θ H , x[j], δx[j] =                N H i=1 θ H k,i ϕ (-2) k,i x k [j] + δx k [j] -ϕ (-2) k,i x k [j] δx k [j] δx k [j] > ϵ, N H i=1 θ H k,i ϕ (-1) k,i x k [j] + δx k [j] 2 otherwise. ( 13
)
Denoting F (resp. F (Θ)) the dim(J ) × n matrix of measured flows (resp. estimated flows) at all n samples, we define the error ϵ(Θ) = F (Θ) -F , and the objective function

E(Θ) = F (Θ) -F 2
, where ∥∥ is the Frobenius norm. Finally we look for the optimal:

Θ ⋆ = arg min E(Θ) subject to θ H ⪰ 0, diag T (θ R ) ⪰ 0. ( 14
)
Here diag T (θ R ) denotes the diagonal elements of T (θ R ).

The sample rate fs is considered high enough so that the numerical error is negligible. Therefore we consider that we identify the continuous system, and the specific contribution of the discretization scheme to the global error [START_REF] Cieśliński | Discrete gradient algorithms of high order for one-dimensional systems[END_REF] is not addressed here.

Constrained optimization

To perform a constrained minimization of E, we rely on the Interior Point Method [START_REF] Boyd | Convex optimization[END_REF]. As this method is well documented, in this section we only provide a basic mathematical layout as a reminder.

We define the loss function L:

L : Θ → E(Θ) - 1 t   Nx k=1 N H i=1 ln θ H k,i + dim(J) k=1 ln θ R k   logarithmic barrier . ( 15
)
The logarithmic barrier [START_REF] Boyd | Convex optimization[END_REF] penalizes the minimization if all coefficients are not strictly positive. The parameter t is set by the user to enforce or relax the constraint 2 . A necessary condition to minimize L is finding a solution to:

F (Θ, µ) = ∂ϵ(Θ) ⊺ ϵ(Θ) -µ µ ⊙ Θ -1/t = 0, ( 16
)
where ∂ϵ is the Jacobian of ϵ, µ is the derivative of the logarithmic barrier w.r.t Θ, and ⊙ denotes the element-wise product. A solution is estimated with a damped Gauss-Newton iteration [START_REF] Deuflhard | Newton methods for nonlinear problems: Affine invariance and adaptive algorithms[END_REF]. Starting from a particular set of parameters Θ 0 meeting the constraints, the set is improved iteratively using:

Θ k+1 , µ k+1 ⊺ = Θ k , µ k ⊺ -α J -1 Θ k , µ k F Θ k , µ k (17)
where α ∈ [0, 1] is a damping coefficient computed with a line search [START_REF] Boyd | Convex optimization[END_REF] and J is the Jacobian of F defined as:

J (Θ, µ) = ∂ϵ(Θ) ⊺ ∂ϵ(Θ) -I diag (µ) diag(Θ) . (18)
Here diag(µ) (resp. diag(Θ)) denotes the square diagonal matrix with the elements of µ (resp. Θ) on its diagonal. The iteration is stopped when the error is sufficiently low, or, since the objective function is non convex, if the error starts increasing.

RESULTS FOR A VIRTUAL PASSIVE PEAKING EQ

Circuit parameters and data generation

We consider a passive peaking EQ [START_REF] Dennis A Bohn | Operator adjustable equalizers: An overview[END_REF] (Fig. 3a). The potentiometer wiper position is parametrized by γ ∈ [0, 1], where γ = 0 corresponds to the lowest position, and γ = 1 to the highest. This parameter determines the shape of the frequency response (Fig. 3b). The potentiometer, resistor and capacitor are all considered to be linear. The inductor is saturating with an effort law of the type [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF]:

iL = I0 ΦL Φsat -tanh ΦL η Φsat , (19)
where I0, Φsat, and η are model parameters (hysteresis is neglected here). Circuit parameters are set so that the center frequency is 50 Hz and the quality factor is 1. They are shown in Table 1. Synthetic measurement data are artificially generated for an input voltage of the form vin = U0 cos (2πf0t). The values of f0 and γ are chosen so that the circuit is maximally resonant. This way, nonlinearities of the inductor can be accurately captured for a plausible U0. The theoretical PH-ODE of the passive peaking EQ is found to be that in Appendix B. The generation is achieved discretizing Eq. ( 2), and performing a standard Newton-Raphson iteration at each sample (a detailed numerical scheme can be found in [START_REF] Najnudel | Simulation of an ondes Martenot circuit[END_REF]). Finally, some noise (SNR from 38 to 98 dB, with a normal distribution) is added to the data to test the robustness of the identification method.

Choice of kernel parameters and initialization

A minima, the approximant must reproduce a linear effort on the measured range of each state. We fix NH = 6. For this NH, we determine (for each component) the smallest ρ so that the relative error (g(x) -g(x))/g(x) × 100 on the test effort g(x) = x stays within some arbitrary bounds, chosen to be 10 %. Here, ρC = 4.10 -5 and ρL = 4.10 -3 yield satisfying results (Figs 4c-4d, for the inductor). Figure 4a (resp. 4b) shows the resulting basis for the approximation of ∇ 2 H (resp. ∇H).

Before performing the optimization procedure, an initial guess Θ 0 has to be estimated. To this end, the problem is linearized around the desired solution. Measurement data are generated with an input voltage small enough to observe a quasi-linear effort, so we have:

F ≈ (J -R) D M X, (20)
where D is a diagonal matrix with positive coefficients, and X is the dim(J )×n matrix of average states

x[i] = (x[i]+x[i+1])/2
and inputs at all n samples. Since n ≫ dim(J ), the matrix M is extracted using M = F X † , where † denotes the pseudo-inverse. Denoting Ms = -(M + M ⊺ )/2 the opposite of the symmetric part of M , the matrix R is initialized to a positive-definite matrix close to Ms. To this end, Ms is decomposed as:

Ms = U ΣV ⊺ . (21)
Then the matrix R defined as: is positive-definite, and R is initialized to:

R = U ΣU ⊺ (22)
R 0 = R/ R , (23)
to account for the (unknown) contribution of D. Finally, all θ J k and θ H k,i are initialized to 1.

Results

The optimization procedure returns a set of estimated parameters after less than 50 iterations. Here, constraint enforcement is privileged over convergence speed as there are no real-time requirements. Still, the estimation is faster compared to deep neural networks methods (mainly because there are far less parameters to estimate). New simulations are computed with the parameters estimated for each SNR. Figures 5a-5c show that the simulated states match very closely with the "measured" states (here with measured state SNR = 50 dB). Figures 5b-5d show the estimated effort laws. The estimated effort law for the capacitor is linear as expected. The saturating behavior of the inductor is accurately captured within the range of measured data. Figure 6 shows the simulated state Normalized Mean Square Error (NMSE = 20 log ∥x -x∥ /∥x∥ ) vs the measured state SNR. The NMSE stays low (around -60 dB) regardless of the SNR.

Finally, to evaluate the robustness of the estimated model, simulations and measurements are also run with different input amplitudes and frequencies than those used for the estimation. Figures 7a-7d show that the simulated states match closely with the measurements in that case also.

CONCLUSION

In this paper, we have presented an identification method to retrieve parameters of a circuit modeled as a port-Hamiltonian system, given measurements of state x, input u, and output y. This Table 1: Data generation parameters for the virtual passive peaking EQ. method allows the joint estimation of constitutive laws of storage components, as well as the interconnection matrix encoding the circuit topology. In turn, the estimated model may be used for passive guaranteed simulations.

Rp (Ω) R (Ω) γ C (F) I0 (A) Φsat (Wb) η U0 (V) f0 (
The method has been tested on a virtual peaking EQ, with accurate results. However, the method needs to be more thoroughly assessed against real measurements. In particular, we should control that the discretization scheme in the simulation error does not introduce too much numerical dispersion, which would alter the optimization process. In that case, the discretization error would have to be taken into account [START_REF] Matsuda | Estimation of ordinary differential equation models with discretization error quantification[END_REF], or the sample rate would have to be increased.

Another immediate perspective for this work is to extend the method in order to include nonlinear dissipative components. This would allow the inclusion of transistors and vacuum tubes, which are an important part of audio circuits.

It would also be interesting to adapt the method to a co-energy variables formulation. This way, measurements would only need to be voltages and currents, instead of charges and magnetic fluxes, which are much more difficult to obtain in practice. Identification from partial measurements (incomplete state, or input and output only) could be studied as well. This will be the object of future work.

f (k) : x → d dx f (k-1) (x). (24)
Similarly, we denote f (-k) a function defined by:

f (-k) : x → x 0 f (1-k) (u)du + C, ( 25
)
where C is a constant. In this work, C is taken so that f (-k) (0) = 0.

B. PH-ODE OF THE PASSIVE PEAKING EQ

With

αp = R Rp + R , Gp = 1 Rp + R , Rp ∥ R = RpR Rp + R , (26)
the PH-ODE of the passive peaking EQ on Fig. 3a is given by:

v C i L v in iout         i C . 1 . . v L -1 -γRp ∥ R Rp R (1 -γ) + 1 γRpGp γRp ∥ R i in . -γRpGp -Gp -αp vout . γRp ∥ R αp -Rp ∥ R . (27)
DAFx.8

Figure 1 :

 1 Figure 1: RLC circuit in series and its PHS representations.

Figure 2 :

 2 Figure 2: Two different approximants of the same function f (x) = sin(2πx) + sin(6πx) + sin(10πx) with noisy data.

2 4

 4 it can be increased dynamically during iteration. DAFx.Proceedings of the 23 rd International Conference on Digital Audio Effects (DAFx2020), Vienna, Austria, September 2020-21 Passive peaking EQ circuit. (b) Frequency response for several wiper positions.

Figure 3 :

 3 Figure 3: Passive peaking EQ and its frequency responses.

  (a) Chosen radial basis for ∇ 2 H. (b) Antiderivatives of the chosen radial basis for ∇ 2 H. (c) Resulting approximation of the effort g(x) = x.(d) Relative error on the effort.

Figure 4 :

 4 Figure 4: Chosen radial basis for the inductor and resulting approximation of the test effort g(x) = x.

  Capacitor simulated state. (b) Capacitor estimated effort. (c) Inductor simulated state. (d) Inductor estimated effort.

Figure 5 :

 5 Figure 5: Estimation results for the virtual passive peaking EQ, with measured state SNR = 50 dB.

Figure 6 :

 6 Figure 6: Simulated state Normalized Mean Square Error with respect to measured state SNR.


Judy Najnudel

Thomas Hélie

David Roze

Rémy Muller

Rémy Müller

From statistical physics to macroscopic port-Hamiltonian Systems: A roadmap

Keywords: physical modeling, statistical physics, thermodynamics, port-Hamiltonian systems

come

From statistical physics to macroscopic port-Hamiltonian Systems: A roadmap

INTRODUCTION

Physical modeling is concerned with the representation of a system governed by laws of physics. An accurate physical model provides insight into the system's behavior, beyond the conditions available to the experimenter.

Several physical modeling techniques have been developped over the years: see e.g. [START_REF] Ersal | A review of proper modeling techniques[END_REF] for an overview. In this paper, we consider in particular the statespace form known as port-Hamiltonian systems (PHS) introduced in [START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators[END_REF]; [START_REF] Van Der Schaft | Porthamiltonian systems theory: An introductory overview[END_REF]. PHS are multi-physical (that is, mechanical, electrical, thermal ... or a combination as well), and modular (a power-conserving or dissipative interconnection assembly of PHS is still a PHS). Most importantly, the PH formalism structurally fulfills the physical power balance. Moreover, numerical methods that preserve this guarantee in the discrete-time domain are available for simulations: [START_REF] Falaize | Passive guaranteed simulation of analog audio circuits: A port-hamiltonian approach[END_REF]. Nevertheless, simulating PHS with numerous degrees of freedom can prove an issue, as it involves very large matrices.

On the other hand, statistical physics is especially designed to describe complex systems with a reduced number of variables. From a collection of interacting elementary particles subject to constraints, statistical physics predicts macroscopic thermodynamic phenomena, such as entropy creation and phase transitions: [START_REF] Stowe | An introduction to thermodynamics and statistical mechanics[END_REF]; [START_REF] Landsberg | Thermodynamics and statistical mechanics[END_REF].

While thermodynamics has been broadly studied in the context of PHS modeling in e.g. [START_REF] Eberard | Port-Hamiltonian systems extended to irreversible systems: The example of the heat conduction[END_REF]; [START_REF] Eberard | An extension of Hamiltonian systems to the thermodynamic phase space: Towards a geometry of nonreversible processes[END_REF]; [START_REF] Ramirez | Irreversible port-Hamiltonian systems: A general formulation of irreversible processes with application to the CSTR[END_REF]; [START_REF] Delvenne | Finite-time thermodynamics of port-Hamiltonian systems[END_REF]; van der Schaft (2020), the proper derivation of thermodynamic variables through statistical physics is seldomly addressed. In this paper, we combine the two approaches and propose a series of steps to systematically construct a simplified yet physicallybased, structured macroscopic PHS, from a complex system described by statistical physics. This paper is structured as follows. Section 2 presents the key principles behind statistical physics, and details the derivation of a macroscopic PHS from a stochastic description. Section 3 applies the method to two illustrative examples: an ideal gas, and a ferromagnet. Finally, section 4 gives some conclusions and perspectives.

THEORETICAL FRAMEWORK

Micro-state of a system

Consider a system composed of particles of matter occupying a volume. Denote M the set of all possible configurations of all (or subsets of) countable particles for all volume values. A particular configuration m ∈ M is called a microstate of a system. For instance, each micro-state of a gas can be described in classical mechanics by a particular set of positions and momenta.

Each micro-state is mapped to a number of particles, a volume, and an energy, corresponding to the following functions:

N : M → N * m → N (m) micro-state number of particles, V : M → R + * m → V(m) micro-state volume, E : M → R m → E(m) micro-state energy,
(1) where the energy function E is assumed to have an inferior bound, i.e., E(m) ≥ c for some c ∈ R.

Accessible micro-states under constraints

According to experimental conditions, any quantity ϕ ∈ F = {N , V} can be fixed while E and the remaining unfixed quantities in F are left free to fluctuate 1 . In this case, the experimental constraints determine the set of accessible micro-states M a . For instance, a closed system (no exchange of particles with the environment) has a fixed N (m) = N 0 ∀m, while V and E can fluctuate. However these fluctuations are coupled through the quantity E(m) + P 0 V(m), where P 0 is the (constant) external pressure 2 . Similarly, an isochoric system has a fixed volume V(m) = V 0 , and the fluctuations of N and E are coupled through the quantity E(m) -µ 0 N (m), where µ 0 is the (constant) external chemical potential. Denote

Q = {E, E + P 0 V, E -µ 0 N , E + P 0 V -µ 0 N }. If the system is thermally insulated, the constraint on Ψ ∈ Q is Ψ (m) = Ψ 0 with a given Ψ 0 ∈ R, hence : M a = {m ∈ M | Ψ (m) = Ψ 0 }.
(2)

The case of a system in thermal contact with its environment is addressed in section 2.4.

Stochastic modeling and measure of uncertainty

The system may be in any accessible micro-state, and goes randomly from one accessible micro-state to another. As it is not possible to predict these fluctuations in a deterministic fashion, statistical physics adopts a stochastic description. This description assigns to each micro-state m a probability p(m) to be the actual state of the system.

Given a probability distribution p, a measure of the uncertainty on the fluctuating micro-state is the statistical entropy [START_REF] Gray | Entropy and information theory[END_REF], defined as (for a discrete distribution 3 ):

S : p -→ -k m∈Ma p(m) ln p(m), (3)
where k is a positive constant. The entropy is indeed a measure of lack of information. For instance, a probability distribution assigning 1 to some micro-state m 0 and 0 to all others would mean that the system is in the micro-state m 0 for certain. From (3), the system entropy would be zero for such a distribution. Conversely, an equiprobable distribution between all micro-states would maximize the uncertainty for an observer, as well as the entropy.

The entropy is a positive quantity, since 0 ≤ p(m) ≤ 1 ∀m. Moreover, the total entropy of two independent systems is the sum of their entropies, since p(m, m ) = p(m)p(m ). It follows that the entropy is an extensive quantity.

Principle of maximum entropy at thermodynamic equilibrium

Thermodynamic equilibrium is reached when the constraints on the system are met. At that point, the system stops evolving and provides a minimal amount of information. Therefore, at equilibrium, the micro-state probability distribution maximizes the statistical entropy, given the constraints on the system. Moreover, at thermodynamic equilibrium, statistical physics assumes the ergodic hypothesis (see e.g. [START_REF] Patrascioiu | The ergodic-hypothesis: a complicated problem in mathematics and physics[END_REF] for a discussion about its validity; here we will admit it as a working hypothesis).

2

The derivation of such coupling quantities is not detailed here. However, they are listed in Table 1, see also [START_REF] Graben | Unified treatment of adiabatic ensembles[END_REF].

3 It is similarly defined for a continuous distribution: the sum is replaced with an integral.

This hypothesis states that at equilibrium, the system visits all accessible micro-states, given a sufficiently long period of time. As a result, the temporal mean Ψ of a fluctuating quantity Ψ ∈ Q coincides with its expectation:

E[Ψ ] := m∈Ma p(m)Ψ (m) = Ψ . (4)
It follows that the micro-state probability distribution at equilibrium p is:

p = arg max p S(p) subject to      E[Ψ ] = Ψ m∈Ma p(m) = 1
(5)

In the following, the constant k in ( 3) is taken as the Boltzmann constant k b = 1.38 × 10 -23 J.K -1 . This ensures that the statistical entropy coincides with the thermodynamic entropy at equilibrium.

To solve (5), we introduce and optimize the Lagrangian:

L : (p, λ Ψ , λ 0 ) → -k b m∈Ma p(m) ln p(m) + λ Ψ E[Ψ ] -Ψ + λ 0   m∈Ma p(m) -1   , (6)
where λ Ψ and λ 0 are Lagrange multipliers. Case 1: system thermally insulated. From ( 2), E[Ψ ] = Ψ 0 . The Lagrangian simplifies, and solving ∂L = 0 yields the equiprobable distribution:

p m | Ψ 0 = 1 /Ω, Ω = card (M a ) .
(7) Case 2: system in thermal contact. Solving ∂L = 0 yields:

p m | Ψ = exp λ Ψ Ψ (m) /k b Z (λ Ψ ) , Z (λ Ψ ) = m∈Ma exp λ Ψ Ψ (m) /k b . (8)
A more detailed derivation can be found in [START_REF] Jaynes | On the rationale of maximum-entropy methods[END_REF].

For systems in thermal contact with their environment, the Lagrange multiplier λ Ψ has a direct physical interpretation. Indeed, from ( 4)-( 8), we deduce :

Ψ = ∂ ∂λ Ψ k b ln Z (λ Ψ ) . (9)
Moreover, reinjecting (8) in (3) yields the thermodynamic entropy S:

S = S(p ) = k b ln Z (λ Ψ ) -λ Ψ Ψ := S Ψ .
(10) It follows that S is a Legendre transform of k b ln Z and that:

-

λ Ψ = ∂ ∂ Ψ S Ψ = - 1 T , ( 11
)
where T is the temperature (both internal and external, at equilibrium).

By applying maximum entropy to each constraint, we systematically derive the micro-state probability and the entropy for the corresponding statistical ensemble (Table 1, see also [START_REF] Ray | Ensembles and computer simulation calculation of response functions[END_REF]).

Table 1. Statistical ensembles and associated constraints.

Ensemble Constraint p (m) Entropy Example Micro-canonical E(m) = E 0 1/Ω k b ln Ω Gas in an isolated tank Isoenthalpic-isobaric E(m) + P 0 V(m) = H 0 1/Ω k b ln Ω

Gas in a closed tank

No with a piston, thermally insulated

thermal contact Unnamed E(m) -µ 0 N (m) = L 0 1/Ω k b ln Ω Gas in a porous tank, thermally insulated Unnamed E(m) + P 0 V(m) 1/Ω k b ln Ω
Gas in a porous tank with a piston,

-µ 0 N (m) = R 0 thermally insulated Thermal contact Canonical E[E] = Ē exp(-E(m)/k b T ) Z(T ) k b ln Z(T ) + Ē/T Gas in a closed tank, in contact with a thermostat Isothermal-isobaric E[E + P 0 V] = H exp(-H(m)/k b T ) Z(T,P 0 ) k b ln Z(T, P 0 ) + H/T
Gas in a closed tank with a piston, in contact with a thermostat

Grand-canonical E[E -µ 0 N ] = L exp(-L(m)/k b T ) Z(T,µ 0 ) k b ln Z(T, µ 0 ) + L/T Gas in a porous tank, in contact with a thermostat Unnamed E[E + P 0 V -µ 0 N ] = R exp -R(m) /k b T R/T
Gas in a porous tank with a piston, in contact with a thermostat 2.5 Macro-state of a system and PHS formulation

Following the ergodic hypothesis, the macroscopic number of particles at equilibrium is

N = E[N ]. Similarly, the macroscopic volume at equilibrium is V = E[V].
A third macroscopic variable is the thermodynamic entropy S. If all external efforts constraining the system (that is, T, µ 0 , P 0 , or a combination of them depending on the experimental conditions) are kept constant in time, there is no dynamics since all macroscopic variables become constant in time. However if these quantities are allowed to vary (slowly, so that the ergodic hypothesis is still verified), one can study the system dynamics between successive equilibrium states.

In order to do that, we rely on port-Hamiltonian systems: [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]; van der Schaft et al. ( 2014), under a differential-algebraic formulation [START_REF] Beattie | Port-Hamiltonian descriptor systems[END_REF]. This formulation allows the representation of a dynamical system as a network of:

(1) storage components of state x and energy E (x), storing the power P stored = ∇E (x) ẋ;

(2) dissipative components of flow w and effort z (w), such as the dissipated power P diss = z (w) w is nonnegative for all w;

(3) connection ports conveying the outgoing power P ext = u y where u are inputs and y are outputs.

The flows f and efforts e of all components are coupled through a skew-symmetric interconnection matrix J = -J (possibly dependent on x):

  ẋ w y   f = J   ∇E (x) z (w) u   e . ( 12
)
For instance, an open system in a thermostat is governed by the conservative PHS in Fig. 1, where δ e S/dt denotes the outgoing entropy flow, int indexes internal flows and efforts, ext indexes external flows and efforts.

∇E(x) u T int µ int -P int δeS /dt Ṅext Vext               Ṡ . . . -

Method recap

To sum up, the derivation of a macroscopic PHS from a stochastic description is performed in 5 steps:

Step 1 Express micro-state m and laws N , V, E.

Step 2 Pick a statistical ensemble of interest, according to experimental hypotheses.

Step 3 Compute the corresponding micro-state probability distribution at equilibrium.

Step 4 Compute the internal energy Ē and the macroscopic state x = [S, N, V ] .

Step 5 Express the internal energy Ē as a function of the macroscopic state to obtain E(S, N, V ).

A PHS is finally derived by adding ports accounting for the environment, as shown for each example in the next section.

Remark on thermodynamic potentials

If only some external efforts constraining the system are kept constant, it is possible to work with an appropriate Legendre transform of the internal energy, in order to reduce the PHS formulation. For instance, the (conservative)

PHS of an open system in a constant thermostat, becomes that of Fig. 2, with F (T, N, V ) = E(S, N, V ) -T S.

Table 2 recaps usual Legendre transforms of the internal energy, also called thermodynamic potentials.

Table 2. Thermodynamic potentials.

Constant Potential

Internal energy E(S, N, V )

T
∇F (x) u µ int -P int Ṅext Vext       ẋ Ṅint . . -1 . Vint . . . -1 y µ ext 1 . . . P ext . 1 . . Figure 2.
Conservative PHS of an open system in a constant thermostat (dots represent zeros).

APPLICATIONS

Ideal gas in a thermostat

Consider N indiscernible, non-interacting atoms in a closed, rigid cube (constant volume V ) in contact with a thermostat at temperature T .

Step 1. A micro-state of an atom is a vector m = [n x , n y , n z ] ∈ N * 3 , where n i is the energy quantization along axis i. Denoting by m the atom mass, a the cube edge length, and h the Planck constant, the Hamiltonian of the atom is given by [START_REF] Davies | The physics of low-dimensional semiconductors: An introduction[END_REF]:

E(m) = n 2 x + n 2 y + n 2 z h 2 /8ma 2 . ( 13
)
Step 2. From the hypotheses, the statistical ensemble of interest is the canonical ensemble (see Table 1).

Step 3. The micro-state probability distribution at equilibrium is p (m) = exp -E(m) /k b T /Z(T ).

Step 4. From (8), the partition function for one atom is:

Z 0 (T ) = 2π T /T0 3/2 , with T 0 = h 2 /ma 2 k b (14) (proof in Appendix A).
For N atoms, the partition function becomes:

Z(T ) = Z0(T ) N /N! = 2π T /T0 3N/2 /N!.
(15) For the canonical ensemble, (9) can be rewritten as:

Ē = k b T 2 ∂ ln Z ∂T (T ), (16)
so that:

Ē = 3 /2 N k b T.
(17) From (10), we deduce:

S = k b ln Z(T ) + Ē/T = S 0 + 3 /2 N k b ln T /T0 , S 0 = 3 /2 N k b ln 2πe / N ! 2 /3N . ( 18
)
Step 5. From (18), we deduce:

T = T 0 exp 2 (S-S0) /3Nk b . (19)
Reinjecting in ( 17), we obtain:

E(S) = 3 /2 N k b T 0 exp 2 (S-S0) /3Nk b .
(20) It is immediately verified that, as expected, ∂E/∂S = T . As N and V are constant, they are removed from the state of the PHS and the only state variable is the entropy S. The constraint due to the thermostat is expressed at the ports of the PHS. Finally, we obtain the following PHS:

∇E(x) u T int δeS /dt ẋ Ṡ . -1 y T ext 1 .
.

Ferromagnetic core in a thermostat

Here, we sum up results of [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF]. The goal is to illustrate the modeling of a more complex system exhibiting phase transitions and dissipation.

Consider a ferromagnetic core with constant number of atoms N , constant volume V , in contact with a thermostat at temperature T .

Step 1. Following the model of [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF] (see also [START_REF] Newell | On the theory of the Ising model of ferromagnetism[END_REF]; [START_REF] Strecka | A brief account of the Ising and Ising-like models: Mean-field, effective-field and exact results[END_REF]), the core is represented as a set of N adimensional magnetic moments, interacting with one another. A micro-state of the core is a particular configuration m ∈ M = {-1, 1} N . The corresponding energy is the Heisenberg Hamiltonian:

E(m) = -1 /2 m J ex m, (21)
where each coefficient J exi,j is the exchange energy between moment i and moment j [START_REF] Liechtenstein | Exchange interactions and spin-wave stiffness in ferromagnetic metals[END_REF]. Assuming isotropic interactions affecting nearest neighbours only, this exchange energy simplifies to:

J exi,j = J i, j nearest neighbours, i = j, 0 otherwise, ( 22
)
where J is a constant energy characterizing the material.

Step 2. From the hypotheses, the statistical ensemble of interest is the canonical ensemble (see Table 1).

Step 3. The micro-state probability distribution at equilibrium is p (m) = exp -E(m) /k b T /Z(T ).

Step 4. Above dimension 2 in space, there is no analytic expression for Z(T ). Assuming small micro-state fluctuations, we rely on the mean-field approximation: Utermohlen (2018), which yields In ( 23), the auxiliary variable m ∈ [-1, 1] can be interpreted as a mean magnetic moment, while q is the (constant) number of nearest neighbours of each moment. Reinjecting ( 23) in ( 16), the mean-field internal energy is found to be:

Z M F (m, T ) = exp -N Jqm 2 /2 k b T 2 cosh Jqm /k b T N . (23)
Ē ≈ E M F (m, T ) = E 0 m 2 /2 -m tanh mTc /T , ( 24
)
where E 0 = N Jq, and T c = Jq/k b is the critical temperature, above which the core becomes paramagnetic. Reinjecting ( 23) in ( 10), we obtain the thermodynamic entropy:

S ≈ S M F (m, T ) = S 0 f mTc /T , (25)
where S 0 = N k b and f (χ) = ln 2 cosh (χ) -χ tanh (χ).

Step 5. Finally, introducing the total magnetic flux B Vcore = m µ 0 M s V , where µ 0 = 4π × 10 -7 H.m -1 is the vacuum magnetic permeability and M s is the core saturation magnetization, we obtain:

E core (S, B Vcore ) = E 0 1 /2 B Vcore /B Vs 2 -B Vcore /B Vs g(S) , (26)
where B Vs = µ 0 M s V and g(S) = tanh f -1 S /S0 . The internal effort is ∇E core (S, B Vcore ) = [T core , H core ] , where T core and H core denote the core temperature and the internal magnetic field, respectively. Figure 3 shows the core going from two meta-stable equilibrium states to one stable equilibrium state, as the entropy increases. This corresponds to a phase transition from ferromagnetic to paramagnetic.

Dissipation. When the core is constrained by an exterior magnetic field H in (created by a coil for instance), jumps between meta-stable states, called Barkhausen jumps, occur: [START_REF] Bertotti | Hysteresis in magnetism: For physicists, materials scientists, and engineers[END_REF]. These jumps are damped (Fig. 4) due to domain structure and non-homogeneities (see [START_REF] Kittel | Physical theory of ferromagnetic domains[END_REF]). This damping can be modeled with a linear magnetic resistor r core connected in series with the core (Fig. 5a). As the difference of energy before and after a jump is entirely dissipated as heat [START_REF] Bertotti | Hysteresis in magnetism: For physicists, materials scientists, and engineers[END_REF], Barkhausen jumps are also responsible for the variation of entropy creation δ i S in the core. The associated thermal power is equal to the magnetic power dissipated through r core , so that:

r core H 2 rcore -T core δiS /dt = 0, (27) where dt stands for an infinitesimal increment of time. We deduce the dissipative flow and effort: w = [-T core , H rcore ] z(w) = rcoreH 2 rcore/T core , r core H rcore = δiS /dt, ḂVr core . (28)

Figure 4. Damped Barkhausen jumps and resulting hysteresis during a cycle. The red curve is the theoretical core internal magnetic field, the blue curve is the real trajectory followed by the total magnetic flux of the core constrained by an external magnetic field.

Finally, the second law of thermodynamics Landsberg (2014) states: Ṡ = δiS /dt -δeS /dt.

(29) Kirchoff's laws in receptor convention, as well as ( 26)-( 28)-( 29), yield the PHS formulation in Fig. 5b. Note that in addition to the power balance, the entropy balance is made explicit in this formulation.

This model has been successfully implemented in several applications simulating audio circuits (see [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF]).

CONCLUSION

In this paper, we presented a method for deriving an ad hoc macroscopic port-Hamiltonian formulation for a constrained system described by statistical physics. This method is based on two key principles of statistical physics (that is, maximum entropy and the ergodic hypothesis), and expresses the internal energy of the system as a function of extensive thermodynamic variables. The exchanges of energy and entropy are made explicit through the PH formulation. Constraints are modeled in ports as inputs, and entropy creation is modeled as dissipation. Alternative formulations based on thermodynamic potentials are also possible, under certain conditions.

The method is applied to model two systems: an ideal gas in a thermostat, and a ferromagnet in a thermostat. In the case of the ferromagnet, the method successfully captures complex phenomena such as phase transitions, meta-stability, dissipation and entropy creation, with only two state variables and one dissipative component.

A perspective of this work is to extend the method to other microscopic quantities and constraints, not restricted to the system number of atoms and volume. Another perspective is to generalize the mean-field approximation (for systems composed of particles interacting with one another), and introduce appropriate auxiliary variables in order to obtain approximations of the partition function with an arbitrary degree of accuracy. This shall be the object of future work.

  Helmholtz free energy F (T, N, V ) = E(S, N, V ) -T S P Enthalpy H(S, N, P ) = E(S, N, V ) + P V T, P Gibbs free energy G(T, N, P ) = H(S, N, P ) -T S µ Hill energy L(S, µ, V ) = E(S, N, V ) -µN T, µ Grand potential J(T, µ, V ) = L(S, µ, V ) -T S P, µ Ray energy R(S, µ, P ) = E(S, N, V ) + P V -µN T, P, µ Guggenheim energy R(S, µ, P ) -T S = 0

Figure 3 .

 3 Figure 3. Energy function of the ferromagnetic core.

  Core with dissipation constrained by an external magnetic field.

Figure 5 .

 5 Figure5. Equivalent circuit and PHS for the core constrained by an external magnetic field with dissipation. arXiv:1705.09081.[START_REF] Bertotti | Hysteresis in magnetism: For physicists, materials scientists, and engineers[END_REF]. Hysteresis in magnetism: For physicists, materials scientists, and engineers. Academic Press.[START_REF] Davies | The physics of low-dimensional semiconductors: An introduction[END_REF]. The physics of low-dimensional semiconductors: An introduction. Cambridge University Press. Delvenne, J.C. and Sandberg, H. (2014). Finite-time thermodynamics of port-Hamiltonian systems. Physica D, 267, 123-132.[START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]. Modeling and control of complex physical systems: the port-Hamiltonian approach. Springer.[START_REF] Eberard | Port-Hamiltonian systems extended to irreversible systems: The example of the heat conduction[END_REF]. Port-Hamiltonian systems extended to irreversible systems: The example of the heat conduction. IFAC Proceedings Volumes, 37(13), 243-248. Eberard, D., Maschke, B., and van der Schaft, A. (2007). An extension of Hamiltonian systems to the thermodynamic phase space: Towards a geometry of nonreversible processes. Rep. Math. Phys., 60(2), 175-198. Ersal, T., Fathy, H.K., Rideout, D.G., Louca, L.S., and Stein, J.L. (2008). A review of proper modeling techniques. J. Dyn. Syst. Meas. Contr., 130(6). Falaize, A. and Hélie, T. (2016). Passive guaranteed simulation of analog audio circuits: A port-hamiltonian approach. Appl. Sci., 6(10), 273. Graben, H. and Ray, J.R. (1991). Unified treatment of adiabatic ensembles. Phys. Rev. A, 43(8), 4100. Gray, R.M. (2011). Entropy and information theory. Springer. Ising, E. (1925). Beitrag zur Theorie des Ferromagnetismus. Zeitschrift für Physik, 31(1), 253-258. Jaynes, E.T. (1982). On the rationale of maximum-entropy methods. Proc. IEEE, 70(9), 939-952. Kittel, C. (1949). Physical theory of ferromagnetic domains. Rev. Mod. Phys., 21(4), 541. Landsberg, P.T. (2014). Thermodynamics and statistical mechanics. Courier Corporation.[START_REF] Liechtenstein | Exchange interactions and spin-wave stiffness in ferromagnetic metals[END_REF]. Exchange interactions and spin-wave stiffness in ferromagnetic metals. J. Phys. F: Met. Phys, 14(7), L125.[START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators[END_REF].An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators. J. Franklin Inst., 923-966. Najnudel, J., Müller, R.,[START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF].A power-balanced dynamic model of ferromagnetic coils. In Proc. of the 23rd Int. Conf. on Digital Audio Effects

Note that F is susceptible to contain other functions, depending on the physics (electric, magnetic, etc) of the considered system.

Appendix A. PARTITION FUNCTION OF AN ATOM IN A BOX

Denoting 0 = h 2 /(8ma 2 ), we have

(A.1)


Laura Thiebaud
email: laura.thiebaud@univ-lorraine.fr

Sophie Legeai
email: sophie.legeai@univ-lorraine.fr

Nicolas Stein
email: nicolas.stein@univ-lorraine.fr

-free electrochemical deposition in an ionic liquid

Keywords: Electrodeposition, Nanowires, Tellurium, Ionic liquid, Hollow nanostructures 1

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

Introduction

Tellurium (Te) and its binary or ternary compounds are used in many domains such as piezoelectric devices [START_REF] Lee | High-power density piezoelectric energy harvesting using radially strained ultrathin trigonal tellurium nanowires assembly[END_REF], gas sensors [START_REF] Tsiulyanu | Concentration induced damping of gas sensitivity in ultrathin tellurium films[END_REF], photoconductors [START_REF] Li | Facile electrochemical synthesis of tellurium nanorods and their photoconductive properties[END_REF], solar cells [START_REF] Ito | Electrochemical deposition of Te and Se on flat TiO2 for solar cell application[END_REF] and thermoelectric generators and coolers [START_REF] Abad | Thermoelectric properties of electrodeposited tellurium films and the sodium lignosulfonate effect[END_REF]. Thermoelectricity, which consists in a reversible conversion of heat into electricity and conversely, allows the harvesting of wasted heat and therefore is particularly interesting for the improvement of energy systems efficiency [START_REF] Martín-González | Nanoengineering thermoelectric for 21 st century: Energy harvesting and other trends in the field[END_REF]. The thermoelectric performances are evaluated by the figure of merit ZT given by the equation 1:

Equation1 ZT= S 2 σ Кe+ Кg T
where S is the Seebeck coefficient, σ is the electrical conductivity, Кe is the electronic thermal conductivity, Кg is the phonon thermal conductivity and T is the temperature.

The performances of thermoelectric compounds can be improved by nanostructuring [START_REF] Martín-González | Nanoengineering thermoelectric for 21 st century: Energy harvesting and other trends in the field[END_REF][START_REF] Borca-Tasciuc | Thermal properties of electrodeposited bismuth telluride nanowires embedded in amorphous alumina[END_REF][START_REF] Dresselhaus | New directions for low-dimensional thermoelectric materials[END_REF][START_REF] Kanatzidis | Nanostructured thermoelectrics: The new paradigm ?[END_REF].

Among nanostructures, one-dimensional nanostructures like nanowires [START_REF] Caballero-Calero | Thermoelectric nanowires: A brief prospective[END_REF][START_REF] Kim | Thermoelectricity in semiconductor nanowires[END_REF][START_REF] Yang | Thermal conductivity of simple and tubular nanowires composites in the longitudinal direction[END_REF][START_REF] Zhang | Nanostructures for thermoelectric applications: synthesis, growth mechanism, and properties studies[END_REF] and core-shell nanowires [START_REF] Zhang | Nanostructures for thermoelectric applications: synthesis, growth mechanism, and properties studies[END_REF][START_REF] Yang | Thermal conductivity modeling of core-shell and tubular nanowires[END_REF][START_REF] Poon | Thermal conductivity of core-shell nanocomposites for enhancing thermoelectric performance[END_REF][START_REF] Yang | The effective thermoelectric properties of core-shell composites[END_REF][START_REF] Kang | Reduction of lattice thermal conductivity in single Bi-Te core/shell nanowires with rough interface[END_REF] are particularly promising due to the high phonon interface scattering and the blocking of the phonon conduction along the wire axis which lead to an increase of the thermoelectric conversion efficiency. Te nanowires can represent a base to create core-shell structures, being the core that could be further capped with a shell to get arrays of Te-based thermoelectric nanostructures like BiTe [START_REF] Kang | Reduction of lattice thermal conductivity in single Bi-Te core/shell nanowires with rough interface[END_REF][START_REF] Zhang | Enhanced thermoelectric properties of core/shell heterostructure nanowires composites[END_REF][START_REF] Zhou | Thermoelectric properties of individual electrodeposited bismuth telluride nanowires[END_REF][START_REF] Fang | Thermoelectric properties of silver telluride-bismuth telluride nanowire heterostructure synthesized by site-selective conversion[END_REF], AgTe [START_REF] Fang | Thermoelectric properties of silver telluride-bismuth telluride nanowire heterostructure synthesized by site-selective conversion[END_REF][START_REF] Gnanadurai | Investigation on the influence of thickness and temperature on the Seebeck coefficient of silver telluride thin films[END_REF][START_REF] Yang | Enhanced thermoelectric properties in bulk nanowires heterostructure-based nanocomposites through minority carrier blocking[END_REF] or LaTe [START_REF] May | Thermoelectric performance of lanthanum telluride produced via mechanical alloying[END_REF].

Various approaches have been developed to synthesize Te nanowires, such as physical and chemical vapor deposition [START_REF] Hawley | Shape-controlled vapor-transport growth of tellurium nanowires[END_REF][START_REF] Mohanty | Large scale synthesis of highly pure single crystalline tellurium nanowires by thermal evaporation method[END_REF], hydrothermal process [START_REF] Park | Aqueous chemical synthesis of tellurium nanowires using a polymeric template for thermoelectric materials[END_REF][START_REF] Liu | Shape-controlled synthesis and growth mechanism of one-dimensional nanostructures of trigonal tellurium[END_REF][START_REF] Liang | Synthesis of tellurium nanowires and their transport property[END_REF], galvanic displacement reaction [START_REF] Wu | Size controlled synthesis of Tellurium nanorices by galvanic displacement reaction of aluminium[END_REF][START_REF] Elazem | Morphology change of galvanically displaced one-dimensional tellurium nanostructures via controlling the microstructure of sacrificial Ni thin films[END_REF][START_REF] Jeong | Template-free synthesis of vertically oriented tellurium nanowires via a galvanic displacement reaction[END_REF] or electrochemical deposition (ECD) [START_REF] She | Template-free electrodeposition of one-dimensional nanostructures of tellurium[END_REF][START_REF] Zhao | Ordered tellurium nanowire arrays and their optical properties[END_REF][START_REF] Kumar | Optical and electrical studies of vertically oriented tellurium nanowires arrays produced by template electrodeposition[END_REF][START_REF] Al-Salman | Template-free electrochemical synthesis of high aspect ratio Sn nanowires in ionic liquids: A general route to large-area metal and semimetal nanowires arrays?[END_REF][START_REF] Szymczak | Template-free electrodeposition of tellurium nanostructures in a room-temperature ionic liquid[END_REF]. Arrays of Te 1D nanostructures have already been obtained by ECD, mainly by template-assisted electrodeposition, using anodic aluminium oxide (AAO) [START_REF] Zhao | Ordered tellurium nanowire arrays and their optical properties[END_REF] or polycarbonate etched membranes templates [START_REF] Kumar | Optical and electrical studies of vertically oriented tellurium nanowires arrays produced by template electrodeposition[END_REF]. However, template-assisted methods involve multiple fabrication steps, like the preparation of the template and its removing which is necessary for the further deposition of a thin upper layer leading to core-shell structures. A single step template-free synthesis route would be more convenient and easier to perform. Room-Temperature Ionic Liquids (RTILs) are promising solvents for template-free synthesis as they can act as capping agents leading to the synthesis of nanostructures. As example, Ma et al. synthesized Te nanowires via a chemical route using triethanolamine tetrafluoroborate ionic liquid as an additive in an aqueous solution [START_REF] Ma | Growth of tellurium nanowire bundles from an ionic liquid precursor[END_REF]. Additionally RTILs present specific advantages for ECD purposes such as a wide electrochemical window, a good electrical conductivity and a high thermal stability which gives access to highly crystalline materials. Template-free electrodeposition of 1D nanostructures in RTILs is little present in the literature. Steichen and Dale showed that selenium nanorods can be grown in an imidazolium-based ionic liquid mixture: 1-ethyl-3-methylimidazolium tetrafluoroborate/chloride [C2mIm][BF4]/[C2mIm]Cl [START_REF] Steichen | Synthesis of trigonal selenium nanorods by electrodeposition from an ionic liquid at high temperature[END_REF]. Al nanowires were obtained by Su et al. in a Lewis acidicionic liquid AlCl3/trimethylamine hydrochloride [START_REF] Su | Electrodeposition of aluminium wires from the Lewis acidic AlCl3/trimethylamine hydrochloride ionic liquid without using a template[END_REF]. With the use of imidazolium based ionic liquids, several works highlight the synthesis of metals (Zn [START_REF] Al-Salman | Template-free electrochemical synthesis of high aspect ratio Sn nanowires in ionic liquids: A general route to large-area metal and semimetal nanowires arrays?[END_REF], Sn [START_REF] Al-Salman | Template-free electrochemical synthesis of high aspect ratio Sn nanowires in ionic liquids: A general route to large-area metal and semimetal nanowires arrays?[END_REF][START_REF] Leong | Electrochemistry on tin in the 1-ethyl-3methylimidazolium dicyanamide room temperature ionic liquid[END_REF], Co [START_REF] Hsieh | Speciation of cobalt-chloride-based ionic liquids and electrodeposition of Co wires[END_REF]) and alloys (NiZn [START_REF] Yang | Single-step large-scale and template-free electrochemical growth of Ni-Zn alloy filament arrays from a zinc chloride based ionic liquid[END_REF], CuZn [START_REF] Hsieh | Electrodeposition of CuZn from chlorozincate ionic liquid : from hollow tubes to segmented nanowires[END_REF], AuZn [START_REF] Hsieh | One-step electrochemical fabrication of nanoporous gold wire arrays from ionic liquid[END_REF] or FeCoZn [START_REF] Yang | Direct electrodeposition of FeCoZn wire arrays from a zinc chloride-based ionic liquid[END_REF]) nanowires. Very few works were published concerning Te nanostructures synthesis in RTILs. Al-Salman and co-workers recently [START_REF] Al-Salman | Template-free electrochemical synthesis of high aspect ratio Sn nanowires in ionic liquids: A general route to large-area metal and semimetal nanowires arrays?[END_REF] synthesized long Te polycrystalline nanowires of several micrometers with diameters between 25 and 35 nm in the 1-butyl-1-methylpyrrolidinium bis(trifluoromethylsulfonyl)imide ionic liquid containing SiCl4. The authors show in this article that the capping agent is not here the ionic liquid but rather SiCl4 that promotes the growth of Te nanowires. However, the use of SiCl4 leads to the presence of an amorphous and passive layer on the surface of the nanowires, composed of Si, C and O, which will alter the electrical properties of these materials.

In a previous work [START_REF] Szymczak | Template-free electrodeposition of tellurium nanostructures in a room-temperature ionic liquid[END_REF], we showed that it is possible to obtain free standing single crystalline Te nanowires at moderate temperature (60-100°C) in a piperidinium-based ionic liquid mixture: EOPipTFSI:EOPipBr (1-ethyl-1-octylpiperidinium cation, bis(trifluoromethylsulfonyl)imide and bromide anions) without the need of any template.

However, the obtained length of the nanowires was limited to 1 µm at best, which is not enough for future applications in thermoelectricity. The morphology, size and organization of the structures are very important characteristics that impact physical and chemical properties of nanostructured materials. The modulation of these parameters can be achieved by varying electrodeposition conditions in order to obtain one-dimensional nanowires with a high length to diameter aspect ratio for optimal thermoelectric properties. In our previous work, only the influence of electrolyte temperature was investigated. It appeared that the best crystallinity of the deposits was obtained at 100°C. Here, we present a systematic study of other growth parameters like the applied fixed potential and the composition of the electrolyte on the deposits morphology. The morphology of Te electrodeposits was systematically investigated by scanning electron microscopy (SEM) and transmission electron microscopy (TEM). This work shows that 1D nanostructures with different morphologies can be obtained resulting from specific mass transport phenomena.

Experimental

TeCl4 (99%, Alfa Aesar) was used as Te precursor. EOPipTFSI and EOPipBr were synthesized according to [START_REF] Triolo | Nanoscale organization in piperidinium-based room temperature ionic liquids[END_REF][START_REF] Traore | New insight into indium electrochemistry in a Tf2N-based room-temperature ionic liquid[END_REF] and were characterized by nuclear magnetic resonance analysis (Bruker Spectrospin 300 MHz) and stored in an argon filled glove box (MBraunLabstar; O2/H2O levels < 1 ppm). A mixture of EOPipTFSI and EOPipBr was used as electrolyte. EOPipBr improves the solubility of Te IV salts in TFSI-based ionic liquids probably by the formation of complex halide species [START_REF] Jeng | Electrochemistry of tellurium(IV) in the basic aluminium chloride-1-methyl-3-ethylimidazolium chloride room temperature molten salts[END_REF]. As the viscosity strongly increases with the EOPipBr content [START_REF] Szymczak | Electrodeposition of stoichiometric bismuth telluride Bi2Te3 using a piperidinium ionic liquid binary mixture[END_REF], the amount of EOPipBr was limited to 5% (mol/mol) that represents the best compromise between a good solubility of the Te IV salt and a low electrolyte viscosity.

Electrochemical experiments were performed in a glove box using a Biologic potentiostat (VSP300). Platinum-coated glass slides (25×10×1 mm, Pt coating thickness 300 nm with a thin Cr sublayer that ensures a good adherence of the Pt layer), supplied by Applications Couches Minces (ACM), were used both as working and counter electrodes and a Pt wire was used as quasi-reference electrode. The volume of the bath was equal to 5 mL leading to an immerged slide surface of 0.5 cm 2 . Some specific electroplating synthesis were performed by using a bulk Te anode, which was beforehand prepared by melting Te shots in a square tin at 300°C. The resulting metallic block was then polished with successively finer grades of SiC Emery paper up to 2400 before use.

Results and Discussion

A typical cathodic linear voltammogram of Te IV in the piperidinium-based ionic liquid mixture is presented in Figure 1.

We can observe two reduction peaks at -0.23V (C1) and -1.2V (C2) vs Pt reference electrode.

The electrochemical behavior is the same as observed by Jeng et al. [START_REF] Jeng | Electrochemistry of tellurium(IV) in the basic aluminium chloride-1-methyl-3-ethylimidazolium chloride room temperature molten salts[END_REF] and Tsai et al. [START_REF] Tsai | Voltammetric study and electrodeposition of tellurium, lead, and lead telluride in room-temperature ionic liquid 1ethyl-3-methylimidazolium tetrafluoroborate[END_REF] in imidazolium based ionic liquids. Consequently, the cathodic peak C1 can be attributed to the reduction of Te IV into Te 0 , this latter being further reduced in Te -II at peak C2, according to [START_REF] Jeng | Electrochemistry of tellurium(IV) in the basic aluminium chloride-1-methyl-3-ethylimidazolium chloride room temperature molten salts[END_REF][START_REF] Tsai | Voltammetric study and electrodeposition of tellurium, lead, and lead telluride in room-temperature ionic liquid 1ethyl-3-methylimidazolium tetrafluoroborate[END_REF]. This latter signal represents the cathodic limit for the deposition of Te 0 films [START_REF] Szymczak | Template-free electrodeposition of tellurium nanostructures in a room-temperature ionic liquid[END_REF]. It should be noticed that the shape of the cathodic linear voltammograms is the same whatever the tellurium salt concentration, exhibiting the C1 and C2 peaks. However the potential values of the peaks can vary in a 50 mV potential range.

Influence of Te IV concentration

The influence of Te IV concentration on Te electrodeposition was first studied at C1 peak potential value (Ep). The results show unambiguously that the increase of the Te precursor concentration has a negative impact on the synthesis of low dimensional structures. For [Te IV ] = 25 mM (Figure 2a), corresponding to the maximal solubility of TeCl4 in EOPipTFSI:EOPipBr 95:5 (mol%), the obtained deposit is dense, composed of nanostructures homogeneous in size of about 150 ± 20 nm in diameter but only 660 ± 50 nm in length . For [Te IV ] = 20 mM (Figure 2b), a compact deposit is obtained, composed of 1.0 ± 0.1 µm long self-standing nanowires with a mean diameter of 140 ± 20 nm, homogeneous in size. By decreasing the concentration at 5 mM (Figure 2c), longer nanowires are obtained, having a mean diameter of 130 ± 20 nm and a mean length of 1.2 ± 0.2 µm corresponding to an aspect ratio of 9.2 ± 3.0, the longest nanowires reaching 1.8 ± 0.1 µm. In these conditions, the nanowires are less homogeneous in size. One can also observe on Figure 2 that for [Te IV ] = 5 mM, the array of nanowires is less dense than for higher concentration values, suggesting that the density of nuclei is smaller. This is in accordance with the electrocrystallization rules proposed by Scharifker and Hills [51].

SAED analyses were performed and typical diffraction patterns are displayed in the Figure 2, a, b and c. All the data exhibit single crystalline phase and the diffraction spots can be indexed to the hexagonal phase of tellurium (JCPDS card No. 36-1452). Whatever the concentration value, the corresponding zone axis is compatible with a growth direction along the c axis.

This result is in good agreement with the shapes of nanowires synthesized for [Te IV ] = 5 mM, exhibiting a well-defined hexagonal shape.

Additional XRD crystallographic analyses of the deposits were performed (see Figure 3). The patterns reveal a preferential orientation along the [001] direction as discussed below.

Except the peaks related to the XRD sample holder and the Pt substrate, all the diffraction peaks can be indexed to the hexagonal phase of tellurium (JCPDS card No. 36-1452), confirming that the deposits are single phase. Sharpness of peaks suggests that the nanowires are well crystallized. In order to obtain a degree of preferential orientation, the Harris' texture coefficient (TC) was systematically calculated according to the following equation: Figure 4 shows that the texture coefficient of the (003) plane is higher than those of other (hkl) planes for Te IV concentrations of 20 mM and 25 mM. These data corroborate the SEM images (Figure 2a andb), proving that the deposits are composed of dense and [001] oriented nanowires. The [001] orientation is due to the highly anisotropic crystal structure of tellurium(space group P3121). In contrast, no marked orientation is observed for the deposits obtained at 5 mM, which can be explained by a random tilt of the nanowires on the substrate (Figure 2c).

Influence of applied potential

The influence of applied potential was studied for two concentration values: [Te IV ] = 5 mM (Figure 5 a-c), leading to the best aspect ratio at peak potential and [Te IV ] = 20 mM (Figure 5 In order to study the formation of HNs, ECD experiments were performed for different coulometric charge values (Q) at low Te IV concentration (5 mM) and high overpotential (Ediff), corresponding to the same conditions as Figure 5c. Figure 6 presents the SEM and TEM images of the obtained deposits. As previously, a statistic study from the analysis of 1D

nanostructures was carried out to determine their mean dimensions (length and diameter) as a function of Q. Corresponding values and calculated aspect ratio are presented in Figures 7a and7b, respectively.

At 100 mC.cm -2 (Figure 6a), the deposit covers the entire working electrode surface and is composed of nanowires. When Q increases up to 200 mC.cm -2 (Figure 6b), some longer 1D nanostructures appear above a sublayer of smaller nanowires. This phenomenon is more and more pronounced as Q increases. One can notice that the nanowires are tip-shaped (inset of Figure 6b) and grow in diameter as fast as in length as shown in Figure 7a. Consequently, the aspect ratio remains constant, 10.0 ± 1.5.

The first HNs appear for a coulometric charge of 400 mC.cm -2 (Figure 6c). By still increasing Q, the number of HNs increases (Figure 6, d-f), some of them being hexagonal shaped.

Lengths and diameters remain roughly proportional leading thus to the same aspect ratio as previously, 10.0 ± 1.5 (Figure 7a).

From the TEM images, it appears then that the nanostructures first grow as nanowires (Figure 6, a and b) and progressively turn into hollow nanostructures (Figure 6, c-f). Consequently, the structures are not hollow from one end to the other (see for example Figure 6, d-e). It seems then that there are several stages of growth: a first stage during which nanowires form, with a center part growing faster than the external parts, leading to tip-shaped structures, as was already observed in literature [START_REF] Kumar | Optical and electrical studies of vertically oriented tellurium nanowires arrays produced by template electrodeposition[END_REF][START_REF] Guo | Butterfly-scale architecture directed electrodeposition of Ag microband arrays for electrochemical detection[END_REF]. In the second part of the growth, 1D hollow nanostructures appear. According to Choi's [START_REF] Siegfried | Directing the architecture of cuprous oxide crystals during electrochemical growth[END_REF][START_REF] Choi | Shape control of inorganic materials via electrodeposition[END_REF] works, this formation might be due to a limitation by the diffusion of the precursor ions towards the center part of the growing structures. This depletion zone formed in the central part of each structure leads to a faster growth of external parts, resulting in hollow structures formation. The formation of HNs was also observed during hydrothermal synthesis of Te nanostructures by Mayers and Xia [START_REF] Mayers | Formation of tellurium nanotubes through concentration depletion at the surfaces of seeds[END_REF].

The authors observed the formation of nanotubes when the supply in Te IV precursor is ratelimiting whereas nanowires were formed otherwise.

The limitation by species supply seems then to favor the formation of HNs as the coulometric charge increases. To justify this assumption, additional experiments were conducted by using a Te soluble anode as counter electrode in the experimental conditions leading to HNs formation, namely a low tellurium salt concentration (5 mM) and a high overpotential (Ediff).

A typical deposit is presented in Figure 8.

One can observe on Figure 8 that only nanowires are formed. As opposite to the experiments conducted in the same growth conditions without the use of a soluble anode (see Figure 6f), no HNs are noticed. Using a soluble anode, the concentration of Te IV in the electrolyte is supposed to be constant during the electrolysis. This hypothesis is based on the presence of the anodic peak corresponding to Te(0) dissolution, observed by Szymczak et al. in the same electrolyte [START_REF] Szymczak | Electrodeposition of stoichiometric bismuth telluride Bi2Te3 using a piperidinium ionic liquid binary mixture[END_REF]. Consequently, the diffusion layer thickness remains constant during the whole experiments leading to constant diffusion flux.

Obtained nanowires are not very homogeneous in size and their density is low. However, the mean length is 5.0 ± 2.0 µm whereas without the use of a soluble anode, it is of 2.7 ± 1.0 µm.

The mean diameter remains the same that previously observed, 250 ± 100 nm. The use of a soluble anode allows then to obtain NWs instead of HNs, with larger length and equivalent diameter. Thus this additional result confirms that the growth of HNs is related to species supply.

The conducted experiments in this work show that the origin of the nanostructuration (NWs and HNs) is governed by the applied overpotential and the cationic concentration. It should be noted that, from the literature data [START_REF] Szymczak | Template-free electrodeposition of tellurium nanostructures in a room-temperature ionic liquid[END_REF][START_REF] Ma | Growth of tellurium nanowire bundles from an ionic liquid precursor[END_REF][START_REF] Steichen | Synthesis of trigonal selenium nanorods by electrodeposition from an ionic liquid at high temperature[END_REF][START_REF] Su | Electrodeposition of aluminium wires from the Lewis acidic AlCl3/trimethylamine hydrochloride ionic liquid without using a template[END_REF][START_REF] Leong | Electrochemistry on tin in the 1-ethyl-3methylimidazolium dicyanamide room temperature ionic liquid[END_REF][START_REF] Hsieh | Speciation of cobalt-chloride-based ionic liquids and electrodeposition of Co wires[END_REF][START_REF] Yang | Single-step large-scale and template-free electrochemical growth of Ni-Zn alloy filament arrays from a zinc chloride based ionic liquid[END_REF][START_REF] Hsieh | Electrodeposition of CuZn from chlorozincate ionic liquid : from hollow tubes to segmented nanowires[END_REF][START_REF] Hsieh | One-step electrochemical fabrication of nanoporous gold wire arrays from ionic liquid[END_REF][START_REF] Yang | Direct electrodeposition of FeCoZn wire arrays from a zinc chloride-based ionic liquid[END_REF], ionic liquids could play the role of capping agents. Besides, halide anions can also be responsible of the growth of nanostructures as observed by Al-Salman et al. [START_REF] Al-Salman | Template-free electrochemical synthesis of high aspect ratio Sn nanowires in ionic liquids: A general route to large-area metal and semimetal nanowires arrays?[END_REF]. In our work, it is difficult to discriminate between halide effect and IL ions effect, as the only Te(IV) salts that are soluble in TFSI-based ionic liquids are halide salts.

Conclusions

In this work, we have shown that it is possible to control the shape and the size of the deposited nanostructures by varying the experimental parameters that influence the deposition kinetics. It appears in particular that the applied potential has a very strong influence on the deposit morphology. Under charge transfer control, micrometer faceted grains are observed.

Under diffusional control, the deposits are composed of nanowires. In high overpotential conditions, nanowires turn progressively into hollow nanostructures due to a limitation by species supply. This phenomenon is more accentuated at low concentration.

This work shows that Te nanowires of high crystallinity can be obtained by electrodeposition in ionic liquids. They could be used as starting materials for the synthesis of core-shell structures. Hollow nanostructures synthesized in specific conditions are also interesting to create high performances thermoelectric materials. Indeed, this geometry could further decrease the lattice part of the thermal diffusion coefficient by phonon confinement.

  The morphology of the deposits was analyzed by SEM (Philips XL30) at 5 kV. TEM images of individual nanostructures were recorded by a Philips CM200 transmission electron microscope operating at an accelerating voltage of 200 kV. A statistic study from TEM images was carried out on a dozen nanostructures per image to evaluate their mean dimensions (length and diameter) for each experimental condition. The growth direction of the nanostructures was characterized by imaging and selected area electron diffraction (SAED). The diffraction spots were indexed according to the JCPDS card No. 36-1452 of the tellurium hexagonal phase. The crystalline structure of the deposits was investigated by X-ray powder diffraction measurements (XRD, Bruker D8Advance) with Cu Kα radiation (λ =1.5406Å) at 40 kV and 40 mA.

Figure 2

 2 presents the SEM images of Te deposits obtained for a coulombic charge of 2 C.cm -² and [Te IV ] = 25 mM (a), 20 mM (b) and 5 mM (c).

  hkl) and I 0 (hkl) are the measured intensity of a generic peak observed in the experimental XRD and the literature value from the data base (JCPDS card No. 36-1452), respectively, and N the number of diffraction peaks considered in the analysis. The TC results are presented in Figure 4.

Figure 1 : 19 Figure 3 :

 1193 Figure 1: Cathodic behaviour of Te IV in EOPipTFSI:EOPipBr 95:5 (mol%). [TeCl4] = 20 mM. Pt-coated glass slide substrate. T = 100°C, scan rate 5 mV.s -1 .

Figure 4 :

 4 Figure 4: Texture coefficients of tellurium nanowires according to (hkl) planes,

Figure 5 :

 5 Figure 5: SEM images of Te nanostructures deposited in EOPipTFSI:EOPipBr 95:5 (mol%) for different applied potential values. T = 100°C, Q = 2 C.cm -2 , Pt-coated glass slide.

Figure 6 :

 6 Figure 6: SEM images of Te nanostructures electrodeposited for different coulometric charge values in EOPipTFSI:EOPipBr 95:5 (mol%), Edep = Ediff, [Te IV ] = 5 mM, (a) 100 mC.cm -2 , (b) 200 mC.cm -2 , (c) 400 mC.cm -2 , (d) 800 mC.cm -2 , (e) 1 C.cm -2 and (f) 2 C.cm -2 , T = 100°C, Pt-coated glass slide

Figure 7 :

 7 Figure 7: Evolution of a) length and diameter and b) aspect ratio of Te one-dimensional

 d-f) for which dense arrays of nanowires were deposited. Three potential zones were tested: the half-wave potential value (E½), corresponding to the potential where the current is half the C1 peak current, the C1 peak potential value (Ep), and a potential value (Ediff) more negative than Ep (Ediff (V) = Ep -0.15). The deposits obtained at peak potential value previously presented in Figure 2 are reminded for easier comparison (Figure 5, b ande).

Whatever the concentration, the deposits obtained at E1/2 are compact, composed of faceted micrometric grains (Figure 5, a andd). SAED analysis shows that these structures are single crystalline and grow preferentially along the [001] direction, as previously observed. One can see on Figure 5 that the diameter of the nanostructures decreases with applied overpotential; this phenomenon is more accentuated comparing the E1/2 deposits and the Ep ones. This is in agreement with electrocrystallization principles, that is, the critical size of nuclei is inversely proportional to the deposition overpotential [START_REF] Paunovic | Fundamentals of electrochemical deposition[END_REF].

When the applied potential is much more negative than C1 peak potential (Figure 5, c and f), the SEM observations of the deposits show of a compact sublayer composed of small grains and a mix of hollow nanostructures (HNs) and nanowires of several micrometers in length. It has already been shown that, using ionic liquids, HNs can been obtained in specific electrochemical conditions (Sn [START_REF] Leong | Electrochemistry on tin in the 1-ethyl-3methylimidazolium dicyanamide room temperature ionic liquid[END_REF], CuZn [START_REF] Hsieh | Electrodeposition of CuZn from chlorozincate ionic liquid : from hollow tubes to segmented nanowires[END_REF] and Cu10Sn3 [START_REF] Hsieh | Direct template-free electrochemical growth of hexagonal CuSn tubes from an ionic liquid[END_REF]).

One can also observe on Figure 5 that the amount of hollow nanostructures increases as the concentration of Te IV decreases. This suggests that their formation results from an insufficient species supply. The formation of hollow nanostructures was already observed by Hsieh et al.

when studying the synthesis of CuZn nanostructures by template-free electrodeposition [START_REF] Hsieh | Electrodeposition of CuZn from chlorozincate ionic liquid : from hollow tubes to segmented nanowires[END_REF], but, in their case, the reverse observation was made : HNs were obtained for low overpotential values under charge transfer control, whereas nanowires were obtained at higher overpotential when mass transport is controlling the growth.


J Vincent
email: vincent@laplace.univ-tlse.fr

P Borderies

J R Poirier

V Gobin

3D-FDTD-based hybridisation technique for wave interaction with scatterers in layered media

A simple method for the computation of 3D scatterers remotely located from the source with layered media is presented. This domain decomposition method is based on the total-field/scattered-field formulation of the 3D finite-difference in time-domain method. The electromagnetic field along the propagation path is not directly involved; it focuses only on the electromagnetic field in the surroundings of the scatterer. The incident electromagnetic field radiated from the source is computed separately with an adaptive algorithm and then applied on the Huygens' surface, leading to the proposed hybridisation technique.

Introduction: Computing the electromagnetic field far from the source with the finite-difference in time-domain (FDTD) method can be hard to achieve because of the numerical dispersion in large-scale problems. In [START_REF] Zhou | LF ground-wave propagation over irregular terrain[END_REF][START_REF] Zhou | Application in low-frequency ground-wave propagation of parallel FDTD based GPU[END_REF], a very large 2D computational domain is developed to evaluate the electric field along a propagation path. However, this technique does not take into account the effects of the 3D topography because of the amounts of space and time needed.

The domain decomposition method based on the total-field/scatteredfield formulation (TF/SF) [3] is a solution which can overcome the large-scale problem. This formulation was first used to illuminate scattering bodies in free space with a plane wave [START_REF] Merewether | On implementing a numeric Huygens' source scheme in a finite difference programme to illuminate scattering bodies[END_REF][START_REF] Umashankar | A novel methjod to analyze electromagnetic scattering of complex objects[END_REF]. The incident electromagnetic field is known and applied on a 3D closed surface named Huygens' surface. Later the ground is included in [START_REF] Huang | Ground-penetrating radar antenna modeling[END_REF] and in [START_REF] Demarest | FDTD modeling of scatterers in stratified media[END_REF] where an incident plane wave creates transmitted and reflected plane waves by hitting the surface.

In this Letter, the TF/SF formalism is used in order to avoid the propagation of the electromagnetic field radiated by a source and to focus on the interaction with 3D scatterers which can be above, on or under the surface of a layered media. The incident electromagnetic field is known everywhere in space and time or can be computed with another method or software and then applied on the six faces of the Huygens' surface.

TF/SF FDTD formulation with a media: From the linearity of Maxwell's equations, one can decompose the physical total electric field E TOT and magnetic field H TOT as

E TOT = E INC + E SCAT , H TOT = H INC + H SCAT ( 1
)
where E INC and H INC are the incident field components. As stated above, their values are known everywhere in the computation space and at each time step. These are the fields that would exist in all concerned media, without any other object in the computation space. E SCAT and H SCAT are the unknown SF components created by the introduction of at least one scatter. Thus, the TF/SF FDTD formalism creates two distinct regions in the computation space bounded by the Huygens' surface (see Fig. 1). Inside the surface, the total electromagnetic field is propagating in the region, and outside the Huygens' surface is the volume where only the SF is taken into account. As an example, Fig. 2 shows the case of a vertical current element and an antenna receiver at the bottom of a valley, computed with an infinite half-space with the two FDTD formalisms. The receiver is located at several wavelengths from the source. Fig. 2a shows the entire domain which is meshed and computed with the usual 3D FDTD algorithm. The TF/SF formulation is represented in Figs. 2b andc As the source is remotely located from the scattering object in the TF/ SF formalism and the propagation path is not computed, the computation time and the memory space needed are obviously reduced. Nevertheless, the values of the incident electromagnetic field applied on all the mesh points of the Huygens' surface at all time steps have to be reasonably accurate.

Incident wave source for TF/SF FDTD formalism with a media: From (1) and from the free-space TF/SF FDTD formulation given in equations (5.48a) to (5.59b) of [3] (pp. 204-210), one can find the assumed known correction terms used in the TF/SF FDTD formulation with a media. Caution should be taken with the coefficients since the dielectric parameters of the media (ε 0 ε r , σ, μ 0 μ r , σ*) are involved in the computations. In this Letter, only media with relative permeability μ r = 1 and without magnetic loss (σ* = 0) are considered.

Obtained from (5.48a) and (5.51a) of [3], ( 2) and (3) show the coefficients which have to be included on one face of the Huygens' surface, in order to consider a media in the TF/SF formulation. The 22 other formulas which correct the electromagnetic field at each halved time step on the whole Huygens' surface can be easily written from the two following equations and from [3]

E x,TOT n + 1 i, j 0 , k = E x,TOT n + 1 i, j 0 , k - 2Dt (21 + sDt)Dy H z,INC n + 1 2 i, j 0 -1/2, k (2)
H z,SCAT n + 1/2 i, j 0 -1/2, k = H z,SCAT n + 1/2 i, j 0 -1/2, k - Dt m 0 Dy E x,INC n i, j 0 , k (3)
In ( 2) and ( 3), E INC and H INC are the known incident fields applied on the whole Huygens' surface at each halved time step. When their values are in the time domain, the pulse is introduced directly in the algorithm.

In frequency domain (e.g. the electromagnetic field radiated at a frequency by an antenna, obtained with the boundary integral equation method), one has to use ( 4) and ( 5) for a time harmonic notation e +jωt :

E FDTD x, y, z, t = < E BIE (x, y, z) • e +jvt • 1 -e -at (4)
H FDTD x, y, z, t = < H BIE (x, y, z) • e +jvt • 1 -e -at (5)
where F FDTD (respectively F BIE ) is the electromagnetic field applied in FDTD (respectively boundary integral equation) and the term (1-e -αt ) is used to prevent the rapid changes of values at the beginning of the algorithm. The parameter α is chosen by the user such that the signal has reached its maximum amplitude after a few initialisation periods. The incident electromagnetic field in both the air and the ground is computed by a numerical integration of Sommerfeld's integrals with an adaptive algorithm [START_REF] Vincent | Modelling the effects of realistic environments in the nearfield and farfield of low frequency antennas with 3D FDTD method[END_REF]. These values are applied on the Huygens' surface surrounding the receiver. The volume inside the Huygens' surfaces is a 60 cells side cube centred in the 3D FDTD domain. The parameter α is chosen such that 99% of the signal is reached after five periods (i.e.

[1-exp(-α × 5/f 0 )] = 0.99). Uniaxial-PMLs are used to terminate the 3D computational domain [START_REF] Gedney | An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices[END_REF] and the thin wire (r = 5 mm) is computed with the wire formalism of Holland and Simpson [START_REF] Holland | Finite-difference analysis of EMP coupling to thin struts and wires[END_REF]. As shown in Fig. 3, the current observed in the middle of the wire receiver needs seven periods to be considered as a sinusoidal waveform. These values are compared with the current obtained by computing the total scene with the boundary equation method of the solver FEKO. The results are in good agreement (relative error of 1.59% in the sinusoidal waveforms). Conclusion: A 3D-FDTD-based hybridisation technique is presented for wave interaction with scatterers in layered media. It rests on the SF-TF formalism with a media. The incident electromagnetic field is known everywhere in space at each time step and applied on the Huygens' surface. This domain decomposition technique is very useful to avoid large domain computations and to focus on the object located upon or inside the layered media.

1 :Fig. 1

 11 Fig. 1 Representation of TF/SF FDTD formalism. Antenna receiver is located inside volume bounded by Huygens' surface (dashed line) in region 2 (TF). Region 1 (SF) is bounded by boundary conditions (UPML)

Fig. 2

 2 Fig. 2 FDTD computations of radiating current element and antenna receiver at bottom of valley, with regular algorithm and with TF/SF formalism. Scattering object is placed inside the Huygens' surface (dashed line) a Entire domain computation with usual FDTD algorithm b Geometry considered in TF/SF formalism, incident electromagnetic field is applied on Huygens' surface c Scattering object considered in TF/SF formalism with media

  Computation of a receiver antenna far from the source: As an example, a 1125 m long vertical wire used as a receiver is located at 375 m over a medium wet flat ground with relative permittivity ε r = 15 and conductivity σ =10 -3 Sm -1 . The source is remotely located at 23.75 km from this receiver and 25 m above the ground. It is a 25 m long infinitesimal vertical current element (I 0 = 325.76 × 10 -6 + j × 51.581 A) radiating at the frequency f 0 = 100 kHz. The FDTD computation space is a 100 cells side cube with the mesh size Δx = Δy = 25 m and Δz = 75 m. The ground surface (z = 0 m) is located at the middle of the z-directed cells.

Fig. 3

 3 Fig. 3 Current in middle of antenna receiver, located at 23.75 km from source

Acknowledgments: This work was supported by Direction Générale de l'Armement (DGA) and Office National d'Études et de Recherches Aérospatiales (ONERA).