# Detection Feasibility of H2 in Ultra-hot Jupiter Atmospheres

Anastasia Morgan  Laboratory for Atmospheric and Space Physics, University of
Colorado, Boulder, CO 80303, USA  P. Wilson Cauley  Laboratory for Atmospheric
and Space Physics, University of Colorado, Boulder, CO 80303, USA  Kevin
France  Laboratory for Atmospheric and Space Physics, University of Colorado,
Boulder, CO 80303, USA  Allison Youngblood  Exoplanets and Stellar
Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Tommi T. Koskinen  Lunar and Planetary Laboratory, University of Arizona,
Tucson, AZ 85721-0092, USA

###### Abstract

Ultra-hot Jupiters (UHJs) have recently been the focus of several atmospheric
studies due to their extreme properties. While molecular hydrogen (H2) plays a
key role in UHJ atmospheres, it has not been directly detected on an
exoplanet. To determine the feasibility of H2 detection via transmission
spectroscopy of the Lyman and Werner bands, we modeled UHJ atmospheres with H2
rotational temperatures varying from 2000 K to 4000 K orbiting A-type stars
ranging from $T_{eff}$ = 8,500 K to $T_{eff}$ = 10,300 K. We present simulated
transmission spectra for each planet-star temperature combination while adding
Poisson noise varying in magnitude from 0.5% to 2.0%. Finally, we cross-
correlated the spectra with expected atmospheric H2 absorption templates for
each temperature combination. Our results suggest that H2 detection with
current facilities, namely the Hubble Space Telescope, is not possible.
However, direct atmospheric transmission spectroscopy of H2 may be viable with
future UV-capable flagship missions.

##  1 Introduction

UHJs are a fascinating population of exoplanets due to their extreme
temperatures and resulting atmospheric chemistry, which more closely resembles
late K-type stars than a typical planetary atmosphere. Molecular hydrogen
plays a key role in heat transport in UHJ atmospheres as it dissociates and
recombines across the terminator (Bell & Cowan, 2018). Despite this and being
the most abundant molecule in gas giant atmospheres (Yelle, 2004), H2 has not
yet been directly detected in an exoplanet atmosphere. The large atmospheric
scale heights of UHJs present optimal observational opportunities to detect
H2, particularly in the far-UV (FUV) where H2 has strong rovibrational
transitions in the Lyman and Werner bands.

##  2 Methods

To assess the potential detectability of H2 in UHJ atmospheres, we simulated
transmission spectra of planets around A-type host stars with a range of
planet and stellar temperatures. Observation of planetary transmission spectra
with H2 transition lines between 1100 Å and 1400 Å requires a strong stellar
FUV flux (the background light source) which can be provided by early-type
host stars. We tested H2 rotational temperatures ranging from 2,000 K to 4,000
K in intervals of 500 K and stellar temperatures from 8,500 K to 10,300 K in
intervals of 300 K.

We fixed the planet radius at Rp = 1.5 RJ and assumed spherical symmetry for
the atmosphere. For all planets we adopted a H2 density profile derived for
KELT-20 b, which was calculated using GGchem (Woitke & Helling, 2021) and
assumed an isothermal atmosphere for each temperature in the explored range.
The isothermal assumption is reasonable given the small radial extent of the
H2 population. After using GGchem to calculate the mixing ratio of H2 and the
atmosphere’s mean molecular weight as a function of pressure, we solved the
hypsometric equation from Koskinen et al. (2022, see their eq. 10) to
determine the atmospheric height at each pressure. We then assumed a base
number density of 1017 cm-3 at $r=1.0R_{p}$ and a pressure of 1 bar, where the
H2 mixing ratio is 0.848, and calculated the H2 number density from $r=1.0\
R_{p}$ out to $r=1.3\ R_{p}$, where it drops to a negligible value of $\approx
10$ cm-3. We note that KELT-20 b is near the lower range of our planetary
temperatures and thus represents an optimistic case for the hotter planets
since H2 dissociates at higher atmospheric pressures as the temperature
increases, decreasing the depth of the H2 transmission spectrum.

We determined H2 rovibrational level populations using Maxwell-Boltzmann
statistics, i.e., we assumed thermal equilibrium (Barthélémy et al., 2007). We
then used the optical depth templates from McCandliss (2003) to calculate the
relative optical depth of each rovibrational transition. Next, we scale the
relative optical depth by the column density from our atmospheric model,
resulting in an annulus of H2 optical depth profiles around the planet.

We generated H2 transmission spectra for each planet and stellar temperature
combination at 10 equidistant transit times, which roughly mimics a 3-transit
HST program, where we determine the orbital parameters by assuming a circular
orbit and solving for the semi-major axis necessary to give the assumed
equilibrium temperature for the given host star $T_{eff}$. Finally, we add
Poisson noise at the level of 0.5%, 1.0%, 1.5%, and 2.0% to each individual
transmission spectrum. We note that we explored larger levels of Poisson noise
but found that anything beyond $\approx 2\%$ required a very large number of
simulated HST exposures to reach a significant detection.

For the cross-correlation analysis we first shifted the transmission spectra
into the rest frame of the planet for assumed orbital velocities between 0 km
s-1 and 500 km s-1. We averaged the resulting velocity-shifted transmission
spectra and cross-correlated this average spectrum with a pre-computed H2
absorption template for system velocities between $\pm$300 km s-1. We
restricted the cross-correlation calculation to between 1300 Å and 1400 Å
where there is sufficient stellar flux and H2 absorption lines to produce a
signal. We repeated this experiment 500 times for each planet-star combination
to test the robustness of the cross-correlation signal to the assumed noise
levels. We took the average cross-correlation significance from these 500
trials as the final detection significance for each case.

##  3 Results

Figure 1: (a) Significance map for stellar and planet temperature combinations
reflecting the average detection significance with 2% transmission spectrum
noise. (b) Noiseless transmission spectrum of H2 example for $T_{planet}$ =
4000 K and $T_{eff}$ = 9100 K. (c) Cross-correlation map example of
$T_{planet}$ = 4000 K and $T_{eff}$ = 9100 K showing a signal detection of
approximately 5$\sigma$.

Our simulations suggest that high signal-to-noise is required for H2
detection. As we show in Figure 1(c), H2 is detectable with a signal
significance of about 5$\sigma$ for a planet with $T_{planet}=4000$ K orbiting
a star with $T_{eff}$ = 9100 K and 2% transmission spectrum noise. The signal
map shown in Figure 1(a) indicates that the greatest detection of H2 is
approximately 5$\sigma$ and most likely to occur in higher temperature planets
orbiting lower temperature stars for the considered temperature ranges. We
note that this trend occurs due to the fixed planetary radius with increasing
stellar radius; a fixed value of $R_{p}/R_{*}$ would likely produce better
cross-correlation signals due to the increased signal-to-noise of the star’s
FUV flux. This trend is consistent with the lower transmission spectrum noise
levels tested, where we measured signal detection significance values of
$\approx 12\sigma$ at 0.5% noise, $\approx 10\sigma$ at 1.0% noise, and
$\approx 8\sigma$ at 1.5%.

##  4 Conclusion

Given our assumed optimistic H2 density profile, the simulated detection
significance values are likely still unachievable with HST. For example, using
the STIS E140M exposure time calculator we estimate that the best signal-to-
noise currently achievable for the FUV transmission spectrum of KELT-9 b, the
hottest known UHJ orbiting a bright A0V star, is $\approx 3\%-5\%$. We
estimate that a $5\sigma$ detection for KELT-9 b would require 6-7 transit
observations of 5 HST orbits each, or 30-35 total orbits. The $\leq 2\%$ noise
levels required for detection via a reasonable amount of observing time will
likely be achievable with next generation UV instruments, such as LUVOIR and
the LUMOS spectrograph (France et al., 2017), although systematics may push
the necessary transmission spectrum precision even lower. Finally, it is worth
noting that our assumption of LTE for the H2 gas in the lower atmosphere needs
to be tested with models which account for the depth at which FUV photons are
deposited in the atmosphere, potentially pumping the H2 into non-thermal
distributions.

## References

  * Barthélémy et al. (2007) Barthélémy, M., Lilensten, J., & Parkinson, C. D. 2007, A&A, 474, 301, doi: 10.1051/0004-6361:20065891
  * Bell & Cowan (2018) Bell, T. J., & Cowan, N. B. 2018, The Astrophysical Journal Letters, 857, L20
  * France et al. (2017) France, K., Fleming, B., West, G., et al. 2017, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10397, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1039713, doi: 10.1117/12.2272025
  * Koskinen et al. (2022) Koskinen, T. T., Lavvas, P., Huang, C., et al. 2022, The Astrophysical Journal, 929, 52
  * McCandliss (2003) McCandliss, S. R. 2003, Publications of the Astronomical Society of the Pacific, 115, 651
  * Woitke & Helling (2021) Woitke, P., & Helling, C. 2021, Astrophysics Source Code Library, ascl
  * Yelle (2004) Yelle, R. V. 2004, Icarus, 170, 167


[a,b]Nobuyuki Matsumoto

# Decimation map in 2D for accelerating HMC

Richard C. Brower  Taku Izubuchi

###### Abstract

To accelerate the HMC with field transformation, we consider a variant of the
trivializing map, the decimation map, which can be regarded as a coarse-
graining transformation. Using the 2D $U(1)$ pure gauge model, combined with
the guided Monte Carlo algorithm, we show that the integrated autocorrelation
time of the topological charge can be exponentially improved in the wall clock
time. Our study indicates that incorporating renormalization group picture is
a powerful and essential ingredient to accelerate the HMC at large $\beta$.

##  1 Introduction

Critical slowing down is intrinsic to Monte Carlo study of lattice quantum
field theories, but is not yet resolved satisfactorily. The state of the art
algorithm for the gauge field generation in QCD is the Hybrid Monte Carlo
(HMC) [1], and there are two major directions to cope with the critical
slowing down in this algorithmic framework (see [2, 3] for reviews). One is to
align the velocity in the molecular dynamics (MD) among all the Fourier modes
[4, 5, 6, 7, 8], and the other is to construct a field transformation such
that the resulting effective action has advantageous sampling properties [10,
11] (see also [12, 13, 14, 15, 16, 17, 18]).

This study follows the latter approach following Lüscher’s seminal work [10],
and we consider trivialization of link variables. A difference, however, is
that we divide the lattice into local blocks while the original work
considered the global trivialization. In fact, our field transformation
corresponds to eliminating links from the theory, namely decimation of the
variables. When a link is trivialized, the effective action for the
surrounding links exactly agrees with that of the integrated theory. Our
approach thus incorporates the idea of renormalization group.

In general, such an algorithm can be extremely complex. We thus first test our
idea with the simplest gauge theory possible: the pure $U(1)$ gauge theory in
two spacetime dimensions. The peculiarity of this model is that there is no
propagation mode; the correlation length defined by the plaquette correlation
function is zero. A nontrivial quantity is the topological charge under the
periodic boundary condition, which counts the winding number of $U(1)$
plaquette angles around the faces, i.e., the number of vortices. The density
of vortices is determined by the susceptibility, which is finite in the
continuum limit. In other words, there is a typical volume scale in which a
vortex appears. Accordingly, we need to move the variables collectively on
fine lattices to induce topological tunneling. In fact, in the HMC, we observe
that the tunneling becomes exponentially hard as in QCD. The goal of this
study is to demonstrate that the decimation map can reduce the autocorrelation
of the topological charge significantly.

The algorithm consists of two parts as in the original trivializing map. We
first construct a series of local trivializing maps, namely the decimation
map, by solving the linear equation for the gradient flow kernel. We here
parameterize the function space with Wilson loops and directly invert the
linear equation with an iterative solver. We then use the guided Monte Carlo
[19] for configuration generation, which is a variant of the HMC replacing the
action in the MD Hamiltonian by an approximate one. The effects of finite flow
step size and the approximated action are under control and the volume scaling
is power-law while keeping the exact detailed balance.

We test our algorithm on a small system of the physical volume $V_{\rm
phys}=6^{2}/g^{2}$, where $g$ is the dimensionful coupling constant. This
volume corresponds to the typical scale of a single vortex. We show that we
obtain $\times 73$ speedup in the integrated autocorrelation time in wall
clock at $\beta=7.1$ ($16\times 16$ lattice), and furthermore, the exponent
towards the continuum limit is decreased with a factor of $0.62$. Our
decimation map exemplifies that the HMC of a gauge theory can be accelerated
with a suitable field transformation.

##  2 Decimation map

The idea of the decimation map can be described simply in the 2D $U(1)$ model
with the Wilson action:

$\displaystyle S(U)\equiv-\beta\sum_{x}\cos\kappa_{x},$ (1)

where $\beta=1/(ag)^{2}$ with the lattice spacing $a$ and the coupling $g$
that determines the scale in the system. $\kappa_{x}$ is the plaquette angle:

$\displaystyle\kappa_{x}\equiv\frac{1}{i}\log(U_{x,0}U_{x+0,1}U_{x+1,0}^{\dagger}U_{x,1}^{\dagger}).$
(2)

As in figure 1, we iteratively choose a set of independent link variables to
trivialize.

Figure 1: The basic idea of the decimation map is a successive local
trivialization of link variables. At each stage, a set of independent link
variables (colored with magenta) are chosen and trivialized.

Each stage of trivialization doubles the size of fundamental Wilson loops in
the theory.

Let us write by $U^{\prime}_{x,\mu}$ the link variables to be trivialized and
by $\tilde{U}_{x,\mu}$ the remaining variables. Once the trivialization is
performed successfully with the map:

$\displaystyle U^{\prime}_{x,\mu}={\cal F}_{x,\mu}(V;\tilde{U}),$ (3)

where $V_{x,\mu}$ are the trivialized variables, the action will be
transformed accordingly from the original one $S=S(U^{\prime},\tilde{U})$ to:

$\displaystyle S_{\rm eff}(\tilde{U})\equiv
S(U^{\prime}(V;\tilde{U}),\tilde{U})-\ln\det{\cal F}^{*}(V;\tilde{U}),$ (4)

where ${\cal F}^{*}$ is the Jacobian of the map ${\cal F}$. Then one can
immediately see that $S_{\rm eff}$ agrees with the action we obtain after
integrating over $U^{\prime}$ from the original action:

$\displaystyle\int(dU^{\prime})(d\tilde{U})e^{-S(U^{\prime},\tilde{U})}=\int(d\tilde{U})e^{-S_{\rm
eff}(\tilde{U})}.$ (5)

We thus see that $S_{\rm eff}(\tilde{U})$ describes a coarse-grained system
and each stage of the map corresponds to the decimation of links.

Such observation has an important implication that the decimation map drives
the system away from the continuum limit, where the HMC becomes inefficient
(see figure 2).

Figure 2: A qualitative picture describing the resulting flow in the action
space. It is likely that the theory approaches towards the renormalized
trajectory (RT), observing that the tunneling rate is not as high as for the
Wilson actions of the larger lattice spacings, at least for the first a few
stages. However, it is uncertain how far this picture holds in the flow and
its precise structure is out of the scope of this work.

One can therefore expect that HMC in the decimated system has smaller
autocorrelation than in the original system. After obtaining configurations in
the decimated system, one can use ${\cal F}$ to supply (or integrate in) the
decimated variables to calculate observables for the original action.

##  3 Obtaining the map numerically

In our simple example, there are only two Wilson loops attached to each link
$U_{x,\mu}$. Let us write by $\kappa_{1,2}$ the two phase angles of the loops.
One can take the basis for the function space to trivialize $U_{x,\mu}$ as:

$\displaystyle
e_{m_{1},m_{2}}\equiv\cos(m_{1}\kappa_{1}+m_{2}\kappa_{2})\quad(m_{1},m_{2}\in\mathbb{Z}).$
(6)

Note that when $m_{1}=m_{2}$, they do not depend on $U_{x,\mu}$, and thus
$e_{m,m}$ $(m\in\mathbb{Z})$ are the zero modes of the differential operators.

As in [10], we take the gradient flow ansatz:

$\displaystyle\dot{U}_{t,x,\mu}U_{t,x,\mu}^{-1}=\partial K_{t},\quad
U_{t,x,\mu}\equiv{\cal F}_{t,x,\mu}(V),\quad{\cal F}_{t=1}={\cal F},$ (7)

where $\partial$ is the derivative with respect to $\phi\equiv\log U_{x,\mu}$.
Then, one can write down the equation for $K_{t}$ by demanding the effective
action to decrease linearly in $t$. For example, for the first stage:

$\displaystyle[-\partial^{2}+\beta
t(\sin\kappa_{1}-\sin\kappa_{2})\partial]K_{t}=\beta(\cos\kappa_{1}+\cos\kappa_{2})+\sum_{m}e_{m,m}v_{t,m},$
(8)

For the higher stages, terms with various $m_{i}$ appear on both hand sides.
By representing the differential operator in the basis (6), one can
numerically solve the linear equation (8) with, e.g., CGNE with
preconditioning by $\partial^{-2}$ in the space without the zero modes. The
zero mode part $\sum_{m}e_{m,m}v_{t,m}$ can be obtained acting the
differential operator to the solution after solving the equation in the
projected space. We find that for $\beta\lesssim 8.9$, the range
$|m_{1,2}|\leq 128$ is sufficient (see figure 5). The inversion takes a few
minutes to about a day using GPU depending on the value of $\beta$ and the
stage of decimation. The details of the inversion will be described in the
subsequent paper.

In practice, we discretize the flow with the Runge-Kutta methods. We here
choose the midpoint integrator. For the map to be bijective, we ensure that
the step size is kept within a bound that can be derived as in [10]. For the
midpoint integrator, given a kernel of the form:

$\displaystyle
K_{t}(\phi)=\sum_{m_{1},m_{2}}e_{m_{1},m_{2}}c_{t,m_{1},m_{2}},$ (9)

the step size $\epsilon$ must satisfy:

$\displaystyle\epsilon B_{t+\epsilon/2}\cdot\left(1+\frac{\epsilon
B_{t}}{2}\right)<1,\quad
B_{t}\equiv\sum_{m_{1},m_{2}}(m_{1}-m_{2})^{2}|c_{t,m_{1},m_{2}}|.$ (10)

##  4 Guided Monte Carlo

The functions $e_{m,m}$ are the Wilson loops of the double size winding $m$
times. They are independent of $U_{x,\mu}$, but have different dependencies on
the surrounding links. They are the basis functions for determining the
weights in the transformed system. In fact, the integral of the zero mode over
$t$ gives the next effective action:

$\displaystyle S_{\rm eff}=\sum_{m}e_{m,m}\int dt\,v_{t,m}.$ (11)

One can use this knowledge to simplify the force calculation with the help of
the guided Monte Carlo algorithm [19] as below.

The key point of the guided Monte Carlo algorithm is that, as far as the time
evolution in the HMC is symplectic, the detailed balance is exact. However,
since the energy conservation is used to have nonvanishing acceptance in the
Metropolis test, it is important that the Hamiltonian associated to the
symplectic integrator is close to the exact one.

We replace the action in the MD Hamiltonian by an approximate one to evade the
evaluation of the gradient of the Jacobian as well as the force propagation.
Since the zero modes $v_{t,m}$ are calculable at each flow time $t$ as a sum
of Wilson loops, we can calculate the effective action directly in the loop
space using eq. (11). We use the Simpson’s formula for approximating the right
hand side, which involves only one additional evaluation of $K_{t}$ at $t=0$
to those calculated for the midpoint integrator. The difference between this
approximate action and the effective action of the discretized flow comes from
the error in the midpoint integrator and is $O(\epsilon^{2})$, where
$\epsilon$ is the step size of the flow. Note that, though for a fixed
$\epsilon$ the acceptance rate becomes exponentially small when enlarging the
volume, this effect can be compensated by decreasing the step size $\epsilon$
of the trivializing flow just as decreasing the MD step size in the
conventional HMC. The cost scaling in enlarging the volume is therefore power
law.

We set the approximate effective action in the trivialized region to be
constant in the MD, and thus we choose the inner links purely randomly. This
corresponds to assigning the zero mass in the MD. Nontrivial updates are
performed for the outer remaining links through the Wilson loops in the large
unit. In the Metropolis test, however, since the replacement involves
discretization error of the flow, we need to include all the updated links in
the weight function.

In practice, we choose the step size adaptively by making it proportional to
the bijection bound (10) while the overall scaling is determined to keep the
acceptance around 0.8. The adaptively chosen step size becomes extremely small
at large $\beta$, however, for the following two reasons. One is that the
number of relevant basis functions grows rapidly as $\beta$ is increased,
which will in turn decrease the bijection bound (10) (see also figure 4 on
this rapid growth). The other is to control the acceptance rate in the guided
Monte Carlo for increasing the lattice sites. In our calculation, the second
point is the bottleneck. This cost due to the latter point may be circumvented
by using a higher-order integrator or by switching to calculating the exact
force of the effective action as in [10]. As will be shown in section 5, the
exponent of the cost towards the continuum limit can be decreased with the
decimation map, and thus one can expect an acceleration with the latter
algorithm. However, this point needs to be verified with an additional study.

The separation of the inner and outer links corresponds to separating the UV
and IR modes in the system. For simplicity, suppose that the decimation map is
obtained without discretization error. One can then choose the inner links
without autocorrelation. However, since the environmental outer variables are
fixed, this fluctuation is only around a fixed background configuration and
depends only on the local links surrounding the trivialized region. It can be
thus interpreted as the UV fluctuation. On the other hand, the update of the
outer links is completely insensible to the inner links, whose dynamics is
determined solely by the coarse-grained action. Thus, the global update of the
outer links corresponds to updating the IR modes. By using this separation, we
can investigate which part of the algorithm is responsible for accelerating
the algorithm.

##  5 Results

To enumerate the effectiveness of our approach, we calculate the integrated
autocorrelation time:

$\displaystyle\tau_{\rm int}({\cal O})\equiv\frac{1}{2}+\sum_{i\geq
1}\frac{\langle{\cal O}_{0}{\cal O}_{i}\rangle}{\langle{\cal O}_{0}{\cal
O}_{0}\rangle},$ (12)

where ${\cal O}_{i}\equiv{\cal O}(U_{i})$ is the value of the observable
${\cal O}$ evaluated with the $i$-th configuration $U_{i}$. Our interest is in
the topological charge [20, 21]:

$\displaystyle Q\equiv-\frac{1}{2\pi}\sum_{x}\kappa_{x},$ (13)

which is integer-valued. It is assumed to take the principal branch of the
logarithm to evaluate $\kappa_{x}$ in eq. (2). One can derive analytically
that the susceptibility approaches in the continuum limit to:

$\displaystyle\chi_{Q}\equiv\frac{\langle Q^{2}\rangle}{V_{\rm
phys}}\to\frac{g^{2}}{(2\pi)^{2}}\quad(a\to 0),$ (14)

which shows that a vortex (an instanton) appears typically in a volume of
$(2\pi)^{2}/g^{2}$.

We fix the physical volume to $V_{\rm phys}=6^{2}/g^{2}$, which is about the
scale of a vortex, and take the continuum limit $a\to 0$ by fixing $g$. The
trajectory length of the MD is fixed to $1.0$ in all simulations. In the
conventional HMC, the MD step size is scaled such that the acceptance rate is
around 0.8. In the guided Hamiltonian simulations, the step size is fixed to
0.05 (i.e., 20 steps); though this is not optimal, this number is not relevant
to the actual cost because the algorithmic overhead is in the field
transformation part. As mentioned in section 4, flow step size of the
decimation map determines the acceptance rate, which is scaled for the
acceptance to be around 0.8.

Figure 3 shows the scaling of $\tau_{\rm int}(Q)$ in the units of Monte Carlo
steps and wall clock time. The number of stages is taken up to four. The
computation is performed on CPU without parallelization.

Figure 3: $\tau_{\rm int}(Q)$ in Monte Carlo steps (left) and in wall clock
time (right). The conventional HMC (blue solid line) is compared to the guided
Monte Carlo with the decimation map up to four stages. In Monte Carlo step,
$\tau_{\rm int}(Q)\simeq 8000$ at $\beta=7.1$ ($16\times 16$ lattice) is
reduced to $\tau_{\rm int}(Q)\simeq 1$. In wall clock, the speedup at
$\beta=7.1$ is x73 and the decrease of the exponent is by a factor 0.62 with
four stages.

As shown in the left panel, for $\beta=7.1$ ($16\times 16$ lattice),
$\tau_{\rm int}(Q)\simeq 8000$ with the conventional HMC is reduced to
$\tau_{\rm int}(Q)\simeq 1$ with the decimation map in the Monte Carlo unit.
In wall clock, we see x73 speedup at $\beta=7.1$, and furthermore, the
exponent has decreased with a factor of 0.62.

Concerning the algorithmic overhead, it is the second stage of decimation that
is the most costly (see figure 4).

Figure 4:  The time consumed to generate a configuration in wall clock. Since
the second stage is the most costly, this part dominates the time even when we
include higher stages.

This is simply because the number of relevant basis functions is the largest.
Figure 5 shows the magnitude of the coefficients $c_{t,m_{1},m_{2}}$ of eq.
(9) for $\beta=8.9$ at $t=1$ (corresponding to the finest lattice during the
transformation). After the first stage, the functions with large $m_{1}-m_{2}$
become relevant as the effective action contains terms of higher
representations, multiply winded loops. At the later stages, the number of
relevant basis functions decreases reflecting that the lattice is becoming
coarse.

Figure 5: Density plot of the magnitude of the coefficient function in the
flow kernel at $t=1$ (corresponding to the finest lattice during the
transformation) for $\beta=8.9$. The first to fourth stages from top left to
right bottom. The relevant basis functions expand in the $m_{1}-m_{2}$
direction after the first stage because the effective action includes multiply
winding Wilson loops. After the second stage, the number of relevant functions
decreases reflecting that the lattice action is becoming coarse.

Figure 6 shows the time series of $Q$ at $\beta=7.1$ with and without the
decimation map. With the conventional HMC, we observe that $Q$ is varying
extremely slowly in large scale; the fast fluctuation only moves $Q$ back and
forth between two nearby sectors. By implementing the decimation map, the
fluctuation becomes centered correctly at $Q=0$ without visible long-range
autocorrelation even though a single update is taking a noticeable time.

Figure 6: Time series of the topological charge in wall clock with the
conventional HMC (left) and with the decimation map with four stages (right)
at $\beta=7.1$. With the conventional HMC, though $Q$ changes back and forth
frequently within two nearby sectors, it takes a long time to cover all the
topological sectors because of the large autocorrelation. With the decimation
map, on the other hand, though each update takes a noticeable time, it spans
all the sectors with a few Monte Carlo steps thanks to the minimal
autocorrelation.

Finally, to study whether it is the UV or IR part of the algorithm that is
generating the acceleration, we separate the inner (UV) and outer (IR) updates
by freezing one of them alternatively, and measure the tunneling rate:

$\displaystyle R\equiv\langle|Q_{i}-Q_{i+1}|\rangle.$ (15)

Figure 7 shows $R$ for each part with various $\beta$.

Figure 7:  Comparison of the tunneling rate $R$ compared between the inner UV
update and the outer IR update with three and four stages of decimation. For
small $\beta$, since a vortex can fit into a local trivialized region easily,
the inner update gives larger $R$. However, as we increase $\beta$, the link
variables start to form a continuum, and consequently they need to be updated
collectively. Accordingly, the outer IR update becomes more effective at large
$\beta$.

For small $\beta$, the volume of the trivialized region is large enough to
create a vortex. Correspondingly, the inner update has a larger tunneling rate
compared to the outer update. However, as we enlarge $\beta$, the chance of
creating a vortex inside the trivialized region becomes small, and the
effectiveness flips.

As discussed in section 4, the outer update is necessary to change the global
configuration. Consequently, even when the topological charge changes in the
UV update, the central value of the fluctuation will not change unless the IR
modes are altered. In large $\beta$ theories, since the problem is to update
the IR modes effectively, mapping the theory to coarse lattice actions by
incorporating renormalization group seems crucial for the speed up.

##  6 Conclusion and Outlook

We considered the decimation map in the 2D $U(1)$ model that can be regarded
as a coarse-graining transformation. With this map, combined with the guided
Monte Carlo algorithm, we showed that the integrated autocorrelation time of
the topological charge can be exponentially improved in the wall clock time.
We observe x73 speed up on the $16\times 16$ lattice at $\beta=7.1$ and the
decrease of the exponent by a factor of 0.62. Our study exemplifies that it is
possible to accelerate the HMC in gauge theories once an appropriate field
transformation is constructed. We further argued that the coarse-graining
picture is crucial for accelerating the HMC at large $\beta$.

Towards the application in QCD, there are three points to be addressed:
generalization to non-Abelian groups, generalization to higher dimensions, and
inclusion of the fermion; and none of them is straightforward. However,
machinery of solving the linear equation itself can be applied to non-Abelian
cases in principle by parameterizing the function space with the Wilson loops.
The nontrivial relations among the Wilson loops called Mandelstam constraints
[22] (see also [18]) will not be an issue in finding a solution with iterative
methods. Generalization to higher dimension needs work even for the $U(1)$
case because the variables correlate through various directions and the
effective action involves loops with complicated shapes. Nevertheless, if we
can trivialize codimension-one hypersurfaces completely, the resulting
effective theory will be described by the double-sized Wilson loops by gauge
invariance. However, aiming further to include fermion, it must be necessary
to find an approximation scheme rather than solving the equation exactly while
retaining the renormalization-group picture. Toy 2D models such as the
$CP^{N-1}$ model and the Schwinger model may be useful for such
investigations. Studies along these lines are in progress and will be reported
elsewhere.

## Acknowledgments

The authors thank Tom Blum, Peter Boyle, Norman Christ, Gerald Dunne, Sam
Foreman, Anna Hasenfratz, Luchang Jin, Xiaoyong Jin, Chulwoo Jung, James
Osborn, Akio Tomiya, Julian Urban and Urs Wenger for valuable discussions.
This work is supported by the U.S. Department of Energy under contract No. DE-
AC02-07CH11359. N.M. was supported by the Special Postdoctoral Researchers
Program of RIKEN.

## References

  * [1] S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, _Hybrid Monte Carlo_ , Phys. Lett. B 195, 216-222 (1987).
  * [2] G. Kanwar, _Flow-based sampling for lattice field theories_ , PoS(LATTICE2023)114.
  * [3] P. Boyle, _Advances in algorithms for solvers and gauge generation_ , PoS(LATTICE2023)122.
  * [4] G. Parisi, _Prolegomena to Any Future Computer Evaluation of the QCD Mass Spectrum_ , Progress in Gauge Field Theory, edited by ’t Hooft et al. (Plenum, New York, 1984), p.351.
  * [5] G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, B. Svetitsky and K. G. Wilson, _Langevin Simulations of Lattice Field Theories_ , Phys. Rev. D 32, 2736 (1985).
  * [6] C. T. H. Davies, G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, K. G. Wilson, P. Rossi and B. Svetitsky, _Fourier Acceleration in Lattice Gauge Theories. 1. Landau Gauge Fixing_ , Phys. Rev. D 37, 1581 (1988).
  * [7] G. Katz, G. Batrouni, C. Davies, A. S. Kronfeld, P. Lepage, P. Rossi, B. Svetitsky and K. Wilson, _Fourier Acceleration. 2. Matrix Inversion and the Quark Propagator_ , Phys. Rev. D 37, 1589 (1988).
  * [8] C. T. H. Davies, G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, P. Rossi, B. Svetitsky and K. G. Wilson, _Fourier Acceleration in Lattice Gauge Theories. 3. Updating Field Configurations_ , Phys. Rev. D 41, 1953 (1990).
  * [9] T. Nguyen, P. Boyle, N. H. Christ, Y. C. Jang and C. Jung, _Riemannian Manifold Hybrid Monte Carlo in Lattice QCD_ , PoS LATTICE2021, 582 (2022) [arXiv:2112.04556 [hep-lat]].
  * [10] M. Luscher, _Trivializing maps, the Wilson flow and the HMC algorithm_ , Commun. Math. Phys. 293, 899-919 (2010) [arXiv:0907.5491 [hep-lat]].
  * [11] G. P. Engel and S. Schaefer, _Testing trivializing maps in the Hybrid Monte Carlo algorithm_ , Comput. Phys. Commun. 182, 2107-2114 (2011) [arXiv:1102.1852 [hep-lat]].
  * [12] M. S. Albergo, G. Kanwar and P. E. Shanahan, _Flow-based generative models for Markov chain Monte Carlo in lattice field theory_ , Phys. Rev. D 100, no.3, 034515 (2019) [arXiv:1904.12072 [hep-lat]].
  * [13] S. Foreman, X. Y. Jin and J. C. Osborn, _Deep Learning Hamiltonian Monte Carlo_ , [arXiv:2105.03418 [hep-lat]].
  * [14] D. Albandea, P. Hernández, A. Ramos and F. Romero-López, _Topological sampling through windings_ , Eur. Phys. J. C 81, no.10, 873 (2021) [arXiv:2106.14234 [hep-lat]].
  * [15] S. Foreman, T. Izubuchi, L. Jin, X. Y. Jin, J. C. Osborn and A. Tomiya, _HMC with Normalizing Flows_ , PoS LATTICE2021, 073 (2022) [arXiv:2112.01586 [cs.LG]].
  * [16] M. Caselle, E. Cellini, A. Nada and M. Panero, _Stochastic normalizing flows as non-equilibrium transformations_ , JHEP 07, 015 (2022) [arXiv:2201.08862 [hep-lat]].
  * [17] S. Bacchio, P. Kessel, S. Schaefer and L. Vaitl, _Learning trivializing gradient flows for lattice gauge theories_ , Phys. Rev. D 107, no.5, L051504 (2023) [arXiv:2212.08469 [hep-lat]].
  * [18] P. Boyle, T. Izubuchi, L. Jin, C. Jung, C. Lehner, N. Matsumoto and A. Tomiya, _Use of Schwinger-Dyson equation in constructing an approximate trivializing map_ , PoS LATTICE2022, 229 (2023) [arXiv:2212.11387 [hep-lat]].
  * [19] A. M. Horowitz, _A Generalized guided Monte Carlo algorithm_ , Phys. Lett. B 268, 247-252 (1991).
  * [20] M. Luscher, _Topology of Lattice Gauge Fields_ , Commun. Math. Phys. 85, 39 (1982).
  * [21] A. Phillips, _Characteristic Numbers of U(1) Valued Lattice Gauge Fields_ , Annals Phys. 161, 399-422 (1985).
  * [22] S. Mandelstam, _Charge - Monopole Duality and the Phases of Nonabelian Gauge Theories_ , Phys. Rev. D 19, 2391 (1979).