text:
\section{Abstract}
 We demonstrate advantages of samples made by mechanical stacking of exfoliated
van der Waals materials for controlling the topological surface state of a 3-dimensional
topological insulator (TI) via interaction with an adjacent magnet layer.  We assemble bilayers with pristine interfaces using exfoliated flakes of the TI BiSbTeSe$_2$ and the magnet Cr$_2$Ge$_2$Te$_6$, thereby avoiding problems caused by interdiffusion that can affect interfaces made by top-down deposition methods.  The samples exhibit an anomalous Hall effect (AHE) with abrupt hysteretic switching.
For the first time in samples composed of a TI and a separate ferromagnetic layer, we
demonstrate that the amplitude of the AHE can be tuned via gate voltage with a strong
peak near the Dirac point. This is the signature expected for the AHE due to Berry curvature
associated with an exchange gap induced by interaction between the topological surface state
and an out-of-plane-oriented magnet.

\section{Background on Topological Insulator/Magnet Heterostructures}

Interactions between three-dimension topological insulators (TIs) and magnets can induce exotic topological phases like the quantum anomalous Hall~\cite{yu2010quantized, chang2013experimental,deng2020quantum,deng2021high}  or axion insulator states~\cite{mogi2017magnetic, liu2020robust}, and might be used to harness the properties of topological-insulator surface states in spintronic devices~\cite{tokura2019magnetic}.  Coupling of a topological-insulator surface state to a magnetization that is perpendicular to the surface is predicted to open an exchange gap at the Dirac point, thereby inducing Berry curvature in the electronic states near this gap~\cite{yu2010quantized, chang2013experimental}.  If the electron chemical potential is swept through the states in the vicinity of the exchange gap, this Berry curvature should cause an anomalous Hall effect (AHE) peaked when the chemical potential is in the gap.  For an ideal scenario at low temperatures, the peak value of the Hall conductivity for an exchange gap on just one surface of a TI thin film should reach a quantized value, $e^2/2h$~\cite{xiao2010berry,mogi2022experimental}.  Previous experiments have explored coupling of magnetism to TIs via two strategies: (1) by making the TI itself magnetic (by introducing magnetic dopants~\cite{chang2013experimental, hor2010development, chang2015high, he2018topological}
or by using an intrinsically magnetic TI material~\cite{ deng2020quantum, liu2020robust}) in which case the maximum temperature to which magnetism persists is generally less than about 40 K~\cite{mogi2017magnetic},
or (2) by proximity coupling of a non-magnetic TI deposited onto a separate magnetic layer, or vice-versa (using such magnets as EuS~\cite{Wei2013, Katmis2016},
Y$_3$Fe$_5$O$_{12}$~\cite{lang2014proximity, jiang2015independent,liu2020changes}, BaFe$_{12}$O$_{19}$~\cite{li2019magnetization},
Cr$_2$Ge$_2$Te$_6$~\cite{Alegria2014, Chong2018, yao2019record, Mogi2019, Mogi2021},
Tm$_3$Fe$_5$O$_{12}$~\cite{tang2017above, chen2019topological},
Zn$_{1-x}$Cr$_x$Te \cite{watanabe2019quantum}, Ga$_{1-x}$Mn$_{x}$As~\cite{Lee2018},
 Eu$_3$Fe$_5$O$_{12}$~\cite{zou2022enormous},  or MgAl$_{0.5}$Fe$_{1.5}$O$_4$ \cite{Riddiford2022})
in which case for many materials the magnetism survives to much higher temperatures.
However, the second strategy -- implementing well-controlled proximity-coupled TI/magnet samples -- has proved to be challenging, in that the resulting AHE signals are often very small and occasionally even non-hysteretic, so that in some cases it can be difficult to distinguish an AHE signal due to an exchange gap from nonlinear Hall responses expected for two-band conduction.
References~\cite{bhattacharyya2021recent, grutter2021magnetic}
present excellent recent reviews.
As of yet no proximity-coupled TI/magnet samples have been fabricated with sufficient control to exhibit the full prototypical response expected from the creation of an exchange-gap – observation of a peak in the AHE as a function of gating the electron chemical potential through the gap.  Materials issues that have made such measurements challenging include bulk conduction in the TI, spatial inhomogeneities due to dopants, difficulty in positioning the energy of the Dirac point to be well-separated from bulk bands, and  disorder at the interface.

Here, we study TI/magnet bilayers using an all-van der Waals (vdW) heterostructure, by stacking together in a glovebox exfoliated flakes of the TI BiSbTeSe$_{2}$ and the insulating magnet Cr$_2$Ge$_2$Te$_6$.  In contrast to previous work on TI/magnet samples grown by top-down deposition methods, or made by assembling exfoliated flakes with air exposure~\cite{Chong2018
, the use of a glovebox-assembled vdW structure ensures a defect- and diffusion-free atomically-flat interface~\cite{geim2013van}. Further, our use of a thin insulating vdW magnet enables both spatially-uniform magnetic coupling to the topological surface state and large tunability of the electron chemical potential using electrostatic gating.  We measure the anomalous Hall response while continuously controlling the contribution of the surface state by gating. The exfoliated samples exhibit an abrupt hysteretic switching signal for measuring the AHE, and for the first time in any TI/magnet bilayer structure we observe a clear peak in the amplitude of the AHE when the chemical potential is gated through the
Dirac point.
This gate dependences is the signature expected from Berry curvature induced by an exchange gap, as distinct from
other possible mechanisms such as a non-topological magnetized conducting layer at the interface~\cite{huang2012transport}
or spin Hall magnetoresistance~\cite{nakayama2013spin}.
Our interpretation is corroborated by theoretical calculations of the AHE produced by an exchange gap.
Our results demonstrate the advantages of a clean all-vdW exfoliated topological insulator-2D magnet system for manipulating topological surface states, and pave the way for improved control of topological magneto-electric effects in TI/magnet heterostructures.

\section{Materials Characterization and Sample Assembly}

In the most-commonly-studied 3D TIs (e.g., Bi$_2$Se$_3$, Bi$_2$Te$_3$), defects and unintentional doping can result in significant bulk conduction, making it difficult to isolate or tune the properties of the surface states. In contrast, for quarternary TIs like Bi$_{2-x}$Sb$_x$Te$_{3-y}$Se$_y$, it has been shown that by optimally tuning the Sb and Se compositions the surface Dirac cone can be engineered to lie isolated inside a well-formed bulk band gap and close to the Fermi energy \cite{Arakane2012,Xu2014, Han2018}.  For this reason, we utilize exfoliated films of quaternary BiSbTeSe$_2$ (BSTS) in our samples, similar to those shown to be of high quality in ref.~\cite{Han2018}.
In Fig.\ 1c, we show the measured longitudinal resistance $R_{xx}$  for an exfoliated 12 nm flake of BSTS (with no magnetic layer) as a function of bottom gate voltage at 2 K. The surface can be gated through the Dirac point as evidenced by the peak in the resistance and change in sign of the Hall coefficient ($R_H$) around $V_{bg} \sim -15$ V. The inset in Fig.\ 1c shows the temperature dependence of $R_{xx}$ in the same sample.  Below $\sim$ 80 K, $R_{xx}$ decreases and then saturates. This indicates that bulk conduction has largely frozen out in this regime and the transport is dominated by the 2D surfaces \cite{ ren2010large, ren2011optimizing, kushwaha2016sn, Xu2014}.

Using similar BSTS material, we fabricate dual-gated vdW heterostructure devices by stacking a few-layer flake of the vdW magnetic insulator Cr$_2$Ge$_2$Te$_6$ (CGT) on top of a BSTS flake (Fig.\ 1a). The samples are sandwiched between hBN layers (typically $<$ 10 nm) at the top and bottom to protect from degradation due to oxygen or moisture. The hBN layers also function as gate dielectrics.
Figure 1d shows a high-angle angular dark field scanning transmission electron microscopy image of a cross section of one of our hBN/BSTS/CGT/hBN structures demonstrating clean, abrupt interfaces,
nawithout the interdiffusion that is often unavoidable in deposited heterostructures (e.g., \cite{Riddiford2022}).
A top metal electrode and a Si back gate are used to tune the carrier densities of the top and bottom surfaces of the TI, respectively. Figure 1b shows an optical image of a representative device (channel length $L=6$ $\mu$m, width $W = 2$ $\mu$m) and the configuration used for transport measurements. Platinum electrical contacts are pre-patterned on the bottom hBN layer using electron beam lithography.

\section{Anomalous Hall Effect}

The out-of-plane magnetization of the CGT is expected to break time-reversal symmetry in the BSTS surface through proximity coupling and result in the opening of an exchange gap ($\Delta$)
in the adjacent Dirac surface state (Fig.\ 2b), as described by the following 2D Dirac Hamiltonian:
\begin{equation}
H(k) =   \hbar v_F(k_x\sigma_y - k_y\sigma_x) + \Delta \sigma_z
\end{equation}
where $\sigma_{x,y,z}$ are the Pauli spin matrices, $k_x$ and $k_y$ are in-plane wave vectors, and $v_F$ is the Fermi velocity.
States in the vicinity of the gap have non-zero Berry curvature, with equal and opposite values on opposite sides of the gap. Therefore when states on opposite sides of the gap have unequal occupations, the result is a nonzero Hall conductance, $\sigma_{xy}$. The peak value of the Hall conductance, when the electron chemical potential lies in the gap, should be $e^2/2h$ in the low temperature limit~\cite{xiao2010berry,mogi2022experimental}. Consequently, a peak of the Hall resistivity ($\rho_{xy} = \sigma_{xy}/(\sigma^2_{xx}+ \sigma^2_{xy}))$, where $\sigma_{xx}$ is the longitudinal conductivity) should be found when the chemical potential is tuned through the gap.

The panels in Fig.\ 2a and Fig.\ 2e show, respectively, the Hall voltage ($V_{xy}$) and longitudinal resistance ($R_{xx}$) measured as a function of an out-of-plane magnetic field for a dual-gated BSTS(10 nm)/CGT(2 nm) device at different fixed top-gate voltages ($V_{tg})$.  The measurements are performed at 4.4 K and under a constant bias current, with the bottom-gate voltage $V_{bg}$ fixed at 0 V. The current flows primarily through the BSTS layers since CGT is insulating at this temperature \cite{wang2018electric, lohmann2019probing, verzhbitskiy2020controlling}. The CGT flake (roughly 9 $\times$ 6 $\mu m$ in size) covers the entire device region, although because the BSTS flake (roughly 14 $\times$ 14 $\mu m$ in size) extends beyond the probed region there are additional current paths not covered by the magnetic layer that lead to shunting in $R_{xy}$. We observe hysteretic step-like changes in $V_{xy}$ corresponding to an anomalous Hall effect (AHE), indicative of a strong perpendicular anisotropy for the magnetism in CGT that is coupled to the BSTS surfaces.
In Fig.\ 2c, we show the anomalous Hall resistance ($R_{xy}$) after subtraction of the ordinary Hall background for the measurement done at $V_{tg} = 0$ V. We find a strong remanence at zero field and coercive fields $\sim 25$ mT.
The AHE magnetic-field sweep does not contain any non-monotonicity or ``hump'' sometimes associated with a topological Hall effect.
Interestingly, the form of this AHE response in the BSTS/CGT heterostructure is quite different from the magnetic properties of standalone few-layer CGT flakes, which typically show negligible remanence at zero field and close-to-zero coercive fields \cite{Gong2017, Idzuchi2019, wang2018electric}. This difference indicates that the proximity effect at the interface of the TI/magnet enhances the perpendicular magnetic anisotropy in our samples.

As the chemical potential of the top surface is tuned by varying the $V_{tg}$, we observe a modulation in the observed AHE signal. The amplitude of the AHE response is maximum at $V_{tg} = 0.55$ V and becomes smaller as the gate voltage is tuned on either side of this maximum. In Fig.\ 2d, we plot the extracted signal size of the anomalous Hall resistance response as a function of $V_{tg}$ (black trace). This trend tracks approximately with the gating behaviour of the longitudinal resistance $R_{xx}$
shown by the dotted blue trace. We therefore identify the maximum in the AHE response as due to tuning of the electron chemical potential within the exchange gap, as expected from the Berry curvature picture.
The association of these changes with surface states is further supported by the observation that the slope of the ordinary Hall background in $V_{xy}$ (Fig.\ 2a) changes sign across the AHE peak (e.g. compare $V_{tg} = -1.2$ V and $V_{tg} = 2.8$ V), indicating a change in surface charge carriers from electron-like to hole-like as the Fermi level is tuned across the gap.
Similar data from additional samples for the gating dependence of the AHE are presented in the supplemental information.

Although both the Hall resistance and the longitudinal resistance are sharply peaked, their peak values occur at slightly different values of gate voltage ($V_{tg} =$ 0.25 V for $R_{xx}$, compared to $V_{tg} =$ 0.55 V for $R_{xy}$). This difference can be explained by considering the contributions to transport from both the top and bottom surfaces of the TI. The Hall resistance should peak when the chemical potential is inside the exchange-induced band gap whereas the longitudinal resistance should peak at the chemical potential with minimal density of states, i.e., at the Dirac point of the bottom surface. The energy offset between the peaks suggest that the Dirac points of the two surfaces have a relative shift, as illustrated schematically in Fig.\ 3c where the bands are calculated based on the tight-binding Hamiltonian of Bernevig-Hughes-Zhang model with the parameters of BSTS (see details in the Supplemental Materials).
 The existence of such a relative shift is not surprising since the bottom surface is interfaced with hBN while the top surface is coupled to the CGT layer. With a finite band gap on top and a potential energy difference between both layers, we employ the Kubo-Streda formula to calculate the resistances as a function of chemical potential as shown in Fig.~\ref{figmuR}d. The predicted trend agrees well with experiments.  Details about the calculation method are included in the supplemental material.

It would be interesting to convert the measured peak in Hall resistance to Hall conductivity $\sigma_{xy}$ to see how the peak value compares the quantized value of $e^2/2h$.  However, the fact that the large TI flake in our device extends beyond the edge of the magnetic layer causes the local resistivity to be spatially inhomogeneous, preventing a determination of $\sigma_{xy}$ within just the portion of the TI under the magnetic layer.  In future work we plan to develop fabrication techniques that will allow more-precise definition the lateral device geometry, along with tests of whether our devices can be tuned into the quantum AHE regime.

\section{Longitudinal Resistance}

We now turn to the magnetic-field dependence of the longitudinal resistance. As shown in Fig.\ 2e, we observe a trend in $R_{xx}$ versus out-of-plane magnetic field, as is expected from magneto-transport in topological insulators~\cite{wang2012room}. Interestingly, for $V_{tg}$ close to the AHE peak (i.e., as the electron chemical potental approaches the exchange gap), a butterfly pattern emerges with small peaks of enhanced conduction near the magnetic-field values where the magnetization reverses. This type of conductance enhancement has been observed previously in samples exhibiting the anomalous Hall effect \cite{Checkelsky2012, Wei2013} and has been associated with transport along magnetic domain walls present in the multidomain state near the coercive field. Although our system is likely not in the quantum anomalous Hall regime, chiral domain wall states should also exist along magnetic domain walls in our samples, leading to an analogous enhancement.  A derivation of the domain wall dispersion is included in the supplemental material. This observation provides further evidence that our measurement is probing the physics of the proximity-coupled topological states at the vdW interface.

\section{Temperature Dependence}

The temperature dependence of $R_{xy}$ for the same BSTS/CGT device as Fig.\ 2  is shown in Fig.\ 4.  The AHE response becomes smaller as the temperature is increased above the base temperature of 4.4 K and is absent above $\sim 68$ K, close to the reported $T_c$ of bulk CGT \cite{Gong2017}. This again confirms that the observed magnetism in the top surface of BSTS is induced by proximity exchange coupling with CGT. We do observe a small  enhancement in  $T_c$ for our samples compared to the reported values of few-layer CGT  (40 - 47 K for 3 - 4 layers \cite{Gong2017} and 53-57 K for 5 layers~\cite{wang2018electric}). The extracted decrease in the AHE response and in the coercive fields ($B_c$) are plotted in Fig.\ 3b.

\section{Conclusions}
We have demonstrated that the use of mechanical assembly of van der Waals materials to form a pristine interface provides a strategy that avoids materials challenges which have inhibited research progress in studying topological insulator/magnet heterostructures grown by molecular beam epitaxy or other deposition techniques. Mechanically-assembled BiSbTeSe$_2$/Cr$_2$Ge$_2$Te$_6$ bilayers exhibit a strong anomalous Hall effect (AHE) with abrupt hysteretic switching.  The amplitude of the AHE is tunable with gate voltage, with a prominent peak when the electron chemical potential is swept through the Dirac point by gating, the first time this behavior has been demonstrated in a bilayer of a topological insulator proximity-coupled to a separate magnetic layer.  This is the signature expected from the Berry curvature of the topological surface states in the presence of an exchange gap due to interaction with the magnetic layer.  The high quality of mechanically-assembled van der Waals structures provides a platform for future studies of the quantum anomalous Hall and axion insulator states, and the rich phenomenology of topological magneto-electric phenomena predicted for these states.

\section{Acknowledgements}
We thank Maciej Olszewski, Bozo Vareskic, Liguo Ma, Avi Shragai, Patrick Hollister and Thow Min Cham for assistance with sample fabrication and measurements, and Nitin Samarth for comments on a draft of the paper.
Research at both Cornell and the Univ.\ of Washington was funded by the AFOSR/MURI project 2DMagic (FA9550-19-1-0390).  Research at Cornell was also supported by the NSF (DMR-2104268) and the Cornell Center for Materials Research (CCMR, supported by the NSF via grant DMR-1719875), and utilized the shared facilities of both the CCMR and the Cornell NanoScale Facility, a member of the National Nanotechnology Coordinated Infrastructure (supported by the
NSF via grant NNCI-2025233). Research at the University of Utah was supported by the National Science Foundation under the Quantum Leap Big Idea Grant No.\ 1936383 as well as the Kuwait Foundation for the Advancement of Sciences (KFAS) under project code "CB20-68EO-01".

\newpage

\section{References}
\bibliographystyle{naturemag}

\section{Supplementary information / Methods}
\section{Fabrication Details}
The dual-gated BSTS/CGT vdW heterostructure devices were fabricated using a layer-by-layer dry transfer technique described in \cite{wang2013one}, inside a glove box with O$_2$ and H$_2$O levels less than 0.5 ppm.

Pure elemental metals Bi, Sb, Te, Se were ordered from Sigma-Aldrich Co., with a purity of 5N grade. Quartz tubes with dimensions of an outer diameter of 14 mm, an inner diameter of 12 mm and a wall thickness of 2 mm were ordered from Technical Glass Products, Inc. The raw materials were prepared with a ratio of 1:1:1:2 of Bi:Sb:Te:Se, respectively. The materials were crushed and mixed for 30 minutes producing a mixed fine powder of weight ~3 g. The powder was then loaded into a carbon-coated quartz tube. The tube is prepared for sealing by flushing with argon gas four times followed by vacuum pumping of pressure ~10–6 torr. Finally, the tube was sealed producing an ampule of length 6–8 cm. The crystal was then grown via the method described in ref. \cite{Han2018}. Bulk crystals of Cr$_2$Ge$_2$Te$_6$ and hBN were commercially obtained (HQ graphene and 2D Semiconductors, respectively).

The individual flakes were exfoliated from the bulk crystals and picked up sequentially using a polymer stamp (Poly(Bisphenol A carbonate)). The stacked flakes were then transferred to a hBN flake which was pre-patterned with 6-nm thick Pt electrodes for making electrical contact to the BSTS and subsequently tip cleaned to remove resist residue. Finally, a metallic Ti/Pt gate was lithographically defined above the top hBN flake.
\section{Data from additional devices}
\renewcommand{\thefigure}{S1}

 \begin{figure}[ht!]
    \includegraphics[width =\linewidth]{dev1_data.eps}
    \caption{Additional data from a BSTS(14 nm)/CGT(4.2 nm) device (channel length L = 6 $\mu$m, width W = 2
 $\mu$m)
  at 4.5 K. (a) Anomalous Hall resistance and (b) longitudinal resistance measured at a function of magnetic field at $V_{tg}$ = 0. (c) Top-gate dependence of the size of the AHE signal (solid black line) and $R_{xx}$ measured at zero magnetic field (dashed blue line).}
\end{figure}
\renewcommand{\thefigure}{S2}
 \begin{figure}[ht!]
    \includegraphics[width =\linewidth]{dev2_data.eps}
    \caption{Additional data from a BSTS(10 nm)/CGT(5.0 nm) device (channel length L = 6 $\mu$m, width W = 2
 $\mu$m) at 4.5 K. (a) Anomalous Hall resistance and (b) longitudinal resistance measured at a function of magnetic field at $V_{tg}$ = 0. (c) Top-gate dependence of the size of the AHE signal (solid black line) and $R_{xx}$ measured at zero magnetic field (blue dashed line).}
\end{figure}

\section{Magneto-transport measurements}
Electrical transport measurements were performed in a vector magnet cryostat (1T, 1T, 7T) with a base temperature of 4.2 K for the devices studied in Fig.~2 and Fig.~3 in the main text. Characterization of BSTS flakes (Fig.~1c in the main text) was done in a Quantum Design Physical Properties Measurement System (PPMS).  A Keithley 6221 source meter was used to apply a constant ac excitation current (100 nA - 5 uA) at 97.7 Hz to the device. The longitudinal and Hall  signals were collected using lock-in techniques (Ametek DSP 7270). The top and bottom gate dc gate voltages were applied using two Keithley 2400 source meters.
\section{Details of theoretical calculations}
\subsection{Model}
For calculating the surface state dispersion, we employ the tight-binding model Hamiltonian
\begin{align}
\label{Ham3DTI}
    H=&M_{\bm{k}} \sigma_x + \lambda \sigma_z (s_y \sin k_x - s_x \sin k_y) + \lambda_z \sigma_y \sin k_z + M_{\bm{k}}' \sigma_0 s_0 \notag \\
    =&\begin{pmatrix}
       M_{\bm{k}}' &  -\lambda (i\sin k_x + \sin k_y) &   M_{\bm{k}}-i\lambda_z\sin k_z &  0 \\
-\lambda (-i\sin k_x + \sin k_y) &   M_{\bm{k}}' &     0    &  M_{\bm{k}}-i\lambda_z\sin k_z \\
   M_{\bm{k}}+i\lambda_z\sin k_z &         0 &   M_{\bm{k}}' &  \lambda (i\sin k_x + \sin k_y) \\
        0 & M_{\bm{k}}+i\lambda_z\sin k_z  & \lambda (-i\sin k_x + \sin k_y) & M_{\bm{k}}'  \\
\end{pmatrix}
\end{align}
where $\bm{\sigma}$ and $\bm{s}$ are Pauli matrices for orbital and spin separately.
The mass term $M_{\bm{k}}=\Delta + 2 \sum_i t_i (1-\cos k_i a_i)$ where $a_i$ is the lattice constant along $i$-th direction with $i=x,y,z$. $\lambda$ and $\lambda_z$ are spin-orbit coupling related hopping energies. $M_{\bm{k}}'=2 \sum_i t_{0i} (1-\cos k_i a_i)$ introduces the particle-hole asymmetry. We set $a_x=a_y=a$, $t_x=t_y=t$, $t_{0x}=t_{0y}=t_0$.
Following ref.\ \cite{https://doi.org/10.48550/arxiv.2203.14301}, we employ the parameters $t_0a^2=7.46 $ eV\AA$^2$, $t_{0z} a_z^2=-1.22$ eV\AA$^2$, $t a^2=24.7 $ eV\AA$^2$, and $t_{z} a_z^2=3.19$ eV\AA$^2$. The Fermi velocity terms are $2\lambda a = 3.23 $ eV\AA~ and $2\lambda_z a_z = 1.03$ eV\AA. The band gap is $\Delta=-0.0553 $eV. We take the lattice constant as $a=4.2 $\AA~and $a_z=1.5a$.

\subsection{Transport calculation}
We employ the Kubo-Streda formula with relaxation time approximation to calculate the conductivity \begin{align}
    \sigma_{\alpha \beta}=\frac{e^2}{h} \frac{\hbar^2}{Na^2}\sum_{\bm{k}} \int_{-\infty}^{+\infty} d\xi n_F(\xi) {\rm Tr}[\hat{v}^\alpha \frac{d G^R}{d \xi} \hat{v}^\beta (G^A-G^R) - \hat{v}^\alpha  (G^A-G^R) \hat{v}^\beta \frac{d G^A}{d \xi} ]
\end{align}
where $e$ is the elementary charge, $h$ is the Planck constant, $a$ is the lattice constant, $N$ is the number of unit cell, $n_F$ is the Fermi-Dirac distribution, $\hat{v}^\alpha=\frac{1}{\hbar}\frac{\partial H}{\partial k_\alpha}$ is the velocity operator, $G^{R/A}(\xi)=1/(\xi-H \pm i \gamma)$ are the retarded and advanced Green's functions. For simplicity, a constant self-energy $\gamma$ is employed.

The longitudinal conductivity can be obtained as
\begin{align}
    \sigma_{\alpha \alpha} =
    \frac{e^2}{h} \frac{\hbar^2}{Na^2}\sum_{\bm{k}} \int_{-\infty}^{+\infty} d\xi \frac{1}{2} \frac{d n_F(\xi)}{d \xi} {\rm Tr}[\hat{v}^\alpha (G^R-G^A) \hat{v}^\alpha (G^R-G^A)].
\end{align}

\subsection{Domain wall state}
Here we analyze the chiral edge states along a domain wall by using a low-energy continuum model, following \cite{semenoff2008domain,ren2016topological}. In this model, the Hamiltonian reads
\begin{eqnarray}
       H= v (k_x \sigma_y - k_y \sigma_x) + \sigma_z m(x,y)
       \label{ham1}
\end{eqnarray}
where $v$ is the Fermi velocity, $\bm{\sigma}$ are Pauli matrices of spin and $m(x,y)$ represents the position dependent exchange field from the ferromagnetic order.
We consider a single domain wall along $x=0$ where the mass term from the magnetic order is positive in one side while is negative in the other side, i.e., $m(x,y)=-\Delta$ for $x>0$ and $m(x,y)=\Delta$ for $x<0$ with $\Delta$ being positive. The mass term opens an energy gap $2\Delta$ in each region. We focus on the interface and look for domain wall states with eigenenergies within the energy gap.

To obtain these domain wall states, one can solve the Dirac equation as follows
\begin{eqnarray}
       \left[\begin{array}{cc}
              m(x,y)  &  v(k_x - ik_y)  \\
             v ( k_x + ik_y) & -m(x,y)
                \end{array}
          \right]  \left[\begin{array}{c}
              u(x,y)  \\
              v(x,y)
         \end{array}
          \right] =   E \left[\begin{array}{c}
              u(x,y)  \\
              v(x,y)
         \end{array}
          \right].
       \label{drceq1}
\end{eqnarray}
The eigenwavefunctions are expected to have the following form. Along $y$-direction, a plane wave is expected due to the transition symmetry. Along $x$ direction, localized states with an exponentially decaying wavefunction is expected since the energies are located within the energy gap. Therefore, the trial wavefunction for this equation is assumed to be, in region A,
\begin{eqnarray}
       u(x,y)=u_0 e^{-\lambda_0 x}e^{ik_y y}, \nonumber \\
	 v(x,y)=v_0 e^{-\lambda_0 x}e^{ik_y y},
       \label{wavfa}
\end{eqnarray}
in region B,
\begin{eqnarray}
       u(x,y)=u_1 e^{\lambda_1 x}e^{ik_y y}, \nonumber \\
	 v(x,y)=v_1 e^{\lambda_1 x}e^{ik_y y}
       \label{wavfcn}
\end{eqnarray}
where $\lambda_i$ ($i=0$,$1$) are positive.
By substituting this trial wavefunction into the Dirac equation, one can obtains
\begin{eqnarray}
	E^2=\Delta^2+ v^2(k_y^2-\lambda_i^2), \label{solt} \\
	\frac{v_0}{u_0}=\frac{\Delta + E}{\lambda-k_y}, \frac{v_1}{u_1}=\frac{\Delta - E}{\lambda+k_y}.
\label{solt1}
\end{eqnarray}
At $x=0$, both wavefunctions must have the same value which give rise to an additional contraint. Combining these equations, one can find that $\lambda=\Delta/v$ and the in-gap eigenstates exhibit the dispersion relation $E=-vk_y$.  Such states will therefore provide enhanced conductivity parallel to the domain wall.

\section{References}


text:
\section{Introduction}
\noindent
The stochastic-alpha-beta-rho (SABR) model~\citep{hagan2002sabr} has been widely adopted in the financial industry for the pricing and risk management of European options, owing to its intuitive and parsimonious parameterization. It has been a standard practice for practitioners to obtain the option price and delta from the asymptotic approximation of implied volatility~\citep{hagan2002sabr}, but the approximation loses accuracy and allows arbitrage as the variance of volatility becomes large. Despite numerous attempts to improve implied volatility approximation~\citep{obloj2007fine,paulot2015asym,lorig2017lsv,yang2017cev-sabr,choi2021sabrcev}, it does not seem possible to obtain an approximation accurate for all parameter ranges. To date, there are several \textit{full-scale} methods for pricing the SABR model: Monte--Carlo simulations~\citep{chen2012low,leitao2017one,leitao2017multi,cai2017sabr,choi2019nsvh,cui2021efficient}, finite difference methods~\citep{park2014sabr,vonsydow2019benchop}, and continuous-time Markov chains~\citep{cui2018ctmc-sabr}. These approaches require heavy computation and complex implementation compared with the approximation approach. Therefore, they are difficult to implement in practice.

We contribute to the literature by proposing a novel and efficient pricing scheme under the normal ($\beta=0$) SABR model. Despite being a special case, the normal SABR model has gained attention for modeling interest rates owing to its flexibility to allow negative asset prices. The normal model based on the arithmetic Brownian motion (BM) has long been used in fixed income markets, as opposed to the Black--Scholes model based on geometric BM.\footnote{See \citet{choi2022bachelier} for a survey of the normal (Bachelier) model.} The normal SABR model is a natural extension of the normal model with stochastic volatility that exhibits a volatility smile. \citet{antonov2015mixing} adopt the normal SABR model as a key component of the mixture approach for modeling a low-interest-rate environment. \citet{choi2019nsvh} propose the hyperbolic normal stochastic volatility (NSVh) model as a broader class of normal stochastic volatility models.

Based on the price transition law of \citet{choi2019nsvh}, we express the option price as a double sum of compound Gaussian quadratures consisting of Gauss--Hermite and Gauss--Laguerre quadratures. Our scheme is highly efficient in the sense that quadrature points with only 49 (i.e., $7\times 7$) nodes can produce very accurate option values. Unlike with the asymptotic implied volatility approach, the resulting prices obtained with our approach are arbitrage-free. Using a similar algorithm, we can also evaluate the option delta (equivalently, the probability distribution). Our study extends the previous literature on the normal SABR model. It compliments \citet{choi2019nsvh} by showing that the price transition law is useful for deterministic pricing as well as Monte--Carlo simulation. It extends the integral representations of the normal SABR model previously observed by \citet{korn2013exact} and \citet{antonov2015mixing}, as we provide an alternative representation and a practical numerical scheme together. The study also extends the stochastic volatility benchmark proposed by \citet{vonsydow2019benchop}. Given the availability of highly accurate option values, our method can serve as a testing benchmark for general SABR pricing methods proposed in the future.

The remainder of this paper is organized as follows. Section~\ref{s:model} reviews the normal SABR model. Section~\ref{s:quad} introduces the proposed quadrature scheme. Section~\ref{s:num} presents the numerical results. Finally, Section~\ref{s:conc} concludes the paper.

\section{Normal SABR model} \label{s:model}
\noindent
The SABR model~\citep{hagan2002sabr} is a stochastic volatility model specified by
\begin{equation}
\frac{dF_t}{F_t^\beta} = \sigma_t \, \rho dW_t \quad \text{and}\quad \frac{d\sigma_t}{\sigma_t} = \nu\, dZ_t,
\end{equation}
where $F_t$ and $\sigma_t$ are the processes for the forward price and volatility, respectively; $\nu$ is the volatility of volatility; $\beta\in[0,1]$ is the elasticity parameter; and
$W_t$ and $Z_t$ are the BMs correlated by $dW_t dZ_t = \rho\, dt$.
To simplify the notations for the remainder of this paper, we also denote
\begin{equation} \label{eq:rhoc}
\rho_\ast = \sqrt{1-\rho^2} \qtext{and} \xi = \frac12 \nu\sqrt{T},
\end{equation}
where $T$ is the option expiry.

The normal SABR model is the case with $\beta=0$. Unlike when $\beta>0$, the asset price $F_t$ can freely become negative in the normal SABR model, requiring no boundary condition at the origin.\footnote{In some studies~\citep{antonov2015mixing}, the model is explicitly referred to as the normal \textit{free-boundary} SABR. When $0<\beta< 1$, the absorbing boundary condition is imposed at the origin, resulting in a mass at zero. See \citet{gulisashvili2018mass,chen2019feeling,choi2021note} for details.}
Here, we are concerned with the undiscounted price and delta of the call option with strike price $K$ and expiry $T$:
$$ C(K) = E([F_T - K]^+) \qtext{and} D(K) = \frac{\partial C(K)}{\partial F_0}.$$
Since $F_0$ and $K$ affect the option value only through moneyness, $F_0-K$, under the normal SABR model, the option delta is also equal to the complementary probability distribution of $F_T$:
$$
D(K) = \frac{\partial C(K)}{\partial F_0} = -\frac{\partial C(K)}{\partial K} = \text{Prob}(F_T\ge K).
$$

\subsection{\citet{hagan2002sabr}'s normal volatility approximation} \label{ss:volapprox}
\noindent
\citet{hagan2002sabr} present the celebrated implied Black--Scholes volatility of the SABR model, from which the option price can be quickly computed using the Black--Scholes formula. In their study, \citet{hagan2002sabr} first derive the implied normal volatility, even for all $\beta\in[0,1]$, because they use perturbation from normal diffusion. They then convert the derived volatility to the equivalent Black--Scholes volatility using another approximation. However, for the normal SABR model, normal volatility is more relevant, as has been discussed. Therefore, the normal volatility approximation in the intermediate step~\citep[(A.59)]{hagan2002sabr} is preferred among practitioners, as it avoids further approximation error. The implied normal volatility for $\beta=0$ is given by
\begin{equation} \label{eq:hagan}
\sigma_\textsc{n}(K) = \sigma_0 \frac{z}{x(z)} \left(1+\frac{2-3\rho^2}{6} \xi^2\right)\\ \qtext{for} z = \frac{\nu}{\sigma_0}(K - F_0)
\end{equation}
where $K$ is the strike price, $\xi$ is given in Eq.~\eqref{eq:rhoc}, and $x(z)$ is defined by\footnote{At and near $z=0$, $z/x(z)$ should be evaluated as 1 from Taylor's expansion near $z=0$:
$$
\frac{x(z)}{z} = 1 - \frac{\rho}{2}\,z + \frac{3\rho^2-1}{6}\,z^2 - \frac{(5\rho^2 - 3)\rho}8\, z^3
+ \;\cdots.
$$}

\begin{equation} \label{eq:Hxv}
x(z) = \log\left( \frac{V(z)+z+\rho}{1+\rho} \right) \qtext{for} V(z) =\sqrt{1+2\rho z+z^2}.
\end{equation}
The option price is obtained by plugging $\sigma_\textsc{n}(K)$ into the Bachelier option formula:
\begin{equation} \label{eq:p_hagan}
C(K) \approx (F_0 - K) N(d_\textsc{n}) + \sigma_\textsc{n}\sqrt{T} \, n(d_\textsc{n})
\qtext{for} d_\textsc{n} = \frac{F_0 - K}{\sigma_\textsc{n}(K) \sqrt{T}},
\end{equation}
where $n(z)$ and $N(z)$ are the probability density function (PDF) and cumulative distribution function (CDF), respectively, of the standard normal distribution.

This volatility approximation, Eq.~\eqref{eq:hagan}, is the leading and first-order terms of the asymptotic expansion around a small $\xi$. Therefore, the accuracy of the approximation noticeably deteriorates when $\xi$ is not small. Despite various attempts to improve approximation~\citep{obloj2007fine,paulot2015asym,lorig2017lsv}, the intrinsic limitation arising from asymptotic approximation is difficult to overcome. Inaccurate option prices can also lead to inaccurate risk management. The option delta $D(K)$ implied from the volatility approximation,
\begin{equation} \label{eq:cdf_hagan}
D(K) = -\frac{\partial C(K)}{\partial K} \approx \frac{C(K-h) - C(K+h)}{2h} \quad\text{for small $h$},
\end{equation}
can deviate from the true delta significantly.

\subsection{Existing option price representation} \label{ss:int}
\noindent
Based on the heat kernel on hyperbolic geometry (i.e., Poincar\'e half-plane), \citet{henry2005general,henry2008analysis} derives a two-dimensional integral representation of the option price under the normal SABR model, and \citet{korn2013exact} correct mistakes in the formula. Their formula, with the simplification of \citet[\S~3.5.2]{antonov2019modern}, is as follows:
\begin{equation} \label{eq:korntang}
C(K) = [F_0-K]^+ + \frac{1}{2\pi\sqrt{2}}\frac{\sigma_0}{\rho_\ast\nu} \int_{s_0}^\infty Q(s) ds,
\end{equation}
where $\rho_\ast$ is given in Eq.~\eqref{eq:rhoc}, and
$$ s_0 = \cosh^{-1} \left( \frac{\sqrt{k^2 + \rho_\ast^2} - \rho k}{\rho_\ast^2}\right)
\qtext{with} k = z + \rho = \frac{\nu}{\sigma_0}(K-F_0) + \rho.
$$
The integrand $Q(s)$ is given by
$$
Q(s) = \left[ e^{-\frac{s}{2}} N\left(-\frac{s}{2\xi} + \xi\right) + e^{\frac{s}{2}}N\left(-\frac{s}{2\xi} - \xi\right) \right]\int_{\alpha_-(s)}^{\alpha_+(s)} \frac{d\alpha}{\sqrt{\cosh s - \cosh d(\alpha)}}
$$
where $d(\alpha)$ and $\alpha_\pm(s)$ are
\begin{gather*}
d(\alpha) = \cosh^{-1} \left( \frac12\left(\alpha + \frac1\alpha\right) + \frac{(k-\rho\alpha)^2}{2\rho_\ast^2 \alpha} \right),  \\
\alpha_\pm(s) = \rho k + \rho_\ast^2 \cosh s \pm \sqrt{\sinh^2 s - (k - \rho \cosh s)^2}.
\end{gather*}

\citet[\S~3.5.1]{antonov2015mixing} also derive a different two-dimensional integral representation:
\begin{equation} \label{eq:antonov}
C(K) = [F_0-K]^+ + \frac{\sigma_0}{\pi\nu} \int_{s_0}^\infty \frac{G(4\xi^2, s)}{\sinh s} \sqrt{\sinh^2 s - (k - \rho \cosh s)^2}\; ds
\end{equation}
where $G(t,s)$ is the CDF of the heat kernel~\citep[\S~3.4.5]{antonov2019modern}:
$$ G(t,s) = \frac{2e^{-t/8}}{t\sqrt{\pi t}} \int_s^\infty du\, u e^{-u^2/2t} \sqrt{\cosh u - \cosh s}.
$$
This representation is also based on the hyperbolic geometry heat kernel; however, the outcome differs from that of \citet{korn2013exact} because the integration is performed in a different order among the variables.

To evaluate these integral representations, Eqs.~\eqref{eq:korntang} and \eqref{eq:antonov}, one must resort to numerical integration. However, the use of generic numerical integrals over two-dimensional indefinite regions can be an extremely slow process because it requires evaluations of dense points. \citet{korn2013exact} present the option prices evaluated for two test cases, but they do not provide the implementation details of the numerical integration.

\section{Gaussian quadrature integration scheme} \label{s:quad}
\noindent
Using the results of \citet{choi2019nsvh}, we present a new integral representation of the option price, which differs from those in Section~\ref{ss:int}. This new representation is evidently easier to evaluate using Gaussian quadratures. First, we review \citet{choi2019nsvh}'s price transition law.

\subsection{NSVh model and the closed-form price transition formula} \label{ss:nsvh}
\noindent
\citet{choi2019nsvh} introduce the NSVh model as
\begin{equation}
dF_t = \sigma_t \left(\rho\,d\Zdrift[\lambda\nu/2]_t + \rho_\ast\,dX_t\right) \quad \text{and}\quad
\frac{d\sigma_t}{\sigma_t} = \nu\; d\Zdrift[\lambda\nu/2]_t,
\end{equation}
where $Z_t$ and $X_t$ are two independent BMs, and $\Zdrift[\mu]_t = Z_t + \mu\, t$ is the BM with a drift $\mu$. By introducing the drift term of $Z_t$, the NSVh model generalizes the normal SABR; the normal SABR model is the NSVh model with $\lambda=0$. Based on Bougerol's identity in hyperbolic geometry~\citep{alili1997full}, the NSVh model admits the following closed-form transition law for $\sigma_T$ and $F_T$~\citep[Corollary~1]{choi2019nsvh}:
\begin{equation} \label{eq:nsvh_price}
 \begin{aligned}
\sigma_T &= \sigma_0 \exp\left(\nu \bar{Z}_T\right),\\
F_T &\,{\buildrel d \over =}\, F_0 + \frac{\sigma_0\rho}{\nu}\,\left( e^{\nu \bar{Z}_T} -  e^{\lambda \nu^2 T/2} \right) + \frac{\sigma_0\rho_\ast}{\nu} \cos\theta\, \phi\hspace{-0.25em}\left(\nu \bar{Z}_T, \nu\sqrt{R_T^2 + \bar{Z}_T^2}\right),\\
\end{aligned} \end{equation}
where
$$\bar{Z}_T =  \Zdrift[(\lambda-1)\nu/2]_T \qtext{and} \phi(Z,D) = e^{Z/2} \sqrt{2\cosh D-2\cosh Z} \quad (Z\le D)$$
and $R_T$ is the two-dimensional squared Bessel process (i.e., $R_T=X_T^2+Y_T^2$ for two independent BMs, $X_T$ and $Y_T$), and $\theta \in [0, \pi]$ is a uniformly distributed random angle. As these random variables can be easily sampled, the transition law serves as a closed-form exact simulation scheme for the normal SABR model, which outperforms the method of \citet{cai2017sabr} for the $\beta=0$ case. See \citet{choi2019nsvh} for further details.

\subsection{New integral representation of the option price}
\noindent
The closed-form transition law in Eq.~\eqref{eq:nsvh_price} also enables an efficient pricing method for vanilla options. The random variables $Z_T$, $R_T$, and $\theta$ follow normal, exponential, and uniform distributions, respectively, whose PDFs are simple.

Let us introduce the integral variables $u$ and $v$, which represent $Z_T$ and $R_T$, respectively:
$$u = \frac{\Zdrift[-\nu/2]_T}{\sqrt{T}} = \frac{Z_T}{\sqrt{T}} - \xi \quad \left( u_\lambda = \frac{\bar{Z}_T}{\sqrt{T}} = u + \lambda\xi\right)\qtext{and} v = \frac{R_{T}^2}{T}.$$
Here, $u_\lambda$ is introduced to maintain the generality of the NSVh model. Under the normal SABR model (i.e., $\lambda=0$), however, $u_\lambda$ is equal to $u$.
The probability densities around variables $u$, $v$, and $\theta$ are respectively given by
$$f_u(u) = e^{-\xi u-\frac{\xi^2}2} n(u),\quad f_v(v) = \frac12 e^{-\frac12 v}, \qtext{and} f_\theta(\theta) = \frac{1}{\pi}.$$
Note that the term $e^{-\xi u-\frac{\xi^2}2}$ is the Radon-Nikodym derivative that arises from our definition of $u$. We define $u$ as such because the term makes the integrand more suitable for numerical integration, as opposed to naively defining $u = Z_T/\sqrt{T}$. For example, $e^{-\xi u}$ offsets $e^{Z/2}$ from $\phi(Z,D)$, mitigating the exponential growth as shown below.

The terminal asset price $F_T$ is expressed as a function of $u$, $v$, and $\theta$,
$$
F_T(u, v, \theta)= F_0 + \frac{\rho\sigma_0}{\nu}e^{2\lambda \xi^2} \left( e^{2\xi u} - 1 \right) + \frac{\rho_\ast\sigma_0}{\nu} \cosh \theta\; \phi\hspace{-0.25em}\left(2\xi u_\lambda, 2\xi \sqrt{v + u_\lambda^2}\right)
$$
Therefore, the undiscounted price of the European option struck at $K$ under the normal SABR model is expressed as an expectation over the three variables:
$$
C(K) = \int_{-\infty}^\infty du\, e^{-\xi u - \frac{\xi^2}2} n(u) \int_0^\infty \frac{dv}{2} e^{-\frac{v}2} \int_0^\pi \frac{d\theta}{\pi} [F_T(u,v,\theta) - K]^+.
$$
We further decompose the payout by introducing
\begin{align}
k(u) &= \frac{\nu}{\sigma_0} e^{-\xi u_\lambda} (K-F_0) - 2\rho e^{\lambda \xi^2}\; \sinh \left(\xi u\right), \label{eq:k} \\
h(u,v) &= \rho_\ast \left( 2\cosh\left(2\xi\sqrt{v + u_\lambda^2}\right) - 2\cosh(2\xi u_\lambda) \right)^{1/2}. \label{eq:h}
\end{align}
We can express the payout and option price as
\begin{gather}
F_T - K = \frac{\sigma_0}{\nu}e^{\xi u_\lambda}\left( h(u,v)\cos\theta - k(u) \right),\\
C(K) = \frac{\sigma_0}{\nu}e^{(\lambda-1/2)\xi^2} \int_{-\infty}^\infty du\, n(u) \int_0^\infty \frac{dv}{2} e^{-\frac{v}2} \int_0^\pi \frac{d\theta}{\pi} [h(u,v)\cos\theta - k(u)]^+.
\end{gather}

Next, we identify the integration region of positive payout. Note that, because $h(u,v)$ is a monotonically increasing function of $v$, taking values from $0$ to $\infty$ for all values of $u$, it is always possible to find $v^*(u)$ such that $h(u,v^*) = |k(u)|$:
\begin{equation} \label{eq:vstar}
v^*(u) = \frac1{4\xi^2}\,\text{acosh}^2 \left( \cosh\left(2\xiu_\lambda\right) + \frac{k^2(u)}{2\rho_\ast^2}  \right) - u_\lambda^2.
\end{equation}
Only when $k(u)=0$, $v^*(u)$ = 0.
As long as $v\ge v^*$, it is also possible to find $\theta^*$ such that $h(u,v)\cos\theta^* = |k(u)|$:
\begin{equation} \label{eq:thetastar}
\theta^*(u,v) = \text{arccos}\left(\frac{|k(u)|}{h(u,v)}\right)\quad \left(0\le \theta^* \le \frac{\pi}{2}\right).
\end{equation}
When $k(u)=0$ or $v\to\infty$, $\theta^*(u,v) = \pi/2$.

Using $v^*(u)$ and $\theta^*(u,v)$, we express the region of $(v,\theta)$, where the payout, $h(u,v)\cos\theta - k(u)$, is positive or negative, depending on the sign of $k(u)$.
When $k(u)\ge 0$, the payout is positive, if
$$ \{(v,\theta): 0\le v \le v^*(u) \qtext{and} 0\le \theta \le \theta^*(u,v) \}.
$$
Therefore,
\begin{align*}
&\int_0^\infty \frac{dv}{2} e^{-\frac{v}2} \int_0^\pi \frac{d\theta}{\pi} [h(u,v)\cos\theta - k(u)]^+ \\
&= \int_{v^*}^\infty \frac{dv}{2\pi}\, e^{-\frac{v}2} \int_0^{\theta^*} d\theta \, (h(u,v)\cos\theta - k(u))\\
&= \int_{v^*}^\infty \frac{dv}{2\pi}\, e^{-\frac{v}2} \left(\sqrt{h^2(u,v) - k^2(u)} - \theta^*(u,v)k(u) \right)
\end{align*}

When $k(u)<0$, it is easier to identify the region of negative payout:
$$ \{(v,\theta): 0\le v \le v^*(u) \qtext{and} \pi - \theta^*(u,v)\le \theta \le \pi \}.
$$
Using the identity $[x]^+ = x - [x]^-$, where $x=h(u,v)\cos\theta - k(u)$,
the integration along $\theta$ and $v$ becomes
\begin{align*}
\int_0^\infty & \frac{dv}{2} e^{-\frac{v}2} \int_0^\pi \frac{d\theta}{\pi} [h(u,v)\cos\theta - k(u)]^+ \\
&= -k(u) - \int_0^\infty \frac{dv}{2} e^{-\frac{v}2} \int_0^\pi \frac{d\theta}{\pi} [h(u,v)\cos\theta - k(u)]^- \\
&= -k(u) - \int_{v^*}^\infty \frac{dv}{2\pi}\, e^{-\frac{v}2} \int^{\pi}_{\pi-\theta^*} d\theta \, (h(u,v)\cos\theta - k(u))\\
&= -k(u) + \int_{v^*}^\infty \frac{dv}{2\pi}\, e^{-\frac{v}2} \left(\sqrt{h^2(u,v) - k^2(u)} + \theta^*(u,v)k(u) \right)
\end{align*}
Combining the two cases, we finally represent the option price as
\begin{equation} \label{eq:call}
\begin{aligned}
C(K)
=& \frac{\sigma_0}{\nu}e^{(\lambda-1/2)\xi^2} \int_{-\infty}^\infty du\, n(u)\left( [-k(u)]^+ + \int_{v^*}^\infty \frac{dv}{2\pi} e^{-\frac{v}2} \left(\sqrt{h^2(u,v) - k^2(u)} - \theta^*(u,v)|k(u)| \right)\right) \\
=& \frac{\sigma_0}{\nu}e^{(\lambda-1/2)\xi^2} \int_{-\infty}^\infty du\, n(u)\\
&\quad \left( [-k(u)]^+ + e^{-\frac{v^*}2} \int_{0}^\infty \frac{dv}{2\pi} e^{-\frac{v}2} \left(\sqrt{h^2(u,v^*+v) - k^2(u)} - \theta^*(u,v^*+v)|k(u)| \right)\right),
\end{aligned}
\end{equation}
where $k(u)$, $h(u,v)$, $v^*(u)$, and $\theta^*(u,v)$ are defined in
Eqs.~\eqref{eq:k}, \eqref{eq:h}, \eqref{eq:vstar}, and \eqref{eq:thetastar}, respectively.

The option delta, $D(K)=\text{Prob}(F_T\ge K)$, can be similarly obtained as
\begin{equation} \label{eq:delta}
\begin{aligned}
D(K) &= \int_{-\infty}^\infty du\, e^{-\xi u - \frac{\xi^2}2} n(u) \int_0^\infty \frac{dv}{2\pi} e^{-\frac{v}2} \int_0^\pi d\theta\, \boldsymbol{1}_{h(u,v)\cos\theta > k(u)} \\
&= \int_{-\infty}^\infty du\, e^{-\xi u - \frac{\xi^2}2} n(u) \left( \boldsymbol{1}_{-k(u)} + \mathrm{sgn}(k(u))\int_{v^*}^\infty \frac{dv}{2\pi} e^{-\frac{v}2}\,\theta^*(u,v)\right)\\
&= \int_{-\infty}^\infty du\, e^{-\xi u - \frac{\xi^2}2} n(u) \left(\boldsymbol{1}_{-k(u)} + \mathrm{sgn}(k(u))e^{-\frac{v^*}2} \int_0^\infty \frac{dv}{2\pi} e^{-\frac{v}2}\,\theta^*(u,v^* + v) \right),
\end{aligned}
\end{equation}
where $\mathrm{sgn}(x)$ is the sign function and $\boldsymbol{1}_x$ is the indicator function (i.e., $\boldsymbol{1}_x = 1$ if $x>0$ or 0 otherwise). From option delta $D(K)$, the CDF at $K$ can also be obtained as $ F(K) = 1 - D(K)$.

\subsection{Integration with Gaussian quadratures}
\noindent
The integral representations in the previous section allow for efficient evaluations because Gaussian quadratures exist for the probability densities of the variables. While we integrate $\theta$ analytically, we use the Gauss-Hermite quadrature for $u$ and the Gauss-Laguerre quadrature for $v$. Let $\{u_i\}$ and $\{w_i\}$ for $i=1,\ldots,N$ be the points and weights, respectively, of the Gauss-Hermite quadrature with respect to the weight function $n(u)$, and let $\{v_j\}$ and $\{\bar{w}_j\}$ for $j=1,\ldots,M$ be those of the Gauss--Laguerre quadrature associated with the weight function $e^{-\frac{v}2}$ ($v\ge 0$). The integrations with respect to $n(u)$ and $e^{-v/2}$ are then accurately approximated by the weighted sums,
\begin{equation*}
\int_{-\infty}^\infty du\; n(u) f(u) \approx \sum_{i=1}^N w_i f(u_i) \qtext{and}\int_{v^*}^\infty dv\; e^{-\frac{v}{2}} g(v) \approx e^{-\frac{v^*}{2}} \sum_{k=1}^M \bar{w}_j g(v^* + v_j),
\end{equation*}
for some functions $f(u)$ and $g(v)$, respectively. The points and weights of the Gaussian quadrature are optimally chosen from orthogonal polynomials, in the sense that the quadrature of size $N$ can precisely evaluate the moments of the probability density up to the order $2N-1$. The quadrature points and weights can be easily generated using public numerical libraries.\footnote{We use the python functions,  \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.eval_hermitenorm.html}{\texttt{scipy.special.eval\_hermitenorm}} and \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.roots_genlaguerre.html}{\texttt{scipy.special.roots\_genlaguerre}}, for Gauss--Hermite and Gauss--Laguerre quadratures, respectively.} Because we perform two-dimensional integration, we construct the compound quadrature $(u_i, v_j)$ with weight $w_i\bar{w}_j$, whose size is $NM$.

Let us define the indexed function values at the quadrature points by
$$
k_i = k(u_i), \quad v^*_i = v^*(u_i), \quad h_{ij} = h(u_i, v^*_i + v_j),\quad \theta^*_{ij}=\theta^*(u_i,v^*_i + v_j).
$$
Then, the option price in Eq.~\eqref{eq:call} is evaluated as a double sum over the composite Gaussian quadratures:
\begin{equation} \label{eq:quad_opt}
C(K) = \frac{\sigma_0}{\nu} e^{(\lambda-1/2)\xi^2} \sum_{i=1}^N  w_i \left([-k_i]^+ + e^{-\frac{v^*_i}{2}} \sum_{j=1}^M \frac{\bar{w}_j}{2\pi} \left(\sqrt{h_{ij}^2-k_i^2} - \theta^*_{ij}\,|k_i|\right) \right).
\end{equation}
The calculation can be performed simply as a vector--matrix operation:
$$
C(K)
= \frac{\sigma_0}{\nu} e^{(\lambda-1/2)\xi^2}
\boldsymbol{w}^T \left(\boldsymbol{b} + \boldsymbol{A} \boldsymbol{\bar{w}}\right),
$$
where $\boldsymbol{w}$ is the size $N$ vector of $w_i$, $\boldsymbol{\bar{w}}$ is the size $M$ vector of $\bar{w}_j$, $\boldsymbol{b}$ is the size $N$ vector of $[-k_i]^+$, and $\boldsymbol{A}$ is the $N\times M$ matrix of
$$ A_{ij} = \frac{1}{2\pi} e^{-\frac{v^*_i}{2}} \left[\sqrt{h_{ij}^2-k_i^2} - \theta^*_{ij}\,k_i\right].
$$
The option delta in Eq.~\eqref{eq:delta} can be similarly evaluated using the Gaussian quadrature:
$$
D(K) = \sum_{i=1}^N \left[w_i e^{-\xi u_i - \frac{\xi^2}2}  \left(\boldsymbol{1}_{-k_i} + \mathrm{sgn}(k_i) e^{-\frac{v^*_i}{2}}\sum_{j=1}^M \frac{\bar{w}_j}{2\pi} \theta^*_{ij}\right) \right].
$$

\section{Numerical Experiments} \label{s:num}
\begin{table}
\caption{Parameter sets used in the simulations} \label{t:param}
\begin{center}
	\begin{tabular}{|c||c|c|c|c|c|} \hline
	Case & $\sigma_0$ & $\nu$ & $\rho$ & $T$ & $F_0$ \\ \hline\hline
1 & 100 & 0.5 & 0 & 30 & 350 \\
2 & 100 & 0.5 & -0.3 & 30 & 350 \\
3 & 100 & 0.5 & -0.6 & 30 & 350 \\
	\hline
	\end{tabular}
\end{center}
\end{table}

\noindent
We test the accuracy of our Gaussian quadrature method for calculating the option price and delta. We use the three parameter sets shown in Table~\ref{t:param}. Case 1 is one of the parameter sets tested by \citet[Tables 3--4 and Figure 4]{korn2013exact}, but scaled by $10^4$. This case is understood as a parameter set for the swaption volatility smile in the unit of basis point. We price this case again because it is a challenging example. Given a large value of $\xi=1.37$, the error from Hagan's normal volatility approximation in Eq.~\eqref{eq:hagan} is significantly large, as shown by \citet{korn2013exact}. From Case 1, we vary correlation $\rho$ from 0\% to $-30\%$ in Case 2 and $-60\%$ in Case 3.

\begin{table}[ht]
\caption{\label{t:case1} The call option price and delta for Case 1 obtained from $N\times M$ Gaussian quadrature points and Hagan's volatility approximation. The exact values are obtained with $90\times 180$ points. The delta value and error are in \% unit.}
\begin{center} \small
\begin{tabular}{c||c|c|c|c|c||c|c|c|c|c} \hline
	Case 1 & \multicolumn{5}{c||}{Option price} & \multicolumn{5}{c}{Option delta (\%)} \\ \hline
	& \multicolumn{4}{c|}{Error from the exact value} & Exact & \multicolumn{4}{c|}{Error from the exact value} & Exact \\ \cline{2-5}\cline{7-10}
	Strike & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value \\ \hline
0 & ~0.16 & 0.14 & 0.07 & 92.28 & 572.02 & -0.19 & -0.10 & -0.06 & ~21.45 & 84.47 \\
100 & ~0.21 & 0.15 & 0.08 & 70.49 & 489.88 & -0.26 & -0.14 & -0.09 & ~21.85 & 79.42 \\
200 & ~0.30 & 0.17 & 0.10 & 49.61 & 414.24 & -0.41 & -0.23 & -0.14 & ~19.03 & 71.16 \\
300 & ~0.45 & 0.27 & 0.16 & 35.07 & 349.19 & -0.64 & -0.44 & -0.31 & ~8.43 & 58.06 \\
350 & -0.01 & 0.00 & 0.00 & 32.92 & 322.16 & -0.00 & -0.00 & ~0.00 & -0.00 & 50.00 \\
400 & ~0.45 & 0.27 & 0.16 & 35.07 & 299.19 & ~0.64 & ~0.44 & ~0.31 & -8.43 & 41.94 \\
500 & ~0.30 & 0.17 & 0.10 & 49.61 & 264.24 & ~0.41 & ~0.23 & ~0.14 & -19.03 & 28.84 \\
600 & ~0.21 & 0.15 & 0.08 & 70.49 & 239.88 & ~0.26 & ~0.14 & ~0.09 & -21.85 & 20.58 \\
700 & ~0.16 & 0.14 & 0.07 & 92.28 & 222.02 & ~0.19 & ~0.10 & ~0.06 & -21.45 & 15.53 \\
	\hline
\end{tabular} \vspace{1em}
\end{center}
\end{table}

\begin{table}[ht]
\caption{\label{t:case2} The call option price and delta for Case 2 obtained from $N\times M$ Gaussian quadrature points and Hagan's volatility approximation. The exact values are obtained with $90\times 180$ points. The delta value and error are in \% unit.}
\begin{center} \small
\begin{tabular}{c||c|c|c|c|c||c|c|c|c|c} \hline
	Case 2 & \multicolumn{5}{c||}{Option price} & \multicolumn{5}{c}{Option delta (\%)} \\ \hline
	& \multicolumn{4}{c|}{Error from the exact value} & Exact & \multicolumn{4}{c|}{Error from the exact value} & Exact \\ \cline{2-5}\cline{7-10}
	Strike & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value \\ \hline
0 & ~0.22 & 0.21 & 0.10 & 105.57 & 580.55 & -0.13 & -0.09 & -0.06 & ~22.48 & 86.50 \\
100 & ~0.28 & 0.19 & 0.11 & 82.14 & 495.84 & -0.20 & -0.09 & -0.07 & ~24.33 & 82.66 \\
200 & ~0.39 & 0.14 & 0.08 & 57.25 & 415.99 & -0.30 & -0.12 & -0.08 & ~25.09 & 76.51 \\
300 & ~0.48 & 0.27 & 0.14 & 33.37 & 344.19 & -0.45 & -0.24 & -0.15 & ~21.48 & 66.23 \\
350 & ~0.46 & 0.36 & 0.19 & 23.83 & 312.82 & -0.53 & -0.36 & -0.24 & ~16.24 & 59.01 \\
400 & ~0.52 & 0.29 & 0.18 & 17.55 & 285.36 & -0.13 & -0.10 & -0.08 & ~8.54 & 50.68 \\
500 & ~0.53 & 0.26 & 0.12 & 17.29 & 243.03 & ~0.56 & ~0.35 & ~0.22 & -7.05 & 34.52 \\
600 & ~0.34 & 0.16 & 0.07 & 28.63 & 214.53 & ~0.31 & ~0.17 & ~0.10 & -14.24 & 23.41 \\
700 & ~0.25 & 0.13 & 0.05 & 43.88 & 194.70 & ~0.20 & ~0.11 & ~0.06 & -15.70 & 16.83 \\
	\hline
\end{tabular} \vspace{1em}
\end{center}
\end{table}

\begin{table}[ht]
\caption{\label{t:case3} The call option price and delta for Case 3 obtained from $N\times M$ Gaussian quadrature points and Hagan's volatility approximation. The exact values are obtained with $90\times 180$ points. The delta value and error are in \% unit.}
\begin{center} \small
\begin{tabular}{c||c|c|c|c|c||c|c|c|c|c} \hline
	Case 3 & \multicolumn{5}{c||}{Option price} & \multicolumn{5}{c}{Option delta (\%)} \\ \hline
	& \multicolumn{4}{c|}{Error from the exact value} & Exact & \multicolumn{4}{c|}{Error from the exact value} & Exact \\ \cline{2-5}\cline{7-10}
	Strike & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value & 7$\times$7 & 10$\times$10 & 14$\times$14 & Hagan & value \\ \hline
0 & -0.04 & 0.04 & 0.03 & 72.67 & 569.45 & -0.03 & -0.09 & -0.07 & ~16.99 & 89.20 \\
100 & -0.04 & 0.13 & 0.01 & 54.56 & 481.52 & -0.13 & -0.04 & -0.04 & ~19.30 & 86.47 \\
200 & ~0.04 & 0.14 & 0.06 & 33.96 & 397.03 & -0.25 & -0.03 & -0.02 & ~21.92 & 82.16 \\
300 & ~0.13 & 0.13 & 0.06 & 10.93 & 318.23 & -0.33 & -0.09 & -0.06 & ~23.80 & 74.72 \\
350 & ~0.05 & 0.14 & 0.07 & -0.91 & 282.24 & -0.32 & -0.18 & -0.11 & ~23.26 & 68.93 \\
400 & ~0.30 & 0.18 & 0.09 & -11.93 & 249.61 & -0.39 & -0.30 & -0.19 & ~20.32 & 61.23 \\
500 & ~0.40 & 0.22 & 0.11 & -25.71 & 198.02 & ~0.67 & ~0.48 & ~0.33 & ~5.97 & 41.57 \\
600 & ~0.14 & 0.05 & 0.02 & -25.13 & 165.13 & ~0.33 & ~0.19 & ~0.11 & -5.36 & 25.56 \\
700 & ~0.07 & 0.02 & 0.01 & -17.86 & 144.45 & ~0.18 & ~0.10 & ~0.06 & -8.32 & 16.74 \\
	\hline
\end{tabular} \vspace{1em}
\end{center}
\end{table}

Tables~\ref{t:case1}--\ref{t:case3} present the results for the three cases. We first obtain the exact values of option price and delta from the Gaussian quadrature integration (Eqs.~\eqref{eq:call} and \eqref{eq:delta}, respectively) with a very dense quadrature set with $N=90$ and $M=180$. The option prices for $K=300$, 350, and 400 in Case 1 is consistent with the values reported by \citet[Table~3]{korn2013exact}. Then, we report the errors of the option price and delta from a small quadrature size and from Hagan's volatility approximation (Eqs.~\eqref{eq:p_hagan} and \eqref{eq:cdf_hagan}). We increase the quadrature size $N=M$ from  7, to 10, to 14, roughly doubling the total number of points from 49, to 100, to 196. The tables show that our quadrature integration accurately evaluates the price and delta. With only $7\times 7$ points, the price error is less than one basis point and the delta error is within 1\%. These errors are negligible for practical purposes. Figure~\ref{f:1} displays the implied normal volatility (left) and delta (right) for the three test cases. The result with $7\times 7$ points is visually indistinguishable from the exact value with $90\times 180$ points.
By contrast, Hagan's volatility approximation shows significant deviations from the exact value. Moreover, the option delta incorrectly goes below 0 or above 100\%. Given that the option delta is equal to the complementary CDF under the normal SABR model, this implies an arbitrage opportunity. It illustrates a possible danger of using a volatility approximation approach.

\begin{figure}[ht!]
	\caption{\label{f:1} The volatility smile (left) and delta (right) as functions of strike price for the three test cases in Table~\ref{t:param}.}
	\includegraphics[width=0.49\textwidth]{Case1-vol.pdf}
	\includegraphics[width=0.49\textwidth]{Case1-delta.pdf}\\ \vspace{1ex}
	\includegraphics[width=0.49\textwidth]{Case2-vol.pdf}
	\includegraphics[width=0.49\textwidth]{Case2-delta.pdf}\\ \vspace{1ex}
	\includegraphics[width=0.49\textwidth]{Case3-vol.pdf}
	\includegraphics[width=0.49\textwidth]{Case3-delta.pdf}\\ \vspace{1ex}
\end{figure}

\section{Conclusion} \label{s:conc}
\noindent
The normal ($\beta=0$) SABR model is an important special case of the popular SABR model for modeling interest rates as it allows negative asset prices. Although Hagan's implied volatility approximation is easy to use, its accuracy deteriorates and arbitrage is possible. This study provides an efficient numerical scheme for option prices and deltas. Based on the new price transition law, our proposed scheme adopts a compound Gaussian quadrature. Numerical tests show that a quadrature with only $7\times 7$ nodes accurately evaluates European options without arbitrage. Given that the SABR model still requires better pricing methods for general cases ($0\le \beta\le 1$), our new scheme will serve to provide the benchmarks for testing the methods proposed in the future.

\section*{Declarations of Interest}
\noindent The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

\begin{singlespace}

redpajama_set_name: RedPajamaArXiv


text:
\section{\label{sec:level1}Introduction\protect}
Quantum annealing (QA) is a metaheuristic, a generic approximate algorithm, to solve combinatorial optimization problems by quantum fluctuations \cite{Kadowaki1998,Kadowaki1998b,Brooke1999,Santoro2002,Santoro2006,Das2008,Tanaka_book2017}. Combinatorial optimization is a class of problems to find the global minimum of a cost function of many discrete variables. The time necessary to solve a combinatorial optimization problem grows typically exponentially as a function of the problem size, a hallmark of difficult problems. There exist various methods developed to solve each of those difficult problems, among which simulated annealing (SA) \cite{Kirkpatrick1983} is a classical metaheuristic that uses classical stochastic processes in analogy with thermal fluctuations in statistical mechanics.  QA has often been compared with SA to investigate if quantum effects may ever be useful to accelerate computation \cite{Ronnow2014}, a recent example of which is a careful study showing evidence of a limited quantum speedup \cite{AlbashLidar2017}.

The cost function of a combinatorial optimization problem can be expressed  as the Hamiltonian of an Ising model with complicated interactions \cite{Lucas2014}. This means that the search for an optimal solution is analogous to the identification of the ground state of an Ising spin glass, one of the most difficult problems in statistical physics \cite{Nishimori_book_sg}. Thus, the study of the performance of QA has important significance from the viewpoint not just of computer science but also of physics.

In the initial stage of QA, one applies a strong transverse field to the system to keep the system state in the trivial ground state with all spins (qubits) pointing to the same direction along the transverse field. One then reduces the transverse field toward zero, and the system evolves following the time-dependent  Schr\"odinger equation.  Since the initial trivial state and final non-trivial state are completely different, a quantum phase transition takes place at a certain strength of the transverse field during the course of QA.
Let us consider a process in which the strength of the transverse field, relative to the magnitude of the coefficients in the Ising model, is reduced from a very large value to zero in a finite amount of time $\tau$ \cite{Farhi2000,Farhi2002}.  The system stays close to the instantaneous ground state if $\tau$ is very large compared to the inverse square of the minimum energy gap between the ground state and the first excited state, $\tau\gg\Delta^{-2}_{min}$, according to the adiabatic theorem of quantum mechanics (see e.g. \cite{AlbashLidar2016,Morita2008}). If there exists a second-order quantum phase transition at a certain value of the transverse field, the energy gap decreases polynomially as a function of the system size according to finite-size scaling \cite{Nishimori_book}. Then, the problem can be solved in polynomial time and is considered easy. In contrast, if the quantum phase transition is of first order, the minimum gap decays exponentially as the system size increases, and the computation time increases exponentially. This implies that a first-order transition sets a serious bottleneck for QA. In this way, the study of quantum phase transition in the transverse field Ising model gives important information on the performance of QA.

The infinite-range Ising model with $p$-body ferromagnetic interactions under a transverse field ($p$-spin model) is a prototype problem that has been studied in the context of QA due to its amenability to analytical investigations \cite{Jorg2010,Bapst2012}.  It is known that this model has a first-order quantum phase transition for $p\ge 3$, and is therefore hard to solve by the conventional QA for the reason described above. Nevertheless, it was shown that this first-order transition can be avoided by the introduction of antiferromagnetic transverse interactions into the Hamiltonian, by which first-order phase transitions are reduced to second order \cite{Seki2012}. A similar phenomenon has been shown to exist in the Hopfield model \cite{Seki2015}, which involves randomness in interactions and is thus much more nontrivial compared with the simple $p$-spin model. Since the antiferromagnetic transverse interactions make the Hamiltonian non-stoquastic \cite{Bravyi2008}, the above example of the reduction of the order of phase transition may be taken as a remarkable example in which a non-stoquastic Hamiltonian shows an exponential speedup over its stoquastic counterpart \cite{NishimoriTakada2017}.  See also Refs. \citen{Farhi2012,Crosson2014,Hormozi2017} for numerical evidence that a non-stoquastic Hamiltonian may sometimes accelerate the performance of QA.

The goal of the present paper is to analyze the properties of first-oder quantum phase transitions in the $p$-spin model with antiferromagnetic transverse interactions.  In particular, we use the semi-classical method \cite{Braun1993,Garg1998,Garg1999,Bapst2012,Kechedzhi2016} to derive the explicit formula of the energy gap as a function of system parameters. Similar semi-classical analyses of the energy gap have already been reported in the context of QA for stoquastic Hamiltonians \cite{Bapst2012,Isakov2016,Jiang2017,Kechedzhi2016}. Our contribution is to generalize these results to the non-stoquastic case, by which we show how the energy gap behaves as the effects of the antiferromagnetic transverse interactions become prominent and quantum fluctuations play important roles \cite{Susa2017}.

We formulate and solve the problem in the next section. The final section concludes the paper.

\section{Semi-classical analysis of the phase diagram and the energy gap}

We first define the problem and then proceed to describe a semi-classical  approach to analyze the phase diagram and the energy gap.

\subsection{Problem formulation}

The problem we study in this paper is described by the Hamiltonian
\begin{align}
\label{eq:hamiltonian1}
\hat{H}(s,\lambda)=s\Big(\lambda\hat{H}_0+(1-\lambda)\hat{V}_{{\rm TI}}\Big)+(1-s)\hat{V}_{{\rm TF}}\quad (0\le s, \lambda \le 1),
\end{align}
where $H_0$ is the $p$-spin model, $\hat{V}_{{\rm TI}}$ is the antiferromagnetic transverse interaction, and $\hat{V}_{{\rm TF}}$ denotes the usual transverse field,
\begin{align}
\hat{H}_0=-N\left(\frac{1}{N}\sum_{i=1}^N\hat{\sigma}_i^z\right)^p,\quad
\hat{V}_{{\rm TI}}=N\left(\frac{1}{N}\sum_{i=1}^N\hat{\sigma}_i^x\right)^2,\quad
\hat{V}_{{\rm TF}}=-\sum_{i=1}^N\hat{\sigma}_i^x.
\end{align}
The initial values of the parameters are $s=0$ and $\lambda$ arbitrary since $\lambda$ disappears when $s=0$. We then let the system evolve adiabatically toward the goal of $s=\lambda =1$, where only the problem Hamiltonian $\hat{H}_0$ remains.  Notice that the Hamiltonian is stoquastic only when $\lambda=1$.  For any $\lambda$ less than 1, i.e. with the antiferromagnetic transverse interactions $\hat{V}_{{\rm TF}}$, the Hamiltonian is non-stoquastic, and the effects of the term $\hat{V}_{{\rm TI}}$ are stronger for smaller $\lambda$.

When $p>2$ and $\lambda=1$, this model is known to go under a first-order transition as a function of $s$ \cite{Jorg2010}. This first-order transition has been shown to be avoided by the introduction of $\hat{V}_{{\rm TI}}$ \cite{Seki2012}. We choose $p$ to be odd to avoid the trivial double degeneracy for even $p$.

It is convenient to rewrite the Hamiltonian as
\begin{align}
\hat{H}=-s\lambda N\left(\frac{2}{N}\hat{S}_z\right)^p+s(1-\lambda)N\left(\frac{2}{N}\hat{S}_x\right)^2-2(1-s)\hat{S}_x.
\end{align}
where
\begin{align}
\hat{S}_{x}=\frac{1}{2}\sum_{i=1}^N\hat{\sigma}_i^x,\;\;\hat{S}_{y}=\frac{1}{2}\sum_{i=1}^N\hat{\sigma}_i^y,\;\;\hat{S}_{z}=\frac{1}{2}\sum_{i=1}^N\hat{\sigma}_i^z.
\end{align}
The total spin operator
\begin{eqnarray}
\hat{S}^2=\hat{S}_x^2+\hat{S}_y^2+\hat{S}_z^2
\end{eqnarray}
commutes with the Hamiltonian, and therefore the total spin $S$ is conserved in the process of time evolution. In the initial ground state, the eigenvalue of $\hat{S}_x$ is maximum, $S=N/2$. The final ground state belongs to the same subspace. Therefore, we only need to consider the subspace with $S=N/2$. Moreover, it is known that the first excited state state also belongs to the same subspace \cite{Bapst2012}.

For later convenience, we normalize spin operators as $\hat{m}_i=\hat{S}_i/S$. The Hamiltonian per spin is
\begin{eqnarray}
\hat{E}=\frac{\hat{H}}{N}=-s\lambda\hat{m}_z^p+s(1-\lambda)\hat{m}_x^2-(1-s)\hat{m}_x.
\label{operator_E}
\end{eqnarray}
Let us write the eigenstate of the normalized spin operator as
\begin{eqnarray}
\hat{m}_z\ket{S,m}=m\ket{S,m},\;\;\left(m=-1,-1+\frac{1}{S},...,1\right).
\end{eqnarray}
We express the eigenstate of $\hat{E}$ as
\begin{eqnarray}
\ket{\Psi}=\sum_mC_m\ket{S,m}.
\end{eqnarray}
The Schr\"odinger equation can be written in terms of the coefficients as
\begin{align}
\nonumber
&EC_m=\Big(-s\lambda m^p+\frac{s(1-\lambda)}{2}(1-m^2+S^{-1})\Big)C_m\\
\nonumber
&-\frac{1-s}{2}\left(\sqrt{(1+m)(1-m+S^{-1})}\,C_{m-1/S}+\sqrt{(1-m)(1+m+S^{-1})}\,C_{m+1/S}\right) \\
\nonumber
&+\frac{s(1-\lambda)}{4}\Big(\sqrt{\left(1-m+2S^{-1}\right)\left(1+m-S^{-1}
\right)\left(1-m+S^{-1}\right)(1+m)}\,C_{m-2/S} \\
&+\sqrt{\left(1+m+2S^{-1}\right)\left(1-m-S^{-1}\right)\left(1+m+S^{-1}
\right)(1-m)}\,C_{m+2/S}\Big).
\end{align}
Since we are interested in the semi-classical limit ($S\gg 1$), we extract the leading terms as
\begin{align}
\label{eq:syure1}
\nonumber
EC_m&=\Big(-s\lambda m^p+\frac{s(1-\lambda)}{2}(1-m^2)\Big)C_m\\
\nonumber
&-\frac{1-s}{2}\sqrt{1-m^2}\left(C_{m-1/S}+C_{m+1/S}\right) \\
&+\frac{s(1-\lambda)}{4}(1-m^2)\left(C_{m-2/S} +C_{m+2/S}\right).
\end{align}

\subsection{Classical potential and the phase diagram}

The asymptotic behavior of the present model for $S\gg 1$ can be analyzed by the semi-classical method \cite{Braun1993,Garg1998,Garg1999,Bapst2012,Kechedzhi2016,Isakov2016,Jiang2017}, i.e. an asymptotic expansion in terms of $1/S$.
Let us write the coefficient $C_m$ as
\begin{eqnarray}
C_m \propto \exp \big(i\Phi(m)\big),
\label{Cm}
\end{eqnarray}
and expand the exponent as
\begin{eqnarray}
\Phi(m)=S\Phi_0(m)+\Phi_1(m)+\mathcal{O}(S^{-1}).
\label{Phim}
\end{eqnarray}
Equation (\ref{eq:syure1}) is then solved to leading order as
\begin{eqnarray}
\label{eq:ekinji}
E=w(m)+t_1(m)\cos\dot{\Phi}_0+t_2(m)\cos(2\dot{\Phi}_0),
\end{eqnarray}
where
\begin{align}
w(m)&=-s\lambda m^p+\frac{s(1-\lambda)}{2}(1-m^2),\\
\label{eq:t1t2}
t_1(m)&=-(1-s)\sqrt{1-m^2},\;\;
t_2(m)=\frac{s(1-\lambda)}{2}(1-m^2).
\end{align}

The starting point of the semi-classical analysis is to identify the effective classical potential, from which the quantum tunneling rate and the energy gap are calculated.  To this end, we first notice that $\dot{\Phi}_0(m)$ in Eq. (\ref{eq:ekinji}) can be identified with the momentum. The operator $\hat{q}$ defined by
\begin{equation}
\hat{q}=-i\frac{d}{dm}
\end{equation}
shifts the argument of $C(m)$ as
\begin{equation}
C(m+ j/S)=e^{ij\hat{q}/S}\, C(m)~(j=\pm 1, \pm 2,\cdots).
\end{equation}
Then, Eq. (\ref{eq:syure1}) can  be written as
\begin{align}
EC(m)=w(m)C(m)+t_1(m) \cos \left(\frac{\hat{q}}{S}\right)\, C(m)
+t_2(m) \cos \left(\frac{2\hat{q}}{S}\right) \, C(m).
\label{ECm0}
\end{align}
Comparison of this equation with Eq. (\ref{eq:ekinji}) reveals that $S\dot{\Phi}_0$ can be identified with the momentum in the classical limit. This means that the classical motion is allowed in the range where $\dot{\Phi}_0$ is real and thus $\cos\dot{\Phi}_0$ lies between $-1$ and 1.

Let us write Eq. (\ref{eq:ekinji}) as
\begin{align}
E=w(m)+2t_2(m)\left(\cos\dot{\Phi}_0(m)+\frac{t_1(m)}{4t_2(m)}\right)^2
-\frac{t_1^2(m)}{8t_2(m)}-t_2(m),
\label{eq:eheihou}
\end{align}
where
\begin{equation}
-\frac{t_1(m)}{4t_2(m)}=\frac{1-s}{2s(1-\lambda)\sqrt{1-m^2}}~(\ge 0).
\label{eq:extremum}
\end{equation}
The classical potential is derived as follows.  For a given fixed value of $m$, the energy $E$ can assume values corresponding to real values of $\dot{\Phi}_0$, the latter condition coming from  real values of the momentum for a classical trajectory.  As depicted in Fig. \ref{fig:potential}, we define the upper bound of the energy as $U_{+}$ (point c in Fig. \ref{fig:potential}) and the lower bound as $U_{-}$ (point b) at a given $m$ as obtained by changing $\cos\dot{\Phi}_0$ in the range $[-1,1]$ for fixed $m$. The energy can take values between these limits.  If we translate this fact into the classically allowed range of $m$ for a given $E$, the parameter $m$ can take values between  $U_{+}$ and  $U_{-}$ for fixed $E$, e.g. between points a and b as well as between d and e in Fig. \ref{fig:potential}.  The region between b and d is forbidden classically since the momentum becomes imaginary.
\begin{figure*}[thb]
	\centering
		\includegraphics[width=6cm,clip]{potential2.pdf}
	\caption{For a fixed value of $m$, the classical energy can take values between b and c.  This fact is translated into the range between a and b for a fixed value of $E$.  Also the range between d and e is allowed classically}
\label{fig:potential}
\end{figure*}
To find an explicit expression of the classical potential, we therefore have to minimize Eq. (\ref{eq:eheihou}) with respect to $\cos\dot{\Phi}_0$ in the range between $-1$ and 1.  It is necessary to consider the cases with $1-s>2s(1-\lambda)$ and $1-s<2s(1-\lambda)$ separately since the right hand side of Eq. (\ref{eq:extremum}) may or may not exceed 1 depending on the direction of the inequality. Before going further, we notice that the condition
\begin{eqnarray}
\label{eq:secondline}
1-s=2s(1-\lambda)
\end{eqnarray}
is known to give the second-order transition line in the $s$-$\lambda$ phase diagram \cite{Seki2012}. The second order phase transition thus corresponds to the border of classification of the classical potential.

When $1-s>2s(1-\lambda)$, the energy $E$ assumes its minimum at $\cos\dot{\Phi}_0=1$ according to Eqs. (\ref{eq:eheihou}) and (\ref{eq:extremum}). Then the classical potential is
\begin{align}
\nonumber
U_1(m)&:=w(m)+t_1(m)+t_2(m)\\
&=-s\lambda m^p-(1-s)\sqrt{1-m^2}+s(1-\lambda)(1-m^2).\label{potential1}
\end{align}
This coincides with the classical potential naively obtained by the replacement of $\hat{m}_z$ by $m$ and $\hat{m}_x$ by $\sqrt{1-m^2}$ in Eq. (\ref{operator_E}) as was done in Refs. \citen{Seoane2012,NishimoriTakada2017}.

The situation is slightly more complicated when $1-s<2s(1-\lambda)$. Suppose that the following relation holds,
\begin{eqnarray}
1\le-\frac{t_1(m)}{4t_2(m)}=\frac{1-s}{2s(1-\lambda)\sin\theta},\label{cond1}
\end{eqnarray}
which is rewritten in terms of $m(=\cos\theta)$ as
\begin{eqnarray}
|m|\ge\sqrt{1-\left(\frac{1-s}{2s(1-\lambda)}\right)^2}=: m_0
\end{eqnarray}
If $|m|>m_0$, i.e. if Eq. (\ref{cond1}) holds, the condition $\cos\dot{\Phi}_0(m)=1$ gives the minimum of the energy $E$ of Eq. (\ref{eq:eheihou}) as a function of $\cos\dot{\Phi}_0(m)$.  The resulting expression is the same $U_1(m)$ of Eq. (\ref{potential1}) as in the case of $1-s>2s(1-\lambda)$. On the other hand, if $|m|<m_0$, the minimum of $E$ is at $\cos\dot{\Phi}_0(m)=-t_1/4t_2$. Therefore the explicit expression of the classical potential for $1-s<2s(1-\lambda)$ is
\begin{align}
U(m)=\begin{cases}
	w(m)+t_1(m)+t_2(m)=U_1(m)  & (|m|\ge m_0) \\
	w(m)-\displaystyle\frac{t_1^2(m)}{8t_2(m)}-t_2(m)=:U_2(m) & (|m|\le m_0),
  \end{cases}
\end{align}
where $U_2(m)$ is expressed explicitly as
\begin{eqnarray}
\label{eq:U2}
U_2(m)=-s\lambda m^p-\frac{(1-s)^2}{4s(1-\lambda)}.
\end{eqnarray}
It is clear that the minimum of $U(m)$ as a function of $m$ is located in the range $|m|\ge m_0$ since $U_2(m)$ is a monotonically decreasing function of $m$. This means that we may consider only $U_1(m)$ as long as we are interested in the minimum of $U(m)$ as a function of $m$ to discuss the ground state. This observation justifies the naive classical analysis using Eq. (\ref{potential1}) to draw the phase diagram in Refs. \citen{Seoane2012,NishimoriTakada2017}.

Analyses of the potential $U_1(m)$ to draw the phase diagram have already been carried out in detail in Ref. \citen{Seoane2012}, and we just quote the result. Examples of the phase diagram are shown in Fig. \ref{fig:phasediagrams} for $p=3, 5$ and 11. The latter two cases have a line of second-order transitions for smaller $\lambda$ that replaces the first-order line for larger $\lambda$.  This means that we can avoid an exponentially small gap at a first-order transition by choosing an appropriate path in the phase diagram from $(s=0, \lambda)$ to $s=\lambda=1$.
\begin{figure*}[thb]
\centering
    \begin{tabular}{c}
      \begin{minipage}{0.33\linewidth}
			\centering
          \includegraphics[width=5cm,clip]{p3souzu2.pdf}
          \hspace{1.6cm} $p=3$
      \end{minipage}
      \begin{minipage}{0.33\linewidth}
        \begin{center}
          \includegraphics[width=5cm,clip]{p5souzu2.pdf}
          \hspace{1.6cm} $p=5$
        \end{center}
      \end{minipage}
      \begin{minipage}{0.33\linewidth}
        \begin{center}
          \includegraphics[width=5cm,clip]{p11souzu.pdf}
          \hspace{1.6cm} $p=11$
        \end{center}
      \end{minipage}
    \end{tabular}
    \caption{Phase diagram on the $s$-$\lambda$ plane for $p=3$ (left), $p=5$ (middle), and $p=11$ (right). The red bold curve is a line of first-order transitions between the quantum paramagnetic phase and the ferromagnetic phase. The blue dotted curve is a line of second order transitions.
}
\label{fig:phasediagrams}
\end{figure*}

\subsection{Exact location of the critical point}

It is possible to derive a closed-form expression of the critical point where the line of first-order phase transitions terminates within the ferromagnetic phase as in Fig. \ref{fig:phasediagrams}, seen most clearly for $p=11$. The Landau theory of phase transitions tells us that the derivatives of the free energy vanish up to third order at such a critical point \cite{Nishimori_book}. Let us denote the magnetization at the critical point as $m_*=\cos\theta_*$.  Then, the critical condition derived from the derivatives of the potential $U_1(\cos\theta)$ is
\begin{align}
\label{eq:firstder}
&ps\lambda\cos^{p-2}\theta_*\sin\theta_*-(1-s)+2(1-\lambda)s\sin\theta_*=0\\
\label{eq:secondder}
&ps\lambda\cos^{p-2}\theta_*\left(1-(p-2)\frac{\sin^2\theta_*}{\cos^2\theta_*}\right)+2s(1-\lambda)=0
\\
\label{eq:thirdder}
&p-7-(p-1)\cos2\theta_*=0.
\end{align}
Equation (\ref{eq:thirdder}) is solved as
\begin{equation}
m_*=\sqrt{\frac{p-4}{p-1}},
\end{equation}
which is the value of the magnetization at the critical point.
Equations (\ref{eq:firstder}) and (\ref{eq:secondder}) then give the critical values of $\lambda$ and $s$ as
\begin{align}
\label{eq:lambdacritical}
\lambda_*&=\frac{1}{1+pm_*^{p-4}}\\
\label{eq:scritical}
s_*&=\frac{1+pm_*^{p-4}}{1+pm_*^{p-4}\left(1+\sqrt{\frac{3}{p-1}}\left(m_*^2+2\right)\right)}.
\end{align}
These equations give the exact location of the critical point, which we have confirmed to agree with the numerical result for several values of $p$.

\subsection{Closed-form expression of the energy gap}

The energy gap at a first-order phase transition is given as the overlap of two wave functions $|\Psi_1\rangle$ and $|\Psi_2\rangle$ localized at two degenerate minima, $m=m_1$ and $m=m_2$, of the classical potential $U_1(m)$, according to the general semi-classical argument as described in Ref.  \citen{Bapst2012}. To calculate the overlap, we need an explicit formula for the wave functions $\phi_1(m)=\langle S, m|\Psi_1\rangle$ and $\phi_2(m)=\langle S, m|\Psi_2\rangle$.  Since $m_1$ and $m_2$ are the two degenerate minima of the classical potential at a first-order transition, the wave functions $\phi_1(m)$ and $\phi_2(m)$ should decay exponentially away from those minima. Also, the momentum $\dot{\Phi}_0$ should be imaginary for quantum tunneling to take place, and so we write $\Phi_0(m)=i\Psi (m)$ with real $\Psi (m)$. We therefore find the following expression of the wave function from Eqs. (\ref{Cm}) and (\ref{Phim}),
\begin{equation}
\phi_j(m)= \exp \left(-S \left|\int_{m_j}^m \dot{\Psi} (m')dm'\right|\right)\quad (j=1, 2).
\end{equation}
The coefficient $\alpha$ of the exponential decay of the gap $e^{-\alpha N}$ is then given in terms of the overlap as
\begin{align}
\alpha &= -\lim_{N\to\infty}\frac{1}{N} \ln \langle \Phi_1|\Phi_2\rangle \nonumber\\
&=-\lim_{N\to\infty}\frac{1}{N}\ln \int_{-1}^1 \phi_1(m) \phi_2(m)\, dm \nonumber\\
&=-\lim_{S\to\infty}\frac{1}{2S}\ln \int_{-1}^1 \exp \left(-S\left|\int_{m_1}^m \dot{\Psi}(m')dm'\right|
-S\left|\int_{m_2}^m \dot{\Psi}(m')dm'\right|\right)\, dm \nonumber\\
&=-\frac{1}{2}\sup_m \left\{- \left|\int_{m_1}^m \dot{\Psi}(m')dm'\right|- \left|\int_{m_2}^m \dot{\Psi}(m')dm'\right| \right\}.
\label{alpha_sup}
\end{align}
The function $\dot{\Psi}(m)$ is defined by Eq. (\ref{eq:ekinji}) with $\dot{\Phi}_0(m) =i\dot{\Psi}(m)$, which changes cosine to hyperbolic cosine in Eq. (\ref{eq:ekinji}). The inverse of hyperbolic cosine has two branches, and we choose the positive one, $\dot{\Psi}(m)\ge 0$. Then the supremum in Eq. (\ref{alpha_sup}) is achieved  when $m$ lies in the range $m_1<m<m_2$, where we assumed $m_1<m_2$. We therefore have an $m$-independent argument of supremum  \cite{Bapst2012},
\begin{equation}
-\int_{m_1}^m \dot{\Psi}(m')dm'-\int_{m}^{m_2} \dot{\Psi} (m')dm' =-\int_{m_1}^{m_2}\dot{\Psi}(m')dm'.
\end{equation}
Other possibilities of $m$ in Eq. (\ref{alpha_sup}) lead to a negative term additional to the above expression such as
\begin{equation}
-2\int_{m_2}^m \dot{\Psi}(m')dm'
\end{equation}
for $m>m_2$ and
\begin{equation}
-2\int_{m}^{m_1} \dot{\Psi}(m')dm'
\end{equation}
for $m<m_1$. Therefore the final expression of the coefficient is
\begin{align}
\alpha&=\frac{1}{2}\int^{m_2}_{m_1}\dot{\Psi}(m')dm',
\nonumber\\
\dot{\Psi}(m)&={\rm arccosh}\,\left(\frac{-t_1(m) -\sqrt{t_1^2(m)-8t_2(m)\left(w(m)-t_2(m)-E\right)}}{4t_2(m)}\right)\nonumber\\
&={\rm arccosh}\left(\frac{(1-s)-\sqrt{(1-s)^2-4s(1-\lambda)(-s\lambda m^p-E)}}{2s(1-\lambda)\sqrt{1-m^2}}\right).
\label{eq:arccosh}
\end{align}
Here, $E$ is the classical energy at the transition point,
\begin{eqnarray}
E=U_1(m_1)=U_1(m_2).
\end{eqnarray}
Equation (\ref{eq:arccosh}) for $\dot{\Psi}(m)$ has been derived from Eq. (\ref{eq:ekinji}) under the replacement $\dot{\Phi}_0=i\dot{\Psi}$, which changes $\cos (x)$ to $\cosh (x)$.  Physically, this means that quantum tunneling takes place in the classically forbidden range of the parameter $m_1<m<m_2$ (between points b and d in Fig. \ref{fig:potential}).

Values of $\alpha$ are plotted in Fig. \ref{fig:alpha-ABCD} fro $p=3, 5, 11,$ and 101.  Except for $p=3$, it is seen that $\alpha$ decreases mildly as $\lambda$ decreases from 1, until the behavior changes at the point where the first-order line meets the second-order line. Beyond this point, $\alpha$ decreases very rapidly toward 0. The dependence of $\alpha$ on $\lambda (>1/2)$ is weak for $p$ large and is close to 1/2, which is the value of $\alpha$ at $\lambda=1$ in the limit $p\to\infty$ \cite{Bapst2012}. The case of $p=3$ is an exception since the transition is always of first order for $0<\lambda \le 1$ \cite{NishimoriTakada2017}.
\begin{figure*}[h]
	\centering
		\begin{tabular}{cc}
			\begin{minipage}[t]{0.5\hsize}
				\centering
				\includegraphics[width=7cm]{p3alpha2.pdf}
				\hspace{1.6cm} $p=3$
			\end{minipage}
			\begin{minipage}[t]{0.5\hsize}
				\centering
				\includegraphics[width=7cm]{p5alpha2.pdf}
				\hspace{1.6cm} $p=5$
			\end{minipage}
\\
			\begin{minipage}[c]{0.5\hsize}
				\centering
				\includegraphics[width=7cm]{p11alpha2.pdf}
				\hspace{1.6cm} $p=11$
			\end{minipage}
			\begin{minipage}[c]{0.5\hsize}
				\centering
				\includegraphics[width=7cm]{p101alpha2.pdf}
				\hspace{1.6cm} $p=101$
			\end{minipage}
		\end{tabular}
\caption{The coefficient $\alpha$ of the energy gap closing rate along the first phase transition line on the $s$-$\lambda$ phase diagram. The horizontal axis represents $\lambda$, the right vertical axis is for $\alpha$ in a logarithmic scale, the left axis is $s$. The red bold line represents first-order phase transitions. The blue dash-dotted line denotes second-order phase transitions, $s=1/(3-2\lambda)$. The green dots are for $\alpha$.}
\label{fig:alpha-ABCD}
\end{figure*}

\section{Conclusion}

We have derived a closed-form expression of the energy gap at first-order quantum phase transitions in the $p$-spin model with antiferromagnetic transverse interactions. This is the first case, as far as we are aware of, that such an expression has been obtained for a non-stoquastic Hamiltonian.  The result reveals that the rate of gap closing depends on the system parameters $s$ and $\lambda$ relatively mildly until the first-order transition line hits the second-order line as the parameter $\lambda$ decreases to about 0.5.  After the two transition lines meet, the coefficient $\alpha$ becomes rapidly decreasing as $\lambda$ further decreases and eventually vanishes at a critical point within the ferromagnetic phase. These behaviors have been anticipated by direct numerical diagonalization of the Hamiltonian \cite{Seki2012}, and our contribution is to have derived the explicit analytical formula, which may lay a foundation for further developments.  The exact location of the critical point within the ferromagnetic phase has also been derived. These results can be taken as the final step to clarify the static properties of the non-stoquastic $p$-spin model from the perspective of quantum annealing.

It would be an interesting future direction to take into account finite temperatures to understand thermal effects on the performance of QA, possibly following the line layed out in Ref. \citen{Kechedzhi2016}.

\begin{acknowledgments}
This work was supported by the JSPS KAKENHI Grant No. 26287086.
\end{acknowledgments}

\bibliographystyle{jpsj.bst}
\input{ref1.bbl}
\end{document}