{'fraction_non_alphanumeric': 0.06845478080040741, 'fraction_numerical': 0.018758222637185417, 'mean_word_length': 3.6476331360946745, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Singular Values of the Attenuated Photoacoustic Imaging Operator

We analyse the ill-posedness of the photoacoustic imaging problem in the case of an attenuating medium. To this end, we introduce an attenuated photoacoustic operator and determine the asymptotic behaviour of its singular values. Dividing the known attenuation models into strong and weak attenuation classes, we show that for strong attenuation, the singular values of the attenuated photoacoustic operator decay exponentially, and in the weak attenuation case the singular values of the attenuated photoacoustic operator decay with the same rate as the singular values of the non-attenuated photoacoustic operator.


Introduction
In standard photoacoustic imaging, see e.g. [20], it is assumed that the medium is non-attenuating, and the imaging problem consists in visualising the spatially, compactly supported absorption density function h : R 3 → R, appearing as a source term in the wave equation ∂ tt p(t, x) − ∆p(t, x) = δ (t)h(x), t ∈ R, x ∈ R 3 , p(t, x) = 0, t < 0, x ∈ R 3 , ((1.1)) from measurements m(t, x) of the pressure p for (t, x) ∈ (0, ∞) × ∂Ω, where ∂Ω is the boundary of a compact, convex set Ω containing the support of h.
In this paper, we consider photoacoustic imaging in attenuating media, where the propagation of the waves is described by the attenuated wave equation where A κ is the pseudo-differential operator defined in frequency domain by: for some attenuation coefficient κ : R → C which admits a solution of (1.2). Heref denotes the one-dimensional inverse Fourier transform of f with respect to time t, that is, for f ∈ L 1 (R): The attenuated photoacoustic imaging problem consists in estimating h from measurements m of p on ∂Ω over time. The formal difference between (1.2) and (1.1) is that the second time derivative operator ∂ tt is replaced by a pseudo-differential operator A κ . We emphasise that standard photoacoustic imaging corresponds to κ 2 (ω) = ω 2 .
We associate with this solution the time-integrated photoacoustic operator in frequency domain: e iκ(ω)|x−y| |x − y| h(y) dy.
One goal of this paper is to characterise the degree of ill-posedness of the problem of inverting the time-integrated photoacoustic operator by estimating the decay rate of its singular values. We mention however, that although the attenuated photoacoustic operator, giving the solutionp, is related to the integrated photoacoustic operator by just time-differentiation, the singular values and functions of the photoacoustic operator have not been characterized so far.
In this paper, we are identifying two classes of attenuation models (classes of functions κ), which correspond to weakly and strongly attenuating media. We prove that for weakly attenuating media the singular values ((λ n (P * κPκ )) 1 2 ) ∞ n=1 decay equivalently to n − 1 3 , as in the standard photoacoustic imaging case, where this result has been proven in [16]. For the strongly attenuating models, the singular values are decaying exponentially, which is proven by using that in this case the operatorP * κPκ is an integral operator with smooth kernel.

The Attenuated Wave Equation
To model the wave propagation in an attenuated medium, we imitate the wave equation for the electric field E : R × R 3 → R 3 in an isotropic linear dielectric medium described by the electric susceptibility χ : R → R (extended by χ(t) = 0 for t < 0 to negative times): or written in terms of the inverse Fourier transformsĚ andχ with respect to the time: − ω 2 c 2 (1 +χ(ω))Ě(ω, x) − ∆Ě(ω, x) = 0. ((2.1)) Analogously, we want to incorporate attenuation by replacing the second time derivatives in our equation (1.1) by a pseudo-differential operator A κ of the form (1.3) for some function κ : R → C (corresponding to ω c 1 +χ(ω) in the electrodynamic model).
We will interpret the equation (1.2) as an equation in the space of tempered distributions S (R × R 3 ) so that the Fourier transform and the δ-distribution are both well-defined. To make sense of A κ as an operator on S (R × R 3 ) and to be able to find a solution of (1.2), we impose the following conditions on the function κ. (iii) we have the symmetry κ(−ω) = −κ(ω) for all ω ∈ R.
The condition (i) in Definition 2.1 ensures that the product κ 2 u of κ 2 with an arbitrary tempered distribution u ∈ S (R) is again in S (R) and therefore, the operator A κ is well-defined.
Remark: We use F when we are talking of the Fourier transform as an operator and use in the calculationŝ φ = Fφ andφ = F −1 φ. We will use the notation F also for the Fourier transform on different spaces (in particular for the three-dimensional Fourier transform on S(R 3 )).
The condition (ii) in Definition 2.1 is motivated by the fact that the functionχ in the electrodynamic model (2.1) is the inverse Fourier transform of a function whose support is inside [0, ∞) and can therefore be holomorphically extended to the upper half plane. We will see later, see Proposition 2.6, that this condition guarantees that the attenuated wave equation (1.2) has a causal solution in S (R × R 3 ), that is a solution whose support is contained in [0, ∞) × R 3 .
Finally, the condition (iii) in Definition 2.1 is required so that the attenuation operator A κ maps real-valued distributions to real-valued distributions: To see this, let u ∈ S (R × R 3 ) be a real-valued distribution. Then, for two real-valued functions φ ∈ S(R) and ψ ∈ S(R 3 ), the relation (2.3) implies By substituting the variable ω by −ω in the Fourier integral below, we get that with κ r given by κ r (ω) = κ(−ω). Thus, the condition A κ u, φ ⊗ ψ S ,S = A κ u, φ ⊗ ψ S ,S is equivalent to κ 2 = κ 2 r . Besides the case of a constant, real function κ (something we are not interested in), this is equivalent to κ = −κ r because of the condition mκ(z) ≥ 0 for all z ∈ H.

Solution of the Attenuated Wave Equation.
In this section, we want to determine the solution p ∈ S (R × R 3 ) of (1.2). To this end, we do a Fourier transform of the wave equation and end up with a Helmholtz equation for each value ω ∈ R, which in the case m κ(ω) > 0 has a unique solution in the space of tempered distributions.

Lemma 2.3
Let κ be a complex number with positive imaginary part, that is κ ∈ H, f ∈ L 2 (R 3 ) with compact essential support. Then, the Helmholtz equation has a unique solution u ∈ S (R 3 ), which is explicitly given by Proof: Writing φ = F −1φ , where F : S(R 3 ) → S(R 3 ) denotes the three-dimensional Fourier transform, we find with the function ψ ∈ S(R 3 ) defined by ψ = κ 2 φ + ∆φ, and thereforeψ(k) = Fψ(k) = (κ 2 − |k| 2 )φ(k), that The inner integral is the inverse Fourier transform of a product and can thus be written as the convolution of two inverse Fourier transforms: ((2.6)) Using spherical coordinates, we obtain by substituting ρ = |k| and cos θ = k,y |k||y| that Extending the integrand on the right hand side to a meromorphic function on the upper half complex plane, we can use the residue theorem to calculate the integral and find by taking into account that κ ∈ H that Inserting (2.7) into (2.6) and further into (2.5), and remarking that ψ is indeed an arbitrary function in S(R 3 ), we end up with (2.4).
To translate the initial condition in (1.2) that the solution p vanishes for negative times into Fourier space, we use that the Fourier transform of such a function can be characterised by being polynomially bounded on the upper half complex plane away from the real axis.
We will briefly summarise the theory as we need it. For a detailed exposition, we refer to [9,Chapter 7.4].
We define the adjoint Fourier-Laplace transformǔ : H → C of u by choosing for every point z ∈ H an arbitrary η z ∈ (0, m z) and by settinǧ Hereẽ z denotes for every z ∈ C with e z < 0 an arbitrary extension of the function e z | [0,∞) to the negative axis such thatẽ z ∈ S(R). Note that if u is a regular distribution: , then this is exactly the holomorphic extension of the inverse Fourier transform of U to the upper half plane: With this construction, we have that the inverse Fourier transform of the tempered distribution u η = e −η u for η > 0 is the regular distribution corresponding toǔ(· + iη), that is, Now, the causality of a distribution u, that is, supp u ⊂ [0, ∞), can be written in the form of a polynomial bound on its Fourier-Laplace transformǔ.
Lemma 2.5 We use again for every z ∈ C the notation e z ∈ C ∞ (R) for the function e z (t) = e izt .
(i) Let u ∈ D (R) be a distribution with supp u ⊂ [0, ∞) and such that e −η u ∈ S (R) for every η > 0.
Then, we find for every η 1 > 0 constants C > 0 and N ∈ N such that the adjoint Fourier-Laplace transformǔ of u fulfils (ii) Conversely, if we have a holomorphic functionǔ : H → C such that there exist for every η 1 > 0 constants C > 0 and N ∈ N with thenǔ coincides with the adjoint Fourier-Laplace transform of a distribution u ∈ D (R) with supp u ⊂ [0, ∞) and the property that e −η u ∈ S (R) for all η > 0. Proof: (i) Let η 1 > 0 and η 0 ∈ (0, η 1 ) be arbitrary. We choose a function ψ ∈ C ∞ (R) with ψ(t) = 1 for t ∈ (−∞, 0] and ψ(t) = 0 for t ∈ [1, ∞). Then, we write the given distribution u in the form u = u 1 + u 2 by setting Since we have by assumption e −η0 u ∈ S (R) and since u 1 has by construction compact support, we get that e −η0 u 2 ∈ S (R). Thus, because of supp Moreover, since e −η0 u 1 has compact support supp(e −η0 u 1 ) ⊂ [0, 1], we find, see for example [9, Theorem 2.3.10], constants A 1 > 0 and N 1 ∈ N so that We now define as in Definition 2.4 for z ∈ C with e z > 0 an extensionẽ z ∈ S(R) of the function e z | [0,∞) and choose for every z ∈ H the function φ =ẽ iz+η0 in (2.11) and (2.12). Then, there exists a constant C > 0 such that with N = max{N 1 , (ii) To construct the distribution u, we define from the given functionǔ for every η > 0 the distribution u η ∈ S (R) via the relation (2.8), so that the inverse Fourier transform of u η is given by the regular distribution corresponding to the function ω →ǔ(ω + iη).
Now, we want to show that e η u η is in fact independent of η, so that there exists a distribution u such that u η = e −η u for every η > 0. To do so, we first remark that the derivative ∂ η u η of u η with respect to η fulfils for every φ ∈ S(R) whereφ = Fφ denotes the Fourier transform of φ and f (t) = −t, so thatφ = iF(f φ). Thus, ∂ η u η = f u η and therefore, ∂ η (e η u η ) = e η (∂ η u η − f u η ) = 0, proving that e η u η is independent of η.
So, the distribution u = e η u η ∈ D (R) is well-defined and fulfils by construction that e −η u ∈ S (R) for every η > 0.
Next, we want to show that supp u ⊂ [0, ∞). Let φ ∈ C ∞ c (R) and write againφ = Fφ. Then, by our construction, we have for every η 1 > 0 that Sinceφ is the Fourier transform of a function with compact support, we can extend it holomorphically to C and get for every N 1 ∈ N 0 a constant C 1 > 0 such that the upper bound holds.
Therefore, we can shift the line of integration in (2.13) by an arbitrary value η > 0 upwards in the upper half plane and get with the upper bound (2.10) that Choosing now φ such that supp φ ⊂ (−∞, 0) and taking the limit η → ∞, the right hand side tends to zero, showing that e −η1 u, φ S ,S = 0 whenever supp φ ⊂ (−∞, 0). Thus, e −η1 u, and therefore also u, has only support on [0, ∞).
Finally, we verify that the Fourier-Laplace transform of u is given byǔ. Indeed, given any ω ∈ R and η > 0, we have by construction for every η 1 ∈ (0, η) and every extensionẽ z ∈ S(R) of e z | [0,∞) for z ∈ C with e z > 0 that Since supp u ⊂ [0, ∞), we know that this expression is independent of the concrete choice of the extensionẽ z . Moreover, both sides are independent of η 1 . Thus, letting on the right hand sideẽ z converge to e z and η 1 to η, Fẽ iω+(η1−η) will tend to √ 2π times the δ-distribution at ω, and we therefore get 1 √ 2π e −η1 u,ẽ i(ω+iη)+η1 S ,S =ǔ(ω + iη).
We now return to the solution of the attenuated wave equation (1.2).
Proposition 2.6 Let κ be an attenuation coefficient and A κ : S (R × R 3 ) → S (R × R 3 ) be the corresponding attenuation operator. Let further h ∈ L 2 (R 3 ) with compact essential support.
Then, the attenuated wave equation where the Laplace operator ∆ : Moreover, p is of the form whereφ denotes the Fourier transform of φ and G denotes the integral kernel Proof: Let p ∈ S (R×R 3 ) be a solution of (2.14) with supp p ⊂ [0, ∞)×R 3 . We evaluate the equation (2.14) for ϑ = φ ⊗ ψ ∈ C ∞ c (R × R 3 ) and write Fφ =φ. It then follows that For arbitrary z ∈ H, we define the adjoint Fourier-Laplace transformp(z) ∈ S (R 3 ) of p by whereẽ z ∈ S(R), z ∈ C with e z < 0, is an arbitrary extension ofẽ z (t) = e zt for t ≥ 0, see Definition 2.4. Then, z → p(z), ψ is holomorphic in the upper half plane H and we have We replace κ in the integrand now by its holomorphic extensionκ : H → H, see (ii) in Definition 2.1, and also extend the Fourier transformφ of the compactly supported function φ holomorphically to C.
Since z → p(z), ψ is the adjoint Fourier-Laplace transform of a distribution with support on [0, ∞), it is polynomially bounded, see Lemma 2.5. Moreover, we have by Definition 2.1 of the attenuation coefficient a polynomial bound onκ and get therefore with the dominated convergence theorem that Since all functions in the integrand are holomorphic in the upper half plane the integral is independent of η and we can therefore remove the limit with respect to η. Using again the dominated convergence theorem, we can evaluate now the limit with respect to ξ and obtain for arbitrary η > 0 the equality Inserting this into the equation (2.17) and arguing in the same way for the two other terms therein, we see thatp(z) ∈ S (R 3 ) solves for every z ∈ H the equatioñ Thus, by Lemma 2.3, we get for every z ∈ H with mκ(z) > 0 that is the only solution. However, sinceκ is holomorphic, its imaginary part cannot vanish in any open set unlessκ were a constant, real function which is excluded by the symmetry condition (iii) in Definition 2.1. Therefore, we can uniquely extend the formula (2.18) for p(z), ψ S ,S by continuity to all z ∈ H.
It remains to verify that supp p ⊂ [0, ∞) × R 3 . To see this, we use that mκ(z) ≥ 0 for every z ∈ H to estimate the integral in (2.18) by p(z), ψ S ,S ≤ C|z| with some constant C > 0. Therefore, by Lemma 2.5, p(z), ψ S ,S is the Fourier-Laplace transform of a distribution with support in [0, ∞).

Finite Propagation Speed.
Seeing the equation (2.14) as a generalisation of the wave equation, it is natural to additionally impose that the solution propagates with finite speed.

Definition 2.7
We say that the solution p ∈ S (R × R 3 ) of the equation (2.14) propagates with finite speed We can give an explicit characterisation of the equations whose solutions propagate with finite speed in terms of the holomorphic extensionκ of the attenuation coefficient κ.
Lemma 2.8 The solution p of the attenuated wave equation (2.14) propagates with finite speed c > 0 if the holomorphic extensionκ of the attenuation coefficient κ fulfils Conversely, if there exists a sequence (z ) ∞ =1 ⊂ H with the properties that • there exists a parameter η 1 > 0 such that m(z ) ≥ η 1 for all ∈ N, • we have |z | → ∞ for → ∞, and • there exists a parameter δ > 0 such that then p propagates faster than with speed c.
Proof: Since the solution p ∈ S (R × R 3 ) is a regular distribution with respect to the second component, see (2.15), having finite propagation speed is equivalent to the condition that the distribution p(x) ∈ S (R), given by Letting h tend to a three dimensional δ-distribution, we see that the distribution g(x) ∈ S (R), defined by , this means that the distribution g τ (x) ∈ S (R), given by to the upper half plane using the explicit formula (2.16) for G κ , we obtain the adjoint Fourier- According to Lemma 2.5, we can therefore equivalently characterise a finite propagation speed in terms of a polynomial bound on the functionǧ τ (·, x).
• If m(κ(z) − z c ) ≥ 0 for every z ∈ H, then the adjoint Fourier-Laplace transform of g τ (x) fulfils that for every x ∈ R 3 there exists a constant C > 0 such that Thus, the condition (2.10) of Lemma 2.5 is satisfied and therefore supp g τ (x) ⊂ [0, ∞), so that p propagates with the finite speed c > 0.
so that condition (2.9) of Lemma 2.5 is violated and therefore the support of g τ (x) cannot be contained in [0, ∞).
where σ : R → R is a monotonically increasing function of bounded variation and A ≥ 0 and B ∈ R are arbitrary parameters, and vice versa, see [3,Formula 3.3].
Then,κ(z) − Az is still of the form (2.19) and therefore is a holomorphic function mapping H to H. In particular, it satisfies m(κ(z) − Az) ≥ 0 for all z ∈ H. Thus, if A > 0, p propagates with the finite speed c = 1 A according to Lemma 2.8. Evaluatingκ along the imaginary axis, we find that asymptotically as ω → ∞ Thus, Moreover for A = 0, we see from (2.20) that for every choice of c > 0, we have the behaviourκ(iω) − iω c = iω(− 1 c + o(1)) and therefore m(κ(iω) − iω c ) ≤ − ω 2c for all ω ≥ ω 0 for a sufficiently large ω 0 . Thus, by Lemma 2.8, p cannot have finite propagation speed for A = 0.

Examples of Attenuation Models
The following examples of attenuation coefficient have been collected in [11], where also references to original papers can be found. In this section, we review them and catalog them into two groups which are characterised by different spectral behaviour.
If the attenuation in the medium increases faster than some power of the frequency, we are in the case of strong attenuation.
A common example, which has the drawback of an infinite propagation speed, is the thermo-viscous model, see Table 1. In [12], the authors modified this model to obtain one with finite propagation speed, see Table 2. Other models, trying to match the heuristic power law behaviour of the attenuation are the Name: Thermo-viscous model, see for example [10,Chapter 8.2] Attenuation coefficient: Propagation speed:

Attenuation type: Strong attenuation coefficient
Indeed a Taylor expansion with respect to 1 ω around 0 yields for ω → ∞: Range ofκ: To see analytically thatκ maps the upper half plane H into itself, we first remark that because of the symmetryκ(−z) = −κ(z), it is enough to show that the first quadrant  Table 3 and Szabo's model, see Table 4, where we chose the modified version introduced in [11] as the original one does not lead to a causal model.

power law in
In these tables and in the following we always use the principal branch of the complex roots, that is, we define for γ ∈ C (re iϕ ) γ = e γ(log(r)+iϕ) for every r > 0, ϕ ∈ (−π, π).
Remark: The attenuation coefficients in Table 2, Table 3, and Table 4 do not fulfil the smoothness assumption κ ∈ C ∞ (R). However, this requirement originates mainly from our choice of solution concept for the attenuated wave equation (1.2) and we may still consider formula (4.1) as definition of the solution p of (1.2) if κ is non-smooth. In particular, the smoothness assumption is not required for the derivation of the decay of the singular values of the integrated photoacoustic operator.
In the case where the attenuation decreases sufficiently fast as the frequency increases, we call the medium weakly attenuating.

Model:
Kowar-Scherzer-Bonnefond model, see [12] Attenuation coefficient: Propagation speed: Attenuation type: Strong attenuation coefficient A Taylor expansion with respect to ω −γ around 0 yields for ω → ∞: Range ofκ: To see whereκ maps the upper half plane, we writeκ in the form And f 2 is a Möbius transform which maps Q +− to the half ballB 1 2 ( 1 2 ) ∩ Q ++ . Thus, since the product of two points in the first quadrant Q ++ is in the upper half plane,κ(Q ++ ) ⊂ H and because of the symmetryκ(−z) = −κ(z), we therefore havẽ κ(H) ⊂ H. for some constants c > 0 and κ ∞ ≥ 0 and a bounded function Clearly, the non-attenuating case where κ(ω) = ω c with c > 0, so that the attenuated wave equation (1.2) reduces to the linear wave equation, falls under this category and we will later on treat the case of a weak attenuation coefficient as a perturbation of this non-attenuated case.
A non-trivial example of a weak attenuation model is the model by Nachman, Smith, and Waag, see Table 5.

The Integrated Photoacoustic Operator
Let us now return to the attenuated photoacoustic imaging problem. Thus, we consider the operator mapping the source term h in the attenuated wave equation (1.2) (interpreted in the sense of (2.14)) Model: Power law, see for example [18] Attenuation coefficient: The second inequality uses Young's inequality to estimate |z| γ ≤ γ|z| + 1 − γ.
Propagation speed: Attenuation type: Strong attenuation coefficient Range ofκ: That the range ofκ is a subset of H follows immediately from to the measurements, which shall correspond to the solution of the attenuated wave equation on the measurement surface ∂Ω measured for all time.
According to Proposition 2.6, the solution p of the attenuated wave equation is given by (2.15). This means that the temporal inverse Fourier transform of p is the regular distribution corresponding to the functioň where G κ is defined by (2.16).
We therefore introduce our measurementsm as the functioň Instead of considering the operator mapping h tom, we will divide the data by −iω, meaning that we consider the map from h to the inverse Fourier transform of the measurements which were integrated over time. Additionally, we want to assume that the measurements are performed outside the support of the source.
Remark: This assumption that the absorption density functions h has compact support in the domain Ω is very common in the theory of photoacoustics, see for instance [2,13].
Model: Modified Szabo model, see [11] and, for the original version, [18] Attenuation coefficient: The second inequality uses Young's inequality to estimate |z| 1+γ 2 Propagation speed: Attenuation type: Strong attenuation coefficient A Taylor expansion with respect to ω γ−1 around 0 yields for ω → ∞: Range ofκ: To determine the range, we writeκ in the form Now, since γ − 1 < 0, f maps the first quadrant Q ++ to a subset of Q ++ . Thus, since the product of two points in the first quadrant is in the upper half plane, we havẽ κ(Q ++ ) ⊂ H and because of the symmetryκ(−z) = −κ(z) thereforeκ(H) ⊂ H. Then, we callP the integrated photoacoustic operator of the attenuation coefficient κ in frequency domain.
Name: Nachman-Smith-Waag model, see [15] Attenuation coefficient: Propagation speed: Attenuation Type: Weak attenuation coefficient Indeed a Taylor expansion with respect to 1 ω around 0 yields for ω → ∞: Range ofκ: for z ∈ Q ++ and the symmetryκ(−z) = −κ(z) then implies thatκ(z) ∈ H for every z ∈ H. (ii) In the case of a weak attenuation coefficient, see Definition 3.2, we split the operator intoP κ = P as the photoacoustic operator with constant attenuation and the perturbatioň Now,P κ h is seen to be the inverse Fourier transform of the function P κ h, defined by can be directly seen to be a bounded linear operator, since we have, recalling that κ ∞ ≥ 0, Thus, combining the two inner integrals to an integral over Ω ε , we find that This is a special case of the more general result in [16,Lemma 4.1].
To obtain the singular values of the operatorP κ , we consider the operatorP * κPκ on L 2 (Ω ε ). It turns out that this operator is a Hilbert-Schmidt integral operator, in particular therefore compact. So, by the singular theorem for compact operators, its spectrum consists of at most countably many positive eigenvalues and the value zero. Proposition 4.3 LetP κ : L 2 (Ω ε ) → L 2 (R × ∂Ω) be the integrated photoacoustic operator of a weak or a strong attenuation coefficient κ for some bounded, convex domain Ω ⊂ R 3 with smooth boundary and some ε > 0. Then,P * κPκ : In particular,P * κPκ is a Hilbert-Schmidt operator and thus compact.
We remark that the convexity and smoothness assumptions on Ω are only needed for the weak attenuation case. For strong attenuation, a Lipschitz domain Ω is sufficient.
(ii) In the case of a weak attenuation coefficient, we write F κ similar to the proof of Lemma 4.2 as the sum of a contribution F κ : • To calculate the double integral in F (0) κ , we first parametrise ∂Ω depending on the values x, y ∈ Ω ε with x = y. We choose a positively oriented, orthonormal basis (e j ) 3 j=1 ⊂ R 3 with e 3 = y−x |y−x| and consider the curve given as the intersection of the boundary ∂Ω and the plane E ϕ through 1 2 (x + y), spanned by the vectors e 3 and cos ϕ e 1 + sin ϕ e 2 . Setting a ϕ = min where we pick for ϕ ∈ (0, π) the function r in such a way that r(ϕ, z) > −r(ϕ − π, z) so that the two maps ψ(ϕ − π, ·) and ψ(ϕ, ·) parametrise together Γ ϕ , see Figure 1.
• If we would formally interchange the order of integration in the definition (4.10) of F (0) κ , the integral over ω would lead to a δ-distribution at the zeros of the exponent 1 c (|ξ − x| − |ξ − y|). We therefore start by analysing this exponent, which is up to the prefactor 1 c given by if we use the parametrisation ψ for integrating over ∂Ω.
The zeros of g are exactly those points (ϕ, z) such that ψ(ϕ, z) is in the bisection plane of x and y. Thus, we have by construction of the parametrisation ψ, see (4.11), that g(ϕ, z) = 0 is equivalent to z = 0. ((4.12)) Furthermore, we can prove that only vanishes at the two points where ψ(ϕ, z) is the intersection point of the line through x and y with ∂Ω: We assume by contradiction that ∂ z g(ϕ, z) = 0 at a point (ϕ, z) ∈ U with ψ(ϕ, z) not lying on the line through x and y. Then, the first vector ψ(ϕ,z)−x |ψ(ϕ,z)−x| − ψ(ϕ,z)−y |ψ(ϕ,z)−y| would be a non-zero vector in E ϕ , and the second vector ∂ z ψ(ϕ, z) is by construction a non-zero tangent vector on Γ ϕ ⊂ E ϕ at ψ(ϕ, z). Thus, ∂ z g(ϕ, z) = 0 would imply that the first vector is a non-trivial multiple of the outer unit normal vector ν(ψ(ϕ, z)) to Γ ϕ at ψ(ϕ, z). However, if w 1 , w 2 ∈ R 2 are two unit vectors with w 1 − w 2 = n, n = 0, then w 1 , n = − w 2 , n , since for given w 1 ∈ S 1 , w 2 is the intersection point of S 1 with the line parallel to n through w 1 . Thus, we would have but, because of the convexity of Ω, the projections ψ−x |ψ−x| , ν • ψ and ψ−y |ψ−y| , ν • ψ have to be both positive, which is a contradiction. Therefore, ∂ z g(ϕ, z) = 0 if and only if ψ(ϕ, z) is on the line through x and y.
Plugging this together with formula (4.15) for µ into (4.16), we finally get for F (0) κ the representation This is an integral over the intersection γ x,y = {ξ ∈ ∂Ω | ξ − 1 2 (x + y), y − x = 0} of ∂Ω and the bisection plane of the points x and y. To express it in a parametrisation invariant form, we write it as a line integral over γ x,y , which means we want to use the volume element |∂ ϕ ψ(ϕ, 0)| dϕ. Calculating the outer unit normal vector explicitly from (4.11), we find Thus, we get the relation With this, we can write (4.17) as the parametrisation free line integral In particular, we have for all x, y ∈ Ω ε with x = y ((4.19)) for some constant A > 0, and therefore F (0) κ ∈ L 2 (Ω ε × Ω ε ). • To estimate the first perturbation F (1) κ , we remark that, according to Definition 3.2, κ * is bounded and square integrable. We can therefore pull in the definition (4.10) of F (1) κ the integration over the variable ω as an inverse Fourier transform inside the surface integral and find that Choosing now a radius R > diam Ω, we get by applying Hölder's inequality and increasing the domain of integration that Ωε Ωε Thus, switching in the two inner integrals to spherical coordinates around the point ξ, we find Ωε Ωε which shows that F κ ∈ L 2 (Ω ε × Ω ε ).
• The second perturbation F (2) κ can be bounded directly via for all x, y ∈ Ω ε . Using that for every bounded set D ⊂ C, there exists a constant C > 0 such that |e z − z − 1| ≤ C|z| 2 for all z ∈ D holds, we find a constantC > 0 so that Plugging this into (4.20), we get that Thus, in particular, we have F We conclude therefore that F κ = F , which shows as in the first part of the proof thatP * κPκ is a Hilbert-Schmidt operator and compact.

Singular Values of the Integrated Photoacoustic Operator
We have seen in Proposition 4.3 that the operatorP * κPκ , given by (4.9), is a compact operator. The inversion of the photoacoustic problem is therefore ill-posed. To quantify the ill-posedness, we want to study the decay of the eigenvalues (λ n (P * κPκ )) n∈N ofP * κPκ , where we enumerate the eigenvalues in decreasing order: 0 ≤ λ n+1 (P * κPκ ) ≤ λ n (P * κPκ ) for all n ∈ N. We differ again between the two cases of a strong and of a weak attenuation coefficient κ.

Strongly Attenuating Media.
To obtain the behaviour of the eigenvalues ofP * κPκ in the case of a strong attenuation coefficient κ, see Definition 3.1, we will use Corollary A.4 which gives a criterion for a general integral operator with smooth kernel to have exponentially fast decaying eigenvalues in terms of an upper bound on the derivatives of the kernel, see (A.24). We therefore only have to check that the kernel (4.8) ofP * κPκ fulfils these estimates. The calculations are straightforward (although a bit tedious) and can be found explicitly in Appendix B.
Then, the kernel F κ ofP * κPκ , explicitly given by (4.8), fulfils the estimate for some constants B, b > 0, where N ∈ N denotes the exponent for = 0 in the condition (2.2) and β ∈ (0, N ] is the exponent in the condition (3.1) for the strong attenuation coefficient κ.
Proof: Putting the derivatives with respect to s inside the integrals in the definition (4.8) of the kernel F κ , we get that where G κ denotes the integral kernel (2.16). We remark that the term 1 ω 2 comes from the fact that we consider the integrated photoacoustic operator instead of the operator which maps directly the measurements to the initial data.
Using Proposition B.3 to estimate the derivative of G κ , we find a constant C > 0 so that From the uniform estimate which is directly obtained from the definition (2.16) of G κ by using that |x − ξ| ≥ ε for all ξ ∈ ∂Ω and x ∈ Ω ε , it then follows that Applying now Lemma B.4 (for the first term in the integrand, we use Lemma B.4 with j = 0), we find constants B, b > 0 so that Combining Proposition 5.1 with Corollary A.4, we obtain the decay of the singular values of the integrated photoacoustic operator.
Then, there exist constants C, c > 0 so that the eigenvalues (λ n (P * κPκ )) n∈N ofP * κPκ in decreasing order fulfil Proof: According to Proposition 5.1, we know that there exist constants B, b > 0 so that the integral kernel F κ of the operatorP * κPκ fulfils the estimate (5.1). Applying thus Corollary A.4 with µ = N β − 1 to the operatorP * κPκ , we obtain the decay rate (5.2).

Weakly Attenuating Media.
To analyse the operatorP * κPκ in the case of a weak attenuation coefficient κ, see Definition 3.2, we splitP κ as in the proof of Lemma 4.2 inP κ =P (0) κ +P (1) κ , see (4.4) and (4.5). We will show that decomposing the operator asP * κPκ =P κ decay faster than those ofP κ so that Q κ does not alter the asymptotic decay rate of the eigenvalues ofP κ corresponds to a constant attenuation and its behaviour was already discussed in [16].

Lemma 5.3
Let κ be a weak attenuation coefficient, Ω ⊂ R 3 be a bounded, convex domain with smooth boundary, and ε > 0. We define the operatorP by (4.4). Then, there exist constants C 1 , C 2 > 0 such that we have for all n ∈ N.

((5.3))
Proof: The idea of the proof is to show that the operatorP κ has the same eigenvalues as an elliptic pseudofifferential operator T : L 2 (M ) → L 2 (M ) of order −2 on a closed manifold M . Then, we can apply the result [17,Theorem 15.2] to obtain the asymptotic behaviour of the eigenvalues.
• First, we want to replace the operatorP (0) κ * P (0) κ by a pseudodifferential operator T on a closed manifold with the same eigenvalues.
We have seen in the proof of Proposition 4.3 that the operatorP κ is an integral operator with integral kernel F (0) κ defined by (4.10). We now generate the closed manifold M by taking two copies of Ω ε and identifying their boundary points: M = (Ω ε × {1, 2})/ ∼ with the equivalence relation (x, a) ∼ (x,ã) if and only if x =x and either a =ã or x ∈ ∂Ω ε . This is called the double of the manifold with boundary Ω ε , see for example [14,Example 9.32]. Then, the operator T : L 2 (M ) → L 2 (M ) given by has the same non-zero eigenvalues asP To write T in the form of a pseudodifferential operator, we extend the kernel F (0) κ to a smooth functionF with φ(y) = 1 for y ∈ Ω ε and supp φ ⊂ Ω and setting Then, defining g up to the normalisation factor (2π) 3 2 as the inverse Fourier transform ofF (0) κ with respect to y: we can write the kernel F (0) κ with the Fourier inversion theorem in the form where g is smooth, since g(x, ·) is the Fourier transform of a function with compact support.
• In the expression (5.5) for g, the integral over ω from the definition (5.4) ofF (0) κ can be seen as a one-dimensional inverse Fourier transform, which allows us to get rid of two one-dimensional integrals.
To this end, we pull the outer integral over R 3 inside both other integrals and write it in spherical coordinates around the point ξ: y = ξ + rθ with r > 0 and θ ∈ S 2 . This gives us For every ξ ∈ ∂Ω and every θ ∈ S 2 , the two inner integrals with respect to r and ω each represent a Fourier transform and we get with ρ(r) = re −κ∞ (r+|ξ−x|) Thus, we find ((5.6)) • We are now interested in the leading order asymptotics of g(x, k) as |k| → ∞. To obtain this, we will apply the stationary phase method, see for example [9,Theorem 7.7.5].
So, let ψ ∈ C ∞ (U ; R 3 ) be a parametrisation of ∂Ω and Θ ∈ C ∞ (V ; R 3 ) be a paramtetrisation of S 2 with some open sets U, V ⊂ R 2 . Then, according to the stationary phase method, the asymptotics is determined by the region around the critical points of the phase function in the integrand in (5.6). The optimality conditions with respect to ϑ imply that k is normal to the tangent space of S 2 in the point Θ(ϑ) at a critical point (η, ϑ) of Φ x,k , that is Θ(ϑ) = ± k |k| . The optimality conditions with respect to η then imply for a critical point (η, ϑ) of Φ x,k that the projections of the two vectors − ψ(η)−x |ψ(η)−x| and Θ(ϑ) = ± k |k| on the tangent space of ∂Ω at ψ(η) coincide. Since both vectors have unit length, this means that up to the sign of the normal component they have to be equal. Additionally, we use that, because of the cut-off term φ(ψ(η) + |ψ(η) − x|Θ(ϑ)) in the integrand, the critical points at which the vector Θ(ϑ) points outwards the domain Ω at ψ(η) do not contribute to the integral. Therefore, a relevant critical point (η, θ) is such that Θ(ϑ) is pointing inwards at ψ(η) and since − ψ(η)−x |ψ(η)−x| is also pointing inwards, we are left with the two critical points (η ( ) , ϑ ( ) ) given by In particular, we have Φ x,k (η ( ) , ϑ ( ) ) = 0.
For the second derivatives of Φ x,k at the critical points, we find For the determinants of the derivatives with respect to η and ϑ, we obtain (this can be readily checked for parametrisations ψ and Θ corresponding to normal coordinates at the points (ψ(η ( ) ), Θ(ϑ ( ) ))) where ν : ∂Ω → S 2 denotes the outer unit normal vector field on ∂Ω.
Therefore, the stationary phase method, see for example [9, Theorem 7.7.5], implies for x ∈ Ω ε , k ∈ R 3 , and µ > 0 that we have asymptotically for µ → ∞ Thus, g is of the form with g −2 (x, ·) being a positive function which is homogeneous of order −2.
Therefore, a parameterix of the pseudodifferential operator T on L 2 (M ) is an elliptic pseudodifferential operator of order 2 and thus has, according to [17,Theorem 15.2], eigenvalues which grow as n 2 3 . Consequently, the eigenvalues of T and thus also those ofP * κPκ decay as n − 2 3 .
To estimate the eigenvalues of the termP (1) κ * P (1) κ inP * κPκ , we show with Mercer's theorem that the operator is trace class.
Taking now an arbitrary sequence (x k , y k ) k∈N ⊂ Ω ε converging to an element (x, y) ∈ Ω ε , we get with the dominated convergence theorem and the continuity of r κ that Thus, R κ is continuous and therefore Mercer's theorem, see for example [6, Chapter III, §5 and §9], implies that ∞ n=1 λ n (P (1) is by definition a decreasing sequence, Abel's theorem, see for example [8, §173], gives us (5.7).
From the decay rates of the singular values ofP Lemma 5.5 Let κ be a weak attenuation coefficient, Ω ⊂ R 3 be a bounded, convex domain with smooth boundary, and ε > 0.
Proof: We start with the positive semi-definite operatorP * κPκ . We splitP κ as in the proof of Lemma 4.2 inP κ =P (0) κ +P (1) κ withP (1) κ given by (4.5). Then, we can write Q κ in the form To estimate the eigenvalues of the operatorP κ , we use that for all m, n ∈ N the inequalities κ )) ∞ n=1 into these inequalities, we find that lim n→∞ n 5 6 |λ n (P (1) Estimating, again with [7, Chapter II.2.3, Corollary 2.2], the eigenvalues of the sum Q κ of the two operatorš P for all m, n ∈ N, which yields with the behaviour (5.7) of the eigenvalues ofP κ the result (5.9).
Combining Lemma 5.3 and Lemma 5.5, we obtain that the singular values of the photoacoustic operatoř P κ decay as in the unperturbed case.

A. Eigenvalues of Integral Operators of Hilbert-Schmidt Type
In this section, we derive estimates for the eigenvalues of operators T of the form on a bounded, open set U ⊂ R m with an Hermitian integral kernel F ∈ C(Ū ×Ū ) (that is, F (x, y) = F (y, x)). In particular, such an operator T is a self-adjoint Hilbert-Schmidt operator, and we want to additionally assume that T is positive semi-definite. So, the eigenvalues (λ n (T )) n∈N of T are non-negative and we enumerate them in decreasing order: To obtain the asymptotic decay rate of the eigenvalues of such an operator T , we proceed as in [4] where a characterisation for a decay rate of the form λ n (T ) = O(n −k ) was presented in terms of an upper estimate on the derivatives of the kernel F . The extension to an exponential decay rate is rather straightforward.
First, we show that when approximating an operator T 1 by a finite rank operator T 2 , we can estimate the eigenvalues above the rank of the finite rank operator T 2 in terms of the supremum norm of the difference of their kernels, see for example [22,Satz II].
Lemma A.1 Let U ⊂ R m be a bounded, open set, T i : L 2 (U ) → L 2 (U ), i = 1, 2, be two integral operators with Hermitian integral kernels F i ∈ C(Ū ×Ū ). Moreover, let T 1 be positive semi-definite and T 2 have finite rank r ∈ N 0 .
Then, the eigenvalues (λ n (T 1 )) n∈N of T 1 (sorted in decreasing order) satisfy Proof: The min-max theorem (see for example [7, Chapter II.2.3]) states that for every self-adjoint, compact operator T and every fixed m ∈ N |λ m (T )| = min where the minimum is taken over all operators A : L 2 (U ) → L 2 (U ) with rank less than or equal to m − 1.
Let us fix n ∈ N now. Applying (A.2) with T = T 1 − T 2 shows that there exists an operator A with rank(A) ≤ n − 1 such that Because rank(T 2 ) ≤ r and rank(A) ≤ n − 1, rank(T 2 + A) ≤ n + r − 1, and we therefore have Taking the sum over all n ∈ N, we get The eigenvalues of T 1 − T 2 do not need to be all non-negative, however, since T 2 has rank at most r, the operator T 1 − T 2 cannot have more than r negative eigenvalues. Moreover, their norm is bounded by Thus, we can estimate the sum in (A.6) by Moreover, choosing an orthonormal eigenbasis (ψ n ) ∞ n=1 ⊂ L 2 (U ) of the compact, self-adjoint operator T 1 − T 2 , we get with Mercer's theorem, see for example [6,Chapter III,§5 and §9], that ( Thus, approximating the kernel in our integral operator by one of its Taylor polynomials, we get a convergence rate for the eigenvalues depending on the approximation error of the Taylor polynomial. To improve this estimate, we first subdivide the domain U in smaller domains, so that the approximation error of the Taylor polynomial is smaller.
Regarding an upper bound for the eigenvalues, it is indeed enough to keep the subdomains along the diagonal of U × U , see [4,Lemma 1].
Then, the integral operator T 2 : L 2 (U ) → L 2 (U ) with the integral kernel F 2 is also positive semi-definite and fulfils Proof: For each ∈ {1, . . . , N }, let P : L 2 (U ) → L 2 (U ), P h = hχ Q be the orthogonal projection onto the subspace L 2 (Q ). In particular, because the sets Q are pairwise disjoint, we have P k P = δ k, P . ((A.11)) With this notation, we can write Now, we first show that T 2 is indeed positive semi-definite. Let us assume by contradiction that this is not the case. Then, there exists a function h ∈ L 2 (U ) so that h, T 2 h < 0. Thus, because of the representation (A.12) of T 2 , we find an index ∈ {1, . . . , N } such that P h, T 1 P h = h, P T 1 P h < 0.
However, this contradicts the fact that T 1 should be positive semi-definite.
To get a relation between the eigenvalues of T 1 and T 2 , we construct a sequence (T (k) ) N k=0 of positive semi-definite operators interpolating between T (0) = T 1 and T (N ) = T 2 . We define recursively for every k ∈ {1, . . . , N } Before continuing, we want to verify that this definition indeed yields T (N ) = T 2 . We first remark that, because of the orthogonality relation (A.11), the equation (A.13) can be written as Thus, we get recursively for every that So, again using the representation (A.12) of T 2 , we see that Inserting the recursive definition (A.13) for T (k) , we get from the subadditivity of the supremum the estimate Since (1 − 2P k ) 2 = 1 and eigenvalues are invariant under conjugation (that is, we have λ n (AT (k−1) A −1 ) = λ n (T (k−1) ) for every invertible operator A), this simplifies to r n=1 λ n (T (k) ) ≤ r n=1 λ n (T (k−1 ).
We therefore get recursively the inequality According to Mercer's theorem, see for example [6, Chapter III, §5 and §9], we therefore get as in (A.9) that F (x, y + sv) , j ∈ N, j ≤ k.
Then, there exist constants A > 0 and a > 0 such that for every n ∈ N λ n (T ) ≤ A min where we take the minimum over all values j in the set Proof: For some δ ∈ (0, 1), we partition the domain U in pairwise disjoint open sets Q , = 1, . . . , N , with diameter not greater than δ such that N =1Q ⊃ U . We remark that there exists a constant a > 0 so that we can find for every δ ∈ (0, 1) such a partition with N sets where According to Lemma A.2, we can now get an upper bound for the behaviour of the lower eigenvalues of T by considering the eigenvalues of the integral operatorT : L 2 (U ) → L 2 (U ) with kernel F N =1 χ Q ×Q or, equivalently, the eigenvalues of the integral operators T : L 2 (Q ) → L 2 (Q ) with the integral kernels F χ Q ×Q .
To obtain an estimate for the eigenvalues of the operators T , we consider instead of T the finite rank operator which we get by approximating the kernel F on Q × Q by a polynomial and then apply Lemma A.1, see [22, §2].
So, we pick in every set Q an arbitrary point z and expand F on Q × Q in a Taylor polynomial of degree j − 1 for some j ≤ k with respect to the second variable around the points z . Then, we get F (x, y) = F j, (x, y) + C j, (x, y), x, y ∈ Q , with the Taylor polynomial F j, explicitly given by and with the remainder term C j, , which can be uniformly estimated by Since F j, is not necessarily Hermitian, we symmetrise it by defining the kernelF j, on Q × Q as Then,F j, is of the formF j, (x, y) = rj i=1 a i, (x)b i, (y) for some functions a i, , b i, ∈ C(Q ) with r j given by two times the number of elements in the set {α ∈ N m 0 | |α| ≤ j − 1}, which is r j = 2 j+m−1 m . Thus, the integral operatorT j, : L 2 (Q ) → L 2 (Q ) with kernelF j, has a finite rank which is not grater than r j . Moreover, we get from (A.18) the uniform estimate For fixed n ∈ N and j ∈ N, we now want to choose the parameter N ∈ N in such a way that N r j < n and that we can make the parameter δ as small as possible. We pick δ = a m r j n , where we assume that j is chosen such that r j < n a m , so that δ < 1 is fulfilled (an upper bound on δ is needed for an estimate of the form (A.17)). Because of r j = 2 j+m−1 m ≤ 2(j + m − 1) m , this condition on r j can be ensured by imposing Evaluating (A.22) at these parameters, we find that Simplifying the expression by estimating r j ≤ 2(j + m − 1) m ≤ 2(j + m) m and 2r j + 1 ≤ 4(j + m) m , we finally get that where we can choose j ∈ {1, . . . , k} arbitrary as long as the condition (A.23) is fulfilled.
In particular, Proposition A.3 includes the trivial case where the kernel F is a polynomial of degree K, in which case M K+1 = 0 and we obtain λ n (T ) = 0 for all n ≥ 2(a(K + m + 1)) m , since then K + 1 ∈ J K+1,n . Moreover, we find that for a general F ∈ C k (Ū ×Ū ), we may always pick j = k for n ≥ 2(a(k + m)) m to obtain that the eigenvalues decay at least as λ n (T ) ≤ Cn − k m for some constant C > 0.
For smooth kernels F , the optimal choice of j depends on the behaviour of the supremum M j of the directional derivative as a function of j. To simplify this, we estimate j µj ≤ (j + m) µ(j+m) and obtain withj = j + m λ n (T ) ≤Ãn miñ Then, by solving the optimality condition 0 = f n (ζ) = ((1 + µ) log ζ + 1 + µ + log α n )f n (ζ), we find that f n attains its minimum at ζ n = e −1 α  Since we only need the asymptotic behaviour for n → ∞, let us pick a value n 0 ∈ N such that ζ n > 1, that is α n < e −(1+µ) < 1, and ζ n < 1 a m n 2 = b α n , that is α µ 1+µ n < be, for all n ≥ n 0 . This can be always achieved since α n → 0 as n → ∞. Now, the minimum in (A.26) is restricted to the set of natural numbersj ≤ 1 a m n 2 so that we cannot simply insert forj the minimum point ζ n of the function f n . Instead, we estimate the minimum from above by the value of f n at the largest integer ζ n below ζ n : λ n (T ) ≤Ãn miñ j≤ 1 a m √ n 2 f n (j) ≤Ãnf n ( ζ n ).
Thus, using the expressions (A.27) and (A.28) for α n and ζ n , we find constants C, c > 0 such that λ n (T ) ≤ Cn m √ n exp −cn

B. Estimating the Kernel of the Integrated Photoacoustic Operator
To be able to use Corollary A.4 to estimate the eigenvalues of the operatorP * κPκ , we need to find an upper bound for the derivatives of the integral kernel F κ of the operatorP * κPκ , which is given by (4.8). Since where G κ denotes the fundamental solution of the Helmholtz equation, given by (2.16), we start with the directional derivatives of the function G κ (ω, ·). Since G κ (ω, ·) is radially symmetric, this means we can write it for arbitrary ω ∈ R and x ∈ R 3 \ {0} in the form G κ (ω, x) = g κ,ω ( 1 2 |x| 2 ) with g κ,ω (ρ) = − iω 4π √ 2π for ρ > 0, ((B.1)) this problem reduces to the calculation of one dimensional derivatives of the function g κ,ω .
Thus, the directional derivatives of G κ (ω, ·) can be calculated from the derivatives of g κ,ω , which we may estimate directly. This formula also extends continuously to the case a = 0.
Proof: By our assumptions on κ, we have that there is a constant C ≥ 0 such that m κ(ω) ≥ κ 0 |ω| β − C for all ω ∈ R. Thus, By Jensen's inequality, applied to the convex function f (r) = r N j , we can estimate Then, we find with the substitution ν = ω β that where Γ(ρ) =