Mössbauer Imaging

In a Mössbauer experiment, if a spatially-extended absorbing sample is rotated relative to a moving γ-ray source, lines of constant γ-ray Doppler shift are generated through the absorber parallel to the motion of the source. As a result, resonant absorption takes place along a series of parallel lines cutting through the absorber, where a particular line is determined by the velocity of the source. The result is a series of measurements of line integrals of the absorption coefficient through the absorber. An image or spatial map of the absorption coefficient distribution may then be reconstructed using tomographic image-reconstruction algorithms. Moreover, when measurements are recorded both as a function of the source velocity and the absorber rotational velocity, spectral information may also be recovered as a function of position. Spatial resolution is proportional to the rate of rotation of the absorber, but is ultimately signal-to-noise limited.


Mossbauer Imaging
In a Miissbauer experiment, if a spa-tiaJly-extended absorbing sample is rotated relative to a moving y-ray source, lines of constant y-ray Doppler shift are generated through the absorber parallel to the motion of the source. As a result, resonant absorption takes place along a series of parallel lines cutting through the absorber, where a particular line is determined by the velocity of the source. The result is a series of measurements of line integrals of the absorption coefficient through the absorber. An image or spatial map of the absorption coefficient distribution may then be reconstructed using tomographic imagereconstruction algorithms. Moreover, September-October 1987 when measurements are recorded both as a function of the source velocity and the absorber rotational velocity, spectral information may also be recovered as a function of position. Spatial resolution is proportional to the rate of rotation of the absorber, but is ultimately signal-tonoise limited.
Accepted: April 30, 1987 The Mossbauer effect is an important spectroscopic tool in materials science [1,2j1. In a typical Mossbauer experiment, 'Y-rays emitted by a radioactive source are allowed to impinge upon a sample of absorbing material; the transmitted or resonantly-scattered 'Y-rays are then subsequently detected ( fig. I). Thus, in a transmission experiment, a detector counts the number of 'Y-rays passing through the absorber, and in a scattering experiment the detector counts the 'Y-rays re-emitted after resonant absorption. The Mossbauer effect, or nuclear 'Y-ray fluorescence, is made possible by the recoilless emission and absorption of the 'Y-rays in nuclei embedded in a solid lattice. As a result, the resonant Iinewidths, relative to the 'Y-ray energy, are characteristically very narrow. Moreover, the resonant absorption can be "tuned" merely by moving the source at a small velocity relative to the absorber, which imparts some kinetic energy, or Doppler shift, to the 'Y-ray. Changes in source-absorber speeds as small as a few millimeters per second are often sufficient to destroy the resonance. In the conventional Mossbauer experiment, an absorption spectrum of the material under study is generated by moving the source relative to the absorber and counting the transmitted 'Y-rays (or resonantly-scattered 'Y-rays in a scattering experiment) as a function of the relative source-absorber velocity. In this paper, we propose the idea of Mossbauer imaging. An image is by definition some quantity or parameter of interest displayed as a function of position, i.e., a picture of the parameter as it is distributed in space. In Mossbauer imaging, the simplest example of such a parameter is the Mossbauer absorption coefficient, and the aim would be to reconstruct and display its spatial distribution. A conventional Mossbauer experiment measures only a bulk average of the resonant absorption coefficient over the absorbing specimen. Thus, spatial inhomogeneities within an extended absorber, due either to variations in the absorption cross section of a pure material or to the mixture of different nuclei within a composite material, are inevitably averaged out in the bulk measurement process. A Mossbauer imaging experiment, however, would permit the reconstruction of a two-dimensional map of the Mossbauer absorption coefficient, that is, a two-dimensional picture of the strength of the resonant absorption as a function of position. In another, somewhat more complex, version of an imaging experiment, spectroscopic information as a function of position should also be recoverable, rather than merely differences in the absorption coefficient as a function of position. As a consequence, in the latter version, true M6ssbauer spectroscopy can be performed in an imaging mode, so that heterogeneous or composite samples may be investigated.
The idea of M6ssbauer imaging was inspired by the success of nuclear-magnetic-resonance (NMR) imaging, since NMR and the M6ssbauer effect share some fundamental characteristics, both being nuclear resonance phenomena. NMR imaging has recently found notable success in diagnostic medicine [3]. While there are no foreseeable applications of M6ssbauer imaging in medicine, applications to materials science are thought to exist. The ability to perform M6ssbauer spectroscopy in an imaging mode, i.e., to do spectroscopy as a function of position within a sample, rather than in bulk, should prove to be of value in the analysis of heterogeneous materials. Although NMR and Mossbauer imaging may have different applica-326 tions, it is instructive to compare the two techniques. In NMR imaging, resonant-frequency information is translated into spatial information by imposing a magnetic field gradient on a system of precessing nuclear spins. In Mossbauer imaging, we will see that by rotating the y-ray absorber, a velocity gradient is imposed along the radial direction of the absorber and the resonance associated with a particular Doppler shift (analogous to a particular NMR precessional frequency) takes place along a family of lines perpendicular to the velocity gradient (analogous to the NMR magnetic field gradient). This process gives rise to a set of measurements of the line integrals of the absorption coefficient, from which an image can be recovered using well-known tomographic image reconstruction algorithms. These algorithms are mathematically similar to those employed in both x-ray and NMR tomography [4,5]. A more detailed description of the process by which the line-integral measurements arise is given in section 3.
The basis of the imaging approach can be most easily described using a classical interpretation of Mossbauer absorption, which will suffice for our purposes. As noted earlier, in this description the relative source-absorber velocity imparts a Doppler shift to the y-ray. For purposes of illustration, consider the simplest case of a material with a single resonance line (neglecting isomer shifts); then resonance will take place at zero Doppler shift, or zero relative velocity between the source and absorber. By rotating an extended absorbing object relative to the source, lines of constant Doppler shift (or equivalently, of constant y-ray energy) are generated across the absorber. As a result, resonant absorption takes place along one line at a time, giving rise to one line-integral measurement of the absorption coefficient. The location of the line is determined by three parameters: the absorber rotational velocity, the source velocity and the instantaneous position of the absorber during the measurement. From a complete set of such lineintegral measurements, a spatial map of the y-ray absorption can, in principle, be tomographically reconstructed. We shall see that the spatial resolution within the reconstructed image is proportional to the ratio of the naturallinewidth of the resonance (in units of velocity) to the rate of rotation of the absorber, but is ultimately signal-to-noise limited. This means that, in principle, arbitrarily high resolution can be achieved, but at the expense of rapidly increasing detector integration times. In particular, the integration time will increase in direct proportion to the number of resolution elements (pixels) in the desired image. Thus, in a two-dimensional image, doubling the linear resolution will square the number of pixels, which in turn will square the required integration time. As a consequence, very high resolution will require exceedingly long integration times unless sources of high intensity become more readily available.
At this stage one might observe that the equivalent line integral measurements could be generated merely by directing a collimated beam of y-rays at the absorber over a wide range of incident angles. This approach is analogous to conventional x-ray computerized tomography in a transmission experiment or to single photon emission tomography in a scattering experiment. The advantage of the proposed method is that an extended source can be used rather than a collimated source, which implies a substantially greater incident y-ray flux and a corresponding enhancement in the signal-to-noise ratio.
In the following three sections, we consider the idealized case of the y-rays incident parallel to the x -axis. This can be achieved by moving the source sufficiently far away so that the source subtends a small angle at the absorber. The assumption of parallel incidence simplifies the analysis considerably and permits an analytical solution to the image reconstruction problem; it also allows the analytical inversion of the complete M6ssbauer spectrum at each point (section 3). Moreover, the parallelincidence case also permits the closed-form derivation of the image point spread function, which is defined as an image of a point object (Dirac delta function) and characterizes the resolving power of the imaging technique (section 4). In section 5 we consider an extended source which gives rise to y-rays incident over a range of directions that are no longer parallel. In the latter case, the image reconstruction cannot be inverted analytically as for parallel incidence, but the reconstruction can nonetheless be performed using iterative algebraic reconstruction techniques [4,5]' Although such iterative reconstruction methods often fail to provide the insight and intuition of analytical solutions (which permit, for example, the closed-form derivation of the image point spread function), such methods are often more flexible in incorporating a priori information into the inversion algorithm. An important example is the incorporation of y-ray attenuation due to a variety of scattering and absorption (resonant and nonresonant) mechanisms. For relatively thick objects, attenuation can no longer be ignored in the inversion. For the sake of tractability, we will, however, ignore attenuation in the next several sections. 327

Standard Experiment
Suppose a y-ray emitted by a moving source is incident on a two-dimensional stationary absorber ( fig. 2). For our purposes, it will suffice to describe the resonant absorption classically by modeling the nuclear resonance phenomenon, within a single absorbing nucleus, as a harmonic oscillator characterized by resonance frequency w. and linewidth T (analogous to the mean lifetime of the y-ray energy level) and driven by the electromagnetic field of the incident y-ray. Letting A (t) denote the response of this elemental oscillator, we have is the y-ray electromagnetic field. In eq (2), v, is the velocity of the moving source, c is the velocity of light and Ao is a constant. Although one might also choose to include a finite linewidth in the y-ray field Ainc(t), we shall, for mathematical simplicity. lump alilinewidth effects into the parameter T in eq (1). Note that wovjc in eq (2) is just the classical Doppler shift due to the motion of the source. In an equivalent notation, wo=EoIft and r=ftIT.
where Eo is the energy of the incident y-ray and r is the (energy) linewidth. Now in eq (I), set where v. is a parameter (in units of velocity) that allows for shifts in the nuclear resonance of a given nucleus. Thus, the "intrinsic velocity," Va, is characteristic of different absorbing nuclei (or their environment). We will later allow Va to vary with position. The solution of the differential eq (1) on substituting eq (2) is where (5) is the Doppler-shifted source frequency. Substituting eqs (3) and (5) into eq (4), and assuming v, and v.< <c, eq (4) becomes Ao e""" Classically, the intensity of the re-emitted y-ray is proportional to IA 1 2 , i.e., or, in an equivalent notation, which is the familiar Lorentzian (or Breit-Wigner) distribution [1,2]. Because the detection of the re-emitted y-ray is an incoherent process, we can regard eq (6) as proportional to the probability of the y-ray re-emission as a function of the velocity To make eq (6) a valid probability density function in v, normalize eq (6) to have unit area when integrated with respect to v, i.e., define (7) where is the linewidth in units of velocity. From eq (7),

328
One can also show that In this and subsequent integrals, we shall assume that the limits of integration are from -00 to + 00 unless otherwise indicated. We can assume this by defining functions such as O'(x, y, vJ to be zero outside their support (Le., beyond the boundaries of the absorber and the domain of Va, and so on). Now, for convenience, suppose we temporarily consider the limit of very small V T (narrow resonant Iinewidth). As V T goes to zero, the probability density function P(T)(Va-v s ) becomes highly peaked about Va = v" and in the limit v T --.(), we can use eq (8) in (9). Thus, defining we have from eqs (8) and (9), which is the expected result for the standard Mossbauer experiment, that is, f(v s ) is the absorption cross section evaluated at Va = Vs averaged over the sample.

Imaging with Parallel y-Ray Incidence
Now let the absorber rotate clockwise about the z-axis at constant angular velocity n. Let (x,y) denote a stationary coordinate system and let (x',y') denote a moving coordinate system embedded in the rotating absorber, as illustrated in figure  3. At 1=0, let the stationary (x,y) and the rotating (x',y') coordinate systems coincide; at this instant denote the absorber's cross-section function by O'o(x,y, va). For 1>0, the absorber has rotated through an angle nt. Define O',(x,y, va) as the ab-sorption cross section relative to the stationary axes (x,y) at time t; then CT,(X,y, Va)=Uo(x',y', va), where x'=x cos Ot-y sin Ot (lla) or, inverting, Again consider a y-ray incident parallel to the will then take place when the total velocity v.+v: of a point on the absorber equals the source velocity v" where Va is the intrinsic material "velocity" and va' is the additional velocity contribution due to the motion of the absorber. Since the y-ray is incident parallel to the x-axis, va' will be the x-component of the absorber velocity; thus, from eq (12),  . 3). Now, in view of eq (13), for the rotating absorber, we can replace v, in eq (9) by v,-Oy. giving (14) Now considering a very narrow resonance, we again take the limit v.-+O. Defining ,.-00 and using, from eq (8), (14), we have This is our fundamental equation from which we shall solve for the function CTJ.,X',y', va) from the dataf,(v" 0) recorded as a function of v" n, and t.
Below, we consider four special cases of increasing generality: 1) one-dimensional imaging with no Va distribution (i.e., assuming a single spectral line at va=O); 2) two-dimensional imaging with no v. distribution; 3) one-dimensional image with an arbitrary and unknown Va distribution; and finally 4) two-dimensional imaging with an arbitrary and unknown Va distribution.
One-Dimensional Imaging with No v. Distribution. We consider the simplest imaging case here.
Suppose the absorber is one-dimensionaJ and rotating. For convenience, assume the amorber is aligned with the (rotating) y'-axis, and suppo~ further that the measurements are made at the in~t3nt the absorber rotates past vertical, i.e .• when the rotating absorber coincides with the y-axis at 1=0.
Thus, the spatial dependence of the one-dimensional cross section <To(y) is given explicitly as a function of vjfi.

Two-Dimensional Imaging (Tomography) with No Va
Distribution. Here we let so that eq (15) becomes Recall that the measurementsf,(v., fi) are made as a function of time as the absorber rotates. To make the t-dependence in eq (17) explicit, rewrite eq (17) as follows.
where eqs (10) and (l2b) were used in the last step. From eq (lS),f,(v" 0) is a line integral through the two-dimensional function <To(x', y') along a line parameterized by the radial distance vjO and the angle Ot from they'-axis. The inversion ofeq (IS) for <To(x',y') can be carried out by standard tomographic techniques. It is easiest to invert eq (IS) using Fourier methods, as follows. First define the one-dimensional Fourier trans- Inserting eq (IS) into (19)  Equation (22) is the well-known central-slice theorem in tomography [4,5]. In the present context it states that l,Ck, 0), evaluated at time t and spatial frequency k, is the two-dimensional spatial Fourier transform, O-o(k" , ky), of uo(x', y') evaluated on a line in the (k", ky ) Fourier plane through the origin at angle Ot with the k,,-axis. Thus, taking many revolutions of data and, during each revolution, stepping the value of v. will generate sufficient data to provide complete coverage of the (k", ky) Fourier plane, whereupon an inverse twodimensional Fourier transform of Cro(k", ky) yields the reconstruction <To(x', y'). Note that the required range of measured v. values is from 0 to fiRtnaJ., where Rmax is the largest radial dimension of the absorber. For the one-dimensional absorber, again assume that the measurement is made as the absorber rotates past vertical at t =0. Then, substituting <T,(x,y, ' Ya} into eq (15) and setting t=O, gives

One-Dimensional Imaging with Unknown
We wish to solve for the two-dimensional function <To(y, val, which has one spatial dimension (y) and one spectral dimension (va>. Rewrite eq (23) as Substituting into eq (19) and interchanging orders of integration, gives  It(k, n)=uo(-kn sin nt, kn cos nt, k). This is the generalization of the central-slice theorem to three-dimensions. Thus, by suitably varying the three parameters v" nand t, full coverage in the three-dimensional Fourier space can be achieved, and an inverse Fourier transform of uo(k", ky, k.) yields the reconstruction CTo(x',y', vJ.

Spatial Resolution and the Image Point Spread Function
In previous sections, we considered the idealized case of an infinitesimally narrow linewidth v r • This simplifies the mathematics, but implies that spatial resolution can be arbitrarily small. In this section, we show that spatial resolution is controlled by the naturallinewidth v r ; in particular, spatial resolution is proportional to the ratio of Vr to the rotational velocity n of the absorber.
Below we derive the point spread function (PSF) of the two-dimensional tomographic reconstruction problem, with no v. dependence for simplicity (the second case above). The PSF is by definition merely the reconstructed image of a twodimensional Dirac delta function, or point object, and its width provides a reasonable measure of image resolution. Alternatively, the reconstructed image may be regarded as the true image convolved with the PSF.
To derive the point spread function in two dimensions, we repeat the above derivation using eq (14) in place of eq (15). Substituting is the Fourier transform of P(T) ( -v.). Equation (29) is merely eq (22) multiplied by P (T)(k).
We can evaluate P(T)(k) by inserting eq (7) into (30); integrating, the result is Thus, the Fourier transform, Cro(k", ky), of O'o(x', y') is weighted by the exponential radiallydependent function exp(kv T ). The point spread function is obtained by setting Cro(k x , k,)= 1 in eq (29); as a result, the PSF is just the two-dimensional inverse Fourier transform of P(T)(k). Since k fi is the radial coordinate in the Fourier plane, P(T)(k) has radial symmetry and the resulting PSF also has radial symmetry. Hence, the inverse twodimensional Fourier transform of p(r)(k) can be evaluated in polar form: where J o (') is the zero-order Bessel function and ,'= v' X,I + y'2. This is the desired point spread function whose width gives the image resolution. Let ro' signify the full width at half maximum of eq (31); then rrJ =0.76 vTIO.
The important thing to note is that '0' is inversely 332 proportional to the rotational velocity fi. Finally, one can show that PSF(,') becomes a two-dimensional Dirac delta function in the limit vT-D, as required. As an example, for typical linewidths, v" of less than a mm/s, this result shows that '0' is on the order of a tenth of a millimeter when the absorber rotates at one revolution per second (0=217') .

Imaging with an Extended Source (Non-Parallel Incidence)
In the latter two sections, we assumed parallel y-ray incidence in the direction of the x-axis. Parallel incidence is easily achieved by placing the source far enough away so that it subtends a small angle at all points on the absorber. This is practical when the source is sufficiently strong to provide a y-ray flux intercepting the absorber that results in an acceptable signal-to-noise ratio. Unfortunately, arbitrarily strong M6ssbauer sources are not available, so that moving the source close to the absorber to maximize the incident flux rnay be important.
In the latter case, the source must be regarded as spatially extended, say, in the y-direction ( fig. 4), and the incident y-rays are no longer parallel. Assume that the extended source is moving in the xdirection at velocity v. and that the absorber rotates with angular velocity fi as before. Also, for simplicity, assume that the intrinsic velocity v.=O. In this geometry, one can show that the component of velocity of a point on the rotating absorber and of a point on the source in the direction of the line joining them (i.e., along the line of flight of a y-ray passing between the two points) are equal to each other (the resonance condition) if and only if the line defined by the two points passes through (x=O,y=v.lfi) in the stationary (x,y) coordinate system. To see this, refer to the line L in figure 5, which we suppose is the path of a y-ray emitted from point A on the source. Suppose this path intersects the y-axis at height R, as shown. Then L is tangent to a circle around the axis of rotation of radius R cos (), where () is the angle between Land the x-axis. From previous arguments, all points in the absorber lying on L have the same component of velocity in the direction of L, namely OR cos (). But the additional component of velocity of the incident y-ray along L due to the motion of the source is v. cos (), as can be seen from figure 5. The resonance condition is then obtained by equating the absorber and source components of velocity along L; thus setting OR cos () = v. cos () results in R =vJ'O. Consequently, only those 'Y-rays moving along paths through (x =0, y =vJ'O) are resonantly absorbed.
The complete set of lines passing through this point and intercepting the source generates a sector, as shown in figure 4, where the boundaries of the sector are the two lines intersecting the end points of the extended source and coming together at the point (0, vJ'O) in the (x,y) system. The density of lines peaks at the point (0, vJ'O) and falls off in proportion to 1I1x I in the x-direction ( fig. 6).
Through the peak in the y-direction, the width of this function is much narrower; in particular, the width in the y-direction can be shown to be vTIO, which is approximately the width of the point spread function derived previously for parallel incidence. Neglecting, for the moment, attenuation of the 'Y-rays through the absorber, the line-density function sketched in figure 6 may be equivalently interpreted as the probability density function of resonant absorption events per unit area, which is just proportional to the inverse of the sector area shown in figure 4.
The highly peaked distribution of the probability of resonance events is significant, since it implies that merely by varying Vs and counting absorption events as a function of time, it should be possible to generate a blurred image even without image reconstruction. Figure 4 shows that, by rotating the absorber, the 'Y-rays are effectively "focused" through the point (0, vJ'O). This focusing-like effect can be used to create a blurred image without further processing. A complete or deblurred image reconstruction cannot in this case be carried out analytically, but the reconstruction nevertheless can be performed using iterative algebraic reconstruction techniques (ART) [4,5]. The ART algorithms assume an initial estimate of the unknown image and then, on the basis of this estimate, recompute the measurements (line integrals of the resonant absorption coefficient). The recomputed measurements are then subtracted from the true measurements, giving rise to a set of error terms, which are then used to systematically adjust the current image estimate to reduce this error. When this iterative procedure is continued, it can be shown that the overall error tends to zero and, in the absence of measurement uncertainty and incomplete data, the reconstructed image converges to the true image. Such iterative approaches also make it relatively easy to incorporate attenuation effects in the algorithm.
In the above discussion, we assumed for simplicity that the extended source was one-dimensional (Le., extended along the y-direction). If the source 333 has some depth as well as height (Le., some extension in the z-direction), the additional degree of freedom arising from variations in the x-component of the 'Y-ray velocity will broaden somewhat -Q _v.  the region of sensitivity beyond the sector illustrated in figure 4. To minimize this broadening effect, the source should be collimated in the z-direction as much as practical with a suitable arrangement of slits in this direction. A three-dimensional absorber could then be reconstructed on a "slice-by-slice" basis, as in two-dimensional x-ray tomography.

Applications
Applications of the technique to interior imaging appear to be limited to lighter materials, or absorbers of small size, that permit sufficient y-ray penetration. One possibility is the imaging of composite materials, in which a Mossbauer element is embedded in a lighter matrix. Another possibility is surface imaging, whereby the surface of a flat absorber is rotated about an axis normal to the surface and a one-dimensional extended source is displaced a small amount from the plane of the surface. If all the impinging 'Y-rays have the same (or nearly the same) z-component of velocity, the region of sensitivity in which resonant absorption takes place will look much the same as the sector illustrated in figure 4. In this case, the value of v, used previously should be replaced by the component of the source velocity in the plane of the surface.
Some potential high-resolution applications in materials science include the imaging of grain boundary segregation. the imaging of residual stress distributions, and the imaging of the distribution of magnetism in ferromagnetic materials due to perturbations in the Mossbauer spectra arising from variations in the magnetic hyperfine field.

Conclusion
In conventional M6ssbauer spectroscopy, the measurement is spatially integrated over the absorber. For parallel 'Y-ray incidence, a rotating absorber generates a line of "sensitivity" Over which resonant absorption takes place. where the location of the line is a function of the source velocity and the absorber rotational velocity. In particular, this line is parallel to the velocity of the source and passes through the point (x =0, y=v.lfl). For an extended source, the line broadens into a sector subtending the source and coming together at the point (x = 0, y = v.ln). For a highly extended source, the result is as if the 'Y-rays were effectively "focused" through the point (x = 0, y = v'/n). This focusing effect can be used to generate a blurred 334 image of the absorption coefficient distribution in the absence of further processing when the measurements are recorded as a function of v, and time (or the instantaneous absorber rotational position). Furthermore, by exploiting the three measurement degrees of freedom, v" 0, and time (or rotational position), complete spectroscopic information can, in principle, be recovered as a function of position. Thus, Mossbauer spectroscopy can be performed in an imaging mode. Spatial resolution is proportional to the ratio of the naturallinewidth V T to the rotational velocity fl, but becomes rapidly signalto-noise limited as the resolution increases.