G Sallen

A Tribu

T Aichele

R André

C Bougerol

S Tatarenko

K Kheng

J Ph

Poizat

Photon correlation spectroscopy on a single quantum dot embedded in a nanowire

Keywords: numbers: 78, 67, Lt, 78, 55, Et

We have observed strong photoluminescence from a single CdSe quantum dot embedded in a ZnSe nanowire. Exciton, biexciton and charged exciton lines have been identified unambiguously using photon correlation spectroscopy. This technique has provided a detailed picture of the dynamics of this new system. This type of semi conducting quantum dot turns out to be a very efficient single photon source in the visible. Its particular growth technique opens new possibilities as compared to the usual self-asssembled quantum dots.

Semiconductor nanowires (NWs) appear as promising building blocks for nanoscale devices and circuits with impressive potential applications including nanoelectronics [1,2,3], optoelectronics (light emitting diodes [4,5], nanolasers [6]), thermoelectrical energy conversion [7], and biological or chemical sensors [8]. Moreover, high quality defect free nanowires can be grown on low-cost, routinely used substrates such as silicon, which means that they could easily be used for fabricating commercial devices and could possibly be integrated with mainstream Si microelectronics devices.

NW growth methods allow for the variation of the chemical composition [9,10] or doping [11] along the longitudinal or radial directions. This enables the fabrication of well controlled 1D nanoscale heterostructures [10]. For example, as shown in this work, it is possible to insert a slice of a low band gap semiconductor within a high bandgap NW and thus realize a light emitting quantum dot (QD) operating as a single photon source [12,13]. So far, work on the light emitting properties of single quantum dots has mainly concerned self-assembled QDs formed by surface forces induced by lattice mismatch between different materials. Such QDs have been widely used in the past decade as single photon sources [14] and for their potential application in quantum information processing (see for example [15,16]). QDs in NWs appear to be an interesting alternative to self-assembled quantum dots. The absence of a wetting layer offers a better confinement which could enable room temperature production of single photons [13]. Radial growth techniques enable engineering of optical guides allowing more efficient light extraction than in bulk materials [17,18]. Furthermore, NW based heterostructures, being much less limited by lattice mismatches, greatly widen the possible materials combinations and enable well controlled stacking of several QDs in a single NW, offering interesting possibilities for quantum information processing [19].

In this letter we present the first detailed optical characterization of excitonic emission in a single CdSe QD embedded in a ZnSe NW. We have already shown that this system is an efficient single photon source operating at temperature as high as 220 K [13]. Single pho- ton emission from NWs has otherwise only been demonstrated at 4 K in InAs QDs embedded in InP NWs [12]. Our system emits light around 550 nm (2.2eV) where silicon avalanche photodiodes (APDs) are very efficient. This has allowed us to perform the first thorough spectroscopic analysis of a QD embedded in a NW by using photon correlation spectroscopy [20]. We have identified unambiguously the exciton, biexciton, and charged exciton lines and obtained information on the charging dynamics of this QD.

The wires are grown by Molecular Beam Epitaxy (MBE) in the Vapour-Liquid-Solid (VLS) growth mode catalysed by gold particles on a Si substrate. In order to make QDs, a small region of CdSe is inserted in the ZnSe NW. This is done by interrupting the ZnSe growth, changing to CdSe growth for a short time and then growing ZnSe again [13]. Details about the growth of the ZnSe NWs can be found in reference [21]. The wire diameter (around 10 nm) is of the order of the bulk exciton Bohr diameter for CdSe (11 nm). This means that the carriers in the CdSe QD are in the strong quantum confinement regime. For the study of single NWs, the sample is sonicated in methanol, causing NWs to break off the substrate into the solution. Droplets of this solution are then deposited on a Si substrate, and a low density of individual NWs is obtained after evaporation. As shown in fig. 1, individual NWs can be isolated, allowing single QD optical spectroscopy.

The experimental apparatus is a standard microphotoluminescence (µPL) set-up. The samples are mounted on a XYZ piezo motor system in a He flow cryostat at a temperature of 4 K. The optical excitation is provided by a 405 nm continuous-wave (CW) diode laser illuminating the sample via a microscope objective of numerical aperture N A = 0.65 located in the cryostat. The NW emission is collected by the same objective and sent to a 50/50 beamsplitter for correlation measurements. In each arm of the beamsplitter, the light is dispersed by a monochromator (1200 grooves/mm grating, 30 cm and 50 cm focal length respectively). Each monochromator has a switchable mirror inside, which can direct the luminescence either onto a charge coupled device (CCD) camera for the measurement of the PL spectrum or through the exit slit towards a low jitter (40 ps), high quantum efficiency APD. The detectors send electrical pulses into a time-correlated single photon module that builts an histogram of the time delays between successive photons. This histogram is proportional to the second order correlation function g (2) (t) [14]. The overall temporal resolution of our set-up is essentially limited by the jitter of the APDs and the dispersion of the monochromator gratings. This time resolution was measured by recording the autocorrelation function of 1 ps pulses from a frequency-doubled Ti:Sapphire laser. A full width at half maximum of 90 ps was obtained for the autocorrelation function peak.

The results presented in this work all come from the same QD. A typical µPL spectrum is shown in the inset of fig. 2 where three lines can be seen. A comparison with relative energy positions of known emission lines in spectra of self-assembled CdSe/ZnSe QDs [22,23] suggests that these lines correspond to the exciton (X), the biexciton (XX) and the charged exciton (CX). Unambiguous proof for the assignement of these lines will be given below using photon correlation spectroscopy. The width of the lines is due to spectral diffusion [24].

One of the characteristic features of such NW QD structures is their polarization properties. As shown in [13] the excitation efficiency and the luminescence are both strongly polarization dependent [25,26,27]. A 90% contrast is obtained for the excitation efficiency depending on the direction of the linear polarization of the pumping laser. The light emission is also 90% linearly polarized in the same direction as the best pumping polarization, independently of the excitation polarization.

It can be seen in fig. 2 that the saturation level of the XX line is more than three times larger than that of the X line. This is due to a strong storage effect of the dark exciton state owing to a rather large dark-bright exciton splitting (∆E = 6 meV) [28]. The dark exciton state reduces the intensity of the X line owing to the leakage from the bright to the dark exciton states caused by spin flip but it remains an efficient intermediate state for populating the XX state [29].

The level scheme used to model our system is shown in figure 3. It is based on models used in references [30,31,32] where the carriers can enter the QD either individually (power dependent rates γ C1 and γ C2 ) or already bound as excitons (power dependent rate r). We have added the dark exciton, which plays a key role in our system. It includes the bright and dark exciton, the charged exciton, and the biexciton. A triexciton level and a charged biexciton level are also included to avoid artificial saturation of the biexciton and of the charged exciton but they are not represented in figure 3.

All of the power independent parameters of the model can be evaluated independently prior to photon correlation experiments by performing lifetime measurements [28]. Their values are listed in table I. These values are compatible with what has already been measured on selfassembled CdSe/ZnSe QDs [23].

We come now to the main results of this work on pho- ton correlation experiments. We present first the data concerning the neutral QD in fig. 4. The autocorrelation of the X line emission is shown in fig. 4 (a) exhibiting a clear antibunching which is characteristic of the statistics of a single photon emitter [14]. Fig 4 (b) shows the cross correlation measurement between the X and the XX line. It displays the typical asymmetric shape with bunching and antibunching features that is the signature for the cascaded emission of a XX photon followed by a X photon [33]. This allows us to identify unambiguously these two lines as exciton and biexciton of the same QD. Note that the narrow bunching peak can only be fitted if the dark exciton is included in the model.

For all the correlation graphs (figs. 4 and 5) the right vertical axes are the raw number of coincidences. The left axes represent the normalized correlation function according to a Poisson statistics where the coincidences involving background photons have been subtracted. The corrected correlation function g (2) is related to the uncorrected one g

(2) u by g (2) -1 = (g

u -1)/ρ 2 , where ρ = S/(S + B) with S and B respectively the number of signal and background photons as measured in the spectrum of fig. 2.

The autocorrelation of the CX line is shown in fig. 5 (a). As for the X line, it exhibits a clear antibunching. Also, on a larger time scale, a bunching effect can be observed. This is due to the hopping from the charged state to the neutral state of the QD as confirmed by the X-CX cross correlation displayed in fig. 5 (b). This proves that the CX line comes from the same QD as the X line. The situation is depicted schematically in fig. 5 (c), which shows that the photons are emitted either from the charged or from the neutral state of the QD. The average time spent by the quantum dot in the charged state is given by the characteristic time of an exponential fitting of the bunching peak in fig. 5 (a) which is 5 ns.

It should be noticed that the antibunching dip in figure 5 (b) is not symmetrical. Negative (positive) time corresponds to the probability of detecting a CX (X) photon after having detected a X (CX) photon. Formation of a CX in an empty QD (t < 0) requires the loading of three charges, whereas the formation of an X in a QD with a single charge (level E c in fig. 3) is faster since it requires only the loading of a single charge (t > 0) [30,31,32]. Note that we do not know whether this charge is positive or negative. We can only say that by comparison with spectra for self-assembled CdSe/ZnSe QDs [22,23] the charge state is probably negative.

As can be seen in figs. 4 and 5 the experimental results are fitted very well by the model shown in fig. 3 taking into account the temporal resolution of our experimental setup (90 ps). Inclusion of the dark exciton turned out to be essential for the modelling of the photon correlation data. Allowing for the coexistence of two excitation mechanisms, namely charge by charge (described by γ C1 and γ C2 ) or directly feeding the QD with an already bound exciton (described by r) has also turned out to be necessary for the fitting. The QD charge hopping time depends on the value of these parameters. Note that XX autocorrelation and the CX-XX cross-correlation (not shown) are also fitted well using the same set of parameters [24]. The model gives also the correct intensities of the spectral lines within 10% [24].

In conclusion, we have used photon correlation spectroscopy to characterize a new type of light emitting quantum dot embedded in a nanowire. We obtained a very good fit to the experimental data with a model based on a standard excitonic level scheme. This allowed us to extract quite complete dynamics of the neutral and charged excitons including charge hopping between these two states of the QD. CdSe/ZnSe QDs in NWs are nanoobjects situated between self-assembled QDs and CdSe based colloidal nanocrystals [34,35]. The latter operate at room temperature but have a blinking problem and a lifetime above 20 ns. On the other hand, self-assembled QDs are non-blinking and feature a sub-nanosecond lifetime allowing GHz repetition rates. QDs embedded in NWs have the potential to combine the best of both worlds by offering non-blinking room temperature [13] single photon sources with a high repetition rate. Furthermore, the versatility of this particular nanostructure growth technique offers interesting perspectives for engineering semi-conducting QDs, such as coupled QDs or waveguide coupled QDs.

Figure 1 :

 1 Figure 1: Scanning electron microscope image of a single CdSe/ZnSe nanowire deposited on a silicon substrate.

Figure 2 :

 2 Figure 2: Inset : µPL spectrum at an excitation power of P = 15µW . The electronic background noise has been subtracted. Main plot : Line intensities as a function of excitation power.

Figure 3 :

 3 Figure 3: Level scheme including the empty dot (E0), the dark exciton (EX D ), the bright exciton (EX B ), the biexciton (EXX), the charged dot (EC), and the charged exciton (ECX ). The various rates between the different level are indicated in the figure.

Figure 4 :

 4 Figure 4: (a) Exciton emission autocorrelation, (b) excitonbiexciton cross-correlation for an excitation power P = 15µW . The left axes are the correlation functions corrected from the background and the right axes are the raw coincidence rates (see text). The fit is performed using the model based on fig. 3. The power dependent parameters used for the fit are r = 0.6ns -1 , γC1 = γC2 = 1ns -1 . The other parameters are given in table I.

Figure 5 :

 5 Figure5: (a) Charged exciton autocorrelation for an excitation power P = 8µW , (b) Exciton-Charged exciton crosscorrelation for an excitation power P = 10µW . The fits are performed using the model based on fig.3. The power dependent parameters used for the fits are respectively r = 0.31ns -1 , γC1 = 0.09ns -1 , γC2 = 0.058ns -1 for the CX autocorrelation (a) and r = 0.37ns -1 , γC1 = 0.25ns -1 , γC2 = 0.28ns -1 for the X-CX cross-correlation (b). The other parameters are given in table I. (c) Representation of the stream of photons coming alternatively from the neutral and the charged QD.

Table I :

 I Table of power independent transition rates in ns -1 .

	γX γXX γCX γSP 1 γSP 2 γNR
	1.4 2.5 1.7 1.4	0	0.2

We thank L. Besombes for many fruitful discussions, F. Donatini for very efficient technical support, and R. Cox for a careful reading of the manuscript. T.A. acknowledges support by Deutscher Akademischer Austauschdienst (DAAD). Part of this work was supported by European project QAP (Contract No. 15848).


Claire Joseph
email: claire.joseph@univ-antilles.fr

Maryse Moutamal
email: maryse.moutamal@aims-cameroon.org

Quasi-boundary method for a fractional ill-posed problem

Keywords: Mathematics Subject Classification. 35K57, 35R25, 26A33 Riemann-Liouville derivative, quasi-boundary method, fractional diffusion equation, illposed problem, inverse-problem

A quasi-boundary method is used to study an ill-posed, time-fractional diffusion equation involving the fractional Riemann-Liouville derivative. In particular, we consider an ill-posed problem for a family of well-posed problems, and prove, by means of eigenfunction expansions, that the solutions of the latter problems converge to the solutions associated with the former problem. The analysis presented includes providing conditions for the rate of the convergence.

Introduction

Let d ∈ N * and Ω be a bounded open subset of R d , for a boundary ∂Ω of class C 2 . For T > 0, we set Q = Ω × (0, T ), Σ = ∂Ω × (0, T ) and consider the fractional diffusion equation:

   D α
RL y(x, t) -∆y(x, t) = f (x, t) in Q, y(σ, t) = 0 on Σ, I 1-α y(x, T ) = y 1 (x)

on Ω, (1.1) FVP where 0 < α < 1, f ∈ L 2 (Q), y 1 ∈ H 1 0 (Ω). The operators I 1-α y and D α RL y are, respectively, the Riemann-Liouville fractional integral of order 1 -α, and the left Riemann-Liouville fractional derivative of order α of y.

The fractional diffusion equation has been of significant interest for many decades. The equation (specifically the time-fractional diffusion equation) is obtained by replacing the first order time derivative of the classical diffusion equation with a time fractional derivative. In comparison with the conventional first order derivative, Left Riemann-Liouville fractional derivatives are characterised by a convolution integral (see Definition 2.6). This shows that the (fractional) derivative depends on the behavior of function y over the interval [0, t]. This is the reason why researchers speak about the memory effect associated with a fractional derivative, and why, in this context fractional derivatives are used in other fields such as Physics, Biology or Economics, where the memory association of a field is mandatory.

The model compounded in Equation (1.1), can be used to investigate environmental phenomenon. In such cases, one might not have all the necessary information on complete the model. In the case considered here, the initial condition is missing so that Equation (1.1) appears as an ill-posed backward time-fractional diffusion equation. The problem compounded in Equation (1.1) does not satisfies the Hadamard conditions. This is because we cannot prove that (1.1) admits a unique solution which depends continuously of y 1 . Nowadays, there are many methods such as inverse methods that may be used to approach such ill-posed problems. In this paper, we use a Quasi-boundary method which was originally introduced by the quasi-reversibility method developed in [START_REF] Lattès | The Method of Quasireversibility : Applications to Partial Differential Equations[END_REF].

The quasi-boundary method is based on perturbing the final condition. Some researches on the topic have shown that this method gives better numerical results than the quasi-reversibility method. For example, in [START_REF] Yang | A quasi-boundary regularization method for identifying the initial value of time-fractional diffusion equation on spherically symmetric domain[END_REF], Yang et al. apply the quasi-boundary method to approximate an inverse problem for identifying the initial data for a time-fractional diffusion equation on a pherically symmetric domain. Jayakumar [START_REF] Jayakumar | Modified quasi-boundary value method for the multidimensional nonhomogeneous backward time fractional diffusion equation[END_REF] use a modified quasi-boundary method to solve a non-homogeneous time fractional diffusion problem involving the left fractional Caputo derivative. More recently, Huynh et al. [START_REF] Huynh | Recovering the space source term for the fractional-diffusion equation with Caputo-Fabrizio derivative[END_REF] applied a modified quasi-boundary method for a fractional diffusion equation involving the Caputo-Fabrizio fractional derivative.

We refer to [START_REF] Karapinar | Identifying the space source term problem for time-space-fractional diffusion equation[END_REF][START_REF] Liu | A backward problem for the time-fractional diffusion equation[END_REF][START_REF] Shi | A fractional-order quasi-reversibility method to a backward problem for the time fractional diffusion equation[END_REF][START_REF] Wang | Quasi-reversibility method to identify a space-dependent source for the time-fractional diffusion equation[END_REF][START_REF] Denche | A modified quasi-boundary value method for ill-posed problems[END_REF][START_REF] Metzler | Boundary value problems for fractional diffusion equations[END_REF][START_REF] Showalter | The final value problem for evolution equations[END_REF][START_REF] Yang | The quasi-reversibility regularization method for identifying the unknown source for time fractional diffusion equation[END_REF] and references therein for more information in regard to the quasi-reversibility method and quasi-boundary method. In this context, the Riemann-Liouville and Caputo based fractional derivatives are closely related. The best of the authors' knowledge, and, judging from the open literature available, there are no studies on the quasi-boundary method for fractional diffusion equations involving the Riemann-Liouville fractional derivative. In this paper, we approach the ill-posed problem compounded in Equation (1.1) through a family of well-posed problems. More precisely, we consider, for any β > 0, the following quasi-boundary value problem:

   D α RL y β (x, t) -∆y β (x, t) = f (x, t) (x, t) ∈ Q, y β (σ, t) = 0 (σ, t) ∈ Σ, I 1-α y β (x, T ) + βI 1-α y β (x, 0 + ) = y 1 (x) x ∈ Ω,
and prove that the family of solutions y β converge to the solution of Equation (1.1) in an appropriate Hilbert space, specifying the rate of convergence. This paper is structured as follows: In Section 2, we provide some definitions on fractional operators, examples of their properties and some preliminary results. In Section 3, we use the spectral method to prove the existence and uniqueness of the solution of the problem. The convergence results are provided in Section 4. def2 Definition 2.2 [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] Let x and y be two complex numbers such that Re(x) > 0 and Re(y) > 0. The Beta function is given by

B(x, y) = 1 0 t x-1 (1 -t) y-1 dt.
Remark 2.2 One can prove [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] that

B(x, y) = Γ(x)Γ(y) Γ(x + y) .
(2.1) beta_formule Definition 2.3 [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] Let α > 0 and β > 0. Then, the two-parameter Mittag-Leffler function is given by

E α,β (z) = ∞ k=0 z k Γ(αk + β) , z ∈ C. (2.2) mittag
Thus, we have

E α,α (z) = ∞ k=0 z k Γ(αk + α) , z ∈ C. (2.3) Eaa
In what follows, we set

E α,1 (t) = E α (t).
majora Theorem 2.1 [START_REF] Podlubny | Fractional Differential Equations[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF] Let 0 < α < 2 and β ∈ R. We consider that µ satisfies πα 2 < µ < min{π, πα}.

In this case, there exists a constant C = C(α, β, µ) > 0, such that

|E α,β (z)| ≤ C 1 + |z| , µ ≤ |arg(z)| ≤ π.
Definition 2.4 [START_REF] Haubold | Mittag-Leffler Functions and Their Applications[END_REF][START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF] Let α, β, ρ ∈ C such that Re(α) > 0 and Re(β) > 0. Then, the generalized Mittag-Leffler function is defined by

E ρ α,β (t) = +∞ n=0 (ρ) n t n Γ(αn + β)n! , for all t ∈ C,
where (ρ) n = ρ(ρ + 1) . . . (ρ + n -1).

Remark 2.3 Note that, when ρ = 1 we obtain

E 1 α,β (t) = E α,β (t), where E α,β is the classical Mittag-Leffler function defined in (2.2).
The following result gives the Laplace transform of the generalized Mittag-Leffler function.

Lemma 2.1 [START_REF] Haubold | Mittag-Leffler Functions and Their Applications[END_REF] Let α, β, ρ be complexes such that Re(α) > 0, Re(ρ) > 0 and Re(β) > 0. Then, we have

L -1 s ρ-1 s α + as β + b ; t = t α-ρ +∞ k=0 (-a) k t (α-β)k E k+1 α,α+(α-β)k-ρ+1 (-bt α ), (2.4) Lpmittag
where |as β /(s α + b)| < 1 and L -1 is the inverse Laplace transform.

Definition 2.5 [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] Let f : R + → R be a continuous function, and α > 0. Then, the expression

IRL I α f (t) = 1 Γ(α) t 0 (t -s) α-1 f (s)ds, t > 0,
is called the Riemann-Liouville integral of order α of the function f . def26 Definition 2.6 [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] Let f : R + → R. The left Riemann-Liouville fractional derivative of order α ∈ (0, 1) of f is defined as

D α RL f (t) = 1 Γ(1 -α) . d dt t 0 (t -s) -α f (s)ds, t > 0,
provided that the integral exists.

Remark 2.4 From Definition 2.5, we see that :

D α RL f (t) = d dt I 1-α f (t). ( 2

.5) DandI

We also with the right Caputo fractional derivative given by: Definition 2.7 [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] Let f : R + → R, 0 < α < 1. The right Caputo fractional derivative of order α of f is defined by

D α C f (t) = -1 Γ(1 -α) T t (s -t) -α f ′ (s)ds, 0 < t < T, (2.6)
provided that the integral exists.

We consider a result obtained through integration by parts, which involves the left Riemann-Liouville fractional derivative and the right Caputo fractional derivative.

integrale1 Lemma 2.2 [11] Let 0 < α < 1, y ∈ C ∞ ( Q) and φ ∈ C ∞ ( Q).
Then, we have,

T 0 Ω (D α RL y(x, t) -∆y(x, t))φ(x, t)dxdt = Ω φ(x, T )I 1-α y(x, T )dx - Ω φ(x, 0)I 1-α y(x, 0)dx + T 0 ∂Ω y(σ, t) ∂φ ∂v (σ, t)dσdt - T 0 ∂Ω ∂y ∂v (σ, t)φ(σ, t)dσdt + Ω T 0 y(x, t)(-D α C φ(x, t) -∆φ(x, t))dxdt, (2.7)
where D α C is the right Caputo fractional of order 0 < α < 1.

integrale0T Corollary 2.1 [START_REF] Mophou | Optimal control of fractional diffusion equation[END_REF] Let D(0, T ) be the set of C ∞ functions on (0, T ) with compact support. Then for all φ ∈ D(0, T ),

T 0 D α RL y(t)φ(t)dt = - T 0 y(t)D α C φ(t)dt,
where D α C is the right fractional Caputo derivative.

article1 Theorem 2.2 [12] Let 1/2 < α < 1, y 0 ∈ H 1 0 (Ω) and f ∈ L 2 (Q). Then, the problem    D α RL y(x, t) -∆y(x, t) = f (x, t) in Q, y(σ, t) = 0 on Σ, I 1-α y(x, 0) = y 0 in Ω.
has a unique solution y ∈ L 2 ((0, T ), H 1 0 (Ω)). Moreover, I 1-α y ∈ C([0, T ], H 1 0 (Ω)), and, there exists a constant C > 0, such that the following estimations hold:

∥y∥ L 2 ((0,T );H 1 0 (Ω)) ≤ ∆ ∥y 0 ∥ H 1 0 (Ω) + ∥f ∥ L 2 (Q) , (2.8
)

∥I 1-α y∥ C([0,T ];H 1 0 (Ω)) ≤ Π ∥y 0 ∥ H 1 0 (Ω) + ∥f ∥ L 2 (Q) , (2.9)
with

∆ = max   C 2T 2α-1 2α -1 , C α 2T λ 1   , and
Π = sup   C √ 2, C 2T 1-α (1 -α)   .
On other hand, it is well-known that (-∆) is a symmetric uniform elliptic operator. Thus, it admits real eigenvalues, 0

< λ 1 ≤ λ 2 ≤ λ 3 ≤ ... with λ k → ∞ when k → ∞.
Moreover, there exists an orthonormal basis {w k } ∞ k=1 of L 2 (Ω) , where w k ∈ H 1 0 (Ω) is an eigenfunction corresponding to λ k : -∆w k = λ k w k . Further, we have,

Ω ∇φ(x) • ∇ψ(x)dx = λ k Ω φ(x)ψ(x)dx, ∀p ∈ H 1 0 (Ω). ( 2

.10) valprop1

In what follows, for all φ, ψ ∈ L 2 (Ω), we denote

(φ, ψ) L 2 (Ω) = Ω φ(x)ψ(x)dx,
as the inner product in L 2 (Ω) and ∥φ∥ L 2 (Ω) as the associated norm. We set a(φ, ψ)

= Ω ∇φ(x) • ∇ψ(x)dx, ∀φ, ψ ∈ H 1 0 (Ω). (2.11) formebi1
Then, the bilinear functional a(., .) defines an inner product on H 1 0 (Ω), and we have

∥φ∥ 2 H 1 0 (Ω) = a(φ, φ), (2.12) equivalence1
which is a norm on H 1 0 (Ω). Since

w k √ λ k ∞ k=1
is an orthonormal basis of H 1 0 (Ω) for the inner product a(., .), we can write

||ϕ|| 2 H 1 0 (Ω) = +∞ i=1 λ i (ϕ, w i ) 2 L 2 (Ω) , ∀ϕ ∈ H 1 0 (Ω).
(2.13) H101

3 Approximate problem approx In this section, using eigenfunctions expansions of the Laplace operator, we prove the existence and uniqueness of solution to the approximate problem given by

   D α RL y β (x, t) -∆y β (x, t) = f (x, t) (x, t) ∈ Q, y β (σ, t) = 0 (σ, t) ∈ Σ, I 1-α y β (x, T ) + βI 1-α y β (x, 0 + ) = y 1 (x) x ∈ Ω, (3.1) eqA1 where 1/2 < α < 1, f ∈ L 2 (Q), y 1 ∈ H 1 0 (Ω) and I 1-α y β (x, 0 + ) = lim t↓0 I 1-α y β (x, t).
Let us assume that (3.1) has a solution y β ∈ C ∞ ( Q). If we multiply the first equation in (3.1) by a function v ∈ H 1 0 (Ω) and integrate by parts over Ω, we obtain

Ω D α RL y β (x, t)v(x)dx + Ω ∇y β (x, t) • ∇v(x)dx = Ω f (x, t)v(x)dx. (3.2) chin Observing that (D α RL y β (t), v) = D α RL (y β (t), v
) and using (2.11), problem (3.1) becomes for all t ∈ (0, T ):

   D α RL (y β (t), v) L 2 (Ω) + a(y β (t), v) = (f (t), v) L 2 (Ω) in Ω, ∀v ∈ H 1 0 (Ω), y β (t) = 0 on ∂Ω, I 1-α y β (x, T ) + βI 1-α y β (x, 0 + ) = y 1 in Ω. (3.3) eqB1
We can then consider the following problem :

Given 1/2 < α < 1, y 1 ∈ H 1 0 (Ω) and f ∈ L 2 (Q), find eq11 y β ∈ L 2 ((0, T ), H 1 0 (Ω)), (3.4a)
I 1-α y β ∈ C([0, T ]; H 1 0 (Ω)), (3.4b) such that diff1 D α RL (y β (t), v) L 2 (Ω) + a(y β (t), v) = (f (t), v) L 2 (Ω) ∀t ∈ (0, T ), ∀v ∈ H 1 0 (Ω), (3.5a) eq21 I 1-α y β (T ) + βI 1-α y β (0 + ) = y 1 in Ω. (3.5b) eq22
In this context, the following existence and uniqueness theorems hold.

existeyb Theorem 3.1 Let 1/2 < α < 1 and a(., .) be the bilinear form defined by (2.11). Then, the approximate problem (3.4)-(3.5) has a unique solution y β ∈ L 2 ((0, T ), H 1 0 (Ω)) given by

y β (t) = +∞ i=1          y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α ) t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i . (3.6) solution11
where λ i is the eigenvalue of the operator -∆ corresponding to the eigenfunction w i . E α,α as given in (2.3),

y 1 i = (y 1 , w i ) and f i (t) = (f (t), w i ) are respectively, the i-th component of y 1 and f (t) in the orthonormal basis {w i } ∞ i=1 of L 2 (Ω). Moreover, I 1-α y β ∈ C([0, T ], H 1 0 (Ω))
and there exists a constant C > 0 such that,

∥y β ∥ L 2 ((0,T ),H 1 0 (Ω)) ≤ Π ∥y 1 ∥ H 1 0 (Ω) + ∥f ∥ L 2 (Q) , (3.7) estimation_yb and I 1-α y β C([0,T ],H 1 0 (Ω)) ≤ Θ ∥y 1 ∥ H 1 0 (Ω) + ∥f ∥ L 2 (Q) , (3.8) estimation_Iyb
where

Π = max   2C β T 2α-1 2α -1 , 2C 2 T α β 2 (1 -α)(2α -1) + 4C 2 T α α -1 2   and Θ = sup   2C β , 2C 4 T 1-α β 2 (1 -α) + 2CT 1-α 1 -α   .
Proof. If we replace v by w i in (3.5a) and use the fact that

a(y β (t), w i ) = λ i (y β (t), w i ) L 2 (Ω) = λ i y β i ,
we deduce from (3.5) that y β i is a solution of the ordinary differential equation

D α RL y β i (t) + λ i y β i (t) = f i (t), t ∈ (0, T ), I 1-α y β i (T ) + βI 1-α y β i (0 + ) = y 1 i , (3.9) edo11
where y 1 i = (y 1 , w i ). Now, using the Laplace transform, we obtain from the first equation of (3.9) that,

Dα RL y β i (s) + λ i ŷ β i (s) = fi (s), (3.10) eqlaplace1 where Dα RL y β i (s) = L(D α RL y β i (t))(s), ŷ β i (s) = L(y β i (t))(s), fi (s) = L(f i (t))(s)
and L denotes the Laplace transform operator. From (2.5), we have Dα RL y β i (s) = -I 1-α y β i (0 + ) + s α ŷ β i (s), which, combining with (3.10), gives

-I 1-α y β i (0 + ) + s α ŷ β i (s) + λ i ŷ β i (s) = fi (s).
Hence,

ŷ β i (s) = I 1-α y β i (0 + ) × 1 s α + λ i + fi (s) × 1 s α + λ i ,
and it follows from (2.4) that

y β i (t) = I 1-α y β i (0 + )t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds.
Therefore,

I 1-α y β i (t) = I 1-α (I 1-α y β i (0)t α-1 E α,α (-λ i t α )) + I 1-α t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds , = A + B where A = I 1-α (I 1-α y β i (0)t α-1 E α,α (-λ i t α )), B = I 1-α t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds .
Let us now compute A and B. We have,

A = I 1-α (I 1-α y β i (0)t α-1 E α,α (-λ i t α )) = I 1-α y β i (0) × I 1-α (t α-1 E α,α (-λ i t α )) = I 1-α y β i (0) × 1 Γ(1 -α) t 0 (t -s) -α s α-1 E α,α (-λ i s α )ds = I 1-α y β i (0) Γ(1 -α) +∞ k=0 (-λ i ) k Γ(αk + α) t 0 (t -s) -α s α-1 s αk ds = I 1-α y β i (0) Γ(1 -α) +∞ k=0 (-λ i ) k t αk-1 Γ(αk + α) 1 0 (1 -u) -α u α-1+αk tdu = I 1-α y β i (0) Γ(1 -α) +∞ k=0 (-λ i ) k t αk Γ(αk + α) B(1 -α, αk + α),
which in view of (2.1) gives

A = I 1-α y β i (0) +∞ k=0 (-λ i t α ) k Γ(αk + 1) = I 1-α y β i (0)E α (-λ i t α ). (3.11) CalA
On the other hand,

B = I 1-α t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds = 1 Γ(1 -α) t 0 (t -s) -α s 0 (s -u) α-1 E α,α (-λ i (s -u) α )f i (u)du ds = 1 Γ(1 -α) t 0 f i (u) t u (t -s) -α (s -u) α-1 E α,α (-λ i (s -u) α )ds du = 1 Γ(1 -α) t 0 f i (u) +∞ k=0 (-λ i ) k Γ(αk + α) t u (t -s) -α (s -u) α-1+αk ds du = 1 Γ(1 -α) t 0 f i (u) +∞ k=0 (-λ i ) k (t -u) αk Γ(αk + α) 1 0 (1 -z) -α z α-1+αk dz du = 1 Γ(1 -α) t 0 f i (u)   +∞ k=0 (-λ i ) k (t -u) αk Γ(αk + α) B(1 -α, αk + α)   du = t 0 f i (u)   +∞ k=0 (-λ i (t -u) α ) k Γ(αk + 1)   du = t 0 f i (u)E α (-λ i (t -u) α )du.
Thus,

B = t 0 f i (u)E α (-λ i (t -u) α )du (3.12) CalB
Finally, adding (3.11) to (3.12), we obtain

I 1-α y β i (t) = I 1-α y β i (0)E α (-λ i t α ) + t 0 f i (u)E α (-λ i (t -u) α )du. (3.13) I1-a
Hence,

I 1-α y β i (T ) = I 1-α y β i (0)E α (-λ i T α ) + T 0 f i (u)E α (-λ i (T -u) α )du.
From (3.9) 2 , we have that,

I 1-α y β i (0)E α (-λ i T α ) + T 0 f i (u)E α (-λ i (T -u) α )du + βI 1-α y β i (0) = y 1 i ,
from which, we deduce that

I 1-α y β i (0) = y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α ) . (3.14) Iyb0
Finally, we obtain

y β i (t) =      y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α )      t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds.
The rest of proof can be done in three steps.

Step 1: We give the formulation of a solution to an approximate problem associated to Equations (3.4)- (3.5). Let V m be a subspace of H 1 0 (Ω) generated by the w 1 , w 2 , . . . , w m . Consider the following approximate problem associated to Equations (3.4) -(3.5): Find y β m : t ∈ (0, T ] → y β m (t) ∈ V m , the solution of

D α RL (y β m (t), v) L 2 (Ω) + a(y β m (t), v) = (f (t), v) L 2 (Ω) , ∀v ∈ V m , (3.15) eq111 I 1-α y β m (T ) + βI 1-α y β m (0) = y 1 m = m i=1 y β 1 i w i . (3.16) eq121
As y β m (t) ∈ V m , we have

y β m (t) = m i=1 (y(t), w i ) L 2 (Ω) w i = m i=1 y β i (t)w i .
Proceeding as per the computation of y β , we show that y β m is a solution of the problem given by Equations (3.15) -(3.16) and obtain,

y β m (t) = m i=1          y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α ) t α-1 E α,α (-λ i t α )          w i + m i=1 t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i .
(3.17) solm1

Step 2: We show that the sequences (y β m ) and (I 1-α y β m ) are respectively, Cauchy sequences in L 2 ((0, T ); H 1 0 (Ω)) and C([0, T ]; H 1 0 (Ω)). Let m and p be two integers such that p > m ≥ 1. Then, from (3.17)

y β p (t) -y β m (t) = p i=m+1 y β i (t)w i . Set b i = y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α )
. Then, we have that,

a(y β p (t) -y β m (t), y β p (t) -y β m (t)) = p i=m+1 λ i [y β i (t)] 2 ≤ 2 p i=m+1 λ i t 2α-2 E 2 α,α (-λ i t α )|b i | 2 + 2 p i=m+1 λ i t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds 2 .
Hence,

∥y β p (t) -y β m (t)∥ 2 L 2 ((0,T );H 1 0 (Ω)) = T 0 a(y β p (t) -y β m (t), y β p (t) -y β m (t))dt ≤ A p + B p , with A p = 2 p i=m+1 λ i |b i | 2 T 0 t 2α-2 E 2 α,α (-λ i t α )dt, B p = 2 p i=m+1 T 0 λ i t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds 2 dt.
Note that from Theorem 2.1, we know that there exists a generic constant C > 0 such that

A p = 2 p i=m+1 λ i |b i | 2 T 0 t 2α-2 E 2 α,α (-λ i t α )dt ≤ 2C 2 p i=m+1 λ i |b i | 2 T 0 t 2α-2 dt ≤ 2C 2 p i=m+1 λ i |b i | 2 t 2α-1 2α -1 T 0 ≤ 2C 2 T 2α-1 2α -1 p i=m+1 λ i |b i | 2 .
Remark 3.1 From the latter estimation, we see that we have to take 1/2 < α < 1 to give a sense to our computation.

Using again Theorem 2.1 and noting that 1 -α ̸ = 0, we obtain

p i=m+1 λ i |b i | 2 = p i=m+1 λ i y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α ) 2 ≤ 1 β 2 p i=m+1 λ i y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du 2 ≤ 2 β 2 p i=m+1 λ i |y 1 i | 2 + 2 β 2 p i=m+1 λ i T 0 E α (-λ i (T -u) α )f i (u)du 2 ≤ 2 β 2 p i=m+1 λ i |y 1 i | 2 + 2C 2 β 2 -(T -u) 1-α 1 -α T 0 p i=m+1 T 0 |f i (u)| 2 du. Consequently, p i=m+1 λ i |b i | 2 ≤ 2 β 2 p i=m+1 λ i |y 1 i | 2 + 2C 2 T 1-α β 2 (1 -α) p i=m+1 T 0 |f i (u)| 2 du. (3.18) sum_bi
and we have that

A p ≤ 4C 2 T 2α-1 β 2 (2α -1) p i=m+1 λ i |y 1 i | 2 + 4C 4 T α β 2 (1 -α)(2α -1)   p i=m+1 T 0 |f i (u)| 2 du   . ( 3

.19) CalAP

On the other hand, using the Cauchy-Schwartz inequality,

B p = 2 p i=m+1 T 0 λ i t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds 2 dt = 2 p i=m+1 T 0 λ i t 0 (t -s) α 2 -1 4 E α,α (-λ i (t -s) α ) (t -s) α 2 -3 4 f i (s) ds 2 dt ≤ 2 p i=m+1 T 0 λ i t 0 (t -s) α-1 2 E 2 α,α (-λ i (t -s) α )ds t 0 (t -s) α-3 2 |f i (s)| 2 ds dt,
which in view of Theorem 2.1 gives

B p ≤ 2C 2 p i=m+1 T 0 t 0 (t -s) -1/2 ds t 0 (t -s) α-3 2 |f i (s)| 2 ds dt ≤ 4C 2 T 1/2 p i=m+1 T 0 t 0 (t -s) α-3 2 |f i (s)| 2 dsdt = 4C 2 T 1/2 p i=m+1 T 0 |f i (s)| 2 T s (t -s) α-3 2 dtds ≤ 4C 2 T α α -1 2 p i=m+1 T 0 |f i (s)| 2 ds.
Thus,

B p ≤ 4C 2 T α α -1 2 p i=m+1 T 0 |f i (s)| 2 ds (3.
∥y β p (t) -y β m (t)∥ 2 L 2 ((0,T );H 1 0 (Ω)) ≤ A p + B p ≤ 4C 2 T 2α-1 β 2 (2α -1) p i=m+1 λ i |y 1 i | 2 + 4C 2 T α β 2 (1 -α)(2α -1) + 4C 2 T α α -1 2   p i=m+1 T 0 |f i (s)| 2 ds   .
Therefore,

∥y β p (t) -y β m (t)∥ L 2 ((0,T );H 1 0 (Ω)) ≤ 2C β T 2α-1 2α -1   p i=m+1 λ i |y 1 i | 2   1/2 + 4C 2 T α β 2 (1 -α)(2α -1) + 4C 2 T α α -1 2   p i=m+1 T 0 |f i (s)| 2 ds   1/2 .
(3.21) normyb In view of Equation (3.13), we have

I 1-α (y β p (t) -y β m (t)) = p i=m+1 |b i |E α,1 (-λ i t α )w i + p i=m+1 t 0 f i (u)E α,1 (-λ i (t -u) α )du w i ,
from which we deduce that,

∥I 1-α (y β p (t) -y β m (t))∥ 2 H 1 0 (Ω) ≤ a(I 1-α (y β p (t) -y β m (t)), I 1-α (y β p (t) -y β m (t))) ≤ 2 p i=m+1 λ i |b i | 2 E 2 α,1 (-λ i t α ) + 2 p i=m+1 λ i t 0 f i (u)E α,1 (-λ i (t -u) α )du 2 .
If we set

C p = 2 p i=m+1 λ i |b i | 2 E 2 α,1 (-λ i t α ), Z p = 2 p i=m+1 λ i t 0 f i (u)E α,1 (-λ i (t -u) α )du 2 ,
we have from Theorem 2.1, (3.18) and the Cauchy-Schwartz inequality that,

C p ≤ 2C 2 p i=m+1 λ i |b i | 2 ≤ 2C 2   2 β 2 p i=m+1 λ i |y 1 i | 2 + 2C 2 T 1-α β 2 (1 -α) p i=m+1 T 0 |f i (u)| 2 du   ≤ 4C 2 β 2 p i=m+1 λ i |y 1 i | 2 + 4C 4 T 1-α β 2 (1 -α)   p i=m+1 T 0 |f i (u)| 2 du   and Z p ≤ 2 p i=m+1 λ i t 0 E 2 α,1 (-λ i (t -u) α )du t 0 |f i (u)| 2 du ≤ 2C p i=m+1 t 0 (t -u) -α du t 0 |f i (u)| 2 du ≤ 2Ct 1-α 1 -α p i=m+1 t 0 |f i (u)| 2 du .
Using the estimations of C p and Z p , we obtain

∥I 1-α (y β p (t) -y β m (t))∥ 2 H 1 0 (Ω) ≤ 4C 2 β 2 p i=m+1 λ i |y 1 i | 2 + 4C 4 T 1-α β 2 (1 -α) + 2Ct 1-α 1 -α   p i=m+1 t 0 |f i (u)| 2 du   . Thus, sup t∈[0,T ] ∥I 1-α (y β p (t) -y β m (t))∥ H 1 0 (Ω) ≤ 2C β   p i=m+1 λ i |y 1 i | 2   1/2 + 4C 4 T 1-α β 2 (1 -α) + 2CT 1-α 1 -α   p i=m+1 t 0 |f i (u)| 2 du   1/2 . (3.22) supI As y 1 ∈ H 1 0 (Ω) and f ∈ L 2 (Q), lim m,p→+∞   p i=m+1 λ i |y 1 i | 2   1/2 = lim m,p→+∞   p i=m+1 t 0 |f i (u)| 2 du   1/2 = 0.
∥I 1-α (y β p (t) -y β m (t))∥ H 1 0 (Ω) = 0.
Consequently, (y β m ) and (I 1-α y β m ) are Cauchy sequences in L 2 ((0, T ); H 1 0 (Ω)) and C([0, T ], H 1 0 (Ω)) respectively. This implies that

y β m → y β in L 2 ((0, T ); H 1 0 (Ω)), (3.23) limy1 and I 1-α y β m → ξ in C([0, T ]; H 1 0 (Ω)). Since y β ∈ L 2 ((0, T ), H 1 0 (Ω)
) and I 1-α y β are continuous, we have ξ = I 1-α y β and

I 1-α y β m → I 1-α y β in C([0, T ]; H 1 0 (Ω)). ( 3

.24) limI1

Step 3: We show that y β satisfies (3.4) -(3.5). Let φ ∈ D(0, T ) and µ ≥ 1 an integer. Then, from (3.15), we have for all m ≥ µ,

T 0 (f (t), v) L 2 (Ω) φ(t)dt = T 0 D α RL (y β m (t), v) L 2 (Ω) φ(t)dt + T 0 a(y β m (t), v)φ(t)dt, ∀v ∈ V µ ,
which according to Corollary 2.1 implies that,

T 0 (f (t), v) L 2 (Ω) φ(t)dt = - T 0 (y β m (t), v) L 2 (Ω) D α C φ(t)dt + T 0 a(y β m (t), v)φ(t)dt, ∀v ∈ V µ .
Therefore, passing to the limit and using (3.23), we obtain

T 0 (f (t), v) L 2 (Ω) φ(t)dt = - T 0 (y β (t), v) L 2 (Ω) D α C φ(t)dt + T 0 a(y β (t), v)φ(t)dt, ∀v ∈ V µ . Since ∪ µ≥1 V µ is dense in H 1 0 (Ω) because (w i ) is a base of H 1 0 (Ω), we have for all v ∈ H 1 0 (Ω) that T 0 (f (t), v) L 2 (Ω) φ(t)dt = - T 0 (y β (t), v) L 2 (Ω) D α C φ(t)dt + T 0 a(y β (t), v)φ(t)dt, ∀v ∈ H 1 0 (Ω).
Using, once again, Corollary 2.1, we can write

T 0 (f (t), v) L 2 (Ω) φ(t)dt = T 0 D α RL (y β (t), v) L 2 (Ω) φ(t)dt + T 0 a(y β (t), v)φ(t)dt, ∀v ∈ H 1 0 (Ω).
This implies that for all v ∈ H 1 0 (Ω), (f (t), v) L 2 (Ω) φ(t) = D α RL (y β (t), v) L 2 (Ω) φ(t) + a(y β (t), v)φ(t), ∀t ∈ (0, T ). From (3.24), we have I 1-α y β m (0) → I 1-α y β (0) in H 1 0 (Ω), and

I 1-α y β m (T ) → I 1-α y β (T ) in H 1 0 (Ω). But I 1-α y β m (T ) + βI 1-α y β m (0) = m i=1 y 1 i w i → +∞ i=1 y 1 i w i = y 1 .
Thus,

I 1-α y β (T ) + βI 1-α y β (0) = y 1 .
To complete the proof of Theorem 3.1, we need to prove Equation (3.7) and Equation (3.8). Since y β is the solution of (3.4)-(3.5), we have

y β (t) = +∞ i=1      y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du β + E α (-λ i T α ) t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i .
Proceeding as above for estimations on y β m , we can prove that there exists a constant C > 0 such that

∥y β (t)∥ L 2 ((0,T );H 1 0 (Ω)) ≤ 2C β T 2α-1 2α -1   +∞ i=1 λ i |y 1 i | 2   1/2 + 4C 2 T α β 2 (1 -α)(2α -1) + 4C 2 T α α -1 2   +∞ i=1 T 0 |f i (s)| 2 ds   1/2 and sup t∈[0,T ] I 1-α y β (t) H 1 0 (Ω) ≤ 2C β   +∞ i=1 λ i |y 1 i | 2   1/2 + 4C 4 T 1-α β 2 (1 -α) + 2CT 1-α 1 -α   +∞ i=1 T 0 |f i (s)| 2 ds   1/2
, from which we deduce, respectively, Equation

.

Convergence results

convergence

In this section we provide some convergence results.

conv_I_y1 Theorem 4.1 For all y 1 ∈ H 1 0 (Ω), we have

lim β→0 ∥I 1-α y β (T ) -y 1 ∥ = 0.
That is I 1-α y β (T ) converges to y 1 in H 1 0 (Ω).

Proof. Since y 1 ∈ H 1 0 (Ω), we know that

∀ϵ > 0, ∃N ϵ ∈ N such that +∞ i=Nϵ+1 λ i |y 1 i | 2 < ϵ 2 . Also, since f ∈ L 2 (Q), we know that ∀ϵ > 0, ∃N ϵ ∈ N such that +∞ i=Nϵ+1 T 0 |f i (s)| 2 ds < ϵ 2 .
Let ϵ > 0 and choose N > 0 such that

+∞ i=N +1 λ i |y 1 i | 2 < ϵ 2 and +∞ i=N +1 T 0 |f i (s)| 2 ds < ϵ 2 .
Then, we have

∥I 1-α y β (T ) -y 1 ∥ 2 H 1 0 (Ω) = a -βI 1-α y β (0), -βI 1-α y β (0) = β 2 +∞ i=1 λ i |b i | 2 ≤ A + B,
where

A = 2β 2 +∞ i=1 λ i |y 1 i | 2 (β + E α (-λ i T α )) 2 . B = 2β 2 +∞ i=1 λ i (β + E α (-λ i T α )) 2 T 0 E α (-λ i (T -u) α )f i (u)du 2 .
We firstly have,

A = 2β 2 N i=1 λ i |y 1 i | 2 (β + E α (-λ i T α )) 2 + 2β 2 +∞ i=N +1 λ i |y 1 i | 2 (β + E α (-λ i T α )) 2 ≤ β 2 N i=1 2λ i |y 1 i | 2 E 2 α (-λ i T α ) + 2 +∞ i=N +1 λ i |y 1 i | 2 ≤ β 2 N i=1 2λ i |y 1 i | 2 E 2 α (-λ i T α ) + ϵ.
Secondly, using the Cauchy-Schwartz inequality, we can write

B = 2β 2 +∞ i=1 λ i (β + E α (-λ i T α )) 2 T 0 E α (-λ i (T -u) α )f i (u)du 2 ≤ 2β 2 +∞ i=1 C 2 T 1-α (1 -α)(β + E α (-λ i T α )) 2 T 0 |f i (u)| 2 du ≤ β 2 N i=1 2C 2 T 1-α (1 -α)E 2 α (-λ i T α ) T 0 |f i (u)| 2 du + C 2 T 1-α (1 -α) ϵ.
Finally, using the estimations of A and B, we have

∥I 1-α y β (T ) -y 1 ∥ 2 H 1 0 (Ω) ≤ β 2   N i=1 2λ i |y 1 i | 2 E 2 α (-λ i T α ) + N i=1 2C 2 T 1-α (1 -α)E 2 α (-λ i T α ) T 0 |f i (u)| 2 du   + 1 + C 2 T 1-α (1 -α) ϵ. Since N i=1 2λ i |y 1 i | 2 E 2 α (-λ i T α ) + N i=1 2C 2 T 1-α (1 -α)E 2 α (-λ i T α ) T 0 |f i (u)| 2 du < ∞,
we choose β such that

β 2 < ϵ   N i=1 2λ i |y 1 i | 2 E 2 α (-λ i T α ) + N i=1 2C 2 T 1-α (1 -α)E 2 α (-λ i T α ) T 0 |f i (u)| 2 du   -1
.

Theorem 4.2 Supposing there exists ϵ ∈ (0, 2) such that

D = 2 +∞ i=1 λ i |y 1 i | 2 E ϵ α (-λ i T α ) + 2C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (s)| 2 ds E ϵ α (-λ i T α )
converges, then ∥I 1-α y β -y 1 ∥ H 1 0 (Ω) converges to zero with order ϵ -2 β ϵ . Proof. Let ϵ ∈ (0, 2) such that D converges and k ∈ (0, 2). We fix a natural integer i, and define

g i (β) = β k [β + E α (-λ i T α )] 2 .
Differentiating g i with respect to β, we obtain

g ′ i (β) = (k -2)β k + kβ k-1 E α (-λ i T α ) [β + E α (-λ i T α )] 3 = β k-1 × (k -2)β + kE α (-λ i T α ) [β + E α (-λ i T α )] 3 . Observing that g ′ i (β) = 0 if β = 0 or (k -2)β + kE α (-λ i T α ) = 0. We have (k -2)β + kE α (-λ i T α ) = 0 ⇔ β = k 2 -k E α (-λ i T α ).
As g i (β) > 0, g i (0) = 0 and lim

β→+∞ g i (β) = 0. Indeed, lim β→+∞ g i (β) = lim β→+∞ β k β 2 = lim β→+∞ 1 β 2-k = 0.
We know that g i achieves its maximum at

β 0 = k 2 -k E α (-λ i T α )
. Hence, we have,

g i (β) ≤ g i (β 0 ) ⇔ g i (β) ≤ (β 0 ) k [β 0 + E α (-λ i T α )] 2 ⇔ g i (β) ≤ k 2 -k k E k α (-λ i T α ) [β 0 + E α (-λ i T α )] 2 ⇔ g i (β) ≤ k 2 -k k E k-2 α (-λ i T α ).
Since we can write

∥I 1-α y β -y 1 ∥ 2 H 1 0 (Ω) ≤ 2β 2 +∞ i=1 λ i |y 1 i | 2 (β + E α (-λ i T α )) 2 + 2β 2 +∞ i=1 λ i (β + E α (-λ i T α )) 2 T 0 E α (-λ i (T -u) α )f i (u)du 2 = 2β 2-k +∞ i=1 λ i |y 1 i | 2 g i (β) + 2β 2-k +∞ i=1 λ i T 0 E α (-λ i (T -u) α )f i (u)du 2 g i (β), it follows that ∥I 1-α y β -y 1 ∥ 2 H 1 0 (Ω) ≤ 2β 2-k +∞ i=1 λ i |y 1 i | 2 g i (β) + 2β 2-k C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (u)| 2 du g i (β) ≤ β 2-k k 2 -k k   2 +∞ i=1 λ i |y 1 i | 2 E k-2 α (-λ i T α ) + 2C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (u)| 2 du E k-2 α (-λ i T α )   .
If we choose, k = 2 -ϵ (then ϵ = 2 -k), we then obtain

∥I 1-α y β -y 1 ∥ 2 H 1 0 (Ω) ≤ β ϵ 2 -ϵ ϵ 2-ϵ   2 +∞ i=1 λ i |y 1 i | 2 E -ϵ α (-λ i T α ) + 2C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (u)| 2 du E -ϵ α (-λ i T α )   ≤ β ϵ 2 ϵ 2   2 +∞ i=1 λ i |y 1 i | 2 E -ϵ α (-λ i T α ) + 2C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (u)| 2 du E -ϵ α (-λ i T α )   .
Since D converges, there exists a constant K > 0 such that

2 +∞ i=1 λ i |y 1 i | 2 E -ϵ α (-λ i T α ) + 2C 2 T 1-α 1 -α +∞ i=1 T 0 |f i (u)| 2 du E -ϵ α (-λ i T α ) < K.
which implies that

∥I 1-α y β (T ) -y 1 ∥ 2 H 1 0 (Ω) ≤ β ϵ 2 ϵ 2 K = ϵ -2 β ϵ (4K) = ϵ -2 β ϵ K ′ .
It then suffices to take K ′ = 4K to achieve the proof.

Theorem 4.3 For all y 1 ∈ H 1 0 (Ω), the problem compounded in Equation (1.1) has a solution y if and only if the sequence I 1-α y β (0 + ) converges in H 1 0 (Ω). Furthermore, we have that y β converges to y as β tends to zero in L 2 ((0, T ); H 1 0 (Ω)).

Proof. We proceed in two steps.

Step 1: We show that if I 1-α y β (0) converges in H 1 0 (Ω), then the problem (1.1) admits a solution. Assume that lim β→0 I 1-α y β (0) = y 0 exists. Since y 0 ∈ H 1 0 (Ω), we can write

y 0 = +∞ i=1 y 0 i w i where y 0 i = (y 0 , w i ).
Let y the solution of the following equation

   D α RL y(x, t) -∆y(x, t) = f (x, t) in Q, y(σ, t) = 0 on Σ, I 1-α y(x, 0) = y 0 in Ω.
where 1/2 < α < 1. Then, from Theorem 2.2, we know that y ∈ L 2 ((0, T ); H 1 0 (Ω)) is given by

y(t) = +∞ i=1 t α-1 E α,α (-λ i t α )y 0 i + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i .
Moreover, I 1-α y ∈ C([0, T ], H 1 0 (Ω)). Thus, I 1-α y(T ) ∈ H 1 0 (Ω) exists. Now, let t ∈ [0, T ], we have

y β (t) -y(t) = +∞ i=1 I 1-α y β i (0)t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i - +∞ i=1 y 0 i t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i = +∞ i=1 I 1-α y β i (0) -y 0 i t α-1 E α,α (-λ i t α )w i .
consequently,

∥y β -y∥ 2 L 2 ((0,T ),H 1 0 (Ω)) = T 0 a(y β (t) -y(t), y β (t) -y(t))dt = T 0   +∞ i=1 λ i I 1-α y β i (0) -y 0 i 2 t α-1 E α,α (-λ i t α ) 2   dt = +∞ i=1 λ i I 1-α y β i (0) -y 0 i 2 T 0 t 2α-2 E 2 α,α (-λ i t α )dt ≤ C 2 +∞ i=1 λ i I 1-α y β i (0) -y 0 i 2 T 0 t 2α-2 dt ≤ C 2 T 2α-1 2α -1 I 1-α y β (0) -y 0 2 H 1 0 (Ω)
.

This implies that y β converges to y in L 2 ((0, T ); H 1 0 (Ω)) because lim β→0 I 1-α y β (0) = y 0 .

On the other hand, we have

I 1-α y i (T ) = y 0 i E α (-λ i T α ) + T 0 E α (-λ i (T -u) α )f i (u)du, and
I 1-α y β i (T ) = I 1-α y β i (0)E α (-λ i T α ) + T 0 E α (-λ i (T -u) α )f i (u)du.
Hence, we obtain

I 1-α y β (T ) -I 1-α y(T ) 2 H 1 0 (Ω) = +∞ i=1 λ i I 1-α y β i (0) -y 0 i 2 E 2 α (-λ i T α ) ≤ C 2 I 1-α y β (0) -y 0 2 H 1 0 (Ω) .
This implies that I 1-α y β (T ) → I 1-α y(T ) strongly in H 1 0 (Ω) and since, from Theorem 4.1 I 1-α y β (T ) → y 1 strongly in H 1 0 (Ω), the uniqueness of the limit allows us to conclude that I 1-α y(T ) = y 1 and y is a solution of the problem compounded in Equation (1.1).

Step 2: We show that if the problem given by Equation (1.1) admits a solution y then I 1-α y β (0) converges in H 1 0 (Ω). Let y be a solution of the problem associated with Equation (1.1), then as in the proof of existence in Theorem 3.1, we know that y i = (y(t), w i ) L 2 (Ω) is a solution of the ordinary differential equation

D α RL y i (t) + λ i y i (t) = f i (t), t ∈ [0, T ], I 1-α y i (T ) = y 1 i . (4.1) edo2
Using the Laplace transform of the first equation in Equation (4.1), we obtain

y i (t) = I 1-α y i (0)t α-1 E α,α (-λ i t α ) + t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds. (4.2) sol_fvp_1
Observing that

I 1-α (t α-1 E α,α (-λ i t α )) = E α (-λ i t α )
and

I 1-α t 0 (t -s) α-1 E α,α (-λ i (t -s) α )f i (s)ds = t 0 f i (u)E α (-λ i (t -u) α )du,
we have I 1-α y i (t) = I 1-α y i (0)E α (-λ i t α ) + t 0 f i (u)E α (-λ i (t -u) α )du.

and because I 1-α y i (T ) = y 1 i , we can write

I 1-α y i (0)E α (-λ i T α ) + T 0 f i (u)E α (-λ i (T -u) α )du = y 1 i ,
from which, we deduce that I 1-α y i (0) =

y 1 i - T 0 f i (u)E α (-λ i (T -u) α )du E α (-λ i T α
) .

Thus, we can write .

y(t) = +∞ i=1          y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du E α (-λ i T α ) t α-1 E α,α (-λ i t α )         
As I 1-α y(0) ∈ H 1 0 (Ω), we choose N > 0 such as

∀ϵ > 0, +∞ i=N +1 λ i |I 1-α y i (0)| 2 < ϵ 2 .
This means that ∀ϵ > 0,

+∞ i=N +1 λ i y 1 i - T 0 f i (u)E α (-λ i (T -u) α )du E α (-λ i T α ) 2 < ϵ 2 ,
and we have I 1-α y β (0) -I 1-α y γ (0)

2 H 1 0 (Ω) = +∞ i=1 λ i (γ -β) y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du βγ + (β + γ)E α (-λ i T α ) + E 2 α (-λ i T α ) 2 . ≤ (γ -β) 2 (βγ) 2 N i=1 λ i y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du 2 . + (γ -β) 2 (β + γ) 2 +∞ i=N +1 λ i y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du 2 E 2 α (-λ i T α ) ≤ (γ -β) 2 βγ 2 N i=1 λ i y 1 i - T 0 E α (-λ i (T -u) α )f i (u)du 2 + (γ -β) 2 (β + γ) 2 ϵ 2 .
I 1-α y β (0) -I 1-α y γ (0)

2 H 1 0 (Ω) ≤ γ -β βγ 2   2 N i=1 λ i |y 1 i | 2 + 2C 2 T 1-α 1 -α N i=1 T 0 |f i | 2 (u)du   + γ -β β + γ 2 ϵ 2 . ≤ 2 β 2 + 2 γ 2   2 N i=1 λ i |y 1 i | 2 + 2C 2 T 1-α 1 -α N i=1 T 0 |f i | 2 (u)du   + 2ϵ. Since   2 N i=1 λ i |y 1 i | 2 + 2C 2 T 1-α 1 -α N i=1 T 0 |f i | 2 (u)du   < ∞ and lim γ,β→∞ 2 β 2 + 2 γ 2 = 0,
we deduce that lim γ,β→∞ I 1-α y β (0) -I 1-α y γ (0)

H 1 0 (Ω) = 0.
This implies that the sequence I 1-α y β (0) is of Cauchy and thus it converges in H 1 0 (Ω).

Conclusion

In this work, we have considered an ill-posed problem associated with a family of well-posed problems and prove, using spectral methods, that the solutions of the latter problems converge to the solution of the former problem in an appropriate Hilbert space. This analysis is useful if we want to control an ill-posed problem which will be the subject of future work. Moreover, the convergence results obtained can be used to find a numerical solution for problem compounded in Equation (1.1).

  20) CalBP Adding (3.19) to (3.20), we obtain

Then 2 H 1 0

 21 (t) -y β m (t)∥ (Ω) dt = 0 and sup t∈[0,T ]

  s) α-1 E α,α (-λ i (t -s) α )f i (s)ds w i

0 E 0 E

 00 α (-λ i (T -u) α )f i (u)du E α (-λ i T α ) E α (-λ i t α ) + t 0 f i (u)E α (-λ i (t -u) α )du I_fvp Let β, γ > 0.Then, from (3.14), we haveI 1-α y β (0) -I 1-α y γ (0) = (γ -β) y 1 i -T α (-λ i (T -u) α )f i (u)du βγ + (β + γ)E α (-λ i T α ) + E 2 α (-λ i T α )


Fabien Sixdenier

Abstract -The automotive industry progresses towards full electrification of cars. Therefore, more and more current sensors will be used in each car. These sensors must consume little electrical energy and be as cheap as possible due to the cost-sensitive automotive market. The latter leads to the use, where possible, of Hall-effect open-loop sensors composed of cheap magnetic materials, such as Fe-Si alloys for the magnetic circuit. Depending on the wide range of current and bandwidth of the different onboard electrical devices, each sensor has a different magnetic circuit that must be designed quickly but still accurately. This paper aims to facilitate the design work of engineers by providing them with a simulation model that can be easily used in a circuit type simulation software such as LTspice. This model can bring high-fidelity results in terms of sensor linearity and frequency bandwidth which allows the optimization of the sensor before proceeding to the prototyping step, reducing the time and costs of the design phase. Index Terms-SPICE, finite element analysis, magnetic flux density, current measurement.

I. INTRODUCTION

18 S IMULATIONS are of high importance for design purposes Hall effect sensors, Finite Element (FE) simulations are of 21 great help to optimize the magnetic circuit size and shape.

22

All sensor dimensions must be optimized depending on the 23 maximum current to be measured, the electromagnetic envi-24 ronment, and the bus-bar size. The air-gap length is one of the 25 most important parameters to set correctly because it defines, 26 among other things, the linearity limit of the sensor. However, 27 especially simulating several current levels to obtain linearity 28 curves is time-consuming and requires high computing power.

AQ:2

Manuscript received February 24, 2022; accepted April 16, 2022. This AQ:3 work was supported by the collaboration between LEM Tech France and the Ampère Laboratory of the Claude Bernard Lyon 1 University. The associate editor coordinating the review of this article and approving it for publication was Prof. Boby George. (Corresponding author: Atef Lekdim.) Atef Lekdim is with LEM Tech France SAS, 69800 Saint-Priest, France (e-mail: ale@lem.com). AQ:4

Fabien Sixdenier and Riccardo Scorretti are with the University AQ:5 of Lyon, 69361 Lyon, France, also with Université Claude Bernard Lyon 1, 69100 Villeurbanne, France, also with the INSA Lyon, 69621 Villeurbann, France, and also with the Ampère Laboratory, CNRS, UMR5005, École Centrale de Lyon, 69622 Villeurbanne, France (e-mail: fabien.sixdenier@univ-lyon1.fr; riccardo.scorretti@ec-lyon.fr). Adnan Grihe was with LEM Tech France SAS, 69800 Saint-Priest, AQ:6

France. He is now with ??? (e-mail: adnan.grihe@gmail.com). Digital Object Identifier 10.1109/JSEN.2022.3169920

Moreover, to obtain the frequency bandwidth, several FE time-domain simulations must be done (one per frequency). All these simulations can take several days or even weeks to complete.

One faster alternative is to take advantage of the electrical and magnetic analogy, and to model the magnetic circuit based on an equivalent magnetic network. The equivalent magnetic network modelling is widely used in electrical machines [START_REF] Guo | Improved equivalent magnetic network modeling for analyzing working points of PMs in interior permanent magnet machine[END_REF]- [START_REF] Cao | A generalized equiva-445 lent magnetic network modeling method for vehicular dual-permanent-446 magnet Vernier machines[END_REF]. The equivalent magnetic network is also used in less complicated systems as shown in [START_REF] Hane | Hysteresis modeling of magnetic devices 449 based on reluctance network analysis[END_REF] and [START_REF] Sixdenier | Current sensor modeling with a 453 FE-tuned MEC: Parameters identification protocol[END_REF]. The authors in [START_REF] Hane | Hysteresis modeling of magnetic devices 449 based on reluctance network analysis[END_REF] have studied an oblong magnetic circuit without airgap excited by one winding. The equivalent magnetic network is composed of 28 magnetic reluctances which is enough to model local dynamic phenomena. Previously, the authors in [START_REF] Sixdenier | Current sensor modeling with a 453 FE-tuned MEC: Parameters identification protocol[END_REF] have studied a circular magnetic circuit with airgap for a closed-loop Hall effect sensor. Their equivalent magnetic circuit is composed of only 5 magnetic reluctances, because dynamic effects are considered differently (averaging homogenization) than in [START_REF] Hane | Hysteresis modeling of magnetic devices 449 based on reluctance network analysis[END_REF]. In this paper, the method developed in [START_REF] Sixdenier | Current sensor modeling with a 453 FE-tuned MEC: Parameters identification protocol[END_REF] is adapted in a more generalized way and applied to two oblong geometries. The goal of this part of the work is to develop a generic magnetic reluctance network topology to guarantee the reproducibility of the method for a large spectrum of magnetic circuit geometry/dimensions.

The main advantage of this approach is that the equivalent magnetic circuit can be easily implemented in an electrical The next section describes the modelling methodology that we used to model the magnetic circuit with sufficient accuracy.

The third section will present some results given by our model compared to experimental measurements. Before the conclusion part, a table is given to compare the proposed model to others published in literature.

II. MODELLING METHODOLOGY

Our modelling methodology is separated into two main parts. First, the geometric parameters of a magnetic equivalent network are identified at different DC current levels.

Second, with the same network topology and the best identified geometric parameters, dynamic flux circuit tubes and specific components able to work in transient simulations are created.

The components are put in place in the model to determine the sensor linearity, frequency bandwidth and response time.

A. Magnetic Core Geometries

To show that the method is generic, it is applied on two oblong magnetic cores with the same mechanical dimensions, but with different air gap localizations. Fig. 1 presents the studied magnetic cores. They are fabricated with silicon steel tape, wound on an oblong mandrel. The figure at the left exhibits the airgap in the small branch (ASB) and at right the air-gap in the long branch (ALB).

B. First Finite Element Stationary

Simulation on COMSOL

The goal of these preliminary FE stationary simulations is to have an idea about the magnetic flux density distribution along the magnetic circuit and in the external environment measuring range (up to 1200 A). Thus, parametric simulation 93 was done to obtain the magnetic flux density distribution along 94 the mean length of the magnetic circuit and along the air 95 gap. This flux density distribution will be used as a reference 96 by the optimization algorithm. For the rest of this section, 97 we present only the simulation results of the ALB geometry. 98 The conclusions for the ASB geometry are very similar. 99 Fig. 2 shows the ALB model and the mean path used to 100 plot the magnetic flux density. The geometric and physical 101 symmetries allow modelling of only the half of the magnetic 102 circuit.

103 Fig. 3 exhibits the stationary simulation results in term of 104 magnetic flux density distribution along the half mean path as 105 presented in Fig. 2. The simulation was done for 12 current 106 levels (100 to 1200 A). 107 Fig. 3 shows that for all values of current (only 4 out of 108 12 represented), the induction starts from its minimum in the 109 air gap to its maximum at the opposite side of the air-gap. 110 Generic reluctance network topology applied to the ALB geometry. This topology will reproduce the flux density distribution inside the core and in the air-gape and the distribution of the fringing field in the air area around the magnetic core.

The reluctances R gp1 , R gp2 , R gp11 and R gp12 reproduce the air-gap magnetic behaviour. The length of each reluctance is calculated according to the dimensions of the magnetic circuit and the fictitious angles α 1 and α 2 .

D. Parameter Identification

An optimization method is used for the identification of the parameters.

1) Geometric Parameter Identification:

The following geometric parameters are all defined relatively as ratios of cross sections.

K 1 = S f ri1 S cm ; K 2 = K 3 = S f ri2 S cm ; K 4 = S f ri4 S cm K gp1 = K gp2 = S gp1 S cm ; K gp11 = S gp11 S cm ; K gp12 = S gp12 S cm .
S cm is the magnetic circuit cross-section. S fri1 to S fri4 define the fringing flux cross-sections of reluctances R fri1 to R fri4 respectively. S gp1 , S gp2 , S gp11 , S gp12 define the surfaces of reluctances R gp1 , R gp2 , R gp11 , R gp12 that represent the fringing fluxes in the air-gap.

For each DC current level, the simulated curve of the induction distribution in the magnetic circuit from COMSOL is compared to the one given by the reluctance network model. The objective function (minimization criteria) is the quadratic error between these two curves (see [START_REF] Sixdenier | Current sensor modeling with a 453 FE-tuned MEC: Parameters identification protocol[END_REF] for details). The reluctance network is modelled using the "Modified Nodal Analysis" which is used in SPICE type software [START_REF] Litovski | VLSI Circuit Simulation and Optimiza-457 tion[END_REF], [START_REF] Mccalla | Fundamentals of Computer-Aided Circuit Simulation[END_REF]. Concretely, the optimization loop of Fig. 6 changes the 168 geometric parameters (the ratios K 1 , K 2 , K 3 , K 4 , K gp1 , K gp11 169 and K gp12 ) so that the difference between flux density obtained 170 by COMSOL simulation "B COMSOL " and that one calculated 171 from the model "B MODEL " is minimized.

172 First, the optimization algorithm found that K 2 and K 3 173 had a very low value. Doing so, the algorithm confirms that 174 reluctances R fri2 and R fri3 are not necessary for this magnetic 175 circuit shape. Therefore, the corresponding values K 2 , K 3 are 176 removed in Fig. 7.

177 Secondly, we notice that the values of K 1 , K 4 , K gp1 , are The use of the parameters of table I ends up with the following flux density distributions shown in Fig. 8. B MODEL curves have step shapes because they correspond to the flux density in each reluctance of the magnetic circuit and air-gap. For each curve B MODEL in Fig. 8 we have 5 steps which corresponds to Rgp11, Rgp1, R6, R4 and R2 from Fig. 5. The inset in Fig. 8 presents the first two steps (zoom) up to halfway in the airgap. Going from the highest level of induction to the lowest one, each step means that some magnetic fringing flux is diverted.

In order to prove that using mean values for the geometric parameters does not change the accuracy of the model, we calculated the difference of the flux density in the air-gap between COMSOL simulations and the model by using the set of parameters in Table I (mean values) and the model without averaging the values (without mean values). We have chosen the flux density in the air-gap because it is the one that openloop Hall-effect current sensors use for sensing. For each case, the difference versus primary current value is plotted in Fig. 9.

Clearly, there is no difference in the air-gap except for the last current value (1200 A) which is not exaggerated as an error. Thus, mean values of the geometric parameters can be used for any simulated primary current. The same conclusions have been found for the ASB magnetic core. Now that the network topology and the geometrical parameters are identified, we can focus on the second main part of the magnetic circuit).

I E E E P r o o f
H dyn × L m = L m × H stat φ S cm + L m × γ S cm × dφ dt (2)
The latter equation is implemented with a behavioural voltage source (BV) and a resistor as shown in Fig. 11 a.

2) Air Flux Tube: With the employed analogy, a simple capacitor (which represents the permeance of an air flux tube) is sufficient to represent flux paths in the air (airgap and fringing field paths). The component is presented on Fig. 11 b.

3) Winding Component: It is a component that allows the passage from electric quantities (expressed in volts and amperes) to magnetic ones (magnetomotive forces and flux derivatives expressed in A and Wb/s respectively). This component is modelled thanks to two behavioural voltage sources and one resistor (electrical resistance of the winding) as shown in Fig. 11 c.

4) Whole Model of the Open Loop Sensor Magnetic Core:

With all the components made so far, it is very easy to build the whole model of the sensor shown in Fig. 5. The built whole sensor model is presented in Fig. 11 d.

III. RESULTS ON THE WHOLE SENSOR MODEL

The model is validated using three types of measurements:

• The linearity measurements,

• The bandwidth measurements,

• The response time measurements (di /dt).

The linearity is validated on the ALB and ASB sensors, but the bandwidth and the time response are validated with a commercialized LEM sensor HAH3 1200 S07/SP3 [10]. All the validation experiments were made at room temperature.

A. Linearity Validation

The measurements were made using a Data Acquisition System which is an assembly of different devices (DC suppliers, Digital Multi-Meters, High precision current sensor, DC voltage suppliers, …) all managed using LabView. Sequencies can be programmed to lunch automated measurements, for example to apply a specific current profile and/or a specific temperature profile.

The linearity measurement principle is illustrated in Fig. 12.

Indeed, the current measured by the tested sensor is compared to the reference current measured by the high precision sensor.

The applied current profile starts with 100 A to 1100 A with a step of 100 A.

The figures Fig. 13 and Fig. 14 show the comparison between the LTspice simulated linearity and the measured linearity for the ALB and ASB, respectively. Two linearity errors are chosen to evaluate the differences between the measurement and the simulation linearity: 0.5 % and 1 %.

Globally, the curves behave the same way. The model overestimates by 65 A at 0.5 % linearity error and by 90 A at 1 % linearity error. This means that the model can predict the current limit where the linearity errors exceed 0.5 % and 1 % with good accuracy (6 % and 7.9 % respectively) for the ALB core sensors.

The model represents the ASB sensor behaviour too. The model also overestimates by 65 A at 0.5 % linearity error and by 73 A at 1 % linearity error. This means that the model can predict the current limit where the linearity errors exceed 0.5 % Fig. 12. Linearity test bench functioning principle. The DAS applies an increasing current profile through the DC suppliers. The applied current is measured by the high precision current sensor and by the tested sensor through its Hall die placed in the centre of the air-gap. The Hall die was calibrated accordingly to the magnetic circuit. The DAS collects the measurements for data post-processing. and 1 % with good accuracy (7.2 % and 7.9 % respectively) 306 for the ASB core sensors.

B. Bandwidth Validation 308

The measurements were made using the experiment shown 309 in Fig. 15. The AC supplier was set as follows: The impedance analyzer calculates the gain and the phase shift by comparing the injected AC current measured with the high bandwidth sensor to the current measured with the 3 phases sensor. Five sensors were submitted to the test. Thus, in total 12 phases (monophase sensors) were measured.

I E E E P r o o f
The mean curves of the 12 sensors are considered as the reference. The LTspice bandwidth simulation is made in the same conditions as the measurements (except that it is done only once). Fig. 16 and Fig. 17 show, respectively, the Bode gain plot and Bode phase plot for the simulation and the measurement.

The model can provide Bode plots with good accuracy.

Quantitatively, in the case of the Bode gain plot, the model accurately fits the targeted values between -1 and -3 dB which corresponds to frequencies approximately between 20 kHz and 75 kHz. Indeed, in this interval, the frequency error is very low (f is less than 2.1 kHz) as shown in the inset of Fig. 16. For the Bode phase plot, the model precisely fits the targeted values of less than 1 • phase shift which

338

It is to notice that the simulations did not include the Hall 339 die frequency behaviour because its bandwidth response is 340 very large compared to the magnetic circuit. Indeed, [START_REF] Crescentini | Hall-effect current sen-468 sors: Principles of operation and implementation techniques[END_REF] 341 shows that the -3 dB gain frequency is around a few hundreds A Scope was used to visualize the shunt voltage which is the image of the injected pulse current, and the Hall voltage measured by the tested 3 phases sensor. The response time is calculated at 90 % of the shunt voltage that corresponds to 1200 A.

As for the bandwidth measurements, 5 three phases' sensors AQ:7 were submitted to the response time test. Thus, in total 12 phases (monophase sensors) were measured.

The table below summarizes the calculated response times for the 12 measured sensors.

The mean value is around 2 μs. with the measurements for the two magnetic cores (ALB and 396 ASB). The current limit of linearity for 0.5 % or 1 % error 397 is predicted by the model with an accuracy below 10 % for both types of sensors. Secondly, the frequency behaviour of the sensor is also predictable by the established model. It has been found that between -0.5 dB and -3 dB, the model predicts the frequencies with an error below 2 kHz in the range of 20 kHz to 70 kHz. The phase shift diagram predicted by the model gave good results too until 1 • . After that, the error between the model and the measurements diverges, as the phase shift is generally much more difficult to predict with accuracy. Finally, the response time of the sensor can also be estimated by our model with an error below 0.5 μs for the worst case. Now that the model of the oblong magnetic circuit with air gap is validated, engineers can, for any design:

• Simulate all types of current waveforms,

• Identify trends by varying the geometric dimensions (in a small range: air gap, cross-section, …),

• Estimate frequency bandwidths,

• Estimate the response time, in a very short time.

A time consumption estimation is done considering LTspice and COMSOL linearity and bandwidth simulations.

Clearly, using the LTspice model will save much time.

As the proposed model (and the associated protocol) can also be used to describe other devices under the condition that the magnetic core geometry remains close to the oblong shape and to the studied dimensions, engineers can also work on the design of electronic circuits for closed-loop Hall effect sensors. If the geometry and/or the dimensions are completely changed, the identification protocol presented in the paper can be easily reused/readapted because the protocol has been fully automated using the COMSOL "LiveLink Matlab" module. Thanks to the model developed in this paper, engineers can give some preliminary answers during the design phase, before the generally more costly prototyping phase.

  of Open-Loop Hall-Effect Current Sensor for Circuit-Type Simulation Software Analyses Atef Lekdim, Fabien Sixdenier, Member, IEEE, Riccardo Scorretti , and Adnan Grihe

AQ: 1 19

 1 to improve electrical devices. In current sensors, especially 20

Fig. 1 .

 1 Fig. 1. Oblong wound magnetic cores fabricated with Fe-Si tape of 0.23 mm thickness.

(

  fringing flux) to propose a representative reluctance network model. The simulations have been done with COMSOL in 3D mode for different primary DC current levels. Indeed, it is important to validate the reluctance model for complete current

Fig. 2 .

 2 Fig. 2. 3D COMSOL model of ALB magnetic core and the tracing path of the magnetic flux density in the xy plane.

Fig. 3 .

 3 Fig. 3. Magnetic flux density distribution over the half of the mean length L m of the ALB geometry for four simulated DC currents.

Fig. 4 .

 4 Fig. 4. Flux density (color scale) and flux lines of the ALB geometry are presented in the xy plane. The simulation was done at 400 A primary current.

  Fig. 5.Generic reluctance network topology applied to the ALB geometry. This topology will reproduce the flux density distribution inside the core and in the air-gape and the distribution of the fringing field in the air area around the magnetic core.

Fig. 6 .

 6 Fig. 6. Geometric parameter optimization diagram. The geometric parameters are modified according to the COMSOL curves until the tolerance is respected in the minimization criteria.

Fig. 7 .

 7 Fig. 7. Variation of the geometric parameters as a function of the applied primary current I p for the ALB magnetic circuit.

Fig. 8 .

 8 Fig.8. Comparison between B MODEL and B COMSOL distribution for several primary currents Ip using a unique set of geometric parameters of the ALB magnetic circuit.

Fig. 9 .

 9 Fig.9. Flux density error in the air-gap as a function of primary current I P . The error is the difference between COMSOL and MODEL result in the air-gap. The curve in blue is regarding the MODEL without the use of average geometric parameters and the curve in red is regarding the MODEL with the use of average geometric parameters.

219 2 )

 2 Magnetic Material Modelling: The dynamic magnetic 220 material model is based on the following equation: 221 H dyn (B) = H stat (B) + γ × d B dt (1) 222 Assumptions, limits, and rules of use of this model are fully 223 detailed in [9]. The dynamic excitation field (H dyn ) is assumed 224 as the sum of a static field (H stat ) that depends on the value of 225 B (flux density) and the value of the derivative dB/dt multiplied 226 by a dynamic parameter γ that must be identified. 227 In our case, H stat (B) is just a quasi-static anhysteretic non-228 linear function as shown in Fig. 10. 229 The value of γ is determined by comparing the simulated 230 and measured hysteresis loops obtained on an oblong Fe-Si 231 closed magnetic core (without air-gap) with a sinusoidal flux 232 density waveform of amplitude B = 1.7 T at the frequency 233 of 100 Hz. AMH-200K-S of LABORATORIO ELECTRO 234 FISICO was used as a hysteresis measurement equipment. 235 After identification, the obtained dynamic parameter γ 236 is 0.0412. Fig. 10 shows the measured and the optimized 237 hysteresis loops. 238 E. Specific Components Implementation in LTspice 239 Three specific components were built, and symbols were 240 created on LTspice software to model the whole sensor more 241 easily. These components are the magnetic circuit flux tube, 242 the air flux tube, and the winding.

243 1 )Fig. 10 .

 110 Fig. 10. BH loop measured at 100 Hz and at imposed sinusoidal flux density (1.7 T) compared to optimized one. The Anhysteretic curve corresponds to the mean curve between the ascendant and descendant parts of the measured BH loop.

Fig. 11 .

 11 Fig. 11. LTspice magnetic components. a: magnetic flux tube; b: air flux tube; c: winding component. They are assembled following the proposed network topology to build the LTspice model of open loop Hall sensor presented in d. All the parameters (including the geometric and dynamic parameters) are introduced in their corresponding magnetic components thanks to ".param" option.

Fig. 13 .

 13 Fig. 13. Linearity for ALB core at I Pn = 600 A. The magnetic circuit of an open loop Hall effect sensor is always designed according to the nominal current I Pn and maximal current I Pmax of a given application.

Fig. 14 .

 14 Fig. 14. Linearity for ASB core at IP n = 600 A. The magnetic circuit of an open loop Hall effect sensor is always designed according to the nominal current I Pn and maximal current I Pmax of a given application.

Fig. 15 .••

 15 Fig. 15. Bandwidth measurement functioning principle. The DC voltage supplier is used to feed up the tested sensor.

  313

Fig. 16 .

 16 Fig.[START_REF] Giraud | Modeling quasi-487 static magnetic hysteresis: A new implementation of the play model 488 based on experimental asymmetrical B(H) loops[END_REF]. Bode gain plot and frequency difference. The inset shows the frequency difference between measurement and simulation as a function of the gain.

Fig. 17 .

 17 Fig. 17. Bode phase plot and frequency difference. The inset shows the frequency difference between measurement and simulation as a function of the phase shift.

Fig. 18 .

 18 Fig. 18. Response time experiment. The DC voltage supplier is used to feed up the tested sensor.

333

  error is less than 250 Hz as shown in the inset of Fig.17

  .

334Beyond 1 •

 1 phase shift, the accuracy of the model decreases 335 (frequency difference increases). The phase is generally much 336 more difficult to predict than the gain, so this result is not 337 surprising.

342 of kHz and the 1 •

 1 phase shift frequency is around a few tens 343 of kHz.

344C.•

  Fig. 19. Current step signal scop screenshot of the sensor 12. The red rectangle shows the time response at 90% of the primary signal.

Fig. 20 .

 20 Fig. 20. Simulated di/dt step signal of the LTspice model. At 90% of 1200 A, the response time is 2.5µs.

Fig. 19

 19 Fig. 19 shows an oscilloscope screenshot of a dv/dt (image of di /dt) measurement of sensor 12. The blue curve is the shunt voltage of the primary current and the cyan curve is the Hall output voltage of the tested sensor. As shown with cursors, the sensor 12 responds with t r = 2.4 μs. The output Hall voltage is proportional to the flux density in the airgap. Neglecting the response time of the Hall probe, the response time of the LTspice model is evaluated by measuring directly the induction in the air-gap. Indeed, the Hall die

387

  the geometrical parameters to identify are independent of the 388 current level. Thus, only one set of parameters is necessary to 389 simulate the sensor's response for any current amplitude and 390 waveform. 391 The added value of the specific components implementation 392 in LTspice (dynamic flux tubes, air tubes and winding compo-393 nents) allows predicting exact results before the prototyping 394 phase. First, the model linearity results show good agreement 395

  Atef Lekdim received the master's and Ph.D. degrees in electrical engineering from Claude Bernard Lyon 1 University in 2013 and 2017, respectively. He is currently working as a Research and Innovation Engineer with LEM Tech France. LEM Company is one of the leaders in electrical sensing devices. His main research interest includes magnetic modeling and magnetic material characterization for current sensor applications. Fabien Sixdenier (Member, IEEE) received the B.S., M.S., and Ph.D. degrees in electrical engineering from Université Claude Bernard Lyon 1, France, in 2000, 2002, and 2005, respectively. Since 2006, he has been an Associate Professor with Université Claude Bernard Lyon 1. He is currently working on characterization and modeling of soft magnetic materials. In particular, he developed several magnetic hysteresis models implemented in finite elements and circuit type (SPICE) software that can take into account temperature. He is also working on middle frequency transformers design methodology. Dr. Sixdenier was a recipient of the Best Paper Award of the Electrical Power and Energy Conference (EPEC) in 2015. AQ:10 Riccardo Scorretti received the M.S. degree in informatics engineering from the University of Florence, Italy, in 1999, and the Ph.D. degree in electrical engineering from the École Centrale de Lyon, France, in 2003. Since 2005, he has been a Researcher with CNRS. He is currently working on numerical modeling applied to numerical dosimetry and magnetic materials. Adnan Grihe received the M.Sc. degree in electrical energy from Franche Comte University, France, in 2020. He is currently working on the development of electrical certification training for engineers and technicians working in the field of hybrid and electric vehicles. He is interested in characterization and modeling of soft magnetic materials and renewable energy.

TABLE II SUMS

 II UP THE RESPONSE TIMES FOR THE 12 SENSORScorresponds approximately to 1.2 kHz. Indeed, the frequency

  μs, the model deviates by only 0.5 μs for the worst case,

		TABLE III
		COMPARISON OF THE PROPOSED MODEL TO PREVIOUSLY PUBLISHED ONES
		TABLE IV
		SIMULATION TIME CONSUMPTION
	371	response time is depending on the spinning frequency of the
	372	bias current [11]-[14] which is greater than 8 MHz.
	373	Fig. 20 shows the LTspice model response to the di /dt step
		signal. At 90% of 1200 A, the response time is 2.5μs.
	376
	377	which is highly precise. The response time accuracy is mainly
	378	related to the dynamic parameter γ . Thus, the identification
	379	must be done with high precision.
	380	Table III exhibits a qualitative comparison between the
	381	proposed model and other magnetic models. Only few models
	382	were found in literature that includes the magnetic circuit
	383	behaviour.
	384	IV. CONCLUSION
	385	In this paper, a generic reluctance network model of the
	386	magnetic core of open-loop Hall effect current sensors with a
		dedicated identification protocol is presented. It is shown that

374

If we consider the measured response time mean value 375 of 2

funding or financial support for the research of this article. AQ:3 = Please provide the expansion of the acronym LEM for your funding agency. Providing the correct acknowledgment will ensure proper credit to the funder.