Inverse problems of combined photoacoustic and optical coherence tomography

Optical coherence tomography (OCT) and photoacoustic tomography are emerging non‐invasive biological and medical imaging techniques. It is a recent trend in experimental science to design experiments that perform photoacoustic tomography and OCT imaging at once. In this paper, we present a mathematical model describing the dual experiment. Because OCT is mathematically modelled by Maxwell's equations or some simplifications of it, whereas the light propagation in quantitative photoacoustics is modelled by (simplifications of) the radiative transfer equation, the first step in the derivation of a mathematical model of the dual experiment is to obtain a unified mathematical description, which in our case are Maxwell's equations. As a by‐product, we therefore derive a new mathematical model of photoacoustic tomography based on Maxwell's equations. It is well known by now that without additional assumptions on the medium, it is not possible to uniquely reconstruct all optical parameters from either one of these modalities alone. We show that in the combined approach, one has additional information, compared with a single modality, and the inverse problem of reconstruction of the optical parameters becomes feasible. © 2016 The Authors. Mathematical Methods in the Applied Sciences Published by John Wiley & Sons Ltd.


Introduction
Only recently, there have been developed experimental setups that can perform photoacoustic and optical coherence tomography experiments in parallel; see [1] and further developments in [2,3]. Combined setups can be used for imaging biological tissue up to a depth of approximately 5 mm. They have been used for in vivo studies of human and mouse skins. In the current state of experiments, the two recorded modalities are visualised by superposition after registration, and we refer the reader to the review paper [4].
In this paper, we derive a mathematical model for quantitative imaging of the multi-modal experiment based on Maxwell's equations. This study characterises the additional information obtainable by the dual experiment. In this paper, we assume that the sample is an inhomogeneous, isotropic, non-magnetic and linear dielectric medium, and the imaging parameters of the medium are the conductivity, the susceptibility and the Grüneisen parameter.
Even though the mathematical modelling of optical coherence tomography based on Maxwell's equations is well established, most of the literature relies on simplified models of the Helmholtz equation [5][6][7][8][9], where, as a consequence, the frequency dependence of the susceptibility is neglected. It is a recent trend to consider the more general case of Maxwell's equations for the mathematical description of the optical coherence tomography (OCT) system [10][11][12][13].
The mathematical modelling of quantitative photoacoustics, on the other hand, is based on the radiative transfer equation or its diffusion approximation [14,15]. Maxwell's equations have only been considered for mathematical models in thermoacoustic tomography where low-frequency radiation is used for illumination [16,17].
To formulate the dual-modal setup, we consider Maxwell's equations as the basic modelling equations because both imaging techniques rely on the same excitation. In particular, this provides a more general model for quantitative photoacoustics based on Maxwell's equations, which is also applicable in the case of high frequency radiation. This paper is organised as follows. In Section 2, we give an overview of the inverse problem in photoacoustics. We present the governing equations and the commonly used simplifications. Next, we formulate mathematically the problem of optical coherence  (5) tomography starting from Maxwell's microscopic equations. For an extensive mathematical modelling of optical coherence tomography, we refer to [11]. The system of equations for the dual experiment is obtained in Section 4. In Section 5, the equation system is transformed to a Fredholm integral equation for the Grüneisen parameter under the far field and Born approximations, respectively. Finally, we discuss the connection between the optical parameters appearing in the radiative transfer equation and Maxwell's equations. To simplify the reading, we summarise the basic notations used in this paper in Table I.

Photoacoustic imaging
Photoacoustic tomography is a hybrid imaging technique, which measures the acoustic response of an object upon illumination with an electromagnetic wave. In the experiments, the object is illuminated with a short laser pulse, typically with light in the near infrared spectrum. This pulse is absorbed by the medium, where the exact amount varies locally depending on the material properties, in particular on the absorption coefficient, and thereby, a fraction is transformed into heat energy. This leads to a local change of pressure in the medium, where the transformation factor is again depending on the local thermodynamic properties of the medium, typically reduced to the single Grüneisen parameter. Finally, this pressure propagates as an ultrasonic wave through the medium and is recorded at every point on a surface around the object as a function of time.
For a mathematical description of this measurement setup, we have to model the three different phenomena: (i) the propagation of the laser light in the medium; (ii) the transformation of the absorbed energy into a pressure distribution; and (iii) the propagation of the pressure wave inside the medium.
For the reader's convenience, we summarise in Table II the variables, which are used throughout this section, as a backup reference. The light propagation is usually modelled via the radiative transfer equation. Assuming that we have a laser pulse of fixed frequency > 0, the photon density # .t, x/ at time t 2 R and coordinate x 2 R 3 moving in the direction # 2 S 2 fulfils the equation where a is the absorption coefficient and s is the scattering coefficient at the radiation frequency ; see for example [18]. Here, the constant c denotes the speed of light, and ‚ is the phase function, that is, ‚.x, Q # , # / is the probability density that a photon heading into a direction Q # 2 S 2 is scattered at the position x 2 R 3 into the direction # 2 S 2 . Because in photoacoustic imaging, the laser excitation occurs to be very short and the absorption of the light happens much faster than the propagation of the acoustic wave, the light distribution as a function of time is not relevant, but only the total energy absorbed at each point matters. Thus, we turn our attention to the total energy fluenceˆ# originating from photons moving in the direction # 2 S 2 as new variable:ˆ# where h denotes the Planck constant. Then,ˆ# solves the equation However, because the phase function ‚ is unknown, often the diffusion approximation div .
with the diffusion coefficient given by is used, which is an equation for the total energy fluence x/, # i for some functions ‚ 1 and 1 , which is a good approximation if we have a strongly scattering medium; see [19,Chapter XXI,§5]. The absorbed energy at a point x is then given by a .x/ N .x/. Moreover, it is common to assume that the produced pressure density p .0/ inside the medium is proportional to the absorbed energy where the proportionality coefficient is the Grüneisen parameter : For the propagation of the acoustic wave, the standard assumption is that we have an elastic medium with mass density %, bulk modulus K and vanishing shear modulus. Then, according to linear elasticity theory, for instance [20], we obtain the equation Moreover, the density is usually considered to be approximatively constant throughout the medium, so that we may simplify this equation to the linear wave equation with the speed of sound c s D q K % . Finally, we acquire the measurements at some detector surface D R 3 . Putting all this together, the inverse problem of photoacoustic imaging is to recover the material parameters a , s (which in this approximation only enters via the diffusion coefficient ; see (4)), and c s from the given measurement data g.
This problem can be solved in two steps. First, we consider the acoustic part, given by the equation system For given, non-trapping speed of sound c s (in practice, a common assumption is that c s is approximatively constant), the measurements g allow to uniquely reconstruct the initial pressure p .0/ ; see [21]. If the detector surface has a simple geometry (the easiest examples are a planar or a spherical detector) and the speed of sound is assumed to be constant, then there are a lot of explicit reconstruction formulas for p .0/ available: [22][23][24][25][26][27][28][29][30][31][32][33].
For unknown speed of sound, it was suggested in [34] to use multiple photoacoustic measurements with a focussed illumination; see also [35,36] for reconstruction formulas for this kind of illumination. For a more extensive review of these works, we refer to [37].
So, let us assume that we can invert the acoustic problem and recover p .0/ from the photoacoustic measurements. Then, as the second step, it remains the problem of reconstructing the material parameters a , s , from the internal data p .0/ via the relations where we assumed that we know the boundary data N .0/ at the detector surface, because we are expected to know all the material parameters outside the object, in particular in the vicinity of the detector, and thus, p .0/ allows us to calculate N around the detector surface. However, in practice, these data are not used, because the absorption is considered very small outside the object of interest. We are only using these data mathematically to be able to define the modelling equations in all R 3 , which is state of the art in this field. In [38], it was shown that even with multiple photoacoustic measurements, it is not possible to reconstruct all three parameters, but from the knowledge of one of the parameters (usually, the Grüneisen parameter is assumed to be approximatively constant, although there is no physical evidence for this), the other two can be explicitly calculated; see also [39] for a review about the reconstruction from internal data.
Thus, the inverse photoacoustic problem is not well posed, as it requires an additional three-dimensional data to allow the reconstruction of all material parameters.

Optical coherence tomography
Opposed to photoacoustic imaging, optical coherence tomography is not a hybrid imaging technique. It visualises the back-scattered light of an object. However, analogously to photoacoustic imaging, the object is illuminated with a short laser pulse in the visible or near infrared spectrum. In practical realisations, the illumination is triggered to produce a plane wave, and the scattered wave is measured at a planar detector parallel to the plane wave illumination.
Because the time resolution of the detectors, which would be required for the desired spatial resolution, is practically hard to achieve, the detection of the reflected wave is accomplished by superimposing it with a reference beam, which is a copy of the wave used for illumination (produced by a beam splitter from the original wave) sent to a mirror instead of to the object. Then, the total intensity of this superposition is measured at every point on the detection plane.
For the reader's convenience, we summarise the variables used in this section, in Table III.

Microscopic Maxwell's equations
The interaction of the laser pulse, which is an electromagnetic wave, and the medium is best described by Maxwell's microscopic equations. These are equations for the time evolution of the microscopic electric and magnetic fields e and b for given microscopic charge density and microscopic electric current j : If the initial data e.0, / and b.0, / are compatible (that is, they are satisfying the Eqns (11a) and (11b) at t D 0) and the continuity equation holds, then the Eqns (11c) and (11d) imply (11a) and (11b). To see this, take the divergence of the two vector Eqns (11c) and (11d). Then, using the identity div x curl x f D 0 for every vector field f , it follows with (12) that Integrating these equations over time gives the Eqns (11a) and (11b).
To obtain a formula for the evolution of the charge density and the electric current j , we consider a medium that consists of N charged particles, at the positions x i : R ! R 3 , i D 1, : : : , N, as functions of time, with masses m i and charge densities q i i , where q i 2 R and i 2 C 1 c .R 3 / are non-negative functions with k i k L 1 D 1. Because these particles have physical radii smaller than 10 14 m, i can be considered an approximation of a ı-distribution.
Then, the charge density and the electric current j are given by the following: Clearly, they fulfil the continuity Eqn (12). For the motion of the particles, we assume that they move according to the Lorentz force: The Eqns (11), (13) and (14)

Macroscopic Maxwell's equations
However, solving the huge system of Eqns (14) is practically impossible. Therefore, we average the fields with respect to the space variable by some non-negative weighting function w 2 C 1 c .R 3 / with support around 0 and kwk L 1 D 1: Then, taking the average of the equation system (11) with respect to the function w, we obtain Maxwell's macroscopic equations for the averaged electromagnetic fields E and B depending on the averaged charge density N and the averaged electric current N j .

Remark 3.2
Because the medium usually consists of molecules, most of the particles are clustered around the centres of these molecules. Let us relabel the particles .
such that x i,k is the ith particle of the kth cluster of particles, i D 1, : : : , I k , k D 1, : : : , K, write q i,k for its electric charge and let N x k be the mass centre of the kth molecule.
Then, we introduce the zeroth-order approximations Q and J of the charge density N and the electric current N j , respectively, which are obtained from the averaged fields N and N j when the molecules are assumed to be point like, that is, setting x i,k D N x k for all i: We note that Q and J also fulfil the continuity equation We now write the deviation N Q of the averaged charge density from its zeroth-order approximation Q as the divergence of some function P : R R 3 ! R 3 , the so-called electric polarisation: According to Helmholtz's theorem, this is always possible but defines P only up to the addition of a divergence-free vector field. Then, we realise that we have with the continuity Eqns (12) and (17) and the definition (18) Doing the Helmholtz decomposition of the vector field N j J @ t P, we therefore find a function M : R R 3 ! R 3 , the magnetic polarisation, (uniquely defined up to a gradient vector field) so that These decompositions (18) and (19) of N and N j into macroscopic ( Q and J ) and polarisation (P and M ) parts are standard in the physics literature to formulate Maxwell's macroscopic equations; see for example [40,Chapter 6.6].
Clearly, the equation system (15) is not enough to determine the averaged fields E and B, because we would still require a solution of the microscopic problem to obtain the averaged quantities N and N j . The common way around this is to impose heuristic relations between these functions and the fields E and B. A common choice for biological tissue is a linear relationship.
for some function 2 C 1 c .R R 3 / with . , x/ D 0 for all < 0, x 2 R 3 a non-magnetic, isotropic and linear dielectric medium. For a non-magnetic, isotropic and linear dielectric medium, we may solve the equation system (15) for given parameter for E and B.

Remark 3.4
In the standard notation introduced in Remark 3.2, in particular (18) and (19), the assumptions corresponding to (20) are such that no magnetic polarisation exists in the medium, that is, M D 0, and that the macroscopic current J and the electric polarisation P are of the form

Problem formulation
Initially, the electromagnetic fields should describe a laser pulse, which did not yet interact with the medium. Thus, we want to look for a solution .E , B/, which for t < 0 is the superposition of a solution E .0/ , B .0/ of the vacuum equations (that is a solution of (15) with N D 0 and N j D 0) with support outside the medium (that is outside the support of . , / for every 2 R) and a stationary solution .E s , B s / produced by the undisturbed medium (we want to assume that without exterior influence, the electric and magnetic fields generated by the medium do not vary in time).

Proposition 3.5 Let
R 3 be an open set. For given functions s 2 C 1 c .R 3 / and j s 2 C 1 c .R 3 ; R 3 / with supp s , supp j s , and div j s D 0, we define the stationary electromagnetic fields E s , B s 2 C 1 .R 3 ; R 3 / by the following: Moreover, let E .0/ be a solution of the wave equation Then, there exists a solution E , B of Maxwell's macroscopic equations (15) for a non-magnetic, isotropic and linear dielectric medium with the parameter , as in Definition 3.3, and with charge density In particular, E is the solution of the differential equation together with the condition (26), which is equivalent to the integral equation where N j is given by (20).

Proof
We first check that (26) and (27) Finally, we obtain Because div x E .0/ D 0, we obtain with the identity curl x curl x F D grad x div x F x F for every function F 2 C 2 .
It remains to check that the current N j , defined by relation (20), indeed fulfils N j .t, x/ D j s .x/ for t Ä 0. Using that .t, x/ D 0 for t Ä 0, condition (24) and (25), we find Z 1 Thus, we have shown that (26) and (27) solve Maxwell's equations (15) for t Ä 0. Combining Maxwell's macroscopic equations (15c) and (15d), it follows that E in particular solves the equation To derive the integral equation, we use that E E s , B B s is a solution of (15) with N replaced by N s and N j by N j j s . Therefore, according to the representation (A.1) of the solution derived in Proposition A.1, we find for E that

Measurements
To simplify the calculations, we assume that the medium satisfies E s D 0 and B s D 0 in a domain . Next, we want to model the measurements of optical coherence tomography. We consider illuminating waves of the form with a function f 2 C 1 c .R/ and a fixed polarisation vector Á 2 R 2 f0g, where we choose f such that E .0/ fulfils the condition (24). Moreover, we place point detectors everywhere on the detector surface with d > 0 sufficiently large, in particular, we require that E .0/ .t, / D 0 for t 0 and 2 D. For more details, we refer to [11,Eqn (29]. Then, the reference wave E z produced by reflecting this incoming wave at a perfect mirror placed in the plane given by the equation x 3 D z is given by the following: The measurements carried out in optical coherence tomography are now the intensities on the detector surface D.
We combine these measurements to where the second term can be obtained by measuring the back-scattered wave without superimposing the reflected wave E z and the last term is explicitly known from the initial laser pulse E .0/ . Equation (34) Inserting the explicit formula (32) for the reflected wave E z , we obtain We use the convention for the Fourier transform of an integrable function F 2 L 1 .R/ with respect to time, and we put a subindex at F if we need to specify the variable with respect to which we do the Fourier transform.
Taking the inverse Fourier transform with respect to z in (35) and using Plancherel's formula we obtain 2 c Provided that the Fourier transform of the function f is nowhere zero and that the polarisation Á fulfils Á 1 ¤ 0 and Á 2 ¤ 0, we obtain from this the data where h can be directly calculated from the measurements Q I j and the knowledge of the initial wave E .0/ : The inverse problem of optical coherence tomography is now to find the material parameter : R R 3 ! R from the function h : R D ! R 2 according to the system of equations At first, it seems that we have enough data to solve this problem: we have the three-dimensional function h for every choice of function f and polarisation Á. However, because the problem is linear in the electric field E , we do not gain any information by changing f or having more than two linearly independent polarisations Á. To see this, we perform a Fourier transform F with respect to time.
Moreover, let E .k/ , k D 1, 2, 3, be the solutions of Then, the measurement data fulfil the relation where the coefficients c k 2 R are determined by Á .3/ D P 2 kD1 c k Á .k/ . Proof As in Proposition 3.5, we consider the equivalent integral equation (29) for E s D 0 Applying a Fourier transform with respect to time to the previous equation and using (20) in the frequency domain, we obtain that Now, the linearity of this equation implies that the solution F t E .3/ of the third equation can be written in the form Restricting this relation to the detector surface D, we obtain (40).
In particular, even if one considers the assumptions of Born and far-field approximations (specified later), we see that the measurements are equivalent to the determination of the four-dimensional Fourier transform of the function on a cone in R 3 ; see [11,Proposition 5.2]. Thus, also for OCT, we do not have sufficient data to uniquely recover the material parameter .
However, under certain simplifications, there exist works that provide reconstructions. A main assumption is that the medium is nondispersive, meaning that the temporal Fourier transform of does not depend on frequency. Under this assumption, Mark et al. [41] proposed algorithms under the Born approximation, and in [42,43], the one-dimensional case was considered. Also in [9,44], the OCT system is described using a single backscattering model.

Combined system
We have seen that the inverse problems of quantitative photoacoustic imaging and OCT both lack data to obtain a unique reconstruction of the material parameters depending on the specific modelling. Because the illumination for photoacoustic imaging and OCT is modelled analogously, they rely on the same electromagnetic material parameters of the medium. Therefore, we are suggesting to use the two measurements in combination to obtain additional informations.
To this end, we want to rewrite the optical problem in photoacoustics in terms of the material parameter (20) used in Maxwell's equations for the modelling of OCT instead of the absorption and scattering coefficients appearing in the radiative transfer equation. Thus, we need an expression of the absorbed energy of an electromagnetic wave in an dielectric medium.
We summarise the variables, which are used throughout this section, as a backup reference in Table IV.

Absorbed energy
We therefore start again with Maxwell's microscopic equations (11), (13) and (14) to describe the medium and define the averaged mechanical energy as the sum of the kinetic energies of the particles and their electrostatic potential energy: Assuming now that no particles enter the support of w. x/ or leave the region where w is constant, the change in energy is given by the following: We consider the limit when the charge distributions i of the single particles tend to ı-distributions, then the Lorentz force (14) is given by the following: Using the identity (43), and the fact that We further assume that the particles are moving relatively slow so that we may approximate the electric field e by the electrostatic approximation, which can be represented via (A.2) by neglecting the contribution of the current j : where e .0/ denotes the solution of the homogeneous wave equation with initial data e .0/ .0, x/ D e.0, x/ and @ t e .0/ .0, x/ D curl x b.0, x/. Here, we have used that for t sufficiently large, B ct .0/ will contain the complete medium and (13a) where i is replaced by a ı-distribution.
Then, we can write the change in energy in the form Until now, we had no particular assumptions on the size of the support of w. We have seen that the particles have diameter of order 10 14 m. The size of molecules in biological tissue is of order 10 10 m; a water molecule for example has a diameter of 0.3 nm. Thus, we specify r to be of order 10 10 m. In what follows, we assume in addition to Practically, the typical wavelength of photoacoustic imaging and OCT is around 1000 nm, which is three orders of magnitude larger than the size of molecules. Then, a ball with radius Q r of order 10 6 m will denote the resolution limitation, and we can assume that Q r r. For the averaging here, we assume that e .0/ is almost constant on the support of w. x/. Then, we may approximate The second term in (44) is nonzero only when x i 2 B r .x/ and x j … B r .x/. We want to assume that locally (related to resolution limit), there is almost no energy exchange: For the case, where x i 2 B r .x/ and x j … B Q r .x/, because Q r r, we can approximate x i x j ' y x j , for every y 2 B r .x/, and we thus obtain Now, because the support of w is small compared with the medium, we may take in the second term the sum over all particles x j without introducing a large error and find that Thus, in a non-magnetic, isotropic and linear dielectric medium, we have that the total absorbed energy is given by the following:

Remark 4.1
A different derivation of the absorbed energy follows from the definition of the total electromagnetic energy density and the energy transport of the electromagnetic wave (Poynting vector) as follows: respectively. Here, denotes the electric displacement. Using the elementary relation div x .F G / D G curl x F F curl x G , the Maxwell's equations (15c) and (15d) follow that Then, from the definition of W, integrating the previous equation over a subvolume of the medium and applying the divergence theorem give where N is the outward pointing unit normal vector. Eqn (50) describes the balance of energy: by the definition of W, the left hand side quantifies the variation of the total energy in the volume. The first term on the right-hand side describes the variation of the energy across the boundary, and the last two terms describe the dissipated power. According to [40], the integrands of the last two terms in the right-hand side integrated over time denote the electric energy density, which reads as follows:

Inverse problem
Using the expression (47) for the absorbed energy as a replacement of a .x/ N .x/ in (6), we obtain by combining (6) and (38) the inverse problem of calculating and from the measurement data p .0/ and h for giving initial pulse f and polarisation Á obeying the equation system Now, although varying the initial pulse f does not give us additional information with respect to the measurements of OCT, Lemma 3.6, its influence in the photoacoustic part is nonlinear and provides us with independent data. Thus, by choosing a one-parameter family of functions f , for example, almost monochromatic waves f centred at a frequency for arbitrary > 0, we obtain for every > 0 the equation system Now, the first two equations uniquely determine E as a function of . Then, depending on , we would have to solve the third equation for the four-dimensional function , which would leave us with two three-dimensional equations for the three-dimensional Grüneisen parameter . Although the dimensions of the data seem to perfectly fit for a reconstruction, it is still an open problem if these equations uniquely determine the material parameters; although there are partial results in this direction, see [17], where the injectivity of a very similar system is discussed.
As a plausibility argument, we want to show in the next section that at least in a simplified case, the parameters can be uniquely recovered.

Weakly scattering medium and multiple pulsed laser illuminations
The full system (51) can be equivalently transformed to a Fredholm integral equation of the second kind under the following assumptions: (i) Born approximation: the medium is weakly scattering, meaning that F t is sufficiently small. (ii) Far-field approximation: typically, in an OCT system, the measurements are performed in a distance much bigger compared with the size of the medium. (iii) Specific illumination: the support of the initial pulses tends to a delta distribution such that the optical parameter can be assumed constant in this spectrum.
In a first step, we derive a Fredholm integral equation of the first kind for the unknown Grüneisen parameter.

Proposition 5.1
The complete system (51) under the previous three assumptions is reduced to the integral equation where KOEQ p. , # ; y/ D We start with (51a) and (51b). We recall Eqn (41) derived in the proof of Lemma 3.6, then we obtain that holds for every frequency > 0. Under the Born approximation, we replace F t E in the integral by the incident field, to obtain E .b/ : We set x D R# , R > 0 and # 2 S 2 . Considering now the far-field approximation, we can replace the previous expression by its asymptotic behaviour for R ! 1, uniformly in # . This gives us, see for example [11,Equation (4.1) in Chapter 4.1]: Then, for every # 2 Q D D fÂ 2 S 2 j 9R > 0 : RÂ 2 Dg and 2 R C , the previous equation using Eqn (51d) and the definition (54) can be rewritten as follows: where we assumed that Ff . / ¤ 0. Similarly, we proceed with the Eqn (51c) of the photoacoustic tomography (PAT) measurement, written as follows: From Plancherel's formula, it follows that For the last equation, we used that p .0/ is real valued. Considering again the Born approximation, we obtain To take advantage of the third assumption, we choose initial pulses so that Ff has only support on a domain f! 2 R j j!j 2 OE ", C "g with a sufficiently small " > 0 such that F t . , x/ can be assumed to be for every x 2 R 3 constant on this support (we remark that because f and are a real-valued functions, the real and imaginary parts of their Fourier transforms have to be even and odd, respectively). If we additionally normalise R 1 1 jFf .!/j 2 d! D 1 2 , we find that Moreover, because is a real-valued function, we know that its Fourier transform fulfils the Kramers-Kronig relation [45] =.F t .!, Thus, combining (56), (57) and (58), we obtain the integral equation (52) for the function 1 where the kernel depends on the data Q p obtained from the PAT measurements.

Remark 5.2
(i) Compared with [11,Formula (20)], we have to mention that the different term in front of the integral in Eqn (55) is due to the different scaling between the coefficients and ; see the definition of P. (ii) In practice, we do not have Eqn (56) for every 2 R C because the band-limited source allows measurements only for a fixed frequency spectrum. (iii) In the most commonly used formulas, the imaginary part of the permittivity or the conductivity appears in (57). The differentiation of our result is due to the definition of , Remark 3.4, and should not confuse the reader.

Analytical formulas for special material
The particular form of the kernel K, Eqn (53), motivated us to transform the integral equation (52) of the first kind to one of the second kind. More precisely, in the special case, where the object mainly consists of a single material with the frequency dependence given by some function˛and a density distributionˇ, that is, if the initial pressure is of the form for some small function ", we may write the integral equation (52) considering (53) as follows: If we assume that A. / ¤ 0, then we may rewrite this as a Fredholm integral equation of the second kind for the Fourier transform y/ .y/ e ihk,yi dy of the functionˇ : where the kernel O K is for all values k 2 R 3 , which are of the form k D c .# C e 3 /, # 2 Q D, defined by the following: KOE". , # ; y/e ihÄ,yi dy and O h is similarly defined by the following: In particular, we see that for a sufficiently small function ", the integral equation for can have at most one solution [46,47].

Relation between the material parameters
In this section, we want to discuss the connection between the parameter , defined in (20), and the absorption and scattering coefficients, a and s , appearing in the radiative transfer equation (3). Recently, there has been an established connection between the radiative transfer equation (vector and scalar) and equations in classical electromagnetism [48,49]. However, even under simplifying assumptions, these results do not lead to a direct connection between the material parameters , a and s as the following remark shows.

Remark 6.1
This remark is based on classical results [50,51] and follows [49]. Let the medium contain N well separated scattering particles. We assume that there exist balls B k , k D 1, ..., N containing only one particle with large radius. We define by N k the unit outward normal vector on the boundary @B k . Let us denote by F C and F the outgoing (scattered) and incoming (incident) field for a given particle, respectively. If we neglect the interference effects between the incoming and outgoing fields, we can set as incident on the subdomain B k the field where S is defined in Remark 4.1. Then, E C k D E E k . We redefine the time averaged Poynting vector (the real part of the complex Poynting vector) S˛. , x/ :D c 8 <e F t E˛. , x/ F t B˛. , x/ Á ,˛D , C, describing the average flux of the incident and the scattered field, respectively. Then, the averaged rate of the energy incident on kth particle is given by S k . / :D 1 Bk R Bk jS . , x/j dx. Then, because every domain contains only one particle, the coefficients obtain the forms a,k . / D If the assumptions of the previous section still hold, then we can approximate the absorbed power as in the derivation of (57) resulting to a,k . / ' Observing the previous relations, we can identify the connection between the real part of the Fourier transform of and the absorption coefficient. But still arises the question if the obtained absorption and scattering coefficients are the ones that appear in the RTE; see Section 2. The fact that they depend on the incident average rate results to an indirect dependence on the initial illumination characterising the coefficients as non material parameters. If we omit the Born approximation, the dependence on the electric fields is stronger. In addition, the assumptions on the inner structure of the medium make the previous formulas applicable only in special cases.

Conclusions
To our knowledge, this paper is the first to use Maxwell's equations to model PAT. So far, the Maxwell's equations were only considered for modelling thermoacoustic tomography where low-frequency radiation is used. However, because recently, the connection between the radiative transfer equation and classical electromagnetism has been established [48,49], there is no need to differentiate from the proposed model considering the different modalities photoacoustic or thermoacoustic. Additionally, in our model, the frequency dependence of the optical parameter is not neglected, and thus, the effect of the different incident frequencies to the medium is included. This modelling allows to develop a unifying model for OCT and PAT, and thus, we can describe the combined setup.
Subtracting the solution e .0/ of the homogeneous wave equation with the initial conditions e .0/ .0, x/ D e.0, x/ and @ t e .0/ .0, x/ D curl x b.0, x/ 4 c j .0, x/, we find that e e .0/ solves the inhomogeneous wave equation with zero initial data. So, by Duhamel's principle, for example, [52,Chapter 2.4,Theorem 4], the solution is given by the following: