I. INTRODUCTION

P HOTOACOUSTIC tomography is a non-ionizing 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 [1]- [4].

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 [5]- [8]. 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 [9] 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 [10] derived Fourier-based reconstruction formulas for the spherical detection geometry from knowledge of Robin boundary data.

Most other reconstruction algorithm assume that waves propagate freely across the detection boundary and that the pressure field (Dirichlet data) is measured exactly. These assumptions are not satisfied in practice. The pressure wave are S. Acosta, Department of Pediatrics, Baylor College of Medicine and Predictive Analytics Laboratory, Texas Children's Hospital, Houston TX, USA. (sebastian.acosta@bcm.edu) Manuscript received date; revised date. This work was partially supported by NSF grant DMS-1712725.

affected by the presence of the sensors and the acoustic sensors 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 [11]- [13], on Fabry-Perot interferometry [14]- [17] or on fiber-optic refractometry [18]- [20]. In [21] 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 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 [12], [13], [22]- [25] 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 [11], [26]- [29]. (c) Although the sensors contain elastic materials that may support shear waves, our analysis is valid for 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 [30]- [32]. 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 [26], [27], [32].

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 Ω p 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 II we derive a model for the transduction from pressure to electrical voltage which is the physical quantity acquired by the piezoelectric sensors. In section III 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 coupled models for piezoelectric measurements and wave propagation, in section IV we define the forward problem associated with photoacoustic imaging. Then in section V we state and prove the solvability of the imaging problem with piezoelectric measurements. A reconstruction algorithm is proposed in section VI where some numerical simulations are presented. The conclusions follow in section VII.

II. 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 3 × 3 symmetric stress tensor σ is related to the 3 × 3 symmetric 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 2-axis 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 3 × 1 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 [31,Ch. 5], [32,Ch. 5] and [27]. Under small perturbations of field conditions, the linearized constitutive relation for the piezoelectric effect is the following

where D is the 3 × 1 electric displacement vector (electric charge per area), E is an externally applied 3 × 1 electric field (voltage per length) and ε is the 3 × 3 dielectric permittivity tensor (capacitance per length). Following [27], it is convenient to express the symmetric stress tensor σ (force per area) as a 6 × 1 vector and the piezoelectric tensor d (electric charge per force) as a 3 × 6 matrix,

In the absence of an external electric field in (6), the normal electric displacement D 3 is given by

As assumed above, the piezoelectric tensor is transversely isotropic, which allows us to simplify the notation as follows

Combining the constitutive relations ( 1) and ( 7) 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 is the normal vector on Γ and ∂ n represents the derivative along the normal direction or 3-axis. Now we take two time-derivatives of (9) and combine with (3)-( 4) to obtain

Express ∆ = ∂ 2 n + ∆ ⊥ where ∆ ⊥ represents the surface Laplacian on the transverse plane, recall that c 2 p = (λ+2µ)/ρ p and solve for ∂ 2 n p p in (5) to plug it into (10) to obtain

By definition of voltage as an electric potential, the voltage V generated across the piezoelectric film by the 3-component of the electric displacement D satisfies

where ε 33 is the dielectric permittivity along the 3-axis.

Integrating (12) across the piezoelectric film and combining the result with (11), 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 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 I 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

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 I.

III. 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 TABLE I: Estimates for the physical parameters of PVDF piezoelectric sensors [11], [26], [27], [30], [32]- [35].

Parameter

Value Units 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,

PVDF

where H is the mean curvature of the surface Γ. See [36], [37] for a derivation.

As in [21], 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 < z < ǫ and defined by Γ z = {y = x + zn(x) : x ∈ Γ} where n(x) is the normal vector at x ∈ Γ. For smooth Γ and sufficiently small ǫ, the surfaces Γ z are welldefined, smooth and mutually disjoint. Each point y ∈ Ω p can be uniquely represented in the form y = x + zn(x) for x ∈ Γ and 0 < z < ǫ. In addition, the normal vector at y ∈ Γ z coincides with the normal vector at x ∈ Γ. See details concerning parallel surfaces in [38,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. For an introduction to this method, see [39]. 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, governed by the wave equation ( 5) and the transmission conditions ( 16)-( 17

b) O(ǫ 1 )-terms:

p is constant as a function of ζ and that the second effective transmission condition is

Combining ( 15) and ( 20)-( 21), we obtain closed-form effective governing equations for the leading order term p 0 of the acoustic field in the domain Ω,

) Similar models for photoacoustics are studied in [8], [21], [40], [41].

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

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 decibels as a function of the incidence angle θ and piezoelectric coefficient κ over a realistic range of values shown in Table I. 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

IV. THE FORWARD PROBLEM

Now we proceed to define the forward 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 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 . The mapping (28), which we seek to invert for photoacoustic imaging, encodes the physical principles for acoustic wave propagation from the unknown pressure profile f to the measured electrical voltage V generated by the piezoelectric sensors.

We work with the standard Sobolev spaces based on squareintegrable 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 wellposedness in Sobolev spaces of the initial value problem ( 27) is a well-established result [42], [43].

V. 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. We work with the manifold Ω endowed with the Riemannian metric c -2 dx 2 . See [44], [45] for details.

Assumption 1 (Nontrapping condition):

Let Ω be a simply connected bounded domain with smooth boundary Γ. Assume there exists T o < ∞ such that any (unit speed) 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.

Theorem 1: Under the Assumption 1 for the manifold (Ω, c -2 dx 2 ) and time T > T o , the forward mapping F :

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 [8], [44], [46] 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 [42]. We also note that Theorem 1 is slightly different from what is presented in [21] 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 [42, §7.2, Thms. 3-5] or [43, Ch. 3, §8, Thm. 8.1] for details. Lemma 1: Let u solve (33).

holds for some constant C > 0.

The definition of the Bochner spaces H 1 ([0, T ]; H 0 (Γ)) and C k ([0, T ]; H 1-k (Γ)) can be found in [42], [43]. Using Lemma 1, 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 .

Proof of Theorem 1: Under Assumption 1 for the manifold (Ω, c -2 dx 2 ) and time T > T o , observability of waves from the boundary [44], [45] yields that

for some constant C = C(Ω, c, T ). Now, from the definition (32) of u and the stability estimate in Lemma 1 for k = 1, we obtain that

Since the norm of C 1 ([0, T ]; H 0 (Γ)) dominates the norm of H 0 ([0, T ]×Γ), combining ( 35) and ( 36) we obtain the desired result (30).

VI. 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 [47,Ch. 6] with Nesterov's acceleration. See Stefanov and Yang [48] for an excellent analysis of the Landweber method for a similar problem. The iterative process is defined in Algorithm 1. In brief, the Landweber method is based on inverting (28) by solving

where the normal operator (F * F) is positive definite and γ > 0 is chosen small enough so that the spectrum of (I -γF * F) is contained in (-1, 1) making it a contraction. The parameter γ is known as the relaxation factor. The parameter µ, known as the momentum factor, allows for the acceleration of the convergence. Unfortunately, it is hard to choose γ and µ optimally. We resort to trial and error to set them satisfactorily. In the presence of noise, the number K of iterations is chosen according to a regularization rule.

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

It can be shown, using straight-forward integration by parts, that indeed F * is the adjoint of F or equivalently that

for all sufficiently smooth f and ψ, where V = F(f ) according to (28) and ∂ t ϕ| t=0 = F * (ψ) according to (39). The brackets •, • X denote the inner product of the Hilbert space H 0 (X) which extends as the duality pairing between the Sobolev functionals and functions.

return u K Fig. 3: A: Coarse mesh for the FEM method. B: Exact profile to be reconstructed showing the brain vasculature imaged with MRI technology [49]. C: Synthetic piezoelectric measurements as modeled by (28).

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 nondimensionalized the physical parameters displayed in Table I in order to have c = 1, diam(Ω) = 2 and ρ = 1. The final time T = 2. The non-dimensional parameters of the piezoelectric film have been chosen as follows ρ p = 1.5, c p = 1.0. The parameters of the backing layer are ρ b = 2.0 and 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 I.

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. In order to visualize the impact of improperly modeling the piezoelectric measurements, we have implemented two reconstructions of the initial acoustic profile. In both case, Fig. 5: Error profiles for synthetic measurements obtained using the piezoelectric model (28). A: Reconstruction using an interpretation of the measurements as piezoelectric data. The relative error is 0.54%. B: Reconstruction using a naive interpretation of the measurements as Dirichlet data. The relative error is 6.43%.

the same synthetic measurements are used. For the first reconstruction, we properly interpret the given measurements as generated by the piezoelectric model (28) and carry out the reconstruction using Algorithm 1. After 50 iterations we obtain a relative error of 0.54%. For the second reconstruction, we improperly interpret the measurements as the Dirichlet data of the pressure field. This is the naive model commonly employed by others but inconsistent with the piezoelectric transduction. The reconstruction is carried out using Algorithm 1 modified by setting κ = 0 which is equivalent to assuming that the measurements are Dirichlet data. After 50 iterations we obtain a relative error of 6.43%. The error profiles for both reconstructions are displayed in Figure 5. Fig. 6: Relative error versus iteration number for white, pink and red noise. In all cases the noise level is at 10%.

The effect of noise added to the piezoelectric measurements is also explored. 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 like ω -1 as ω → ∞. Red or Brownian noise has a PSD decaying like ω -2 . Figure 6 displays the relative error as a function of the iteration number for the three types of noise at a 10% level. We observe that the reconstruction method can handle best the white noise. The error for the pink noise is greater. And the reconstruction error for the red noise is the greatest of the three types of noise. This phenomenon could be explained by studying the spectral characteristics of the discrete normal operators (F * F) or (FF * ). The iterations from the Landweber algorithm correspond to truncated Neumann series [48],

where binomial formula yields

Therefore, the Landweber iterate u K could be expressed as

for some coefficients c kj . This last expression motivates the study of how the discrete normal operator (FF * ) responds to the noise contained in the piezoelectric measurements. Figure 7 displays the average power spectral response to time-and space-tracings of white noise as the input to (FF * ). We observe that the operator (FF * ) suppresses the high frequency components. This explains why the reconstruction algorithm can handle white noise better than pink or red noise. The latter noise types have a larger portion of their power residing over low frequencies. Therefore, the reconstruction algorithm does not suppress those types of noise.

VII. DISCUSSION AND CONCLUSIONS

We have proposed a mathematical model for the acoustic measurements transduced by piezoelectric sensors. This model, taking the form (28), 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 III) and the solvability of the associated photoacoustic problem (see sections IV-V).

The directional response for plane waves derived in section III displayed in Figure 2 matches the experimental measurements carried out by other researchers [11]- [13], [30]. Design considerations for the mechanical and electrical properties of the piezoelectric film (condensed in the parameter κ defined in ( 14)) play a role in the appearance of critical angles where the sensitivity of the sensor vanishes. The incorporation of this type of sensor response into reconstruction algorithms has been highlighted as one of the challenges associated with photoacoustic imaging [7], [24], [50]- [52] and partially investigated in [6], [8]- [10], [21]. Our present paper represents a novel contribution to this research effort.

We note that our mathematical model for the piezoelectric sensor is qualitatively similar to the model for the Fabry-Perot transducer proposed in [21] in spite of the completely different physical principles from which they are derived. Hence, a common mathematical framework can be used to analyze both types of sensors. 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.

We have also proposed an iterative reconstruction Algorithm 1 based on an accelerated Landweber method. We implemented proof-of-concept synthetic simulations to highlight the importance of incorporating the proper modeling of the piezoelectric transducer. We carried out reconstructions with and without the proper piezoelectric model. See the error profiles shown in Figure 5. We observe that certain features, whose wavefronts may have reached the measurement boundary at a non-normal incidence, are missed by the naive reconstruction. We also investigated the effect of noise of different spectral characteristics. The reconstruction method can handle white noise better than pink or red noise. See Figures 6. This is due to suppression of high-frequency components in the measurements. See Figure 7.

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 fact that the measurements were synthetically generated on a mesh different from the reconstruction mesh. 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 ill-conditioned) and the reconstructed images suffer from aberration. Remedies for this phenomenon have been investigated, including regularization, two-grid methods and numerical schemes with lower dispersion. See [45, Sect. 6.8-6.10], [53] and references therein. However, these improvements fall outside of the scope of this paper.

As future research, it would be very important to explicitly model the influence of matching layers commonly employed in the design of piezoelectric sensors [30], [32]. Matching layers play an important role in optimizing the transmission of high-frequency acoustic energy into the sensor. In this paper, this transmission is modeled by (27b) which for a flat surface simplifies to

Here Z b = ρ b c b is the acoustic impedance of the thick backing substrate which usually does not match the acoustic impedance Z = ρc of the fluid or the acoustic impedance Z p = ρ p c p of the piezoelectric film. Matching layers are designed to have an impedance Z ml such that Z < Z ml < Z p Z b so that the wave field experiences a more gradual change of media across the fluid-sensor interface. The validity of the proposed asymptotic 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 I, the wave speed in PVDF materials ranges from 1300 -2300 m/s. Hence, we expect our model to be valid for frequencies 12 MHz. For higher frequencies, the wave field is affected by the presence of the piezoelectric film and the thin matching layer [30]. In such a scenario, our transmission model would be inappropriate. Hence, there remains a need for a more complete model that can accurately handle multiple frequency scales.