# Propagation effects in resonant high-harmonic generation
and high-order frequency mixing in laser plasma

Vasily V. Strelkov<EMAIL_ADDRESS>Prokhorov General Physics Institute
of the Russian Academy of Sciences, Vavilova Street 38, Moscow 119991, Russia
Margarita A. Khokhlova  King’s College London, Strand Campus, London WC2R 2LS,
UK

###### Abstract

We study the phase matching of resonant high-harmonic generation (HHG) and
high-order frequency mixing (HFM) in plasma. We solve numerically the
propagation equations coupled with the time-dependent Schrödinger equation for
the nonlinear polarisation. The macroscopic harmonic signal is enhanced in the
vicinity of a multiphoton resonance with the transition between the ground and
autoionising states of the generating ion. We show that narrow and strong
resonances (as for gallium and indium ions) provide compensation of the plasma
dispersion in a spectral region above the exact resonance improving the phase
matching and leading to a high macroscopic signal. The compensation does not
take place for a wider resonance (as for manganese ions), instead the phase
matching is achieved in the HFM process. Comparing the XUV generated in
manganese plasma and in neon gas, we show that the resonant HHG in plasma is
an order of magnitude more effective than in the gas, moreover another order
of magnitude can be gained from the propagation using HFM in plasma.

††preprint: APS/123-QED

Attosecond pulses have proven to be a unique instrument to observe and control
ultrafast electron dynamics in different kinds of matter [1, 2, 3], and the
improvement of attopulse brightness is of key fundamental and technological
importance.

The production of attopulses largely relies on the process of high-harmonic
generation (HHG), where many photons of the strong laser driver are converted
into a comb of high harmonics (HH) through a highly nonlinear nonperturbative
response of the medium. The total efficiency of this process is composed of
the microscopic, single-particle, response and of the macroscopic response
describing the propagation of the light through the medium. Commonly, one of
these responses is improved to reach the higher efficiency of HHG, and thus
brighter attopulses: microscopic, by choosing the generating medium smartly,
or macroscopic, by optimising the propagation conditions such as the density
of the medium and the laser parameters.

One way to boost the microscopic HHG response is to use a resonance of HHs and
a transition between the states of the generating system. In the case of the
transition between the ground state and an autoionising state (AIS), the
resonant HHs are generated up to orders of magnitude more efficiently than the
nonresonant ones [4]. Resonant HHG enhancement has been observed in plasma
plumes [5, 6, 7, 8, 4], as well as for a ‘giant’ resonance in Xe [9].

For the HHG propagation, it has been shown [10, 11] that there are three
mechanisms limiting the coherent growth of the HH energy with propagation
length: absorption, phase matching and blue shift. However, there is other
fundamentally-close high-order process, taking place when the strong driver is
accompanied by another generating field — high-order frequency mixing (HFM)
[12, 13, 14] — which experiences the phase-matching and blue-shift limitations
to lesser extent [15, 16, 17, 18, 11, 19, 20].

In this Letter we study resonant HHG, as a process resulting in a boosted
microscopic response, considering theoretically the macroscopic properties of
this process. We also investigate the resonant HFM, optimizing the propagation
by choosing parameters of the additional weak generating field.

We start from the analysis of the propagation properties of the HHG process
along with the microscopic resonant response. For HHG the coherence length is
defined as

$L_{\mathrm{coh}}=\frac{\pi}{|\Delta k|}\,,$ (1)

where the wavevector mismatch is $\Delta k=k_{q}-qk_{0}$, here $q$ is a HH
order, $k_{0}$ and $k_{q}$ are the driver and HH wavevectors, correspondingly.
The wavevectors can be written as

$k_{0}=\frac{\omega_{0}}{c}\left(1+\Delta n_{0}^{\mathrm{(ion)}}+\Delta
n_{0}^{\mathrm{(el)}}\right)+\delta k_{0}^{\mathrm{(geom)}}\,,$ (2)
$k_{q}=\frac{q\omega_{0}}{c}\left(1+\Delta n_{q}^{\mathrm{(ion)}}+\Delta
n_{q}^{\mathrm{(el)}}\right)\,,$ (3)

where $\omega_{0}$ is the driver frequency, $\Delta n_{0}^{\mathrm{(el)}}$ is
the contribution to the refractive index of the driver due to the free
electrons written (atomic units are used throughout if not specified
otherwise) as

$\Delta n_{0}^{\mathrm{(el)}}=-\frac{2\pi N}{\omega_{0}^{2}}\,,$ (4)

with the electronic density $N$ (in usual units the latter equation is written
as $\Delta n_{0}^{\mathrm{(el)}}=-4.5\times
10^{-22}N[$cm${}^{-3}](\lambda_{0}[\mu$m$])^{2}$); $\Delta
n_{q}^{\mathrm{(ion)}}$ is the additive to the refractive index for the
harmonic due to the ions. Note that the common assumption for the HHG phase
matching is that the HH refractive index is unity. While it is a very reliable
assumption for nonresonant HHG, here we show that this is not the case for the
resonant one. Namely, the resonant term for $\Delta n_{q}^{\mathrm{(ion)}}$ is
written as

$\Delta n_{q}^{\mathrm{(ion)}}=-\frac{\pi
Nf}{\Omega}\quad\frac{\Delta\omega}{\Delta\omega^{2}+\Gamma^{2}/4}\,,$ (5)

where $\Omega$ and $\Gamma$ are the resonant frequency and width,
$\Delta\omega=q\omega_{0}-\Omega$ is the detuning from the resonance, $f$ is
the oscillator strength of the resonant transition, and we assume that the
ionic and the electronic densities are equal.

The other terms contributing to the medium refractive index in Eqs. (2) and
(3) are negligible: $\Delta n_{q}^{\mathrm{(el)}}$ is rather small because the
HH frequency is high, while $\Delta n_{0}^{\mathrm{(ion)}}$ is low because the
fundamental frequency is far from ionic resonances. We also neglect the
geometric dispersion $\delta k_{0}^{\mathrm{(geom)}}$, which can be taken into
account similarly to the $\Delta n_{0}^{\mathrm{(el)}}$ contribution if
necessary.

One can see from Eqs. (4) and (5) that above the resonance the contribution of
the ions to the detuning can compensate the contribution of the free
electrons. The maximum value of $|\Delta n_{q}^{\mathrm{(ion)}}|$ is achieved
near the resonance for $|\Delta\omega|=\Gamma/2$. For $\Delta\omega=\Gamma/2$:

$\Delta n_{q}^{\mathrm{(ion,max)}}/\Delta
n_{0}^{\mathrm{(el)}}=f\omega_{0}^{2}/(2\Gamma\Omega)\,.$ (6)

The absorption of the generated radiation near the resonance is also
important. The absorption length is $L_{\mathrm{abs}}=1/(N\sigma)$ with the
photoionisation cross-section $\sigma$, which can be written using the
parameters of the resonance (see, for instance, [4]) as

$L_{\mathrm{abs}}=\left[\frac{\pi
Nf}{c}\quad\frac{\Gamma}{\Delta\omega^{2}+\Gamma^{2}/4}\right]^{-1}\,.$ (7)

From Eqs. (7) and (1)-(5) for $\Delta\omega\ll\Omega$, we can write

$\frac{L_{\mathrm{abs}}}{L_{\mathrm{coh}}}=\left|\frac{2\Omega(\Delta\omega^{2}+\Gamma^{2}/4)}{\pi\Gamma\omega_{0}^{2}f}-\frac{\Delta\omega}{\pi\Gamma}\right|\,.$
(8)

In particular, near the resonance ($\Delta\omega=\Gamma/2$) we have

$\frac{L_{\mathrm{abs}}}{L_{\mathrm{coh}}}=\frac{1}{2\pi}\left|1-\frac{2\Gamma\Omega}{f\omega_{0}^{2}}\right|=\frac{1}{2\pi}\left|1-\frac{\Delta
n_{0}^{\mathrm{(el)}}}{\Delta n_{q}^{\mathrm{(ion,max)}}}\right|\,.$ (9)

From the latter equation we see that when the electronic contribution to the
dispersion dominates, HHG is phase-matching limited ($L_{\mathrm{coh}}\ll
L_{\mathrm{abs}}$); this is the case for high $\Gamma$, $\Omega$ and low $f$,
$\omega_{0}$. For the opposite case, free-electron and ion contributions can
compensate each other (such compensation was discussed in [21]), so the
macroscopic HHG is absorption-limited.

We study several resonant HHG conditions. We consider the transition
$3d^{10}4s^{2}\quad{}^{1}S_{0}\xrightarrow{}3d^{9}4s^{2}4p\quad^{1}P_{1}$ in
GaII (the transition energy is 21.88 eV, see [22, 23, 4] for other
parameters), the transition
$4d^{10}5s^{2}\quad{}^{1}S_{0}\xrightarrow{}4d^{9}5s^{2}5p(^{2}D)\quad^{1}P_{1}$
in InII (the transition energy is 19.92 eV, see [24] for other parameters) and
the ‘giant’ $3p\xrightarrow{}3d$ resonance in MnII (the transition energy is
51 eV, see [25] for other parameters).

We find that for the IR driver for all three cases the ratio (6) is below
unity, thus HHG is phase-matching limited. However, when we switch to a
visible driver, the situation becomes complicated. For the resonant H7
generation in GaII (driver wavelength is $\sim$ 400 nm111It has been shown
that HHG with 400 nm driver in plasma plumes [4, 8] achieves a much higher
efficiency than using an 800 nm one., which is the second harmonic of the
Ti:Sapph laser), we find that ratio (6) is 0.96, and for the resonant H7
generation in InII (driver wavelength is $\sim$ 435 nm), this ratio is as high
as 2.1. Thus, above the exact resonance the contribution of the free electrons
can be compensated by the ionic contribution in Ga as well as in In plasma. In
the vicinity of the frequencies, for which the compensation takes place, the
HHG is absorption-limited, while outside this bandgap it is phase-matching
limited, see Eq. (8). For the resonant H17 generation in Mn plasma (driver
wavelength is $\sim$ 420 nm) the ratio (6) is 0.034, thus the HHG is phase-
matching limited even using a visible generating field.

The wavevector mismatch, which is a key factor limiting the HHG macroscopic
efficiency in this case, can be essentially reduced for the process of HFM.
Let us consider the HFM process involving $Q$ quanta from the strong
fundamental and $m$ quanta from the additional weak field. The HFM process
leads to generation of frequency components

$\omega_{Q,m}=Q\omega_{0}+m\omega_{1}\,,$ (10)

where $Q$, $m$ and $Q+m$ are natural, integer and odd numbers, respectively.
For a weak additional field in comparison to the driver, HFM leads to a
notable efficiency only for processes involving few weak-field quanta: $|m|\ll
Q$. If the plasma contribution dominates in the dispersion, the phase mismatch
neglects [15, 16, 17, 18, 11, 19] for this process for $m<0$ and for the
frequency of the weak field

$\omega_{1}=\frac{|m|}{Q}\omega_{0}\,,$ (11)

which means that the weak-field frequency is much lower than the one of the
driver. Note that for orders $q$ that neighbour the ideally phase-matched $Q$,
the phase mismatch is rather small.

One more phenomenon that influences the macroscopic efficiency of HHG and HFM
is the plasma-induced blue shift. The ionisation of the medium by the
generating fields leads to frequency shifts of these fields. In spite of being
small in comparison to the fields’ frequencies, these shifts are multiplied by
a high order of the process that leads to a noticeable shift of the generated
frequencies. The latter shift can limit the macroscopic efficiency of the HHG
and HFM processes [11, 27]. It was shown [20] that for HFM the compensation of
the plasma-induced blue shifts is similar to the compensation of the phase
mismatch. Thus, the HFM process described by (10) and (11) for $q\approx Q$
almost does not suffer from the blue shift of the generating fields.

Below we study numerically the HHG and HFM macroscopic signals taking into
account the propagation of both the generating and generated fields in Ga and
Mn plasma. We assume that initially the plasma consists of singly-charged ions
and free electrons. The 1D propagation equations for the fields are integrated
as described in [27, 28] using single-ion response calculated via numerical
solution of the 3D time-dependent Schrödinger equation (TDSE) at every
propagation step. The TDSE is solved for the ionic model potential (similar to
the one used in [29, 4]) with the parameters chosen to reproduce the
properties of the real ion: the ionisation energy, the AIS energy and its
width.

To understand the propagation of the XUV field in plasma, we numerically
calculate the refractive index and the absorption coefficient in absence of
the laser field, that defines the actual HH refraction and absorption only
approximately because the laser field modifies them.

First, we study the effect of the detuning from the exact resonance between
the harmonic energy and the ground-AIS transition energy on the propagation of
the resonant HHG emission. Fig. 1 shows the calculated H7 intensity in Ga
plasma for three driving wavelengths in the vicinity of the 7-photon resonance
with the ground-AIS transition. When the driving wavelength is such that H7 is
in the exact resonance, as well as when H7 is above the exact resonance, the
harmonic signal monotonically saturates with the propagation length. This
means that the generation is absorption-limited in agreement with our
analytical model above. When the fundamental frequency is such that ‘resonant’
H7 is below the resonance, the initial growth of the macroscopic signal is
stronger, but it stops growing at a much shorter distance. This behaviour is
followed by oscillations of the harmonic signal, which is typical for phase-
matching limited generation. Thus, even a minor variation of the driving
wavelength (resulting in the detuning of H7 near the resonance) changes the
behaviour of its macroscopic signal.

Figure 1: H7 macroscopic signals generated in Ga plasma as a function of
propagation length. The drivers’ wavelengths are 396 nm (above resonance), 400
nm (in resonance), 404 nm (below resonance). The laser pulse intensity is
$3\times 10^{14}$ W/cm2, its duration is 60 fs, and the plasma density is
$10^{18}$cm-3.

Figure 2: XUV absorption line and intensities of microscopic and macroscopic
(a) H7 responses of Ga and (b) H17 responses of Mn plasma for different
fundamental wavelengths. Other laser pulse parameters for (a) are the same as
in Fig. 1, and for (b) the driving pulse intensity is $9\times 10^{14}$ W/cm2,
its duration is 40 fs. The upper panels at each graph show two main additives
to the refractive indexes of the driver and of the resonant HH, defining phase
matching of the HHG.

To understand this feature in more detail, we perform the calculations for a
number of driving wavelengths, see the results in Fig. 2. The macroscopic
signal presented in the figure is the maximal HH macroscopic intensity
achieved during propagation in the target. We see that the absorption line and
the resonant enhancement of the microscopic HHG response do not coincide. For
Ga plasma (Fig. 2a), at first glance, the generation conditions below the
absorption line are favourable due to the high microscopic response and low
absorption. However, Fig. 1 as well as the macroscopic response in Fig. 2a
show higher signal above the resonance.

To clarify this, Fig. 2 shows also the numerically calculated two main
contributions to the refractive indexes, which define the wavevector mismatch
for HHG in plasma, see Eqs. (2) and (3). One can see that above the resonance
this two terms compensate each other improving phase matching for HHG. As
discussed above, such compensation is possible for strong and narrow
resonances (such as ones in GaII and InII), and only using visible or UV
fundamental for which the contribution $\Delta n_{0}^{(el)}$ is moderate.

When H7 in Ga plasma is below the resonance, the microscopic response is high,
but the phase matching is weak and HHG is limited by both absorption and phase
matching, providing the weakest signal shown in Fig. 1. Above the resonance
the microscopic response is weaker but the phase matching is better, providing
higher HHG signal. The XUV absorption plays a key role in the saturation of
this signal, as well as for the resonant case.

For the nonresonant harmonics, the plasma contribution to dispersion
dominates, leading to a poor phase matching for their generation. Thus, the
resonant harmonic is enhanced with respect to the neighboring ones not only
due to microscopic, but also due to macroscopic response properties. For the
resonant HHG in Ga plasma [4] the observed enhancement of the macroscopic
signal was higher than the calculated in the same paper enhancement of the
microscopic signal. This can be attributed to the described above favorable
phase matching of the resonant HHG.

We also study HHG in Mn plasma in the vicinity of the resonance. Fig. 2b
presents the absorption coefficient as well as the microscopic and macroscopic
responses. One can see that the peaks of the three lines coincide, but their
widths are different. The figure also shows the main contributions to the
refraction coefficients for H17 and for the driver. We see that the phase-
matching conditions qualitatively differ from the ones in Ga plasma: in Mn
plasma the resonant refraction coefficient is too low to compensate the plasma
dispersion (in agreement with our analytical model), so HHG is phase-matching
limited in the resonant case, as well as below and above it. Fig. 3 shows the
resonant H17 intensity with propagation length. One can see deep oscillations
of the macroscopic signal, which are typical for phase-matching-limited
generation.

To characterise quantitatively the HHG efficiency in plasma, we calculate the
HHG macroscopic signal in (initially neutral) Ne. Note that a rather high
driving intensity is necessary to achieve a cutoff higher than H17 of the 400
nm fundamental. We choose Ne because HHG in Ne is compatible with such a high
fundamental intensity. In Fig. 3 one can see that the H17 generation
efficiency in Mn plasma can be an order of magnitude higher than in Ne. This
is achieved due to the very high resonant microscopic response of MnII and in
spite of the weak phase matching of HHG in Mn plasma.

Figure 3: Macroscopic HHG and HFM signals in Mn plasma with propagation. The
XUV intensities are divided by the maximal XUV intensity from HHG in Ne. The
driver wavelength for HHG is 420 nm and for HFM is 415 nm, the weak-field
wavelength is 3800 nm. The driving pulse intensity and duration are as in Fig.
2b, the weak-field intensity is 0.5% of the driver. The inset shows the sketch
of the HHG and HFM processes.

As we discussed above, the efficiency of phase-matching-limited HHG can be
improved by realising HFM. The phase mismatch vanishes for harmonics with
orders $q\approx Q$, see Eq. (11). The highest conversion efficiency is
achieved for $q$ slightly lower than $Q$ [27]; so to optimise HFM with $q=17$,
we choose $Q=19$ and $m=-2$ and from Eq. (11) we find that the low-frequency
field wavelength is 3800 nm. In Fig. 3 we see that the improved phase matching
for HFM allows a growth of the signal up to a much longer propagation length
than for HHG. This leads to a higher macroscopic response which exceeds the
one for HHG in Ne by two orders of magnitude.

To summarise, we study theoretically the phase matching of resonant HHG and
HFM in plasma. We integrate the propagation equations using the nonlinear
polarisation found via numerical solution of the TDSE at every propagation
step. The macroscopic HH signal is enhanced in the vicinity of multiphoton
resonances with the transition between the ground state and AIS of the
generating ion. We show that narrow and strong resonances in Ga II (same as in
In II) ions provide compensation of the plasma dispersion in the region above
the exact resonance. As a result, the improved phase matching leads to a high
macroscopic signal in this spectral region. However, the compensation does not
take place for wider resonances, as in Mn plasma, where the phase matching can
be achieved instead in the HFM process. We compare the HHG efficiency in Ne
gas and Mn plasma and show that the resonant HHG in plasma is an order of
magnitude more effective than in the gas, while another order of magnitude can
be gained switching from HHG to HFM in plasma.

Funding: This study was funded by RSF (Grant No. 24-22-00108).

Acknowledgments: MK acknowledges Royal Society funding under URF\R1\231460.

## References

  * Villeneuve [2018] D. M. Villeneuve, Attosecond science, Contemp. Phys. 59, 47 (2018).
  * Ryabikin  _et al._ [2023] M. Y. Ryabikin, M. Y. Emelin, and V. V. Strelkov, Attosecond electromagnetic pulses: generation, measurement, and application. attosecond metrology and spectroscopy, Phys. Usp. 66, 360 (2023).
  * Khokhlova  _et al._ [2023] M. Khokhlova, E. Pisanty, and A. Zaïr, Shining the shortest flashes of light on the secret life of electrons, in  _Advanced Photonics, Vol. 5, Issue 6_ (2023) p. 060501.
  * Singh  _et al._ [2021] M. Singh, M. A. Fareed, V. Strelkov, A. N. Grum-Grzhimailo, A. Magunov, A. Laramée, F. Légaré, and T. Ozaki, Intense quasi-monochromatic resonant harmonic generation in the multiphoton ionization regime, Optica 8, 1122 (2021).
  * Singh  _et al._ [2024] M. Singh, M. A. Fareed, R. G. Shirinabadi, R. Marcelino, F. Zhu, F. Légaré, and T. Ozaki, Recent advances in high-order harmonic generation from laser-ablated plumes at the advanced laser light source laboratory, Fundamental Plasma Physics , 100043 (2024).
  * Ganeev  _et al._ [2016] R. A. Ganeev, Z. Wang, P. Lan, P. Lu, M. Suzuki, and H. Kuroda, Indium plasma in single- and two-color mid-infrared fields: Enhancement of tunable harmonics, Phys. Rev. A 93, 043848 (2016).
  * Fu  _et al._ [2022] W. Fu, Y. H. Lai, J. Liang, and W. Li, Investigation of high-harmonic cutoff of metal ions driven by near-infrared laser, Opt. Express 30, 23090 (2022).
  * Fu  _et al._ [2023] W. Fu, J. Wang, J. Yu, and W. Li, Extension of high-order harmonic generation cutoff from laser-ablated tin plasma plumes, Opt. Express 31, 15553 (2023).
  * Shiner  _et al._ [2011] A. D. Shiner, B. E. Schmidt, C. Trallero-Herrero, H. J. Wörner, S. Patchkovskii, P. B. Corkum, J.-C. Kieffer, F. Legare, and D. M. Villeneuve, Probing collective multi-electron dynamics in xenon with high-harmonic spectroscopy, Nat. Phys. 7, 464 (2011).
  * Constant  _et al._ [1999] E. Constant, D. Garzella, P. Breger, E. Mével, C. Dorrer, C. Le Blanc, F. Salin, and P. Agostini, Optimizing high harmonic generation in absorbing gases: Model and experiment, Phys. Rev. Lett. 82, 1668 (1999).
  * Khokhlova and Strelkov [2020] M. A. Khokhlova and V. V. Strelkov, Highly efficient XUV generation via high-order frequency mixing, New J. Phys. 22, 093030 (2020).
  * Eichmann  _et al._ [1995] H. Eichmann, A. Egbert, S. Nolte, C. Momma, B. Wellegehausen, W. Becker, S. Long, and J. K. McIver, Polarization-dependent high-order two-color mixing, Phys. Rev. A 51, R3414 (1995).
  * Gaarde  _et al._ [1996] M. B. Gaarde, A. L’Huillier, and M. Lewenstein, Theory of high-order sum and difference frequency mixing in a strong bichromatic laser field, Phys. Rev. A 54, 4236 (1996).
  * Balcou  _et al._ [1999] P. Balcou, A. S. Dederichs, M. B. Gaarde, and A. L’Huillier, Quantum-path analysis and phase matching of high-order harmonic generation and high-order frequency mixing processes in strong laser fields, J. Phys. B: At. Mol. Opt. Phys. 32, 2973 (1999).
  * Shkolnikov  _et al._ [1993] P. L. Shkolnikov, A. E. Kaplan, and A. Lago, Phase matching for large-scale frequency upconversion in plasma, Opt. Lett. 18, 1700 (1993).
  * Milchberg  _et al._ [1995] H. M. Milchberg, C. G. Durfee III, and T. J. McIlrath, High-order frequency conversion in the plasma waveguide, Phys. Rev. Lett. 75, 2494 (1995).
  * Meyer  _et al._ [1996] S. Meyer, H. Eichmann, T. Menzel, S. Nolte, B. Wellegehausen, B. N. Chichkov, and C. Momma, Phase-matched high-order difference-frequency mixing in plasmas, Phys. Rev. Lett. 76, 3336 (1996).
  * Birulin  _et al._ [1996] A. V. Birulin, V. T. Platonenko, and V. V. Strelkov, High-harmonic generation in interfering waves, J. Exp. Theor. Phys. 83, 33 (1996).
  * Cohen  _et al._ [2007] O. Cohen, T. Popmintchev, D. M. Gaudiosi, M. M. Murnane, and H. C. Kapteyn, Unified microscopic-macroscopic formulation of high-order difference-frequency mixing in plasmas, Phys. Rev. Lett. 98, 043903 (2007).
  * Hort  _et al._ [2021] O. Hort, A. Dubrouil, M. A. Khokhlova, D. Descamps, S. Petit, F. Burgy, E. Mével, E. Constant, and V. V. Strelkov, High-order parametric generation of coherent XUV radiation, Opt. Express 29, 5982 (2021).
  * Ganeev  _et al._ [2013] R. A. Ganeev, T. Witting, C. Hutchison, V. V. Strelkov, F. Frank, M. Castillejo, I. Lopez-Quintas, Z. Abdelrahman, J. W. G. Tisch, and J. P. Marangos, Comparative studies of resonance enhancement of harmonic radiation in indium plasma using multicycle and few-cycle pulses, Phys. Rev. A 88, 033838 (2013).
  * Pindzola  _et al._ [1982] M. S. Pindzola, D. C. Griffin, and C. Bottcher, Electron-impact excitation autoionization of ga ii, Phys. Rev. A 25, 211 (1982).
  * Peart  _et al._ [1987] B. Peart, I. C. Lyon, and K. Dolder, Measurements of absolute photoionisation cross sections of ga+ and zn+ ions, Journal of Physics B: Atomic and Molecular Physics 20, 5403 (1987).
  * Duffy and Dunne [2001] G. Duffy and P. Dunne, The photoabsorption spectrum of an indium laser produced plasma, Journal of Physics B: Atomic, Molecular and Optical Physics 34, L173 (2001).
  * Kjeldsen  _et al._ [2004] H. Kjeldsen, F. Folkmann, B. Kristensen, J. B. West, and J. E. Hansen, Absolute cross section for photoionization of mn+ in the 3p region, Journal of Physics B: Atomic, Molecular and Optical Physics 37, 1321 (2004).
  * Note [1] It has been shown that HHG with 400 nm driver in plasma plumes [4, 8] achieves a much higher efficiency than using an 800 nm one.
  * Birulia  _et al._ [2024] V. A. Birulia, M. A. Khokhlova, and V. V. Strelkov, Macroscopic effects in generation of attosecond xuv pulses via high-order frequency mixing in gases and plasma, New Journal of Physics 26, 023005 (2024).
  * Strelkov and Khokhlova [2023] V. V. Strelkov and M. A. Khokhlova, Phase-matching gating for isolated attosecond pulse generation, Photonics 10, 10.3390/photonics10101122 (2023).
  * Strelkov [2010] V. Strelkov, Role of autoionizing state in resonant high-order harmonic generation and attosecond pulse production, Phys. Rev. Lett. 104, 123901 (2010).


# Inference with Hamiltonian Sequential Monte Carlo Simulators

Remi Daviet Daviet: Wharton Marketing Department, University of Pennsylvania

###### Abstract

The paper proposes a new Monte-Carlo simulator combining the advantages of
Sequential Monte Carlo simulators and Hamiltonian Monte Carlo simulators. The
result is a method that is robust to multimodality and complex shapes to use
for inference in presence of difficult likelihoods or target functions.
Several examples are provided.

Keywords: Sequential Monte Carlo, Hamiltonian Monte Carlo, Laplace Estimators

##  1 Introduction

Econometrics has traditionally relied on methods based on the optimization of
a target function to perform inference. Such methods include the widespread
Maximum Likelihood Estimators (MLE), Generalized Method of Moments (GMM) or
the Least Squares Estimators (LSE). These methods are especially effective
when the target function is unimodal and is shown to converge rapidly to a
point.

Recent developments in economic modeling have lead to the emergence of more
complicated target functions where multimodality, complex shapes and slow
convergence make traditional inference approaches ineffective. Numerical
maximization methods such as Nelder-Mead, Newton-Raphson or Expectation
Maximization do not guarantee the convergence toward a global maximum
[Gourieroux and Monfort, 1995]. In addition, these methods provide a point-
estimate and do not give information about the shape of the target function.

These more complicated target functions are found in almost every subfield of
economics: Multimodal likelihoods for instance can be the result of using DSGE
models in Macroeconomics [Herbst and Schorfheide, 2014], GARCH models in
Finance [Doornik et al., 2000], BLP models in industrial organization [Bajari,
2003], Spatial linear models in Urban Economics [Mardia and Watkins, 1989]. In
fact, with many non-linear models, likelihood functions are non-smooth and
multimodal [Koop and Potter, 1999]. A similar observation can be made for
models with structural breaks and outliers Koop and Potter [2000]. Various
attempts have been made to mitigate the problem, generally relying on Bayesian
or Quasi-Bayesian Methods using Markov-Chain Monte Carlo (MCMC) or Sequential
Monte Carlo (SMC) simulations.

MCMC simulation has the advantage of being able to recover more information
about a target function than optimization algorithms [Chernozhukov and Hong,
2003]. The only restriction is that the target function has to be positive. In
the case of a non-positive target function, the exponential of the function
can be taken. This exponentiated function can be treated as a quasi-
likelihood. The draws from the simulation provide a Monte Carlo approximation
of the distribution of interest, facilitating inference. Application of this
method to economics can be found in several papers across all subfields. For
instance, Herbst has proposed several Monte Carlo methods to solve the
estimation problems with Macroeconomic DSGE models [Herbst and Schorfheide,
2014]. We find similar approaches in industrial organization for demand
estimation models [Jiang et al., 2009], or in time series analysis [Burda,
2015].

The key problem with traditional MCMC simulators using Metropolis-Hastings or
Gibbs Sampling methods is that they do not behave necessarily well under
multimodality, concentrated mass or complex shapes. The Markov Chain used to
simulate the distribution can get trapped in one of the modes if it is not
close enough to the other modes and there is no probability mass between them.
When facing complex concentrated shapes, these simulators are not performing
well exploring the space of interest due to a high level of rejections when
trying to move away from the current point in the chain. Some specific
simulators have shown interesting properties concerning these problems. The
families of Population Monte Carlo (PMC) and more generally SMC simulators
solve the problem of multimodality by not having a unique Markov Chain but
rather a large set of particles exploring the parameter space [Cappe et al.,
2004, Del Moral et al., 2006, Durham and Geweke, 2013]. The sequential
approach of the SMC also partly solves the problem of concentrated mass by
allowing the target function to progressively converge toward its final form
[Chopin, 2002]. The MCMC simulators using Hamiltonian dynamics, sometimes
refereed to as Hamiltonian Monte Carlo (HMC) simulators, have shown to be very
effective in exploring the parameter space when the target function has a
complex elongated shape or isolated concentrated mass.

Our contribution is to provide a method that is robust to most types of
multimodality and complex shapes by combining the advantages of the SMC and
the HMC. Moreover, we implement a kernel based resampling method to improve
robustness and efficiency. While most econometrics papers describing MCMC or
SMC methods are put in a Bayesian framework, we keep the description in a
general framework. Using our method, practitioners should be able to perform
inference without having to worry about the shape or complexity of their
target function.

The advantages of recovering the full shape of the target function are
multiple. Multimodality can easily be identified, concentration can be
measured and counter-factual checking or prediction can be done by integration
over the parameter space. Moreover, the Monte Carlo approach allows for easy
variable transformation without the need of derivation of a complicated
Jacobian.

Our method can also be used for maximization of a complicated function $f(x)$
by simulated annealing: the target quasi-likelihood is of the form
$\exp(f(x))^{\gamma}$, where $\gamma$ is chosen by the practitioner. The
quasi-likelihood will become concentrated on the maximizers of $f(x)$ as
$\gamma\to\infty$ [Hwang, 1980]. The use of SMC methods for optimization via
simulated annealing has already been detailed and proven effective [Zhou and
Chen, 2013].

We will first present a brief review of the various uses of statistical
simulation. We will then describe three common simulations methods
(Metropolis-Hastings, HMC, and SMC). We will compare their performance using
simple examples. We will finally introduce our HSMC method, its properties and
present several simple examples as applications.

##  2 Statistical simulation

Statistical simulation generates samples $\\{\theta_{n}\\}_{n=1}^{N}$
following a distribution with density $f(\theta_{n})$ to get an approximation
of various quantities depending on that distribution. These quantities are
usually moments, modes or quantiles.

To compute population moments, we can use the sample moments which converge
under a Law of Large Number to the desired quantity under weak regularity
assumptions. For some continuous function $h(\theta_{n})$ and a Euclidean
vector $\theta_{n}$, we have:

$\frac{1}{N}\sum_{n}h(\theta_{n})\overset{p}{\to}\int
h(\theta)f(\theta)d\theta$

There are several ways to test for the existence of several modes, count them
and find their location. A review of standard kernel based methods can be
found in Vieu [1996]. Most popular methods based on critical kernel bandwidth
are presented in Silverman [1981] and Minnotte [1997]. Another popular
approach called excess mass approach can be found in several other papers
[Müller and Sawitzki, 1991, Butucea et al., 2007]. The goal of this paper is
not to treat mode estimation methods in details and we will leave further
readings to the interested reader.

Transformations of the target function can be computed directly by evaluating
the transformations of the resulting draws, without the need for analytical
derivation or approximation of the Jacobian term which is in many cases very
complicated.

##  3 Common simulation methods

As it is often impossible to draw directly from a distribution, several
methods have been developed to draw a sample that approximately follows the
target distribution. The researcher only needs to be able to evaluate a kernel
of the target distribution density function.

The Metropolis-Hastings (MH) algorithm and the HMC are both Markov Chain Monte
Carlo (MCMC) methods. Consider a set of particles $\\{\theta_{t}\\}$ whose
distribution seeks to approximate an underlying target function of interest
$f(x)$. The generic principle of an MH or HMC algorithm is as follows. A
particle $\theta_{t}$ randomly moves to a new position $\theta^{\prime}_{t}$
in the space $\Theta$. After each move, with some probability the new position
is accepted and $\theta_{t+1}=\theta^{\prime}_{t}$. If the move is rejected,
$\theta_{t+1}=\theta_{t}$. The passage from $\theta_{t}$ to $\theta_{t+1}$ is
called an MCMC step. After performing a series of initial steps called burn-in
period the Markov Chain should reach its equilibrium distribution which
corresponds to the target distribution. We can then use the history of the
positions of $\theta_{t}$ as a sample approximating the target distribution.
The difference between MH and HMC is in the way the random move is done, and
the formulas to compute acceptance probabilities at the end of the move.

With SMC methods, $N$ particles $\\{\theta_{nt}\\}_{n=1}^{N}$ are moved
simultaneously using MCMC steps and resampled at each iteration. These methods
have the advantage of being able to recover the shape of multimodal
distributions with separate areas of concentrated mass. These distributions
usually cannot be recovered properly using standard MCMC methods as the Markov
Chains get trapped under one of the modes without being able to move to the
other mode if no probability mass connects them.

###  3.1 Metropolis-Hastings

The Metropolis-Hastings algorithm [Metropolis et al., 1953, Hastings, 1970] is
the most common accept-reject MCMC method to approximate a target density
proportional to $f(\theta)$. Starting from an arbitrary point $\theta_{0}$,
for each iteration a new particle is drawn from a distribution with density
$q(\theta^{\prime}|\theta_{t})$. For instance, this distribution could be a
$\mathcal{N}(\theta_{t},\Sigma)$. Then, we set $\theta_{t+1}=\theta^{\prime}$
with probability $\min(1,\frac{f(\theta^{\prime})\cdot
q(\theta_{t}|\theta^{\prime})}{f(\theta_{t})\cdot
q(\theta^{\prime}|\theta_{t})})$, and $\theta_{t+1}=\theta_{t}$ otherwise.

After an initial burn-in period of $b$, the values
$\\{\theta_{t}\\}_{t=b}^{T}$ should be approximately distributed folowing the
normalized target density.

###  3.2 Hamiltonian Monte Carlo

As with MH, the purpose of the HMC method is to formulate a Markov chain for
which, under certain conditions, a multiple of $f(\theta_{t})$ is the density
of the stationary distribution. It relies on Hamiltonian dynamics to move a
particle $\theta_{t}$ to a new point $\theta^{\prime}$ in the $\Theta$ space.
This new point, called proposal will then be accepted or rejected as the new
value for $\theta_{t+1}$ following a method similar to Metropolis-Hastings
rejection step. The movement in the $\Theta$ space can be constrained and we
then refer to the method as a Constrained Hamiltonian Monte Carlo (CHMC).

To describe the HMC step, we need to define the function $U(\theta_{t})$ with
gradient denoted $\nabla U(\theta_{t})$:

$U(\theta_{t})=-\log(f(\theta_{t}))$

We also need to define an auxiliary vector $p$ of dimension $\dim(\theta_{t})$
that will represent the momentum of the particle when moving following
Hamiltonian dynamics. Intuitively, the particle will move on a surface where
the potential energy depending on the altitude is represented by
$U_{t}(\theta_{n})$. When the particle goes up a slope, it will slow down or
even turn back. The continuous movement is approximated by a series of $L$
steps of size $\epsilon$. At the end of the last step, the momentum is
reversed to make the proposal symmetric. If we start at the final position
with the final reversed momentum, we will find the particle going back to the
original position after $L$ steps. This ensures reversibility and facilitate
the computation of the acceptance probability.

The HMC step proceeds as follow:

  1. 1.

define the starting position of the proposal
$\theta^{\prime}=\hat{\theta}_{t}$

  2. 2.

draw an initial momentum vector $p$ from a multivariate Gaussian distribution:
$p\sim\mathcal{N}(0,M)$

  3. 3.

update the momentum vector by half a step taking the gradient into account:
$p^{\prime}=p-\frac{\epsilon}{2}\cdot\nabla U(\theta^{\prime})$

  4. 4.

Repeat for $l=1,\ldots,L$

    1. (a)

update the position by a full step:
$\theta^{\prime}=\theta^{\prime}+\epsilon\cdot p^{\prime}$

    2. (b)

update the momentum by a full step, except at the end of the trajectory: if
$(l\neq L)$, then $p^{\prime}=p^{\prime}-\epsilon\cdot\nabla
U(\theta^{\prime})$

  5. 5.

update the momentum vector by half a step:
$p^{\prime}=p^{\prime}-\frac{\epsilon}{2}\cdot\nabla U(\theta^{\prime})$

  6. 6.

negate the momentum vector: $p^{\prime}=-p^{\prime}$

  7. 7.

compute the acceptance probability:
$a=\min\left(1,\exp\left[U(\theta_{t})-U(\theta^{\prime})+\frac{\sum
p^{2}}{2}-\frac{\sum{p^{\prime}}^{2}}{2}\right]\right)$

  8. 8.

set $\theta_{t+1}=\theta^{\prime}$ with probability $a$, and
$\theta_{t+1}=\theta_{t}$ otherwise

The HMC mutation step needs to be tuned by choosing appropriately the
quantities $(M,L,\epsilon)$. To learn more about how to choose these
quantities, see the review paper by Neal [2011].

When we want to put constraints on some of the dimensions of $\theta$, we can
modify the HSMC method to have the particles to bounce off the constraints as
if they were walls. Most of the time, the constraints are of the form
$\theta_{dt}\leq u_{d}$ or $\theta_{dt}\geq l_{d}$ for some dimensions $d$ of
$\theta_{t}$. The position updating step 4a is then replaced by:

  1. 1.

for each dimension $d$ of $\theta^{\prime}$:

    1. (a)

update position in dimension $d$:
$\theta^{\prime}_{d}=\theta^{\prime}_{d}+\epsilon\cdot p^{\prime}_{d}$

    2. (b)

if $\theta^{\prime}_{d}$ is constrained, repeat the following until
$\theta^{\prime}_{d}$ satisfies all constraints:

      1. i.

if $(\theta^{\prime}_{d}>u_{d})$, then
$\theta^{\prime}_{d}=u_{d}-(\theta^{\prime}_{d}-u_{d})$ and
$p^{\prime}_{d}=-p^{\prime}_{d}$

      2. ii.

if $(\theta^{\prime}_{d}<l_{d})$, then
$\theta^{\prime}_{d}=l_{d}+(l_{d}-\theta^{\prime}_{d})$ and
$p^{\prime}_{d}=-p^{\prime}_{d}$

With this modified approach, if the particle passes a constraint ”wall” during
a position update, the symmetric of the particle $\theta^{\prime}$ relative to
the wall in the $\Theta$ space is taken, and the momentum in the constrained
dimension is reversed. This approach simulates the particle bouncing off the
wall and preserves reversibility. A similar approach can be used for more
complex constraints of the form $G(\theta)\geq 0$.

###  3.3 Sequential Monte Carlo

As mentioned before, SMC methods use multiple particles moving in parallel
[Del Moral et al., 2006]. While there are many different types of SMC, we are
going to describe one of the most popular approaches alternating Importance
Resampling and MH steps. The goal is to obtain a sample
$\\{\theta_{n}\\}_{n=1}^{N}$ from a sequence of distributions with densities
$f_{1}(\theta_{n}),\ldots,f_{T}(\theta_{n})$. The simulator we propose works
best when target densities in the sequence are smoother at the beginning and
progressively converging toward the final sharper target density. We also
require an initial distribution with density $f_{0}(\theta_{n})$ that is easy
to sample from and covers well the mass of the first distribution in the
sequence. An SMC simluator provides us for each step $t=1,\cdots,T$ with a set
of values $\\{\theta_{n}\\}_{n=1}^{N}$ that approximately follow the
distribution $f_{t}(\theta_{n})$.

A sequence of distributions can be found in many applications relevant to
economics. In frequentist econometrics, the sequence can be the likelihood or
quasi-likelihood for data collected until time $t$, e.g.
$f_{t}(\theta_{n})\propto L(\theta_{n};y_{1},\ldots,y_{t})$. The quasi-
likelihood can be built from any estimator maximizing or minimizing a target
function such as the Generalized Method of Moments or a Least Squares
Estimator [Chernozhukov and Hong, 2003]. The counterpart in the Bayesian
framework would be the posterior distribution of the parameter $\theta$ given
the data until time $T$. Compared to a standard MCMC approach which requires
an evaluation of the target function with every observation at each step, the
SMC approach is less computationally intensive. Moreover, adding the
observations one or a few at a time creates a desirable tempering effect
[Chopin, 2002]. This approach is particularly efficient in large datasets
where new observations come regularly as updating the estimator can be done in
one step.

A second application is kernel density estimation when observations are added
a few at a time and bandwidth is progressively shrunk. An application of this
approach can be found later in this paper.

Another application proposed by Neal [2001] shows the benefits of moving
progressively from a tractable distribution $f_{1}(\theta_{n})$ to a target
distribution $f(\theta_{n})$ by geometrically reweighting them :
$f_{t}(\theta_{n})\propto
f(\theta_{n})^{\phi_{t}}f_{1}(\theta_{n})^{1-\phi_{t}}$ with
$0\leq\phi_{1}<\ldots<\phi_{T}=1$.

Finally, SMC can be used for maximization using a simulated annealing
approach. Our sequence of distributions will then become
$f_{t}(\theta_{n})\propto f(\theta_{n})^{\gamma_{t}}$ with $\gamma_{t}$
increasing to high values as $t$ increases.

The algorithm is as follow:

  1. 1.

Initialization: Draw $N$ particles $\\{\theta_{n}^{(0)}\\}_{n=1}^{N}$ from
$f_{0}(\theta_{n})$

  2. 2.

Repeat for $t=1,\ldots,T$

    1. (a)

Correction: assign weight $w_{n}^{(t)}=f_{t}(\theta_{n})/f_{t-1}(\theta_{n})$
to each of the particles $\\{\theta_{n}^{(t-1)}\\}_{n=1}^{N}$

    2. (b)

Selection: draw $N$ new particles $\\{\hat{\theta}_{n}^{(t)}\\}_{n=1}^{N}$
with replacement from the current sample of particles using weights
$w_{n}^{(t)}$. Give the new particles a weight of $1$.

    3. (c)

Mutation: For each particle, perform a MH step as described in section 3.1 to
obtain a new sample of particles $\\{\theta_{n}^{(t)}\\}_{n=1}^{N}$.

###  3.4 Examples

We are now going to illustrate the advantages of various algorithms with
simple examples. The first example is a Rosenbrock’s banana function and
illustrates the problems classical MH face when the function has an elongated
shape. The second example is a 6-dimensional normal distribution. This example
is designed to show the difference in convergence speeds when dimensionality
increases.

####  3.4.1 Banana function

The Rosenbrock’s banana function (Figure 1) is defined as follow:

$f(\theta)\propto
e^{\frac{1}{8}\left(-5\left(y-x^{2}\right)^{2}-x^{2}\right)}$

Figure 1: Contour plot of the Rosenbrock function

Using the MH algorithm, we run two trials with different $\sigma$ for the
$\mathcal{N}(\theta,\sigma I)$ proposal. When $\sigma=0.2$, the acceptance
rate is $66.1\%$ but the chain fails to cover the distribution after 1000
iterations. When $\sigma=1$, the coverage improves but the acceptance rate
decreases to $39.6\%$. A visual representation of the paths is provided in
Figure 2.

Using the HMC approach, with $\epsilon=0.05$ and $L=20$ steps, the target
density is well covered and the acceptance rate is $99.8\%$. A visual
representation of the paths is provided in Figure 3.

Figure 2:  Metropolis-Hastings paths for a Rosenbrock function

Figure 3:  Hamiltonian Monte-Carlo path for a Rosenbrock function

####  3.4.2 6-dimensional normal distribution

Using a $N(\mu,I)$ as the target, where $\mu=(10,10,10,-10,-10,-10)$, we
compare the performance of the MH and HMC approach starting the chain at the
point $(-15,-15,-15,15,15,15)$. The $\sigma$ for the MH method is chosen to
get at least a $60\%$ acceptance rate.

In the case of the MH approach, the chain requires a burn in of about 800
iterations to reach the mass of the target density. With the HMC approach,
only 50 iterations are required for the chain to converge. A visual
representation of the paths is provided in Figure 4.

Figure 4: Comparison of Metropolis-Hastings path (Top) and Hamiltonian Monte-
Carlo path (Bottom) for a 6-dimensional normal distribution.

##  4 Hamiltonian Sequential Monte Carlo

In this section, we are going to describe our method and its properties. The
HSMC method replaces the MH step of standard SMC by a Hamiltonian step. The
re-sampling is also done using a leave-one-out approximation
$\hat{f}_{t-1}(\theta_{n})$ of the observed distribution of the particles
instead of the theoretical distribution $f_{t-1}(\theta_{n})$ as in specific
cases the particles do not have the time to converge to their stationary
distribution in one step. This method shows good convergence rates when
applied to importance resampling [Delyon et al., 2016].

###  4.1 Algorithm

  1. 1.

Initialization: Draw $N$ particles $\\{\theta_{n}^{(0)}\\}_{n=1}^{N}$ from
$f_{0}(\theta_{n})$

  2. 2.

Repeat for $t=1,\ldots,T$

    1. (a)

Correction: assign weight
$w_{n}^{(t)}=f_{t}(\theta_{n})/\hat{f}_{t-1}(\theta_{n})$ to each of the
particles $\\{\theta_{n}^{(t-1)}\\}_{n=1}^{N}$, where
$\hat{f}_{t-1}(\theta_{n})$ is a ”leave-one-out” kernel density estimate.

    2. (b)

Selection: draw $N$ new particles $\\{\hat{\theta}_{n}^{(t)}\\}_{n=1}^{N}$
with replacement from the current sample of particles using weights
$w_{n}^{(t)}$. Give the new particles a weight of $1$.

    3. (c)

Mutation: For each particle, perform a Hamiltonian step as described in
section 3.2 to obtain a new sample of particles
$\\{\theta_{n}^{(t)}\\}_{n=1}^{N}$.

In the initialization phase we obtain a sample of particles distributed
according to $f_{0}(\theta_{n})$ distribution. For the HSMC method to perform
well we need a distribution that covers well the whole $\Theta$ space and has
mass where the other distributions $f_{t}(\theta_{n})$ in the sequence have
mass too.

In the loop, before the correction phase, we have particles all weighted to
$1$ that provides a Monte-Carlo simulation of the distribution
$f_{t-1}(\theta_{n})$. In the correction phase, we reweight them to obtain an
approximation of the distribution $f_{t}(\theta_{n})$ by importance sampling.

In the selection phase, we perform sampling importance resampling to obtain an
approximation of the distribution $f_{t}(\theta_{n})$ using particles
$\\{\hat{\theta}_{n}^{(t)}\\}_{n=1}^{N}$ with equal weights. Note that at the
end of the selection phase, we can expect to have several particles sharing
the same value.

Finally, in the mutation phase, we explore the $\Theta$ space by moving the
particles using a Hamiltonian Monte Carlo (HMC) approach. Since
$f_{t}(\theta_{n})$ is the stationary distribution of our HMC, the particles
both before and after the HMC step should approximate $f_{t}(\theta_{n})$.
However, in the case $\\{\hat{\theta}_{n}^{(t)}\\}_{n=1}^{N}$ do not follow
exactly the distribution $f_{t}(\theta_{n})$, performing a HMC step should
improve the approximation by the convergence properties of HMC [Neal, 2011].

The main contribution of the paper lies in the use of Hamiltonian dynamics for
the mutation phase in a SMC method. In the literature, SMC algorithms are
generally found to use a standard Metropolis-Hastings step in their mutation
phase. Conversely, HMC methods using multiple particles do not have the
resampling phases 2a/2b.

###  4.2 Algorithm’s properties

Our HSMC method fits into the SMC framework described in Chopin [2004] and his
central limit theorem can be applied to our simulator. Provided that our
Hamiltonian mutation step preserves the distribution $f_{t}(\theta_{n})$, the
following convergences hold almost surely as $N\to\infty$ for any measurable
function $\phi$ such that the expectations below exists:

$N^{-1}\sum_{n=1}^{N}\phi(\theta_{n}^{(t)})\to\mathbb{E}_{f_{t}}[\phi(\theta_{n})]$
$\frac{\sum_{n=1}^{N}w_{n}^{(t)}\phi(\theta_{n}^{(t-1)})}{\sum_{n=1}^{N}w_{n}^{(t)}}\to\mathbb{E}_{f_{t}}[\phi(\theta_{n})]$
$N^{-1}\sum_{n=1}^{N}\phi(\hat{\theta}_{n}^{(t)})\to\mathbb{E}_{f_{t}}[\phi(\theta_{n})]$

The proof that our Hamiltonian transition kernel satisfies the conditions to
have $f_{t}(\theta_{n})$ as a stationary distribution can be found in the
review paper on Hamiltonian Monte Carlo by Neal [2011].

##  5 Method variations

While the Hamiltonian step is designed to be applied on continuous
distributions, it is easy to extend our method for spaces with discrete
dimensions. We can split our space $\Theta$ in two blocks
$\\{\Theta_{c},\Theta_{d}\\}$ where $\Theta_{c}$ includes the dimensions where
$\Theta$ is continuous and $\Theta_{d}$ includes the dimensions where $\Theta$
is discrete. From this separation into blocks, a standard Metropolis within
Gibbs [Gilks et al., 1995] step can be used with the Hamiltonian step being
used in the continuous block.

Another variation of the method consists in running the algorithm in parallel
for $J$ groups of $N$ particles. At the end, the sample properties of the $J$
groups can be compared. If the properties differ too much from each other, we
can suspect a convergence problem. This is the approach taken by Durham and
Geweke [2013] in their adaptive SMC simulator.

Finally, since the Hamiltonian step preserves the distribution of interest,
multiple steps can be made at each mutation phase. This solution increases the
performance of the method when particles are not exploring the $\Theta$ space
fast enough.

##  6 Examples

We used the HSMC method on kernel density estimates of two functions known to
challenge classical optimizers and MCMC simulators. The first function is
created for this paper and called the smiley function. It is a mixture of 3
Rosenbrock smile functions often found as an example to show the limitations
of MCMC algorithms. The second function is a dropwave function. Both functions
present multimodality and follow complex shapes. For ease of visualization, we
kept the functions two-dimentional.

We chose to simulate Gaussian kernel density estimates as they can mimic the
progressive convergence of several target functions when data are added
progressively. In our case, data are added by blocks of $100$ points and the
bandwidth of the kernel density estimate is $n^{-1/5}$. This rate has been
chosen as it is the order of the bandwidth reduction rate when using optimal
bandwidth for most distributions [Yatchew, 1998].

We also used the algorithm on a non-linear logit model where the log-
likelihood shows both multimodality and complex shapes.

###  6.1 Smiley kernel density estimate

We generated a sample of 2048 data points with coordinates $(x,y)$ using a
density proportional to the following function:

$\displaystyle g(x,y)=$
$\displaystyle\exp\left(\frac{1}{5}\left(-6\left(-(2.5\,-x)^{2}-1.5y+38\right)^{2}-(2.5\,-x)^{2}\right)\right)$
$\displaystyle+\exp\left(\frac{1}{5}\left(-6\left(-(x+2.5)^{2}-1.5y+38\right)^{2}-(x+2.5)^{2}\right)\right)$
$\displaystyle+\exp\left(\frac{1}{5}\left(-5\left(y-x^{2}\right)^{2}-x^{2}\right)\right)$

The contour plot and 3D plot of the function used can be found in Figure 5.

Figure 5: contour plot and 3D plot of the smiley function

We can easily see the multimodality and the elongated shapes on the function
plots. The data generated match the shape of the smiley function and have been
represented in Figure 6.

Figure 6: Generated data

The target function we want to simulate is a kernel density estimate using
these data points. We used 4 independent groups of 512 particles to simulate
the kernel density for a total of 2048 particles. The data points are
partitioned in 20 blocks of 100 points and 1 block of 48 points, for a total
of 21 blocks. Consequently, the sequence will include $T=21$ target functions
with the addition of the initial density. The initial density
$f_{0}(\theta_{n})$ used is a bivariate normal distribution with parameters
$\\{\mu_{1}=0,\mu_{2}=10,\sigma_{1}=10,\sigma_{2}=20,\rho=0\\}$. The HMC
tunning parameters $(M,L,\epsilon)$ used are $(I_{2},20,0.05)$, where $I_{2}$
is the $2\times 2$ identity matrix. The results of the HSMC simulation as well
as a smoothed histogram of the simulated points have been represented in
Figure 7. We can see that the HSMC method provided satisfying results and
successfully converged with only 21 iterations. The lowest mutations
acceptance rate we observed across several runs was $2043/2048$.

Figure 7: HSMC Simulation results

To illustrate the performance of our HSMC method, we also tried to simulate
the same sequence of functions using a parallel HMC approach with 21
iterations on the same kernel density. To do so we performed the mutation
phase 21 times on $f_{21}(\theta_{n})$. The result represented in Figure 8
show that several particles failed to converge. Moreover, there is no
guarantee that the particles around each mode are distributed according to the
mass around these modes in the target function. As particles tend to converge
to the closest mode, it is possible to have a first mode with twice the mass
of a second mode but only half of the particles around it.

Figure 8: Parallel HMC Simulation results

###  6.2 Constrained dropwave kernel density estimate

We generated a sample of 4096 data points with coordinates $(x,y)$ using a
density proportional to the following function defined on
$[-2.5;2.5]\times[-2.5;2.5]$:

$\displaystyle g(x,y)=$
$\displaystyle\exp\left(\frac{\cos\left(5\sqrt{x^{2}+y^{2}}\right)+1}{x^{2}+y^{2}+2}\right)$

The contour plot and 3D plot of the function used can be found in Figure 9.

Figure 9: contour plot and 3D plot of the dropwave function

The data generated match the shape of the dropwave function and have been
represented in Figure 10.

Figure 10: Generated data

The target function we want to simulate is a constrained kernel density
estimate using these data points. The Gaussian kernel is defined on
$\mathbb{R}^{2}$ but we want to limit the domain to
$[-2.5;2.5]\times[-2.5;2.5]$. To do so we keep the Gaussian kernel as is but
use the constrained HMC method during our mutation phase. We used 4
independent groups of 512 particles to simulate the kernel density for a total
of 2048 particles. The data points are partitioned in 40 blocks of 100 points
and 1 block of 96 points, for a total of 41 blocks. Consequently, the sequence
will include $T=41$ target functions with the addition of the initial density.
The initial density $f_{0}(\theta_{n})$ used is a bivariate normal with
parameters $\\{\mu_{1}=0,\mu_{2}=0,\sigma_{1}=10,\sigma_{2}=10,\rho=0\\}$. The
HMC tunning parameters $(M,L,\epsilon)$ used are $(I_{2},20,0.05)$, where
$I_{2}$ is the $2\times 2$ identity matrix. The results of the HSMC simulation
as well as a smoothed histogram of the simulated points are presented in
Figure 11. We can see that the constrained HSMC method also provided
satisfying results and successfully converged with 41 iterations. The lowest
mutations acceptance rate we observed across several runs was $2023/2048$.

Figure 11: CHSMC Simulation results

###  6.3 Non-linear logit

We use a Logit discrete choice model [Train, 2009] where the deterministic
part of the utility function is non-linear. In this hypothetical experiment,
an individual who has been given a black t-shirt is presented with another
t-shirt of a different color $x_{t}$. The individual can exchange it against
his own t-shirt $x_{0}$ or keep his current t-shirt. The individual chooses
$x_{t}$ or $x_{0}$ in order to maximize utility represented by:

$U(x_{t})=V(x_{t})+e_{1t}$ $U(x_{0})=e_{0t}$

where $e_{1t}$ and $e_{0t}$ are the standard extreme value errors. The
deterministic part of the utility function is represented by the fuction over
the color palette $x\in[-2;8]$:

$V(x)=2\frac{\sin(\beta_{2}\cdot x)}{1+0.5\cdot(\beta_{1}-x)^{2}}$

This function is represented on Figure 12 for the values $beta_{1}=3$ and
$beta_{2}=3$.

Figure 12: $V(x)$, $x\in[-2;8]$

The choice probabilities for $x_{t}$ and $x_{0}$ can be shown to equal:

$\frac{e^{V(x_{t})}}{1+e^{V(x_{t})}}\qquad,\qquad
1-\frac{e^{V(x_{t})}}{1+e^{V(x_{t})}}$

We can note that as $\beta_{1}\to\infty$, the choice probabilities converge to
$0.5$. We simulated $400$ experiments by drawing randomly an $x_{t}$ uniformly
on $[-2;8]$.

Figure 13: Theoretical Likelihood Figure 14: Theoretical Loglikelihood

A first look at the shape of the likelihood function represented in Figure 13
seems to indicate that the likelihood behaves nicely with a global maximum
around the true value. However, examining the log-likelihood on Figure 14
reveals the multimodality of the function.

The MH and HMC methods both fail ton converge after $5000$ iterrations and do
not get close to the global-maximum on numerous trials.

For the HMC method, the sequence of functions is created by adding the
observations progressively in the likelihood function, adding $50$
observations for each iteration. We compare the result to a standard SMC
algorithm.

Figure 15: Smoothed HSMC simulation

Figure 16: Smoothed SMC simulation

We can see that the HSMC method (Figure 16) could converge in only 4
iterations, whereas the SMC method (Figure 16) couldn’t converge completely.
First, in the SMC case, $4$ iterations might constitute a too short sequence
to obtain convergence. Adding more intermediate functions in the sequence
might improve convergence. Second, we notice that some of the particles in the
SMC case were trapped in a mode around the point $(4,-2)$ and these particles
are unlikely to be moved to the main area of mass concentration with more
iterations. With the HSMC algorithm, if some particles are trapped in a
different mode and misrepresent the mass in the underlying function, they get
relocated during the resampling phase. This is the consequence of using
$f_{t}(\theta_{n})/\hat{f}_{t-1}(\theta_{n})$ as a weight instead of the
theoretical $f_{t}(\theta_{n})/f_{t-1}(\theta_{n})$ used in standard SMC.

##  7 Future work and conclusion

We have shown that the HSMC algorithm is able to approximate by simulation
many functions with irregularities. The encouraging performance of the
algorithm finds direct potential applications in empirical work where standard
MCMC methods have shown limitations. One of this potential application is the
use of optimal instruments for non-linear BLP models in industrial
organization [Reynaert and Verboven, 2014].

## References

  * Bajari [2003] P. Bajari.  Discussion of Allenby, Chen and Yang.  _Quantitative Marketing and Economics_ , 1(3):277–283, 2003.
  * Burda [2015] M. Burda.  Constrained hamiltonian monte carlo in BEKK GARCH with targeting.  _Journal of Time Series Econometrics_ , 7(1):95–113, 2015.
  * Butucea et al. [2007] C. Butucea, M. Mougeot, K. Tribouley, et al.  Functional approach for excess mass estimation in the density model.  _Electronic Journal of Statistics_ , 1:449–472, 2007.
  * Cappe et al. [2004] O. Cappe, A. Guillin, J.-M. Marin, and C. P. Robert.  Population monte carlo.  _Journal of Computational and Graphical Sstatistics_ , 13:907–929, 2004.
  * Chernozhukov and Hong [2003] V. Chernozhukov and H. Hong.  An MCMC approach to classical estimation.  _Journal of Econometrics_ , 115(2):293–346, 2003\.
  * Chopin [2002] N. Chopin.  A sequential particle filter method for static models.  _Biometrika_ , 89(3):539–552, 2002.
  * Chopin [2004] N. Chopin.  Central limit theorem for sequential Monte Carlo methods and its application to bayesian inference.  _Annals of statistics_ , pages 2385–2411, 2004.
  * Del Moral et al. [2006] P. Del Moral, A. Doucet, and A. Jasra.  Sequential Monte Carlo samplers.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(3):411–436, 2006.
  * Delyon et al. [2016] B. Delyon, F. Portier, et al.  Integral approximation by kernel smoothing.  _Bernoulli_ , 22(4):2177–2208, 2016.
  * Doornik et al. [2000] J. A. Doornik, M. Ooms, et al.  Multimodality and the GARCH likelihood.  In _World Congress of the Econometric Society, Seattle, August_ , 2000\.
  * Durham and Geweke [2013] G. Durham and J. Geweke.  Adaptive sequential posterior simulators for massively parallel computing environments.  _Bayesian Model Comparison (Advances in Econometrics); Jeliazkov, I., Poirier, DJ, Eds_ , pages 1–44, 2013.
  * Gilks et al. [1995] W. R. Gilks, N. Best, and K. Tan.  Adaptive rejection metropolis sampling within gibbs sampling.  _Applied Statistics_ , pages 455–472, 1995.
  * Gourieroux and Monfort [1995] C. Gourieroux and A. Monfort.  _Statistics and econometric models_ , volume 1.  Cambridge University Press, 1995.
  * Hastings [1970] W. K. Hastings.  Monte carlo sampling methods using markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Herbst and Schorfheide [2014] E. Herbst and F. Schorfheide.  Sequential Monte Carlo sampling for DSGE models.  _Journal of Applied Econometrics_ , 29(7):1073–1098, 2014.
  * Hwang [1980] C.-R. Hwang.  Laplace’s method revisited: weak convergence of probability measures.  _The Annals of Probability_ , pages 1177–1182, 1980.
  * Jiang et al. [2009] R. Jiang, P. Manchanda, and P. E. Rossi.  Bayesian analysis of random coefficient logit models using aggregate data.  _Journal of Econometrics_ , 149(2):136–148, 2009\.
  * Koop and Potter [1999] G. Koop and S. M. Potter.  Bayes factors and nonlinearity: evidence from economic time series.  _Journal of Econometrics_ , 88(2):251–281, 1999\.
  * Koop and Potter [2000] G. Koop and S. M. Potter.  Nonlinearity, structural breaks or outliers in economic time series.  _Nonlinear Econometric Modeling in Time Series Analysis_ , pages 61–78, 2000.
  * Mardia and Watkins [1989] K. Mardia and A. Watkins.  On multimodality of the likelihood in the spatial linear model.  _Biometrika_ , 76(2):289–295, 1989.
  * Metropolis et al. [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.
  * Minnotte [1997] M. C. Minnotte.  Nonparametric testing of the existence of modes.  _The Annals of Statistics_ , pages 1646–1660, 1997.
  * Müller and Sawitzki [1991] D. W. Müller and G. Sawitzki.  Excess mass estimates and tests for multimodality.  _Journal of the American Statistical Association_ , 86(415):738–746, 1991.
  * Neal [2001] R. M. Neal.  Annealed importance sampling.  _Statistics and Computing_ , 11(2):125–139, 2001\.
  * Neal [2011] R. M. Neal.  MCMC using hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , 2:113–162, 2011\.
  * Reynaert and Verboven [2014] M. Reynaert and F. Verboven.  Improving the performance of random coefficients demand models: the role of optimal instruments.  _Journal of Econometrics_ , 179(1):83–98, 2014\.
  * Silverman [1981] B. W. Silverman.  Using kernel density estimates to investigate multimodality.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 97–99, 1981.
  * Train [2009] K. E. Train.  _Discrete choice methods with simulation_.  Cambridge university press, 2009.
  * Vieu [1996] P. Vieu.  A note on density mode estimation.  _Statistics & probability letters_, 26(4):297–307, 1996.
  * Yatchew [1998] A. Yatchew.  Nonparametric regression techniques in economics.  _Journal of Economic Literature_ , 36(2):669–721, 1998.
  * Zhou and Chen [2013] E. Zhou and X. Chen.  Sequential Monte Carlo simulated annealing.  _Journal of Global Optimization_ , 55(1):101–124, 2013.