Kernel well-posedness and computation by power series in backstepping output feedback for radially-dependent reaction–diffusion PDEs on multidimensional balls

Systems & Control Letters 177 (2023) 105538 Contents lists available at ScienceDirect Systems & Control Letters journal homepage: www.elsevier.com/locate/sysconle Kernel well-posedness and computation by power series in backstepping output feedback for radially-dependent reaction–diffusion PDEs on multidimensional balls✩ Rafael Vazquez a,∗, Jing Zhang b, Jie Qi b, Miroslav Krstic c aDepartamento de Ingenieria Aeroespacial y Mecanica de Fluidos, Universidad de Sevilla, Camino de los Descubrimiento s.n., 41092 Sevilla, Spain bCollege of Information Science and Technology, Donghua University, Shanghai 201620, China cDepartment of Mechanical and Aerospace Engineering, University of California San Diego, La Jolla, CA 92093-0411, USA article info Article history: Received 10 July 2022 Received in revised form 20 February 2023 Accepted 20 April 2023 Available online 12 May 2023 Keywords: Partial differential equations Spherical harmonics Infinite-dimensional systems Backstepping Parabolic systems abstract Recently, the problem of boundary stabilization and estimation for unstable linear constant-coefficient reaction–diffusion equation on n-balls (in particular, disks and spheres) has been solved by means of the backstepping method. However, the extension of this result to spatially-varying coefficients is far from trivial. Some early success has been achieved under simplifying conditions, such as radiallyvarying reaction coefficients under revolution symmetry, on a disk or a sphere. These particular cases notwithstanding, the problem remains open. The main issue is that the equations become singular in the radius; when applying the backstepping method, the same type of singularity appears in the kernel equations. Traditionally, well-posedness of these equations has been proved by transforming them into integral equations and then applying the method of successive approximations. In this case, with the resulting integral equation becoming singular, successive approximations do not easily apply. This paper takes a different route and directly addresses the kernel equations via a power series approach (in the spirit of the method of Frobenius for ordinary differential equations), finding in the process the required conditions for the radially-varying reaction (namely, analyticity and evenness) and showing the existence and convergence of the series solution. This approach provides a direct numerical method that can be readily applied, despite singularities, to both control and observer boundary design problems. ©2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Prologue It is with admiration that we contribute this article to celebrate the 80th birthday of Professor Art Krener in a special issue in Systems & Control Letters—a journal in which he has published some of his most impactful papers, such as his 1983 introduction of nonlinear observers, jointly with Isidori, and his 1984 singleauthored introduction of approximate feedback linearization. Art’s contributions are not only deep but also remarkably broad, including nonlinear controllability/observability, nonlinear ✩This work was supported in part by the National Natural Science Foundation of China (62173084), the Fundamental Research Funds for the Central Universities, China and Graduate Student Innovation Fund of Donghua University, China (CUSF-DH-D-2019089) and the scholarship from China Scholarship Council (CSC201806630010). Rafael Vazquez acknowledges financial support of the Spanish Ministerio de Ciencia, Innovación y Universidades, Spain under grant PGC2018-100680-B-C21. ∗Corresponding author. E-mail addresses: [email protected] (R. Vazquez), [email protected] (J. Zhang), [email protected] (J. Qi), [email protected] (M. Krstic). representation theory, feedback linearization, observer design, control of bifurcations, optimal/bang–bang control, nonlinear H-infinity control, solutions to Hamilton–Jacobi PDEs, and estimation and analysis of stochastic systems. His career in control theory started with his work on the maximum principle for problems where the first order conditions are inconclusive. In order to generate conclusive results with the help of higher order terms, he turned to Lie algebra and differential geometry tools and employed his extension of Chow’s theorem whose answer to the questions of accessibility and reachable states leads to higher order necessary conditions for optimality. Art then transitioned to representation theory, to answer when a control-affine state–space model can be represented as an input–output model via Volterra series and, somewhat conversely, when a nonlinear input–output map can be represented as a bilinear state–space model. His work was a significant inspiration for our spatial Volterra series approach to control of nonlinear PDEs [1]. Art’s most impactful contribution resulted when he turned his attention to controllability in his 1977 IEEE Transactions on Automatic Control paper with mathematical physicist Hermann. https://doi.org/10.1016/j.sysconle.2023.105538 0167-6911/©2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 Krener provided a comprehensive set of results on controllability, observability, and minimality for nonlinear systems, tackling them analogously to the linear case—using rank conditions and leading to decompositions into controllable and uncontrollable subsystems. The initial idea for the celebrated feedback linearization originated in Krener’s necessary and sufficient conditions for the existence of a change of state variables that yields a transformed system that is linear. Brockett, Hunt and Su, and Jacubczyk and Respondek extended this idea by adding feedback to the state transformation, producing the general framework now known as feedback linearization. (One could argue that ‘‘PDE backstepping’’ [2] is a PDE realization of the feedback linearization idea, with a PDE target system extending the finite-dimensional Brunovsky form.) Recognizing himself that feedback linearizability is not satisfied for some important physical systems, Krener then went on to develop a theory of ‘‘approximate feedback linearization’’, published in Systems & Control Letters [3]. Pursuing a dual of the feedback linearization approach, for the problem of state estimation, Krener developed an approach to the design of observers for nonlinear systems, in collaboration with Isidori in Systems & Control Letters [4], and later with Respondek, which employs a state transformation and the injection of measured output. What we do in Section 6of the present paper is a PDE analog of Krener’s approach. Art pursued nonlinear control for systems with disturbances in several waves—first in a paper with Isidori, Gori-Giorgi and Monaco in an award-winning paper that gives existence conditions and constructions of decoupling and noninteracting control laws, and then several years later in the context of nonlinear H-infinity control and estimation. In the early and mid-1990s, inspired by ‘‘bifurcation control’’, wildly popular in the physics community at that time, and linked to interesting applications, Krener and Kang developed control designs for normal forms that include locally quadratic dependencies, which remain the definitive results of the art on this subject. Art Krener cast an eye on PDE control systems at least as early as the mid-1990s, in the framework of a US Air Force funded project on nonlinear control of jet engine instabilities. While other researchers focused on the first-order Galerkin approximations of the models of rotating stall instabilities in axial flow compressors in jet engines, Art formulated and studied generalized higher-order Moore–Greitzer nonlinear models [5]. It is a delight to see Art Krener take on PDE control as the preoccupation for the present stage of his career [6,7]. Going beyond the conventional development of operator Riccati equations, and the limitation to the study of their well posedness, in his trademark fashion, Art is producing computable approximate solutions to Riccati PDEs using Al’brekht’s approach, which he has already brought to its state-of-the art form for nonlinear and stochastic ODEs. 2. Introduction In this paper we introduce an explicit boundary outputfeedback control law to stabilize an unstable linear radiallydependent reaction–diffusion equation on an n-ball (which in 2-D is a disk and in 3-D a sphere). This paper extends the spherical harmonics [8] approach of [9], which assumed constant coefficients, using some of the ideas of [10]. For a finite number of harmonics, we design boundary feedback laws and output injection gains using the backstepping method [2] (with kernels computed using a power series approach) which allows us to obtain exponential stability of the origin in the L2norm. Higher harmonics will be naturally open-loop stable. The required conditions for the radially-varying coefficients are found in the analysis of the numerical method and are non-obvious (evenness of the reaction coefficient). The idea of using a power series to compute backstepping kernels was first seen in [11] (without much analysis of the method itself, but rather numerically optimizing the approximation) and later in [12], where piecewise-smooth kernels require the use of several series. Here, we prove that the method provides a unique converging solution, in the spirit of the method of Frobenius for ordinary differential equations. Some partial results towards the solution of this problem were obtained in [13] and [14] for the disk and sphere, respectively; however they required symmetry conditions. Older results in this spirit were obtained in [15] and [16]. This paper extends and completes our conference contribution [17] where the ideas where initially presented (without proof). To the best of our knowledge, this paper presents the first rigorous proof of convergence of a power series solution for the backstepping kernel equations. Thus, this work consolidates the method as a valid alternative to more traditional numerical approaches, which include finite difference approximations of the kernel equations [2,18–20], the use of symbolical successive approximation series [21], or the numerical solution of the integral version of the kernel equations [22,23]. The main advantages of the method are its simplicity (it does not require the sometimes cumbersome conversion to integral equations thus preventing mistakes or any consideration about discrete meshes), speed (modern computing systems can reach high orders of the series in seconds), precision (one reaches a simple polynomial in one variable for the gain at the boundary that does not require interpolation), adaptability (it can be adapted to settings with discontinuous kernels by breaking the domain in pieces, see [12]), and capacity to produce kernels depending on parameters (by symbolically solving the kernel equations). The main drawback is the analyticity requirement of the system coefficients, even though most physical systems and examples seen in backstepping papers indeed possess analytic coefficients, and possibly a slow convergence rate of the series in some cases. Previous results and applications in multi-dimensional domains include multi-agent deployment in 3-D space [24] (by combining the ideas of [9] and [25]), convection problem on annular domains [21], PDEs with boundary conditions governed by lower-dimensional PDEs [10,26], multi-dimensional cuboid domains [27]. The backstepping method has proved itself to be an ubiquitous method for PDE control, with many other applications including, among others, flow control [28,29], nonlinear PDEs [1], hyperbolic 1-D systems [30–32], or delays [33]. Nevertheless, other design methods are also applicable to the geometry considered in this paper (see for instance [34] or [35]). The structure of the paper is as follows. In Section 3we introduce the problem. In Section 4we state our stability result. We study the well-posedness of the kernels in Section 5, which is the main result of the paper, proving existence of the kernels and providing means for their computation; interestingly, odd and even dimensions require a slightly different approach. We briefly talk next about the observer in Section 6, but skip most details based on its duality with respect to the controller. Then, we give some simulation results in Section 7. We finally conclude the paper with some remarks in Section 8. 3. n-D reaction–diffusion system on an n-ball Consider the following constant-coefficient reaction–diffusion system in an n-dimensional ball of radius R: ∂u ∂t=ϵ(∂2u ∂x2 1+∂2u ∂x2 2+. . . +∂2u ∂x2 n)+λu=ϵ△nu+λ( x)u,(1) 2 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 where u=u(t, x), with  x= [x1,x2,...,xn]T, is the state variable, evolving for t>0 in the n-ball Bn(R) defined as Bn(R)={ x∈Rn: ∥ x∥ ≤ R},(2) with boundary conditions on the boundary of Bn(R), which is the (n−1)-sphere Sn−1(R) defined as Sn−1(R)={ x∈Rn: ∥ x∥ = R}.(3) The boundary condition is assumed to be of Dirichlet type, u(t, x)⏐⏐⏐ x∈Sn−1(R)=U(t, x),(4) where U(t, x) is the actuation. On the other hand the measurement y(t, x) is defined as y(t, x)=∂ru(t, x)⏐⏐⏐ x∈Sn−1(R),(5) where ∂rdenotes the derivative in the radial direction (normal to the (n−1)-sphere), which would be defined as ∂ru(t, x)= ∇u· x ∥ x∥. Differently from [9], we consider non-constant λ( x) verifying the following assumption. Assumption 3.1. The coefficient λ( x) is an analytic function of  x and depends exclusively on the radius r= ∥ x∥. Following [9], both the state and the actuation variable can be written in n-dimensional spherical coordinates, also known as ultraspherical coordinates (see [8], p. 93), which consist of one radial coordinate rand n−1 angular coordinates  θ. Then, using a (complex-valued) Fourier–Laplace series of Spherical Harmonics1 to handle the angular dependencies, defined as u(t,r, θ)= l=∞ ∑ l=0 m=N(l,n) ∑ m=0 um l(r,t)Yn lm( θ),(6) U(t, θ)= l=∞ ∑ l=0 m=N(l,n) ∑ m=0 Um l(t)Yn lm( θ),(7) where N(l,n) is the number of (linearly independent) n-dimensional spherical harmonics of degree l, given by N(0,n)= 1 (representing the mean value over the n-ball) and, for l>0, N(l,n)=2l+n−2 l(l+n−3 l−1),(8) with Yn lm being the mth n-dimensional spherical harmonic of degree l. The coefficients in (6)–(7) are possibly complex-valued. Following [9] and using (6)–(7) one reaches the following independent complex-valued 1-D reaction–diffusion equation for each spherical harmonic coefficient: ∂tum l=ϵ rn−1∂r(rn−1∂rum l)−l(l+n−2) ϵ r2um l+λ(r)um l,(9) evolving in r∈(0,R],t>0, with boundary conditions um l(t,R)=Um l(t),(10) In these equations, we have considered Dirichlet boundary conditions. The measurement would be the flux at the boundary, namely ∂rum l(t,R). Note that following [36, p. 640], a second boundary condition, reflecting the second-order character of (9) and the need to avoid singular behaviors, can be expressed as: |um l(t,0)|<∞.(11) 1Spherical harmonics were introduced by Laplace to solve the homonymous equation and have been widely used since, particularly in geodesics, electromagnetism and computer graphics. A very complete treatment on the subject can be found in [8]. In the above equations, the integers mand lstand for the order and degree of the harmonic, respectively. Note that the higher the degree (corresponding to high frequencies), the more ‘‘naturally’’ s Eqs. (9)–(10) is, as seen next. Define the L2norm ∥f∥L2=√∫R 0|f(r)|2rn−1dr.(12) and the associated L2space as usual, where |f|2=ff ∗, being f∗ the complex conjugate of f. Lemma 3.1. Given λ(r)and R, there exists L ∈Nsuch that, for all l >L, the equilibrium um l≡0of system (9)–(10) is open loop exponentially stable, namely, for Um l=0in (10) there exists a positive constant D1, such that for all t ∥um l(t,·)∥L2≤e−D1t∥um l(0,·)∥L2.(13) D1is independent of l, and only depends on n, λ(r),ϵ, and R, and can be chosen as large as desired just by increasing the values of L. The proof is skipped as it mimics [10] just by using the L2norm as a Lyapunov function and Poincare’s inequality. Thus one only needs to stabilize the unstable mode with l<L. Since the different modes are not coupled, it allows us to stabilize them separately and re-assembling them. Moreover since only a finite number of harmonics is stabilized, there is no need to worry about the convergence of the control law as in [9], with its Spherical Harmonics series being just a finite sum. Our objective can now be stated as follows. Considering only the unstable modes, design an output-feedback control law for Um l using, for each mode, only the measurement of ∂rum l(t,R). Our design procedure is established in the next section along with our main stability result. 4. Stability of controlled harmonics Next, for the unstable modes we design the output-feedback law. The observer and controller are designed separately using the backstepping method, by following [9]; in this reference it is shown that both the feedback and the output injection gains can be found by solving a certain kernel PDE equation, which is essentially the same for both the controller and the observer. Thus, for the sake of brevity and to avoid repetitive material, we only show how to obtain the (full-state) control law, giving the basic observer design and some additional remarks later in Section 6. 4.1. Design of a full-state feedback control law for unstable modes Based on the backstepping method [2], our idea is utilizing an invertible Volterra integral transformation wm l(t,r)=um l(t,r)−∫r 0 Kn lm(r, ρ)um l(t, ρ)dρ, (14) where the kernel Kn lm(r, ρ) is to be determined, which defined on the domain Tk= {(r, ρ)∈R2;0≤ρ≤r≤R}to convert the unstable system (9)–(10) into an exponentially target system: ∂twm l=ϵ∂r(rn−1∂rwm l) rn−1−ϵl(l+n−2)wm l r2−cwm l,(15) wm l(t,R)=0,(16) where the constant c>0 is an adjustable convergence rate. From (14) and (16), let r=R, we obtain the boundary control as the following full-state law Um l(t)=∫R 0 Kn lm(R, ρ)um l(t, ρ)dρ. (17) 3 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 Following closely the steps of [9] to find conditions for the kernels, and defining Kn lm(r, ρ)=Gn lm(r, ρ)ρ(ρ r)l+n−2, we finally reach a PDE that the G-kernels need to verify: λ(ρ)+c ϵGn lm =∂rr Gn lm+(3−n−2l)∂rGn lm r−∂ρρ Gn lm+(1−n−2l)∂ρGn lm ρ, (18) with only one boundary condition: Gn lm(r,r)= −∫r 0(λ(σ)+c)dσ 2rϵ.(19) We assume as usual that these kernel equations are well-posed and the resulting kernel is bounded in T; this will be analyzed later in Section 5, providing also a numerical method for its computation. 4.2. Closed-loop stability analysis of unstable modes To obtain the stability result of closed-loop system, we need three elements. We begin by stating the stability result for the target system. We follow by obtaining the existence of an inverse transformation that allows us to recover our original variable from the transformed variable. Then we relate the L2norm with spherical harmonics. With these elements, we construct the proof of stability mapping the result for the target system to the original system. This is done by showing that the transformation is an invertible map from L2into L2. We first discuss the stability of the target system, having the following lemma: Lemma 4.1. For all l ∈N, and for c ≥0, the equilibrium wm l≡0of system (15)–(16) is exponentially stable, i.e., there exists a positive constant D2such that for all t, ∥wm l(t,·)∥L2≤e−D2t∥wm l(0,·)∥L2,(20) where the constant D2is independent of n, l or m, and only depends on c, ϵ, and R; it can be chosen as large as desired just by increasing the value of c. Proof. Consider the Lyapunov function: V2(t)=1 2∥wm l(t,·)∥2 L2,(21) then, taking its time derivative, we obtain ˙ V2=∫R 0 ¯wm l∂twm l+wm l∂t¯wm l 2rn−1dr ≤ −(ϵ 4R2+c)∥wm l∥2 L2 (22) choosing c=D2−ϵ 4R2(23) we then obtain, independent of the value of n, ˙ V2≤ −2D2V2,(24) thus proving the result. □ Lemma 4.2. For |l| ≤ L, let c be chosen as in Lemma 4.1, and assume that the kernel Kn lm(r, ρ)is bounded and integrable. The system (9) with boundary control (17) is closed-loop exponentially stable, namely there exists positive constants C and D2such that ∥um l(t,·)∥L2≤Ce−D2t∥um l(0,·)∥L2,(25) C and D2are independent of m or l, and only depend on n, L, λ(r), ϵand R. Proof. The proof consists of two parts, one is existence of an inverse transformation, and then showing the equivalence of norms of the variables un lm and wn lm; the result then follows from the stability of the target system. As shown in [9], when Kn(r, ρ) is bounded and integrable, the map (14) is invertible and its inverse transformation is um l(t,r)=wm l(t,r)+∫r 0 Ln lm(r, ρ)wm l(t, ρ)dρ, (26) which is also bounded and integrable. Call now ¯ Kand ¯ Lthe maximum of the bounds of the function ˇ Kn lm and ˇ Ln lm for a given n and all l≤Lin their respective domains. It is easy to get ∥wm l(t,·)∥2 L2≤M1∥um l(t,·)∥2 L2,(27) ∥um l(t,·)∥2 L2≤M2∥wm l(t,·)∥2 L2.(28) where M1=2+R4¯ K/(2n) and M2=2+R4¯ L/(2n). Combining then Lemma 4.1 with the norm equivalence between um land wm l system stated as in (27) and (28), it is easy to obtain ∥um l(t,·)∥L2≤√M2∥wm l(t,·)∥L2≤√M2e−D2t∥wm l(0,·)∥L2 ≤√M1M2e−D2t∥um l(0,·)∥L2.(29) Let C=√M1M2, the result then follows. □ Note that combining Lemmas 3.1 and 4.2 and taking D= min{D1,D2}, we get the following stability result for all spherical harmonics and thus the full physical system. Theorem 1. Under the assumption that the kernel Kn lm(r, ρ)is bounded and integrable, the equilibrium um l≡0of system (9)–(10) under control law (17) is closed-loop exponentially stable, namely, there exists a positive constant D, such that for all t ∥um l(t,·)∥L2≤Ce−Dt ∥um l(0,·)∥L2.(30) where D can be chosen as large as desired just by increasing the value of L and c in the control design process. 5. Well-posedness of the kernel equations Next, we state the main result of the paper, which was in part assumed in Theorem 1, also giving the requirements for λ(r). In addition the proof of the result also provides a numerical method to compute the kernels, which is an alternative to successive approximations which do not work in this case (due to the singularities at the origin; see for instance [13] to see the resulting singular integral equation that needs to be solved). Theorem 2. Under the assumption that λ(r)is an even real analytic function in [0,R], then for a given n >1and all values of l ∈N, there is a unique power series solution Gn lm(r, ρ)for (18)–(19), even in its two variables in the domain T, which is real analytic in the domain. In addition, if λ(r)is analytic, but not even, then there is no power series solution to (18)–(19) for most values l ∈N. The requirement of evenness for λ(r) might seem unusual. However, if we carefully consider Assumption 3.1, and since r= ∥ x∥ = √x2 1+x2 2+. . . +x2 n, in physical space λ( x) will be nonanalytic, unless it is even. Thus, while solutions to the kernel equations might exist for non-even λ(r), we cannot expect them to be analytic. This result notwithstanding, if one is interested in controlling only very low-order harmonics, kernels do exist without this requirement, as shown in [13,14], which only consider the 0-th order harmonic (the mean) respectively for a disk and a sphere, and only require boundedness of λ(r). 4 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 5.1. Proof of Theorem 2 We start by giving out an algorithmic method to compute the power series for Gn lm(r, ρ), which will allow us to prove Theorem 2 as well as numerically approximating the kernels. First of all, we show that the evenness of λ(r) is a necessary condition to find an analytic solution. Next, it is possible to establish that the series for Gn lm(r, ρ) only has even powers. Exploiting this property to suitably express (18)–(19), we finally show the existence of the power series and thus Theorem 2 follows. Convergence and related issues (radius of convergence) is studied towards the end, finishing the proof. 5.1.1. Computing a power series solution for the kernels Starting from the most basic assumption of Theorem 2, we consider that λ(r) is analytic in [0,R], therefore it can be written as a convergent series (encompassing cand ϵfor notational convenience): λ(r)+c ϵ=∞ ∑ i=0 λiri,(31) which, by the evenness of λ, may only contain even powers,2this is, λi=0 if iis odd. We then, in the spirit of the method of Frobenius for ordinary differential equations, seek for a solution of (18)–(19) of the form: Gn lm(r, ρ)=∞ ∑ i=0⎛ ⎝ i ∑ j=0 Cijrjρi−j⎞ ⎠,(32) where the dependence on n,land mhas been omitted for simplicity (the solution will depend on these values). The series in (32) collects together (in the parenthesis) all the polynomial terms with the same degree. It is easy to see that the boundary condition (19) implies: ∀i, i ∑ j=0 Cij = − λi 2(i+1),(33) which in particular implies C00 = −λ0 2ϵ. On the other hand, the left-hand side of (18) becomes λ(ρ)+c ϵGn lm =⎡ ⎣ ∞ ∑ i=0⎛ ⎝ i ∑ j=0 Cijrjρi−j⎞ ⎠⎤ ⎦ ∞ ∑ i=0 λiρi =∞ ∑ i=0⎛ ⎝ i ∑ j=0 Bijrjρi−j⎞ ⎠,(34) where we have defined Bij = i ∑ k=j Ckjλi−k.(35) Finally, to express the right-hand side of (18), denote γ=n+ 2l−2≥0 and define the operators D1=∂rr +(1 −γ)1/r∂rand D2= −∂ρρ +(−1−γ)1/ρ∂ρ. Then D1Gn lm =∞ ∑ i=1⎛ ⎝ i ∑ j=1 j(j−γ)Cijrj−2ρi−j⎞ ⎠,(36) 2It is a known fact of analysis that even functions (respectively, odd functions) contain only even powers (respectively, odd powers) in their Taylor series. This fact has an easy proof by substituting the series in the definition of evenness f(r)=f(−r) (respectively, oddness f(r)= −f(−r)) and checking the conditions verified by the coefficients. D2Gn lm =∞ ∑ i=1⎛ ⎝ i−1 ∑ j=0 (i−j)(j−i−γ)Cijrjρi−j−2⎞ ⎠,(37) and thus, rewriting the sum to be homogeneous with (34),we find (D1+D2)Gn lm =∞ ∑ i=−1⎛ ⎝ i+1 ∑ j=−1 Dijrjρi−j⎞ ⎠,(38) where, (assuming Cij =0 if i,j<0 or j>i), Dij =(j+2)(j+2−γ)C(i+2)(j+2) −(i−j+2)(i−j+2+γ)C(i+2)j,(39) Equating (38) and (34), we obtain a system of equations: ∀i≥ −1,Di(i+1) =Di(−1) =0,(40) ∀i≥0,0≤j≤i,(j+2)(j+2−γ)C(i+2)(j+2) −(i−j+2) (i−j+2+γ)C(i+2)j=Bij.(41) With λ(r) and nare fixed, we want to show that the kernel equations are solvable for all values of l∈N. Thus, γtakes increasing values. In addition we can assume γ= 1, since the case n=3, l=0 was already addressed in [14] showing that it reduces to the usual 1-D kernel equations for parabolic systems [2], which admits a power series solution according to [11]. The first two equalities, if γ= 1, imply ∀i≥1,Ci1=Ci(i−1) =0,(42) and, in particular, C10 =C01 =0, whereas the second equality results in a system of equations that needs to be solved recursively, starting at i=0. It can be rewritten as follows to start at i=2 (since C00,C10 and C01 are already determined). ∀i≥2,0≤j≤i−2, (j+2)(j+2−γ)Ci(j+2) −(i−j)(j−i−γ)Cij = i−2 ∑ k=j Ckjλi−2−k. (43) Note that for each i≥2, there are i+1 coefficients in (32) but i+2 relations: one from (33), two from (42) and i−1 from (43). Thus, it would seem that (33)–(42)–(43) is in general an incompatible system. This is indeed the case if λ(r) is not even, i.e., if the series (31) contains odd powers, as shown in the next section. 5.1.2. Evenness requirement of λ(r) We start with the following result. Lemma 5.1. If λ(r)is not even, then there are values of l ∈Nfor which there is no solution to (18)–(19) in the form of (32). Proof. We show that, if there exists iodd such that λi= 0, then there is no solution in the form of a power series. First, if λ1= 0, then from (33) we know that C01 +C10 = −λ1 4ϵ, however since form (42) one has C01 =C10 =0, this cannot hold. Consider now there is indeed a value i>1 for which a coefficient λiis distinct from zero and let us show the result by contradiction. Consider the first such i. Now, since the right-hand side of (43) depends on C(i−2)j, one gets that for all odd i′<i Ci′jmust zero from (33)– (42)–(43) all having a zero right-hand side (this can be formalized with an induction argument; we skip the details). Thus, at i, the following system of equations has to be verified: Ci1=Ci(i−1) =0,(44) i ∑ j=0 Cij = − λi 2ϵ(i+1),(45) 5 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 and for 0 ≤j≤i−2, (j+2)(j+2−γ)Ci(j+2) −(i−j)(i−j+γ)Cij =0,(46) Let us consider lsufficiently large such that γ > i, so that the coefficient (j+2−γ)in(46) is distinct from zero in the full range of j, namely 0 ≤j≤i−2. Then none of the coefficients in (46) is zero. Therefore, combining (44) with (46), from Ci1we can find Ci3, then Ci5, and so on. Similarly, from Ci(i−1) we can find Ci(i−3), Ci(i−5) and so on. These two sequences do not overlap because iis odd and therefore, one finds Cij =0 for all 0 ≤j≤iwhich is not compatible with (45) unless λi=0, which contradicts our initial assumption. □ Next we show that evenness of λimplies evenness of the kernels. Lemma 5.2. If λ(r)is even, then, a solution to (18)–(19) in the form of (32) only has even powers. Proof. We need to prove that Cij =0 if either ior jis odd. From the proof of Lemma 5.1, we directly know that for odd i one has Cij =0. Fix, then, ieven and consider jodd; for i=2, the result is obvious. Assuming Ci′j=0 for all even numbers i′<iand jodd, let us prove the result by induction on the first coefficient. As before, we would need to solve (45)–(43). The right-hand side Bij =∑i−2 k=jCkjλi−2−kof (43) is zero as in (46) by the induction hypothesis (if keven) or directly zero if kodd. Then, following again the proof of Lemma 5.1, we have the same system of Eqs. (45)–(46) for our even iand odd j’s. Now: Cij =(j+2)(j+2−γ) (i−j)(i−j+γ)Ci(j+2), so starting from Ci(i−1) =0 we find Ci(i−3) =0, then Ci(i−5), and so on; however, with ibeing even, this sequence ends now in Ci1(thus, the proof of Lemma 5.1 does not apply because the sequences starting at Ci1and Ci(i−1) overlap). Thus, one finds Cij = 0 for all odd values of jbetween 1 and i−1. □ 5.1.3. Well-posedness of the coefficient system Next, we show that the coefficients of the power series can always be found, which by the previous lemmas only requires studying the even coefficients. For simplification, we redefine (31) and (32) as: λ(r)+c ϵ=∞ ∑ i=0 λir2i,Gn lm(r, ρ)=∞ ∑ i=0⎛ ⎝ i ∑ j=0 Cijr2jρ2(i−j)⎞ ⎠,(47) without bothering to redefine the coefficients (note that (35) does not require any change). Defining as well γ′=γ 2=n 2+l−1≥0, the new system of equations to be solved is ∀i, i ∑ j=0 Cij = − λi 2(2i+1),(48) and ∀i≥1,0≤j≤i−1,(j+1)(j+1−γ′)Ci(j+1) +(i−j)(j−i−γ′)Cij = i−1 ∑ k=j Ckjλi−1−k=B(i−1)j.(49) Let us outline the solution procedure, and later derive some conclusions. Solving in (49) every Cij as a function of Ci(j+1) we get: Cij =(j+1)(j+1−γ′)Ci(j+1) +B(i−1)j (i−j)(i−j+γ′),(50) which can be written more briefly if we define, for i>0 and 0≤j<i, aij(γ′)=(j+1)(j+1−γ′) (i−j)(i−j+γ′),(51) as Cij =aij(γ′)Ci(j+1) +B(i−1)j (i−j)(i−j+γ′).(52) To be able to simplify a bit the equation, redefine ˆ B(i−1)j=B(i−1)j (i−j)(i−j+γ′)(53) then, Cij =aij(γ′)Ci(j+1) +ˆ B(i−1)j.(54) and iterating this equality until reaching Cii, we get Cij =⎡ ⎣ k=i−1 ∏ k=j aik(γ′)⎤ ⎦Cii +ˆ B(i−1)j+ i−1 ∑ r=j+1 k=i−1 ∏ k=r aik(γ′)ˆ B(i−1)r, and inserting this into (48), we reach an equation for Cii, namely Cii = − 1 κ(i, γ ′)[λi 2ϵ(2i+1) +Hi],(55) where κ(i, γ ′)=1+ i−1 ∑ j=0 k=i−1 ∏ k=j aik(γ′),(56) and Hi= i−1 ∑ j=0ˆ B(i−1)j+ i−2 ∑ j=0 i−1 ∑ r=j+1 k=i−1 ∏ k=r aik(γ′)ˆ B(i−1)r = i−1 ∑ j=0ˆ B(i−1)j+ i−1 ∑ r=1 r−1 ∑ j=0 k=i−1 ∏ k=r aik(γ′)ˆ B(i−1)r = i−1 ∑ j=0ˆ B(i−1)j+ i−1 ∑ r=1 r k=i−1 ∏ k=r aik(γ′)ˆ B(i−1)r = i−1 ∑ j=0⎛ ⎝1+j k=i−1 ∏ k=j aik(γ′)⎞ ⎠ˆ B(i−1)j(57) It is quite clear that these κ(i, γ ′) will play an important role; in particular, if they are non-zero, one can always find a unique solution for the coefficients Cij. Thus one needs to show that κ(i, γ ′)= 0 for any possible ior γ′. The following lemma shows this is indeed the case, by exploiting a connection of the aij coefficients with Gauss’ hypergeometric functions. Lemma 5.3. Let i be a positive integer and γ′≥0a real number. Then, it holds that κ(i, γ ′)=2i! i! Γ(γ′+1) Γ(i+γ′+1) >0,(58) where Γdenotes the Gamma function [37, p. 255]. Proof. Recalling from (56) and (51) the definitions of κ(i, γ ′) and aij, respectively, one has κ(i, γ ′)=1+ i−1 ∑ j=0 k=i−1 ∏ k=j (k+1)(k+1−γ′) (i−k)(i−k+γ′)(59) 6 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 which can be rewritten in terms of binomial numbers and rising/falling factorials3[38] as κ(i, γ ′)= i ∑ j=0(i j)(i−γ′)i−j (1 +γ′)i−j(60) and reordering the sum and using (i j)=(i i−j), κ(i, γ ′)= i ∑ j=0(i j)(i−γ′)j (1 +γ′)j= i ∑ j=0 (−1)j(i j)(γ′−i)j (1 +γ′)j,(61) where the obvious fact (x)j=(−1)j(−x)jhas been used. Consider now the finite polynomial pi(x, γ ′) defined as pi(x, γ )= i ∑ j=0 (−1)j(i j)(γ′−i)j (1 +γ′)jxj(62) From the definition of Gauss’ hypergeometric function [37, p. 561], denoted as 2F1(a,b;c;x), in the polynomial case (aor bnonpositive integer) and noting (−1)j(i j)=(−i)j j!, it is immediate that pi(x, γ )=2F1(−i, γ ′−i;1+γ′;x) (63) and therefore, from Gauss’ summation theorem [37, p. 556], which is applicable in this case since 1 +2i>0, κ(i, γ ′)=pi(1, γ )=2F1(−i, γ ′−i;1+γ′;1) =Γ(1 +γ′) Γ(1 +γ′+i) Γ(1 +2i) Γ(1 +i)=2i! i! Γ(1 +γ′) Γ(1 +γ′+i),(64) finishing the proof. □ The next result is an immediate conclusion of the positivity of κ(i, γ ′): Lemma 5.4. If λ(r)is even, then, for all values of l ∈R, the coefficients in (47) that solve (18)–(19) can be uniquely found up to any order i. To conclude the proof of Theorem 2, we need to prove analyticity of the series (47). This step, however, requires splitting the problem in two possible cases: odd dimension (thus, γ′= n/2+2l−1 is not an integer) and even dimension (γ′integer). 5.1.4. Proof of analyticity for odd dimension In the odd-dimension case, define the following coefficients: Li0=1, Lij =(i j)(i+γ′)(i−1+γ′). . . (i−j+γ′+1) (1 −γ′)(2 −γ′). . . (j−γ′),j>0 (65) with Li0defined as 1; these are well-defined given that γ′is non-integer. They can also be expressed as Lij =(i j)Γ(1 −γ′)Γ(i+1+γ′) Γ(j+1−γ′)Γ(i−j+1+γ′) Now, in (48)–(49), denote Cij =Lij ˇ Cij. Replacing in the recurrence we get B(i−1)j=(j+1)(j+1−γ′)(i j+1) 3Rising factorials (x)nare sometimes expressed using the Pochhammer’s symbol, with a slightly different notation, namely (x)n. Γ(1 −γ′)Γ(i+1+γ′) Γ(j+2−γ′)Γ(i−j+γ′)ˇ Ci(j+1) −(i−j)(i−j+γ′)(i j) Γ(1 −γ′)Γ(i+1+γ′) Γ(j+1−γ′)Γ(i−j+1+γ′)ˇ Cij =(i−j)(i−j+γ′)Lij (ˇ Ci(j+1) −ˇ Cij)(66) Define now ˇ B(i−1)j=B(i−1)j (i−j)(i−j+γ′)Lij (67) and the new set of recurrence equations for ˇ Cij becomes rather simple: ∀i, i ∑ j=0 Lij ˇ Cij = − λi 2(2i+1),(68) ∀i≥1,0≤j≤i−1,ˇ Ci(j+1) −ˇ Cij =ˇ B(i−1)j,(69) and the recurrence is easily solvable in terms of one element; for instance, ˇ Cii: ˇ Cij =ˇ Cii − r=i−1 ∑ r=j ˇ B(i−1)r(70) for j=1,...,i. Replacing in (68) we reach i ∑ j=0 Lij ˇ Cii − i−1 ∑ j=0 Lij r=i−1 ∑ r=j ˇ B(i−1)r= − λi 2(2i+1) Thus: ˇ Cii =−λi 2(2i+1) +∑k=i−1 k=0∑r=i−1 r=kLik ˇ B(i−1)r ∑k=i k=0Lik =−λi 2(2i+1) +∑r=i−1 r=0ˇ B(i−1)r(∑k=r k=0Lik) ∑k=i k=0Lik (71) Call Rij = k=j ∑ k=0 Lik,Ri= k=i ∑ k=0 Lik =Rii.(72) Then ˇ Cii =−λi 2(2i+1) +∑r=i−1 r=0ˇ B(i−1)rRir Ri (73) Now, solving for the remaining coefficients from (70): ˇ Cij =−λi 2(2i+1) +∑r=i−1 r=0ˇ B(i−1)rRir Ri− r=i−1 ∑ r=j ˇ B(i−1)r =−λi 2(2i+1) +∑r=i−1 r=0ˇ B(i−1)rRir −∑r=i−1 r=jˇ B(i−1)rRi Ri =−λi 2(2i+1) +∑r=j−1 r=0ˇ B(i−1)rRir −∑r=i−1 r=jˇ B(i−1)r(Ri−Rir ) Ri (74) Finally, recovering the coefficients Cij from Cij =Lij ˇ Cij and using (67): Cij = −λi 2(2i+1) Lij +∑r=j−1 r=0B(i−1)r LijRir (i−r)(i−r+γ′)Lir −∑r=i−1 r=jB(i−1)r Lij(Ri−Rir ) (i−r)(i−r+γ′)Lir Ri , (75) which is quite explicit. 7 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 Notice that since ρ≤r, |Gn lm(r, ρ)| ≤ ∞ ∑ i=0⎛ ⎝ i ∑ j=0|Cij|r2jρ2(i−j)⎞ ⎠≤∞ ∑ i=0 r2i⎛ ⎝ i ∑ j=0|Cij|⎞ ⎠, thus, defining αi=∑i j=0|Cij|, if we can prove that ∑∞ i=0αir2i converges for a certain radius of convergence R, so does Gn lm(r, ρ) for ρ≤r≤R, and thus we obtain the required analyticity. Now: αi= i ∑ j=0|Cij| ≤|λi| 2(2i+1) ∑i j=0|Lij| |Ri|+∑i j=0∑r=j−1 r=0|B(i−1)r|| LijRir (i−r)(i−r+γ′)Lir | |Ri| +∑i j=0∑r=i−1 r=j|B(i−1)r|| Lij(Ri−Rir ) (i−r)(i−r+γ′)Lir | |Ri| =|λi| 2(2i+1) ∑i j=0|Lij| |Ri|+∑i−1 r=0∑i j=r+1|B(i−1)r|| LijRir (i−r)(i−r+γ′)Lir | |Ri| +∑r=i−1 r=0∑r j=0|B(i−1)r|| Lij(Ri−Rir ) (i−r)(i−r+γ′)Lir | |Ri| ≤|λi| 2(2i+1) ∑i j=0|Lij| |Ri| + i−1 ∑ r=0|B(i−1)r|(|Rir |∑i j=r+1|Lij|+|Ri−Rir |∑r j=0|Lij| (i−r)(i−r+γ′)|Lir ||Ri|),(76) To prove the convergence of the power series ∑∞ i=0αir2iconsider the following lemma, inspired by [39]: Lemma 5.5. Consider g(x)=∑∞ i=0gix2iand h(x)=∑∞ i=0hix2i analytic functions, both with radius of convergence R. Let i0be a nonnegative integer, let (ai)∞ i=0be a sequence of real numbers, and define f (x)=∑∞ i=0aix2i, where aiverify, for i >i0 ai≤bi|gi|+ci i−1 ∑ j=0 aj|hi−1−j|(77) where the sequences bi,ci≥0are decreasing for i >i0, with ci also verifying limi→∞ ci=0. Then, f (x)is analytic with radius of convergence at least R. Proof. Since gand hanalytic with radius of convergence Rwe can write |gi|,|hi| ≤ MR−2i, where the definition as power series of squares has been taken into account. Thus: ai≤biMR−2i+ci i−1 ∑ j=0 ajMR−2i+1+2j(78) Define ˇ ai=aifor i≤i0and, for i>i0,ˇ ai=biMR−2i+ ci∑i−1 j=0ˇ ajMR−2i+1+2j. Obviously ai≤ˇ aiand therefore the radius of convergence of f(x) would be at least the radius of convergence of ˇ f(x)=∑∞ i=0ˇ aix2i. Now: ˇ ai+1=bi+1MR−2i−2+ci+1 i ∑ j=0ˇ ajMR−2i+2j−1 =bi+1MR−2i−2+ci+1Mˇ ai+ci+1R−2 i−1 ∑ j=0ˇ ajMR−2i2+j+1(79) It is sufficient to compute lim i→∞ ˇ ai+1 ˇ ai =lim i→∞ bi+1MR−2i−2+ci+1Mˇ ai+ci+1R−2∑i−1 j=0ˇ ajMR−2i+2j+1 ˇ ai =lim i→∞Mci+1 +lim i→∞ bi+1MR−2i−2+ci+1R−2∑i−1 j=0ˇ ajMR−2i+2j+1 biMR−2i+ci∑i−1 j=0ˇ ajMR−2i+1+2j =R−2lim i→∞ bi+1+ci+1∑i−1 j=0ˇ ajR2j+1 bi+ci∑i−1 j=0ˇ ajR1+2j ≤R−2lim i→∞ bi+ci∑i−1 j=0ˇ ajR2j+1 bi+ci∑i−1 j=0ˇ ajR1+2j =R−2,(80) where the inequality holds for sufficiently large i>i0and thus in the limit, therefore proving the lemma. □ To apply Lemma 5.5 to (76) we need to bound some of the terms. In particular, if we are able to find biand cisuch that ∑i j=0|Lij| 2(2i+1)|Ri|≤bi, max r∈{0,...,i−1}(|Rir |∑i j=r+1|Lij|+|Ri−Rir |∑r j=0|Lij| (i−r)(i−r+γ′)|Lir ||Ri|)≤ci,(81) we get αi≤bi|λi|+ci i−1 ∑ r=0|B(i−1)r| =bi|λi|+ci i−1 ∑ r=0 i−1 ∑ k=r|Ckr ||λi−1−k| =bi|λi|+ci i−1 ∑ k=0 k ∑ r=0|Ckr ||λi−1−k| =bi|λi|+ci i−1 ∑ k=0 αk|λi−1−k|.(82) Thus, assuming that biand civerify the conditions given in Lemma 5.5, and given that λ(x) has a radius of convergence of at least R, we see that Gn lm(r, ρ) converges and defines an analytic function for ρ≤r≤R, thus proving Theorem 2 for the odd-dimension case. It remains to find such biand ci. Proceeding exactly as in Lemma 5.3 with a slight modification, we directly find Ri=2F1(−i,−(γ′+i);1−γ′;1) =Γ(1 −γ′)Γ(2i+1) Γ(−γ′+1+i)Γ(i+1) =2i! i! Γ(1 −γ′) Γ(1 −γ′+i)(83) Now, let N=γ′−1/2 and i>2N. One can see that for 1≤j≤N, |Lij| = (i j)(i+γ′)(i−1+γ′). . . (i−j+γ′+1) (γ′−1)(γ′−2) . . . (γ′−j) ≤ |LiN | =(i N)(i+γ′)(i−1+γ′). . . (i+3/2) (γ′−1)(γ′−2) . . . (1/2) ≤i! N!(i−N)! (i+N+1)! (i+1)! 1 (N−1)!(84) and |Ri| = 2i! i!⏐⏐⏐⏐ Γ(γ′−1) Γ(i+1−γ′)⏐⏐⏐⏐≥2i! i!(i−N)!N!.(85) 8 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 Thus, for i>2N, calling dithe following sequence di=∑N j=0|Lij| |Ri|≤i!(i+N+1)! (i+1)(N−1)!2i!,(86) it is clear that diis a decreasing sequence, since from the ratio test r=limi→∞ di+1 di=1/4. Now, set i0=2N. For i>i0, ∑i j=0|Lij| 2ϵ(2i+1)|Ri|=∑N j=0|Lij|+∑i j=N+1|Lij| 2(2i+1)|Ri|≤di+1 2(2i+1) =bi(87) It is obvious that biis decreasing, since diis decreasing. Now we need to find a sequence cifor the second term in (81). First of all, |Rir |∑i j=r+1|Lij|+|Ri−Rir |∑r j=0|Lij| (i−r)(i−r+γ′)|Lir ||Ri| ≤2∑i j=r+1|Lij|∑r j=0|Lij| (i−r)(i−r+γ′)|Lir ||Ri|(88) The following lemmas are needed to find a bound to (88). Lemma 5.6. Let N =γ′−1/2and i >2N+1. Then, define j∗= ⌊i−1+N 2⌋. It holds that |Lij|≤|Lij∗|. Proof. Consider the ratio |Lij+1| |Lij|. It is easy to see that |Lij+1| |Lij|=(i−j+γ′)(i−j) |j+1−γ′|(j+1) (89) Now, if j≤N<i/2−1/2, then |j+1−γ′| = γ′−j−1> γ′−1=N+1/2. Then, we have |Lij+1| |Lij|=(i−j+γ′)(i−j) (γ′−j−1)(j+1) >i/2(i/2+γ′) (N+1/2)(N+1) ≥1.(90) Thus, the sequence always increases as long as j≤N, and we can look for a maximum j∗>N. Then, for j>N, denote the ratio of (89) by f: f=|Lij+1| |Lij|=(i−j+γ′)(i−j) (j+1−γ′)(j+1) (91) Now, f≤1 implies (i−j)2+γ′(i−j)≤(j+1)2−γ′(j+1). Thus, (i−j)2−(j+1)2+γ′(i+1) ≤0. Manipulating the expression, we find (i2−1) −2j(i+1) +γ′(i+1) ≤0 and canceling the term (i+1) the following inequality for jis reached: j≤i−1+γ′ 2(92) Therefore, if (and only if) the bound given by (92) is verified, |Lij+1| |Lij|≤1. Therefore, we conclude that the maximum of the sequence |Lij|is reached at j=j∗= ⌊i−1+γ′ 2⌋=⌊i−1+N 2+1/4⌋=⌊i−1+N 2⌋(93) thus finishing the proof. □ Lemma 5.7. Let N =γ′−1/2and i >2N+1. Then we have ∑i j=r+1|Lij|∑r j=0|Lij| (i−r)(i−r+γ′)|Lir ≤2Lij∗, where j∗is defined in (93). Proof. Now, to bound the term ∑i j=r+1|Lij|∑r j=0|Lij, consider two possibilities and use Lemma 5.6. If r<j∗, then ∑i j=r+1|Lij| ∑r j=0|Lij ≤(i−r)(r+1)|Lij∗||Lir |. On the other hand, if r≥j∗, then ∑i j=r+1|Lij|∑r j=0|Lij ≤(i−r)(r+1)|Lir+1||Lij∗|. Therefore, if r<j∗: ∑i j=r+1|Lij|∑r j=0|Lij| (i−r)(i−r+γ′)|Lir |≤(r+1)|Lij∗| i−r+γ′≤(j∗+1)|Lij∗| i−j∗+γ′ else, if r≥j∗, ∑i j=r+1|Lij|∑r j=0|Lij| (i−r)(i−r+γ′)|Lir |≤(r+1)|Lir+1||Lij∗| (i−r+γ′)|Lij| and using (91), ∑i j=r+1|Lij|∑r j=0|Lij| (i−r)(i−r+γ′)|Lir |≤(i−r)|Lij∗| (r+1−γ′)≤(i−j∗)|Lij∗| (j∗+1−γ′) Now, since i−1+N 2≤j∗≤i+N 2, one has that j∗+1 i−j∗+γ′≤ i+N 2+1 i−i+N 2+N+1/2=i+N+2 i+N+1<2, and similarly, i−j∗ j∗+1−γ′≤i−i−1+N 2 i−1+N 2+1/2−N=i+1−N i−N<2, thus concluding the proof. □ Thus, we are left with showing that |Lij∗| |Ri|is decreasing, which is expected, since Lij∗is one of the elements that appear in the sum Ri. Using the expression (83) and the formula for Lij that involves the Gamma function, we obtain the following: |Lij∗| |Ri|=i!2Γ(1 −γ′)Γ(i+1+γ′)Γ(i+1−γ′) 2i!j∗!(i−j∗)!Γ(γ′−1)Γ(j∗+1−γ′)Γ(i−j∗+1+γ′). (94) Now, the decreasing character of the sequence is established as follows (consider the case where i+Nis odd, so that j∗=i−1+N 2; the even case is analogous). Consider Stirling’s approximation to the factorial, namely n! ≈ √2πn(n e)n. Then: i!2 2i!j∗!(i−j∗)!=i!2 2i!(i−1+N 2)!(i+1−N 2)! ≈√i π(i2−(N−1)2) 1 22i ei (i−1+N 2)i−1+N 2(i+1−N 2)i+1−N 2 .(95) On the other hand, Stirling’s approximation to the Gamma function [37, p. 257] reads Γ(z)≈√2π z(z e)z, thus Γ(i+1−γ′)Γ(i+1+γ′) Γ(j∗+1−γ′)Γ(i−j∗+1+γ′) =Γ(i+1/2−N)Γ(i+3/2+N) Γ(i−N 2)Γ(i+N 2+2) ≈1 ei√(i−N 2)(i+N 2+2) (i+1/2−N)(i+3/2+N) ×(i+1/2−N)i+1/2−N(i+3/2+N)i+3/2+N (i−N 2)i−N 2(i+N 2+2)i+N 2+2.(96) Putting together (95)–(96), we obtain the following. |Lij∗| |Ri|≈iΓ(1 −γ′) √πΓ(γ′−1)√f1(i)f2(i)if3(i),(97) where we have broken the approximation into three functions: f1(i)=i π(i2−(N−1)2) (i−N 2)(i+N 2+2) (i+1/2−N)(i+3/2+N),(98) 9 R. Vazquez, J. Zhang, J. Qi et al. Systems & Control Letters 177 (2023) 105538 on Control Applications, 1997, pp. 651–656. [6] A.J. Krener, Optimal boundary control of a nonlinear reaction diffusion equation via completing the square and Al’brekht’s method, IEEE Trans. Automat. Control 67 (12) (2021) 6698–6709. [7] A.J. Krener, Boundary control of the beam equation by linear quadratic regulation, Systems Control Lett. 153 (2021) 104949. [8] K. Atkinson, W. Han, Spherical Harmonics and Approximations on the Unit Sphere: An Introduction, Springer, 2012. [9] R. Vazquez, M. Krstic, Boundary control of reaction–diffusion PDEs on balls in spaces of arbitrary dimensions, ESAIM:Control Optim. Calculus Variations 22 (4) (2016) 1078–1096. [10] R. Vazquez, M. Krstic, J. Zhang, J. Qi, Stabilization of a 2-D reaction– diffusion equation with a coupled PDE evolving on its boundary, in: Proceedings of the 2019 IEEE Conference on Decision and Control, CDC, 2019. [11] P. Ascencio, A. Astolfi, T. Parisini, Backstepping PDE design: A convex optimization approach, IEEE Trans. Automat. Control 63 (2018) 1943–1958. [12] L. Camacho-Solorio, R. Vazquez, M. Krstic, Boundary observers for coupled diffusion-reaction systems with prescribed convergence rate, Systems Control Lett. 134 (2020) 104586. [13] R. Vazquez, M. Krstic, Boundary control of a singular reaction–diffusion equation on a disk, in: CPDE 2016 (2nd IFAC Workshop on Control of Systems Governed By Partial Differential Equations), 2016. [14] R. Vazquez, M. Krstic, Boundary control and estimation of reaction– diffusion equations on the sphere under revolution symmetry conditions, Internat. J. Control 92 (1) (2019) 2–11. [15] F. Bribiesca Argomedo, C. Prieur, E. Witrant, S. Bremond, A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients, IEEE Trans. Automat. Control 58 (2013) 290–303. [16] S.J. Moura, N.A. Chaturvedi, M. Krstic, PDE estimation techniques for advanced battery management systems—Part I: SOC estimation, in: Proceedings of the 2012 American Control Conference, 2012. [17] R. Vazquez, M. Krstic, J. Zhang, J. Qi, Output feedback control of radiallydependent reaction-diffusion PDEs on balls of arbitrary dimensions, IFAC-PapersOnLine 53 (2) (2020) 7635–7640. [18] S. Kerschbaum, J. Deutscher, Backstepping control of coupled linear parabolic PDEs with space and time dependent coefficients, IEEE Trans. Automat. Control 65 (7) (2019) 3060–3067. [19] J. Deutscher, J. Gabriel, A backstepping approach to output regulation for coupled linear wave-ODE systems, Automatica 123 (2021) 109338. [20] H. Anfinsen, O.M. Aamo, Adaptive Control of Hyperbolic PDEs, Springer, New York, 2019. [21] R. Vazquez, M. Krstic, Boundary observer for output-feedback stabilization of thermal convection loop, IEEE Trans. Control Syst. Technol. 18 (2010) 789–797. [22] L. Jadachowsi, T. Meurer, A. Kugi, An efficient implementation of backstepping observers for time-varying parabolic PDEs, IFAC Proc. Vol. 45 (2) (2012) 798–803. [23] J. Auriol, D. Bresch-Pietri, Robust state-feedback stabilization of an underactuated network of interconnected n+ m hyperbolic PDE systems, Automatica 136 (2022) 110040. [24] J. Qi, R. Vazquez, M. Krstic, Multi-Agent deployment in 3-D via PDE control, IEEE Trans. Automat. Control 60 (4) (2015) 891–906. [25] T. Meurer, M. Krstic, Finite-time multi-agent deployment: A nonlinear PDE motion planning approach, Automatica 47 (2011) 2534–2542. [26] J. Qi, M. Krstic, S. Wang, Stabilization of reaction–diffusion PDE distributed actuation and input delay, in: Proceedings of the 2018 IEEE Conference on Decision and Control, CDC, 2018. [27] T. Meurer, Control of Higher-Dimensional PDEs: Flatness and Backstepping Designs, Springer, 2013. [28] R. Vazquez, M. Krstic, Control of Turbulent and Magnetohydrodynamic Channel Flow, Birkhauser, 2008. [29] R. Vazquez, E. Trelat, J.-M. Coron, Control for fast and stable laminar-to- high-Reynolds-numbers transfer in a 2D Navier-Stokes channel flow, Disc. Cont. Dyn. Syst. Ser. B 10 (2008) 925–956. [30] L. Hu, R. Vazquez, F. Di Meglio, M. Krstic, Boundary exponential stabilization of 1-D inhomogeneous quasilinear hyperbolic systems, SIAM J. Control Optim. 57 (2) (2019) 963–998. [31] J. Auriol, F. Di Meglio, Minimum time control of heterodirectional linear coupled hyperbolic PDEs, Automatica 71 (2016) 300–307. [32] G. Andrade, R. Vazquez, D. Pagano, Backstepping stabilization of a linearized ODE-PDE Rijke tube model, Automatica 96 (2018) 98–109. [33] M. Krstic, Delay Compensation for Nonlinear, Adaptive, and PDE Systems, Birkhauser, 2009. [34] R. Triggiani, Boundary feedback stabilization of parabolic equations, Appl. Math. Optimiz. 6 (1980) 201–220. [35] V. Barbu, Boundary stabilization of equilibrium solutions to parabolic equations, IEEE Trans. Automat. Control 58 (2013) 2416–2420. [36] R. Haberman, Elementary Applied Partial Differential Equations, second ed., Prentice-Hall, 1987. [37] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, ninth ed., Dover, 1965. [38] D.E. Knuth, Two notes on notation, Amer. Math. Monthly 99 (5) (1992) 403–422. [39] V. Leon, B. Scardua, On singular Frobenius for second order linear partial differential equations, 2019, preprint downloaded from arXiv, https://arxiv. org/abs/1907.02620. [40] K.B. Howell, Ordinary Differential Equations: An Introduction to the Fundamentals, second ed., CRC Press, 2019. [41] E. Butkov, Mathematical Physics, Addison-Wesley, 1995. 16