# Higher order harmonic generation and strong field ionization with Bessel-
Gauss beams in a thin jet geometry

Michael Davino  1 Adam Summers  2,3 Tobias Saule  1 Jan Tross  2 Edward
McManus  1 Brandin Davis  1 and Carlos Trallero-Herrero1,2  1Department of
Physics, The University of Connecticut, 196 Auditorium Road, Storrs, CT 06269,
USA  2J. R. Macdonald Lab., Kansas State University, 116 Cardwell Hall,
Manhattan, KS 66506, USA  3CFO, Institut de Ciencies Fotoniques, Mediterranean
Technology Park, Av. Carl Friedrich Gauss 3, 08860, Castelldefels, Barcelona,
Spain<EMAIL_ADDRESS>
###### Abstract

A promising alternative to Gaussian beams for use in strong field science is
Bessel-Gauss (BG or Bessel-like) laser beams as they are easily produced with
readily available optics and provide more flexibility of the spot size and
working distances. Here we use BG beams produced with a lens-axicon optical
system for higher order harmonic generation (HHG) in a thin gas jet. The
finite size of the interaction region allows for scans of the HHG yield along
the propagation axis. Further, by measuring the ionization yield in unison
with the extreme ultraviolet (XUV) we are able to distinguish regions of
maximum ionization from regions of optimum XUV generation. This distinction is
of great importance for BG fields as the generation of BG beams with axicons
often leads to oscillations of the on-axis intensity,which can be exploited
for extended phase matching conditions. We observed such oscillations in the
ionization and XUV flux along the propagation axis for the first time. As it
is the case for Gaussian modes, the harmonic yield is not maximum at the point
of highest ionization. Finally, despite Bessel beams having a hole in the
center in the far field, the XUV beam is well collimated making BG modes a
great alternative when spatial filtering of the fundamental is desired.

††journal: osajournal††articletype: Research Article

##  1 Introduction

With the discovery of high harmonic generation (HHG) in the late 1980s and
1990s [1, 2, 3, 4], and the proceeding evolution of HHG based light sources
[5, 6], a variety of new fields of physics were developed. HHG has allowed for
the ability to produce pulses of coherent radiation with pulse durations of
less than 100 attoseconds [7] and with photon energies exceeding 1keV [8].
These light sources unlocked the ability to probe ultrafast processes, such as
photoionization [9, 10, 11] chemical bond rearrangement [12, 13, 14] and
electronic dynamics in solid-state materials [13, 15], on their natural time
scales. Additionally, the multi-octave spectral bandwidth inherently present
in XUV and X-ray attosecond pulses has given the ability to probe matter with
elemental specificity across a wide range of absorption edges [16]. As such,
the development of HHG light sources with higher flux, higher photon energy
and shorter pulse duration has received immense and continued interest [17].
A multitude of approaches to produce more advanced HHG light sources have been
pursued. These have included novel driving laser sources, target geometries,
and generation media [6, 5, 8, 18]. However, comparatively little work has
been undertaken to explore the advantages of using non-Gaussian mode profiles
in the driving laser source.
In this work we explore the use of non-Gaussian beams (i.e. non-Gauss-Hermite)
in HHG. Here we specifically investigate using Bessel-Gauss (BG) beams [19,
20], also known as Bessel-like [21, 22], for the generation of the XUV
radiation. In the past we have characterized the properties of such beams and
their relevance for strong field science [23].
BG beams have a series of enticing properties [24, 25]. Among the most notable
is the ability to overcome limitations of the Rayleigh range [26, 27]. Bessel-
like beams maintain focused, on-axis peak intensity for much longer distances
than focused Gaussian beams [23]. This allows for the use of an extended
generation target geometry as well as potentially more favorable phase
matching conditions. Finding and establishing the best conditions for phase
matching has been one of the most active areas of high harmonic and attosecond
science [28]. The potential phase matching benefits from Bessel-like beams
primarily arise from the difference in on-axis phase front velocity as well as
the extended intensity profile that can be achieved using BG beams for HHG as
opposed to conventional Gaussian focal profiles [18]. Further, recent
theoretical work on two-color gating with BG beams shows promise for the
generation of attosecond pulses of high quality beam profiles [29, 30, 31] or
with extended phase matching conditions [31].
Experimental studies have previously been presented exploiting the generation
of XUV pulses also with a lens-axicon optical system in a semi-infinite gas
cell [32]. It was demonstrated that BG beams are not only capable of producing
XUV pulses but also that they show an increase in XUV flux for some harmonics
as compared to standard Gaussian beams which can be attributed to an
improvement in phase matching conditions [33].
An additional advantage of Bessel-like beams for HHG is self-filtering. In
most experiments, utilizing the radiation produced by HHG requires the
elimination of the fundamental beam. This is typically accomplished by
employing a series of very thin filters. While these filters block
significantly more IR than XUV or X-ray radiation, they still reduce the total
flux of the HHG beam. Using transmissive filters for filtering the driving
radiation from the HHG radiation also presents a potential problem as HHG
pulse durations approach the zeptosecond regime. The limitation with this
technique is clear with the realization that the optical path of 100 zs is 30
pm, or roughly two atomic layers of $\alpha$-quartz. Such variation could
happen in thin metallic filters due to strain in vacuum or heat because of
laser interaction. While some promising work in self-filtering schemes has
been presented in using annular beams [34], these beams lack the substantial
tunability of BG beams produced by a lens-axicon system which includes an
extended effective Rayleigh range, focal distance, and spot size [23].
Further, the focusing of annular beams for HHG results in a decrease in
harmonic yield relative to a purely Gaussian beam whereas BG have shown an
increase in harmonic yield in a semi-infinite geometry [32].
In this paper we show the first HHG obtained from a capillary (nozzle) where
the fundamental is focused using a lens-axicon optical system. The capillary
produces a thin gas target region and allows us to explicitly study how HHG
flux evolves along the focal profile of the BG beam. This is in contrast to
extended gas cell type geometries that do not allow for position sensitive
yield measurements. Measuring the local yield for a more complete
understanding the interplay between laser intensity and HHG phase matching.
This is vital information when working with the complex focal profiles present
in Bessel-Gauss beams. The total ionization yield is simultaneously measured
with the harmonics, giving the ability to disentangle the focal intensity
profile from local phase matching conditions. We have shown in the past that
such parallel measurements provide an unambiguous determination of the maximum
peak intensity and the relative position of the optical beam and gas jet.
Mapping the HHG yield from such beams is an important step in implementing BG
beams for high-flux HHG sources.

##  2 Experimental methods and results

Figure 1: Experimental setup. The output of a 1 kHz Ti:Sapph ($\lambda$ = 795
nm, $\tau$ = 30 fs) laser system is incident on a tip block, lens, and axicon
producing a BG beam to be focused into an Argon (Ar) gas jet. The XUV beam
produced from the gas jet is inherently isolated from the fundamental light as
the on-axis BG beam vanishes in the far field leaving only the XUV to enter
the spectrometer comprising of a slit, focusing grating, MCP and phosphor
screen. The phosphor is imaged by an external camera. To attain a profile of
the XUV beam the grating and slit are translated perpendicular to the beam.
Also shown is a simulated intensity profile of a femtosecond BG beam.

A schematic of the experimental setup is shown in Figure 1. The laser source
used for this study is a 1 kHz Titanium:Sapphire system producing 2 mJ pulses
with a duration of 30 fs and a central wavelength of 795 nm. Next, to prepare
the Bessel-Gauss beam for use in generating XUV, the laser output is sent into
an optical setup using a lens and axicon in combination [23, 32] with a 250 mm
focusing lens set a variable distance L before a one-degree base angle axicon.
This versatile geometry allows for studying the effects of different BG
focusing conditions on the process of HHG. A 2 mm diameter tip block is placed
before the lens and axicon to mitigate the imperfections of the tip of the
axicon (as done by S. Akturk et al. [35]). Both the lens and axicon are
mounted on a computer-controlled stage to allow the focal profile (see inset
of Figure 1) to be scanned through the gas jet. The stage has a 30 mm travel
range allowing the distance from the axicon tip to the gas jet to be varied
between 95 mm and 125 mm. For larger values of L this distance is further
limited by the spot size of the beam on the entrance window to the HHG
chamber, which induces nonlinear effects as the peak intensity increases.
Argon gas is injected into the source chamber as a thin gas jet produced by a
hollow core (250 µm core diameter) capillary mounted on a 3D manipulator
approximately 5 cm from the entrance window of the chamber. The core diameter
of the capillary was chosen to be sufficiently smaller than the length scale
on which oscillations in on-axis intensity occur in BG beams [23] (see Figure
1 inset). The Argon backing pressure for this experimental study was
approximately 2-3 bar. To keep the gas flow low a needle / leak valve is
placed in the gas line, limiting the gas flow to 200 mbar L/sec.
A channel electron multiplier (CEM, channeltron) is placed several centimeters
behind the gas jet and slightly offset to allow for simultaneous measurement
of ionization yield and HHG signal. This technique was used to demonstrate
phase matching conditions in the past [36, 37]. The source chamber is kept at
6.66 x 10-5 mbar during operation by a 2000 L/s magnetic levitation turbo pump
to prevent reabsorption in the ambient backing pressure [38].
An HHG spectrometer (same as in references [39, 40]) is placed approximately 1
m behind the generation region. Differential pumping tubes between the source
chamber and the spectrometer allow for two orders of magnitude pressure
difference. The spectrometer consists of a slit followed by a curved,
variable-spacing, holographic grating (groove density 1,200 grooves/mm), that
disperses the harmonics across a Z-stack MCP (Photonis, 105 mm by 35 mm) with
a P46 phosphor screen on the output. The slit and grating are mounted on a
translational stage allowing a beam profile of the XUV to be taken (see inset
of Figure 1). A Hamamatsu Orca Flash 2.8, low-noise, 12-bit resolution camera
is used to capture the harmonic spectrum. An example of a raw XUV spectrum is
shown in Figure 2. The integration time was 300 ms and harmonics were observed
up to the 27th order.
The ionization yield and the XUV spectra were taken simultaneously while
varying the nozzle position in a 1D scan along the beam propagation with
respect to the lens and axicon. Figure 3 shows two scans for different
experimental conditions with L = 13 mm (a) and L = 12.3 mm (b). The top plot
is a trace of the ionization yield with the HHG yield attained from
integrating the shown 2D map of the harmonic spectrum in Figure 3 (bottom).
Varying the distance, L, between the lens and axicon alters the focusing
conditions of the BG beam. This can be seen in both the ionization traces as
well as the HHG yield dependence as a function of nozzle position.

Figure 2: Typical spectrum from imaging the MCP/Phosphor screen. The
integration time was 300 ms. Harmonics up to the 27th order were observed
corresponding to 41 eV. Figure 3: From top to bottom; traces of ionization and
HHG yield as well as 2D maps of HHG yield for L = 13 mm (a) and L = 12.3 mm
(b). Images were taken at 37 (a) and 53 (b) different positions along z, the
axis of propagation. At each position the HHG output spectrum was recorded to
produce the 2D map. This map was then integrated to produce the trace of the
full HHG yield over all detected harmonic orders. Vertical and horizontal red
dashed lines in both (a) and (b) indicate the position of maximum HHG yield
and the ionization yield at that position respectively. Black circles indicate
other positions at which the ionization level was the same as at the position
of max XUV flux.

Several striking features are immediately visible. The first is the long range
over which strong-field ionization takes place (> 2 cm in Figure 2 (b)). The
second major feature in the ionization yield is the oscillations as the beam
is scanned through the jet. Measurement of the ionization yield gives us
direct access to the number of ions produced in the strong field. Much like
the intensity profile shown in the Figure 1 inset, oscillations are observed
in both the HHG and ionization yields along the axis of the BG beam. Both the
extended ionization range and the oscillations are features absent from
strong-field ionization using Gaussian beams.
Oscillations are also present in the HHG data. However, for both the L = 13 mm
and L = 12.3 mm cases the HHG yield does not follow with a one-to-one
correspondence to the ionization. In both Figure 3 (a) and Figure 3 (b) we
observe that the ionization yield drops to approximately half the peak value
in the second maximum while the HHG yield drops by an order of magnitude
respectively. Further, horizontal red dashed lines in Figure 3 (a) and (b)
indicate the level of ionization at the focal position of max XUV flux. There
are multiple positions at which this same ionization yield is found (marked by
black circles), but each with significantly less XUV yield. These
discrepancies in correlation between the harmonic yield and ionization yield
point to the presence of different phase matching conditions through different
regions, that are not only a function of the laser intensity or ionization
rate. Considering the expected on-axis intensity oscillations for BG beams,
one might also expect to observe the resulting oscillations in both the XUV
and ionization yield. However, the results shown in Figure 3 indicate an added
degree of complexity to understanding the full picture for HHG.

Figure 4: Simulated intensity profiles of a BG beam from a lens-axicon system
with parameters: 1 mm tip block, L = 12 mm, and a 4 mm FWHM spot size incident
on the axicon that closely resembles the experimental conditions. (a) shows
intensity as a function of the distance from the tip of the axicon (z,
increasing from right to left) and the radial distance. Two line profiles
derived from the distribution in (a) are also shown. Black is an integration
of the intensity in the central spot of the beam profile (i.e. from radius = 0
$\mu$m to the first radial zero of the intensity), and red is the on-axis
lineout (radius = 0 $\mu$m), normalized to their max values. Also shown are
two transverse beam profiles taken at z = 70 mm (b), and z = 190mm (c). These
positions are marked with dashed black lines in (a). The profile at z = 70 mm
corresponds to the peak of the first oscillation in the BG focal region from
the axicon tip where maximum XUV and ionization yield is seen. The profile at
z = 190 mm shows the vanishing on-axis intensity of a BG beam.

To better understand the oscillations in ionization and XUV yield, we have
performed simulations for the three-dimensional propagation of the BG beams
using a quasi-discrete Hankel transform method [41]. The simulation uses
parameters comparable to those used in the experiment. Figure 4 shows the
simulated intensity profile as the beam propagates from right to left in the
figure. The figure also shows the on-axis intensity as a function of
propagation distance (red) as well as the integrated intensity profile to the
first radial zero of the BG transverse profile (black). As expected, both the
integrated and peak intensities exhibit oscillations typical of Bessel-like
beams [35] that closely follow the measured ionization oscillations, however,
the maximum of both curves occur at different points. The integrated intensity
is maximized at the start of the focal region (z $\approx$ 10 mm in Figure 4),
closer to the axicon tip, and degrades rapidly moving away from the axicon. In
contrast, the on-axis intensity is more constant across the focal region and
is maximized towards the end of the region (z $\approx$ 110 mm in Figure 4),
further from the axicon. Comparing experiment to simulation, the HHG and
ionization yields presented in Figure 3 mirror the integrated intensity
profile in Figure 4 rather than the on-axis (peak) intensity profile.
Therefore, we argue that the maximum ionization yield occurs at the first
intensity oscillation in the BG focus from the axicon due to volume averaging
effects between the gas jet and focal profile. The XUV yield follows this as
HHG is macroscopic, coherent process which scales with the number of emitters
squared [37, 36].
With the simulation presented in Figure 4 some light is also shed on the
ionization data previously touched on. Horizontal red dashed lines in Figure 3
mark the ionization level equivalent to the value seen at the position along
the focal region of maximum XUV yield and black circles mark other positions
where there is similar ionization yield. Though there are multiple positions
where this same ionization yield occurs, they vary drastically in XUV yield,
particularly in decreasing at positions further from the axicon tip. To
understand this observation we need to remind the reader the experimental
limitation of our ionization measurements. By using a channeltron we are not
able to measure the ionization fraction. Instead we can only conclude that to
achieve a similar current in the channeltron output but from a smaller
interaction region, a larger fraction of ions need to be generated. Figure 4
shows that the size of the central spot in the BG focal region shrinks moving
away from the axicon while the peak intensity remains largely constant with
each intensity oscillation. This indicates that the ionization fraction
increases with each oscillation in intensity as the peak intensity remains
similar while the spot size shrinks for a constant gas target. Thus, HHG phase
matching is degraded with each oscillation due to the presence of a higher
concentration of ions. Strong field ionization and macroscopic harmonic
simulations are needed to ratify our experimental findings.

##  3 Concluding remarks

We measured higher-order harmonics and ionization yields in unison from
Bessel-like beams in a thin Ar gas jet, which is scanned along the propagation
axis. The BG beams are created with a lens-axicon optical system that provides
flexibility of the spot size, working distance and focal range when compared
with Gaussian optics. The ionization follows the intensity oscillations in BG
beams created by non-perfect axicons and subsequently there are multiple
regions along these oscillations of favorable conditions for HHG. Harmonics
and ionization are both maximum around the first oscillatory point where there
is a lower on-axis peak intensity. We attribute the larger yields to a larger
volume integral at this spatial point. This surprising result is complimented
by numerical simulations of BG beams produced by a comparable lens and axicon
system. The numerical simulations show that the integrated field strength over
the center part of the spatial distribution is indeed larger for the initial
oscillation compared to the spatial oscillation at the highest peak intensity.
Finally, despite the driving field having a donut-like far field beam profile,
the HHG beam is well collimated in the far field. Thus, generation of
attosecond pulses with BG beams presents a novel spatially tunable, self-
filtering approach for the strong field needed for the HHG. Besides removing
the need of spatial filters, the presence of multiple regions for XUV
generation also opens the door for multiple interacting regions as an approach
for the generation of attosecond pulses.

##  4 Acknowledgements

This work was funded by the Air Force Office of Scientific Research, Grant No.
FA9550-17-1-0369. T.S. was partially supported by the Chemical Sciences,
Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office
of Science, U.S. Department of Energy (DOE) Grant No. DE-SC0019098.

##  5 Disclosures

The authors declare no conflicts of interest.

## References

  * [1] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K. Rhodes, “Studies of multiphoton production of vacuum-ultraviolet radiation in the rare gases,” J. Opt. Soc. Am. B 4, 595 (1987).
  * [2] A. L’Huillier and P. Balcou, “High-order harmonic generation in rare gases with a 1-ps 1053-nm laser,” Physical Review Letters 70, 774–777 (1993).
  * [3] J. L. Krause, K. J. Schafer, and K. C. Kulander, “High-order harmonic generation from atoms and ions in the high intensity regime,” Physical Review Letters 68, 3535–3538 (1992).
  * [4] P. B. Corkum, “Plasma perspective on strong field multiphoton ionization,” Phys. Rev. Lett. 71, 1994–1997 (1993).
  * [5] P. Corkum and F. Krausz, “Attosecond science,” Nature physics 3, 381–387 (2007).
  * [6] F. Krausz and M. Ivanov, “Attosecond physics,” Rev. Mod. Phys. 81, 163–234 (2009).
  * [7] T. Gaumnitz, A. Jain, Y. Pertot, M. Huppert, I. Jordan, F. Ardana-Lamas, and H. J. Wörner, “Streaking of 43-attosecond soft-X-ray pulses generated by a passively CEP-stable mid-infrared driver,” Opt. Express 25, 27506 (2017).
  * [8] T. Popmintchev, M.-C. Chen, P. Arpin, M. M. Murnane, and H. C. Kapteyn, “The attosecond nonlinear optics of bright coherent X-ray generation,” Nat. Photonics 4, 822–832 (2010).
  * [9] E. Goulielmakis, Z. H. Loh, A. Wirth, R. Santra, N. Rohringer, V. S. Yakovlev, S. Zherebtsov, T. Pfeifer, A. M. Azzeer, M. F. Kling, S. R. Leone, and F. Krausz, “Real-time observation of valence electron motion,” Nature 466, 739–743 (2010).
  * [10] M. Uiberacker, T. Uphues, M. Schultze, A. J. Verhoef, V. Yakovlev, M. F. Kling, J. Rauschenberger, N. M. Kabachnik, H. Schröder, M. Lezius, K. L. Kompa, H. G. Muller, M. J. Vrakking, S. Hendel, U. Kleineberg, U. Heinzmann, M. Drescher, and F. Krausz, “Attosecond real-time observation of electron tunnelling in atoms,” Nature 446, 627–632 (2007).
  * [11] M. Schultze, M. Fieß, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A. L. Cavalieri, Y. Komninos, T. Mercouris, C. A. Nicolaides, R. Pazourek, S. Nagele, J. Feist, J. Burgdörfer, A. M. Azzeer, R. Ernstorfer, R. Kienberger, U. Kleineberg, E. Goulielmakis, F. Krausz, and V. S. Yakovlev, “Delay in photoemission,” Science 328, 1658–1662 (2010).
  * [12] A. R. Attar, A. Bhattacherjee, C. D. Pemmaraju, K. Schnorr, K. D. Closser, D. Prendergast, and S. R. Leone, “Femtosecond x-ray spectroscopy of an electrocyclic ring-opening reaction,” Science 356, 54–59 (2017).
  * [13] R. Geneaux, H. J. Marroux, A. Guggenmos, D. M. Neumark, and S. R. Leone, “Transient absorption spectroscopy using high harmonic generation: a review of ultrafast x-ray dynamics in molecules and solids,” Philosophical Transactions of the Royal Society A 377, 20170463 (2019).
  * [14] P. M. Kraus, M. Zürch, S. K. Cushing, D. M. Neumark, and S. R. Leone, “The ultrafast X-ray spectroscopic revolution in chemical dynamics,” Nat. Rev. Chem. 2, 82–94 (2018).
  * [15] A. L. Cavalieri, N. Müller, T. Uphues, V. S. Yakovlev, A. Baltuška, B. Horvath, B. Schmidt, L. Blümel, R. Holzwarth, S. Hendel, M. Drescher, U. Kleineberg, P. M. Echenique, R. Kienberger, F. Krausz, and U. Heinzmann, “Attosecond spectroscopy in condensed matter,” Nature 449, 1029–1032 (2007).
  * [16] F. Silva, S. M. Teichmann, S. L. Cousin, M. Hemmer, and J. Biegert, “Spatiotemporal isolation of attosecond soft X-ray pulses in the water window,” Nat. Commun. 6, 6611 (2015).
  * [17] L. Young, K. Ueda, M. Gühr, P. H. Bucksbaum, M. Simon, S. Mukamel, N. Rohringer, K. C. Prince, C. Masciovecchio, M. Meyer _et al._ , “Roadmap of ultrafast x-ray atomic and molecular physics,” Journal of Physics B: Atomic, Molecular and Optical Physics 51, 032003 (2018).
  * [18] G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, “Linking high harmonics from gases and solids,” Nature 522, 462–464 (2015).
  * [19] J. H. McLeod, “The Axicon: A New Type of Optical Element,” J. Opt. Soc. Am. 44, 592 (1954).
  * [20] J. H. McLeod, “Axicons and Their Uses,” J. Opt. Soc. Am. 50, 166 (1960).
  * [21] J. Durnin, J. Miceli, and J. H. Eberly, “Diffraction-free beams,” Phys. Rev. Lett. 58, 1499–1501 (1987).
  * [22] J. Durnin, J. H. Eberly, and J. J. Miceli, “Comparison of Bessel and Gaussian beams,” Opt. Lett. 13, 79 (1988).
  * [23] A. M. Summers, X. Yu, X. Wang, M. Raoul, J. Nelson, D. Todd, S. Zigo, S. Lei, and C. A. Trallero-Herrero, “Spatial characterization of Bessel-like beams for strong-field physics,” Opt. Express 25, 1646 (2017).
  * [24] S. Akturk, “Tailored-beam ultrashort laser pulses,” Quantum. Phys. Lett. 1, 97–112 (2012).
  * [25] D. Mcgloin and K. Dholakia, “Bessel beams: Diffraction in a new light,” Contemp. Phys. 46, 15–28 (2005).
  * [26] V. Jarutis, R. Paškauskas, and A. Stabinis, “Focusing of Laguerre-Gaussian beams by axicon,” Opt. Commun. 184, 105–112 (2000).
  * [27] R. Grunwald, M. Bock, V. Kebbel, S. Huferath, U. Neumann, G. Steinmeyer, G. Stibenz, J.-L. Néron, and M. Piché, “Ultrashort-pulsed truncated polychromatic Bessel-Gauss beams,” Opt. Express 16, 1077 (2008).
  * [28] C. D. Lin, A.-T. Le, C. Jin, and H. Wei, _Attosecond and Strong-Field Physics: Principles and Applications_ (Cambridge University Press, 2018).
  * [29] J. Luo, Q. Cheng, and D. Xu, “Generation of intense isolated attosecond pulses with high spatiotemporal quality by two-color polarization-gating Bessel-Gauss beams,” Opt. Commun. 339, 247–253 (2015).
  * [30] Z. Wang, W. Hong, Q. Zhang, S. Wang, and P. Lu, “Efficient generation of isolated attosecond pulses with high beam quality by two-color Bessel–Gauss beams,” Opt. Lett. 37, 238 (2012).
  * [31] Y. Li, Q. Zhang, W. Hong, S. Wang, Z. Wang, and P. Lu, “Efficient generation of high beam-quality attosecond pulse with polarization-gating Bessel-Gauss beam from highly-ionized media,” Opt. Express 20, 15427 (2012).
  * [32] L. Van Dao, K. B. Dinh, and P. Hannaford, “Generation of extreme ultraviolet radiation with a Bessel-Gaussian beam,” Appl. Phys. Lett. 95 (2009).
  * [33] A. Averchi, D. Faccio, R. Berlasso, M. Kolesik, J. V. Moloney, A. Couairon, and P. Di Trapani, “Phase matching with pulsed Bessel beams for high-order harmonic generation,” Phys. Rev. A 77, 021802 (2008).
  * [34] R. Klas, A. Kirsche, M. Tschernajew, J. Rothhardt, and J. Limpert, “Annular beam driven high harmonic generation for high flux coherent XUV and soft X-ray radiation,” Opt. Express 26, 19318 (2018).
  * [35] S. Akturk, B. Zhou, B. Pasquiou, M. Franco, and A. Mysyrowicz, “Intensity distribution around the focal regions of real axicons,” Opt. Commun. 281, 4240–4244 (2008).
  * [36] A. D. Shiner, C. Trallero-Herrero, N. Kajumba, B. E. Schmidt, J. B. Bertrand, K. T. Kim, H.-C. Bandulet, D. Comtois, J.-C. Kieffer, D. M. Rayner, P. B. Corkum, F. Légaré, and D. M. Villeneuve, “High harmonic cutoff energy scaling and laser intensity measurement with a 1.8 $\mu$m laser source,” Journal of Modern Optics .
  * [37] A. D. Shiner, C. Trallero-Herrero, N. Kajumba, H. C. Bandulet, D. Comtois, F. Légaré, M. Giguère, J. C. Kieffer, P. B. Corkum, and D. M. Villeneuve, “Wavelength scaling of high harmonic generation efficiency,” Phys. Rev. Lett. 103 (2009).
  * [38] B. L. Henke, E. M. Gullikson, and J. C. Davis, “X-ray interactions: Photoabsorption, scattering, transmission, and reflection at E = 50-30, 000 eV, Z = 1-92,” At. Data Nucl. Data Tables 54, 181–342 (1993).
  * [39] J. Tross, G. Kolliopoulos, and C. A. Trallero-Herrero, “A self-referencing attosecond interferometer,” APS 2017, Q1–139 (2017).
  * [40] J. Troß and C. A. Trallero-Herrero, “High harmonic generation spectroscopy via orbital angular momentum,” J. Chem. Phys. 151, 84308 (2019).
  * [41] M. Guizar-Sicairos and J. C. Gutiérrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A 21, 53 (2004).


Alternative proofs of Shafer’s inequality for inverse hyperbolic tangent

Yogesh J. Bagul1,∗, Ramkrishna M. Dhaigude2

1,∗Department of Mathematics,

K. K. M. College Manwath,

Dist: Parbhani(M.S.) - 431505, India.

Email<EMAIL_ADDRESS>
2Department of Mathematics,

Government Vidarbha Institute of Science,

and Humanities, Amravati(M. S.)-444604, India

Email<EMAIL_ADDRESS>
Abstract. We point out that a concise proof of Theorem 2 in the article, ’On a
quadratic estimate of Shafer’ by L. Zhu contains a small mistake. Correcting
this mistake and giving alternative proofs of Theorem 2 is the main aim of
this note.

††footnotetext:  ∗Corresponding author 2010  Mathematics Subject
Classification: 26A09, 26D05, 26D15. Keywords: Shafer’s inequality, inverse
hyperbolic tangent, increasing-decreasing function.

##  1 Introduction and Correction

In 2008, L. Zhu [6] published a new proof of the following theorem:

######  Theorem 1.

Let $0<x<\sqrt{15}/4.$ Then

$\displaystyle\frac{\operatorname{\tanh^{-1}}x}{x}<\frac{8}{3+\sqrt{25-\frac{80}{3}x^{2}}}.$
(1.1)

The inequality (1.1) was originally established by R. E. Shafer [3, 4, 5] and
its alternative proof is given in [6] in a concise way. Though the proof of
Theorem 1 is given in a simple way in [6], it contains a small mistake which
can be explained as follows:

While giving the proof of Theorem 1, it is shown in [6] that the function

$H(x)=\frac{25-\left(\frac{8x}{\operatorname{\tanh^{-1}}x}-3\right)^{2}}{x^{2}}$

is decreasing on $(0,\sqrt{15}/4).$ This is accomplished by showing

$I(t)=\frac{-4\sinh^{2}t+3t\sinh t\cosh
t+t^{2}\cosh^{2}t}{t^{4}\cosh^{2}t}=\frac{A(t)}{B(t)}$

to be decreasing on $(0,\operatorname{\tanh^{-1}}\sqrt{15}/4)$ due to the
transformation $H(x)=16I(t),$ where $\operatorname{\tanh^{-1}}x=t.$ A careful
observation shows that the denominator $B(t)$ of $I(t)$ is mistaken as
$t^{4}\cosh^{2}t$ instead of $t^{2}\sinh^{2}t.$ Fortunately, the function
$I(t)$ remains decreasing for either expression for $B(t)$ and the final
conclusion is unaffected. For final conclusion, the following lemma is used.

######  Lemma 1.

([2]) Let $A(x)=\sum_{n=0}^{\infty}a_{n}x^{n}$ and
$B(x)=\sum_{n=0}^{\infty}b_{n}x^{n}$ be convergent for $|x|<R,$ where $a_{n}$
and $b_{n}$ are real numbers for $n=0,1,2,\cdots$ such that $b_{n}>0.$ If the
sequence ${a_{n}/b_{n}}$ is strictly increasing(or decreasing), then the
function $A(x)/B(x)$ is also strictly increasing(or decreasing) on $(0,R).$

However, it is necessary to show that how $I(t)$ is decreasing on
$(0,\operatorname{\tanh^{-1}}\sqrt{15}/4)$ with $B(t)=t^{2}\sinh^{2}t.$ In
fact, with this $B(t)$ the proof becomes more clear and convincing. Here we
present the proof.

Corrected proof of Theorem 1. As in the concise proof of Theorem 2 in [6], we
have

$A(t)=\sum_{n=1}^{\infty}a_{n}t^{2n+2},$

where

$a_{n}=\frac{-4\cdot 2^{2n+2}+3(2n+2)\cdot 2^{2n+1}+(2n+1)(2n+2)\cdot
2^{2n}}{2(2n+2)!}.$

Now

$\displaystyle B(t)$ $\displaystyle=t^{2}\sinh^{2}t$
$\displaystyle=\frac{t^{2}}{2}\left(\cosh 2t-1\right)$
$\displaystyle=-\frac{1}{2}t^{2}+\frac{1}{2}t^{2}\sum_{n=0}^{\infty}\frac{2^{2n}}{(2n)!}t^{2n}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{2^{2n-1}}{(2n)!}t^{2n+2}=\sum_{n=1}^{\infty}b_{n}t^{2n+2}$

where $b_{n}=\frac{2^{2n-1}}{(2n)!}=\frac{(n+1)(2n+1)\cdot 2^{2n}}{(2n+2)!}.$
Then we write

$\displaystyle\frac{a_{n}}{b_{n}}$
$\displaystyle=\frac{(n+1)(2n+1)+6(n+1)-8}{(n+1)(2n+1)}$
$\displaystyle=\frac{2n^{2}+9n-1}{2n^{2}+3n+1}$
$\displaystyle=1+\frac{2(3n-1)}{2n^{2}+3n+1}:=1+2c_{n}$

where

$c_{n}=\frac{3n-1}{2n^{2}+3n+1}\,\,\text{and}\,\,c_{n+1}=\frac{3n+2}{2n^{2}+7n+6},\,\,n=1,2,3,\cdots$

We claim that

$c_{n}\geq c_{n+1},\,\,n\geq 1.$

Equivalently,

$\frac{3n-1}{2n^{2}+3n+1}\geq\frac{3n+2}{2n^{2}+7n+6}$

or

$19n^{2}+11n-6\geq 13n^{2}+9n+2.$

i.e., $6n^{2}+2n\geq 8$ which is true for $n\geq 1.$ Therefore a sequence
$\left\\{\frac{a_{n}}{b_{n}}\right\\}$ is decreasing for $n=1,2,3,\cdots.$
Hence by Lemma 1, $I(t)$ is also decreasing on
$(0,\operatorname{\tanh^{-1}}\sqrt{15}/4).$ ∎

Next, it is interesting to see other simple proofs of Theorem 1.

##  2 Alternative simple proofs

We give two alternative simple proofs of Theorem 1. The first proof, is very
elementary and uses basic calculus only.

First simple proof of Theorem 1. If we let $\operatorname{\tanh^{-1}}x=t,$
then it suffices to prove that

$\frac{t}{\tanh t}<\frac{8}{3+\sqrt{25-\frac{80}{3}\tanh^{2}t}},$

for $t\in(0,\operatorname{\tanh^{-1}}\sqrt{15}/4).$ Equivalently we want

$\left(8\frac{\sinh t}{t}-3\cosh
t\right)^{2}>25\cosh^{2}t-\frac{80}{3}\sinh^{2}t.$

i.e.

$64\sinh^{2}t-48t\sinh t\cosh
t>16t^{2}\cosh^{2}t-\frac{80}{3}t^{2}\sinh^{2}t.$

Or

$192\sinh^{2}t-144t\sinh t\cosh t>48t^{2}\cosh^{2}t-80t^{2}\sinh^{2}t.$

i.e.

$12\sinh^{2}t-9t\sinh t\cosh t>3t^{2}-2t^{2}\sinh^{2}t.$

Now suppose,

$f(t)=12\sinh^{2}t+2t^{2}\sinh^{2}t-9t\sinh t\cosh t-3t^{2}.$

Successive differentiations with respect to $t$ give

$\displaystyle f^{\prime}(t)$ $\displaystyle=15\sinh t\cosh
t-5t\sinh^{2}t+4t^{2}\sinh t\cosh t-9t\cosh^{2}t-6t,$ $\displaystyle
f^{\prime\prime}(t)$ $\displaystyle=10\sinh^{2}t+15\cosh^{2}t-20t\sinh t\cosh
t+4t^{2}\sinh^{2}t+4t^{2}\cosh^{2}t$ $\displaystyle-9\cosh^{2}t-6,$
$\displaystyle f^{\prime\prime\prime}(t)$ $\displaystyle=12\sinh t\cosh
t-12t\sinh^{2}t-12t\cosh^{2}t+16t^{2}\sinh t\cosh t,$ $\displaystyle
f^{iv}(t)$ $\displaystyle=16t^{2}\sinh^{2}t+16t^{2}\cosh^{2}t-16t\sinh t\cosh
t$ $\displaystyle=16t^{2}\sinh^{2}t+16t\cosh t(t\cosh t-\sinh t)>0$

due to well-known inequality $\frac{\sinh t}{t}<\cosh t,\,\,t>0.$ This implies
that $f^{\prime\prime\prime}(t)$ is strictly increasing for $t>0$ and hence
$f^{\prime\prime\prime}(t)>f^{\prime\prime\prime}(0).$ Since,
$f^{\prime\prime\prime}(0)=f^{\prime\prime}(0)=f^{\prime}(0)=f(0),$ we
continue the argument and conclude that $f(t)>f(0)=0.$ This completes the
proof.

Second simple proof of Theorem 1. Since

$I(t)=\frac{t^{2}\cosh^{2}t+3t\sinh t\cosh t-4\sinh^{2}t}{t^{2}\sinh^{2}t}$

we have

$I^{\prime}(t)=-\frac{1}{4t^{3}\sinh^{3}t}i(t),$

where

$i(t)=\left(24\sinh t-8\sinh 3t-3t\cosh t+3t\cosh 3t+8t^{3}\cosh
t+12t^{2}\sinh t\right).$

Substituting the two formulas

$\cosh kt=\sum_{n=0}^{\infty}\frac{k^{2n}}{(2n)!}t^{2n}\,\,\text{and}\,\,\sinh
kt=\sum_{n=0}^{\infty}\frac{k^{2n+1}}{(2n+1)!}t^{2n+1}$

into the previous formula to get

$\displaystyle i(t)$
$\displaystyle=24\sum_{n=0}^{\infty}\frac{1}{(2n+1)!}t^{2n+1}-8\sum_{n=0}^{\infty}\frac{3^{2n+1}}{(2n+1)!}t^{2n+1}-3t\sum_{n=0}^{\infty}\frac{1}{(2n)!}t^{2n}$
$\displaystyle+3t\sum_{n=0}^{\infty}\frac{3^{2n}}{(2n)!}t^{2n}+8t^{3}\sum_{n=0}^{\infty}\frac{1}{(2n)!}t^{2n}+12t^{2}\sum_{n=0}^{\infty}\frac{1}{(2n+1)!}t^{2n+1}$
$\displaystyle=24\sum_{n=0}^{\infty}\frac{1}{(2n+1)!}t^{2n+1}-8\sum_{n=0}^{\infty}\frac{3^{2n+1}}{(2n+1)!}t^{2n+1}-3\sum_{n=0}^{\infty}\frac{1}{(2n)!}t^{2n+1}$
$\displaystyle+\sum_{n=0}^{\infty}\frac{3^{2n+1}}{(2n)!}t^{2n+1}+8\sum_{n=0}^{\infty}\frac{1}{(2n)!}t^{2n+3}+12\sum_{n=0}^{\infty}\frac{1}{(2n+1)!}t^{2n+3}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{24}{(2n+1)!}t^{2n+1}-\sum_{n=1}^{\infty}\frac{8\times
3^{2n+1}}{(2n+1)!}t^{2n+1}-\sum_{n=1}^{\infty}\frac{3}{(2n)!}t^{2n+1}$
$\displaystyle+\sum_{n=1}^{\infty}\frac{3^{2n+1}}{(2n)!}t^{2n+1}+\sum_{n=1}^{\infty}\frac{8}{(2n-2)!}t^{2n+1}+\sum_{n=1}^{\infty}\frac{12}{(2n-1)!}t^{2n+1}$
$\displaystyle=\sum_{n=4}^{\infty}\frac{d_{n}}{(2n+1)!}t^{2n+1}$

where

$\displaystyle d_{n}$ $\displaystyle=24-8\times
3^{2n+1}-3(2n+1)+3^{2n+1}(2n+1)+8(2n-1)(2n)(2n+1)$
$\displaystyle+12(2n)(2n+1)$ $\displaystyle=24-8\times
3^{2n+1}-6n-3+3^{2n+1}\cdot 2n+3^{2n+1}+16n(4n^{2}-1)$
$\displaystyle+24n(2n+1)$
$\displaystyle=3^{2n+1}\cdot(2n-7)+64n^{3}+48n^{2}+2n+21>0.$

So $i(t)>0$ holds for all $t>0$ giving us $I^{\prime}(t)<0.$ Thus $I(t)$ is
decreasing on $(0,\operatorname{\tanh^{-1}}\sqrt{15}/4)$ and so is $H(x)$ on
$(0,\sqrt{15}/4).$ Consequently,

$H(0+)>H(x)$

and with $H(0+)=80/3$ we get the inequality (1.1). ∎

Competing Interests. Authors would like to state that they do not have any
competing interest.

Author’s contributions. Both the authors contributed equally for this paper.

Acknowledgement. The authors would like to sincerely thank Prof. L. Zhu for
checking the manuscript before submission and suggesting second simple proof.

## References

  * [1] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products, Elsevier, edn. 2007.
  * [2] V. Heikkala, M. K. Vamanamurthy and M. Vuorinen, _Generalized elliptic integrals,_ Comput. Methods Funct. Theory, Vol. 9, No. 1, pp. 75-109, 2009.
  * [3] R. E. Shafer, On quadratic approximation, SIAM J. Numerical Analysis, Vol. 11, No. 2, pp. 447-460, 1974.
  * [4] R. E. Shafer, Analytic inequalities obtained by quadratic approximation, Univ. Beograd. Publ. Elektrotehn. Fak. Ser. Mat. Fiz., No. 577-No. 598(1977), pp. 96-97, 1977.
  * [5] R. E. Shafer, On quadratic approximation II, Univ. Beograd. Publ. Elektrotehn. Fak. Ser. Mat. Fiz., No. 602-No. 603(1978), pp. 163-170, 1978.
  * [6] L. Zhu, On a quadratic estimate of Shafer, J. Math. Ineqal., Vol. 2, No. 4, pp. 571-574, 2008. Doi: 10.7153/jmi-02-51


\[ 2  \left[ \max_k P_0\left(  \phi_{n,k,t,a_0} - \phi_{\infty, t,a_0}\right)^2 \right]^{1/2}.\]

By Lemma <ref>, <ref> implies that this is bounded by $
2C \left(\max_k C_{1k}\right)^{1/2}+2C\left(\max_k C_{2k}\right)^{1/2}+2C\left(\max_k C_{3k}\right)^{1/2}$
for $C$ depending only on $\eta$, and $C_{1k}$, $C_{2k}$ and $C_{3k}$ defined as in the proof of Lemma 6. By <ref>, this upper bound is $\fasterthan(1)$, which implies that $\left|\theta_{n}(t,a_0) -\theta_0(t,a_0) \right| = \fasterthan(1)$.

For uniform consistency, Lemma <ref> and the triangle inequality, <ref> implies that

\begin{align*}
\sup_{u\in[0,t]} \left|\theta_{n}(u,a_0) -\theta_0(u,a_0) \right| &\leq \sup_{u\in[0,t]}  \left| \d{P}_n \phi_{\infty,u,a_0}^* \right|+ \sup_{u\in[0,t]}\left|\frac{1}{K} \sum_{k=1}^K \frac{K n_k^{1/2}}{n} \d{G}_{n}^k \left(\phi_{n,k,u,a_0} - \phi_{\infty,u,a_0}\right)\right| \\
&\qquad+ \sup_{u\in[0,t]}\left|\frac{1}{K} \sum_{k=1}^K \frac{K n_k}{n} P_0 \left(\phi_{n,k,u,a_0} - \phi_{\infty, u, a_0}\right) \right|.
\end{align*}

Lemma <ref> implies that $\{\phi_{\infty, u, a_0} : u \in [0,t]\}$ is a $P_0$-Donsker class, so the first term on the right-hand side of the inequality above is $\bounded(n^{-1/2})$. By Lemma <ref>, <ref> and <ref> imply that the second term  is $\fasterthan(n^{-1/2})$. As above, the third term is $\fasterthan(1)$ by Lemma <ref>. We thus find that $\sup_{u \in [0,t]}\left|\theta_{n}(u,a_0) -\theta_0(u,a_0) \right| = \fasterthan(1)$.

Since <ref> holds automatically when $S_\infty = S_0$, $G_\infty = G_0$ and $\pi_\infty = \pi_0$, Lemma <ref> implies that

\begin{align}
\theta_{n}(t,a_0) -\theta_0(t,a_0) &= \d{P}_n \phi_{0,t,a_0}^* + \frac{1}{K} \sum_{k=1}^K \frac{K n_k^{1/2}}{n} \d{G}_{n}^k \left(\phi_{n,k,t,a_0} - \phi_{0,t,a_0}\right) + \frac{1}{K} \sum_{k=1}^K \frac{K n_k}{n} \left[P_0 \phi_{n,k,t,a_0} - \theta_0(t,a_0)\right]\label{eq:theta_decomp_asy_lin}.
\end{align}

Since <ref> and <ref> hold by assumption, the second summand on the right-hand side is $\fasterthan(n^{-1/2})$ by Lemma <ref>, where we replace the symbol $\infty$ by $0$ throughout.  By Lemma <ref>, $P_0 \phi_{n,k,t,a_0} - \theta_0(t,a_0)$ equals

\[E_{0} \left[ S_{n,k}(t \,|\, a_0, W) \int_0^t \frac{S_0(u\m \,|\, a_0, W)}{S_{n,k}(u \,|\, a_0, W)}\left\{\frac{\pi_0(a_0 \,|\, W) G_0(u \,|\, a_0, W)}{\pi_{n,k}(a_0 \,|\, W)G_{n,k}(u \,|\, a_0, W)} - 1\right\} (\Lambda_{n,k} - \Lambda_0)(du \,|\, a_0, W)   \right].\]

By the Duhamel equation (Theorem 6 of 14), we have that

\[\frac{S_0(u\m \,|\, a_0, W)}{S_{n,k}(u \,|\, a_0, W)}(\Lambda_{n,k} - \Lambda_0)(du \,|\, a_0, W) =\left( \frac{S_0}{S_{n,k}} - 1 \right)(du \,|\, a_0, W),\]

and so the above equals

\begin{align*}
&E_{0} \left[ S_{n,k}(t \,|\, a_0, W) \int_0^t \left\{\frac{\pi_0(a_0 \,|\, W) G_0(u \,|\, a_0, W)}{\pi_{n,k}(a_0 \,|\, W)G_{n,k}(u \,|\, a_0, W)} - 1\right\} \left(\frac{S_0}{S_{n,k}} - 1 \right)(du \,|\, a_0, W)   \right] \\
&\qquad=E_{0} \left[ S_{n,k}(t \,|\, a_0, W) \int_0^t \left\{ \frac{\pi_0(a_0 \,|\, W) }{\pi_{n,k}(a_0 \,|\, W)} - 1\right\} \left(\frac{S_0}{S_{n,k}} - 1 \right)(du \,|\, a_0, W)   \right]\\
&\qquad\qquad + E_{0} \left[ S_{n,k}(t \,|\, a_0, W) \int_0^t \frac{\pi_0(a_0 \,|\, W) }{\pi_{n,k}(a_0 \,|\, W)} \left\{ \frac{G_0(u \,|\, a_0, W)}{G_{n,k}(u \,|\, a_0, W)}  - 1 \right\} \left(\frac{S_0}{S_{n,k}} - 1 \right)(du \,|\, a_0, W)   \right]\\
&\qquad=E_{0} \left[\left\{ \frac{\pi_0(a_0 \,|\, W) }{\pi_{n,k}(a_0 \,|\, W)} - 1\right\} S_{n,k}(t \,|\, a_0, W) \int_0^t\frac{S_0}{S_{n,k}}(du \,|\, a_0, W) \right] \\
&\qquad\qquad+E_0\left[\frac{\pi_0(a_0 \,|\, W) }{\pi_{n,k}(a_0 \,|\, W)}S_{n,k}(t \,|\, a_0, W) \int_0^t \left\{ \frac{G_0(u \,|\, a_0, W)}{G_{n,k}(u \,|\, a_0, W)}  - 1 \right\} \left( \frac{S_0}{S_{n,k}} - 1\right)(du \,|\, a_0, W)   \right] \\
&\qquad=E_{0} \left[ \frac{\left\{ \pi_{n,k}(a_0 \,|\, W) -  \pi_0(a_0 \,|\, W) \right\} \left\{ S_{n,k}(t \,|\, a_0, W) - S_0(t \,|\, a_0, W)\right\}}{\pi_{n,k}(a_0 \,|\, W)}\right] \\
&\qquad\qquad+E_{0} \left[ \frac{\pi_0(a_0 \,|\, W) }{\pi_{n,k}(a_0 \,|\, W)}S_{n,k}(t \,|\, a_0, W) \int_0^t \left\{ \frac{G_0(u \,|\, a_0, W)}{G_{n,k}(u \,|\, a_0, W)}  - 1 \right\} \left( \frac{S_0}{S_{n,k}} - 1\right)(du \,|\, a_0, W)   \right].
\end{align*}

In view of <ref>, the absolute value of the right-hand side of the last equality is bounded above by

\begin{align*}
&\eta E_{0} \left[ \left| \pi_{n,k}(a_0 \,|\, W) -  \pi_0(a_0 \,|\, W) \right| \left| S_{n,k}(t \,|\, a_0, W) - S_0(t \,|\, a_0, W)\right|\right] \\
&\qquad+\eta E_0 \left[ \left| S_{n,k}(t \,|\, a_0, W) \int_0^t \left[ \frac{G_0(u\,|\, a_0, W)}{G_{n,k}(u \,|\, a_0, W)}  - 1 \right]  \left(\frac{S_0}{S_{n,k}} -1 \right)(du \,|\, a_0, W)  \right| \right].
\end{align*}

The maximum over $k$ of these expressions equals $\eta (r_{n,t,a,1} + r_{n,t,a,2})$. Therefore, since $\frac{1}{K} \sum_{k=1}^K \frac{K n_k}{n} \leq 2$ (as discussed in the proof of Theorem <ref> above), we have

\[ \left|\frac{1}{K} \sum_{k=1}^K \frac{K n_k}{n} \left[P_0 \phi_{n,k,t,a_0} - \theta_0(t,a_0)\right]  \right| \leq 2\eta(r_{n,t,a_0,1} + r_{n,t,a_0,2})\ ,\]
is $\fasterthan(n^{-1/2})$ by <ref>. This establishes that $\theta_n(t,a_0) = \theta_0(t,a_0) + \d{P}_n \phi_{0,t,a_0}^* + \fasterthan(n^{-1/2})$, as claimed in Theorem <ref>. Since $\phi_{0,t,a_0}^*$ is uniformly bounded, $P_0\phi^{*2}_{0,t,a_0}<\infty$, and since $P_0\phi^*_{0,t,a_0}=0$, it follows, as claimed, that

\[n^{1/2} \d{P}_n \phi_{0,t,a_0}^* \indist N\left( 0, P_0\phi_{0,t,a_0}^{*2}\right).\]

For the uniform statements, we use the same decomposition as in (<ref>). Since <ref>, <ref> and <ref> hold by assumption,

\[\sup_{u \in [0,t]} \left|  \frac{1}{K} \sum_{k=1}^K \frac{K n_k^{1/2}}{n} \d{G}_{n}^k \left(\phi_{n,k,u,a_0} - \phi_{0,u,a_0}\right)\right| = \fasterthan(n^{-1/2})\]

by Lemma <ref>, where we replace the symbol $\infty$ by $0$ throughout. Following the derivation above, we have that

\[ \sup_{u \in [0,t]} \left| \frac{1}{K} \sum_{k=1}^K \frac{K n_k}{n} \left[P_0 \phi_{n,k,u,a_0} - \theta_0(u,a_0)\right]\right| \leq \sup_{u \in [0,t]} 2\eta(r_{n,u,a_0,1} + r_{n,u,a_0,2}),\]

which is $\fasterthan(n^{-1/2})$ by <ref>. Therefore, $ \sup_{u \in [0,t]}\left| \theta_n(u,a_0) - \theta_0(u,a_0) - \d{P}_n \phi_{0,u,a_0}^*\right| = \fasterthan(n^{-1/2})$
as claimed. Since $\{\phi_{0,u,a_0}^* : u \in [0,t]\}$ is a uniformly bounded $P_0$-Donsker class by Lemma <ref>, $\left\{n^{1/2}  \d{P}_n \phi_{0,u,a_0}^* : u \in [0, t]\right\}$ converges weakly to a tight mean-zero Gaussian process with covariance $(u,v) \mapsto P_0\left( \phi_{0,u,a_0}^*\phi_{0,v,a_0}^*\right)$.

We first note that for $P_0$-almost every $(a,w)$ and all $t \leq \tau$, conditions (A1)–(A5) imply that

\begin{align*}
E_0\left[ \frac{\Delta I(Y \leq t)}{G_0(Y \,|\, a, w)}\, \middle|\, A = a, W =w\right] &= -\int_0^\infty \frac{I(y \leq t) }{G_0(y \,|\, a, w)}G_0(y  \,|\, a, w) S_0(dy \,|\, a, w) = -\int_0^{t}  S_0(dy \,|\, a, w) \\
&= 1 - S_0(t  \,|\, a, w) = P_{0,F}(T \leq t \, | \, A = a, W=w)\ .
\end{align*}

Therefore, for any such $(t, a,w)$

\[ s \mapsto  E_0 \left[ s \left[s - 2 \left\{ 1 - \frac{\Delta I(Y \leq t)}{G_0(Y \,|\, A, W)}\right\}\right]\,\middle|\, A = a, W =w\right] =s \left\{s - 2 S_0(t \, | \, a,w)\right\} \]

is uniquely minimized by $s = S_0(t \, | \, a,w)$. Hence, any minimizer $S^*$ of

\[ S \mapsto E_0\left[ \int_0^\tau E_0 \left[ S(t \, | \, A,W) \left[S(t \, | \, A,W) - 2 \left\{ 1 - \frac{\Delta I(Y \leq t)}{G_0(Y \,|\, A, W)}\right\}\right]\,\middle|\, A , W\right] dt\right]\]

satisfies $S^*(t \, | \, a,w) = S_0(t \, | \, a,w)$ for $P_0$-almost every $(a,w)$ and all $t \leq \tau$. If the integral with respect to $t$ can be exchanged with the conditional expectation with respect to $(Y,\Delta)$ given $A=a$ and $W=w$, then the result follows. To justify exchanging these integrals, we use Fubini's theorem, which requires demonstrating that

\[ \int_0^\tau E_0\left\{ \left| S(t \, | \, a,w) \left[S(t \, | \, a,w) - 2 \left\{ 1 - \frac{\Delta I(Y \leq t)}{G_0(Y \,|\, a, w)}\right\}\right] \right| \, \middle|\,A = a, W= w\right\}dt< \infty\]

for $P_0$-almost every $(a,w)$. Since $S(t \, | \, a,w)$ and $G_0(t \, | \, a,w)$ are both contained in $[0,1]$, it holds that

\begin{align*}
&\left| S(t \, | \, a,w) - 2 \left\{ 1 - \frac{\delta I(y \leq t)}{G_0(y\,|\, a, w)}\right\}\right|  \\
&=  \delta I(y\leq t) \left| S(t \, | \, a,w) - 2\left\{1 - \frac{1}{G_0(y\,|\, a, w)}\right\} \right| + [1-\delta I(y \leq t)] | S(t \, | \, a,w) - 2| \\
&= \delta I(y\leq t) \left\{ S(t \, | \, a,w) - 2 + \frac{2}{G_0(y\,|\, a, w)}\right\} + [1-\delta I(y \leq t)] [ 2- S(t \, | \, a,w)]\leq 2\left\{ 1+  \frac{\delta I(y\leq \tau)}{G_0(y\,|\, a, w)}\right\}\ .
\end{align*}

Thus, we find that

\begin{align*}
& \int_0^\tau E_0\left[ \left| S(t \, | \, a,w) \left[S(t \, | \, a,w) - 2 \left\{ 1 - \frac{\Delta I(Y \leq t)}{G_0(Y \,|\, a, w)}\right\}\right] \right| \, \middle|\, A = a, W = w\right] dt \\
&\qquad\leq 2\tau + 2 \int_0^\tau E_0\left[\frac{\Delta I(Y\leq \tau)}{ G_0(Y\,|\, a, w)} \,\middle|\,A = a, W = w\right] dt = 2\tau  - 2\tau \int_0^\tau S_0(dy \, | \, a, w) \leq 4\tau
\end{align*}

for $P_0$-almost every $(a,w)$. Hence, the order of integration may be exchanged.

Similarly, for any $(a, w)$ in the support of $(A,W)$ and $t \leq \tau$ such that $S_0(t\m \, | \, a, w) > 0$, we note that

\begin{align*}
E_0\left[ \frac{(1-\Delta)I(Y < t) }{S_0(Y \,|\, A, W)}\,\middle|\, A = a, W =w\right] &= -\int_0^\infty \frac{I(y < t)}{S_0(y \,|\, a, w)} S_0(y  \,|\, a, w) G_0^+(dy \,|\, a, w) \\
&= -\int_{[0,t)}  G_0^+(dy \,|\, a, w) = 1 - G_0(t  \,|\, a, w)\\
&= P_{0,F}(C < t \, | \, A = a, W=w)\ ,
\end{align*}

where $G_0^+$ is the right-continuous version of $G_0$ (which is left-continuous by definition). Therefore, for any such $(t, a,w)$, it follows that

\[ g \mapsto  E_0 \left[ g \left[g - 2 \left\{ 1 - \frac{(1-\Delta) I(Y < t)}{S_0(Y \,|\, A, W)}\right\}\right]\,\middle|\, A = a, W =w\right] =g \left\{g - 2 G_0(t \, | \, a,w)\right\} \]

is uniquely minimized by $g = G_0(t \, | \, a,w)$. For $t \leq \tau$ such that $S_0(t\m \, | \, a, w)= 0$, we have that

\begin{align*}
E_0\left[ \frac{(1-\Delta)I(Y < t) }{S_0(Y \,|\, A, W)}\,\middle|\, A = a, W =w\right] &= -\int_0^\infty \frac{I(y < t)}{S_0(y \,|\, a, w)} S_0(y  \,|\, a, w)\,G_0^+(dy \,|\, a, w) \\
&= -\int_{[0,t^+(a,w))}  G_0^+(dy \,|\, a, w) = 1 - G_0(t^+(a,w)  \,|\, a, w)\ ,
\end{align*}

where $t^+(a,w) := \inf\{t : S_0(t\, | \, a, w)= 0\}$. Hence, it follows that

\[ G \mapsto E_0\left[ \int_0^\tau E_0 \left\{ G(t \, | \, A,W) \left[G(t \, | \, A,W) - 2 \left\{ 1 - \frac{(1-\Delta) I(Y < t)}{S_0(Y \,|\, A, W)}\right\}\right]\,\middle|\, A , W\right\} \,dt \right] \]

is minimized by $G^*$, where $G^*(t \, | \, a,w) = G_0(t \, | \, a,w)$ if $S_0(t\m \, | \, a, w) > 0$  and $G^*(t \, | \, a,w) = G_0(t^+(a,w) \, | \, a,w)$ otherwise for $P_0$-almost every $(a,w)$. Similarly as was done above, it is straightforward to show that the inner function is integrable, thus justifying the use of Fubini's theorem for exchanging the order of integration. This completes the proof.


# Molecular Simulations of Liquid Jet Explosions and Shock Waves Induced by
X-Ray Free-Electron Lasers

Leonie Chatzimagas  Jochen S. Hub<EMAIL_ADDRESS>Theoretical
Physics and Center for Biophysics, Saarland University, Saarbrücken 66123,
Germany

###### Abstract

X-ray free-electron lasers (XFELs) produce X-ray pulses with high brilliance
and short pulse duration. These properties enable structural investigations of
biomolecular nanocrystals, and they allow resolving the dynamics of
biomolecules down to the femtosecond timescale. Liquid jets are widely used to
deliver samples into the XFEL beam. The impact of the X-ray pulse leads to
vaporization and explosion of the liquid jet, while the expanding gas triggers
the formation of shock wave trains traveling along the jet, which may affect
biomolecular samples before they have been probed. Here, we used molecular
dynamics simulations to reveal the structural dynamics of shock waves after an
X-ray impact. Analysis of the density in the jet revealed shock waves that
form close to the explosion center, travel along the jet with supersonic
velocities and decay exponentially with an attenuation length proportional to
the jet diameter. A trailing shock wave formed after the first shock wave,
similar to the shock wave trains in experiments. Although using purely
classical models in the simulations, the resulting explosion geometry and
shock wave dynamics closely resemble experimental findings, and they highlight
the importance of atomistic details for modeling shock wave attenuation.

## Introduction.

X-ray free-electron lasers (XFELs) are a source of X-ray radiation that enable
novel experiments in the field of structural biology. The high peak
brilliance, short pulse duration and high repetition rates enable resolving
the structure of biomolecular nanocrystals, also in a time-resolved manner [1,
2, 3, 4, 5, 1, 6, 7]. Liquid jets have been used to deliver the samples
rapidly into the beam [8, 9, 10]. Upon impact of the X-ray beam, not only the
sample is destroyed, but also the segment of the liquid jet exposed to the
X-ray pulse is vaporized. The expanding vapor drives the explosion of the
liquid jet.

In recent years, several studies analyzed XFEL-induced explosions of liquid
droplets or liquid jets [11, 12, 13, 14, 15, 16] as well as the relevance of
explosions for the design of crystallographic studies. Stan et al. used time-
resolved imaging to study explosions in droplets and jets, revealing that the
expanding vapor launches shock waves traveling across the drops or along the
jet [12, 11]. In jets, the shock front may split, leading to a sequence of
succeeding pressure and density oscillations. These findings have implications
for the design of experiments performed at high X-ray repetition rate, for two
reasons [9, 13]: First, the gap in the jet formed by the explosion should be
replenished before the arrival of the next X-ray pulse, requiring either
increased jet velocities or reduced X-ray pulse repetition rates. Second,
shock waves traveling backwards along the jet may pass the samples before they
have been probed. The pressure and density oscillations may perturb the
samples leading to lower crystallographic resolution.

The influence of shock waves on biomolecules has been studied both
computationally and experimentally. Molecular dynamics (MD) simulations
modeling laser-induced shock waves revealed that hemoglobin is compressed
during the passage of a shock wave, but the tetrameric structure remained
intact [17, 18]. Experimentally, comparing crystallographic data of hen egg-
white lysozyme (HEWL) microcrystals collected from two succesive pulses with
$1.1\text{\,}\mathrm{MHz}\text{/}$ repetition rate revealed no perturbation of
the HEWL microcrystals [9, 19, 20]. However, since the jet diameter was
smaller or approximately equal to the focal spot of the X-ray pulse, no shock
waves might have been created according to descriptions by Blaj et al. [13].
In contrast, two recent studies found a degradation of diffraction data
quality of HEWL and hemoglobin microcrystals as well as structural changes in
the hemoglobin microcrystals when probing with an effective pulse repetition
rate of $4.5\text{\,}\mathrm{MHz}\text{/}$ [10, 21]. Additional experimental
and computational studies are required to clarify how the experimental design
controls the formation and propagation of shock waves, and wether the shock
waves affect the biological samples. Here, we used MD simulations to model jet
explosions as well as the formation, dynamics, and attenuation of shock waves.

## Methods.

We developed a purely classical MD model of the liquid jet and its exposure to
an XFEL pulse. We modeled the liquid jet as a water cylinder in vacuum using
the all-atom SPC/E water model [22], if not stated otherwise. To test the
effect of the water model and the importance of atomic details, additional
simulations were carried out using the all-atom TIP4P/2005 water model [23] or
the coarse-grained MARTINI water model [24]. We modeled the explosion of the
liquid jet via instantaneous heating of the central jet segment with a
Gaussian-shaped temperature distribution by assigning new velocities to the
water molecules drawn from the Maxwell-Boltzmann distribution. Since we
focussed on the shock wave dynamics, inside the liquid water, and not on the
plasma dynamics, the photoelectric effect and Auger electrons were neglected.
We set up simulation systems with jet diameters between
$10\text{\,}\mathrm{n}\mathrm{m}$ and $80\text{\,}\mathrm{n}\mathrm{m}$ hit by
X-ray beams with intensity profiles with FWHM between
$1.5\text{\,}\mathrm{n}\mathrm{m}$ and $12\text{\,}\mathrm{n}\mathrm{m}$
(Table LABEL:S-tab:param), where the jet diameter was always larger than the
FWHM. For each system, we carried out 19 to 100 independent simulations in the
microcanonical ensemble (Table LABEL:S-tab:param). Details on the model and
simulation parameters are presented in the Supplementary Material.

## Liquid jet explosion and shock waves.

Figure 1: Simulation of jet explosion with $20\text{\,}\mathrm{nm}\text{/}$
diameter induced by a modeled X-ray pulse with $3\text{\,}\mathrm{nm}\text{/}$
FWHM. (a) Snapshots of a single simulation. (b) Radial density at
$0.1\text{\,}\mathrm{ps}\text{/}$ to $15.1\text{\,}\mathrm{ps}\text{/}$ after
the X-ray impact averaged over 50 independent simulations.

In the simulations, the impact of the X-ray pulse triggered an explosion of
the central jet segment leading to the formation of a gap and, thereby,
splitting the jet into two segments (Figs. 1 (a) and LABEL:S-jetcombiR5 (a)).
As the gap growed over tens of picoseconds, thin water films formed at the
ends of the segments, which later adopted a conical shape, and finally, after
several tens of picoseconds, folded back to the jet. The simulated explosions
developed on a by far smaller time scale compared to the experimental jet
explosion observed by Stan et al. [11] that developed in the microsecond
range. The time scale difference is likely a result of the larger experimental
jet diameter of $20\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$ diameter,
which is approximately three orders of magnitude larger compared to our MD
model. Despite these different length and time scales, the jet explosion
dynamics in simulations resemble dynamics observed in the experiments.

To analyze shock waves traveling along the jet, we computed the mass density
as a function of the axial and radial direction (Figs. 1 (b) and
LABEL:S-jetcombiR5 (b), Supplementary Material Methods). The data revealed the
formation of two density peaks positioned symmetrically around the jet center
and close to the explosion site. The peaks traveled along the jet and decayed
to density values slightly below the density of an equilibrated water jet,
revealing the rarefaction caused by the shock wave. Notably, the density peak
was maximized near the jet axis (Fig. 1 (b), $r=0$) but small or event absent
at the jet surface ($r\approx 9$ nm), suggesting that the shock wave energy is
dissipated at the jet surface.

Figure 2: (a) Time evolution of the density in a water jet with
$20\text{\,}\mathrm{nm}\text{/}$ diameter after explosion induced by a modeled
X-ray pulse with $3\text{\,}\mathrm{n}\mathrm{m}$ FWHM, averaged over 50
simulations. Lines represent the density at time delays of 0.1 to
$20.1\text{\,}\mathrm{ps}\text{/}$ after X-ray impact (see color bar). (b)
Density versus time at fixed distances of $26\text{\,}\mathrm{nm}\text{/}$
(blue) or $36\text{\,}\mathrm{nm}\text{/}$ (orange) from the jet center.

Averaging the two-dimensional densities along the radial direction yields the
time-dependent density as function of the axial direction $z$. This analysis
revealed a second density peak succeeding the first density peak (Fig. 2 (a)),
as evident from the density evolution at fixed distances from the jet center
(Fig. 2 (b)). All the simulations with different combinations of jet diameter
and FWHM of the X-rays pulse showed an explosion of the jet and a similar
evolution of the densities (Fig LABEL:S-densall). Hence, despite (i) the
approximations underlying the classical models and (ii) the smaller time and
length scales, the MD simulations reproduce the sequential shock waves
observed in previous experiments [11, 13].

## Attenuation of the shock front.

Figure 3: Shock wave attenuation and propagation. (a) Density peak attenuation
of the shock front scaled by the water jet diameter. The relative height of
the density peaks, $\Delta$ Density, attenuates with the propagation distance
to the jet center. The linear attenuation on the semi-log plot indicates an
exponential density decay. The attenuation approximately scales with the jet
diameter. (b) Propagation of the shock front. The shock propagates with
constant, supersonic velocity after few picosecond. This velocity is
approximately the same for different simulation sets and in average
$\sim$2350\text{\,}\mathrm{m}\text{/}\mathrm{s}$$ (see Table
LABEL:S-tab:param).

To analyze the attenuation of the shock front, we determined the height of the
first density peak relative to the density of the equilibrated jet. Figure 3
(a) presents the height of the peaks for different jet diameters and X-ray
pulse widths, plotted versus the propagation distance from the jet center. The
attenuation of the density peaks followed an exponential decay, as evident
from the linear decay on the semi-logarithmic scale. To test the influence of
the jet diameter, we scaled the propagated distance of the density peaks by
the jet diameter. A linear least-squares fit on the data revealed similar
slopes for the peak decay in jets with different diameter, demonstrating that
the decay length is proportional to the jet diameter or, equivalently, to the
jet circumference (see Table LABEL:S-tab:param). These findings suggest that
the shock wave energy is dissipated due to surface effects, compatible with
the reduced shock wave densities at the surface (Figs. 1 (b)).

These results are further compatible with experimental observations
considering (i) that, in experiments, the jets are three orders of magnitude
larger and the shock front travels further into the jet compared to the
simulations, reflecting that decay lengths increase with jet diameter; (ii) as
pressure and density are positively correlated, the density evolution of the
leading peak agrees with the pressure peak decay calculated from experiments
by Blaj et al. [13], who suggested that the pressure peaks decay exponentially
and that the decay length is approximately proportional to the jet diameter.

## Shock wave velocity.

Figure 3 (b) presents the peak propagation with time, demonstrating that,
after a deceleration in the first $1\text{\,}\mathrm{ps}\text{/}$ to
$4\text{\,}\mathrm{ps}\text{/}$, the peaks propagate at a constant velocity
after a few picoseconds. The shock wave velocities for all simulation systems
is summarized in Table LABEL:S-tab:param, revealing velocities between
$2200\text{\,}\mathrm{m}\text{/}\mathrm{s}$ and
$2600\text{\,}\mathrm{m}\text{/}\mathrm{s}$, where the velocities increase
with smaller jets or with increased deposited X-ray energy. Hence, the density
peaks travel with supersonic velocities, as expected for shock waves.

## Evolution of the gap size.

Figure 4: Evolution of the gap size. The gap was defined as the segment with a
density below $400\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$. The resulting
gap size was scaled by the jet diameter and previously proposed models were
fitted [11, 13].

To analyze the gap growth induced by the jet explosion, we defined the gap as
the axial segment with a density below
$400\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$, and we scaled the gap size
by the jet diameter (Fig. 4). Notably, the gap size evolution in simulation
resembled the evolution in experiments by Stan et al. [11]. However, in the
simulations, the rise of the gap size was delayed by few picoseconds, as
required for clearing the central segment after X-ray impact.

Two models have been put forward to describe the gap growth: a logarithmic
growth used by Stan et al. [11] and a power law proposed by Gañán-Calvo [15].
To test wether these models are compatible with our simulations, we fitted the
models to the simulated gap growth, augmented by a delay time before the
initial rise of the gap (Supplementary Material Methods). As shown in Fig. 4,
both models are in excellent agreement with the data. However, since the
simulation times-scales cover less then two orders of magnitude, the data is
insuffcient to unambigiously determine all parameters in the power law model
(see Table LABEL:S-tab:paramganan), or to decide whether one of the models is
preferable for describing the gap growth. Taken together, the gap growth in
simulations agrees with previous experiments and is compatible with the models
by Stan et al. [11] and Gañán-Calvo [15].

## Influence of water model and atomic details.

To test the influence of the selected water model, we carried out additional
simulations using the four-site TIP4P/2005 model. The results obtained using
TIP4P/2005 are in good agreement with results described above obtained with
the three-site SPC/E water model (Fig. LABEL:S-peakveldiff,
LABEL:S-gapsizediff), demonstrating that the choice of the atomistic water
model is not critical for modeling shock wave propagation and gap growth.

To test the importance of atomic details for modeling jet explosion and shock
wave propagation, we repeated the simulations with the coarse-grained MARTINI
water model [24]. MARTINI describes liquid water as a Lennard-Jones fluid,
modeling four molecules by one Lennard-Jones bead. Hence, the MARTINI water
model is greatly simplified as it lacks atomic details such as hydrogen bonds,
leading to a smoothed potential energy landscape. Despite these
simplifications, we observed qualitatively similar explosion and shock wave
dynamics compared to simulations using the SPC/E water model (Fig. 5, Fig.
LABEL:S-peakvelmartini and LABEL:S-gapgrowthmartini): the density peaks
attenuate exponentially and propagate with a constant velocity after a few
picoseconds; Further, the gap size evolution shows similar trends compared to
experimental estimates [11]. However, the simulations reveal a quantitative
difference of the shock wave attenuation dynamics compared to SPC/E
simulations: While the exponential decay length is still approximately
proportional to the jet diameter, the decay length scaled by the jet diameter
is larger with MARTINI as compared to SPC/E simulations (see Fig 5 and Table
LABEL:S-tab:param). We hypothesize that the slower decay is a consequence of
the smoothed potential energy landscape of the MARTINI model, which may lead
to reduced internal friction and, thereby, to slower dissipation of the shock
wave energy. Hence, atomistic details are relevant for quantitatively
describing the attenuation dynamics of the density peaks.

Figure 5: Influence of atomistic details on the density peak attenuation. The
density decays over shorter distances when using the atomistic SPC/E model
(blue) as compared to using the coarse grained MARTINI water model (orange).

## Conclusion.

We carried out a large set of MD simulations of water jet explosions after
impact of an XFEL laser pulse, modeled by a temperature jump of water
molecules at the jet center. Despite the approximations underlying our
classical simulations, we found good agreement with previous experiments in
terms of explosion geometry, shock wave dynamics, and gap growth. These
results suggests that photoelectric and Auger effects, which are certainly
critical for modeling the explosion at the core of X-ray impact, are less
important for shock wave formation and propagation or for the qualitative
dynamics of the gap growth. In the simulations, jet explosion triggered the
formation of leading and trailing shock fronts that traveled with supersonic
velocities along the jets. The shock waves attenuated exponentially with a
decay length proportional to the jet diameter, suggesting a role of surface
effects in dissipation of shock wave energy. Modeling the jet without atomic
details led to slower decay of the shock wave, likely due to an overly
smoothed potential energy landscape of water–water interactions. We expect
these insights to be useful for designing experiments at XFELs with high
repetition rate.

###### Acknowledgements.

This project was supported by the Deutsche Forschungsgemeinschaft (HU
1971/3-1).

## References

  * Chapman [2019] H. N. Chapman, Annu. Rev. Biochem. (2019).
  * Pande  _et al._ [2016] K. Pande, C. D. Hutchison, G. Groenhof, A. Aquila, J. S. Robinson, J. Tenboer, S. Basu, S. Boutet, D. P. DePonte, M. Liang, _et al._ , Science 352, 725 (2016).
  * Nogly  _et al._ [2018] P. Nogly, T. Weinert, D. James, S. Carbajo, D. Ozerov, A. Furrer, D. Gashi, V. Borin, P. Skopintsev, K. Jaeger, _et al._ , Science 361, eaat0094 (2018).
  * Pandey  _et al._ [2020] S. Pandey, R. Bean, T. Sato, I. Poudyal, J. Bielecki, J. C. Villarreal, O. Yefanov, V. Mariani, T. A. White, C. Kupitz, _et al._ , Nat. Methods 17, 73 (2020).
  * Neutze  _et al._ [2000] R. Neutze, R. Wouts, D. van der Spoel, E. Weckert, and J. Hajdu, Nature 406, 752 (2000).
  * Roedig  _et al._ [2017] P. Roedig, H. M. Ginn, T. Pakendorf, G. Sutton, K. Harlos, T. S. Walter, J. Meyer, P. Fischer, R. Duman, I. Vartiainen, _et al._ , Nat. Methods 14, 805 (2017).
  * Gisriel  _et al._ [2019] C. Gisriel, J. Coe, R. Letrun, O. M. Yefanov, C. Luna-Chavez, N. E. Stander, S. Lisova, V. Mariani, M. Kuhn, S. Aplin, _et al._ , Nat. Commun. 10, 1 (2019).
  * Chapman  _et al._ [2011] H. N. Chapman, P. Fromme, A. Barty, T. A. White, R. A. Kirian, A. Aquila, M. S. Hunter, J. Schulz, D. P. DePonte, U. Weierstall, _et al._ , Nature 470, 73 (2011).
  * Grünbein  _et al._ [2018] M. L. Grünbein, J. Bielecki, A. Gorel, M. Stricker, R. Bean, M. Cammarata, K. Dörner, L. Fröhlich, E. Hartmann, S. Hauf, _et al._ , Nat. Commun. 9, 3487 (2018).
  * Grünbein  _et al._ [2021] M. L. Grünbein, L. Foucar, A. Gorel, M. Hilpert, M. Kloos, K. Nass, G. N. Kovacs, C. M. Roome, R. L. Shoeman, M. Stricker, S. Carbajo, W. Colocho, S. Gilevich, M. Hunter, J. Lewandowski, A. Lutman, J. E. Koglin, T. J. Lane, T. van Driel, J. Sheppard, S. L. Vetter, J. L. Turner, R. B. Doak, T. R. M. Barends, S. Boutet, A. L. Aquila, F.-J. Decker, I. Schlichting, and C. A. Stan, Phys. Rev. Research 3, 013046 (2021).
  * Stan  _et al._ [2016a] C. A. Stan, D. Milathianaki, H. Laksmono, R. G. Sierra, T. A. McQueen, M. Messerschmidt, G. J. Williams, J. E. Koglin, T. J. Lane, M. J. Hayes, _et al._ , Nat. Phys. 12, 966 (2016a).
  * Stan  _et al._ [2016b] C. A. Stan, P. R. Willmott, H. A. Stone, J. E. Koglin, M. Liang, A. L. Aquila, J. S. Robinson, K. L. Gumerlock, G. Blaj, R. G. Sierra, _et al._ , J. Phys. Chem. Lett. 7, 2055 (2016b).
  * Blaj  _et al._ [2019] G. Blaj, M. Liang, A. L. Aquila, P. R. Willmott, J. E. Koglin, R. G. Sierra, J. S. Robinson, S. Boutet, and C. A. Stan, Phys. Rev. Fluids 4, 043401 (2019).
  * Ursescu  _et al._ [2020] D. Ursescu, V. Aleksandrov, D. Matei, I. Dancus, M. D. de Almeida, and C. A. Stan, Phys. Rev. Fluids 5, 123402 (2020).
  * Gañán-Calvo [2019] A. M. Gañán-Calvo, Phys. Rev. Lett. 123, 064501 (2019).
  * Beyerlein  _et al._ [2018] K. R. Beyerlein, H. O. Jönsson, R. Alonso-Mori, A. Aquila, S. Bajt, A. Barty, R. Bean, J. E. Koglin, M. Messerschmidt, D. Ragazzon, _et al._ , Proc. Natl. Acad. Sci. U. S. A. 115, 5652 (2018).
  * Wiederschein [2009] F. Wiederschein,  _Investigation of Laser-Induced-Liquid-Beam-Ion-Desorption (LILBID) with Molecular Dynamics Simulations_ , Ph.D. thesis, Georg-August-Universität Göttingen (2009).
  * Vöhringer-Martinez [2008] E. Vöhringer-Martinez,  _Dynamics, ionization and charge separation in superheated metastable water_ , Ph.D. thesis, Georg-August-Universität Göttingen (2008).
  * Wiedorn  _et al._ [2018] M. O. Wiedorn, D. Oberthür, R. Bean, R. Schubert, N. Werner, B. Abbey, M. Aepfelbacher, L. Adriano, A. Allahgholi, N. Al-Qudami, _et al._ , Nat. Commun. 9, 1 (2018).
  * Yefanov  _et al._ [2019] O. Yefanov, D. Oberthür, R. Bean, M. O. Wiedorn, J. Knoska, G. Pena, S. Awel, L. Gumprecht, M. Domaracky, I. Sarrou, _et al._ , Struct. Dyn. 6, 064702 (2019).
  * Grünbein  _et al._ [2021] M. L. Grünbein, A. Gorel, L. Foucar, S. Carbajo, W. Colocho, S. Gilevich, E. Hartmann, M. Hilpert, M. Hunter, M. Kloos, _et al._ , Nat. Commun. 12, 1 (2021).
  * Chatterjee  _et al._ [2008] S. Chatterjee, P. G. Debenedetti, F. H. Stillinger, and R. M. Lynden-Bell, J. Chem. Phys. 128, 124511 (2008).
  * Abascal and Vega [2005] J. L. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
  * Marrink  _et al._ [2007] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. De Vries, J. Phys. Chem. B 111, 7812 (2007).