*2.3. Thermoacoustic Equation*

Once an initial pressure is found, via either optical or radio frequency electromagnetic radiation, the pressure must then propagate to ultrasound transducers located outside the tissue boundary. Extremely short (nanosecond) pulses of energy are assumed to be irradiating the medium. Therefore the initial pressures are assumed to be delta functions, and only present for the initial conditions when formulating the equations. The pressure wave in photoacoustics is usually modelled using a standard scalar wave equation, with initial conditions:

$$\begin{aligned} \frac{\partial^2}{\partial t^2} p(\mathbf{x}, t) &= c^2(\mathbf{x}) \nabla^2 p(\mathbf{x}, t) \cdot \mathbf{x} \in \Omega\\ p(\mathbf{x}, 0) &= f(\mathbf{x})\\ \frac{\partial}{\partial t} p(\mathbf{x}, 0) &= 0, \end{aligned} \tag{7}$$

where, *f*(*x*) corresponds to the initial pressure found from the initial pulsed excitation, and *p*(*<sup>x</sup>*, *t*) is the pressure wave that propagates to the transducers. Similar to the electromagnetic wave equation, the pressure propagation also requires a perfectly matched layer (PML) to prevent reflections from the boundary of the numerical domain [36,37]. Therefore, the pressure is received at transducers along a curve (surface in 3D) *γ*, which is not the boundary of the domain *∂*Ω.

#### *2.4. Time Reversal Reconstruction Algorithm*

The pressure received at the transducers can be represented as *g*(*y*, *t*), with *y*∈*γ*. The goal of the reconstruction problem is to reconstruct the tissue parameters, *μa* in the optical excitation or *σ* in the RF excitation case, given *g*(*y*, *t*). One way of recovering the parameters is to try and find the initial pressure *f*(*x*) from the pressure measurements *g*(*y*, *t*), and then divide the pressure by the assumed known energy distribution to recover the parameters. One of the more common methods to recover the initial pressure is using the time-reversal technique. For the time-reversal to rigorously reconstruct *f*(*x*), requires Huygen's principle to be valid. Unfortunately, Huygen's principle does not hold when the speed of sound is not constant, and more problematically, it does not hold in two dimensions [38]. Therefore, in most domains of interest the time-reversal method recovers an approximation of the initial pressure *f*(*x*).

As the name suggests, time reversal entails taking the received pressure *g*(*y*, *t*), and using it as a source on *γ*, with time moving in reverse. If the forward simulation was run until a stoppage time *T*, then when simulated using the time reversal method *p*(*<sup>x</sup>*, *T*) ≈ *f*(*x*). There are several other methods that can be used to recover *f*(*x*) besides time-reversal, but there is no universally accepted reconstruction algorithm that works in all cases of interest. Most of these methods require pressure information at the sensors as functions of time. For this reason, the computationally expensive time-domain wave equation was used for solving the pressure wave equation instead of a Helmholtz equation similar to Equation (5). Since the initial pressure acts as a delta function, the received signal is inherently broadband. A time domain simulation is more efficient for measuring broadband response, and provides a better representation of what would be measured during a physical experiment.

## *2.5. Phantom Geometry*

Two phantoms are used to demonstrate the effectiveness of the single simulation tool used in this study. The first homogeneous phantom consists of a circular region of interest with radius 20 mm. Within this region, two objects are placed, a circle of radius 1 mm placed approximately 4 mm to the

left of center, and an ellipse with major axis 2 mm and minor axis 1 mm placed 4 mm to the right of center. The second phantom consists of an approximation to a human breast. The breast region of interest is represented by a circle of radius 50 mm. Glandular tissue rendered as a circle of radius 10 mm surrounds an elliptical tumor with major axis 7 mm and minor axis 3.5 mm. The tumor is located 32.5 mm deep from the right side of the phantom. The optical sources and acoustic detectors are placed around the region of interest in a continuous fashion.

For creating the finite element meshes, a characteristic length of 2 mm was used for efficiency reasons, though near the absorbing objects, the size of the elements smoothly decreased in order to properly model the objects with sufficient resolution. For the homogeneous phantom mesh, the total number of nodes in the mesh was 3635 with 7268 elements. The mesh for the breast phantom had 12,166 nodes and 24,330 elements. For the time domain simulation of the acoustic propagation, a time step of 30 nanoseconds was used in a Newmark numerical integration method. The homogeneous phantom ran for 700 time steps, while the breast phantom ran for 1700 time steps. A constant speed of sound of 1.5 mm/μs was used for both the forward and time-reversed simulations.

The time estimates for the breast phantom are: 0.6 s for the optical simulation, and 1.1 s for the RF simulation. The acoustic simulation took 18.1 s in the optical source case, and 17.3 s in the RF source case. The reconstruction using time reversal took 385.8 s and 408.0 s for the optical and RF source cases, respectively. The homogeneous phantom timings are: 0.16 s and 0.315 s for the optical and RF simulations; 2.1 s and 2.0 s for the forward acoustic simulation; and the time reversal took 38.4 s and 39.4 s for the optical and RF cases, respectively. The timings were done on a Dell Precision 5820 desktop PC, on a single thread. ONELAB has the functionality to run on multiple threads, as well as a graphical processing unit (GPU), but these options were not used in measuring these time estimates.

## *2.6. Phantom Parameters*

The absorption coefficient *μa* of the homogeneous phantom for the background was set to 0.001 mm<sup>−</sup>1. The circle and ellipse, which represent absorbers, have an absorption coefficient of 0.425 mm<sup>−</sup>1; the average absorption coefficient of blood at 800 nm. As is common in human tissue, the background is assumed to have a higher reduced scattering coefficient *μ s* than the absorbers. The reduced scattering coefficient for the background is set to 1 mm<sup>−</sup>1, while the absorbers *μ s* is set to zero for this phantom. In terms of electrical properties, the human body in general does not have a significant magnetic response at radio frequencies, and so the relative permeability *μr* is set to 1 for all objects. A background relative permittivity of *εr* = 5 is similar to human tissue at 434 MHz. Unlike in the PACT case, objects such as tumors have a higher scattering coefficient as well as absorption in the radio frequency regime. Therefore, the relative permittivity of the circle and ellipse were set to *εr* = 25. The conductivity, which governs the amount of absorption, was set to 0.1 S/m for the background, approximating general tissue. The circle and ellipse are given values similar to that of a tumor, 10 S/m [39].

General breast tissue has an absorption coefficient *μa* = 0.005 mm<sup>−</sup>1, with a reduced scattering coefficient of 1.52 mm<sup>−</sup><sup>1</sup> when an 800 nm source is used. The optical properties of the glandular tissue and tumor depend on the assumed material composition of the tissues. The assumed amount of hemoglobin and percentage of water can have a dramatic effect on the properties at any given wavelength. In this study, the same composition as Reference [18] is used for the glandular and tumor tissue. The spectral characteristics of hemoglobin, deoxyhemoglobin, and water at 800 nm are taken from References [40,41]. The specific optical properties used for this breast phantom are given in Table 1.

**Table 1.** Optical properties (*μa* the absorption coefficient, and *μ s* the reduced scattering coefficient) of different breast tissue at 800 nm.


Human tissue is dispersive at radio frequencies, and so the electrical properties vary over frequency, though not as dramatically as the optical parameters. The properties for the generic and glandular tissue at 434 MHz were taken from the ITIS database provided by ETH-Zurich [42], which references a technical report compiled by the United States Air Force [43,44]. The tumor properties were extrapolated from data provided at 100 MHz in Reference [39] using a Cole–Cole dispersion model. The specific electromagnetic parameters are given in Table 2. The relative magnetic permeability is again set to unity, since the body does not exhibit strong magnetic response at 434 MHz.

**Table 2.** Electrical properties (*<sup>ε</sup>r* the relative permittivity, and *σ* the conductivity) of different breast tissue at 434 MHz.

