research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

3D reconstruction of magnetization from dichroic soft X-ray transmission tomography

aALBA Synchrotron, Cerdanyola del Vallès, Barcelona 08290, Spain, bAdvanced Photon Source, Argonne National Laboratory, Argonne, IL 60439, USA, cDepartment of Electrical Engineering and Computer Science, Northwestern University, Evanston, IL 60208, USA, dMaterials Science Division, Argonne National Laboratory, Argonne, IL 60439, USA, eDepartamento de Física, Universidad de Oviedo, Oviedo 33007, Spain, and fCentro de Investigación en Nanomateriales y Nanotecnología, CINN (CSIC – Universidad de Oviedo), El Entrego 33940, Spain
*Correspondence e-mail: ferrer@cells.es

Edited by P. A. Pianetta, SLAC National Accelerator Laboratory, USA (Received 3 January 2018; accepted 14 April 2018; online 30 May 2018)

The development of magnetic nanostructures for applications in spintronics requires methods capable of visualizing their magnetization. Soft X-ray magnetic imaging combined with circular magnetic dichroism allows nano­structures up to 100–300 nm in thickness to be probed with resolutions of 20–40 nm. Here a new iterative tomographic reconstruction method to extract the three-dimensional magnetization configuration from tomographic projections is presented. The vector field is reconstructed by using a modified algebraic reconstruction approach based on solving a set of linear equations in an iterative manner. The application of this method is illustrated with two examples (magnetic nano-disc and micro-square heterostructure) along with comparison of error in reconstructions, and convergence of the algorithm.

1. Introduction

Advances in nanomagnetism towards applications in spintronics involve magnetic heterogeneous systems with increasing complexity including multiple materials, and complex geometries. Spectroscopic methods and imaging tools are required to characterize and visualize the local magnetic properties in different regions throughout the heterostructures. Several tools for magnetic imaging have been developed in the past years. Using visible-light photons, Kerr microscopy and vector magnetometry provide, via magneto-optical effects, mappings of the magnetization with lateral resolutions around ∼1 µm limited by the photon wavelengths (Hubert & Schaëfer, 1998[Hubert, A. & Schaëfer, R. (1998). Magnetic Domains: The Analysis of Magnetic Microstructures. Berlin/Heidelberg: Springer.]; Berger & Pufall, 1999[Berger, A. & Pufall, M. R. (1999). J. Appl. Phys. 85, 4583-4585.]; Teixeira et al., 2011[Teixeira, J. M., Lusche, R., Ventura, J., Fermento, R., Carpinteiro, F., Araujo, J. P., Sousa, J. B., Cardoso, S. & Freitas, P. P. (2011). Rev. Sci. Instrum. 82, 043209.]). Their probing depth is a few nanometres in metallic samples. Higher lateral resolution is achieved with electron-based techniques by using secondary or photo-emitted electrons. Scanning electron microscopy with polarization analysis and X-ray photoemission electron microscopy provide very good lateral resolution (∼10 nm) but they also have shallow probing depths (1–2 nm) due to strong inelastic scattering of low-energy electrons (Scheinfein et al., 1990[Scheinfein, M. R., Unguris, J., Kelley, M. H., Pierce, D. T. & Celotta, R. J. (1990). Rev. Sci. Instrum. 61, 2501-2527.]; Unguris, 2001[Unguris, J. (2001). Magnetic Imaging and its Application to Materials, edited by M. De Graef & Y. Zhu, Vol. 36, pp. 167-194. San Diego: Academic Press.]; Locatelli & Bauer, 2008[Locatelli, A. & Bauer, E. (2008). J. Phys. Condens. Matter, 20, 093002.]; Foerster et al., 2016[Foerster, M., Prat, J., Massana, V., Gonzalez, N., Fontsere, A., Molas, B., Matilla, O., Pellegrin, E. & Aballe, L. (2016). Ultramicroscopy, 171, 63-69.]). This limitation is circumvented for high-energy electrons that traverse samples. Lorentz transmission electron microscopy (Lorentz TEM) (De Graef, 2001[De Graef, M. (2001). Magnetic Imaging and its Application to Materials, edited by M. De Graef & Y. Zhu, Vol. 36, ch. 2. San Diego: Academic Press.]; McVitie et al., 2015[McVitie, S., McGrouther, D., McFadzean, S., MacLaren, D. A., O'Shea, J. J. & Benitez, M. J. (2015). Ultramicroscopy, 152, 57-62.]) has an excellent lateral resolution (∼1 nm) and is sensitive to the magnetization of the whole sample in systems up to ∼100 nm in thickness. Because the technique probes the magnetization perpendicular to the direction of the propagation of electrons, it is highly sensitive to in-plane magnetization. A further development of Lorentz TEM, called vector field electron tomography, is especially interesting (Phatak et al., 2008[Phatak, C., Beleggia, M. & De Graef, M. (2008). Ultramicroscopy, 108, 503-513.]; Prabhat et al., 2017[Prabhat, K. C., Aditya Mohan, K., Phatak, C., Bouman, C. & De Graef, M. (2017). Ultramicroscopy, 182, 131-144.]). It combines both the excellent lateral resolution and the in-depth magnetic sensitivity of Lorentz TEM with the three-dimensional (3D) volume reconstruction capabilities of tomography to obtain the potential vector A and, from it, the full magnetization configuration in 3D magnetic systems. Besides these methods, two X-ray-based techniques are being used to probe magnetization with high sensitivity and spatial resolution. Soft X-ray transmission microscopy using circularly polarized X-rays of energies typically below 1 keV with a 100–300 nm penetration depth has been used to image buried magnetic layers by exploiting the large magnetic dichroism occurring at specific electronic transitions in magnetic atoms (Blanco-Roldán et al., 2015[Blanco-Roldán, C., Quirós, C., Sorrentino, A., Hierro-Rodríguez, A., Álvarez-Prado, L. M., Valcárcel, R., Duch, M., Torras, N., Esteve, J., Martín, J. I., Vélez, M., Alameda, J. M., Pereiro, E. & Ferrer, S. (2015). Nat. Commun. 6, 8196.]; Fernández-Pacheco et al., 2017[Fernández-Pacheco, A., Streubel, R., Fruchart, O., Hertel, R., Fischer, P. & Cowburn, R. P. (2017). Nat. Commun. 8, 15756.]). More recently, hard X-rays (photon energies of about 6–10 keV), that have larger penetration depths (∼µm), have been successfully used to reconstruct 3D magnetization. In this case, as the magnetic dichroic absorption is very small compared with soft X-rays, the method uses transversally coherent X-rays to exploit diffraction and phase contrast in ptychography mode, which involves the acquisition of several thousands (∼105) of diffraction images, allowing to resolve the magnetization of the sample with lateral resolution around 100 nm typically (Donnelly et al., 2017[Donnelly, C., Guizar-Sicairos, M., Scagnoli, V., Gliga, S., Holler, M., Raabe, J. & Heyderman, L. J. (2017). Nature (London), 547, 328-331.]).

Here we present the development of a new iterative algorithm to obtain quantitative 3D vector magnetization reconstruction using soft X-ray microscopy (Blanco-Roldán et al., 2015[Blanco-Roldán, C., Quirós, C., Sorrentino, A., Hierro-Rodríguez, A., Álvarez-Prado, L. M., Valcárcel, R., Duch, M., Torras, N., Esteve, J., Martín, J. I., Vélez, M., Alameda, J. M., Pereiro, E. & Ferrer, S. (2015). Nat. Commun. 6, 8196.]; Hierro-Rodriguez et al., 2017a[Hierro-Rodriguez, A., Quirós, C., Sorrentino, A., Blanco-Roldán, C., Álvarez-Prado, L. M., Martín, J. I., Alameda, J. M., Pereiro, E., Vélez, M. & Ferrer, S. (2017a). Phys. Rev. B, 95, 014430.],b[Hierro-Rodriguez, A., Quirós, C., Sorrentino, A., Valcárcel, R., Estébanez, I., Álvarez-Prado, L. M., Martín, J. I., Alameda, J. M., Pereiro, E., Vélez, M. & Ferrer, S. (2017b). Appl. Phys. Lett. 110, 262402.]). This is achieved by taking X-ray transmission images with opposite dichroism (positive and negative) for two tomogram series by rotating the object around two orthogonal tilt axes. The total number of images required is around 500 or less. The reconstructions are based on a joint processing of both tomogram series for obtaining the reconstructed magnetic configuration. In §2[link] the forward problem of magnetic soft X-ray transmission microscopy is presented by introducing the equations that describe the projected images. In §3[link], we analyze scalar and vector field reconstruction problems solving them by using a modified iterative algebraic reconstruction technique. The algorithm will be made openly available under TomoPy, which is a library for tomographic image reconstruction (Gürsoy et al., 2014[Gürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188-1193.]), and the simulated tomograms will be accessible through TomoBank (De Carlo et al., 2018[De Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, J., Ludwing, W., Mancini, L., Marone, F., Mokso, R., Pelt, D., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004.]). §4[link] illustrates the application of the method by reconstructing two simulated magnetic microstructures and evaluating its accuracy. Finally, the conclusions of the work are presented in §5[link].

2. Magnetic soft X-ray transmission microscopy

In an X-ray transmission microscope, the transmitted X-ray intensity through the sample under investigation is recorded at each pixel of a two-dimensional detector forming a transmission image. In the following, we assume a simplified geometry with incoming parallel beam, although the condenser optics of the soft X-ray microscope focuses the beam onto the sample, and the objective Fresnel zone plate lens (FZP) collects the transmitted beam producing a magnified image (×1500) at the charged-coupled detector (CCD) (Pereiro et al., 2009[Pereiro, E., Nicolás, J., Ferrer, S. & Howells, M. R. (2009). J. Synchrotron Rad. 16, 505-512.]). The FZP has a limited depth of field that will affect the projections while rotating, if the sample size exceeds it. This has to be taken into account for samples with relatively large lateral dimensions. By exploiting the broad photon energy spectrum emitted by synchrotron light-sources, it is possible to tune the X-ray wavelength in order to match atom-specific absorption edges leading also to resonant atom-specific images (Pereiro et al., 2009[Pereiro, E., Nicolás, J., Ferrer, S. & Howells, M. R. (2009). J. Synchrotron Rad. 16, 505-512.]; Olivares-Marín et al., 2015[Olivares-Marín, M., Sorrentino, A., Lee, R. C., Pereiro, E., Wu, N. L. & Tonti, D. (2015). Nano Lett. 15, 6932-6938.]). Moreover, if the polarization of the incident beam is circular (right- or left-handed), atom-specific magnetic images can be recorded by taking advantage of magnetic dichroism effect (Olivares-Marín et al., 2015[Olivares-Marín, M., Sorrentino, A., Lee, R. C., Pereiro, E., Wu, N. L. & Tonti, D. (2015). Nano Lett. 15, 6932-6938.]; Sorrentino et al., 2015[Sorrentino, A., Nicolás, J., Valcárcel, R., Chichón, F. J., Rosanes, M., Avila, J., Tkachuk, A., Irwin, J., Ferrer, S. & Pereiro, E. (2015). J. Synchrotron Rad. 22, 1112-1117.]). In this framework, the X-ray intensity after passing through the sample can be described as follows,

[I= I_0\exp \left\{ \textstyle\int L^{-1} \left[ 1+\delta \left( {\bf{k}}\cdot{\bf{m}} \right) \right]\,{\rm{d}}t \right\}. \eqno(1)]

In equation (1)[link], I and I0 are the transmitted and incident intensities of the X-ray beam, respectively. L−1 is the inverse of the attenuation length for the X-rays, and it depends on the photon energy and on the electron density of the sample which will be variable in heterogeneous systems. δ is the dichroic coefficient for the magnetic material under analysis; it is the scale factor of the magnetic sensitivity and depends on the electronic levels of the absorbing atoms. [{\bf{k}}\cdot{\bf{m}}] is the dot product of the X-ray wavevector and magnetization; it provides the sensitivity of the dichroism at different relative orientations of the sample and the photon beam. m is the reduced magnetization vector (m = M/MS, with M the magnetization vector and MS the saturation magnetization). dt is the elementary path along the X-ray linear trajectory spanned by the line integral. The latter runs along the entire beam path, from its source to the detector, passing through the sample space. Thus, it is clear that L−1 and m are explicit functions of t as they are sampled by the X-ray beam. The transmittance (T = I/I0) is separated into two terms,

[T=\exp \left\{ \textstyle\int L(t)^{-1}\,{\rm{d}}t + \textstyle\int L(t)^{-1} \delta(t) \left[ {\bf{k\cdot{m}}}(t) \right] \,{\rm{d}}t \right\}. \eqno(2)]

The first integral does not depend on the magnetism but only on the charge distribution in the sample whereas the second one includes the magnetic contributions. For practical convenience and to separate the magnetic and nonmagnetic parts we take the logarithms of the transmittance for right-handed [+δ, equation 3(a)[link]] and left-handed [−δ, equation (3b)[link]] polarizations,

[\eqalignno{ \log\left[T_{+\delta}\right] & =\textstyle\int L(t)^{-1}\,{\rm{d}}t \,+\, \textstyle\int L(t)^{-1} \delta(t) \left[ {\bf{k\cdot{m}}}(t) \right] \,{\rm{d}}t, &(3a) \cr \log\left[T_{-\delta}\right] & =\textstyle\int L(t)^{-1}\,{\rm{d}}t \,-\, \textstyle\int L(t)^{-1} \delta(t) \left[ {\bf{k\cdot{m}}}(t) \right] \,{\rm{d}}t. &(3b) }]

Hence, by simply adding and subtracting equations (3a) and (3b), separate expressions for the non-magnetic [equation (4a)[link]] and magnetic [equation (4b)[link]] contributions are obtained,

[\eqalignno{ \log\left[T_{+\delta}\right] + \log\left[T_{-\delta}\right] & = 2\textstyle\int L(t)^{-1}\,{\rm{d}}t, & (4a) \cr \log\left[T_{+\delta}\right] - \log\left[T_{-\delta}\right] & = 2\textstyle\int L(t)^{-1} \delta(t)\left[ {\bf{k\cdot{m}}}(t) \right] \,{\rm{d}}t. & (4b) }]

Equation (4a)[link] will be used to obtain the values of the attenuation length (L), which is a scalar field, while equation (4b)[link] will allow the magnetization configuration (m) of the system to be extracted.

3. Scalar and vector field tomographic reconstruction

In a general tomography problem, a certain property of an unknown sample is volume reconstructed by the analysis of its different projections (Kak & Slaney, 1988[Kak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographic Imaging. New York: IEEE Press.]). These are acquired by taking images of the sample at different rotation angles around a certain axis or axes.

Fig. 1(a)[link] shows an image of the general problem for transmitted X-rays (red arrows) and two orthogonal rotation axes. We indicate as Xtilt and Ytilt the projection series acquired by rotating the system around the X and Y axes, respectively. The sample is located between the X-ray source and the detector. Both rotation axes pass through the centre of the sample taken as the origin of the reference frame. In the sample space, a volume model formed by voxels is created (grey cube around the sample in Fig. 1a[link]) with its centre located at the origin.

[Figure 1]
Figure 1
(a) Scheme of the tomography problem showing the X-ray beam, the sample, the rotation axes, the detector and a volume model for the reconstruction. Two projections at 0° (b) and 60° (c) rotated around the Y axis (plane XZ) on a one-dimensional detector are also presented.

Let us call x the property that we want to reconstruct, thus each voxel of the model will contain a value for this parameter. In order to find the x field by using the volume model, the images recorded with the detector at different projection angles (measured data) are compared with the mathematic projection of the model at the same tilt angles. In this way, the problem of reconstructing the x field can be written as a system of linear equations in the following form,

[{y^\phi} - {A^\phi}x = 0, \eqno(5)]

where

[{y^\phi} = \left [{\matrix{ {y_{1,1}^\phi} \cr \vdots \cr {y_{n,m}^\phi} \cr \vdots \cr {y_{N,M}^\phi} \cr } } \right], \quad x = \left [{\matrix{ {{x_{1,1,1}}} \cr \vdots \cr {{x_{i,j,k}}} \cr \vdots \cr {{x_{I,J,K}}} \cr } } \right],]

[{A^\phi} = \left [{\matrix{ {l_{1,1,1}^{\,\phi, 1,1}} & \cdots & {l_{i,j,k}^{\,\phi, 1,1}} & \cdots & {l_{I,J,K}^{\,\phi, 1,1}} \cr \vdots & \ddots & {} & {} & \vdots \cr {l_{1,1,1}^{\,\phi, n,m}} & {} & {l_{i,j,k}^{\,\phi, n,m}} & {} & {l_{I,J,K}^{\,\phi, n,m}} \cr \vdots & {} & {} & \ddots & \vdots \cr {l_{1,1,1}^{\,\phi, N,M}} & \cdots & {l_{i,j,k}^{\,\phi, N,M}} & \cdots & {l_{I,J,K}^{\,\phi, N,M}} \cr } } \right].]

The column vector yϕ represents the transmittance at each pixel of the detector arranged in raster order (concatenated stacking of all the rows of the detector) for a certain projection ϕ. The detector has N × M pixels. The column vector x is composed of the values for the x field stored in the volume model arranged also in raster order [concatenated stacking of volume model data ordered by running the rows (j index), the columns (i index) and finally the layers (k index)]. The general volume model has I × J × K voxels, indicating the number of rows, columns and layers, respectively. The matrix Aϕ is the projection matrix which allows the values of each detector pixel yϕ to be obtained as a function of the model parameters x for a certain projection angle ϕ. This matrix has size NM × IJK and is sparse. Its elements ([l_{i,j,k}^{\,\phi, n,m}]) are indexed indicating to which detector pixel (n,m) and at which tilt angle (ϕ) the volume is being projected, and also which cell of the model (i,j,k) is involved. The main point here is that, for different projections, different linear combinations of the voxels contribute to the integrated intensity in the same detector pixel. Thus, the value of each element of Aϕ is calculated as the length of a specific ray through each voxel at a projection angle ϕ. These lengths are the `weights' of the physical property enclosed in the voxels (x) to form the linear combination determined by the beam. It is important to mention here that an efficient implementation of the calculation of these elements is crucial for the performance of the final reconstruction algorithm (Siddon, 1985[Siddon, R. L. (1985). Med. Phys. 12, 252-255.]; Jacobs et al., 1998[Jacobs, F., Sundermann, E., De Sutter, B., Christiaens, M. & Lemahieu, I. (1998). J. Comput. Information Technol. 6, 89-94.]).

To further clarify this point, we analyse a simple example of two different projections in the Ytilt configuration at 0° (Fig. 1b[link]) and 60° (Fig. 1c[link]). The presented situation is reduced to a single row (one-dimensional) of the detector with three pixels (y1, y2, y3) and a two-dimensional slice (x1,1,…, xi,j,…, x3,3) in the XZ plane (Ytilt) of the volume model. The length of the X-ray beam inside each cell ([l_{i,\,j}^{\,\phi,n}]) is indexed using the previous convention. Thus, for instance the linear combination resulting in the signal integrated by pixel 2 for 0° and 60° projections can be written as

[y_2^0 = l_{1,2}^{\,0,2}{x_{1,2}} + l_{2,2}^{\,0,2}{x_{2,2}} + l_{3,2}^{\,0,2}{x_{3,2}}]

and

[y_2^{60} = l_{1,3}^{\,60,2}{x_{1,3}} + l_{2,3}^{\,60,2}{x_{2,3}} + l_{2,2}^{\,60,2}{x_{2,2}} + l_{2,1}^{\,60,2}{x_{2,1}} + l_{3,1}^{\,60,2}{x_{3,1}}.]

It is clear now that the line integral appearing in equation (4a)[link] is numerically reproduced by this linear combination along the X-ray path. The L−1 values are considered homogeneous inside each voxel and the infinitesimal line integral element (dt) is substituted by the length of the ray through the involved model cell. The resulting equations system has as many equations as the number of detector pixels times the number of different projection angles. For instance, a 256 × 256 pixel detector and 100 projections leads to a system with more than 6.5 million equations. To solve this problem we use the algebraic reconstruction technique (ART) (Kak & Slaney, 1988[Kak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographic Imaging. New York: IEEE Press.]). First, an initial value (it can be 0) is assigned to the parameters contained in the model voxels. After this, the volume model is projected into the detector space by using Aϕ for the initial tilt angle of the tomogram. By calculating the difference between the experimental data (yϕ) and the numerically projected one, an error vector (eϕ) is obtained as depicted in equation (6)[link],

[{e^\phi} = {y^\phi} - {A^\phi}x. \eqno(6)]

The model parameters are updated as indicated in equation (7)[link]. Vectors xnew and xold represent the model parameters after and before the update, respectively. Cϕ and Rϕ are diagonal matrices where each of their elements is calculated as the inverse of the sum of all column and row elements in matrix Aϕ, respectively. For their calculation, indices `c' and `r' indicate the column and row index of the projection matrix Aϕ, respectively, thus they run from 1 to IJK for index c and from 1 to NM for index r. These matrices are included in order to compensate for the number of beams interacting with the same voxel, and for the number of pixels which are hit by the same ray, thus they prevent overweighting. Finally, [Aϕ]T represents the transposed Aϕ matrix,

[x^{\rm{new}} = x^{\rm{old}} + \left[ C^{\,\phi}\left(A^{\phi}\right)^{\rm{T}}\,R^{\,\phi}\right]e^{\phi},\eqno(7)]

where

[{C^{\,\phi}} = \left [{\matrix{ {c_{1,1}^\phi} & 0 & \cdots & {} & 0 \cr 0 & \ddots & {} & {} & {} \cr \vdots & {} & {c_{c,c}^\phi} & {} & \vdots \cr {} & {} & {} & \ddots & 0 \cr 0 & {} & \cdots & 0 & {c_{IJK,IJK}^\phi} \cr } } \right], \quad c_{c,c}^\phi = {1 \over {\sum\limits_r {l_{r,c}^{\,\phi}} }},]

[R^{\,\phi} = \left [{\matrix{ {r_{1,1}^\phi} & 0 & \cdots & {} & 0 \cr 0 & \ddots & {} & {} & {} \cr \vdots & {} & {r_{r,r}^\phi} & {} & \vdots \cr {} & {} & {} & \ddots & 0 \cr 0 & {} & \cdots & 0 & {r_{NM,NM}^\phi} \cr } } \right], \quad r_{r,r}^\phi = {1 \over {\sum\limits_c {l_{r,c}^{\,\phi}} }}.]

Note here that an iteration is completed when all the recorded projections are taken into account to update the model; thus to complete an iteration by using the scheme shown in equation (7)[link] we apply the following protocol: (1) calculate the Aϕ matrix for the first projection angle and calculate the error vector eϕ for that angle, (2) update the model using equation (7)[link]. After this, go back to step (1) but using the second projection angle instead of the first one and calculate the error using the previously updated model, update the model again with the new error and continue repeating the protocol until the last projection is taken into account. By performing several iterations the solution for the reconstructed field converges, which is common to ART. We have chosen ART for the iterative reconstruction due to the fast convergence (8–10 iterations usually lead to convergence) and because the method does not assume any a priori information related to the noise or object model.

This method can be directly applied in order to reconstruct scalar fields for the situation shown in equation (4a)[link] to reconstruct the attenuation length in 3D. In order to apply the protocol to reconstruct the vector field as in the case of equation (4b)[link], it is only necessary to calculate a different projection matrix which takes into account the signal recorded with magnetic contrast due to the dot product. To do this we express the X-ray wavevector in spherical coordinates (kx = sinθcosφ, ky = sinθsinφ, kz = cosθ) and perform the dot product [equation (8)[link]],

[\eqalignno{ \log\left[T_{+\delta}\right]-\log\left[T_{-\delta}\right] = {}& 2\textstyle\int\limits L(t)^{-1}\delta(t) \Big[\sin\theta\cos\varphi\,m_x(t) \cr& \!+\sin\theta\sin\varphi\,m_y(t) + \cos\theta\,m_z(t)\big]\,{\rm{d}}t. &(8)}]

Note here that the previously used projection angle defined as ϕ is the same as θ in spherical coordinates. In the reference frame as sketched in Fig. 1[link], the Ytilt series implies a rotation of an angle θ with a fixed angle φ = 0° (defining the XZ plane). In the case of the Xtilt series, the fixed φ angle is 90° and the rotation angle is equal to θ (defining the YZ plane). This means that with only one tilt series it is not possible to reconstruct all the components of the magnetization vector; thus to obtain the necessary information we acquire two tomogram series: one around the Y axis (Ytilt) and the other around the X axis (Xtilt). The first one will give information about the mx and mz components [equation (9a)[link]], while the second one contains information from my and mz [equation (9b)[link]],

[\eqalignno{ & Y{\rm{tilt}} \,\,\longrightarrow \log\left[T_{+\delta}\right]-\log\left[T_{-\delta}\right] = &(9a) \cr& \qquad\qquad\quad 2\textstyle\int\limits L(t)^{-1}\delta(t)\Big[\sin\theta\,m_x(t) +\cos\theta\,m_z(t)\Big]\,{\rm{d}}t, \cr & X{\rm{tilt}} \,\,\longrightarrow \log\left[T_{+\delta}\right]-\log\left[T_{-\delta}\right] = &(9b) \cr& \qquad\qquad\quad 2\textstyle\int\limits L(t)^{-1}\delta(t)\Big[\sin\theta\,m_y(t) +\cos\theta\,m_z(t)\Big]\,{\rm{d}}t. }]

The vector field reconstruction projection matrix will need to take into account now not only the length of the ray through each voxel but also the projection angle sine and cosine due to the magnetic contrast. As we are working now with two different tilt series, for the same projection angle θ we have two different images acquired (one from Ytilt and the other from Xtilt). We can arrange the data in the form of a column vector as in the scalar case, but now concatenating the Xtilt data after the Ytilt data. Also, it is necessary to create a volume model where now each voxel contains three parameters which are the three magnetization vector components. In this way, the linear equations system for the vector field case can be written as follows,

[y^{\varphi,\theta} - A^{\varphi,\theta}x = 0, \eqno(10)]

where

[{y^{\varphi,\theta}} = \left [{\matrix{ {y_{1,1}^{0,\theta }} \cr \vdots \cr {y_{n,m}^{0,\theta }} \cr \vdots \cr {y_{N,M}^{0,\theta }} \cr {y_{1,1}^{90,\theta }} \cr \vdots \cr {y_{n,m}^{90,\theta }} \cr \vdots \cr {y_{N,M}^{90,\theta }} \cr } } \right], \quad x = \left [{\matrix{ {x_{1,1,1}^X} \cr \vdots \cr {x_{i,j,k}^X} \cr \vdots \cr {x_{I,J,K}^X} \cr {x_{1,1,1}^Y} \cr \vdots \cr {x_{i,j,k}^Y} \cr \vdots \cr {x_{I,J,K}^Y} \cr {x_{1,1,1}^Z} \cr \vdots \cr {x_{i,j,k}^Z} \cr \vdots \cr {x_{I,J,K}^Z} \cr } } \right],]

[{A^{\varphi, \theta }} = \left [{\matrix{ {\sin \theta \,{B^{0,\theta }}} & 0 & {\cos \theta \,{B^{0,\theta }}} \cr 0 & {\sin \theta \,{B^{90,\theta }}} & {\cos \theta \,{B^{90,\theta }}} \cr } } \right],]

[{B^{0,\theta }} = \left [{\matrix{ {l_{1,1,1}^{\,0,\theta, 1,1}} & \cdots & {l_{i,j,k}^{\,0,\theta, 1,1}} & \cdots & {l_{I,J,K}^{\,0,\theta, 1,1}} \cr \vdots & {} & \vdots & {} & \vdots \cr {l_{1,1,1}^{\,0,\theta, n,m}} & \cdots & {l_{i,j,k}^{\,0,\theta, n,m}} & \cdots & {l_{I,J,K}^{\,0,\theta, n,m}} \cr \vdots & {} & \vdots & {} & \vdots \cr {l_{1,1,1}^{\,0,\theta, N,M}} & \cdots & {l_{i,j,k}^{\,0,\theta, N,M}} & \cdots & {l_{I,J,K}^{\,0,\theta, N,M}} \cr } } \right],]

[{B^{90,\theta }} = \left [{\matrix{ {l_{1,1,1}^{\,90,\theta, 1,1}} & \cdots & {l_{i,j,k}^{\,90,\theta, 1,1}} & \cdots & {l_{I,J,K}^{\,90,\theta, 1,1}} \cr \vdots & {} & \vdots & {} & \vdots \cr {l_{1,1,1}^{\,90,\theta, n,m}} & \cdots & {l_{i,j,k}^{\,90,\theta, n,m}} & \cdots & {l_{I,J,K}^{\,90,\theta, n,m}} \cr \vdots & {} & \vdots & {} & \vdots \cr {l_{1,1,1}^{\,90,\theta, N,M}} & \cdots & {l_{i,j,k}^{\,90,\theta, N,M}} & \cdots & {l_{I,J,K}^{\,90,\theta, N,M}} \cr } } \right].]

The elements of the yφ,θ column vector are indexed as [y_{n,m}^{\varphi,\theta}] where φ and θ indicate whether the data are related to the Ytilt (φ = 0) or Xtilt (φ = 90) series, and the value of the projection angle, respectively. The sub-indices n and m specify the pixel in the detector as in the scalar case. Now, this column vector has 2(N × M) elements. The column vector containing the volume model parameters has 3(I × J × K) elements and is arranged by concatenating the X, Y and Z components of the vector property to be reconstructed. Their elements are labelled indicating the voxel index (i,j,k) and the vector component (X, Y or Z). Due to these new sizes for detector and volume model vectors the projection matrix Aφ,θ has 2(NM) × 3(IJK) elements, and due to the selected arrangement of yφ,θ and x vectors its arrangement is different from the scalar case. The matrix can be separated into six different blocks arranged in two rows and three columns. Each block is a sub-matrix of size NM × IJK and allows for the projection of a different vector component in the Ytilt or Xtilt situation. The first row is referred to the Ytilt projection and the second one to the Xtilt. The three columns are related by X, Y and Z vector field components, respectively. Thus, by using equations (9a) and (9b)[link], the first row supports the reconstruction of X and Z components of the vector field, while the second one deals with Y and Z. This means that first-row second-column and second-row first-column blocks are zeros. First-row first-column and second-row second-column blocks are multiplied by sin θ to project mx and my, respectively, and both row block-elements of the third-column project mz and are multiplied by cosθ. The base sub-matrices are indicated as B0,θ and B90,θ. They are calculated as the projection matrix in the scalar case and contain the lengths of the analysed ray passing through the model voxels in the Ytilt and Xtilt configurations, respectively. The elements of matrices B0,θ and B90,θ are labelled as [l_{i,j,k}^{\,0,\theta,n,m}] and [l_{i,j,k}^{\,90,\theta,n,m}] indicating the Ytilt (φ = 0) or Xtilt (φ = 90) configuration, the detector pixel (n,m) and the voxel index in the volume model (i,j,k).

ART can be directly applied to the equation system described in equation (10)[link] as it was applied for the reconstruction of a scalar field; the only difference is that matrices C and R must be calculated with Aφ,θ without multiplying its sub-matrix elements (B0,θ and B90,θ) by sinθ and cosθ. This is because the C and R matrices only compensate for the number of beams interacting with the same voxel and for the number of pixels which are hit by the same ray. Finally, it is important to note that we are not directly reconstructing the reduced magnetization vector; we are reconstructing 2L(t)−1δ(t)m(t) [equation (4b)[link]]. The contribution of the attenuation length can be easily accounted for by using the scalar field reconstruction using equation (4a)[link], and then using those values in the model to isolate δ(t)m(t). The latter is proportional to the magnetization configuration.

4. Reconstruction of magnetic micro/nanoparticles

Two different magnetic particles have been simulated in order to test the capabilities of the aforementioned reconstruction approach: a Permalloy (Ni80Fe20; Py) nano-disc and a Py/Al/Py micro-square heterostructructure (Fig. 2[link]).

[Figure 2]
Figure 2
Schematics of the magnetic nano-disc (a) and micro-square (b) systems simulated to test the vector reconstruction algorithm. The relaxed magnetic configuration of both structures is presented showing simulation volume slices in the disc (c) and square (d) cases.

The micromagnetic simulations have been performed by using the Mumax3 code (Vansteenkiste et al., 2014[Vansteenkiste, A., Leliaert, J., Dvornik, M., Helsen, M., Garcia-Sanchez, F. & Waeyenberge, B. V. (2014). AIP Adv. 4, 107133.]). The systems have been simulated using 5 nm × 5 nm × 5 nm cells with MS of Py = 810 × 103A m−1, exchange stiffness Aexch = 1.3 × 10−11 J m−1 and an out-of-plane uniaxial anisotropy energy density KU = 6 × 104 J m−3. The latter is usually present in Py structures and supports the formation of stripe domains in thicker samples (Voltan et al., 2016[Voltan, S., Cirillo, C., Snijders, H. J., Lahabi, K., García-Santiago, A., Hernández, J. M., Attanasio, C. & Aarts, J. (2016). Phys. Rev. B, 94, 094406.]), but it is not generally considered in thin film simulations due to its negligible contribution in comparison with the Py thin film shape anisotropy term (which is two orders of magnitude larger). In our case, as we are going to simulate a rather thick magnetic heterostructure, it needs to be taken into account. The Py disc was 290 nm in diameter with a thickness of 40 nm. It presents a clear vortex configuration at remanence as represented from the central slice of the simulation (Fig. 2c[link]). The colour scale for the magnetization orientation is represented by small vortices with opposite polarities (black → mz negative; white → mz positive).

The micro-particle consists of a 1 µm-side square with 100 nm/40 nm/40 nm thickness structure of Py/Al/Py. Different slices of the simulation showing the magnetization configurations at remanence are presented for the thick Py layer (I–III, Fig. 2d[link]) and for the thin one (IV, Fig. 2d[link]). The system supports stripe domains in the thick Py layer and mainly in-plane magnetization is present in the thin one. This decoupling is mediated by the non-magnetic Al spacer. The output of the micromagnetic simulation encloses the reduced magnetization vector at each simulation cell. This is combined with a model containing the attenuation length of each material for the Fe L3 energy edge (Py → L−1 = 5.13 × 106 m−1, Al → L−1 = 8.24 × 105 m−1), and a dichroic factor of 0.22 only where the Py is present. By projecting these models using equation (4b)[link], we simulate the X-ray transmission tomography data which will be used to test the reconstruction algorithm.

In order to completely validate the reconstruction method, we investigate tomograms with data from the full projection (FP) range (−90 to 90°, 1° step), and tilt series with a limited number of projections (missing wedge, MW). We have chosen for the latter a range from −60 to 60° with 1° step. This limitation is typical of X-ray transmission tomography setups (Pereiro et al., 2009[Pereiro, E., Nicolás, J., Ferrer, S. & Howells, M. R. (2009). J. Synchrotron Rad. 16, 505-512.]; Olivares-Marín et al., 2015[Olivares-Marín, M., Sorrentino, A., Lee, R. C., Pereiro, E., Wu, N. L. & Tonti, D. (2015). Nano Lett. 15, 6932-6938.]). The nano-disc and micro-square particles have been reconstructed taking 64 and 46 Z layers in the volume model, respectively. Fig. 3[link] presents a comparison between ground truth and the reconstructed data for the nano-disc particle in the 32nd Z layer of the volume. Components X, Y and Z of the reconstructed vector field are shown for the FP and MW situations. Moreover, in order to quantify the quality of the reconstructions, the normalized root mean square error (NRMSE) of the reconstruction compared with the ground truth is presented. This parameter is calculated as follows,

[{\rm{NRMSE}} = {{ 1 }\over{ X_{\max}-X_{\min} }} \left[ {{1}\over{N}}\textstyle\sum\limits_{i\,=\,1}^N \left(X_{\rm{GT}}-X\right)^2\right]^{1/2}. \eqno(11)]

X represents the data contained in the reconstructed model slice to be analysed; XGT is the ground truth data; Xmin and Xmax are the minimum and maximum values of the pixels in the analysed slice, respectively; and N is the total number of pixels involved. The high quality of the reconstructions is directly observed for the MW and FP situations. However, the latter presents an excellent agreement with the ground truth data, while the value of all vector components in the MW case is smaller than the original one (5–20%).

[Figure 3]
Figure 3
Ground truth [mx (a), my (b), mz (c)], full projection [mx (d), my (e), mz (f)] and missing wedge [mx (g), my (h), mz (i)] reconstructions of the disc particle. Images correspond to Z layer 32 of the volume model. The normalized root mean square error is indicated for the specific images.

We have analysed also the NRMSE for all the Z slices of the reconstruction volume to observe the MW effects along the thickness in the reconstructed solution (Fig. 4[link]). Vertical dashed lines have been superimposed in the graphs where the edges of the disc are present. The results for the estimated error of X (Fig. 4a[link]) and Y (Fig. 4b[link]) components are almost equal, while that for Z (Fig. 4a[link]) presents a different behaviour. This occurs because the information for the reconstruction of the Z component is present in the Xtilt and Ytilt series, while the other components are reconstructed from the information of individual series. The evolution of the NRMSE calculated for the whole reconstructed volume model instead of by the Z layer is presented for the FP (Fig. 3d[link]) and MW (Fig. 3e[link]) situations. The MW case presents in general a larger error than the full projection one. The reconstruction of the nano-disc top and bottom surfaces is also affected by the MW configuration leading to an ambiguity at the borders. It is important to mention here that the increased error associated with the MW series is mainly associated with the X and Y components. The out-of-plane component presents almost the same error in the FP (NRMSE = 0.02) and MW (NRMSE = 0.03) situations; the only difference here is a small oscillation in the NRMSE versus iterations. In any case, both reconstructions clearly allow the disc structure at the centre of the model volume and its magnetic configuration to be identified.

[Figure 4]
Figure 4
Normalized root mean square error from the reconstructed mx (a), my (b) and mz (c) of the nano-disc as a function of the Z layer. Full projection (FP, black squares) and missing wedge (MW, red dots) situations are analysed. Dashed vertical lines indicate where the disc is located. NRMSE for the whole volume model as a function of iterations in the FP (d) and MW (e) situations. Components mx (black squares) and my (red dots) present almost the same behaviour and appear superimposed.

In the case of the magnetic micro-square heterostructure, a direct comparison between ground truth and reconstructions is presented in Figs. 5[link] and 6[link]. The first one shows the 16th Z slice (thick Py region) while the second represents the 38th Z slice (thin Py region).

[Figure 5]
Figure 5
Ground truth [mx (a), my (b), mz (c)], full projection [mx (d), my (e), mz (f)] and missing wedge [mx (g), my (h), mz (i)] reconstructions of the square particle. Images correspond to Z layer 16 of the volume model. NRMSE is indicated for the specific images.
[Figure 6]
Figure 6
Ground truth [mx (a), my (b), mz (c)], full projection [mx (d), my (e), mz (f)] and missing wedge [mx (g), my (h), mz (i)] reconstructions of the square particle. Images correspond to Z layer 38 of the volume model. NRMSE is indicated for the specific images.

Again, the FP and MW situations are studied. The latter presents a smaller intensity in all the reconstructed components and the agreement of the FP range reconstruction is much better than the MW one. It is important to mention here that the NRMSE for the Z component of the vector presents almost the same value for the MW and FP reconstructions in the out-of-plane dominated area of the magnetic structure.

The aforementioned effect in the error is clearly observed in the NRMSE representation as a function of the Z slice of the volume model for all the vector components (Fig. 7[link]). The Z component presents a maximum in the error for the region where the out-of-plane magnetization dominates (thick Py region), and in the range dominated by mx and my the NRMSE decreases. This implies that it is more difficult to reconstruct the Z component despite the redundancy in the tomographic data (Z component reconstructed from both tilt series). This can also be observed in the convergence plots showing the NRMSE calculated for the whole reconstructed model as a function of the iteration number (Figs. 7d and 7e[link]). Again, the main difference between MW and FP is an important increase in the error associated with the in-plane components while the value for the Z component is almost the same. It is also observed that the MW induces artifacts at the interfaces inside the object between different layers of the heterostructure. These artifacts are due to the missing data and conventional tomography also experiences such artefacts, especially for in situ imaging studies where the sample is in a chamber or cell and not all views are accessible with X-rays. However, the trend of the reconstructed vector field is qualitatively and in some cases quantitatively in agreement with the original magnetization configuration, indicating that the algorithm is capable of successfully reconstructing the 3D magnetization from magnetic soft X-ray transmission tomograms.

[Figure 7]
Figure 7
Normalized root mean square error from the reconstructed mx (a), my (b) and mz (c) of the micro-square as a function of the Z layer. Full projection (FP, black squares) and missing wedge (MW, red dots) situations are analysed. Dashed vertical lines indicate where different materials are located. NRMSE for the whole volume model as a function of iterations in the FP (d) and MW (e) situations. Components mx (black squares) and my (red dots) present almost the same behaviour and appear superimposed.

5. Conclusions

A new method to reconstruct the magnetization vector field of arbitrary magnetic systems using soft X-ray transmission tomography has been described. The method takes advantage of the natural high dichroic contrast of magnetic materials at soft X-ray energies which leads in practice to acquisition times of only a few hours to achieve expected resolutions of around 40 nm or better. The technique is useful to characterize magnetic samples with thicknesses up to ∼300 nm and up to several micrometres of lateral dimensions. Both scalar and vector reconstruction problems have been analysed in detail and solved by using ART. The vector case requires two differently oriented tilt series to obtain the three components of the magnetization. To test the method, two different magnetic particles have been simulated, and their respective tomograms calculated. We have studied both full projections and also incomplete series due to missing wedges to mimic actual experimental limitations. The results for the full projections are always better than those for the missing wedge as expected; however, both approaches provide qualitative and even quantitative descriptions of the magnetic structures. The method is well suited for providing detailed information of the magnetization of buried magnetic structures or interfaces, and consequently appears to be a valid characterization technique of 3D magnetism in spintronic devices.

Footnotes

Present address: SUPA, University of Glasgow, School of Physics and Astronomy, G12 8QQ Glasgow, UK.

Funding information

The following funding is acknowledged: Spanish MINECO (grant No. FIS2013-45469; grant No. FIS2016-76058 (AEI/FEDER, EU); contract No. FIS2016-76058 (AEI/FEDER, EU) to AHR); FICYT-Asturias (grant No. FC-GRUPIN14-040); US Department of Energy (DOE) Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory (contract No. DE-AC02-06CH11357 to DG); US Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division (contract to CP). AHR acknowledges support from European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Action (reference H2020-MSCA-IF-2016-746958).

References

First citationBerger, A. & Pufall, M. R. (1999). J. Appl. Phys. 85, 4583–4585.  Web of Science CrossRef Google Scholar
First citationBlanco-Roldán, C., Quirós, C., Sorrentino, A., Hierro-Rodríguez, A., Álvarez-Prado, L. M., Valcárcel, R., Duch, M., Torras, N., Esteve, J., Martín, J. I., Vélez, M., Alameda, J. M., Pereiro, E. & Ferrer, S. (2015). Nat. Commun. 6, 8196.  Google Scholar
First citationDe Carlo, F., Gürsoy, D., Ching, D. J., Batenburg, J., Ludwing, W., Mancini, L., Marone, F., Mokso, R., Pelt, D., Sijbers, J. & Rivers, M. (2018). Meas. Sci. Technol. 29, 034004.  Web of Science CrossRef Google Scholar
First citationDe Graef, M. (2001). Magnetic Imaging and its Application to Materials, edited by M. De Graef & Y. Zhu, Vol. 36, ch. 2. San Diego: Academic Press.  Google Scholar
First citationDonnelly, C., Guizar-Sicairos, M., Scagnoli, V., Gliga, S., Holler, M., Raabe, J. & Heyderman, L. J. (2017). Nature (London), 547, 328–331.  Web of Science CrossRef Google Scholar
First citationFernández-Pacheco, A., Streubel, R., Fruchart, O., Hertel, R., Fischer, P. & Cowburn, R. P. (2017). Nat. Commun. 8, 15756.  Google Scholar
First citationFoerster, M., Prat, J., Massana, V., Gonzalez, N., Fontsere, A., Molas, B., Matilla, O., Pellegrin, E. & Aballe, L. (2016). Ultramicroscopy, 171, 63–69.  Web of Science CrossRef Google Scholar
First citationGürsoy, D., De Carlo, F., Xiao, X. & Jacobsen, C. (2014). J. Synchrotron Rad. 21, 1188–1193.  Web of Science CrossRef IUCr Journals Google Scholar
First citationHierro-Rodriguez, A., Quirós, C., Sorrentino, A., Blanco-Roldán, C., Álvarez-Prado, L. M., Martín, J. I., Alameda, J. M., Pereiro, E., Vélez, M. & Ferrer, S. (2017a). Phys. Rev. B, 95, 014430.  Google Scholar
First citationHierro-Rodriguez, A., Quirós, C., Sorrentino, A., Valcárcel, R., Estébanez, I., Álvarez-Prado, L. M., Martín, J. I., Alameda, J. M., Pereiro, E., Vélez, M. & Ferrer, S. (2017b). Appl. Phys. Lett. 110, 262402.  Google Scholar
First citationHubert, A. & Schaëfer, R. (1998). Magnetic Domains: The Analysis of Magnetic Microstructures. Berlin/Heidelberg: Springer.  Google Scholar
First citationJacobs, F., Sundermann, E., De Sutter, B., Christiaens, M. & Lemahieu, I. (1998). J. Comput. Information Technol. 6, 89–94.  Google Scholar
First citationKak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographic Imaging. New York: IEEE Press.  Google Scholar
First citationLocatelli, A. & Bauer, E. (2008). J. Phys. Condens. Matter, 20, 093002.  Web of Science CrossRef Google Scholar
First citationMcVitie, S., McGrouther, D., McFadzean, S., MacLaren, D. A., O'Shea, J. J. & Benitez, M. J. (2015). Ultramicroscopy, 152, 57–62.  Web of Science CrossRef Google Scholar
First citationOlivares-Marín, M., Sorrentino, A., Lee, R. C., Pereiro, E., Wu, N. L. & Tonti, D. (2015). Nano Lett. 15, 6932–6938.  Google Scholar
First citationPereiro, E., Nicolás, J., Ferrer, S. & Howells, M. R. (2009). J. Synchrotron Rad. 16, 505–512.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationPhatak, C., Beleggia, M. & De Graef, M. (2008). Ultramicroscopy, 108, 503–513.  Web of Science CrossRef Google Scholar
First citationPrabhat, K. C., Aditya Mohan, K., Phatak, C., Bouman, C. & De Graef, M. (2017). Ultramicroscopy, 182, 131–144.  Web of Science CrossRef Google Scholar
First citationScheinfein, M. R., Unguris, J., Kelley, M. H., Pierce, D. T. & Celotta, R. J. (1990). Rev. Sci. Instrum. 61, 2501–2527.  CrossRef Web of Science Google Scholar
First citationSiddon, R. L. (1985). Med. Phys. 12, 252–255.  CrossRef CAS PubMed Web of Science Google Scholar
First citationSorrentino, A., Nicolás, J., Valcárcel, R., Chichón, F. J., Rosanes, M., Avila, J., Tkachuk, A., Irwin, J., Ferrer, S. & Pereiro, E. (2015). J. Synchrotron Rad. 22, 1112–1117.  Web of Science CrossRef IUCr Journals Google Scholar
First citationTeixeira, J. M., Lusche, R., Ventura, J., Fermento, R., Carpinteiro, F., Araujo, J. P., Sousa, J. B., Cardoso, S. & Freitas, P. P. (2011). Rev. Sci. Instrum. 82, 043209.  Web of Science CrossRef Google Scholar
First citationUnguris, J. (2001). Magnetic Imaging and its Application to Materials, edited by M. De Graef & Y. Zhu, Vol. 36, pp. 167–194. San Diego: Academic Press.  Google Scholar
First citationVansteenkiste, A., Leliaert, J., Dvornik, M., Helsen, M., Garcia-Sanchez, F. & Waeyenberge, B. V. (2014). AIP Adv. 4, 107133.  Google Scholar
First citationVoltan, S., Cirillo, C., Snijders, H. J., Lahabi, K., García-Santiago, A., Hernández, J. M., Attanasio, C. & Aarts, J. (2016). Phys. Rev. B, 94, 094406.  Web of Science CrossRef Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds