{'fraction_non_alphanumeric': 0.03754693366708386, 'fraction_numerical': 0.01523643685041084, 'mean_word_length': 4.640883977900552, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photoacoustic imaging taking into account thermodynamic attenuation

In this paper we consider a mathematical model for photoacoustic imaging which takes into account attenuation due to thermodynamic dissipation. The propagation of acoustic (compressional) waves is governed by a scalar wave equation coupled to the heat equation for the excess temperature. We seek to recover the initial acoustic profile from knowledge of acoustic measurements at the boundary. We recognize that this inverse problem is a special case of boundary observability for a thermoelastic system. This leads to the use of control/observability tools to prove the unique and stable recovery of the initial acoustic profile in the weak thermoelastic coupling regime. This approach is constructive, yielding a solvable equation for the unknown acoustic profile. Moreover, the solution to this reconstruction equation can be approximated numerically using the conjugate gradient method. If certain geometrical conditions for the wave speed are satisfied, this approach is well--suited for variable media and for measurements on a subset of the boundary. We also present a numerical implementation of the proposed reconstruction algorithm.

There have been recent efforts to incorporate acoustic attenuation in the modeling of photoacoustic tomography (PAT) and thermoacoustic tomography (TAT). See [29,30,31,32,33,34,35,36,37,38,18,39] and references therein. Most of these results aim at modeling attenuation in the frequency domain to account for dissipation and dispersion. In this paper, however, we adopt a model where the propagation of acoustic waves is thermodynamically coupled to the diffusion of heat. The photoacoustic effect, on which PAT is physically based, consists of two transformations of energy. First, electromagnetic energy is absorbed and transformed into heat. Second, there is a conversion of heat into mechanical energy due to thermal expansivity of the tissues. Concerning this second step, due to the thermodynamic interaction between temperature and pressure, the reverse transformation of energy also takes place. Since heat diffuses, this process attenuates the mechanical energy of the acoustic waves. See [40,Ch. 8] for an introduction to thermoelasticity in biomechanics. We claim that this type of acoustic attenuation should be naturally considered in PAT because PAT itself is based on the thermo-elastic interaction.
We realize that, in mathematical terms, the PAT problem coincides with a problem of boundary observability, which is one of the central concepts of control theory for partial differential equations [41,42,43,44]. We have already employed similar tools to address the PAT problem in an enclosure [26] and other related problems [45,46]. The objective of this paper is to constructively employ the tools of observability for hyperbolic equations together with certain regularity properties of parabolic equations to solve the PAT problem in the presence of thermodynamic attenuation. For the thermoelastic system, there is a series of works on establishing exact, approximate and null controllability or observability estimates [47,48,49,50,51,52,53,54,43]. In these works, however, either the boundary condition or the distribution of control/observation is not of the type we need to model the PAT problem. Therefore, we modify some ideas provided by these references to seek a solution for PAT in the weak coupling regime. Although it might be possible to use Carleman estimates from [53,54,43] to treat the strongly coupled system, we refrain from doing so because the thermoelastic coupling in PAT is known to be relatively weak. Hence, we claim that the results of this paper are sufficient for the nature of PAT in biological tissues.

Mathematical Formulation and Main Results
In this paper we study the photoacoustic tomography problem in the presence thermodynamic attenuation. This is modeled by the linear equations of elasticity coupled with thermal diffusivity [40]. Let Ω ⊂ R n be a smooth bounded connected domain with boundary ∂Ω where n ≥ 2. The propagation of thermoelastic waves in isotropic media is governed by the following system [40,Ch. 8], for the displacement u and where θ denotes the deviation from the reference temperature θ ref . Also, ρ is the mass density, λ and µ are the Lamé coefficients, α denotes the thermal diffusivity and c p is the specific heat at constant pressure. The thermoelastic coupling is given by βK where β is the coefficient of volumetric thermal expansion and K is a constant reference bulk modulus. In soft biological tissue, λ µ which implies that K/(λ + 2µ) ≈ 1. For the moment let us assume that λ, µ and ρ are constants in Ω. Later, we will drop this assumption.
Since photoacoustic imaging is primarily concerned with the compressional waves, we define the pressure p = −(λ + 2µ) div u and the square of the compressional wave speed c 2 = (λ + 2µ)/ρ and proceed to obtain a scalar model for the thermoacoustic waves. Simultaneously, we seek to reveal the strength of the thermoelastic coupling by writing the governing equations in unitless form. Let L be a characteristic length scale of the domain Ω (such as its diameter). Let the characteristic time scale be given by is the square of a reference wave speed. We define the following unitless variables and parameters: Pressurep = p/K and temperatureθ = θ/θ ref .
The unitless coupling parameter > 0 is introduced to analyze the case where is sufficiently small. This is valid for small coefficient of thermal expansion β. The unitless product G = σ is known as the Grüneisen coefficient. Table 1 displays rough estimates for the values of these physical parameters for soft biological tissues. We obtain that 0.05 0.1 and 0.5 σ 10. For notational convenience, we assume that σ = 1. This presents no impediment to the theory as we could easily treat the case σ > 0. At this point, to alleviate the notation, we also drop the caret to denote the unitless quantities. The unitless scalar governing system then becomes, We take (1)-(2) as the starting point for the modeling of this problem. While the equivalence of the thermoelastic system and equations (1)-(2) only holds when λ, µ and ρ are constants, we will consider c(x) to be a sufficiently smooth function in Ω. This follows the common practice of considering a variable wave speed to model heterogeneous media, even when the wave equation is not in divergence form. The system (1)- (2) is augmented by the following initial and boundary conditions, Here, γ : ∂Ω → [0, ∞) denotes the acoustic impedance coefficient at the boundary ∂Ω. Physically, γ = 0 models an acoustically hard surface (such as reflectors) and γ → ∞ approximates an acoustically soft boundary. In general we allow γ to vary on the boundary ∂Ω to model the inhomogeneous nature of an enclosing surface and the interface with sensors or air. The length of the observation window of time is given by τ < ∞ which is defined below. In (4), the symbol ∂ ν denotes the outward normal derivative at the boundary ∂Ω. All the coefficients in (1)-(4) are assumed to be sufficiently smooth.
Concerning the initial conditions, it is common in the modeling of photoacoustic tomography to assume the following. These two expressions are respectively valid when pressure relaxation and thermal diffusion are negligible in the very short lapse of heat deposition from the optical source. These conditions can be achieved in biological tissues by using nanosecond optical pulses [1,2,3,12,5,6,4,55,56]. It turns out that the assumption p 1 = 0 is not needed from the mathematical point of view because it is possible to recover both terms of the initial Cauchy data for the acoustic field. On the other hand, the assumption θ 0 = p 0 is mathematically crucial because it removes an important degree of freedom in the analysis. It would not be possible to recover, in stable manner, an independent initial condition for the thermal field. This is a well-known consequence of the smoothing effect of the heat equation.
In order to consider partial measurements and the presence of acoustic reflectors, we divide the boundary as the disjoint union ∂Ω = Γ ∪ (∂Ω \ Γ) where Γ is the portion where we make observations of the acoustic field. As reviewed in the next section, the forward problem (1)-(4) has a unique solution, and we can define the measurement map given by where p is the solution of (1)-(4) with initial conditions satisfying Assumption 2.1. The goal of the photoacoustic tomography problem is to find the initial profile p 0 from knowledge of Mp 0 . This is a challenging problem with intricate dependencies between the domain Ω, the partial boundary Γ, the wave speed c and the time interval (0, τ ). The admissible dependencies are made precise by a sophisticated assumption of geometric character. Following Bardos, Lebeau and Rauch [42], we assume that our problem enjoys the geometric control condition for the Riemannian manifold (Ω, c −2 δ ij ) with only a portion Γ of the boundary ∂Ω being accessible for observation. We assume that Γ is a smooth open domain relative to ∂Ω. This geometric assumption is the following.
Assumption 2.2 (Geometric Condition). There exists τ < ∞ such that any geodesic ray, originating from any point in Ω at t = 0, eventually reaches Γ in a non-diffractive manner (after possible geometrical reflections on ∂Ω \ Γ) before time t = τ . We refer to τ as the controllability or observability time.
The main theoretical result of this paper is the following theorem.
In addition to this theoretical result, we also propose a convergent iterative reconstruction algorithm. This algorithm is described in Section 4.

Proof of the Main Result
In order to properly analyze the inverse problem we must first state some mathematical properties of the initial boundary value problem (1)-(4). Our guiding references are [57,58,59,60]. We use the following definition of energy for the thermoacoustic system, It is convenient to seek a space that is complete under the energy norm. Since p = const. and θ = const. is a non-trivial solution to the governing equations, the sought space should rule out nonzero constant solutions. This is achieved by considering only solutions such that This is motivated by the fact that the left-hand sides of (6)- (7) are indeed independent of time provided that p and θ solve the governing equations. Therefore, it is only needed to require (6)-(7) at time t = 0. Then the appropriate energy space is given by Notice that H is a closed subspace of H 1 (Ω) × H 0 (Ω) × H 1 (Ω). So it is complete under the norm of H 1 (Ω) × H 0 (Ω) × H 1 (Ω) as well as under the energy norm. The problem (1)-(4) for (p, ∂ t p, θ) is well-posed on the space H and the energy is non-increasing. We state this in the form of a lemma. The proof follows from the standard analysis of partial differential equations and semigroup theory [57,58,59,60].
Now we proceed to prove the main result.

Reconstruction Algorithm
In this section we explicitly recover the initial acoustic profile p 0 in terms of the boundary measurements Mp 0 = p| (0,τ )×Γ . This is accomplished by using the observability properties obtained in the previous section leading to the invertibility of the normal operator (M * M).
In order to obtain an applicable expression for the operator M * , in this section we state the dual or adjoint problem associated with (1)-(4). This is equivalent to constructing the well-known Hilbert Uniqueness Method (HUM) for control of partial differential equations. See [61,41] for details. Throughout, we assume that > 0 is sufficiently small for Theorem 2.3 to apply. This adjoint problem is to find a weak solution (ψ, ξ) for the following IBVP, for a given η ∈ H 0 ((0, τ ) × Γ) (extended as zero on (0, τ ) × ∂Ω \ Γ). Notice that this problem is solved backwards in time with vanishing Cauchy data at time t = τ . Since this problem is well-posed [57], we obtain the following definition.
Now we proceed to describe how the action of C * can be approximated using the conjugate gradient method. See [62,Sect. 4.6] for a standard description of the conjugate gradient method in a Hilbert space setting. For sake of completeness, we describe the inversion of a generic equation (M * M)φ = ζ. Let φ 0 be an initial guess for the true solution φ * . Define r 0 = ζ − (M * M)φ 0 as the initial residue and s 0 = r 0 . For k ≥ 0, define .
Since the operator (M * M) : H 1 (Ω) → H 1 (Ω) is bounded and coercive, then there are positive constants m and M such that . The conjugate gradient iterates can be shown to converge as follows (see [62,Sect. 4.6] and references therein), Notice that at each iteration, one must apply the operator (M * M) which amounts to solve the problem (1)-(4) (under Assumption 2.1) followed by solving the adjoint problem (11)- (14). In practice, this can be approximated using numerical methods for PDEs. However, depending on the method of choice, there are intrinsic complications that may prevent a convergence estimate such as (16) from being satisfied in the limit as the discretization is refined. We shall not elaborate any further on these complications as they lie outside of the scope of this paper. For details on these numerical issues we refer to [41,63,64,65,66] and references therein. In this paper, we adopted the two-grid approach described in [65,41] using second order finite difference methods.

Numerical Results
Now we present some numerical results to illustrate the performance of the reconstruction algorithm described in Section 4. We implemented a numerical solver for the governing system (1)-(4) and its adjoint (11)-(14) based on second order finite differences. To avoid spurious numerical instabilities, we adopted the two-grid approach described in [65,41]. We worked in R 2 where the domain Ω was taken as the unit-square. The initial profile p 0 corresponds the Shepp-Logan phantom.
We present two examples. One with constant wave speed c(x) ≡ 1, and the other with variable wave speed c = c(x) defined below. In both cases, we used the following parameters: impedance γ(x) = c −1 (x) over the boundary of Ω, thermal diffusivity α = 0.01, and coupling parameter = 0.1. The observability time was chosen to be τ = 2 which is enough for more than 99% of the energy contained in the initial profile to dissipate or leave the domain through the boundary when the wave speed is constant.
We shall compare the results from the proposed algorithm against the results from purely acoustic time-reversal. The latter is accomplished by producing measurements using the thermoacoustic forward solver M, and then back-propagating the boundary measurements in a purely acoustic medium ( = 0), that is, by ignoring the thermodynamic attenuation. See details in [26,16] for the purely acoustic time-reversal approach. The acoustic time-reversal is approximated using the same finite difference method. The initial guess for the conjugate gradient algorithm is the approximate solution obtained from the purely acoustic time-reversal algorithm. Although the proposed reconstruction algorithm has been described in the H 1 (Ω) setting, a similar study could be performed in the H 0 (Ω) setting where the inner-products in the conjugate gradient algorithm would need to be understood appropriately. In this section, we present results from the implementation both in the H 1 (Ω) and H 0 (Ω) formulations.

Constant wave speed
For the first example where c ≡ 1, Figure 1 displays the exact initial profile and the reconstructions. The relative errors in the H 1 (Ω) and H 0 (Ω) formulations are reported in Table 2 for the first few iterations of the conjugate gradient algorithm. We notice that by ignoring the thermodynamic attenuation in the purely acoustic timereversal reconstruction (Iter = 0), the edges in the Shepp-Logan phantom are blurred considerably. Some of the sharpness is recovered by accounting for the attenuation in the proposed algorithm even after a single iteration. Table 2.
Constant wave speed example. Relative error at each iteration of the conjugate gradient method described in Section 4. Iter = 0 corresponds to the initial guess given by a purely acoustic time-reversal algorithm.

Variable wave speed
For the second example we have selected a variable wave speed defined as a layer of higher speed surrounding the smaller ellipses in the Shepp-Logan phantom. The actual profile is illustrated in the top-right panel of Figure 2. The relative errors in the H 1 (Ω) and H 0 (Ω) formulations are reported in Table 3 for the first few iterations of the conjugate gradient algorithm. Again, as shown in the lower panels of Figure 2, we see great improvements over the purely acoustic time-reversal reconstruction. We highlight the ability in capturing the jump discontinuities and the reduction of the artifacts introduced by ignoring the attenuation.

Conclusion
We have presented a PAT/TAT model based on thermoelasticity. The thermoelastic coupling accounts for how pressure changes can induce temperture changes in a body and vice versa. The coupling between temperture and deformation is a fundamental feature of PAT/TAT. The current literature dealing with PAT/TAT only considers one side of the thermoelastic interaction (the photoacoustic effect). By considering both effects simultaneously we account for a natural attenuation phenomenon of the acoustic signal. We related the thermoelastic model of PAT/TAT with boundary observavility for the thermoelastic system. We showed uniqueness and stability of recovering the initial profile of the pressure from boundary data. The initial profile was reconstructed using a conjugate gradient method. The recovery analysis of the inital wave profile works under a geometric control assumption on the wave speed. The performance of the reconstruction method is very similar in both the constant and variable wave speed scenarios.
The coupling uniteless parameter of the thermoelastic model is approximately between 0.05 and 0.1 as obtained from Table 1. Theorem 2.3 requires to be sufficiently small. Given that in PAT/TAT the thermodynamic interaction is small, such a condition on is reasonable. Nonethelesss, it might be possible to remove this condition by using Carleman estimates for the coupled thermoelastic system (e.g., [53,54,43]).