Pyroelectric ultrasound sensor model: directional response

Ultrasound is typically measured using phase-sensitive piezoelectric sensors. Interest in phase-insensitive sensors has grown recently, with proposed applications including ultrasound attenuation tomography of the breast and acoustic power measurement. One advantage of phase-insensitive detectors, in contrast to conventional phase-sensitive detectors, is that they do not suffer from a narrow directional response at high frequencies due to phase cancellation. A numerical model of a phase-insensitive pyroelectric ultrasound sensor is presented. The model consists of three coupled components run in sequence: acoustic, thermal, and electrical. The acoustic simulation models the propagation and absorption of the incident ultrasound wave. The absorbed acoustic power density is used as a heat source in the thermal simulation of the time-evolution of the temperature in the sensor. Both the acoustic and thermal simulations are performed using the k-Wave MATLAB toolbox with an assumption that shear waves are not supported in the medium. The final component of the model is a pyroelectric circuit model which outputs the sensor response based on the temperature change in the sensor. The modelled pyroelectric sensor response and directional dependence are compared to empirical data.


Introduction
The pyroelectic effect is a property of certain dielectric materials which have a spontaneous electric polarization, where a change in temperature of the material results in a charge 4  separation, and hence a measurable current or electric potential difference [1][2][3]. It has been used in many applications, including infrared sensors, pyroelectric thermometers, laser power meters, [1,4] and, more recently, for ultrasound power measurement [5,6] and ultrasound tomography of the breast [7]. In order to understand the behaviour of pyroelectric ultrasound sensors and to optimise their design for particular applications, an accurate model is needed.
A typical design for a pyroelectric ultrasound sensor is similar to a PVDF (polyvinylidene difluoride) membrane hydrophone in that it consists of a layer of pyroelectric sensor material sandwiched between electrodes which measure the voltage response. The size of the electrodes determine the size of the sensitive region of the sensor, since only the charge separation between the electrodes is measured. Since the sensor responds to heating, a backing layer is added, the purpose of which is to absorb incoming ultrasound as effectively as possible so that the heat is deposited close to the sensitive membrane. A schematic of a generic pyroelectric sensor is shown in figure 1.
Most conventional ultrasound sensors are made from piezoelectric materials which have an electrical response to mechanical deformation of the material. To a first approximation, piezoelectric sensors can be thought of as measuring the acoustic pressure integrated across the surface of the sensor, and can have a sufficiently fast response time to measure the phase of an acoustic wave even at high ultrasonic frequencies. For normally incident waves, a larger piezoelectric sensor will have a higher sensitivity, but for non-normally incidence waves, the variation in phase across the sensor surface leads to cancellation and a lower overall sensitivity. Reducing the size of the piezoelectric sensor reduces this phasecancellation, and when the sensitive region is much smaller than the wavelength an omni-directional response results. However, this also reduces the sensor's sensitivity, and thus the size of a piezoelectric sensor is a trade-off between directionality and sensitivity. Pyroelectric sensors, on the other hand, are inherently phase-insensitive-they respond to the energy deposited in the sensor due to absorption of the ultrasound wave-and so are not subject to these directionality effects due to phase-cancellation, and the sensitivity is only weakly affected by angle of incidence. Furthermore, pyroelectric sensors typically have a much slower temporal response than piezoelectric sensors, since pyroelectric sensors respond to temperature change, whereas piezoelectric sensors respond to deformation. In order to improve the temporal response of a pyroelectric sensor, either the power of the heat source needs to be increased or the sensor must be made to respond to smaller changes in temperature. The former option is not viable for medical imaging, as the heating quickly becomes damaging to tissue, and the latter option makes the sensor sensitive to thermal noise, reducing the signal-to-noise ratio. The directional response of pyroelectric ultrasound sensors has been measured [7], but the reason for the shape of the directional response curve is not well understood, which is one of the principal problems we seek to elucidate in this paper.
An accurate simulation of the pyroelectric effect requires the electric, thermal, and mechanical fields within the pyroelectric material to be modelled. Previous approaches for modelling pyroelectric ultrasound sensors have either relied on simple analytical approximations for the heating and for the acoustic field [5], or have been specialized to particular geometries which may not generalize easily [6]. In the case of pyroelectric infrared sensors, some modelling approaches have used finite element methods [8,9] which can be computationally costly, especially for the acoustic simulations.
The model presented in this paper simulates, in sequence, the acoustic field within the pyroelectric sensor due to an external ultrasound source, the temperature change in the sensor caused by the absorption of acoustic energy from the ultrasound wave, and the pyroelectric voltage response of the sensor. Section 2 gives a brief overview of the theoretical aspects of the pyroelectric effect, and the assumptions that have been made for modelling the pyroelectric sensor The pyroelectric element is parametrized by a set of coordinates Ω, with cross-sectional area A, depth h, and a volume of Vp = |Ω| = A · h. The spontaneous polarization density P S and the pyroelectric coefficient γ are shown to be orthogonal to the sensor surface, and we assume that the electric displacement D and its time derivative, the displacement current density J D , are also orthogonal to the sensor surface. Differential surface element, dA = dx dy, on the pyroelectric sensor material. The volume below the surface element forms a column of height h. The change in surface charge on dA induced by a change in spontaneous polarization due to a change in temperature is assumed to depend only on the average rate of change in temperature in the column.
response. The governing equations for the sensor response are given and an expression for the measured voltage is derived. In section 3, the three components of the model are outlined: acoustic, thermal, and voltage response simulations. In section 4 it is shown how the modelled voltage responses and directional responses vary with changes in certain parameters of interest. We focus especially on the directional response of the sensor, as it is a key feature of pyroelectric sensors that has not previously been modelled. In section 5 the results of the model are compared to measured data. Finally, in section 6 we discuss the results and suggest potential applications of this pyroelectric ultrasound model, as well as aspects which can be further worked on.

Pyroelectric effect
Crystal and semi-crystalline polymer structures which exhibit the pyroelectric effect are a subset of piezoelectric materials, which in turn are a subset of all dielectric materials [1,3]. Pyroelectric materials can broadly be categorized as those materials which exhibit a change in surface charge in response to a change in temperature. More precisely, this means that the electric displacement field D = ε 0 E + P inside the material has a temperature dependence, where E is the electric field, P is the polarization density, and ε 0 is the vacuum permittivity. Pyroelectric materials possess a permanent polarization P S , called the spontaneous polarization, which persists even without an external electric field. The total polarization density P is the sum of the spontaneous polarization density P S and the induced polarization density P ind which arises from an electric field. The spontaneous polarization density is a function of temperature, whereas the induced polarization density is a function of the electric field. PVDF is a non-linear dielectric, but we assume that the external electric field is negligible [10], so that the electric field inside the sensor arises from the spontaneous polarization only. This electric field is relatively weak, and hence we can make a first order approximation of the D-E dependence so that the sensor material is modelled as a linear dielectric, where the induced polarization is proportional to the electric field, P ind = ε 0 (ε r − 1)E, and the electric displacement is then given by D = εE + P S , where ε = ε 0 ε r , and ε r is the relative permittivity. The spontaneous polarization reacts to changes in temperature due to the shifting of the electric dipole moments in the material, resulting in a change in the surface charge of the crystal, which can be measured.

Pyroelectric coefficient
The pyroelectric effect can be described as a coupling of electric, thermal, and mechanical effects in materials, as shown in figure 2. Here we may assume that magnetic effects are negligible, since in the absence of external magnetic fields the only magnetic effects would arise from an inhomogeneous electric field (in time or space), and the only electric fields in our model are weak ones arising from the spontaneous polarization density of the pyroelectric material. The measured pyroelectric effect consists primarily of two separate phenomena: the primary pyroelectric effect, which is the change in electric displacement or polarization due to a change in temperature, and the secondary pyroelectric effect, which is the indirect change in the polarization due to a thermal expansiondriven strain. These effects are highlighted in figure 2.
Using equilibrium thermodynamics, relationships between thermal, mechanical, and electric properties which give rise to the pyroelectric effect can be expressed up to first order by the coupled differential forms for D = D(σ, E, T) and S = S(σ, E, T), given by [1,2,12] where σ is the stress tensor, T is the temperature, and the superscripts indicate variables which are held constant.
(A similar expression could also be written for entropy, χ, but it does not affect the electric displacement and hence does not contribute to the pyroelectric effect; instead it is involved in the converse effect, the electrocaloric effect.) The pyroelectric coefficient describes how the spontaneous polarization of a material changes with temperature and is given by [1,12,13] Although the pyroelectric coefficient is a vector quantity, most pyroelectric crystals are only polarizable in a single direction, and in practice the pyroelectric crystal is poled orthogonal to a flat surface, so that with respect to the surface the pyroelectric coefficient can be treated as a scalar quantity [14]. The two primary components contributing to the pyroelectric coefficient in equation (3) can be expressed separately using equations (1) and (2), detailed below.
Equations (1) and (2) may be written in more compact forms by replacing the partial derivatives with the tensor coefficients they quantify, where we have used the Einstein summation convention, and we have used the fact that the coefficients which quantify the direct piezoelectric effect are equal to the coefficients quantifying the converse piezoelectric effect [2]. d ijk are the coefficients of the rank-3 piezoelectric tensor, ε ij are the coefficients of the rank-2 permittivity tensor, γ i are the coefficients of the rank-1 pyroelectric tensor, s ijkℓ are the coefficients of the rank-4 elastic compliance tensor, and α ij are the coefficients of the rank-2 thermal expansion tensor. Assuming a constant electric field, dE = 0, and constant strain, dS = 0, equations (4) and (5) can be rearranged to isolate the stress differential terms, dσ ij , and then equated to get the expression for the total pyroelectric effect, by making use of the tensor inverse identities (d −1 ) imn d ijk = δ mj δ nk and (s −1 ) ijmn s ijkℓ = δ mk δ nℓ , where δ ij is the Kronecker delta function. For many pyroelectric materials, the secondary pyroelectric effect has a non-negligible contribution to the total effect, and in the case of PVDF, a polymer widely used for ultrasound detection, the primary and secondary effects have nearly equal Diagram depicting the relationships between electrical, mechanical, and thermal effects in crystals. Intensive variables comprise the outer triangle, and extensive variables the inner triangle. The ranks of the tensors corresponding to the variables in the circles are shown in round brackets, while the ranks of tensors corresponding to properties between variables are shown in square brackets. The primary pyroelectric effect is highlighted in blue, while the secondary pyroelectric effect is highlighted in red. Diagram based on version seen in J F Nye (1985) and S B Lang (1974) [1,2]. Reproduced with permission from [2]. The diagram was originally conceived by G Heckmann (1925) [11]. Reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature, Springer Berlin Heidelberg.
contribution [10]. We will therefore consider the total pyroelectric coefficient in equation (6) as the pyroelectric coefficient of interest. There also exists a tertiary pyroelectric effect which may arise in piezoelectric materials from non-uniform heating, leading to non-uniform stresses from thermal pressure in the material [1,15]. We assume that the tertiary effect is negligible, so that pyroelectric response dependence on thermal gradients may be ignored.

Pyroelectric current response
The electric displacement field in the pyroelectric material induces a surface charge density on the faces of the material normal to the electric displacement field. The time rate of change of the electric displacement then leads to a change in this surface charge, which can be measured as a current by attaching electrodes to the surfaces of the pyroelectric material.
Several simplifying assumptions are made for describing the pyroelectric response of the sensor in the model presented here. It is assumed that for the purposes of modelling the pyroelectric current response to heating from an ultrasound source, the surface charge on any particular area element dA (refer to figure 1) of the pyroelectric sensor depends only on the change in temperature below that area, and not to changes in the temperature beneath adjacent area elements. This assumption is justified by a combination of 1-dimensional models [5] and by the thermal diffusion in the sensor propagating primarily in the direction normal to the surface of the sensor, which in turn is justified by the heat equation where heat flow follows the gradient of temperature. The temperature gradient along the surface of the sensor is relatively uniform compared to the gradient normal to the surface, since the time scale of thermal diffusion is on the order of t diff ≈ L 2 /4κ, where L is the length scale of thermal diffusion and κ is the thermal diffusivity, whereas the time scale it takes for the full acoustic  wave to reach the surface of the sensor has the upper bound t delay ≈ d sep /c, where d sep is the difference between the longest and shortest distances travelled by the wave-front of the acoustic wave to reach the sensor surface (for a plane wave at normal incidence, this distance is zero), and c is the sound speed in the background medium. For a sensor with a width of d = 50 mm and a water background, we can estimate that d sep = d, so that t delay = 50 mm/c water = 33.4 µs.
For typical pyroelectric materials such as polyvinylidene difluoride (PVDF) and lead zirconate titanate (PZT) the thermal diffusion is on the order of κ ≈ 10 −7 m 2 s −1 , so the length scale of thermal diffusion over a time span of t delay is L ≈ 3.7 µm, which is negligible compared to the width of the sensor. This assumption does neglect edge effects in cases where the width of the incoming ultrasound wave is smaller than the width of the sensor, in which case we expect heat to spread towards the edges of the sensor as well as along the depth, but this should be a small effect, since even for long diffusion times on the order of a second, the thermal diffusion length scale is about 0.63 mm, which is still considerably smaller than the width of a typical sensor.
It is also assumed that the spontaneous polarization P S within the pyroeletric material is uniform at a given temperature T, and is normal to the surface of the sensor for all temperatures, hence implying that the pyroelectric coefficient γ σ,E is also uniform and normal to the sensor surface. The uniformity of the spontaneous polarization along the depth of the pyroelectric material and the temperature stability of the pyroelectric coefficient are supported fairly well by experimental results for PVDF, for temperatures below 40 • C [16][17][18][19].
indicating that we may expect the pyroelectric current to have a proportional dependence on the pyroelectric coefficient and the time rate of change of temperature. According to previous works, in the 1-dimensional case where the temperature only varies along the depth of the pyroelectric crystal and is uniform across the slices parallel to the crystal surface, the pyroelectric current response is given by the average rate of change in temperature along the depth of the crystal (along the direction of polarization) [20][21][22][23][24]. We will follow the arguments presented in [21], and extend the result for temperatures which vary throughout the crystal by applying the method along columns below surface elements dA, assuming that the change in charge induced on the surface element dA is proportional only to the average rate of change in temperature directly below it in the sensor. By doing this, we can get the pyroelectric current density at the surface of the sensor.
We partition the sensor material into N slabs of thickness ∆z i for i ∈ {1, . . . , N}, where we consider each slab being sufficiently thin so that it can be considered to have a constant spontaneous polarization density, P S,i = P S,iẑ . If we consider the volumes beneath area elements dA of the sensor element as parallel plate capacitors with capacitance where ε = ε r ε 0 is the permittivity of the pyroelectric material, ε 0 is the permittivity of the vacuum, and ε r is the relative permittivity, then we can relate the surface charge element to the potential difference between the top and bottom surfaces of the pyroelectric sensor by The potential difference V h 0 is the line integral of the electric field between the top and bottom electrodes of the sensor, where the electric field in this case is induced by the internal surface charges caused by the jump discontinuities in the polarization between adjacent slabs. We can relate the electric field to the polarization density field through the field boundary conditions between adjacent layers, where ς f,i and ς b,i are the free and bound surface charge densities between layers i and i + 1, respectively. We assume that the sensor is a perfect insulator and hence has no free charges, so ς f,i = 0 for i = 1, . . . , N. Within the electrodes at i = 0 and i = N + 1, the polarization density and electric fields are zero since the electrodes are conductors, and we assume that the momentary electric field within the electrodes caused by changes in polarization density in the sensor are short-lived compared to the thermal timescales causing the changes, and can hence be neglected. Combining these, we get that We will assume that the dielectric constant ε i is constant throughout the sensor material, so we can let ε i = ε. Since the voltage V h 0 is the line integral of the electric field across the thickness of the sensor, it follows from equation (10) that the surface charge element can be expressed as By expressing the differential charge as the change in charge from an initial value, dq h = dq h (0) + dq h (t) and linearly expanding the spontaneous polarization density with respect to temperature as P S,i ≈ P S,i | T0 + ∂PS,i ∂T T0 T i (t), we can express the differential change in charge in the limit as N → ∞ by where γ is again the pyroelectric coefficient. We are interested in the magnitude of the change in charge, so we will drop the sign from now on. We assume that the pyroelectric coefficient is uniform within the sensor, so by taking the time derivative of equation (15) and differentiating both sides with respect to the area, we get the pyroelectric current density, Integrating equation (16) over the surface yields the total pyroelectric current, which is proportional to the pyroelectric coefficient and the surface integral of the depth-averaged rate of change in temperature, as shown in the first equality in (17). Alternatively, we can split the proportionality to the volume-averaged rate of change in temperature and the surface area of the sensor, shown in the second equality of (17).

Ultrasound sensor model
The model described here consists of three coupled components, which together simulate the pyroelectric response of an ultrasound sensor to an ultrasound wave.

Acoustic simulation
The acoustic simulation in this model is done using the k-Wave MATLAB package, assuming linear acoustic propagation with absorption in a fluid medium, meaning that shear waves are assumed to be negligible. The dynamics of these acoustic waves are governed by the equations: Sensor Capacitance, C (nF) 0.1

Time [s]
Normalised voltage response where d and u are the particle displacement and velocity, ρ and p are the acoustic density and pressure, c 0 and ρ 0 are the equilibrium sound speed and density, S M is a mass source term, and L is a linear loss operator. The operatorL describes absorption and dispersion that follows a frequency power law of the form α = α 0 ω y , where y is the power law exponent. It is defined using the fractional Laplacian and has the form, where α 1 (x) = −2α 0 (x)c 0 (x) y−1 is the absorption proportionality coefficient, α 2 (x) = 2α 0 (x)c 0 (x) y tan(πy/2) is the dispersion proportionality coefficient [25]. The volume rate of heat deposition, Q, is determined by a combination of the absorption characteristics of the medium as well as the form of the acoustic wave applied to the medium. In general, we can describe the heat deposition term as a moving mean of dissipation, averaged over a characteristic period τ of the acoustic wave, where D is a term quantifying the dissipation of acoustic energy, and ⟨·⟩ τ denotes the temporal moving mean over a window of length τ , which can be expressed as The characteristic period τ represents a time-scale over which the small-scale periodic behaviour of the acoustic wave can be averaged, so that only the macroscopic behaviour of the wave is taken into account. The form of the dissipation term D depends on how the absorption is modelled, but it can be related to acoustic energy density and acoustic intensity by treating it as an energy loss term in the acoustic wave, so that it obeys the energy conservation corollary [26] ∂w(x, t) ∂t where w is the acoustic energy density and I is the acoustic intensity. If we were to be more precise, in the case of acoustic absorption due to thermal conduction and viscous effects, w and I depend on the temperature of the medium, and I further depends on the stress tensor σ. However, we will assume uniform stress dσ = 0 throughout the material and a sufficiently well-behaved applied acoustic wave such that I = pu is a reasonable approximation for the acoustic intensity. In the steadystate we can say that the first term in equation (24) is zero, and hence the acoustic intensity and dissipation are left with trivial time dependencies. Hence, we can take the moving time average of the acoustic intensity over an arbitrary window of length τ that is an integer multiple of the period of the ultrasound wave. Thus, the heating term due to dissipation is given by which is then used in the thermal simulation.

Thermal simulation
The thermal simulation can be greatly simplified if we can neglect diffusion, which is valid if the timescale of thermal diffusion is much greater than the heating duration. The heat diffusion timescale can be characterized by the thermal relaxation time, t thermal , which is the time that it takes the temperature of a heated region to drop to 1/e of its initial temperature. An approximation for t thermal for acoustic plane wave heating is given by the equation where κ is the thermal diffusivity of the material, and d is a characteristic length. For the PVDF layer d is the thickness h of the membrane, whereas for the backing layer the characteristic length is the acoustic penetration depth of energy, 1/2α, where α is the acoustic attenuation. Assuming a PVDF membrane with thickness 100 µm along with a backing layer with attenuation coefficient α = 760 NPm −1 at 1 MHz. For the PVDF pyroelectric sensor with material properties given in table 2, the thermal relaxation time in the PVDF layer is 13 ms, whereas for the backing layer, assuming a 1 MHz ultrasound source, the thermal relaxation time is 564 ms. At higher ultrasound frequencies the thermal relaxation time in the backing layer would be even shorter. These times are comparable to the heating durations considered in this paper, and hence we need to use a thermal model that accounts for diffusion. The thermal dynamics of the sensor are also modelled in k-Wave, which uses the heat equation: where c p is the specific heat capacity, K is the thermal conductivity, and Q is the volume rate of heat deposition as defined by equation (25). This part of the simulation gives us the temperature field, T(x, t), within the sensor region Ω. The spatial mean of the temperature field within this region, ⟨∆T(x, t)⟩ Ω , gives us the pyroelectric current response through equation (17), which is used as an input in the final part of the model-the voltage response.

Voltage response
An electrical circuit model approach for the pyroelectric sensor considers the sensor as a current generator connected in parallel with a capacitor and a resistor, which correspond to the capacitance and resistance of the material, respectively [1,5,27]. In the version of the circuit model seen in figure 3.3, presented in [5], the pyroelectric circuit component is connected to a unit gain amplifier with an input resistance R a , which in turn is connected to a shunt resistor of resistance R L . The current i(t) measured across the electrodes attached to the sensor surfaces is given by equation (17), and so the voltage response can be attained using Kirchhoff's laws to yield [1,5,27] where R eq is the parallel combination of R p and R a , and the values for capacitance C and resistance R p of the sensor correspond to their total values across the electrodes. The thermal simulation gives us the change in temperature field ∆T(x, t) in the sensor, from which we can numerically solve for the pyroelectric sensor response, V(t), using equation (28).
The full pyroelectric signal is the voltage as a function of time, V(t). In this paper we take the absolute maximum of V(t) to reduce the signal to a single quantity that can be compared between measurements. More technical approaches can also be considered, such as the correlation scheme outlined in [28], where the cooling time of the sensor in physical measurements is accounted for.

Sensor voltage and directional responses
Using the model from the previous section, we investigate how various sensor parameters affect the pyroelectric sensor response. The model contains many parameters that can be independently varied over large ranges, but we focus on a few key physical parameters of particular interest: sensitive region diameter (d s , the diameter of the electrode, as shown in figure 6), PVDF thickness, source-sensor separation, thermal diffusivity, and the electrical parameters R eq and C. For these simulations, we will use the parameters outlined in table 1 as the control, which are chosen to be close to typical values found in pyroelectric ultrasound sensors in the literature. We will additionally use the acoustic and thermal properties for the PVDF and backing layer as listed in table 2, and the water background is assumed to be non-attenuating.
The results presented in this paper were computed using a 2D model for faster computation times. The computational grid used in these simulations had a grid stepsize of 100 µm per pixel. The voltage response and directional response were tested for convergence by checking that increasing the grid resolution to a stepsize of 50 µm produced voltage response curves with relative root mean squared difference of less than 3%, and directional response curves with a difference of less than 1%.

Sensor voltage response
To model the sensor voltage response, we consider situations where the sensor is undriven (there is no ultrasound) for the initial t in = 2 s, is then exposed to constant heating due to the absorption of a normally-incident wave for t Q = 2 s, and then is allowed to cool down for a further t cool = 6 s. Figure 4 shows how the diameter of the sensitive region affects the voltage curve at various angles of incidence for the ultrasound wave. Increasing the sensor size results in a greater voltage response, since they capture a greater portion of the charge separation arising from the pyroelectric effect. Increasing the sensor size also causes the voltage curves at different angles of incidence θ to converge to one-another, indicating a flatter directional response. The directional response aspect is investigated more in section 4.2. Figure 5 shows how varying the electrical parameters R eq and C affects the voltage response curves. There appears to be a symmetry with respect to changes in R eq or C, where lower values of R eq C result in fast voltage responses with sharp peaks, whereas larger values of R eq C result in slow voltage responses with smooth peaks. From equation (28) we can see that for large values of R eq C, the voltage response becomes proportional to the temperature change and for small R eq C it becomes proportional to the rate of temperature change. This is discussed further in section 6.2.

Directional response
To model the sensor directional response, we consider situations where the sensor is exposed to constant heating due to ultrasound absorption for t Q = 0.15 seconds, and then is allowed to cool down for t cool = 0.05 s. We used (27) to obtain ∂T ∂t and use this as the left hand side in (28); the latter equation was then solved using Matlab functions ode45 and ode15s. As we are taking the maximum voltage response for each angle, the cooling portion does not actually affect the directional response as the peak voltage amplitude is reached within the first 150 ms. A directional response measurement is done by varying the angle of a beam-like source with respect to the surface of the sensor, while maintaining the same separation between the source and sensor, as depicted in figure 6. We therefore define our directional response as (29) Figures 7, 8, and 9 show how changes in the geometry of the source-sensor arrangement affect the directional response. The separation between the source and sensor does not have much effect beyond a very short separation of about 1 cm, whereas the sensitive region diameter and PVDF thickness seem to affect the directionality significantly. The sourcesensor separation should not affect the directional response unless there are significant near-field effects, and a sensor with a sensitive region diameter of 1 mm is likely to be less affected by these unless the separation is very small, as we see in figure 9. The PVDF thickness and sensor width, however, directly affect the heating profile of the PVDF layer, so they should be expected to affect the directional response.
A more direct way to affect the heating of the sensor is by adjusting the thermal conductivity, K, of the PVDF and backing layers. In figure 10 we see how increasing the thermal conductivity affects the directional response. With higher thermal conductivity, the sensor response seems to become more stable near normal incidence, and diminish less overall for all angles from 0 • to 40 • . This could be because with higher thermal conductivity, the loss in the power density of the incoming ultrasound with increased angle is compensated for by the fast transfer of heat in the highly thermally conductive medium.
Finally, we again look at how changes in the total resistance R eq and sensor capacitance C affect our modelled results, this time for directional response. We can see in figure 11 how for the same acoustic and thermal simulation, variations in the electrical parameters affect the directivity. Aside from differences in smoothness, no significant differences are apparent in the qualitative nature of the response curves.

Comparison to experimental measurements
The pyroelectric ultrasound sensor model outlined in section 3 was tested against experimental data described in [5,7], with experimental parameters outlined in tables 2 and 3. Both the voltage response and directional response were tested with appropriate parameters for comparison. For the purpose of saving on computation time, the simulations were done with two spatial dimensions instead of three.

Voltage response
The most fundamental aspect that the model must replicate is the voltage response of the pyroelectric sensor. Using the parameters from tables 2 and 3, chosen to match the parameters used in empirical tests for pyroelectric ultrasound sensors at NPL, we can simulate the same setup where a 1 MHz ultrasound piston source is driven for 2 s in front of a 50 mm sensitive region diameter pyroelectric sensor, then turned off, after which the sensor is allowed to cool down for another 6 s. A comparison of the measured and modelled normalised voltage responses is shown in figure 12. The normalised simulated voltage curve has a relative root mean squared error of 26.4% compared to the normalised experimental voltage curve.

Directional response
Using the parameters from tables 2 and 3, we simulated a directional response measurement for a pyroelectric ultrasound sensor in the farfield of a 2.03 MHz source, driven long enough so that a steady-state is attained inside the sensor volume. In the 2D simulation, the ultrasound source is approximated by a line source, which is rotated with respect to the center of the sensor surface, in such a way that the center of the sensor surface and the center of the source maintain a constant distance of R for all angles, as shown in figure 6.
The pyroelectric ultrasound sensor simulation is used to calculate the directional response in the case of a 1 mm sensitive region diameter sensor. There is a symmetry with respect to the normal incidence, as the directional responses for angles θ and −θ are identical, so the simulation only needs to be run for positive angles, and then mirrored to get the full directional response. The empirically measured directional response and the corresponding simulated directional response are shown in figure 13. The phase-sensitive directional responses for the same sensitive region diameters are shown for comparison. The phase-sensitive sensor was modelled by recording the peak of the average acoustic pressure over the sensor surface, V PS ≡ max t |⟨ p(x, t⟩ ∂Ωt |, where ∂Ω t denotes the top surface of the sensor volume Ω. The shape of the phase-sensitive directionality curve is very responsive to the diameter of the sensitive region of the sensor, so an adjusted diameter of 0.925 mm was used to produce a better fit. This is reasonable, as there is often some sensitivity outside the geometrically defined sensitive region due to inherent limitations on how constrained the electric field can be during poling. The modelled pyroelectric response in figure 13 captures the slope of the directional response. The central peak of the experimental pyroelectric directional response is believed to be an artifact arising from reflections between the sensor and the ultrasound transducer, which were not included in the model [7].

Shape of the directional response
We argued in section 3.2 that diffusion is not negligible for heating durations above 10 ms For the experimental directionality measurements considered in section 5.2, the acoustic source was driven for 0.15 s, so diffusion should be expected to contribute to the pyroelectric signal.
We can compare how much of the pyroelectric signal comes from direct heating of the PVDF layer and how much comes from diffusion from the backing layer by making modified heat source terms, Q PVDF , and Q ABS for some original heat source Q. We define Q PVDF to be equal to Q within the PVDF layer and zero everywhere else, and Q ABS to be equal to Q in the backing layer and zero everywhere else. By using these heat sources in the thermal simulation, we get the directional response comparison in figure 14. The pyroelectric response due to diffusion from the backing layer is on average 3.3 times greater than the pyroelectric response from the direct heating of the PVDF layer, indicating that for heating a duration of 0.15 seconds with a 2 MHz ultrasound beam, the backing layer is the primary contributor to the sensor response. If we make a modified Q ABS2 by setting the rate of heating within Q ABS to zero for depths beyond 0.5 mm away from the top surface of the backing layer, we find that the resulting pyroelectric response from the Q ABS2 heating term has a relative root mean squared error of 0.68 % from the Q ABS pyroelectric response. Thus, the majority of the directional response arises from heating in the top 0.5 mm of the backing layer, with approximately 25 % of the response arising from the direct heating of the PVDF membrane. The distribution of the rate of heat deposition in the crosssections of the sensor maintains a fairly constant shape in the top layers of the sensor. Additionally, since the heating of the  sensor along the lateral direction is almost uniform, the heat flow from the backing layer into the PVDF is nearly linear, resulting in an enhancement of the direct heating caused directional response. Because of this, we can expect the pyroelectric directional response to match closely with the directionality of the total rate of heat deposition in the PVDF layer, Ω Q dV. In figure 15 we compare the directional responses of the 1 mm sensor to the total power in the PVDF layer and see that there is a strong agreement. In cases where thermal diffusion can be neglected, these curves should be expected to match exactly. However we can also see from figure 10 that in situations where thermal diffusion is very high, this equivalence between the directionality in V and´Ω Q dV no longer holds.  Modelled directional response compared to measured directional response from [7]. Included are the directional responses of 1 mm phase insensitive (PI) pyroelectric sensors and phase sensitive (PS) piezoelectric sensors. The modelled piezoelectric sensor was given a sensitive region diameter of 0.925 mm for a better fit with the experimental data, where the nominal sensitive region diameter is 1 mm. The modelled pyroelectric directional response was normalized to have a peak value of 1.075.

Voltage response curve
The voltage response V(t) in our model is the solution of the first order ODE in equation (28), which can equivalently be expressed in the integral form: where we have assumed the initial condition V(0) = 0, and have defined ∆T(s) = ⟨∆T(x, s)⟩ Ω . For short heating durations the change in temperature in the PVDF during heating is almost linear, so we can approximate that d∆T(s) ds = d∆T dt = O (1), which allows us to compute an analytical expression for the voltage response  The appearance of the product R eq C within the exponential term in equation (31) explains why the normalized voltage responses and directional responses along the bottom-left to top-right diagonals in figures 5 and 11 are identical to one another. The smaller the duration of R eq C, the sooner the voltage reaches its peak value of γAR eq d∆T dt . From this observation we can make a rough evaluation of the feasibility of using pyroelectric ultrasound sensors to measure the time of arrival of an ultrasound beam. At R eq C = 0.1 µs, V a (t) achieves over 99% of its asymptotic peak value in less than 0.5 µs, which is the period of a 2 MHz wave, meaning that we would need R eq C less than this in order to be able to measure the phase of the incoming ultrasound wave. If we assume that the pyroelectric sensor can be treated like a parallel plate capacitor with a dielectric medium, then we can also write C = εA/h. If we furthermore treat the sensor as a thick wire, then we may write its resistance R p = ξh/A, where ξ is the resistivity of the pyroelectric material, and R eq C = R p C/(1 + R p /R a ) = ξε/(1 + R p /R a ). Lowering the amplifier resistance R a so that it approaches zero would also make R eq C arbitrarily small, but in practice the resistance of the amplifier would eventually be so low that it would be dominated by the resistance of the connecting wires. Assuming copper wire which has a resistivity of approximately ξ a = 1.68 × 10 −8 Ωm, with length 10 cm and cross-sectional area of 1 mm 2 , we would have that R a = 1.68 mΩ. Pyroelectric materials are insulators, so the total resistance would be nearly equal to the resistance of the amplifier, R a . The capacitance required to achieve R eq C < 0.1 µs is then C < 59.5 µF, which is easily achievable. However, lowering the amplifier resistance to this extent would likely lower the sensor signal response far below the noise threshold. Faster response times can also be achieved by increasing the pyroelectric sensor resistance, R p . Without significantly altering the membrane geometry, this would require the use of a novel material with a significantly higher resistivity.

Validity of model assumptions
One of the major assumptions made in this model is that the spontaneous polarization density P S , and electric displacement field D, are orthogonal to the sensor surface for all temperatures. This assumption allows us to use the slab assumption to argue that the current response is proportional to the volume averaged temperature. In cases where this assumption fails to be accurate, the voltage across the pyroelectric material would need to be considered from more general assumptions, such as having an internal bound charge density η b = −∇ · P S .
For large temperature gradients a tertiary pyroelectric effect may also contribute to the voltage response, which our model does not account for. However, for thin membranes, the only situation where large temperature gradients may be expected are when the sensitive region diameter is larger than the cross-section of the incident ultrasound field, so that there exists a boundary between the heated and unheated portions of the sensor. In this paper we have only considered sensitive region diameters that are at most as large as the source beam diameter.

Conclusion
This paper has described and demonstrated a numerical model for pyroelectric ultrasound sensors. The model was shown to agree well with experimental data for the voltage response and directional response of a pyroelectric ultrasound sensor up to a normalization factor.
The modelling was also used to give insights into the shapes of the responses. By varying the parameters of the model, significant changes in the voltage response and directional response were observed. The observation that larger sensors and thinner PVDF membranes appear to produce a flatter directional responses could serve to steer the direction of sensor design in the context of ultrasound tomography.
Potential future work could investigate the extent to which shear waves affect the model results, and in what circumstances. This model can also be used to simulate pyroelectric UST data, and opens up the possibility of performing model based inversion of such data for imaging.