Local well-posedness of a coupled Westervelt-Pennes model of nonlinear ultrasonic heating

High-Intensity Focused Ultrasound (HIFU) waves are known to induce localized heat to a targeted area during medical treatments. In turn, the rise in temperature influences their speed of propagation. This coupling affects the position of the focal region as well as the achieved pressure and temperature values. In this work, we investigate a mathematical model of nonlinear ultrasonic heating based on the Westervelt wave equation coupled to the Pennes bioheat equation that captures this so-called thermal lensing effect. We prove that this quasi-linear model is well-posed locally in time and does not degenerate under a smallness assumption on the pressure data.


Introduction
High-Intensity Focused Ultrasound (HIFU) is an innovative medical tool that relies on focused sound waves to induce localized heating to the targeted tissue [34]. Due to its non-invasive nature and relatively brief treatment time, it has excellent potential to be used in the therapy of various benign and malignant tumors; see, e.g., [10,13,21,23,36]. The ability to accurately determine the properties of the pressure and temperature field in the focal region is crucial in these procedures and motivates the research into the validity of the corresponding mathematical models.
It is well-known that the heating of tissue influences the speed of propagation of sound waves and, in turn, the position of the focal region; this effect is commonly referred to as thermal lensing [7,11,12]. In this work, we analyze a mathematical model of nonlinear ultrasonic heating that captures this effect. More precisely, we study a coupled problem consisting of the Westervelt wave equation of nonlinear acoustics [35]: (1.1) p tt − c 2 (Θ)∆p − b∆p t = k(Θ) p 2 tt and the Pennes bioheat equation [27]: where p is the acoustic pressure, Θ the temperature, and Q(p t ) is the absorbed acoustic energy. We refer to Section 2 below for further details on the modeling and the involved material parameters.
To the best of our knowledge, this is the first work dealing with a rigorous mathematical analysis of a coupled Westervelt-Pennes model. Westervelt's equation is a quasilinear strongly damped (for b > 0) wave equation that has been extensively studied by now in various settings with constant material parameters; see, e.g., [14,15,17,18,24] and the references given therein, where results concerning local well-posedness, global well-posedness, and asymptotic behavior of the solution have been established. The results on the well-posedness of the Westervelt equation with an additional strong nonlinear damping and with L ∞ (Ω) varying coefficients have been obtained in [4,25]. We mention that this wave equation can also be rigorously recovered in the limit of a thirdorder nonlinear acoustic equation for vanishing thermal relaxation time; see the analysis in [3,16].
A distinguishing feature of the present quasilinear thermo-acoustic problem is the dependence of propagation speed on the temperature, which we assume to be polynomial (in accordance with the real-world setting) and non-degenerate. Having temperaturedependent medium parameters presents a challenge due to the higher regularity the pressure should possess to tackle the nonlinearities. To resolve this issue, our theoretical approach relies on higher-order energy analysis of a suitable linearization combined with a fixed-point argument, under the assumption of smooth and small (with respect to pressure) data. Although the heat equation (1.2) has regularizing properties, it does not seem feasible to transfer these to the pressure equation (1.1) and make use of the damping property of heat conduction, as in the classical thermo-elastic systems; see, e.g., [19,20,29] and the references given therein. This issue arises due to the very weak nonlinear coupling in the present model.
A critical step in the analysis is handling the higher-order time-derivative of the pressure in the nonlinear term; that is, k(Θ) p 2 tt . Due to the temperature-dependent coefficients, we have to rely on higher-order energies compared to the analysis of Westervelt equation in homogeneous media in [14] and assume (p, p t ) |t=0 = (p 0 , p 1 ) ∈ H 3 (Ω) × H 2 (Ω).
Another well-known difficulty in the analysis of the Westervelt equation is related to the possible degeneracy for large values of p since the factor in front of p tt in (1.1) can be written as 1 − 2k(Θ)p. This invokes the condition 1 − 2k(Θ)p > 0 almost everywhere, which in turn requires p L ∞ to remain small enough in time. The issue is commonly resolved by using a Sobolev embedding under the assumption of small pressure data, e.g., H 2 (Ω) ֒→ L ∞ (Ω), as in [14].
We organize the rest of our exposition as follows. We provide more detailed insight into mathematical bio-acoustic modeling in Section 2. Section 3 focuses on the energy analysis of a (partially) linearized uncoupled problem. In Section 4, we present the study of the coupled nonlinear model by relying on the result from the previous section and Banach's fixed-point theorem. Our main well-posedness result is contained in Theorem 4.1. We conclude the paper with a discussion and an outlook on future work.

A model of ultrasonic heating based on the Westervelt equation
Volume coupling of the acoustic pressure p to the temperature field Θ is achieved via appropriate source terms and the use of temperature-dependent acoustic material parameters; see, e.g., [6,7,11,26,32]. In this work, we model the soft tissue as a thermoviscous fluid and consider the Westervelt equation in pressure form with the temperature-dependent speed of sound c = c(Θ): , where the right-hand side coefficient is given by Here, ρ is the medium density and β acou the acoustic coefficient of nonlinearity. The third term in (2.1) models thermoviscous losses. The parameter b is known as the sound diffusivity [22]. Assuming harmonic excitation with angular frequency ω, it is connected to the absorption coefficient α via where c a is the ambient speed of sound (in the tissue); cf. [26]. Equation (2.1) is obtained from the model considered in [7] under the assumption of constant medium density and upon approximating b c 2 p ttt by b∆p t . Note that if the attenuation obeys a frequency power law, equation (2.1) generalizes to involve a fractional damping term; see, e.g., [26]. This case is thus of interest for future analysis as well, but outside the scope of the current work.
The temperature distribution in the tissue is modeled by the Pennes bioheat equation [27] with a nonlinear source term: The function Q = Q(p t ) represents the acoustic energy absorbed by the tissue at any given point. The term ρ b C b W (Θ−Θ a ) models the removal of heat by blood circulation. Here, ρ b and C b are the density and specific heat capacity of blood, respectively, and W is the volumetric perfusion rate of the tissue measured in milliliters of blood per milliliter of tissue per second. The values of typical material properties in the human tissue can be found, for example, in [7, Table 3]. The coefficients ρ a and κ a denote the ambient density and thermal conductivity (i.e., the tissue density and thermal conductivity). C a is the ambient heat capacity and Θ a is the ambient temperature. In the body, the latter is usually taken to be 37°C; see [7].
Altogether, we end up with the following coupled problem: where we have introduced the function For simplicity, we consider (2.3a) together with homogeneous Dirichlet boundary conditions (2.3b) p| ∂Ω = 0, Θ| ∂Ω = 0, and the initial data  [2]. We thus make the following assumptions on the function q in our analysis. Note that throughout the paper, we use x y to denote x ≤ Cy, where C > 0 is a generic constant that may depend on Ω, the final time T , and medium parameters.
We assume that there exists q 0 > 0, such that Furthermore, there exist γ 1 ≥ 0 and C 1 > 0, such that By these assumptions and Taylor's formula, it further follows that The function k is assumed to be related to q via (2.2) throughout this work. Therefore, we have we conclude that there exists γ 2 > 0, such that Modeling the absorbed acoustic energy. The acoustic energy absorbed by the tissue is represented by the source term Q = Q(p t ) in the heat equation. We will make the following general assumptions concerning its properties in our analysis, which allow us to cover important particular cases from the literature.
The mapping Q is Lipschitz continuous on bounded subsets of the space L ∞ (0, T ; L ∞ (Ω)) with values in L 2 (0, T ; L 2 (Ω)), that is, and such that Q(0) = 0. Additionally, Note that by plugging in v = 0 above, these assumptions further imply that In [26,28], the absorption term is modeled as which clearly satisfies our assumptions if p t ∈ L ∞ (0, T ; L ∞ (Ω)) and p tt ∈ L ∞ (0, T ; L 2 (Ω)). More commonly, the absorption term appears in the literature averaged over a certain time interval. In, e.g, [7, §2.2], the absorbed energy is given by Here j is a positive integer, τ is the period of ultrasound excitation and t ′ is a sufficient time from the start of the simulation so that a steady-state has been reached. In [11], the absorbed energy is averaged over the whole time interval Both of these functionals satisfy Assumption 2. In case of (2.9), for example, we note that for all t ∈ [0, T ], and by using Minkowski's inequality (see [ In case of a time-averaged absorbed energy, we have Auxiliary results. We collect here several useful inequalities that are repeatedly used in the analysis below. We assume throughout that Ω ⊂ R d , where d ∈ {1, 2, 3}, is an open, bounded, and sufficiently smooth set. We will often rely on the Ladyzhenskaya inequality for u ∈ H 1 (Ω): H 1 . By using (2.10) together with Young's inequality, we further find that for u ∈ H 1 0 (Ω) and any ε > 0 (2.11) with ε = Cε 4/d . This estimate can also be obtained (on bounded domains) by employing Ehrling's lemma; see [30,Lemma 8.2]. Further, given a ∈ H −1 (Ω) and b ∈ W 1,3 (Ω) ∩ L ∞ (Ω), the following bound holds: To keep the presentation self-contained, we also state here the version of Gronwall's inequality that will be employed in the proofs. Let v be non-negative and integrable. Suppose that u : I → R is in C 1 (I) and satisfies: Then it holds that Proof. The inequality follows by combining the arguments of [5, Appendix B] and [9, Lemma 3.1].

Analysis of a linearized problem
We first analyze a de-coupled linearization of (2.3a), given by in Ω × (0, T ), and supplemented by the boundary (2.3b) and initial (2.3c) conditions. To facilitate the analysis, we make the following regularity and non-degeneracy assumptions on the involved coefficients and source terms.
Assumption 3. Given T > 0, the variable coefficients and the source terms satisfy the following assumptions.
(F) Let f 1 ∈ L 2 (0, T ; H 1 0 (Ω)), ∂ t f 1 ∈ L 2 (0, T ; H −1 (Ω)), and f 2 ∈ H 1 (0, T ; L 2 (Ω)). From the last assumption, by [31,Theorem 7.22], we have f 1 ∈ C([0, T ]; L 2 (Ω)) and Energies. To accommodate the energy analysis, we introduce the following lower and higher-order acoustic energies: In the analysis, we will also use the combined acoustic energy with the associated dissipation rate The initial acoustic energy is set to Further, the heat energy is given by with the associated dissipation Solution spaces. To formulate the existence result, we also introduce the following solutions spaces for the pressure: and the temperature: with the short-hand notation : tr ∂Ω u = 0, tr ∂Ω ∆u = 0 . We claim that the linearized problem is well-posed under the above-made assumptions. Proposition 3.1. Let T > 0 and let Assumption 3 hold. Further, assume that . Then there exists a unique solution (p, Θ) ∈ X p ×X Θ of (3.1). Furthermore, the acoustic pressure satisfies a.e. in time, with Proof. Since the system is de-coupled, we can analyze the equations in (3.1) sequentially.
Analysis of the pressure equation. The analysis of the pressure equation can be rigorously conducted by employing a Galerkin discretization in space based on the smooth eigenfunctions of the Dirichlet-Laplacian; see. e.g., [8,Ch. 7]. We focus here on presenting the energy analysis.
Energy analysis. Testing the (semi-discrete) pressure equation with p t , integrating over Ω, and using integration by parts yields the following identity: a.e. in time. From here, by Hölder's and Young's inequalities, we have On account of Assumption 3, we know that and thus for any ε > 0, it holds that where we have applied Poincare's inequality together with Young's ε-inequality in the estimate of the last term. Note that by fixing ε > 0 small enough, we can absorb the last term in (3.7) by the dissipative term on the left. By using the embedding H 1 (Ω) ֒→ L 4 (Ω) together with the Poincaré inequality, the first term on the right-hand side of (3.7) can be absorbed by the dissipative term √ b∇p t (t) 2 L 2 as well if we assume the norm α t /b L ∞ (L 2 ) to be small. However, to avoid this smallness assumption, we use inequality (2.11) instead and split this term into two parts: an energy term and a dissipation term with an arbitrary small factor ε > 0. This idea will be used repeatedly in the proof below. Indeed, by using inequality Consequently, by recalling Assumption 3 and fixing ε > 0 small enough, so that where C is the hidden constant in (3.7), we obtain where we have also used again the uniform bound on α given in Assumption 3. Estimate (3.8) indicates that further testing is needed to absorb the energy E 1 on the right. Thus, we test the first (semi-discrete) equation in (3.1) with −∆p t and integrate in space, which yields Clearly, by selecting ε > 0 small enough in the above estimate, the second term on the right-hand side will be absorbed by the dissipation on the left. Hence, by choosing ε > 0 as small as needed, keeping in mind that ∆p = 0 on ∂Ω, and using Poincaré's inequality, we obtain To retrieve the energy E 1 on the left, we will next work with the time-differentiated pressure equation. Indeed, on account of the regularity assumptions on the coefficients and source term, we can differentiate the semi-discrete pressure equation with respect to t: Multiplying (3.10) by p tt , integrating over Ω, and using integration by parts with respect to time in the first term, we obtain The first two r terms on the right can be estimated as follows: for some ε > 0, where we have relied on the embedding H 1 (Ω) ֒→ L 4 (Ω). By applying estimate (2.11), we can further bound the last term: where we have also utilized elliptic regularity (since ∂Ω is smooth): The first and the fifth term on the right-hand side of (3.11) can be estimated as follows: We then further estimate the last term above using again inequality (2.11): Keeping in mind Assumption (3), and plugging (3.13) into (3.12), we have By using Young's inequality together with the Poincaré's inequality, we find that Recalling that ∆p = 0 on ∂Ω, we can estimate the last term on the right-hand side of (3.11) as follows: We see that the first term on the right can be absorbed by the dissipation in (3.11) and the last one is an energy term. By collecting the above estimates with ε > 0 small enough, we arrive at (3.14) Adding (3.14) to (3.9), exploiting Assumption 3, using Poincaré's inequality, and possibly reducing ε, so that the ε terms can be absorbed by the left side, we obtain To be able to absorb the term √ b∇∆p(t) 2 L 2 on the right, we should additionally test the pressure equation with ∆ 2 p: Integrating by parts and using the fact that p tt = ∆p = ∆p t = 0 on the boundary for smooth Galerkin approximations, as well as that Recalling how the energy E 2 is defined in (3.3), from here we obtain By Hölder's inequality, we further have Using Young's and Poincaré's inequalities, and the embedding H 1 (Ω) ֒→ L 6 (Ω) yields By adding inequalities (3.15) and (3.16), and selecting ε > 0 small enough, we have By collecting the above estimates, we arrive at a bound that involves the combined acoustic energy: where Λ(t) and F(t) are defined in (3.5) and (3.6), respectively. By Gronwall's inequality, we then immediately have Additional bootstrap arguments. We can obtain more information on the pressure field by relying on the (semi-discrete) PDE. Indeed, by the acoustic PDE we have We can then further estimate the right-hand side of (3.19) by employing the acoustic energy: where we have also used estimate (3.2) to bound the f 1 (t) 2 L 2 term. Adding this bound to (3.18) yields Similarly, Adding γ·(3.20) to (3.17) with small enough γ > 0 yields on which we can apply Gronwall's inequality.
Additionally, from the time-differentiated equation (3.10), standard arguments (see, e.g., [8,Ch. 7,p. 383]) give the following bound in the dual space H −1 (Ω): where we have used the embedding L 6/5 (Ω) ֒→ H −1 (Ω) together with Hölder's inequality to get Thus, we have with a uniform bound By using again the embedding L 6/5 (Ω) ֒→ H −1 (Ω) and Hölder's inequality, we have, similarly to before, and thus Then adding γ·(3.21) to (3.17) with γ > 0 small enough, and using Gronwall's inequality yields Combining the three derived estimates yields (3.4), as first in a semi-discrete setting. The obtained uniform bound allows us to employ standard compactness arguments and prove existence of a solution p ∈ X p to the pressure equation; see, e.g., [8,Ch. 7] for similar arguments. By the weak/weak-⋆ lower semi-continuity of norms, p satisfies the same energy bound (3.4). Note that p ∈ X p implies Uniqueness. Uniqueness in the pressure equation follows by showing that the only solution of the homogeneous problem is zero. To this end, let p ∈ X p solve We can repeat our previous energy analysis up to (3.17), where instead of testing with ∆ 2 p (which is not a valid test function), we take the gradient of the equation and test with ∇∆p ∈ L ∞ (L 2 (Ω)). In this manner, from (3.4) we obtain E[p](t) = 0, which immediately yields p = 0.
Analysis of the heat equation. We next rewrite the heat equation as According to, e.g., [37, Ch. 1, Theorem 1.3.2], the unique solution Θ ∈ X Θ of this problem satisfies for all t ∈ [0, T ]; see also [33,Ch. 2,Theorem 3.2]. Thanks to the assumed properties of the mapping Q, we have .
Thus, by the embedding H 1 (0, T ) ֒→ C[0, T ], from (3.22) we have where B will be a suitably chosen ball in X T , with the solution (p, Θ) ∈ X p × X Θ of (4.1) ( with the boundary (2.3b) and initial (2.3c) conditions. Our main results reads as follows.
then there exist a unique solution (p, Θ) of (2.3) in X T . Furthermore, the solution depends continuously on the data with respect to · X T .
Proof. As already announced, we intend to rely on Banach's fixed-point theorem to arrive at the claim. To facilitate the fixed-point argument, we define the pressure and temperature norms: and We can then also define the combined norm as follows: To have an equivalence between this norm and the energies, we introduce the total pressure energy E[p] as and the associated dissipation rate as Then on account of Assumption 3, there exist positive constants C 1 , . . . , C 4 , such that

E[p](T ) + D[p](T )) and
We next introduce a ball in X T : where the radii R 1 > 0 and R 2 > 0 are to be determined by the proof. The constant k 1 > 0 is such that |k(Θ)| ≤ k 1 ; cf. assumption (2.5). In the course of the proof we will impose a smallness condition on the pressure, but not on the temperature data, which is why we have introduced two different radii here.
Note that the solution of the linear problem with α = r = 1 and f 1 = f 2 = 0, belongs to this ball if δ > 0 is small enough and R 2 large enough, so that (Ω) + δ 2 + 1), so this set is non-empty. We consider the ball to be equipped with the distance Then (B, d) is a complete metric space. We first prove that T is a self-mapping. Proof. We wish to rely on the well-posedness result from the previous section. To this end, we set α(x, t) = 1 − 2k(Θ * )p * , r(x, t) = q(Θ * ), f 1 (x, t) = 2k(Θ * )p 2 * t , f 2 (x, t) = 0. to fit problem (4.1) into the framework of Proposition 3.1. We next verify Assumption 3 on these functions. Since and so the non-degeneracy condition is fulfilled. Further, by the embeddings H 1 (Ω) ֒→ L 3 (Ω) and H 2 (Ω) ֒→ L ∞ (Ω), we have From here and properties (2.6) of the function k, it follows that Again by the embedding H 1 (Ω) ֒→ L 3 (Ω) and properties of the function k, it holds that We can analogously estimate the function r: and thus We can further estimate the source term in the pressure equation as follows: . By using the embedding L 6/5 (Ω) ֒→ H −1 (Ω) and the inequality .
On account of Proposition 3.1, the mapping T is well-defined, and, furthermore, a.e. in time, with Λ(t) and F(t) defined in (3.5) and (3.6), respectively; that is, H −1 . By our calculations above, we immediately have is a positive constant that depends on T, R 1 , and R 2 . Furthermore, by relying on (4.4), we obtain Altogether, from (4.5) and the above bounds, we have . Thus, from (4.6), by decreasing R 1 and δ, we can achieve that Xp ≤ R 2 1 . Further, by the embedding H 2 (Ω) ֒→ L ∞ (Ω), we know that Xp , which we can then bound by γ ∈ (0, 1/(2k)) by possibly additionally reducing δ and R 1 . It remains to show that Θ X Θ ≤ R 2 . Proposition 3.1 with f 2 = 0 implies that With the equivalence of the temperature norm and energy (4.3), we have Thus, if we additionally choose R 2 large enough, so that Lemma 4.2. For sufficiently small R 1 and δ, the mapping T is strictly contractive in the topology induced by · X T .
Proof. To prove contractivity, take any (p (2) * ). We introduce the differences p = p (1) − p (2) , Our goal now is to prove that where C is a positive constant that depends on T, R 1 , and R 2 . Observe that (p, Θ) solves the following problem: with the right-hand sides . We can rearrange the acoustic source term f 1 as follows: (2) * t ) := f 11 + f 12 + f 13 and next wish to show that it satisfies Assumption 3.
Further, we know that and By keeping in mind properties (2.5) and (2.6) of the function k, this implies that and thus To obtain a bound on ∇f 11 , we note that Xp .
Using properties (2.6) of the function k, we can bound the first term on the right: Further, we have By using the fact that |k(s)| 1 q 0 , we find Hence, Consequently, from the derived bounds we infer By collecting the derived estimates of separate contributions to f 1 , we arrive at The estimate of ∂ t f 1 L 2 (H −1 ) . Our next task is to estimate ∂ t f 1 L 2 (H −1 ) . As above, we estimate the contributions ∂ t f 1j L 2 (H −1 ) for j = 1, 2, 3 separately. We start by noting that By employing the H −1 estimate stated in (2.12), we then find that (4.16) Hence, we obtain from above .
Finally, we estimate the last term on the right-hand side of (4.18) as (4.21) Using the embedding H 1 (0, t) ֒→ C[0, t], we find that Consequently, we obtain from (4.21), Collecting (4.19), (4.20), and (4.22) results in Finally, by collecting the bounds of separate contributions, we infer that The estimate of f 2 H 1 (L 2 ) . We can bound the source term in the heat equation as follows: (4.24) Since p (j) * t ∈ B for j = 1, 2, we have by the Sobolev embedding Hence, in view of the assumption (2.7), this yields Similarly, using (2.8), we have Plugging (4.25) and (4.26) into (4.24), we obtain The energy bound for the difference equations. Now we can apply the energy results of Proposition 3.1 to system (4.7) by setting Adding the energy estimate for the pressure to the energy bound (3.22) for the temperature (where nowf = f 2 = f 2 ), we obtain (p, Θ) 2 with Λ = Λ(t) defined in (3.5). We have (4.28) We estimate the terms on the right-hand side of (4.28) as follows: using (2.4), we have Further, by using assumption (2.6) we have ).
Next we find that where we have used (2.5) in the last estimate. Using the bound ∇Θ .
Also, we have as above ).
Further, we have the estimate ).
Collecting the above estimates leads to Λ L 1 (0,T ) ≤ C(T, R 1 , R 2 ), where C = C(T, R 1 , R 2 ) is a positive constant that depends on T, R 1 , and R 2 . Finally, taking into account (4.28) and recalling (4.15), (4.23), (4.25), and (4.27) we obtain T (p Thus, by selecting the radius R 1 > 0 sufficiently small, we can guarantee that T is a strict contraction in B. On account of Lemmas 4.1 and 4.2, an application of the contraction mapping theorem implies that there exists a unique (p, Θ) = T (p, Θ) in B which solves the coupled problem.
An application of Gronwall's inequality leads to (p (1) − p (2) , Θ (1) − Θ (2) ) 2 This last inequality yields the desired result, which also implies uniqueness in X T by taking the data to be the same.

Conclusion and Outlook
In this work, we have analyzed the coupled Westervelt-Pennes model of HIFUinduced heating. By relying on the energy analysis of a linearized problem and a subsequent fixed-point argument, we proved the local-in-time well-posedness of this model under the assumption of smooth and (with respect to pressure) small data.
Although our well-posedness result does not require any smallness assumption on T , it is not a global existence result due to the possible dependence of data on the final time. A global existence result would require to show that there is a universal neighborhood of the origin in the topology induced by the norm of the initial conditions for which the solution exists and is bounded uniformly in time. To achieve this, we would need to refine our energy estimates by constructing compensating functionals to capture the dissipation of appropriate components of the norm of the solution. This analysis will be the subject of future research.
We note further that in the energy estimates in Section 3, b must be a positive constant, independent of Θ. To permit more realistic modeling scenarios, future analysis will involve studying the case b = b(Θ) and allowing for (time-or space-) fractional damping in the model.