{'fraction_non_alphanumeric': 0.043677220859134634, 'fraction_numerical': 0.013163531565214005, 'mean_word_length': 3.786559791220058, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 37, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coherent quantum tomography

We discuss a quantum mechanical indirect measurement method to recover a position dependent Hamilton matrix from time evolution of coherent quantum mechanical states through an object. A mathematical formulation of this inverse problem leads to weighted X-ray transforms where the weight is a matrix. We show that such X-ray transforms are injective with very rough weights. Consequently, we can solve our quantum mechanical inverse problem in several settings, but many physically relevant problems we pose also remain open. We discuss the physical background of the proposed imaging method in detail. We give a rigorous mathematical treatment of a neutrino tomography method that has been previously described in the physical literature.


Introduction
We present a new imaging modality which we call coherent quantum tomography. The reason for this name will be discussed in section 3. Our goal is to model and solve the following problem: An object Ω ⊂ R n is to be imaged nondestructively. Properties of the object are described by a (hermitean) matrix valued function H on Ω, a quantum mechanical Hamiltonian. Particles are fired through Ω and they experience time evolution due to the Hamiltonian H and the final state of the particle is measured. How much about the function H can we infer from this data?
Of course, much depends on what initial states in the state space are possible and what kinds of measurements are made. We will mainly focus on two kinds of situations, a general case with ideal data and a more realistic one with neutrino oscillations, but the discussion of section 2 is valid for any data. These problems will be described in more detail below.
The imaging modality we discuss is a quantum mechanical one and relies on coherence between states in a multiple dimensional state space. We restrict our attention to finite dimensional state spaces, but many practically interesting systems are indeed finite dimensional. The practical setting resembles that of X-ray tomography (sending rays through an object and measuring what happens to them) but the physical phenomena behind the measurement and the mathematical model are rather different. Our analysis does, however, lead to a certain type of X-ray transform.
The application we have in mind is imaging the Earth with neutrinos. The idea is to produce neutrinos somewhere at the surface (in a nuclear reactor or a particle accelerator) and measure neutrinos of different flavours coming from this source. Neutrinos come in three flavours, electron, muon and tau neutrinos (ν e , ν µ , ν τ ) and typical nuclear reactions produce electron neutrinos or electron antineutrinos. These neutrino flavours will oscillate -turn into each other while propagating -and details of this oscillation depend on the medium they traverse. The problem is then to recover properties of the medium from neutrino oscillation data. For basics of neutrino physics, we refer to the book [15] and section 3.1.
Neutrino tomography of the Earth is not a new idea. Tomography of our planet using neutrino absorption has been discussed in [46,22,45] and using matter effects in neutrino oscillations in [35,28,25]. A discussion and comparison of these two methods can be found in [66]. It has also been proposed that geoneutrinos produced in nuclear decays in the Earth could be used for imaging the Earth [8] and that neutrino tomography could be used as a means of earthquake prediction [64]. A spectrometric approach to using neutrinos in imaging the Earth was recently proposed in [47]. We aim to bridge the gap between the physical literature of the idea of neutrino tomography and the mathematical literature of inverse problems and imaging by providing a physically meaningful and mathematically rigorous discussion of the neutrino oscillation tomography problem.
Due to the quantum mechanical nature of the problem, our data is phaseless. For some other results concerning phaseless data, we refer to [38,65,21].
Problems considered in this paper can also be posed on manifolds, re-placingΩ by a compact Riemannian manifold with boundary and lines by geodesics. Euclidean geometry is most relevant for many physical applications, so for the sake of simplicity we restrict our attention to it. A reader familiar with differential geometry will observe that many arguments carry over to Riemannian manifolds. We end up introducing some new ray transform problems (see problems 5.1 and 3.1), which can also be asked on Riemannian manifolds. We provide a partial solution in the Euclidean setting in theorems 7.2 and 7.3. See section 7 for more details on ray transforms.
1.1. Outline. We will introduce the quantum mechanical problem in detail in section 2, and then we will have the language to state our main results in section 2.1. We will overview the physical aspects of the proposed indirect measurement in section 3. Section 4 is devoted to converting our phase free data to a more useful form. In section 5 we will prove some of our main results using a pseudolinearization argument. This will involve some weighted X-ray transforms where the weight is a matrix; these are introduced in section 5.1. In section 5.2 we will reduce our quantum mechanical problem to a ray transform problem. We will also discuss the linearized problem in section 5.3. Our quantum mechanical problem involves a time ordered exponential, and we investigate the same problem without time ordering in section 6. In section 7 we will show that some of our new weighted X-ray transforms are injective with very rough weights.

The quantum mechanical problem
We will use natural units: c = = 1. Consider a point particle moving in n-dimensional Euclidean space, with position at a time t given by γ(t) = x 0 + tv ∈ R n , where |v| = 1. The particle is ultrarelativistic -its speed is practically the speed of light -and it is assumed not to scatter or be absorbed.
Besides moving in space, the particle has a state in a quantum mechanical state space C N where time evolution is given by the Schrödinger equation. Vectors in the state space C N will be referred to as states.
Note that the dimension n of the ambient space and the dimension of the state space N need not have anything to do with each other. The dimension N need not be high to accommodate interesting physical systems. The case N = 1 is in a sense empty, but two state systems (such as spin) provide an example where N = 2 and for neutrino oscillations N = 3 is relevant (excluding possible sterile neutrinos).
Suppose that the Hamiltonian governing the time evolution in the quantum mechanical state space C N is a hermitean N ×N matrix H(x) depending on the position of the particle in R n . In this case the Schrödinger equation for the state Ψ(t) ∈ C N is i∂ t Ψ(t) = H(γ(t))Ψ(t).
For notational simplicity, let us drop the line γ and consider H a function of time. Let U H (t 2 , t 1 ) be the time evolution operator (solution operator) for which Ψ(t) = U H (t, t 0 )Ψ 0 is the unique solution to Formally, we may write U H as a time ordered exponential 1 1 For a discussion of the time evolution operator as a time ordered exponential, see textbooks in quantum mechanics, such as [49]. The time ordered exponential has a series representation known as the Dyson series.
If [H(t), H(s)] = 0 for all t, s ∈ [t 1 , t 2 ], then this is just the matrix exponential. We will consider unordered time evolution, where T exp is replaced by exp, in section 6.
Consider the time evolution from creation of a particle at t = 0 to observation at t = T . Suppose there is a set I ⊂ C N of possible initial states Ψ(0) and a set F ⊂ C N of observable final states. We assume that we can prepare the particle to any initial state in I at will. If the final state Ψ(T ) = U H (T, 0)Ψ(0) is measured via an operator A, we can observe (from the statistics of a large number of measurements) the norms of the projections of Ψ(T ) to eigenspaces of A. If the spectrum of A is nondegenerate, we can measure |a * Ψ(T )| 2 for all eigenvectors (eigenstates) a of A. The set F consists of all nondegenerate eigenvectors of all observables that can be used for measurements.
Consider a domain Ω ⊂ R n which we cannot enter but through which we can fire particles. Our aim is to obtain information about the Hamiltonian H in Ω from measurements of particles sent through it. This is a nondestructive measurement problem or an inverse problem.
We assumed the particle to be ultrarelativistic (|v| = 1) and also scattering and absorption to be negligible. This means essentially that the medium only affects the particle in its time evolution in the state space C N . As mentioned in section 1, the main example we have in mind is neutrinos travelling through the Earth, the objective being to recover the inner structure of the planet from neutrino oscillation measurements. Information from such measurements should, of course, be combined with similar measurements by other methods (most importantly seismic imaging), but here we focus on the quantum mechanical tomography problem. We adopt the following notation: For a line segment γ : [0, T ] →Ω we denote U γ H := U γ H (T, 0), the time evolution operator over the entire length of the line.
The problem is now this: We will assume throughout that the Hamiltonian only depends on the position. It is also possible (and physically relevant, see section 3.2) that the Hamiltonian depends on the direction v. If the dependence on v is polynomial, the X-ray tomography problems that arise below will have to be generalized to tensor tomography problems. Recovering a tensor field from its integrals over all lines (or geodesics) can only be possible up to potential tensors. This natural gauge condition is the only obstruction in many situations; see for example [53,54,40,41]. Definition 2.2. We say that the pair of subsets I , F ⊂ C N gives ideal data if for any two unitary matrices U and V the condition |a * U b| = |a * V b| for all a ∈ F and b ∈ I implies U = e iϕ V for some ϕ ∈ R.
We also assume perfect data in the sense that the measurements are done over every line through the domain. Partial data problems with restricted sets of lines are also physically relevant, but we do not pursue them here.
The first problem is inferring as much as possible about U γ H from the data. It follows from lemma 4.1 below that in the case of ideal data we can recover the matrix U γ H up to a phase factor. This implies that if we add a position dependent scalar multiple of the identity matrix, f (x)I, to H(x), the data is not changed. We are therefore completely blind to multiples of the identity in H, and the true problem is whether this is the only obstruction to recovering H.

2.1.
Results. We will show in theorem 5.2 that ideal data (in the sense of definition 2.2 above) is enough to determine the trace free part of the Hamiltonian matrix H everywhere, provided that certain weighted X-ray transforms are injective.
Theorem 5.2: Let Ω ⊂ R n , n ≥ 2, be a convex bounded smooth domain and write M =Ω. Suppose the sets I , F ⊂ C N , N ≥ 2, give ideal data. Let H : M → C N ×N be continuous and pointwise hermitean. Assume that we know |a * U γ H b| for all a ∈ F , b ∈ I and every line γ : [0, T ] →Ω through Ω. If I W is injective for any continuous weight W : SM → SU (N 2 ), then the data uniquely determines the trace free part of H. Here I W is the X-ray transform with weight W (see section 5.1 for a definition) and SM = M × S n−1 denotes the sphere bundle of M .
In theorem 5.4 we will consider the special case where H(x) = H 0 (x)+f (x)G(x), where H 0 and G are known matrix functions and f is an unknown scalar function; the result is that ideal data determines f , provided that G is nowhere a multiple of the identity.
for a function f ∈ C 1,α (M ). Assume that we know |a * U γ H b| for all a ∈ F , b ∈ I and every line γ : [0, T ] →Ω through Ω. If H 0 and G are known and G is nowhere a multiple of the identity, then the data uniquely determines the function f . We will discuss the linearized problem in section 5.3.
In section 6 we will introduce an unordered version of our main problem 2.1. We emphasize that the unordered problem 6.1 is neither a generalization nor a specialization of problem 2.1. For this problem we will only be able to solve a linearized version; the solution of the unordered problem linearized at a constant background is given in theorem 6.13. Theorem 6.13: Let Ω ⊂ R n , n ≥ 2, be a bounded domain. Let H ∈ C(Ω, C N ×N ) take values in hermitean matrices and let H 0 be a fixed hermitean matrix. Suppose that for every line γ : [0, T ] →Ω with endpoints on ∂Ω we know up to an unknown multiple of the identity matrix. Then we can recover the trace free part of H(x) for every x ∈ Ω. Finally, we will analyze the X-ray transform problem with matrix weights that arises in theorem 5.2. We will show in theorem 7.2 that in Euclidean domains of dimension three or higher a continuous C N -valued function is determined by its weighted integrals over all lines when the weight is an invertible C N ×N matrix depending on position and direction in a sufficiently continuous way. In fact, we will prove a local injectivity result for this weighted transform; this is stated separately in theorem 7.3.
Theorem 7.2: Let M be the closure of a bounded convex domain in R n , n ≥ 3, and let β > 0. Let W ∈ C β (M, Lip(S n−1 , C N ×N )) be invertible at every point on SM = M × S n−1 . (This regularity assumption is true, in particular, if W ∈ C 1,β (SM, C N ×N ).) If f ∈ C(M, C N ) satisfies I W f = 0, then f = 0. Combining theorems 5.2 and 7.2, we obtain the following theorem which is our main answer to problem 2.1.
Let Ω be a bounded convex domain in R n , n ≥ 3, and let β > 0. Suppose H,H ∈ C 1,β (Ω, C N ×N ). Suppose the sets I , F ⊂ C N , N ≥ 2, give ideal data. The two Hamiltonians H andH give the same data in the sense of problem 2.1 if and only ifH = H + ϕI for some scalar function ϕ : Ω → R.

Physical discussion
More detailed physical remarks and background in addition to section 2 are in order, and we give them in this section. A reader who is interested in solving the mathematical problem posed above instead of justifying it physically may wish to skip this section.
In our model a quantum mechanical particle has well defined position and momentum at all times. This situation is in fact forbidden by the uncertainty principle but it is an excellent approximation. In a more careful model particles would have to be modeled by wave packets which are nontrivially distributed in position and momentum. We have (implicitly) assumed that uncertainty in momentum is very small in comparison to the mean value of momentum. Therefore momentum of a particle as a single vector is a relatively well defined concept, and if momentum is large, even very small relative uncertainty allows for substantial absolute uncertainty. The necessary amount of uncertainty in position may be so small that a particle can essentially be considered point-like. Well defined position and momentum (which could also be informally described as a plane wave supported at a point) are a fairly good approximation if momentum is large and the background medium does not have features on shorter length scales than position uncertainty. When the Earth is imaged by neutrinos, these assumptions are easily met. (We will also discuss momentum dependent Hamiltonians in section 3.2.) Our mathematical model is only relevant in physical settings where the particles are ultrarelativistic and experience very little scattering or absorption. Such particles can be reasonably approximated by pointlike particles with a quantum mechanical state as modeled in section 2. This is a common leading order approximation in particle physics.
The original motivation behind this article -to which the model is not restricted -is to image the Earth by sending rays through it. Most particles experience so much absorption or scattering that practically no signal makes it through the planet. For neutrinos, however, absorption and scattering are very weak and quite difficult to measure. Intensity of a neutrino beam is not very sensitive to medium properties, so X-ray type imaging is difficult. On the other hand, a quantum mechanical phenomenon called neutrino oscillation is sensitive to electron density in the medium. This calls for quantum mechanical measurements of the type described in section 2. The fact that quantum mechanical evolution is sensitive to something but intensity is not shows that the underlying physical phenomenon is quantum mechanical interference. Such interference is only possible for a coherent superposition of quantum states, whence the term "coherent quantum tomography". Neutrino oscillation will be discussed in much more detail in section 3.1 below.
The observation that we could -in principle at least -see more things with neutrinos than other particles due to the weakness of their interactions is not a new one. The Universe became transparent to photons at the age of about 380000 years and optical observations do not allow to see a younger Universe directly. For neutrinos the corresponding age of the Universe is about one second, whence neutrinos "see" significantly older structures than photons do. The problem is that neutrinos are much, much more difficult to measure than photons; more information is available but it is in a form difficult to measure. For more about neutrino astrophysics, see [15,26,61]. This is the case for imaging the Earth with neutrinos, too: current technology is insufficient for very precise measurements but there is no fundamental physical obstruction to seeing new things with neutrinos.
The usual computerized X-ray tomography is based on measuring attenuation of photon intensity. Photons can be replaced with neutrinos as mentioned in the introduction, and computerized tomography with muons has been proposed in [58]. We emphasize that neutrino absorption tomography (as discussed in [46,22,45,66]) is significantly different from neutrino oscillation tomography discussed in this article. For physical articles on neutrino oscillation tomography, we refer to [35,28,25,66]. Neutrino tomography could complement existing geophysical imaging methods as discussed in the cited articles.
Of other existing imaging methods based on coherence we mention optical coherence tomography [51] and quantum optical coherence tomography [60]. There are also other kinds of inherently quantum mechanical imaging modalities which often go by the name quantum tomography (see eg. [52,16,30]).
Our imaging method is also related -from a physical point of view -to inverse scattering problems for the Schrödinger equation (see eg. [34,9]). We study quantum mechanical particles that obey the Schrödinger equation, but scattering is not a very significant phenomenon for the present imaging problem. Moreover, in our problem quantum mechanical time evolution takes place in a finite dimensional state space instead of a Hilbert space of functions.

Neutrino oscillation.
For an introduction to neutrino physics and in neutrino oscillation, we refer to [15,14,23,20,59,1,29,2]. We restrict our attention mostly to a quantum mechanical description of neutrino oscillation; quantum field theoretical descriptions can be found, for example, in [15,1].
There are three kinds of neutrinos (ν e , ν µ , ν τ ) and the state space (excluding motion) of a neutrino is thus a three dimensional complex space. We think of this space as C 3 and identify the aforementioned neutrino types with the standard basis vectors e 1 , e 2 , e 3 . This is the so-called flavour basis, and the corresponding basis states (ν e , ν µ , ν τ ) are known as flavour neutrinos. A neutrino state is a superposition of different flavour neutrinos.
Neutrinos are typically produced as flavour neutrinos, so in our quantum mechanical model I ⊂ {e 1 , e 2 , e 3 }. Neutrinos are typically also observed as flavour neutrinos; if all flavours can be measured, we have F = {e 1 , e 2 , e 3 }. Flavour basis is the natural choice if we only consider creation and annihilation in weak charged current reactions (mediated by W ± bosons); if weak neutral current processes are included, the picture is more complicated.
It is sometimes convenient to study neutrino physics in a two-flavour model and it is possible that there are more than three neutrino flavours (the additional ones are known as sterile neutrinos). If there are sterile neutrinos, they cannot be produced or observed directly, so the basis vectors e 4 , e 5 , . . . will not appear in I and F .
Let us then discuss the Hamiltonian related to neutrinos. There is a hermitean 3 × 3 matrix M describing neutrino masses. There is a matrix U diagonalizing M , so that where m 1 , m 2 , m 3 are the neutrino masses. The eigenstates of the matrix M are known as mass states. It is of crucial importance in neutrino physics that the mass states are not the same as the flavour states. The flavour neutrinos do not have well-defined masses, but suitable linear combinations of them (the mass states) do. The relation of flavour and mass states is expressed by the lepton mixing matrix, the Pontecorvo-Maki-Nakagawa-Sakata matrix U . Suppose the neutrino is created at a momentum 2 p > 0 which is much larger than any of the masses. If the neutrino had mass m ∈ R, then its energy would be E = p 2 + m 2 ≈ p + m 2 2p . Thus in a good approximation we may write the Hamiltonian matrix (recall that the Hamiltonian is the energy operator in quantum mechanics) as H = pI + 1 2p M 2 , where M is the mass matrix introduced above. Since multiples of the identity are irrelevant for the time evolution (they only introduce a total phase), we can consider the Hamiltonian This simple Hamiltonian gives rise to a remarkable physical phenomenon: neutrino oscillation. If a neutrino is produced in a flavour state, the probability that it is observed in another flavour state (known as transition probability) oscillates in time. It is evident from equation (1) that the key parameters for neutrino oscillations are the PMNS matrix U and the differences of the squared neutrino masses. The momentum p in (1) can be replaced with the energy E, as these are almost the same for ultrarelativistic particles.
This oscillation is produced by quantum mechanical interference of different neutrino mass states. Interference is only possible if the different mass states are coherent and in the same place. Different masses give rise to slightly different speeds, and different mass states drift apart slowly, whence coherence is lost and neutrino oscillation comes to an end. Transition probabilities stop oscillating, but the asymptotic probabilities are not trivial (not only zero or one). For coherence and decoherence in neutrino oscillations, we refer to [20,Section 4]. (We wish to point out that lost coherence can be recovered at detection; see [1,Section 5.4].) Neutrino oscillation is a kinematic effect; it requires no interactions. The neutrino oscillation described hitherto is the kind that takes place in vacuum. Interactions with the medium change the Hamiltonian of (1) and also the oscillation. It is via these changes that we hope to determine properties of the medium using neutrino oscillation.
If the neutrinos propagate in matter with electron number density N e , we have to add the potential term (with a minus sign for antineutrinos) to the Hamiltonian of (1), neglecting again multiples of the identity matrix. Here E is the neutrino energy and G F is the Fermi constant. This particular form of the potential is due to only electron neutrinos interacting with electrons in the medium via weak charged current; the neutral current potentials of electrons and protons cancel each other and the one of neutrons is a multiple of identity. For details of neutrino oscillation in matter, we refer to [15,Section 9.2].
Because of the form of the potential in (2), neutrino oscillations provide a means of measuring electron density as a function of position. Replacing momentum with energy, we arrive at the total Hamiltonian where f is a scaled electron density. The problem is to recover the function f from the data described in section 2, possibly with measurements for several energies E. We emphasize that ideal data is not available, so theorem 5.4 is not applicable.
From the point of view of integral geometry, this corresponds to replacing scalar fields by (sums of) tensor fields, the order of the polynomial being the order of the tensor field. A particularly simple case to analyze is a first order polynomial (corresponding to a sum of a function and a one-form). For more details on tensor tomography, we refer to [40,53].
The polynomial assumption is also natural physically, and we illustrate this with examples related to particles in an electromagnetic field. The Dirac equation in R 3 for a fermion of charge q and mass m can be written as i∂ t Ψ = HΨ with the Hamiltonian 3 , the 4 × 4 matrices α 1 , α 2 , α 3 are related to the Dirac gamma matrices by α i = γ 0 γ i , A is the magnetic vector potential and ϕ is the electric scalar potential. The Hamiltonian depends on the momentum p, which can be turned into dependence on direction and energy using the identity p = Ev. The nonrelativistic counterpart of the Dirac equation is the Pauli equation with the Hamiltonian acting on a spinor Ψ ∈ C 2 , where B = ∇ × A is the magnetic field, B · σ = B 1 σ 1 + B 2 σ 2 + B 3 σ 3 and the σs are the Pauli matrices. There are several possible choices for the matrices appearing in (3) and (4). For completeness, we give one here: where I 2 is the 2 × 2 identity matrix. The mathematically oriented reader can neglect positioning of indices (upper or lower); we used some upper indices with the Dirac equation and Dirac matrices to be consistent with notation in particle physics. The natural problem related to these equations is to recover the electromagnetic potentials ϕ and A from the data described in section 2. The weighted ray transform described in section 5.1 has to be generalized further: we have to consider matrix valued tensor fields, or functions SM → C N ×N that are polynomial in the fiber variable v. As with all tensor tomography problems, we expect that there is some gauge freedom. But this is to be expected on physical grounds as well, since one can only ever measure a (magnetic or other) potential up to gauge.
We state this tensor tomography problem separately:

Gauge conditions for measurements
Given I , F ⊂ C N , how much about U ∈ U (N ) can be inferred from the knowledge of |a * U b| for all a ∈ F and b ∈ I ? The best we can hope for is knowing U up to a factor e iϕ , ϕ ∈ R, which means ideal data. Proof. We need to show that if for two matrices U, V ∈ U (N ) we have |a * U b| = |a * V b| for all a, b ∈ C N with |a| = |b| = 1, then U = λV for some λ ∈ C, |λ| = 1. By a scaling argument we can drop the assumption |a| = |b| = 1.
It follows from the condition that for any a ∈ C N a vector is orthogonal to U a if and only if it is orthogonal to V a. Therefore there is λ a ∈ C so that V a = λ a U a. We wish to show that all the constants λ a , a ∈ C N \ {0}, are in fact equal.
If a, b ∈ C N are linearly dependent, then by linearity λ a = λ b . If a, b ∈ C N are linearly independent, then Since U a and U b are linearly independent, we have λ a = λ a+b = λ b . Problem 4.2. How large do the sets I , F ⊂ C N have to be to give ideal data? The assumption of the previous lemma looks like an overkill.
The following lemma shows that if we assume the Hamiltonian to be trace free (the trace is invisible from the data anyway), then in the case of ideal data we can recover the time evolution operator U γ H for every line γ without any gauge ambiguity. After this reduction we may restrict ourselves to trace free Hamiltonians, special unitary time evolution and most importantly an inverse problem with no obvious gauge invariance.
Since we have ideal data, we know U γ H up to phase (by lemma 4.1), and since A is traceless, the matrix U γ A is special unitary. Therefore we know the matrix U γ A up to multiplication by an element of the group G N ⊂ U (1) of N th roots of unity.
We know that for very short lines γ the matrix U γ A is very close to the identity, so by finiteness of G N we can recover the matrix U γ A itself (without any gauge ambiguity) for sufficiently short lines. By convexity of the domain every line with endpoints on the boundary can be continuously deformed into an arbitrarily short one; it then follows from continuity of H and finiteness of G N that we can in fact recover U γ A for all γ. By assumption, this information is enough to determine A, the trace free part of H.
The argument in the proof of the lemma above should also work with piecewise continuity, but we refrain from pursuing optimal regularity.
It is a simple exercise to show that if the sets I , F ⊂ C N give ideal data and U, V ∈ U (N ), then also the pair (U I , V F ) gives ideal data. From this it follows that if the tomography problem can be solved on a domain ω ⊂ R n and the trace free part of the Hamiltonian is known in Ω \ ω for some larger set Ω ⊂ R n , then the problem can also be solved in Ω. If the particle travels on a concatenation of several curves, the total time evolution operator is simply the composition of time evolution operators over each segment. Without ideal data there is more gauge freedom than a simple phase, and it is less obvious to find such layer stripping results. (1) |a * U b| = |a * V b| for all a ∈ F and b ∈ I .
1 =⇒ 2: Suppose |a * U e 1 | = |a * V e 1 | for both a ∈ {e 1 , e 2 }. Since the columns of U and those of V are orthonormal and the dimension is two, we have also |a * U e 2 | = |a * V e 2 | for both a ∈ {e 1 , e 2 }. Let It is easy to verify that this is indeed possible due to the constraint λµ = κν. We may linearize (w.r.t. the phase of each matrix element) to find the dimension of the tangent space of the full symmetry group and thus the dimension of that group. To simplify matters, we assume each matrix element to be nonzero. Experiments show that the group of symmetries preserving the norms of elements of a generic 3 × 3 unitary matrix has dimension 6, not 5. We have been unable to identify the group itself in a spirit similar to condition 2 of lemma 4.4.
The following problem needs to be solved before one can study our main problem for N ≥ 3 with I = F = {e 1 , . . . , e N } the standard basis of C N . This problem with N = 3 is relevant for neutrino oscillations, see section 3.1. Is there a description as simple as that of lemma 4.4 for N = 2?

Pseudolinearization
Here we introduce a so-called pseudolinearization that reduces our nonlinear problem to a linear one with unknown parameters. It is enough to know some general characteristics of these parameters (certain weight functions) to show injectivity of the linear problem. It also gives rise to an iterative reconstruction scheme which we shall discuss after theorem 5.2. A similar idea has been previously used for rigidity problems [56,57].
We will phrase our results in the setting of continuous unknown functions. If we assume more regularity, we can assume more regularity in the weights of the corresponding ray transform problems. As it will shortly turn out, pseudolinearization leads to this integral transform. Weighted and especially attenuated X-ray transforms have been studied before (see eg. [12,27,3,50,48,4,31,32,63]), but our results about X-ray transforms with matrix weights are new to the best of our knowledge. The matrix weight is, in a sense, an exponential, and could therefore be regarded as an attenuation (see remark 5.3).
The attenuated X-ray transform with matrix attenuation is related to recovering a matrix valued connection from is parallel transport; for results about connection problems, see [10,17,39,36].
The (scalar) weights that arise in the proof of theorem 5.4 are not attenuations. We have results for more general weights in section 7 and we do not wish to pursue optimal regularity at the expense of clarity, so we will not study the attenuated case further here.
For This gives rise to an improved guess H 1 for H ∞ , and one may proceed iteratively to find successive approximations H 2 , H 3 , H 4 , . . . using the same method. Then hopefully H k → H ∞ as k → ∞ in a suitable sense. This idea been used to solve inverse problems numerically via pseudolinearization in [6,5].
Remark 5.3. The weight that arose in the proof of theorem 5.2 is of a special kind. Identifying C N ×N with C N ⊗ C N , the weight is W (x, v) = A(x, v) ⊗ B(x, v) * . Let us define the differential operator X = v · ∇ x acting on functions on SM ; this operator is known as the geodesic vector field and it generates the geodesic flow. Comparing equations (7) and (5), we observe that XA = iAH and XB = −iBH, whence XW = iW Q, where Q = H ⊗ I + I ⊗H (which is hermitean). The weight can therefore be regarded as an attenuation. We will prove an injectivity result for weighted X-ray transforms even for very rough weights in section 7, and we need not assume that the weight is an attenuation. A more detailed discussion of X-ray transforms is given in section 7.
The following result is based on weighted X-ray transforms with scalar weights. We present the theorem in the amount of generality allowed for by our result on local weighted X-ray transforms, theorem 7.3. Results for such transforms are known also on manifolds [13,57], so the result can be easily generalized. Suppose two functions f,f ∈ C 1,α (M ) give the same data. As in (6), we get for all lines γ : [0, T ] → M with endpoints on ∂M . If the line goes through x ∈ M and has direction v ∈ S n−1 at time t, let us denote By (8) know that the scalar function g = f −f has vanishing weighted X-ray transform for each of the weights w ij . Since the matrix whose coefficients are w ij is nonzero, the weights cannot vanish simultaneously. The weights are also C 1,α on SM , so every point (x, v) ∈ SM has a neighborhood where one of the weights w ij is bounded away from zero. In particular, if x ∈ ∂M and v is tangent to ∂M , there is such a neighborhood. Furthermore, since g is real valued, its weighted X-ray transform vanishes also with the weight Re(λw ij ) for any λ ∈ C; we can thus always locally find a strictly positive real weight for the X-ray transform. Now, using theorem 7.3 with N = 1, we observe that f −f = 0 in some neighborhood of x. We may iterate this argument (a similar layer stripping method was presented in [62]) to show that f −f = 0 in all of M . The method presented in [62] requires a foliation condition that is satisfied by a strictly convex Euclidean domain. (Instead of iteration, one can also argue by contradiction as in the proof of theorem 7.2.) 5.3. Linearization. Since we have been able to give some answers to the full nonlinear problem 2.1, we will only briefly discuss its linearization.
The linearized problem will lead to a similar weighted X-ray transform problem that arose in the proof of theorem 5.2. The difference is that the weight will be known (dependent on the reference Hamiltonian at which we linearize) and therefore easier to invert. For completeness, we prove below that this is indeed the correct linearization. Linearization may be a relevant approach when the Hamiltonian is known to be close to a known Hamiltonian or when one wishes to track (slow) time dependence of the Hamiltonian.

Unordered time evolution
In section 2 we introduced the time evolution operator and mentioned that it could be expressed as a time ordered exponential. We will now investigate what happens if we replace the time ordered exponential with the usual matrix exponential. This problem is no longer quantum mechanically relevant, but it makes an interesting mathematical inverse problem all the same.
If the Hamiltonian H satisfies [H(x), H(y)] = 0 for all x, y ∈ Ω, then the time ordered and the usual unordered exponentials of all integrals of H coincide. We do not assume this commutator property in this section, and therefore the time ordered and regular matrix exponentials are different.
We have a set of initial states I ⊂ C N and a set of final states F ⊂ C N and for all lines γ : (0, T ) → Ω with endpoints on ∂Ω we know the inner products for all Φ ∈ F and Ψ ∈ F . The only difference to problem 2.1 is that T exp is replaced with exp.
In the case of ideal data this means that the matrix can be recovered up to a phase. For ideal data, the unordered problem is now this: This problem quickly leads to this auxiliary problem: Problem 6.2. Let A, B be two hermitean N × N matrices. When do we have e iA = e iB ? Can we classify the pairs of matrices that have the same exponential?
We give a solution to this problem in proposition 6.3 below, but, unfortunately, this result is difficult to work with. Therefore instead of the full problem we treat the linearization. Dropping time ordering from the exponential makes a huge difference in the linearization, as we shall see below. The answer to problem 6.1, linearized at any constant Hamiltonian, is given in theorem 6.13 below.
6.1. When do two matrices have the same exponential? some integer N ≥ 1. For an N × N matrix A, let us denote its eigenspace corresponding to the eigenvalue λ ∈ C by E A λ = ker(A − λI). Let us also define the "periodic eigenspaces" With these spaces we can formulate a simple answer to problem 6.2.
Proposition 6.3. Let N ≥ 1 be an integer and A and B two hermitean N × N matrices. Then the following are equivalent: (1) e iA = e iB , (2) F A λ = F B λ for all λ ∈ R. Proof. 2 =⇒ 1: Let 0 ≤ λ 1 < λ 2 < · · · < λ n < 2π be the numbers for which dim(F A λ k ) > 0 and denote the corresponding spaces by S k = F A λ k = F B λ k . Now C N is an orthogonal direct sum of the spaces S 1 , . . . , S n and the matrices A and B are block diagonal in this decomposition (that is, A(S k ), B(S k ) ⊂ S k for all k). Let A k , 1 ≤ k ≤ n, denote the corresponding block of A, and similarly B k , U k and V k for B, U = e iA and V = e iB .
Considering kth block only, we have e iA k = U k and e iB k = V k . All eigenvalues of A k are in λ k + 2πZ, so the only eigenvalue of U k is e iλ k . Similarly, the only eigenvalue of V k is also e iλ k , so U k = e iλ k I = V k . Since U k = V k holds for each k, we have arrived at U = V .
1 =⇒ 2: Let us show that F A λ = E e iA e iλ for all λ ∈ R. The desired conclusion follows from this claim.
The inclusion F A λ ⊂ E e iA e iλ is elementary and holds also in the case E e iA e iλ = 0. But . This concludes the proof of the claim. Let Using this power series we can also define the operator ϕ(ad A ) : SH → SH.
The solution of problem 6.2 given by proposition 6.3 is difficult to work with, so we linearize the problem. For this we will need the derivative of the exponential. The following two lemmas are well known, but we record them explicitly for the sake of an easy reference. The rest of this section is devoted to lemmas which we will need in section 6.3.
Lemma 6.5. The derivative of the exponential map is d exp A = e A ϕ(ad A ). Lemma 6.6. Ad exp(A) = e ad A for any A ∈ SH. Lemma 6.7. For any A ∈ SH we have SH = ker ad A ⊕ im ad A .
Proof. We can assume A to be diagonal. We fix a basis so that we can write A in the block form where λ 1 < λ 2 < · · · < λ m are the distinct eigenvalues of −iA and I k is the k × k identity matrix. The kernel ker ad A consists of block diagonal matrices whose blocks are of size k 1 × k 1 , k 2 × k 2 , . . . , k m × k m (ie. with block structure corresponding to that of A). For a complex k 1 × k 2 matrix C, let us denote (12) B It now suffices to show that the matrix B C , for any C, is in the image of ad A ; all other matrices with zero block diagonal (corresponding to the block structure of A) can be expressed as sums of matrices of this form. Now a simple calculation shows that [A, We always have ker ad A ∩ im F A = {0}. Furthermore, SH = ker ad A ⊕ im F A if and only if no two eigenvalues of −iA differ by a number in 2πZ\{0}.
Proof. Let us write A again in the block form (11). It is clear that ker ad A consists of block diagonal matrices (with block structure corresponding to that of A). Let us see why such matrices are not contained in the image of F A . Let B C denote the off-diagonal matrix of (12) as in lemma 6.7. (Again we can without loss of generality only consider the first two eigenvalues; other pairs can be dealt with similarly.) Now a calculation shows that (13) F Hence im F A contains only matrices with zero diagonal, and ker ad A ∩ im F A = {0}.
Let us then suppose that the condition on eigenvalues of −iA holds.
This implies that SH = ker ad A ⊕ im F A .
Suppose then that SH = ker ad A ⊕ im F A . This implies that for each k 1 × k 2 matrix D there is a matrix C of the same dimensions so that F A B C = B D . Due to the identity (13) this is requires that e iλ 2 −iλ 1 = 1, so λ 1 − λ 2 / ∈ 2πZ. Lemma 6.9. We have for any A ∈ SH, B ∈ ker ad A and C = ad A D ∈ im ad A .
Proof. Using lemma 6.5 we have Now ad A B = 0 so ϕ(ad A )B = B. By lemma 6.6 we have Combining these gives the claimed equation.
Lemma 6.10. For any A ∈ SH the following are equivalent: (1) The derivative d exp A is a bijection.
(2) All matrices that commute with e A commute with A.
Proof. 1 ⇐⇒ 3: By lemma 6.7 we have SH = ker ad A ⊕ im ad A . Let B ∈ ker ad A and C ∈ SH. By lemma 6.9 e −A d exp A (B + ad A C) = B + C − e −A Ce A . Thus by lemma 6.8 d exp A is bijective if and only if condition 3 holds. 2 ⇐⇒ 3: Let us change the basis so that A comes to the block form where λ 1 < λ 2 < · · · < λ m are the eigenvalues of −iA and I k is the k×k identity matrix. Then Now it is clear that there are matrices that commute with e A but not A if and only if e iλ i = e iλ j for some i = j. But this is equivalent with λ i − λ j ∈ 2πZ.
Lemma 6.11. For any A ∈ SH the set S A = {t ∈ R; d exp tA is not bijective} is countable and discrete.
Proof. Let σ ⊂ R be the set of eigenvalues of iA. Then by lemma 6.10 This set is clearly countable and discrete.
6.3. Solution of the linearized problem. We are now ready to solve a linearized version of our problem. Physically linearization amounts to assuming the Hamiltonian to be very close to some known Hamiltonian H 0 . We make the additional assumption that H 0 is a constant matrix, independent of the position x ∈ Ω. We believe that the linearized problem can be solved in a similar way even for nonconstant H 0 (at least generically if not for every H 0 in every domain).
Denoting by H the small deviation from H 0 , we see from equation (10) that the linearized problem is to recover H from the knowledge of up to a multiple of e −iT H 0 . We will need a little geometrical lemma to solve this problem. Lemma 6.12. Let Ω ⊂ R n be a domain. Then for every r ∈ R the set A r = {(x, y) ∈ ∂Ω × ∂Ω; d(x, y) = r} ⊂ ∂Ω × ∂Ω is closed and has no interior points.
Proof. Closedness follows from the continuity of the distance function d : ∂Ω × ∂Ω → R. Since A r = ∅ for r < 0 and A 0 is the (closed) diagonal of ∂Ω × ∂Ω, it remains to consider r > 0.
Suppose (x, y) is an interior point of A r for some r > 0. Then there is δ > 0 so that d(z, w) = r for all z ∈ B(x, δ) ∩ ∂Ω and w ∈ B(y, δ) ∩ ∂Ω. Fixing z = x shows that near y the boundary ∂Ω is a sphere centered at x; and similarly near x it is a sphere centered at y. Let L be the line containing x and y. Let L be a line parallel to L with a distance less than δ to it. Let z be the unique point in the intersection L ∩ B(x, δ) ∩ ∂Ω and similarly w ∈ L ∩ B(y, δ) ∩ ∂Ω. Then clearly d(z, w) < r, a contradiction. Proof. Let us denote by A the space of trace free hermitean N ×N matrices. We can then decompose the unknown matrix function uniquely as H(x) = A(x) + h(x)I where A(x) ∈ A and h(x) ∈ R. Both A and h are continuous functions.
By lemmas 6.11 and 6.12 we know that d exp d(x,y)H 0 is injective for a dense set of pairs (x, y) ∈ ∂Ω × ∂Ω. By continuity and density, we can recover the integral T 0 A(γ(t))dt for every line γ : [0, T ] →Ω. (Note that by lemma 6.9 e iT H 0 d exp −iT H 0 maps multiples of identity to multiples of identity.) By injectivity of the X-ray transform in Euclidean domains (see eg. [19,33,44,7]), we can recover A everywhere in Ω.
There may be lines for which the corresponding derivative of the exponential is not injective. The set of singular lines has no interior points so we were able to use a density argument, but we suspect that this renders recovery somewhat unstable. A stability result for the linearized problem would give a local uniqueness result for the full problem.

X-ray transforms with matrix weights
Let M be the closure of a compact, smooth domain in R n . Let SM = M × S n−1 denote the unit sphere bundle on M and L the manifold of maximal directed lines in M , which is the quotient of SM by the geodesic flow. Let π M : SM → M and π L : SM → L be the projection and the quotient map.
Let W : SM → C N ×N be a continuous matrix function. We define the weighted X-ray transform I W f : L → C N of a continuous function f : M → C N by We will first restrict our attention to the case n = 2 and later use results from this case for higher dimensions.
For any p ∈ M such that W (p, v) is invertible for all v ∈ S n−1 there is a continuous function Q p W : L → C N ×N so that (π * L Q p W )W = I on the fiber S p M = {p} × S n−1 . This will allow us to make the weight trivial at any single chosen point p ∈ M in all directions.
Let I α denote the Riesz potential The size of the constant c α is irrelevant for us, but we remark that c 1 = 1/2π. Let P be the back projection operator defined on functions g : L → C N by P g : M → C N , We will invert I W for suitable weights W by freezing the weight at a point near the boundary and using the identity P I I f = I 1 f and injectivity of the Riesz potential in suitable spaces. This gives a proof of injectivity in two dimensions in a sufficiently small domain (where the weight is close enough to the identity matrix), which in turn gives a more general injectivity results in dimension three or higher. (Injectivity in small domains was proved for the attenuated X-ray transform in [11] with smooth scalar attenuations.) Our proof is based on the ideas of [27], but we improve it in two important ways. We allow rougher weights; in particular, we need not assume that the weight has even a first derivative. Most importantly -from the point of view of problem 5.1 at least -we allow the weight to be any invertible matrix. A similar problem for X-ray transforms with matrix weights has been studied before by several authors [36,18,24,37,63]. There is also a recent preprint [42] with very general results for matrix weights on manifolds.
Injectivity of I W requires that the weight W is sufficiently regular. Our regularity assumption is quite relaxed, as we do not require differentiability. Assumptions on the weight are less stringent in dimension three and higher because the problem is formally overdetermined. There are counterexamples to injectivity if the weight is not regular enough; see [27,3].
Take any β > 0 and consider a weight W which satisfies W ∈ C β (M, Lip(S 1 , C N ×N )) and W (0, v) = 0 for all v ∈ S 1 . There are p ∈ (1, ∞) and ε > 0 so that for any δ > 0 we have for any f ∈ L p (B(0, δ) ∩ M, C N ) where C p is a constant depending only on p.
Proof. Let us denote by inequalities including constants that only depend on the exponents. Let us assume that f is supported in B(0, δ).
Let i ∈ {1, 2} and j = 3 − i. A calculation gives (We may use a simple approximation argument to justify the above calculation with the assumed regularity of W and f .) Therefore, if K(x + y) 1 |y| |f (x + y)| dy.
We therefore want to have K ∈ L s (M ) for some s > 2.
We are now ready to give an answer to problem 5.1. Using theorem 5.2, we can determine the trace free part of the Hamiltonian H from our quantum measurements, provided that H ∈ C 1,β for some β > 0. This was stated as theorem 2.3. Proof. Suppose f = 0 but I W f = 0. Take p in the intersection of spt(f ) and the boundary of its convex hull so that there is a subtangent plane of the convex hull at p touching the support only at p. We shrink M to be this convex hull. Let ν ∈ S n−1 a vector pointing inward from p.
For t > 0 denote x t = p + tν and suppose t is so small that x t ∈ int M . Let G be the set of two dimensional subspaces of the chosen tangent plane at p. Denote M T t = (x t + T ) ∩ M for t > 0 and T ∈ G. Let δ t be the maximum of diameters of M T t , T ∈ G. We will consider the weighted X-ray transform of f on each M T t , and for the sake of simplicity we will consider T ∈ G fixed and suppress it from the notation. To that end, let us denote the restrictions of f and W to M t by f t and W t , respectively, and the set of lines on M t by L t .
Define Q t : S 1 → C N ×N by Q t (v) = W t (x t , v) −1 . This is a Lipschitz function, and we may define another Lipschitz functionQ t : L t → C N ×N by lettingQ t ( ) = Q t (v ), where v ∈ S 1 is the direction of the directed line . The Lipschitz constant ofQ t is uniformly bounded in t.
Since I Wt f t = 0, we have I (π * L tQ t)Wt f t =Q t I Wt f t = 0. Let us denotẽ W t = (π * LtQ t )W t . The function π * LtQ t : M t × S 1 → C N ×N is Lipschitz (uniformly in t) and only depends on the direction, so we have = ∇P IW t f t + ∇P I I−Wt f t p ≤ ∇P IW t f t p + ∇P I I−Wt f t p ≤ ∇P IW t f t p + cδ ε t f t p for a sufficiently large fixed p and some constants C, c, ε > 0. Note that this estimate is uniform in T ∈ G. Since δ t → 0 as t → 0, there is t 0 > 0 is such that cδ ε t < C/2 for all t < t 0 . For t < t 0 we get thus But IW t f t = 0 so f t = 0 for t < t 0 . This means that f = 0 in 0≤t<t 0 T ∈G M T t , which is a neighborhood of p. Due to our choice of p, this is impossible. Therefore I W f = 0 is only possible if f = 0.
It is not actually necessary that W (x, v) is invertible for all (x, v) ∈ SM . The proof above only used invertible for v in a subspace homeomorphic to S n−2 of S n−1 , given x. If one wants to recover a function from its integrals over some set of lines L ⊂ L, it can only be done stably if the normal bundle of L covers SM (see eg. [43,55,13]).
Inspecting the proof of theorem 7.2, we see that we have in fact also proven the following local result. A similar result on manifolds (with n ≥ 3, N = 1 and W ∈ C ∞ ) was shown by Stefanov, Uhlmann and Vasy [57, Corollary 3.2]. As mentioned above, an earlier Euclidean version can be found in [27] (for C 2 scalar weights) and in [55] (for C ∞ scalar weights using microlocal tools). Theorem 7.3. Let M be the closure of a bounded convex C 1 domain in R n , n ≥ 3, and let β > 0. Let W ∈ C β (M, Lip(S n−1 , C N ×N )) be a weight. Fix p ∈ ∂M . Suppose p has a neighborhood U so that W (x, v) is invertible for all x ∈ U and v in a neighborhood of tangential directions at p. Suppose for f ∈ C(U, C N ) we have I W f (γ) = 0 for all lines γ that stay in U . Then there is a possibly smaller neighborhood U ⊂ U of p so that f | U = 0.
The proof also gives a stability estimate for the local problem; cf. (15).