Orbit tomography in constants-of-motion phase-space

Tomographic reconstructions of a 3D fast-ion constants-of-motion phase-space distribution function are computed by inverting synthetic signals based on projected velocities of the fast ions along the diagnostic lines of sight. A spectrum of projected velocities is a key element of the spectrum formation in fast-ion D-alpha spectroscopy, collective Thomson scattering, and gamma-ray and neutron emission spectroscopy, and it can hence serve as a proxy for any of these. The fast-ion distribution functions are parameterised by three constants of motion, the kinetic energy, the magnetic moment and the toroidal canonical angular momentum. The reconstructions are computed using both zeroth-order and first-order Tikhonov regularisation expressed in terms of Bayesian inference to allow uncertainty quantification. In addition to this, a discontinuity appears to be present in the solution across the trapped-passing boundary surface in the three-dimensional phase space due to a singularity in the Jacobian of the transformation from position and velocity space to phase space. A method to allow for this apparent discontinuity while simultaneously penalising large gradients in the solution is demonstrated. Finally, we use our new methods to optimise the diagnostic performance of a set of six fans of sightlines by finding where the detectors contribute most complementary diagnostic information for the future COMPASS-Upgrade tokamak.

Tomographic reconstructions of a 3D fast-ion constants-of-motion phase-space distribution function are computed by inverting synthetic signals based on projected velocities of the fast ions along the diagnostic lines of sight.A spectrum of projected velocities is a key element of the spectrum formation in fast-ion D-alpha spectroscopy, collective Thomson scattering, and gamma-ray and neutron emission spectroscopy, and it can hence serve as a proxy for any of these.The fast-ion distribution functions are parameterised by three constants of motion, the kinetic energy, the magnetic moment and the toroidal canonical angular momentum.The reconstructions are computed using both zeroth-order and first-order Tikhonov regularisation expressed in terms of Bayesian inference to allow uncertainty quantification.In addition to this, a discontinuity appears to be present in the solution across the trapped-passing boundary surface in the three-dimensional phase space due to a singularity in the Jacobian of the transformation from position and velocity space to phase space.A method to allow for this apparent discontinuity while simultaneously penalising large gradients in the solution is demonstrated.Finally, we use our new methods to optimise the diagnostic performance of a set of six fans of sightlines by finding where the detectors contribute most complementary diagnostic information for the future COMPASS-Upgrade tokamak.

Introduction
In a fusion plasma, fast ions are energetic particles much more energetic than the bulk plasma.Fast ions are both required and inevitable.They are born as fusion products and they are generated by heating mechanisms, such as neutral beam injection (NBI) and ion-cyclotron range-of-frequencies heating.They are needed to sustain the high temperatures required for nuclear fusion to take place, but can also drive instabilities in the plasma [1][2][3][4][5][6].Fast ions can themselves be transported away from the plasma core by magnetohydrodynamic modes [7][8][9][10][11][12], thus rendering the study of fast ions vital in the development of magnetic confinement fusion devices.
The motion of fast ions is conveniently described by the three constants of motion in the guiding-centre approximation: the kinetic energy E, the magnetic moment µ and the toroidal canonical angular momentum P ϕ , defining 3D trajectories of the ions, called orbits.Confined orbits are closed in the poloidal projection, see e.g.[13,14], whereas lost orbits hit the wall.The fast-ion distribution function can therefore be expressed in the phase space spanned by the constants of motion E, µ and P ϕ with the addition of a binary index σ.
The diagnostic signals and the fast-ion distribution function are often related via so-called weight functions, which describe the sensitivity of the diagnostics in a given phasespace region.The diagnostic signals used in this work are all synthetic, i.e. they are calculated numerically.The projected velocity of the fast ions along the sightlines is used as a proxy for different diagnostic signals, such as fast-ion D α (FIDA) spectroscopy [35][36][37][38][39][40], collective Thomson scattering [41], gamma-ray spectroscopy (GRS) [42,43] and neutron emission spectroscopy (NES) [19].This allows an analytical study of the orbit sensitivity [44].
Weight functions in constants-of-motion phase space have also been calculated numerically for FIDA in [56].However, the work in the present paper builds upon the analytical work on weight functions in constants-of-motion phase space developed in [44].The fast-ion species is here taken to be deuterium.Due to a singularity [57,58] in the Jacobian of the transformation from Cartesian position-and velocitycoordinates, steep gradients will appear in the fast-ion distribution function across the trapped-passing boundary of phase space, if the Jacobian is absorbed in the distribution function, as is often the case [57,58].A key objective of the present paper is to compute the reconstructions in constants-of-motion space directly, since this space allows an analytical description of the weight functions [44], and it is used in stability codes [59,60].Another key point is to demonstrate the ability to reconstruct the expected steep gradient across the trappedpassing boundary of phase-space, while simultaneously promoting a smooth reconstruction in the remainder of phase space, as is expected from the Fokker-Planck equation governing the evolution of the distribution function, due to collisions.We return to this point in section 5.
The remainder of this paper is organised as follows.In section 2, we introduce the forward problem and discretise it to write it as a matrix-vector equation.This section ends with a short presentation on the diagnostic setup assumed in this study.In section 3, we formulate the inverse problem and explain the equivalence between Tikhonov regularisation and Bayesian inference, which allows us to make estimates on the uncertainty quantification.Reconstructions of fast-ion distribution functions in the future COMPASS-Upgrade tokamak using zeroth-order Tikhonov regularisation, penalising large peaks in the solution, are performed in section 4.These reconstructions are followed in section 5 by reconstructions using first-order Tikhonov regularisation, penalising large gradients in the solution, where we also modify the regularisation to allow for expected steep gradients at the trapped-passing boundary.As a practical application of the developed orbit tomography in constants-of-motion phase space, we perform an optimisation study in section 6 on the diagnostic FIDA setup of the future COMPASS-Upgrade tokamak.We conclude and give prospects in section 7.

Forward problem and symmetry constraints
In this section, we consider the general setup of the forward problem and how the symmetries of the system and the dynamics of fast ions allow us to reduce the phase-space dimensions.
M. Rud et al

Symmetry constraints
In the forward problem, a signal is calculated from the fastion distribution via the, in general six-dimensional, linear Fredholm integral equation of the first kind, by integrating over three position coordinates x and three velocity coordinates v.As we take the projected velocity to be a proxy for the signal, the units of the signal is 1/(m s −1 ) [41].
A diagnostic is imagined to measure projected velocities u, as this is a fundamental quantity in several fast-ion diagnostics, such as FIDA spectroscopy, CTS, GRS and NES, as mentioned in the introduction.The weight function is the kernel w (6D) with the unit signal per ion, which quantifies the diagnostic sensitivity to each point in phase space, and the distribution function is measured in ions per phase-space volume.
The upper and lower limit of an interval in the measured spectrum, u 1 and u 2 , reflect the finite spectral resolution of energetic particle diagnostics.The viewing angle between the diagnostic sightline and the magnetic field vector at the spatial location of the measurement is denoted by φ.This angle is essential information, as it may lead to either a blue-or redshift of the signal, depending on the direction of the fast ion, which can lead to asymmetric signals for oblique detectors.The fundamental ingredient in the shift of the signal is the projected velocity u of the fast-ion along the diagnostic sightline [41], where γ is the gyro angle of the fast-ion motion around the magnetic field vector.The velocity component parallel to the magnetic field is denoted by v ∥ , and the velocity component perpendicular to the magnetic field is v ⊥ .The linear relation in equation ( 1) is a good model for FIDA and CTS as well as NES and GRS dominated by beam-target interactions, i.e. when considering the interactions between fast ions and thermal ions in the plasma.Beam-beam reactions are non-linear but often contribute a fairly small amount to the signal [19,61].
Equation ( 1) can be simplified by reducing the dimensionality of the problem, by taking advantage of ignorable coordinates.A detailed derivation of this can be found in [32].Expressing the weight functions and the distribution function in terms of the constants of motion in the guiding-centre approximation, (E, µ, P ϕ ), given by [2,13,14] equation (1) then becomes COM (E, µ, P ϕ ; σ) dP ϕ dµdE.(6) Note that ϕ denotes the toroidal angle coordinate, whereas φ denotes the angle between the diagnostic sightline and the magnetic field vector at the measurement location.The gyro angle γ of the fast ion is assumed to be uniformly distributed in γ ∈ [0, 2π].The fast gyration means that a single point along the fast-ion guiding-centre orbit generates an entire spectrum of signals in the detector, provided that the specific point in phase space is observed by the given detector.The binary index σ is needed in the distribution function, since there is an ambiguity in the direction of the particles relative to the magnetic field.A passing particle moving in the same direction as the magnetic field can have the same (E, µ, P ϕ )-triplet as a particle moving in the opposite direction as the magnetic field.Thus, co-going and trapped orbits are collected in the σ = 1-part of phase space, and will be denoted with positive energies, and counter-going orbits are collected in the σ = −1part of phase space, and will be denoted by negative energies.By expanding the energy range in this way, we ensure that a single point in phase space corresponds to only a single orbit [44,57].
In effect, we have reduced the six-dimensional forward problem (1) to a three-dimensional one by integrating over the ignorable coordinates using the constraints imposed by the symmetries of the system.The Jacobian of the transformation from the Cartesian coordinates to the constantsof-motion coordinates in (6) is absorbed in f (3D) COM .As noted in the introduction, the absorption of the Jacobian introduces a singularity along the trapped-passing boundary, leading to a steep gradient in f (3D) COM on a discrete grid.Alternatively, the Jacobian could be absorbed in the weight functions changing their units, but leaving the distribution function smooth.
Discretisation of (6), allows us to write it in the form of a matrix-vector equation where s ∈ R N×1 is a vector containing all measurements, f ∈ R M×1 is the vectorised fast-ion distribution function, such that each element in the vector corresponds to a specific point in phase space, and W ∈ R N×M is a matrix, where each row is a vectorised weight function.This enables us to conveniently employ different regularisation techniques in the tomographic inversion of the signals to reconstruct the fast-ion distribution function.The number of columns M in W equals the number of rows in f, which is the number of points in phase space, in which the fast-ion distribution function is resolved.The number of rows N in W equals the number of rows in s, which is the sum of the number of bins in which each spectrum is collected.E.g. 40 spectra with 100 bins in each would give N = 4000.Usually, we choose N ̸ = M,  which makes W a real rectangular matrix, which cannot be inverted.Before we turn to an inspection of the weight functions in constants-of-motion phase space, we briefly review the diagnostic setup used for the future COMPASS-Upgrade tokamak.

Diagnostic setup
Figure 1 shows a CAD (computer-aided design) drawing of the sightlines considered in this work, denoted in magenta, observing measurement volumes along an NBI drawn in orange.Five fans with 10 sightlines each are included.The diagnostic ports are located in section 3, 7, 8 and 9, as indicated in figure 1(b).The diagnostic NBI, as well as two diagnostic ports are located in section 8.One port above the NBI, and one below.The fans originating from these ports will be denoted as 8u and 8l, respectively.The radial coordinates of the measurement volumes as well as the angles between all sightlines and the magnetic field in each measurement volume are collected in table 1.The z-coordinate of all measurement volumes is z = −0.01m, as they are located along the NBI beam line.
The diagnostic system could be a FIDA-system, but as clarified above, by using only projected velocities, we keep the investigation general, without referring to a specific type of diagnostic.

Weight functions and spectra
The calculation of the weight functions in constants-of-motion phase space is explained in detail in [44].Here, we present a subset of the weight functions used in the calculation of the synthetic signals.Here, only the sensitivity in the σ = −1 part of phase space is shown, since no co-going orbits of these energies give rise to The sensitivity to trapped orbits is collected in the σ = 1 part of phase space which is not shown here.
such negative projected velocities in this sightline.The sightline is sensitive to some trapped orbits with this u, but they are collected in the σ = 1 part of phase space, and thus not shown in this figure.We refer to both co-passing and co-stagnation orbits as co-going orbits, and to both counter-passing and counter-stagnation orbits as counter-going orbits.
In figure 3 we fix the fast-ion energy to E = 65 keV and compare the sensitivity to projected velocities u ∈   The magnetic equilibrium used in this work is the same as that used in [44], calculated by the FIESTA code [62], with all profiles referenced in [63].The plasma current of 1.6 MA is going clockwise, when viewed from above, parallel to the toroidal magnetic field in an ITER-like configuration.The magnetic field on axis is 4.3 T. The fast-ion distribution that we will take as the true distribution function is calculated using the EBdyna orbit solver [64] corresponding to 2MW of NBI power injected in COMPASS-Upgrade at tangency radius 0.45 m.
The synthetic spectra s are calculated as where e ∼ N (0, C s ) is assumed to be additive Gaussian noise with zero mean and covariance matrix C s , which is here taken to be diagonal, where s i is the ith element of s clean .Note from (9) that if s i is too small, white Gaussian noise is added with a standard deviation of 10 −3 max {s clean }.In figure 4 we show example spectra obtained both with and without noise.The spectrum without noise obtained from the seventh sightline from port 9 is shown by the full black line in figure 4(a).The same spectrum including 10% Gaussian noise is shown by the blue dashed line.All spectra as obtained by all sightlines are shown in figure 4(b).Note that we have removed the signal in the interval u ∈ [−0.7, 0.7] • 10 6 m s −1 in all spectra, as a corresponding experimental signal will not be observable due to the much stronger signal due to the thermal ions and can therefore not be used.We mimic this by removing the thermal part of the spectrum.It is evident that the signal in figure 4(a) is asymmetric about u = 0 m s −1 , which is expected, since the sightline is oriented horizontally with an oblique angle into the plasma.The seventh sightline from port 7 thus obtains a spectrum that is approximately the mirror image about u = 0 m s −1 of this.
In the next section, we invert these spectra to infer the fastion distribution function.

The inverse problem: inferring the distribution function
In this section, the inverse problem is formulated as well as solved to obtain tomographic reconstructions of a fast-ion distribution in the future COMPASS-Upgrade tokamak.The outset is equation (7) in which we know W and s and wish to recover f.Since W can be rectangular, it is not invertible (this is only possible if W is square).The least squares solution to the inversion is where ||•|| 2 is the Euclidean two-norm.However, due to the large condition number of the weight matrix, small changes in the signal s, due to inevitable noise in the data collection (8), can lead to large errors in the solution f LS .Therefore the problem is said to be ill-posed.cf Hadamard's criteria [65].One way to address this issue is to adopt a probabilistic approach based on Bayesian inference.Another way is to regularise the problem and penalise undesired properties of the solution in order to force the solution to change continuously and stably with the signal.In this work we will utilise the well-known fact [66] that for a Gaussian process, with the prior covariance matrix given by the Tikhonov regularisation matrix, the maximum a posteriori (MAP) estimate of the distribution function is exactly the solution found by Tikhonov regularisation, and thus that the two approaches are equivalent.Tikhonov regularisation allows us to penalise properties based on prior information, and the Gaussian process provides quantification of uncertainties in the solution.
We do not use the same model, i.e. discretisation of the distribution function, in the forward and the inverse problem.Doing so is in some mathematical communities referred to as committing an inverse crime [67].This would be equivalent to assuming that our physical model represents reality in an experiment exactly, and that would be a wrong assumption, which we therefore avoid.Since this study is based on synthetic signals, in order to not commit an inverse crime, we interpolate the weight matrix W → W onto a coarser grid, on which we solve the inverse problem.This interpolation introduces interpolation noise in the model, in addition to the random noise put on the signal in (8), thus challenging the inference further.

Bayesian inference and Tikhonov regularisation
The goal is to infer the most probable distribution function f given the noisy measurement s, where the forward problem is formulated in (8), in which e, and thus also s clean , is unknown.For a regularisation matrix L the regularised Tikhonov solution is where the hat notation indicates that the signal s and the weight matrix W have been weighted by the noise level, so we are solving a weighted least squares problem.The regularisation matrix L can be related to the prior covariance matrix of f in a Bayesian perspective.Since the prior covariance matrix is positive semi-definite, we can decompose it as C −1 f = λ 2 L T L, such that finding the MAP estimate of f is equivalent to finding the solution to (11).To make the dependence of the prior covariance matrix on λ explicit, we will from now on denote it by C f λ .Choosing the prior covariance matrix C f λ of the distribution function to be proportional to the identity matrix, C −1 f λ = λ 2 1 M , we find that f λ is the solution to the zeroth-order Tikhonov regularisation, which is one of the methods adopted in this work in section 4.This will penalise large values in the solution.In section 5, we apply first-order Tikhonov regularisation which penalises large gradients in the solution, thus promoting smooth solutions.
For a Gaussian likelihood and a Gaussian prior distribution the mean of the posterior distribution is found at the fλ that maximises the posterior distribution, i.e. the MAP estimate.Thus, fλ is given by [66] and the posterior covariance matrix is The diagonal elements of Cf λ are the variances of the solution fλ at each point in phase-space, and they provide us a quantification of the uncertainties in the solution.
The next problem is to determine the required regularisation strength λ to optimise the solution f λ .Several methods to find λ exist, such as the L-curve [68] and the generalised cross-validation [69].In this work we solve the inverse problem with several different λ's and choose the optimal regularisation strength as the one that minimises the regularisation error ε, given by i.e. the two-norm of the difference between the reconstruction and the true solution, such that λ opt = argmin λ ε(λ).Here, the true solution f true is the true distribution used in the forward problem interpolated onto the coarser grid in order to avoid an inverse crime.This is only possible if we know the true solution, which we do in the numerical simulations, but not in actual experiments.The focus of the present work is not to provide a method to choose λ, but rather to demonstrate the reconstructions in constants-of-motion phase space.As a test of the orbit tomography in constants-of-motion space, we reconstruct a 3D fast-ion distribution corresponding to 2MW of NBI power injected in COMPASS-Upgrade at tangency radius 0.45 m.

Reconstruction of fast-ion distribution function using zeroth-order Tikhonov regularisation
The reconstructions of the distribution function in COMPASS-Upgrade are calculated using all five ports with a fan of 10 sightlines in each.In figure 5 we show how the reconstruction error (14) decreases when additional fans of sightlines are added.In this case we used a clean signal, i.e. with no noise, zeroth-order Tikhonov regularisation and we used the same grid for the forward and inverse problem (inverse crime), such that the only variation is the number of sightlines.The figure should be understood in the way that the first data point with a reconstruction error ∼0.7 is obtained by including only the upper fan from port 8.The second data point is obtained by including both the upper fan from port 8 and the oblique fan from port 7, and so on.Since the upper fan has predominantly perpendicular sightlines, the reconstruction error is greatly reduced by including the fan from port 7, which has oblique sightlines.It is evident that including additional sightlines with different viewing angles improves the reconstruction.All reconstructions using zeroth-order Tikhonov regularisation calculated in this work are calculated using (12).
We now turn to reconstructions using noisy signals.See figures 6(a) and 7(a) for the simulated distribution, which we take to be the true distribution, and figures 6(b) and 7(b) for the reconstruction of the positive and negative energy slices respectively, using zeroth-order Tikhonov regularisation.The reconstruction is calculated using the detectors in port 3, 7, 8u, 8l and 9.In general the reconstruction looks similar to the true distribution, in the sense that high density is found at the correct energies and in the trapped region for σ = 1.Since we remove the projected velocity interval u ∈ [−0.7, 0.7] • 10 6 m s −1 , we are not able to reconstruct the lowest energies, which is expected.An unsatisfactory feature, however, is that relatively high density is reconstructed along the upper boundary of the trapped region, as opposed to the centre of the trapped region.The high density in the trapped region for E = 18 keV is also not reconstructed.This is further illustrated in figure 8(a) where the difference between the reconstruction and the true distribution is shown for σ = 1 and in 8(b), where the standard deviation of the posterior distribution is shown, also for σ = 1.Note that since the 5 keV energy plane is hidden by the thermal feature, the 5 keV energy plane in figure 8(a) is merely showing the negative of the true distribution.The standard deviation is the square root of the diagonal elements of Cf λ , which is a measure of the uncertainty in the amount of ions per phase-space volume.It is evident that large differences from the true distribution occur in places with large uncertainties, i.e. especially in the deeply trapped region across all energies.
The singularity along the trapped-passing boundary challenges further the reconstruction there.Additional challenges are investigated in the following section.

Singular value decomposition (SVD) of the weighted weight matrix
In this section, we utilise the fact that the zeroth-order Tikhonov regularisation can be conveniently expressed in terms of the SVD of the weighted weight matrix W. This allows us to identify which components in the solution are damped by the regularisation.Writing the SVD of the signalto-noise-ratio-weighted W and signal, as and writing the identity matrix as 1 N = VV T , the solution fλ from ( 12) can be rewritten as The columns of V are the right singular vectors of W, the columns of U are the left singular vectors of W, and Σ T Σ is a square diagonal matrix, with the square of the singular values of W along the diagonal.Therefore ( 16) can be written as where M is the number of points in phase space f ∈ R M×1 .σ i is the i'th singular value, v i is the i'th right singular vector and u i is the i'th left singular vector.The singular values should not be confused with the σ distinguishing positive and negative energies in the phase space.From (17), we see that the reconstruction is a linear combination of the right singular vectors, with the reconstruction coefficients depending on e.g. the singular values and the regularisation parameter λ.When λ ≫ σ i the reconstruction coefficient scales as ∼ σ i /λ 2 and the component of the reconstruction in the v i 'th direction will thus be heavily damped.This fact is summarised in the filter factors To investigate the difficulty in the reconstruction, focusing on the energy slice E = 18 keV, we can project the true solution f true onto an increasing number of right singular vectors v, to see how many are needed in order to obtain a good approximation of the true solution.At the same time, we can plot the filter factors with the optimal parameter λ from the MAP estimate (14), to see how many singular vectors actually contribute significantly to the reconstruction.When λ becomes larger than the singular values, the filter factor gets below 1/2.The projection of f true onto the right singular vectors is given by In figure 9, we show the expansion coefficients α i , the singular values σ i and the filter factors.From the expansion coefficients in figure 9(a), we see that the projection of f true has approximately the same magnitude along first 2844 singular vectors.These singular vectors correspond to the ones with associated singular values much larger than machine precision, as observed in figure 9(b).This is why only the first 3000 coefficients are shown.However, from figure 9(c) it is evident that only the first 760 singular values are larger than the optimal regularisation parameter.For additional terms in (17) the reconstruction coefficients will scale increasingly as ∼ σ i /λ 2 → 0 as σ i → 0. The optimal regularisation parameter suppresses the smaller singular values, and reflects the 10% noise and the interpolation noise.In figure 10 we show the projection f v onto the first 2844 and 760 singular vectors respectively to study the difference.We plot only the σ = 1 case, i.e. positive energies.We notice that the projection onto the first 2844 singular vectors in figure 10  reason for this is found in the geometrical setup of the diagnostics.The trapped orbits spend most of their time around the banana-tips of the orbit, so to obtain as much information as possible about the trapped orbits, it is desirable to have the banana-tips inside the measurement volumes.Note that the reconstruction is reasonably good along the upper topological boundary of the trapped region, since these trapped orbits can be so small that they are completely contained within the measurement volumes.The measurement volumes in the case of COMPASS-Upgrade are located around the midplane, mainly on the low-field side, extending only slightly into the high-field side due to the attenuation of the neutral beam in the plasma.Thus the measurement volumes only contain the banana-tips of the marginally trapped orbits with tips close to the magnetic axis, and very deeply trapped orbits close to the upper topological boundary.The rest of the trapped orbits are observed on the low-field side, where they do not spend much time, and not much information about them is collected.In order to better reconstruct the deeply trapped region of phase space, we therefore propose that it would be beneficial to place measurement volumes off the midplane.This is however a challenge for a FIDA diagnostic, since the measurement volumes are restricted to the NBI beam path through the plasma.
In the next section, we incorporate physics-based prior information into the problem, and thus reconstruct the fastion distribution function while penalising large gradients in constants-of-motion space by employing first-order Tikhonov regularisation.

First-order Tikhonov regularisation in 3D constants-of-motion space
It is known a priori that the fast-ion distribution function must be smooth due to the physics of the interactions, which are governed by the Fokker-Planck equation.For local Coulomb collisions, the Fokker-Planck equation governing the evolution of the distribution function in velocity-space is given by [70] where : denotes the double-dot product.If the gradient in Cartesian velocity-space is large, the collision operator blows up and ensures smoothness.
To penalise large gradients in the solution, we choose L to be the gradient operator approximated by a finite-difference operator L = L 1 , with the subscript denoting that we approximate the first-order finite difference operator.
In the general is given by where M = m 1 m 2 m 3 is the total number of phase-space points, D mi is the m i × m i finite difference matrix and the Kronecker product ⊗ is defined as In the numerical implementation, the first direction is taken to be the µ-direction, the second direction to be the P ϕ -direction, and the third to be the E-direction.For the inverse problem in the three-dimensional constants-of-motion phase space, we impose the smoothness property such that the finite differences in each direction h i are used as the regularisation parameters, From this we can already see the satisfying interpretation that needing less regularising, i.e. smaller λ i 's, corresponds to larger distances h i 's in the i'th direction.For an extended explanation of how (20) and h i are defined, see appendix.In addition to collisional physics, we also utilise that we know that the fast-ion distribution function is non-negative.The first-order Tikhonov regularisation problem is therefore solved with a non-negativity constraint, where and ⊕ is the direct sum of matrices.
As mentioned previously the Jacobian of the transformation to constants-of-motion phase space is singular at the trappedpassing boundary, since the parallel velocity of marginally trapped orbits goes to zero.Thus, absorbing the Jacobian in the distribution function, in order to preserve the units of the weight functions, which is done in this work, introduces steep gradients on a discrete grid.These have to be taken into account when regularising, since we must allow for discontinuities in phase space.This is done by penalising large gradients along the trapped-passing boundary, but not across it.In particular, this is realised by putting the direction across the boundary in the null space of the regularisation matrix L (3D) 1 .At all points in phase space the gradient is calculated, where f m is the mth element of f, and m is an index running over all points of phase space, not to be confused with the mass of the ion.Now define the set of vectors ξm for m = 1, . . ., M, at all points in phase space, as being the unit vector normal to the trapped-passing boundary at all points along the boundary, and the zero vector everywhere else in phase space.The modified gradient allowing for discontinuities across the trappedpassing boundary is then At all points the projection matrix ξm ξT m projects the solution onto the direction across the boundary at the specific phasespace point indexed by m.The matrix 1 3 − ξm ξT m is thus the projection onto the remaining directions.The regularisation matrix allowing steep gradients across the trapped-passing boundary will be referred to as the modified regularisation matrix L(3D) 1 with a hat-notation.For simplicity and transparency it is desirable to express this in terms of the original first-order Tikhonov regularisation matrix L (3D) 1 , as a matrix product The task is now to find the matrix X, which allows the steep gradients.To construct the modified regularisation matrix, one has to consider the bookkeeping of matrix elements, when comparing (20) and (25).The modified regularisation matrix eventually becomes where is the direct sum of matrices and K (M,3) is the commutation matrix with the property of commuting the Kronecker product [71], for every m × n matrix A and v × w matrix B. Also, the commutation matrix is orthogonal with To showcase the difference between L (3D) 1 and L(3D) 1 , we compare the reconstructions of a test distribution in figure 11.We choose the test distribution to have constant densities in the trapped region and in the co-passing region, respectively, with no variation for different energies.Thus only a single plane of constant energy is shown.The test distribution function is reconstructed using a signal with no added noise.We use different grids for the forward problem and the inverse problem, so no inverse crime is committed.It is clearly seen that the discontinuity along the trapped-passing boundary is much better reconstructed using L(3D)
In figure 12 we compare the reconstruction of the σ = 1 part of phase space using plain first-order Tikhonov regularisation with the regularisation matrix L (3D) 1 , with the reconstruction using the modified first-order Tikhonov regularisation, with L(3D)

1
, where large gradients are allowed across the trapped-passing boundary.It is evident that penalising large gradients provides a smooth solution, also across the trapped-passing boundary, as expected.Even though the general features of the true solution (high density in the trapped region) are recovered, the large jumps across the trappedpassing boundary are clearly missing.Using the modified firstorder Tikhonov regularisation results in steeper gradients in the solution across the trapped-passing boundary, as shown in figure 12(b).For both the plain first-order Tikhonov problem and the modified first-order Tikhonov problem, we see that providing the additional physics-based prior information has improved the reconstruction substantially compared to using zeroth-order Tikhonov regularisation, as in the previous section.Specifically, we observe that the peak in the deeply trapped region is well reconstructed.
In figure 13, we compare the same cases, this time in the σ = −1 part of phase space.These reconstructions are almost identical, since there are no expected steep gradients in the valid regions in the σ = −1 phase space.The only differing trend we observe is that using the modified firstorder Tikhonov regularisation seems to concentrate the solution along the trapped-passing boundary, more so than when using the plain first-order Tikhonov regularisation.

Application: optimisation study
In this section, we investigate the optimal location of an additional sixth detector starting from port 3, where a detector is already placed, and rotating the location of the sixth detector 180 • toroidally around the tokamak, always keeping the detector in the midplane.A synthetic signal and reconstruction is calculated for every 10 • , irrespective of whether a port will be available at those positions in the future COMPASS-Upgrade tokamak or not.This is to investigate which angle with respect to the magnetic field will be observed depending on the location of the detector, and subsequently, which locations will yield the best reconstructions.Two additional cases are included at 112.5 • and 135 • as detectors are already placed at these locations.We wish to check that these locations are definitely not the optimal location of an additional detector, since the detector would be redundant.In order to eliminate all variables but the location of the sixth detector, the first five detector locations are fixed, no noise is added on the synthetic signals, and an inverse crime is committed.This means that the forward and the inverse problem are calculated on the same phase-space grid, such that no interpolation noise is introduced.The regularisation technique used is that of zerothorder Tikhonov, to keep the problem and the implementation simple, such that we are sure that the only thing varying from case to case is the location of the sixth detector.Since we commit an inverse crime, we simply calculate the signal via the forward problem using the weight matrix W, and invert the exact same signal using the exact same weight matrix W. This problem is still challenging, since there are 5940 measurements and 19004 unknowns.The problem is thus heavily underdetermined and ill-posed, thus requiring regularisation.The resulting reconstruction errors from using all six fans of sightlines, calculated from (14), are shown in figure 14 with the abscissa denoting how much the location of the sixth detector has been rotated toroidally around the tokamak starting from port 3, and the ordinate showing the reconstruction error.It is clearly evident that when the sixth detector is located on top of

Conclusion
In this paper, we have performed tomographic reconstructions of a 3D fast-ion distribution function in constants-ofmotion phase-space.This is done by inverting synthetic signals based on projected velocities of the fast ions along the diagnostic lines of sight.In this work, the FIDA setup of the future COMPASS-Upgrade was used.The weight functions, based on projected velocities only, are used to calculate synthetic signals in the forward problem, which are then inverted to reconstruct fast-ion distribution functions in COMPASS-Upgrade using both zeroth-order and first-order Tikhonov regularisation.The formulation of Tikhonov regularisation within Bayesian inference allows us to calculate the variance of the reconstruction, providing a measure of uncertainty.From this it is evident that large differences between the reconstruction and the true distribution occur in regions of  phase space with large uncertainties.The challenges regarding the inversion are understood both mathematically, by analysing the SVD of the weight matrix in the case of zerothorder Tikhonov, and physically, by considering the diagnostic setup.Our interpretation of the study is that the large regularisation needed in the inversion dampens singular vectors containing information about the trapped region of phase space containing orbits, which have their banana-tips outside the measurement volumes, which renders us unable to do a good reconstruction in that region using zeroth-order Tikhonov regularisation.This happens because the trapped orbits spend most of their time at the banana-tips.Trapped orbits, along the upper topological boundary of the trapped region, can be so small they are completely contained within the measurement volumes, and can therefore be reconstructed.Since the NBIs for FIDA measurements are usually in the midplane of the tokamak, FIDA diagnostics will often have more difficulties to obtain diagnostic information about trapped particles as compared to passing particles.Adding physics-based prior information, such as collisional physics and non-negativity constraints, enables us to obtain better reconstructions in this region.
We also succeed in modifying the penalty matrix of the first-order Tikhonov regularisation, such that we allow for large gradients across the trapped-passing boundary, which is expected when the phase-space Jacobian is absorbed in the fast-ion distribution function, rather than in the weight functions.
Finally, we applied the developed methods to perform an optimisation study on the positioning of an extra sixth fan of sightlines.This clearly showed that placing an extra fan at the same toroidal position as existing ones does not add extra information.It also identified a specific location, where most new information could be obtained.Specifically this optimal location on the midplane is at a 120 • toroidal angle with respect to port 3. We interpret this location as the detector location complementing the existing detector locations best, regardless of whether one actually does a tomographic reconstruction or not.
This work allows further study in orbit tomography in constants-of-motion phase space using experimental data.Reconstructions of the fast-ion distribution function in this space can probe effects of particle-wave interactions in the plasma, such as fast ions being transported by Alfvén eigenmodes or fast ions driving instabilities.related to the finite differences by (22).To begin with, we consider the continuous case in order to understand which differential operators we are approximating with matrices in the discrete case.Consider some continuous function γ of two variables x and y, such that γ = γ(x, y) and the gradient of γ is ∇γ = ∂γ ∂x , ∂γ ∂y . (A.1) In first-order Tikhonov regularisation we penalise the square of the two-norm of the gradient of the distribution function, and thus we take the two-norm of the gradient of γ, which is a linear operation, and we are thus able to approximate this using matrix-vector products.This shows that we need to be aware of which differential operator we are approximating, and which behaviour we are penalising in the regularisation.
In the discrete case, the partial differential operators in (A.2) are approximated by finite-difference operators, which will be explained in the following.
Before dealing with the 3D case we consider the 2D case in some unspecified 2D space.This could e.g.be the (v ∥ , v ⊥ ) velocity space.Consider a discrete function G ∈ R 2×2 taking values in a 2D space, which we take to be two-by-two in this example.G is written as We take G to be the discrete version of γ.Taking the partial derivatives with respect to x and y in the continuous case will now correspond to calculating the finite differences along columns and rows respectively.In the two-by-two case the finite difference operator with zero Neumann boundary conditions is given by where h i is the finite distance between two neighbouring pixels in the i'th direction in the discrete space, within which we wish to promote smoothing.Thus, if the inverse problem is calculated in the (E, p) coordinates in the two-dimensional case, but the smoothness property is imposed in the (v ∥ , v ⊥ ) coordinate system, then h i is calculated from the Jacobian of the coordinate transformation between the two systems [21].

Figure 1 .
Figure 1.CAD drawing of the future COMPASS-Upgrade tokamak and the diagnostic sightlines used in this work.The sightlines are shown as magenta lines, and the NBI's are shown as orange cylinders.Three horizontal oblique fans of sightlines are located in segment 3, 7 and 9, respectively.Two perpendicular fans of sightlines are located in segment 8. One fan observing from above, which we will denote as port 8u, and one fan observing from below, which will be denoted as port 8l. Figure reprinted with permission from [44].(a) and (b) Reproduced from [44].© 2024 The Author(s).Published by IOP Publishing Ltd on behalf of the IAEA.All rights reserved CC BY 4.0.
Two weight functions are shown in figure2for the two projected velocity intervals u ∈ ±[1.53, 1.59]  • 10 6 m s −1 centred at u = ±1.56• 10 6 m s −1 respectively, corresponding to an energy interval of E u = mu 2 /2 ∈[24.4,26.5]  keV, for the seventh sightline (when counting the measurement volumes from the high-field side towards the low-field side) from port 9. Note that E u varies strongly during a gyration of an ion.It gives a minimum required fast-ion energy to be observable at that point in the spectrum.The measurement volume is centred at (R, z) = (1.06,−0.01) m and the angle of the sightline with respect to the magnetic field is φ = 16.0 • .We can thus see how the sensitivity changes for different fast-ion energies.The weight function for positive u ∈ [1.53, 1.59] • 10 6 m s −1 is shown in figure2(a).Note that only the sensitivity in the σ = 1 part of phase space is shown, since counter-going orbits of these energies do not give rise to such large positive projected velocities in this sightline.This is contrary to figure2(b), which shows the weight function for u ∈ −[1.59, 1.53] • 10 6 m s −1 .

Figure 2 .
Figure 2. Weight functions quantifying the sensitivity of the seventh sightline from port 9 in COMPASS-Upgrade, with φ = 16 • , to projected velocities u =∈ ±[1.53, 1.59] • 10 6 m s −1 .(a) Weight function for u ∈ [1.53, 1.59] • 10 6 m s −1 .Only the σ = 1 part of phase space, which includes co-going and trapped orbits, is shown, since no counter-going orbits give rise to u > 0 in this sightline.(b) Weight function for u ∈ −[1.59, 1.53] • 10 6 m s −1 .Only the σ = −1 part of phase space is shown, since no co-going orbits give rise to u < 0 in this sightline.The sensitivity to trapped orbits is collected in the σ = 1 part of phase space which is not shown here.

[ 1 .
53, 1.59]  • 10 6 m s −1 across different sightlines from the same fan, i.e. observing different measurement volumes along the NBI.Apart from the fact that each sightline observes the plasma at different radial coordinates, they also observe the plasma with different angles φ, with respect to the magnetic field in the centre of the measurement volumes.From figures 3(a) to 3(d) the angle φ drops from 35.5 • to 3.7 • .

Figure 4 .
Figure 4. (a) Example spectrum with and without 10% noise as obtained by the seventh sightline (when counting the measurement volumes from the high-field side towards the low-field side) from port 9 in COMPASS-Upgrade.(b) An example of all spectra incl.noise obtained by the entire FIDA setup of COMPASS-Upgrade.Note that in all spectra, the thermal spectrum u ∈ [−0.7, 0.7] • 10 6 m s −1 , has been removed.

Figure 5 .
Figure 5. Reconstruction error decreasing as additional fans of sightlines are included.The first data point shows the reconstruction error calculated using only data from the upper fan from port 8.The next data point shows the error calculated using data from both the upper fan from port 8 and the fan from port 7, and so on adding extra fans of sightlines for each data point.

Figure 6 .
Figure 6.True and reconstructed fast-ion distribution at σ = 1 for 2 MW NBI injected into COMPASS-Upgrade at tangency radius 0.45 m.
(a) is similar to the true distribution in figure 6(a), whereas the projection onto the first 760 singular vectors in figure 10(b) is similar to the reconstruction in figure 6(b).The interpretation of this is that the large regularisation parameter dampens the singular vectors that contain information about the deeply trapped region, due to them having small singular values.The physical

Figure 8 .
Figure 8.(a) Difference between reconstruction and true distribution fλ − ftrue.(b) Standard deviation of the posterior distribution fλ , i.e. the square root of the diagonal elements of Cf λ .

Figure 10 .
Figure 10.Projection of the true distribution function onto (a) 2844 right singular vectors and (b) 760 singular vectors given by equation (18).

Figure 11 .
Figure 11.Reconstruction of test distribution, showing only a single energy plane E = keV, using the standard regularisation matrix L (3D) 1 in (b) and the modified regularisation matrix L(3D) 1 allowing steep gradients across the trapped-passing boundary in (c).The test distribution is reconstructed using a signal including interpolation noise.

Figure 12 .
Figure 12.Reconstructions of the fast-ion distribution function for σ = 1 using (a) plain first-order Tikhonov regularisation and (b) the modified first-order Tikhonov regularisation allowing steep gradients across the trapped-passing boundary.

Figure 15
shows the angle φ for all 60 sightlines in this case, where angles larger than 90 • are mirrored around 90• , φ ′ = π − φ.It is evident that the extra sixth fan observes a large number of angles not already covered by the original five fans.

Figure 13 .
Figure 13.Reconstructions of the fast-ion distribution function for σ = −1 using (a) plain first-order Tikhonov regularisation and (b) the modified first-order Tikhonov regularisation allowing steep gradients across the trapped-passing boundary.

Figure 14 .
Figure 14.Reconstruction error ε(λ), equation (14), for varying position of an extra sixth fan of sightlines.The reconstructions are calculated using zeroth-order Tikhonov regularisation by committing an inverse crime and without noise on the signal.The dashed lines indicate the toroidal locations of the existing five detectors.

Figure 15 .
Figure 15.Angle φ, mirrored around 90 • , between the sightline and the magnetic field in each measurement volume along the NBI for all fans of sightlines.The angles of the extra fan, denoted by the '+' markers, are shown for the case when the extra fan has been rotated 120 • with respect to port 3.

Table 1 .
Table showing basic data of the diagnostic sightlines.The first row denotes the radial position of the 10 measurement volumes, while the remaining rows denote the angles between the sightlines from each port and the magnetic field in the measurement volumes.