Real-Time Tomography of Gas-Jets with a Wollaston Interferometer

Abstract

A tomographic gas-density diagnostic using a Single-Beam Wollaston Interferometer able to characterize non-symmetric density distributions in gas jets is presented. A real-time tomographic algorithm is able to reconstruct three-dimensional density distributions. A Maximum Likelihood-Expectation Maximization algorithm, an iterative method with good convergence properties compared to simple back projection, is used. With the use of graphical processing units, real-time computation and high resolution are achieved. Two different gas jets are characterized: a kHz, piezo-driven jet for lower densities and a solenoid valve-based jet producing higher densities. While the first jet is used for free electron laser photon beam characterization, the second jet is used in laser wake field acceleration experiments. In this latter application, well-tailored and non-symmetric density distributions produced by a supersonic shock front generated by a razor blade inserted laterally to the gas flow, which breaks cylindrical symmetry, need to be characterized.

1. Introduction

The aim of future advanced accelerator concepts is to use accelerating structures that can sustain higher electric field strengths to downsize the structure compared to conventional microwave cavities. Laser Wake Field Acceleration (LWFA) is one such promising technology, which is being investigated by several research groups around the world. The idea of LWFA is to excite a wake field in an underdense plasma by a short energetic laser pulse, the so-called driver. The plasma can be generated by ionizing a gas jet using the driver pulse itself. In this plasma, the driver pulse excites electron density oscillations, which are co-propagating with the driver almost at the speed of light forming a plasma wave, which is moving with a relativistic phase velocity. Due to the local charge separation (the ions remain stationary on the corresponding time scales), longitudinal electric fields with amplitudes in excess of 100 GV/m are generated, orders of magnitude higher than in conventional accelerator structures. Electrons can be accelerated in this moving field structure when they are trapped in the wave with the right phase, i.e., at the correct position and with the correct velocity. An inherently synchronized method for trapping electrons is down-ramp injection, which can be induced when the plasma wave propagates across a sudden drop in density, a so-called down-ramp [1,2]. Here, the plasma wave, which has been generated in a first, high-density region, enters a second region with lower density. The electrons in the high-density region oscillate faster since the plasma frequency ${\omega }_p=\sqrt{n_e{e}^2/\left({\epsilon }_0{m}_e\right)}$ depends on the electron density $n_e$; $m_e$ and e are the mass and charge of the electron, respectively. Therefore, the wavelength of the plasma wave, ${\lambda }_p\approx c/{\omega }_p$, is shorter than in the subsequent low-density region. When this density transition occurs over a distance of the order or shorter than the plasma wavelength, some of the electrons forming a density peak in the plasma wave in the high-density region will then no longer be associated with the density peak of the low-density plasma wave, but they will find themselves in the accelerating phase of the wave. It has been shown that a density change by a factor of two to three can be achieved by a shock front in supersonic gas jets [3]. At such a sharp transition, many electrons are loaded into the same phase-space volume of the plasma wave, which will result in a narrow energy spread of the accelerated electrons. The parameters of the generated electron beam (energy spectrum, charge and duration) critically depend on the plasma wave and its evolution, which in turn is very sensitive to the plasma density distribution. With these facts in mind, a precise and fast real-time density measurement is required for controlling the down-ramp process [4] and (when shot-to-shot-fluctuations in the gas jet cannot be sufficiently suppressed) for increasing the shot-to-shot stability of the generated electron beam. The view of high quality beams from LWFA, a fast and accurate characterization of the gas jet is of great importance. A comprehensive overview on density measurements of gas jets can be found in [5].

While cylindrically-symmetric gas jet distributions can easily be characterized from a single interferometric measurement taken in one direction, e.g., using a Michelson-type interferometer or a Nomarski-type interferometer using a Wollaston prism, tomographic methods are necessary when non-symmetric distributions, e.g., caused by a density jump as mentioned above, are required. The latter methods, which include interferometric measurements of the gas density taken along many directions, are more elaborate and sometimes time-consuming to analyze, when high accuracy and spatial resolution are needed [6,7,8]. Therefore, analysis methods that can be applied in quasi-real time are highly advantageous. Here, we present the physics setup and the computational method, which fulfill the above demands.

The paper is organized as follows: Section 2 introduces the principle of the density measurement using a Wollaston prism. The measurement is then put in the context of LWFA application. In Section 3, the details of the optical setup, as well as the setup for tomography of a non-symmetric density distribution are introduced. Section 4 presents experimental results from the characterization of different gas jets, while Section 5 summarizes the paper.

2. The Single Beam Wollaston Interferometer Setup

2.1. Theory

A Wollaston prism (Figure 1B) is a combination of two prisms made of a birefringent crystal (e.g., quartz), which has a single optical axis. The two prisms are cut and combined such that the optical axes of the prisms are perpendicular to each other. The main feature of such a birefringent crystal is that the refractive index for light polarized parallel to the optical axis (${\eta }_e$) is different from that for light polarized perpendicular to the optical axis (${\eta }_o$). Hence, a light ray polarized at ${45}^{\circ }$ with respect to the optical axis is split up into two equally-intense beams with perpendicular polarization, which experience different refractive indices and hence propagate with different phase velocities. In a Wollaston prism, one of the two beams experiences a higher and the second a lower refractive index in the first birefringent prism, while the situation is reversed in the second prism. Hence, the two rays, which still propagate collinearly in the first prism, will be refracted differently at the interface between the two prisms and then propagate under an angle in the second prism. After exiting the second prism at its back surface, the two rays will diverge with an angle $ϵ$ (dispersion angle), which depends on the geometric and optical properties of the crystals forming the Wollaston prism. As shown by Small [9],

$$ϵ=2\beta ({\eta }_e-{\eta }_o),$$

where $\beta$ is the prism angle of both prisms (see Figure 1A). The Wollaston prism used in this experiment was manufactured by Societe d’Optique de Précision Fichou (Figure 1B) and has a dispersion angle of ${20}^{″}$ = 5.8 mrad. If two light rays with initial relative angle $ϵ$ pass through the prism, one polarization component of the first ray will coincide with the polarization component of the second ray (compare Figure 1A). Before passing through the lens, these two rays are separated by a distance $d=ϵf_{\circ }$. For the Wollaston prism and a lens with focal length $f_{\circ }=30$cm, $d=1.74$ mm. If a second polarizer at $-{45}^{\circ }$ is placed after the prism, these two rays can interfere. This concept is illustrated in Figure 1A.

Depending on the position of the prism relative to the focal point of the lens, the interfering rays will pick up a phase difference. If the prism is centered at the focal point, every optical path runs equal distances through the two halves of the prism, i.e., no relative phase difference is generated. This scenario is called normal mode or the Infinite Fringe Width (IFW) setup. Displacing the prism by a distance b from the focal point results in different path lengths of interfering rays, as they travel different distances in each half of the prism. These phase shifts lead to regular interference patterns on the screen. This setup ($b\ne 0$) is called differential mode or the Finite Fringe Width (FFW) setup [10]. The FFW setup is used in this work and is further described in Section 2.3.1. The spacing S between the undisturbed fringes is given by:

$$S=\frac{\lambda }{ϵ}\frac{p}{b},$$

(1)

where p is the distance from the Wollaston prism to the screen [11]. The spacing can be decreased by increasing b, i.e., by changing the position of the Wollaston prism with respect to the focal plane. If the position of the screen remains fixed, which is the case when the lens images a certain plane on the screen, p changes accordingly. Placing a medium with refractive index $\eta \ne 1$ that covers only parts of the laser beam will result in a shift of the fringe spacing S, since rays passing through this medium will pick up an additional phase shift with respect to the unperturbed rays.

2.2. Estimation of Phase Shift Due to a Gas Jet Density Distribution

A local density gradient imposes a varying refractive index $\eta$, resulting in a phase shift $\Delta \varphi$ of rays passing through that particular region in comparison to rays that would have propagated through a vacuum. Therefore, the undisturbed fringe spacing is locally shifted by $\Delta S$. The relationship between the atomic density $\rho$ and refractive index $\eta$ of a gas is given by the Lorentz–Lorenz equation [12]:

$$\rho =\frac{{\eta }^2-1}{{\eta }^2+2}\frac{N_A}{A}\approx (\eta -1)\frac{2}{3}\frac{N_A}{A},$$

(2)

where $N_A$ is Avogadro’s number and A is the molar refractivity, e.g., $A_{Ar}\approx 4.20\times {10}^{-6}{\mathrm{m}}^3/\mathrm{mol}$. The approximation is valid for $\eta \approx 1$, which is fulfilled for a gas with densities $n<{10}^{19}/{\mathrm{cm}}^3$. Inverting Equation (2) yields the dependency of the refractive index on the density. The phase shift between interfering rays is given by $\Delta \varphi =\varphi \left(y\right)-\varphi (y-d)$, where $\varphi$ is the phase imposed by the neutral gas atoms and y is the coordinate along the direction perpendicular to the interference fringes. Assuming the linear approximation (Equation (2)), a homogeneous density within a gas jet of diameter l yields:

$$\varphi \approx l(\eta -1)\frac{2\pi }{\lambda }=l\frac{3\pi }{\lambda }\frac{A}{N_A}\rho .$$

(3)

This expression yields an upper limit for the phase shift, since the gas jet has a circular throat and, hence, produces a cylindrically-symmetric gas flow, i.e., the path length of a ray propagating through the gas is at most l. As described above, the fringe spacing is affected by the phase shift in the FFW setup; the fractional fringe distance shift is given by:

$$\frac{\Delta S}{S}=\frac{\Delta \varphi }{2\pi }.$$

(4)

For comparison purposes, and because higher pressures than 15 .... 20 bar are possible, the gas density using the solenoid valve (Parker Miniature High-Speed Valve) will be estimated. The pressure of the gas before leaving the nozzle, the backing pressure $P_b$, as well as the nozzle design play crucial roles with respect to the gas density distribution in the jet. The gas jet used in the present experiment allows for backing pressures up to 80 bar and has a conical nozzle with a half opening angle $\theta ={45}^{\circ }$ (see Figure 2). As described by Chen [13], the on-axis particle density at a height x of a conical gas jet is approximately given by:

$$\frac{\rho }{{\rho }_0}\approx 0.15{\left(\frac{0.74d}{xtan\theta }\right)}^2,$$

where x is the distance to the throat of the nozzle (see Figure 2), ${\rho }_0$ is the atomic number density of the gas before leaving the nozzle and d is the throat diameter of the valve. The gas jet used in this study has a throat diameter $d=500\mathsf{\mu }$m and $L=250\mathsf{\mu }$m. Using the ideal gas law, which is a good approximation for the monoatomic gas argon at room temperature, ${\rho }_0$ is determined by the backing pressure $P_b$, which refers to the pressure in the chamber before the nozzle, as well as the temperature of the nozzle $T_0$ via: ${\rho }_0=P_b/\left(k_B{T}_0\right)$, where $k_B$ is the Boltzmann constant. Figure 3 shows the on-axis density $\rho$ with respect to h (distance from the end of the nozzle) for different backing pressures from 10–40 bar for the dimensions of the Parker solenoid valve.

A practical example for gas jets to produce a plasma for LWFA is as follows: in order to excite the plasma wake resonantly with a laser pulse of a 15-fs duration, an electron density in the order of $1.75\times {10}^{18}{\mathrm{cm}}^{-3}$ is necessary [14]. The required density (assuming argon being ionized eight-fold) at $h=2.5$mm (distance from the nozzle) is achieved with backing pressures of 30 bar. According to Equation (3), the phase shift of a gas jet with an average density ranging from $1\times {10}^{18}$–$1\times {10}^{19}{\mathrm{cm}}^{-3}$ lies between 0.1 rad and 1 rad. A phase shift of this order is expected to result in a clearly visible change of the interference pattern described by Equation (4).

2.3. Experimental Setup for Gas Density Measurements

The Wollaston interferometer to measure the gas density is set up on a single breadboard (75 cm × 125 cm), such that the optical setup and vacuum chamber with the gas jet form a compact and transportable unit.

2.3.1. Wollaston Interferometer

The interferometer is depicted schematically in Figure 4 and Figure 5. A linearly-polarized continuous He:Ne laser with a wavelength ${\lambda }_I=632.8$ nm and an output power of 21 mW is used as a coherent light source for the interferometer. To prevent up-streaming air due to the heat of the laser from causing unwanted phase shifts, the laser is placed outside of the black cardboard box that contains the interferometer. The box reduces noise as it minimizes the airflow in the whole experiment, as well as external photons hitting the charge-coupled device (CCD) sensor. The initial laser diameter is 0.7 mm (1/$e^2$ width). The telescope (20×) attached to the laser provides a beam with a diameter of 14 mm, which is suitable to backlight a gas jet with a diameter of a few millimeters. The noise of the acquired images is expected to be lowest in the center of the beam due to the higher intensity of the laser in that region. The He:Ne laser enters and leaves the vacuum chamber through two windows, and the gas flows in perpendicular direction, downwards towards the vacuum pump. The Parker solenoid valve (Figure 6A) is operated at 0.5–3 Hz and with opening times $T\le 12$ ms to reduce the gas load on the pump. The Wollaston prism is placed between two crossed polarizers, such that the interference fringes are parallel to the jet. The interference pattern is captured with a CCD camera and an f = 200 mm camera lens. The minimal exposure time of the camera is 18 $\mathsf{\mu }$s. The exposure time is set to 30 $\mathsf{\mu }$s in order to obtain an image with good contrast without reaching saturation of the sensor.

Figure 6B shows a typical interference pattern of the undisturbed gas jet. The change of the fringes due to the gas distribution is clearly observable when the Parker solenoid valve is operated at 35 bar backing pressure.

2.4. Experimental Set-Up for Non-Rotational Measurements

Figure 7 depicts the experimental setup to obtain tomographic projection data of the solenoid gas jet with a supersonic shock front generated by inserting a razor blade into the gas flow; in order to keep the gas jet and the imaging system fixed, but to position the blade around the gas jet with two degrees of freedom. This is achieved with a rotational piezo stage and a linear piezo positioner, which is attached to the rotational stage. The rotational stage allows one to measure phase projections along different directions with respect to the orientation of the shock front. The resulting density gradients are shown in Section 4.2. The radial positioning can then be used to vary the position of the blade with respect to the center of the gas jet. This allows us to change the properties of the shock front. Furthermore, the radial degree of freedom can be used to correct the eccentricity between the center of the gas jet and the axis of the rotational stage due to mechanical imperfections. The vertical distance between the razor blade and the nozzle of the gas jet is fixed in this setup to 1.6 mm.

3. Data Analysis for Interferometry and Real-Time Tomography

The Wollaston interferometry setup introduced in the previous section can be used to measure phase projections along the propagation direction of the interferometry laser. An optically-transmitting object that has a refractive index close to one (e.g., a gas jet) can be characterized. If the studied object can be assumed to be rotationally symmetric around a central axis, an Abel inversion yields the 3D density distribution from a single phase projection, i.e., a 2D projection obtained in one single observation direction, which is perpendicular to the axis of symmetry. In order to obtain the 3D density distribution of an object without such symmetry (e.g., a gas jet with a shock front), several projections along different angles are needed. The 3D distribution can then be reconstructed with tomographic algorithms, for instance the Maximum Likelihood Expectation Maximization (ML-EM), cf. Section 3.2.1. The images showing the interference fringes obtained with the Wollaston interferometer (Section 2.3.1) need to be evaluated numerically in order to obtain quantitative information about the phase projections and hence to be able to reconstruct the asymmetric, spatial density distribution.

3.1. Numerical Tools

The undisturbed fringes are mathematically described by a harmonic oscillation with a fixed frequency multiplied by the Gaussian amplitude of the laser beam profile. The wave period of the oscillation is given by the undisturbed fringe width S (Equation (1)), i.e., the spacing between the fringes. A gas jet will locally change the fringe width according to the phase introduced to the different light rays in the laser beam by the gas flow. The following phase extraction problem arises [15]: given noisy discrete values of a function of the form:

$$g\left(y\right)=b\left(y\right)·e^{\mathrm{i}[2\pi fy+\alpha (y\left)\right]},$$

where $\alpha \left(y\right)$ represents the phase shift induced by the gas jet. The term $b\left(y\right)$ takes global intensity changes due to the Gaussian profile of the interferometry laser into account; y is the coordinate perpendicular to the gas flow and interferometry laser; and f is the oscillation frequency of the undisturbed fringes. This problem can be solved by fitting the measured data to an ansatz for $\alpha \left(y\right)$. However, this is not appropriate, since one has to make (possibly incorrect) assumptions about the form of $\alpha \left(y\right)$. This is particularly true when a shock front has to be characterized by $\alpha \left(y\right)$. A more elegant way to directly unwrap the phase from the noisy data makes use of a Fourier transformation:

$$\mathcal{F}\left(g\right)\left(\omega \right)={\int }_{-\infty }^{+\infty }{e}^{-\mathrm{i}2\pi \omega y}·b\left(y\right)·e^{\mathrm{i}[2\pi fy+\alpha (y\left)\right]}\mathrm{d}y.$$

Since the acquired signal is real, the full information is contained in the positive frequency domain where the spectrum has two significant peaks: one at the fringe frequency f and another around zero due to the slowly varying Gaussian intensity profile of the beam. A Gaussian window is applied to the spectrum to cut out the low frequency peak. In order to eliminate the $2\pi fy$ phase factor, we apply a rotation of the Fourier-transformed signal by:

$$R_f\mathcal{F}\left(g\right)\left(\omega \right)=\mathcal{F}\left(g\right)(\omega +f)={\int }_{-\infty }^{+\infty }{e}^{-\mathrm{i}2\pi \omega y}·b^{\prime }\left(y\right)·e^{\mathrm{i}\alpha \left(y\right)}\mathrm{d}y,$$

with the result that the rotation by $R_f$ shifts the peak at frequency f to zero.

In order to obtain the sought phase shift $\alpha$, we apply the inverse Fourier-transform:

$${\mathcal{F}}^{-1}\left(R_f\mathcal{F}\left(g\right)\right)\left(y\right)=b^{\prime }\left(y\right)·e^{\mathrm{i}\alpha \left(y\right)}=:A\left(y\right).$$

As a result, A contains the information about the phase shift $\alpha$, as well as remaining effects of b denoted by $b^{\prime }$. These are minor effects due to imperfections of the mirrors, lenses and polarizers or inhomogeneities of the two windows through which the laser passes. These effects are eliminated by the point-wise product $A·A_{ref}^{-1}$, where $A_{\mathrm{ref}}$ is the signal obtained by the same transformation applied to the reference data, i.e., without gas jet. The phase of the remaining term corresponds to the phase caused by the gas jet alone.

The term $e^{\mathrm{i}2\pi fy}$ would arise in $A_{\mathrm{ref}}$, as well, and is therefore canceled in the last step, even if the rotation $R_f$ is omitted. As explained in Section 2, the Wollaston interferometer measures the phase difference between rays separated by a distance d in the object plane. Therefore, phase shifts obtained at two points in the image plane that are separated by a distance corresponding to d have to be added up. We obtain the additional phase $\varphi$ by $\varphi \left(y\right)=\frac{2\pi }{\lambda }\int (\eta (y,z)-1)\mathrm{d}z.$ This is the phase that a ray has picked up when passing through the gas jet, in comparison to a ray that has passed through a vacuum only. The integration is taken along the ray’s path at position y; $\lambda$ is the wavelength of the laser; and $\eta (y,z)$ refers to the distribution of the refractive index in the gas jet. As long as cylindrical symmetry of the gas distribution is assumed, obtaining the radial density distribution $\rho \left(r\right)$ from $\varphi \left(y\right)$ is possible. The corresponding integral relation is called inverse Abel transformation or Abel inversion. The phase shift $\varphi$, extracted with the method explained before, is proportional to a projection of the gas jet along the propagation direction of the laser. Here, cylindrical symmetry of the gas distribution is assumed. This assumption enables one to reconstruct the radial density distribution from $\varphi$. This is achieved by the inverse Abel transform [16,17], which reads:

$$\varphi \left(y\right)=\frac{\Delta s}{\lambda }2\pi =\frac{2\pi }{\lambda }2{\int }_0^{\infty }\left[\eta (y,z)-1\right]\mathrm{d}z=\frac{4\pi }{\lambda }{\int }_y^{\infty }\left[\eta \left(r\right)-1\right]\frac{r}{\sqrt{r^2-y^2}}\mathrm{d}r,$$

$$F\left(y\right)=\varphi \left(y\right)·\frac{\lambda }{2\pi }=2{\int }_y^{\infty }\left[\eta \left(r\right)-1\right]\frac{r}{\sqrt{r^2-y^2}}\mathrm{d}r,$$

$$\stackrel{\mathrm{inv}.\mathrm{Abel}}{\Rightarrow }\eta \left(r\right)-1=-\frac{1}{\pi }{\int }_r^{\infty }\frac{\mathrm{d}F\left(y\right)}{\mathrm{d}y}\frac{\mathrm{d}y}{\sqrt{y^2-r^2}}.$$

(5)

The derivative and the integral occurring in Equation (5) are calculated following [18]. The key points are summarized here. It is not suitable to approximate the derivative of F by the discrete differential quotient due to noise, since this would amplify noise drastically. A more sophisticated method based on Gaussian filters is needed. Consider the following Fourier identity for f, g functions with compact support.

$$\begin{array}{cc}\hfill F\left(f\right)·\mathcal{F}\left(\frac{\mathrm{d}g}{dx}\right)& ={\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }f\left(x\right)·{\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }\frac{\mathrm{d}g\left(x\right)}{dx}\hfill \\ & ={\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }f\left(x\right)·{\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }\left(\mathrm{i}2\pi \omega \right)g\left(x\right)\hfill \\ & ={\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }\left(\mathrm{i}2\pi \omega \right)f\left(x\right)·{\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }g\left(x\right)\hfill \\ & ={\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }\frac{\mathrm{d}f\left(x\right)}{dx}·{\int }_{-\infty }^{+\infty }\mathrm{d}x{e}^{-\mathrm{i}2\pi x\omega }g\left(x\right)\hfill \\ & & =F\left(\frac{\mathrm{d}f}{dx}\right)·\mathcal{F}\left(g\right)\hfill \end{array}$$

(6)

When g is set to a Gaussian distribution, rearranging this identity yields the smoothed derivative of a noisy signal f by using only the derivative of g instead of f (Equation (6)). The derivative of the Gaussian function g can be easily evaluated from its analytical expression. Another issue is the singularity inside the integral of Equation (5), when the integral is approximated on a discrete domain. This is solved by setting the first value of the integral ($y=r$) to the second value ($y=r+\Delta y$). For the limit of many points (i.e., infinitesimal grid spacing $\Delta y$), the numerical value of the integral will converge to the analytical value [18].

3.2. Tomographic Reconstruction

The goal of tomographic reconstruction algorithms is to estimate a 3D density distribution based on its measured projections along directions with different angles. If rotational symmetry can be assumed, an Abel inversion provides an analytic method to reconstruct the 3D density distribution from a single projection. Problems where rotational symmetry cannot be assumed demand for another reconstruction method. These kinds of problems arise frequently in medical physics when a density image of tissue is desired, but only projections from X-ray scans or PET (Positron Emission Tomography) data are available. An efficient algorithm is the so-called back projection. As computation power increased drastically, more sophisticated algorithms were developed. One of these is the Maximum Likelihood-Expectation Maximization (ML-EM) algorithm [19], which is an iterative method with good convergence properties compared to the back projection, method but, on the other hand, having a higher computational cost.

The idea of a simple back projection is to redistribute (back-project) the measured projections homogeneously along the projection lines. This procedure can be implemented efficiently, but it can only provide rough information about the distribution. The reconstructed images become blurred, albeit infinite projections are available [19].

3.2.1. Maximum Likelihood-Expectation Maximization

The basic principle of this algorithm is to advance an initial guess of the distribution iteratively by comparing the forward projected data of the current guess with the measured data from all angles. The n-th estimate of the i-th voxel is named $x_i^n$. Often, a homogeneous distribution is chosen as the initial guess $x_i^0$. The ML-EM step that advances the guess is computed by the following equation [19]:

$$x_i^{n+1}=x_i^n\frac{1}{{\sum }_j{A}_{ij}}\sum _j{A}_{ij}{R}_j^n,$$

where $R_j^n$ is the ratio of the value of measurement pixel $y_j$ and the forward projected data of the n-th estimate:

$$R_j^n=\frac{y_j}{{\sum }_k{A}_{kj}{x}_k^n}.$$

The matrix A is the system matrix, which accounts for how much each voxel contributes to each measurement. By this procedure, the advanced guess will approach a distribution that is compatible with all measurements. More precisely, the computed guess converges to a distribution that has a high probability (maximum likelihood) to be the original distribution, given the measured data. First, the ML-EM algorithm is applied to a single phase projection from the undisturbed gas jet. The reconstructed radial density distribution agrees well with the result from Abel inversion (Figure 8). The Abel inversion shows non-smooth behavior around $r=0$ due to the singularity. This problem does not arise with the ML-EM algorithm.

To validate the performance of the algorithm on non-rotational symmetric distributions, tests with known distributions are carried out. A 2D Gaussian multiplied by a step function is chosen, as the shock wave by the razor blade is expected to have a similar shape. Figure 9 shows the generic distribution, as well as the reconstructed image from $N_a=7$ projections after 15 iterations. Random noise distributed according to a Gaussian distribution with a standard deviation of $\sigma =0.01$ is added to the projection data before running the ML-EM reconstruction, and qualitative reconstruction of the original distribution without rotational symmetry is achieved.

3.2.2. Convergence and Error Studies of ML-EM

The convergence properties of the ML-EM algorithm are summarized in Figure 10, which shows the $L_1$ norm for the first 15 iterations for different numbers of measurements (number of projection angles $N_a$). Gaussian noise (standard deviation $\sigma =$ 0.01–0.05) is added to the measurement. The reconstructed distribution matches the original distribution best if noise is lowest, and more interestingly, the error is not strongly correlated to the number of projections. The downside of many projection angles is that more noise is picked up. It turns out that $N_a$ around seven achieves the best convergence properties for the distribution and a noise level of $\sigma =0.01$.

3.3. Real-Time Computation of the Density Reconstruction

Real-time tomographic reconstruction requires large computational resources. Since the ML-EM algorithm can be parallelized, real-time reconstruction can be achieved using either multi-core Central Processing Units (CPUs) or Graphical Processing Units (GPUs).

The measurement can be divided into three different parts as sketched in Figure 11. In the first step, the actual image is taken. The blade is held at a fixed position, and the gas jet opens. In this experiment, the gas jet’s open-close frequency is limited by the vacuum pump to about 2 Hz. Since we average over 10 shots of the gas jet, to improve the statistics, the time needed for each step in angle is fixed to approximately 5 s. Between two projection measurements at different angles, the blade has to be moved (second step). The maximal angular velocity is 1 deg/s; hence, for a complete measurement (90 deg), the total driving time is about 90 s.

The last step of the measurement is the tomographic reconstruction of the density distribution. Since driving the blade does not require any computational resources, this process can be done in parallel. As soon as the stage is moving, the computation of the density distribution can be started. After approximately 10 s, the first estimate of the 2D density distribution is available and improves consecutively with the number projections from different angles. The time required to reconstruct the density for $N_a=7$ projections is below 2 s.

In order to create the parallel implementation of the algorithm, the TomoPy package for Python [20] together with the Astra toolbox [21] was used. TomoPy is a Python-based framework for tomographic image reconstruction and data processing tasks developed at the Advanced Photon Source of Argonne National Laboratory. TomoPy includes many functions to perform pre-processing and image reconstruction using different algorithms [20]. The Astra toolbox is an open-source project developed at the University of Antwerp. The toolbox provides tomographic image reconstruction of 2D and 3D datasets. The toolbox uses Compute Unified Device Architecture (CUDA) to offload the reconstruction algorithms to the GPUs, but most of the 2D algorithms can also be executed on the CPU [21].

Using the TomoPy toolbox, the algorithm described in Section 3.2.1 is implemented in Python to perform the image reconstruction. Additionally, the Astra toolbox is used in order to target the GPU platform. The GPU platform is mostly targeted in order to evaluate the potential to perform real-time 3D image reconstruction, as for the 2D reconstruction, the power of the CPU is mostly sufficient.

Reconstruction times of the algorithm are shown in Figure 12. The figure shows reconstruction times using 1, 4 and 8 CPU cores, as well as reconstruction offloaded to the GPU. The hardware used in the reconstruction was 2× Intel Xeon E5-2609 v2 CPUs and 1× NVIDIA Tesla K40c. The time represented in the figure shows only the reconstruction time without input and output operations, which are constant for all the implementations. The benchmark was done for a 3D reconstruction using 100 slices.

With these technologies, the bottle neck of the experiment is clearly not the computational part anymore. Even when, for example, a stage like the U-651 from Physik Instrumente (PI) with a maximal speed of 540 deg/s were been used, the time needed to drive the blade could be pushed below one second, still on the same order as the computation.

The achieved speed-ups on multi-core CPUs and GPU will allow performing reconstruction in real time.

4. Examples of Gas Jet Density Measurements

In this section, two different argon gas jets (piezo-driven or solenoid-based valves) are studied under various conditions. A shock front arising from a razor blade inserted into the gas flow of the solenoid jet is evaluated via tomography.

4.1. Piezo Gas Jet for Free Electron Laser Beam Instrumentation

Gas-based monitors are attractive for beam instrumentation, as they can be used to build quasi-noninvasive diagnostics in free electron lasers. The low density of the active medium allows the electrons and photons, respectively, to pass the monitor with only minimal impact on the beam. The gas atoms are ionized by the beam, and the photoelectrons and/or ions can be characterized by appropriate detectors. Piezo gas jets can deliver these atoms directly into the beam line vacuum. When combined with a turbo-molecular vacuum pump and used with microsecond opening times [22], these systems can be operated without the need for beam windows.

SwissFEL is equipped with a Terahertz streak camera for measurements of the X-ray pulse length and arrival time [23]. In this device, the X-ray pulses ionize xenon atoms, and the photoelectrons are streaked by an externally-applied terahertz field [24]. For the X-ray pulse length measurement of SwissFEL, the knowledge of the gas jet diameter determines the maximum interaction length, which is used to assess the Guoy phase shift of the terahertz wave [25].

The interaction point of the X-ray pulse and the exit nozzle of the gas jet is given by the size of the ion or electron spectrometers, respectively. We have thus set up a measurement of a piezo gas jet (Amsterdam Piezovalve) with an opening time of 7 $\mathsf{\mu }$s. We present here measurements performed 30 mm away from the nozzle. In order to measure the gas density up to the relevant distance, the piezo valve is mounted on a linear vacuum stage. The opening time of the valve is set to 55 $\mathsf{\mu }$s to measure the static gas flow without opening and closing effects of the valve (exposure time of CCD: 30 $\mathsf{\mu }$s). Figure 13 shows the radial density distribution, at distances of 5 mm–29 mm from the piezo valve, at 6 bar backing pressure. In the same figure, the full width at half maximum of the radial distribution is shown. The gas expands linearly with a half opening angle of 9${}^{\circ }$ in the considered region. Figure 14 depicts the central density of the gas flow. After expanding over a distance of almost 30 mm, the density is reduced by a factor of 20 down to $1.5\times {10}^{15}{\mathrm{cm}}^{-3}$.

Similar methods may be used to study the electron beam that is used to generate the X-rays in a free electron laser. The tunnel ionization rate of the gas atoms depends on the peak current in the beam, which is a key parameter for the amplification in the FEL. The knowledge of the gas density, and its distribution, are important for the understanding the data.

4.2. Solenoid Gas Jet for LWFA

A disadvantage of piezo valves is that the valve closing mechanism can only sustain backing pressures $\le 10$ bar. Therefore, such valves are not suitable for high density applications. For the purpose of LWFA characterization, high-Z neutral atom densities are required, which are larger by one order of magnitude than the peak density, which the piezo valve can produce. The Parker solenoid valve (Figure 6A) is expected to provide densities exceeding $1\times {10}^{18}{\mathrm{cm}}^{-3}$ at a distance of 2.6 mm from the nozzle, when operated with backing pressures around 35 bar (Figure 3).

Due to geometric limitations in our setup (blade holder), the focused laser beam in the LWFA setup cannot interact with the gas closer to the nozzle than 2.6 mm. The plasma frequency ${\omega }_p$, which is determined by the plasma density, has to be matched to the laser pulse duration in order to excite the plasma wake efficiently. In the experiment, the plasma frequency can be controlled by the density and ultimately by the backing pressure applied to the valve. The central density measurement of the solenoid gas jet (without blade) for argon is summarized in Figure 15. From this plot, one can infer the required backing pressure to create the resonant electron density at a certain distance from the nozzle. For a laser pulse of $\tau =15$ fs duration (FWHM), the required argon density is $1.75\times {10}^{18}{\mathrm{cm}}^{-3}$, assuming that the gas is ionized eight-fold (all valence electrons) [14]. Figure 15 also contains the density estimates from Section 2.2. Good qualitative agreement between the model and the measurement is observed. However, the measured densities are systematically higher than the estimate. A reason may be the fact that the estimate is based on the straight streamline model, which assumes that the expansion angle of the gas flow is equal to the opening angle of the nozzle [13]. However, the measured expansion angle, defined by a linear fit to the FWHM of the gas distribution, is around 52${}^{\circ }$, which is significantly smaller than the opening angle of the conical nozzle (90${}^{\circ }$). This would result in a higher particle density and could explain the discrepancy between the used model and the measurement results. Figure 15 (right) depicts the argon density with respect to the backing pressure applied to the valve at fixed vertical distances h. It is observed that the argon density is increasing for backing pressures up to 37.5 bar. For even higher pressures, the data indicate a stagnation and even a regression of the density. This can be explained by the working principle of the solenoid valve. When operating at higher backing pressures, a larger force is needed to lift the puppet out of the seal. Therefore, the valve may not open properly with too high backing pressures. This can then result in a smaller density. To overcome the force due to the backing pressure, a high-voltage pulse (burst) is applied to the valve, the duration of which can be set internally. For the measurements shown here, the burst duration was set to 220 $\mathsf{\mu }$s.

4.3. Shock Front Characterization in an LWFA

For the purpose of density down-ramp injection in the linear regime of LWFA, a razor blade is inserted laterally into the gas jet to create a supersonic shock front [4]. A typical phase projection of the shock front is given in Figure 16. Two images next to each other are shown, one with positive and one with negative signal values. This is an artifact due to the working principle of the Wollaston prism, since it produces two parts in the interferograms; one represents the reference (unperturbed) and the other one the signal (cf. Section 2.1). From the projection, it can already be noted that the density gradient is decreasing with increasing vertical distance above the nozzle. The tomographic data obtained with the rotational setup, cf. Figure 7, are analyzed via the ML-EM reconstruction algorithm explained in Section 3.2.1. Figure 17 depicts the reconstructed density distribution of the shock front at 35 bar backing pressure. To study and optimize the ramp properties, this measurement is carried out for different backing pressures $P_b$ and blade positions $L_b$ from the center up to the edge of the gas flow (edge: $L_b=-3.2$ mm). The ML-EM tomography is computed at several distances from the nozzle h. The range of these parameters is summarized in

Table 1. Parameters for shock front measurements.

	${\mathit{P}}_{\mathit{b}}$ (bar)	${\mathit{L}}_{\mathit{b}}$ (mm)	h (mm)	h (pixel)
min, max	30.0, 40.0	$-3.2,-2.5$	1.8, 3.0	1200, 1050
step	2.5	0.1	0.2	25

. In order to evaluate the characteristics of the shock front numerically, the following quantities are defined:

• Height $h_s$ of the shock front: density difference between ramped and undisturbed distribution at the ramp,
• Ramp factor r: ramped peak density divided by undisturbed peak density,
• $w_1$: half-width (left) defined by the ramp peak density,
• $w_2$: half-width (right) defined by the height $h_s$.

For a better understanding, the parameters $w_1$ and $w_2$ are indicated in Figure 18, showing the density of a shock front, as well as the undisturbed distribution along the z-direction. The main effect of the backing pressure is the central density, which is plotted in Figure 15. The ramp characteristics are governed by the parameters $L_b$, h and y, and their respective effects are summarized in Figure 19 at a fixed backing pressure of 35 bar. A sharp density down-ramp is desirable to achieve a small electron energy spread in an LWFA ([2,26,27,28]). The length of the down-ramp is quantified by $w_2$, which is increasing with distance from the blade. A slight dependency of $w_1$ on the blade position $L_b$ is observed with a local minimum at $L_b$ = −2.8 mm. Another important quantity is the ramp factor r, which is mainly determined by $L_b$ and y. It is larger when the blade is positioned closer to the center of the gas flow, as the ramp is created in a region of higher undisturbed density. The measurements indicate that the ramp factor also depends on the horizontal position y. In particular, r has a local maximum 0.5–0.6 mm displaced from the center. This is understandable by looking at the profiles of a two-dimensional Gaussian. When the profile is off-centered, the distribution is flatter. This means that for a centered profile, the density is rising more after the ramp, which results in a lower ramp factor.

4.4. Error and Stability Analysis

For the phase measurements of the piezo valve, an average of 1000 images is taken in order to reduce noise. This leads to an error of below 1 mrad (standard deviation $\sigma$), which is needed in order to measure low densities down to $1.5\times {10}^{15}{\mathrm{cm}}^{-3}$ with an uncertainty of $0.4\times {10}^{15}{\mathrm{cm}}^{-3}$.

In the case of the solenoid valve, the purpose is to create and measure a shock front with densities up to three magnitudes larger than in the scenario described above. Combined with the fact that tomography requires more data (phase projections from different angles), it is decided to reduce the number of images per measurement drastically. For instance, the peak density in the shock front at 35 bar, 1 mm away from the nozzle, is 1.5 rad and has a standard deviation of 0.05 rad. After averaging, the phase signal has noise in the order of 0.01 rad. As the maximum of the phase signal is around 1 rad, the ML-EM test (Section 3.2.2) is realistic with a noise level of $\sigma =$ 0.01 (the generic distribution for this test is normalized, such that the maximum of the projections is equal to one). In agreement with the ML-EM convergence test, the number of iterations for the 3D density reconstruction is set to 10. Figure 20 shows the density distribution of the shock front at 35 bar after 10 iterations and the noise of the reconstructed image. The noise is calculated as follows. The data are fitted with a Savitzgy–Golay filter, which smooths noisy data by interpolation within a symmetric window around each data point (window size: 51 pixels, polynomial order: three) [29,30]. The fit is subtracted from the noisy data, then divided into segments of 30 pixels (0.2 mm), over which the standard deviation is calculated. This provides a local error estimate of the reconstructed density. The noise ($\sigma <1.4\times {10}^{17}{\mathrm{cm}}^{-3}$) in the region of the shock front is a factor of 20 lower than the peak density $2.8\times {10}^{18}$, i.e., the density change in the shock front is reconstructed with sufficient precision.

5. Conclusions

A simple single-beam Wollaston interferometer has been described. The interferometer was used to measure gas densities in electron beam monitors for the free electron laser and to characterize shock fronts in an LWFA. Convergence and error studies of the used ML-EM algorithm show adequate accuracy of the presented problems. The use of parallel computation and GPU technology reduces the computational time below the data taking time. In that sense, this represents real-time density computation.

The presented setup is limited by the slow rotational stage and the low frequency of the gas jet due to the weak vacuum pump. These technical limits can be easily overcome; hence, the time for a full 3D density reconstruction is estimated to be less than 10 s.

We plan to explore the possibility of interleaving the measurement with the operation of the LWFA. We aim at periodic measurements of the gas-jet density and potentially integrate the obtained results, in a feedback loop, with the overall goal to improve the quality and stability of the electron beam.

Further studies will also include the investigation of better initial conditions w.r.t. the convergence of the ML-EM algorithm and the understanding of the artefacts such as the small peaks at the center, visible in the 2D density reconstructions.