# Efficient Solutions of Fermionic Systems using Artificial Neural Networks

Even M. Nordhagen  Department of Physics and Njord Center, University of Oslo,
N-0316 Oslo, Norway  Jane M. Kim  Department of Physics and Astronomy and
Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI
48824, USA  Bryce Fore  Physics Division, Argonne National Laboratory,
Argonne, IL 60439, USA  Alessandro Lovato  Physics Division, Argonne National
Laboratory, Argonne, IL 60439, USA  Morten Hjorth-Jensen  Department of
Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State
University, East Lansing, MI 48824, USA  Department of Physics and Center for
Computing in Science Education, University of Oslo, N-0316 Oslo, Norway

###### Abstract

We discuss differences and similarities between variational Monte Carlo
approaches that use conventional and artificial neural network
parameterizations of the ground-state wave function for systems of fermions.
We focus on a relatively shallow neural-network architectures, the so called
restricted Boltzmann machine, and discuss unsupervised learning algorithms
that are suitable to model complicated many-body correlations. We analyze the
strengths and weaknesses of conventional and neural-network wave functions by
solving various circular quantum-dots systems. Results for up to 90 electrons
are presented and particular emphasis is placed on how to efficiently
implement these methods on homogeneous and heterogeneous high-performance
computing facilities.

##  I Introduction

Solving the Schrödinger equation for systems of many interacting bosons or
fermions is classified as an NP-hard problem due to the complexity of the
required many-dimensional wave function, resulting in an exponential growth of
degrees of freedom. Reducing the dimensionalities of quantum mechanical many-
body systems is an important aspect of modern physics, ranging from the
development of efficient algorithms for studying many-body systems to
exploiting the increase in computing power. To write software that can fully
utilize the available resources has long been known to be an important aspect
of these endeavors. Despite tremendous progress has been made in this
direction, traditional many-particle methods, either quantum mechanical or
classical ones, face huge dimensionality problems when applied to studies of
systems with many interacting particles.

Over the last two decades, quantum computing and machine learning have emerged
as some of the most promising approaches for studying complex physical systems
where several length and energy scales are involved. Machine learning
techniques and in particular neural-network quantum states Goodfellow  _et
al._ (2016) have recently been applied to studies of many-body systems, see
for example Refs. Carleo and Troyer (2017); Carrasquilla and Torlai (2021);
Pfau  _et al._ (2020); Calcavecchia  _et al._ (2014); Carleo  _et al._ (2019);
Boehnlein  _et al._ (2022); Adams  _et al._ (2021); Lovato  _et al._ (2022),
in various fields of physics and quantum chemistry, with very promising
results. In many of these studies, one has obtained results that align well
with exact analytical solutions or are in close agreement with state-of-the-
art quantum Monte Carlo calculations.

The variational and diffusion Monte Carlo algorithms are among the most
popular and successful methods available for ground-state studies of quantum
mechanical systems. They both rely on a suitable ansatz for the ground-state
of the system, often dubbed trial wave function, which is defined in terms of
a set of variational parameters whose optimal values are found by minimizing
the total energy of the system. Devising flexible and accurate functional
forms for the trial wave functions requires prior knowledge and physical
intuition about the system under investigation. However, for many systems we
do not have this intuition, and as a result it is often difficult to define a
good ansatz for the state function.

According to the universal approximation theorem, a deep neural network can
represent any continuous function within a certain error Hornik  _et al._
(1989) — see also Refs. Murphy (2012); Hastie  _et al._ (2009); Bishop (2006);
Goodfellow  _et al._ (2016) for further discussions of deep leaning methods.
Since the variational state wave function in principle can take any functional
form, it is natural to replace the trial wave function with a neural network
and treat it as a machine learning problem. This approach has been
successfully implemented in recent works, see for example Refs. Pfau  _et al._
(2020); Carleo and Troyer (2017); Cassella  _et al._ (2022); Adams  _et al._
(2021); Lovato  _et al._ (2022), and forms the motivation for the present
study. Here, the neural network of choice was derived from so-called
restricted Boltzmann machines, much inspired by the recent contributions by
Carleo et al., see for example Refs. Carleo and Troyer (2017); Carleo  _et
al._ (2019). Note that neural-networks representations of variational states
are more general, as they do not in principle require prior knowledge on the
ground-state wave function, thereby opening the door to systems that have yet
to be solved. Particular attention however has to be devoted to the symmetries
of the problem, whose inclusion is critical to achieve accurate results .

In this work, we will focus on systems of electrons confined to move in two-
dimensional harmonic oscillator systems, so-called quantum dots. These are
strongly confined electrons and offer a wide variety of complex and subtle
phenomena which pose severe challenges to existing many-body methods. Due to
their small size, quantum dots are characterized by discrete quantum levels.
For instance, the ground states of circular dots show similar shell structures
and magic numbers as seen for atoms and nuclei. These structures are
particularly evident in measurements of the change in electrochemical
potential due to the addition of one extra electron. Here these systems will
serve as our test of the applicability of artificial neural network
variational states, including restricted Boltzmann Machines.

The theoretical foundation and the methodology are explained in section II.
The subsequent sections present our results with an analysis of computational
methods and resources. In the last section we present our conclusions and
perspectives for future work.

##  II Method

For any Hamiltonian $\hat{\mathcal{H}}$ and trial wave function $\psi_{T}$,
the variational principle guarantees that the expectation value of the energy
$E_{T}$ is greater than or equal to the true ground state energy $E_{0}$,

$E_{0}\leq
E_{T}=\frac{\langle\psi_{T}|\hat{\mathcal{H}}|\psi_{T}\rangle}{\langle\psi_{T}|\psi_{T}\rangle}.$
(1)

Thus approximate solutions to the time-independent Schrödinger equation can be
obtained by choosing a careful parameterization of the wave function and
minimizing the energy $E_{T}$ with respect to the parameters. Since the
integrals representing $E_{T}$ are normally high dimensional, it is most
efficient to evaluate them by means of Monte Carlo methods

$E_{T}\approx\langle
E_{L}\rangle=\frac{1}{n}\sum_{i=1}^{n}E_{L}(\bm{{R}}_{i}),\ \
\bm{{R}}_{i}\sim|\psi_{T}(\bm{{R}})|^{2}.$ (2)

This involves collecting $n$ samples of configurations and averaging over the
so-called local energies

$E_{L}(\bm{{R}})=\frac{1}{\psi_{T}(\bm{{R}})}\hat{\mathcal{H}}\psi_{T}(\bm{{R}}).$
(3)

We apply the variational Monte Carlo (VMC) method to various circular quantum
dots systems. These are systems of interacting electrons confined to move in a
two-dimensional harmonic oscillator well. The (scaled)111Natural units are
used with energy given in units of $\hbar$ and length given in units of
$\sqrt{\hbar/m}$ Hamiltonian is given by

$\hat{\mathcal{H}}=\frac{1}{2}\sum_{i}\left[-\nabla_{i}^{2}+\omega^{2}r_{i}^{2}+\sum_{j\neq
i}\frac{1}{r_{ij}}\right],$ (4)

where $\omega$ is the oscillator frequency, $r_{i}$ is the distance between
electron $i$ and the origin, and $r_{ij}$ is the distance between electrons
$i$ and $j$. We will henceforth assume the total number of electrons $N$ to be
even and the total spin of the system to be zero.

A simple ansatz can be built starting from the analytical solutions to the
non-interacting case. The harmonic oscillator eigenfunctions are given by

$\phi_{m,n}(x,y)\propto
e^{-\omega(x^{2}+y^{2})}H_{m}(\sqrt{\omega}x)H_{n}(\sqrt{\omega}y),$ (5)

where $H_{n}$ are the Hermite polynomials of degree $n$. To constrain the
antisymmetry of the many-body wave function, products of the lowest $N/2$
spatial states and the two spin states $\xi_{\pm}(\sigma)$ are used as a basis
for a Slater determinant

$\psi_{\text{SD}}(\bm{{R}})=\det\Big{[}\Big{\\{}\phi_{m,n}(x_{i},y_{i})\xi_{k}(\sigma_{i})\Big{\\}}\Big{]},$

where $m,n,k$ label the single-particle state, $i$ labels the particle, and
$\bm{{R}}$ contains all coordinates of the $N$ particles. As an aside, we do
not include the spin projections $\sigma_{i}$ as explicit inputs to the wave
function as we will describe how to treat them separately in Section II.B. We
then define a reference state by pulling the common exponential term out of
the determinant and inserting a single variational parameter $\alpha$

$\psi_{\text{Ref}}(\bm{R};\alpha)=e^{-\alpha\omega\sum_{i}(x_{i}^{2}+y_{i}^{2})}\det\Big{[}\Big{\\{}H_{m}(\sqrt{\omega}x_{i})H_{n}(\sqrt{\omega}y_{i})\xi_{k}(\sigma_{i})\Big{\\}}\Big{]}.$
(6)

Correlations among electrons can be handled by a Padé-Jastrow factor Drummond
_et al._ (2004),

$g(\bm{{R}};\beta)=\exp\left(\sum_{i=1}^{N}\sum_{j>i}^{N}\frac{a_{ij}r_{ij}}{1+\beta
r_{ij}}\right),$ (7)

where $\beta$ is a variational parameter and

$a_{ij}=\begin{cases}1/3\quad&\text{if }\sigma_{i}=\sigma_{j}\\\
1\quad&\text{if }\sigma_{i}\neq\sigma_{j}\end{cases},$

in order for the Kato cusp condition to be satisfied Huang  _et al._ (1998).
The product of the Slater determinant and the Padé-Jastrow factor is commonly
named the Slater-Jastrow ansatz,

$\psi_{\text{Slater-
Jastrow}}(\bm{R};\alpha,\beta)=\psi_{\text{Ref}}(\bm{R};\alpha)\times
g(\bm{{R}};\beta).$ (8)

###  II.1 Gaussian-binary restricted Boltzmann machine

There are many possible choices for a machine learning inspired wave function,
but using an artificial neural network is natural. Inspired by Ref. Carleo and
Troyer (2017), our choice is to start from a restricted Boltzmann machine
(RBM) configured for continuous inputs, illustrated in Fig. 1. The inputs
$\bm{{x}}\in\mathbb{R}^{2N}$ are the flattened particle positions and
interactions between the particles are mediated by $H$ hidden binary nodes.
After summing over all the possible values of the hidden nodes, the marginal
distribution of the inputs to the Gaussian-binary RBM takes the form

$\displaystyle P(\bm{{R}};\bm{{a}},\bm{{b}},\bm{{w}})=$
$\displaystyle\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma_{i}^{2}}\right)\prod_{j=1}^{H}\left[1+\exp\left(b_{j}+\sum_{i=1}^{2N}\frac{x_{i}w_{ij}}{\sigma_{i}^{2}}\right)\right].$
(9)

Here, $\bm{{a}}\in\mathbb{R}^{2N}$ and $\bm{{b}}\in\mathbb{R}^{H}$ are the
bias parameters of the input and hidden nodes, respectively. The weights
between the input and hidden nodes are $\bm{{w}}\in\mathbb{R}^{2N\times H}$,
while $\bm{{\sigma}}\in\mathbb{R}^{2N}$ are the widths of the Gaussian input
nodes (not to be confused with the spin projections). It is possible to train
these widths by reparameterizing them as $\sigma_{i}=\exp(s_{i})$, but in this
work all of the widths were fixed to $\sigma=1/\sqrt{\omega}$ and only the
biases and weights are treated as variational parameters. See Appendix A.1 for
the derivation of the marginal probability.

Figure 1: Architecture of a restricted Boltzmann machine. Inter-layer
connections between the visible and the hidden layer are represented by the
black lines, where, for instance, the line connecting $x_{1}$ to $h_{1}$
represents the weight $w_{11}$. The red lines represent the visible biases,
where the line going from the bias unit to the visible unit $x_{3}$ represents
the bias weight $a_{3}$. The purple lines represent the hidden biases, where
the line going from the bias unit to the hidden unit $h_{3}$ represents the
bias weight $b_{3}$.

Notice how the marginal distribution in Eq. (9) mimics the Gaussian parts of
our aforementioned ansätze in Eqs. (6) and (8). Based on such observations,
our next step to is construct two corresponding ansätze

$\psi_{\text{RBM}}(\bm{{R}};\bm{{a}},\bm{{b}},\bm{{w}})=P(\bm{{R}};\bm{{a}},\bm{{b}},\bm{{w}})\times\det\Big{[}\Big{\\{}H_{m}(\sqrt{\omega}x_{i})H_{n}(\sqrt{\omega}y_{i})\xi_{k}(\sigma_{i})\Big{\\}}\Big{]},$
(10)

and

$\psi_{\text{RBM+PJ}}(\bm{{R}};\bm{{a}},\bm{{b}},\bm{{w}},\beta)=P(\bm{{R}};\bm{{a}},\bm{{b}},\bm{{w}})\times
g(\bm{{R}};\beta)\times\det\Big{[}\Big{\\{}H_{m}(\sqrt{\omega}x_{i})H_{n}(\sqrt{\omega}y_{i})\xi_{k}(\sigma_{i})\Big{\\}}\Big{]}.$
(11)

The two trial wave functions above apply different levels of physical
intuition. While $\psi_{\text{RBM}}$ does not contain specific information
about the electron-electron interactions, $\psi_{\text{RBM+PJ}}$ contains a
correlation factor that explicitly upholds the cusp condition. Both ansätze
contain knowledge about the required antisymmetry and the Gaussians in the
marginal distribution help localize the wave functions to satisfy the boundary
conditions far from the oscillator well. Also, as the marginal distribution is
positive definite, these ansätze will never collapse into the bosonic state
even if the marginal distribution is not symmetric.

###  II.2 Code optimization

Parallel computing is an important part of our efforts for developing an
efficient VMC solver. However, increasing the available computational
resources alone is often not sufficient. One should also consider developing
sophisticated algorithms that deliberately minimize the number of floating
point operations, cache misses, and communication between parallel processes.

The kinetic energy term of the Schrödinger equation is usually one of the more
computationally expensive parts to compute in terms of computing cycles. It
includes computing, amongst other elements, the Laplacian of the wave
function. The Laplacian term in the expression for the local energy can be
written as

$\frac{\nabla_{i}^{2}\psi_{T}}{\psi_{T}}=\nabla_{i}^{2}\ln\psi_{T}+\left(\nabla_{i}\ln\psi_{T}\right)^{2}.$

This way of writing the kinetic energy term is beneficial for two reasons:
First, our trial wave function has an exponential shape, which is taken care
of by the log-function. This is often the case for many other ansatzes for the
trial functions. Second, this form allows for separating various elements of
the trial wave function. By writing the trial wave function as a product of
various terms, here $\psi_{T}=\prod_{j}\psi_{j}$, the kinetic energy terms
from each particle $i$ can be written as a sum of their corresponding
Laplacians and gradients

$\frac{\nabla_{i}^{2}\psi_{T}}{\psi_{T}}=\sum_{j}\nabla_{i}^{2}\ln\psi_{j}+\left(\sum_{j}\nabla_{i}\ln\psi_{j}\right)^{2}.$

Obtaining analytical expressions of the gradient and Laplacian for all the
wave function elements is usually computationally advantageous. However, in
many Monte Carlo studies they are normally evaluated numerically using
automatic differentiation Neidinger (2010); Baydin  _et al._ (2018). Nowadays,
automatic differentiation algorithms are employed routinely in VMC
calculations.

The computational complexity of calculating the determinant in the näive way
is proportional to $\mathcal{O}(n^{3})$. This calls for a reduction of
dimensionality as well as efficient evaluations of the Slater determinant. In
this work we do not consider open systems and assume that all single-particle
states up to the Fermi level are filled up. We can then split the Slater
determinant in a spin-up and a spin-down part Pfau  _et al._ (2020) without
affecting the expectation value of the energy, that is

$\psi_{\text{RBM}}=\det\big{[}\\{\phi_{nm}(\bm{{r}}_{\uparrow})\xi(\sigma_{\uparrow})\\}\big{]}\times\det\big{[}\\{\phi_{nm}(\bm{{r}}_{\downarrow})\xi(\sigma_{\downarrow})\\}\big{]}.$

Figure 2: Flow chart of the solver with emphasis on the parallel processing.
The sampling is parallelized across $m$ processes, where the trial wave
function (WF) is broadcast from process 0. Then $E_{L}$,
$\nabla_{\theta}\ln\psi_{T}$ and $E_{L}\cdot\nabla_{\theta}\ln\psi_{T}$ are
sampled independently on each process, and an average is taken after all
processes are done with sampling.

The gradient and Laplacian of the logarithm of a determinant with respect to
particle $i$ are given by

$\displaystyle\nabla_{i}\ln\det=$
$\displaystyle\sum_{j}\nabla_{i}d_{ji}d_{ij}^{-1},$
$\displaystyle\nabla_{i}^{2}\ln\det=$
$\displaystyle\sum_{j}\nabla_{i}^{2}d_{ji}d_{ij}^{-1},$

where $d_{ji}$ is element $(j,i)$ of the matrix and $d_{ij}^{-1}$ is the
corresponding element of the inverse matrix Hammond  _et al._ (1994). The
general solution requires inversion of the matrix, which is known to be
costly. Fortunately, we can use a trick to reduce the cost: If we move one
particle at the time in our sampling over configurations, this means that we
change either the elements of one row or alternatively one column of the
Slater matrix. In this case, there is a simple relation between the old and
the new inverse matrix

$d_{kj}^{-1}=\begin{cases}\frac{1}{R_{i}}d_{kj}^{-1}\quad&\text{if }j=i\\\
d_{kj}^{-1}-\frac{S_{ij}}{R_{i}}d_{ki}^{-1}\quad&\text{if }j\neq
i\end{cases},$

such that the new inverse Slater matrix can be found by a few operations when
the previous is known. Here, $R_{i}$ is the ratio between the new and the old
determinant and $S_{ij}$ is the cross product between columns in the new rows
and the old matrix,

$\displaystyle R_{i}=$ $\displaystyle\sum_{j}d_{ij}d_{ji}^{-1},$
$\displaystyle S_{ij}=$ $\displaystyle\sum_{l}d_{il}d_{lj}^{-1}.$

By using these expressions, the entire Slater determinant matrix is inverted
only once per simulation. We avoid also including spin flips in the
simulations. We limit ourselves in this work to systems where we can use these
approximations. For systems where all single-particle states up to the Fermi
level are filled, the above serves as a useful approximation if the
Hamiltonian does not contain spin-dependent terms. If not, every suggested
move should include possible spin flips as well.

Like most other Monte Carlo schemes, the algorithm can be split into smaller
individual parts and run efficiently in parallel. For each optimization step,
the system is sampled independently in several processes, and the results from
all the processes are averaged before performing the parameter optimization.
In this way, we achieve near perfect parallelization with message parsing
interface (MPI). A flow chart of the simulation code can be found in Fig. 2.
Notice that the parameters are updated with respect to the gradients of the
expectation value of the local energy. The latter is given by

$\nabla_{\theta}\langle E_{L}\rangle=2\left(\langle
E_{L}\nabla_{\theta}\ln\psi_{T}(\bm{{\theta}})\rangle-\langle
E_{L}\rangle\langle\nabla_{\theta}\ln\psi_{T}(\bm{{\theta}})\rangle\right),$

as discussed by Ref. Umrigar and Filippi (2005), among others. When the
simulations are run with a burn-in period, each process should have a burn-in
time equal to the burn-in time for a single process. The theoretical
parallelization efficiency is then given by $(t_{\text{burn-
in}}+t_{\text{sample}})/(mt_{\text{burn-in}}+t_{\text{sample}})$ where $m$ is
the number of parallel processes, $t_{\text{burn-in}}$ is the burn-in time and
$t_{\text{sample}}$ is the total sampling time. Additionally, the weight
optimization can not be parallelized, but has a negligible computational cost
compared to the sampling. The communication can also be neglected even with
low communication speed.

##  III Results and discussions

In this section we compare the computational complexity, ground state energy,
energy convergence, contribution from the different Hamiltonian terms and
electron densities for various trial wave function ansätze. The RBM ansatz
consists of a Slater determinant with Hermite polynomials as the basis,
multiplied with the RBM marginal distribution, as presented in Eq. 10. We have
not made any attempt to include optimized single-particle state functions
through mean-field optimizations like Hartree-Fock theory. The second ansatz
is the RBM ansatz with a correlation factor described in Eq. 7. We have also
include results obtained using a traditional Slater-Jastrow ansatz (Eq. 8),
and as reference we use a plain Slater determinant (without a correlation
factor, Eq. 6). The diffusion Monte Carlo (DMC) results obtained by Høgberget
(2013) are our main reference for the ground state energy. For quantum dots
with two electrons we compare our results with corresponding analytical ones
from Ref. Taut (1994).

In Fig. 3 the CPU time per iteration is plotted as a function of system size
for the various ansätze. The time per iteration differs strongly for the
various ansätze, with the RBM+PJ and Slater-Jastrow ansätze as the most
computationally intensive ones. This is caused by the Padé-Jastrow factor,
which requires all the electron-electron distances at all states. The
execution time is not only longer per iteration, but large systems require in
general lower learning rates, and thus significantly more iterations to
converge.

Figure 3: CPU time per iteration as a function of the number of electron $N$
in a quantum dot well. Here we have employed $10^{4}$ Monte Carlo cycles per
iteration (with $10^{3}$ iterations in total to reach an acceptable
statistical error). The number of hidden nodes in the Boltzmann machine was
set to $H=6$.

The computational overhead associated with parallel processing is negligible
compared to other parts of the code, when ignoring the burn-in steps. In fact,
the variation in CPU time as a result of noise is much more important than the
variation due to parallel processing.

Table 1: The ground state energy of two-dimensional quantum dots with $N$
electrons and frequency $\omega$. Other references include diffusion Monte
Carlo results taken from Høgberget (2013) and semi-analytical results obtained
by Taut (1994). The energy is given in Hartree units, and the numbers in
parenthesis are the statistical uncertainties in the last digit. Bold values
correspond to the lowest ground-state energy obtained in this work. For
abbreviations see the text.
$N$ | $\omega$ | RBM | RBM+PJ | Slater-Jastrow | Ref. |  DMC (Ref. Høgberget (2013)) |  Semi-analytical (Ref. Taut (1994))
---|---|---|---|---|---|---|---
2 | 0.1 | 0.46774(5) | 0.440975(8) | 0.44129(1) | 0.5279(1) | 0.44079(1) |
| 1/6 | 0.70389(7) | 0.666665(1) | 0.66710(1) | 0.7693(1) |  | 2/3
| 0.28 | 1.07100(6) | 1.021668(7) | 1.02192(1) | 1.1388(1) | 1.02164(1) |
| 0.5 | 1.72343(7) | 1.659637(6) | 1.65974(1) | 1.7983(3) | 1.65977(1) |
| 1.0 | 3.0789(1) | 2.999587(5) | 2.99936(1) | 3.1484(3) | 3.00000(1) | 3.0
6 | 0.1 | 3.6971(1) | 3.5700(2) | 3.5695(1) | 3.8552(5) | 3.55385(5) |
| 0.28 | 7.9318(3) | 7.6203(2) | 7.6219(1) | 8.0517(9) | 7.60019(6) |
| 0.5 | 12.2640(6) | 11.80494(7) | 11.8104(2) | 12.2799(9) | 11.78484(6) |
| 1.0 | 20.5635(6) | 20.1773(1) | 20.1918(2) | 20.697(1) | 20.15932(8) |
12 | 0.1 | 12.6772(4) | 12.3416(4) | 12.29962(9) | 12.9742(9) | 12.26984(8) |
| 0.28 | 26.389(2) | 25.7266(2) | 25.7049(4) | 26.625(2) | 25.63577(9) |
| 0.5 | 40.375(1) | 39.2348(2) | 39.2421(5) | 40.227(2) | 39.1596(1) |
| 1.0 | 67.620(3) | 65.7911(7) | 65.7026(4) | 66.744(3) | 65.7001(1) |
20 | 0.1 | 30.7906(8) | 30.1444(2) | 30.0403(2) | 31.253(2) | 29.9779(1) |
| 0.28 | 63.592(3) | 62.1445(5) | 62.0755(7) | 63.681(3) | 61.9268(1) |
| 0.5 | 96.356(5) | 94.101(1) | 94.0433(9) | 95.755(4) | 93.8752(1) |
| 1.0 | 159.428(3) | 156.104(1) | 155.8900(4) | 157.904(6) | 155.8822(1) |
30 | 0.1 | 61.853(2) | 60.774(2) | 60.585(1) | 62.449(4) | 60.4205(2) |
| 0.28 | 126.891(8) | 124.437(2) | 124.195(2) | 126.717(5) | 123.9683(2) |
| 0.5 | 191.455(9) | 187.488(2) | 187.325(3) | 189.977(6) | 187.0426(2) |
| 1.0 | 315.364(8) | 308.989(2) | 308.576(1) | 311.70(2) | 308.5627(2) |
42 | 0.1 | 109.767(7) | 108.128(2) | 107.928(2) | 110.630(7) | 107.6389(2) |
| 0.28 | 224.257(9) | 220.588(3) | 220.224(2) | 223.837(8) | 219.8426(2) |
| 0.5 | 337.43(1) | 331.410(3) | 331.276(3) | 335.18(1) | 330.6306(2) |
| 1.0 | 553.40(1) | 543.746(3) | 542.977(2) | 548.07(2) | 542.9428(8) |
56 | 0.1 | 179.035(8) | 176.659(2) | 176.221(1) | 180.08(1) | 175.9553(7) |
| 0.28 | 364.52(2) | 359.456(6) | 358.470(2) | 363.81(1) | 358.145(2) |
| 0.5 | 547.20(3) | 538.666(5) | 537.841(4) | 544.12(3) | 537.353(2) |
| 1.0 | 894.12(2) | 881.010(5) | 879.514(3) | 887.20(5) | 879.3986(6) |
72 | 0.1 | 274.12(1) | 270.870(3) | 270.296(3) | 275.34(2) |  |
| 0.28 | 556.63(2) | 549.899(8) | 548.315(4) | 555.45(2) |  |
| 0.5 | 833.85(3) | 822.78(2) | 821.089(6) | 829.31(3) |  |
| 1.0 | 1355.37(2) | 1341.54(2) | 1339.85(1) | 1349.65(6) |  |
90 | 0.1 | 399.84(1) | 395.486(4) | 394.621(4) | 401.19(2) |  |
| 0.28 | 809.99(2) | 800.504(6) | 799.187(5) | 808.35(2) |  |
| 0.5 | 1211.92(5) | 1198.12(3) | 1195.025(9) | 1205.65(4) |  |
| 1.0 | 1973.95(5) | 1948.75(2) | 1946.27(1) | 1958.58(5) |  |

The ground state energies of two-dimensional quantum dots of various sizes and
frequencies are presented in table 1. The RBM+PJ and Slater-Jastrow ansätze
provide the lowest ground state energies as expected as the correlation factor
does explicit satisfy the Kato’s cusp conditions. The relative errors with
respect to diffusion Monte Carlo (DMC) calculations tend to be less than 0.2%
for both methods. The various results obtained with the RBM+PJ ansatz show
that this ansatz dominates for small quantum dots, but is outperformed by the
Slater-Jastrow ansatz for larger systems. We suspect this is due to the fact
that the former ansatz is more complex and contains significantly more
parameters than the latter, and has therefore a hard time finding the global
minimum. The optimization could be improved with a more sophisticated
optimization algorithm.

For low frequency dots ($\omega<0.28$), the RBM produces ground state energies
lower than the reference energy, but fails for large high frequency dots
($N>6$, $\omega>0.28$). For low frequency dots the interaction energy
dominates, and the RBM manages to learn the correlations. For high frequency
dots, the interaction energy is less important, and the neural network gets
redundant. If not all weights are set to 0, the neural network will be more a
burden for this case. Additionally, the Gaussian part in the RBM is not
identical to the one in the Hermite functions in the sense that not all
coordinates are multiplied with the same parameter, but rather subtracted with
individual parameters. Since the interaction part is less important for these
systems, the wave function is closer to the non-interacting one. We also
observe that the reference values are similar to corresponding Hartree-Fock
energies up to 30 particles Mariadason (2018), which is expected as both
approaches try to find the optimal single Slater determinant without a
correlation factor.

Some simulations were also performed for the RBM and RBM+PJ ansätze with
sorted network inputs to enforce anti-symmetry under exchange of two
particles, as suggested by Ref. Saito (2018), among others. Sorted inputs
showed promising results, where the ground state energy typically dropped in
the fourth or fifth digit. For example, the RBM+PJ ansatz found the ground
state energy of a system with $N=30$ electrons and $\omega=0.5$ to be
187.311(1) Hartree. This is lower than the Slater-Jastrow energy. We will
investigate this further in a follow-up paper Kim  _et al._ .

Figure 4: Energy convergence for quantum dots with $N=2$ and $\omega=1/6$. The
figure shows how the ground state energy approaches the analytical value (Ref.
Taut (1994)) for various ansätze. Figure 5: Energy convergence for quantum
dots with $N=56$ and $\omega=0.1$. The reference value was obtained by Ref.
Høgberget (2013) using diffusion Monte Carlo. The labeling of the various
results is the same as in the previous figure.

For the traditional ansätze that do not contain artificial neural networks,
the learning rate determines to a large extent how fast the trial wave
function converges towards the true ground state wave function. However, for
the RBM and the RBM+PJ approaches, the results are sensitive to the chosen
learning rate values. A too large learning rate can easily cause exploding
energies, and with a too small learning rate the obtained energies with a
given trial function might not converge at all. Our strategy is to find the
highest learning rate that does not lead to exploding energies. In general
this requires efficient grid searching methods. Here, when the energies have
converged, we decrease the learning rate by a factor of 10 and let it run
until it converges again. Also, the RBMs tend to converge step-wise, making it
hard to know whether or not they have converged.

In Fig. 4, the convergence of the local energy is plotted for our four ansätze
for $N=2$ electrons. Here the results are compared with analytical
calculations for $N=2$. A similar behavior is shown for $N=56$ electrons and
$\omega=0.1$ in Fig. 5. For the $N=2$ electrons case, RBM+PJ turns out to be
the most accurate ansatz with small absolute errors. During training, the
Padé-Jastrow parameter is rather constant, but the adjustment of weights seems
important for the model to reach the energy minimum. For larger systems, the
Slater-Jastrow ansatz provides a slightly lower energy than the RBM+PJ ansatz.
Since the number of inputs for the RBM increases with number of particles, the
network contains additional training parameters as we increase the number of
electrons. Therefore, the networks get more complex for large quantum dots,
and they are naturally harder to train.

Figure 6: Electron density profiles, $\rho(x,y)$, for two-dimensional quantum
dots with frequency $\omega=0.5$ and $N=2$, 6 and 12 electrons seen from the
top. The surface plot and the contour plot on the $xy$-plane illustrate the
density, and the graph on the $yz$-plane represents the cross-section through
$x=0$. They were obtained using RBM (left column) and RBM+PJ (right column),
with $M=2^{30}$ Monte Carlo cycles after convergence. The plots are noise-
reduced using a Savitzky-Golay filter. For abbreviations and description of
the natural units used, see the main text for more details. Figure 7: Electron
density profiles, $\rho(x,y)$, for two-dimensional quantum dots with frequency
$\omega=0.5$ and $N=20$, 30 and 42 electrons seen from the top. The surface
plot and the contour plot on the xy-plane illustrate the density, and the
graph on the yz-plane represents the cross-section through $x=0$. They were
obtained using RBM (left column) and RBM+PJ (right column), with $M=2^{30}$
Monte Carlo cycles after convergence. The plots are noise-reduced using a
Savitzky-Golay filter. For abbreviations and description of the natural units
used, see main text.

The spatial one-body density plots for $\omega=1.0$ and the RBM and RBM+PJ
ansätze are presented in Figs. 6,7. The electron densities have a wave shape
for both ansätze, with two nodes for $N=2$, three nodes for $N=6$ and so on
like observed by Refs. Ghosal  _et al._ (2007); Høgberget (2013). The RBM
seems to exaggerate the states with more distinct peaks (higher peaks and
lower wave valleys) compared to the RBM+PJ results. This can be explained by
more localized electrons, which becomes even more apparent in low frequency
dots (see Fig. 8), as the interactions are modelled differently. The RBM would
hardly be able to model the correct electron-electron distances, as the
network itself is purely linear. The electron density was found to be shape-
invariant for high frequency dots ($\omega>0.28$), with decreasing spatial
range as the frequency increases.

Figure 8: One-body density profile, $\rho(x,y)$, of two-dimensional quantum
dots with frequency $\omega=0.1$ and $N=2$, 6 and 12 electrons seen from left
to right. The ansätze used are RBM (upper panels) and RBM+PJ (lower panels).
The surface plot and the contour plot on the xy-plane illustrate the density,
and the graph on the yz-plane represents the cross-section through $x=0$. The
surface plots are noise-reduced using a Savitzky-Golay filter. For
abbreviations and description of the natural units used, see main text. Figure
9: One-body density profile, $\rho(x,y)$, of two-dimensional quantum dots with
frequency $\omega=0.01$ and $N=2$ electrons. The ansätze used are RBM (left)
and RBM+PJ (right). The surface plot and the contour plot on the xy-plane
illustrate the density, and the graph on the yz-plane represents the cross-
section through $x=0$. The surface plots are noise-reduced using a Savitzky-
Golay filter. For abbreviations and description of the natural units used, see
main text.

When reducing the frequency further, the interaction energy dominates over the
kinetic energy and harmonic oscillator energy (see Figs. 10 and 11), and the
electrons naturally become more localized. The Padé-Jastrow factor solves this
by radially localizing the electrons and conserving the circular symmetry,
such that the electrons are confined to specific orbitals. This can be seen
from Fig. 8, where the spatial one-body density is plotted for low frequency
dots ($\omega=0.1$) and system sizes $N=2$, 6 and 12. Radial localization is
also what we would expect from the Hamiltonian, which is strictly circular
symmetric. On the other hand, the RBM seems to localize the electrons both in
radial and angular direction, with the number of electrons corresponding to
the number of peaks. This is a nonphysical solution to the problem, and shows
that the RBM ansatz breaks down for low frequencies. The RBM+PJ ansatz, on the
other hand, confines the electrons in orbitals. Distinct peaks in radial
direction is what is expected from Wigner crystallization, which we might see
indications of with density parameters $r_{s}\approx 6.7$, $r_{s}\approx 1.2$
and $r_{s}=0.3$ respectively for the three system sizes with the RBM+PJ
ansatz. Notice for instance the small “pit” on top of the $N=2$ plot for
RBM+PJ, which is not seen for higher oscillator frequencies. In Fig. 9, we
have decreased the one-body density even further to $\omega=0.01$ for $N=2$.
As expected, the electrons become even more localized and are clearly showing
Wigner crystallization effects with a density parameter of $r_{s}\approx 29$.
For the RBM ansatz, the electrons are strongly localized (in all directions)
and the electron densities barely overlap. For the RBM+PJ ansatz, we observe
strong orbital confinement where the Wigner crystallization was not the target
of this study, but the framework seems capable of a profound study of this
phenomenon.

Figure 10: Energy distribution for $N=2$ electrons for the RBM+PJ ansatz, with
frequencies ranging from $\omega=0.01$ to $10.0$. Figure 11: Energy
distribution for $N=20$ electrons for the RBM+PJ ansatz, with frequencies
ranging from $\omega=0.01$ to $10.0$.

To understand the behavior of the one-body density for various frequencies, we
investigate the expectation values of the kinetic energy
($\langle\hat{T}\rangle$), the harmonic oscillator potential energy
($\langle\hat{V}_{\text{ext}}\rangle$) and the two-body interaction energy
($\langle\hat{V}_{\text{ext}}\rangle$). The energy distribution is plotted for
the RBM+PJ ansätze for $N=2$ and $N=20$ (Figs. 10 and 11) for frequencies
$\omega\in\\{0.01$, 0.1, 0.28, 0.5, 1.0, 2.0, 3.0, 5.0 and $10.0\\}$. For
large values of the oscillator frequency it is the one-body part of the
Hamiltonian which dominates in absolute value (kinetic energy and harmonic
oscillator potential energy) compared with the expectation value of the two-
body interaction. For such frequencies we notice also that the results can
almost be interpreted in terms of the virial theorem. This theorem provides a
useful relation between the kinetic and potential energy Fock (1930). For
circular quantum dots without the two-electron interaction, the theorem reads
$2\langle\hat{T}\rangle=2\langle\hat{V}_{\text{ext}}\rangle$. From our results
we notice that, due to the interaction energy, the ratio between kinetic and
harmonic oscillator energies are not exactly equal two for large frequencies.
However, for such frequencies the interaction energy plays a less prominent
role, resulting thus in a ratio which is close to the ideal value. When we
decrease the frequency however, and the system becomes more dilute with an
increase of the mean distance between particles, one notices a somewhat
counter intuitive behavior. The expectation value of the kinetic energy and
the harmonic oscillator potential energy decrease (and the electrons become
more localized as seen in the one-body densities above). However, many-body
correlations increase in importance with decreasing frequency. This is
reflected in the increased role of the expectation value of the two-body
Coulomb interaction, an effect which is simply due to the infinite range of
the Coulomb interaction. If we were to multiply the Coulomb interaction with a
finite range factor, this effect would disappear. Figures 10 and 11 show this
behavior rather clearly.

##  IV Computational details

The quantum dot systems are studied using a general framework for variational
Monte Carlo simulations (The code is available on github.com/evenmn/VMaChine).
Importance sampling was used to accelerate the simulations Metropolis  _et
al._ (1953). To minimize the local energy, we applied the Adam optimizer with
$\beta_{1}=0.9$ and $\beta_{2}=0.999$ as suggested by Ref. Kingma and Ba
(2014). We use adaptive number of cycles, starting from $2^{20}$ and then
increased to $2^{24}$ for the ten last iterations after the energy had
converged, and further to $2^{30}$ for the very last iteration to reduce the
statistical uncertainty and noise in electron density. The particle step
length was chosen to get an acceptance ratio close to 99.5%, spawning from
$10^{1}$ for the smallest and weakest systems to $10^{-3}$ for the large and
narrow oscillators. The optimal learning rate was found by grid searches, and
vary from $10^{1}$ to $10^{-5}$. Both the step length and the learning rate
depend strongly on the system and the trial wave function.

For the RBM and RBM+PJ ansätze, only the raw particle positions were input to
the RBM. This choice was made for performance reasons, as inputting processed
positions like the electron-electron distances would lead to significantly
more computations as it would require the processing and increase network
complexity.

The Gaussian parameter was initialized to $\alpha=1.0$, which is the
analytical optimum without interaction. We use Xavier initialization for the
RBM weights, putting they all close to zero Glorot and Bengio (2010). The
special case with all weights set to zero corresponds to our reference ansatz
with $\alpha=1.0$. For all the RBMs but the most narrow ones ($\omega=1.0$),
the number of hidden nodes was set to 6, giving 40 to 1272 parameters for the
various system sizes. For $\omega=1.0$, we used $H=12$ hidden nodes to achieve
a lower energy.

All the simulations were run on CPUs. In total the computational cost of this
project was of the order of $10^{6}$ CPU hours with largest amount of cycles
spent, by obvious reasons, on the largest systems.

##  V Conclusions and Perspectives

We found that the RBM ansatz gives a significantly lower ground state energy
than the reference ansatz for low frequency quantum dots. This may indicate
that the RBM manages to capture some of the electron-electron correlations.
Based on the one-body density plots, the RBM found the electrons to be
localized both angularly and radially compared to the ansätze containing a
Padé-Jastrow factor, which confine electrons radially only. For high frequency
dots, the RBM fails in the sense that the obtained ground state energy is
larger than the reference energy. This can be explained by the fact that when
the interactions get less important, the reference ansatz is a good guess. The
RBM+PJ ansatz gives energies close to the DMC energies for small quantum dots.
The ansatz performance needs to be seen in light of the computational cost, as
the ansätze containing a Padé-Jastrow factor are far more computationally
intensive.

From the one-body density and energy distribution plots, it is apparent that
the RBM ansatz is not able to capture the correlations at the same level as
the Padé-Jastrow factor. Because of the linearity of the network, it is
impossible for it to compute the distance between the particles, which is
crucial to model the interactions correctly. One solution to this could be to
input the electron-electron distance into the network, as discussed by Pfau
_et al._ (2020), but this would increase the computational intensity. Also,
despite including a Slater determinant, the trial wave function is not
necessarily anti-symmetric when including the an RBM. We performed some
simulations where anti-symmetry was forced by sorting the network inputs,
which showed promising results in terms of the ground state energy. However,
although promising results can be obtained, RBMs are less flexible than
general neural networks that make fewer assumptions about the specific
mathematical forms of the trial functions. We have encoded explicitly the
anti-symmetry via a Slater determinant. Furthermore, two-body correlations are
constructed using a Jastrow factor. An RBM with Gaussian distributions is
capable of capturing the one-body part of the problem, but is less flexible in
finding two-body or more complicated many-body correlations. Although the RBM
results reported here are promising compared with existing VMC calculations,
recent results with neural networks like those presented in for example Refs.
Pfau  _et al._ (2020); Cassella  _et al._ (2022); Lovato  _et al._ (2022);
Adams  _et al._ (2021); Carrasquilla and Torlai (2021) offer much more
flexible and promising research venues for deep learning methods applied to
many-body problems. Results obtained with deep neural networks for these
systems will be presented in a future work Kim  _et al._ .

### Acknowledgements

AL and BF are supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under contracts DE-AC02-06CH11357, by the NUCLEI
SciDAC program, and the DOE Early Career Research Program. MHJ is supported by
the U.S. Department of Energy, Office of Science, office of Nuclear Physics
under grant No. DE-SC0021152 and U.S. National Science Foundation Grants No.
PHY-1404159 and PHY-2013047. JK is supported by the U.S. National Science
Foundation Grants No. PHY-1404159 and PHY-2013047. EMN is supported by the
Norwegian Research Council under grant 287084.

##  Appendix A Derivation of RBM distributions

In this appendix, we will derive the marginal and conditional distributions of
a Gaussian-binary restricted Boltzmann machine with the system energy

$E(\bm{{x}},\bm{{h}})=\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma_{i}^{2}}-\sum_{j=1}^{H}b_{j}h_{j}-\sum_{i=1}^{2N}\sum_{j=1}^{H}\frac{x_{i}w_{ij}h_{j}}{\sigma_{i}^{2}}.$

There are $2N$ visible units $x_{i}$ with related bias weights $a_{i}$ and $H$
hidden units $h_{j}$ with related bias weights $b_{j}$. $w_{ij}$ are the
weights connecting the visible units to the hidden units. The joint
probability distribution is given by the Boltzmann distribution

$P(\bm{{x}},\bm{{h}})=\frac{1}{Z}\exp(-\beta E(\bm{{x}},\bm{{h}})),$

where $Z$ is the partition function,

$Z=\iint d\bm{{x}}d\bm{{h}}P(\bm{{x}},\bm{{h}}),$

and $\beta=1/k_{B}T$ is the reciprocal temperature that will be fixed to 1. As
the marginal and conditional distributions are closely related both for the
visible and hidden layer, we present the distributions in sections respective
for the two layers.

###  A.1 Distributions of visible units

The distributions of the visible units are used to find properties related to
the visible units. If we recall a restricted Boltzmann machine, the
transformation between the visible units and the hidden units is
$f_{j}(\bm{{x}};\bm{{\theta}})=b_{j}+\sum_{i=1}^{2N}w_{ij}x_{i}/\sigma_{i}^{2}$.
By this expression, we can express the joint probability distribution as

$\displaystyle P(\bm{{x,\bm{{h}}}})$
$\displaystyle=\frac{1}{Z}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma_{i}^{2}}+\sum_{j=1}^{H}b_{j}h_{j}+\sum_{i=1}^{2N}\sum_{j=1}^{H}\frac{x_{i}w_{ij}h_{j}}{\sigma_{i}^{2}}\right),$
(12)
$\displaystyle=\frac{1}{Z}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma^{2}}+\sum_{j=1}^{H}h_{j}f_{j}(\bm{{x}};\bm{{\theta}})\right).$

The marginal distribution of the visible units is given by the sum over all
possible hidden states, $\\{\bm{h}\\}\in\\{0,1\\}$:

$P(\bm{{x}})=\sum_{\\{\bm{{h}}\\}}P(\bm{{x}},\bm{{h}}),$

as the hidden units can take binary values only. By inserting the expression
of the joint probability distribution from equation (12), we obtain

$\displaystyle P(\bm{{x}})$
$\displaystyle=\frac{1}{Z}\sum_{\\{\bm{{h}}\\}}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma^{2}}+\sum_{j=1}^{H}h_{j}f_{j}(\bm{{x}};\bm{{\theta}})\right),$
$\displaystyle=\frac{1}{Z}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma^{2}}\right)\times\left(\sum_{\\{\bm{{h}}\\}}\prod_{j=1}^{H}\exp(h_{j}f_{j})\right),$
$\displaystyle=\frac{1}{Z}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma^{2}}\right)\prod_{j=1}^{H}\sum_{h_{j}=0}^{1}\exp(h_{j}f_{j}),$
$\displaystyle=\frac{1}{Z}\exp\left(-\sum_{i=1}^{2N}\frac{(x_{i}-a_{i})^{2}}{2\sigma^{2}}\right)\prod_{j=1}^{H}\left[1+\exp(f_{j}(\bm{{x}};\bm{{\theta}}))\right].$

This is what we will use as the marginal distribution of the visible units.

## References

  * Goodfellow  _et al._ (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (The MIT Press, Cambridge, Massachusetts, 2016).
  * Carleo and Troyer (2017) G. Carleo and M. Troyer, Science 355, 602 (2017).
  * Carrasquilla and Torlai (2021) J. Carrasquilla and G. Torlai,  (2021), 10.48550/arxiv.2101.11099.
  * Pfau  _et al._ (2020) D. Pfau, J. S. Spencer, Alexander G. D. G. Matthews, and W. M. C. Foulkes, Physical Review Research 2, 033429 (2020).
  * Calcavecchia  _et al._ (2014) F. Calcavecchia, F. Pederiva, M. H. Kalos, and T. D. Kühne, Physical Review E 90, 053304 (2014).
  * Carleo  _et al._ (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Reviews of Modern Physics 91, 045002 (2019).
  * Boehnlein  _et al._ (2022) A. Boehnlein, M. Diefenthaler, N. Sato, M. Schram, V. Ziegler, C. Fanelli, M. Hjorth-Jensen, T. Horn, M. P. Kuchera, D. Lee, W. Nazarewicz, P. Ostroumov, K. Orginos, A. Poon, X.-N. Wang, A. Scheinker, M. S. Smith, and L.-G. Pang, Reviews of Moddern Physics 94, 031003 (2022).
  * Adams  _et al._ (2021) C. Adams, G. Carleo, A. Lovato, and N. Rocco, Physical Review Letters 127, 022502 (2021).
  * Lovato  _et al._ (2022) A. Lovato, C. Adams, G. Carleo, and N. Rocco,  (2022), 10.48550/arxiv.2206.10021.
  * Hornik  _et al._ (1989) K. Hornik, M. Stinchcombe, and H. White, Neural Networks 2, 359 (1989).
  * Murphy (2012) K. P. Murphy,  _Machine Learning: A Probabilistic Perspective_ (The MIT Press, Cambdridge, Massachusetts, 2012).
  * Hastie  _et al._ (2009) T. Hastie, R. Tibshirani, and J. Friedman,  _The Elements of Statistical Learning: Data Mining, Inference and Prediction_ (Springer Verlag, Berlin, 2009).
  * Bishop (2006) C. M. Bishop,  _Pattern Recognition and Machine Learning_ (Springer Verlag, Berlin, 2006).
  * Cassella  _et al._ (2022) G. Cassella, H. Sutterud, S. Azadi, N. D. Drummond, D. Pfau, J. S. Spencer, and W. M. C. Foulkes,  (2022), 10.48550/arxiv.2202.05183.
  * Note (1) Natural units are used with energy given in units of $\hbar$ and length given in units of $\sqrt{\hbar/m}$.
  * Drummond  _et al._ (2004) N. D. Drummond, M. D. Towler, and R. J. Needs, Physical Review B 70, 235119 (2004).
  * Huang  _et al._ (1998) C.-J. Huang, C. Filippi, and C. J. Umrigar, The Journal of Chemical Physics 108, 8838 (1998).
  * Neidinger (2010) R. D. Neidinger, SIAM Review 52, 545 (2010).
  * Baydin  _et al._ (2018) A. G. Baydin, B. A. Pearlmutter, A. Andreyevich Radul, and J. M. Siskind, Journal of Machine Learning Research 18, 1 (2018).
  * Hammond  _et al._ (1994) B L Hammond, W A Lester, and P J Reynolds,  _Monte Carlo Methods in Ab Initio Quantum Chemistry_ (World Scientific, Singapore, 1994).
  * Umrigar and Filippi (2005) C. J. Umrigar and C. Filippi, Physical Review Letters 94, 150201 (2005).
  * Høgberget (2013) J. Høgberget,  _Quantum Monte Carlo Studies of Generalized Many-body Systems_ , Master’s thesis, University of Oslo (2013).
  * Taut (1994) M. Taut, Journal of Physics A , 1045 (1994).
  * Mariadason (2018) A. A. Mariadason,  _Quantum Many-Body Simulations of Double Dot System_ , Master’s thesis, University of Oslo (2018).
  * Saito (2018) H. Saito, Journal of the Physical Society of Japan 87, 074002 (2018).
  * (26) J.M. Kim, B. Fore, E. M. Nordhagen, A. Lovato, and M. Hjorth-Jensen, “Deep learning and confined electrons in two dimensions,” .
  * Ghosal  _et al._ (2007) A. Ghosal, A. D. Güçlü, C. J. Umrigar, D. Ullmo, and H. U. Baranger, Physical Review B 76, 085341 (2007).
  * Fock (1930) V. Fock, Zeitschrift für Physik 63, 855 (1930).
  * Metropolis  _et al._ (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba,  (2014), 10.48550/arxiv.1412.6980.
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , Proceedings of Machine Learning Research, Vol. 9, edited by Y. W. Teh and M. Titterington (PMLR, 2010) p. 249.


OCU-PHYS 537

NITEP 97

Nonvanishing Finite Scalar Mass
in Flux Compactification

Takuya Hirosea and Nobuhito Marua,b,

aDepartment of Mathematics and Physics, Osaka City University,
Osaka 558-8585, Japan

bNambu Yoichiro Institute of Theoretical and Experimental Physics (NITEP),
Osaka City University, Osaka 558-8585, Japan

We study possibilities to realize a nonvanishing finite Wilson line (WL)
scalar mass in flux compactification. Generalizing loop integrals in the
quantum correction to WL mass at one-loop, we derive the conditions for the
loop integrals and mode sums in one-loop corrections to WL scalar mass to be
finite. We further guess and classify the four-point and three-point
interaction terms satisfying these conditions. As an illustration, the
nonvanishing finite WL scalar mass is explicitly shown in a six dimensional
scalar QED by diagrammatic computation and effective potential analysis. This
is the first example of finite WL scalar mass in flux compactification.

##  1 Introduction

It has been considered that the hierarchy problem is one of the guiding
principles to search for the physics beyond the Standard Model (SM) of
particle physics. In the SM, the quantum correction to the mass of Higgs field
is sensitive to the square of the ultraviolet cutoff scale of the theory (for
example, Planck scale or the scale of grand unified theory). Since an
experimental values of the Higgs mass is 125 GeV, the solution of the
hierarchy problem requires an unnatural fine-tuning of parameters or exploring
a new physics beyond the SM at order of TeV scale. Although the latter
approaches have been mainly studied so far, no signature of new physics has
been found, which is likely to increase the new physics scale, namely Higgs
mass. Therefore, the solution seems to be desirable such that the Higgs mass
is zero at the classical level and is generated by the quantum effects.

As one of the approaches to the hierarchy problem, higher dimensional theory
with magnetic flux compactification has been studied. Magnetic flux
compactification has been originally studied in string theory [1, 2] Even in
the field theories, flux compactification has many attractive properties:
attempt to explain the number of the generations of the SM fermion [3],
computation of Yukawa coupling [4, 5, 6]. Recently, it has been shown that the
quantum corrections to the masses of zero-mode of the scalar field induced
from extra component of higher dimensional gauge field (called Wilson-line
(WL) scalar field) are canceled [7, 8, 9, 10, 11, 12]. The physical reason of
the cancellation is that the shift symmetry from translation in extra spaces
forbids the mass term of scalar field. In that situation, the zero-mode of the
scalar field can be identified with Nambu-Goldstone (NG) boson of
spontaneously broken translational symmetry. It is not possible for these
results to apply to the hierarchy problem as it stands since the scalar field
is also massless at quantum level. If we identify Wilson-line scalar field
with Higgs field, we need some mechanism to generate an explicit breaking term
of the translational symmetry in compactified space and the scalar field must
be a pseudo NG boson such as pion.

In this paper, we study the possibility to realize nonvanishing finite WL
scalar mass in flux compactification. First, we generalize loop integrals in
the quantum correction to WL scalar mass at one-loop. Then, the conditions for
the loop integral and mode sum to be finite are derived. We further guess and
classify the four-point and three-point interaction terms generating the
finite one-loop quantum correction to WL scalar mass. Of these interaction
terms, we focus on a simplest interaction term and illustrate the finite
quantum correction to the WL scalar mass in a six dimensional scalar QED in
two ways: diagrammatic computation or effective potential analysis. This is
the first example of finite WL scalar mass in flux compactification.

This paper is organized as follows. We introduce a six-dimensional theory with
flux compactification and derive Kaluza-Klein mass spectrum of scalar field,
fermion field and SU(2) gauge field in section 2. In section 3, we generalize
the loop integrals of quantum correction to the masses of Wilson-line scalar
field. After deriving the conditions for the quantum correction to be finite,
the interaction terms providing finite quantum corrections are classified. In
section 4, we focus on an interaction term of all interaction terms classified
in section 3 and calculate finite quantum correction to WL scalar mass in a
six dimensional scalar QED. In the last section, we summarize our conclusion.
The property of Hurwitz zeta function is summarized in appendix A.

##  2 Flux Compactification and Kaluza-Klein Mass

In this section, we introduce our setup and summarize the Kaluza-Klein mass
spectrum of various fields, which are required for calculating the quantum
correction to WL scalar mass.

###  2.1 Flux compactification

Let us first consider a six-dimensional U(1) gauge theory with a constant
magnetic flux. The six-dimensional spacetime is a product of four-dimensional
Minkowski spacetime $M^{4}$ and two-dimensional torus $T^{2}$. For later
discussion, let us consider the following Lagrangian

$\displaystyle\mathcal{L}_{gauge}$ $\displaystyle=-\frac{1}{4}F_{MN}F^{MN},$
(1) $\displaystyle\mathcal{L}_{scalar}$
$\displaystyle=-D_{M}\overline{\Phi}D^{M}\Phi,$ (2)
$\displaystyle\mathcal{L}_{fermion}$
$\displaystyle=i\overline{\Psi}\Gamma^{M}D_{M}\Psi,$ (3)

where the spacetime index is given by
$M,N=0,1,\cdots,6,~{}\mu,\nu=0,1,2,3,~{}m,n,=5,6$ respectively and we follow
the metric convention as $\eta_{MN}=(-1,+1,\cdots+1)$. The field strength and
the covariant derivative of U(1) gauge field $A_{M}$ are defined by
$F_{MN}=\partial_{M}A_{N}-\partial_{N}A_{M}$, $D_{M}=\partial_{M}-igA_{M}$
with a gauge coupling constant $g$. $\Phi$ is a bulk scalar field.
$\Gamma^{M}$ are six-dimensional gamma-matrices. 111For more detail
convention, see [9, 13]. $\Psi$ is a six-dimensional Weyl spinor with two-
component Weyl spinors, which satisfies $\Gamma^{7}\Psi=-\Psi$ in terms of
$\Gamma^{7}\equiv-\Gamma^{0}\Gamma^{1}\Gamma^{2}\Gamma^{3}\Gamma^{5}\Gamma^{6}$.
$\Psi$ has $\psi$ with a charge $-g$ and $\chi$ with a charge $+g$ and is
expressed by

$\displaystyle\Psi=\left(\begin{array}[]{l}\psi_{L}\\\
\psi_{R}\end{array}\right),~{}~{}\psi_{L}=\left(\begin{array}[]{l}\psi\\\
0\end{array}\right),~{}~{}\psi_{R}=\left(\begin{array}[]{l}0\\\
\bar{\chi}\end{array}\right).$ (10)

We introduce the magnetic flux in our model. The magnetic flux is given by the
nontrivial background (or vacuum expectation value (VEV)) of the fifth and the
sixth component of the gauge field $A_{5,6}$. This background must satisfy
their classical equation of motion $\partial^{m}\braket{F_{mn}}=0$. In flux
compactification, the background of $A_{5,6}$ is chosen as

$\displaystyle\braket{A_{5}}=-\frac{1}{2}fx_{6},~{}~{}~{}\braket{A_{6}}=\frac{1}{2}fx_{5},$
(11)

which introduces a constant magnetic flux density $\braket{F_{56}}=f$ with a
real number $f$. Note that this solution breaks an extra-dimensional
translational invariance spontaneously. Integrating over $T^{2}$, the magnetic
flux is quantized as follows

$\displaystyle\frac{g}{2\pi}\int_{T^{2}}dx_{5}dx_{6}\braket{F_{56}}=\frac{g}{2\pi}L^{2}f=N\in\mathbb{Z},$
(12)

where $L^{2}$ is an area of two-dimensional torus. In the following, we set
$L=1$ for simplicity.

It is useful to define $\partial,z$, and $\phi$ as

$\displaystyle\partial\equiv\partial_{z}=\partial_{5}-i\partial_{6},\quad
z\equiv\frac{1}{2}\left(x_{5}+ix_{6}\right),\quad\phi=\frac{1}{\sqrt{2}}\left(A_{6}+iA_{5}\right).$
(13)

In terms of these complex coordinates and variables, the VEV of $\phi$ is
given by $\braket{\phi}=f\bar{z}/\sqrt{2}$. We expand $\phi$ around the flux
background $\phi=\braket{\phi}+\varphi$, where $\varphi$ is a quantum
fluctuation. To distinguish $\varphi$ from an introduced bulk scalar $\Phi$,
we call $\varphi$ Wilson line (WL) scalar field. Defining the covariant
derivatives in the complex coordinates is also useful to obtain Kaluza-Klein
(KK) masses later, which are defined as

$\displaystyle D$
$\displaystyle=D_{5}-iD_{6}=\partial-\sqrt{2}g\phi=\mathcal{D}-\sqrt{2}g\varphi,$
(14) $\displaystyle\bar{D}$
$\displaystyle=D_{5}+iD_{6}=\bar{\partial}+\sqrt{2}g\bar{\phi}=\bar{\mathcal{D}}+\sqrt{2}g\bar{\varphi},$
(15) $\displaystyle\mathcal{D}$
$\displaystyle=\mathcal{D}_{5}-i\mathcal{D}_{6}=\partial-\sqrt{2}g\braket{\phi},$
(16) $\displaystyle\bar{\mathcal{D}}$
$\displaystyle=\bar{\mathcal{D}}_{5}+i\bar{\mathcal{D}}_{6}=\bar{\partial}+\sqrt{2}g\braket{\bar{\phi}}.$
(17)

Note that $\mathcal{D}_{m}$ means the covariant derivative with VEV.

Finally, we consider an SU(2) Yang-Mills theory with a constant magnetic flux.
The Lagrangian of Yang-Mills theory is given by

$\displaystyle\mathcal{L}_{YM}=-\frac{1}{4}F^{a}_{MN}F^{aMN},$ (18)

where $a=1,2,3$ are gauge indices. The field strength and the covariant
derivative of SU(2) gauge field $A^{a}_{M}$ are defined by
$F^{a}_{MN}=\partial_{M}A^{a}_{N}-\partial_{N}A^{a}_{M}-ig[A_{M},A_{N}]^{a}$,
$D_{M}A^{a}_{N}=\partial_{M}A^{a}_{N}+g\epsilon^{abc}A^{b}_{M}A^{c}_{N}$.
$\epsilon^{abc}$ is a totally anti-symmetric tensor of SU(2). In the case of
Yang-Mills theory, we introduce a flux background as

$\displaystyle\braket{A_{5}^{1}}=-\frac{1}{2}fx_{6},\quad\braket{A_{6}^{1}}=\frac{1}{2}fx_{5},\quad\braket{A_{5}^{2,3}}=\braket{A_{6}^{2,3}}=0.$
(19)

For later convenience, we define the covariant derivatives in the complex
coordinates as follows

$\displaystyle DX^{a}$ $\displaystyle=(D_{5}-iD_{6})X^{a}=\partial
X^{a}-\sqrt{2}g[\phi,X]^{a}=\mathcal{D}X^{a}-\sqrt{2}g[\varphi,X]^{a},$ (20)
$\displaystyle\bar{D}X^{a}$
$\displaystyle=(D_{5}+iD_{6})X^{a}=\bar{\partial}X^{a}+\sqrt{2}g[\bar{\phi},X]^{a}=\bar{\mathcal{D}}X^{a}+\sqrt{2}g[\bar{\varphi},X]^{a},$
(21) $\displaystyle\mathcal{D}X^{a}$
$\displaystyle=(\mathcal{D}_{5}-i\mathcal{D}_{6})X^{a}=\partial
X^{a}-\sqrt{2}g[\braket{\phi},X]^{a},$ (22)
$\displaystyle\bar{\mathcal{D}}X^{a}$
$\displaystyle=(\bar{\mathcal{D}}_{5}+i\bar{\mathcal{D}}_{6})X^{a}=\bar{\partial}X^{a}+\sqrt{2}g[\braket{\bar{\phi}},X]^{a}.$
(23)

###  2.2 Kaluza-Klein mass spectrum

To compute one-loop correction to WL scalar mass in flux compactification, we
need to derive mass eigenvalues for fields propagating in a loop. In analogy
to the quantum mechanics in magnetic field, we regard the covariant derivative
$\mathcal{D}$ and $\bar{\mathcal{D}}$ as creation and annihilation operators
by

$\displaystyle
a=\frac{1}{\sqrt{2gf}}i\bar{\mathcal{D}},~{}~{}~{}a^{\dagger}=\frac{1}{\sqrt{2gf}}i\mathcal{D},$
(24)

which satisfy the commutation relation $[a,a^{\dagger}]=1$. Hereafter, we
denote $\alpha=2gf$.

The ground state mode function is determined by
$a\xi_{0,j}=0,a^{\dagger}\bar{\xi}_{0,j}=0$, where $j=0,\cdots,|N|-1$ accounts
for the degeneracy of the ground state. Creation operator and annihilation
operator acts on mode functions as

$\displaystyle
a^{\dagger}\xi_{n,j}=\sqrt{n+1}\xi_{n+1,j},~{}~{}~{}a\xi_{n,j}=\sqrt{n}\xi_{n-1,j},$
(25)

and we can construct the higher mode function $\xi_{n,j}$ in the same way as
the harmonic oscillator (in detail, see [14])

$\displaystyle\xi_{n,j}=\frac{1}{\sqrt{n!}}(a^{\dagger})^{n}\xi_{0,j},~{}~{}~{}\bar{\xi}_{n,j}=\frac{1}{\sqrt{n!}}(a)^{n}\bar{\xi}_{0,j},$
(26)

where $n=0,1,2\cdots$ is Landau level. The higher mode function satisfies an
orthonormality condition

$\displaystyle\int_{T^{2}}dx^{2}\bar{\xi}_{n^{\prime},j^{\prime}}\xi_{n,j}=\delta_{n,n^{\prime}}\delta_{j,j^{\prime}}.$
(27)

####  2.2.1 Scalar field

We decompose (2) into a four-dimensional part and an extra-dimensional part
(see [7])

$\displaystyle\mathcal{L}_{scalar}=-D_{\mu}\overline{\Phi}D^{\mu}\Phi-
D_{m}\overline{\Phi}D^{m}\Phi.$ (28)

Now, we focus on the second term in (28) and extract mass term

$\displaystyle\mathcal{L}_{scalar~{}mass}$
$\displaystyle=-\mathcal{D}_{m}\overline{\Phi}\mathcal{D}^{m}\Phi$
$\displaystyle=-\overline{\Phi}\alpha\left(a^{\dagger}a+\frac{1}{2}\right)\Phi.$
(29)

Then, the KK mass of scalar field is obtained by

$\displaystyle m^{2}_{scalar}=\alpha\left(n+\frac{1}{2}\right),$ (30)

where the fact that $a^{\dagger}a$ is a number operator is used. Note that the
KK mass of WL scalar $\varphi^{2,3}$ induced from SU(2) Yang-Mills theory
corresponds to (30) in Feynman gauge222If you have a concern about gauge-
fixing term and ghost, see [10].. Although SU(2) Yang-Mills theory involves a
ghost field, the KK mass of the ghost field also agrees with (30) in Feynman
gauge.

####  2.2.2 Fermion field

We decompose (3) as in the case of scalar field (see [9])

$\displaystyle\mathcal{L}_{fermion}=i\overline{\Psi}\Gamma^{\mu}D_{\mu}\Psi+i\overline{\Psi}\Gamma^{5}D_{5}\Psi+i\overline{\Psi}\Gamma^{6}D_{6}\Psi,$
(31)

and we focus on the second and the third terms in (31). Decomposing these
terms in terms of two-component Weyl spinors $\psi$ and $\chi$, the mass terms
of fermion field are expressed by

$\displaystyle\mathcal{L}_{fermion~{}mass}$ $\displaystyle=-\chi(\partial-
gf\bar{z})\psi-\bar{\chi}(\bar{\partial}-gfz)\bar{\psi}.$ (32)

In this case, there are two pairs of annihilation and creation operators

$\displaystyle a_{-}$ $\displaystyle=\frac{i}{\sqrt{\alpha}}(\partial-
gf\bar{z}),~{}~{}~{}a^{\dagger}_{-}=\frac{i}{\sqrt{\alpha}}(\bar{\partial}+gf\bar{z}),$
(33) $\displaystyle a_{+}$
$\displaystyle=\frac{i}{\sqrt{\alpha}}(\bar{\partial}-gf\bar{z}),~{}~{}~{}a^{\dagger}_{+}=\frac{i}{\sqrt{\alpha}}(\partial+gf\bar{z}),$
(34)

where $a_{-},a^{\dagger}_{-}$ act on $\psi$ and $a_{+},a^{\dagger}_{+}$ act on
$\chi$. Using these annihilation and creation operators, we obtain the mass-
squared operators for $\psi$ and $\chi$

$\displaystyle m^{2}_{\psi}=\alpha
a^{\dagger}_{-}a_{-},~{}~{}~{}m^{2}_{\chi}=\alpha(a^{\dagger}_{+}a_{+}+1).$
(35)

Note that (35) means an existence of chiral fermion in flux compactification
since zero-mode of $\psi$ is massless but zero-mode of $\chi$ is massive. We
can rewrite (31) in terms of Dirac fermion $\psi_{Lj}$ and $\Psi_{n,j}$,

$\displaystyle\psi_{Lj}=\left(\begin{array}[]{cc}\psi_{0,j}\\\
0\end{array}\right),~{}~{}~{}\Psi_{n,j}=\left(\begin{array}[]{cc}\psi_{n+1,j}\\\
\bar{\chi}_{n,j}\end{array}\right),$ (40)

and obtain

$\displaystyle\mathcal{L}_{fermion~{}mass}=\sqrt{\alpha(n+1)}\overline{\Psi}_{n,j}\Psi_{n,j}.$
(41)

We conclude that $\psi_{Lj}$ is massless and the KK mass of fermion
$\Psi_{n,j}$ is given by

$\displaystyle m^{2}_{fermion}=\alpha(n+1).$ (42)

####  2.2.3 SU(2) gauge field

Decomposing (18) and focusing on $F^{a}_{\mu 5}F^{a\mu 5}+F^{a}_{\mu 6}F^{a\mu
6}$ terms (see[10]), we obtain a mass term for an SU(2) gauge field,

$\displaystyle\mathcal{L}_{mass}$ $\displaystyle=-\frac{1}{2}\partial
A_{\mu}^{a}\bar{\partial}A^{a\mu}+g^{2}[A_{\mu},\braket{\phi}]^{a}[A^{\mu},\braket{\bar{\phi}}]^{a}-\frac{g}{\sqrt{2}}\Big{\\{}-\partial
A_{\mu}^{a}[A^{\mu},\braket{\bar{\phi}}]^{a}+\bar{\partial}A^{a}_{\mu}[A^{\mu},\braket{\phi}]^{a}\Big{\\}}$
$\displaystyle=-\frac{1}{2}A^{a}_{\mu}[-\mathcal{D}\bar{\mathcal{D}}]A^{a\mu}.$
(43)

Diagonalizing the covariant derivatives, we find the KK mass of the SU(2)
gauge field

$\displaystyle m^{2}_{YM}=\alpha\left(\begin{array}[]{ccc}n_{1}&0&0\\\
0&n_{2}&0\\\ 0&0&n_{3}+1\\\ \end{array}\right).$ (47)

Note that Abelian gauge field is not expressed such as (43) since the part of
commutator $[A_{M},A_{N}]^{a}$ is absent in the case of Abelian gauge theory.
Therefore, the KK spectrum of Abelian gauge field is nothing but an ordinary
KK mass spectrum $m^{2}_{U(1)}=(n/R)^{2}+(m/R)^{2}$ ($n,m$ is integer and
$R=L/(2\pi)$).

##  3 Analysis on the divergence structure of loop integral and classification
of interaction terms

In this section, we systematically analyze the divergence structure of the
quantum corrections to WL scalar mass and classify possible interactions
providing a finite mass.

###  3.1 The divergence structure of loop integral: part 1

Figure 1: Example of quantum corrections to the WL scalar mass at one-loop.
The fields in the loop are not specified in this figure.

In this subsection, we investigate the divergence structures for quantum
correction to the WL scalar mass at one-loop. In general, there are two types
of Feynman diagrams in figure 1. From these diagrams and the results of
section 2.2, the general form of loop integral in the quantum correction can
be written as

$\displaystyle I(x;a,b)$
$\displaystyle=\sum_{n=0}^{\infty}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{k^{2a}}{\left(k^{2}+\alpha\left(n+x\right)\right)^{b}}$
$\displaystyle=\frac{1}{\alpha^{b-a}}\left(\frac{4\pi}{\alpha}\right)^{\epsilon-2}\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon+b-a-2\right)}{\Gamma(b)\Gamma\left(2-\epsilon\right)}\zeta[\epsilon+b-a-2,x],$
(48)

where the dimensional regularization was employed for loop integral in the
second line. $\Gamma(z)$ is a gamma function, $\zeta[s,a]$ is Hurwitz zeta
function which is summarized in appendix A, and $d=4-2\epsilon$ dimensions.
$x$ is the part of KK mass characterized by the field running in the loop.
$x=1/2$ corresponds to the KK mass of scalar field (30). $x=1$ mainly
corresponds to the KK mass of fermion field (42). $x=0$ mainly corresponds to
the KK mass of SU(2) gauge field (47). $a$ denotes the number of derivatives
acting on the single field and $b$ corresponds to the number of the
propagator. From our interest, we focus on $b=1$ or $b=2$ since we consider
one-loop correction to WL scalar mass.

Since the WL scalar cannot have a bare mass term, the loop integral and mode
sum for one-loop correction to WL scalar mass must be finite to realize
nonvanishing finite WL scalar mass. To clarify this point, we investigate

$\displaystyle
J(x;a,b)\equiv\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon+b-a-2\right)}{\Gamma\left(2-\epsilon\right)}\zeta[\epsilon+b-a-2,x]$
(49)

in (48). For one-loop corrections, it is enough to consider the case $b=1$ or
$b=2$. In the case of $b=1$, the Gamma function part of (49) is expressed by

$\displaystyle\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon-a-1\right)}{\Gamma\left(2-\epsilon\right)}$
$\displaystyle=(a+2-\epsilon-1)(a+2-\epsilon-2)\cdots(2-\epsilon)\Gamma(\epsilon-a-1)$
$\displaystyle=(-1)^{a}\Gamma(\epsilon-1).$ (50)

Thus, $J(x;a,1)$ becomes

$\displaystyle J(x;a,1)=(-1)^{a}\Gamma(\epsilon-1)\zeta[\epsilon-a-1,x].$ (51)

In the case of $b=2$, the same part of (49) is expressed by

$\displaystyle\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon-a\right)}{\Gamma\left(2-\epsilon\right)}$
$\displaystyle=(a+2-\epsilon-1)(a+2-\epsilon-2)\cdots(2-\epsilon)\Gamma(\epsilon-a)$
$\displaystyle=(-1)^{a}(\epsilon-a-1)\Gamma(\epsilon-1).$ (52)

Thus, $J(x;a,2)$ becomes

$\displaystyle
J(x;a,2)=(-1)^{a}(\epsilon-a-1)\Gamma(\epsilon-1)\zeta[\epsilon-a,x].$ (53)

Here, Gamma function and Hurwitz zeta function can be expanded in $\epsilon$

$\displaystyle\Gamma(\epsilon-1)$
$\displaystyle=\frac{\Gamma(\epsilon)}{\epsilon-1}=-\left(\frac{1}{\epsilon}-\gamma_{E}+1+\mathcal{O}(\epsilon)\right),$
(54) $\displaystyle\zeta[\epsilon-p,x]$
$\displaystyle=\zeta[-p,x]+\zeta^{(1,0)}[-p,x]\epsilon+\mathcal{O}(\epsilon^{2}),$
(55)

where $\gamma_{E}=0.5772\cdots$ is the Euler-Mascheroni constant, $p$ is an
arbitrary positive integer, and
$\zeta^{(1,0)}[s,a]=\partial\zeta[s,a]/\partial s$. Using (54), (55) and (115)
in appendix A, the condition for $\Gamma(\epsilon-1)\zeta[\epsilon-p,x]$ being
finite is that $p$ is even,

$\displaystyle\Gamma(\epsilon-1)\zeta[\epsilon-p,x]=\mathrm{finite},~{}~{}~{}\mathrm{if}~{}p=even.$
(56)

Applying this result to (51) and (53), $J(x;a,1)$ takes finite value at odd
$a$, $J(x;a,2)$ does at even $a$.

###  3.2 Classification of interaction terms: part 1

We classify the interaction terms providing finite correction to WL scalar
mass at one-loop. We consider interaction terms which has no derivatives
acting on $\varphi$ or $\bar{\varphi}$ because we consider one-loop
corrections to WL scalar mass.

####  3.2.1 Four-point interaction

Four-point interaction term generates a correction to WL scalar mass of the
left one in figure 1. The diagram corresponds to $J(x;a,1)$ ($a$: odd), from
which we can guess the four-point interaction terms as follows,

  * •

scalar field loop

$\displaystyle
J(1/2;a,1)\rightarrow\bar{\varphi}\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a}}\overline{\Phi}\partial^{\mu_{1}}\cdots\partial^{\mu_{a}}\Phi,$
(57)

  * •

fermion field loop

$\displaystyle
J(1;a,1)\rightarrow\bar{\varphi}\varphi\bar{\psi}({\ooalign{\hfil/\hfil\crcr$\partial$}})^{2a-1}\psi,$
(60)

  * •

SU(2) gauge field loop

$\displaystyle
J(0;a,1)\rightarrow\bar{\varphi}\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a}}A^{a}_{\nu}\partial^{\mu_{1}}\cdots\partial^{\mu_{a}}A^{a\nu}.$
(61)

We did not consider a four-point interaction with such as
$\bar{\varphi}\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a}}\bar{\psi}\partial^{\mu_{1}}\cdots\partial^{\mu_{a}}\psi$
since the fermion mass $m_{fermion}=\sqrt{\alpha(n+1)}$ is emerged from a
numerator in the fermion propagator and then the form of Hurwitz zeta function
is complicated. On the other hand, $({\ooalign{\hfil/\hfil\crcr$k$}})^{2a-1}$
is obtained by (60). Computing quantum correction, the trace of  / $k$ from a
numerator in the propagator of fermion multiplied by
$({\ooalign{\hfil/\hfil\crcr$k$}})^{2a-1}$ is given by $k^{2a}$. If $a$ is
odd, this term contributes to quantum correction to WL scalar mass.

####  3.2.2 Three-point inteaction

Three-point interaction term generates a correction of the right one in figure
1. The diagram corresponds to $J(x;a,2)$ ($a$: even), from which we can guess
the three-point interaction terms as follows,

  * •

scalar field loop

$\displaystyle J(1/2;0,2)$
$\displaystyle\rightarrow\bar{\varphi}\overline{\Phi}\Phi+\varphi\overline{\Phi}\Phi,$
(62) $\displaystyle J(1/2;a,2)$
$\displaystyle\rightarrow\bar{\varphi}\partial_{\mu_{1}}\cdots\partial_{\mu_{a/2}}\overline{\Phi}\partial^{\mu_{1}}\cdots\partial^{\mu_{a/2}}\Phi+\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a/2}}\overline{\Phi}\partial^{\mu_{1}}\cdots\partial^{\mu_{a/2}}\Phi,$
(63)

  * •

fermion field loop

$\displaystyle J(1;a,2)$
$\displaystyle\rightarrow\bar{\varphi}\bar{\psi}({\ooalign{\hfil/\hfil\crcr$\partial$}})^{a-1}\psi+\varphi\bar{\psi}({\ooalign{\hfil/\hfil\crcr$\partial$}})^{a-1}\psi,$
(68)

  * •

SU(2) gauge field loop

$\displaystyle
J(0;a,2)\rightarrow\bar{\varphi}\partial_{\mu_{1}}\cdots\partial_{\mu_{a/2}}A^{a}_{\nu}\partial^{\mu_{1}}\cdots\partial^{\mu_{a/2}}A^{a\nu}+\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a/2}}A^{a}_{\nu}\partial^{\mu_{1}}\cdots\partial^{\mu_{a/2}}A^{a\nu},$
(69)

where $J(1/2;0,2)$ is allowed because of $\zeta[0,1/2]=0$. We are particularly
interested in the interaction term in (62), therefore we will discuss in
section 4.

###  3.3 The divergence structure of loop integral: part 2

In more general, we can consider the interaction term with coefficient
depending on KK mode. The more general form of loop integral in the quantum
correction is given by

$\displaystyle I^{\prime}(x;a,b)$
$\displaystyle=\sum_{n=0}^{\infty}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{k^{2a}f(n)}{\left(k^{2}+\alpha\left(n+x\right)\right)^{b}}$
$\displaystyle=\frac{1}{(4\pi)^{d/2}}\frac{\Gamma\left(a+\frac{d}{2}\right)\Gamma\left(b-a-\frac{d}{2}\right)}{\Gamma(b)\Gamma\left(\frac{d}{2}\right)}\sum_{n=0}^{\infty}\frac{f(n)}{(\alpha(n+x))^{b-a-\frac{d}{2}}},$
(70)

where $f(n)$ is a coefficient generated by an interaction term depending on KK
mode $n$. Since the more complicated the form of $f(n)$ is, the more difficult
we express as Hurwitz zeta function, the discussion on the finiteness of loop
integral becomes hard to proceed. As a candidate, $f(n)=((\alpha(n+q))^{c}$
($q$ and $c$ are real numbers) is considered. In this section, we assume
$f(n)=\alpha(n+q)$ for simplicity. Thus, $I^{\prime}(x;a,b)$ is expressed by

$\displaystyle I^{\prime}(x;a,b)$
$\displaystyle=\frac{1}{\alpha^{b-a-1}}\left(\frac{4\pi}{\alpha}\right)^{\epsilon-2}\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon+b-a-2\right)}{\Gamma(b)\Gamma\left(2-\epsilon\right)}$
$\displaystyle\quad\quad\quad\quad\times\Big{(}\zeta[\epsilon+b-a-3,x]+(q-x)\zeta[\epsilon+b-a-2,x]\Big{)}.$
(71)

If $q\neq x$, the divergence will inevitably appears from either
$\zeta[\epsilon+b-a-3,x]$ or $\zeta[\epsilon+b-a-2,x]$. To avoid the
divergence and see whether the quantum correction is finite, we need to choose
$q=x$ (equivalent to the choice $f(n)=$ KK mass), and then investigate

$\displaystyle
K(x;a,b)\equiv\frac{\Gamma\left(a+2-\epsilon\right)\Gamma\left(\epsilon+b-a-2\right)}{\Gamma\left(2-\epsilon\right)}\zeta[\epsilon+b-a-3,x],$
(72)

in (71). Substituting $b=1$ or $b=2$ in (72) and using (50) or (52), we obtain

$\displaystyle K(x;a,1)$
$\displaystyle=(-1)^{a}\Gamma(\epsilon-1)\zeta[\epsilon-a-2,x],$ (73)
$\displaystyle K(x;a,2)$
$\displaystyle=(-1)^{a}(\epsilon-a-1)\Gamma(\epsilon-1)\zeta[\epsilon-a-1,x].$
(74)

Applying the result (56) to (73) and (74), $K(x;a,1)$ takes finite value at
even $a$, $K(x;a,2)$ does at odd $a$.

###  3.4 Classification of interaction terms: part 2

We consider the case of four-point interaction terms ($a$: even) and guess
their form providing finite quantum corrections to WL scalar mass,

  * •

scalar field loop

$\displaystyle K(1/2;0,1)$
$\displaystyle\rightarrow\bar{\varphi}\varphi\overline{\Phi}\left(a^{\dagger}a+\frac{1}{2}\right)\Phi,$
(75) $\displaystyle K(1/2;a,1)$
$\displaystyle\rightarrow\bar{\varphi}\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a}}\overline{\Phi}\left(a^{\dagger}a+\frac{1}{2}\right)\partial^{\mu_{1}}\cdots\partial^{\mu_{a}}\Phi,$
(76)

  * •

fermion field loop

$\displaystyle
K(1;a,1)\rightarrow\bar{\varphi}\varphi\bar{\psi}({\ooalign{\hfil/\hfil\crcr$\partial$}})^{2a-1}(a^{\dagger}a+1)\psi$
(79)

  * •

SU(2) gauge field loop

$\displaystyle
K(0;a,1)\rightarrow\bar{\varphi}\varphi\partial_{\mu_{1}}\cdots\partial_{\mu_{a}}A^{a}_{\nu}(a^{\dagger}a)\partial^{\mu_{1}}\cdots\partial^{\mu_{a}}A^{a\nu}.$
(80)

The case of three-point interaction term is hard to guess because the three-
point interaction term cannot be expressed in terms of a mass-squared
operator. Thus, we do not consider the three-point interaction terms in this
section. If we compute loop integral by using $f(n)=(\alpha(n+x))^{c}$, the
meaning of $c$ is the number of mass-squared operators.

One might think that it would be possible that finite quantum corrections from
(62) and (75) are cancelled if both interactions (62) and (75) are present in
a theory. However, we need not worry such a cancellation since the structure
of $J(1/2;0,2)=\Gamma(\epsilon)\zeta[\epsilon,1/2]$ is different from
$K(1/2;0,1)=\Gamma(\epsilon-1)\zeta[\epsilon-2,1/2]$.

###  3.5 Comments on interactions between the field with different KK mode
indices

Due to the presence of annihilation and creation operators, there are
interactions between the field with different KK mode indices. For example,
they correspond to $\varphi\chi_{n,j}\psi_{n+1,j}$ (see [9]) and $\varphi
A^{a}_{\mu,n,j}A^{a\mu}_{n+1,j}$ (see [10]). The generalization of loop
integral with these interaction terms is very complicated because we must use
two propagators in different KK modes and cannot systematically identify the
coefficient depending on KK mode $n$.

##  4 Illustration of nonvanishing finite WL scalar mass

In section 3.2, we have classified the interaction terms generating finite
quantum correction at one-loop. Now we focus on (62) since it has no
derivatives and is the simplest interaction term of all interaction terms in
section 3.2. Therefore, we explicitly calculate finite quantum corrections to
WL scalar mass from (62) by diagrammatic calculation and effective potential
analysis. This is the first example of WL scalar mass in flux
compactification.

###  4.1 Set up

The Lagrangian we consider is given by (1), (2) and (62),

$\displaystyle\mathcal{L}=-\frac{1}{4}F_{MN}F^{MN}-D_{M}\overline{\Phi}D^{M}\Phi+\kappa(\bar{\phi}\overline{\Phi}\Phi+\phi\overline{\Phi}\Phi),$
(81)

where $\kappa$ is a dimensionless coupling constant. Using the expansion of
$\phi=\braket{\phi}+\varphi$, the Lagrangian is deformed as

$\displaystyle\mathcal{L}$
$\displaystyle\supset-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-D_{\mu}\overline{\Phi}D^{\mu}\Phi-m^{2}_{scalar}\bar{\Phi}\Phi$
$\displaystyle\quad-
ig\sqrt{2\alpha}\bar{\varphi}\bar{\Phi}a^{\dagger}\Phi+ig\sqrt{2\alpha}\varphi\bar{\Phi}a\Phi-2g^{2}\bar{\varphi}\varphi\overline{\Phi}\Phi$
$\displaystyle\quad+\kappa(\bar{\varphi}\overline{\Phi}\Phi+\varphi\overline{\Phi}\Phi)+\kappa(\braket{\bar{\phi}}\overline{\Phi}\Phi+\braket{\phi}\overline{\Phi}\Phi),$
(82)

where we note that the unnecessary terms are omitted. To derive a four-
dimensional effective Lagrangian by KK reduction, we need to expand $\Phi$ in
terms of mode functions $\xi_{n,j}$

$\displaystyle\Phi=\sum_{n,j}\Phi_{n,j}\xi_{n,j}.$ (83)

Integrating over $T^{2}$, the four-dimensional effective Lagrangian is
obtained by

$\displaystyle\mathcal{L}_{4D}$
$\displaystyle=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-\partial^{\mu}\bar{\varphi}\partial_{\mu}\varphi$
$\displaystyle\quad+\sum_{n,j}\left(-D_{\mu}\overline{\Phi}_{n,j}D^{\mu}\Phi_{n,j}-\alpha\left(n+\frac{1}{2}\right)\overline{\Phi}_{n,j}\Phi_{n,j}\right.$
$\displaystyle\quad-
ig\sqrt{2\alpha(n+1)}\overline{\varphi}\overline{\Phi}_{n+1,j}\Phi_{n,j}+ig\sqrt{2\alpha(n+1)}\varphi\overline{\Phi}_{n,j}\Phi_{n+1,j}-2g^{2}\bar{\varphi}\varphi\overline{\Phi}_{n,j}\Phi_{n,j}$
$\displaystyle\quad+\kappa\bar{\varphi}\overline{\Phi}_{n,j}\Phi_{n,j}+\kappa\varphi\overline{\Phi}_{n,j}\Phi_{n,j}+\kappa\braket{\phi}_{I}\overline{\Phi}_{n,j}\Phi_{n,j}+\kappa\braket{\bar{\phi}}_{I}\overline{\Phi}_{n,j}\Phi_{n,j}\Big{)},$
(84)

where $\braket{\phi}_{I}$ and $\braket{\bar{\phi}}_{I}$ are expressed by

$\displaystyle\braket{\phi}_{I}=\int_{T^{2}}dx^{2}\braket{\phi}\bar{\xi}_{n,j}\xi_{n^{\prime},j^{\prime}},~{}~{}~{}\braket{\bar{\phi}}_{I}=\int_{T^{2}}dx^{2}\braket{\bar{\phi}}\bar{\xi}_{n,j}\xi_{n^{\prime},j^{\prime}}.$
(85)

When $\braket{\phi}=f\bar{z}/\sqrt{2}$, $\braket{\phi}_{I}$ and
$\braket{\bar{\phi}}_{I}$ lead to zero because of odd function with respect to
integral variables $z$ or $\bar{z}$.

###  4.2 Diagrammatic computation

Before computing a finite quantum correction, we first review that the quantum
correction to WL scalar mass is cancelled at one-loop in the case of
$\kappa=0$ (see [7]). Computation of Feynman diagram is expressed as

$\displaystyle\mathcal{I}_{4pt}$
$\displaystyle=-i2g^{2}\sum_{n,j}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{k^{2}+\alpha\left(n+\frac{1}{2}\right)},$
(86) $\displaystyle\mathcal{I}_{3pt}$
$\displaystyle=+i2g^{2}\sum_{n,j}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{\alpha(n+1)}{\left(k^{2}+\alpha\left(n+\frac{1}{2}\right)\right)\left(k^{2}+\alpha\left(n+\frac{3}{2}\right)\right)}.$
(87)

Note that $\mathcal{I}_{3pt}$ is the example of section 3.5. The sum of
$\mathcal{I}_{4pt}$ and $\mathcal{I}_{3pt}$ is obtained by

$\displaystyle\mathcal{I}_{4pt}+\mathcal{I}_{3pt}=-2ig^{2}|N|\sum_{n}\int\frac{d^{4}k}{(2\pi)^{4}}\left(\frac{n+1}{k^{2}+\alpha\left(n+\frac{3}{2}\right)}-\frac{n}{k^{2}+\alpha\left(n+\frac{1}{2}\right)}\right).$
(88)

By the shift $n\rightarrow n+1$ in the second term, (88) becomes zero.

For $\kappa\neq 0$, we get a new quantum correction to WL scalar mass from the
right diagram in figure 1. Computing the diagram results in

$\displaystyle\mathcal{I}$
$\displaystyle=+i\kappa^{2}\sum_{n,j}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{\left(k^{2}+\alpha\left(n+\frac{1}{2}\right)\right)^{2}}=\frac{i\kappa^{2}|N|}{\alpha^{2}}\left(\frac{4\pi}{\alpha}\right)^{\epsilon-2}J(1/2;0,2).$
(89)

Applying $J(1/2;0,2)$ to (53), we obtain
$J(1/2;0,2)=\Gamma(\epsilon)\zeta[\epsilon,1/2]$ and then

$\displaystyle\mathcal{I}$
$\displaystyle=\frac{i\kappa^{2}|N|}{\alpha^{2}}\left(\frac{4\pi}{\alpha}\right)^{\epsilon-2}\Gamma(\epsilon)\zeta[\epsilon,1/2]$
$\displaystyle=-i\frac{\kappa^{2}|N|\ln
2}{32\pi^{2}}\left(\frac{4\pi}{\alpha}\right)^{\epsilon}+\mathcal{O}(\epsilon),$
(90)

where

$\displaystyle\Gamma(\epsilon)$
$\displaystyle=\frac{1}{\epsilon}-\gamma+\mathcal{O}(\epsilon),$ (91)
$\displaystyle\zeta[\epsilon,1/2]$ $\displaystyle=0-\epsilon\frac{\ln 2}{2}$
(92)

are used in the second line of (90). This correction is finite in
$\epsilon\rightarrow 0$ limit. Thus, the quantum correction to WL scalar mass
at one-loop is given by

$\displaystyle\delta m^{2}=i\mathcal{I}=\frac{|N|\ln
2}{32\pi^{2}}\frac{\kappa^{2}}{L^{2}}.$ (93)

Note that we introduced a factor of torus area $L^{2}$, which comes from the
normalization factors for KK mode function. Obviously, $\delta m^{2}=0$ is
reproduced for $\kappa=0$ in six-dimensional scalar QED (see [7]). One of the
interesting phenomenological applications is that the quantum correction
$\delta m^{2}$ to WL scalar mass can be interpreted as Higgs mass. This idea
is based on gauge-Higgs unification, namely a zero-mode of $\varphi$ is
regarded as Higgs field. Even if the compactification scale is Planck scale
$1/L\sim\mathcal{O}(M_{Planck})$, Higgs mass could be realized by the
interaction term (62) generated by some dynamics at $\mathcal{O}$(1) TeV
scale. This is analogous to the mass of pion as a pseudo NG boson for chiral
symmetry. The reason why the pion mass is not Planck scale is that chiral
symmetry is dynamically broken at extremely lower energy scale comparing to
the Planck scale, namely, QCD scale. Note that the similar discussion cannot
be applied to the ordinary gauge-Higgs unification since $\kappa$ is replaced
by an $SU(2)_{L}$ gauge coupling in this scenario.

###  4.3 Effective potential analysis

Next, we consider the quantum correction to WL scalar mass in terms of
effective potential. In our setup (84), we read the KK mass spectrum of $\Phi$
to be $\alpha(n+1/2)-\kappa\braket{\phi}_{I}-\kappa\braket{\bar{\phi}}_{I}$.
Thus, the four-dimensional effective potential is given by

$\displaystyle
V=\sum_{n=0}^{\infty}\int\frac{d^{4}k}{(2\pi)^{4}}\ln\left(k^{2}+\alpha\left(n+\frac{1}{2}\right)-\kappa\braket{\phi}_{I}-\kappa\braket{\bar{\phi}}_{I}\right),$
(94)

where we take into account a degree of freedom of complex scalar field $\Phi$.
To obtain quantum correction to WL scalar mass from four-dimensional effective
potential, we differentiate effective potential with respect to
$\braket{\phi}_{I}$ and $\braket{\bar{\phi}}_{I}$. Thus, $\delta m^{2}$ is
obtained as

$\displaystyle\delta m^{2}$
$\displaystyle=\left.\frac{\partial^{2}V}{\partial\braket{\phi}_{I}\partial\braket{\bar{\phi}}_{I}}\right|_{\braket{\phi}_{I}=0}$
$\displaystyle=-\kappa^{2}\sum_{n=0}^{\infty}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{\left(k^{2}+\alpha\left(n+\frac{1}{2}\right)\right)^{2}}=i\mathcal{I}.$
(95)

This result (95) agrees with (89) or (93).

###  4.4 Explicit breaking of translational invariance in extra space

Before confirming explicit breaking of translational invariance in extra
space, we review NG boson of translational invariance in extra space in the
case of $\kappa=0$. If $\kappa=0$, the quantum correction to WL scalar mass
becomes zero. This physical reason can be understood from the fact that the
zero-mode of $\varphi$ is a NG boson of translational invariance in extra
space. The breaking of translational invariance by the background
$\braket{\phi}$ can be replaced for a shift of $\varphi$,

$\displaystyle\delta_{T}\varphi=(\epsilon\partial+\bar{\epsilon}\bar{\partial})\varphi+\frac{\bar{\epsilon}}{\sqrt{2}}f,$
(96)

where $\epsilon\equiv(\epsilon_{5}+i\epsilon_{6})/2$ and
$\epsilon_{5},~{}\epsilon_{6}$ means constant parameters of translations in
extra spaces. Focusing on the zero-mode of $\varphi$ and noticing
$\partial\varphi=\bar{\partial}\varphi=0$, we obtain

$\displaystyle\delta_{T}\varphi=\frac{f}{\sqrt{2}}\bar{\epsilon}.$ (97)

This shows that the zero-mode of $\varphi$ is identified with a NG boson under
the translation in extra spaces.

If $\kappa\neq 0$, $\varphi\overline{\Phi}_{n,j}\Phi_{n,j}$ (or
$\bar{\varphi}\overline{\Phi}_{n,j}\Phi_{n,j}$) in (84) is expected to break
the translational invariance. To confirm it, we consider the following local
six-dimensional transformation [9]

$\displaystyle\varphi^{\prime}=\varphi-\frac{1}{\sqrt{2}}\partial\Lambda,~{}~{}~{}\Phi^{\prime}=e^{g\Lambda}\Phi,~{}~{}~{}\overline{\Phi}^{\prime}=e^{-g\Lambda}\overline{\Phi},$
(98)

where $\Lambda=f(\epsilon\bar{z}-\bar{\epsilon}z)$. Infinitesimal
transformation of $\epsilon,\bar{\epsilon}$ is expressed as

$\displaystyle\delta_{\Lambda}\varphi=-\frac{1}{\sqrt{2}}\partial\Lambda,~{}~{}~{}\delta_{\Lambda}\Phi=g\Lambda\Phi,~{}~{}~{}\delta_{\Lambda}\overline{\Phi}=-g\Lambda\overline{\Phi}.$
(99)

Transformations of $\varphi$ and $\Phi$ are the combination of translation
$\delta_{T}$ and infinitesimal transformation $\delta_{\Lambda}$,

$\displaystyle\delta\varphi$
$\displaystyle=(\delta_{T}+\delta_{\Lambda})\varphi=\sqrt{2}f\bar{\epsilon},$
(100) $\displaystyle\delta\Phi$
$\displaystyle=(\delta_{T}+\delta_{\Lambda})\Phi=-i\sqrt{\alpha}(\epsilon
a^{\dagger}+\bar{\epsilon}a)\Phi.$ (101)

Using (25) and (83), we obtain

$\displaystyle\delta\Phi$
$\displaystyle=-i\sqrt{\alpha}\sum_{n,j}\Phi_{n,j}(\epsilon
a^{\dagger}+\bar{\epsilon}a)\xi_{n,j}=\sum_{n,j}\delta\Phi_{n,j}\xi_{n,j},$
(102) $\displaystyle\delta\Phi_{n,j}$
$\displaystyle=-i\sqrt{\alpha}(\epsilon\sqrt{n+1}\Phi_{n+1,j}+\bar{\epsilon}\sqrt{n}\Phi_{n-1,j}).$
(103)

For $\delta\overline{\Phi}_{n,j}$, it is given by complex conjugate of (103),

$\displaystyle\delta\overline{\Phi}_{n,j}$
$\displaystyle=+i\sqrt{\alpha}(\bar{\epsilon}\sqrt{n+1}\overline{\Phi}_{n+1,j}+\epsilon\sqrt{n}\overline{\Phi}_{n-1,j}).$
(104)

Let us confirm the explicit breaking of translational invariance of the
interaction term $\varphi\overline{\Phi}_{n,j}\Phi_{n,j}$. First, a
transformation of $\overline{\Phi}_{n,j}\Phi_{n,j}$ is

$\displaystyle\delta\left(\sum_{n,j}\overline{\Phi}_{n,j}\Phi_{n,j}\right)$
$\displaystyle=i\sqrt{\alpha}\sum_{n,j}\Big{(}\bar{\epsilon}\sqrt{n+1}\overline{\Phi}_{n+1,j}\Phi_{n,j}+\epsilon\sqrt{n}\overline{\Phi}_{n-1,j}\Phi_{n,j}$
$\displaystyle\quad\quad\quad\quad\quad-\epsilon\sqrt{n+1}\overline{\Phi}_{n,j}\Phi_{n+1,j}-\bar{\epsilon}\sqrt{n}\overline{\Phi}_{n,j}\Phi_{n-1,j}\Big{)}$
$\displaystyle=0.$ (105)

Thus, the mass term of $\Phi_{n,j}$ is invariant. For
$\varphi\overline{\Phi}_{n,j}\Phi_{n,j}$, a transformation is

$\displaystyle\delta\left(\sum_{n,j}\varphi\overline{\Phi}_{n,j}\Phi_{n,j}\right)$
$\displaystyle=(\delta\varphi)\sum_{n,j}\overline{\Phi}_{n,j}\Phi_{n,j}+\varphi\delta\left(\sum_{n,j}\overline{\Phi}_{n,j}\Phi_{n,j}\right)$
$\displaystyle=\sqrt{2}f\bar{\epsilon}\sum_{n,j}\overline{\Phi}_{n,j}\Phi_{n,j}\neq
0.$ (106)

This result means the explicit breaking of translational invariance in extra
spaces. For $\kappa\neq 0$, the zero-mode of $\varphi$ is identified with a
pseudo NG boson of translational invariance in extra spaces.

One might claim that the interaction terms (62) are not gauge invariant since
$\varphi$ or $\overline{\varphi}$ transforms under the gauge symmetry as (98).
In ordet to overcome such a claim, $\varphi$ or $\overline{\varphi}$ should be
expressed by a gauge invariant non-local Wilson line operator and the
interaciton terms (62) should be regarded as one of the terms of expanding the
Wilson line operators in small $\varphi$ or $\overline{\varphi}$. Noting that
the Wilson line operators333In non-Abelian case, the path ordered operation
must be taken into account, $U_{5,6}=P\exp[ig\oint A_{5,6}dx^{5,6}]$.

$\displaystyle U_{5}=\exp\left[ig\oint A_{5}dx^{5}\right],\quad
U_{6}=\exp\left[ig\oint A_{6}dx^{6}\right]$ (107)

can be written in terms of $\varphi,\overline{\varphi}$ and $z,\bar{z}$ as

$\displaystyle U_{5}=\exp\left[\frac{g}{\sqrt{2}}\oint(\varphi dz+\varphi
d\bar{z}-\overline{\varphi}dz-\overline{\varphi}d\bar{z})\right],$ (108)
$\displaystyle U_{6}=\exp\left[\frac{g}{\sqrt{2}}\oint(\varphi dz-\varphi
d\bar{z}+\overline{\varphi}dz-\overline{\varphi}d\bar{z})\right],$ (109)

we find that the cubic terms introduced in this paper can be expressed by the
non-local Wilson line operators

$\displaystyle
i(U_{5}-U_{5}^{\dagger})\overline{\Phi}\Phi-i(U_{6}-U_{6}^{\dagger})\overline{\Phi}\Phi$
$\displaystyle\supset 2\sqrt{2}ig\oint\varphi
d\bar{z}\overline{\Phi}\Phi-2\sqrt{2}ig\oint\overline{\varphi}dz\overline{\Phi}\Phi$
$\displaystyle=2\sqrt{2}ig_{4}\varphi\overline{\Phi}\Phi-2\sqrt{2}ig_{4}\overline{\varphi}\overline{\Phi}\Phi$
(110)

where $g_{4}$ is a four dimensional gauge coupling constant. Note that the
$\overline{\Phi}\Phi$ term cannot be included in (110). If this term is
allowed, the WL scalar mass would be divergent.

We comment on how the finite WL scalar mass can be expressed in terms of the
Wilson line operators. If the WL scalar mass is generated in the broken phase,
where the VEV of the WL scalar field is non-zero, it is straightforward to
express the WL scalar mass by the Wilson line operators as in the gauge-Higgs
unification. As for the WL scalar field mass in the present paper, the mass is
generated in the unbroken phase and is independent of the VEV of the WL scalar
field. Therefore, we cannot express the WL scalar field mass by the Wilson
line operators explicitly.

Under the constant shift of $A_{5}\to A_{5}-f\epsilon_{6}/2,A_{6}\to
A_{6}+f\epsilon_{5}/2$, the operators

$\displaystyle U_{5}-U_{5}^{\dagger}=2i\sin\left[g\oint
A_{5}dx^{5}\right],\quad U_{6}-U_{6}^{\dagger}=2i\sin\left[g\oint
A_{6}dx^{6}\right]$ (111)

are not obviously invariant, which means that the interaction terms (110)
explicitly break the shift symmetry. Clarifying the origin of the interaction
terms (110) is not easy and beyond the scope of this paper. We expect that the
origin of the interaction terms would be connected to the quantum gravity
effects, nontrivial backgrounds such as a vortex, or some non-perturbative
dynamics. We leave this issue for our future work.

##  5 Conclusion and Discussion

We have studied the possibility to realize nonvanishing WL scalar mass in flux
compactification. Using KK mass of various fields in the bulk, we have
generalized loop integrals in the quantum correction to the WL scalar mass and
have systematically analyzed their structure of divergence. The conditions for
the loop integral and the mode sum to be finite were derived. Then, we have
classified four-point and three-point interaction terms providing finite
quantum corrections to WL scalar mass at one-loop.

Of these interaction terms, we focused on
$\bar{\varphi}\overline{\Phi}\Phi+\varphi\overline{\Phi}\Phi$ type interaction
since these interaction terms have no derivatives and are the simplest of all
interaction terms. Using these interactions in a six dimensional scalar QED,
we have explicitly illustrated the nonvanishing finite quantum correction to
WL scalar mass in two ways: diagrammatic computation or effective potential
analysis, which is the first example of finite WL scalar mass in flux
compactification. This result is easy to understand, which analogous to the
case of pion. WL scalar was originally a NG boson of translational symmetry in
extra space and therefore massless. Introducing an interaction term
$\varphi\overline{\Phi}\Phi$ of explicit breaking of translational symmetry
made WL scalar into a pseudo NG boson and WL scalar obtain a finite mass.

Phenomenologically, our result can be applied to Higgs physics. If we regard
(93) as Higgs mass, $\kappa/L\approx\mathcal{O}(1)$ TeV. If the interactions
explicitly breaking translational invariance are generated around TeV scale by
some dynamics, Higgs mass will be expected to be obtained even if the
compactification scale is an order of the Planck scale. Extension to a
realistic model of electroweak symmetry breaking is very interesting to
consider in future study.

## Acknowledgments

This work is supported in part by JSPS KAKENHI Grant Number JP17K05420 (N.M.).

##  Appendix A The property of Hurwitz zeta function

Hurwitz zeta function is given by

$\displaystyle\zeta[s,a]=\sum_{n=0}^{\infty}\frac{1}{(n+a)^{s}}.$ (112)

It is known that Hurwitz zeta function is related to Riemann zeta function by
the following identical equation

$\displaystyle\zeta[s,1]$ $\displaystyle=\zeta(s),$ (113)
$\displaystyle\zeta[s,1/2]$ $\displaystyle=(2^{s}-1)\zeta(s).$ (114)

Since Riemann zeta function satisfies $\zeta(-2n)=0$ ($n$ is a positive
integer), Hurwitz zeta function also satisfies

$\displaystyle\zeta[-2n,1]=0,~{}~{}~{}\zeta[-2n,1/2]=0.$ (115)

In particular, $\zeta[s,1/2]$ is satisfied by $\zeta[0,1/2]=0$.

## References

  * [1] R. Blumenhagen, B. Kors, D. Lust and S. Stieberger, “Four-dimensional String Compactifications with D-Branes, Orientifolds and Fluxes,” Phys. Rept. 445, 1 (2007) [hep-th/0610327].
  * [2] L. E. Ibanez and A. M. Uranga, “String theory and particle physics: An introduction to string theory,” Cambridge Univ. Press (2012).
  * [3] E. Witten, “Some Properties of O(32) Superstrings,” Phys. Lett. 149B, 351 (1984).
  * [4] D. Cremades, L. E. Ibanez and F. Marchesano, “Computing Yukawa couplings from magnetized extra dimensions,” JHEP 0405 079 (2004) [hep-th/0404229].
  * [5] H. Abe, K. -S. Choi, T. Kobayashi and H. Ohki, “Higher Order Couplings in Magnetized Brane Models,” JHEP 0906 080 (2009) [hep-th/0903.3800].
  * [6] Y. Matsumoto and Y. Sakamura, “Yukawa couplings in 6D gauge-Higgs unification on $T^{2}/Z_{N}$ with magnetic fluxes,” PTEP 2016 (2016) 5, 053B06 [hep-th/1602.01994].
  * [7] W. Buchmuller, M. Dierigl, E. Dudas and J. Schweizer, “Effective field theory for magnetic compactifications,” JHEP 1706 (2017) 039 [hep-th/1611.03798].
  * [8] D. M. Ghilencea and H. M. Lee, “Wilson lines and UV sensitivity in magnetic compactifications,” JHEP 1706 039 (2017) [hep-th/1703.10418].
  * [9] W. Buchmuller, M. Dierigl, E. Dudas, “Flux compactifications and naturalness,” JHEP 1808 151 (2018) [hep-th/1804.07497].
  * [10] T. Hirose and N. Maru, “Cancellation of One-loop Corrections to Scalar Masses in Yang-Mills Theory with Flux Compactification,” JHEP 1908 (2019) 054 [hep-th/1904.06028].
  * [11] M. Honda and T. Shibasaki, “Wilson-line Scalar as a Nambu-Goldstone Boson in Flux Compactifications and Higher-loop Corrections,” JHEP 03 (202) 031 [hep-th/1912.04581].
  * [12] T. Hirose and N. Maru, “Cancellation of One-loop Corrections to Scalar Masses in Flux Compactification with Higher Dimensional Operators,” [hep-th/2012.03494].
  * [13] J. Wess and J. Bagger, “Supersymmetry and supergravity,” Princeton University Press, Princeton U.S.A. (1992).
  * [14] Y. Hamada and T. Kobayashi, “Massive Modes in Magnetized Brane Models,” Prog. Theor. Phys. 128 (2012) 903 [hep-th/1207.6867].