Direct inversion of the Longitudinal ray transform for 2D residual elastic strain fields

We examine the problem of Bragg-edge elastic strain tomography from energy resolved neutron transmission imaging. A new approach is developed for two-dimensional plane-stress and plane-strain systems whereby elastic strain can be reconstructed from its Longitudinal ray transform (LRT) as two parts of a Helmholtz decomposition based on the concept of an Airy stress potential. The solenoidal component of this decomposition is reconstructed using an inversion formula based on a tensor filtered back projection (FBP) algorithm whereas the potential part can be recovered using either Hooke’s law or a finite element model of the elastic system. The technique is demonstrated for two-dimensional plane-stress systems in both simulation, and on real experimental data. We also demonstrate that application of the standard scalar FBP algorithm to the LRT in these systems recovers the trace of the solenoidal component of strain and we provide physical meaning for this quantity in the case of 2D plane-stress and plane-strain systems.


Introduction and context
Elastic strain imaging via energy-resolved neutron transmission measurement (also known as 'Bragg-edge imaging') forms a natural tensor-tomography problem aimed at reconstructing the full triaxial elastic strain field within a physical sample from a set of lower-dimensional scalar images.
The full solution to this tomography problem will have a key impact in a number of areas in science and engineering focused on the study of residual stress in materials.An important topical example includes the development of additive manufacturing techniques for metallic components where residual stresses generated by the thermo-mechanics of deposition are a significant and ever present concern.Tomographic techniques for strain have the potential to provide a unique insight in this area.
The Bragg-edge strain tomography problem has been studied for more than a decade, with various experimental demonstrations on special cases (e.g.axisymmetric systems and in situ applied loads) (e.g.[1,2,3]), and, more recently, solutions for general systems using Bayesian and least-squares techniques constrained by equilibrium (e.g.[4,5]).In this paper we examine this problem from the perspective of developing a direct inversion algorithm.
With reference to Figure 1, strain images of this type refer to projections of the average of elastic strain, ϵ, along straight-line ray paths through a sample Ω of the form where L is the path-length associated with a ray passing through the point x 0 ∈ Ω, travelling in the direction ξ, and, as in the rest of the paper, we use the summation convention for repeated indices.For convenience, strain outside of the boundary of the sample is assigned a value of zero.From many measurements of this form, we wish to reconstruct the original strain field.Bragg-edge strain measurements are naturally related to the Longitudinal Ray Transform (LRT), I, which can be written for suitable with the extension to all of L 2 (S m ; R n ) achieved in the usual way (see below for definitions and notation).Unfortunately, the LRT has a large null space that creates a well-known issue with direct tomographic reconstruction of strain from Bragg-edge imaging and LRT measurements in general [6].For f ∈ L 2 (S m ; R n ), this null space consists of potential fields of the form du for any u that vanishes at infinity.
The structure of this null space, particularly in the case of bounded support, is important for reconstruction.In this context, we will explore the mechanics of linear elastic systems in the context of tensor decompositions and inversion formulas related to the LRT.Through this process, we will provide two direct inversion techniques for the LRT for two-dimensional elastic strain fields that satisfy mechanical equilibrium on bounded domains in the absence of externally applied traction forces.While the detailed context and precise definitions will follow, along the way we will demonstrate; 1.In the case of two-dimensional elastic strain fields, the assumption of zero boundary traction, a condition on stress, implies that the Helmholtz decomposition of strain on a bounded and unbounded domain are equivalent (up to extension by zero).We also demonstrate more generally in R n that the Helmholtz decomposition of a symmetric tensor field on a bounded or unbounded domain are equivalent if and only if the harmonic component of the Helmholtz decomposition on the bounded domain is zero (Lemma 1).
2. The general inversion formula for the solenoidal component of the LRT on m-rank tensor fields over the entirety of R n due to Sharafutdinov [7] is equivalent to that of Louis [8] and Derevtsov et al. [9] for m = n = 2.The latter inversion formula was previously restricted to fields with zero harmonic component on the unit ball -we extend its use to all of R 2 (Lemma 2).
3. In R 2 , the Helmholtz decomposition of any elastic strain field can be specified directly through the concept of an Airy stress potential and Hooke's law (Proposition 1).
4. The application of standard scalar filtered back projection to the LRT recovers the trace of the solenoidal component, which in the case of elastic strain in R 2 is proportional to the hydrostatic component of stress (see section 7).
We begin by introducing the notation used throughout the paper.

Notation and definitions
First, Ω will be an open subset of R n with Lipschitz boundary possibly equal to R n in the following definitions, and we will write S n−1 for the unit sphere in R n .Given a vector v ∈ R 2 , we write v ⊥ for the anti-clockwise rotation of v by 90 degrees.
; Ω) will be the space of smooth m-rank symmetric tensor fields on Ω with continuous derivatives of all orders and C ∞ c (S m ; Ω) the subspace of C ∞ (S m ; Ω) comprising fields with compact support in Ω.
We use the following differential operators: will be the symmetric derivative defined in [7].This coincides with the gradient when m = 0 and for , where ⊗ refers to dyadic product and (•) T refers to the transpose operation; where e ijk is the usual Levi-Civita permutation symbol.Equivalently This is the contraction of the gradient of a tensor field and for the general formula see [7].For u ∈ C ∞ (S 1 ; R n ), Div(u) is the standard divergence of u; Div ⊥ -The perpendicular divergence which is the formal adjoint of −d ⊥ and maps C ∞ (S m+1 ; R 2 ) → C ∞ (S m ; R 2 ).This is the same as the operator δ ⊥ in [9].
We additionally say that a tensor field is divergence-free or solenoidal if its divergence is zero.The differential operators are initially defined on smooth tensor fields, but can be extended to fields with distributional coefficients.
For function spaces we use: L 2 (S m ; Ω) -The space of square-integrable m-rank symmetric tensor fields on Ω with norm ∥u∥ L 2 (S m ;Ω) .
H k (S m ; Ω) -The Sobolev space of square-integrable m-rank symmetric tensor fields on Ω whose weak derivatives up to order k are also squareintegrable.
Ḣ1 0 (S m ; Ω) -The homogeneous Sobolev space which is the closure of The homogeneous Sobolev spaces are equivalent to the standard Sobolev spaces of fields with trace zero when Ω is bounded, but different for unbounded Ω.We will mostly be concerned with tensors of rank either m = 1 or 2 and use the standard notations f : g for contraction of 2-rank tensors and f • g for multiplication of a 2-rank tensor with a 1-rank tensor, or the dot product of 1-rank tensors.
We now return to the topic and begin with a review of Helmholtz decomposition and inversion of the LRT, both in general, and in the context of elastic strain in R 2 .

Helmholtz decompositions and LRT inversion formulas
As per [7] and others, the null space of the LRT forms part of the orthogonal Helmholtz decomposition in R n of symmetric tensor fields of the form where Here the differential operators are understood to act in the sense of distributions on R n , and for given f the decomposition (3) is unique.
Using the fundamental theorem of calculus, it is easy to check that If = I s f , and so at best we can hope to recover s f from the LRT of f .In fact, such recovery is possible as demonstrated by Sharafutdinov [7] (see (6) below).However, an interesting practical problem exists when applying this to real systems; even if f is compactly supported, u and s f in (3) may have unbounded support, and for practical computation it is usually necessary to consider a bounded domain.In this light, let us introduce solenoidal decomposition on a bounded domain.
Let Ω ⊂ R n be a bounded domain with Lipschitz boundary and outward surface normal n on ∂Ω.Similar to (3), there is a unique decomposition of f ∈ L 2 (S m ; R n ), m ≥ 1, restricted to this set of the form (see [10] for the case of vector fields) where u Ω ∈ Ḣ1 0 (S m−1 ; Ω), h Ω ∈ H 1 (S m−1 ; Ω), known as the 'harmonic part', satisfies Div(dh Ω ) = 0 on Ω, and s f Ω ∈ L 2 (S m ; Ω) satisfies the weak equation It is clear from (5) that s f Ω extended by zero to R n is divergence-free.A key point for our result is that, for fields where the boundary trace makes sense, this extension by zero is only divergence-free when the boundary condition ( s f Ω ) i 1 ... im n im | ∂Ω = 0 holds.For an in depth discussion of weak formulation of the Helmholtz decomposition in the case of L 2 vector fields, see [10].
To relate reconstruction formulae for the LRT on R n to formulae on a bounded set, we must consider the relationship between the decompositions (3) and (4).Indeed, when the harmonic part vanishes in (4), the solenoidal decomposition on the bounded set Ω is related to the one on R n as in the following lemma.
Lemma 1. Suppose that Ω contains the support of f .If h Ω in (4) is zero, then s f and u in (3) are equal to the extension by zero of s f Ω and u Ω to R n .Conversely, if s f and u in (3) are supported in Ω, then h Ω = 0.
Proof.Assume that decomposition (4) holds with h Ω = 0 and extend u Ω and s f Ω to v and g on R n by setting them equal to zero outside of Ω.By (5), g is then divergence-free on R n .And since u Ω ∈ Ḣ1 0 (S m−1 ; Ω) then v ∈ Ḣ1 0 (S m−1 ; R n ) and dv is du Ω extended by zero to R n .By uniqueness of the decomposition in (3), u = v and s f = g.
Conversely, suppose that s f and u in (3) are supported in Ω and define ũΩ and s fΩ by restricting their domain to Ω.Then, since u ∈ Ḣ1 0 (S 1 ; R n ) with support contained in Ω, its restriction ũΩ is in Ḣ1 0 (S 1 ; Ω).Additionally, we can see that (5) holds for s fΩ because the same must hold for s f on R n for any φ ∈ H 1 (S m−1 ; R n ).By uniqueness of the decomposition we see that (4) holds with u Ω = ũΩ , s f Ω = s fΩ and h Ω = 0 on Ω, as claimed.Now let us turn to inversion of the LRT.Various inversion formulas exist that can uniquely recover s f from If (e.g.[7,9,8]).Sharafutdinov [7] provides the general result for f ∈ L 2 (S m ; R n ) as where c k are specified scalar coefficients, powers of the Laplacian (−∆) 1/2 and (−∆) −1 are defined via the Fourier transform, the operators i and j respectively refer to product and contraction with the Kronecker tensor, and µ m is the formal adjoint of I when the measure on S n−1 is normalised to one.
In practical terms, µ m is related to the adjoint of the X-ray transform 1 (i.e.scalar back-projection), R * , acting component-wise with back-projections weighted by the diadic product of ξ with itself m-times; Note that the constant factor is present because of the normalisation of the measure on S n−1 in [7].
For 2D elastic strain ϵ ∈ L 2 (S 2 ; R 2 ), (6) simplifies to where c 0 = 3/4, c 1 = −1/4, tr is the trace operator, I is the 2-rank identity and In comparison, Derevtsov and Svetov [9] and Louis [8] consider recovery when Ω is the unit ball in R 2 , implicitly assuming also that the harmonic part of the field is equal to zero so that s ϵ = s ϵ Ω by Lemma 1.
In this context, [8] provides a much simpler inversion formula of the form while Derevtsov and Svetov [9] provide the same formula (9) but, due to a typographical error, multiplied by a factor of 2 on the right side.We now show in Lemma 2 that (8) and ( 9) are indeed equivalent.This extends the inversion results of [8,9] from the unit ball to R 2 , and handles the case of non-vanishing harmonic part, which was not considered in [8,9].Lemma 2. For any ϵ ∈ L 2 (S 2 ; R 2 ), the right hand sides of (8) and (9) are equal and hence (9) can be used on all of R 2 regardless of any harmonic component.
Proof.Taking the component-wise Fourier transform with spatial frequency vector κ, (8) can be written 1 Equivalent to the Radon transform in 2D.
where g = I * Iϵ.Since s ϵ is solenoidal s εκ = 0 and we can write s ε = ακ ⊥ (κ ⊥ ) T for some α ∈ L 2 (R 2 ).Hence (10) becomes Multiplying by κ ⊥ (κ ⊥ ) T and rearranging; Now g is also solenoidal and hence can also be written ĝ In the spatial domain this implies over all of R 2 : which is identical to ( 9) but on all of R 2 .
Given Lemma 2, we use only ( 9) which provides a component-wise approach to reconstruction of the solenoidal component of strain in R 2 of the form where Λ is the Ram-Lak filter (or similar) used in standard scalar Filtered Back Projection (FBP).
Because of Lemma 2, we know that this inversion formula recovers the solenoidal part on all of R 2 with potentially unbounded support regardless of the finite nature of the sample.By Lemma 1, the solenoidal component of ϵ will have support contained in a bounded domain only if its harmonic part vanishes, and so it is important to know when this will occur in the context of strain.
Before we address this, we first provide a brief review of the mechanics of stress and strain on the plane in the context of this work.

Elasticity theory and residual stress
Consider a sample consisting of an elastic body in R 3 represented by the bounded domain Ω with outward surface normal n.Within Ω we can decompose the total strain at each point, ϵ T , into an elastic component, ϵ and an 'eigenstrain', ϵ * (e.g.permanent strain introduced by plasticity, phase change, thermal expansion, etc.) [11,12] The elastic component of strain is related to stress, σ, through Hooke's law, which in its most general form, can be written in terms of a 4-rank stiffness tensor; σ ij = C ijkl ϵ kl .In the isotropic case with Young's modulus E and Poisson's ratio ν Governing equations can be assembled for this system on the basis of equilibrium, compatibility of strain and boundary conditions.In the absence of body forces (gravity, magnetism, etc.) mechanical equilibrium holds that The total strain physically originates as the symmetric gradient of a displacement field (i.e. is potential) and can be expressed as ϵ T = du for some u, where, in general, u ̸ = 0 on ∂Ω.This condition is known as strain 'compatibility' which for a simply connected domain can be expressed as a vanishing Saint-Venant operator 2 , W (ϵ T ) = 0, or The final ingredient is to specify boundary conditions experienced by the sample.These can vary, but in the case of 'residual stress' problems, the surface of the sample is typically free of any traction Equations ( 14), ( 15) and ( 16) together form an elliptic boundary value problem for ϵ based on a known eigen-strain ϵ * .While σ and ϵ are inherently three-dimensional in nature, there are two typical limiting assumptions on the plane that have practical utility [14]: Plane-strain is a limiting case for thick prismatic samples, while plane-stress relates to thin two-dimensional samples where 'thick' and 'thin' refer to dimensions in the x 3 direction.The above analysis applies directly to both cases where Ω ∈ R 2 with the exception that, in the plane-stress case, the isotropic elasticity tensor becomes 2 The Saint-Venant operator is defined by In R 3 , this simplifies to six unique components specified by the 2-rank symmetric incom- kl where e ijk is the Levi-Civita permutation symbol.In a simply connected domain in R n , W (f ) = 0 if and only if f = du for some u.On a multiply connected domain with k holes, n(n + 1)k/2 additional integral constraints are required along with W (f ) = 0 to imply f = du (see [13, Proposition 2.8]).

Problem statement
We are now in a position to state precisely the inverse problem we seek to solve in this work.
The rest of the paper is focused on developing a solution to this problem and demonstrating its numerical implementation.

Helmholtz decomposition of strain in R 2
We begin by connecting the stress σ and strain ϵ initially defined only on the bounded set Ω to the solenoidal decomposition (3) on all of R 2 .Given that the stress σ satisfies ( 14) in the classical sense (i.e. is twice differentiable) on Ω and satisfies the traction-free boundary condition (16), in fact σ extended as zero outside of Ω is divergence free in the distributional sense and is therefore its own solenoidal part with no potential part if decomposed according to (3).Our goal in this section is to use this fact, together with (20) or (19) to find the solenoidal decomposition of ϵ.
This can be achieved through the concept of an Airy stress function.In both the plane-stress and plane-strain cases, it is possible to write σ in terms of a scalar Airy stress potential, ψ ∈ H 2 (Ω) in such a way that it automatically satisfies equilibrium: When combined with Hooke's law (i.e. ( 13) or ( 17)), it follows that strain can also be written in terms of this same potential as for plane-strain conditions, or in the case of plane-stress.
Both (20) and ( 19) already appear to be in the form of Helmholtz decompositions, however the issue is that the Airy stress potential appearing in (18) may not satisfy equilibrium in a distributional sense when extended as zero to R 2 .The next lemma shows that when the traction-free boundary condition ( 16) is satisfied, in fact there is an Airy stress potential which extends as zero.
Proposition 1. Suppose that σ ∈ L 2 (S 2 ; R 2 ) has support contained in a bounded and simply connected set Ω and satisfies (14) in the distributional sense on R 2 .Then there exists unique for a constant M > 0 which depends on Ω but not σ.
Proof.First consider the case when σ ∈ C ∞ c (S 2 ; R 2 ) satisfies ( 14) and has support contained in Ω which is itself inside an open ball B R of radius R centred at the origin.The two columns of σ, σ i1 and σ i2 , are divergence free vector fields on R 2 and so the path integrals of e ik3 σ ij dx k between any two points are independent of path due to Green's theorem.For x 0 ∈ ∂B R and any x ∈ R 2 , we define new functions via the path integrals in which the path is left unspecified.Defining the vector field ϕ = (ϕ 1 , ϕ 2 ) it follows, due to path independence and the fundamental theorem of calculus, that Additionally, since Ω is simply connected, for any x ∈ R 2 \ Ω we can choose a path from x 0 to x outside of Ω and by its path integral definition (23), we have ϕ(x) = 0. Thus, we conclude that ϕ is also supported in Ω.
Next, from (24) we obtain This implies as before that line integrals of e ik3 ϕ i dx k between two points are independent of path, and we define Also as before, this implies that ψ is supported in Ω and Putting together the previous construction and using path independence we see that ψ is directly related to σ by the formula Since the support σ is bounded we can restrict the area of integration in the previous integrals to bounded rectangles, and then use the Cauchy-Schwartz inequality to prove (22) where the constant M depends only on the size of Ω.
We have now proved the proposition for the case when σ is smooth.For σ ∈ L 2 (S 2 ; R 2 ) we approximate by a sequence σ j ∈ C ∞ c (S 2 ; R 2 ) of divergence free fields such that σ j → σ in L 2 (S 2 ; R 2 ) and each σ j is supported within a domain with its boundary within a distance of 2 −j from ∂Ω.By (22) the corresponding potentials ψ j also converge in H 2 (R 2 ) to a function ψ and by continuity of the derivatives from H 2 to L 2 we see that (21) also holds.The supports of the potentials will also shrink to Ω and so we see that the support of ψ is contained in Ω.
Finally, note that from (21) the potential ψ ∈ H 2 (R 2 ) satisfies the biharmonic equation This equation has a unique solution in H 2 (R 2 ) and so the proof is complete.
From Lemma 1, we can conclude the following: If a two dimensional residual elastic strain field on the simply connected bounded domain Ω exists in the absence of boundary traction, its extension by zero to all of R 2 has a unique Helmholtz decomposition of the form where s ϵ and dω are compactly supported within Ω.Note that we only assume that the support of ϵ is contained within the simply connected set Ω, not that the support of ϵ is itself simply connected.By uniqueness and comparison to (19) and (20), this decomposition can be written in terms of the Airy stress potential as in the case of plane-strain, or for plane-stress.Note that in each case s ϵ is proportional to σ.
From this decomposition and the inversion formula for s ϵ we now seek to recover the full elastic strain tensor over a sample.Before we approach this task, we provide a brief comment on recent experimental work in this area.

Isotropic strain and scalar Filtered Back Projection
Some recent work in Bragg-edge strain tomography has approached this problem through an assumption that strain is isotropic at all points within the sample; i.e. ϵ = ε I for some scalar mean strain ε.This assumption is plainly false in almost all cases; the only hydrostatic stress field (and hence strain field) that satisfies equilibrium is constant for all x.However, the assumption does allow for a direct means of reconstruction by standard scalar FBP since Iϵ = Rε for this case.
For example, in Busi et al [15] the authors perform a slice-by-slice FBP to recover an assumed isotropic strain within an additively manufactured stainless steel cube from a set of 19 Bragg-edge strain images.Similarly, Zhu et al [16] recover an assumed scalar isotropic strain in a laser welded steel sample using a similar technique.
Clearly the assumption of isotropic strain was invalid in both cases, however the question remains: What has been recovered?How does the scalar FBP of the LRT relate to the strain field within the sample?
To answer this question, we examine the trace of the solenoidal component of elastic strain in (11) to obtain the following (note that |ξ| = 1); Hence the recovered scalar field stemming from an isotropic assumption is precisely the trace of the (in-plane) solenoidal component, and in general there are no further conclusions that can be made.However, if the strain field is inherently two-dimensional, we can extend this result by considering stress in terms of the Airy potential.As before, under plane-stress or plane-strain conditions, s ϵ can be interpreted through the natural Helmholtz decompositions ( 27) and (29).From this perspective, it follows that for plane-strain and for plane-stress

Recovery of ϵ from s ϵ
We now turn our attention to the problem of recovering ϵ from s ϵ using the constraints provided by elasticity theory.To this end, we present three approaches to the solution of Problem 1.

Recovery of ϵ from compatibility
Applying the Saint-Venant operator to (25) implies W (ϵ) = −W (ϵ * ) = W ( s ϵ) and we can replace the compatibility relation (15) to form a boundary value problem for ϵ; Under two-dimensional plane-stress or plane-strain conditions we can satisfy equilibrium via (19) or (20), and the compatibility condition becomes a non-homogeneous bi-harmonic equation subject to the boundary condition Potentially this provides a direct approach to recover ψ and hence ϵ through numerical solution.However, it should be recognised that computing the right hand side of (33) involves taking second order numerical derivatives.
In the presence of experimental uncertainty, this is likely to be a very unstable process.

Recovery of the potential component
An alternate approach involves the recovery of the potential part of ϵ using equilibrium.From ( 25) and (32), the equilibrium of the system implies This is in the form of a standard structural elasticity problem for ω as a displacement field resulting from a distributed body force and trivial Dirichlet boundary condition.For 2D plane-stress conditions In contrast to the previous approach, calculation of b only involves computing first derivatives, and hence is potentially a much more stable process.

Recovery of ϵ from Hooke's law
By far the most direct means for recovering ϵ from s ϵ is through Hooke's law.From Proposition 1, we have that σ = (d ⊥ ) 2 ψ for some Airy stress function ψ, and therefore using ( 27) and (29) together with (19) and ( 20), we can write for plane-strain, or for plane-stress conditions.9. Numerical demonstration: Simulated data

Strain fields
Numerical demonstrations of the above process were performed on three synthetic two-dimensional plane-stress strain fields.The first of these fields was generated over the unit disk from an Airy stress potential of the form with α = 15, and elastic properties E = 1 and ν = 0.34.The three independent components of this strain field are shown in Figure 2a.
The second and third fields corresponded to finite element simulations of physical samples that were the focus of prior experimental work [4].All relevant details can be found in the reference, however a brief description of each sample is as follows; 1. Crushed Ring: A sample formed by plastically deforming an initially stress-free steel ring along its diameter.The geometry of the sample and applied deformation is shown in Figure 3a.The residual strain field in this sample originates from a distributed eigen-strain related to plastic deformation (see Figure 4a) 2. Offset Ring-and-Plug: A cylindrical steel sample constructed by shrinkfitting an oversize cylindrical 'plug' into an undersize hole that is offset from the centreline (see Figure 3b).The strain field within this sample originates from the interference between the offset ring and the plug (see Figure 5a).In the context of ( 12), the interference imposes a discrete eigen-strain with localised support on the interface.
Both samples were 14mm thick and were simulated as steel with E = 209GPa, ν = 0.34 and a yield stress of 650MPa.The finite element model for the first sample required a non-linear solve based on an elasto-plastic material model, while the second sample was modelled using linear-elasticity.Both models were built and solved in the software package PTC/Creo.
All three strain fields were represented as three scalar components mapped to regular two-dimensional grids.The size and resolutions of these grids were as follows: Airy -400 × 400, spacing 0.006, Crushed Ring -500 × 500, spacing 48µm, Ring and Plug -521 × 521, spacing 50µm.In each case, all three strain components were extended by zero outside the sample boundaries.What follows is a demonstration of the reconstruction of these fields from synthetic LRT data.

Procedure
The demonstrations were was carried out with the help of the Matlab 'radon' and 'iradon' functions.In this context, the implementation was as defined in the following process: Matlab functions.These transformed components are then multiplied by appropriate κ-space filters corresponding to ∂/∂x 1 and ∂/∂x 2 before transforming back to the real domain using 'fftshift' and 'ifft2' 5. From these derivatives, calculate the two components of the vector b using (38) and (39).
6. Using the Matlab PDE solver, calculate a finite element solution for the displacement field ω satisfying (36) and (37) subject to the calculated vector field b.
7. Calculate a second reconstruction for ϵ as the sum ϵ = s ϵ+dω, where dω is computed from the shape functions within the finite element solution.The target element size for the finite element model in step 6 was set to be 0.5% of the maximum sample dimensions.This was conservatively chosen through a standard mesh-independence investigation.

Results
In all three cases the reconstructions based on Hooke's law and the finite element recovery of the potential component were visually indistinguishable from each other.However, the reconstruction based on Hooke's law was slightly more accurate in terms of a root-mean-square error.It was also interesting to note that, in each case, the reconstructed solenoidal component was approximately zero outside the sample boundary (as expected from Lemma 1).This is examined further in Section 9.5 below.The difference between the reconstructions and the original field was small; typically around 1-5% of the maximum value of the original components.However, it was observed that this did not significantly decrease along with the number of projections.The source of this persistent discrepancy was discretisation error related to minor deviations from the equilibrium relation introduced by various interpolations onto the regular grid.This is examined further in the following section.
Figure 6 shows the computed Saint-Venant incompatibility of the recon-structed solenoid compared to the original for all three fields.These images were calculated using a similar transform-filter-transform approach in the Fourier domain.
The Airy stress field shows incompatibility distributed over the sample domain, whereas the other two samples show more localised support.In the case of the crushed-ring, this is likely to have originated from localised plastic shear within the elasto-plastic finite element model, while the offset ring-and-plug indicates a clear dipole around the circumference of the plug corresponding to the interference.
As expected, the incompatibility of the reconstructed solenoidal components are identical to that of the original fields within a small amount of numerical noise.

Reconstruction in the presence of measurement uncertainty
A further set of simulations was carried out in order to examine the behaviour of reconstructions in the presence of Gaussian noise.In this respect both approaches were found to be quite stable and converged to the original field with an increasing number of projections (notwithstanding the discretisation error identified earlier).
Although not strictly necessary, slight improvement was found by limiting the order of terms in the numerical derivatives used to compute b.This was achieved by cutting-off the κ-space filters for frequencies above a certain threshold.A cut-off frequency equal to 0.7 times the maximum magnitude provided a good compromise between noise and fidelity.
For the Airy stress field, Figure 7 shows the convergence of the reconstructed fields along with the number of projections in the presence of Gaussian random noise with a standard deviation of 10% of the maximum LRT value.Results from three systems are shown corresponding to different spatial resolutions (i.e.grid size).In each case, convergence of the relative error to zero is observed to occur at O(n −1/2 ) until the lower limit corresponding to the discretisation error is reached.
Generally speaking, the reconstruction based on Hooke's law had a lower persistent error and the size of the persistent error was observed to be directly related to the resolution of the grid.
It should be noted that, in the presence of noise the calculation of the Saint-Venant operator was found to be inherently unstable regardless of any reasonable cut-off frequency used in the relevant filters.The overall error in the reconstruction of the Airy stress field in the presence of 10% Gaussian measurement noise as a function of the number of projections.The relative error is computed as the root-mean-square of the residual divided by the root-mean-square of the original strain field summed over all components.Dotted lines show the minimum error possible for the given mesh density (calculated using 50,000 projections with no added noise).

Boundary traction and compact support
In order to examine the effect of the boundary conditions, a further set of simulations were carried out on the strain field specified in Appendix A of Gregg et al [17] with e 0 = R = 1 (see Figure 8a).This is an axi-symmetric 'plane-stress' strain field on the unit disk originating from the hydrostatic eigen-strain ϵ * rr = ϵ * θθ = (1 − r) 2 , and subject to a zero traction boundary condition (i.e.σ rr (1) = 0).In polar coordinates it has the form A simulated reconstruction based on 1000 equally spaced LRT projections from a 400 × 400 Cartesian grid is shown in Figure 8b.As expected, the reconstructed strain matches the original field accurately and the support of the reconstruction is contained within the boundary of the sample.Outside of the boundary, the reconstructed solenoidal component was around three orders-of-magnitude smaller than the original field.Figure 8c shows the residual between the LRT of the original field and the reconstruction.
Figure 8d shows the same field with the addition of a constant hydrostatic strain of magnitude ε = 0.2.Like the original field, this altered version sat-isfies equilibrium at all points within the sample, however it clearly violates the traction-free boundary condition since |σ • n| = 0.2 on ∂Ω.
An attempted reconstruction of this field based on the same process is shown in Figure 8e.A visual inspection of the result clearly indicates the reconstruction has failed to reproduce the original field.
It is also interesting to note that the reconstructed field is far from zero outside the boundary of the sample.This observation, together with Lemma 1 suggests that the apparent support of s ϵ reconstructed from data gives a reliable indicator of the existence of a harmonic potential component, and hence the appropriateness of the traction-free assumption for a given experimental system.
It is also clear that the LRT of the reconstructed solenoid does not match that of the original field.Figure 8f shows the difference between these two sinograms computed with s ϵ masked to zero outside the boundary.The residual is of a significant magnitude and appears to correspond directly to the added hydrostatic/harmonic component.This poses an interesting question: Given the harmonic component is compatible, can it be recovered through reconstruction of a non-zero boundary condition similar to the process carried out by Hendriks et al [1]?This question will form the focus of future work in this area.

Numerical demonstration: Experimental data
As a final demonstration, the reconstruction approach was applied to experimental data measured from physical samples using the RADEN energy resolved imaging instrument within the Materials and Life Sciences institute at the J-PARC spallation neutron source in Japan [18].All relevant details of this experiment are described in Gregg et al [4].The outcome of this experiment was measured strain-sinograms from the crushed-ring and offset ring-and-plug samples corresponding to a set of 50 golden-angle projections.As per (1), these measurements correspond to average strain along ray-paths, which require multiplication by appropriate values of L to compute the LRT (see Figure 9a and 9b).
Figure 9d and 9f show the results of the reconstruction based on Hooke's law compared to traditional neutron diffraction based strain measurements from the KOWARI engineering diffractometer at the Australian Centre for Neutron Scattering within the Australian Nuclear Science and Technology Organisation [20].This reference data (Figure 9c and 9e) is in the form of  interpolated/inferred fields computed from scattered measurements using a technique that guarantees equilibrium is satisfied at each point [19].
Overall the reconstruction has performed well in terms of overall mag-nitude and distribution within the limits of resolution.In particular, the reconstructions show remarkable similarity to that of previous work from the same data by Gregg et al [4] based on constrained least squares optimisation of Fourier basis functions.

Conclusion
A direct link has been established between the concept of Airy stress potentials in two-dimensional elastic systems and the standard Helmholtz decomposition at the heart of the LRT and its null space.For homogeneous, isotropic materials under plane-strain or plane-stress conditions, when the stress field satisfies equilibrium and has zero boundary traction, then the Helmholtz decomposition of the strain field can be written in terms of an Airy stress potential allowing for identification of the solenoidal and potential parts, which will have compact support.
Through this lens, direct approaches for the reconstruction of two-dimensional elastic strain fields from LRT data have been developed and demonstrated.We show that a tensorial version of standard FBP recovers the solenoidal (divergence free) component of the strain field, which can then be used to determine the original field through the application of Hooke's law or a process involving the numerical solution of a standard elasticity problem.In simulation, both approaches were found to be robust to measurement noise.Both approaches also performed well on real experimental data.
From this perspective, it was also possible to identify the result of standard scalar FBP when applied to LRT measurement as the trace of the solenoidal component.In some situations (e.g.plane-stress or plane-strain) this can be related to the trace of the stress tensor, however in general, more information is required to bring meaning to such a reconstruction in a three-dimensional system.

Acknowledgements
This work is supported by the Australian Research Council through a Discovery Project Grant (DP170102324).Access to the RADEN and KOWARI instruments was made possible through the respective user access programs of J-PARC and ANSTO (J-PARC Long Term Proposal 2017L0101 and ANSTO Program Proposal PP6050).

Figure 1 :
Figure 1: Geometry of the Longitudinal Ray Transform and Bragg-edge strain measurements.

Figure 2 :
Figure 2: A reconstruction of a synthetic strain field computed from an Airy stress field.(a) The original strain field.(b) A reconstruction of the solenoidal component of this field from a simulated LRT consisting of 200 equally spaced projections over 360 • .(c) The recovered potential component from elastic finite element modelling.(d) The reconstructed strain field formed by the sum of the solenoidal and potential components.

Figure 3 :
Figure 3: Two samples representing strain fields used to perform numerical demonstrations of the reconstruction algorithm.(a) A crushed steel ring containing a distributed eigenstrain field.(b) An offset ring and plug system containing a discrete eigen-strain field generated through mechanical interference.

Figure 4 :
Figure 4: A reconstruction of a synthetic strain field computed from an elasto-plastic finite element model of the crushed ring.(a) The original strain field.(b) A reconstructed of the solenoidal component of this field from a simulated LRT consisting of 200 equally spaced projections over 360 • .(c) The recovered potential component from elastic finite element modelling.(d) The reconstructed strain field formed by the sum of the solenoidal and potential components.

Figure 5 :
Figure 5: A reconstruction of a synthetic strain field computed from an linear-elastic finite element model of the offset ring and plug system.(a) The original strain field.(b) A reconstructed solenoidal component of this field from a simulated LRT consisting of 200 equally spaced projections over 360 • .(c) The recovered potential component from elastic finite element modelling.(d) The reconstructed strain field formed by the sum of the solenoidal and potential components.

Figures 2 ,
4 and 5 show the results of this process based on simulated LRT data from 200 equally spaced angular projections over 360 • .Each figure shows the original strain field together with the reconstructed solenoidal component, the recovered potential component, and the final reconstruction based on the sum of the two.

Figure 6 :
Figure 6: The Saint-Venant operator as applied to the reconstructed solenoidal components compared to the same for the original strain fields.

Figure 7 :
Figure7: The overall error in the reconstruction of the Airy stress field in the presence of 10% Gaussian measurement noise as a function of the number of projections.The relative error is computed as the root-mean-square of the residual divided by the root-mean-square of the original strain field summed over all components.Dotted lines show the minimum error possible for the given mesh density (calculated using 50,000 projections with no added noise).

Figure 8 :
Figure 8: Demonstration of the effect of the no-traction boundary condition on reconstruction.(a) An axisymmetric residual 'plane-stress' strain field that satisfies the no-traction boundary condition (see [17]), along with (b) its reconstruction, and (c) the residual LRT between these fields.(d) The same field with an additional hydrostatic component that violates the no-traction condition, along with (e) a failed reconstruction and (f) the nonzero residual.

Figure 9 :
Figure 9: Reconstruction of residual strain fields from real data.(a) and (b) MeasuredLRT data from the crushed-ring and offset ring-and-plug samples using Bragg-edge strain imaging on the RADEN energy-resolved neutron imaging instrument[4]  (c) and (e) Reference measurements from each sample taken using traditional neutron diffraction based strain measurement techniques on the KOWARI engineering diffractometer (see[19]).(d) and (f) reconstructed strain fields formed by the sum of the reconstructed solenoidal and recovered potential components.