Diffractive small angle X-ray scattering imaging for anisotropic structures

Insights into the micro- and nano-architecture of materials is crucial for understanding and predicting their macroscopic behaviour. In particular, for emerging applications such as meta-materials, the micrometer scale becomes highly relevant. The micro-architecture of such materials can be tailored to exhibit specific mechanical, optical or electromagnetic behaviours. Consequently, quality control at micrometer scale must be guaranteed over extended areas. Mesoscale investigations over millimetre sized areas can be performed by scanning small angle X-ray scattering methods (SAXS). However, due to their long measurement times, real time or operando investigations are hindered. Here we present a method based on X-ray diffractive optics that enables the acquisition of SAXS signals in a single shot (few milliseconds) over extended areas. This method is applicable to a wide range of X-ray sources with varying levels of spatial coherence and monochromaticity, as demonstrated from the experimental results. This enables a scalable solution of spatially resolved SAXS.

T he macroscopic mechanical behaviour of materials is directly linked to the arrangement of the individual building elements at a micro or nanoscopic level. Examples of such materials occurring in nature include bone and wood. In both cases the arrangement of fibres provide the mechanical stability or compliance of the structure 1,2 . Bioinspired artificial structures [3][4][5] with superior properties such as ultra low weight to strength ratio can be utilised in a wide range of fields such as biomedical engineering 6,7 , aerospace engineering 8 and the automotive industry 9 . Such structures can be realised with technologies like additive manufacturing and injection moulding where fibrous structures are either mimicked or directly introduced in the raw material. The further development and optimisation of such structures requires accurate understanding of the local fibre orientation which eventually determines the macroscopic mechanical properties. Current technologies capable of probing the fibre orientation at a microscopic level include visible light microscopy, high resolution X-ray computed tomography 10 and scanning small angle X-ray scattering (SAXS) 11 .
In particular the X-ray based approaches have the benefit to penetrate matter and therefore can image inner or concealed areas. However, both high resolution X-ray computed tomography and spatially resolved SAXS are limited to sample sizes of a few millimetres. In addition, the scan time even with highly brilliant synchrotron beams is in the range of tens of hours 12 , therefore hindering in situ and eventually forbidding real-time or operando investigations. SAXS signals can be partially recorded in the real space over larger areas (centimetre range) by utilising interferometric methods such as X-ray grating interferometry (XGI) 13 . XGI utilises periodic phase modulating structures combined with coherent illumination in order to create an interference fringe at specific distances downstream (Talbot effect). Scattering introduced by the sample causes broadening of the incoming beam and consecutively reduction of the contrast of the interference fringe. The relationship between the reduction in visibility and the small angle scattering properties of the sample has been well developed and understood [14][15][16][17] . Conventionally, the interference pattern is limited to only one direction, meaning that only a radial component of the scattering intensity is sensed. For isotropically scattering samples this is not a limitation, however, for the investigation of samples with anisotropic scattering properties rotation of either the sample or the imaging modality is necessary 18,19 , leading to prolonged acquisition times and thus inhibiting the observation of time evolving effects. Recent developments in XGI have allowed an omnidirectional sensitivity of the scattering signal over the imaging plane in a single shot 20 . The method relies on the direct recording of the interference fringe of an array of circular gratings, each creating a self-image at a specific distance downstream the grating which is scaling with the square of the grating period. In order to observe this effect, sufficient coherence, approximately equal to the grating period at the grating, and high enough resolution, at least half of the interference fringe period, are required. Taking into account the low refraction angles of X-rays we can conclude that the existing method is only applicable on highly coherent sources, such as synchrotrons, in combination with micrometre resolution detectors. Therefore, the applicability of the method is constrained to at best millimetre sized samples that can only be investigated at large scale facilities.
Here we propose a general methodology that provides omnidirectional scattering sensitivity, does not stringently require coherent or monochromatic illumination, and finally does not rely on high resolution detectors. The key advancement is the introduction of an optical element capable of creating an illumination field, by incoherent superposition, that locally encodes the scattering properties with omnidirectional sensitivity. We demonstrate the compatibility of the method with synchrotron sources and conventional micro-and macrofocal X-ray tubes. Enabling real-time investigations and examinations of larger objects. Finally, the method is applied to the investigation of the fibre orientations in reinforced polymers manufactured by injection moulding.

Results
Demonstration of working principle. The optical element is composed of unit cells that cover the full field of view (FoV) to be examined. Each unit cell has a width of W and is an ensemble of equally spaced concentric annuli with a period P, with each annulus containing a periodic diffractive substructure with a period p 1 significantly smaller than the width of the annuli. A schematic representation of a unit cell with two diffractive annuli is shown in Fig. 1a. The width of the annulus can be chosen freely, however, for simplicity we assume that it is half of the coarse period P. The diffraction from the substructure results in a conical beam splitting. Provided sufficient propagation distance after the optical element, a modulation pattern with a circular symmetry will start to appear as illustrated in Fig. 1b for the case of a single unit cell with only one diffractive annulus. For clarity, only a cross section of the ± 1 diffraction orders is shown while neglecting the zeroth order. The projected unit cell size W 0 defines the pixel size in the retrieved images and is refereed to as a macropixel. For a period of p 1 of the substructure, the diffraction angle of the ± 1 orders for perpendicularly incident radiation of wavelength λ is ± sin À1 λ p 1 . In the hard X-ray regime where λ $ 10 À11 m and for micrometre sized structures p 1 $ 10 À6 m the diffraction angles can be approximated by ± λ p 1 . The distance at which the ± 1 orders have separated completely and overlap with the transmitted X-rays from the non diffractive areas, and subsequently the intensity of the modulation maximises, is given by D ¼ 1 2 Pp 1 λ for a parallel beam geometry. A theoretical derivation of the conditions for maximising the visibility is presented in the Methods section. If a sample that exhibits small angle scattering is placed before the optical element, the visibility of the recorded pattern will be reduced. In the Methods section, it is analytically derived that the visibility reduction can be extracted under all angles on the imaging plane and corresponds to a ring of the real space correlation of the sample at a length scale ξ ¼ zλ P . Therefore, this type of measurement has the capability to extract the orientation of the underlying scattering structures. For full SAXS analysis measurements at additional correlation lengths would be required. In Fig. 1c a scanning electron microscopy (SEM) image of fabricated optics containing a single diffractive annulus with a fine diffractive period of 1 μm, coarse period of 35.75 μm, and design energy of 17 keV are shown. The measured average visibility illustrated in Fig. 1d at 17 keV for distances between 16 and 75 cm between the optical element and the detector, shows a clear maximum of 52% at the theoretical position. The insets depict the intensity pattern of a single unit cell at the marked distances.
Imaging of injection moulded fibre reinforced polymers. To demonstrate the capability of the method to extract the fibre orientation in relevant samples we imaged a glass fibre reinforced injection moulded tensile sample shown in Fig. 2a. The injection procedure took place from the two inlets annotated with red circles, thus creating a weld line near the centre of the sample. In Fig. 2b the retrieved orientation and degree of anisotropy of the fibres are shown. The colour corresponds to the orientation according to the colour wheel, and the degree of anisotropy to the intensity. The experimental details and retrieval method are given in the Methods section. In Fig. 2c the corresponding orientation The details of the simulation are given in the Methods. In general, there is a strong agreement between the orientation predicted from the proposed method and the state of the art simulation method. However, the most prominent and crucial difference can be found at the weld line. The simulation fails to predict the accurate location of the weld line, potentially due to the assumption that the injection is happening with exactly the same pressure from both inlets. In addition, our measurements show an asymmetry of the weld line in the horizontal axis, which the simulation failed to predict.
We further investigated the weld line area by comparing the additional contrasts that are available through the proposed method. In the conventional absorption image in Fig. 3a no visible sign of the weld line or any micro-structural difference can be inferred since a uniform absorption signal is observed. On the other hand, looking at the average scattering image shown in Fig. 3b, it becomes evident that some non uniformity is present. Finally, clear evidence is provided by extracting the degree of anisotropy image shown in Fig. 3c. The weld line becomes highly visible together with the neighbouring affected zones, where a significantly lower fibre alignment is inferred. In Fig. 3d the local fibre orientation and degree of anisotropy is plotted with a vector plot, the direction corresponds to the estimated fibre direction and the colour to the degree of anisotropy.
In order to further support the observed difference in the weld line we performed a micrometre computed tomography (μCT) scan of the affected area. The Experimental details of the μCT are given in the Methods section. We chose to show three slices, one in the xy plane going through the middle of the sample at the location of the weld line, and two in the yz plane, one directly at the weld line and one 1 cm further down in the unaffected area. In Fig. 3e the illustrated tomographic slices are marked on a sketch of the tensile sample. The sagittal slice from the unaffected area shows the fibres to be oriented along the x direction (perpendicular to the observation plane). But as we move the observation plane closer to the weld line an increasing number of in plane fibres start to show up. Finally, once the centre of the weld line has been reached all fibres seem to be located in plane as shown in the second sagittal plane in Fig. 3g. From the coronal slice shown in Fig. 3f a deeper understanding of the observed degree of anisotropy shown in Fig. 3c can be accomplished. The zoom-ins correspond to areas of high and low degree or anisotropy marked by the blue and green boxes, respectively. As expected the blue box shows fibres highly aligned while the green fibres more randomly distributed. The observed asymmetry of the weld line can also be justified from the coronal slice, as a higher amount of aligned fibres can be observed at the bottom of the weld line. Finally, the exact apparent shape of the weld line might not match completely the one extracted from the slice due to the fact that we are looking at a projection and not a tomographic reconstruction of the fibre orientation.
Demonstration of real time fibre orientation extraction. The combination of transparent optics, single shot acquisition mode, and high flux beams from synchrotron sources allow the imaging procedure to be executed in a few milliseconds, thus enabling real-time imaging. As a demonstrator we chose to image a carbon fibre overhand knot in real time (25 frames per second) while it was being tightened. The tensile experiment is sketched in Fig. 4a. The one end of the knot was pulled while the other was anchored to a fixed point. In Fig. 4b and c, the retrieved orientation maps of the carbon fibre knot at time points 0.0 and 11.0 s are shown. The lines correspond to the retrieved direction of the carbon fibres and the colour to the normalised degree of anisotropy (by 0.27) translating in this case directly to the local fibre density variations 18 . In the Supplementary Movie 1 the orientation of the underlying fibres and the local fibre density of the knot throughout the whole experiment are shown. Since only one end of the knot is being pulled a sliding motion in the direction of pull is observed. This causes the free part of the knot to develop an asymmetric pattern, with the left side containing fibres more densely packed than the right. This is observed by the increase in the degree of anisotropy. The intertwined part of the knot shows a stronger signal due to the overlapping fibres of the strands.
Imaging on conventional X-ray tubes. For the translation of the method from synchrotron facilities to laboratory X-ray sources two major challenges need to be considered: source size and the polychromatic spectrum of such sources. Intrinsically, the method does not rely on coherent illumination, however, the projected source size should be small enough in order not to completely blur out the generated pattern. Practically, this can be avoided if the projected source size is kept smaller than half of the projected pattern's period P 0 , leading to the condition: where L is the source to detector distance and σ the width of the X-ray source. Given that condition (1) is not usually fulfilled, an aperture array can be introduced in front of the source to mitigate this effect. The opening of the apertures should be small enough in order to fulfil condition (1), and their spacing such that a superimposition of the individual generated patterns corresponding to each aperture is achieved (Lau condition) 21 . In Supplementary Note 1, we present proof of concept images recorded with a macrofocal X-ray tube (1 mm) in combination with an aperture array. In addition, discussion regarding the flux efficiency of such an approach is presented. Regarding the polychromaticity it is demonstrated in the Methods section that theoretically a spectral acceptance of 90% can be achieved, therefore confirming the achromatic character of the method and its compatibility with a wide range of X-ray sources. In Fig. 5a, a recorded intensity modulation with a visibility of 20% is shown resulting from a microfocal tube operated at 70 kVp. As a demonstrator we imaged a carbon fibre loop with predictable fibre orientation. As shown in Fig. 5b, a reliable reconstruction of the fibre orientation was achieved. Fig. 5c presents the fibre orientation map of a carbon fibre reinforced injection moulded component and well illustrates the application of our method to an industrially relevant sample. Four injection locations were used, indicated by the round head arrows. The top part of the sample was not imaged due to geometric constrains of our setup. The pointy arrows mark the resulting weld lines.

Discussion
In this work we introduced a generalised imaging method for obtaining scattering signals with an omnidirectional sensitivity at a defined correlation length. Due to the lack of any strong coherence or monochromaticity requirements the method is compatible with a wide range of X-ray sources and X-ray detectors as demonstrated in the experimental section. The presented experiments demonstrate the potential of the method to address current challenges in the field of fibre reinforced structures where both fibre orientation and density play a crucial role in the observed macroscopic mechanical properties. It could do so by being integrated in the development, quality control, and or testing cycles. Thus overcoming the limitation of the current state of the art simulative methodology. Nevertheless, the applicability of the method is not only limited to fibre reinforced polymer structures, but can be extended to biological specimens such as bone 22 or brain tissue 23 . On a more speculative basis, the compatibility with high power and energy X-ray tubes opens up the possibility to image hard condensed systems such as alloys. Potentially, lamellar orientation could be extracted which again is crucial for the mechanical properties 24 of the alloys. On the soft matter side, applications could include imaging of three dimensional printed structures 25 such as biologically inspired composites, shape-morphing systems, and soft sensors.
The length scale sensitivity of our method is equal to half of the period of the diffractive structure; with the current fabrication methods of X-ray optics periods spanning between a few micrometres down to hundreds of nanometres 26 can be implemented. In addition, the method optimally (in terms of fringe visibility and photon utilisation) requires phase shifting gratings, by utilising heavier materials such as Ir or Au optics for high energies can be fabricated and thus targeting either samples of industrial scale or hard condensed matter.
The degree of coherence of the X-ray sources plays an integral role in determining the final fringe visibility since geometric blurring needs to be taken into account. Therefore, highly coherent sources will allow for an increased flat visibility. The noise level of the retrieved scattering signal is directly linked to the nominal visibility 27 . Thus, we foresee that an increase in coherence, which is expected with the upcoming machine upgrades of light sources, can benefit the method in two ways: either by allowing to reduce the acquisition time in order to observer faster dynamic effects, or allow for an increase in sensitivity for weakly scattering materials such as biological samples.
The main limitation of the current method is the fixed autocorrelation length. Meaning that in oder to ripe the benefits of full SAXS quantification methodologies 28 additional measurements at different autocorrelation lengths are required. In order to avoid this, a different unit cell design could be implemented allowing multiple autocorrelation lengths in a single shot (M. Kagias, Z. Wang, G. Lovric, K. Jefimovs and M. Stampanoni, manuscript in preparation).
The single shot aspect of the method unlocks the potential for either time resolved studies with millisecond temporal resolution over centimetre sized areas, previously impossible with conventional SAXS methods, or accelerated tensor tomographic imaging with acquisition times in the range of a few seconds compared with hours with the current methods. Finally, the method can also be extended to other forms of radiation from low coherence and brilliance sources such as neutrons.

Methods
Visibility calculation and design optimisation. We perform a theoretical analysis of the visibility estimation assuming that the diffractive substructure is composed of a binary grating with complex transmissions of T h and T l , duty cycle d c ¼ with p h and p l being the lengths of the two states, respectively, and p 1 ¼ p h þ p l . Our goal is to extract optimum parameters for maximising the fringe visibility. The diffraction efficiencies can be calculated from the coefficients of a Fourier series expansion of the diffractive substructure. The efficiencies for the jnj ! 1 orders are given by The transmitted order is given from For simplicity we will consider only a radial cross section of a single unit cell of the optical element. After a propagation distance z the recorded intensity IðxÞ can be considered as the sum of the non diffracted I 0 ðxÞ and diffracted I n ðxÞ orders n¼ ± 1; ± 2;::: The visibility is estimated as twice the ratio of the first and zeroth Fourier components of the recorded intensity. Due to the linearity of Fourier series, we can calculate the zeroth and first Fourier components for each diffraction order separately, then perform the addition in the Fourier space. Each intensity I n ðxÞ can be considered as a pulse train. The zeroth order pulse train will have levels of 1 and A 0 , in the case of a generic order n, the pulse train will have levels of 0 and A n a b c 0 1 Normalised degree of anisotropy z. Taking into account the above assumption we conclude that the zeroth order Fourier component will be given by And the first by n¼1;2;::: Resulting in a visibility of VðzÞ ¼ 4 π 1 À A 0 À 2 P 1 n¼1;2;::: A n cosð2πn 2 λ Our goal is to estimate the distance z at which the visibility will maximise. The positions where the visibility maximises or minimises are given from the solutions of the following equation which can be rewritten as X 1 n¼1;2;::: A set of solutions of Eq. (9) can be given when the sinusoidal term is equal to 0 for every n. This can be achieved if In order to find which of these solutions represent maxima we calculate the second order partial derivative of the visibility which should be positive which results in the following set of locations where the visibility maximises in the case of monochromatic radiation These solutions are general and do not depend on the type (phase, attenuation or mixed) of the diffractive structure used but only on the wavelength and periods. Further optimisation can be performed at these locations to maximise the observed visibility. Specifically, we want to optimise for the duty cycle and transmission levels of the diffractive substructures. For simplicity we will focus on the case where the diffractive structure is only inducing a phase shift of ΔΦ. In this case the diffraction efficiencies of formulas (2) and (3) become and In addition, the denominator of the second fraction of the visibility term according to Eq. (7) will converge to 2 since no intensity is lost in the grating structure. From the above we can straightforwardly conclude that for pure phase shifting optics the visibility maximises when the phase shift ΔΦ is an odd multiple of π and the duty cycle is 0.5. In the case of mixed gratings optimum visibility is expected for different duty cycle and phase/transmission values.
Spectral acceptance. In the previous section we derived the conditions for maximising the visibility in the case of monochromatic radiation. Here we will demonstrate that fringe formation process still takes place with a broad spectrum. For the sake of simplicity we will assume a purely phase shifting element, with a duty cycle of 0.5, with energy dependent refractive index decrement of δðλÞ, designed to introduce a phase shift of π at a given wavelength λ d or energy E d . In addition, we will take into account only the ± 1 diffraction orders. At the nominal design distance of z 1 ¼ 1 2 Pp 1 λ d from Eqs. (7), (13), and (14) the visibility becomes Under the assumption that the refractive index decrement is inversely proportional to the square of the photon energy δðEÞ / 1=E 2 (which holds true in the absence of absorption edges etc.) Eq. (15) becomes The spectral acceptance ΔE=E d is calculated as the full width half maximum (FWHM) of the lobe occurring around the design energy E d , which is calculated to be~90% from Eq. (16). Indicating that our method is highly compatible with polychromatic sources. Nonetheless, the spectral acceptance rapidly reduces once higher order design distances are considered. For example, for m ¼ 3 and m ¼ 5 the spectral acceptance becomes ΔE=E d % 34% and ΔE=E d % 17%, respectively.
Modelling of interaction between sample and optical element. Let us assume that the sample exhibits uniform scattering properties within the area of a single unit cell. This scattering behaviour is described by the symmetric two dimensional scattering function Sðq x ; q y Þ with q x and q y being the scattering vectors in the x and y directions, respectively. For convenience we will utilise the polar representation of the scattering function Sðq; θÞ. If the sample is placed right before the optical element and scatters with a vector q under angle θ, the diffraction orders of the radial cross section of the optical element under angle ψ will propagate with an angle given by By following similar calculations as in the previous subsection it can be shown that the visibility at angle ϕ will be scaled by a factor cos ξq cosðϕ À θÞ ½ ; ð18Þ where ξ ¼ zλ P . To retrieve the total visibility reduction under angle ϕ noted by V r ðξ; ϕÞ we need to integrate over all scattering vectors By changing to the cartesian coordinates q x , q y , ξ x and ξ y the above integral can be written as which implies that the visibility reduction is the Fourier cosine transform of the scattering function. Nonetheless, it is known that the real space correlation function Gðξ x ; ξ y Þ and the scattering function form a pair of the Fourier cosine transform therefore it can be deducted that the visibility reduction corresponds to the real space correlation of the sample. The distance z at which the detector is placed defines the length scales that are being probed. In the case the measurement is carried out at the optimum distance described in the previous subsection, we obtain a ring from the real space correlation function at a correlation length equal to half of the period of the fine diffractive structure.
Retrieval of the directional visibility reduction. The recorded pattern for each unit cell in the absence of a sample can be approximated as a radial cosine where A f corresponds to the average intensity of the flat pattern and V f ðξ; ϕÞ the radial flat visibility. Under the assumption that the investigated sample has uniform scattering properties in the range of one unitcell, the recorded intensity at each unit cell once the sample has been introduced will be For each recorded unit cell a two dimensional discrete Fourier transform is performed. The transmission of the sample for each unit cell is calculated as the ratio of the zeroth orders of the flat and sample measurements T m;n ¼Î s ðq m;n ¼ 0Þ=Î f ðq m;n ¼ 0Þ. The visibility reduction at each angle ϕ is then given by where N is the ratio of the unit cell size W and the coarse period P. The angularly dependent visibility reduction can be modelled as following: where α 0 corresponds to the average scattering intensity, α 1 the oriented component of the scattering, and ϕ d the direction of the underlying structure. The ratio α 1 =α 0 is the degree of anisotropy. These three parameters can be extracted by Fourier analysis. transferred to the Cr mask by reactive ion etching of Cl and CO 2 plasm. Finally, the Si was etched to the specified depth depending on the design photon energy.
Omnidirectional imaging on synchrotron X-ray sources. In both cases monochromatic X-rays, selected by a multilayer monochromator (2% bandwidth) were used. After passing through the sample and the optical element the X-rays were converted to visible light by a 300-μm thick scintillator and finally recorded with a pco.edge 5.5 CMOS camera equipped a pixel size of 6.5 μm. For the tensile sample, the energy was set to 25 keV. The optical element had a unit cell size of 71.5 μm, fine period of 1.2 μm, and a coarse period of 35.75 μm. Each unit cell contained only one annulus. Due to the parallel beam geometry the macropixel size of the retrieved images correspond to the unit cell size of 71.5 μm. The field of view was 1:5 0:4 cm 2 and was limited by the beam size. For this reason multiple exposures were performed in order to image the full sample. For the time resolved experiment, optics designed for 17 keV, with a unit cell size of 84.5 μm, fine and coarse periods of 1 and 42.25 μm, respectively. Again due to the parallel beam geometry a final macropixel size of 84.5 μm was achieved in the retrieved images. The exposure time for each frame was set to 40 ms giving a frame rate of 25 frames per second which is considered real time.
Omnidirectional imaging on X-ray tubes. The carbon fibre loop was measured with the HAMAMATSU L10101 X-ray tube, operated at 70 kVp and 100 μA, resulting in a source size of 10 μm. The optical element had a diffractive period p 1 ¼ 1 μm, annular period P ¼ 35:75 μm and unit cell size W ¼ 71:5 μm. The total system had a length of 1 m with the optical element placed in the middle. The sample was placed right before the optical element. Due to the symmetric configuration the macropixel size was also 71:5 μm. The detector in use was the PI-SCX:4300 equipped with a (Gadox) scintillator and a pixel size of 24 μm. The exposure time for the flat and sample images was 60 sec each. For the carbon fibre reinforced component a similar configuration was used. Instead of a phase shifting optical element we utilised an absorbing structure in order to mitigate beam hardening effects. A mosaic scan was performed to cover the whole structure, a total of 47 images where acquired with at least 50% overlap and exposure time of 10 min each.
Simulation of fibre orientation in injection moulded sample. The injection moulding simulation was performed in Moldex3D 29 and the resulting fibre orientation tensors were exported in a structural finite element mesh for structural simulations in Ansys Workbench. The simulations were performed in a grid with approximately 1 1 0:5 mm 3 element sizes. The injection simulation was performed in PA66 GF30 glass-filled polymer with a filling time of 0.6 s, melt temperature of 295 C, mould temperature of 100 C and a maximum filling pressure of 500 MPa. The resulting 3D orientation tensor information was projected onto ellipses in the xy-plane and averaged over the z-direction. The orientation and the degree of orientation for the fibre orientation plot were calculated as the direction of the major-axis and the ratio between major-and minor-axis, respectively.
Micro CT of weld line. The measurement was done using a ZEISS Xradia 520 Versa X-ray microscope with a Flat Panel (FPX) detector. The polychromatic micro-focus source was operated at 80 kV, 7 W and with a low energy filter (LE2) to optimise transmission. In total, 3201 radiographs were acquired over 360 sample rotation with an exposure time of 1 s per radiograph, resulting in a total scan length of 1 h 30 min. The transmission images were reconstructed using the commercial software from ZEISS to an isotropic voxel size of 12:5 μm.

Data availability
The data sets generated and analysed during this study are available from the corresponding author upon reasonable request.

Code availability
The code used for data analysis as well as for display of the data is available from the corresponding author upon reasonable request.