Member, IEEE Alexander Medvedev

Student member, IEEE Viktor Bro

Member, IEEE Rosane Ushirobira

Linear time-invariant discrete delay systems in Laguerre domain

Keywords: Delay systems, linear systems, modeling

The paper provides the formulae connecting the Laguerre spectrum of the output signal of a linear discretetime time-invariant system to the Laguerre spectrum of its input signal and presents a standard state-space system description. The Laguerre domain system representation is meaningful when the input signal is square-summable and the system is stable. A certain type of polynomials arising in the evaluation of the output spectrum due to the presence of time delay is defined and key properties of these are investigated. The polynomials are characterized by a three-term recurrence relation that also facilitates their numerically reliable calculation.

I. INTRODUCTION

The problem of representing a signal as a point in a functional space lies in the heart of mathematical modeling. The coordinates of the point correspond then to the coefficients of a series approximating the signal. The ambition to keep the dimensionality of the functional space parsimonious immediately leads to the use of orthonormal functional bases and extensive tools have been developed in this classical field of mathematics during the last two centuries [START_REF] Szegő | Orthogonal Polynomials[END_REF].

Since a linear time-invariant dynamical system is completely characterized by its impulse response, functional orthonormal bases can be equally applied to system modeling [START_REF] Heuberger | Modelling and Identification with Rational Orthogonal Basis Functions[END_REF]. When a basis is spanned by functions with parameters, a priori information about the system dynamics can be utilized for selecting values of the functional basis parameters so that the system behavior is captured with few terms, i.e. the resulting approximation is sparse.

Often, mathematical modeling is performed to design control or estimation algorithms by solving, generally speaking, an optimization problem. The numerical properties of the model are crucial for obtaining a robust design solution.

Model balancing involves an orthogonalization procedure and constitutes an effective way of improving numerical properties of state-space realizations as systems of difference or differential equations [START_REF] Laub | Computation of system balancing transformations -and other applications of simultaneous diagonalization algorithms[END_REF]. Yet other operators than difference and differential can be employed to capture the system dynamics. The quest for an alternative parametrization to improve numerical properties of models has resulted in e.g. δ -operator [START_REF] Middleton | Digital Control and Estimation: A Unifying Approach[END_REF] and γ-operator [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Problems: Accuracy Aspects[END_REF]. The former stems from the concept of VB and AM were supported in part by the Swedish Research Council, under the grant 2019-04451.

VB and AM are with Information Technology, Uppsala University, Box 337, Uppsala SE-751 05, Sweden (e-mail: {viktor.bro,alexander.medvedev}@it.uu.se).

RU is with Inria, Univ. Lille, CNRS, UMR 9189 -CRIStAL, F-59000 Lille, France (e-mail: rosane.ushirobira@inria.fr) unifying the discrete-and continuous-time system frameworks, in particular by obtaining sampled-data models that converge to their continuous counterparts when the sampling time tends to zero. The latter operator has been originally introduced as a means of mitigating finite word length effects in digital implementation but has a wider field of potential applications.

As already mentioned, orthonormal functional bases with parameters offer a mathematically strict way of approaching the matter of favorable numerical properties in modeling. A sound choice for a basis representing linear time-invariant (LTI) systems in discrete and continuous time is the set of Laguerre functions [START_REF] Clement | Laguerre functions in signal analysis and parameter identification[END_REF] since they are orthogonalized polynomially weighted exponents [START_REF] Wahlberg | System identification using Laguerre models[END_REF], [START_REF] Mäkilä | Laguerre and Kautz shift approximations of delay systems[END_REF]. By mimicking the general form of a solution to an LTI system, Laguerre functions allow adjusting the effective time interval captured by the model through the value of a free parameter. The Laguerre series can be used to represent systems, as approximations of the transfer functions as operators, but also to describe input and output signals that decay at infinity. When the input-output signal pair is given in terms of Laguerre series, one speaks of a Laguerre domain system description.

The introduction of a continuous time-delay operator in a mathematical model results in an infinite-dimensional dynamics. Conventionally, time-delay systems are represented by differential-difference equations [START_REF] Fridman | Discrete-time delay systems[END_REF]. In the Laguerre domain, a time-delay system is given by a difference equation, but it is still infinite-dimensional [START_REF] Hidayat | Laguerre domain identification of continuous linear time delay systems from impulse response data[END_REF]. In discrete time, a time delay increases the order of the state-space realization as well as imposes a specific block structure on the system matrix [START_REF] Fridman | Discrete-time delay systems[END_REF]. Since difference equations are inherently formulated in terms of the discrete-time shift operator, pure time delays in such a model are a matter of parametrization. Nevertheless, the numerical properties of a model are essential and tools for improving them are desired.

Continuous LTI systems with input or output delay were previously treated in Laguerre domain [START_REF] Hidayat | Laguerre domain identification of continuous linear time delay systems from impulse response data[END_REF], [START_REF] Fischer | L 2 time delay estimation by means of Laguerre functions[END_REF]. The corresponding problem in discrete time has been addressed in full only recently [START_REF] Bro | Laguerre-domain modelling and identification of linear discrete-time delay systems[END_REF]. Both in continuous and discrete time, certain types of polynomials arise in Laguerre domain description of time-delay systems. In continuous time, the polynomials are readily recognized as associated Laguerre polynomials [START_REF] Szegő | Orthogonal Polynomials[END_REF], whereas the polynomials appearing in the discrete-time system description are not, to the best knowledge of the authors, covered in the literature and, therefore, are studied in the present contribution.

The main result of this paper is in formulating the distinctive properties of the polynomials arising in the Laguerre domain representation of discrete LTI systems with delay and completely characterizing the polynomial class at hand with a three-term relationship. A binomial identity that can also be of use in other areas (see Lemma 1) is proven as a partial step to the main result.

II. THE DISCRETE LAGUERRE FUNCTIONS

In the Z-domain, the discrete Laguerre functions can be expressed as

L j (z) = √ 1 -p z - √ p T j (z), T (z) 1 - √ pz z - √ p , (1)
for all j ∈ N, where the constant 0 < p < 1 is the discrete Laguerre parameter. Let H 2 d be the Hardy space of analytic functions on the complement of the unit disc, that are square-integrable on the unit circle. The functions L j ( j ∈ N) constitute an orthonormal complete basis in H 2 d with respect to the inner product

W,V 1 2πi D W (z)V (z) dz z , (2)
where V (z) = V (z -1 ) and the contour integral is taken over the unit circle D.

From (1), the complete set of Laguerre functions can be obtained in terms of the discrete-time Laguerre shift operator T by the recursion formula:

L j+1 (z) = T (z)L j (z), ∀ j ∈ N.
Clearly, the orthogonality property of the Laguerre functions with respect to (2) is invariant to (forward or backward) time shift, so the basis

L j (z; d) = z d L 0 (z)T j (z), j = 0, 1, . . . , d ∈ Z, (3)
is also orthonormal. This alternative definition (with d = 1) is used for instance in [START_REF] Fischer | Laguerre shift identification of a pressurized process[END_REF]. Furthermore, the j-th Laguerre coefficient of W ∈ H 2 d is evaluated as the projection of W onto L j

w j = W, L j , ∀ j ∈ N, (4)
and the set {w j } j∈N is referred to as the Laguerre spectrum of W . This parallels the notion of the Fourier spectrum that is obtained by projecting a function onto a set of harmonic functions, while the Laguerre spectrum is calculated by projecting it on a set of weighted (real) exponentials. Then, W is completely defined by its Laguerre series or, in Laguerre domain, by the sequence of the Laguerre coefficients, i.e.

W (z) = ∞ ∑ k=0 w k L k , W (z) → w 0 w 1 w 2 . . . .
The Laguerre shift operator T , defined in (1), acts as the right-shift operator q on the sequence of Laguerre coefficients {w k } k∈N , i.e. the Laguerre spectrum

T (z)W (z) = ∞ ∑ k=0 w k T (z)L k , T (z)W (z) → 0 w 0 w 1 . . . , or q w 0 w 1 w 2 . . . = 0 w 0 w 1 . . . .
The time domain representations of the Laguerre functions ), the space of square-summable sequences defined for non-negative arguments, where Z -1 denotes the inverse Ztransform.

j = Z -1 L j (z) ( j ∈ N) yield an orthonormal basis in l 2 = l 2 [0, ∞

III. LAGUERRE DOMAIN SYSTEM DESCRIPTION

A. Problem formulation Consider the LTI system with an input delay

x(t + 1) = Ax(t) + Bu(t -τ), ∀t ∈ R + (5) y(t) = Cx(t),
where x : R → R n , A, B, C are real matrices of suitable dimensions, x(0) = 0, and τ ∈ N. The problem at hand is to calculate the Laguerre spectrum of the output {y j } j∈N in terms of A, B, C, and τ, given the Laguerre spectrum of the input {u j } j∈N .

Proposition 1. Consider the discrete linear time-invariant system in (5), with y, u ∈ l 2 . For the input Laguerre spectrum {u j } j∈N , the coefficients of the output Laguerre spectrum {y j } j∈N are given by

y j = (1 -p)C(I - √ pA) -2 j-1 ∑ k=0 T (A) j-k-1 Bv k + + √ pC(I - √ pA) -1 Bv j , (6)
where

v j = (1 -p)(-1) j-τ j-1 ∑ k=0 L (τ) j-k ( √ p)u k + √ p τ u j , (7)
L (τ) m ( √ p) (- √ p) m-τ τ-1 ∑ n=0 m + n n m -1 τ -n -1 (-p) n ,
and it is agreed that n k = 0 for k > n by definition. Proof. The proof is obtained by observing that ( 5) is composed of two blocks in series, one constituting an LTI system without explicit time delay and another that is a pure delay. The Laguerre-domain state-space equations for (5) are wellknown, see [START_REF] Nurges | Laguerre state equations for multivariable discrete systems[END_REF] and [START_REF] Fischer | Laguerre shift identification of a pressurized process[END_REF]. The Laguerre-domain description in [START_REF] Wahlberg | System identification using Laguerre models[END_REF] was recently obtained in [START_REF] Bro | Laguerre-domain modelling and identification of linear discrete-time delay systems[END_REF].

The sum in ( 7) is a (discrete) convolution of the polynomials L (τ) m with the Laguerre spectrum of the input signal. In the convolution operation, each Laguerre coefficient of the input u k , k = 0, . . . , j -1 is weighted by the polynomial L (τ) j-k and, therefore, all the input coefficients contribute to the value of v j . This is in contrast with the finite memory of the delay operator in time domain.

Instead of projecting the input and output signal onto the polynomials L j as in [START_REF] Szegő | Orthogonal Polynomials[END_REF], one can consider using the basis L j (z; d) as defined in [START_REF] Laub | Computation of system balancing transformations -and other applications of simultaneous diagonalization algorithms[END_REF]. Then the basis functions can be adapted to the properties of the modeled plant, i.e. (5), when there is a priori information about the delay value. Yet, when the plant ( 5) is subject to initial conditions, a delayed basis might have difficulty approximating the transient response in the output. Moreover, this approach is not feasible for systems with large delays.

IV. COMPARISON TO CONTINUOUS CASE

It is instructive to compare Laguerre-domain descriptions of the time delay operators in continuous and discrete time. A comparison between Laguerre-domain LTI models based on a unifying shift-operator description is provided in [START_REF] Fischer | Laguerre shift identification of a pressurized process[END_REF]. The Laplace transform of the k-th continuous Laguerre function is given by

k (s) = √ 2p c s + p c T k c (s), T c (s) s -p c s + p c ,
for k ∈ N, where p c > 0 represents the Laguerre parameter, and T c is the continuous Laguerre shift operator. Let H 2 be the Hardy space of functions analytic in the open left half-plane. The set { k } k∈N constitutes an orthonormal complete basis in H 2 with respect to the inner product

W,V c 1 2πi ∞ -∞ W (s)V (-s) ds. ( 8
)
The well-known family of associated Laguerre polynomials [START_REF] Szegő | Orthogonal Polynomials[END_REF] is explicitly given by

L m (ξ ; α) = m ∑ n=0 1 n! m + α m -n (-ξ ) n , ∀m ∈ N, ξ ∈ R. (9) Denote L m (ξ ) = L m (ξ ; α)| α=-1
and consider the continuous (pure) delay system

y(t) = u(t -τ c ), τ c ∈ R + , ∀t ∈ R + . (10)
For the system (10), the following relationship holds between the Laguerre coefficients of the input {u j } j∈N and those of the output {y j } j∈N , all evaluated in sense of (8) [START_REF] Hidayat | Laguerre domain identification of continuous linear time delay systems from impulse response data[END_REF], [START_REF] Fischer | L 2 time delay estimation by means of Laguerre functions[END_REF]:

y j = e -κ 2 j ∑ k=1 L k (κ)u j-k + u j , ∀ j ∈ N, where κ = 2pτ c . (11)
For both [START_REF] Wahlberg | System identification using Laguerre models[END_REF] and [START_REF] Fischer | L 2 time delay estimation by means of Laguerre functions[END_REF], the Laguerre coefficients of the output are given by a discrete convolution between the Laguerre coefficients of the input and a polynomial kernel, where the number of terms in the convolution sum is equal to the Laguerre order of the output coefficient. Further, there is a direct term, i.e. y j that statically depends on u j , this property is expected since the signals are related to each other as a whole and not their points or samples. Notably, the convolution operators are also "causal", in the sense that u k , k > j are not involved in right-hand side of ( 7) and [START_REF] Fischer | L 2 time delay estimation by means of Laguerre functions[END_REF].

There are also significant differences between [START_REF] Wahlberg | System identification using Laguerre models[END_REF] and [START_REF] Fischer | L 2 time delay estimation by means of Laguerre functions[END_REF]. In continuous time, the argument of the kernel function incorporates a product of the Laguerre parameter p c and the delay value τ c , thus highlighting the time-scaling effect of p c . In discrete time, the argument of the kernel function is the Laguerre parameter itself and the delay value is not involved. Therefore, the set of the kernel functions is completely defined and can be calculated once and for all feasible delays τ as soon as the value of the Laguerre parameter p is specified.

Another principal difference is in the number of terms in the discrete and continuous case polynomials. In L (τ) m , the number of terms is at most equal to the time delay, whereas the number of terms in L m is defined by the order of the polynomial.

This difference is intimately related to the matter of delay operator dimensionality. The number of terms in the discrete case polynomial cannot exceed the time delay value since this is the dimension of the underlying dynamical system. In continuous time, the delay operator is infinite-dimensional and there is no limit to the number of terms in L m .

Notice that, both the continuous and discrete time, Laguerre shift operators are special cases of the symmetric bilinear transformation that is widely used in systems theory

s → as + b cs -a , a 2 + bc > 0, a, b, c ∈ R.
Bilinear transformations can be introduced formally, e.g. to improve numerical properties of a system identification algorithm [START_REF] Verhaegen | The use of a bilinear transformation of the shift operator in subspace model identification[END_REF] or with motivation from approximation theory and numerical integration, as for δ -operator [START_REF] Middleton | Digital Control and Estimation: A Unifying Approach[END_REF] and γ-operator [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Problems: Accuracy Aspects[END_REF] defined by

δ (z) = z -1 T , γ(z) = 2 T z -1 z + 1 , T > 0.
Obviously, γ-operator is closely related Tustin's discretization method

z = e T s = 1 + T s 2 1 -T s 2 , s = 2 T 1 -z -1 1 + z -1 ,
where T is interpreted as the sampling time. An important observation that can be made from ( 6) is that the convolution kernel involves a matrix function of the state matrix A, namely

T (z)| z=A = T (A).
Similarly, T c (A), δ (A), and γ(A) appear in the corresponding state-space representations, see [START_REF] Fischer | Laguerre shift identification of a pressurized process[END_REF], [START_REF] Middleton | Digital Control and Estimation: A Unifying Approach[END_REF], [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Problems: Accuracy Aspects[END_REF].

V. PROPERTIES OF THE POLYNOMIALS

Some examples of the polynomials L (τ)

m for different values of τ and j are provided below for reference.

For τ = 1, the polynomial expression simplifies to

L (1) m = (- √ p) m-1 .
For τ = 2, the first five polynomials are

L (2) 1 = 2 √ p, L (2)
2 = 1 -3 √ p 2 , L (2)
3 = -2 √ p + 4 √ p 3 , L (2)
4 = 3 √ p 2 -5 √ p 4 , L (2)
5 = -4 √ p 3 + 6 √ p 5 .
For τ = 3, the first five polynomials are

L (3) 1 = 3 √ p 2 , L (3)
2 = 3 √ p -6 √ p 3 , L (3)
3 = 1 -8 √ p 2 + 10 √ p 4 , L (3)
4 = -3 √ p + 15 √ p 3 -15 √ p 5 , L (3)
5 = 6 √ p 2 -24 √ p 4 + 21 √ p 6 .

The polynomials L (τ)

m possess interesting properties some of which can be established by inspection.

m+1 (ξ ) = a (τ) m (ξ )L (τ) m (ξ ) + b m L (τ) m-1 (ξ ) (12)
holds with a

(τ) m (ξ ) a (τ) m,1 ξ + a (τ) m,2 ξ -1 a (τ) m,1 = - m + τ m + 1 , a (τ) m,2 = - m -τ m + 1 , b m = - m -1 m + 1 . (13)
Proof. See Appendix B.

For τ = 2 a

(2) m being monic can be settled through normalization, yet ( 13) cannot be reduced to a first-order polynomial. Therefore, the orthogonality of L (τ)

2 (ξ ) = - 4 3 ξ , b 2 = - 1 3 a (2) 3 (ξ ) = - 5 4 ξ - 1 4 ξ -1 , b 3 = - 1 2 a (2) 4 (ξ ) = - 6 5 ξ - 2 5 ξ -1 , b 4 = - 3 5 For τ = 3 a (3) 2 (ξ ) = - 5 3 ξ + 1 3 ξ -1 , b 2 = - 1 3 a (3) 3 (ξ ) = - 3 2 ξ , b 3 = - 1 2 a (3) 4 (ξ ) = - 7 5 ξ - 1 5 ξ -1 , b 4 = - 3
m remains an open question. In continuous time, the associated Laguerre polynomials obey the following three-term relationship

L m+1 (ξ ) = 1 m + 1 (ξ + 2m)L m (ξ ) - m -1 m + 1 L m-1 (ξ ), (14)
and the orthogonality follows immediately. The presence of ξ -1 in ( 13) is a consequence of P2 since the lower powers of ξ have to be canceled when the order of L (τ) m (ξ ) is increased. This follows from that fact that the discrete-time delay operator has finite dimension equal to τ. In the continuous-time polynomials L m (ξ ), the number of terms grows with the polynomial order.

VI. DISCUSSION

A system with a time delay in the input or the output signal constitutes a cascade coupling of two dynamical subsystems with dissimilar properties. When both subsystems are unknown and must be identified from input-output data, the parameter estimation problem becomes challenging [START_REF] Belkoura | Parameters estimation of systems with delayed and structured entries[END_REF].

Transforming a continuous LTI delay system into the Laguerre domain results in an equivalent discrete state-space description of the dynamics without resorting to a sampling procedure. Besides, the model acquires an extra degree of freedom, the Laguerre parameter, that can be utilized to obtain beneficial numerical properties. In the Laguerre domain, the continuous delay operator turns into the discrete convolution operator in [START_REF] Hidayat | Laguerre domain identification of continuous linear time delay systems from impulse response data[END_REF]. Since the convolution kernel is polynomial and depends on the unknown delay value, the system identification problem becomes nonlinear. Despite the deceivingly simple discrete form, the operator is still infinite-dimensional and cannot be expressed in a state-space form with a finite number of states.

In discrete time, the effect of transforming an LTI delay system into the Laguerre domain is twofold. On the one hand, similarly to continuous time, the model can be regularized by employing the Laguerre parameter. On the other hand, the convolution operator in [START_REF] Wahlberg | System identification using Laguerre models[END_REF] does not include the delay value in the argument of the kernel but rather, as the sharp bound on the number of the terms in the kernel polynomials. The convolution operator is finite-dimensional, as follows from the finite number of the polynomial terms, but of unknown order since the time delay is unknown.

In the system identification context, a practical (but probably not the best) way of identifying time-delay systems in the Laguerre domain is gridding over an interval of feasible time delay values and solving an LTI identification problem for each point of the grid. Indeed, since time delay is a non-negative and bounded parameter, as well as the interval of suitable values is usually a priori known in biomedical applications, this approach is demonstrated to produce good results [START_REF] Hidayat | Laguerre domain identification of continuous linear time delay systems from impulse response data[END_REF], [START_REF] Bro | Laguerre-domain modelling and identification of linear discrete-time delay systems[END_REF].

Laguerre function basis has for a long time been used in nonlinear modeling to parametrize kernels in Volterra (convolution) models [START_REF] Wiener | Nonlinear problems in random theory[END_REF], thus leading to so-called Volterra-Laguerre (VL) models. The structure of a Volterra model readily allows for a time delay in the input or the output signal, albeit implicitly. When an estimate of the delay is sought, the parameter has to appear explicitly in the model. This is achieved by attributing the delay to the kernels of the Volterra model and projecting them onto the Laguerre basis. Assuming the kernels are separable and approximating them by truncated Laguerre series conveniently justifies their interpretation as impulse responses of LTI systems, cf. [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Problems: Accuracy Aspects[END_REF]. Continuous VL models with explicit time delay are dealt with in [START_REF] Bro | Identification of continuous Volterra models with explicit time delay through series of Laguerre functions[END_REF] and the discrete counterpart is covered in [START_REF]Continuous and discrete Volterra-Laguerre models with delay for modeling of smooth pursuit eye movements[END_REF], applied to mathematical modeling of the human smooth pursuit system.

VII. CONCLUSION

Laguerre domain descriptions of dynamical systems provide a common structure for modeling of continuous and discrete systems along with enforcing beneficial numerical properties of the models. The utility of this framework for time-delay systems was highlighted and mathematical properties of the resulting models were analyzed, in particular with respect to the special types of polynomials arising in mapping system descriptions from the time domain to the Laguerre domain.

APPENDIX A PROOF OF PROPOSITION 2

To prove the sought result, an auxiliary lemma is required. Lemma 1. The following identity holds for all m, τ ∈ N

(-1) τ τ ∑ n=0 (-1) n m + n n m τ -n = 1. ( 15
)
Proof. First, notice that the terms in the sum in [START_REF] Verhaegen | The use of a bilinear transformation of the shift operator in subspace model identification[END_REF] are all polynomials in m of order τ. Using the factorial formula for the binomial coefficients gives

m + n n m τ -n = (m + n) τ n!(τ -n)! = 1 n!(τ -n)! τ ∑ k=0 s(τ, k)(m + n) k ,
where s(•, •) are the Stirling numbers of the first kind and (x) n is the falling factorial, i.e.

(m + n) τ = (m + n)(m + n -1) . . . (m + n -τ + 1).
Then, for the result to hold, the sum of the coefficients by each non-zero power of m in (15) should be zero, and the free term has to be equal to one. The latter property is checked by calculating

m = 0 : (-1) τ τ ∑ n=0 (-1) n n n 0 τ -n = (-1) 2τ = 1.
The leading coefficient (by m τ )

(-1) τ τ ∑ k=0 τ ∑ n=0 (-1) n n!(τ -n)! s(τ, k) = (-1) τ τ ∑ k=0 s(τ, k) τ ∑ n=0 (-1) n n!(τ -n)! = (-1) τ τ! τ ∑ k=0 s(τ, k) τ ∑ n=0 (-1) n τ n = 0, in virtue of the identity τ ∑ n=0 (-1) n τ n = 0.
The sum of the coefficients by a certain power of m can be obtained by derivation of the left-hand side of (15) with respect to m

d j dm j (-1) τ τ ∑ n=0 (-1) n n!(τ -n)! τ ∑ k=0 s(τ, k)(m + n) k = (-1) τ τ ∑ n=0 (-1) n n!(τ -n)! τ ∑ k=0 s(τ, k) d j dm j (m + n) k = (-1) τ τ ∑ n=0 (-1) n n!(τ -n)! τ ∑ k=0 s(τ, k)(k) j (m + n) k-j = (-1) τ τ! τ ∑ k=0 s(τ, k)(k) j τ ∑ n=0 (-1) n τ n (m + n) k-j . Now, from Corollary 2 in [20], for 1 ≤ j ≤ τ τ ∑ n=0 (-1) n τ n (m + n) k-j = 0,
which completes the proof. Now the claim of Proposition 2 can be proven.

Proof. The result is proven by induction with respect to τ. For τ = 1, it is easily checked that

S (1) m = (-1) 0 m 0 m -1 0 = 1, ∀m.
Since m does not actually appear in the previous expression, it is instructive to evaluate By using the identity

S (2) m = (-1) 1 m 0 m -1 1 + (-1) 0 m + 1 1 m -1 0 = -(m -1) + (m + 1) = 2, ∀m,
m -1 τ -n = m τ -n - m -1 τ -n -1 , one has S (τ+1) m = (-1) τ τ ∑ n=0 (-1) n m + n n m τ -n + (-1) τ-1 τ-1 ∑ n=0 (-1) n m + n n m -1 τ -n -1 . Finally, S (τ+1)
m = S (τ) m + (-1) τ τ ∑ n=0 (-1) n
where the last identity is due to Lemma 1.

APPENDIX B PROOF OF PROPOSITION 3

According to P2, the polynomials L The next coefficient (multiplying ξ m-τ ) is zero in both polynomials because of P1. Further, for recursion [START_REF] Bro | Laguerre-domain modelling and identification of linear discrete-time delay systems[END_REF] to hold, the least power term of a

where m-1 are involved in the right-hand side of [START_REF] Bro | Laguerre-domain modelling and identification of linear discrete-time delay systems[END_REF].

s 1 (k) = - m + τ m + 1 m + k k m -1 τ -1 -k (-1) m-τ+k , s 2 (k) = - m -τ m + 1 m + k + 1 k + 1 m -1 τ -2 -k (-1) m-τ+k+1 , s 3 (k) = - m -1 m + 1 m + k k + 1 m -2 τ -2 -k ( -
Consider first the sum

s 1 (k) + s 3 (k) = (-1) m-τ+k+1 m + 1 m + k k m -1 τ -1 -k (m + τ) +(m -1) m + k k + 1 m -2 τ -2 -k .
By applying the identities

m -1 τ -1 -k = m -1 τ -1 -k m -2 τ -2 -k , m + k k + 1 = m k + 1 m + k k ,
it can be concluded that Making use of the identities Finally, by observing that

m -1 τ -k -1 = m -τ + k + 1 τ -k -1 m -1 τ -k -2 , m + k + 1 k + 1 = m + k + 1 m m + k k + 1 ,
m τ -k -1 = m m -τ + k + 1 m -1 τ -k -1 , m + k + 1 k + 1 = m + 1 k + 1 m + k + 1 k ,
the equality in ( 16) is proven.

P1:m

  All the powers in L (τ) m are either odd or even, which implies that the polynomials themselves are odd or even functions. P2: There are at most τ terms in L (τ) m . For m ≥ τ, all L (τ) m have exactly τ terms. P3: The lowest order of a term in L (τ) m is |m -τ|, the highest is τ + m -2. Therefore, only the polynomials L = τ, for all m ∈ N. Proof. See Appendix A. The polynomials L (τ) m can be obtained in a recursively. Proposition 3. The following three-term recurrence relation L (τ)

5

 5 Notice that Favard's theorem securing orthogonality of a polynomial set through a three-term recurrence relation demands that a(τ) m (ξ ) = ξ + c mfor some c m , is satisfied. While the matter of a (τ)

1 τ

 1 to reveal the mechanism by which the powers of m cancel each othern -1 .

m

  cannot have more than τ terms for any m. The polynomials L (τ) m+1 and ξ L (τ) m comprise terms of power from mτ + 1 to m + τ -1, while the terms in the polynomials L (τ) m-1 and ξ -1 L (τ) m are of power from mτ -1 to m + τ -3. The leading coefficient of L (τ) m+1 (the coefficient of ξ m-τ+1 ) coincides with that of ξ a

  again, the coefficients multiplying ξ m+τ-2 are zero, both in ξ -1 L (τ) m and b m L (τ) m-1 , due to P1. Now, it has to be shown that the relationm + k + 1 k m τ -1k (-1) m-τ+k = s 1 (k) + s 2 (k) + s 3 (k),

  1) m-τ+k , holds for k = 0, . . . , τ -2 (i.e. the powers mτ + 1 to m + τ -3), where the non-zero coefficients of all three polynomials ξ L

s 1 (k)+s 3 ( 1 m - 1 τ - 1 -k m + k k + 1 τ 1 m - 1 τ - 2

 131111112 Now consider the whole right-hand side expression in (12)s 1 (k) + s 2 (k) + s 3 (k) = (-1) m-τ+k+1 m + (k + m + 1) m -(mτ) m + k + 1 k +k .

one arrives to s 1 1 τ - 1

 111 (k) + s 2 (k) + s 3 (k) = (-1) m-τ+k+1 m(k + 1) (m + 1)(mτ + k + 1) m -k m + k + 1 k + 1 .


Rosane Ushirobira

Senior Member, IEEE Denis Efimov

A Simple Frequency Estimator For Power Systems

Fast and accurate frequency estimation of an unbalanced three-phase power system in the presence of measurement offset is a challenging task. A simple to implement and tune method is proposed in this work through a delay-based linearregression framework. Using two challenging scenarios, we verify the performance of our procedure in the presence of severe grid abnormalities, including voltage unbalance, harmonics, frequency jump, and DC offset. Comparative hardware-in-the-loop experimental results demonstrate the suitability of the proposed approach.

I. INTRODUCTION

Frequency plays an important role in monitoring the electric power grid. Various abnormalities such as unbalance, DCoffset, and voltage sag/swell often make the frequency estimation challenging. Moreover, fast estimation is also required to take any potential preventive action. Numerous techniques are available in the literature that deal with fast and accurate frequency estimation. Some popular tools are: phase-locked loop (PLL) ( [START_REF] Mellouli | A third-order MAF based QT1-PLL that is robust against harmonically distorted grid voltage with frequency deviation[END_REF], [START_REF] Rodríguez | Decoupled double synchronous reference frame PLL for power converters control[END_REF]) frequency-locked loop (FLL) [START_REF] Rodríguez | Multiresonant frequency-locked loop for grid synchronization of power converters under distorted grid conditions[END_REF], demodulation [START_REF] Verma | A robust Lyapunov demodulator for tracking of single-/three-phase grid voltage variables[END_REF], least-square [START_REF] Achlerkar | Recursive least squares based adaptive parameter estimation scheme for signal transformation and grid synchronization[END_REF], Kalman filter [START_REF] Amirat | Phasor estimation for grid power monitoring: Least square vs. linear Kalman filter[END_REF], neural network [START_REF] Rahoui | Neural networks based frequency-locked loop for grid synchronization under unbalanced and distorted conditions[END_REF] and so on.

PLLs [START_REF] Rodríguez | Decoupled double synchronous reference frame PLL for power converters control[END_REF] dynamic performance depends on the loop-filter (proportional-integral-type) tuning parameters. A fast dynamic response comes at the cost of reduced harmonics robustness. Total quasi type-1 PLL (TQT1-PLL) [START_REF] Mellouli | A third-order MAF based QT1-PLL that is robust against harmonically distorted grid voltage with frequency deviation[END_REF] uses a proportional controller-based loop-filter, increasing the dynamic response without sacrificing much the harmonic robustness. FLL [START_REF] Rodríguez | Multiresonant frequency-locked loop for grid synchronization of power converters under distorted grid conditions[END_REF] and its many variants [START_REF] Rahoui | Neural networks based frequency-locked loop for grid synchronization under unbalanced and distorted conditions[END_REF] use an integral controller to estimate the frequency. This approach is sensitive to the grid harmonics components. So, additional filters are required. Demodulation [START_REF] Verma | A robust Lyapunov demodulator for tracking of single-/three-phase grid voltage variables[END_REF] and least-square [START_REF] Achlerkar | Recursive least squares based adaptive parameter estimation scheme for signal transformation and grid synchronization[END_REF] use an initial phase-angle-based technique for frequency estimation. The demodulation technique applies additional moving average filters, and adaptive tuning of the filter's window length can be challenging and computationally demanding. Kalman filter [START_REF] Amirat | Phasor estimation for grid power monitoring: Least square vs. linear Kalman filter[END_REF] can provide good robustness to measurement noise. However, it requires accurate information of the process and measurement noise covariance matrices. Moreover, results in [START_REF] Amirat | Phasor estimation for grid power monitoring: Least square vs. linear Kalman filter[END_REF] show that the least-square technique has a better statistical performance than the Kalman filter.

From the brief literature review, it can be found that there is a demand for a simple frequency estimator that can handle Hafiz Ahmed is with the Nuclear Futures Institute, Bangor University, Bangor LL57 1UT, UK (e-mail: hafiz.h.ahmed@ieee.org). Rosane Ushirobira and Denis Efimov are with Inria, Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France (e-mail: {rosane.ushirobira,denis.efimov}@inria.fr). Denis Efimov is also with ITMO University, 49 av. Kronverkskiy, 197101 Saint Petersburg, Russia. H. Ahmed is funded through the Sêr Cymru programme by Welsh European Funding Office (WEFO) under the European Regional Development Fund (ERDF).

various grid abnormalities without requiring additional preand or post-filtering stages. One such simple solution has been proposed in [START_REF] Vediakova | Finite time frequency estimation for multi-sinusoidal signals[END_REF], where the unknown frequency estimation problem has been solved through linear regression. However, it works only for single-phase cases and can not handle DC offset. In this work, we have applied the idea of [START_REF] Vediakova | Finite time frequency estimation for multi-sinusoidal signals[END_REF] to estimate the unknown frequency of an unbalanced three-phase grid and our proposed approach can handle DC offset without any additional filtering stages.

The rest of this letter is organized as follows: the linearregression-based model of the three-phase unbalanced voltages and the design of the unknown frequency estimator are given in Section II, comparative experimental results are in Section III, and finally, Section IV concludes this letter.

II. SIGNAL MODELING AND ESTIMATOR DESIGN

Unbalanced three-phase grid voltages in the stationary reference are given by:

v α (t) = V α0 + V p cos(ωt + φ p ) + V n cos(ωt + φ n ), (1) v β (t) = V β0 + V p sin(ωt + φ p ) -V n sin(ωt + φ n ), (2)
where the subscript p and n denote the positive and negative sequence for i ∈ {p, n}, θ i = ωt + φ i is the total phase, V i , ω, φ i , V α0 , and V β0 represent the amplitude, angular frequency, phase angle, and DC offsets in phases α and β, respectively. Let us first consider the signal v α and its delayed versions v τ α (t) = v α (t -τ ) with τ > 0 as given by:

v τ α (t) = V α0 + V p cos(ω(t -τ ) + φ p ) + V n cos(ω(t -τ ) + φ n ), = V α0 + V p cos(θ p ) cos(ωτ ) + V p sin(θ p ) sin(ωτ ) + V n cos(θ n ) cos(ωτ ) + V n sin(θ n ) sin(ωτ ). (3)
Similarly, the delayed signals

v 2τ α (t) = V α (t -2τ ) and v 3τ α (t) = V α (t -3τ ) are given by: v 2τ α (t) = V α0 + V p cos(θ p ) cos(2ωτ ) + V p sin(θ p ) sin(2ωτ ) + V n cos(θ n ) cos(2ωτ ) + V n sin(θ n ) sin(2ωτ ). (4)
v 3τ α (t) = V α0 + V p cos(θ p ) cos(3ωτ ) + V p sin(θ p ) sin(3ωτ ) + V n cos(θ n ) cos(3ωτ ) + V n sin(θ n ) sin(3ωτ ). (5)
Using ( 1) and ( 3)-( 5), one can obtain that:

v 2τ α (1 + 2 cos(ωτ )) -v τ α (1 + 2 cos(ωτ )) = -v α + v 3τ
α . This leads to the following relationship:

v α -v τ α + v 2τ α -v 3τ α -2 cos(ωτ )(v τ α -v 2τ α ) = 0. (6)
The above calculations are equally applicable to the signal v β and one can write that: 6) and ( 7), the following first-order linear regression model can be obtained:

v β -v τ β + v 2τ β -v 3τ β -2 cos(ωτ )(v τ β -v 2τ β ) = 0. (7) By combining (
y = ϕΩ, (8)
where

y = v α + v β -v τ α -v τ β + v 2τ α + v 2τ β -v 3τ α -v 3τ β , ϕ = 2(v τ α + v τ β -v 2τ α -v 2τ β )
, and Ω = cos(ωτ ). As a result, the unknown grid frequency can be easily estimated from the parameter Ω in eq. ( 8). The unknown parameter Ω in ( 8) can be estimated by applying the gradient method [START_REF] Ioannou | Robust adaptive control[END_REF] and is given by:

Ω = εϕ(t)(y(t) -ϕ(t) Ω), (9)
where ε > 0 is the tuning parameter. From Ω, the unknown frequency can be evaluated as:

ω = arccos( Ω)/τ.
It is well-known that the regressor ϕ is persistently exciting if τ < π/ω [START_REF] Ioannou | Robust adaptive control[END_REF]. As a result, exponential convergence of the estimator ( 9), and consequently, the frequency ω can easily be established.

Note that instead of combining ( 6) and ( 7), two separate frequency estimators can be designed by writing these equations into the linear regression form [START_REF] Vediakova | Finite time frequency estimation for multi-sinusoidal signals[END_REF]. Then, by taking an average of the estimated frequencies, the effect of noise can be reduced. However, this will increase the computational complexity as an additional integrator will be required.

The estimator [START_REF] Ioannou | Robust adaptive control[END_REF] has two parameters to tune: the delay τ and the gain ε. The delay can be selected as τ = T /4, where T is the fundamental signal period. The tuning gain ε has to be selected as a trade-off between fast convergence speed and disturbance rejection property. Through an extensive numerical simulation, it has been found that ε = 10 can be a good choice. III. RESULTS AND DISCUSSIONS dSPACE 1104 board-based hardware-in-the-loop experimental studies are considered in this section. Details of the experimental setup are given in [START_REF] Ahmed | On the enhancement of generalized integrator-based adaptive filter dynamic tuning range[END_REF]. As a comparison technique, TQT1-PLL [START_REF] Mellouli | A third-order MAF based QT1-PLL that is robust against harmonically distorted grid voltage with frequency deviation[END_REF] is used. Parameters of the TQT1-PLL are selected the same as in [START_REF] Mellouli | A third-order MAF based QT1-PLL that is robust against harmonically distorted grid voltage with frequency deviation[END_REF]. Both techniques are discretized using the explicit Euler method with a sampling frequency of 10 kHz.

In Test-I, we consider 2 Hz frequency step change and a sudden addition of 0.05 p.u., -0.06 p.u., and 0.07 p.u. DC offsets in phases a, b, and c, respectively. Experimental results are given in Fig. 1

→ V +1 = 1∠0 • , - → V -5 = 0.04∠45 • , - → V +7 = 0.03∠60 • , - → V -11 = 0.02∠90 • , and - → V +19 = 0.01∠90 • [p.u.
]. This corresponds to a total harmonic distortion of ≈ 0.5%. Experimental results in Fig. 1 (c) and (d) show that the developed method has faster convergence (less than 1 cycle) to the actual value with an excellent harmonic robustness property. Results in Fig. 1 show that the developed method is very suitable as a real-time frequency monitoring tool for smart grid. In addition, it can also be used as the frequency estimator for stationary reference frame-based grid-synchronization methods [START_REF] Rodríguez | Multiresonant frequency-locked loop for grid synchronization of power converters under distorted grid conditions[END_REF], [START_REF] Rahoui | Neural networks based frequency-locked loop for grid synchronization under unbalanced and distorted conditions[END_REF].

IV. CONCLUSION

This letter studied the frequency estimation problem of an unbalanced three-phase electric power grid with DC offset. First, using consecutive delays, the unknown frequency estimation problem has been converted into a problem of a delaybased linear regression. Then, a simple gradient estimator is used to estimate the unknown frequency in real-time. Since the proposed method is linear, it is very simple to implement and tune. Comparative experimental studies demonstrated the suitability of the proposed method in estimating the unknown grid frequency from unbalanced and distorted grid voltage measurements.

  (a) and (b). The proposed technique converged in ≈ 20 msec. while TQT1-PLL took ≈ 40 msec. In Test-II, suddenly the grid become distorted with -

Figure 1 .

 1 Figure 1. Comparative experimental results with TQT1-PLL [1].


Mateo Uldemolins

Andrej Mesaros

Pascal Simon

Effect of Van Hove singularities on Shiba states in two-dimensional s-wave superconductors

Magnetic impurities in a superconductor induce Yu-Shiba-Rusinov (YSR) states inside the superconducting gap, whose energy depends on the strength of the coupling to the impurity and on the density of states (DOS) at the Fermi level. We consider DOS exhibiting a logarithmic or a power-law divergence at the Fermi level due to Van Hove singularities (VHS) and high-order Van Hove singularities (HOVHS), respectively. We find that the energy of the YSR states has the same functional form as in the constant DOS scenario, with the effect of the singularity being an enhancement of the effective coupling constants. In particular, the critical magnetic coupling strength at which the Shiba transition occurs is always lowered by a factor 1/ρ(∆/Ec), where ∆ is the superconducting gap, Ec is the bandwidth, and ρ(E) is the factor in DOS which diverges at E = 0 for a VHS or HOVHS. Further, since the critical magnetic coupling is significantly reduced, a new regime becomes accessible where the transition point is controlled by the non-magnetic coupling constant. Interestingly, the slope of the Shiba energy curve at the Shiba transition is independent of impurity parameters and purely reflects the band structure. Additionally, we find that our main conclusions remain valid even when the Fermi level is not precisely tuned to the Van Hove singularity, but instead lies on an energy range of order the superconducting gap. Our results show that tuning a superconducting material towards a VHS or HOVHS enhances the possibilities for engineering YSR states, and for characterizing the superconductor itself.

I. INTRODUCTION

The appearance of VHS in the DOS as the system is tuned through a Lifshitz transition enhances electronic correlations which in turn can trigger quantum emergent phenomena such as superconductivity [1,2] or charge-density-wave instability [3]. Recent experimental progress showcasing the ability to bring a VHS arbitrarily close to the Fermi level in twisted bilayer graphene (TBG) [4,5] and in heterostructures [6] has reawakened the interest in the field. On the theory side, the notion of VHS was recently extended to HOVHS where the DOS is enhanced from logarithmic to power-law divergence and they were predicted to be realizable by tuning a single parameter in Moiré superlattices [7,8]. Although some progress has been made to characterize HOVHS in interacting systems [9], the interplay of HOVHS and magnetic impurities in conventional superconductors remains unknown.

The classical spin of a magnetic impurity embedded in a superconductor is an example of a pair-breaking defect which leads to the formation of localized excitations inside the superconducting energy gap known as Yu-Shiba-Rusinov (YSR) states [10][11][12][13][14]. Ever since their first observation by Yazdani et al. [15], the advancements in scanning tunneling spectroscopy (STS) and atomic manipulation techniques have prompted a growing interest in the subject (see [16] for a recent review). YSR states have proved to be a useful tool to probe the orbital properties of the magnetic adsorbate [17][18][19][20], and recently, the advent of STS functionalized tips has opened a new * pascal.simon@universite-paris-saclay.fr avenue to study spin-dependent phenomena [21][START_REF] Huang | Spin-dependent tunneling between individual superconducting bound states[END_REF][START_REF] Schneider | [END_REF][24]. Simultaneously many works have proposed YSR as building blocks to engineer exotic topological states of matter [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41] with experimental spectral signatures consistent with topological superconductivity displaying Majorana zero modes [42][43][44][45][46][47] (see [48] for a recent review and refs. therein). Both endeavors require the YSR states to have a large spatial extent as well as a strong coupling to the substrate. While past observations indicate that quasi two-dimensional materials constitute an ideal playground to achieve the former [49,50], a quantitative experimental connection between the YSR state energy and the impurity-substrate hybridization has been only analyzed very recently [51,52] and deserves further study.

In the present work, motivated by the recent discovery of superconductivity in TBG [53], the quest to find YSR states in graphene [54][START_REF] Del Río | Observation of yushiba-rusinov states in superconducting graphene[END_REF][START_REF] Lopez-Bezanilla | [END_REF], and the promising pathway to VHS tunability in these materials, we explore the influence of conventional and high-order VHS on YSR states. We find that tuning the Fermi level to a Van Hove singularity is a powerful way to enhance the coupling of YSR states to the substrate, while at the same time, it provides a possibility to extract information from the DOS singularity itself.

Our paper is organized as follows. In Sec. II, we present our model Hamiltonian for the Shiba impurity coupled to a superconductor with a Fermi level tuned near a VHS or HOVHS. Sec. III contains our analytical and numerical results for both types of Van Hove singularities. Finally, in Sec. IV we provide a short conclusion and discussion of our results. The effect of a small perturbation in the chemical potential away from the singularity is discussed in App. A, and some technical details are presented in App. B.

II. MODEL HAMILTONIAN

We consider a pointlike, isotropic, magnetic impurity on a two-dimensional, s-wave superconductor. The Bogoliubov-de Gennes Hamiltonian (BdG) of the system in the Nambu basis Ψ

= (ψ ↑ , ψ ↓ , ψ † ↓ , -ψ † ↑ ) T reads H = ξ k τ z + ∆τ x + (Kτ z -Jσ z )δ(r -r imp ), (1)
where k and r designate the electron's momentum and position, ξ k is the energy dispersion of the electrons in the superconducting substrate, ∆ is the superconducting gap, K is the amplitude of the non-magnetic scattering potential and J ≡ JS/2 denotes the coupling strength between the electrons and the magnetic impurity with classical spin S at r imp , and Pauli matrices τ i and σ i act on particle-hole and spin space, respectively. In the present study we disregard any quantum effects associated with the magnetic impurity (e.g. Kondo screening), and comment on the validity of the model in the discussion section.

In this work we focus on three expressions for the energy dispersion ξ k defined in the continuum, since the information relevant for the Shiba state comes from the vicinity of the Fermi momentum, and thereby the three energy dispersions are representative of three broad classes of systems: ones with constant DOS at the Fermi energy, ones with a classical VHS, and ones belonging to a family of HOVHS labeled by a parameter ν, as we detail now. Firstly, owing to the toroidal topology of the Brillouin zone, the energy dispersion of two-dimensional systems is always endowed with at least two saddle points, defined as

∇ k ξ k = 0, det D < 0, (2)
where D the Hessian matrix of ξ k . First discussed in the context of phonons [57] and later applied to electronic systems, it is well known that such saddle points yield a logarithmic singularity in the DOS at the corresponding energy, henceforth denoted conventional Van Hove singularity:

ρ(ξ) = 1 2E c log E c ξ , (3)
where ξ is the energy measured from the Fermi level.

Here the prefactor stems from imposing the normalization condition 1 = Ec -Ec ρ(ξ)dξ, and E c is an energy cutoff which is introduced to delimit the energy range over which Eq. ( 3) is a good description of the system's DOS. In our calculations we always assume that E c is the largest energy scale, of the order of the bandwidth. Next, for systems with constant density of states at the Fermi energy, we use the constant

ρ(ξ) = 1 2E c , (4)
which is conveniently normalized with the use of the same bandwidth. Finally, in recent works the notion of Van Hove singularity was extended to high order Van Hove singularities [8], where degenerate saddle points, i.e. fulfilling det D = 0, lead to a power-law divergent DOS,

ρ(ξ) = ν + 1 η + 1 1 E ν+1 c |ξ| ν η if ξ < 0, 1 if ξ > 0. (5)
Here -1 < ν < 0, ν ∈ Q, is the parameter labeling the HOVHS, the E c has the same meaning as above, and η ≡ ρ(-|ξ|) ρ(|ξ|) is the particle-hole asymmetry ratio. Note that we include the case of particle-hole symmetric HOVHS at η = 1, and that the previous two dispersions are also symmetric. The exponent ν can take an infinite number of rational values in its range, each value having a corresponding η value; a set of possibilities is tabulated in Ref. [8].

Having defined the models, we now calculate the energy of the Shiba states looking for in-gap solutions of the Schrödinger equation in particle-hole space.

III. RESULTS AND DISCUSSION

A. Conventional Van Hove singularity

A logarithmic divergence of the DOS at the Fermi level as in Eq. ( 3), yields the following self-consistent expression for the Shiba energy (see Appendix B for details),

E S = ±∆ 1 -J 2 + K 2 4 J 2 + 1 -J 2 + K 2 2 , (6)
where the existence of two symmetric eigenvalues stems from the particle-hole symmetry constraint of the BdG Hamiltonian, and J = f vh (∆, E S )J and K = f vh (∆, E S )K are effective coupling energies with

f vh (∆, E S ) = π 2E c log E c ∆ 2 -E 2 S , (7)
a renormalizing factor proportional to the density of states in Eq. ( 3). An expression for E S of the classical VHS in the strong-coupling limit E S ∆, and for vanishing K was first obtained in Ref. [58]. Remarkably our expression for the Shiba energy in Eq. ( 6) has the same functional form as the one for a constant DOS, namely taking the dispersion in Eq. ( 4) one obtains the standard Shiba energy [12], i.e. Eq. ( 6) with

f const = π 2E c . ( 8
)
An important consequence is that a Van Hove singularity in the DOS at the Fermi level enhances the effective coupling energies in the Shiba problem, since f const < f vh .

In the presence of a magnetic impurity, localizing an electron entails an energy trade-off due to the exchange interaction, therefore the energy of the in-gap excitation E S FIG.

1. Positive branch of the Shiba energy, Eqs. ( 6) and (11), as a function of the magnetic coupling strength J for a constant (green), logarithmically-divergent (red), power-law divergent with ν = -1/3, η = 2 (orange) and power-law divergent with ν = -1/2, η = 1 (blue) DOS. The black dot indicates the value of K/Ec. In the small K regime (a), Jc depends strongly on the underlying DOS, while for large K regime (b), Jc ∼ K. The solid lines indicate the numerical solution of the self-consistent equations ( 6) and ( 11) while the dashed lines represent the zeroth-order approximation in ES.

decreases as J increases (see Fig. 1). This competition triggers the Shiba quantum phase transition at E S = 0, where the in-gap state becomes occupied and it turns into the new superconducting ground state [13,14]. The critical magnetic coupling J c prompting the phase transition reads

J j c = K 2 + 1 f j (∆, 0) 2 , (9)
where j ∈ {const, vh} labels the chosen dispersion. In all our calculations we assume K E c , therefore the relation f const (∆, 0) < f vh (∆, 0) implies that a Van Hove singularity in the DOS at the Fermi level reduces J c with respect to its value with a constant DOS. This remains true even when Fermi level is not exactly tuned to the Van Hove singularity but within an interval of the order of the superconducting gap instead, although the expression for J vh c becomes slightly more complicated (see App. A for details).

The DOS enhancement at the VHS enables one to access a new regime of the impurity system. Namely, the theory is always applicable in a "small K" regime defined by K 1 fj (∆,0) , which in the case of the constant DOS reiterates the theoretical condition K E c . This regime in the VHS case leads to the critical magnetic coupling J c reduction due to the larger f j (∆, 0) [Fig. 1 (a)]. However, for a VHS system, we can consider a "large K" regime, where K 1 f vh (∆,0) without contradicting any assumption of the model, and then the J j c ∼ K would become independent of the DOS [see Fig. 1 (b)].

We note that the VHS introduces a dependence of the J c on the superconducting energy gap ∆, in contrast to the standard Shiba impurity for constant DOS. Specifically it implies that lowering the superconducting gap enhances the impurity coupling, which may motivate a search for systems with optimal size of superconducting gap given the experimental energy resolution.

The general form of the critical magnetic coupling J j c reveals a quantity that is independent on the impurity parameters J, K and instead characterizes the DOS itself. Namely, the slope of the E S (J) curve at the parityswitching point only depends on the DOS and ∆, dE S dJ J=Jc = -∆f j (∆, 0). (10) This is also the case when the chemical potential µ lies on a [-∆, ∆] interval around the singularity. As shown in App. A, a small displacement of the chemical potential yields a µ-dependent slope but its relative change as µ varies along the aforementioned interval is of the order of 1%. In Fig. 2, we illustrate the unchanging nature of the slope which could be measured in experiments by varying the coupling constant across the transition as it was recently reported [51,59,60]. It could be possible to obtain an estimate of the renormalization parameter, and subsequently to extract information about the nature of the DOS divergence.

B. Higher order Van Hove singularity

In the case of a power-law divergent DOS at the Fermi level as in Eq. ( 5), the self-consistent expression for the Shiba energy reads,

E S = ±∆ 1 -2f 2,ν K + (f 2 1,ν + f 2 2,ν )(K 2 -J 2 ) 4f 2 1,ν J 2 + (1 -2f 2,ν K + (f 2 1,ν + f 2 2,ν )(K 2 -J 2 )) 2 , ( 11
)
where the renormalization parameters

f 1,ν (∆, E S ) = π 2E c 1 + ν cos π 2 ν ∆ 2 + E 2 S E c ν , (12a)
f 2,ν (∆, E S ) = π 2E c 1 -η 1 + η 1 + ν sin π 2 ν ∆ 2 + E 2 S E c ν , (12b)
also inherit the structure of the DOS. The seemingly more complicated appearance of the Shiba-energy equation stems from the particle-hole asymmetry of the DOS. Indeed, if the DOS is symmetric around the Fermi level, i.e. η = 1, the factor f 2,ν vanishes and the expression for the Shiba energy simplifies to Eq. ( 6) where now J = f 1,ν (∆, E S )J and K = f 1,ν (∆, E S )K. Since -1 < ν < 0, for a sufficiently large energy cutoff E c we have f 1,ν > f vh > f const and we obtain the general result:

The effect of a particle-hole symmetric singularity (either a VHS or a HOVHS) in the DOS at the Fermi level is to enhance the effective coupling between the impurity and the electrons in the host.

The critical magnetic coupling for the Shiba transition in case of HOVHS takes a slightly more complicated form,

J hovh,ν c = K 2 + 1 -2f 2,ν (∆, 0)K f 1,ν (∆, 0) 2 + f 2,ν (∆, 0) 2 . ( 13
)
For the particle-hole symmetric case, η = 1, this expression has the same form as for the VHS, Eq. ( 9). Therefore, regarding the effect of the non-magnetic scattering K for the particle-hole symmetric cases of HOVHS we conclude the same about two possible regimes (small vs. large K) as discussed for the VHS after Eq. ( 9). For the particle-hole asymmetric case, η = 1, we note that for reasonable values of the non-magnetic scattering potential, i.e. K E c , the divergence in the DOS always reduces the critical coupling J c with respect to its value for a constant DOS. As shown in Fig. 1 (a), in the small K regime, the stronger the divergence of the DOS, the smaller J c . We note that the f 2,ν (∆, 0)K term stemming from the DOS asymmetry would only become relevant if

|K| 1 f2,ν (∆,0) ∼ E 1+ν c
, i.e. in the large K regime, but then the K 2 term prevails, therefore J c is governed by the non-magnetic scattering amplitude K and it does not depend on the details of the underlying DOS [Fig. 1 (b)]. We remark that when the Fermi level is tuned within an interval of the order of ∆ around the Van Hove singularity, particle-hole symmetric DOSs (η = 1) also lead to a Shiba energy and critical J expressions of the form (11) and ( 13) respectively, where the f 2,ν factor stems from the asymmetry induced by the displacement of the chem- 6) and (11), as function of the magnetic coupling strength J zoomed around the Shiba transition for logarithmicallydivergent (red), power-law divergent with ν = -1/3, η = 2 (orange), and power-law divergent with ν = -1/2, η = 1 (blue) DOS. The tangent, dashed lines indicate the slope at the crossing point which remains constant for increasing values of K denoted in the color bars and which becomes steeper for a stronger divergence. Note that if the DOS is asymmetric, the minimum of Jc is not at K = 0 [cf. blue and orange curves, and Eq. ( 13) with and without vanishing f2,ν (∆, 0) factor]. ical potential. Nevertheless, the main features presented here remain the same (see App. A for details).

Interestingly, we again find a quantity that is independent of the impurity parameters, namely the slope of the E S (J) curve at the Shiba transition point:

dE S dJ J=Jc = -∆ f 1,ν (∆, 0) 2 + f 2,ν (∆, 0) 2 f 1,ν (∆, 0) , (14)
implying also that the particle-hole asymmetry of the HOVHS divergence would be reflected in Shiba states (see Fig. 2).

IV. CONCLUSIONS

In summary, we found that tuning the Fermi level to a conventional VHS or HOVHS enhances the coupling between the superconducting substrate and the magnetic impurity, thereby offering a pathway to improve YSR states engineering. These results become particularly significant in the advent of twisted transition metal dichalcogenides and graphene heterostructures which allow to tune the system over large regions of the parameter space. From a practical standpoint one might inquire about the validity of the results in the vicinity of the Van Hove singularity. We explore this matter in detail in Appendix A. The key takeaway is that our main conclusions remain valid even when the Fermi level lies on an interval of the order of the smallest energy scale in the problem, i.e. µ ∈ [-∆, ∆]. If the chemical potential deviates from the Van Hove singularity by an energy of the order of the superconducting gap, essentially the same electronic states are mixed when superconductivity is switched on, hence a similar behavior of the Shiba states is expectable. The robustness against perturbations in the chemical potential strengthens the experimental relevance of these results; nevertheless, it is fair to remark here we only considered an s-wave superconductor while electronic instabilities due to a divergent DOS may yield more exotic pairings [61] or even destroy superconductivity [62,63]. Moreover, in this work we treated the spin impurity as a classical degree of freedom. For the conventional case of constant density of states, this approach can be theoretically supported by three arguments: by invoking the limit of a large impurity spin in the Kondo model, by the finding that the Kondo Hamiltonian with strong on-site anisotropy typical of impurity adsorption on a surface yields a sub-gap structure similar to that obtained in the classical limit [64], or by the expectation that typical measurement temperatures are sufficiently higher than the Kondo temperature yet sufficiently lower than the superconducting transition temperature so that the quantum effects of impurity are not pronounced. In our case of a diverging density of states, we expect that the first two arguments still hold, however, theory predicts an enhancement of the Kondo temperature [65] which might lead to a reduction of the range of experimental validity of the Shiba model compared to the case of constant density of states. An interesting open challenge is to consider the effect of (HO)VHS on a quantum impurity, in a Kondo model or an even more general Anderson impurity model.

Additionally, we showed that it is possible to extract an impurity-independent quantity from the E S -J curves, namely the slope of the curve at the Shiba transition. State-of-the-art experiments have shown that it is possible to continuously tune the exchange coupling constant between magnetic molecules [51,60] or magnetic adatoms [59] to the substrate by varying the distance between the impurity and the STM tip, thereby unambiguously identifying the transition point. We propose that this technique could be employed to compare the strength of the divergence in the DOS of different compounds regardless of the nature of the impurities, and even provide an estimate of the divergence law. In this Appendix we asses the robustness of the previous results to a deviation in the chemical potential away from the Van Hove singularity of the order of the superconducting gap ∆. To that purpose, we assume that the expressions (3) and ( 5) introduced in the main text remain a good description of the system's DOS and replace the energy dispersion ξ k in the Hamiltonian (1) by ξ kµ, where µ indicates the chemical potential. To be consistent with the previous assumptions, E c must stay the largest energy scale of the problem, therefore we can only address the situations where the chemical potential lies on a comparably small range, i.e. µ ∈ [-∆, ∆].

We find that our main conclusions remain valid even when the Fermi level is not precisely tuned to the Van Hove singularity. The enhancement of the effective coupling constants with respect to the constant DOS scenario continues to exist, we can still access a "large K" regime where the non-magnetic scattering potential controls the critical J, and the slope of the E S (J) curves at the Shiba transition remains independent of the impurity parameters. Further, we note that the relative change of the meaningful observables, namely J c and the slope at the transition point, are of the order of few percent for our range of µ.

In the following we present the results for the two broad classes of DOS singularities.

Conventional Van Hove singularity

The Shiba energy now fulfills a self-consistent equation analogous to that discussed in Section III B in the context of asymmetric HOVHS:

E S = ±∆ 1 -2f 2,vh K + (f 2 1,vh + f 2 2,vh )(K 2 -J 2 ) 4f 2 1,vh J 2 + (1 -2f 2,vh K + (f 2 1,vh + f 2 2,vh )(K 2 -J 2 )) 2 , ( A1
)
where the f 2,vh factor stems from the asymmetry induced by the perturbation in the chemical potential. We recall that when the DOS is symmetric around the Fermi level f 2 = 0 for any class of Van Hove singularities and the Shiba energy equation simplifies to the standard form [Eq. ( 6)]. The renormalizing factors f 1,vh and f 2,vh are obtained by evaluating integrals (B10a) and (B10b), and while it is possible to obtain a closed form in terms of the dilogarithm, the full expressions are too cumbersome to be presented here. Instead, we present the numerical solution of the self-consistent Shiba equation (A1) in Fig. 3, which clearly shows the smallness of the effect. Nevertheless, in the limit µ ω 1, i.e. for a small perturbation in the chemical potential in the strong-coupling limit where ω 2 = ∆ 2 -E 2 S ∼ ∆ 2 , the renormalizing factors have the following approximate compact form (see App. B for details):

f 1,vh (∆, E S , µ) ≈ π 2E c log E c µ 2 + ω 2 , (A2a)
f 2,vh (∆, E S , µ) ≈ π 2E c µ ω . (A2b)
These expressions justify the insignificant variation of the E S (J) curves: the factor f 1,vh inherits the logarithmic form from the DOS as in the µ = 0 case while the correction characterizing the newly introduced asymmetry is such that f 2,vh /f 1,vh ∼ µ ∆ 1 log(Ec/∆)

1.

The approximate factors become most accurate at the Shiba transition where E S strictly vanishes, therefore by expanding them to lowest order in µ ∆ it is possible to capture the behavior of the critical magnetic coupling:

J c J c (µ = 0) = 1 -sgn(K)|α vh | µ ∆ + |β vh | µ ∆ 2 , (A3)
where

|α vh | ≡ |α vh (K, ∆, E c )| < 1 | log(∆/Ec)| and |β vh | ≡ |β vh (K, ∆, E c )|, with |β(K = 0)| ∼ 1 | log(∆/Ec)| .
Thus, it becomes clear that the relative corrections to J c are small assuming a realistic bandwidth, i.e. ∆ E c . Interestingly, in the absence of a non-magnetic potential scattering (K = 0) the system is particle-hole symmetric, and therefore electron or hole doping has the same effect on J c [inset (a) in Fig. 3]. In a likewise manner, the slope of the E S (J) curve at the transition point is independent of the impurity parameters, and hence an even function of µ as well [inset (b) in Fig. 3]:

Slp c ≡ dE S dJ J=Jc = -∆ f 1,vh (∆, 0, µ) 2 + f 2,vh (∆, 0, µ) 2 f 1,vh (∆, 0, µ) , Slp c Slp c (µ = 0) ≈ 1 -|γ vh | µ ∆ 2 , (A4)
where

|γ vh | ≡ |γ vh (∆, E c )| ∼ 1 | log(∆/Ec)| .

Higher order Van Hove singularity

Analogously, if the Fermi level is tuned to the vicinity of a HOVHS, the self-consistent Shiba energy equation also takes the form of Eq. (11), where the renormalization parameters now read

f 1,ν (∆, E S , µ) = π 2E c 1 + ν cos π 2 ν µ 2 + ω 2 E c ν sin[(π -ϕ)ν] + η sin(ϕν) (1 + η) sin π 2 ν , (A5a)
f 2,ν (∆, E S , µ) = π 2E c 1 + ν sin π 2 ν µ 2 + ω 2 E c ν µ 2 + ω 2 ω sin[(π -ϕ)(ν + 1)] -η sin[ϕ(ν + 1)] (1 + η) cos π 2 ν + µ ω f 1,ν (∆, E S , µ), (A5b)
with ϕ = arg(µ + iω). We find that the renormalizing parameters also inherit the power-law dependence from the DOS [cf. µ = 0 case, Eq. ( 12)], with µ entering as a correction to ω as it was the case in the conventional VHS. Therefore, the energy curves when the Fermi level is tuned to the vicinity of a HOVHS maintain the same structure that their µ = 0 counterparts [see Fig. The energy of the in-gap Shiba states E S is obtained by evaluating the Fourier transform of (B1) at r = 0 and solving

det 1 - dk (2π) 2 E S + ξ k τ z + ∆τ x E 2 S -ξ 2 k -∆ 2 (Kτ z ∓ J) = 0. ( B2
) We solve Eq. (B2) by performing the substitution dk (2π) 2 = ρ(ξ)dξ and evaluating

I 1 (0) = - Ec -Ec dξ ρ(ξ) 1 ξ 2 + ω 2 , (B3a)
I 2 (0) = - Ec -Ec dξ ρ(ξ) ξ ξ 2 + ω 2 , (B3b)
where

ω 2 = ∆ 2 -E 2 S .
Note that the same Eq. (B2) also follows from calculating the poles of the transfer matrix in the usual T -matrix method.

a. Logarithmic DOS (conventional Van Hove singularity)

Here

ρ(ξ) = 1 2E c log E c ξ . (B4)
In the limit ω Ec → 0, Here

I 1 (0) = 1 E 2 c 1 0 dx log(x) x 2 + (ω/E c ) 2 ∼ 1 E 2 c ∞ 0 dx log(x) x 2 + (ω/E c ) 2 = - 1 2E c π ω log E c ω , (
ρ(ξ) = 1 1 + η ν + 1 E 1+ν c |ξ| ν η if ξ < 0, 1 if ξ > 0. (B7)
with -1 < ν < 0, ν ∈ Q . To calculate the integrals we use result 3.251.10 from Ref. [START_REF] Gradshteyn | Table of integrals, series, and products[END_REF]:

1 0 dx x p-1 (1-x q ) -p q = π q 1 sin p q π , q > p > 0. (B8)

Taking an infinite energy cutoff, E c → ∞, we have

I 1 (0) = - ν + 1 E 1+ν c ω ν-1 ∞ 0 dx x ν x 2 + 1 = - ν + 1 2E 1+ν c ω ν-1 1 0 dt t -1+ν 2 (1 -t) ν-1 2 = -π ν + 1 2E 1+ν c cos π 2 ν ω ν-1 , (B9a)
I 2 (0) = - ν + 1 E 1+ν c 1 -η 1 + η ω ν ∞ 0 dx x ν+1 x 2 + 1 = ν + 1 2E 1+ν c η -1 1 + η ω ν 1 0 dt t -1-ν 2 (1 -t) ν 2 = π 1 -η 1 + η ν + 1 2E 1+ν c sin π 2 ν ω ν , (B9b)
where we performed the changes of variables x = ξ ω and t = 1

x 2 +1 , and we identified q = 1 and p = 1-ν 2 (B9a) and p = -ν 2 (B9b).

Fermi level tuned to an interval [-∆, ∆] around the Van Hove singularity

We now allow for a finite chemical potential and the relevant integrals read

I µ 1 (0) = - Ec -Ec dξ ρ(ξ) 1 (ξ -µ) 2 + ω 2 , (B10a)
I µ 2 (0) = - Ec -Ec dξ ρ(ξ) ξ -µ (ξ -µ) 2 + ω 2 , (B10b)
where

ω 2 = ∆ 2 -E 2 S .
a. Logarithmic DOS (conventional Van Hove singularity Substituting Eq. (B4) and in the limit ω Ec → 0,

I µ 1 (0) ∼ 1 2E 2 c ∞ 0 dx log(x) [x -(µ/E c )] 2 + (ω/E c ) 2 + log(x) [x + (µ/E c )] 2 + (ω/E c ) 2 = -π 2E c ω log E c µ 2 + ω 2 , (B11) I µ 2 (0) ∼ - 1 2E c ∞ 0 dx x log(x) 1 [x + (µ/E c )] 2 + (ω/E c ) 2 - 1 [x + (µ/E c )] 2 + (ω/E c ) 2 -µ I µ 1 (0) ∼ π 2E c µ ω , (B12)
where in integral (B12) we further took the limit µ ω and approximated the term in brackets by µ Ec 4x (x 2 +(ω/Ec) 2 ) 2 .

b. Power-law DOS (higher order Van Hove singularity) Substituting Eq. (B7), the integrals (B10) are readily solved using the standard keyhole contour of infinite radius with a branch cut on the positive real axis: where Re iϕ = µ + iω.

I µ 1 (0) ∼ -π(1 + ν) R ν E 1+ν

FIG. 2 .

 2 FIG.2. Positive branch of the Shiba energy, Eqs.(6) and(11), as function of the magnetic coupling strength J zoomed around the Shiba transition for logarithmicallydivergent (red), power-law divergent with ν = -1/3, η = 2 (orange), and power-law divergent with ν = -1/2, η = 1 (blue) DOS. The tangent, dashed lines indicate the slope at the crossing point which remains constant for increasing values of K denoted in the color bars and which becomes steeper for a stronger divergence. Note that if the DOS is asymmetric, the minimum of Jc is not at K = 0 [cf. blue and orange curves, and Eq. (13) with and without vanishing f2,ν (∆, 0) factor].

  Appendix A: Perturbation in the chemical potential away from the Van Hove singularity

FIG. 3 .

 3 FIG. 3. Variation of the positive branch of the Shiba energy, Eq. (A1) as a function of the magnetic coupling strength when the Fermi level tuned to a [-∆, ∆] interval around a conventional VHS. The red (green) area corresponds to a vanishing (large) non-magnetic scattering potential K. The solid line indicates the µ = 0 curve while the dashed-dotted lines correspond to µ = ±∆. Inset (a) shows the normalized critical J as a function of the chemical potential for different values of K. The markers indicate the corresponding curves in the main plot. Inset (b) depicts the normalized slope of the ES(J) curve at the transition point as a function the chemical potential. In both insets the solid line represents the exact numerical solution and the dashed line the series expansion to lowest order in µ/∆ [Eq. (A3)]. Plot parameters: Ec = 1000∆.

FIG. 4 .

 4 FIG. 4. Variation of the positive branch of the Shiba energy as a function of the magnetic coupling strength when the Fermi level tuned to a [-∆, ∆] interval around different HOVHS as indicated in the legend. The solid lines indicate the µ = 0 curves and the dashed-dotted lines correspond to µ = ±∆ (omitted in the bottom panel to improve readability). Top panel: Zero non-magnetic scattering potential. Bottom panel: Large non-magnetic scattering potential. Note that an asymmetric DOS (η = 1, in orange) leads to different curves for ±K but whose slope at ES = 0 is the same. Insets show the normalized critical J as a function of the chemical potential for different HOVHS. The relative change in Jc at K = 0 [inset (b)] is one order magnitude larger than in the other considered scenarios because the actual values of Jc are close to 0. The markers indicate the corresponding curves in the main plot. The dashed lines in (b) correspond to the series expansion to lowest order [Eq. (A6)] which was omitted in (d) to improve readability. Note that the ±K curves in the inset are symmetric upon electron or hole-doping only if the DOS is symmetric around the singularity (η = 1). Further, relative change of the absolute value of Jc is symmetric around µ = 0 only if the DOS is symmetric around the singularity. Plot parameters: Ec = 1000∆.

2 +

 2 (ω/E c ) 2 = 0. (B6)b. Power-law DOS (higher order Van Hove singularity)

  c ω sin[(πϕ)ν] + η sin(ϕν) (1 + η) sin(πν) , πϕ)(ν + 1)]η sin[ϕ(ν + 1)] (1 + η) sin[π(ν + 1)] -

ACKNOWLEDGMENTS

We would like to thank Tristan Cren for fruitful discussions.

series expansion in powers of µ ∆ captures the behavior of the critical J at small doping. We find

where

As discussed in the previous Section, in the absence of non-magnetic scattering potential we obtain an even function of µ, which interestingly only depends on the power-law exponent to lowest order. Further, for a finite K the linear term does not vanish, and its sign depends on K in the same manner as in the conventional VHS. The slope of the E S (J) curve at the Shiba transition remains independent of the impurity parameters [compare ±K curves at the same µ in Fig. 4 (c)] and it takes the form of Eq. ( 14) with the newly-introduced renormalization parameters [Eq. (A5)]. At low doping,

where

The slope at the transition is only sensitive to the bulk parameters, therefore, if the DOS is symmetric around the singularity (η = 1), we obtain an even function in µ as discussed in the previous Subsection (see Fig. 5). To calculate the energy of the in-gap YSR states E S we introduce the spinors φ + = (ψ ↑ , ψ † ↓ ) T and φ -= (ψ ↓ , -ψ † ↑ ) T [27]. We choose the impurity to be at the origin, r imp = 0, and we write the eigenvalue equation restricted to particle-hole space, (ξ k τ z + ∆τ x -E S )φ ± (k) = -(Kτ z ∓ Jσ z )φ ± (0). (B1)