{'fraction_non_alphanumeric': 0.03610393339075509, 'fraction_numerical': 0.01661642262417456, 'mean_word_length': 4.614547652629458, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solvability for Photoacoustic Imaging with Idealized Piezoelectric Sensors

Most reconstruction algorithms for photoacoustic imaging assume that the pressure field is measured by ultrasound sensors placed on a detection surface. However, such sensors do not measure pressure exactly due to their non-uniform directional and frequency responses, and resolution limitations. This is the case for piezoelectric sensors that are commonly employed for acoustics-based biomedical imaging. In this paper, using the method of matched asymptotic expansions and the basic constitutive relations for piezoelectricity, we propose a simple mathematical model for piezoelectric transducers. The approach simultaneously models how the pressure waves induce the piezoelectric measurements and how the presence of the sensors affects the pressure waves. Using this model, we analyze whether the data gathered by piezoelectric sensors leads to the solvability of the photoacoustic imaging problem. We conclude that this imaging problem is well-posed in certain normed spaces and under a geometric assumption. We also propose an iterative reconstruction algorithm that incorporates the model for piezoelectric measurements. Numerical implementation of the reconstruction algorithm is presented.


Introduction
Photoacoustic or thermoacoustic tomography is a nonionizing imaging modality designed to advantageously combine the high contrast of optical absorption with the high resolution from broadband ultrasound waves. The imaging of optical absorption reveals important functional and pathological information about biological tissues [39,8,35,37].
One of the open challenges concerning photoacoustic inversion is the incorporation of realistic models for acoustic measurements. This need for modeling the physics of ultrasound sensors has been recognized in [24,20,38,4]. It has been claimed that ultrasound measurements can be described as a linear combination of the pressure field and its normal derivative at the boundary. With that motivation, Dreier and Haltmeier [15] recently established explicit formulas for the inversion of the two-dimensional wave equation from Neumann boundary data for circular and elliptical domains. In a related effort, Zangerl, Moon and Haltmeier [48] derived Fourier-based reconstruction formulas for the spherical detection geometry from knowledge of Robin boundary data.
Most other reconstruction algorithms assume that the pressure field (Dirichlet data) is measured on the detection boundary. This assumption is not satisfied in practice because acoustic transducers do not measure pressure directly. They measure a surrogate for pressure that depends on the actual transducer mechanism. The most common mechanisms for ultrasound applications are based on the piezoelectric effect [44,40,41], on Fabry-Perot interferometry [9,14,50,22] or on fiber-optic refractometry [43,42,45]. In [3] we formulated a model specifically tailored to the Fabry-Perot sensor design. In this paper, we derive a similar model for piezoelectric sensors and analyze the well-posedness arXiv:2002.09929v1 [math.AP] 23 Feb 2020 of photoacoustic imaging with such measurements. In order to attain a balance between accuracy and simplicity, the model we develop here is based on the following underlying idealizations: (a) The sensors are treated as point-like detectors. Hence, we do not account for resolution limitations due to the finite size of sensing elements. See [47,23,36,29,40,41] for investigations concerning this issue. We also assume that the domain to be imaged is fully enclosed by the detection surface. (b) We assume that in each detector, the piezoelectric film is flat and its thickness is small in comparison to the wavelengths under consideration. In practice, industrial processes can manufacture piezoelectric films with thickness 30 -100 µm approximately [19,30,32,44,28]. (c) Although the sensors contain elastic materials that may support shear waves, our analysis is entirely based on compressional waves governed by the scalar wave equation. (d) For the piezoelectric film, the poling direction is along its thickness, and the piezoelectric properties are transversely isotropic in the plane perpendicular to the poling direction. The sensing film is mechanically isotropic. (e) Sensors may have complex structures, including a casing for structural integrity, electrodes, bonding layers and multiple paddings designed to match the mechanical impedance of the acoustic medium [11, 34,33]. However, we assume a simple design consisting of the piezoelectric film, sandwiched by electrode foils of negligible thickness, mounted on a much thicker backing layer. This follows models described in [19,30,33].
An illustration of the idealized setup is shown in Figure 1. The acoustic domain, denoted by Ω, contains soft tissue with variable density ρ and variable wave speed c. The piezoelectric film Ω s has a small uniform thickness > 0, constant density ρ p and constant wave speed c p . The thick backing layer Ω b has constant density ρ b and constant wave speed c b . The interface between the acoustic domain and the piezoelectric film is denoted Γ. The interface between the piezoelectric film and the backing layer is denoted Γ .
In section 2 we derive a model for the transduction from pressure to electrical voltage which is the physical quantity acquired by the piezoelectric sensors. In section 3 we derive an effective boundary condition for the transmission of waves from the acoustic medium of interest into the piezoelectric sensor. This effective transmission condition accounts for the influence that the sensor exerts on the acoustic waves. Using the cou- Fig. 1: Acoustic domain Ω with density ρ and wave speed c. Piezoelectric material Ω p of uniform thickness , density ρ p and wave speed c p . The thick backing layer Ω b has density ρ b and wave speed c b . The interface between the acoustic domain and the piezoelectric film is denoted by Γ. The interface between the piezoelectric film and the backing layer is denoted by Γ .
pled models for piezoelectric measurements and wave propagation, in section 4 we define the forward problem associated with photoacoustic imaging. Then in section 5 we state and prove the solvability of the imaging problem with piezoelectric measurements. A reconstruction algorithm is proposed in section 6 where some numerical simulations are presented. The conclusions follow in section 7.

Model for Piezoelectric Measurements
We start the modeling of the piezoelectric measurements from the basic constitutive relations for both piezoelectric and mechanical variables. Since the sensing material is mechanically isotropic, the stress tensor σ is related to the strain tensor s as follows, where the direction along the thickness of the piezoelectric film is denoted as the 3-axis, and the 1-axis and 2axis are the transverse plane. Here δ ij is the Kronecker delta, and λ and µ are the first and second Lamé coefficients. The equation of mechanical motion is where u is the material displacement vector. For irrotational deformations, i.e. in the absence of shear stress, the above equation can be simplified in order to relate the particle displacement u to the pressure p p in the piezoelectric film, where the pressure p p is defined as Combining (3) and (4), we find that the pressure field p p satisfies the wave equation, where the wave speed c p is defined by c 2 p = (λ + 2µ) /ρ p . The piezoelectric transducer measures the electrical voltage V across the piezoelectric film generated by the mechanical deformation due to the transmitted acoustic waves. We proceed to derive the mathematical relationship between the voltage V and the pressure p p in the piezoelectric material. Our guiding references are [34,Ch. 5], [33,Ch. 5] and [30]. Under small perturbations of field conditions, the linearized constitutive relation for the piezoelectric effect is the following where D is the electric displacement (electric charge per area), E is an externally applied electric field (voltage per length), σ is the stress tensor (force per area). The piezoelectric properties are the dielectric permittivity tensor ε (capacitance per length) and the piezoelectric tensor d (electric charge per force). In the absence of shear stress and of an external electric field, the normal electric displacement D 3 is given by As assumed above, the piezoelectric tensor in transversely isotropic, which allows us to simplify the notation as follows d ⊥ = d 31 = d 32 , d = d 33 , and D = D 3 . Combining the constitutive relations (7) and (1) we obtain the electric displacement in terms of the strain, where e ⊥ = 2d ⊥ (λ + µ)+dλ and e = d (λ + 2µ)+2d ⊥ λ. As a consequence, using the definition of strain s in terms of the displacement u, we obtain where ∂ n represents the derivative along the 3-axis, i.e., normal to the piezoelectric film. Now we take two timederivatives of (9) and combine with (3)-(4) to obtain Plugging (5) into (10) we obtain We have expressed ∆ = ∂ 2 n + ∆ ⊥ where ∆ ⊥ represents the surface Laplacian on the transverse plane.
The voltage V generated across the film by the electric displacement D, which is obtained from Gauss law where is the thickness of the piezoelectric film, ε is its dielectric permittivity, and dn is the differential across the film. Combining (11) and (12), we obtain our model for the piezoelectric measurements with vanishing initial state. Here the symbol ∝ denotes equality up to a multiplicative constant (which is typically estimated through experimental calibration). In (13), it is assumed that the pressure field is constant across the piezoelectric film. This assumption is rigorously justified in the next section. The model (13) for the piezoelectric sensor is qualitatively similar to the model for the Fabry-Perot transducer proposed in [3] in spite of the completely different physical principles from which they are derived. The symbol κ appearing in the model (13) is a unitless coefficient defined by the elastic and piezoelectric properties of the sensing film where we have expressed λ = ρ p c 2 p ν/(1 − ν) and µ = ρ p c 2 p (1−2ν)/(2−2ν) in terms of Poisson's ratio ν to obtain the last equality. Common values for all these physical parameters are shown in Table 1 for polyvinylidene fluoride (PVDF) piezoelectric sensors.
We note from (13) that a theoretically perfect transduction from pressure to voltage would be attained if the coefficient κ = 0. However, due to the nature of the poling processes employed to manufacture these piezoelectric materials, the coefficients d and d ⊥ have opposite signs and generally |d|> |d ⊥ |. This implies that 1 < (1 − d ⊥ /d) < 2. Hence, in order for κ = 0, the Poisson's ratio would have to be ν = 0.5 which requires the piezoelectric material to be incompressible. In practice, Poisson's ratio for PVDF films ranges from 0.2 to 0.4 approximately. We note that κ ranges from 0.3 to 1.5, for the realistic range of values for the Poisson's ratio ν and the piezoelectric ratio d ⊥ /d displayed in Table 1.  [19,30,11,49,44,10,33,12].

Effective Model for Wave Propagation
Typically, the piezoelectric film and the backing layer are acoustically more rigid and heavier than the biological medium of interest. Therefore, the presence of the sensors induces partial reflections of the waves.
Here we seek to model how the sensors exert influence on the pressure waves. This model takes the form of an effective impedance boundary condition that replaces the more involved transmission process for waves traveling from the acoustic domain Ω, through the piezoelectric film Ω p and into the backing layer Ω b . We assume that the pressure field p b in the backing layer is outgoing which translates into satisfying a radiation condition of the form, where H is the mean curvature of the surface Γ. See [6,1] for a derivation. As in [3], we make some geometric assumptions about the domain Ω p occupied by the piezoelectric film. We let Ω p = {y ∈ Ω c : 0 < dist(y, Γ) < }. For sufficiently small , the domain Ω p can be expressed as a family of parallel surfaces parametrized by 0 < s < and defined by Γ s = {y = x + sn(x) : x ∈ Γ} where n(x) is the normal vector at x ∈ Γ. For smooth Γ and sufficiently small , the surfaces Γ s are well-defined, smooth and mutually disjoint. Each point y ∈ Ω p can be uniquely represented in the form y = x + sn(x) for x ∈ Γ and 0 < s < . In addition, the normal vector at y ∈ Γ s coincides with the normal vector at x ∈ Γ. See details concerning parallel surfaces in [26,Sect. 6.2].
The transmission of the pressure field from the acoustic domain into the piezoelectric film is governed by the following transmission conditions at the interface Γ, where p and p p are the pressure in the acoustic medium and piezoelectric film, respectively. The first condition in (16) ensures the continuity of the pressure field. The second condition in (16) ensures the continuity of particle motion in the normal direction. Similar transmission conditions hold at the interface Γ , where p b is the pressure in the backing layer. The pressures p, p p and p b satisfy the wave equation with respective wave speeds c, c p and c b . Now, we proceed to use the method of matched asymptotic expansions to derive an effective model for the interplay between the pressure fields and the piezoelectric sensor. First we consider the formal asymptotic expansions for the pressure fields, and introduce a change of variable in order to extract the effect of the piezoelectric film thickness , The boundary value problem for the pressure field p p in the piezoelectric film, is governed by the wave equation (5) and the transmission conditions (16)- (17) can be recast in terms of ζ and terms with same powers of are gathered to obtain the following cases.
O( 0 )-terms: for ζ ∈ (0, 1), p 0 = p 0 p and ∂ ζ p 0 p = 0, at ζ = 0, p 0 p = p 0 b and ∂ ζ p 0 p = 0, at ζ = 1, which imply that p 0 p is constant as a function of ζ and that the first effective transmission condition is that O( 1 )-terms: for ζ ∈ (0, 1), which imply that ∂ ζ p 1 p is constant as a function of ζ, and that the second effective transmission condition is Combining (15) and (20)- (21), we obtain closedform effective governing equations for the leading order term p 0 of the acoustic field in the domain Ω, Similar models for photoacoustics are studied in [4,5,3,2] As an example for the response of the piezoelectric sensor design, we can analyze its behavior for plane waves and for a flat boundary Γ. Both the boundary value problem (22)-(23) and the model for the measurements (13) play an important role in this analysis. A plane wave of the form p inc = e i(x·k−ωt) propagating in the direction of k, induces a reflection governed by (23). The total pressure field p is the superposition of the incident and reflected wave, where R is the reflection coefficient, k r is the reflection wavenumber satisfying |k|= |k r |= ω/c and n · k r = −n · k, where n is the outward normal on Γ. We can write n · k = |k|cos θ where θ is the angle of incidence. Plugging (24) into (23) and neglecting the O( ) terms, we find that the reflection coefficient satisfies After plugging (24)-(25) into the model (13) and evaluating at the origin x = 0, we find that the piezoelectric measurements satisfy the following directivity pattern where κ is given by (14). Figure 2 displays the directional response (26) in dB as a function of the incidence angle θ and piezoelectric coefficient κ over a realistic range of values shown in Table 1. We observe that for values of κ > c 2 /c 2 p , a critical angle appears. For incidence at this critical angle, vanishing measurements are obtained by the piezoelectric sensor design. This critical angle is given by θ cr = arcsin c c −1 p κ −1/2 .

The Forward Problem
Now we proceed to define the inverse problem of photoacoustic imaging in terms of the wave propagation model (22)-(23) and the model for piezoelectric measurements (13). We neglect higher order terms O( ) studied in the previous section, so that the pressure field p is assumed to satisfy the following initial value problem, where 0 < T < ∞ is the measurement time to be determined later. Recall that the underlying assumption concerning media properties are that c is bounded from below and above, and is smooth in Ω, and that c p , c b , ρ p , ρ b , ρ and κ are constants. The forward mapping, which we seek to invert for photoacoustic imaging, is given by where, according to the piezoelectric model (13), the measured electric voltage V satisfies for the pressure field p evolving according to (27) from the initial condition f . We work with the standard Sobolev spaces based on square-integrable functions defined on the domain Ω or the boundary (0, T ) × Γ. The associated inner product extends as the duality pairing between functionals and functions. For the Sobolev space H 0 (Ω), the inner product is weighted by c −2 so that the differential operator c 2 ∆ is formally self-adjoint. The well-posedness in Sobolev spaces of the initial value problem (27) is a well-established result [18,27].

The Inverse Problem
The inverse problem associated with photoacoustic imaging is the following: Given the voltage measurements V modeled by (29) on Γ × (0, T ), induced by the pressure field p satisfying (27), find the unknown initial condition f . The solvability of this inverse problem depends on the geometry of the domain Ω, the profile of the wave speed c and the time T < ∞. These conditions are made precise in the following assumption, known as the geometric control condition or a nontrapping condition for the geodesic flow [7,21].

Assumption 1 (Nontrapping condition)
Let Ω be a simply connected bounded domain with smooth boundary Γ. Assume there exists T o < ∞ such that any geodesic ray of the manifold (Ω, c −2 dx 2 ), originating from any point in Ω at time t = 0, reaches the boundary Γ at a nondiffractive point before t = T o .
With this assumption in place, we can state the main result of the paper in the form of a theorem. ) is injective, that is, the photoacoustic imaging problem is uniquely solvable. Moreover, the following stability estimate, holds for some constant C > 0.
We wish to make some comments before we proceed with the proof. Notice in (30) that we are only able to dominate f in the norm of H 0 (Ω) (rather than in the norm of its stated space H 1 0 (Ω)) with the measured data V in the norm of H 1 ([0, T ]; H 0 (Γ)). By contrast, when the Dirichlet data is measured on [0, T ] × Γ, then the imaging operator (left inverse of F) enjoys stability estimates as a mapping from H 0 ([0, T ] × Γ) to H 0 (Ω), or from H 1 ([0, T ] × Γ) to H 1 0 (Ω). See [7,31,4] for details. Hence, there is an apparent loss of stability due to the nature of the piezoelectric measurement model (29). The double time-integration needed to invert the left-hand side of (29a) does not fully restore the regularity lost by the application of the hyperbolic differential operator on the right-hand side of (29a). This is a well-known property concerning regularity of hyperbolic equations [18]. We also note that Theorem 1 is slightly different from what is presented in [3] where the imaging operator was shown to satisfy a stability estimate as a mapping from H 0 ([0, T ]; H 1 (Γ)) to H 0 (Ω). Hence, formally, there is a mild loss of stability of one degree either in space or in time, but not both. Now, it is convenient to define the following operation where p and V satisfy (29) and p(0) = f ∈ H 1 0 (Ω). Then it follows that u solves the following initial value problem, The following lemma is a well-established result. See holds for some constant C > 0.
Using this lemma, we proceed to prove the main theoretical result of the paper. In what follows, the generic constant C > 0 changes from inequality to inequality, but it does not depend on f , p or V .

Numerical Simulations
Now we propose and numerically implement a reconstruction algorithm to solve the PAT problem at the discrete level. The reconstructions presented here are based on the Landweber iterative method [17, Ch. 6] accelerated using Nesterov's approach. The iterative process is defined in Algorithm 1. In this algorithms, the parameter γ is known as the relaxation factor which guarantees the stability and convergence of the iterations. The parameter µ, known as the momentum factor, allows for the acceleration of the convergence. In the presence of noise, the maximum number K of iterations is chosen according to some stopping or regularization rule such as a discrepancy principle. For the Landweber method, it is necessary to evaluate the adjoint F * of the forward operator F. This evaluation amounts to solve the following final boundary value problem, where η solves in order to define the adjoint mapping as Algorithm 1 Accelerated Landweber iteration Set K, 0 < γ < 2 F −2 and 0 ≤ µ < 1. Initial guesses u 0 = F * V and v 0 = u 0 for k = 1, 2, ..., Both the forward map F and its adjoint F * are discretized using a piecewise linear finite element method (FEM) in space and the explicit Newmark method for the time stepping. The discretization parameters are chosen to satisfy the CFL stability condition. The FEM is implemented on triangulations of the domain Ω. For the numerical simulations, we have non-dimensionalized the physical parameters displayed in Table 1  Coarse mesh for numerical approximation of the F and F * mappings. B: Exact profile to be reconstructed showing the brain vasculature imaged with MRI technology [25]. C: Synthetic piezoelectric measurements as modeled by (29). c b = 1.0. The piezoelectric-elastic coupling coefficient κ = 0.9. These non-dimensional parameters are consistent with the ranges of their dimensional counterparts listed in Table 1. Figure 3 displays a coarse mesh used for the FEM, the exact pressure profile to be reconstructed and the boundary measurements. These measurements were synthetically generated by applying the discrete version of the forward operator F using the aforementioned numerical method for the wave equation. The FEM for the reconstruction procedure has 109,762 degrees of freedom and 6,400 time steps covered the time window for T = 2. The mesh employed to generate the measurements was more refined, with mesh size approximately half of the mesh size employed in the reconstruction steps, and the data was down-sampled to the reconstruction mesh using linear interpolation. The performance of the Algorithm 1 for various values of the momentum factor µ is displayed in Figure  4. Significant improvements are observed for increasing values of µ. For instance, in order to reach below 1% relative error, the original Landweber method (µ = 0.0) takes 36 iterations, whereas the accelerated method with µ = 0.6 takes 12 iterations. For values µ > 0.7, some instability is observed. The error profiles for the initial guess (back-projection) and for the 12 th iterate are shown in Figure 5. We observe that after a few iterations, the algorithm is able to recover most of the sharp features of the image. The effect of noise added to the piezoelectric measurements is also explored here. We have considered three types of noise: white, pink and red. These are characterized by increasing autocorrelation distances or equivalently faster decay of their spectral power. White or Gaussian noise has equal power spectral density (PSD) at different frequencies. Pink noise has a PSD decaying as ω −1 as ω → ∞. Red or Brownian noise has a PSD decaying as ω −2 . Figure 6 displays samples from these three types of noise. Figure 6 also displays the relative error as a function of the iteration number and increasing levels of noise for the three types of noise. We observe that the implemented reconstruction method can handle best the white noise. The error for the pink noise case is greater. And the reconstruction error for the red noise case is the greatest of the three types of noise.

Conclusions
We have proposed a mathematical model for the acoustic measurements transduced by piezoelectric sensors. This model, taking the form (13), is highly idealized as described in the Introduction. This idealization allows us to attain a balance between correctness and simplicity in order to analyze the properties of the model (see section 3) and the solvability of the associated photoacoustic problem (see sections [4][5]. The validity of the model is limited to small values for the thickness of the piezoelectric film with respect to the wavelength of the acoustic waves. Advanced industrial processes are able to manufacture piezoelectric films with thickness in the range 30 -100 µm approximately. As shown in Table  1, the wave speed in PVDF materials ranges from 1300 -2300 m/s. Hence, we expect our model to be valid for frequencies 13 MHz.
The directional response for plane waves derived in section 3 and displayed in Figure 2 matches the experimental measurements carried out by other researchers [44,40,41,11]. Specifically, design considerations concerning the mechanical and electrical properties of the piezoelectric film and backing layer play a role in the appearance of critical angles where the sensitivity of the sensor vanishes. See Figure 2. The incorporation of sensor response into reconstruction algorithms has been highlighted as one of the challenges associated with photoacoustic imaging [38,13,16,36,46] and partially investigated in [20,4,48,15,3]. Our present paper represents a novel contribution, tailored to piezoelectric sensors, to this research effort.
From the theoretical perspective as stated in Theorem 1, we conclude that the photoacoustic imaging problem is solvable for piezoelectric measurements. However, some stability is lost compared to measuring directly the Dirichlet data. Based on these theoretical properties of the inverse problem, we have also proposed a reconstruction iterative Algorithm 1 based on an accelerated Landweber method. We implemented some proof-of-concept synthetic simulations with and without noise. Robustness in the presence of noise of different spectral characteristics is observed.
One limitation of the proposed method relates to the numerical characteristics of the FEM and Newmark time-stepping. Notice in Figure 4 that the error initially decays exponentially as the theory for the Landweber method predicts. However, as the iterations continue, the error eventually stagnates. This stagnation may be attributed to the numerical error introduced by the scheme employed to approximate the action of the forward operator F and its adjoint F * . Among other issues, the discrete version of F * may not be an exact adjoint for the discrete version of F. Also, the numerical scheme to solve the wave equation suffers from numerical dispersion. Different frequency components of the initial pressure profile travel at different group velocity towards the detection boundary. As a consequence, the discrete version of F * F loses coercivity (becomes illconditioned) and the reconstructed images suffer from aberration. These complications are not mitigated by refinement of the FEM mesh. Remedies for this phenomenon have been investigated, including regulariza-tion, two-grid methods and numerical schemes for the wave equation with lower dispersion. See [21, Sect. 6.8-6.10], [51] and references therein. However, these improvements fall outside of the scope of this paper.
Other possible areas for future research include efforts to relax the idealizations listed in the Introduction. From the mathematical point of view, it seems plausible to model the resolution limits due to the finite size of detectors, the effect of partial measurements, and the contribution from shear waves supported by the piezoelectric film and backing layer.