Quantitative phase contrast imaging of a shock-wave with a laser-plasma based X-ray source

X-ray phase contrast imaging (XPCI) is more sensitive to density variations than X-ray absorption radiography, which is a crucial advantage when imaging weakly-absorbing, low-Z materials, or steep density gradients in matter under extreme conditions. Here, we describe the application of a polychromatic X-ray laser-plasma source (duration ~0.5 ps, photon energy >1 keV) to the study of a laser-driven shock travelling in plastic material. The XPCI technique allows for a clear identification of the shock front as well as of small-scale features present during the interaction. Quantitative analysis of the compressed object is achieved using a density map reconstructed from the experimental data.

freely in the vacuum after they have passed through the object of interest. The distance between the object and the detector determines the degree of interference between phase-shifted and unperturbed rays.
Our X-ray source was produced by irradiating a solid target at intensities larger than 10 17 W/cm 2 , generating hot electrons that cross the backlighter target. These hot electrons can interact with bound electrons in the atoms, generating vacancies that are immediately filled by electrons from higher shells. The photons that are emitted during this bound-bound transition have wavelengths dependent on the transition energy e.g. L-shell to K-shell transitions result in Kα emission. These sources are characterized by subpicosecond time duration 27 (shorter than the hydrodynamic time) and a spectrum with few energetic lines (Kα, Kβ etc.) sitting on Bremsstrahlung continuum 28 . The energy conversion efficiency, from laser to X-rays, for on-target laser intensities above 10 18 W/ cm 2 oscillates between 29 1 × 10 −5 and 3 × 10 −4 . As an example, using 1 J laser pulse, the X-ray photon flux for the copper Kα transition (8.047 keV) is 8 × 10 10 photons/4π sr/pulse. These characteristics make laser-driven X-ray sources attractive candidates for X-ray imaging in experiments at high-power laser facilities.
In our experiment, we used a 5 m diameter tungsten wire as a backlighter. Using a wire instead of a flat foil limits the spread of electrons in the target and hence limits the source size 29 . Tungsten wire was used in previous experiment at this facility 30 , providing good performance in terms of photon flux, energy and, above all, source size. A small source size makes XPCI possible in such a configuration. Tungsten can be manufactured easily down to 5 m diameter wire. Mass limited target results in a better flux compared to pin-holes which limit the flux. The wire was illuminated by a laser pulse with the following parameters: 0.5 ps time duration, maximum pulse energy of 25 J and wavelength of 1.06 m. The focusing optics used for the backlighter pulse was a 45 deg f = 400 mm off-axis parabola. The FWHM of the focal spot was ~5 m, the energy encircled in this diameter is ~6 J (1/4 of the pulse energy), leading to on-target intensities of around 6 × 10 19 W/cm 2 for the backlighter beam.
The X-ray source was used to image a shock launched inside a polystyrene cylinder (300 m diameter and 300 m long) using a second laser pulse with a wavelength of 1.06 m, energy of 25 J and adjustable time duration between 1 and 10 ns. The focusing optic was a 90 deg f = 1500 mm off-axis parabola leading to a focal spot diameter of 50 m, which corresponds to an intensity on the target of 3 × 10 14 W/cm 2 . We did not use any smoothing techniques (random phase plate etc.). The interaction face of the plastic cylinder was coated with a 100 nm-thick layer of aluminium in order to increase laser absorption and avoid laser shine-through at early times. The time delay between the two pulses was varied between 1 and 10 ns. We chose a plastic object for two reasons. First, the transmission of polystyrene in the X-ray photon range (6-9 keV) is above 70%, allowing for a clear demonstration of the performance of XPCI. Second, we have checked that the EOS used by DUED reproduces with high accuracy the experimental data 31 for plastic at pressures in the range 1-10 Mbars.

Result
A real image is given by the convolution of several contributions, summarized by the following formula 32-34 : where I real is the "real" observed image, I ideal the image obtained assuming an ideal point-like source and perfect detector, PSF det is the detector point spread function (PSF), and σ s is the source size. The last parameter has to be rescaled by the ratio R 1 /R 0 (where R 0 is the source-object distance and R 1 is the object-detector distance). The blurring effect induced by the source is proportional to R 1 , meaning that the smaller the value of R 1 , the lower the negative effects due to the finite source size. To help to choose the correct distances, we developed a simulation code based on the Fresnel-Kirchoff formalism 35 (cf. Image Simulation in Methods), allowing generating a synthetic XPC-image starting from a density map generated by hydrodynamic codes (in particular we used the code DUED 36 ). The parameters used to simulate the shock are the same as those detailed above (laser and object), while the X-ray source was modelled as a Gaussian beam with a FWHM of 5 m and 7.5 keV photon energy. The X-ray energy is close to the average one for a tungsten Bremsstrahlung backlighter.  For the evaluation of Contrast = (I max − I min )/(I max + I min ), the maximum (I max ) and minimum (I min ) light intensity in the neighbourhood of the shock-front were considered. From the plot in Fig. 2, it is clear that increasing the Contrast at magnifications higher than 4 requires the distance R 0 to be increased as well. This behaviour can be explained by the sensitivity of the propagation-based geometry to source lateral coherence L ⊥ = λR 0 /σ s , where λ and σ s are the source wavelength and size respectively. The spatial coherence requirements for propagation-based XPCI have been studied by Wu et al. [38][39][40] . Our experimental configuration is represented in Fig. 2 by a red cross, where R 0 = 25 cm and R 1 = 205 cm. This operating point was chosen as a compromise between the geometrical constraints imposed by the vacuum chamber, the layout of the laser optics and the detector resolution. Our detector was an SR-type imaging plate (IP-SR), with 100 m intrinsic resolution 41,42 . The IP was installed in air, after a 0.5 mm-thick PMMA window and filtered with 13 μm-thick aluminium foil ensuring that all X-ray photons with energies below 5 keV are suppressed.
On each shot, the backlighter spectrum was monitored with a HOPG crystal spectrometer 43 and a gold knife edge was used to measure the source size. The knife edge (see Fig. 1) was placed 55 mm from the source and 152 cm from the detector (IP-SR). The high (26x) magnification used for the ESF ensured that the source-extension contribution dominated over the detector PSF. In Fig. 3, the source Edge Spread Functions (ESF) along the X and Y axes are shown for a single shot. The experimental curve is fitted with a function obtained by convolving a Gaussian with a Heaviside step function. The source widths are different along the two axes, due to the geometry used to illuminate the tungsten wire and the processes involved in laser-matter interaction. As observed by Park et al. 29 , the source size is bigger than the laser focal spot due to lateral spreading of hot electrons  www.nature.com/scientificreports www.nature.com/scientificreports/ in cold matter. In our case, the width along the X-axis was constrained by the wire diameter (5 um), which means that the resolution on this axis is limited by the intrinsic detector resolution, corresponding to 10.6 um on the object plane. On the Y-axis, no physical constraints were present besides the wire length. Hence, due to an average source size of 30 um, the resolution was limited to 27 um. The average photon energy was measured to be around 7.5 keV. The source size so measured is later used to simulate the XPC-image, after being rescaled by the source magnification (R 1 /R 0 ) as reported in Eq. 1.
In order to compare experimental and measured images, both are normalized using a flat-field image without the object. During the experiment, it was not possible to simultaneously acquire an object image and a flat-field image. Moreover, the reproducibility of the X-ray source (in terms of flux, source size and spectrum) was not sufficient to allow using a single flat-field measurement for all data. Therefore starting from the acquired image, we estimated the flat-field performing a polynomial fit to each column within the image, excluding any pixels containing the object.
XPC-images of the object without a laser-driven shock are shown in Fig. 4. Experimental images are on the left-hand side (shot #14881) and the simulated versions are positioned on the right (labelled as Sim.). The synthetic image was generated using a density map of the unperturbed object, the measured spectrum and source size (5 × 30 μm) as source parameters. Profiles taken along the object axis for both images are shown in Fig. 4c. In order to improve the signal-to-noise ratio we performed an average over a 10 pixel width area (shown in red in Fig. 4). The agreement between the two profiles, confirms that the code is able to reproduce both phase-contrast and absorption. The first one is clearly visible as contrast enhancement at object boundary. Pure absorption contribution is visible in the middle of the object.
In Fig. 5, we present some raw phase contrast images that clearly display the matter-vacuum and compressed-uncompressed matter interfaces in laser-driven shock-waves. Note that the pump laser was not perfectly aligned with the object on the last two shots (#14902, #14896). In Fig. 6, we show the same images but filtered with the so-called "wrong filter" (in Method). As it will be explained in the Method section, this filter allows reducing noise while enhancing the edge-contrast. We have added a dashed line to highlight the position of the shock front and used white arrows to indicate the laser interaction point. The shock front is clearly visible in all the images. In addition other, more subtle features can also be discerned. In shot #14902, for example, one can see a plastic shrapnel  www.nature.com/scientificreports www.nature.com/scientificreports/ blown off during the interaction (see Fig. 6). In shot #14896 instead, the laser illuminated the metal support (white arrow), resulting in a shock-wave driven by the emitted radiation, which propagates parallel to the object axis. In Fig. 5, the images #14886 and #14890 show two perfectly-aligned shots. The images were obtained using identical laser pulse parameters: energy, temporal length and delay between the pump and the probe (8.3 ns). In shot #14886, however, the phase contrast is lower because the X-ray source intensity was 30% lower than for shot #14890, which means a lower signal-to-noise ratio. Moreover the source dimensions were larger on shot #14886 (10% and 20% larger along the X-axis and Y-axis respectively). The finite-source size acts as frequency filter, smoothing and removing the oscillations due to phase-contrast. In Fig. 7c, we compare the profiles taken over an area along the object axis (red box in Fig. 5) with the profiles from the simulated images in Fig. 7a,b. Figure 7a is generated taking in account only the absorption contribution, while in Fig. 7b both effects (absorption and phase-contrast) are taken in account. The Sim. XPC + Abs. simulation reproduces the object boundary (arrow 4 in Fig. 7), which is still intact 8.3 ns after the laser interaction. The shock front (arrow 1 in Fig. 7) is also well reproduced: the discrepancy of 13 μm between shot #14886 and the simulation is commensurate with the imaging resolution (10 μm) in the object plane. Instead, the XPCI code fails to reproduce structural features between the shock front and the object boundary (between arrows 2 and 3 in Fig. 7). This behaviour can be explained by the presence, of a strong density gradient behind the shock-front, feature not shown by the hydrodynamic simulations. In ref. 18 we argued that the discrepancy between the arrows 2 and 3 as a consequence of high intensity spikes inside the laser focal spot (no smoothing techniques were used). To mimic this behaviour we performed several hydrodynamic simulations where the central spike inside the laser spot was reduced from 50 m down to 5 m. In this way, we proved that a spike in laser intensity can effect the resulting XPC-image. The double border observable in Fig. 7b, indicated by the black arrows 3 and 4, is caused by the geometry of the imaging system. The source is not aligned with the object vertical face, producing a double-layer effect. This was also seen for shot #14890 in Figs. 5 and 6.  Before one can extract quantitative information from a phase-contrast image, target homogeneity must also be taken into account. The phase-retrieval algorithm requires a homogeneous object, which means that the ratio δ/β should be constant over the entire object for all wavelengths. To prove that this assumption is justified, a simulated shock-wave was used as a test and the ratio δ/β was calculated across the object for all relevant wavelengths. We find that the requirement is fulfilled since emission lines from the source were far from the absorption edges of the elements composing our object 44 .
In Fig. 8, we show the density map generated by the hydrodynamic code (blurred according to the source size measured in the experiment) on the left-hand side, and the two maps retrieved from experimental images (shots #14886 and #14890) on the right. The measured density maps are the result of a phase-retrieval code plus a Feldkamp-Davis-Kress (FDK) tomographic reconstruction algorithm 45 was used. FDK is an algorithm suitable for point-projection tomography. A tomographic algorithm was chosen instead of Abel inversion for the following. First, the idea is to implement the diagnostic also at facilities where it is possible to acquire multiple projections at different angles. Second, one of the main problems of Abel inversion is its instability with noise fluctuations. The FDK algorithm is more stable because it works in the Fourier space, with a consequent smoothing of the noise fluctuation. Of course, an external Fourier filter could be applied before proceeding with an Abel inversion. The problem would be optimizing the filter magnitude in order to smooth noise without compromising the data.
The reconstructed slices along the object axis (#14886 and #14890 in Fig. 7) show the advancement of the shock-wave inside the unperturbed object. In Fig. 9a, we compare the density profile along the Y-axis (in front of the shock-front) with the theoretical prediction. We see that the experimental curve follows the theoretical one, noting that modulations on top of the experimental line are artefacts of the reconstruction code. In Fig. 9, dashed lines represent the density profile convolved with a Gaussian function in order to reproduce the imaging system resolution (the values used are 10.6 μm × 27 μm as estimated above). For the Y-profile, there is a good agreement between the experimental profile and the blurred simulated profile on both shots. Figure 9b instead, shows a significant discrepancy between the simulated density profiles along the object axis and the reconstructed ones. Here, the hydrodynamic code estimates a shock front density of 2.3 g/cm 3 , while the density inferred from the experimental images is 1.3 g/cm 3 . The dashed line is the simulated profile which is treated with the same Gaussian filter as before. The shape of the shock is roughly reproduced by the simulation, however, the peak density of 1.7 g/cm 3 is still higher than the measured value. Furthermore, the behaviour of the density behind the shock is quite different in the simulation and the experiment. The reconstructed density profiles show a faster decay than the simulated one behind the shock. This discrepancy can be explained by the difficulty of accurately modelling the laser-target alignment and the detailed intensity distribution of laser focal spot in the hydrodynamic simulator DUED 36 . It was not possible to acquire an image of the focal spot and pin down exactly its position on the object face during the experiment (especially along the horizontal axis). Every shot required to repeat the laser alignment and focusing procedure for both backlighter and object in order to compensate for the differences between each target.

Discussion
The use of XPCI for shock-wave imaging and characterisation has been reported in several previous works 16,[18][19][20][21] , which demonstrated how XPCI is superior to X-ray absorption imaging in shock detection. However, only recently has XPCI was demonstrated to work also with for broadband laser-plasma X-ray source 18,19 , in a configuration similar to the one employed in ICF facilities. www.nature.com/scientificreports www.nature.com/scientificreports/ In this paper, we have presented an analysis of laser-driven shock-waves using a laser-plasma X-ray backlighter. XPCI allowed obtaining higher object definition than absorption radiography [6][7][8] . We verified that XPCI can reveal structural features that are undetectable using x-ray absorption techniques 3,5 (see Fig. 7). Enhanced contrast at density interfaces, that occurs in XPCI, allows for clear identification of the shock front.
XPCI, as absorption radiography 6 , is also suitable for quantitative analysis of the acquired image. In this work, we have employed two approaches. In the first one we have tried to reproduce the experimental images starting from hydrodynamic simulations. In the second, we have used the acquired images to generate density maps and comparing them to the hydrodynamic simulations. Both ways have shown discrepancies between simulated and experimental data. As explained above, these discrepancies can be attributed to the difficulty of identifying the complex experimental parameters (e.g. laser alignment and laser spot morphology) that are used as inputs to the hydrodynamic simulations. where r is the length of the vector from the source S to the point P 1 on the detector plane, t(P) is the object transfer-matrix, k is the wave number, R 0 and R 1 are respectively the distance source-object and object-detector, r 0 and r 1 are the lengths of the vectors SP and PP 1 respectively (see Fig. 1 for details). Equation 2 is sufficient to generate a phase-contrast image. Nevertheless, two considerations are necessary. First, the computation of a triple integral is required (t(P) is the result of an integral). To implement Eq. 2 in a code at least six loops (for a 2D image) have to be used, resulting in a large computational time. Second, here the source has been considered point-like. In a real experiment, the source has a finite size producing a blurring effect on the detected image. However, if the waves propagate a long distance after the object and the transverse dimension of the sample itself is small compared to the propagation distance, the Fresnel approximation (or paraxial R 0 /r 0 ≃ R 1 /r 1 ≃ 1) can be used. This requires the inequality N eff > 1 to be fulfilled where N F,eff is the so-called effective Fresnel 4 number equal to: where R = R 1 + R 0 , λ is the weighted sum of the source spectrum wavelength (Eq. 10), M = (R 0 + R 1 )/R 0 is the magnification, σ obj is the width of the smallest object feature, σ src is the source width and PSF det is the detector PSF.
Considering the parameters of the employed imaging system, an average photon energy around 7.5 keV, the average source and detector PSF ranging respectively between 5 ÷ 50 μm and 100 ÷ 200 μm, N eff > 1 is fulfilled. Equation 2 can then be reduced to a convolution integral and solved using Fourier transform. Applying the Fresnel approximation, Eq. 2 in the Fourier space becomes: www.nature.com/scientificreports www.nature.com/scientificreports/ where u and v are the spatial frequencies corresponding to the variables x 1 and y 1 , and T(Mu, Mv) is the transform of the object transfer function t(P). This last quantity is a function of (Mu, Mv) because the transform operator is applied with respect to x 1 . In this form, the equation is suitable for numerical evaluation using the Fast Fourier Transform (FFT) algorithm. The sample is described in the above formalism with the transfer function t(P): where the two quantities in the argument of the exponential are obtained from the following integrals: The two integrals are calculated along the path from the source point S to the detector point P 1 . Equations 5 and 6 are only valid in the thin sample approximation. A sample can be considered thin as long as the smallest features to be imagined are larger than the quantity λd , where d is the sample thickness 40 . At this level, the above equations are coupled with the hydrodynamic code. In particular, the two-dimensional density map generated by a hydrodynamic code in cylindrical geometry is used as an input to Eq. 6. In order to increase the calculation speed, the code uses a dedicated FORTRAN library, written for such a purpose, that performs all the calculations in parallel mode via OpenMP library. The transfer matrix is later used to generate the synthetic image considering all the parameters of the real imaging system, namely polychromaticity, source and detector PSF. Regarding the polychromaticity, for a polychromatic source the final image is a weighted sum of the monochromatic ones 32,40 : poly mono image analysis. The condition required to use the Eq. 4 together with that of object homogeneity, and an X-ray photon energy below 50 keV allow phase retrieving from the measured image, by using the Transport Intensity Equation (TIE) 46,47 . In particular, the object project phase is: where − F 1 and F are the inverse and direct Fourier transform and R eff = R 1 /M is the effective propagation distance. Equation 8 has been implemented in a numerical code via FFT.
As explained before, as far as the inequality N F,eff > 1 is fulfilled, the above equations can also be applied in the case of a polychromatic source 48 . The quantities δ , µ , λ in Eq. 8 are then averaged over the source spectrum and the detector function response D(λ): Prior to proceeding with phase retrieval, the image is processed with a Wiener filter 49 in the Fourier space: 2 where PSF is the total point spread function (PSF src *PSF det ) and K is a constant that takes noise variance in account. The filter has two effects: the first one is to partially remove blurring induced by the imaging system and the second is to normalize the noise. One consideration can be done on Eq. 8. The term 4π 2 on the denominator is required by the FFT used for the numerical evaluation of I rel . Removing such term from the filter will result in a non-physical behaviour at the object border, and in general at all the interfaces. The result is an enhancement of the contrast at the border. We found useful applying such "wrong filter" on the measured images for qualitative analysis.
The calculated projected phase map can be used to evaluate a three-dimensional map. Since only one projection is available, this requires that the object is cylindrical. The reconstruction is performed with a Feldkamp-Davis-Kress (FDK) algorithm 45,50,51 . It allows performing a tomographic reconstruction starting from a set of cone-beam projections of the object. In our case, the phase-map projection ϕ(P) at the object plane, obtained from Eq. 8, is used to generate a three-dimensional spectral averaged 50 where r e is the classical electron radius, N A is the Avogadro number, ρ is the mass density of the compound, Z j and A j are the atomic number and the atomic weight of jth element of the compound, ′ f j is the real part of the dispersion correction factor.
The parameters in Eq. 14 are evaluated using the library Xray-lib described in ref. 53 .