Laurent Gosse
email: l.gosse@ba.iac.cnr.it

REDHEFFER PRODUCTS AND NUMERICAL APPROXIMATION OF CURRENTS IN ONE-DIMENSIONAL SEMICONDUCTOR KINETIC MODELS

Keywords: Discrete ordinates, Case's elementary solutions, Hamiltonian-Preserving scheme, Kinetic model of semiconductors, Redheffer star-product, Scattering matrix, Well-Balanced scheme AMS subject classifications. 65M06, 82D37

Redheffer products and numerical approximation of currents in one-dimensional semiconductor kinetic models

1. Introduction and preliminaries. Our main goal is to compute a reliable numerical kinetic density of mobile charge carriers (either electrons or holes) in the vicinity of a p|n junction represented on an interface of the 1D computational grid. The following simplifying assumptions are made throughout this text:

1. carriers recombination is not modeled, we treat only one equation rendering for the dynamics of the majority carriers (i.e. electrons), 2. the dielectric value is depends on impurity concentration (i.e. doping), 3. the carriers mobility can depend on both the doping and the electric field. Numerical issues come from the subtle balance occurring around a p|n junction:

• conduction electrons from donors (resp. holes for acceptors) first are attracted by opposite ions, hence diffuse out of the n (resp. p) area, • until the negative ionization of acceptors (resp. positive ionization of donors) creates a potential barrier stopping this motion. • This "spiky" electric field is supported inside a limited region, located around the junction, where electrons (holes) aren't available for current conduction. It is called a "depletion layer", its length depends on impurity concentration and applied voltage: this area can be identified as an intrinsic semiconductor. The channel of a n + nn + device gets filled with carriers until depletion layers corresponding to both the forward bias and the specificities of the junctions are correctly created and merged.

1.1. The relaxation-time approximation (RTA) kinetic model. A derivation of the "RTA approximation" of the low-density quadratic collision term of Vlasov-Boltzmann (also called Boltzmann-Poisson) equation is presented in [START_REF] Jüngel | Transport Equations for Semiconductors[END_REF], pages 85-86. Here, following [START_REF] Carrillo | Computational macroscopic approximations to the onedimensional relaxation-time kinetic system for semiconductors[END_REF], we simplify further this multi-D weakly nonlinear kinetic equation by restricting our study to simplified devices in one dimension by assuming that most of the salient features of the electronic transport are given in the direction parallel to the force field induced by Poisson's potential. Hence we are led to considering the following kinetic equation, already scrutinized in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapters 11 and 15:

∂ t f + v∂ x f + E∂ v f = 1 τ (x) M θ (v) R f (t, x, v ′ )dv ′ -f , t, v ∈ R + * × R, (1.1)
where f (t, x, v) stands for the kinetic density depending on a time, space and velocity variables, τ is a relaxation time and E(t, x) = -∂ x ϕ(t, x) is the electric field given by:

-λ 2 (x)∂ xx ϕ = ρ -ρ D (x), ρ(t, x) = R f (t, x, v)dv. (1.2)
The parameter λ 2 > 0 is the scaled Debye length which renders the electrons screening distance in the material and ρ D (x) > 0 is the doping profile, i.e. the concentration of donor impurities molten in the semiconductor crystal. Generally, the equations (1.1), (1.2) are set in a bounded space domain, x ∈ (-1, 1) and boundary conditions read:

f (t, x = ∓1, ±|v|) = ρ D (∓1)M θ (v), ϕ(t, x = -1) = 0, ϕ(t, x = 1) = -V, (1.3) M θ (v) = 1 √ 2πθ exp(- v 2 2θ ), v ∈ R,
being the Maxwellian distribution at the lattice temperature θ > 0 and V > 0 a voltage bias. Such inflow boundary data render ohmic contacts at each device's edge. It's essentially the model advocated in e.g. [START_REF] Ben Abdallah | The Milne problem for high field kinetic equations[END_REF][START_REF] Carrillo | Computational macroscopic approximations to the onedimensional relaxation-time kinetic system for semiconductors[END_REF] and studied theoretically in [START_REF] Ben Abdallah | Diffusion Approximation for the one dimensional Boltzmann-Poisson system[END_REF][START_REF] Fuchs | Asymptotic analysis of the degenerate Boltzmann-Poisson system for semiconductors[END_REF].

1.2. Overview of the "elementary solutions" literature. According to e.g. [START_REF] Ballestra | Numerical problems in semiconductor simulation using the hydrodynamic model: a second-order finite difference scheme[END_REF], a main goal when it comes to performing numerical simulations of the aforementioned weakly nonlinear system is the determination of the so-called "Current-Voltage relation", that is to say the (constant) macroscopic flow J(t, x) = R vf (t, x, v)dv at steady-state expressed as a function of the biasing value V . Clearly, this asks for a numerical method which is able to stabilize in large times onto a stationary kinetic density f endowed with a constant moment of order 1: see [START_REF] Ben Abdallah | Relative entropies for kinetic equations in bounded domains (irreversibility, stationary solutions, uniqueness)[END_REF] for theoretical details. Such a requirement is reminiscent of the so-called "well-balanced property" for shallow water equations, for which one seeks numerical stationary flows endowed with a constant discharge, even in presence of a strongly-varying topography term [START_REF] Bouchut | A subsonic-well-balanced reconstruction scheme for shallow water flows[END_REF].

This leads to the question about extending the well-balanced framework from quasi-linear hyperbolic systems of balance laws toward collisional kinetic equations: instead of localizing a topography term at each interface of a Cartesian computational grid, one is now led to (formally) concentrate, by means of Dirac masses [START_REF] Glimm | An S-matrix theory for classical nonlinear physics[END_REF], both the collision process, and the Vlasov acceleration term as well. The handling of linear collisional processes in a kinetic time-dependent setup is described in [START_REF] Gosse | Transient radiative transfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF][START_REF] Gosse | Well-balanced schemes using elementary solutions for linear models of the Boltzmann equation in one space dimension[END_REF][START_REF] Gosse | A well-balanced scheme for kinetic models of chemotaxis derived from onedimensional local forward-backward problems[END_REF] and reviewed in Part II of [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]: it consists in building at each interface a "scattering matrix" which relates the incoming states f n j-1 (|v|), f n j (-|v|) to a set of outgoing ones, f n j-1 2 (±|v|), cf. [START_REF] Glimm | The interaction of nonlinear hyperbolic waves[END_REF] and Fig. 2.1. This scattering matrix is usually obtained from the exact integration of the stationary equations (more astute manners exist, see [START_REF] Bellman | An introduction to invariant imbedding[END_REF]).

In the context of 1+1-dimensional linear kinetic models, this immediately raises the issue of deriving exact solutions of boundary-value problems for the corresponding stationary equations: such a program was thoroughly studied after Case's influential paper [START_REF] Case | Elementary solutions of the transport equation and their applications[END_REF] on "elementary solutions", see e.g. [START_REF] Aamodt | Useful identities for half-space problems in linear transport theory[END_REF][START_REF] Beals | An abstract treatment of some forward-backward problems of transport and scattering[END_REF][START_REF] Case | Linear transport theory[END_REF][START_REF] Cercignani | The method of elementary solutions for kinetic models with velocity dependent collision frequency[END_REF][START_REF] Kaper | Spectral Methods in Linear Transport Theory[END_REF] which can be presented as "spectral Green's functions" [START_REF] Barros | A Numerical Method for One-Group Slab-Geometry Discrete Ordinates Problems[END_REF]. Thanks to Siewert's works [START_REF] Barichello | Unified Solutions to Classical Flow Problems Based on the BGK Model[END_REF][START_REF] Barichello | A discrete-ordinates solution for a non-grey model with complete frequency redistribution[END_REF][START_REF] Siewert | Efficient eigenvalue calculations in radiative transfer[END_REF], this was converted into a numerical algorithm, originally designed for steady-state computations involving least-square procedures, referred to as "Analytical Discrete-Ordinates". The extension to unsteady equations was first performed in [START_REF] Gosse | Transient radiative transfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF] just by observing that the well-balanced approach is equivalent to concentrating the scattering process of collisions at each interface of the grid: this led to high-quality numerical results for various models [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]. One shortcoming at this stage was the difficulty of including a Vlasov-type acceleration term: indeed, several authors tried to derive a set of "elementary solutions" for simple models of collisional Vlasov equations, see [START_REF] Ch | Half-space problem of the Boltzmann equation for charged particles[END_REF][START_REF] Paveri-Fontana | A Neutral Gas Model for Electron Swarms[END_REF][START_REF] Toepffer | Analytical results for the Boltzmann equation[END_REF]. Unfortunately, left apart the special case of Vlasov-Fokker-Planck models (see [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapters 12-13), these spectral solutions were far too intricate to be used for producing an efficient numerical routine. An alternative construction was set up in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapter 11, and led to an efficient algorithm partially built on a high-field asymptotic distribution [START_REF] Frosali | Conditions for runaway phenomena in the kinetic theory of swarms[END_REF][START_REF] Ben Abdallah | The Milne problem for high field kinetic equations[END_REF]. However, this derivation may be felt as somewhat rigid because such an asymptotic density is rarely exactly computable. A remedy may consist in "marrying" both the exact resolutions of a Vlasov equation by means of the Hamiltonian-Preserving method, see [START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF] also [START_REF] Abdallah | A Hybrid Kinetic-Quantum Model for Stationary Electron Transport[END_REF][START_REF] Greengard | A boundary value problem for the stationary Vlasov-Poisson equations: The plane diode[END_REF], which can be rephrased as the action of a local (stochastic) scattering matrix and a linear collisional one. An issue is that one doesn't obtain a scattering matrix out of the matrix product of two of them; the key is to consider instead the Redheffer's star-product, [START_REF] Redheffer | On the relation of transmission-line theory to scattering and transfer[END_REF][START_REF] Timotin | Redheffer products and characteristic functions[END_REF].

1.3. Organization of the paper. Our main objective is to describe a practical algorithm which combines the nice properties of both the Hamiltonian-Preserving (HP) treatment of the Vlasov acceleration term [START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF][START_REF] Wen | The L 1 -Error Estimates for a Hamiltonian-Preserving Scheme for the Liouville Equation with Piecewise Constant Potentials[END_REF] and the scattering matrices obtained through the so-called Analytical Discrete-Ordinates (ADO) method [START_REF] Barichello | Unified Solutions to Classical Flow Problems Based on the BGK Model[END_REF][START_REF] Barichello | A discrete-ordinates solution for a non-grey model with complete frequency redistribution[END_REF][START_REF] Siewert | Efficient eigenvalue calculations in radiative transfer[END_REF], thanks to the Redheffer product of scattering matrices, in the formulation of [START_REF] Timotin | Redheffer products and characteristic functions[END_REF]. The paper is organized as follows: §2 is devoted to the derivation of the scattering matrices corresponding to the linear collision process (see §2.1) and the electric acceleration (see §2.2), respectively. The Redheffer product is presented in §3, together with a simple fix to ensure current preservation at the interfaces (see §3.2). A maximum principle for the resulting Godunov scheme is deduced in Proposition 1 (see §3.3). The drifted-Maxwellian simplification for handling source and drain heavily doped regions is recalled in §4, together with a correct manner to implement inflow boundary conditions. Then, in §5, some numerical results are provided at numerical steady-state in presence of discontinuous parameters. In §6, a more involved situation, where the relaxation parameter τ depends also on the electric field in order to render pinching effects and mobility saturation is studied. After concluding remarks delivered in §7, an Appendix A displays a short numerical study of weakly collisional Landau damping, following [START_REF] Degond | Particle Simulations of the Semiconductor Boltzmann Equation for One Dimensional Inhomogeneous Structures[END_REF] (and [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], pages 225-227).

2. Scattering matrices for acceleration and collisions. Before entering the matters core, let us state a few relevant notions of linear algebra: Definition 1. A N × N real matrix M is called stochastic (or "Markov matrix") if all its entries are non-negative and the sum of each line's components equals 1:

∀1 ≤ i ≤ N, N j=1 M i,j = 1.
Lemma 2.1. Let the N × N real matrix M be stochastic, then:

1. it preserves non-negativity, 2. 1 is an eigenvalue associated to

1 := (1, 1, ..., 1) ∈ R N , 3. M ∞ = 1.
Proof. Consider any X ∈ R N + , the i th component of M X is N j=1 M i,j X j which is non-negative as a sum of non-negative terms (any line of M realizes a convex combination). If X = 1, any component of M X reads N j=1 M i,j = 1 hence 1 is the eigenvector associated to the eigenvalue 1. Lastly, for any X ∈ R N with

X ∞ = 1, the i th component of M X satisfies N j=1 M i,j X j ≤ N j=1 M i,j |X j | ≤ 1.
Then M ∞ ≤ 1 and since 1 is eigenvalue associated to 1, M ∞ = 1.

2.1. Scattering matrix for linear (RTA) collision term. Hereafter, we are given a Cartesian computational grid in t, x determined by the positive parameters ∆t, ∆x; the shorthand notation is used, x j = j∆x, t n = n∆t, the control cell located around x j is C j = (x j-1 2 , x j+ 1 2 ). Moreover, a maximal value of the velocity modulus is fixed, v M AX > 0 and an even (Gaussian) quadrature is set on the interval (-1, 1). As it is symmetric and 0 doesn't belong to it, it can safely be multiplied by v M AX ,

{v k , ω k } k=1,...,2K , v k = -v 2K+1-k = 0, ω k = ω 2K+1-k , K ∈ N. (2.1)
The shorthand notation will be used frequently:

V = {v k } k=K+1,...,2K ∈ (0, v M AX ) K Ω = {ω k } k=K+1,...,2K ∈ R K + .
The numerical schemes advocated in Part II of [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF] hinge on deriving a Godunov scheme for the "localized equation" which reads, when only collisions are included:

∂ t f + v∂ x f = ∆x τ j∈Z M θ (v) R f (t, x, v ′ )dv ′ -f δ(x -x j-1 2 ).
This is essentially the approach suggested in [START_REF] Glimm | An S-matrix theory for classical nonlinear physics[END_REF]: the presence of this "localized source term" induces a static discontinuity at each interface separating 2 control cells C j-1 and C j and implies that the Godunov scheme reads, for each 0

< v k ∈ V ⊂ R K + : f n+1 j (v k ) = f n j (v k ) -∆t ∆x v k f n j (v k ) -f n j-1 2 (v k ) , f n+1 j (-v k ) = f n j (-v k ) + ∆t ∆x v k f n j+ 1 2 (-v k ) -f n j (-v k ) .
Let S τ stand for the scattering matrix which relates incoming/outgoing states at each interface x j-1 2 : the former numerical scheme rewrites in a less usual way,

f n+1 j (V) f n+1 j-1 (-V) = 1 -V ∆t ∆x f n j (V) f n j-1 (-V) + V ∆t ∆x S τ f n j-1 (V) f n j (-V) . (2.
2)

The presentation is now centered on each interface x j-1 2 of the grid (like in [START_REF] Delarue | Probabilistic analysis of the upwind scheme for transport equations Arch[END_REF]), instead of the former one which is centered on each cell C j . Observe that, since 0 < v k ∈ V, the states f n j-1 (v k ), f n j (-v k ) approach in the sense of Glimm [START_REF] Glimm | The interaction of nonlinear hyperbolic waves[END_REF] the Dirac mass located at x j-1 2 . Previous experience implies that the entries of the scattering matrix are known from the solving of the (forward/backward) boundaryvalue problem with inflow data for the stationary equation of the original problem,

v∂ x f = 1 τ M θ (v) R f (x, v ′ )dv ′ -f , x ∈ (0, ∆x). Now, since M θ L 1 (R) = 1
, it is known from [START_REF] Aamodt | Useful identities for half-space problems in linear transport theory[END_REF][START_REF] Cercignani | The method of elementary solutions for kinetic models with velocity dependent collision frequency[END_REF] that if ḡ(x, v) solves a similar equation, namely

v∂ x ḡ = R M θ (v ′ )ḡ(x, v ′ )dv ′ -ḡ, then f (x, v) = M θ (v)ḡ( x τ , v
). Solving the forward/backward problem for this last equation is easy according to Case's theory of "elementary solutions" that we don't recall; we refer instead to [START_REF] Bart | Linear integral equations of the third kind[END_REF][START_REF] Beals | An abstract treatment of some forward-backward problems of transport and scattering[END_REF][START_REF] Case | Elementary solutions of the transport equation and their applications[END_REF][START_REF] Case | Linear transport theory[END_REF][START_REF] Kaper | Spectral Methods in Linear Transport Theory[END_REF]. According to [START_REF] Barichello | Unified Solutions to Classical Flow Problems Based on the BGK Model[END_REF][START_REF] Barichello | A discrete-ordinates solution for a non-grey model with complete frequency redistribution[END_REF][START_REF] Siewert | Efficient eigenvalue calculations in radiative transfer[END_REF], one can derive the 2K × 2K matrix S τ by: 1. determining the eigenmodes ν k (sometimes called "constants of separation") by considering the approximate stationary equation in the following form,

v∂ x f (x, v) + f (x, v) = vMAX 0 M θ (v ′ ) f (x, v ′ ) + f (x, -v ′ ) dv ′ . (2.3)
Analogously with the continuous case, the separation variable ν is introduced:

f (x, v) = ϕ(ν, v) exp(-x/ν).
Plugging into (2.3) and taking the quadrature rule into account yields:

1 ∓ v k ν k ϕ(ν, ±v k ) = 2K ℓ=1+K ω ℓ M θ (v ℓ ) ϕ(ν, v ℓ )+ϕ(ν, -v ℓ ) , k ∈ {1, ..., K}.
We denote Φ ± (ν) = ϕ(ν, ±v k ) k∈{1+K,...,2K} thus an eigenvalue problem arises:

± V ν Φ ± (ν) = Φ ± (ν) -(ΩM θ (V)) ⊗ Φ ± (ν) -(ΩM θ (V)) ⊗ Φ ∓ (ν). (2.4)
In [START_REF] Barichello | A discrete-ordinates solution for a non-grey model with complete frequency redistribution[END_REF][START_REF] Siewert | Efficient eigenvalue calculations in radiative transfer[END_REF], Barichello, Siewert and Wright present an astute manner to recast this eigenvalue problem into a simpler formulation, which allows to prove the "interlacing property": 0

< v 1 < ν 1 < v 2 < ν 2 < ... < v K < ν K . 2.
assembling the matrices M, M based on the "elementary solutions". Since ḡ(x, v) = 1, x -v are obviously solutions associated to the zero-eigenvalue, one deduces that f (x, v) = M θ (v), M θ (v)( x τ -v) are too (they are sometimes called "linear diffusion solutions" [START_REF] Fisch | Separating variables in two-way diffusion equations[END_REF]). The other solutions are associated with finite eigenmodes: ḡ

(x, v) = exp(-x/ν) 1-v/ν , thus f (x, v) = M θ (v) 1-v/ν exp(-x τ ν ). f (x, ±v k ) ≃ M θ (v k ) α + β x τ ∓ v k + E x τ , ±v k , ν , v k ∈ V.
E is a finite superposition of damped modes involving ν, the "Knudsen layers":

E(x, v, ν) = K-1 ℓ=1 A ℓ 1 -v/ν ℓ exp(-x/ν ℓ ) + B ℓ 1 + v/ν ℓ exp(x/ν ℓ ) . (2.5)
In order to compute the 2K outgoing states f (∆x, v k ), f (0, -v k ), v k ∈ V, we must first deduce the set of 2K coefficients α, β, A, B from the supplied boundary data f (0, v k ), f (∆x, -v k ), and then exploit them to derive the sought values. This leads to inverting a 2K × 2K matrix M depending on ν ∈ R K-1

+ M = 1 -V ⊗ ν -1 -1 1 R K 1 + V ⊗ ν -1 -1 exp(-∆x τ ν ) -V 1 + V ⊗ ν -1 -1 exp(-∆x τ ν ) 1 R K 1 -V ⊗ ν -1 -1 ∆x τ + V ,
and to defining the complementary matrix:

M = 1 -V ⊗ ν -1 -1 exp(-∆x τ ν ) 1 R K 1 + V ⊗ ν -1 -1 ∆x τ -V 1 + V ⊗ ν -1 -1 1 R K 1 -V ⊗ ν -1 -1 exp(-∆x τ ν ) V .
3. if D stands for the 2K×2K diagonal matrix which entries are {M θ (V), M θ (V)},

S τ = (D M ) (DM ) -1 (2.6)
is the expression of the desired scattering S-matrix used in (2.2), see Fig. 2.1.

f n j-1 (|v|) fj-1 2 (-|v|) fj-1 2 (|v|) f n j (-|v|) f n j (|v|) f n j-1 (-|v|) fj-1 2 (±|v|) = S τ f n j-1 (|v|), f n j (-|v|) f n j-1 (|v|) f n j (-|v|)
(t, x, v) := f (t,x,v) M θ (v) , (2.2)-(2.6) becomes g n+1 j (V) g n+1 j-1 (-V) = 1 -V ∆t ∆x g n j (V) g n j-1 (-V) +V ∆t ∆x M M -1 g n j-1 (V) g n j (-V) (2.7)
as one can easily check on a computer that for any τ, K > 0, the product M M -1 is stochastic. However, proving rigorously that this property holds uniformly in τ, K is an open problem, except for the simple case K = 1 (see [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapter 9, for details).

2.2. Scattering matrix for the Vlasov acceleration term. The derivation of an approximation of the scattering matrix S E for the linear Vlasov equation,

∂ t f + v∂ x f + E∂ v f = 0, E ∈ R, (2.8)
follows exactly from the same ideas as before, except that it appears simpler as one doesn't need "elementary solutions". The exact solution of a inflow boundary-value (forward/backward) problem for the stationary equation, v∂ x f + E∂ v f = 0 is given simply by "shifting velocities" as it is explained in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF][START_REF] Greengard | A boundary value problem for the stationary Vlasov-Poisson equations: The plane diode[END_REF][START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF][START_REF] Wen | The L 1 -Error Estimates for a Hamiltonian-Preserving Scheme for the Liouville Equation with Piecewise Constant Potentials[END_REF]. This construction (called "Scheme I" in [START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF]) furnishes a scattering matrix for Hamiltonian-Preserving by means of linear interpolation: numerical difficulties can arise because of the possibly weak variance of the probability density f , see e.g. [START_REF] Salton | The Use of the Central Limit Theorem for Interpolating in Tables of Probability Distribution Functions[END_REF]. Once the 2K × 2K matrix S E is known, the numerical treatment of (2.8) simply consists in iterating (2.2) with S E instead of S τ . Hereafter, we recall how to construct S E :

1. one first observes that it is enough to consider E < 0, because the stationary equation

v∂ x f + E∂ v f = 0 is invariant if both v, E change signs. This means: f n j-1 2 (-V) f n j-1 2 (V) = S -E f n j (-V) f n j-1 (V)
.

2. assuming that E < 0, one loops on the indexes k ∈ {1, 2, ..., K}:

• Let v = (v k ) 2 -2E∆x: k = max ℓ {v ℓ ≤ v} and α k = v-vk vk +1 -vk , (S E ) k, k = 1 -α k (S E ) k, k+1 = α k . • If (v k ) 2 + 2E∆x ≤ 0, then (S E ) K+k,k = 1. • Let v = (v k ) 2 + 2E∆x, k = max ℓ {v ℓ ≤ v}, and α -k = v-vk -1 vk-vk -1 : (S E ) K+k,K+ k-1 = 1 -α -k , (S E ) K+k,K+ k = α -k .
3. The former linear interpolation doesn't systematically preserve neither the total mass nor the current. Another choice is to replace it by:

α k = ω k v k -ω kv k ω k+1 v k+1 -ω kv k (S E ) k, k = 1 -α k (S E ) k, k+1 = α k .
The resulting matrix isn't always stochastic: it can contain negative entries.

In the paper [START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF], the so-called "Scheme II" appears to be built in such a way so as to preserve the currents across the scattering process (see Fig 2 .2). Its implementation seems quite heavy (see pages 294-296) though. fj-1 2 (|v|) (outgoing state)

f n j (-|v|) (incoming state) f n j-1 (|v|) (incoming state) fj-1 2 (-|v|) (outgoing state) Fig. 2.2.
Illustration of the scattering process for Vlasov with E < 0 2.3. Legendre polynomials and vanishing moments. Let us fix an interface x j-1 2 ∈ R, and following e.g. [START_REF] Beals | An abstract treatment of some forward-backward problems of transport and scattering[END_REF], consider the incoming/outgoing states:

f in (v) = f n j (-|v|) f n j-1 (|v|) , f out (v) = fj-1 2 (-|v|) fj-1 2 (|v|)
.

By consistency with the stationary kinetic problem, we want mass conservation,

2K k=1 ω k f in (v k ) -f out (v k ) = 0, (2.9)
(even if that is somewhat contradictory with working out Gaussian distributions in a bounded domain of the velocity variable) together with current preservation:

2K k=1+K ω k v k [f n j-1 (v k ) -fj-1 2 (-v k )] -[ fj-1 2 (v k ) -f n j (-v k )] , ⇔ 2K k=1 ω k |v k | f in (v k ) -f out (v k ) = 0. (2.10) Lemma 2.2. Let V ′ := V -vMAX 2
and Id be the identity, the modified S-matrix,

SE = S E + 1 2v M AX 1 ⊗ Ω Ω + 12 v 2 M AX V ′ V ′ ⊗ ΩV ′ ΩV ′ (Id -S E ),
enforces both mass and current preservations: (2.9), (2.10) hold for

f out = SE f in .
Proof. The equalities (2.9), (2.10) mean that f in -f out has 2 vanishing moments against 1 and |v|. By linearity, such a property extends to any test function α + β|v|, α, β ∈ R 2 , too. Hence f in -f out is orthogonal for the L 2 (-v M AX , v M AX ) scalar product to the plane spanned by {1, |v|}. Now, it is easy to see that,

1 √ 2v M AX , √ 6(|v| -v M AX ) v 3 2

M AX

, constitute an orthonormal basis of this plane. The least-squares solution,

f out + 1 2v M AX < f in -f out , 1 > L 2 + 12v ′ v 2 M AX < f in -f out , v ′ > L 2 , for v ′ := |v| -vMAX 2
has both the desired vanishing moments and is in Range( SE ). Remark 2. The S-matrix formalism, suggested in [START_REF] Glimm | The interaction of nonlinear hyperbolic waves[END_REF][START_REF] Glimm | An S-matrix theory for classical nonlinear physics[END_REF], is a manner to compute wave interactions in a very linearly degenerate framework: indeed, the Godunov scheme (2.2) inside a bounded domain can be seen as approximating a "lifted system",

∂ t f + v∂ x f - 1 τ (x) M θ (v) R f (t, x, v ′ )dv ′ -f ∂ x a(x) = 0, ∂ t a(x) = 0,
where a(x) is a staircase function constant on C j 's, and jumping of ∆x at each interface x j-1 2 (this jump yields a Dirac mass). When setting up the discrete-ordinates (2.1) in the v variable, we get a strictly hyperbolic (2K + 1) × (2K + 1) system where 2K characteristic fields are linear, and the last one is only linearly degenerate: the S-matrix gives the jump relations when the linear waves interact with the LD field.

Combining acceleration and collisions with Redheffer product.

L L 1 f in (|v|) f in (-|v|) f out (-|v|) f out (|v|) L ⋆ L 1
∂ t f + v∂ x f + E n j-1 2 ∂ v f = ∆x τ M θ (v) R f (t, x, v ′ )dv ′ -f δ(x -x j-1 2 ), with E n j-1 2 ∈ R, τ ∈ R + , by computing S τ,E n j-1 2 = S τ ⋆S E n j-1 2
the "Redheffer product" [START_REF] Redheffer | On the relation of transmission-line theory to scattering and transfer[END_REF] of the aforementioned matrices S τ and S E n j-1 2

. Roughly speaking, it all amounts to a "space-splitting" algorithm, in which any two different scattering mechanisms are flanked infinitesimally close to one another. Following [START_REF] Timotin | Redheffer products and characteristic functions[END_REF], we recall a definition of this "⋆" product of 2 partitioned 2K × 2K matrices L, L 1 (see Fig. 3.1):

L = A B C D , L 1 = A 1 B 1 C 1 D 1 ,
where both A, B, C, D and A 1 , B 1 , C 1 , D 1 are sets of K × K matrices [START_REF] Lu | Inverses of 2 × 2 block matrices[END_REF]. It reads:

L ⋆ L 1 = A(Id -B 1 C) -1 A 1 B + A(Id -B 1 C) -1 B 1 D C 1 + D 1 C(Id -B 1 C) -1 A 1 D 1 (Id -CB 1 ) -1 D , (3.1)
for Id standing for the K × K identity matrix. Obviously, for numerical efficiency, one must never implement the inverse of these K × K matrices, but instead perform simpler inversions of the resulting 2 linear systems. The Godunov scheme we propose for solving (1.1) is thus a slight modification of (2.2) which reads precisely:

f n+1 j (V) f n+1 j-1 (-V) = 1 -V ∆t ∆x f n j (V) f n j-1 (-V) + V ∆t ∆x S n j-1 2 f n j-1 (V) f n j (-V) , (3.2)
with the scattering matrix, partially deduced from a standard approximation of (1.2),

S n j-1 2 := S τ (x j-1 2 ) ⋆ S E n j-1 2 , E n j-1 2 = - ϕ n j -ϕ n j-1 ∆x . (3.3)
We stress that, since the kinetic density f n j (±V) is defined at each x j , its moments are too, and consequently also the electric potential ϕ thanks to a standard secondorder finite-difference algorithm. Hence the electric field, a first-order derivative of ϕ, is naturally computed at each interface x j- 1 2 , and this is exactly what we want within a "scattering perspective". Of course, the choice (3.3) is not unique: by analogy with second-order Strang's splitting in time, one may want to define instead

S n j-1 2 := S E n j-1 2 /2 ⋆ S τ (x j-1 2 ) ⋆ S E n j-1 2 /2 , (3.4)
thanks to the associativity of the ⋆-product, [START_REF] Timotin | Redheffer products and characteristic functions[END_REF] (for a higher CPU cost, though). Such a symmetric S-matrix reveals itself useful for instance in computing accurately a n + nn + device with zero bias imposed at both its edges. Since the scattering matrix S E doesn't automatically generate outgoing states for which the macroscopic current is a constant, the scheme (3.2)-(3.3) may not stabilize correctly in large times: possible remedies are either the procedure presented in §2.3 or an even simpler fix below.

3.2.

A simple fix to enforce current conservation. We explain now how one can proceed in order to recover a constant current across any interface x j- 1 2 in a cheap and robust manner. The Godunov scheme reads for K = 1 + K, ..., 2K:

f n+1 j (v k ) = f n j (v k ) -∆t ∆x v k f n j (v k ) -f n j-1 2 (v k ) , f n+1 j-1 (-v k ) = f n j-1 (-v k ) + ∆t ∆x v k f n j-1 2 (-v k ) -f n j-1 (-v k ) ,
where, as shown in Fig. 2.

1, f n j-1 2 (±V) = S n j-1 2 (f n j-1 (-V), f n j (-V))
. Let δ ∈ R stand for the possible macroscopic current defect, that is,

δ := 0<v k ∈V ω k v k f j-1 (v k ) + f j (-v k ) -fj-1 2 (v k ) + fj-1 2 (-v k ) . (3.5)
One must distribute this defect on both sides, fj-1 2 (±V) in order to ensure that the macroscopic flow is constant across each interface: let

J + = v k ∈V ω k v k f j (v k ) and J -= v k ∈V ω k v k f j-1 (-v k ), 0 ≤ α ≤ 1 then by modifying outgoing states like, f * j-1 2 (v k ) = fj-1 2 (v k ) + δα J + f j (v k ), f * j-1 2 (-v k ) = fj-1 2 (-v k ) - δ(1 -α) J - f j-1 (-v k ),
one recovers two outgoing states f * j-1 2 (±v k ) such that the macroscopic current is constant across the discontinuity induced by the Dirac mass (which is essential for treating correctly the mass continuity equation). In practice, bypassing this slight correction of the scattered states of the kinetic equation prevents the numerical scheme to stabilize in large times onto an approximate kinetic density f n j (±V) endowed with a flat macroscopic current (the flux term in the mass continuity equation).

We are left with a free parameter α, which value is fixed by mass conservation:

0 = 2K k=1+K ω k f j-1 (v k ) + f j (-v k ) -f * j-1 2 (v k ) + f * j-1 2 (-v k ) .
This constraint yields the following value: if |δ| > 0,

α = 2K k=1+K ω k f j-1 (v k ) -fj-1 2 (v k ) + f j (-v k ) -fj-1 2 (-v k ) δ + f j-1 (-v k ) J - 2K k=1+K ω k f j (v k ) J + + f j-1 (-v k ) J - .
However, in the cases where, by construction, S n j-1 2 preserves the total mass,

2K k=1+K ω k f j-1 (v k ) -fj-1 2 (v k ) + f j (-v k ) -fj-1 2 (-v k ) = 0,
then the expression simplifies and finally one finds:

α = 2K k=1+K ω k f j-1 (-v k ) J - 2K k=1+K ω k f j (v k ) J + + f j-1 (-v k ) J - ∈ [0, 1]. (3.6)
Notice that at numerical steady-state, velocity distributions match each other:

f * j-1 2 (|v k |) = f n j (|v k |), f * j-1 2 (-|v k |) = f n j-1 (-|v k |).
In practice, the fix (3.6) can be set up even when the local scattering matrix S n j-1 2 doesn't exactly preserve the total mass, which is a rather stringent condition because one never has an infinite support in the v variable thus "tails" |v k | ≃ v M AX are usually (slightly) flawed. Results of the forthcoming sections were obtained this way.

3.3. Stability of the resulting Godunov scheme. An interesting property of this method based on scattering S-matrices is its weak CFL criterion.

Proposition 1. Let 0 ≤ f (t = 0, x, v) ∈ L ∞ (R 2
) and assume both that,

• the mild CFL restriction ∆t ≤ ∆x/v M AX holds,

• the S-matrix S n j-1 2
is stochastic for all j, n ∈ Z × N.

Then the Godunov scheme (3.2) is positivity-preserving and obeys to:

∀j, n ∈ Z × N, 0 ≤ f n j (v k ) ≤ f (t = 0, •, •) L ∞ . (3.7)
Proof. We assume that initial data are sampled in such a manner that the L ∞ norm is preserved at time n = 0. From Lemma 2.1, any matrix S n j-1 2 preserves the L ∞ norm and non-negativity: thus for any j, k ∈ Z × {1, ..., 2K}, f n+1 j (v k ) ≥ 0 as a sum of non-negative terms. Now, by using the presentation (2.2), it comes

f n+1 j (V) f n+1 j-1 (-V) ∞ = max k |f n+1 j (v k )| |f n+1 j-1 (-v k )| = max k (1 -v k ∆t ∆x )|f n j (v k )| + v k ∆t ∆x | f n j-1 2 (v k )| (1 -v k ∆t ∆x )|f n j-1 (-v k )| + v k ∆t ∆x | f n j-1 2 (-v k )| ≤ (1 -vk ∆t ∆x ) f n j (V) f n j-1 (-V) ∞ + vk ∆t ∆x S n j-1 2 ∞ f n j-1 (V) f n j (-V) ∞ ≤ (1 -vk ∆t ∆x ) sup j ′ ,k |f n j ′ (v k )| + vk ∆t ∆x S n j-1 2 ∞ sup j ′ ,k |f n j ′ (v k )| ≤ sup j ′ ,k |f n j ′ (v k )|.
The index k is the one for which the "max" is reached. This is enough to ensure that for all j ∈ Z, there holds

sup 1≤k≤2K |f n+1 j (v k )| ≤ sup j ′ ∈Z, 1≤k≤2K |f n j ′ (v k )|,
and (3.7) follows by induction. For instance, Proposition 1 directly applies when setting up the matrix S E derived in §2.2 in the context of Hamiltonian-Preserving approximations of linear Vlasov equations (which solutions are obtained by the method of characteristics, hence the maximum principle clearly holds). Such a uniform bound was obtained in [START_REF] Shi | Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials[END_REF] too, but here we show that it holds under a weaker CFL restriction. For kinetic equations endowed with linear relaxation, Proposition 1 applies to the modified scheme (2.7) in Remark 1: the obstacle is to prove rigorously that M M -1 is stochastic for K > 1.

Remark 3. The bound (3.7) yields compactness in the weak-⋆ topology. In the special case where the S-matrix is stochastic and depends only of time (not of space),

∀j ∈ Z, S n j-1 2 ≡ S n ,
then a L ∞ dissipation estimate holds for divided differences approximating ∂ p f ∂x p (t, •, •), p ∈ N because the linear schemes (3.2), or (2.7), become invariant by x-translations.

4. Drifted-Maxwellian treatment of source and drain areas.

4.1. Computation of an approximate scattering matrix. For completeness, we hereafter recall quickly the algorithm presented in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], pages 222-224. One argues that, since the hypothesis f (t, x, v) ≃ M θ (v)ρ(t, x) is used in deriving (1.1), see [START_REF] Jüngel | Transport Equations for Semiconductors[END_REF] page 84, the same procedure may allow to simplify the acceleration term in strongly doped areas: by approaching it with E(t, x)∂ v (ρ(t, x)M θ (v)), it comes,

∂ t f + v∂ x f = 1 τ (x) M θ (v)(1 + v τ E θ ) R f (t, x, v ′ )dv ′ -f , (4.1)
supplemented by (1.2) [START_REF] Paveri-Fontana | A Neutral Gas Model for Electron Swarms[END_REF][START_REF] Frosali | Conditions for runaway phenomena in the kinetic theory of swarms[END_REF]. The terminology "drifted Maxwellian" comes from:

η(v|E) := 1 + v τ E θ M θ (v) ≃ exp - v 2 -2vτ E 2θ ≃ M θ (v -τ E),
for τ |E| ≪ 1. It corresponds to a modeling where the motion of electrons is governed by the host medium's convection and scattering [START_REF] Paveri-Fontana | A Neutral Gas Model for Electron Swarms[END_REF]. If τ |E| is small and the density of electrons has a velocity profile close to M θ (v -τ E), a displaced Maxwellian determined by the temperature θ and the local drift velocity of the background medium, both models can be expected to yield similar results (and not far from the ones of the corresponding diffusion approximation) because even if the two kinetic equations describe two distinct physical processes, they share the same diffusion approximation. At a given interface x j-1 2 , let's assume we know an average value ρ of the macroscopic density: the method of elementary solutions suggests an ansatz in the form,

f (x, ±v k ) ≃ M θ (±v k ) α + β( x τ ∓ v k ) + E( x τ , ±v k , ν) + ρ Ex θ , v k ∈ V, where E = E n j-1 2
is the electric field, ρ = ρ(f n j-1 , f n j ) is a local approximation of the macroscopic density, for instance, one may detect any spurious negative value with

ρ(f n j-1 , f n j ) = R f n j-1 (v)f n j (v)dv,
E stands again for the damped modes (2.5). For a convenient choice of the coefficients, the expression is meant to match the incoming states:

M θ (v k )[α -βv k + E(0, v k , ν)] = f n j-1 (v k ), M θ (-v k ) α + β( ∆x τ + v k ) + E( ∆x τ , -v k , ν) + E∆x θ ρ = f n j (-v k ),
and with the same coefficients, the outgoing states fj-1 2 (±V) read: fj-

1 2 (v k ) = M θ (v k ) α + β( ∆x τ -v k ) + E( ∆x τ , v k , ν) + E∆x θ ρ , fj-1 2 (-v k ) = M θ (-v k ) [α + βv k + E(0, -v k , ν)] .
Just recalling the matrix S τ from (2.6), it boils down to:

fj-1 2 (V) -ρM θ (V) E∆x θ fj-1 2 (-V) = S τ f n j (V) f n j+1 (-V) -ρM θ (V) E∆x θ . (4.2)
The Godunov scheme (2.2)-(4.2) supplemented by a standard finite-differences Poisson solver for (1.2) is the algorithm which is used in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF] to obtain the numerical results of Fig. 11.5 (see page 229). Its large-time behavior is consistent with [START_REF] Ben Abdallah | Relative entropies for kinetic equations in bounded domains (irreversibility, stationary solutions, uniqueness)[END_REF] but one drawback is that the electrons temperature remains too close to the lattice's θ.

Implementation of inflow Maxwellian boundary conditions.

Realistic computations are carried out in a computational domain which is bounded in the space variable. In order to implement Maxwellian injection on the left/right boundaries within a well-balanced framework, one discretizes the left state (handling the right one is done similarly) as

f lef t (V) = ρ D (x lef t )M θ (V) ∈ R K + .
Without loss of generality, assume that x lef t = 0, C 0 is the "ghost cell" centered in j = 0 so the interface with the first computational cell

C 1 is x 1 2 = ∆x 2 .
To advance in time the Godunov scheme (3.2), one needs "scattered boundary states" f 1 2 (V), solutions of

  f 1 2 (V) - E n 1 2 ∆x θ ρn 1 2 M θ (V) f 1 2 (-V)   = S τ f lef t (V) f n 1 (-V) - E n 1 2 ∆x θ ρn 1 2 M θ (-V) ,
where

E n 1 2 = 1 ∆x (ϕ lef t -ϕ n 1 )
. The K values f 1 2 (-V) are useless. These boundary conditions are implemented in heavily doped regions (the source and the drain for a one-dimensional n + nn + device) where the deviation from the normalized Maxwellian can be fairly considered small. Hence temperature jumps at borders remain small too. The uniform L ∞ bound in Proposition 1 cannot hold here because there is a net mass flow through both the boundaries until the macroscopic current becomes a constant in the whole computational domain, at steady-state [START_REF] Ben Abdallah | Relative entropies for kinetic equations in bounded domains (irreversibility, stationary solutions, uniqueness)[END_REF].

5. Equilibration of a n + nn + one-dimensional transistor. We display hereafter the behavior of these "scattering schemes" on elementary models of n + nn + devices: the temperature is fixed to θ = 1 2 . The computational domain is x ∈ (-1, 1) (the depleted channel lies in the interval x ∈ (-1 2 , 1 2 )) with inflow Maxwellian boundary conditions; it is coarsely gridded with 2 6 points so ∆x = 2 -5 and the CFL meets with Proposition 1: ∆t = 0.95∆x/v M AX . A double Legendre Gaussian quadrature rule is applied for griding symmetrically the velocity variable with v M AX = 2.8. The numerical process we set up for these computations reads:

• in both highly doped zones (source and drain), the scheme (2.2)-(4.2),

• in the low-doped channel, the scheme (3.2)-(3.4) with the current fix (3.6). We aim at checking the positivity-preserving properties and large-time stabilization arising from the standard choice of initial data induced by the doping profile:

f (t = 0, x, v) = ρ D (x)M θ (v),
x ∈ (-1, 1).

The Figs. 5.1. Well-balancing without collisions in the channel. We first want to display the behavior of the most elementary model where the linear Vlasov equation is set up in the device's channel, which is equivalent to let the collision rate diverge τ → +∞. Vlasov equation describes ballistic charge transport, so such a model can be the one of a "metal-semiconductor" device, a Schottky junction. Hence, in the interval x ∈ (-1 2 , 1 2 ), we simply use the scheme (3.2) with either the scattering matrix SE (see Lemma 2.2), or the Hamiltonian-Preserving one S E with the current fix (3.6). Results at time t = 100 are displayed in Fig. 5.1. One can see that the whole picture is completely symmetric around x = 0, the macroscopic current is of the order to 10 -8 , constant in the doped zones (source and drain) and monotone decreasing in the channel. The temperature is around 0.51 in the channel, 2% higher than the lattice's value θ. There is no collision process to tame the electric acceleration for x ∈ (-1 2 , 1 2 ) in the presence of the other discontinuous parameters (see

Steady-state balance with zero bias.

Here we use a constant (non-zero) collision rate in the channel (see Table 5.2) and keep on working out numerically the case where a null bias is applied at the edges of the device. We must deal with another type of difficulty because, even if we have now a collision process in the channel to tame the Vlasov acceleration term, the simple Redheffer product (3.3) is by construction asymmetrical and leads, in large time, to the onset of a spurious current in the whole device. Here, one aims at keeping as much as possible the symmetry of the initial data thus it is necessary to use the "Strang-type" star-product (3.4). Results at numerical steady-state, involving the current fix (3.6), are displayed on Fig. 5.2, where a small (spurious) current of the order of 10 -9 appears again. The fact that these spurious currents are identical with and without collisions in the channel shows that, as soon as the star-product (3.4) is used, no worsening of large-time behavior appears despite handling a higher model complexity. The parameters τ, λ 2 , ρ D (x) are the same as the ones written in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Table 11.2; numerical results can be compared to what is displayed there, on Fig. 15.4, page 303. All considered quantities are very similar to the "high field scheme", except for the temperature, which appears much lower with the present algorithm. We don't have any rigorous explanation for this difference, especially because 3 schemes (well-balanced N -scheme [START_REF] Roe | Optimum positive linear schemes for advection in two and three dimensions[END_REF], high-field WB of Chapter 15 [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], and the present one) yield 3 different temperatures.

We now define the so-called "Kullback-Leibler information",

K(t, x) := R f (t, x, v) log f (t, x, v) ρ(t, x)M θ (v) dv, x ∈ (-1, 1),
which measures the loss of information when the distribution f is replaced by the simpler one ρ(t, x)M θ (v). The evolution in time of K(t, •) L 1 (-1,1) corresponds to the red curve in the graphic, bottom right of Fig. 5.2, which is shown together with the L 2 residual on f (black curve), the

H 1 semi-norm of J (green curve, ∂ x J(t, •) L 2 (-1,1)
) and the mass growth ( ρ(t, •) L 1 (-1,1) blue curve). Residuals on f decrease monotonically, and total mass increases too. As expected, ∂ x J(t, •) L 2 (-1,1) decreases but slower than residuals (a plateau seems to exist around 10 -2 ). The L 1 norm of K(t, •) first grows as carriers begin to fill the channel, then decreases until it reaches a stable plateau: this indicates that close to each junction, the kinetic density f is quite different from a Maxwellian distribution, as signaled by the (small) temperature spikes. This feature was overlooked by a scheme proposed in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]: see Figs. 11.5 and 11.7.

5.3.

Flat currents with moderate bias. Now we aim at mimicking the test displayed in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Fig. 15.8 on page 308. A moderate bias V = -1 2 is imposed on the right side x = 1, and the device is expected to stabilize in large time onto a stationary regime endowed with a constant macroscopic current (which allows to stabilize the increase of the total mass contained in the device ρ(t, •) L 1 ). Results at t = 100 are displayed on Fig. 5.4, andagain, what shows up first is the electrons temperature which is clearly lower than the ones obtained by means of both the "Nscheme" and the "high field scheme" which were set up in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF]. Besides, the (constant)

direct reverse n + n + n depletion exiting current ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ Fig. 5.3.
Schematic view of a transistor as 2 p|n diodes: one directly polarized, one reversed.

value of the current is very close to the most flat one in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], and apart from the 2 small discontinuities at the junctions which appear with the "high field scheme", the macroscopic densities appear identical (and the electric fields too).

According to [START_REF] Jerome | Analysis of Charge Transport: Mathematical Theory of Semiconductor Devices[END_REF] (page 18) or [START_REF] Markowich | Semiconductor equations[END_REF], the length of a depletion layer around a p|n junction scales with λ, the square root of the scaled Debye length. The parameters given in Table 5.2 indicate that the size of the depletion layer located around x ≃ ∓ 1 2 , should be larger in the channel than in the doped areas (source or drain). This is well rendered as one can see looking at the electric field on Fig. 5.4: both spikes signal the electrostatic barriers arising from the equilibration of p|n junctions, and they are supported in the depletion layers. For instance, the depletion layer of the directly polarized junction (see Fig. 5.3) around x = -1 2 starts at x ≃ -0.9 and ends at x ≃ 0.2: we have a depletion length of 0.4 in the source, and 0.7 in the channel. The one around x = 1 2 runs from x ≃ 0.3 to x = 1. The resistance of a depletion layer grows with its length, hence we have a big current generated in the source, which carriers cannot recombine correctly in the channel. This is exactly what we want unless the transistor reduces just to an elementary p|n diode with a lot of "exiting current". In order to recover a big current in the drain, one reverse-polarizes the second p|n junction and sets up a strong potential drop (the forward bias -|V |): this reversed diode acts just as an anode which captures all the electrons, not recombined in the channel. A depletion layer can be thought of as an intrinsic semiconductor, hence ballistic transport of carriers is possible (thanks to Bloch theory) up to lattice vibrations which increase quickly with the temperature θ. Hence it would make sense to let the relaxation time τ , at least in the channel, grow with θ, too.

It is remarkable that the elementary kinetic model (1.1), (1.2) which doesn't contain any recombination/depletion information, succeeds in reproducing such a complex behavior, just relying on the self-consistent Poisson equation (1.2). The dynamics of this model (henceforth of the numerical scheme) are governed by:

• particles enter the computational domain through the ohmic contacts (1.3),

• they cannot settle in the doped areas because λ is low (quasi-neutrality),

• they accelerate and slowly fill in the channel (where λ is bigger) because they feel the anode's attraction (rendered by means of the forward bias) until a stable balance is reached with Poisson's self-consistent repulsive forces. Hence the balancing effect of the scattering matrices is of critical importance as it is the only effect still present in the kinetic model which is able to mimic the depletion layers. The net gain is that it suffices to run a simple model instead of simulating a more involved two-species equation endowed with a recombination process.

6. Velocity saturation with field-dependent relaxation time. It is stipulated, e.g. in [START_REF] Jerome | Analysis of Charge Transport: Mathematical Theory of Semiconductor Devices[END_REF] (see page [START_REF] Dimarco | Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors[END_REF], that in order to render a saturation effect, namely that the drift velocity and the collision rate must be bounded functions of the modulus of the electric field |E| in order to prevent infinite values if |E| becomes unbounded in the vicinity of transition points (junctions). This modification, leading to the well-known phenomenon of velocity saturation has 2 antagonistic outcomes:

• it tames even more the electric acceleration near the boundaries of the channel by strongly increasing the collision rate. Consequently, the possible "temperature spikes" are affected there too, • it worsens both the nonlinearity and the stiffness of the kinetic equation, thus numerically achieving a stationary balance between all these nonlinear terms becomes even more delicate. Corresponding parameters are shown in Table 6.1 where the (discontinuous) Debye length was also decreased in order to decrease the possible mass inflow in the device before reaching steady-state. This choice increases the electric field around the junctions, thus makes the effect of the field-dependent collision rate more noticeable. 6.1. Zero bias with a regime close to quasi-neutrality. With respect to the values displayed in Figs. 5.1 and 5.2, one sees at once on Fig. 6.1 that now, the amplitude of the electric field's modulus has nearly been multiplied by 4. The greater difficulty in achieving a global balancing among all these nonlinear terms (acceleration and collisions) appears as the resulting numerical stationary current ceases to be monotone: however, the saturation process reduces the remaining spurious currents to an order of 10 -9 , nearly identical than the ones appearing in the (less nonlinear) former cases. The steady-state macroscopic density is noticeably lower too, as an effect of the shorter Debye lengths (the electrons screening distance is bigger, hence selfconsistent repulsive effects are stronger). Without any forward bias, the temperature in the channel is now 3% over the lattice's θ. 6.2. Applying a rather strong forward bias. Since we are closer to a regime of quasi-neutrality, it is necessary to apply stronger biases in order to create macroscopic currents in the device. Here we select V = -3 2 , a value triple with respect to what was applied for generating the results of Fig. 5.4. The numerical stationary macroscopic current presents a deviation with respect to its arithmetical mean of the order of 10 -8 , but it is strongly non-monotonic. Its value is around 0.0113, which is nearly the triple of the one showing up on Fig. 5.4, ≃ 0.004 (in spite of a higher value of ρ). This comes from the triple bias yielding local stiffness in the electric field. Temperatures in the channel are higher, too, but significantly lower compared to the ones computed in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], Chapter 15. Clearly, being a moment of order 2, the temperature is a quantity sensitive to truncation errors (especially for high-velocity particles) as its computation also asks for both the macroscopic density and momentum,

T (t, x) = R v - J(t, x) ρ(t, x) 2 f (t, x, v)dv, J(t, x) = R vf (t, x, v)dv := ρu(t, x),
with u(t, x) the macroscopic velocity. Our WB schemes based on scattering S-matrices produce moments of order 0 and 1 which are quite insensitive to ∆x, the grid size, but order 2 moments may still be considered as "fragile". Remark 4. A field-dependent carriers mobility is meant to render the "pinchoff " of a depleted channel submitted to an excessive potential drop. Numerically, it is possible to see the saturating (stationary) currents as a function of the applied forward bias: see Fig. 6.3. The straight line corresponds to the linear regression of currents generated with biases lower than 4, for which the device works purely as a resistor.

Conclusion.

Finite volumes method is built to enforce mass conservation property (in opposition to finite elements). Well-balanced extensions aim at furnishing moreover qualitatively correct stationary mass flows too. Observe that an oscillating current arising close to the junctions indicates a poor resolution of the depletion layer, and is likely to produce a discrepancy on the resulting temperature as well.

Besides, it may be possible to perform an integration of the present methods, where numerical fluxes are based on S-matrices, inside the nowadays popular Discontinuous Galerkin algorithms [START_REF] Cheng | Positivity-preserving discontinuous Galerkin schemes for linear Vlasov-Boltzmann transport equations[END_REF] in order to raise up the overall accuracy. More involved kinetic models, where the inflow boundary value problem for the stationary equations cannot be easily solved, may seem difficult to handle by means of our present approach: in such cases, Bellman's Invariant imbedding theory [START_REF] Bellman | An introduction to invariant imbedding[END_REF], where scattering coefficients are derived by means of Ricatti differential equations, may help.

ω 2 P k 2 R max ′ (v) v-z dv + i τ k R max(v) v-z dv ≃ 1 -i τ kv th v th z -1 k 2 λ 2 0 v th z 2 , v th z → 0
where λ 0 = v th ωP . It corresponds to a damped oscillator which exhibits 2 different regimes depending on the location of τ is located with respect to the critical value:

2ω P τ c = 1 ⇔ τ c = 1 2 √ ρ 0 .
1. For τ < τ c , the regime is strongly collisional and exponential decay of the perturbations occurs, quickly damping the kinetic density f . The asymptotic decay rate, as ω P τ → 0, is given by the dielectric frequency ω D . 2. For τ > τ c , the regime is weakly collisional and damped oscillations are observed on the average value (i.e. the Fourier component in ξ = 0) of the electric field. The complex frequency satisfying D(ω, k; τ ) ≃ 0 reads:

ω = ω r + iω i = ±ω P 1 - 1 (2ω P τ ) 2 - i 2τ
.

Thus one can check numerically the half-period of the oscillations and the speed of decay, which should read T τ = π ωr and exp(-t/2τ ) respectively. This is what is displayed on Fig. A.1 for various values of the relaxation parameter τ : the numerical oscillations are superimposed with the theoretical half-periods and decay rates and a good match appears.

Fig. 2 . 1 . 1 .

 211 Fig. 2.1. Riemann problem in the discrete ordinates and scattering matrix Sτ Remark 1. By defining another density g(t, x, v) := f (t,x,v) M θ (v) , (2.2)-(2.6) becomes

Fig. 3 . 1 .

 31 Fig. 3.1. Illustration of Redheffer star-product

Fig. 5 . 2 .

 52 Fig. 5.2. Steady-state at zero bias (t = 100) with constant collision rate in the channel.

Fig. 5 . 4 .

 54 Fig. 5.4. Steady-state at bias V = -1 2 (t = 100) with constant collision rate in the channel.

Fig. 6 . 1 .

 61 Fig. 6.1. Steady-state at zero bias (t = 100) with field-dependent collision rate.

Fig. 6 . 3 .

 63 Fig. 6.3. Current-Voltage relation showing linear resistor and saturating behavior.

  0

Fig. A. 1 .

 1 Fig. A.1. Numerical decay w.r.t time of the average electric field for τ = 0.7, 0.3, 1 10 .

Table 5 . 1

 51 5.1, 5.2, 5.4, 6.1, 6.2 display (left to right, top to bottom) the macroscopic density ρ, current J, temperature T , electric field E, kinetic density f in the x, v plane, and either various residuals decay (if zero bias is applied), or the macroscopic velocity u at numerical steady-state (t ≃ 100). Doping, relaxation time, scaled Debye length and lattice temperature for the "Schottky diode".

	Parameter Source Channel Drain
	ρ D (x)	1	0.02	1
	τ (x)	0.01	+∞	0.01
	λ 2 (x)	0.15	0.5	0.15
	θ	1 2	1 2	1 2

Table 5 .

 5 [START_REF] Aamodt | Useful identities for half-space problems in linear transport theory[END_REF] hence this test-case is delicate to stabilize, despite its apparent simplicity.

	Parameter Source Channel Drain
	ρ D (x)	1	0.02	1
	τ (x)	0.01	1	0.01
	λ 2 (x)	0.05	0.5	0.05

Table 5 . 2

 52 Doping values, relaxation time, scaled Debye length for a standard diode.

Table 6 . 1

 61 Doping values, relaxation time, scaled Debye length for the saturation diode.

	Parameter Source Channel Drain
	ρ D (x)	1	0.02	1
	τ (x, E) λ 2 (x)	0.01 0.015	1 1+2|E| 0.05	0.01 0.015

Acknowledgement. The author thanks Dr. Ch. Buet for having pointed out the relevance of stochastic matrices in the context of numerical "elementary solutions".

Appendix A. Weakly collisional Landau damping. For the sake of completeness, we set up hereafter a benchmark already considered in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws[END_REF], pages 225-227 and inspired by the theoretical calculations published in [START_REF] Degond | Particle Simulations of the Semiconductor Boltzmann Equation for One Dimensional Inhomogeneous Structures[END_REF]. The computational domain is x ∈ (-1, 1) with periodic boundary conditions; it is gridded with 2 7 points thus ∆x = 2 -6 . The initial data reads

When set up with the ballistic Vlasov-Poisson equation, those initial data lead to the well-known phenomenon of Landau damping, that is to say a stable exponential decay of the electric field endowed with periodic oscillations. It turns out that this behavior


Mathieu Hébert

Yule-Nielsen effect in halftone prints: Graphical analysis method and improvement of the Yule-Nielsen transform

Keywords: Color reproduction, spectral model, halftone color, printing, reflectance, transmittance, Yule-Nielsen model

The well-known Yule-Nielsen modified spectral Neugebauer model is one of the most accurate predictive models for the spectral reflectance of printed halftone colors which expresses the spectral reflectance of halftones raised to the power 1/n as a linear combination of the spectral reflectance of the fulltone colors (Neugebauer primaries) also raised to the power 1/n, where n is a tunable real number. The power 1/n transform, characteristic of the Yule-Nielsen transform, empirically models the nonlinear relationship between the spectral reflectances of halftones and fulltones due to the internal propagation of light by scattering into the printing support, a phenomenon known as "optical dot gain" or "Yule-Nielsen effect". In this paper, we propose a graphical method permitting to observe this non-linear relationship in the case of single-ink halftones and to experimentally check the capacity of the Yule-Nielsen model to predict it accurately. In the case where the Yule-Nielsen transform is not well adapted to the considered type of prints, we propose alternative transforms in order to improve the prediction accuracy.

INTRODUCTION

Digital printing technologies have considerably increased the number of applications and printable materials in a context where designers and even customers themselves are more and more used to select the printing materials and parameters in order to precisely obtain the expected appearance. Precise color management needs a calibration of the printing system for the selected printing support, inks and halftoning technique [START_REF] Hersch | Interaction between light, paper and color halftones: challenges and modelization approaches[END_REF]. The classical calibration method is based on a color profile established by measuring the color of about one thousand of halftones [START_REF] Kipphan | Handbook of Print Media[END_REF], a procedure that must be repeated each time the printing support, the ink set or the halftoning technique is modified. The use of a color prediction model can considerably ease the color calibration since only a few tens of colors need to be printed and measured. One of the most accurate models available today is the Yule-Nielsen modified spectral Neugebauer model [START_REF] Viggiano | The Color of Halftone Tints[END_REF], or simply "spectral Yule-Nielsen model", which is widely used for the calibration of printing systems [START_REF] Balasubramanian | Optimization of the spectral Neugebauer model for printer characterization[END_REF][START_REF] Iino | Building color management modules using linear optimization I. Desktop[END_REF][START_REF] Iino | Building color management modules using linear optimization II. Prepress system for offset printing[END_REF][START_REF] Xia | End to end color printer calibration by total least squares regression[END_REF] and can be improved by the method by Crété & al. for estimating the effective surface coverage of the inks [START_REF] Hersch | Improving the Yule-Nielsen modified spectral Neugebauer model by dot surface coverages depending on the ink superposition conditions[END_REF].

The spectral Yule-Nielsen model predicts the spectral reflectance of halftone raised to the power 1/n as a linear combination of the spectral reflectances of the fulltone colors (i.e. colors made by covering the whole surface with one or several inks, also called "colorants" or "Neugebauer primaries") also raised to the power 1/n, where n is a tunable parameter. The weights attributed to the different fulltone colors are the effective surface coverages of these colors in the halftone. The power 1/n transform empirically models the nonlinear relationship between the spectral reflectances of halftones and fulltones due to the internal propagation of light by scattering into the printing support, a well-known phenomenon also called "optical dot gain" or "Yule-Nielsen effect" [START_REF] Yule | The penetration of light into paper and its effect on halftone reproduction[END_REF]. In this paper, we propose a graphical method showing this non-linear relationship in an intelligible way for single-ink halftones. The graph represents the reflectance of a single-ink halftone as a function of the reflectance of this fulltoned ink. If the spectral measurement contains k wavebands, then k points can be displayed in the graph and it is easy to check how well the Yule-Nielsen function matches them. The measured points generally follow a curved line which fairly coincides with the plot of the Yule-Nielsen transform computed with the appropriate n value and effective ink surface coverage. In some cases however, the coincidence is not perfect and alternative transforms, based on other functions than the power 1/n function, may be found. We propose a few ones which provide modest but appreciable gain in prediction accuracy, this latter being measured in terms of average CIELAB ΔE1994 values over several tens of printed colors. The graphical method is also very interesting to display effects that are rather difficult to detect by looking directly at the spectral reflectances.

SPECTRAL REFLECTANCE AND TRANSMITTANCE OF HALFTONES COLORS

Before analyzing the relevance of the spectral Yule-Nielsen model for a given type of print, we propose to recall in this section the whole prediction method, including the ink spreading assessment proposed by Hersch and Crété [START_REF] Hersch | Improving the Yule-Nielsen modified spectral Neugebauer model by dot surface coverages depending on the ink superposition conditions[END_REF]. The model is presented for reflectance of CMY halftones but it similarly applies for transmittance and for any ink colors.

Halftone colors

Halftoning is a color reproduction technique based on a finite set of coloring substances, the inks, which are deposited according to binary grids of small patterns (inks dots). The halftone surface is a mosaic of colored areas corresponding to the partial overlap of the ink dots. The areas with no ink, those with a single ink layer, and those with two or three superposed ink layers are each one considered as distinct colorants (also called Neugebauer primaries). In CMY printing, based on cyan, magenta and yellow inks, one obtains eight colorants corresponding to the areas with no ink, cyan alone, magenta alone, yellow alone, red (magenta and yellow), green (cyan and yellow), blue (cyan and magenta) and black (cyan, magenta and yellow). In classical clustered-dot or error diffusion prints, the fractional area occupied by each colorant can be deduced from the surface coverages of the primary inks according to Demichel's equations [12]. Denoting as c, m, and y the surface coverages of the cyan, magenta and yellow inks, the surface coverages of the eight colorants are:

( )( )( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )
+ + + ++ = - = -- - = -- = - = -- = - = -- = (1)

The Yule-Nielsen model

The first color prediction model for halftone prints, proposed by Neugebauer [START_REF] Neugebauer | Die theoretischen Grundlagen des Mehrfarbendrucks[END_REF], assumes that the color of the halftone is a linear combination of the colors of the colorants contained in it. The spectral version of this model can be written:

( ) ( ) 8 1 λλ ii i R aR = = ∑ (2)
where R(λ) denotes the spectral reflectance of the printed halftone, ( )

λ i R
the spectral reflectances of the colorants i printed each one on a large area to be measured, and a i the respective surface coverages of the colorants in the halftone.

However, this linear equation (2) does not predict correctly the color of the halftone due to the scattering of light within the paper bulk which induces a color halo around the ink dots. This phenomenon, known as the "Yule-Nielsen effect", is even amplified by the multiple reflections between the paper bulk and the print-air interface [START_REF] Hersch | Interaction between light, paper and color halftones: challenges and modelization approaches[END_REF]. In order to account for this effect, Vigiano [START_REF] Viggiano | The Color of Halftone Tints[END_REF] proposed to apply the Yule-Nielsen transform [START_REF] Yule | The penetration of light into paper and its effect on halftone reproduction[END_REF] on the spectral Neugebauer equation:

( ) ( ) 8 1/ 1 λλ n n ii i R aR =  =  ∑ ( 3
)
where n is real number, generally higher than 1 (when n = 1, the spectral Yule-Nielsen equation ( 3) is equivalent to the Neugebauer equation), to be fitted. It is generally related to the strength of scattering in the printed support [START_REF] Arney | A Probability Description of the Yule-Nielsen Effect I[END_REF], on the halftone screen frequency and on the ink-paper interaction [START_REF] Hébert | Analyzing halftone dot blurring by extended spectral prediction models[END_REF].

Ink spreading assessment

At printing time, the inks may spread on the surface in different amounts according to the presence or not of other inks on the surface [START_REF] Yule | The penetration of light into paper and its effect on halftone reproduction[END_REF]. As a consequence, the effective surface coverages of the inks, thereby those of the colorants in the halftone, are not known. A solution to obtain them is to establish experimentally the correspondence between nominal and effective surface coverages through so-called ink spreading functions. Let us present the method to compute them.

( ) ( ) ( ) ( ) 1/ 1/ / /& λ1 λ λ n nn uv v uv u v R a R aR  ′′ = -+  (4)
The effective surface coverage / uv a′ is fitted by minimizing the sum of square differences between the measured spectral reflectance of the halftone, M(λ), and the spectral reflectance predicted by Eq. ( 4):

( ) ( ) ( ) ( ) { } 2 730nm 1/ 1/ /& 01 λ 380nm arg minλ1 λ λ n nn uv v u v x a M x R xR ≤≤ =  ′ = -- +  ∑ (5)
Repeating this procedure for the 36 halftones yields 12 lists of 3 effective surface coverages / uv a′ associated with the nominal surface coverages {0.25, 0.5, 0.75}. It is assumed that the effective surface coverage is 0, respectively 1, when the nominal surface coverage is 0 (no ink), respectively 1 (full coverage). By linear interpolation, one obtains continuous curves ( )

/ uv af a ′ =
, where a and a′ denote the nominal, respectively the effective surface coverages of the halftoned ink (see examples in Figure 3).

With the fitted surface coverages, the spectral reflectances predicted by Eq. (4) are as close as possible to the measured spectral reflectances but not equal to them. In order to assess the difference between predicted and measured spectra, we use a visual metric, the CIELAB ΔE 94 color distance, calculated by converting the two spectra first into CIE1931-XYZ tristimulus values by selecting the D65 illuminant, and then into CIELAB color coordinates using as white reference the spectral reflectance of the unprinted paper illuminated with the D65 illuminant [START_REF] Wyszecki | Color Science: Concepts and Methods, Quantitative Data and Formulae[END_REF].

The optimal n value can be determined by trying successively different values, e.g. n = 1 to 10 in steps of 0.1, and keep the n value for which the average ΔE 94 value computed over the 36 halftones is minimal. Figure 4 shows two plots of the average ΔE 94 value as a function of the n value for CMY colors printed in inkjet with the Canon PixmaPro 9500 Mark II printer (a) at 150 lpi on Canon PP201 260 g/m 2 glossy coated paper (spectral measurements based on the 45°:0° geometry using the X-rite i1 spectrophotometer) and (b) at 100 lpi on common 80 g/m 2 office paper (spectral measurements based on the diffuse:8° geometry using the X-rite Color i7 spectrophotometer). A minimum is clearly displayed in the case of the coated paper at n = 3.5, whereas the curve asymptotically decreases in the case of the office paper, which indicated that the optimal n tends to infinity (in practice, we can set n to 20 or 100).

Prediction

Once the ink spreading functions have been obtained, the spectral reflectance of any halftone can be predicted. Let us consider a halftone where the cyan, magenta and yellow ink are printed at the nominal surface coverages c 0 , m 0 , and y 0 . These surface coverages are converted into effective ink surface coverages c, m and y by a weighted average of the ink spreading functions in order to account for the superposition-dependent ink spreading [START_REF] Hersch | Improving the Yule-Nielsen modified spectral Neugebauer model by dot surface coverages depending on the ink superposition conditions[END_REF]. The weights are expressed by the surface coverages of the respective colorants on which the ink halftone is superposed. For example, the weight of the ink spreading function c f (cyan halftone over white colorant) is ( )( )

11 my --. In the case of three halftoned inks, effective surface coverages are obtained by performing a few iterations with the following equations:

( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) /0 / 0 /0 / 0 / 0 /0 / 0 / 0 / 0 /0 / 0 / 0 1 1 1 1 1 1 1 1 1 1 1 1 cw cm c y cm y mw mc m y mc y yw yc ym yc m c m y f c m y f c m yf c myf c m c y f m c y f m c yf m cyf m y c m f y c m f y c mf y cmf y + + + =- - + - +- + =- - + - +- + =- - + - +- + (6)
For the first iteration, 0 = cc , 0 = mm and 0 = yy are taken as initial values on the right side of the equations. The obtained values of c, m and y are then inserted again into the right side of the equations, which gives new values of c, m, y and so on, until the values of c, m, y stabilize. They are then plugged into Eq. ( 1) in order to obtain the effective surface coverages of the eight colorants, which are finally used in Eq. ( 3) to predict the spectral reflectance of the halftone.

Spectral transmittance of halftone colors

It has been recently shown that the spectral Yule-Nielsen transform applies to the spectral transmittance of halftone colors, either printed on paper [START_REF] Hébert | Yule-Nielsen based recto-verso color halftone transmittance prediction model[END_REF] or on transparency films [START_REF] Machizaud | Spectral reflectance and transmittance prediction model for stacked transparency and paper both printed with halftone colors[END_REF], in a similar way as for reflectances:

( ) ( )

8 1/ 1 λλ =  =  ∑ n n kk k T aT ( 7
)
where T(λ) denotes the spectral transmittance of the halftone, a i the surface coverage in the halftone of each colorant i whose spectral reflectances ( )

λ i T
is measured using a spectrophotometer in transmittance mode, e.g. the X-Rite Color i7 instrument in total transmittance mode (diffuse:0° geometry).

Note that for a same halftone sample, the reflectance and transmittance equations may involve different n values because scattering phenomena are different in reflection and transmission. The ink spreading functions may also be different and it is recommended to compute them from reflectance measurements for reflectance predictions, and from transmittance measurements for transmittance predictions.

GRAPHICAL ANALYSIS OF THE YULE-NIELSEN TRANSFORM

Although the Yule-Nielsen equation in known to be empirical and its physical justification has not been completely clarified despite several attempts [START_REF] Ruckdeschel | Yule-Nielsen effect in printing: a physical analysis[END_REF][START_REF] Arney | The Physics Behind the Yule-Nielsen Equation[END_REF][START_REF] Lewandowski | Applying the Yule-Nielsen equation with negative n[END_REF], one generally observes good agreement between predicted and measured spectra for a wide range of printing techniques. However, it is rather difficult to interpret precisely the deviations that are sometimes observed between them. We thus propose a graphical method permitting the comparison in a more intelligible way in the case of single-ink halftones.

Yule-Nielsen function and Fulltone-to-halftone reflectance diagram

Let us denote as 0 R , a R and 1 R the spectral reflectances of patches where the ink has the respective surface coverages of 0, a and 1. The Yule-Nielsen equation for a single-ink halftone, already given in Eq. ( 4), is written:

( ) ( ) ( ) ( ) 1/ 1/ 01 λ1 λ λ n nn a R a R aR  ′′ = -+  (8)
where a′ denotes the effective surface coverage of the ink corresponding to the nominal surface coverage a. Using the notations 0 / a yRR = and 10 / xRR = , Eq. ( 8) can also be written, for each wavelength:

1/ 1 n n y a ax  ′′ = -+  (9)
We will call "Yule-Nielsen function" the function 1/

(1 )

nn x a ax ′′ →-+
. By plotting it for given values of n and a, one obtains a curve showing graphically the relationship between fulltone and halftone reflectances normalized by the unprinted support reflectance (Figure 5). When the ink is nonabsorbent at a given wavelength, we have This relationship can be verified experimentally from the measurements of ( )

0 λ R , ( )
λ a R and
( )

1 λ R
, which provides as many points as values contained in the measured spectra. We obtain a "fulltone-to-halftone diagram" as the ones showed in Figures 6, 7 and8. The spectral reflectances needed to plot these diagrams are systematically measured for the calibration of the ink spreading functions.

Eq. ( 9) can also be written

1/ 1/ 1 nn y a ax ′′ = -+ , ( 10
)
which means that the Yule-Nielsen model assumes a linear relationship between . This can also be experimentally verified in a "transformed halftone-to-fulltone diagram" (see Figures 6 and7).

Physical meaning of the Yule-Nielsen function

The ratio 10 / xRR = can be interpreted as the effective transmittance of the ink layer for the average light path of light in the fulltone print. In a halftone, we may consider the layer as composed of absorbing areas with effective transmittance x and non-absorbing areas with effective transmittance 1. In absence of scattering, the effective transmittance of the halftone layer would be 1 y a ax ′′ = -+ , which is the Neugebauer function. But since the layer is generally scattering, we can consider that after transmission through a sublayer of relative thickness 1/n, therefore of transmittance 1/ 1, n t a ax =-+ light is mixed before entering the next sublayer. Since the light crosses n sublayers of relative thickness 1/n, the transmittance of the whole layer is n t . This yields the Yule-Nielsen function defined by Eq. [START_REF] Yule | The penetration of light into paper and its effect on halftone reproduction[END_REF].

Although this interpretation is not realistic if one imagines light trajectories in the printing material, it has at least a physical meaning in term of statistical light transfers. A light scattering model would be necessary to establish the equivalence between this simple model and a more rigorous description of the propagation of light in the scattering halftone. But the experience shows that this equivalence must exist when the Yule-Nielsen model is capable of accurate prediction of colored halftones. Note that this interpretation of the Yule-Nielsen transform has been recently used to define a hybrid additive-subtractive color mixing system [START_REF] Simonot | Between additive and subtractive color mixings: intermediate mixing models[END_REF].

Illustration of the graphical method with halftones printed in inkjet

Figures 6 illustrates the graphical method through the example of a halftone printed with either cyan, magenta or yellow ink at nominal surface coverage 0.5 by using a periodical clustered dot halftoning at 150 lpi, and printed on Canon PP201 260 g/m 2 glossy paper (same paper as for the experiment presented in Figure 4) with the Canon Pixma Pro9500 inkjet printer. The spectral reflectances of the unprinted paper, ( )

0 λ R , of the halftones, ( ) λ a R
, and of the fulltones, ( )

1 λ R ,
were measured from 380 nm to 730 nm in steps of 10 nm (36 wavebands) with the X-rite i1 spectrophotometer based on the 45°:0° geometry. These spectra are plotted in the top row of the figure. In dashed line are plotted the spectral reflectances of the halftones predicted by the spectral Yule-Nielsen model with the optimal n = 3.5 fitted in the calibration step (see. Section 2.3) and the effective surface coverage a′ obtained from equation ( 5) and specified on top of each column. The fulltone-to-halftone diagrams, in the middle row, show the 36 points ( 10RR; 0 a RR ) computed from the measured spectral reflectances, as well as the Yule-Nielsen function given by Eq. ( 9) with n = 3.5 (solid line) and the Neugebauer function given by Eq. ( 9) with n = 1 (dashed line) for the specified effective surface coverage a′ . At the bottom of the figure are shown the transformed fulltone-to-halftone diagrams.

We see from the fulltone-to-halftone diagrams that the Yule-Nielsen function matches well the experimental points, which indicates that the Yule-Nielsen transform correctly models the nonlinear relationship between the reflectances of the halftone, the fulltone, and the unprinted paper. This is confirmed by the low ΔE 94 values computed from the measured and predicted spectra: 0.20 units for the cyan halftone, 0.55 units for the magenta halftone and 0.09 for the yellow halftone. This is also confirmed by the fact that the points ( ) In order to verify the prediction accuracy of the model over a larger set of colors, we printed 125 halftones with the same printer, same inks, same paper and same halftoning technique as the colors used for Figure 6 and measured them using the same spectrophotometer. These 125 colors correspond to all the combinations of cyan, magenta and yellow inks printed at the nominal surface coverages 0, 0.25, 0.5, 0.75 and 1. The average ΔE 94 value over these 125 colors is 0.51 units, and the 95%-quantile of the 125 ΔE 94 values is 1.12 units. We can therefore conclude that the Yule-Nielsen model is well performing for this type of print.

;

Limitations of the Yule-Nielsen function

Let us now study halftones printed with the same printer, the same inks and the same halftoning technique on a nonfluorescent, non-calendered and noticeably porous 80 g/m 2 paper called Biotop. Figure 7 presents similar graphs as those presented in Figure 6, corresponding to the cyan ink halftones at 0.25, 0.5 and 0.75 nominal surface coverages. The Yule-Nielsen functions where calculated with an optimal n value of 100. This example shows that good agreement between prediction and measurement cannot be achieved for all types of prints: Despite the optimized n and a' values, the Yule-Nielsen functions do not match the experimental points satisfyingly, and the points ( )

; xy ′′ , with 1/ 10 () n x RR ′ = and 1/ 0 () n a y RR ′ =
, are not aligned. Consequently, the prediction accuracy of the spectral Yule-Nielsen model is poorer: the ΔE 94 values for the halftones printed at 0.25 and 0.5 surface coverage are beyond the value of the just noticeable color difference threshold (ΔE 94 = 1). On a set of 64 patches designed for verification (all combinations of the cyan, magenta and yellow inks printed at the nominal surface coverages 0.13, 0.38, 0.63 and 0.88), we obtained an average ΔE 94 value of 1.23 units, and a 95%-quantile of 2.36 units. Figure 8 shows three other examples of cyan ink halftone printed at nominal surface coverage 0.5. On the left and the center of the figure, the studied sample is the same as in the left column of Figure 6 but instead of measuring its spectral reflectance with the X-rite i1 spectrophotometer based on the 45°:0° geometry, we used the X-rite Color i7 spectrophotometer based on the diffuse:8° geometry in either specular reflection included mode (left column) or specular reflection excluded mode (central column). The difference that we observe between the two measurement modes is due to the fact that the cyan ink displays a reddish sheen in the specular direction which is more pronounced in the fulltone than in halftones. This sheen is responsible for the deviation between measured and predicted spectra beyond 570 nm in the specular reflection included mode. In the fulltone-to-halftone diagram, the left-most experimental points in the graph follow "u"-shaped line and are thus located on either side of the Yule-Nielsen function, which means that y have different values for a same value of x. No function can reproduce such an incurved line since y can have only one value for one x value. However, the global prediction accuracy over the 125 halftone colors described in the previous section is still acceptable because the colored sheen is visible in a limited number of halftones, i.e. those where the cyan ink has a high surface coverage and is not superposed with other inks. The average ΔE 94 value obtained with the optimal n = 5.5 was 0.91 units, therefore below the just noticeable color difference, and the 95%-quantile 1.95 units. The prediction accuracy is nevertheless poorer than in the specular reflection excluded mode where the reddish sheen is not observed: for the same set of 125 halftone colors, the average ΔE 94 value obtained with the optimal n = 5.6 was 0.62 units and the 95%-quantile 1.19 units, therefore comparable with the prediction accuracy achieved with measurements based on the 45°:0° geometry.

The graphs on the right of Figure 8 show the case of the same cyan ink printed at 120 lpi on Canson 90 g/m 2 tracing paper. Since the tracing paper is semi-transparent, the reflectance is much lower in comparison to the other papers (see the ordinate scale of the plotted reflectances). For this sample, the experimental points are dispersed on either side of the Yule-Nielsen function. The prediction accuracy achieved on a set of 40 CMY halftones is very poor: we obtained an average ΔE 94 value of 2.99 units and a 95%-quantile of 4.58 units.

ATTEMPTS TO IMPROVE THE YULE-NIELSEN MODEL

Rewriting the Yule-Nielsen equation (3) as ( ) ( )

8 1 1 λλ ii i R f af R - =  =    ∑ , ( 11
)
with the function 1/

:

n fx x →
, and the reverse function

1 1/ : n f yy - →
can be a source of inspiration to find other functions f in order to improve the spectral Yule-Nielsen model when it fails to provide accurate prediction accuracy. The proposition by Lewandowski et al. [START_REF] Lewandowski | Applying the Yule-Nielsen equation with negative n[END_REF] to select a negative n value in the Yule-Nielsen model, while n is rather supposed higher than 1, can be considered as a first example of extended Yule-Nielsen function.

The new functions f i should be defined and continuous on [0, 1] and satisfy the condition ( )

11 i f = .
The reflectance a R of a single-ink halftone of surface coverage a′ is written:

( ) ( ) ( ) ( ) 1 11 1 i a ii R f afR a fR - ′′ = -+ , ( 12
)
where 0 R and 1 R denote the reflectances of the printing support without, respectively with ink. It can also be written

( ) ( ) ( ) ( ) 1 00 1 ia i ii fR fR aa fR fR ′′ = -+ (13)
but not ( ) ( )

0 10 1 ia i fRR a a fRR ′′ = -+ (14)
except for the functions of the form

α : fx x →
where α is a real number (i.e. those on which is based the classical Yule-Nielsen model). The fulltone-to-halftone diagram, based on the values [START_REF] Arney | A Probability Description of the Yule-Nielsen Effect I[END_REF] RR and 0 a RR , is therefore not helpful to optimize the new function, but we can check graphically in the transformed fulltone-to-halftone diagram whether, for the selected function, the points

( ) ( ) ( ) ( ) ( ) 01 0
;

ia i i i
fR fR fR fR are aligned.

The searching of alternative functions is empirical and may be tedious. A first idea to generate new functions is to make simple transforms on the classical function itself, by introducing for example a shift: ; :

( ) 1/
n n f x xs f y y s - →- → + ( 15
)
where s is a real number that somehow represents surface reflection, not subject to optical dot gain. Note that since this function does not satisfy the condition ( )

11 f
= , the s value should be rather small.

The classical function f can also be modified in the following way:

( ) ( ) ( ) ( ) 2 1/ 21 22 41 : 1 ; : 21 n n q q qy f x qx q x f y q -  - +-   → -- → - (16)
and function f 1 and f 2 can be combined:

( )( ) ( ) ( ) ( ) 2 1/ 2 1 33 41 : 1 ; : 21 n n q q qy f x qx q x s f y s q - - +- → -- - → + - (17)
Another transform of the function 1 f is

( ) ( ) ( ) 1/ 1 44 : ln 1 / ln ; : 1 n n x f x xs b b f y b b s -  → - -+ → +- +  (18)
where "ln" denotes the natural logarithm function. More complicated functions can be found, for example

( ) ( ) ( ) ( )
where "arcsh" denotes the inverse of the hyperbolic sine function, and ( ) ω arcsh ln 2 n = .

These new functions f 1 to f 5 were tested on the 64 halftones printed on Biotop paper already presented in Section 3.2 and Figure 7. For the model calibration and the predictions, similar procedure is followed as the one detailed in Section 2 by simply replacing the classical function f with the considered function f i . The prediction accuracy yielded by each function for the 64 halftones is presented in Table 1, as well as the values of the tunable parameters which were fitted in order to obtain the lowest average ΔE 94 value. The five functions yield very similar accuracy, which is slightly better than with the classical function f although the average ΔE 94 value remains beyond the just noticeable color difference (ΔE 94 = 1). The improvement, modest but appreciable, is partly explained by the introduction of additional tunable parameters, but also by the fact that the experimental points are better aligned in the transformed fulltone-to-halftone diagram, as shown in Figure 9 through the example of the cyan halftone printed at the surface coverage 0.5 on the Biotop paper (same sample as the one studied in the central column of Figure 7).

) n f x qx q x → -n = 100; q = 0.12 1.11 2.12 ( ) iii i fR fR fR fR issued from the measured spectral reflectances of the Biotop paper with no ink (reflectance R 0 ), with full coverage ink (reflectances R 1 ) and with halftoned ink at 0.5 nominal surface coverage (reflectances R 0.5 ), for the classical function f and the functions f 1 and f 2 computed with the parameters presented in Table 1. The lines of equation:

2 1/ 3 : ( 1 ( )) n f x qx q x s → -- - n = 50; s = 0.06; q = -0.
1 y a ax ′ ′ ′′ = -+
, where a′ is the effective surface coverage of the halftone fitted for each function, are plotted in dashed line.

The three diagrams in Figure 9 have been generated by the classical function f and the functions f 1 and f 2 (those obtained with the functions f 3 to f 5 , quasi-identical to these latters, are not reproduced). The effective surface coverages a′ specified on top of the graphs were fitted for each function. In this example, functions f 1 and f 2 yield the same a′ -value, lower than the one fitted with f. This indicates that f 1 and f 2 model stronger optical dot gain than f. We also see that the experimental points are closer to the straight line of equation: 1 y a ax ′ ′ ′′ = -+ with the new functions than for the classical one. Is seems that the improvement of accuracy is due to the fact that the line passes through the points corresponding to the highest absorbance of the ink (left-most points in the diagram), which is due to the fact that the ΔE 94 value is more sensible to spectral deviations in wavelengths domains where the spectrum is low.

CONCLUSIONS

Despite its physical base is not completely established yet, the spectral Yule-Nielsen model improved by Hersch and Crété's ink spreading assessment method provides a high prediction accuracy for most types of halftone prints. However, it is less performing for certain types of printed surfaces, especially when the printing support is lowly reflective like tracing paper, or when the surface of the print has a colored sheen in the specular direction. The graphic method proposed in this paper enables verifying the relevance of the Yule-Nielsen transform for the considered printing setup, by plotting for single-ink halftones the Yule-Nielsen function ( ) It may happen that the 10 0 (; ) a RRR R points do not draw a well defined curve in the fulltone-to-halftone diagram, i.e., that at wavelength where the unprinted paper and the fulltone have respectively same reflectance the halftone has different reflectances, due for example to optical phenomena such as fluorescence, too strong absorption or color sheen at the surface. Except these optical phenomena are explicitly modeled, the spectral Yule-Nielsen model is incapable of predicting these differences of reflectance and the prediction accuracy in this case is generally poorer.

In the case where the spectral Yule-Nielsen model is judged to be not accurate enough but the 10 0 (; ) a RRR R points draw a well defined line, one can try to find another function than the 1/n x function in order to improve the prediction accuracy. A few ones have been proposed in this paper, but other ones can be invented. The improvement of prediction accuracy may be modest, as in the examples presented here, but still appreciable.

Figure 1 .

 1 Figure 1. Detail of a halftone color with cyan, magenta and yellow inks, showing the ink dots and the mosaic of 8 colorants.

Figure 2 .

 2 Figure 2. Halftones used for the ink spreading assessment. The set of 36 halftones shown in Figure 2 are printed. They correspond to the combinations of one ink u at nominal surface coverage a = 0.25, 0.5 or 0.75, and one solid colorant v obtained by printing each other ink at a surface coverage 0 or 1. Each of these halftones contains two colorants corresponding to the areas where ink u is superposed to colorant v [effective surface coverage / uv a′ , spectral reflectance

Figure 3 .

 3 Figure 3. Example of ink spreading functions giving the effective surface coverage a' of each ink superposed with each colorant made of the other two inks at the surface coverages 1 or 0, as functions of the nominal surface coverage a.

Figure 4 .

 4 Figure 4. Variation as a function of the n value of the average ΔE 94 computed from the predicted and measured spectral reflectances of the 36 single-ink halftones displayed in Figure 2, in the case of (a) inkjet CMY colors printed at 150 lpi on glossy coated paper and (b) in the case of CMY colors printed at 100 lpi on common office paper.

  . The Yule-Nielsen function thus satisfies ( )11 f= and the curve passes through the point (1, 1) in the diagram. When the ink is fully absorbing, we have 1 case where n = 1, the Yule-Nielsen equation becomes the Neugebauer equation (2) and the Yule-Nielsen function is a straight line that intercepts the ordinate axis in 1 ya ′ = -. It is tangent to the Yule-Nielsen function in x = 1, i.e. the relationship between halftone and fulltone reflectance is almost linear when the fulltone has nearly the same reflectance as the unprinted support, therefore when the ink is almost not absorbing.

Figure 5 .

 5 Figure 5. Plots of the Yule-Nielsen function defined by Eq. (9) for various n-values and three different surface coverages a.

  for the fitted values of n and a′ .

Figure 6 .

 6 Figure 6. Top row: Measured spectral reflectances of cyan (left column), magenta (central column), and yellow (right column) halftones printed at the surface coverages 0.5 (reflectances R 0.5 ) and 1 (reflectances R 1 ) on glossy coated paper (reflectance R 0 ). Plots in dashed line represent the spectral reflectances of the halftones predicted by the spectral Yule-Nielsen model for the values of a′ specified on top of each column. The CIELAB ΔE 94 values assess the deviation between predictions and measurements. Central row: Fulltone-to-halftone diagrams showing the points ( 10 RR; 0 a RR) deduced from the measured spectra in the different wavebands, as well as the Yule-Nielsen function given by Eq. (9) with n = 3.5 (solid line) and the Neugebauer function, for the specified value of a′ . Bottom row: Transformed fulltone-to-halftone diagrams showing the experimental points ( 1/ 10 () n RR ;

Figure 7 .

 7 Figure 7. Similar graphs as in Figure6for cyan ink halftones printed at 0.25, 0.5 and 0.75 nominal surface coverages on Biotop non-calendered paper. The Yule-Nielsen functions were calculated with n = 100.

Figure 8 .

 8 Figure 8. Spectral reflectances and fulltone-to-halftone diagrams for cyan ink halftones printed at 0.5 nominal surface coverage on the PP201 glossy coated paper (left and center) and on tracing paper (right). The spectral measurements are based on the diffuse:8° geometry with specular component included (left and right) and excluded (center).

Figure 9 .

 9 Figure 9. Transformed fulltone-to-halftone diagrams representing the points ( ) 1 0 0.5 0 () () ;( ) ()

  with the optimal n and a′ values computed in the calibration step and comparing it with the points 10 0 (/; /) a R RR R computed in the different wavebands from the measured spectral reflectances of the paper,

Table 1 .

 1 Performance of the spectral Yule-Nielsen model based on the classical function f and on transformed functions f i on 64 halftones covering the printer's gamut printed on Biotop non-calendered paper.

	Function				Parameters	Average ΔE 94	95-quantile ΔE 94
	: fx	→	1/ x	n			n = 100	1.23	2.36
	1 : f x	( →-xs	) 1/	n	n = 20; s = 0.09	1.12	2.1 0
	2 :		(			2 1/

ACKNOWLEDGEMENT

The author would like to thank Baptiste Paul, Serge Mazauric and Jacques Machizaud for their contribution to this study.