The inverse transport theory has been explored for decades, and its applications are found in nuclear engineering, atmospheric radiation, remote sensing, astrophysics, biomedical imaging, etc. Mathematically, the most common measurement assumption is the albedo operator; the full albedo operator naturally implies a Schwartz kernel decomposition concerning the singularity strengths. There are still a few open problems with the radiative transfer model.
For the albedo operator, when the measurement is angularly resolved or the input source is angularly dependent, the reconstructions are stable for the ballistic contribution, namely the total absorption coefficient. Still, the angularly averaged measurement with an isotropic source will not have the so-called Schwartz kernel decomposition. The most critical question is whether there are any general uniqueness results or stability estimates for the reconstructions. A more intriguing question is: what kind of mathematical tool could be the key to this problem? There are a few efforts in this direction already [Bal09, BVLM11, BLM08, BJ09, Lan08, ZZ18], and some related works [LLU18, CLL17], but the problem remains largely unresolved.
In quantitative photoacoustic tomography (qPAT), the radiative transfer model is assumed to be accurate in many biomedical environments, such as breast tissues or bones. The theory is already discussed in [BJJ11, MR12]. In practice, the coefficients are assumed to be isotropic for simplicity, see the settings in [MR12]. The scattering-free case is relatively straightforward, but the scattering medium is only studied with the complete Albedo-type measurement, namely, there is angular dependence in the input data, which is unrealistic. With small scattering, the problem is solved through a fixed-point iteration, but the general case is still open. Furthermore, if the source is isotropic, is it possible to come up with a similar theory, such as diffusion theory in the qPAT? [partially solved in RZ25 if the Gruneisen and scattering coefficients are known, it is still unknown whether a large scattering coefficient can be recovered yet. ]
In terms of the forward model, the numerical schemes to solve the radiative transport equation are studied in extensive literature. The SPn and Pn methods were popular in the early days, which used high-order diffusion approximation to solve the spherical harmonic moments of the solution. However, the approximation error is only rigorously proved in the diffusion limit. In the diffusion regime, there are many linear algebraic fast algorithms to deal with the H-matrices, but for the transport regime, the intrinsic complexity is still unknown.
The SPn method and Pn method are locally equivalent for a 1D-like region. It is shown that these two settings are different, and one does not even possess convergence as the order approaches infinity. All the evidence shows that the zeroth moment/flux should live in a relatively low-dimensional space rather than solving the complete phase space equation. If this is not true, then what kind of condition should the source possess?
For qPAT, even in the "almost" diffusion regime, the diffusion models sometimes appear inadequate for low-order expansion; analysis for higher-order expansion is only conducted in the SP1 and SP3 formulations. In the general formulation, the uniqueness and stability will deteriorate, but the analysis is still absent. [Solved. The SPn method's stability decays algebraically.]
As an extension of the intrinsic complexity for the radiative transfer equation, since the flux is more important in practice, and enjoys better regularity, it should have better computational complexity if we make certain analogs to the diffusion model. Is there any fast algorithm for the anisotropic radiative transfer equation? The brute-forced compression is already done, but it isn't easy to generalize into high-order expansion.
Is there any compression like the isotropic radiative transfer (Peierls integral equation) for applications like black-body radiation or other models? How to compress the boundary conditions?
Numerically solving the Peierls integral equation will involve the weak singularity of order (1-d), which implies the boundary layer effect, meaning the regularity will deteriorate quickly as one approaches the boundary. Is there a way to get rid of the fine discretization on the boundary layer but achieve high accuracy in the flux? [It seems not possible to avoid high resolution in the boundary layer unless scattering smoothly varies to zero at the boundary, which makes the problem globally defined. ]
In the time domain, the Peierls volume integral will become a conical volume integral. This cubature could be solved with interpolation, but the interpolation will be the most expensive step in the algorithm. Is there an easier way to get the cubature calculated? [Solved. The complexity would match the frequency-marching, but frequency-marching will lose the causality.]
In the frequency domain, the volume integral should have a natural butterfly-decomposition structure like the Helmholtz equation. Is it possible to apply such a method to compute the optical tomography in the transport regime? The tricky part will be the inversion process; the fast-forwarding operator does not relieve the burden of inversion.
In the setting of Problem 1, the numerical reconstruction is complex as well for the scattering coefficient (the absorption coefficient is easier); however, numerical evidence shows that the back-scattering could play a role in getting the singularity (jump) in the scattering coefficient through the assumption that the multiple-scattering part is much smoother than the single-scattering and scattering-free parts. Such evidence could help reduce the number of unknowns if the scattering contrast is locally supported (i.e., it could be scattered). However, it is also essential to recognize that the problem is ill-posed, and the forward solver is not low-dimensional. In other words, the solution to data mapping is the main compression, and the solution space is not the correct space to adopt the so-called subspace methods. Another thing is to combine with Problem 10 in a continuous-wave setting. At the final stage, one cannot avoid solving the transport equation with optimization if a faster convergence rate is desired. An efficient way to update the operator is the most challenging part.
As a side product of Problem 2, when the boundary source is isotropic, the angular-averaged solution likely has a unique continuation property. A straightforward case might be on the unit disk with coefficients as constants. The general case is challenging to address because there are not many tools for the unique continuation property. The scattering-free case is especially a boundary integral problem, which might be solved with the doubling inequality technique. [For constant-coefficient and scattering-free, solved with Beurling's theorem. ]
For the CGO solution with zero Neumann/Dirichlet data on the partial boundary, can we obtain a higher regularity statement instead of those in L2/H1? Without a boundary, there is uniqueness, so differentiating the equation automatically adds regularity to the estimate.
As an extension of 12. Do we have a unique continuation for boundary integral operators, or what kind of boundary integral possesses such a property? I have asked such a question on mathoverflow.net. It seems a quasi-analytic case can be proved if the derivatives of the kernel grow at most poly-exponentially. [partially solved.]
In 2D, given the boundary incoming source function with finite Fourier modes, do the solution's corresponding Fourier modes on the boundary recover the solution? This is true for the isotropic case. The potential theory approach seems to only work for even modes. For reference, see [Da03].
It is known that the unique continuation property can uniquely determine point sources (or singular sources); however, the usual minimization approach is not adequate to guarantee convergence to the target sources due to the non-convex landscape. Is it possible to derive a strategy when meeting a local minimum? For known media, this problem does not seem complicated, especially for homogeneous cases, but for variable media or practically unknown media, the reconstruction is far more difficult. [This is initially for the elliptic equation, but can be generalized with the transport model, see 20.]
The angularly averaged non-scattering RTE involves a boundary integral that is similar to the Neumann-Poincaré operator. Is it possible to adjust the theories of the NP operator to the RTE case (e.g., spectral properties)? The integral can be reformulated with a complex variable in 2D in an exponentially weighted Cauchy integral form. There is very little information on the integral.
The angular-averaged collision term of RTE used in [LYZ19] is not fully understood, although there is a proximal way to perform the reconstruction, which requires a specific boundary source. Is it possible to reconstruct a single homogeneous source? This problem is also considered in [LUZZ22], but with infinite data (linearization technique).
As an analog to 13, is it possible to show a similar result of the power density matrix non-degenerated (inner product of gradients) for averaged solutions of RTE?
As an extension of 16, is the point source argument applied to the transport model with the averaged setting?
As a follow-up to 5, the mass matrix part of the weak formulation should be taken care of, and the eigenvalue decay should be the same as the stiffness part. We should be able to get to a result independent of order. [Wrong. The norm of the solution is changed to "standard L^2 norm", but the constant dependence is at O(N). We need to know whether the diffusion approximation dominates other moments.]
For the integral formulation, the inverse problem is computationally expensive due to the re-computation of the kernel. It will be a good idea to defer the updating by introducing a two-stage update scheme from the hybrid formulation. [There is a hybrid formulation combined with Jacobi iteration to mitigate the issue.]
Given a positive definite 4-tensor A (in the sense of a vectorized matrix), can we find an easy way to compute the left and right bases that are rank-1 and orthogonal to A? Without the rank-1 requirements, then two disjoint sets of eigenvectors will be enough. This relates to a large class of inverse problems.
The forward problem of the nonlinear absorption RTE is only known to be well-posed for small sources; it remains to be seen if this assumption can be lifted. [solved.]
Can we remove the constraints in the coefficients for PAT data when only reconstructing the absorption coefficient? [solved. We can do more complicated cases. ]