Weighted Radon transforms of vector fields, with applications to magnetoacoustoelectric tomography

Currently, theory of ray transforms of vector and tensor fields is well developed, but the Radon transforms of such fields have not been fully analyzed. We thus consider linearly weighted and unweighted longitudinal and transversal Radon transforms of vector fields. As usual, we use the standard Helmholtz decomposition of smooth and fast decreasing vector fields over the whole space. We show that such a decomposition produces potential and solenoidal components decreasing at infinity fast enough to guarantee the existence of the unweighted longitudinal and transversal Radon transforms of these components. It is known that reconstruction of an arbitrary vector field from only longitudinal or only transversal transforms is impossible. However, for the cases when both linearly weighted and unweighted transforms of either one of the types are known, we derive explicit inversion formulas for the full reconstruction of the field. Our interest in the inversion of such transforms stems from a certain inverse problem arising in magnetoacoustoelectric tomography (MAET). The connection between the weighted Radon transforms and MAET is exhibited in the paper. Finally, we demonstrate performance and noise sensitivity of the new inversion formulas in numerical simulations.


Introduction
In this paper we study unweighted and linearly weighted Radon transforms of vector fields. There is a significant body of work on ray transforms (that involve integration over straight lines) of vector and tensor fields [1][2][3][4][5]. In particular, exponential and attenuated ray transforms were studied in [6][7][8][9], and momentum ray transforms were investigated in [10,11]. However, when it comes to the Radon transforms of vector fields (with integration over hyperplanes), there are very few publications [12,13]; moreover, the consideration is usually restricted to unweighted transforms of potential fields with finitely supported potentials. In the present paper we consider general vector fields (i.e. not purely potential or solenoidal), and we study both unweighted and linearly weighted Radon transforms.
As in the case of ray transforms, when studying the Radon transforms one finds it convenient to use the Helmholtz decomposition. In other words, one splits a general vector field F into the potential and solenoidal parts F p and F s , and considers transversal and longitudinal Radon transforms of both F p and F s . However, even for a finitely supported field F components F p and F s are defined in the whole space R d and they are known to have only a polynomial decay at infinity. Thus, in order to analyze the Radon transforms of F p and F s one first needs to prove that such transforms do exist (i.e. integrals over hyperplanes in R d converge). This is not completely trivial. In particular, the estimate given in the foundational book [5] on ray transforms does not guarantee the convergence of the Radon transforms. Thus, first we obtain an improved estimate for the rate of decay at infinity of the potential and solenoidal parts F p and F s of a fast decaying field F. This estimate guarantees the existence of the unweighted longitudinal and transversal Radon transforms of F p and F s .
Similarly to the case of ray transforms, the longitudinal Radon transforms of a potential field vanish. The same is true for the transversal transform of a solenoidal field. Therefore, reconstructing a general vector field from only the longitudinal or only the transversal transform(s) is not possible. However, it is not unusual in practice [6] that one of the transform types (either longitudinal or transversal) cannot be measured. In order to replace missing information one may consider measuring weighted transforms of the available type. For example, our interest in this problem stems from a certain measurement scheme in the magnetoacoustoelectric tomography (MAET). This scheme does not permit measuring a transversal transform of a certain vector field, but, in addition to longitudinal transforms one can measure linearly weighted longitudinal transforms of that field.
Below we present explicit formulas for solving two distinct problems. The first problem is that of reconstructing a general vector field from known values of its transversal transform, and from d − 1 weighted transversal transforms with various linear weights. The second problem (motivated by MAET) is the reconstruction of a general vector field from d − 1 of its longitudinal transforms and one weighted longitudinal transform (again, with a linear weight). The reader may want to compare our solutions of these problems to the results of [11], where a full vector field is reconstructed from a ray transform and a first-moment ray transform.
The rest of the paper is organized as follows. We define all the needed transforms in Section 2.1 below, and we present explicit solutions to the above two problems in Section 2.2. In Sections 3 and 4 we provide proofs of the theorems formulated in Section 2. Section 5 exhibits a potential application of the Radon transforms of vector fields to a problem arising in MAET. We further validate our theoretical results by numerical simulations, see Section 6. Finally, the proof of Theorem 1 (on the rates of decay of F p and F s ) is relegated into the Appendix.

Definitions and technical estimates.
Consider a continuous function f (x) defined in R d , subject to the condition f (x) = O |x| −d at infinity. Define a hyperplane Π(ω, p) by the equation ω · x = p, where S d−1 is the unit sphere in R d , and (ω, p) ∈ S d−1 × R. The Radon transform Rf is defined as the set of integrals of f over all the hyperplanes: where dA Π (x) is the standard area element on Π(ω, p). Properties of the Radon transform are traditionally studied for functions f (x) from the Schwartz class S(R d ). We recall that this class consists of all C ∞ (R d ) functions f (x) whose derivatives decay at infinity faster than any rational function: where α and β are multiindeces, α = (α 1 , ..., α d ), β = (β 1 , ..., β d ), α j 's and β i 's are non-negative A function f (x) ∈ S(R d ) can be reconstructed from its projections g = Rf using the well known filtered backprojection inversion formula [14]: where R # is the dual Radon transform that acts on a function g(ω, p) defined on S d−1 ×R according to the formula and where the Riesz potential I α f of a function f is expressed through the direct and inverse Fourier transforms F and F −1 as follows Let us consider now a continuous vector field F (x) = (F 1 (x), ..., F d (x)) defined on R d , d ≥ 2, whose components decay fast enough for the existence of integrals over each hyperplane (e.g., |F (x)| = O |x| −d ). Below we define several types of Radon transforms of such a field.
The componentwise Radon transform RF of F is defined in the obvious way: The transversal Radon transform D ⊥ F is the Radon transform of the projection of F onto the normal ω to the plane Π(ω, p): For each fixed direction ω ∈ S d−1 , let us arbitrarily extend ω to an orthonormal basis B = (ω, To simplify the notation, below we will suppress the dependence of ω j 's on ω. Define the longitudinal Radon transforms D k F of F, k = 1, ...d − 1, as follows: For a faster decaying vector field F (x) (e.g. satisfying |F (x)| = O |x| −d−1 ), one can define the weighted transversal transforms W ⊥ k and longitudinal transforms W k with linear weights ω k · x, k = 1, ...d − 1, by the following expressions: The present definitions of the unweighted longitudinal and transversal Radon transforms coincide with those given in [4,6] (where they are mentioned under the names of "probe" and "normal" transforms, respectively). Our definitions of the weighted transforms appear to be new; they naturally extend the notion of "moments ray transforms" [6,11] to the case of Radon transforms.
It is well known that the Radon transform of a scalar function considered on S d−1 ×R is redundant. Indeed, since Π(ω, p) = Π(−ω, −p), one concludes that [Rf ] (ω, p) = [Rf ] (−ω, −p). Similarly, by inspecting equation (3) one can see that D ⊥ F (ω, p) = − D ⊥ F (−ω, −p), where the change of sign occurs due to the factor ω· under the integral. The definitions of transforms D k , W ⊥ k , and W k depend on two vectors, ω and ω k . In general, our definition of basis B permits a significant freedom in choosing the dependence ω k = ω k (ω). However, if we restrict consideration to the case ω k (ω) = −ω k (−ω), the following redundancies will arise Such redundancies can be exploited in practice, to reduce the number of required measurements and to halve the number of floating point operations when implementing inversion formulas, both known and the ones presented below. (For example, operator R # in (2) can be computed by integration over a half of a sphere.) However, since the focus of this paper is mostly theoretical, for simplicity of presentation we will work with projections defined on S d−1 × R.
For the future reference we note the obvious relations Let us now consider a smooth and fast decaying vector field F (x) such that each component F m (x) of F (x) is a function from the Schwartz space S(R d ). We define the potential ϕ as the convolution of the divergence Φ of F with the fundamental solution G of the Laplace equation in R d : where explicit expressions for G(x) are well known: Now the potential part F p of the field F is the gradient of ϕ: and the solenoidal part F s is just the difference The following theorem is a technical result that is an important tool in our investigation.
Theorem 1. Suppose that each component F k (x), k = 1, ..., d of a vector field F (x) is a function from the Schwartz class S(R d ). Then potential ϕ and fields F p and F s given by equations (8)- (10) have the following decay rates at infinity The estimates (11)- (14) are a refinement of the well known estimate on the rate of decay of F p and F s given by Theorem 2.6.2 of [5]: with the similar bound on F p . The importance of estimates (11)- (14) for the present work is in that they guarantee existence of the transversal, longitudinal, and component-wise Radon transforms of F p and F s , so that and W ⊥ k F s cannot be defined, in general. Indeed, according to definitions (5) and (6), such transforms would require integration of fields F s and F p multiplied by linear functions in x, over hyperplanes in R d . Such products decay at infinity at the rate O(|x| 1−d ). Such decay is not sufficient for the existence of the integrals.
2.2. Main theorems. The main results of this paper are the following two theorems: If an infinitely differentiable vector field F (x) = (F 1 (x), ..., F d (x)) satisfies decay conditions (1), its divergence Φ can be reconstructed from the transversal transform D ⊥ F by applying the inversion formula (2) as follows Further, the componentwise Radon transform of F can be reconstructed from D ⊥ F and weighted transversal transforms with linear weights W ⊥ k F , k = 1, .., d − 1, as follows: Finally, field F can be recovered by inverting RF componentwise: where vectors e 1, e 2 , ..., e d form the canonical orthonormal basis in R d , and where R −1 is understood as the filtration/backprojection formula (2).
In order to formulate the next theorem, let us denote by Ψ the componentwise Laplacian Ψ of the solenoidal part of the field F s : Theorem 3. If an infinitely differentiable vector field F (x) = (F 1 (x), ..., F d (x)) satisfies decay conditions (1), the componentwise Laplacian Ψ of its solenoidal part F s and the Radon transform of Ψ can be reconstructed from longitudinal transforms D j F , j = 1, ..., d − 1, using the following formulas: Further, the divergence Φ of the field can be reconstructed from the linearly weighted longitudinal transform W 1 F and previously found Ψ as follows: where R −1 is understood as the filtration/backprojection formula (2). Finally, filed F is reconstructed from Φ and Ψ by convolving these functions with G and its gradient: We provide the proofs of theorems (2) and (3) in Sections 3 and 4, respectively. The proof of theorem (1) can be found in the Appendix.
3. Properties of the transversal transforms and proof of Theorem 2 3.1. Reconstructing the potential part of the field. Most of the material reviewed in the present section 3.1 is known. However, to make the presentation self-contained, we provide elementary proofs.
Then the transversal Radon transform D ⊥ F s of F s vanishes: Proof. Fix an arbitrary pair (ω, p) ∈ S d−1 × R and the corresponding hyperplane Π(ω, p). Consider a sphere S(0, R) of radius R centered at the origin. Further, consider the region Υ(R, p) bounded by a part of S(0, R) and Π(ω, p), and such that the interior normal to the boundary of Υ(R, p) on Π(ω, p) coincides with ω. Let us denote by ∂Υ 1 (R, p) the spherical part of the boundary Υ(R, p), i.e. ∂Υ 1 (R, p) ≡ ∂Υ(R, p) ∩ S(0, R). Since div F s = 0, the following integrals are equal where dS(x) is the standard area element on S(0, R) and n(x) is the exterior normal to the sphere. Now, let us take the limit R → ∞. Due to the fast decrease of F s (x) at infinity, the right hand side in the above equation converges to 0. The left hand side converges to Π(ω,p) F (x) · ω dA Π (x), proving that this integral is equal to 0. Since this is true for arbitrary (ω, p), equation (23) follows.
The corollary below follows immediately from Proposition 4 and Theorem 1.
is a C ∞ vector field defined on R d and decaying at infinity at rates given by equation (1), and F p + F s are defined by equations (8)- (10). Then Suppose h is a Radon integrable function with a Radon integrable derivative ∂ ∂x k h. Then the following relation holds [15]: This leads to the following Lemma. Then Proof. The divergence div H(x) has the rate of decay O |x| −d , justifying the following: where equation (7) is used on the second line of equalities.
In particular for a C ∞ field F satisfying the rates of decay (1) we obtain Since Φ is a function from the Schwartz class S(R d ) it can be reconstructed from projections using the filtered backprojection formula (2), which yields equation (17). The potential part of the field F p (x) can now be computed by combining (8) and (9):

3.2.
Reconstructing the whole filed. Due to Proposition 4, the solenoidal part F s of the field F lies in the null space of the transversal Radon transform D ⊥ , and therefore, cannot be reconstructed from the knowledge of D ⊥ F. Thus, in addition to D ⊥ F , in this section we assume the knowledge of the transversal weighted transforms W ⊥ k F, k = 1, ..., d − 1, defined by (5). This information will allow us to reconstruct the whole field F and thus to complete the proof of theorem 2 First, for the future use we would like to find projections of RF on the vectors of the basis B. By combining equations (7) and (24) one observes: Let us find projections of RF on vectors ω 1 , ..., ω d−1 of the basis B. Note that, due to (4) We start with ω 1 · RF p : for any (ω, p) ∈ S d−1 × R. Since the numbering of vectors ω 1 , ..., ω d−1 in the basis B is arbitrary, we conclude that In other words, a longitudinal transform of a potential field vanishes. Since basis B is orthonormal, by combining (28) with (27) one obtains the following formula: Thus, the componentwise Radon transform of the potential part F p of the field F can be easily recovered from the transversal transform ωD ⊥ F .
Let us find what information can be extracted from the weighted transversal transforms W ⊥ k F. It follows from the definition (5) Due to (28) term ω k · RF p vanishes, and one obtains These equations combined with (23) determine projections of vector-valued function RF s onto the vectors of the orthonormal basis B, leading to the following result: By combining the latter formula with equation (29) we arrive at the formula (18) that gives an explicit expression for [RF ](ω, p). Since field components F j (x) are functions from the Schwartz space, formula (2) can be used to reconstruct F j 's from components of the vector-valued [RF ](ω, p), thus yielding equation (19). The proof of Theorem 2 is complete.

Properties of longitudinal transforms and proof of Theorem 3
In this section we assume that only longitudinal transforms D j F, j = 1, ..., d − 1, and one of the weighted longitudinal transforms (e.g., W 1 F ) are known. Our goal is to reconstruct field F from these data.

4.1.
Reconstructing the solenoidal part of the field. Proposition 7. Suppose F is a smooth vector field satisfying the decay conditions (1), and F p , F s are its potential and solenoidal parts, respectively. Then longitudinal transforms D j F p of F p vanish, j = 1, ..., d − 1, and the Radon transform of the solenoidal part can be expressed through D j F as follows: Proof. Using (7) and (29) one obtains which implies that all longitudinal transforms of the potential part of a field vanish: Further, by expanding vector RF in basis B and using (7) again see that so that (30) holds.
Equation (30) shows that the longitudinal transforms D j F, j = 1, 2, ..., d − 1, contain enough information to obtain the componentwise Radon transform of the solenoidal part of the field. However, a straightforward componentwise application of the inversion formula (2) is not justified in general, since components of the field F p are not in the Schwartz space. It is known that formula (2) remains valid for slower decaying functions (see Chapter 1 of [15]). However, reconstruction of functions decaying at the rate (12) still, in general, cannot be guaranteed. While we conjecture that inversion formula (2) can be used for componentwise inversion of (30), we will not prove this statement here. Instead, we notice that by computing the second derivative of equation (30) in p one obtains the Radon transform of the componentwise Laplacian Ψ of F s : Let us find out the rate of decay of components of Ψ at infinity. Using (10) one obtains Since each component of field F belongs to the Schwartz space, so does Ψ k (x), k = 1, ..., d. Therefore, equation (31) can be inverted componentwise using formula (2), thus proving formula (20). Knowing Ψ, the solenoidal part F s of the field can be recovered as the following convolution:

4.2.
Reconstructing the whole field. In this section we will show that, assuming that Ψ is known (for example, reconstructed using formula (20)), the divergence Φ of the field can be reconstructed from the weighted longitudinal transform W 1 (F ) using formula (21), and the whole field F can be obtained as convolutions (22). As before, we will try to differentiate the weighted transform W 1 F. More precisely, let us evaluate the following expression: The second term in the right hand side of (32) can be transformed as follows: The first term in the right hand side of (33) can be seen to be equal to (e k · ω 1 ) D 1 (F p ); it vanishes as a longitudinal transform of a potential field. The remaining second term in (33) can be simplified further: where integration by parts was performed with respect to y 1 . By combining (32), (33), and (34) we thus obtain Now, let us apply the operator (e k · ω) ∂ ∂p again, this time to equation (35): By summing the above formula in k from 1 to d one obtains This allows one to simplify the first term in the right hand side of (36) as follows: which holds since the longitudinal transform of a potential field D 1 (∇(ω 1 · F s (x)) vanishes. By combining (36) and (37) we arrive at the following formula (38) R(Φ) = R ((x · ω 1 )(ω 1 · Ψ(x))) − ∂ 2 ∂p 2 W 1 (F ) Now the Laplacian Φ of the potential ϕ can be reconstructed by inverting the Radon transform in (38), yielding the whole field can be reconstructed by computing convolutions (22). This completes the proof of theorem 3.

Vector fields in magnetoacoustoelectric tomography
Our interest in the Radon transforms of vector field is motivated, in part, by an inverse problem arising in magnetoacoustoelectric tomography (MAET). This imaging modality is a novel coupledphysics technique designed to image the electrical conductivity of biological objects. It is based on measurements of electric potential arising in conductive tissues when they move in a magnetic field. In detail, one places the object of interest in a strong constant magnetic field and illuminates it with ultrasound pulses [16][17][18][19][20]. Frequently this is done with the object immersed in conductive saline, which provides good acoustic coupling and facilitates the measurements of the arising electric potential with the use of electrodes immersed in the liquid. The said potential results from the interaction of the vibrational motion of electrons and ions contained in a conductive tissue, with magnetic field. This generates the Lorentz forces that separate the particles of opposite polarities and, in turn, results in Ohmic current flowing through the object and the saline. The electric potential associated with this current is then measured outside of the object, providing the data for the future MAET reconstruction.

5.1.
A traditional data acquisition scheme. In the remaining part of the Section 5 and in Section 6 we work with the three-dimensional space.
It has been shown ( [20]) that when the tissue with conductivity σ(x) moves with velocity V(t, x) within magnetic field B, the arising Lorentz force will generate Lorentz currents J L (t, x) given by the formula The vibrational velocity V(t, x) of the tissues arising due to the ultrasound excitation is governed by the standard wave equation with the speed of sound that can be assumed constant within soft tissues. Without loss of generality the speed of sound can be set to 1. Then V(t, x) and the acoustic pressure p(t, x) can be related to the velocity potential ζ(t, x) by equations Here the density ρ is assumed to be constant within soft tissues and equal to the density of water. The scalar velocity potential ζ(t, x) itself also satisfies the wave equation in the whole space R 3 : . The time scales of this model are such that the electromagnetic effects are much faster than the mechanic motion of the liquid [19]. Therefore, the currents in the system can be considered stationary, corresponding to velocity V(t, x) at the given time t. Then, it can be shown that the difference of potentials M (t) measured by a pair of electrodes can be expressed as follows [21] (40) where the lead current I(x) is the current that would flow through the object in the absence of the magnetic and acoustic excitation, if a unit potential difference were applied to the electrode pair. This quantity appears in (40) because I(x) also describes the sensitivity of the measuring system to a dipole placed at the point x. Finally, the domain Ω in the above equation is the volume occupied by the saline and by the object immersed in it. Below, it will be convenient for us to consider a model where Ω is large and can modeled by the whole space R 3 . A measurement corresponding to a given acoustic wave ζ(t, x) is, according to (40), a function of one variable. The goal of MAET is, by using a sufficiently rich set of excitations ζ(t, x), to collect enough information for reconstruction of the conductivity σ(x) of the tissues.
In the early mathematical work on MAET [21,22] mathematicians would assume that the object and the electrodes remain fixed and the transducer is moved around the object providing a large family of excitations ζ(t, x). Then the inverse problem of MAET naturally decouples into two steps. Since curl C(x) is independent from ζ(t, x), one considers (40) as values of projections of the quantity B · C(x) on the complete set of excitations ζ(t, x), and reconstructs B · C(x). Then, the second step is to reconstruct the conductivity σ(x) from B · C(x), possibly from measurements repeated with two or three different orientations of B. Depending on the waveforms ζ(t, x), the  Figure 1. The novel MAET scheme; the electrode/transducer assembly rotates around the object first step frequently can be reduced to one of the known tomography problems. For example, if one illuminates the object by ideal plane waves with various directions ω ∈ S 2 , the resulting measurements can be expressed the Radon transform of B · C(x), that can be easily inverted. Similarly, if one assumes an ideal point-like transducer that produces spherical outgoing waves, the problem reduces to the inverse source problem of thermoand photoacoustic tomography, whose solution is well known by now (see, e.g. [23,24]).
The measuring scheme described above is easy to analyze. However, it does not work well in practice. Indeed, if electrodes and the object are held in a fixed position, there are very few directions from which the transducer can send sound waves into the object without illuminating the electrodes, which generates strong spurious elecrtic pulses that overwhelm the usefull signal. Thus, researchers are investigating a different approach to data acquisition [25,26], which assumes that object is rotated while the electrodes are kept stationary. Equivalently, one can keep the object fixed, and rotate the electrodes and transducer(s). In both cases, the curl C(x) becomes a function of the object (or electrodes') position, and the traditional two step reconstruction procedure described above is not applicable anymore.

5.2.
Acquisition scheme with a rotated object. We thus consider here the novel acquisition scheme for MAET, with a rotating electrode/transducer assembly, as shown in Figure 1. The object under investigation is immersed in a conductive saline, and the assembly rotates around it. For simplicity, we model the propagation of currents in this scheme assuming that the electrodes are placed far away from the object. Here, the conductive medium is presumed to occupy all of R 3 , with the conductivity σ(x) being constant and known outside of the support Ω 0 of the inhomogeneity, i.e. σ(x) = σ 0 for x ∈ R 3 \Ω 0 . Then the lead current I is a function of x and the orientation ν of the electrodes, i.e. I ≡ I ν (x). We assume that, in the absence of the inhomogeneity, the electrodes generate field E 0 ν = ν. In the presence of inhomogeneity, additional potential u ν (x) will arise, so that the current can be expressed as subject to the following condition at infinity: Due to the absence of sinks and sources of charges in the medium, current I ν (x) is solenoidal. By setting to zero the divergence of (41) we find that potential u ν (x) solves the divergence equation, subject to the decay at infinity The above simplified model will allow us to express potential u ν and current I ν for an arbitrary orientation ν through three "basis" solutions. Indeed, let us consider the solutions u (j) (x), j = 1, 2, 3 of (42), (43) corresponding to directions ν = e 1 , e 2 , e 3 , where e j 's are the canonical vectors in R 3 : The corresponding total currents and their curls will be denoted by I (j) (x) and C (j) (x) respectively: Due to the linearity of the problem (42), (43) with respect to the right hand side of (42), for an arbitrary direction ν the potential u ν (x) and current I ν (x) can be represented as the following linear combinations: Let us denote by C ν (x) the curl of the three-dimensional field I ν (x): Recall that MAET measurements are directly related to B · C ν (x) (see equation (40)). Let us assume for now that the transducer is oriented along the vector ω perpendicular to ν, and is producing ideal plane waves. Then, the corresponding measurements M (ν, ω, t) can be expressed as Here the integration is restricted to Ω 0 since the curls C (j) (x) of currents I (j) (x) vanish within any region with constant conductivity, i.e. outside of Ω 0 . By combining equations (44) and (45) one obtains where we introduced the vector field C(x) defined as follows We thus recognize M (ν, ω, t) as a longitudinal Radon transform of the vector field C(x). If one directs vector ν to be parallel to one of the vectors ω 1 or ω 2 orthogonal to ω, measurements M (ν, ω, t) coincide with the longitudinal transforms D k C (ω, p) defined by equation (4): If one manages to reconstruct from MAET measurements field C(x), projections of curls B · C (j) are easily found: B · C (j) (x) = ρe j · C(x), j = 1, 2, 3.
Then, the measurements can be repeated with alternatively directed B, until C (j) 's can be determined. After that, currents I (j) and conductivity σ(x) can be reconstructed, following the techniques presented in [21,22]. In a simplified two-dimensional setting (as in [25]), curls C (j) , j = 1, 2, are oriented orthogonally to the plane in which currents are flowing, and magnetic induction B is parallel to C (j) 's. Additional directions of B are not needed in this case. However, analysis presented in the previous sections of this paper shows that only a solenoidal part of a vector field C(x) can be reconstructed from known longitudinal transforms D k C, k = 1, 2. In general, there is no reason to expect that field C(x) is solenoidal. As a way to remedy this situation, we propose to conduct additional measurements, by illuminating the object with linearly modulated acoustic waves in the form with directions ω varying over S 2 , and ω 1 aligned with the electrode directions. Such measurements N (ω 1 , ω, t) are described by the formula they can be expressed as the weighted longitudinal transform W 1 C: Theorem 3 states that the vector field C(x) can be reconstructed from its longitudinal transforms D 1 C and D 2 C and weighted longitudinal transform W 1 C using formulas (20)- (22). MAET measurements using linearly modulated waves (46) have not been implemented previously, in part because the benefit of such measurements have not been previously discussed in the literature. However, there is no physical obstacles for conducting such an experiment. Indeed, functions in the form (46) are easily seen to satisfy the wave equation. They can be generated in a number of ways. For example, if a transducer array is used for sound generation (as in [26]), such waves can be obtained by scaling linearly the excitation voltage along the transducer elements. If a synthetic flat detector is utilized (as in [25]), one obtains the desired result by a weighted averaging of individual measurements. Such sound waves can also be excited using optically generated ultrasound [27,28], by using optical excitation with linearly varying intensity.
We will not attempt to simulate a full MAET experiment with linearly modulated sound waves in this paper, leaving it to the future work. Below we present numerical simulations of reconstruction of a 3D vector field from its longitudinal transforms D 1 F , D 2 F and the weighted longitudinal transform W 1 F .

Numerical simulations
The goal of this section is to demonstrate the validity of the exact reconstruction formulas (20)- (22) in a numerical experiment. To this end we picked a smooth phantom F (x) defined in the unit ball B(1, 0) in R 3 . Each component F j (x), j = 1, 2, 3 is a linear combination of a rather arbitrary collection of shifted radially symmetric functions ("bumps") k,j were chosen to lie in one of the planes (B(1, 0)) function. The phantom is shown in Figure 2, and the values of constants x k,j , R k,j , and a k,j used in our simulations can be found in Table 1. In addition, M 1 = 5, M 2 = 5, M 3 = 8.
Formulas (5)- (7) show that one can find values of longitudinal transforms D 1 F , D 2 F , and W 1 F by computing the standard and the linearly weighted Radon transforms of each of the F j , j = 1, 2, 3.
The latter transforms of the radial bump functions f x − x (c) k,j , R k,j are given by the following k,j , R k,j (ω, p). While our formulas are valid for any choice of orthonormal basis vectors ω 1 (ω) and ω 2 (ω), for numerical simulations we defined these vectors as follows. Vector ω 2 was chosen to lie in the horizontal plane spanned by canonical vectors e 1 and e 2 ; it was computed as follows: where ω * 2 (ω) = (−(ω · e 2 ), (ω · e 1 ), 0).
The following grid in the variables (ω, p) was used to compute the Radon transforms. Variable p was discretized using a uniform gird with 257 nodes in the interval [−1, 1]. Vector ω(θ, ϕ) = (sin ϕ cos θ, sin ϕ sin θ, cos ϕ) was discretized using a product grid on [0, 2π] × [0, π], with 513 uniformly spaced nodes in the variable θ and 256 Gaussian nodes in the variable t = cos ϕ. For simplicity of presentation we did not use the redundancy in the Radon transform to reduce the required data and the computational complexity. However, in practice it is sufficient to vary ω over half a sphere and multiply the result by the factor of 2.
The inversion of the classical Radon transform required by equations (20) and (21) was implemented by discretizing the 3D version of the formula (2), with α = 0: The inversion was computed in the nodes of 257 × 257 × 257 Cartesian grid in x, for |x| ≤ 1 only. For our first simulation, the derivatives in p in (21) and in (47) were computed by a spectrally accurate algorithm, using the Fast Fourier transform (FFT), in order to achieve high accuracy when processing theoretically exact data. The components of the reconstructed fields F s (x) and F p (x) are shown in Figure 3 (the gray scale used in the images is the same as in Figure 2). When added together, these fields produce an accurate approximation to the exact F (x). When plotted in a grey scale figure (not shown here) the reconstructed F (x) is indistinguishable from the exact field presented in Figure 2. In this case, the relative L 2 error of the reconstruction is 0.09% and the relative L ∞ error does not exceed 0.3%. This is consistent with the exactness of our reconstruction formulas.
Our second numerical simulation aims to demonstrate the noise sensitivity of formulas (20)-(22). In the above mentioned equations, functions Φ and Ψ j , j = 1, 2, ..., d are reconstructed from the second derivatives of the data in p. This is followed by convolutions with smoothing kernels in (22). However, the solenoidal part F s of the field is obtained by convolutions with the fundamental solution (Ψ j * G), j = 1, ..., d, whereas the potential part is computed by convolution of Φ with the gradient ∇G of the fundamental solution. This additional differentiation implies that the potential part should be more sensitive to high spatial frequencies of the noise.
In order to test this conclusion we added to the data D 1 F , D 2 F , and W 1 F a small normally distributed spatially uncorrelated noise with relative intensity 0.1% in L 2 norm. Spectral differentiation in p in (21) and in (47) was replaced by the standard second order symmetric finite difference formula. This has a mild regularizing effect compared with the spectral differentiation. The fields F s (x) and F p (x) reconstructed from the noisy data are shown in Figure 4 (the gray scale used in this figure is the same as in Figure 2). Comparison with the Figure 3 shows that the solenoidal part F s (x) is little affected by this mild noise, while reconstructed F p (x) contains much stronger high frequency artifacts (the reader may want to magnify the figure to see this clearly). Indeed, a quantitative comparison reveals that the relative error in F s (x) is 1.1% in L 2 norm and 1.3% in L ∞ norm. On the other hand, the relative error in F p (x) is 63% in L 2 norm and 74% in L ∞ norm.
The total reconstructed field is the sum of F s (x) and F p (x). It is depicted in Figure 5 (the gray scale is the same as in Figure 2). Due to the high level of artifacts in F p (x), the total field also contains significant error, with the relative error equal to 36% in L 2 norm and 41% in L ∞ norm. It should be noted that the high error in F p (x) is a manifestation of the poor conditioning of the problem of reconstructing the potential part of the field from a linearly weighted longitudinal transform W 1 F . Indeed, formula (30) shows that the Radon transform of F s is expressed as a linear combination of data D 2 F . Thus, the conditioning of finding F s is similar to conditioning of inverting the standard scalar Radon transform. On the other hand, in the equation (35) the Radon transform R F p k is expressed through the derivative of the data ∂ ∂p W 1 (F ). This additional differentiation of data makes the problem of reconstructing F p significantly more ill-posed than that of inverting the regular Radon transform. This leads to the appearance of strong high frequency artifacts in the reconstructed F p .  Figure 5. Field F reconstructed from noisy data In order to convince the reader that this is indeed a high-frequency phenomenon, we applied a low-pass linear filter to the total reconstructed field F (x), obtaining a smoothed field F smooth (x). In detail, each component F smooth k (x) of F smooth (x) was obtained by applying filter η(ξ) in the Fourier domain: ](x), k = 1, 2, 3. where F and F −1 are the forward and inverse Fourier transforms, and filter η(ξ) was given by the formula η(ξ) = 0.5 1 + cos π|ξ| 0.4f Nyquist for |ξ| < 0.4f Nyquist , 0 otherwise, where f Nyquist is the Nyquist frequency of the spatial discretization in x. The relative errors in the so found approximation F smooth (x) where 12% in L 2 norm and 19% in L ∞ norm.
We would like to stress that the reconstruction algorithm presented here, based on direct discretization of our inversion formulas, is meant only to illustrate the exactness of these formulas (when applied to accurate data), and to demonstrate the increased sensitivity of these formulas to noise (in comparison to the standard Radon inversion). The development of a more practical, efficient and robust algorithm is a matter of the future work. Such an algorithm would require a prudent choice of a regularization technique, to reduce the noise sensitivity. An optimal choice of such technique depends heavily on the parameters of a particular application, such as the signal-tonoise ratio, spectral content of the noise, desired resolution, etc. For a general overview of classical regularization methods we refer the reader to the book [29] and article [30]. The regularization methods used recently in vector tomography include the singular value decomposition [31], the method of approximate inverse [32], and an expansion in a series of orthogonal polynomials [12]. These topics, however, are outside of the scope of the present paper.
where C d is the volume of the unit ball in R d . Obviously, for any y ∈ B(R), |y| ≤ R. Since |x| = 2R, |x − y| ≥ R, and |g(x − y)| ≤ C g (1 + |x − y|) M ≤ C g (1 + R) M = 2 M C g (2 + |x|) M . Therefore, I B can be bounded as follows On the other hand, for y ∈ R d \B(R), |f (y)| can be bounded by C f (1+R) K so that the following inequality holds Finally, by combining inequalities (50) and (51) one proves Lemma 8.
We are ready to prove Theorem 1.
Proof. First, we establish the rate of decay at infinity of the potential ϕ(x) given by the convolution (8). We note that divergence Φ(x) belongs to the Schwartz space S(R d ) and, therefore, for any l ≥ 0, there is a constant C l such that |Φ(x)| ≤ C l /(1 + |x| l ). On the other hand, the derivatives of the fundamental solution G(x) decay as follows: Let us introduce an infinitely smooth nonnegative cut-off function η(t), t ∈ R, with η(t) = 1 for every t ∈ (−1/2, 1/2) and η(t) = 0 for |t| ≥ 1. Convolution (8) can be re-written as ϕ(x) = I 1 (x) + I 2 (x), Then |I 1 (x)| is bounded by C G C l /(1 + (|x| − 1) l ) for any l ≥ 0. The second term is the following convolution The latter sum is the sum of convolutions of functions satisfying conditions of Lemma 8, where the role of f is played by F j with M ≥ 2d − 1 (since F j 's are Schwartz functions), and the role of g is played by ∂ ∂x j [G(x)(1 − η(|x|))] with K = d − 1. Therefore, I 2 (x) has the desired rate of decay O |x| 1−d . This term dominates the sum I 1 (x) + I 2 (x) at infinity. This proves equation (11).
The estimate for the derivatives of ϕ(x) can be obtained in a similar way. Indeed where where C G still given by (53). Since second derivatives of F are Schwartz functions, |I 3 (x)| decays faster than any power of 1/|x|. For the term I 4 (x) we observe: The rate of decay of derivatives ∂ 2 ∂x k ∂x j [G(x − y)(1 − η(|x − y|))] coincides with the decay rate of ∂ 2 ∂x k ∂x j G(x); it is given by (52). Now, the application of Lemma 8 establishes that |I 4 (x)| = O |x| −d . This proves (12) for F p . The similar estimate for F s comes from F s (x) = F (x) − F p (x), where the second term dominates at infinity.