CHAOS FOR THE HYPERBOLIC BIOHEAT EQUATION

The Hyperbolic Heat Transfer Equation describes heat processes in which extremely short periods of time or extreme temperature gradients are involved. It is already known that there are solutions of this equation which exhibit a chaotic behaviour, in the sense of Devaney, on certain spaces of analytic functions with certain growth control. We show that this chaotic behaviour still appears when we add a source term to this equation, i.e. in the Hyperbolic Bioheat Equation. These results can also be applied for the Wave Equation and for a higher order version of the Hyperbolic Bioheat Equation.


Preliminaries.
1.1. Introduction. Nowadays, surgery uses high temperature ablative techniques, such as laser, radiofrequency, microwave, or ultrasound energy to heat biological tissues to over 50 • C in a localized and safe way. Theoretical modeling can provide information about the biophysics of these techniques quickly and cheaply. Specifically, the thermal problem is modeled using the Bioheat Equation (BE) as the governing equation [25] − u xx + 1 where u represents the temperature and α and k the thermal diffusivity and conductivity of the material. The source term g refers to internal heat sources and represents different contributions for the heat sources in a biological tissue: The subscript s denotes a surgical heat source (e.g. laser or radiofrequency treatment), p refers to blood perfusion, and m to any source related with metabolic activity. The presence of g in equation (1) differences BE from the classic heat equation.
BE, as the classic heat equation, is based on the Fourier Theory which assumes an infinite thermal energy propagation speed. Although this theory might be suitable for modeling most ablative procedures, there are other surgical procedures in which extremely short periods of time or extreme temperature gradients are involved and it is necessary to consider a non-Fourier model: The Hyperbolic Bioheat Equation (HBE) [24], where τ is the thermal relaxation time, which depends on the material and represents the time that goes by since the temperature gradient is imposed until the heat flux is produced. Again, the presence of internal heat sources differences HBE from the hyperbolic heat equation. The expression of the internal heat sources varies according to the ablative technique employed. Such expression is a function of spatial and temporal variables. The spatial dependence is related with the type of ablative technique employed (radiofrequency, laser, microwaves, or ultrasound). The time-dependence refers to the energy delivering mode (e. g. continuous or pulsed). Therefore, equation (3) can present many different formulations. In [27], the laser heating heat source term in the one-dimensional case for a pulsed protocol was: where M and b are physical parameters, H(t) is the Heaviside function and ∆t is the time that the laser pulse is applied to the tissue. Figure 1 shows an schematic representation of this kind of source. In [28] the authors presented the expression for a spherical laser source and for an infinitely long cylindrical laser source being N , d, P and q biophysical parameters and R, the radial coordinate. ω(t) refers to any type of time-dependent protocol for energy delivering. In [21] there CHAOS FOR THE HYPERBOLIC BIOHEAT EQUATION   3 is another example of a source coming from radiofrequency, where the heat of a continuous source in the spherical coordinates was being Q a physical parameter. And in [20] we found the expression for a spherical microwave or ultrasound source and for an infinitely long cylindrical microwave or ultrasound source Z, m, V and w being biophysical parameters. The kind of tissue determines the ablation technique and the consideration or not of g p and g m in equation (2). In this sense, g p = 0 in non-perfused organs like the cornea, but g p is a temperature-dependent expression in well perfused organs like the liver. In the case of g m only in organs which can generate metabolic heat is taking into account usually with a constant value.
The dynamical behaviour presented by the solutions of the heat equation, when there is no heat source, has been studied on certain spaces of analytic functions with certain growth control [18]. Similar results were also obtained for the Hyperbolic Heat Transfer Equation (HHTE) in the absence of heat sources [9,17]. Certainly, it is also interesting to know some aspects of the dynamical behaviour of the HBE. In this sense the aim of this note is to study the chaotic asymptotic behaviour of certain solutions of a Cauchy problem in which the governing equation is the HBE, that is where g(t, x) represents the corresponding term to the heat source in equation (3), ϕ 1 (x) represents the initial temperature and, and ϕ 2 (x), the initial variation of temperature. We consider a general case, but with g p = 0, since we consider a linear problem. To develop this study, we will represent the solutions of (10) in terms of the solutions of the HHTE which are expressed by a C 0 -semigroup generated by certain first order equation. The solutions of (10), and the associated C 0 -semigroup, will be considered on the product of a certain function space X of analytic functions with certain growth control with itself, i.e. X ⊕ X. The general treatment of the problem shows that our results could be meaningful for all types of sources used in ablative therapies.
1.2. C 0 -semigroups. A family {T t } t≥0 of linear and continuous operators on a Banach space X is said to be a C 0 -semigroup if T 0 = Id, T t T s = T t+s for all t, s ≥ 0, and lim t→s T t x = T s x for all x ∈ X and s ≥ 0.
Let {T t } t≥0 be an arbitrary C 0 -semigroup on X. It can be shown that an operator defined by Ax := lim t→0 1 t (T t x − x) exists on a dense subspace of X; denoted by D(A). Then A, or rather (A, D(A)), is called the (infinitesimal) generator of the semigroup. It can also be shown that the infinitesimal generator determines the semigroup uniquely. If the generator A is defined on X (D(A) = X), the semigroup is expressed as {T t } t≥0 = {e tA } t≥0 [14].
The link between semigroups and differential equations is via the infinitesimal generator. The unique solution of the abstract Cauchy problem where A is a linear operator defined on X, is given by u(t, x) = e tA ϕ(x). In that sense, u(t, x) is called a classical solution of the abstract Cauchy problem (11) and the semigroup {T t } t≥0 = {e tA } t≥0 is called the solution semigroup of (11), whose infinitesimal generator is A.
In the case of a non-homogeneous Cauchy problem of the form where g(t, x) is a source term, we have the (classical) unique solution given by the following expression: 1.3. Linear dynamics of C 0 -semigroups. Given a family of operators {T t } t≥0 , we say that this family of operators is transitive if for every pair of non-void open sets U, V ⊂ X there exists some t > 0 such that T t (U ) ∩ V = ∅. Furthermore, if there is some t 0 such that the condition T t (U ) ∩ V = ∅ holds for every t ≥ t 0 we say that it is topologically mixing. A family of operators {T t } t≥0 is said to be universal if there exists some x ∈ X such that {T t x : t ≥ 0} is dense in X. When {T t } t≥0 is a C 0 -semigroup we refer to it as hypercyclic instead of universal. In this setting, transitivity coincides with universality, but it is strictly weaker than topologically mixing [6].
In addition, two notions of chaos are introduced: Devaney chaos and distributional chaos. First, we recall that an element x ∈ X is said to be a periodic point On the one hand, a family of operators {T t } t≥0 is said to be chaotic in the sense of Devaney if it is hypercyclic (universal) and there exists a dense set of periodic points in X. On the other hand, it is distributionally chaotic if there are an uncountable set S ⊂ X and δ > 0, so that for each ε > 0 and each pair x, y ∈ S of distinct points we have with µ standing for the Lebesgue measure on R + 0 . A vector x ∈ X is said to be distributionally irregular for the C 0 -semigroup {T t } t≥0 if for every δ > 0 we have Dens{s ≥ 0 : ||T s x|| ≥ δ} = 1 and CHAOS FOR THE HYPERBOLIC BIOHEAT EQUATION 5 Dens{s ≥ 0 : ||T s x|| < δ} = 1.
Such vectors were considered in [7] so as to get a further insight into the phenomenon of distributional chaos, showing the equivalence between a distributionally chaotic operator and an operator having a distributionally irregular vector. This equivalence has been shown for C 0 -semigroups in [1].
A criterion for Devaney chaos in terms of the abundance of eigenvectors of the infinitesimal generator of a C 0 -semigroup was stated in [13] by Desch, Schappacher, and Webb. Since then, this criterion has been reformulated and applied to several examples of C 0 -semigroups which are solution of certain partial differential equations, see for instance [2,10,17]. The following version can be found in [17,Th. 7.30].
Assume that there exists a nonempty open connected subset U of C and weakly holomorphic functions f j : U → X, j ∈ J, such that Then {T t } t≥0 is topologically mixing and Devaney chaotic.
The third condition in this result is used in order to prove the density of the span of certain sets of eigenvectors associated to eigenvalues of A with real part greater, equal, and smaller than 0. A criterion stated in these terms was firstly stated for operators by Godefroy and Shapiro in [15]. Theorem 1.2. Eigenvalue criterion for chaos. Let X be a complex separable Banach space and {T t } t≥0 a C 0 -semigroup on X. Suppose that the sets X 0 := span{x ∈ X : ∃λ > 0, T t x = e λt x, ∀t ≥ 0}, In the proof of Theorem 1.1, Condition (3) of its statement is used to satisfy the hypothesis of Theorem 1.2, and hence to prove that the C 0 -semigroup is Devaney chaotic. As a result, we can replace Condition (3) in Theorem 1.1 with the hypothesis of Theorem 1.2.
In addition, there exist several criteria for distributional chaos [1,7]. Nevertheless, either the Desch-Schappacher-Webb criterion or the Eigenvalue criterion for chaos imply distributional chaos, [4,Rem. 3.8], see also [5,Cor. 31]. Moreover, in these cases we can affirm that there is a dense distributionally irregular manifold, that is a dense manifold of distributionally irregular vectors.
More information on sufficient conditions for hypercyclicity and chaos for C 0semigroups and operators can be found in [2,6,12,13,17,19].

The Hyperbolic Heat Transfer Equation.
The chaotic behaviour of the solutions of an abstract Cauchy problem (10) that is given by the Hyperbolic Heat Transfer Equation in the absence of internal heat sources was analyzed on certain spaces of analytic functions with certain growth control in [9], see also [17]. This can be done if we express this second-order PDE as a first-order equation by representing it as a C 0 -semigroup on the product of a certain function space with itself. To do this we set u 1 = u and u 2 = ∂u ∂t . Then the associated first-order equation is We fix ρ > 0 and consider the space endowed with the norm ||f || = sup n≥0 |a n |, where c 0 is the Banach space of complex sequences tending to 0. Then X ρ is a Banach space of analytic functions with a certain growth control. By its definition it is isometrically isomorphic to c 0 . This type of spaces were already used in [18]. X ρ is a Banach space of analytic functions that is densely embedded in C(R) with the topology of uniform convergence on compact sets of R, since it contains all polynomials. Essentially, X ρ is a space of analytic functions with certain increasing control at infinity. In fact, (X ρ , || · ||) is isometrically isomorphic to (c 0 (N 0 ), || · || ∞ ). Since is a linear and continuous operator on X ρ ⊕ X ρ , we have that {e tA } t≥0 , with is well defined on X ρ ⊕ X ρ , and we have that (e tA ) t≥0 is a C 0 -semigroup on X ρ ⊕ X ρ (even uniformly continuous), which is the solution semigroup of (14)  For simplicity we will denote by || · || ρ,ρ the norm in X ρ ⊕ X ρ given by the norm || · || ρ on each copy of the space X ρ . Then the solution semigroup {e tA } t≥0 of (14) is topologically mixing and Devaney chaotic on X ρ ⊕ X ρ .
The proof of Theorem 2.1 is included here because we will use the same notation in many of the upcoming results. This proof consists on checking the Eigenvalue criterion for chaos stated in Theorem 1.2. For this purpose, we take the holomorphic functions of the form where λ ∈ C, z 0 , z 1 ∈ R, and R λ = (τ λ 2 + λ)/α. If we take λ ∈ V ⊂ C, with V the open disk of radius r = αρ 2 /2τ > 0 centered at zero, then ϕ λ,z0,z1 ∈ X ρ . Then, the functions φ z0,z1 : V → X, z 0 , z 1 ∈ R, given by satisfy that φ λ,z0,z1 ∈ ker(λI −A) for every λ ∈ V, z 0 , z 1 ∈ R. Finally, the Eigenvalue criterion is applied taking the sets X 0 , X 1 , X p mentioned in Theorem 1.2 as With a similar approach, one can obtain the chaotic behaviour for the solutions of the Wave Equation: where α > 0 is the square of the speed of wave propagation. Again, this equation can be reformulated as a first order equation taking the following operator B instead of A.
In this case, the solution C 0 -semigroup to the abstract Cauchy Problem in (23) presents a Devaney chaotic behaviour on the space X ρ ⊕ X ρ for any ρ > 0 [9, Th. 2.3] and [17, Ex. 7.5.3].

The Hyperbolic Bioheat Equation.
Let us consider the HBE given in (10). As in the previous case, we express this second-order equation as a first-order equation by setting u 1 = u and u 2 = ∂u ∂t and taking A as it has been already defined in (16). In this way, the associated first-order equation is formulated as where the function g(t, x) represents the internal heat sources. Comparing this expression with the formulation of the HBE in (3), we see that g(t, x) stands for τ k (g + τ g t ).
The unique (classical) solution of the HBE in (25) is given by: where we have used the following notation

Dynamics of the Hyperbolic Bioheat Equation.
Since we know that the C 0 -semigroup {e tA } t≥0 is chaotic on the Banach space X ρ ⊕ X ρ , in order to study the asymptotic behaviour of the solutions of the HBE we have to analyze the asymptotic properties the second term in (26): In addition, as the solution u(t, x) is expected to be an element of X ρ ⊕ X ρ , it is reasonable to consider the source term Ψ(t, x) in this same space. Later, we will see that in fact we are still considering some type of HBE. Moreover, for many applications it is possible to consider the source term g(t, x) as time-independent. So that, we will firstly consider Ψ(t, x) as a time-independent function of the form φ λ,z0,z1 (x), for certain values of λ with negative real part.
On the one hand, we will see that in Theorems 2.2 -2.5 we consider which has introduced a variation in the relationship between u 1 (t, x) and u 2 (t, x) in equation (25): we have (u 1 ) t = u 2 + ϕ λ,z0,z1 instead of just (u 1 ) t = u 2 . At first sight, it may seem that the resulting differential equation is no longer the HBE. Fortunately, that is not the case: Let us express the heat source term in the form: where we assume that j 1 , its time derivative (j 1 ) t , and j 2 are small. Then, from the equation (25), it is just an exercise to check that u 1 (t, x) = u(t, x) satisfies the HBE: which corresponds to a small perturbation of the heat source term g(t, x) in equation (10). On the other hand, taking a term Ψ(t, x) as in (28) is not a great restriction. This is due to the fact that we will consider eigenfunctions of A, φ λ,z0,z1 , in a set whose span is dense in X ρ ⊕ X ρ . Therefore, given any source function g(t, x), the initial term 0 g(t, x) can be approximated by a linear combination of those eigenfunctions multiplied by a certain function depending of t. Then, for this new source term, we will be able to find an initial condition whose orbit was dense under the action of the operators in the family u(t, x). Before going further, we point out that the notation of the statement of the following theorem and in the rest of the results of this section is the same used in the proof of Theorem 2.1.
Proof. Let us take two non-empty open sets U 1 , U 2 ⊆ X ρ ⊕X ρ . Let W ⊂ X ρ ⊕X ρ be an open 0-neighborhood and let U 2 be an open set in X ρ ⊕X ρ such that W +U 2 ⊂ U 2 .
Consider an arbitrary λ ∈ V with (λ) < 0. On the one hand, since the C 0semigroup {e tA } t≥0 is known to be topologically mixing on X ρ ⊕ X ρ , then there exists some t 0 > 0 such that for all On the other hand, since φ λ,z0,z1 (x) does not depend on the time, then we can easily analyze the integral h(t, x): Clearly, h(t, x) ∈ X ρ ⊕ X ρ for every t ≥ 0 and, since and (λ) < 0, we can affirm that there exists t 1 > 0 such that h(t, ·)+ 1 λ φ λ,z0,z1 ∈ W for all t ≥ t 1 . Therefore, for t ≥ max{t 0 , t 1 }, we have that u(t, ·) acting on the initial condition Φ t ∈ U 1 yields which fulfills the definition of topologically mixing for the solution family {u(t, ·)} t≥0 . on X ρ ⊕ X ρ .
Remark 1. The condition that λ belongs to V is necessary to be sure that φ λ,z0,z1 (x) belongs to X ρ ⊕ X ρ . This can be replaced by the more general one |τ λ 2 + λ| < αρ 2 .

Remark 2.
Under the hypothesis of Theorem 2.2 we can get a little more than topologically mixing: since the asymptotic behaviour of the orbits by the C 0 -semigroup {e tA } t≥0 coincides with the asymptotic behaviour of the orbits under the family of operators {e tA + h(t, ·)} t≥0 except by a constant, then we can affirm that every hypercyclic function for the C 0 -semigroup {e tA } t≥0 is universal for the family of operators {e tA + h(t, ·)} t≥0 . Furthermore, by [11,Th. 2.3] we have that a hypercyclic/universal behaviour is presented on every nontrivial autonomous discretization of {e tA + h(t, ·)} t≥0 , i.e. for the sequence of operators {e kt0A + h(kt 0 , ·)} k∈N for every t 0 > 0. Furthermore, the set of hypercyclic/universal functions is shared by the family of operators itself and by every nontrivial autonomous discretization. Even more, any nontrivial single operator is also topologically mixing c.f. [6, Th.

3.5].
Remark 3. One can also consider the cases Ψ(t, x) = f (t)φ λ,z0,z1 (x), with λ ∈ V and (λ) < 0, with a similar proof to the one of Theorem 2.2: This can be done by analyzing the asymptotic behaviour of h(t, x) when t tends to ∞. i) f (t) ∈ L 1 (R + 0 ): We will see that lim t→∞ h(t, x) = 0 for every x ∈ R. Let us fix x ∈ R and ε > 0. Since f (t) ∈ L 1 (R + 0 ), we can find some t > 0 such that ∞ t |f (s)|ds < ε 2|φ λ,z 0 ,z 1 (x)| . Now, as (λ) < 0, we have that there is some t 0 > t such that for all t ≥ t 0 we have that e (λ)(t−t ) < ε 2||f ||1|φ λ,z 0 ,z 1 (x)| . Taking this into account we have for all t ≥ t 0 . ii) f (t) ∈ L p (R + 0 ) for every 1 < p < ∞: We will also see that lim t→∞ h(t, x) = 0. Again, let us fix x ∈ R and ε > 0. As before, there is some t > 0 such that for all t ≥ t 0 . Applying Hlder inequality we get iii) f (t) a bounded locally integrable function that asymptotically tends to a constant L 0 : Fix x ∈ R and ε > 0. We will see that lim t→∞ h(t, x) = −L0 λ φ λ,z0,z1 . Since lim t→∞ f (t) = L 0 there exists some t > 0 such that |f (s) − L 0 | < ε|λ| 8|φ λ,z 0 ,z 1 (x)| for all s ≥ t . Then there exists some t 0 > t such that for all t ≥ t 0 we have We have to estimate which is smaller or equal than The first integral is smaller or equal than The second one can be bounded by To sum up we have h(t, x) + L0 λ φ λ,z0,z1 < ε for all t ≥ t 0 .
Example 1. Let us consider a time dependent source term Ψ(t, x) = f (t)φ λ,z0,z1 (x) such that the time dependent function f (t) is obtained from the laser heating source term g s (t, x) given in (4) and the corresponding right hand side of (3), τ k (g + τ g t ). Thus, for all x ∈ R and t ≥ 0. In this case, the integral h(t, x) involves the distribution δ, and If (λ) < 0, lim t→∞ h(t, x) = 0, and the asymptotic behaviour of the solution family {u(t, ·)} t≥0 of (35) is the same that the solution semigroup {e tA } t≥0 of (14).
Example 2. Analogous calculations can be done for the heat source term coming from radiofrequency given in (7), then and, the h(t, x) integral is given by: Thus, the dynamical properties of this example are the same that of the time independent case.
If we consider the abstract Cauchy problem given by the Wave Equation where we have included a source term g(t, x), one can also obtain similar results to Theorems 2.2 and 2.3 using the results about the chaotic behaviour of the solutions of the wave equation in (23) on the spaces X ρ ⊕ X ρ . One just have to take the operator B in (24) instead of A and to consider its corresponding eigenfunctions, which are of the form φ λ,z0,z1 with R λ = λ 2 α .
Theorem 2.5. Let ρ > 0, W the open disk of radius αρ 2 centered at 0, and Ψ(t, x) = φ λ,z0,z1 (x) with λ ∈ W ∩ iR. The solution family {u(t, ·)} t≥0 of (35) is Devaney chaotic on X ρ ⊕ X ρ . The results obtained for the HBE can be also applied to the relativistic heat equation (RHE). The formulation of the RHE only differs from the HBE formulation in the meaning of the coefficients that accompanied second derivatives in equation (10) (see [22]). The RHE is a hyperbolic-like equation, whose theoretical model is based on the theory of relativity and which was designed to overcome the possible conflict between the HHTE and the second law of thermodynamics.
We conclude this section with a comment regarding the stability of the solutions. We recall that a C 0 -semigroup of the form {e tA } t≥0 defined on a Banach space X is (uniformly) exponentially stable, [14, p. 296], if there exists ε > 0 such that lim t→∞ e εt ||e tA || = 0, This situation is not fulfilled in our case. Nevertheless, a weaker version of this condition can also be considered. We say that {e tA } t≥0 is exponentially stable on a subspace Y ⊂ X if there exists ε > 0 such that for any y ∈ Y we have lim t→∞ e εt ||e tA y|| = 0, Such study is sometimes considered when analyzing the chaotic behaviour of C 0semigroups, see for instance [8,3]. The HHTE can be seen to be exponentially stable on the subspaces X 1,δ := span{φ λ,z0,z1 : λ ∈ V with (λ) < δ, z 0 , z 1 ∈ R}, for every − αρ 2 /2τ < δ < 0. For the case of λ ∈ V with (λ) > 0 we have that the behaviour of the solutions for an initial condition on X 1 escapes to ∞ in norm.
3. An extension: dynamics of solutions of the higher order HBE. An extension of the problem considered above is the dynamics of the solutions of higher order linear heat transfer or diffusion equations. Many interesting physical problems are modeled by a fourth, or even higher, diffusion equation. These equations appear in a wide range of areas, including fluid dynamics, electromagnetism and semiconductors, optical tomography, image processing, etc (as an example, see [26,16] and references therein). In particular, an interesting class of problems described by nonlinear fourth order diffusion equations are thin fluid film flows where surface tension, due to gravity or temperature gradient, is a driving mechanism [23]. The simplest example is the one-dimensional linear fourth order diffusion equation, u t = −αu xxxx u(0, x) = ϕ(x), x ∈ R , where α is the diffusivity. The linear operator A = −α∂ xxxx is the infinitesimal generator of the solution semigroup {T t } t≥0 = {e tA } t≥0 .
are holomorphic and satisfy that φ λ,Z ∈ ker(λI − A m ) for every m ∈ N + . Furthermore, for a fixed ρ > 0 such that α m τ m ρ 2(m+1) > 2 and for λ ∈ U , being U the open disk of radius α m ρ 2(m+1) /2τ m > 0 centered at 0, functions φ λ,Z belong to the Banach space X ρ ⊕ X ρ . Therefore, using the chaos criterion (Theorem 1.1), the C 0 -semigroups {e tAm } t≥0 are mixing a chaotic on X ρ ⊕ X ρ for all m ∈ N + . Discussion for the part of the solution due to the source terms, i.e. the integral t 0 e (t−s)Am Ψ m (s, x)ds, is analogous of the discussion of the HBE equation. As a result, the solution family {u m (t, ·)} t≥0 (44) of the 2(m+1)-order hyperbolic bioheat equation (42) is topologically mixing for λ values with (λ) < 0 and Devaney chaotic with λ ∈ iR on X ρ ⊕ X ρ for each m ∈ N + (under some additional restrictions analogous to the hypothesis of Theorem 2.2 and Theorem 2.3).

Conclusions.
In this paper we have studied the linear dynamics of the solutions of the HBE in certain spaces of analytic functions that are densely embedded in C(R) with the topology of uniform convergence on compact subsets of R. Specifically, we have shown the chaotic behaviour for the solutions of a Cauchy problem in which the governing equation is given by the HBE with a time-independent heat source term that is certain eigenfunction for the differential operator A. Since these eigenfunctions belong to a set whose span is dense in X ρ ⊕ X ρ , see Theorem 1.2, given any initial source term g(t, x), we can find another source term of the form Ψ(t, x) = n i=1 φ λi,zi,0,zi,1 (x) with λ i ∈ V and (R) ≤ 0 for all 1 ≤ i ≤ n such that n i=1 ϕ λi,zi,0,zi,1 (x) n i=1 λ i ϕ λi,zi,0,zi,1 (x) can be found as close as we want of 0 g(t, x) and for this new source term we can find an initial condition with dense orbit under the action of the operators in the family u(t, x).
The general treatment of the problem allows us to apply the results to all types of ablative therapies (laser, radiofrequency, microwave, and ultrasound), and with different modes of energy delivering without making any additional computation.
An extension of these results for Hyperbolic Bioheat Equations of higher order is also provided.