Inverse problems of damped wave equations with Robin boundary conditions: an application to blood perfusion

Knowledge of the properties of biological tissues is essential in monitoring any abnormalities that may be forming and have a major impact on organs malfunctioning. Therefore, these disorders must be detected and treated early to save lives and improve the general health. Within the framework of thermal therapies, e.g. hyperthermia or cryoablation, the knowledge of the tissue temperature and of the blood perfusion rate are of utmost importance. Therefore, motivated by such a significant biomedical application, this paper investigates, for the first time, the uniqueness and stable reconstruction of the space-dependent (heterogeneous) perfusion coefficient in the thermal-wave hyperbolic model of bio-heat transfer from Cauchy boundary data using the powerful technique of Carleman estimates. Additional novelties consist in the consideration of Robin boundary conditions, as well as developing a mathematical analysis that leads to stronger stability estimates valid over a shorter time interval than usually reported in the literature of coefficient identification problems for hyperbolic partial differential equations. Numerically, the inverse coefficient problem is recast as a nonlinear least-squares minimization that is solved using the conjugate gradient method (CGM). Both exact and noisy data are inverted. To achieve stability, the CGM is stopped according to the discrepancy principle. Numerical results for a physical example are presented and discussed, showing the convergence, accuracy and stability of the inversion procedure.


Physical background
Understanding the heat transfer in biological tissues is of crucial importance in biomedical applications such as hyperthermia, thermal ablation or microwave heating, which are typical examples of thermal therapies used to treat cancer, menorrhagia and benign prostate hyperplasia [1]. When mathematically modelling such applications, care should be taken to include the underlying processes that are taking place such as heat conduction, blood perfusion and heat generation due to metabolism. One formulation that takes into account these mechanisms is based on the much celebrated Pennes' parabolic reaction-diffusion equation (obtained by taking τ = 0 in equation (1) below), which was proposed to model the temperature evolution during cancer hyperthermia treatment [51], the thermal radiation from cellular phones [55] and the ablation of afflicted tissues [28], among others. However, although still widely used, the Pennes parabolic model of heat transfer implies infinite speed of heat propagation. This characteristic becomes practically unrealistic when modelling heat propagation in biological tissues for which a non-negligible relaxation time-lag τ (typically between 15-30 seconds) occurs [49]. Therefore, a more appropriate model taking into account that in biological tissues thermal waves propagate with a finite speed is given by the hyperbolic equation [46], where Ω is the spatial domain occupied by the tissue, T is a final time (in s) of interest, u is the tissue temperature (in • C), u b is the (arterial) blood temperature (in • C), κ and C tissue are the thermal conductivity (in W (m • C) −1 ) and heat capacity (in J (kg • C) −1 ) of the tissue, respectively, C b is the heat capacity (in J (kg • C) −1 ) of the blood, w b is the bood perfusion rate (in s −1 ) and Q contains the heat generations (in W m −3 ) due to metabolism and external heating. Equation (1) is supplied with the initial conditions u = Φ 0 , u t = 0 in Ω at t = 0 and the Robin boundary condition where Φ 0 is prescribed initial temperature (in • C), ν is the outward unit normal to the boundary ∂Ω and h is the heat transfer coefficient (in W (m 3 • C) −1 ) between the tissue and the surrounding ambient having a temperature u amb (in • C). In case h → 0, equation (3) becomes the homogeneous Neumann adiabatic boundary condition whilst in case |h| → ∞, equation (3) becomes the Dirichlet boundary condition u = u amb on (0, T) × ∂Ω. (5) In practical applications concerned with hyperthermia treatment or assessing burn injuries, evaluating the blood perfusion rate w b is of main interest. In this paper, we consider the case when the unknown coefficient w b = w b (x) is space-dependent across the tissue and we aim to recontruct it from non-intrusive temperature measurements over the whole boundary ∂Ω or on a portion of it Γ 0 ⊂ ∂Ω, namely, Before performing any analysis, equations (1)- (6) have to be non-dimensionalised. Assuming that u b is a non-zero constant, we introduce the following dimensionless variables: With this non-dimensionalisation, equations (1)- (3) and (6) recast as (denoting F := Q + Q t ) where, without any confusion, the dimensionless ν, Γ 0 , Ω and ∂Ω have been denoted by the same symbols as their dimensional counterparts. Equations (4) and (5) also take the nondimensional form u = u amb on (0, T) × ∂Ω.
Certain thermal experiments [7] may be designed so that the ambient temperature u amb is equal to the arterial blood temperature u b . In such a situation, from (7) it follows that u amb = 0 and equations (10) and (13) simplify, respectively, as ∂ ν u + hu = 0 on (0, T) × ∂Ω. (14) and

Mathematical background
In the mathematical analysis we consider the exposition in a more general setup than the particular physical model in equations (8) and (9). To start with, assume Ω ⊂ R d , d ∈ N * , is an open and bounded C 2 -domain. We mention that a one-dimensional related, but different inverse problem has recently been analysed in [52].
Consider the (homogeneous) Robin problem related to (8), (9) and (14) given by where the functions Q 0 and Q 1 are in L ∞ (Ω). Moreover, assume ϕ 0 ∈ H 1 (Ω), ϕ 1 ∈ L 2 (Ω), F ∈ L 2 (Ω T ) and h = h(x) ∈ L ∞ (∂Ω). In most of the mathematical analysis presented below, Q 0 and Q 1 are independent of each other, but we shall also consider the physical case of equation (1) in which Q 0 (x) = w(x) and Q 1 (x) = 1 + w(x) in theorems 4.4 and 4.5. In the limit that h → 0 we obtain the (homogeneous) Neumann problem in which the Robin boundary condition in (16) is being replaced by (12). This paper is organized as follows. In section 2, the mathematical formulation of the direct problem and its well-posedness are presented. In section 3, the poweful technique of Carleman estimates is discussed and the main theorem 3.3 is established. Section 4 is fully devoted to theoretically investigating the inverse coefficient identification problems of interest and establishing conditional stability estimates. We mention herein that Lipschitz stability estimates, similar to those established in theorems 4.1, 4.3 and 4.4, but for the wave equation c(x)u tt = ∆u with unknown coefficient c(x) were obtained in [36] and [35, theorem 3.6]. In section 5, the inverse problem is recast as a variational problem. The resulting nonlinear least-squares minimizing functional is proved to be Fréchet differentiable and an explicit expression for its gradient is derived. Numerical results obtained using the conjugate gradient method (CGM) for a physical example are presented and discussed in section 6 to confirm the proposed inversion algorithm's convergence, accuracy and stability. Finally, conclusions are provided in section 7.

Comparisons with other works
Our approach to the well-posedness of direct problem is based on Evans' book [17] with adaptions to Robin boundary condition settings and some basic spectral theory of the Laplacian [5,6]. The method of Carleman estimates is inspired by approach developed in the book of Bellassoued and Yamamoto [10], as well as the one of Lerner [44]. Theorem 3.3 extends the Carleman estimate in [10, chapter 4], which can also be viewed as an extension of [44] in the presence of boundary. There are several improvements obtained in this paper when comparing with other literature. One improvement compared with [10] is to include all other boundary contributions into the Carleman estimate (see remark 3.4). This can also be compared with the work of Lasiecka et al in [40], where authors used a different approach to obtain the Carleman estimate and their result looks different from ours. As an application of our Carleman estimate, the stability estimates presented in section 4 are sharper than the works of [10,23,40,47] (see remark 4.2). This mathematical refinement is physically meaningful, as it results in a shorter time of physical measurements on the boundary that is required to determine the solution in the interior of the domain.

Direct problem
The well-posedness of the direct problem (16) is established in the literature, (see, for instance, [37, chapter 4], [17, chapter 7] and [45]). First, we define the time-dependent bilinear form
and Q 1 ∈ L ∞ (Ω) be given input data for the direct problem (16). Then there exists a unique weak solution u ∈ Λ h of the direct problem (16). Moreover, this solution satisfies the stability estimate (18) where the positive constant C depends on the L ∞ -norms of Q 0 and Q 1 .
We defer the proof to appendix A since methods of proving theorem 2.1 can be found in various literature with some minor modifications. For instance, in [37, chapter IV, theorem 5.1], one can find a similar statement for the well-posedness of the Robin problem (16) but with F ∈ L 1 (0, T; L 2 (Ω)) instead of F ∈ L 2 (Ω T ), as in our theorem.
The rest of the paper concerns the analysis of inverse coefficient problems. The main ingedient that is used for proving stability estimates of these inverse problems is based on Carleman estimates, which are introduced and discussed in the next section.

Carleman estimates
A brief review of literature. Applications of Carleman-type estimates for obtaining various observability and controllability results in inverse problems were initiated by the groundbreaking publication of Bukhgeim and Klibanov [12] and have subsequently been tremendously successful in the last few decades [30,31,[33][34][35]57]. For instance, a Carleman-type estimate can be obtained as in [42] for a parabolic setting or as in [43] for a hyperbolic setting. A systematic study of various Carleman estimates without boundary conditions, but under a general pseudo-convexity assumption, can be found in the book by Lerner [44]. As we will be focusing on hyperbolic inverse problems, we will mainly employ tools developed by Bellassoued and Yamamoto in [10], which can also be found in [24,25]. Even though the book [10] focuses on hyperbolic Carleman estimates under Dirichlet boundary conditions, the tool can be extended to cope with Neumann conditions and, more importantly, Robin boundary conditions in our setting modelling the heat transfer with the surrounding ambient environment, as given by equation (10). Note that Carleman estimates with boundary conditions can be traced back to the work of Isakov in [27], Lavrent'ev et al in [41], as well as Tataru's works in [53,54].
Notations. The approach presented in this section can be adapted to more general manifolds. For the sake of simplicity, we will assume that our n-dimensional time-space manifold, denoted by X, is Minkowski, where the metric is given by η = dt 2 − g and g is the Euclidean metric. Therefore, we can identify vectors with covectors. That is, for vector fields V = (V 0 ,Ṽ) and . We will also denote η(V, V) by η(V) and g(Ṽ,Ṽ) by |Ṽ| 2 . For a given function f defined in the time-space, we will use df to denote its exterior derivative, which means it is an n-component-wise vector. We will use ∇f to denote the vector containing the spatial components of df. In other words, ∇f is a d-component-wise vector, where d = n − 1. The Hessian with respect to g of a function, f, is denoted by ∇ 2 f and its action on vectors is given by

Carleman estimate: a control from interior
In this subsection, we will review a derivation of a classical Carleman estimate via the method developed in [22,44]. This will serve us as a guidance in the next subsection, when we extend the method to cope with boundary conditions. Let m ∈ N * , a typical Carleman estimate (inequality) for an m-th order differential operator P acting on R n and a suitable function φ is given by where wherev denotes the Fourier transform of the function v. Note that the main role in (25) is represented by the principal part of the operator P because contributions from lower-order terms can be absorbed into λ∥v∥ 2 H m−1 λ (R n ) . When working with λ-weighted Sobolev spaces and λ-dependent operators P λ , we will treat λ the same as ξ, e.g. the symbol of P λ will be a polynomial in (ξ, λ).
A strategy to prove the celebrated inequality (25) via pseudo-differential operators is explored in great details in [44] and it is summarised below.

and obtain
Cλ∥v∥ 2 which is essentially the estimate (25) for sufficiently large λ.
Remark 3.1. We will briefly explain the reasons behind each step.
(1) The reason of the splitting in equation (26) is to involve the third term, which is one order lower than p 1 in equation (27) form = 0. Moreover, symbolic calculus of pseudodifferential operators can be applied easily to the commutator. (2) Whenm = 0, inequality (28) is the same as the sufficient condition of Carleman estimate proved in theorem 28.2.3 of [22]. (3) The advantage of using the Fefferman-Phong inequality is that (29) holds for principally normal operators, which include operators with complex-valued principal symbols. For operators with real-valued principal symbols, one observes that (29) can be derived via for somem ⩽ 0, which can then be achieved by applying Gårding's inequality to (28) This leads us to introduce the concept of strong pseudo-convexity.

Pseudo-convex functions.
As the pseudo-convexity of a function is important in establishing Carleman estimates, we will now give its definition.

Definition 3.2.
Let X be a n-dimensional manifold and P be an m-th order differential operator with C 1 (X) principal coefficients and L ∞ loc (X) complex-valued lower-order terms. A function φ is called strongly pseudo-convex at ζ 0 ∈ X with respect to P if dφ (ζ 0 ) ̸ = 0 and for all (ξ, λ) ∈ (R n × R + )\{0}, we have that where ξ λ,ϵ = ξ + i(λ + )dφ(ζ 0 ). Moreover, we say that φ is strongly pseudo-convex with respect to P on a set U ⊂ X, if it is strongly pseudo-convex for each ζ 0 ∈ U.
The operator in (16) defined by satisfies the criteria of the operator in definition 3.2 for n = 1 + d. Moreover, P is a secondorder (m = 2) differential operator and its principal symbol is given by are the dual variables of (t, x). As quadratic polynomials can only have complex roots in conjugate pairs, (30) is void for λ ̸ = 0. Therefore, it simplifies the pseudoconvexity criterion for second-order differential operators to and H p2 is the Hamiltonian vector field generated by p 2 .

Carleman estimate with boundary conditions
If we can write a second-order differential operator as P = P 2 + P 1 , where P 1 is a first-order differential operator with L ∞ (Ω T ) coefficients, then ∥P 1 u∥ L 2 (ΩT) can be absorbed by choosing sufficiently large constant C in (25). In other words, for the case of P being (31), it suffices to establish a Carleman estimate for P = P 2 = □. Therefore, we will first establish a Carleman estimate for the conjugate operator of □, i.e. P ϕ,λ = e λϕ □e −λϕ , where φ is a function in such that dψ ̸ = 0 in Ω T and the Hessian of ψ 0 with respect to g is positive in Ω, Here, the small positive parameter β and the large positive parameter γ will be chosen so that Hörmander's pseudo-convexity property (30) is satisfied. It turns out that, for sufficiently large γ and sufficiently small β, a function φ ∈ P is strongly pseudo-convex with respect to P on Ω T , see appendix B for the justification. A sub-class of P is given by for some positive constant δ. We will use this sub-class of functions to obtain the Carleman estimate (35) below. For where Note that equation (34) is an analogue of (26) up to a change of lower-order terms. In the following theorem, we extend the results in [10, chapter 4] to include the lower (t = 0) and upper (t = T) surface contributions. Theorem 3.3. Let u ∈ H 2 (Ω T ) and φ be a function in the set P 0 . Then there exists a positive constant γ 0 (β, δ, ρ), such that, for any γ > γ 0 , there exist positive constants C(α, β, γ, δ, ρ, Ω T ) and λ 0 (γ), for which the following estimate hold: where and The same estimate holds for P being replaced by P + P 1 , where P 1 is a first-order differential operator with L ∞ (Ω T ) coefficients. The only difference is that λ 0 will also be dependent on the L ∞ (Ω T )-norm of the coefficients.
Remark 3.4. In the special case ψ 1 (t) = (t − t 0 ) 2 for some 0 < t 0 < T, then the boundary term derived in [10, chapter 4] is different from B 1 only by the first term in equation (36). The second term B 2 in (37) is a new boundary contribution that is absent under the vanishing conditions of u at t = 0 and t = T, which was imposed in [10, chapter 4].  (32), as in [35]. However, we have stated (and proved) theorem 3.3 for any function φ ∈ P 0 for generality and in order to be consistent with Hörmander's approach to Carleman estimates (see [22]).
Proof of theorem 3.3. We will focus on P + ϕ,λ u, P − ϕ,λ u in (34). As in [10, chapter 4], one has Applying integration by parts to eliminate all the second-order derivatives of u in I k , we have Summing the above expression for I 1 , . . ., I 5 and I 6 , we have where B 0 is the summation of boundary terms in I k and it is given by Now, B 0 becomes Since ∇ 2 ψ 0 (x ′ , x ′ ) > αρ|x ′ | 2 for all x ′ ∈ R d , together with equation (38), this shows that ∥∇u∥ L 2 (ΩT) is controlled by P + ϕ,λ u, P − ϕ,λ u , and we are left to control the time-derivative term of u, i.e. −2βγλ´Ω T φ ψ ′ ′ 1 |u t | 2 d xdt. This means that another positive term alike´Ω T φ|u t | 2 d xdt is needed to overcome this latter term. This can be achieved by considering the following term: Again, using integration by parts to remove all the second-order derivatives of u, one deduceŝ where B 1 is a boundary term given by Now, by Cauchy-Schwarz and Young's inequalities, we have for any 0 < < e cγ 2γ(a+bγ) , where a := sup ΩT |□ψ|, b := sup ΩT |η(dψ)| and c := inf ΩT ψ. Note that there exits d = d(β, γ, Ω T ) such that |□ 2 φ(t, x)| ⩽ dφ(t, x) for all (t, x) ∈ Ω T . Therefore, applying equation (39) and to equation (38), we obtain Since |∇ψ 0 | > 2δ, we have that H(ψ) > 2α(ρ − β)δ 2 for γ > α(ρ+3β) 2 16(ρ−β)δ 2 . Therefore, . This shows that where B 1 , defined by (36), is the contribution from time-like surfaces, and B 2 (0) and B 2 (T) are the contributions from the initial and final surfaces, respectively. Finally, using the fact that φ is positive and bounded in Ω T , we have This completes the proof of theorem 3.3.
Proof of theorem 4.1. Following the standard technique, we consider the difference of u (1) and u (2) . Define U := u (1) − u (2) Now, take a function χ ∈ C ∞ [−T, T] that satisfy where ε is a small positive number. Then w := χ U t satisfies Moreover, w vanishes in {|t ± T| < ε} × Ω. In this proof, the letter C will be used to denote different constants. By setting P := ∂ 2 t − ∆ + Q 1 ∂ t + Q 0 I, u := e λϕ w and φ(t, x) = e γ(ψ0(x)−βt 2 ) in the framework of theorem 3.3, one concludes that where B(χ, U) := χ ′ 2U tt + Q 1 U t + χ ′ ′ U t and C depends on β, ψ 0 and Ω T , whilst λ depends on ψ 0 and the L ∞ (Ω)-norms of Q 0 and Q (1) 1 . Note that assumption (ii) in theorem 4.1 and proposition D.1 of appendix D imply that R and R t ∈ L ∞ ((−T, T) × Ω)). As the support of χ ′ and χ ′ ′ is contained in |t ± T| < 2ε, we have, by the estimate (18), for T > r+−r−(ε) 1 , R and R t . Using the stationary phase method [21, section 7.7], we have Using the inequalities (44)-(46), we obtain Now consider where A is a first-order operator acting on u with L ∞ (X) coefficients b 0 ∈ L ∞ (X) and b 1 ∈ L ∞ (X). Integrating the left-hand-side of equation (48) from 0 to T and using the boundary conditions in (43), we havê This meanŝ Applying estimates (45) and (47) to the above inequality, we concludê Note that ∥ fe λϕ ∥ L 2 ((−T,T)×Ω) can be estimated by the same factor as in (46). This reduces the estimate toˆΩ As assumption (i) in theorem 4.1 says that |R(0, Together with estimate (49) this yields Since r − > r − (ε), we can choose λ sufficiently large to absorb the first term in the right-hand side of (50). Finally, we let ε ↘ 0, which completes the proof.  (20). For i = 1, 2, let u (i) ∈ V be the solution to (19) . Choose a ψ 0 that satisfies condition (40). Assume also that the following conditions hold: Denote by q (i) := ∂ ν u (i) | (0,T)×∂Ω . Then, for T > T 0 := r+−r− β 1 2 , we have the following stability estimate: where C depends on β, Ω T , ψ 0 , ϕ 0 , ϕ 1 and the L ∞ -norms of Q 0 , Q (i) 1 and u (i) for i = 1, 2. Proof. We will use the same strategy as that employed in proving theorem 4.1.
For higher regularity data, we have the following proposition related to theorems 2.3 and 4.5.  (22) and (53).
Proof of theorem 4.5. As the proof is similar to the one of theorems 4.1 and 4.3, we will only highlight the main differences. Let U := u (1) − u (2) , f(x) = Q (2) 0 (x) − Q (1) 0 (x) and R = au (2) t + u (2) , so that the form of (52) changes to ⩽ e λϕ f R t + χ ′ 2U tt + Q (1) where |∂ T w| 2 is the tangential part of |∇w| 2 and G(h, φ) is some smooth function of h and φ. Note that, the boundary control in (54) is over the whole boundary ∂Ω, whilst for the Dirichlet boundary condition (13) the control in (44) is over the portion Γ + of the boundary ∂Ω. Since the right-hand-side of the above inequality shall be the replacement for the boundary term in (47). The rest is the same as in the proof of theorem 4.1.
It is not hard to see that the initial data at t = 0 in (43) is key to obtain the stability control of Q 0 or Q 1 . This inspires us to study, for example, v = χ U tt or higher-order time derivatives of U, which, in principle, will give us some estimates of Q 0 and Q 1 . However, in contrast to (43), Pv fails to be in L 2 ((−T, T) × Ω) (see appendix C). Therefore, the Carleman estimate (35) will no longer be useful. As shown in appendix C, the obstruction is exactly caused by extensions of U, Q 0 , Q 1 and R from (0, T) to (−T, T). However, if we would like to obtain a control of both Q 0 and Q 1 at an intermediate time, say t = T/2 as in [47], then we do not need to worry about singularities caused by extensions. Alternatively, one could use the method mentioned in remark 3.5 to avoid the symmetric extension. Although not detailed herein, it is worth remarking that by employing the microlocal analysis of Lasiecka and Triggiani [39] to study the property of the operator P (defined by (31)) around the boundary in the Robin boundary condition setting, one can improve the stability estimate (53) to read as where the first term ∥p t ∥ L 2 ((0,T);H 1 (Γ)) in the right-hand side of (53) has been replaced by the sharper term ∥p (55). Section 4 was devoted to establish the uniqueness and several conditional stability estimates for the inverse coefficient problems. However, the stability results expressed by the inequalities (41), (51), (53) and (55) involve the time derivatives of the measured boundary data, which, when polluted by noise, give rise in itself to an ill-posed problem of numerical differentiation that needs to be regularized. As such, the next section presents a variational formulation of the inverse problem (8), (9), (11) and (14), which enables an iterative regularizing numerical development based on the CGM.

Variational formulation
In order to solve the inverse problem (8), (9), (11) and (14) for the reconstruction of the unknown dimensionless blood perfusion coefficient w(x) along with the dimensionless temperature u(x, t), we minimize the nonlinear least-squares objective functional J : L 2 (Ω) → R + defined by where u(·, ·; w) denotes the solution of the direct problem (8), (9) and (14) for a given w ∈ L 2 (Ω). In the next theorem, we prove that the objective functional (56) is Fréchet differentiable and derive an expression for its gradient.
Theorem 5.1. The objective functional J defined in (56) is Fréchet differentiable and its gradient is given by where v(x, t) is the solution of the following adjoint problem: where δ is the Dirac delta function, which in numerical computation is approximated by where σ is a small positive constant, typically 10 −3 .
Proof. Taking a small variation δw ∈ L 2 (Ω) of w, we have where δu is the solution of the sensitivity problem: For the sensitivity problem (60), we wish to establish that For this, we start with By applying theorem 2.1 to δu satisfying the problem (60), we obtain Now, from corollary 2.4 and Sobolev's inequality, we have that for any m ⩾ max{ d 2 , 1}. Finally, combining equations (62)-(64), we conclude that where the positive constant C depends on the norms of ϕ 0 and F, as in equation (24), and of the L ∞ (Ω)-norm of w. Then, the inequality (65) implies the required statement (61). Now, multiplying the PDE in (58) by δu(x, t) and integrating over the solution domain (0, T) × Ω, using the terminal and boundary conditions in (58) and also using (60), yield From (59), (61) and (66) we obtain that J is Fréchet differentiable, and its gradient at w is given by (57).
In the following subsection, the CGM is described for the minimization of the objective functional (56).

Iterative procedure
The CGM employed to minimize the objective functional J for the reconstruction of the spacedependent perfusion coefficient w(x) is based on the recursive relation: where the direction of descent R n is given by the Fletcher-Reeves conjugate coefficient β n is given by and the search step size ζ n is computed as the minimizer To evaluate ζ n from (70), we have We set δw n = R n and linearize u(x, t; w n − ζR n ) by a first-order Taylor series expression to obtain u(·, ·; w n − ζR n ) ≈ u(·, ·; w n ) − ζR n ∂u ∂w n (·, ·; w n ) ≈ u(·, ·; w n ) − ζδu(·, ·; w n ), where δu(x, t; w n ) is obtained by solving the sensitivity problem (60) with δw n = R n . Then, differentiating J(w n − ζR n ) with respect to ζ and making it vanish yield ζ n = ´Γ 0´T 0 (u(t, y; w n ) − θ(t, y))δu(t, y; w n )dtdy ∥δu(·, ·; w n )∥ . (71)

Stopping criterion
For ensuring stability, we need to stop the iterations according to the discrepancy principle at the first iteration n * for which where¯ is a small positive value, e.g.¯ = 10 −5 , for exact data or for noisy data θ ϵ .

Algorithm
The CGM proceeds as described in the following steps: (1) Set n = 0 and select an initial guess w 0 ∈ L 2 (Ω).
(3) Stop if the criterion (72) is satisfied. Else go to step 4. (4) Solve the adjoint problem given by (58) to find v(·, ·; w n ). Compute the gradient J ′ (w n ) from equation (57), the conjugate coefficient β n from equation (69), and the direction of descent R n from equation (68). (5) Solve the sensitivity problem given by (60) to obtain δu(·, ·; w n ) by taking δw n = R n and compute the search step size ζ n from equation (71). (6) Update w n from equation (67), set n = n + 1 and go to step 2.
The dimensional temperature u(t, x) and space-dependent perfusion coefficient w b (x) can be obtained, via (7), after u(t, x) and w(x) have been reconstructed.
In the next section, the proposed inversion algorithm's convergence, accuracy and stability are illustrated and discussed on a physical example.

Numerical results and discussion
We consider a physical example concerning the recovery of the blood perfusion rate w b (x) of a one-dimensional, one-layered biological skin tissue Ω = (0, L), which undergoes a laser irradiation of the form, [13], where µ = 700 m −1 is the extinction coefficient of the tissue, I 0 = 500 W m −2 denotes the intensity of the laser and, for simplicity, the heat generation due to metabolism has been neglected.
The above dimensional quantities transform, via (7), into the following dimensionless input data: We wish to recover the dimensionless solution w(x) and u(t, x), allowing us to obtain, via (7), the dimensional blood perfusion rate w b (x) and the tissue temperature u(t, x).
The direct, adjoint and sensitivity problems present in the CGM described in section 5.1 are solved using the Crank-Nicolson finite-difference method (FDM) [14] in one-dimension (d = 1) with a uniform mesh size L/M and time step T/N. The two-point first-order backward finite difference formula is used to approximate the time-derivative u t (t, x) in (57). The trapezoidal rule is used for discretizing all the integrals present. The accuracy error, as a function of the number of iterations n, is defined as where w n stands for the numerical result obtained by the CGM at the iteration number n and w denotes the true dimensionless blood perfusion coefficient (if available).
In the absence of an analytical solution for u(t, x) being available, we generate the input measured data θ in (11) numerically by solving first the direct problem given by using the FDM, with the input data (73) and the true coefficient w(x) = 0.9822 assumed known. We then consider only half of the boundary data for each u(t, 0) = θ(t, 0) and u(t, L) = θ(t, L) obtained from solving the direct problem with M = N = 200 as our input data θ in (11), and solve the inverse problem with a coarser mesh of M = N = 100 in order to avoid committing an inverse crime. The data θ is further perturbed by additive random noise, which is numerically simulated as θ ϵ = θ + p∥θ∥ L ∞ ((0,T)×Γ0) , where p represents the percentage of noise and ε are random variables generated from a Gaussian normal distribution with mean 0 and variance 1.
We run the CGM, based on minimizing the least-squares functional (56) with Ω = (0, L) and Γ 0 = ∂Ω = {0, L} (full boundary temperature measurements) or {0} = Γ 0 ⊂ ∂Ω (partial boundary temperature measurements) from the initial guess w 0 (x) = 0.5, which is reasonably far from the true solution w(x) = 0.9822. Of course, since the inverse problem under consideration is non-linear, the least-squares objective functional (56) is non-convex and a good initial guess is required; otherwise the iterative minimization algorithm can get stuck in a local minimum. Alternatively, one may employ the globally convexification method for coefficient inverse problems originated in [32], further stated in [35, chapters 5-11] and [36, formula (4.1)], and its second generation developed in [8,9]. This method does not sufer from the phenomenon of local minima and works for both the case of vanishing initial conditions [35, chapters 7, 8, 10, 11], which are not in the framework of the Bukhgeim-Klibanov method, as well as for the cases that are in that framework, i.e. non-vanishing initial conditions, such as the one of this paper, [35, chapter 9], [8,9,36].
Figures 1(a) and 2(a) depict the monotonic decreasing convergence of the objective functional that is minimized for full Dirichlet data measured over the whole boundary Γ 0 = ∂Ω = {0, L} and for partial data measured over the portion Γ 0 = {0}, respectively, as a function of the number of iterations n, for p ∈ {0, 1, 3}% noise. For exact data, i.e. p = 0, the objective functionals in figures 1(a) and 2(a) rapidly attain very low values of 1.3 × 10 −11 and 6.9 × 10 −11 , respectively, in 20 iterations. For noisy data p ∈ {1, 3}%, the stopping iteration numbers n * ∈ {8, 7} and n * ∈ {7, 4} are generated according to the discrepancy principle (72) for Γ 0 = {0, L} and Γ 0 = {0}, respectively. Figures 1(b) and 2(b) depict the accuracy error (74), as a function of the number of iterations n, and the optimal iteration numbers n opt ∈ {8, 6} and n opt ∈ {7, 5} for Γ 0 = {0, L} and Γ 0 = {0}, respectively, can be inferred. It can be also observed that, the stopping iteration numbers n * generated according to the discrepancy principle (72) are the same as the optimal ones n opt for p = 1%, while there is only one iteration difference between n * and n opt for p = 3%. Of course, in practice only the values of n * can be obtained by the discrepancy principle (72). The corresponding numerical solutions for the dimensional space-dependent perfusion function w b (x) obtained via (7) are presented in figures 1(c) and 2(c) for Γ 0 = {0, L} and Γ 0 = {0}, respectively. First, it can be seen that in the case of noiseless data, the numerical solutions for the perfusion coefficient w b (x) are very and reasonably accurate for Γ 0 = {0, L} and Γ 0 = {0}, respectively. Second, in the case of noisy data, the numerical solutions are reasonably stable and become more accurate, as the percentage of noise p decreases from 3% to 1% and then to zero. The accuracy errors E(w n * ) ∈ {0.006, 0.058, 0.095} for p ∈ {0, 1%, 3%} noise, obtained from figure 1(b), are smaller than E(w n * ) ∈ {0.033, 0.10, 0.11} for p ∈ {0, 1%, 3%} noise, obtained from figure 2(b), indicating that, as expected, the numerical solution for the case of full boundary data Γ 0 = {0, L} being inverted is more accurate than the numerical solution obtained by inverting limited partial data over Γ 0 = {0} only. In closing, we mention that although this section has illustrated the performance of the CGM only for a physical one-dimensional biological skin tissue, the CGM can also be numerically implemented in higher dimensions, as described recently in [3] for related inverse source problem in the thermal-wave model of bio-heat transfer.

Conclusions
The inverse coefficient problem of recovering the unknown space-dependent blood perfusion coefficient in the hyperbolic thermal-wave model of bio-heat transfer from boundary temperature measurements has been investigated. Uniqueness and conditional Lipschitz stability have been established using the technique of Carleman estimates. These have been found valid over a time interval (0, T) with T > ( r+−r− β ) 1/2 that, for β ∈ (1 − r− r+ , 1), is shorter than the usual T > r 1/2 + previously reported in the literature [23,40]. Further, the problem has been reformulated as a nonlinear least-squares minimization problem and has been numerically solved using the CGM combined with the discrepancy principle for achieving stability. Numerical results associated with a physical example have been presented and discussed. Accurate and stable solutions have been obtained for both exact and random noisy data using the proposed iterative CGM.
Further work associated with the thermal-wave model of bio-heat transfer is possible, e.g. determining the shape, size and location of tumours within biological tissues using nonintrusive boundary observations.

Data availability statement
No new data were created or analysed in this study. In this appendix, we modify the method given in [45] and [17, chapter 7] for proving theorem 2.1.

Acknowledgments
Proof of theorem 2.1. The proof of existence is similar to the construction of weak solutions via Galerkin's method in [17, chapter 7] or, in a more abstract setting, in [45]. Below we give a sketch of the proof. As we will focus on the case when d = 2 or 3, by Sobolev embeddings, we have which is a compact operator. As Q 0 ∈ L ∞ (Ω) can be thought as a bounded operator from L ∞ (Ω) to L 2 (Ω), it follows that Q 0 (−∆ + iI) −1 is a compact operator and hence −∆ + Q 0 is essentially self-adjoint and has a discrete spectrum. This way of establishing spectral properties of self-adjoint operators uses the idea of relative compactness, which can be founded in [15, chapter 4] or [29, section 1, chapter 4]. The idea of using compact embeddings to establish relative compactness of Q 0 can be found in [11]. Therefore, we can choose an orthogonal basis, denoted by {w k } ∞ k=1 , of H 1 (Ω; h) and L 2 (Ω). For instance, we can choose {w k } ∞ k=1 to be the eigenfunctions of −∆ + Q 0 subject to the Robin boundary condition in (16). Standard Galerkin's method gives that for a given K ∈ N * , the function u K (t, x) := K k=1 d k K (t)w k (x) is an approximation solution to u(t, x) that yields a weak solution to (16), where d k K is determined by satisfying ⟨(u K ) tt , w k ⟩ + B[u K , w k ; t] = ⟨F, w k ⟩ for k = 1, . . ., K. The stability estimate (18) of u K can be established via energy estimates [17, chapter 7]. The energy function of u m associated with (16) is given by It is clear that E K (t) is differentiable in t, as the elliptic regularity theorem warrants the smoothness of w k for k = 1, . . ., K. Differentiating (75) yields Then, by the Cauchy-Schwarz and Sobolev trace inequalities, we have E K (t) = E K (0) +ˆt 0 E ′ K (s)ds ⩽ C 1 E K (0) +ˆt 0 ∥F(s, ·)∥ 2 L 2 (Ω) ds + C 2ˆt 0 E K (s)ds.
By Gronwall's inequality, we obtain where C(t) = C 1 (1 + C 2 te C2t ). Inequality (76) gives the stability estimate (18) for u K . Now, Banach-Alaoglu theorem gives the existence of a weak solution to (16) defined by u := lim K l →∞ u K l . This proves the existence of a weak solution. The uniqueness of u follows similarly, whilst the estimate (18) results immediately by taking the limit K l → ∞ in (76). This also implies the continuous dependence on the data F ∈ L 2 (Ω T ), ϕ 0 ∈ H 1 (Ω) and ϕ 1 ∈ L 2 (Ω).