# Direct characterisation of tuneable few-femtosecond dispersive-wave pulses
in the deep UV

Christian Brahms Blackett Laboratory, Imperial College London, Prince Consort
Road, London SW7 2AZ, UK School of Engineering and Physical Sciences, Heriot-
Watt University, Edinburgh, EH14 4AS, UK  Corresponding author:
<EMAIL_ADDRESS>Dane R. Austin  Blackett Laboratory, Imperial College
London, Prince Consort Road, London SW7 2AZ, UK Francesco Tani  Max Planck
Institute for the Science of Light, Staudtstr. 2, 91058 Erlangen, Germany
Allan S. Johnson  Blackett Laboratory, Imperial College London, Prince Consort
Road, London SW7 2AZ, UK Douglas Garratt  Blackett Laboratory, Imperial
College London, Prince Consort Road, London SW7 2AZ, UK John C. Travers
School of Engineering and Physical Sciences, Heriot-Watt University,
Edinburgh, EH14 4AS, UK  Max Planck Institute for the Science of Light,
Staudtstr. 2, 91058 Erlangen, Germany  John W.G. Tisch  Blackett Laboratory,
Imperial College London, Prince Consort Road, London SW7 2AZ, UK Philip St.J.
Russell  Max Planck Institute for the Science of Light, Staudtstr. 2, 91058
Erlangen, Germany  Jon P. Marangos  Blackett Laboratory, Imperial College
London, Prince Consort Road, London SW7 2AZ, UK

Dispersive wave emission (DWE) in gas-filled hollow-core dielectric waveguides
is a promising source of tuneable coherent and broadband radiation, but so far
the generation of few-femtosecond pulses using this technique has not been
demonstrated. Using in-vacuum frequency-resolved optical gating, we directly
characterise tuneable $3\text{\,}\mathrm{fs}$ pulses in the deep ultraviolet
generated via DWE. Through numerical simulations, we identify that the use of
a pressure gradient in the waveguide is critical for the generation of short
pulses.

Sources of very short laser pulses in the deep ultraviolet (DUV,
$200\text{\,}\mathrm{nm}300\text{\,}\mathrm{nm}$) are a key enabling
technology for many areas of physics, such as ultrafast spectroscopy [1, 2].
New approaches are made necessary by the limitations of sources based on
harmonic generation, most importantly the lack of direct spectral tuneability,
low conversion efficiency to the DUV, and the fact that only pulses of similar
duration to the driving field can be generated [3]. One promising avenue is
the use of dispersive wave emission (DWE), a phenomenon observable during
optical soliton [4] and filamentation [5] dynamics among other effects. The
extreme spectral broadening and nonlinear phase evolution experienced by a
pulse undergoing soliton self-compression allows for phase-matching and
coherent energy transfer to a secondary pulse at a different frequency, known
as a dispersive wave. Crucially, in a gas-filled waveguide the wavelength at
which phase-matching occurs is continuously tuneable by way of the gas
pressure, enabling the generation of broadband pulses from the visible to the
vacuum ultraviolet spectral range [6, 7].

Numerical studies suggest that the dispersive wave is generated as a near
transform-limited pulse with a duration of only a few femtoseconds,
significantly shorter than the driving pulse [8]. Measurements of pulses which
were strongly chirped by propagation through windows and ambient air before
characterisation were consistent with few-fs pulse durations at generation,
providing some support for this idea [9, 10]. However, to date, no experiment
has successfully measured a few-fs dispersive wave in the DUV or provided a
route towards delivering such pulses to experiments without distortions.

Here we present the full characterisation of $3\text{\,}\mathrm{fs}$
dispersive-wave pulses with central wavelength tuneable from
$225\text{\,}\mathrm{nm}$ to $300\text{\,}\mathrm{nm}$ using cross-correlation
frequency-resolved optical gating (XFROG). Our measurements represent the
first characterisation of a few-fs pulse generated via DWE and the shortest
tuneable laser pulses in this spectral region to date, demonstrating the power
of DWE-based sources for ultrafast science. The pulses are delivered to the
characterisation set-up in vacuum; the measurements therefore faithfully
reflect the pulses as generated in the waveguide. By numerically simulating
the generation process, we further establish that the use of a pressure
gradient along the waveguide is critical to the generation of few-fs pulses.

Figure 1: Layout of the experimental apparatus for (a) pre-compression, (b)
DUV generation and (c) DUV pulse characterisation. BPF: bandpass filter; BS:
beam-splitter; CM: chirped mirror; HWP: half-wave plate; TFP: thin-film
polariser; $\tau$: delay stage, $\lambda$: spectrometer, A: aperture. (d)
Pressure distribution along the $7.5\text{\,}\mathrm{cm}$ long AR-PCF for the
input fill pressures used in the experiment. (e) Phase-matching efficiency for
DFG between a DUV test pulse and a narrowband gate pulse centred at
$790\text{\,}\mathrm{nm}$ for a $5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ BBO
crystal for three different angles between the crystal axis and the beam
propagation direction, taking into account the change in effective thickness
as the crystal, cut for an angle of $None$, is rotated away from normal
incidence.

The experimental layout is shown in Fig. 1(a)-(c) and consists of three parts:
pre-compression, DUV generation, and pulse characterisation. A titanium-doped
sapphire laser amplifier delivers pulses of $30\text{\,}\mathrm{fs}$ duration
and $150\text{\,}\mathrm{\SIUnitSymbolMicro J}$ energy at a repetition rate of
$1\text{\,}\mathrm{kHz}$. The pulses are spectrally broadened in a
$1\text{\,}\mathrm{m}$ long hollow glass capillary with a
$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$ core diameter filled with a
positive pressure gradient from vacuum to $250\text{\,}\mathrm{mbar}$ of argon
and pass through a chirped-mirror compressor and an attenuator. To generate
DUV pulses, the pre-compressed driving pulses are then coupled into a
$7.5\text{\,}\mathrm{cm}$ long kagomé-type anti-resonant photonic crystal
fibre (AR-PCF) with a core diameter of $33\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ by focusing the beam to a $22\text{\,}\mathrm{\SIUnitSymbolMicro m}$
diameter ($1/\mathrm{e}^{2}$) spot. A cross-section micrograph of the PCF is
shown in the inset in Fig. 1(b). A spatial filter consisting of a
$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter pinhole placed in the
focal plane of a telescope removes distortions from the beam profile to
optimise the coupling. The inset in Fig. 1(a) shows the pulse profile at the
entrance of the AR-PCF as measured using a second-harmonic generation (SHG)
FROG device. Pre-compression of the driving pulse is not required but leads to
higher-quality soliton self-compression and less structure in the DWE spectrum
[8].

The AR-PCF is sealed into the connection between a gas cell and a vacuum
chamber such that its entrance can be pressurised while keeping the exit under
vacuum. This is critical to avoid stretching of the DUV pulse in an exit gas
cell or a window. The shape of the resulting pressure gradient is given by
$P(z)=P_{0}\sqrt{(L-z)/L}$, where $z$ is the distance along the waveguide, $L$
is the length of the waveguide and $P_{0}$ is the fill pressure at the
entrance [7]. Fig. 1(d) shows $P(z)$ in the AR-PCF used in the experiment for
three different fill pressures.

Due to the constraints of the vacuum system, DUV pulse energy measurement and
pulse characterisation cannot be carried out simultaneously. An energy
measurement with the gas cell at $7.5\text{\,}\mathrm{bar}$ pressure of argon
and the vacuum chamber at atmospheric pressure showed a DUV pulse energy of
$30\text{\,}\mathrm{nJ}$ with $1\text{\,}\mathrm{\SIUnitSymbolMicro J}$
driving pulse energy, corresponding to a conversion efficiency of
$3\text{\,}\mathrm{\char 37\relax}$. This is in line with previous results [7]
but can likely be optimised further. Higher energies can be achieved by moving
to a larger core size [11].

The pulse characterisation is based on in-vacuum XFROG between a DUV pulse
(the test pulse) and a narrowband infrared pulse (the gate) in a beta barium
borate (BBO) crystal (see Fig. 1(c)) [12]. In contrast to gas-based
autocorrelation measurements [3], XFROG allows for the full characterisation
of the temporal intensity and phase of the pulses. We use difference-frequency
generation (DFG) as the nonlinear interaction since BBO is opaque below
$190\text{\,}\mathrm{nm}$, precluding the use of sum-frequency generation
(SFG). Compared to third-order processes, DFG has the advantage of higher
sensitivity as well as generating a signal at a different wavelength to the
test pulse; this removes a potential source of background from the measurement
[12]. As shown in Fig. 1(e), using a type I (o-o-e) scheme in a
$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick BBO crystal, broadband phase-
matching for DFG across the DUV can be achieved by tuning the crystal angle.

A portion of the pulse is split off after the hollow capillary and bandpass-
filtered to create the gate pulse; its spectrum is shown in the inset in Fig.
1(a). It enters the vacuum chamber through a separate window after traversing
a motorised delay line. An aperture in the test pulse arm passes the DUV part
of the spectrum undisturbed while removing over $70\text{\,}\mathrm{\char
37\relax}$ of the driving pulse energy. This prevents nonlinear interactions
between the driving pulse and the DUV pulse. The gate and test pulses are
brought to a common focus in the BBO crystal, and the three beams exiting the
crystal are collimated and steered out of the chamber for analysis by a
spectrometer. The signal spot is imaged onto the spectrometer entrance slit
with sufficient magnification to reduce the temporal blurring introduced by
the non-collinear geometry to much less than the delay step in the XFROG
trace. A moveable block allows for the measurement of the test pulse spectrum
in the same position.

Pulse retrieval from XFROG measurements has traditionally relied upon the
accurate characterisation of the gate pulse [12]. Recently, the use of
ptychographic techniques has enabled the retrieval of pulses from XFROG traces
without a known gate pulse [13]. Furthermore, it was demonstrated that
ptychographic XFROG is capable of measuring pulses many times shorter than the
gate [14]. We apply this method to DFG XFROG for the first time. We use the
regularised iterative ptychographic engine (rPIE) outlined in [15] with the
addition of coefficient randomisation and soft thresholding described in [16]
to aid convergence. To transfer the algorithm from SFG to DFG it has to be
adapted, because the signal field in DFG XFROG is given by

$\psi(t)=E(t)G^{*}(t-\tau)\,,$ (1)

where $E(t)$ is the field of the test pulse and $G(t-\tau)$ is the field of
the gate pulse delayed by $\tau$, whereas in the case of SFG XFROG, the gate
pulse appears without the complex conjugate. The adapted algorithm is thus
derived by replacing the gate with its complex conjugate in the SFG algorithm.
The rPIE update function to obtain a new guess for the test pulse,
$E^{\prime}(t)$, is given by

$E^{\prime}(t)=E(t)+\frac{G(t-\tau_{j})}{(1-\alpha)\absolutevalue*{G(t-\tau_{j})}^{2}+\alpha\absolutevalue*{G(t-\tau_{j})}^{2}_{\mathrm{max}}}\Delta\psi(t)\,,$
(2)

where $\alpha$ is a randomly chosen coefficient of order unity, max denotes
the maximum value, and $\Delta\psi$ is given by
$\Delta\psi=\psi^{\prime}(t)-\psi(t)$; $\psi$ is the signal field as
calculated with Eq. 1 using the current guess for $E(t)$ and $G(t)$ and
$\psi^{\prime}$ is obtained by replacing the magnitude of $\psi$ with that of
the measured XFROG trace:

$\psi^{\prime}(t)=\mathcal{F}_{\omega\rightarrow
t}^{-1}\left[\frac{\mathcal{F}_{t\rightarrow\omega}\left[\psi(t)\right]}{\absolutevalue{\mathcal{F}_{t\rightarrow\omega}\left[\psi(t)\right]}}\sqrt{S(\omega,\tau_{j})}\right]\,.$
(3)

Here $\mathcal{F}_{t\rightarrow\omega}$ and $\mathcal{F}_{\omega\rightarrow
t}^{-1}$ denote the forward and inverse Fourier transform over time,
respectively, and $S(\omega,\tau_{j})$ is the slice of the trace measured at
delay $\tau_{j}$. The gate pulse, on the other hand, is updated by

$G^{\prime}(t-\tau_{j})=G(t-\tau_{j})+\frac{E(t)}{(1-\beta)\absolutevalue{E(t)}^{2}+\beta\absolutevalue{E(t)}^{2}_{\mathrm{max}}}\Delta\psi^{*}(t)\,,$
(4)

where $\beta$ is a second coefficient, and optionally projected onto the
separately measured power spectrum of the gate pulse,
$I_{\textsc{g}}(\omega)$:

$G^{\prime\prime}(t)=\mathcal{F}_{\omega\rightarrow
t}^{-1}\left[\frac{\mathcal{F}_{t\rightarrow\omega}\left[G^{\prime}(t)\right]}{\absolutevalue{\mathcal{F}_{t\rightarrow\omega}\left[G^{\prime}(t)\right]}}\sqrt{I_{\textsc{g}}(\omega)}\right]\,.$
(5)

Figure 2: Measured (top row) and reconstructed (bottom row) DFG XFROG traces
on a logarithmic colour scale for three different argon gas pressures at the
PCF entrance. The measured traces have been corrected for the calculated
phase-matching efficiency. Figure 3: Spectral (top row) and time-domain
(bottom row) representation of the DUV pulses as retrieved from the
measurements shown in Fig. 2. In the frequency domain plots, the black line
shows the intensity, the purple line shows the spectral phase, and the red
dashed line shows the separately measured power spectrum of the test pulse.

XFROG measurements of DUV pulses generated by coupling driving pulses with
between $0.65\text{\,}\mathrm{\SIUnitSymbolMicro J}$ and
$1.15\text{\,}\mathrm{\SIUnitSymbolMicro J}$ of energy into the AR-PCF filled
with three different pressures of argon at the entrance are shown in Fig. 2.
The traces were corrected for the calculated phase-matching efficiency to
account for the fact that, while broadband, it is not uniform across the
bandwidth of the DUV pulses. Even before retrieval, it can be seen that the
pulses are close to transform-limited by the absence of any tilt in the XFROG
traces. The retrieval algorithm was run for 400 iterations with the spectral
projection of the gate pulse applied for the first 50. The excellent agreement
between the measured and retrieved traces shows that the retrieval was
successful. The retrieved pulses are shown in Fig. 3 along with the spectrum
of the pulses as measured by partially removing the block for the test and
gate pulses such that the test pulse enters the spectrometer (see Fig. 1(c)).
The good agreement between the spectra confirms that the pulses were retrieved
accurately.

Figure 4: Simulated propagation of the measured input pulse with
$1.15\text{\,}\mathrm{\SIUnitSymbolMicro J}$ energy in a
$7.5\text{\,}\mathrm{cm}$ length, $33\text{\,}\mathrm{\SIUnitSymbolMicro m}$
core diameter AR-PCF filled with argon with a negative pressure gradient from
$6\text{\,}\mathrm{bar}$ to vacuum (identical parameters to the experiment).
(a) The output DUV spectral intensity (black) and phase (purple). (b) The
spectral evolution of the pulse along the fibre on a logarithmic colour scale.
The green line shows the zero-dispersion wavelength. (c) The temporal
evolution of the field in a co-moving frame on a linear colour scale. The
inset shows the evolution of the DUV part of the field as obtained by
bandpass-filtering; the edges of the filter are shown as black dashed lines in
(b). The green line follows the peak of the DUV pulse.

All of the measured pulses consist of a main peak with a duration of around
$3\text{\,}\mathrm{fs}$. To our knowledge, these are the shortest DUV pulses
ever measured, with the exception of (non-tuneable and low-efficiency) third-
harmonic generation using a few-cycle driving pulse [3]. At
$250\text{\,}\mathrm{nm}$ and $300\text{\,}\mathrm{nm}$ central wavelength,
the pulse profile is exceptionally clean without post-pulses or significant
pedestals. This is reflected in the flat spectral phase of these pulses, which
confirms that the pulses are indeed very short at the point of generation and
that our apparatus can deliver them to an experiment without significant
dispersion. In contrast, the pulse at $225\text{\,}\mathrm{nm}$ exhibits a
pedestal after the main peak. This is caused by the more complicated spectral
intensity and phase structure, with a significant difference in arrival time
between parts of the spectrum, as demonstrated by the strong phase gradient
above $235\text{\,}\mathrm{nm}$.

To better understand the conditions necessary for the generation of short DUV
pulses, we have simulated the pulse propagation in the AR-PCF using the
single-mode forward Maxwell equation [17]. The model includes the waveguide
[18] and gas [19] dispersion, the Kerr nonlinearity [20], and the effect of
photoionisation [21] using the Perelomov-Popov-Terent’ev ionisation rate [22].
The initial pulse energy was determined from the experimental energy and the
transmission of the AR-PCF.

The DUV spectrum obtained by propagating the measured input pulse through a
waveguide with the same parameters as used in the experiment at
$6\text{\,}\mathrm{bar}$ argon pressure is shown in Fig. 4(a). The shape of
the spectrum as well as the phase are in good agreement with the XFROG
measurement, including a distinct difference in arrival time between the main
peak and a long-wavelength pedestal. As shown in Fig. 4(b), the dispersive
wave is emitted around $5\text{\,}\mathrm{cm}$ into the fibre. Due to the
pressure gradient, the phase-matched wavelength changes over the course of the
emission process, with lower pressures leading to phase-matching at shorter
wavelengths [7]. As a consequence the DUV pulse shifts and broadens spectrally
as it is being generated, decreasing the transform-limited pulse duration of
the dispersive wave. The lower-frequency pedestal appears further into the
propagation and is displaced from the main pulse due to group-velocity
mismatch. A result of complex generation dynamics, the pulse shape of the
dispersive wave cannot be determined from its power spectrum alone. Accurate
pulse characterisation is thus very important for the use of DWE as a
practical source for ultrafast science.

In the time-domain representation, shown in Fig. 4(c), the self-compression of
the driving pulse is clearly visible; DWE occurs at the point of maximum
compression. The curved trajectory of the DUV pulse, revealed by bandpass
filtering, shows that it begins to catch up to the driving pulse towards the
end of the propagation after falling behind initially. This is due to the
pressure gradient. As demonstrated by the shifting zero-dispersion wavelength,
shown in Fig. 4(b), the dispersion landscape changes; the dispersive wave is
accelerated relative to the driving pulse as a consequence. Importantly,
despite several centimetres of propagation, the duration of the main DUV pulse
remains well below $4\text{\,}\mathrm{fs}$—the gas pressure is low and the DUV
pulse experiences negligible dispersion. This is in contrast to the case of
uniform gas pressure, where the dispersive wave can be stretched to several
times its initial duration unless the point of maximum compression is very
close to the end of the waveguide [8], making the generation process more
sensitive to fluctuations of the input pulse energy. Since it also removes the
need for an exit window as well as any propagation through a pressurised gas
cell, the use of a pressure gradient is critical for the use of DWE as a
source of few-fs pulses.

In summary, we have demonstrated the generation, delivery and characterisation
of compressed, tuneable few-femtosecond pulses from DWE for the first time in
what is also the first application of ptychographic retrieval to non-SFG
XFROG. Our measurements as well as numerical simulations show that to obtain
such short pulses, a pressure gradient in the waveguide has to be used. With
vacuum propagation enabling the dispersion-free delivery of few-femtosecond
DUV pulses, the apparatus is directly compatible with application in pump-
probe schemes. Our work fully establishes DWE as a powerful light source for
ultrafast measurements with unprecedented spectral coverage and temporal
resolution.

## Acknowledgements

We acknowledge valuable assistance from T. Barnard, A. Ermolov, A. Gregory, S.
Jarosch, E. Larsen, C. O’Donovan, S. Parker, C. Strüber, and P. Ye, as well as
funding from the following sources: Engineering and Physical Sciences Research
Council (EPSRC) (EP/I032517/1, EP/R019509/1); European Research Council (ASTEX
290467, HISOL 679649); EPSRC/Defence Science and Technology Laboratory
Multidisciplinary University Research Initiative (EP/N018680/1); EPSRC
Doctoral Training Account

## References

  * [1] Albert Stolow, Arthur E. Bragg, and Daniel M. Neumark.  Femtosecond Time-Resolved Photoelectron Spectroscopy.  Chemical Reviews, 104(4):1719–1758, April 2004.
  * [2] Paul Hockett, Christer Z. Bisgaard, Owen J. Clarkin, and Albert Stolow.  Time-resolved imaging of purely valence-electron dynamics during a chemical reaction.  Nature Physics, 7(8):612–615, August 2011.
  * [3] Florentin Reiter, Ulrich Graf, Martin Schultze, Wolfgang Schweinberger, Hartmut Schröder, Nicholas Karpowicz, Abdallah Mohammed Azzeer, Reinhard Kienberger, Ferenc Krausz, and Eleftherios Goulielmakis.  Generation of sub-3 fs pulses in the deep ultraviolet.  Optics Letters, 35(13):2248, 2010.
  * [4] P. K. A. Wai, C. R. Menyuk, Y. C. Lee, and H. H. Chen.  Nonlinear pulse propagation in the neighborhood of the zero-dispersion wavelength of monomode optical fibers.  Optics Letters, 11(7):464, July 1986.
  * [5] Cristian Redondo Lourés, Thomas Roger, Daniele Faccio, and Fabio Biancalana.  Superresonant Radiation Stimulated by Higher Harmonics.  Physical Review Letters, 118(4):043902, January 2017.
  * [6] N. Y. Joly, J. Nold, W. Chang, P. Hölzer, A. Nazarkin, G. K L Wong, F. Biancalana, and P. St J Russell.  Bright spatially coherent wavelength-tunable deep-UV laser source using an Ar-filled photonic crystal fiber.  Physical Review Letters, 106(20):203901, 2011.
  * [7] Ka Fai Mak, John C Travers, Philipp Hölzer, Nicolas Y Joly, and Philip St J Russell.  Tunable vacuum-UV to visible ultrafast pulse source based on gas-filled Kagome-PCF.  Optics Express, 21(9):10942–53, 2013.
  * [8] John C Travers, Wonkeun Chang, Johannes Nold, Nicolas Y Joly, and Philip St. J. Russell.  Ultrafast nonlinear optics in gas-filled hollow-core photonic crystal fibers [Invited].  Journal of the Optical Society of America B, 28(12):A11–A26, 2011\.
  * [9] Alexey Ermolov, Heli Valtna-Lukner, John Travers, and Philip St.J. Russell.  Characterization of few-fs deep-UV dispersive waves by ultra-broadband transient-grating XFROG.  Optics Letters, 41(23):5535, 2016.
  * [10] Felix Köttig, Francesco Tani, Patrick Uebel, Philip St J. Russell, and John C. Travers.  In CLEO Europe 2015, paper PD_A_7. Optical Society of America, June 2015.
  * [11] Felix Köttig, Francesco Tani, Christian Martens Biersach, John C. Travers, and Philip St.J. Russell.  Generation of microjoule pulses in the deep ultraviolet at megahertz repetition rates.  Optica, 4(10):1272, October 2017.
  * [12] S. Linden, J. Kuhl, and H. Giessen.  Amplitude and phase characterization of weak blue ultrashort pulses by downconversion.  Optics Letters, 24(8):569–571, 1999.
  * [13] T Witting, D Greening, D Walke, P Matia-Hernando, T Barillot, J P Marangos, and J W G Tisch.  Time-domain ptychography of over-octave-spanning laser pulses in the single-cycle regime.  Opt. Lett., 41(18):4218–4221, 2016.
  * [14] Dirk-Mathys Spangenberg, Michael Brügmann, Erich Rohwer, and Thomas Feurer.  All-optical implementation of a time-domain ptychographic pulse reconstruction setup.  Applied Optics, 55(21):5771, 2016.
  * [15] Andrew Maiden, Daniel Johnson, and Peng Li.  Further improvements to the ptychographical iterative engine.  Optica, 4(7):736, July 2017.
  * [16] Pavel Sidorenko, Oren Lahav, Zohar Avnat, and Oren Cohen.  Ptychographic reconstruction algorithm for frequency-resolved optical gating: super-resolution and supreme robustness.  Optica, 3(12):1320, December 2016.
  * [17] A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Ramírez-Góngora, and M. Kolesik.  Practitioner’s guide to laser pulse propagation models and simulation.  The European Physical Journal Special Topics, 199(1):5–76, 2011\.
  * [18] E. A. J. Marcatili and R. A. Schmeltzer.  Hollow Metallic and Dielectric Waveguides for Long Distance Optical Transmission and Lasers.  Bell System Technical Journal, 43(4):1783–1809, July 1964.
  * [19] A. Börzsönyi, Z. Heiner, M. P. Kalashnikov, A. P. Kovács, and K. Osvay.  Dispersion measurement of inert gases and gas mixtures at 800 nm.  Applied Optics, 47(27):4856, September 2008.
  * [20] David P. Shelton and Julia E. Rice.  Measurements and Calculations of the Hyperpolarizabilities of Atoms and Small Molecules in the Gas Phase.  Chemical Reviews, 94(1):3–29, 1994.
  * [21] M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec.  Light Propagation in Field-Ionizing Media: Extreme Nonlinear Optics.  Physical Review Letters, 83(15):2930–2933, October 1999.
  * [22] A M Perelomov, V S Popov, and M V Terent ’ev.  Ionization of atoms in an alternating electric field.  Soviet Physics JETP, 23(50):1393–1409, 1966.


Preheating with deep learning

Jong-Hyun Yoon, Simon Cléry, Mathieu Gross and Yann Mambrini

Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

Abstract

We apply deep learning techniques to the late-time turbulent regime in a post-
inflationary model where a real scalar inflaton field and the standard model
Higgs doublet interact with renormalizable couplings between them. After
inflation, the inflaton decays into the Higgs through a trilinear coupling and
the Higgs field subsequently thermalizes with gauge bosons via its
$SU(2)\times U(1)$ gauge interaction. Depending on the strength of the
trilinear interaction and the Higgs self-coupling, the effective mass squared
of Higgs can become negative, leading to the tachyonic production of Higgs
particles. These produced Higgs particles would then share their energy with
gauge bosons, potentially indicating thermalization. Since the model entails
different non-perturbative effects, it is necessary to resort to numerical and
semi-classical techniques. However, simulations require significant costs in
terms of time and computational resources depending on the model used.
Particularly, when $SU(2)$ gauge interactions are introduced, this becomes
evident as the gauge field redistributes particle energies through
rescattering processes, leading to an abundance of UV modes that disrupt
simulation stability. This necessitates very small lattice spacings, resulting
in exceedingly long simulation runtimes. Furthermore, the late-time behavior
of preheating dynamics exhibits a universal form by wave kinetic theory.
Therefore, we analyze patterns in the flow of particle numbers and predict
future behavior using CNN-LSTM (Convolutional Neural Network combined with
Long Short-Term Memory) time series analysis. In this way, we can reduce our
dependence on simulations by orders of magnitude in terms of time and
computational resources.

###### Contents

  1. 1 Introduction
  2. 2 Tachyonic production of Higgs
  3. 3 Post-inflationary weak interactions
  4. 4 CNN-LSTM time series analysis
  5. 5 Conclusion
  6. A $SU(2)\times U(1)$ on the lattice
  7. B Hyperparameters

##  1 Introduction

Cosmic inflation has been demonstrated to play a pivotal role in establishing
the initial conditions of the early universe within modern cosmology [1, 2, 3,
4]. These characteristics can be realized by a real scalar field, known as the
inflaton. It allows for the exploration of various intriguing phenomenologies
due to the interplay between the inflaton and other fields [5, 6, 7]. However,
in any post-inflationary model, we must eventually explain a thermal universe
where the constituents of thermal plasma possess distinctive properties such
as thermal mass, necessitating a deeper comprehension of thermal field theory.
During the initial stages of thermalization, defining thermal quantities
becomes challenging, particularly with a low number of plasma constituents
[8]. In the initially under-occupied systems, low-momentum particles form a
bath, thermalize, and then dominate the system by absorbing the energy from
the hard excitations. This process is dominated by number-changing
interactions, which can interfere and therefore it is important to take into
account modification of the interaction rate so-called the
Landau–Pomeranchuk–Migdal (LPM) effect [9, 10, 11]. When the rate of plasma
constituent production is sufficiently low, interaction dynamics can be
treated perturbatively [12, 13], but when the rate is large, the system falls
into the non-perturbative regime, where fields redistribute their energies by
turbulent rescattering. This late-time phenomenon is known as free turbulence,
where the dynamics are described by kinetic theories [14, 15]. Given the
inherent challenges of demonstrating free turbulence and tracking
thermalization using lattice simulations, it is tempting to infer patterns in
the occupation number of fields, represented as self-similarity implied by the
kinetic theories and analyze them using deep learning methods for its
proficiency in pattern recognition.

Machine learning and deep learning have revolutionized various fields,
including physics, by providing powerful tools to analyze complex systems and
phenomena [16, 17, 18]. Machine learning involves algorithms that enable
computers to learn from data and make predictions or decisions without being
explicitly programmed, while deep learning utilizes neural networks with
multiple layers to learn complex patterns and representations from data.
Convolutional Neural Networks (CNNs) are particularly effective for processing
spatial data, such as images, while Long Short-Term Memory (LSTM) networks
specialize in capturing temporal dependencies in sequential data [19].

Therefore, combining Convolutional 1D (Conv1D) layers with LSTM networks
offers a powerful solution for analyzing time series data in physics. Conv1D
layers effectively capture local patterns in sequential data, while LSTM
networks model long-term dependencies and temporal dynamics, enabling accurate
analysis of time-dependent datasets across various domains.

In this study, we apply deep learning techniques to the following example
model. We consider the post-inflationary dynamics of standard model Higgs
particles abundantly produced by the inflaton, where the Higgs subsequently
decays into gauge bosons through $SU(2)\times U(1)$ interactions, leading to
the formation of a thermal sector. During the inflaton oscillation period,
numerical simulations are often necessary since the production of scalar or
gauge bosons after inflation involves non-perturbative dynamics (e.g.
preheating)[20, 21, 22, 23, 24, 25]. Since lattice simulations designed to
analyze the preheating system have limitations in observing late-time dynamics
up to reheating, we attempt to overcome these limitations by employing deep
learning techniques.

Our paper is structured as follows: In Section 2, we examine the tachyonic
production mechanism of Higgs via its trilinear interaction with the inflaton.
In Appendix A, we investigate how post-inflationary $SU(2)\times U(1)$ gauge
interactions are implemented on the lattice, and the simulation results will
be presented in Section 3. We discuss the deep learning technique that we
applied to our results in Section 4 and draw conclusions in Section 5.
Appendix B provides an introduction to the hyperparameters used in deep
learning.

##  2 Tachyonic production of Higgs

The minimal interactions between the inflaton $\phi$ and the Higgs $H$ are of
the form [32]

$\Delta
V=\frac{1}{2}{m^{2}_{\phi}}\,\phi^{2}+\frac{1}{4}{\lambda_{\phi}}\,\phi^{4}+\frac{1}{2}{\lambda_{\phi
h}}\,\phi^{2}\,H^{\dagger}H+{\sigma_{\phi
h}}\,\phi\,H^{\dagger}H-{m^{2}_{h}}\,H^{\dagger}H+{\lambda_{h}}\,(H^{\dagger}H)^{2}\;,$
(1)

where we assume $\lambda_{\phi}=0,\,\lambda_{\phi h}=0,\,$ and $m_{h}=0$ for
simplicity.111In the case of the quartic inflaton potential, to elevate the
global minimum of the potential from negative values to zero, one should
introduce an undesirable offset term [27]. In unitary gauge

$H=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\ h\end{pmatrix}\,,$ (2)

the potential simplifies to

$V=\frac{1}{2}{m^{2}_{\phi}}\,\phi^{2}+\frac{\sigma_{\phi
h}}{2}\,\phi\,h^{2}+\frac{\lambda_{h}}{4}\,h^{4}\,.$ (3)

By rewriting Eq. (3) in the form

$V=\frac{1}{2}\left({m_{\phi}}\,\phi+\frac{\sigma_{\phi
h}}{2m_{\phi}}h^{2}\right)^{2}+\frac{1}{4}\,\left(\lambda_{h}-\frac{\sigma^{2}_{\phi
h}}{2m^{2}_{\phi}}\right)h^{4}\,,$ (4)

one finds that $\lambda_{h}$ must be larger than $\sigma^{2}_{\phi
h}/2m^{2}_{\phi}$ to make the potential bounded from below [26]. The equation
of motion of $h$ is

$\Box h+V_{,h}=0$ (5)

Quantizing $h$ with the eigenmodes $h_{k}(t)e^{i\mathbb{k}\vec{x}}$ and
$\mathbb{k}$ the comoving momentum, one may derive linear mode equations in
momentum space222The modes are rescaled to factor out the scale factor,
$h_{k}\rightarrow a^{3/2}h_{k}$.

$\ddot{h}_{k}+\omega^{2}_{k}h_{k}=0,$ (6)

where

$\omega^{2}_{k}=\frac{k^{2}}{a^{2}}+\sigma_{\phi h}\phi(t)+3\lambda_{h}\langle
h^{2}\rangle-3\dot{a}^{2}/4a^{2}-3\ddot{a}/2a\,.$ (7)

Given that negative values of the inflaton background can make the entire
frequency term negative, we expect a tachyonic growth of the Higgs mode during
the inflaton oscillation period. While one can theoretically examine the
mathematical structure of the mode equations and quantitatively determine the
number of produced particles, practical implementation proves to be
challenging for various reasons [27]. Each mode in the equation experiences
movement across the instability band with varying Floquet quotient due to the
expansion of space. Eq. (6) becomes linearized in terms of $k$ only when we
disregard the rescattering of modes with different $k$. This is not the case
in the presence of the strong self-interaction and gauge interaction of Higgs,
even at the onset of the inflaton oscillation period. To deal with the
turbulent particle interactions, we opt to solve the equations of motion
directly in spacetime through a semi-classical approach (see [28] for studies
on gravitational wave production in a similar framework).

##  3 Post-inflationary weak interactions

Given the gauge interactions in the standard model, Higgs particles will
interact with gauge bosons upon production. Previous studies have explored
scenarios where considerable Higgs condensates generate gauge bosons via
resonant channels [24, 25]. We examine a system in which field constituents,
particularly focusing on the Higgs directly produced by the inflaton, interact
with other fields in a turbulent manner.

Figure 1:  Energy densities of interacting fields during preheating. Left:
energy density of fields normalized by $m_{\phi}^{4}$. The $x$-axis represents
time rescaled by the inflaton mass, $m_{\phi}$. Right: ratio of the energy
density of each field to the total energy over the scale factor.

Using a publicly available code, CosmoLattice [29, 30], we conduct lattice
simulations to model post-inflationary weak interactions of the Higgs (see A).
We set $g_{2}=g_{1}=0.6$ and $\lambda_{h}=10^{-2}$ for the coupling constants
at the inflationary physics scale to be not too different from those of the
standard model. As in [28], the inflaton field starts to roll down to the
minimum of potential with the bare mass $m_{\phi}=1.6\times 10^{13}$ GeV,
$\phi_{init}=1\,m_{p}$ and $\dot{\phi}_{init}=-0.71\,m_{\phi}m_{p}$, where the
Planck mass $m_{p}\approx 2.435\times 10^{18}\text{ GeV}$. Here, $init$ refers
to the initial values input into the lattice simulations. While we do not
specify an inflationary model, it is assumed to converge to a quadratic
inflaton potential some time after the end of inflation.333Around the end of
inflation where $\phi>0$, the potential of the form Eq. (4) reaches its
minimum for $h=0$. Parameters may vary within certain limits to accommodate
running coupling constants or experimental constraints.

The energy densities of fields are defined

$\displaystyle\rho_{\phi}$ $\displaystyle\equiv$ $\displaystyle
K_{\phi}+G_{\phi}+V(\phi)\,$ (8) $\displaystyle\rho_{H}$ $\displaystyle\equiv$
$\displaystyle K_{H}+G_{H}+V(H)\,$ (9) $\displaystyle\rho_{B}$
$\displaystyle\equiv$ $\displaystyle K_{U(1)}+G_{U(1)}\,$ (10)
$\displaystyle\rho_{W}$ $\displaystyle\equiv$ $\displaystyle
K_{SU(2)}+G_{SU(2)}\,,$ (11)

where $V(\phi)$ and $V(H)$ represent the potential of the inflaton and Higgs
doublet, respectively. Fig. 1 shows the energy densities of fields over time
for the case where $\sigma_{\phi h}=10^{9}\,\text{GeV}$. The total energy
density of the system is still dominated by the energy of the inflaton field,
while the Higgs contributes only a small fraction of the total energy. Due to
the strong self-interaction of the Higgs ($\lambda_{h}=10^{-2}$), tachyonic
production is highly suppressed but its interaction with the inflaton is
sufficiently strong to withstand the dilution by the quadratic inflaton
potential ($\rho_{R}\propto a^{-4}$ for radiation while $\rho_{\phi}\propto
a^{-3}$ for the massive inflaton field).

For a qualitative analysis, we illustrate the energy ratio of the Higgs and
the $SU(2)$ gauge field to the $U(1)$ gauge field on the left of Fig. 2. We
anticipate observing the equipartition of energy in the gauge sector as an
indication of thermalization, in a ratio of 4:6:2 (or equivalently 2:3:1)
corresponding to the degrees of freedom of the Higgs, $SU(2)$ gauge fields,
and $U(1)$ gauge boson, respectively. We have captured the moment of energy
equipartition shortly after the end of inflation ($a\sim 60$ or $m_{\phi}t\sim
700$).

Figure 2:  Indication of a quasi-equilibrium state in the gauge sector. Left:
energy ratio of Higgs and $SU(2)$ gauge field to $U(1)$ gauge field. The
dashed grey horizontal line represents 3 which corresponds to the ratio of
degrees of freedom for the $SU(2)$ gauge field to the one for the $U(1)$ gauge
field ($\rho_{W}/\rho_{B}=3$). Likewise, the red dashed horizontal line
represents the same quantity but for the Higgs doublet
($\rho_{H}/\rho_{B}=2$). Right: occupation number for the Higgs doublet (red
to blue over time). The comoving momentum $k$ is normalized by the inflaton
mass. The mean particle distribution shifts from IR to UV as time progresses.

We excited initial particle modes for the comoving momentum in the range
$5\leq k/m_{\phi}\leq 10$ and evolved the system in the extended range $5\leq
k/m_{\phi}\lesssim 690$. It should be noted that our focus is on the tail end
of spectra in the UV, which could potentially underestimate Higgs production
in the IR. Expanding the momentum window in the IR direction leads to
increased energy exchange due to the increased scattering of particles,
necessitating further expansion in the UV direction as well. This results in a
significant increase in simulation runtime and it is one of the key
motivations for us to leverage deep learning techniques. For the subsequent
discussion, we define the occupation number for the Higgs doublet as444Zero
initial occupation number for scalar fields in theory corresponds to
$n_{k,init}=1/2$ in our definition. For the Higgs doublet, we take into
account 4 Higgs d.o.f and 1/2 for the normalization factor, resulting in
$n^{H}_{k,init}=1/2\times 4\times 1/2=1$.

$\displaystyle
n^{H}_{k}=\sum_{i}n^{i}_{k}\equiv\sum_{i}\frac{\omega^{i}_{k}}{2}\left(|{(h^{i}_{k})^{\prime}}|^{2}/({\omega^{i}_{k}})^{2}+|h^{i}_{k}|^{2}\right)\,,$
(12)

where $i$ represents the Higgs $d.o.f$ and the frequency is

$\displaystyle\omega^{i}_{k}\equiv\sqrt{k^{2}+a^{2}\,|d^{2}V(H)/d{h^{i}}^{2}|}\,.$
(13)

The prime denotes the conformal time derivative $a(d/dt)$. For gauge bosons,
since the occupation number is gauge-dependent, gauge fixing must be performed
at each measurement [25, 33, 34, 35]. Due to its technical difficulty, we
leave this aspect for future work.

The late-time classical reheating process of self-interacting fields or non-
abelian fields exhibits a turbulent and self-similar evolution of distribution
functions towards equilibrium [15, 36]. The shape of the spectra, along with
the self-similar dynamics, can be comprehended within the framework of wave
kinetic theory. This suggests the presence of a consistent trend in the
evolution of the occupation number, offering a foundation for exploring the
occupation number using deep learning methods.

##  4 CNN-LSTM time series analysis

The CNN-LSTM model integrates two fundamental components: the Convolutional
Neural Network (CNN) and the Long Short-Term Memory (LSTM) network. Each
component plays a distinct role in processing the input data to extract
features and capture temporal dependencies.

The CNN component is responsible for extracting spatial features from the
input data. It operates by applying convolutional filters over the input data
to capture spatial patterns, such as edges, textures, and shapes. This process
is facilitated by multiple convolutional layers followed by pooling layers,
which help reduce spatial dimensions while retaining essential information.

On the other hand, the LSTM component specializes in modeling temporal
dependencies in sequential data. It consists of stacked LSTM layers, each
equipped with memory cells capable of retaining information over extended time
periods. These memory cells are governed by gates that control the flow of
information, enabling the model to selectively update, forget, or output
information based on the input sequence.

In the parallel CNN-LSTM model, the input data are processed separately by the
CNN and LSTM components. Once the CNN and LSTM components have processed their
respective data streams, the output features are concatenated and fed into a
linear layer (or layers) for further processing. This linear layer serves as
the final step in the model, where the extracted features are transformed and
mapped to the desired output space.

In our CNN-LSTM time series analysis, the input data have the shape
(batch_size, seq_length, input_size), which corresponds to (the number of time
points, the length of sequence, the number of $k$ modes) for the preheating
model that we discuss. Both the input and output contain $n_{k}(t)$ of the
Higgs, but the model output predicts values occurring after one sequence.555In
a time series context, the term ‘time’ typically denotes the sequential order
of data points rather than implying that the input variables must directly
represent time intervals. We implement the model using PyTorch [37], a popular
deep learning library in Python. The model architecture consists of a Conv1D
layer followed by an LSTM layer. The Conv1D layer captures local patterns in
the sequential data, while the LSTM layer captures long-term dependencies.

The model is trained using the training data, and the training process
involves optimizing a loss function (mean absolute error) using the Adam
optimizer. The training loop iterates over epochs, updating the model
parameters to minimize the loss. Predictions are made for future time steps
beyond the training data and compared with the actual test data. The test data
are completely separated from the model training and forecasting process.

The free turbulence starts approximately around $m_{\phi}t\approx 400$, so we
divided the data into a training set for $400<m_{\phi}t<1000$ and a test set
for $1000<m_{\phi}t<1500$. We make predictions up to $m_{\phi}t\approx 2000$.

Figure 3:  Occupation number for the Higgs field as a function of $t$ with a
logarithmic scale for two $k$ modes. The gray dashed vertical line
distinguishes between training data (orange dotted curve) and time series
predictions (green solid curve). The trained predictions (blue solid curve)
are obtained by feeding the input data into the model, while the time series
predictions are recursively generated by applying the model to the output.
Test data (red dotted curve) are used to calculate the validation loss. Left:
$n_{k}(t)$ for $k\approx 5m_{\phi}$. Since it corresponds to the region
interacting with the inflaton, the pattern appears to be more complex. Right:
$n_{k}(t)$ for $k\approx 400m_{\phi}$.

The model’s predictions are visualized alongside the actual data as in Fig. 3.
It represents the actual data used for training, the trained predictions, the
test data not used for training, and the predictions for future time steps for
two momentum modes among 15 modes ($5\leq k/m_{\phi}\lesssim 400$). It is
remarkable that our deep learning model managed to learn patterns and trends
in the limited amount of data and make predictions similar to actual
simulation results.

The predictions for each mode are collected and reshaped in Fig. 4. At the end
of our predictions, around $m_{\phi}t\approx 2000$, it appears that we have
already reached a point where the system departs from the classical regime and
quantum effects become significant, indicating proximity to thermalization. By
further developing the deep learning model, we can implement a more realistic
cosmological model that correctly incorporates the production of the Higgs,
which has been underestimated in the IR regime.

Figure 4:  Predicted spectra of the Higgs occupation number. Left: Occupation
numbers computed for different momentum modes (15 points in total) are
collected and reshaped. The darker magenta color indicates a later spectrum in
time ($500\lesssim m_{\phi}t\lesssim 2000$) and the cyan solid (dashed) curve
represents a spectrum at $m_{\phi}t\approx 1000$ ($m_{\phi}t\approx 1500$).
Right: Self-similar form of the spectra corresponding to the left-hand side of
Eq. (14). Only predicted values outside the data range were used to fit $p$
and $q$ ($1500\lesssim m_{\phi}t\lesssim 2000$).

Still, we may use the predicted values to verify the self-similar evolution.
Rewriting Eq. (4) from [15] yields

$\displaystyle\tilde{k}^{4}n(k=\tilde{k}\tau^{p},t)\tau^{q}=\tilde{k}^{4}n(k=\tilde{k},t=t_{0})~{},$
(14)

where $\tilde{k}\equiv k\tau^{-p}$, $\tau\equiv t/t_{0}$, and $t_{0}$ denotes
a fixed point in time. $p$ and $q$ represent the fitting exponents, and $p$
determines how quickly the system approaches equilibrium. Since the right-hand
side of Eq. (14) is stationary, it was used as data for fitting, resulting in
a good fit when $p\approx 0.16$ and $q\approx 0.63$, as illustrated in Fig 4.
The exponents may vary depending on the details of the preheating model and
the fitting conditions.

##  5 Conclusion

We consider a simple post-inflationary scenario where after inflation, the
Standard Model Higgs decays into gauge bosons via weak interactions while
forming a quasi-thermal state in the gauge sector. The non-perturbative nature
of the system is mainly due to the $SU(2)$ gauge interaction ($g=0.6$) and
Higgs self-coupling ($\lambda_{h}=10^{-2}$) along with the trilinear
interaction ($\sigma_{\phi h}\phi H^{\dagger}H$ where $\sigma_{\phi
h}=10^{9}\,$GeV) between Higgs and the inflaton. Analyzing the production and
decay of Higgs in this model appears challenging without relying on numerical
and semi-classical approaches. Demonstrating thermalization through simulation
is also very challenging, primarily because simulations are costly in
practical terms and quantum effects become significant in later stages. We
introduced deep learning techniques as an effort to overcome the limitations
of classical lattice simulations and ultimately reproduce the thermal
universe. The later stage of preheating is referred to as the regime of free
turbulence, where the flow of particle numbers can be analyzed through kinetic
theory, making it a suitable material for applying deep learning techniques.
LSTM networks are commonly used in time series analysis, but they have
limitations in predicting beyond the trained range. And since the momentum
modes of the fields in the free turbulence regime have localized influence
rather than behaving independently, local patterns in the momentum space can
play an important role. CNNs are useful for capturing localized patterns
observed in images and similar data. Therefore, we utilized a deep learning
model that combines LSTM with CNN, allowing a fully connected linear layer to
learn the outputs of both networks simultaneously. We applied this model to
the realistic post-inflationary preheating scenario and predicted the spectra
of the Higgs field beyond simulation results. While still in its early stages,
this study has demonstrated its potential for advancement.

Acknowledgements. J.Y. would like to thank Oleg Lebedev for the useful
discussion. This work was performed using HPC resources from the
‘‘Mésocentre’’ computing center of CentraleSupélec, École Normale Supérieure
Paris-Saclay and Université Paris-Saclay supported by CNRS and Région Île-de-
France (https://mesocentre.universite-paris-scalay.fr/). This project has
received support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
860881-HIDDeN, and the CNRS-IRP UCMN.

##  Appendix A $SU(2)\times U(1)$ on the lattice

Let us start with the $SU(2)\times U(1)$ gauge Lagrangian

$\mathcal{L}_{\textrm{g}auge}=-\frac{1}{4}{(W^{a}_{\mu\nu})^{2}}\,-\frac{1}{4}{B^{2}_{\mu\nu}}\,+(D_{\mu}H)^{\dagger}(D_{\mu}H)+{m^{2}_{h}}\,H^{\dagger}H-{\lambda_{h}}\,(H^{\dagger}H)^{2}\;,$
(15)

where $B_{\mu}$ is the $U(1)$ gauge boson, with
$B_{\mu\nu}=\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$, and $W^{a}_{\mu}$
are the $SU(2)$ gauge bosons, with $W^{a}_{\mu\nu}$ their field strengths. The
covariant derivative of Higgs doublet $H$ is

$D_{\mu}H=\partial_{\mu}H-ig_{2}W^{a}_{\mu}\tau^{a}H-\frac{1}{2}ig_{1}B_{\mu}H,$
(16)

where $g_{2}$ and $g_{1}$ are the $SU(2)$ and $U(1)$ couplings, respectively.

To present the equations of motion in an expanding universe, we extract useful
expressions from [29]. First, we assume a flat Friedmann-Lemaître-Robertson-
Walker (FLRW) metric, characterized by the line element:

$\text{d}s^{2}=g_{\mu\nu}\text{d}x^{\mu}\text{d}x^{\nu}=-a(\eta)^{2\alpha}\text{d}\eta^{2}+a(\eta)^{2}\delta_{ij}\text{d}x^{i}\text{d}x^{j}\
,$ (17)

Here, $a(\eta)$ represents the scale factor, $\delta_{ij}$ denotes the
Euclidean metric, and $\alpha$ is a constant parameter chosen conveniently
depending on the scenario.

The EOMs in an expanding universe are given by

$\displaystyle\phi^{\prime\prime}-a^{-2(1-\alpha)}{\vv\nabla}^{\,2}\phi+(3-\alpha)\frac{{a^{\prime}}}{a}{\phi^{\prime}}$
$\displaystyle=$ $\displaystyle-a^{2\alpha}V_{,\phi}\ ,$ (18) $\displaystyle
H^{\prime\prime}-a^{-2(1-\alpha)}{\vv
D}^{\,2}H+(3-\alpha)\frac{{a^{\prime}}}{a}{H^{\prime}}$ $\displaystyle=$
$\displaystyle-\frac{a^{2\alpha}V_{,|H|}}{2}\frac{H}{|H|}\ ,$ (19)
$\displaystyle\partial_{0}F_{0i}-a^{-2(1-\alpha)}\partial_{j}F_{ji}+(1-\alpha)\frac{{a^{\prime}}}{a}F_{0i}$
$\displaystyle=$ $\displaystyle a^{2\alpha}J^{B}_{i}\ ,$ (20)
$\displaystyle(\mathcal{D}_{0})_{ab}(G_{0i})^{b}-a^{-2(1-\alpha)}(\mathcal{D}_{j})_{ab}(G_{ji})^{b}+(1-\alpha)\frac{{a^{\prime}}}{a}(G_{0i})^{b}$
$\displaystyle=$ $\displaystyle a^{2\alpha}(J_{i})_{a}\ ,$ (21)
$\displaystyle\partial_{i}F_{0i}$ $\displaystyle=$ $\displaystyle
a^{2}J^{B}_{0}\ ,$ (22) $\displaystyle(\mathcal{D}_{i})_{ab}(G_{0i})^{b}$
$\displaystyle=$ $\displaystyle a^{2}(J_{0})_{a}\,,$ (23)

where the prime ′ denotes $\alpha$-time derivative $d/d\eta$ and
$(\mathcal{D}_{\nu}O)_{a}=(\mathcal{D}_{\nu})_{ab}O_{b}\equiv(\delta_{ab}\partial_{\nu}-f_{abc}W_{\nu}^{c})O_{b}$.
The field strength tensors are defined as

$\displaystyle F_{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}\ ,$ (24)
$\displaystyle G_{\mu\nu}$ $\displaystyle\equiv$
$\displaystyle\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}-ig_{2}[W_{\mu},W_{\nu}]\,,$
(25)

and the currents are

$\displaystyle\hskip 51.21504pt(J^{B})^{\mu}$ $\displaystyle\equiv$
$\displaystyle g_{1}\mathcal{I}m[H^{\dagger}(D^{\mu}H)]\,,$ (26)
$\displaystyle\hskip 51.21504ptJ_{a}^{\mu}$ $\displaystyle\equiv$
$\displaystyle 2g_{2}\mathcal{I}m[H^{{\dagger}}T_{a}(D^{\mu}H)]\,.$ (27)

Eqs. (22) and (23) are the $U(1)$ and $SU(2)$ Gauss constraints in an
expanding background that can be used to check the consistency of numerical
simulations. The evolution of the scale factor obeys the Friedmann equations

$\displaystyle\left({a^{\prime}\over a}\right)^{2}$ $\displaystyle=$
$\displaystyle\frac{a^{2\alpha}}{3m_{p}^{2}}\left\langle{K}_{\phi}+{K}_{H}+{G}_{\phi}+{G}_{H}+{K}_{U(1)}+{G}_{U(1)}+{K}_{SU(2)}+{G}_{SU(2)}+{V}\right\rangle\,,$
(28) $\displaystyle{a^{\prime\prime}\over a}$ $\displaystyle=$
$\displaystyle\frac{a^{2\alpha}}{3m_{p}^{2}}\left\langle(\alpha-2)({K}_{\phi}+{K}_{H})+\alpha({G}_{\phi}+{G}_{H})+(\alpha+1)V\right.$
$\displaystyle\hskip
56.9055pt\left.+~{}(\alpha-1)({K}_{U(1)}+{G}_{U(1)}+{K}_{SU(2)}+{G}_{SU(2)})\right\rangle\,,$

where $\langle...\rangle$ represents values averaged over space. One of the
two Friedmann equations is redundant and used to check consistency as Gauss
constraints. The kinetic ($K$) and gradient ($G$) energy densities are

$\displaystyle\begin{array}[]{lcl}{K}_{\phi}&=&\frac{1}{2a^{2\alpha}}{\phi^{\prime}}^{2}\vspace{0.1cm}\vspace{0.1cm}\\\
{K}_{H}&=&\frac{1}{a^{2\alpha}}(D_{0}H)^{\dagger}(D_{0}H)\vspace{0.1cm}\\\
{K}_{U(1)}&=&\frac{1}{2a^{2+2\alpha}}\sum_{i}F_{0i}^{2}\vspace{0.1cm}\\\
{K}_{SU(2)}&=&\frac{1}{2a^{2+2\alpha}}\sum_{a,i}(G_{0i}^{a})^{2}\vspace{0.1cm}\\\
\end{array}\hskip 2.84544pt;\hskip
11.38092pt\begin{array}[]{lcl}{G}_{\phi}&=&\frac{1}{2a^{2}}\sum_{i}(\partial_{i}\phi)^{2}\vspace{0.1cm}\\\
{G}_{H}&=&\frac{1}{a^{2}}\sum_{i}(D_{i}H)^{\dagger}(D_{i}H)\vspace{0.1cm}\\\
{G}_{U(1)}&=&\frac{1}{2a^{4}}\sum_{i,j<i}F_{ij}^{2}\vspace{0.1cm}\\\
{G}_{SU(2)}&=&\frac{1}{2a^{4}}\sum_{a,i,j<i}(G_{ij}^{a})^{2}\,.\vspace*{0.2cm}\\\
\end{array}$ (38)

One can solve the above equations in the temporal gauge with different
equation-solving algorithms. We implement our model on the lattice using a
publicly available code called CosmoLattice [29, 30].

Program variables used for our physics model are rescaled as follows:

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle
0~{}~{},~{}~{}f_{*}=1.6\cdot
10^{13}\,\text{GeV}~{}~{},~{}~{}\omega_{*}=1.6\cdot
10^{13}\,\text{GeV}~{}~{},~{}~{}d\tau=a^{-\alpha}\omega_{*}dt~{}~{},$ (39)
$\displaystyle\tilde{f}$ $\displaystyle\equiv$ $\displaystyle
f/f_{*}~{}~{},~{}~{}\tilde{V}\equiv V/(f^{2}_{*}\omega^{2}_{*})~{}~{},$ (40)

where variables with tildes represent dimensionless program variables.
Important input parameters used for our model are set as follows

evolver=VV2 , N=160 , dt=0.005 , kIR = 5 , kCutOff = 10.0 . (41)

We evolve the field equations using the second kind of Velocity-Verlet
algorithm on 160 lattice points per dimension. dt represents the time step in
program unit.

##  Appendix B Hyperparameters

Figure 5:  Train and validation losses versus epoch. Early stopping was
performed at the epoch marked by the gray dashed vertical line.

Hyperparameters are parameters whose values are set before the learning
process begins. They control aspects of the learning process itself and are
not learned from the data. These parameters are external to the model and must
be set based on heuristics, prior experience, or trial and error. Common
hyperparameters include the learning rate, batch size, number of hidden
layers, number of neurons in each layer, dropout rate, regularization
strength, and optimizer choice. Proper selection and tuning of hyperparameters
are crucial for achieving good performance and generalization of machine
learning models. If hyperparameters are not tuned, the results obtained
through the model can either underfit or overfit the data. Therefore, it is
important to compute and monitor the train loss and validation loss.
Typically, the train loss decreases over epochs, so by solely examining the
validation loss, we can discern whether the model is underfitting or
overfitting. If the validation loss continues to decrease, it indicates
underfitting, suggesting room for improvement through further training.
However, if the validation loss does not decrease, it implies overfitting,
where the model merely memorizes the data without understanding beyond it.
Thus, it is advisable to perform early stopping at this intermediate stage, as
illustrated in Fig 5.

## References

  * [1] A. A. Starobinsky, Phys. Lett. B 91 (1980) 99-102.
  * [2] A. H. Guth, Phys. Rev. D 23 (1981) 347-356.
  * [3] A. D. Linde, Phys. Lett. B 108 (1982) 389-393; Phys. Lett. B 129 (1983), 177-181.
  * [4] V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33, 532-535 (1981).
  * [5] O. Lebedev and J. H. Yoon, Phys. Lett. B 821, 136614 (2021).
  * [6] O. Lebedev, T. Nerdi, T. Solomko and J. H. Yoon, Phys. Rev. D 106, no.4, 043537 (2022).
  * [7] O. Lebedev, T. Solomko and J. H. Yoon, JCAP 02, 035 (2023).
  * [8] A. Kurkela and G. D. Moore, JHEP 12, 044 (2011) doi:10.1007/JHEP12(2011)044 [arXiv:1107.5050 [hep-ph]].
  * [9] L. D. Landau and I. Pomeranchuk, Dokl. Akad. Nauk Ser. Fiz. 92, 535-536 (1953)
  * [10] A. B. Migdal, Phys. Rev. 103, 1811-1820 (1956) doi:10.1103/PhysRev.103.1811
  * [11] M. Gyulassy and X. n. Wang, Nucl. Phys. B 420, 583-614 (1994) doi:10.1016/0550-3213(94)90079-5 [arXiv:nucl-th/9306003 [nucl-th]].
  * [12] K. Harigaya and K. Mukaida, JHEP 05, 006 (2014) doi:10.1007/JHEP05(2014)006 [arXiv:1312.3097 [hep-ph]].
  * [13] K. Mukaida and M. Yamada, [arXiv:2402.14054 [hep-ph]].
  * [14] R. Micha and I. I. Tkachev, Phys. Rev. Lett. 90, 121301 (2003) doi:10.1103/PhysRevLett.90.121301 [arXiv:hep-ph/0210202 [hep-ph]].
  * [15] R. Micha and I. I. Tkachev, Phys. Rev. D 70, 043538 (2004) doi:10.1103/PhysRevD.70.043538 [arXiv:hep-ph/0403101 [hep-ph]].
  * [16] Y. LeCun, Y. Bengio and G. Hinton, Nature 521, 436-444 (2015) doi:10.1038/nature14539
  * [17] J. Schmidhuber, Neural Networks 61, 85-117 (2015) doi:10.1016/j.neunet.2014.09.003
  * [18] T. Plehn, A. Butter, B. Dillon, T. Heimel, C. Krause and R. Winterhalder, [arXiv:2211.01421 [hep-ph]].
  * [19] S. Hochreiter and J. Schmidhuber, Neural Comput. 9, no.8, 1735-1780 (1997) doi:10.1162/neco.1997.9.8.1735
  * [20] A. Rajantie, P. M. Saffin and E. J. Copeland, Phys. Rev. D 63, 123512 (2001) doi:10.1103/PhysRevD.63.123512 [arXiv:hep-ph/0012097 [hep-ph]].
  * [21] J. I. Skullerud, J. Smit and A. Tranberg, JHEP 08, 045 (2003) doi:10.1088/1126-6708/2003/08/045 [arXiv:hep-ph/0307094 [hep-ph]].
  * [22] R. Micha and I. I. Tkachev, doi:10.1142/9789812704498_0020 [arXiv:hep-ph/0301249 [hep-ph]].
  * [23] J. Berges, S. Schlichting and D. Sexty, Phys. Rev. D 86, 074006 (2012) doi:10.1103/PhysRevD.86.074006 [arXiv:1203.4646 [hep-ph]].
  * [24] D. G. Figueroa, J. Garcia-Bellido and F. Torrenti, Phys. Rev. D 92, no.8, 083511 (2015) doi:10.1103/PhysRevD.92.083511 [arXiv:1504.04600 [astro-ph.CO]].
  * [25] K. Enqvist, S. Nurmi, S. Rusak and D. Weir, JCAP 02, 057 (2016) doi:10.1088/1475-7516/2016/02/057 [arXiv:1506.06895 [astro-ph.CO]].
  * [26] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso and D. Podolsky, JCAP 07, 006 (2006) doi:10.1088/1475-7516/2006/07/006 [arXiv:hep-ph/0602144 [hep-ph]].
  * [27] A. A. Abolhasani, H. Firouzjahi and M. M. Sheikh-Jabbari, Phys. Rev. D 81, 043524 (2010) doi:10.1103/PhysRevD.81.043524 [arXiv:0912.1021 [hep-th]].
  * [28] C. Cosme, D. G. Figueroa and N. Loayza, JCAP 05, 023 (2023) doi:10.1088/1475-7516/2023/05/023 [arXiv:2206.14721 [astro-ph.CO]].
  * [29] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, JCAP 04, 035 (2021).
  * [30] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, [arXiv:2102.01031 [astro-ph.CO]].
  * [31] P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 6175-6192 (1997).
  * [32] O. Lebedev, Prog. Part. Nucl. Phys. 120, 103881 (2021).
  * [33] L. Giusti, M. L. Paciello, C. Parrinello, S. Petrarca and B. Taglienti, Int. J. Mod. Phys. A 16, 3487-3534 (2001) doi:10.1142/S0217751X01004281 [arXiv:hep-lat/0104012 [hep-lat]].
  * [34] G. S. Bali, V. Bornyakov, M. Muller-Preussker and F. Pahl, Nucl. Phys. B Proc. Suppl. 42, 852-854 (1995) doi:10.1016/0920-5632(95)00401-T [arXiv:hep-lat/9412027 [hep-lat]].
  * [35] M. Schröck and H. Vogt, Comput. Phys. Commun. 184, 1907-1919 (2013) doi:10.1016/j.cpc.2013.03.021 [arXiv:1212.5221 [hep-lat]].
  * [36] J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, Phys. Rev. D 89, no.11, 114007 (2014) doi:10.1103/PhysRevD.89.114007 [arXiv:1311.3005 [hep-ph]].
  * [37] Adam Paszke et al., CoRR, vol. abs/1912.01703, 2019. [Online]. Available: http://arxiv.org/abs/1912.01703