Fluid dynamics analysis of a gas attenuator for X-ray FELs under high-repetition-rate operation

Newtonian fluid dynamics simulations were performed using the Navier–Stokes–Fourier formulations to elucidate the short time-scale (µs and longer) evolution of the density and temperature distributions in an argon-gas-filled attenuator for an X-ray free-electron laser under high-repetition-rate operation. Both hydro­dynamic motions of the gas molecules and thermal conductions were included in a finite-volume calculation. It was found that the hydro­dynamic wave motions play the primary role in creating a density depression (also known as a filament) by advectively transporting gas particles away from the X-ray laser–gas interaction region, where large pressure and temperature gradients have been built upon the initial energy deposition via X-ray photoelectric absorption and subsequent thermalization. Concurrent outward heat conduction tends to reduce the pressure in the filament core region, generating a counter gas flow to backfill the filament, but on an initially slower time scale. If the inter-pulse separation is sufficiently short so the filament cannot recover, the depth of the filament progressively increases as the trailing pulses remove additional gas particles. Since the rate of hydro­dynamic removal decreases while the rate of heat conduction back flow increases as time elapses, the two competing mechanisms ultimately reach a dynamic balance, establishing a repeating pattern for each pulse cycle. By performing simulations at higher repetition rates but lower per pulse energies while maintaining a constant time-averaged power, the amplitude of the hydro­dynamic motion per pulse becomes smaller, and the evolution of the temperature and density distributions approach asymptotically towards, as expected, those calculated for a continuous-wave input of the equivalent power.

1. Introduction

Currently there are many short-wavelength free-electron lasers (FELs) that have been successfully built and are in user operation, including the Free electron LASer at Hamburg (FLASH) VUV FEL (Ackermann et al., 2007), the Linac Coherent Light Source (LCLS) (Emma et al., 2010), the SPring-8 Ångstrom Compact Free Electron Laser (SACLA) (Ishikawa et al., 2012) and the FERMI@Elettra free-electron laser (Allaria et al., 2012). Many more facilities are either near completion or being commissioned, including the European XFEL in Hamburg, Germany, the SwissFEL X-ray free-electron laser, in Villigen, Switzerland, and the PAL-XFEL in Pohang, South Korea. In addition, an on-going LCLS-II project will upgrade the LCLS to a continuous-wave (CW) superconducting linac, similar to the European XFEL, to produce high-repetition-rate FEL beams up to 1 MHz. The construction of existing and new facilities as well as the upgrade activities were, and continue to be, driven by the great scientific potential of FELs in the frontier research of physics, chemistry, life science, material and energy sciences because of FELs' extraordinary properties, including ultra-short pulse duration, ultra-high peak brightness and nearly full transverse coherence. To maximize the scientific impact of FELs, various photon diagnostic, optics and beam conditioning devices are required to best utilize the beam capabilities.

Inline photon diagnostics such as pulse energy monitors, spectrometers and timing tools are essential in providing important beam properties to the users for properly analyzing experimental data and the operators for FEL tuning and optimization. Similarly, the ability to reduce the FEL pulse energy in a controlled manner is of great importance for helping users to perform alignment, evaluating the threshold for sample damage and measuring nonlinear effects. A FEL monitor or attenuator is required to not only sustain the extremely high peak fluence of single FEL pulses without damage but also to be able to dissipate the absorbed energy sufficiently quickly before the arrival of the subsequent pulse. Especially for soft X-ray FEL beams, the former requirement effectively eliminated all techniques other than using a gas medium. Gas-based intensity monitors and attenuators have been designed and implemented at various facilities including FLASH (Richter et al., 2003; Hahn & Tiedtke, 2007) and LCLS (Hau-Riege et al., 2008; Ryutov et al., 2009). The fast energy dissipation has always been implicitly assumed, but starts to break down as the time separation between pulses decreases for high-repetition-rate FELs.

The impact of the much reduced inter-pulse separation on the performance of these gas devices has been studied recently (Feng et al., 2015a,b, 2016a), revealing nonlinear effects of the precision of intensity measurements for gas monitors as well as the effective attenuation for gas attenuators arising from the FEL-induced density depression or filamentation phenomenon. These simulations, however, only considered heat conduction without including hydro­dynamic motions, and thus are in principle only valid for long time scales determined by the gas thermal diffusivity and the typical length scale of the device, which are typically of the order of tens to thousands of microseconds. In the current study reported here, we have used the Navier–Stokes–Fourier (NSF) formulation summarized by Sone (2007) to simultaneously treat the macroscopic hydro­dynamic motions and the thermal conductions by the gas molecules to elucidate the short time-scale (tens of nanoseconds to many microseconds) evolution of the density and temperature distributions in an argon-gas filled attenuator under high-repetition-rate operation. The rapid formation of the density filament was clearly shown, in a single-pulse input, to be driven by hydro­dynamic motions, with its time scale set by the speed of sound, and its subsequent slow decay facilitated by thermal transport. In the limit of a CW input, the hydro­dynamic motion becomes negligible and the steady-state density and temperature distribution agree quite well with those obtained by thermal conduction calculations (Feng et al., 2015a,b, 2016a). For a large-aspect-ratio device such as a long attenuator and when energy deposition is mostly concentrated at the entrance, our full three-dimensional simulation also reveals appreciable gas advection and thermal conduction in the longitudinal direction and at even longer time scales of tens of milliseconds.

2. Previous thermal conduction studies of gas devices

In previous thermal conduction studies, the performances of both gas attenuators and detectors were studied. For gas energy detectors, the pulse energy measurement was found to be impacted substantially by the total energy absorbed in the interaction region, especially when there are shot-to-shot fluctuations in the input FEL per pulse energy (Feng et al., 2015b). Even for a monitor with a perfect detection scheme, there is an intrinsic measurement uncertainty arising from the finite FEL beam power that must be absorbed for making the measurement. The standard deviation of the measured pulse energy increases nonlinearly with the absorbed beam power. For the detector geometry, beam parameters and cooling scheme that were considered, the absorbed power must be kept sufficiently low to achieve the required precision, typically at a 1% level. In the case of gas attenuators, the actual achieved attenuation at a given cell pressure deviates significantly from that expected for a low-repetition-rate operation, and cannot be easily predetermined (Feng et al., 2016a). This deviation is nonlinear and depends on many parameters, including the average input power, desired attenuation, repetition rate, geometries of both the beam and cell, and gas properties including the thermal diffusivity and mass attenuation coefficient. Operationally, the actual attenuation must be tabulated first as a lookup table with a given set of beam parameters for subsequent use, rather than being simply dialled up by a trivial calculation from the equilibrium gas pressure. Moreover, the shot-to-shot intensity fluctuations of the FEL beam induce additional uncertainty in the actual attenuation by as much as 75% for the cell geometry and beam parameters considered (Feng et al., 2015a).

The performance dependence on the FEL beam power in both gas monitors and attenuators stems from the so-called `gas filamentation effect' previously observed in an optical femtosecond laser experiment (Cheng et al., 2013), whereby the proceeding pulses deposit energy into the beam path and surrounding volume, which in turn causes the gas to heat up and expand, leaving behind a density depression that the trailing pulses will then experience. The magnitude and spatial extent of the density depression evolve continuously and will eventually approach a steady-state distribution that repeats itself after each pulse, if the pulses have a constant energy and are equally spaced in a full duty-cycle mode (Feng et al., 2016a). If, however, the pulse energy fluctuates as in a FEL beam, leading to variations in the deposited energy, this in turn would make the density depression not repeat itself but rather `fluctuate' about the steady-state distribution.

The underlying mechanism of the gas filamentation effect lies in the finite thermal diffusivity of the gas medium; any thermal disturbances in the gas will take a finite time for it to relax back to its starting conditions. For a low-repetition-rate FEL such as the LCLS at 120 Hz, the inter-pulse time of 8.3 ms turned out to be sufficiently long that, at typical operating conditions for both the gas attenuators and monitors, the filamentation effect did not manifest experimentally and consequently was never investigated until recently (Feng et al., 2016a). For the LCLS-II soft and hard X-ray FEL beams, the minimal inter-pulse time of 10 µs at 100 kHz repetition rate was clearly shown to be too short for the gas devices to thermally recover (Feng et al., 2015a,b). If the inter-pulse time is reduced even further, whether as a result of operating at a partial duty cycle or simply increasing to a higher repetition rate, the filamentation effect is expected to be more severe as the on-cycle running average of the beam power becomes higher (Feng et al., 2016b).

In the current study reported here, a more complete Newtonian fluid dynamics simulation was performed in order to understand the short time-scale (microseconds and longer) evolution of the density and temperature distributions, by including both hydro­dynamic motions (flows) of the gas molecules and (diffusive) thermal conductions in a finite-volume calculation, similar to what was done previously considering a fully absorbing boundary condition and taking into account no viscosity effect for a single-pulse femtosecond laser filamentation experiment (Cheng et al., 2013). After introducing the differential equations based on the NSF formulation of the fluid dynamics and establishing the proper boundary conditions for our specific gas attenuator geometry in §3, single-pulse calculations are presented in §4.1. The initial density depression (or filament) within a few microseconds was clearly shown to emanate from outward shock wave motions driven by the pressure gradient, not as an ad hoc consequence of the pressure equilibration in previous studies; the filament formation is accompanied by a precipitous temperature drop, also in contrast to the assumption of an isothermal process made previously [see, for example, §2 of Feng et al. (2016a)]. However, the ensuing temperature evolution and the filament recovery on time scales of many microseconds or longer is predominantly determined by thermal conduction and the associated back flow, with the additional complications stemming from wave reflections from the cell walls. If the inter-pulse separation is sufficiently short such that the filament cannot fully recover, the filament deepens progressively as the trailing pulses remove additional gas particles as demonstrated in the multi-pulse calculations in §4.3. Since the rate of the shock wave removal decreases while that of the (thermal conduction) back flow increases as time elapses, the two competing mechanisms ultimately bring about a dynamic balance within each pulse cycle and allow the density and temperature evolution to establish asymptotically a repeating pattern, which oscillates about the steady-state distributions presented in §4.2 for the case of a CW input of an equivalent average power in which the effects of both the shock waves and thermal back flow become infinitesimally small and the entire gas medium becomes quiescent.

3. NSF formulation of fluid dynamics in gas devices

Our current simulation applies the NSF formulation (Sone, 2007) to a gas attenuator for an X-ray FEL beam under high-repetition-rate operation. The construction and operation of the attenuator is shown schematically in Fig. 1, where a main gas tube of length L_p and diameter 2R_p (<< L_p) is filled with the gas of choice (e.g. argon) through the gas inlet and regulated to an equilibrium pressure P before the X-ray FEL beam is switched on. The upstream and downstream differential pumping stages are used to maintain the required ultra-high vacuum level (1 × 10⁻⁹ Torr) outside of the cell. The outer wall and two end walls are controlled to a constant temperature, i.e. T₀ = 300 K. The FEL beam of size 2R_b in FWHM enters the cell in a pulsed pattern evenly separated by Δt, ranging from one to many microseconds, or equivalently at 1 MHz or lower repetition rate. The X-ray photon energy is 230 eV, and the pulse energy is assumed to be constant at 2 mJ per pulse in the following single-pulse and 100 kHz multi-pulse simulations and will be scaled accordingly to keep the same average power in additional multi-pulse simulations at different frequencies including a CW case.

Figure 1 Schematic of a typical gas device, consisting of a gas tube of length Lp and diameter 2Rp (<< Lp), which is differentially pumped, and filled with a gas and regulated to pressure P. The cell walls are temperature controlled to T0 = 300 K. The FEL beam of size 2Rb in FWHM enters from the left in a pulsed pattern.

The hydro­dynamic motion and thermal conduction of the gas are assumed to be the dominant processes occurring on time scales of tens of nanoseconds, microseconds and longer, which are much greater than those of the FEL–gas interactions including energy deposition via photoionization (< a few to tens of femtoseconds), and the subsequent thermalization from the electronic to the kinetic degrees of freedom of the gas particles (< tens to hundreds of picoseconds) (Cheng et al., 2013). Furthermore, after many nanoseconds the gas is assumed to be fully recombined and in local thermodynamic equilibrium. The energy deposition profile is presumed to be spatially confined to the region of the photoionization events and is defined by the transverse profile of the FEL beam, largely based on the assumption that the photoelectron mean-free path is much shorter than the beam diameter. The amount of energy deposition in a given volume is determined by the beam intensity and local gas density. The energy dissipation in the gas is assumed to come entirely from thermal exchange at the cell walls, and any radiative processes such as X-ray and other fluorescence are neglected at present but could possibly be included phenomenologically if necessary.

The conservation laws of mass, (linear) momentum and energy in a small volume provide the following differential equations in the absence of external forces:

where n is the number density of gas molecules, v is the hydro­dynamic velocity, p (≡ mv) and ∊ are the (average) momentum and energy per molecule, respectively, m is the molecular mass, σ is the stress tensor, Q is the heat source (nonzero in the case of a CW input, otherwise not present as the energy deposition is treated as the initial condition in a pulsed input, cf. §4.1), j_n is the mass flux, j_k is the energy transfer due to thermal conduction, and t is the time. The constitutive laws for the fluxes and stress are given by

where η and λ are the first and second coefficients of viscosity, P is the thermodynamic pressure, I is the identity matrix, and the superscript in the last term denotes the transpose of a matrix. For simplicity, the Stokes hypothesis is used, λ = 2η/3. P is assumed to follow the ideal gas law, P = nk_BT, with T being the temperature and k_B the Boltzmann constant, which is a valid approximation for this work of a dilute gas. Equation (6) is the Fourier law with k being the thermal conductivity. The total energy is partitioned into thermal and kinetic energy as ∊ ≡ 3k_BT/2 + mv²/2 for a monoatomic gas such as argon. An axisymmetry is also assumed to be in the direction of the beam propagation along the longitudinal axis of the tube, and equations (1)–(6) can be reduced and expressed in cylindrical coordinates (z, r) for simplifying numerical simulations.

The conservation and constitutive equations will be solved by specifying the initial condition such that the entire cell has uniform gas density, pressure, temperature and zero velocity fields. At time t = 0, an incoming FEL beam is impinged at the tube entrance z = 0 and propagates across the tube axis. The beam interacts with the gas molecules and deposits energy primarily via photoabsorption. The beam is assumed initially to have a Gaussian radial profile as given by i(z = 0, r) = i₀exp(−2r²/ω²), where i₀ is the peak fluence and ω is the beam waist size. The initial pulse energy is I₀ = πω²i₀/2. The amount of energy per unit area absorbed by the gas in a thin slice of width Δz in z is given by

where μ is the absorption length calculated for the gas having a number density n₀. The relationship is used to add energy pointwise in the system for the first pulse, and for the trailing pulses as well using an updated number density field n(z, r).

For the boundary conditions, the tube is assumed to be a fully closed system and any exchange of gas molecules through the inlet or the two differential pumping apertures is neglected, certainly for the duration of the simulation time. At all boundary points, the normal component of the velocity is set to zero [automatically the normal component of the mass flux becomes zero due to equation (4)]. The temperature on the tube outer boundary r = R_p is controlled to a constant at 300 K, i.e. a fully absorbing thermal boundary condition. At the tube ends, a thermal insulation condition is assumed such that the normal component of the temperature gradient vanishes. The constant temperature condition on the tube boundary is justified by the consideration that the metal tube wall possesses a much larger heat capacity compared with the amount of energy carried by the gas. Even though the thermal diffusivity of the metal, say, copper, is much lower than that of the dilute gas considered in the simulation, the change in temperature at the tube boundary should still be negligible. In contrast, the insulating condition at the ends may require re-examination. During the simulation, it became apparent that a fully absorbing boundary condition at the tube ends would lead to excessive gas densification, which caused a numerical overflow. A possible reason for this result is the coarse mesh size in z necessitated by the extremely high aspect ratio of the cell. A sufficiently fine mesh at the expense of long computational time should potentially resolve this, which, however, must be deferred to future studies. Another future consideration for improvement is to include as part of the simulation the end apertures with pressure regulation as well as the gas inlet.

The above initial-boundary value problem is then solved numerically by applying a finite volume method (LeVeque, 2002; Versteeg & Malalasekera, 2007). The domain is discretized in a rectangular mesh of finite volumes. The nodes are assigned to the centers of the finite volumes. The fluxes of mass, momentum and energy at the boundaries of the finite volumes are evaluated by finite difference approximations based on the nodal values of density, velocity and energy (temperature). The conservation laws are enforced with these fluxes in a local integral form. By the finite volume method, the conservation laws are ensured between finite volumes, even though the time derivatives and the fluxes are approximate. For the sake of brevity, the standard procedure of the method is not described herein.

4. Simulation results and findings

Our numerical simulations were carried out for argon gas filling the tube to an initial pressure of 2.5 Torr (= 0.333 kPa) and an initial temperature of 300 K, having an initial density of 8.045 × 10⁴ molecules µm⁻³. The viscosity and thermal conductivity of argon are given by η = η₀(T/T₀)^1/2 and κ = κ₀(T/T₀)^1/2, with η₀ = 2.23 × 10⁻⁵ Pa s, κ₀ = 0.0173 W m⁻¹ K⁻¹, T₀ = 293 K. The tube was set to be 4 m in length and 10 mm in radius. The FEL beam has a Gaussian transverse profile with a FWHM of 2 mm. The specific X-ray absorption length of 400 mm was used, corresponding to that of argon at 2 Torr and 293 K and for a photon energy of ∼230 eV. The two-dimensional domain (due to the axisymmetry) was discretized with an adaptive, rectangular mesh. Along the radial direction, a total of 125 nodes were used, with 25 uniformly distributed over the beam half-width of 1 mm, and 100 nodes in a parabolic distribution over the rest of the 9 mm area (more concentrated towards the center and the boundary). Along the longitudinal direction, 40 nodes were used, in a linear distribution with a grading factor of 1.1 from the entrance to the other end. The mesh was sufficient to capture most of the physical processes of the present interest, except at the two ends of the tube, which will be the subject of our future study as discussed briefly above. The results of a single pulse, a CW input and multiple pulses at three selected frequencies of 100 kHz, 200 kHz and 1 MHz are presented, and are followed by discussions. The size of the time increment used was 20 ns for the single pulse and 1 MHz multiple pulse simulations, and 40 ns for the 100 and 200 kHz cases, and 100 ns for the CW case.

4.1. Results of a single-pulse input

We first describe the results from the simulation of a single-pulse input with a pulse energy of I₀ = 2 mJ. This is motivated by the goal of understanding the evolution of the thermodynamic quantities such as pressure, temperature and density at both short (<1 µs) and long (>1 ms) time scales. The energy deposition by the FEL beam via the photoionization process is considered instantaneous for the sake of this work as the density of the gas medium remains essentially unchanged from its initial distribution on the time scale of the FEL pulse duration, and can be calculated from the z-dependence of the pulse energy shown in Fig. 2. The beam is essentially completely attenuated at about z = 1.5 m. What happens next to the gas medium on the intermediate time scales is illustrated by Figs. 3–6. In Fig. 3, three snapshots were taken, immediately after the energy deposition t = 0⁺ (= 20 ns) and at 20 µs and at 40 µs afterwards. In the left-hand column, the profiles of the density distribution in (r, z) and evolution in t are shown, having been color-coded by the temperature; whereas in the right-hand column, the profiles of the corresponding velocity component v_r is coded by the radial stress component, −σ_r. The time evolution of the density (coded by temperature) of the gas tube (for z up to 2 m) is given in Fig. 4 for points only on the beam axis. The time evolution of the density and the radial component of the velocity field v_r for points in the radial direction r at a given z = 36 mm are depicted in Fig. 5, coded by temperature and radial stress component, −σ_r, respectively. In Fig. 6 the variation of density, temperature, radial velocity component v_r, longitudinal velocity component v_z and radial stress component −σ_r, in time at z = 36 mm and selected r = 0, 1, 5 and 10 mm are shown in semi-logarithmic plots. The variation of the longitudinal stress component −σ_z is very similar in magnitude to that of −σ_r and hence is not plotted.

Figure 2 Pulse energy as a function of the location z at the start of the simulation t = 0^+, exhibiting the exponential decay commonly assumed and allowing the energy deposited into gas in a finite volume at a given location to be calculated. This is the case of the only pulse in a single-pulse simulation or of the first pulse in a multi-pulse simulation with an initially uniform cold gas in the cell.

Figure 3 Three snapshots taken (a) immediately after the energy deposition t = 0^+, and (b) at 20 µs and (c) at 40 µs afterwards. In the left-hand column, the profiles of the density distribution in (r, z) and evolution in t are shown, being color-coded by the temperature in K. In the right-hand column, the profiles of the corresponding radial velocity component vr is coded by the radial pressure component (−σr) in kPa.

Figure 4 The intermediate to long time scale (up to 400 µs) evolution of the gas density on the beam axis (r = 0) at different z locations as color-coded by the temperature in K. The modulation of the density superposing on the gradual recovering towards the asymptote is the result of the shock waves transporting gas particles away from and back to the beam axis.

Figure 5 Time evolution of the (a) density and (b) radial velocity component vr at z = 36 mm from the entrance as color-coded by the temperature in K and the radial pressure component (−σr) in kPa, respectively. The density modulation on top of a gradual recovery is again seen clearly at r = 0 at different time points up to 400 µs.

Figure 6 Time evolution (0 to 1000 µs) of the (a) density, (b) temperature, (c) radial velocity component vr, (d) longitudinal velocity component vz and (e) radial pressure component (−σr) at z = 36 mm from the entrance of the gas cell at different radial points at r = 0, 1, 5 and 10 mm from the axis. The longitudinal component of the velocity is much smaller in magnitude than its radial counterpart. Both the temperature and density decay exponentially with time as indicated by the roughly straight lines after the 3 µs mark on the semi-log plots if the shock-wave-related wiggles are ignored.

Immediately after the passage of the FEL pulse, a high-temperature and high-pressure condition within the beam path is built up as shown in Fig. 3(a), which launches a shock wave propagating outwards from the tube center as indicated by the positive value of the radial component of the velocity v_r in Fig. 3(b). The shock wave is reflected at the tube wall, and travels back towards the center as shown by the negative value of v_r in Fig. 3(c). The driving force of the shock wave is the excess pressure or non-vanishing stress, −σ_r, whose maximum magnitude coincides with the shock front and remains positive irrespective of the direction of the propagation. The speed of the shock wave is estimated to be approximately 300 m s⁻¹, and is highest near the entrance (z = 0) where the stress (pressure gradient) is the greatest, and decreases for larger z locations and tends to the limit of that of the (acoustic) sound for a uniform gas. At a longer time scale, as shown in Fig. 4, the shock wave bounces back and forth, transporting energy and gas particles as shown by the modulating on-axis density while it is gradually recovering. This is due to the fact that the tube wall is thermally conductive; the shock wave intensity is much reduced upon reflection. However, the shock wave can last for several additional cycles and can backfill the beam center as it propagates inward. While the shock waves do transport energy and gas particles, the density depression or filament remains and the temperature at the core stays relatively high, indicative of an equilibrium condition of the pressure throughout the gas volume, which in turn would sustain a continuing process of thermal conduction to be discussed later.

In Figs. 5 and 6 the same process is elucidated at a different position z = 36 mm from the entrance of the gas cell, and the time evolution (0 to 1000 µs) of the density, temperature, velocity components v_r and v_z and radial stress component −σ_r at different radial points at r = 0, 1, 5 and 10 mm from the axis are shown. The shock wave emission finishes in about 3 µs from t = 0⁺ and, during its initial launch from the center, the on-axis temperature drops precipitously from an initial ∼2350 K (upon adiabatic FEL heating) to ∼900 K as required by the conversion of thermal energy to hydro­dynamic motions. This is precisely when the pressure simultaneously drops from 2.28 kPa (again upon adiabatic FEL heating) to effectively the equilibrium starting value of 0.333 kPa before the FEL arrival. The gas filament or density depression is formed at the completion of the first wave as it carries away the gas particles from the center. The filament is the deepest with a density only about a quarter of the initial value at z = 36 mm, and decreases at larger z positions. The maximum magnitude of the radial component v_r of the velocity reaches above 200 m s⁻¹, gained from the thermal energy. The axial component v_z of the velocity is much smaller in magnitude than its radial counterpart, indicating a long time scale for the equilibration along the longitudinal axis.

What proceeds from this point forward in time (t > 3 µs) is dominated by thermal conduction, and hydro­dynamic motions become secondary as radial pressure becomes nearly equilibrated. The time evolution of density and temperature can be approximately calculated by using the Fourier law as described previously (e.g. Cheng et al., 2013; Feng et al., 2016a), however, starting with a smaller temperature gradient due to the shock wave emission.¹ Thermal conduction transfers thermal energy outwards towards the cooler wall, lowering the temperature, which in turn reduces the pressure at the core. This negative pressure gradient then drives the gas particles backwards towards the center, allowing the filament or density depression to recover. It will take about 1 ms for this whole recovery process to effectively complete, with the time scale dictated by the overall dimensions of the gas tube and thermal conductivity of the gas species. As can be seen in Fig. 6, both the temperature and density decay approximately exponentially with time indicated by the roughly straight lines on the semi-log plots if the shock-wave-related wiggles are ignored.

4.2. Results of a CW input from cold condition

Before studying the effect of the shock wave under the irradiation of repeated FEL pulses, simulations were performed at an increasingly high repetition rate with diminishingly small pulse energy while maintaining a constant time-averaged power. As such, the amplitude of the shock wave per pulse becomes infinitely small, and the evolution of the temperature and density distributions approach asymptotically towards, as expected, those calculated for a CW input of the equivalent power. The case of a CW input at 200 W starting from a cold condition was considered, as the limit to later simulations of the same time-averaged power but at finite repetition rates. The energy is added to the system continuously through the last Q term in equation (3).

Some results are plotted in Figs. 7–11. Fig. 7 shows the FEL power attenuation and the corresponding gas density profiles along the z-axis at t = 0⁺ (= 100 ns), 0.1 ms and 40 ms. The power decay behavior changes from being exponential, characteristic of a cold gas, to being nearly linear at later times as the gas filament (void) develops. Fig. 8 shows a 2D snapshot of the density distribution color-coded by temperature at t = 20 ms. The filament or a region of density depression is deepest at the center of the entrance, measuring only about one-fifth of the asymptotic value near the wall and towards the tube end. Fig. 9 shows the time evolution of the on-axis density profile color-coded by temperature and of the longitudinal velocity component v_z profile color-coded by the corresponding pressure. Fig. 10 shows the time evolution of the on-axis density profile color-coded by temperature and of the radial velocity component v_r profile color-coded by radial pressure at z = 36 mm, and Fig. 11 shows the time evolution of the density, temperature, both velocity components v_r and v_z, and radial pressure at z = 36 mm and at different radial points at r = 0, 1, 5 and 10 mm.

Figure 7 FEL power as a function of the location z at different time points t = 0^+, 0.1 ms and 40 ms, with the decay behavior changing from being exponential at the start to being nearly linear at longer times. The on-axis density as a function of z is shown by the dashed lines at corresponding times.

Figure 8 Density distribution in the gas cell color-coded by the temperature in K at t = 20 ms. The filament or a region of density depression is deepest at the center of the entrance, measuring only about one-fifth of the asymptotic value near the wall and towards the tube end.

Figure 9 On-axis density profile color-coded by the temperature in K as a function of time up to 10 ms, and (b) the longitudinal velocity component vz color-coded by the longitudinal pressure component (−σz) in kPa, indicating a continuing push of gas particles towards the tube end.

Figure 10 Time evolution of the (a) density and (b) radial velocity component vr at z = 36 mm from the entrance as color-coded by the temperature in K and the radial pressure component (−σr) in kPa, respectively.

Figure 11 Time evolution (0^+ to 10^5 µs) of the (a) density, (b) temperature, (c) radial velocity component vr, (d) longitudinal velocity component vz and (e) radial pressure component (−σr) at z = 36 mm from the entrance of the gas cell and different radial points at r = 0, 1, 5 and 10 mm from the axis.

It can be seen from Fig. 10 that the CW beam generates an initial shock wave, similarly to that in the above single-pulse case, but at much less intensity. The shock wave quickly fades away, leaving behind only small ripples bouncing back and forth within the tube. The filament is formed in about 10 µs, with a central depth that is one-fifth of the initial density (at location z = 36 mm) as shown in Fig. 11. The filament is deeper but it takes a longer time to form compared with the single-pulse case. After the initial depletion, the density at the core is stable, while the gas particles are being continuously pushed from inside out and piled up against the tube wall. This entire process takes about 3 ms, after which the density near the wall will eventually drop due to the transport towards the cold tube end driven by the longitudinal pressure gradient and heat conduction (Fig. 9).

The temperature gradually increases from the initial cold condition from the continuing FEL heating, and reaches up to nearly 2300 K at the center entrance after about 3 ms, and then starts to decrease due to the axial process. The pressure variation at longer time scales also indicates the same characteristic behavior as those of the density and temperature, which is not surprising considering that they are related through the ideal gas law when a significant velocity field is no longer present. The pressure and velocity evolutions indicate that the dynamic wave event is present only from the start to about 10 µs, a bit longer than 3 µs of the single-pulse case due to energy being slowly trickled in rather than at once. Furthermore, it can be seen from Fig. 7 how the gas filamentation affects the attenuation as it develops. Lastly, there is a substantial longitudinal velocity component v_z deeper inside the tube, which can reach up to about 70 m s⁻¹ at the 8 ms mark and then gradually decreases indicating a nontrivial and slower but active transport process along the z-axis.

These CW results clearly demonstrate three distinct processes on separate time scales: (i) the central filament formation by shock waves occurring at a few microseconds; (ii) the radial transport process establishing density temperature and pressure equilibrium through thermal conduction from 4 to 8 ms; (iii) the axial transport process through thermal conduction and convection approaching the final distribution at about 40 ms. The CW results are largely consistent with those presented in the previous works (e.g. Feng et al., 2016a, and references therein) simulated using only the thermal conduction without hydro­dynamic motions. Only the current study, however, can provide the real mechanism of how the gas filament is formed and on what time scales.

4.3. Results of a multi-pulse input

Finally, we turn to simulate the case of a gas device under repeated FEL pulse irradiation to further elucidate the effect of the shock waves, especially when the pulse repetition rate and shock wave time scales are comparable. Simulations were carried out at three different frequencies, 100 kHz, 200 kHz and 1 MHz. The per-pulse energies were 2 mJ, 1 mJ and 0.2 mJ, respectively, corresponding to the same time-average power of 200 W, which is the same for the CW simulation in §4.2. Some selected results are plotted in Figs. 12–16. Fig. 12 shows the effective FEL attenuation of the gas cell following the sequence of pulses for different pulse frequencies at 100 kHz, 200 kHz and 1 MHz, exhibiting the exponential decay commonly expected for a cold gas at 0⁺ (= 40 ns for the 100 kHz and 200 kHz cases, and 20 ns for the 1 MHz case), but deviating significantly for subsequent pulses and approaching an asymptote. Fig. 13 shows the time evolution of the on-axis density profile color-coded by temperature and of the radial velocity component v_r profile color-coded by radial pressure at z = 36 mm for a repetition rate of 100 kHz. Figs. 14, 15 and 16 show the time evolution of the density, temperature, both velocity components v_r and v_z, and radial pressure at z = 36 mm and different radial points at r = 0, 1, 5 and 10 mm for the three repetition rates 100 kHz, 200 kHz and 1 MHz, respectively.

Figure 12 Pulse energy as a function of the location z at the start of the simulation t = 0^+, and following repeated pulses at three different frequencies, (a) 100 kHz, (b) 200 kHz and (c) 1 MHz, exhibiting the exponential decay commonly assumed for a cold gas at 0^+, but deviating significantly for subsequent pulses and approaching an asymptote.

Figure 13 Time evolution of (a) the on-axis density profile color-coded by the temperature in K, and (b) the corresponding radial velocity component vr profile color-coded by the radial pressure component (−σr) in kPa at z = 36 mm at a repetition rate of 100 kHz.

Figure 14 Time evolution (0^+ to 2000 µs) of the (a) density, (b) temperature, (c) radial velocity component vr, (d) longitudinal velocity component vz and (e) radial pressure component (−σr) at z = 36 mm from the entrance of the gas cell and different radial points at r = 0 mm (red), 1 mm (cyan), 5 mm (orange) and 10 mm (pink) from the axis for a repetition rate of 100 kHz. The thick dark lines are the CW results for comparison.

Figure 15 Time evolution (0^+ to 2000 µs) of the (a) density, (b) temperature, (c) radial velocity component vr, (d) longitudinal velocity component vz and (e) radial pressure component (−σr) at z = 36 mm from the entrance of the gas cell and different radial points at r = 0 mm (red), 1 mm (cyan), 5 mm (orange) and 10 mm (pink) from the axis for a repetition rate of 200 kHz. The thick dark lines are the CW results for comparison.

Figure 16 Time evolution (0^+ to 2000 µs) of the (a) density, (b) temperature, (c) radial velocity component vr, (d) longitudinal velocity component vz and (e) radial pressure component (−σr) at z = 36 mm from the entrance of the gas cell and different radial points at r = 0 mm (red), 1 mm (cyan), 5 mm (orange) and 10 mm (pink) from the axis for a repetition rate of 1 MHz. The thick dark lines are the CW results for comparison.

In these cases of different pulse frequencies but having the same average power, the amount of energy deposited on the tube axis per pulse is quite different. It results in very different spikes and fluctuations, especially of temperature, radial pressure component and radial velocity component. Overall the fields being established by the consecutive pulses may be described as an evolving baseline superimposed by fluctuations due to the discrete pulses. The baseline trend largely follows that of the CW results as the limiting behavior. The amplitudes of the fluctuations/spikes are determined by the amount of energy each pulse deposits. When the frequency increases and when the energy of each pulse diminishes (but having a fixed time-averaged power), the spikes disappear and the multi-pulse solutions approach those of the true CW input.

Since the baseline temperature in the CW case can reach up to 2300 K and a single pulse of energy 2 mJ can raise the temperature by about 2000 K, the maximum temperature spike at the center of the tube entrance (z, r = 0) can reach as high as 3700 K (at the t = 2 ms mark, and even higher at later times) despite the release of energy by the shock waves, in the case of the 100 kHz frequency as shown in Fig. 14. When the pulse energy is reduced to 1 mJ in the 200 kHz case, the temperature jump by a single pulse is reduced by a factor of two, as shown in Fig. 15. At the highest frequency of 1 MHz simulated, the spikes nearly all vanish as shown in Fig. 16. The fluctuations in the baselines for the CW case are due to the cylindrical shock waves bouncing back and forth inside the tube. The lifetime of these rippling waves is determined by the momentum diffusivity, η/nm, where η is the first coefficient of viscosity, n is the number density and m is the molecular mass.

The filament forms in a few microseconds, mainly by the first shock wave carrying particles away. The subsequent pulses continue to carry more particles away, but at a much lower rate, simply because the density is much lower than initially, as shown in Figs. 13–16. As can also be seen in Fig. 12, it takes more pulses to reach the steady-state attenuation at higher frequencies, as more pulses of lower energy are required to heat up the filament region and move the particles out. By comparing these three cases of different frequencies and the CW case, it is evident that the radial component v_r of the velocity is purely a result of the shock wave. In contrast, the other state variables, the density, temperature and pressure, are more correlated, both in the initial stage of filament formation and further development at later times.

To further illustrate the primary role that the hydro­dynamic processes play in creating the density depression, the time evolution of the density, temperature and radial pressure at intermediate time scales from the first time step 0⁺ to 50 µs or so for different repetition rates of 100 kHz, 200 kHz, 1 MHz and the CW case is overlaid in Fig. 17. After each pulse, both temperature and pressure go up abruptly and the density depression is driven by the pressure gradient but occurs gradually as dictated by the speed of sound; in contrast, in the previous thermodynamic simulations (Feng et al., 2015a,b, 2016a), the density depression was assumed to happen instantaneously by virtue of the assumption of local thermodynamic equilibrium and the lack of a pressure gradient. It must be pointed out that the current hydro­dynamic CW simulation is only equivalent to the previous thermodynamic simulation after reaching steady state when the pressure gradient ceases to change. After the onset of the FEL, a pressure gradient still builds up in addition to the temperature rise as shown in Fig. 17(c) for t < 10 µs, which consequently drives the density lower, albeit without any distinctive wave-like motions.

Figure 17 Time evolution of intermediate time scale (0^+ to 50 µs) of the (a) density, (b) temperature, (c) radial pressure component (−σr) at the center r = 0 and z = 36 mm from the entrance of the gas cell at different repetition rates of 100 kHz (red), 200 kHz (cyan), 1 MHz (blue) and CW (black).

5. Summary and discussions

In summary, by using the classical NSF formulations, we were able to investigate the hydro­dynamic behavior of a gas device under repeated irradiation of ultrafast X-ray FEL pulses at high frequencies. It was found that the shock waves play the primary role in creating a density depression (or filament) by transporting gas particles away from the X-ray laser–gas interaction region, where large pressure and temperature gradients have been built upon the initial energy deposition via X-ray photoelectric absorption and subsequent thermalization. The higher the per pulse energy, the faster the formation of the filament, but eventually reaching a steady state whose distributions of state variables resemble very closely those of the CW input as expected, with the additional ripples stemming from the discrete pulses. In the meantime, there is a counter flow of particles driven by the pressure gradient created in the wake of outward thermal conduction, which itself originates from the high-temperature condition in the filament region even after the shock wave has moved out. The development of the filament will cause the attenuation of the FEL beam to deviate from what is expected for a cold gas medium and to depend on the sequence of the pulse arrival and per pulse energy, which is an operating difficulty that must be mitigated either by precisely knowing the energies of all preceding pulses or by deploying an entirely different device that is less affected by this filamentation phenomenon. Measuring precisely the pulse energy using a gas device creates a dichotomy for it too may suffer the similar FEL induced uncertainty.

It must also be noted that the current study and prior works (e.g. Feng et al., 2016a, and references therein) are all based on the assumption that the energy deposition profile by the FEL pulses is limited to the spatial extent of the beam, to which further studies focusing on the FEL–gas interaction at much shorter time scales and in the plasma regime may suggest otherwise. The authors are actively engaging in such investigations both in simulations and in conducting experiments.

Footnotes

¹Surprisingly, because of the temperature-dependent thermal conductivity of the gas, the decay of the temperature gradient in the current study is slower than that calculated previously with a much higher starting temperature, thus arriving at the same asymptotic temperature at approximately the same time.