# Bandwidth controlled quantum phase transition between an easy-plane quantum
spin Hall state and an s-wave superconductor

Disha Hou  Department of Physics, Beijing Normal University, Beijing 100875,
China Yuhai Liu  Beijing Computational Science Research Center, 10 East
Xibeiwang Road, Beijing 100193, China School of Science, Beijing University of
Posts and Telecommunications, Beijing 100876, China Toshihiro Sato  Institut
für Theoretische Physik und Astrophysik, Universität Würzburg, 97074 Würzburg,
Germany Wenan Guo<EMAIL_ADDRESS>Department of Physics, Beijing Normal
University, Beijing 100875, China Beijing Computational Science Research
Center, 10 East Xibeiwang Road, Beijing 100193, China Fakher F. Assaad
<EMAIL_ADDRESS>Institut für Theoretische Physik und
Astrophysik, Universität Würzburg, 97074 Würzburg, Germany Würzburg-Dresden
Cluster of Excellence ct.qmat, Am Hubland, 97074 Würzburg, Germany Zhenjiu
Wang<EMAIL_ADDRESS>Max-Planck-Institut für Physik komplexer Systeme,
Dresden 01187, Germany

###### Abstract

The quantum spin Hall state can be understood in terms of spontaneous O(3)
symmetry breaking. Topological skyrmion configurations of the O(3) order
parameter vector carry a charge 2e, and as shown previously, when they
condense, a superconducting state is generated. We show that this topological
route to superconductivity survives easy-plane anisotropy. Upon reducing the
O(3) symmetry to O(2)$\times$ Z2, skyrmions give way to merons that carry a
unit charge. On the basis of large-scale auxiliary field quantum Monte Carlo
simulations, we show that at the particle-hole symmetric point, we can trigger
a continuous and direct transition between the quantum spin Hall state and
s-wave superconductor by condensing pairs of merons. This statement is valid
in both strong and weak anisotropy limits. Our results can be interpreted in
terms of an easy-plane deconfined quantum critical point. However, in contrast
to the previous studies in quantum spin models, our realization of this
quantum critical point conserves $U(1)$ charge, such that skyrmions are
conserved.

##  I Introduction

Topology is a key factor for understanding phase transitions. In the
Kosterlitz-Thouless theory Kosterlitz and Thouless (1973), an O(2) local order
parameter in two-dimensional space allows for the definition of the vortex,
the proliferation of which drives the transition. Let us stay in two-
dimensional space, $\bm{x}=(x,y)$, but consider an O(3) local order parameter
$\bm{n}(\bm{x})$ with the unit norm. This combination of space and order
parameter defines a winding number, Fradkin (2013a)

$\frac{1}{4\pi}\int
d^{2}\bm{x}\,\bm{n}\cdot\partial_{x}\bm{n}\times\partial_{y}\bm{n}.$ (1)

For smooth configurations, this quantity is quantized and counts the winding
of the unit vector on the unit sphere: a skyrmion. We can now reduce the O(3)
symmetry to O(2) $\times$ Z2. In the context of spin systems, this would
correspond to restricting the O(3) symmetry to O(2) transformations around,
say, the z-axis and a change of sign of the third component of the
$\bm{n}$-vector. Assuming that energetics favor the $\bm{n}$-vector to be in-
plane (i.e., vanishing z-component) then the topological excitations will
correspond to vortices in the x-y plane. Due to the normalization condition,
the $\bm{n}$-vector at the core of the vortex will have to point along the
z-direction. Since the $\bm{n}$-vector lies in the x-y plane at infinity, the
integrand of Eq. 1 vanishes at infinity and the integral takes half-integer
values: a meron.

The above considerations acquire different interpretations depending upon the
specifics of the local order parameter. In this paper, it corresponds to the
order-parameter of the quantum-spin Hall state. In particular, let
$\hat{H}_{0}=-v_{F}\sum_{\bm{p},\sigma,i=1,2}\hat{\bm{\Psi}}^{\dagger}_{\sigma}(\bm{p})ip_{i}\gamma_{0}\gamma_{i}\hat{\bm{\Psi}}_{\sigma}(\bm{p})$
be the Hamiltonian akin to graphene in the absence of interactions Neto  _et
al._ (2009), using the notation of Ref. Herbut  _et al._ , 2009. Inclusion of
the quantum spin Hall mass term leads to:

$\hat{H}=\hat{H}_{0}+\int_{V}d\bm{x}\bm{N}(\bm{x})\cdot\hat{\bm{\Psi}}^{\dagger}_{\sigma}(\bm{x})i\gamma_{0}\gamma_{3}\gamma_{5}\bm{\sigma}_{\sigma,\sigma^{\prime}}\hat{\bm{\Psi}}_{\sigma^{\prime}}(\bm{x}).$
(2)

The order parameter $\bm{N}$ can be normalized to unity
$\bm{n}=\bm{N}/|\bm{N}|$ if, as will be the case in our model, the single-
particle gap, that is proportional to the length of $\bm{N}$, does not vanish.
Furthermore, $\bm{n}$ is odd under charge conjugation such that the expression
in Eq. (1) carries the same quantum numbers as the charge density $\rho$
measured with respect to half-filling. In particular, in Ref. Grover and
Senthil, 2008 it is shown that

$\rho(\bm{x})=\frac{2e}{4\pi}\bm{n}\cdot\partial_{x}\bm{n}\times\partial_{y}\bm{n}$
(3)

such that the skyrmion (meron pairs) carry charge 2e. Since topological
defects of one phase carry the charge of the other, their proliferation will
lead to a symmetry-broken state. In the above discussion, we have considered
the quantum spin Hall (QSH) state where the skyrmions (or pairs of merons)
carry charge 2e and their proliferation leads to a superconducting state.
Alternatively, one could consider the three spin-density wave mass terms. In
this case, the skyrmion creation/annihilation process would acquire a phase
under the spatial rotation of the lattice. This proliferation of skyrmions is
the essence of so-called ‘deconfined’ quantum critical points (DQCPs) Senthil
_et al._ (2004a, b); Wang  _et al._ (2017).

The above discussion takes place in the continuum, and some type of
regularization is needed to carry out numerical simulations. Starting with
magnetic phases, where the topological defects carry a U(1) rotational charge,
a lattice regularization leads to additional symmetry breaking terms that may
lead to subtleties since they do not exist in the expected IR field theory
Sandvik (2007); Shao  _et al._ (2016). For instance for a square lattice
regularization scheme, U(1) rotational symmetry gives way to a
$\mathbb{Z}_{4}$ invariance of the valence bond solid (VBS) state. This
regularization-induced symmetry reduction introduces novel operators that have
to be argued to be irrelevant at the critical point. In particular, in the
framework of the CP1 field theory of spinon coupled to a U(1) gauge field
Senthil  _et al._ (2004a), this symmetry reduction allows for the creation of
quadruple monopoles of the gauge field.

Although the $Z_{4}$ lattice symmetry-breaking field is relevant in the VBS
state, a necessary condition for the continuous nature of the transition into
an xy-antiferromagnetic (AFM) phase is that this symmetry-breaking field is
irrelevant at the critical point. Numerically, AFM-VBS phase transitions in
the easy-plane case have clear first-order signatures in most cases Kuklov
_et al._ (2006); D’Emidio and Kaul (2016); Sandvik  _et al._ (2002); Kragset
_et al._ (2006); Sen  _et al._ (2007); D’Emidio and Kaul (2017); Desai and
Kaul (2020). Among numerical works, Desai et. al. Desai and Kaul (2020), in p
emphasize the absence of the continuous transition in any easy-plane spin
system: the authors generally claimed the absence of the easy-plane deconfined
fixed point without considering the effect of quadruple monopoles. However,
could it be that the $Z_{4}$ symmetry-breaking field introduces a runaway
flow, leading to a first-order transition?

Instead of encoding the U(1) symmetry as a rotational invariance, one can
encode it in terms of charge conservation. Importantly, charge conservation
will not be broken by lattice regularization. Following the work of Ref. Liu
_et al._ , 2019; Wang  _et al._ , 2021 we set up a set of designer
Hamiltonians which have the potential to realize an easy-plane DQCP without
quadruple monopoles. A dynamically generated QSH insulating state which
spontaneously breaks the O(2) spin rotational symmetry emerges from a Dirac
semi-metal via a spin-orbital interaction. Our particular interest lines in
the phase transition between the QSH and s-wave superconducting (SSC) states.

The fermion basis introduces a simple but much more provoking picture: meron
defects of the QSH order parameter which carry a unit of electron charge are
the fundamental excitation at the critical point; on the other side of the
transition, the condensation of meron-pair creation/annihilation operators
forms the superconducting state. Importantly, the U(1) charge conservation
broken by the SSC phase is an exact symmetry of our lattice Hamiltonian,
meaning that quadruple monopoles are absent by definition.

The aim of this work is to systematically search for the existence of an easy-
plane DQCP without monopoles. The continuity of phase transitions does not
only depend on symmetry. Domany  _et al._ (1984); Blöte  _et al._ (2002) In
our lattice model, a model parameter that is related to the strength of the
easy-plane anisotropy continuously tunes the energy gap of meron
configurations of a QSH order parameter. Regardless of the strength of the
anisotropy and in contrast to lattice spin models, our model shows no obvious
signs of first-order transitions. We argue that this transition flows to the
easy-plane DQCP.

The paper is organized as follows. In Sec. II we introduce our lattice
Hamiltonian, the quantum Monte Carlo algorithm, as well as basic observables.
The numerical results are shown in Sec. III, beginning with the ground state
phase diagram and followed by a detailed investigation of the nature of the
phase transitions. Finally, we draw conclusions and give an outlook in Sec.
IV.

##  II Model and Methods

We consider a model of Dirac fermions in $2+1$ dimensions on the honeycomb
lattice with Hamiltonian

$\displaystyle\hat{H}_{t}=-t\sum_{\langle\bm{i},\bm{j}\rangle}(\hat{\bm{c}}^{\dagger}_{\bm{i}}\hat{\bm{c}}^{\phantom{\dagger}}_{\bm{j}}+H.c.).$
(4)

Here, the spinor
$\hat{\bm{c}}^{{\dagger}}_{\bm{i}}=\big{(}\hat{c}^{{\dagger}}_{\bm{i},\uparrow},\hat{c}^{{\dagger}}_{\bm{i},\downarrow}\big{)}$
where $\hat{c}^{{\dagger}}_{\bm{i},\sigma}$ creates an electron in a Wannier
state centered around lattice site $\bm{i}$ with $z$-component of spin
$\sigma$. This term accounts for nearest-neighbor hopping. The interaction
term that we consider reads:

$\displaystyle\hat{H}_{\lambda}=$
$\displaystyle-\lambda\sum_{\varhexagon}\left[\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{x}_{\bm{i},\bm{j}}\right)^{2}+\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{y}_{\bm{i},\bm{j}}\right)^{2}\right.$
(5)
$\displaystyle\phantom{=\;\;}\left.+\Delta\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{z}_{\bm{i},\bm{j}}\right)^{2}\right]$

where $\hat{\bm{J}}_{\bm{i},\bm{j}}\equiv
i\nu_{\bm{i}\bm{j}}\hat{\bm{c}}^{\dagger}_{\bm{i}}\bm{\sigma}\hat{\bm{c}}^{\phantom{\dagger}}_{\bm{j}}+H.c.$
The components of $\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ are the
Pauli spin-1/2 matrices. This term is a plaquette interaction involving next-
nearest-neighbor pairs $\langle\langle\bm{i}\bm{j}\rangle\rangle$ of sites and
phase factors $\nu_{\bm{ij}}=\pm 1$ identical to the Kane-Mele model Kane and
Mele (2005), see also Ref. Liu  _et al._ , 2019.

The Hamiltonian $\hat{H}=\hat{H}_{t}+\hat{H}_{\lambda}$ with $\Delta=1$ has
been studied in Ref. Liu  _et al._ , 2019. A dynamically generated QSH
insulator that breaks $SU(2)$ spin rotational symmetry spontaneously is found
at intermediate interacting strength ($\lambda$), separating a Dirac semi-
metal(DSM) state at small $\lambda$ and an SSC state at large $\lambda$. The
DSM-QSH transition belongs to the Gross-Neveu Heisenberg universality class
Gross and Neveu (1974) whereas the QSH-SSC transition falls into the class of
DQCP. In the current work, we focus on the case of $\Delta\in[0,1)$ where the
$SU(2)$ spin rotational symmetry is reduced to $U(1)\times Z_{2}$.

We used the ALF (Algorithms for Lattice Fermions) implementation Assaad  _et
al._ (2022) of the well-established auxiliary-field quantum Monte Carlo (QMC)
method Blankenbecler  _et al._ (1981); White  _et al._ (1989); Assaad and
Evertz (2008). Because $\lambda>0$ and $\Delta>0$, we can use a real Hubbard-
Stratonovich decomposition for the perfect square term. We set the imaginary
time interval $\Delta\tau=0.2$ and choose a symmetric Trotter decomposition to
ensure the hermiticity of the imaginary time propagation in the Monte Carlo
simulations.Liu  _et al._ (2019) Additionally, a checkerboard decomposition is
applied to the exponential of hopping matrix $\hat{H}_{t}$. For each field
configuration, time-reversal symmetry and charge are conserved. Hence the
eigenvalues of the fermion determinant occur in complex conjugate pairs, and
we do not suffer from the negative sign problem Wu and Zhang (2005). We
simulated lattices with $L\times L$ unit cells (each containing two Dirac
fermions) and periodic boundary conditions. Following our previous work Wang
_et al._ (2021), we used a projective version of the algorithm (PQMC) Sugiyama
and Koonin (1986); Sorella  _et al._ (1989); Assaad and Evertz (2008). This
algorithm is based on the form:

$\langle\hat{O}\rangle=\lim_{\theta\rightarrow\infty}\frac{\left<\Psi_{T}|e^{-\theta\hat{H}}\hat{O}e^{-\theta\hat{H}}|\Psi_{T}\right>}{\left<\Psi_{T}|e^{-2\theta\hat{H}}|\Psi_{T}\right>}.$
(6)

Provided that the trial wave function $\left.|\Psi_{T}\right>$ is not
orthogonal to the ground state, $\langle\hat{O}\rangle$ corresponds to the
ground state expectation value of the observable $\hat{O}$. To avoid the
negative sign problem, we consider the same time-reversal symmetric trial wave
function as the one used in Ref. Wang  _et al._ , 2021. We explicitly checked
the projection convergence to the ground state at each system size and each
$\Delta$: simulations are performed at $L=6,9,12,15,18$ and the values of
projection length $\theta$ for ground state calculation are listed in Tab. 1.

Table 1:  Projection length $\theta$ at different values of $\Delta$ and $L$.  lattice size | anisotropy strength
---|---
$\Delta=0.1$ |  $\Delta=0.5$ |  $\Delta=0.75$
$L=6$ | 15 | 15 | 15
$L=9$ | 21 | 21 | 21
$L=12$ | 24 | 24 | 24
$L=15$ | 42 | 36 | 36
$L=18$ | 42 | 36 | 36

The basic measurements in our QMC simulations are equal time correlation
functions in real space:

$\displaystyle
S^{O}(\bm{r})\equiv\frac{1}{L^{2}}\sum_{n}\sum_{\bm{r}^{\prime}}\langle\hat{O}^{\dagger}_{\bm{r},n}\hat{O}^{\phantom{\dagger}}_{\bm{r}+\bm{r}^{\prime},n}\rangle,$
(7)

and the structure factor:

$\displaystyle
S^{O}_{m,n}(\bm{q})\equiv\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}\langle\hat{O}^{\dagger}_{\bm{r},m}\hat{O}^{\phantom{\dagger}}_{\bm{r}^{\prime},n}\rangle,$
(8)

where $\hat{O}_{\bm{r},n}$ is a local operator with $\bm{r}$ denoting the unit
cell and $n$ denoting the intra unit-cell dependence that we will refer to as
orbital.

For instance, the spin-orbit coupling operators correspond to
$\hat{\bm{O}}_{\bm{r},n}=\hat{\bm{J}}_{\bm{r}+\bm{\delta}_{n},\bm{r}+\bm{\eta}_{n}}$.
Here $n$ runs over the six next-nearest neighbor bonds of the corresponding
hexagon with legs $\bm{r}+\bm{\delta}_{n}$ and $\bm{r}+\bm{\eta}_{n}$
($n=1,2,...6$), see Fig. 1. To detect QSH ordering which breaks the $O(2)$
spin rotational symmetry, we calculate the structure factor matrix associated
with the $X$ and $Y$ components of $\hat{\bm{O}}_{\bm{r},n}$

$\displaystyle
S^{\text{QSH}}_{m,n}(\bm{q})\equiv\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}\langle\hat{O}^{X}_{\bm{r},m}\hat{O}^{X}_{\bm{r^{\prime}},n}+\hat{O}^{Y}_{\bm{r},m}\hat{O}^{Y}_{\bm{r^{\prime}},n}\rangle,$
(9)

with $m,n=1,2,...6$. We also consider the structure factor matrix associated
with the $Z$ component of $\hat{\bm{O}}_{\bm{r},n}$

$\displaystyle
S^{\text{QSH}_{Z}}_{m,n}(\bm{q})\equiv\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}\langle\hat{O}^{Z}_{\bm{r},m}\hat{O}^{Z}_{\bm{r^{\prime}},n}\rangle.$
(10)

The physical meaning of this quantity will be discussed in Sec. III.

Figure 1:  Honeycomb plaquette illustrating (a) the spin-orbit coupling
operator $\hat{\bm{J}}_{\bm{r}+\bm{\delta}_{n},\bm{r}+\bm{\eta}_{n}}$ and (b)
the pairing operator $\hat{\eta}^{+}_{\bm{r},\bm{\tilde{\delta}_{a}}}$.

To detect SSC ordering which breaks $U(1)$ charge conservation, we consider
the following structure factor matrix:

$\displaystyle
S^{\text{SSC}}_{a,b}(\bm{q})\equiv\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}[\langle\hat{\eta}^{+}_{\bm{r},\bm{\widetilde{\delta}}_{a}}\hat{\eta}^{-}_{\bm{r^{\prime}},\bm{\widetilde{\delta}}_{b}}\rangle+\langle\hat{\eta}^{-}_{\bm{r},\bm{\widetilde{\delta}}_{a}}\hat{\eta}^{+}_{\bm{r^{\prime}},\bm{\widetilde{\delta}}_{b}}\rangle],$
(11)

with $a,b=1,2$, denoting the A(B) sublattice, where the s-wave pairing
operator is defined as

$\hat{\eta}^{+}_{\bm{r},\bm{\tilde{\delta}_{a}}}=\hat{c}^{\dagger}_{\bm{r}+\bm{\tilde{\delta}_{a}},\uparrow}\hat{c}^{\dagger}_{\bm{r}+\bm{\tilde{\delta}_{a}},\downarrow}.$
(12)

Here $\tilde{\delta}_{a}$ runs over the two orbitals in unit cell $\bm{r}$,
see Fig. 1.

We use Eq. 9 and Eq.11 to calculate the order parameter:

$m^{O}=\sqrt{\frac{\Lambda_{1}(S^{O}(\bm{Q}))}{L^{2}}}$ (13)

Here, $\Lambda_{1}()$ indicates the largest eigenvalue of the corresponding
matrix in orbital space, $O$ denotes QSH and SSC order parameters, and
$\bm{Q}=(0,0)$. The corresponding eigenvector will determine the orbital
structure. This is of particular importance for the QSH state since we expect
it to reflect the sign structure $\nu_{\bm{i}\bm{j}}$ of the Kane-Mele model.

To locate the critical points and study the critical properties, after
diagonalizing the corresponding structure factors, we calculated the
renormalization-group invariant correlation ratio

$R^{O}\equiv 1-\frac{S^{O}(\bm{Q}+\Delta\bm{q})}{S^{O}(\bm{Q})}$ (14)

using the largest eigenvalue $S^{O}$ ($O$ = QSH, QSHz, SSC); $\bm{Q}=(0,0)$ is
the ordering wave vector and $\bm{Q}+\Delta\bm{q}$ is a neighboring wave
vector with $|\Delta\bm{q}|=\frac{4\pi}{\sqrt{3}L}$. By definition, $R^{O}\to
1$ for $L\to\infty$ in the corresponding ordered state, whereas $R^{O}\to 0$
in the disordered state. At the critical point, $R^{O}$ is scale-invariant for
sufficiently large $L$ so that the results for different system sizes cross.

##  III Quantum Monte Carlo results

In this section we will first provide the ground state phase diagram and then
will proceed to investigate the nature of the phase transitions.

###  III.1 Ground state phase diagram

As mentioned previously, we are interested in the parameter range of
$\Delta\in[0,1)$ where the spin rotational symmetry of the Hamiltonian
$\hat{H}=\hat{H}_{t}+\hat{H}_{\lambda}$ is lowered to $U(1)\times Z_{2}$. The
DSM and SSC states found in the $SU(2)$ symmetric case Liu  _et al._ (2019) (
$\Delta=1$ ) are naturally stable against weak easy-plane anisotropy
$\Delta\approx 1$ since both states are spin rotational invariant.
Furthermore, since time reversal symmetry is not broken by our symmetry
reduction, we expect the QSH phase to be equally stable. To confirm the above,
we can use the mean-field approach introduced in Ref. Wang  _et al._ , 2021
that carries over to the anisotropic case. Due to the Dirac nature of the
kinetic term in Eq. 4, we foresee the robustness of the DSM phase in the
weakly interacting case. On the other hand, the attractive nature of
$\hat{H}_{\lambda}$ term ( for $\Delta\geq 0$, $\lambda>0$ ) suggests that the
mean-field picture in the large $\lambda$ case Wang  _et al._ (2021) will
still favor an SSC instability. Finally, the dynamically generated QSH state
at intermediate values of $\lambda$ will be restricted to the $U(1)$ plane in
the current case. We present the mean-field phase diagram in Fig. 2. The
details of the calculations are summarized in Appendix A. It is worth
mentioning that, in the mean-field analysis, the QSH and SSC orderings coexist
in a large $\lambda$ region of the phase diagram, see Fig. 2. This is a
natural consequence of the anti-commuting nature of the two Dirac masses. In
this case, the fermion band gap is given by the norm of the four-component
order parameter accounting for the QSH ($m_{\text{QSH}}$) and SSC
($m_{\text{SSC}}$) orders:

$\Delta_{\text{MF}}\propto\sqrt{m^{2}_{\text{QSH}}+m^{2}_{\text{ SSC}}}.$ (15)

Hence, developing superconducting ordering in the background of a QSH state
can simply minimize the mean-field free energy Wang  _et al._ (2021) in the
large  $\lambda$ limit.

Figure 2:  Mean-field phase diagram at zero temperature as a function of
$\Delta$ and $\lambda$.  Figure 3:  Ground state phase diagram in the
$(\Delta,\lambda)$ plane which covers the range of $\Delta\in[0.1,1]$. The
DSM-QSH and QSH-SSC phase boundaries at $\Delta=0.1,0.5$, and $0.75$ are
estimated from PQMC simulations in this paper. The phase boundaries of $SU(2)$
symmetric Hamiltonian (the overline of $\Delta=1$ ) are from Ref Liu  _et al._
(2019).

We summarize the exact ground state phase diagram in $(\Delta,\lambda)$ plane
based on QMC results in Fig. 3, where the overline $\Delta=1$ corresponds to
the $SU(2)$ symmetric model that was studied in Ref.Liu  _et al._ , 2019.
Generally speaking, we found a DSM state at weak interaction (small $\lambda$
region), an SSC state at strong interaction (large $\lambda$ region), as well
as a $U(1)$ broken QSH state at an intermediate region.

The numerical simulations that we performed cover three horizontal lines as a
function of $\lambda$: $\Delta=0.1$, $0.5$ and $0.75$. As shown in Fig. 4(b),
Fig. 5(b), and Fig. 6(b), the QSH correlation ratio $R^{QSH}$ increases toward
$1$ at intermediate values of $\lambda$, indicating a robust QSH state for all
three cases. We have checked explicitly that the eigenvector corresponding to
the largest eigenvalue matches the sign structure $\nu_{\bm{i}\bm{j}}$ of the
Kane-Mele model. Furthermore, in Appendix D we have used the flux-insertion
scheme presented in Ref. Assaad  _et al._ , 2013 to probe for the topological
invariant. The DSM-QSH and QSH-SSC phase boundaries are estimated by fitting
the equal-time correlation ratios of the two largest lattice sizes $L=15$ and
$L=18$ according to the following function

$f(L,\lambda)=R_{c}+a_{1}(\lambda-\lambda_{c})L^{1/\nu}+a_{2}(\lambda-\lambda_{c})^{2}L^{2/\nu}.$
(16)

The details of the fit are listed in Tab. 2 and Tab. 3. Remarkably, the
critical values of $\lambda$ where superconducting order develops, as shown by
the crossing points of $R^{SSC}$ in Fig. 4(a), Fig. 5(a), and Fig. 6(a), match
the values of $\lambda$ where the QSH order vanishes. This indicates that,
regardless of the strength of anisotropy, direct phase transitions exist
between the QSH and SSC phases. The order parameters give consistent results,
as shown in Fig. 7.

The comparison of the mean-field, Fig. 2, and QMC, Fig. 3, phase diagrams are
very instructive. The transition from the DSM to QSH insulator at $\Delta<1$
belongs to the U(1) Gross Neveu universality class. The essence of this
transition, a symmetry-breaking induced electronic mass generation, is
captured at the mean-field level. In fact, an $\epsilon$-expansion around the
upper critical dimension accounts rather well for this transition for the
SO(3) Assaad and Herbut (2013); Liu  _et al._ (2021) and U(1) cases Otsuka
_et al._ (2018). It is hence not unexpected that the comparison between the
mean-field, Fig. 2, and QMC, Fig. 3, phase diagrams is good for this
transition. In contrast, the competition and interplay between the QSH and SSC
phases are radically different at the mean-field and QMC levels. We interpret
this mismatch as a hint that topology – not accounted for at the mean-field
level – is crucial for the understanding of the intertwinement of the QSH and
SSC phases. Of particular importance is that the QMC phase diagram does not
show a coexistence of the QSH and SSC phases. The nature of the transition
will be discussed in the next section.

Table 2:  DSM-QSH crossing points $\lambda_{c}$. anisotropy | $\lambda_{c}$ | $\chi_{r}^{2}$ | $O$
---|---|---|---
$\Delta=0.1$ | 0.0243(2) | 0.14 | QSH
$\Delta=0.5$ | 0.0225(2) | 2.7 | QSH
$\Delta=0.75$ | 0.0208(3) | 0.81 | QSH
Table 3:  QSH-SSC crossing points $\lambda_{c}$. anisotropy | $\lambda_{c}$ | $\chi_{r}^{2}$ | $O$
---|---|---|---
$\Delta=0.1$ | 0.06006(8) | 0.08 | QSH
| 0.05988(2) | 2.41 | SSC
$\Delta=0.5$ | 0.0448(2) | 19.87 | QSH
| 0.04444(4) | 3.77 | SSC
$\Delta=0.75$ | 0.03829(4) | 1.67 | QSH
| 0.03788(3) | 3.44 | SSC

Our QSH insulator at zero temperature is also a gapless phase, reflecting the
emergence of Goldstone modes upon breaking the global XY symmetry. Therefore,
both the spin current operators $\hat{J}^{x,y}$ and the corresponding angular
momentum operator $\hat{S}^{z}$ reveal gapless excitations around zero
momentum. On the other hand, merons in this phase are another low-energy
excitation with a finite gap. Roughly speaking, the binding energy of pairs of
merons, is higher in the case of strong anisotropy. The numerical evidence
that merons bind can be deduced from the spectral functions presented in
Appendix C. A comparison between the single particle and superconducting
spectral functions shows that the cost of adding a pair is less than twice the
single-particle gap. This binding energy can be tuned to zero by increasing
the interaction, thus triggering a direct transition to a superconducting
state.

Figure 4:  Equal-time correlation ratio $R^{\rm SSC}$ (a) and $R^{\rm QSH}$
(b) as a function of $\lambda$ for $\Delta=0.1$.  Figure 5:  Same as Fig. 4
for $\Delta=0.5$.  Figure 6:  Same as Fig. 4 for $\Delta=0.75$.  Figure 7: QSH
and SSC order parameter as a function of $\lambda$, for $\Delta=0.1$ ((a) and
(b)), $0.5$ ((c) and (d)) and $0.75$ ((e) and (f)).

###  III.2 Nature of the QSH-SSC phase transition

Our most important result is the seemingly continuous nature of the QSH-SSC
transition.

Figure 8:  Single particle gap $\Delta_{sp}$ across the QSH-SSC transition for
$\Delta=0.1$ (a), $\Delta=0.5$ (b), and $\Delta=0.75$ (c), respectively.

The QSH-SSC transition has only bosonic excitations at low energy. To
characterize this, we extrapolate the fermionic single particle gap
$\Delta_{sp}$ from Green’s function:

$\displaystyle\sum_{\sigma}\langle
c_{\bm{k},\sigma}(\tau)c^{\dagger}_{\bm{k},\sigma}(0)\rangle\propto
e^{-\Delta_{\text{sp}}\tau},$ (17)

at $\bm{k}=M\equiv(0,\frac{2\pi}{\sqrt{3}}).$ Figure. 8 demonstrates that
$\Delta_{sp}$ remains nonzero across the QSH-SSC transition at $\lambda_{c}$
for all three considered values of the anisotropy.

In the following step, we inquire whether the QSH-SSC transition is continuous
or not. Considering that the computational cost of the AFQMC scales as $\theta
L^{3}$ with $\theta$ being the projection length, we do not apply the commonly
used method for detecting first-order transitions, such as analyzing the
finite-size behavior of the histogram of order parameters or the behavior of
the Binder cumulant. Instead, we study the correlation length to reflect the
nature of the phase transitions. A continuous phase transition is
characterized by a diverging correlation length in the thermodynamic limit:
$\xi\propto|\lambda-\lambda_{c}|^{-\nu}$. On the other hand, for a first-order
transition, the correlation length saturates to a finite value. We use the
real-space, equal-time correlation functions $S^{O}(\bm{r})$ of the order
parameter to define the correlation lengthParisen Toldin  _et al._ (2015)

${\xi}^{O}\equiv\sqrt{\frac{\sum_{\bm{r}}|\bm{r}|^{2}S^{O}(\bm{r})}{\sum_{\bm{r}}S^{O}(\bm{r})}}$
(18)

where $S^{O}(\bm{r})$ is defined in Eq. 7. For continuous symmetry breaking,
an issue with this definition is that it picks up the correlation lengths
along both the longitudinal and the transverse directions. Hence without a
specific symmetry-breaking pinning field to resolve the longitudinal
direction, the correlation length $\xi$ in Eq. (18) is well defined only in
the disordered state or at the critical point. Therefore, we discard the data
for the QSH (SSC) correlation length in the QSH (SSC) state.

As depicted in Fig. 9(a), Fig. 10(a) and Fig. 11(a), around transition points
at different values of $\Delta$, the correlation length $\xi$ of both the QSH
and SSC order parameters grow with system size $L$, without any tendency of
saturation. We define the scaled correlation length $\xi/L$ as the ratio
between the correlation length and system size. As shown in Fig. 9(b), Fig.
10(b), and Fig. 11(b), for both the QSH and SSC correlation lengths, the
ratios for different system sizes cross at the same point, suggesting that the
correlation length diverges with $L$. The divergence of correlation lengths
indicates the continuous nature of the phase transition.

It is worth mentioning that we observe an amazing match of the value of the
scaled correlation length $\xi/L$ for different anisotropy strengths $\Delta$.
The fact that the same value of $\xi/L$ at the three transition points implies
that all three QSH-SSC phase transition points correspond to the same fixed
point. $\xi/L$ is a dimensionless quantity which is a renormalization-group
invariant at the critical point. This number is claimed to be universal in a
$(2+1)$ dimensional system with conformal invariance. On a lattice system,
this universal number is not only pinned by the scaling dimension of the order
parameter but also by the microscopic couplings in different directions, the
boundary conditions, and the shape of the system (e.g. the aspect ratio).
Kamieniarz and Blote (1993); Qian  _et al._ (2016) In our case, the only
difference between the three different values of $\Delta$ is the intrinsic
spin anisotropy which is not related to the lattice geometry, such that
$\xi/L$ should be universal if all three transitions at $\Delta=0.1,0.5$ and
$0.75$ belong to the same universality class. We also observe the interesting
behavior that QSH and SSC operators cross at the same value of $\xi/L$ at the
transition point. This also indicates the identical value of the anomalous
dimension $\eta$ between the two order parameters, which is significant
evidence of an emergent O(4) symmetry Wang  _et al._ (2017); Sato  _et al._
(2017).

Figure 9:  Correlation length (a) and scaled correlation length (b) in
disordered phase for $\Delta=0.1$.  Figure 10:  Same as Fig. 9 for
$\Delta=0.5$.  Figure 11:  Same as Fig. 9 for $\Delta=0.75$.

We also calculate the first-order partial derivative of the free energy
density with respect to the coupling strength $\lambda$,

$\frac{\partial
f}{\partial\lambda}=\frac{1}{L^{2}}\frac{1}{\lambda}\langle\hat{H}_{\lambda}\rangle,$
(19)

to study the nature of the phase transition. This approach requires no
information of the order parameter (and the associated symmetry breaking). In
the case of a first order transition, when the system size is much larger than
correlation length $\xi$, one expects a discontinuity in this derivative at
the transition point.

Figure 12:  Free-energy derivative ${\partial f/\partial\lambda}$ as a
function of $\lambda$ in the vicinity of the QSH-SSC transition point, for
$\Delta=0.1$ (a), $\Delta=0.5$ (b), and $\Delta=0.75$ (c), respectively.

Figure 12 shows $\partial f/\partial\lambda$ as a function of $\lambda$ in the
vicinity of the QSH-SSC transition point for $\Delta=0.1$, $0.5$ and $0.75$.
Our data reveal no clear signs of a jump for the accessible system sizes in
this study. This is consistent with the aforementioned correlation length
analyses.

In the strong anisotropic case $\Delta=0.1$, the slope of the curve scales up
moving towards the transition point when increasing the system size. This
result may be interpreted as the signature of a ‘weakly first-order’
transition. Around a continuous transition point, the derivative of the free
energy scales as

$\frac{\partial f}{\partial\lambda}\propto|\lambda-\lambda_{c}|^{(d+z)\nu-1}$
(20)

in the thermodynamic limit. For a ‘weakly first-order’ transition, one expects
a ‘pseudo critical’ phenomenon where the $\nu(L)$ estimated from finite sizes
would approach $1/(d+z)$ as $L$ reaches $\xi$,Iino  _et al._ (2019) such that
$\partial f/\partial\lambda$ asymptotically shows a jump. However, we won’t be
able to conclude the nature of transition at $\Delta=0.1$ since it is not
clear whether the finite-size slope diverges or saturates upon approaching the
thermodynamic limit. On the other hand, the robustness of the slope in the
case of $\Delta=0.5$ and $\Delta=0.75$ as shown in Fig. 12(b) and (c),
indicates clear continuous phase transitions, unless there exists a non-
diverging $\xi$ that is significantly larger than $L$.

Figure 13:  Equal-time correlation ratio $R^{QSHz}$ for $\Delta=0.1$ (a),
$\Delta=0.5$ (b), and $\Delta=0.75$ (c), respectively.

We now consider the Z-component of the QSH correlations, Eq. 10. In the DSM
and at the DSM to QSH transition the single particle gap vanishes such that
this quantity is expected to decay as a power-law. In particular, in the DSM
the scaling dimension of the fermion operator is given by $\frac{d}{2}$ ($d=2$
is the dimensionality) such that the Z-QSH as well as XY-QSH correlation
functions are expected to decay as $r^{-4}$ with $r$ being the distance. At
the DSM to QSH transition the scaling dimensions of both quantities will
differ. For equal time correlations, and in $d=2$, power-law decay leads to a
divergence in the structure factor provided that the scaling dimension of the
operator is smaller that unity. In the DSM the scaling dimension of the QSH
operators is two, and we do not pick up any a signal: as apparent from Fig.
13, the correlation ratio decays as a function of system size. However at the
DSM-QSH transition it is worth noting a distinct cusp in the
$S^{\text{QSH}_{Z}}$ correlation ratio especially at $\Delta=0.75$ (see Fig.
13(c)). In the ordered XY-QSH phase, the anisotropy opens a gap in the Z-QSH
spectrum, leading to a reduction in the Z-QSH correlation ratio as a function
of system size (see Fig. 13). In the SSC phase the spin degrees of freedom are
gapped, such that the Z-QSH correlations decay exponentially. Again, in this
phase, the Z-QSH correlation ratio decays as a function of system size (see
Fig. 13).

It hence comes as a surprise that at the QSH-SSC transition, we see a distinct
cusp in the Z-QSH correlation ratio. This suggests that at this transition

$S^{\text{QSH}}(\bm{Q})\propto L^{1-\eta_{XY}}$ (21)

and

$S^{\text{QSH}_{Z}}(\bm{Q})\propto L^{1-\eta_{Z}}$ (22)

albeit with different scaling dimensions. Here, $\bm{Q}\equiv(0,0)$. Assuming
the same criticality as in Ref. Qin  _et al._ , 2017, $\eta_{XY}\approx 0.13$
and $\eta_{Z}\approx 0.91$. Given the large anomalous dimension, $\eta_{Z}$,
it is more advantageous to consider the susceptibility as defined in Eq. 36 of
Appendix B. For Lorentz invariant systems, this quantity scales as
$L^{2-\eta_{Z}}$ and suffers less from background effects. Fig. 19 plots the
correlation ratio as obtained from the susceptibility. While we can observe
clear cusps at the QSH-SSC transition we cannot unambiguously claim that this
quantity scales to a finite value in the thermodynamic limit.

Figure 14:  QSH (a), SSC (b), and QSHz (c) dynamical spectrum at the QSH-SSC
critical point ($\lambda=0.038$) for $\Delta=0.75$. Here, $L=18$. Blue lines
in (a) and (b) are the momentum dependence of the extrapolated excitation gap
of $\hat{J}^{XY}$ and $\hat{\eta}^{+}(\hat{\eta}^{-})$ operators. The blue
line in (c) is copied from the SSC dispersion relation in (b) to guide the
eye.

In Fig. 14 we show the spin-orbit, $\hat{J}^{XY}$, and the pairing,
$\hat{\eta}^{+}$($\hat{\eta}^{-}$), dynamical correlation functions at the
QSH-SSC critical point at $\Delta=0.75$. For the definition of the spectral
functions as well as for further data, we refer the reader to Appendix C. The
spectral functions display gapless excitations with the very same velocity.
This stands in accord with emergent Lorentz invariance. The $\hat{J}^{Z}$
spectral function, Fig. 14 (c), allows for an interpretation in terms of
gapless excitations at the $\Gamma$ point albeit with small spectral weight at
low energies.

##  IV Discussions and outlook

Our model realizes an easy-plane quantum spin Hall insulating state that
emerges from spontaneous U(1) broken spin symmetry, with a unnormalized three
component order parameter $\bm{N}$ defined in space-time. The norm of the
order parameter defines the single-particle gap. By tuning a single parameter
$\lambda$ at a given anisotropy strength, we observe DSM-QSH as well as QSH-
SSC transitions. At the DSM-QSH transition, the amplitude of the order
parameter vanishes and the single-particle gap closes. We understand this
transition in terms of Gross-Neveu-XY universality, the exponents of which
should be equivalent to those computed in Ref. Otsuka  _et al._ , 2018.

The focus of the paper, is on the QSH-SSC transition. Here, the QMC results
show that i) the single particle gap does not close at the transition and ii)
in contrast to mean-field calculations, our results on lattice sizes up to
$L=18$ support a continuous and direct transition for all considered values of
the anisotropy. We can safely omit the scenario of fine tuning, where
accidentally, the two transitions occur at the same value of $\lambda$. This
statement is based on the numerical observation that we consistently see a
direct transition for all values of the anisotropy. Furthermore, the
energetics are affected by the value of the imaginary time step we use in our
simulations. Were we at a fine tuning point, then we would have noticed
substantial changes in our results when varying the imaginary time step.

Since the single particle gap remains finite we can normalize the order
parameter vector $\bm{n}=\bm{N}/|\bm{N}|$, and attempt to understand the
transition in terms of fluctuations of $\bm{n}$. In this context, and as
discussed in the introduction, skyrmions or pairs of merons carry charge 2e
and condense at the transition. This suggests that our QSH-SSC transitions
flow to the easy-plane deconfined quantum critical point irrespective of the
anisotropy parameter. This is also supported numerically by the observed
universal value of the scaled correlation length upon changing the anisotropy.
The scaled correlation length takes the same value for the QSH and SSC
fluctuations at the easy-plane DQCP which supports an emergent O(4) symmetry
at least on intermediate length scales. Our model summarizes the very first
calculations of this critical point in a lattice Hamiltonian where the lattice
regularization does not break the IR symmetries of the putative field theory.
As a consequence, our lattice regularization does not introduce quadruple
monopoles as in the easy-plane JQ model. We refer the reader to Appendix E for
a detailed discussion of this point.

One of the characteristics of the easy-plane DQCP, is the emergence of
deconfined spinons at the critical point. In particular, one can adopt a CP1
representation of the order parameter

$\bm{n}\equiv z^{\dagger}\bm{\sigma}z\quad\quad\quad z^{\dagger}z=1$ (23)

(see Appendix E). The claim of DQC Senthil  _et al._ (2004b) is that the CP1
theory supports a deconfined phase at the critical point. The questions then
becomes how we can provide numerical evidence for this. Fractionalized spinons
are not directly measurable at the DQCP point since they do not directly
correspond to any local second quantized operators. However, the existence of
deconfined spinons suggests that at criticality the Z-component of the
correlation functions of the QSH order parameter shows power-law decay since
it just corresponds to spinon correlation functions. Our numerics support this
point of view.

To place our results in a broader perspective, we can ask the question of
whether designer Hamiltonians with higher symmetry can impact criticality.
Although it is well known that phase transitions numerically observed in easy-
plane lattice spin models have a higher tendency to be discontinuous in the
anisotropic case, two possible underlying physical interpretations remain
possible. First, the relevance of the $Z_{4}$ symmetry-breaking perturbation
at the easy-plane deconfined fixed point generally leads to a runaway flow,
explaining the first order nature of transition. Second, even if symmetry-
breaking terms due to lattice regularization (e.g., $O(4)$ down to $U(1)\times
U(1)$ or $U(1)\times Z_{4}$) are imposed to be zero, the easy-plane DQCP may
not even exist in any unitary conformal field theory. In this context our
results may be understood in terms of proximity to a DQCP that is not
accessible to our simulation space. This could correspond to a DQCP in the
complex plane Wang  _et al._ (2017) or in dimensions close to d=2 Nahum
(2020); Ma and Wang (2020).

###  ACKNOWLEDGMENTS

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V.
for funding this project by providing computing time on the GCS Supercomputer
SUPERMUC-NG at Leibniz Supercomputing Centre. F.F.A. acknowledges the DFG for
funding via Würzburg-Dresden Cluster of Excellence on Complexity and Topology
in Quantum Matter ct.qmat (EXC 2147, Project ID 390858490) as well as the
SFB1170 on Topological and Correlated Electronics at Surfaces and Interfaces.
T.S. acknowledges funding from the Deutsche Forschungsgemeinschaft under Grant
No. SA 3986/1-1. Y.L. was supported by the China Postdoctoral Science
Foundation under Grants No. 2019M660432 and No. 2020T130046 as well as the
National Natural Science Foundation of China under Grants No. 11947232 and
No.U1930402. D.H. and W.G. were supported by the National Natural Science
Foundation of China under Grants No. 12175015 and No. 11734002.

## References

  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Fradkin (2013a) E. Fradkin,  _Field Theories of Condensed Matter Physics_, 2nd ed. (Cambridge University Press, 2013).
  * Neto  _et al._ (2009) A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Herbut  _et al._ (2009) I. F. Herbut, V. Juričić, and B. Roy, Phys. Rev. B 79, 085116 (2009).
  * Grover and Senthil (2008) T. Grover and T. Senthil, Phys. Rev. Lett. 100, 156804 (2008).
  * Senthil  _et al._ (2004a) T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. Fisher, Phys. Rev. B 70, 144407 (2004a).
  * Senthil  _et al._ (2004b) T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004b).
  * Wang  _et al._ (2017) C. Wang, A. Nahum, M. A. Metlitski, C. Xu, and T. Senthil, Phys. Rev. X 7, 031051 (2017).
  * Sandvik (2007) A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
  * Shao  _et al._ (2016) H. Shao, W. Guo, and A. W. Sandvik, Science 352, 213 (2016).
  * Kuklov  _et al._ (2006) A. Kuklov, N. Prokof’ev, B. Svistunov, and M. Troyer, Annals of Physics 321, 1602 (2006), july 2006 Special Issue.
  * D’Emidio and Kaul (2016) J. D’Emidio and R. K. Kaul, Phys. Rev. B 93, 054406 (2016).
  * Sandvik  _et al._ (2002) A. W. Sandvik, S. Daul, R. R. P. Singh, and D. J. Scalapino, Phys. Rev. Lett. 89, 247201 (2002).
  * Kragset  _et al._ (2006) S. Kragset, E. Smørgrav, J. Hove, F. S. Nogueira, and A. Sudbø, Phys. Rev. Lett. 97, 247201 (2006).
  * Sen  _et al._ (2007) A. Sen, K. Damle, and T. Senthil, Phys. Rev. B 76, 235107 (2007).
  * D’Emidio and Kaul (2017) J. D’Emidio and R. K. Kaul, Phys. Rev. Lett. 118, 187202 (2017).
  * Desai and Kaul (2020) N. Desai and R. K. Kaul, Phys. Rev. B 102, 195135 (2020).
  * Liu  _et al._ (2019) Y. Liu, Z. Wang, T. Sato, M. Hohenadler, C. Wang, W. Guo, and F. F. Assaad, Nature Communications 10, 2658 (2019).
  * Wang  _et al._ (2021) Z. Wang, Y. Liu, T. Sato, M. Hohenadler, C. Wang, W. Guo, and F. F. Assaad, Phys. Rev. Lett. 126, 205701 (2021).
  * Domany  _et al._ (1984) E. Domany, M. Schick, and R. H. Swendsen, Phys. Rev. Lett. 52, 1535 (1984).
  * Blöte  _et al._ (2002) H. W. J. Blöte, W. Guo, and H. J. Hilhorst, Phys. Rev. Lett. 88, 047203 (2002).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * Gross and Neveu (1974) D. J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974).
  * Assaad  _et al._ (2022) F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. P. Toldin, J. S. E. Portela, and J. Schwab, SciPost Phys. Codebases , 1 (2022).
  * Blankenbecler  _et al._ (1981) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * White  _et al._ (1989) S. White, D. Scalapino, R. Sugar, E. Loh, J. Gubernatis, and R. Scalettar, Phys. Rev. B 40, 506 (1989).
  * Assaad and Evertz (2008) F. Assaad and H. Evertz, in  _Computational Many-Particle Physics_, Lecture Notes in Physics, Vol. 739, edited by H. Fehske, R. Schneider, and A. Weiße (Springer, Berlin Heidelberg, 2008) pp. 277–356.
  * Wu and Zhang (2005) C. Wu and S.-C. Zhang, Phys. Rev. B 71, 155115 (2005).
  * Sugiyama and Koonin (1986) G. Sugiyama and S. Koonin, Annals of Physics 168, 1 (1986).
  * Sorella  _et al._ (1989) S. Sorella, S. Baroni, R. Car, and M. Parrinello, EPL (Europhysics Letters) 8, 663 (1989).
  * Assaad  _et al._ (2013) F. F. Assaad, M. Bercx, and M. Hohenadler, Phys. Rev. X 3, 011015 (2013).
  * Assaad and Herbut (2013) F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
  * Liu  _et al._ (2021) Y. Liu, Z. Wang, T. Sato, W. Guo, and F. F. Assaad, Phys. Rev. B 104, 035107 (2021).
  * Otsuka  _et al._ (2018) Y. Otsuka, K. Seki, S. Sorella, and S. Yunoki, Phys. Rev. B 98, 035126 (2018).
  * Parisen Toldin  _et al._ (2015) F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Phys. Rev. B 91, 165108 (2015).
  * Kamieniarz and Blote (1993) G. Kamieniarz and H. W. J. Blote, Journal of Physics A: Mathematical and General 26, 201 (1993).
  * Qian  _et al._ (2016) X. Qian, Y. Deng, Y. Liu, W. Guo, and H. W. J. Blöte, Phys. Rev. E 94, 052103 (2016).
  * Sato  _et al._ (2017) T. Sato, M. Hohenadler, and F. F. Assaad, Phys. Rev. Lett. 119, 197203 (2017).
  * Iino  _et al._ (2019) S. Iino, S. Morita, N. Kawashima, and A. W. Sandvik, Journal of the Physical Society of Japan 88, 034006 (2019), https://doi.org/10.7566/JPSJ.88.034006 .
  * Qin  _et al._ (2017) Y. Q. Qin, Y.-Y. He, Y.-Z. You, Z.-Y. Lu, A. Sen, A. W. Sandvik, C. Xu, and Z. Y. Meng, Phys. Rev. X 7, 031052 (2017).
  * Nahum (2020) A. Nahum, Phys. Rev. B 102, 201116 (2020).
  * Ma and Wang (2020) R. Ma and C. Wang, Phys. Rev. B 102, 020407 (2020).
  * Francesco  _et al._ (1997) P. D. Francesco, P. Mathieu, and D. Senechal,  _Conformal Field Theory_ (Springer Verlag, Berlin, Heidelberg, New York, 1997) p. Good for Lie Algebras. At least very short and compact.
  * Beach (2004) K. S. D. Beach, (2004), cond-mat/0403055 .
  * Qi and Zhang (2008) X.-L. Qi and S.-C. Zhang, Phys. Rev. Lett. 101, 086802 (2008).
  * Fradkin (2013b) E. Fradkin,  _Field Theories of Condensed Matter Physics_, 2nd ed. (Cambridge University Press, 2013).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 1029 (1988).

##  Appendix A Mean field calculation

In this appendix, we present our mean-field calculation. Expanding Eq. 5 of
the main text as

$\displaystyle H_{\lambda}=$
$\displaystyle-\lambda\sum_{\varhexagon}\left[\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{x}_{\bm{i},\bm{j}}\right)^{2}+\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{y}_{\bm{i},\bm{j}}\right)^{2}\right.$
(24)
$\displaystyle\phantom{=\;\;}\left.+\Delta\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{z}_{\bm{i},\bm{j}}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle-\lambda\sum_{\varhexagon}\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle}\sum_{\langle\langle\bm{i^{\prime}}\bm{j^{\prime}}\rangle\rangle\neq\langle\langle\bm{i}\bm{j}\rangle\rangle}\bm{\hat{J}^{x}_{\langle\langle
i,j\rangle\rangle}\cdot\hat{J}^{x}_{\langle\langle
i^{\prime},j^{\prime}\rangle\rangle}}$
$\displaystyle+\bm{\hat{J}^{y}_{\langle\langle
i,j\rangle\rangle}\cdot\hat{J}^{y}_{\langle\langle
i^{\prime},j^{\prime}\rangle\rangle}}+\Delta\bm{\hat{J}^{z}_{\langle\langle
i,j\rangle\rangle}\cdot\hat{J}^{z}_{\langle\langle
i^{\prime},j^{\prime}\rangle\rangle}}$
$\displaystyle-\lambda\sum_{\varhexagon}\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle}[+(4+2\Delta)\hat{\eta}^{\dagger}_{\bm{i}}\hat{\eta}_{\bm{j}}+h.c+...]$

where

$\displaystyle\bm{\hat{J}}^{\alpha}_{\langle\langle i,j\rangle\rangle}\equiv
i\nu_{\bm{i}\bm{j}}\hat{c}^{\dagger}_{\bm{i}}\bm{\sigma}^{\alpha}\hat{c}_{\bm{j}}+H.c.,$
(25)
$\displaystyle\hat{\eta}_{\bm{i}}\equiv\hat{c}_{\bm{i}\downarrow}\hat{c}_{\bm{i}\uparrow},\
\ \ \
\hat{\eta}^{\dagger}_{\bm{i}}\equiv\hat{c}^{\dagger}_{\bm{i}\uparrow}\hat{c}^{\dagger}_{\bm{i}\downarrow}.$

The ellipsis denotes terms that do not contribute to the SSC or QSH ordering
within the mean-field decomposition.

The mean-field calculation involves selecting a polarization direction for the
two components of the QSH and SSC order parameters. The calculation is done by
numerically minimizing the free energy in the space of the two order
parameters.

The two order parameters as a function of $\lambda$ for the half-filled case
are shown in Fig. 15. For all three different values of anisotropy, we observe
a Dirac semi-metal ($\phi_{QSH}=\phi_{SSC}=0$), a pure QSH state
($\phi_{QSH}\neq 0,\phi_{SSC}=0$) as well as coexistence of QSH and SSC phases
($\phi_{QSH}\neq 0,\phi_{SSC}\neq 0$).

The mean-field phase diagram in Fig. 2 shows greater stability of the QSH
phase at stronger anisotropy. The reason for this becomes transparent when
taking a glimpse at Eq. 24. Here $\Delta$ modulates the magnitude of the pair-
hopping ($\hat{\eta}^{\dagger}_{\bm{i}}\hat{\eta}_{\bm{j}}+h.c.$.) but not of
the in-plane spin-orbit interactions ( $\hat{J}^{x}_{\langle\langle
i,j\rangle\rangle}\cdot\hat{J}^{x}_{\langle\langle
i^{\prime},j^{\prime}\rangle\rangle}+x\leftrightarrow y$ ).

Figure 15:  Mean-field solution for the QSH and SSC order parameters as a
function of $\lambda$ for $\Delta=0.1$ (a), $\Delta=0.5$ (b), and
$\Delta=0.75$ (c), respectively.

##  Appendix B Time displaced observables

To define susceptibilities in the realm of the zero temperature projective QMC
algorithm used in this paper, we distinguish between observables
$\hat{O}_{\bm{q}}=\frac{1}{\sqrt{N}}\sum_{\bm{r}}e^{i\bm{q}\cdot\bm{r}}\hat{O}_{\bm{r}}$
that commute or not with the Hamiltonian. Here, $N=L^{2}$. The key point is
that for $\left[\hat{O}_{\bm{q}},\hat{H}\right]=0$,

$\lim_{L\rightarrow\infty}\lim_{\beta\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)\neq\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)$
(26)

where

$\chi^{FT}(\bm{q},\beta,L)=\int_{0}^{\beta}d\tau\left(\langle\hat{O}_{\bm{q}}(\tau)\hat{O}_{-\bm{q}}\rangle_{T}-\langle\hat{O}_{\bm{q}}\rangle_{T}\langle\hat{O}_{-\bm{q}}\rangle_{T}\right).$
(27)

In the above,
$\langle\bullet\rangle_{T}=\frac{1}{Z}\text{Tr}\left(e^{-\beta\hat{H}}\bullet\right)$
and $\hat{O}_{\bm{q}}(\tau)=e^{\tau\hat{H}}\hat{O}e^{-\tau\hat{H}}$. In
particular, if the ground state is not degenerate then for
$\left[\hat{O}_{\bm{q}},\hat{H}\right]=0$,
$\lim_{L\rightarrow\infty}\lim_{\beta\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)=0$,
while
$\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)$
does not necessarily vanish. An example is the spin susceptibility for a
tight-binding model when the boundary conditions are chosen to ensure that the
ground state is non-degenerate.

Let us now consider the case of finite momentum $\bm{q}\neq 0$ for a momentum-
conserving Hamiltonian. Hence, $\left[\hat{O}_{\bm{q}},\hat{H}\right]\neq 0$.
Provided that the ground state is unique, we will show that the limits can be
interchanged. Our starting point is the Lehmann representation:

$\displaystyle\chi^{FT}(\bm{q},\beta,L)$
$\displaystyle=\frac{\beta}{Z}\sum_{n}e^{-\beta E_{n}}\left|\langle
n|\hat{O}_{\bm{q}}|n\rangle\right|^{2}+$ $\displaystyle\frac{1}{Z}\sum_{n\neq
m}\frac{e^{-\beta E_{m}}-e^{-\beta E_{n}}}{E_{m}-E_{n}}\left|\langle
n|\hat{O}_{\bm{q}}|m\rangle\right|^{2}$

with $\hat{H}|n\rangle=E_{n}|n\rangle$ and $n\in\mathbb{N}$. Since we have
assumed that $\bm{q}\neq 0$, the first term of the right-hand side of the
above equation vanishes. Defining the density of state as,
$N(E)=\lim_{L\rightarrow\infty}\sum_{n}\delta(E_{n}-E)$ we obtain:

$\displaystyle\lim_{L\rightarrow\infty}\lim_{\beta\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)$
$\displaystyle=\int dEN(E){\cal P}\frac{1}{E-E_{0}}\times$
$\displaystyle\left(\left|\langle
E|\hat{O}_{\bm{q}}|E_{0}\rangle\right|^{2}+\left|\langle
E_{0}|\hat{O}_{\bm{q}}|E\rangle\right|^{2}\right)$

and

$\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)=N(E_{0})\lim_{L\rightarrow\infty}\lim_{\beta\rightarrow\infty}\chi^{FT}(\bm{q},\beta,L)$
(30)

For a unique ground state, $N(E_{0})=1$. Hence, under the aforementioned
assumptions we can interchange the limits and it makes sense to define
susceptibilities within the ground state algorithm, where we first take the
limit of zero temperature and then consider larger and larger lattices.

For practical purposes, we compute:

$\chi(\bm{q})=\int_{0}^{\beta}\text{d}\tau\langle\hat{O}_{\bm{q}}(\tau)\hat{O}_{-\bm{q}}(0)\rangle\\\
-\langle\hat{O}_{\bm{q}}(\tau)\rangle\langle\hat{O}_{-\bm{q}}(0)\rangle$ (31)

where

$\langle\hat{O}_{\bm{q}}(\tau)\hat{O}_{-\bm{q}}(0)\rangle\equiv\frac{\langle\Psi_{T}|e^{-\theta\hat{H}}e^{-(\beta-\tau)\hat{H}}\hat{O}_{\bm{q}}e^{-\tau\hat{H}}\hat{O}_{-\bm{q}}e^{-\theta\hat{H}}|\Psi_{T}\rangle}{\langle\Psi_{T}|e^{-(2\theta+\beta)\hat{H}}|\Psi_{T}\rangle}$
(32)

and $|\Psi_{T}\rangle$ is the trial wave function. We consider $\beta=L$ for
all three values of $\Delta$. As shown in the main text, we implicitly checked
that for the considered size L, $\theta$ is chosen to be large enough to
converge to the ground state.

Since the zero temperature approach to susceptibilities matches the result
obtained with the traditional calculations, the scaling behaviors are
identical. In particular in the vicinity of a Lorentz invariant (z=1) critical
point we expect:

$\chi(\bm{q})\propto\xi^{-2\Delta+d+1}$ (33)

where the relationship between the scaling and anomalous dimensions reads
$2\Delta=d+z-2+\eta$ and $\xi$ is the diverging length scale in space and
time. Replacing the length scale with the linear size of our system yields the
desired result:

$\chi(\bm{q})\propto L^{2-\eta}.$ (34)

Hence, as for the finite temperature case, we expect at the critical point
that $\chi$ suppresses background contributions of the non-singular part of
free energy by an additional power (of the dynamical exponent). Liu  _et al._
(2019); Francesco  _et al._ (1997) Assuming Lorentz invariance at the easy-
plane DQCP, this additional power is unity.

We define the susceptibilities of the QSH, SSC, and QSHz in 35, 36 and 37.

$\displaystyle\chi^{\text{QSH}}_{m,n}(\bm{q})\equiv$
$\displaystyle\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}\int_{0}^{\beta}\text{d}\tau
e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}$ (35)
$\displaystyle\langle\hat{O}^{X}_{\bm{r},m}(\tau)\hat{O}^{X}_{\bm{r^{\prime}},n}(0)+\hat{O}^{Y}_{\bm{r},m}(\tau)\hat{O}^{Y}_{\bm{r^{\prime}},n}(0)\rangle.$

$\displaystyle\chi^{\text{QSH}_{Z}}_{m,n}(\bm{q})\equiv\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}\int_{0}^{\beta}\text{d}\tau
e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}\langle\hat{O}^{Z}_{\bm{r},m}(\tau)\hat{O}^{Z}_{\bm{r^{\prime}},n}(0)\rangle.$
(36)

$\displaystyle\chi^{\text{SSC}}_{a,b}(\bm{q})\equiv$
$\displaystyle\frac{1}{L^{2}}\sum_{\bm{r},\bm{r^{\prime}}}\int_{0}^{\beta}\text{d}\tau
e^{\mathrm{i}\bm{q}\cdot(\bm{r}-\bm{r}^{\prime})}$ (37)
$\displaystyle[\langle\hat{\eta}^{+}_{\bm{r},\bm{\widetilde{\delta}}_{a}}(\tau)\hat{\eta}^{-}_{\bm{r^{\prime}},\bm{\widetilde{\delta}}_{b}}(0)\rangle+\langle\hat{\eta}^{-}_{\bm{r},\bm{\widetilde{\delta}}_{a}}(\tau)\hat{\eta}^{+}_{\bm{r^{\prime}},\bm{\widetilde{\delta}}_{b}}(0)\rangle],$

After diagonalizing the corresponding three susceptibilities, we calculated
the renormalization-group invariant correlation ratio:

$R^{O}_{\chi}\equiv 1-\frac{\chi^{O}(\bm{Q}+\Delta\bm{q})}{\chi^{O}(\bm{Q})}$
(38)

with using the largest eigenvalue $\chi^{O}(\bm{q})$ with $O$ referring to
SSC, QSH, and QSHz, respectively. The ordering wave vector is $\bm{Q}=(0,0)$,
and $|\Delta\bm{q}|=\frac{4\pi}{\sqrt{3}L}$. We note that since the
superconducting order parameter breaks U(1) global charge symmetry and the QSH
order parameter breaks inversion symmetry so that the conditions for inter-
changing the limits of zero temperature and infinite size are satisfied.

Figure 16:  Time displaced correlation ratio $R^{SSC}_{\chi}$ (a) and
$R^{QSH}_{\chi}$ (b) as a function of $\lambda$ for $\Delta=0.1$.  Figure 17:
Same as Fig. 16 for $\Delta=0.5$.  Figure 18:  Same as Fig. 16 for
$\Delta=0.75$.  Figure 19:  Time displaced correlation ratio $R^{QSHz}_{\chi}$
for $\Delta=0.1$ (a), $\Delta=0.5$ (b), and $\Delta=0.75$ (c), respectively.

The finite size behavior of $R_{\chi}^{\text{QSH}}$ and
$R_{\chi}^{\text{SSC}}$ show no significant differences compared to the
corresponding equal time correlations considered in the main text. As shown in
Figs. 16, 17 and 18, the crossing points are consistent with the estimation of
critical points from the main text. On the other hand, the scaling of
$R_{\chi}^{\text{QSH}_{Z}}$ remains ambiguous: at the QSH-SSC transition
point, this quantity decays toward zero for strong anisotropy ($\Delta=0.1$).
In the case of $\Delta=0.5$ and $0.75$ it could converge to a finite constant
in thermodynamic limit; upon increasing the system sizes, its decreasing
tendency is similar to the one from equal time correlation ratio in the main
text.

##  Appendix C Spectrum

The spectral function for a given operator at zero temperature reads:

$A(\bm{k},\omega)=\pi\sum_{n}|\langle
n|\hat{O}|0\rangle|^{2}\delta(E_{n}-E_{0}-\omega).$ (39)

Here, $|0\rangle$ in Eq. 39 denotes the ground state, and $|n\rangle$ denotes
the eigenstates of the Hamiltonian with energy $E_{n}$. Given the imaginary-
time Monte Carlo data, the spectrum is obtained using the stochastic
analytical continuation approach Beach (2004) to solve for $A(\bm{k},\omega)$,
given $\langle 0|\hat{O}^{\dagger}(\tau)\hat{O}(0)|0\rangle$:

$\langle 0|\hat{O}^{\dagger}(\tau)\hat{O}(0)|0\rangle=\int d\omega
e^{-\tau\omega}A(\bm{k},\omega).$ (40)

Here, $\hat{O}$ represents the momentum space operators defined in the main
text: ${\hat{J}}^{XY}$, $\hat{\eta}$, $\hat{J}^{Z}$ and $\hat{c}$.
$A(\bm{k},\omega)$ is the corresponding spectral function, denoted as
$A^{\text{QSH}}$, $A^{\text{SSC}}$, $A^{\text{QSHz}}$ and $A^{sp}$,
respectively.

Figure 20:  QSH, SSC, QSHz, and single particle spectrum inside two (QSH, SSC)
phases and near the critical point, for $\Delta=0.1$. Blue lines are the
momentum dependence of the extrapolated excitation gap of $\hat{J}$ and
$\hat{\eta}^{+}(\hat{\eta}^{-})$ operators. We took $L=18$.  Figure 21:  Same
as Fig. 20 for $\Delta=0.5$.  Figure 22:  Same as Fig. 20 for $\Delta=0.75$.

As expected, the single-particle spectrum $A^{sp}$ is clearly gapped in both
the QSH and SSC ordered states, as well as across the QSH-SSC transition
points. As one can observe from Figs. 20, 21 and 22, $A^{sp}$ shows no
fundamental differences among three considered values of $\Delta$.

Deep inside the ordered QSH and SSC phases, and irrespective of the anisotropy
parameter, $\Delta$, the order parameter excitations exhibit Goldstone modes
stemming from global U(1) symmetry breaking. In particular, in Fig. 20, for
the case $\Delta=0.1$, a linear mode is observed for QSH operator at
$\lambda=0.04$ and for SSC operator at $\lambda=0.065$. The same behavior is
visible at $\Delta=0.5$ and $0.75$ in Fig. 21 and 22.

On the other hand, the excitation of both QSH and SSC order parameters at the
critical point shows a linear dispersion relation. Near the $\Gamma$ point,
the Goldstone mode is expected to give rise to a branch cut reflecting the
anomalous dimension:

$\displaystyle A_{\text{QSH}}(\bm{k},\omega)$
$\displaystyle\propto(v^{2}|\bm{k}|^{2}-\omega^{2})^{1-\frac{\eta_{\text{QSH}}}{2}}$
(41) $\displaystyle A_{\text{SSC}}(\bm{k},\omega)$
$\displaystyle\propto(v^{2}|\bm{k}|^{2}-\omega^{2})^{1-\frac{\eta_{\text{SSC}}}{2}}.$

Although the anomalous dimension of the two order parameters can in general be
different, the velocity $v$ is uniquely defined at a Lorentz invariant
critical point. We mark the velocity of these two excitations in Fig. 20, 21
and 22.

The spectrum of the $Z$ component of the spin current operator ($\hat{J}^{Z}$)
is controversial at the transition point. $A^{\text{QSH}_{Z}}(\bm{k},\omega)$
shows a clear gap around the $\Gamma$ point in the cases of $\Delta=0.1$ and
$0.5$. On the other hand, we observe that the gap decreases upon reducing the
anisotropy. As shown in Figs. 22 and 14 for $\Delta=0.75$, the value of gap at
the $\Gamma$ point is comparable to the finite size gap of the superconducting
fluctuations. A consistent picture in terms of easy-plane DQCP with deconfined
spinons at criticality requires that the excitation gap of $\hat{J}^{Z}$ to
scale to zero in the thermodynamic limit.

##  Appendix D Local detection of $Z_{2}$ topology using $\pi$ flux insertion

To detect the topology of our QSH insulator, we employ the magnetic flux
insertion approach Qi and Zhang (2008) that has successfully been tested in
Ref. Assaad  _et al._ , 2013. When $\pi$ fluxes are pumped locally into the
QSH insulator, mid-gap states carrying nontrivial spin quantum numbers are
exponentially localized around the flux insertion points. This approach
directly probes the $Z_{2}$ topological invariant, and we refer the reader to
Ref. Assaad  _et al._ , 2013 for a detailed discussion.

Consider the following kinetic Hamiltonian:

$\displaystyle\hat{H}_{t}=-t\sum_{\langle\bm{i},\bm{j}\rangle}(\hat{\bm{c}}^{\dagger}_{\bm{i}}\hat{\bm{c}}^{\phantom{\dagger}}_{\bm{j}}e^{iA_{\bm{i}\bm{j}}}+H.c.).$
(42)

and the interaction term:

$\displaystyle\hat{H}_{\lambda}=$
$\displaystyle-\lambda\sum_{\varhexagon}\left[\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{x}_{\bm{i},\bm{j}}e^{iA_{\bm{i}\bm{j}}}\right)^{2}+\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{y}_{\bm{i},\bm{j}}e^{iA_{\bm{i}\bm{j}}}\right)^{2}\right.$
(43)
$\displaystyle\phantom{=\;\;}\left.+\Delta\left(\sum_{\langle\langle\bm{i}\bm{j}\rangle\rangle\in\varhexagon}\hat{J}^{z}_{\bm{i},\bm{j}}e^{iA_{\bm{i}\bm{j}}}\right)^{2}\right]$

where $\bm{A}_{\bm{i}\bm{j}}$ is the vector potential that accounts for the
pair of $\pi$-fluxes. To practically insert a $\pi$-flux into our system, we
consider an arbitrary string connecting the centers of the two hexagons. Each
time an electron crosses this string, it acquires an $e^{i\pi}$ phase factor.

Figure 23:  $S^{xy}_{\Omega}(\bm{i})$ (a) and
$S^{\text{charge}}_{\Omega}(\bm{i})$ (b) distributions on a $L=12$ honeycomb
lattice. The simulation is performed deep inside the QSH state
($\Delta=0.1,\lambda=0.045$). Here, $\beta=24$.

Due to the easy-plane anisotropy, the dynamical generation of the QSH
insulator is associated with a long-range order of spin currents in the $U(1)$
plane. Thus, the mid-gap objects localized around the $\pi$ fluxes are Kramers
pairs of spin ‘up’ and ‘down’ states rotating in the $x-y$ plane, as well as
doublets of charge fluxons. The presence of localized spin and charge fluxons
can be captured by the low-energy spectral weight of
$c^{\dagger}_{\bm{i}}\sigma^{x}c^{\phantom{\dagger}}_{\bm{i}}$
($c^{\dagger}_{\bm{i}}\sigma^{y}c^{\phantom{\dagger}}_{\bm{i}}$) and
$c^{\dagger}_{\bm{i}}c^{\phantom{\dagger}}_{\bm{i}}$ operator:

$\displaystyle S^{xy}_{\Omega}(\bm{i})\equiv\int_{0}^{\Omega}d\omega
S^{xy}(\bm{i},\omega)$ (44) $\displaystyle
S^{\text{charge}}_{\Omega}(\bm{i})\equiv\int_{0}^{\Omega}d\omega
S^{\text{charge}}(\bm{i},\omega)\ \ \quad\Omega=0.25$

where

$\displaystyle S^{xy}(\bm{i},\omega)=\pi\sum_{n}|\langle
n|c^{\dagger}_{\bm{i}}\sigma^{x}c^{\phantom{\dagger}}_{\bm{i}}|0\rangle|^{2}\delta(\omega-
E_{n}-E_{0})$ (45) $\displaystyle
S^{\text{charge}}(\bm{i},\omega)=\pi\sum_{n}|\langle
n|c^{\dagger}_{\bm{i}}c^{\phantom{\dagger}}_{\bm{i}}|0\rangle|^{2}\delta(\omega-
E_{n}-E_{0})$

where $|n\rangle$ represents an energy eigenstate with energy $E_{n}$. The
energy window of $\Omega=0.25$ is well below twice single particle gap
($\Delta\approx 1.0$).

The enhanced spectral weight in $S^{xy}_{\Omega}(\bm{i})$ and
$S^{\text{charge}}_{\Omega}(\bm{i})$ around $\pi$ fluxes, as depicted in Fig.
23, clearly demonstrates the existence of spin and charge fluxons. This
numerically proves that the insulating state that we observe at intermediate
values of $\lambda$ is $Z_{2}$ nontrivial.

It is worth noting that our model does not exhibit quantized spin Hall
conductivity due to the absence of $U(1)$ spin conservation at low energies.
Therefore, the topological invariant characterizing our system is the $Z_{2}$
index, rather than the so-called spin Chern number.

##  Appendix E Absence of monopoles

The aim of this section is to show the absence of monopoles in our system. A
‘monopole’ corresponds to a singularity of the $U(1)$ gauge field in CP1
representation. It couples to the physical electric-magnetic vector potential
such that a $U(1)$ electric-magnetic gauge invariance is broken in the
presence of ‘monopoles’ Grover and Senthil (2008). Therefore, the statement of
‘monopole-free’ is a natural consequence of the exact $U(1)$ charge
conservation in our physical system.

A continuum description of our system in terms of fluctuating QSH order
parameter in $S^{2}$ space is:

$\displaystyle S=$ $\displaystyle\int
d^{2}\bm{x}d\tau\psi^{\dagger}[\gamma_{0}\gamma_{\mu}(\partial_{\mu}+ieA_{\mu})$
(46) $\displaystyle+$ $\displaystyle
im\bm{n}\cdot\gamma_{0}\gamma_{3}\gamma_{5}\vec{\sigma}]\psi+S_{\text{anisotropy}}+...$

where $\psi(\bm{x},\tau)$ is the Grassmann variable in the space of
$\mathds{R}^{2}\otimes\mathds{C}^{2}_{valley}\otimes\mathds{C}^{2}_{orbital}\otimes\mathds{C}^{2}_{spin}$.
$\bm{n}$ is the three-component order parameter of the QSH state and
$S_{\text{anisotropy}}$ describes the anisotropic term which breaks SU(2) spin
rotational symmetry down to U(1)$\times{\mathbb{Z}}_{2}$.
$A_{\mu}(\bm{x},\tau)$ is the physical electric-magnetic field. Note that our
lattice model doesn’t break particle number conservation, such that local
$U(1)$ gauge symmetry is satisfied:

$\displaystyle\psi(\bm{x},\tau)\rightarrow
e^{ie\theta(\bm{x},\tau)}\psi(\bm{x},\tau),$ (47) $\displaystyle
A_{\mu}(\bm{x},\tau)\rightarrow
A_{\mu}(\bm{x},\tau)-\partial_{\mu}\theta(\bm{x},\tau).$

Since $\bm{n}$ is normalized to unity, the single particle gap does not close,
and we can integrate out the fermions to obtain the effective action:

$S=S_{g}+S_{c}$ (48)

where

$S_{g}=\int
d^{2}\bm{x}d\tau\frac{1}{m}[(\partial_{\mu}\bm{n})\cdot(\partial_{\mu}\bm{n})]$
(49)

and the lowest-order electric-magnetic response is

$\displaystyle S_{c}=i\int d^{2}\bm{x}\int d\tau
A^{\mu}(\bm{x},\tau)J^{\mu}(\bm{x},\tau)$ (50)

where

$J^{\mu}(\bm{x},\tau)=2e\frac{1}{8\pi}\epsilon^{\mu\alpha\delta}\bm{n}\cdot(\partial_{\alpha}\bm{n}\times\partial_{\delta}\bm{n})$
(51)

Now reformulate the action into a gauge redundant (CP1) representation:

$\vec{N}\equiv z^{\dagger}\vec{\sigma}z\quad\quad\quad z^{\dagger}z=1$ (52)

such that $S_{g}$ is described by
$\frac{1}{g}|(\partial_{\mu}-ia_{\mu})z|^{2}$. A local $U(1)$ gauge redundancy
comes naturally from the $CP^{1}$ reformulation:

$z\longrightarrow e^{i\chi}z,\ \ \ \ a_{\mu}\longrightarrow
a_{\mu}+\partial_{\mu}\chi$ (53)

Crucially the following identity holds:

$\displaystyle\frac{1}{8\pi}\epsilon^{\mu\nu\lambda}\bm{n}\cdot(\partial_{\nu}\bm{n}\times\partial_{\lambda}\bm{n})=\frac{1}{4\pi}\epsilon^{\mu\nu\lambda}\partial_{\nu}a_{\lambda}.$
(54)

This follows from the the saddle point results Fradkin (2013b):

$a_{\mu}=iz^{\dagger}\partial_{\mu}z.$ (55)

To describe skyrmions one has to allow for compactness of $a_{\mu}$ :

$\oint\bm{a}\cdot d\bm{l}=n2\pi\,\,\text{ with }\,\,n\in\mathbb{Z}$ (56)

for a closed loop. Using Stokes theorem, integration at a given time slice
hence gives: $\int
d^{2}\bm{x}\epsilon^{0\mu\nu}(\partial_{\mu}a_{\nu}-\partial_{\nu}a_{\mu})=4\pi
n$ and counts the number of ‘skyrmions’ in the time slice. A magnetic
‘monopole’ corresponds to field configurations, $a_{\mu}$, that have a non-
vanishing ‘magnetic flux’, $\epsilon^{\delta\mu\nu}\partial_{\mu}a_{\nu}$,
through a closed surface in space-time.

With Eq. 54 and Eq. 51, the Chern-Simons part of action reads

$S_{c}=\int
d^{2}\bm{x}d\tau\frac{ie}{2\pi}\epsilon^{\mu\nu\kappa}A_{\mu}\partial_{\nu}a_{\kappa}$
(57)

To show the absence of monopoles, we consider a uniform gauge transformation
inside a sphere:

$A_{\mu}\longrightarrow A_{\mu}+\partial_{\mu}\Theta$ (58)

with

$\Theta(\tau,x,y)=\left\\{\begin{array}[]{ll}\Theta_{0}&\text{ for }r<R\\\
0&\text{ for }r\geq R\end{array}\right..$ (59)

In the above $r=\sqrt{\tau^{2}+x^{2}+y^{2}}$. As a consequence,
$\partial_{\mu}\Theta$ does not vanish only on sphere of radius R, $S_{R}$. In
particular:

$\displaystyle\delta S=$ $\displaystyle\frac{ie}{2\pi}\int
d^{3}\bm{x}\epsilon^{\mu\nu\kappa}\partial_{\mu}\Theta\,\partial_{\nu}\,a_{\kappa}$
(60) $\displaystyle=$
$\displaystyle\frac{ie}{2\pi}\Theta_{0}\int_{S_{R}}d\bm{s}_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}\,a_{\kappa}$
$\displaystyle\equiv$ $\displaystyle\frac{ie}{2\pi}\Theta_{0}Q.$

In the above, $d\bm{s}$ defines a surface element of $S_{R}$ and $Q$ the
number of monopoles within $S_{R}$. If $Q\neq 0$ then U(1) local gauge
invariance (see Eq. 47 ) is not satisfied. Hence, charge conservation, or
equivalently U(1) local gauge invariance, requires $Q=0$ and the absence of
monopoles.

It is important to emphasize that, based on spin coherent path integral
calculation by Haldane, Haldane (1988) the monopole configurations of anti-
ferromagnetically coupled spin $1/2$ system on the square lattice also carry a
nontrivial phase factor upon a $U(1)$ gauge transformation. However, in this
case, instead of the physical electromagnetic field, this $U(1)$
transformation corresponds to the lattice rotation. This symmetry is broken
upon lattice regularization, and only a $Z_{4}$ subgroup of it is conserved.
Under $\pi/2$ lattice rotation, a single monopole contributes by

$\delta S=i\frac{\pi}{2}$ (61)

to the action. Hence, quadruple monopoles configurations are allowed.


# An extension of an asymptotic result of Tricomi concerning a definite
integral

R. B. Paris
Division of Computing and Mathematics,
Abertay University, Dundee DD1 1HG, UK

###### Abstract

We consider the expansion of an integral considered by F.G. Tricomi given by

$\int_{-\infty}^{\infty}xe^{-x^{2}}(\mbox{${\textstyle\frac{1}{2}}$}+\mbox{${\textstyle\frac{1}{2}}$}\mbox{erf}\,x)^{m}dx$

as $m\to\infty$. The procedure involves a suitable change of variable and the
inversion of the complementary error function $\mbox{erfc}\,x$. Numerical
results are presented to demonstrate the accuracy of the expansion.

A second part examines an extension of an integral arising in airfoil theory.

Mathematics subject classification (2010): 33B20, 33E20, 34E05, 41A60

Keywords: Tricomi integral, asymptotic expansion, inversion, airfoil integral

1\. Introduction

In a short note published in the Accademia dei Lincei [4], F.G. Tricomi
obtained the leading asymptotic behaviour of the integral arising in
probability theory

$\int_{-\infty}^{\infty}xe^{-x^{2}}\biggl{(}\frac{1+\mbox{erf}\,x}{2}\biggr{)}^{\\!m}dx$

as $m\to\infty$, where $\mbox{erf}\,x$ is the error function defined by

$\mbox{erf}\,x=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt.$

The factor in brackets in the above integral satisfies
$0<\mbox{${\textstyle\frac{1}{2}}$}+\mbox{${\textstyle\frac{1}{2}}$}\mbox{erf}\,x<1$
for $x\in(-\infty,\infty)$, with the upper and lower bounds approached as
$x\to\pm\infty$, respectively. By a suitable transformation of the integration
variable, Tricomi showed that the leading behaviour of this integral is given
by

$\hskip 28.45274pt\frac{\sqrt{\pi\log\,m}}{m}\qquad(m\to\infty).$ (1.1)

In the first part of this paper we extend the result in (1.1) by obtaining
higher-order terms in the large-$m$ expansion of the integrals

$I_{n,m}=\int_{-\infty}^{\infty}x^{n}e^{-x^{2}}\biggl{(}\frac{1+\mbox{erf}\,x}{2}\biggr{)}^{\\!m}dx\qquad(n=0,1,2),$
(1.2)

the case $n=0$ being trivial. The second part considers the evaluation of an
extension of an integral arising in airfoil theory.

2\. Extension of Tricomi’s analysis

Following Tricomi [4], we make the change of variable

$e^{-t}=\frac{1+\mbox{erf}\,x}{2},\qquad\frac{e^{-x^{2}}}{\sqrt{\pi}}\,\frac{dx}{dt}=-e^{-t}$
(2.1)

to yield

$I_{n,m}=-\int_{0}^{\infty}x^{n}e^{-x^{2}}\frac{dx}{dt}\,e^{-mt}\,dt=\sqrt{\pi}\int_{0}^{\infty}x^{n}e^{-st}\,dt,\qquad
s:=m+1.$ (2.2)

When $n=0$, we then have the trivial result

$I_{0,m}=\frac{\sqrt{\pi}}{m+1}.$

To deal with the cases $n=1$ and $n=2$, we require the inversion of the first
relation in (2.1) written in the form

$\mbox{erfc}\,x=2-2e^{-t}$ (2.3)

to yield $x\equiv x(t)$, which is finite and continuous for $0<t<\infty$. From
the well-known asymptotic behaviour [2, p. 164]

$\mbox{erfc}\,x\sim\frac{e^{-x^{2}}}{\sqrt{\pi}x},\qquad\mbox{erfc}(-x)\sim
2-\frac{e^{-x^{2}}}{\sqrt{\pi}x}\qquad(x\to+\infty),$

it is seen that $x(t)\sim\sqrt{-\log\,t}$ as $t\to 0$ and $x(t)\sim-\sqrt{t}$
as $t\to+\infty$, with $x(t)=0$ when $t=\log\,2$.

To proceed we split the integral into two parts: $[0,\alpha]$ and
$(\alpha,\infty)$, where $\alpha$ is finite and bounded away from zero. Then,
with the change of variable $t=u/s$, we have

$I_{n,m}=\frac{\sqrt{\pi}}{s}\int_{0}^{\alpha
s}x^{n}e^{-u}du+\frac{\sqrt{\pi}}{s}\int_{\alpha s}^{\infty}x^{n}e^{-u}du$
(2.4)

where $x\equiv x(u/s)$. Since $x(t)>-K\sqrt{t}$ for some suitable positive
constant $K$, the second integral is bounded by

$\frac{1}{s}\int_{\alpha
s}^{\infty}|x(u/s)|^{n}e^{-u}du<\frac{K^{n}}{s^{n/2+1}}\int_{\alpha
s}^{\infty}u^{n/2}e^{-u}du=\frac{K^{n}}{s^{n/2+1}}\Gamma(\mbox{${\textstyle\frac{1}{2}}$}n+1,\alpha
s)=O(s^{-1}e^{-\alpha s}),$

as $s\to\infty$, where we have employed the result for the (upper) incomplete
gamma function $\Gamma(a,z)\sim z^{a-1}e^{-z}$ as $z\to+\infty$.

In the first integral in (2.4) the dominant contribution will arise from the
neighbourhood of $u=0$ as $m\to\infty$. We therefore require the inversion of
(2.3) as $t\to 0$ expressed in the form

$\mbox{erfc}\,x=\frac{2u}{s}\biggl{(}1-\frac{u}{2s}+\cdots\biggr{)}\qquad(u\to
0).$ (2.5)

In the appendix it is shown that (see (A.2))

$x^{2}=\log\,s\biggl{(}1+\frac{A(u)}{L}+\frac{B(u)}{2L^{2}}+\frac{C(u)}{4L^{3}}+O(L^{-4},(sL)^{-1})\biggr{)},$
(2.6)

where we have defined

$L=\log\,s,\quad L_{1}=\log\sqrt{\log\,s},\quad a=2\sqrt{\pi},$

and the first few coefficients are

$A(u)=-\log\,au-L_{1},\quad B(u)=-A(u)-1,\quad
C(u)=A^{3}(u)+3A(u)+\mbox{${\textstyle\frac{7}{2}}$}.$

We define the integrals

$\lambda_{n}:=\int_{0}^{\infty}(\log\,u)^{n}e^{-u}du,\qquad(n=0,1,2,\ldots),$

so that $\lambda_{0}=1$ and

$\lambda_{1}=-\gamma,\quad\lambda_{2}=\gamma^{2}+\frac{\pi^{2}}{6},\quad\lambda_{3}=-\gamma^{2}-\frac{\pi^{2}\gamma}{2}-2\zeta(3),$

where $\gamma=0.57721...$ is the Euler-Mascheroni constant and $\zeta$ is the
Riemann zeta function.

Then, in the case $n=2$, we have upon extending the upper limit of integration
to $\infty$ (thereby introducing an exponentially small error of O$(e^{-\alpha
s})$)

$\displaystyle I_{2,m}$ $\displaystyle\sim$
$\displaystyle\frac{\sqrt{\pi}\log\,s}{s}\int_{0}^{\infty}\biggl{(}1+\frac{A(u)}{L}+\frac{B(u)}{2L^{2}}+\frac{C(u)}{4L^{3}}+\cdots\biggr{)}e^{-u}du$
(2.7) $\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}\log\,s}{s}\biggl{\\{}1+\frac{{\cal
A}_{2}}{L}+\frac{{\cal B}_{2}}{2L^{2}}+\frac{{\cal
C}_{2}}{4L^{3}}+\cdots\biggr{\\}}\qquad(m\to\infty),$

where

${\cal A}_{2}=G-L_{1},\qquad{\cal B}_{2}=-G+L_{1}-1,$ ${\cal
C}_{2}=G^{2}+(3-2L_{1})G+L_{1}^{2}-3L_{1}+\frac{\pi^{2}}{6}+\frac{7}{2}$

with $G:=\gamma-\log\,a$.

In Tricomi’s case with $n=1$, we find from (2.6) that

$x=\sqrt{\log\,s}\biggl{(}1+\frac{A(u)}{2L}+\frac{2B(u)-A^{2}(u)}{8L^{2}}+\frac{2C(u)-2A(u)B(u)+A^{3}(u)}{16L^{3}}+O(L^{-4},(sL)^{-1})\biggr{)}.$

Hence we obtain the expansion

$I_{1,m}\sim\frac{\sqrt{\pi\log\,s}}{s}\biggl{\\{}1+\frac{{\cal
A}_{1}}{2L}-\frac{{\cal B}_{1}}{8L^{2}}+\frac{{\cal
C}_{1}}{16L^{3}}+\cdots\biggr{\\}}\qquad(m\to\infty),$ (2.8)

where

${\cal A}_{1}=G-L_{1},\qquad{\cal
B}_{1}=G^{2}+2(1-L_{1})G+L_{1}^{2}-2L_{1}+\frac{\pi^{2}}{6}+2,$ ${\cal
C}_{1}=G^{3}+(4-3L_{1})G^{2}+(3L_{1}^{2}-8L_{1}+8+\frac{\pi^{2}}{6})G+4L_{1}^{2}-8L_{1}-L_{1}^{3}$
$+\frac{\pi^{2}}{6}(4-3L_{1})+2\zeta(3)+7$

and we recall that $s=m+1$. It is seen that the leading term in this expansion
agrees with the result stated in (1.1).

As a numerical verification of the expansions (2.7) and (2.8) we present111In
the tables we write $x\times 10^{y}$ as $x(y)$. in Table 1 the asymptotic
values of $I_{n,m}$ ($n=1,2$) for different $m$ compared with the values
obtained by high-precision evaluation of (1.2). In Table 2 we show the
absolute relative error in the computation of $I_{1,m}$ for different $m$ as a
function of the truncation index $k$ in the expansion (2.8).

Table 1: The values of $I_{n,m}$ ($n=1,2$) compared with the expansions (2.7) and (2.8) for different $m$. $m$ | $I_{1,m}$ | (2.7) | $I_{2,m}$ | (2.6)
---|---|---|---|---
$10^{2}$ | $3.116097(-2)$ | $3.119672(-2)$ | $5.694564(-2)$ | $5.695077(-2)$
$10^{3}$ | $4.058838(-3)$ | $4.060226(-3)$ | $9.413132(-3)$ | $9.414250(-3)$
$10^{4}$ | $4.826833(-4)$ | $4.827422(-4)$ | $1.322800(-3)$ | $1.32285?(-3)$
$10^{5}$ | $5.494877(-5)$ | $5.495164(-5)$ | $1.710063(-4)$ | $1.710084(-4)$
$10^{6}$ | $6.094732(-6)$ | $6.094889(-6)$ | $2.101178(-5)$ | $2.101184(-5)$
Table 2: The absolute relative error in $I_{1,m}$ using the expansion (2.8) for different truncation index $k$. $k$ | $m=10^{4}$ | $m=10^{5}$ | $m=10^{6}$
---|---|---|---
0 | $1.143(-1)$ | $9.447(-2)$ | $8.094(-2)$
1 | $5.526(-3)$ | $3.686(-3)$ | $2.656(-3)$
2 | $1.361(-5)$ | $1.013(-4)$ | $7.440(-5)$
3 | $1.220(-5)$ | $5.214(-5)$ | $2.581(-5)$

3\. An extension of an integral arising in airfoil theory

In airfoil theory, the expression for the downward velocity $w(a)$ at a point
$x=a$ over the wingspan induced by a single vortex filament is [3, p. 201]

$w(a)=-\frac{1}{4\pi}\int_{-1}^{1}\frac{1}{x-a}\,\frac{d\Gamma}{dx}\,dx\qquad(-1<a<1),$

where $\Gamma$ is the circulation and the $x$-coordinate is normalised to the
wingspan so that $-1\leq x\leq 1$. For an elliptical distribution
$\Gamma=\Gamma_{0}(1-x^{2})^{1/2}$, where $\Gamma_{0}$ is the circulation in
the middle of the airfoil, we have

$w(a)=\frac{\Gamma_{0}}{4\pi}\int_{-1}^{1}\frac{x}{x-a}\,\frac{dx}{\sqrt{1-x^{2}}}=\frac{\Gamma_{0}}{4\pi}\biggl{\\{}\int_{-1}^{1}\frac{dx}{\sqrt{1-x^{2}}}+a\int_{-1}^{1}\frac{1}{x-a}\,\frac{dx}{\sqrt{1-x^{2}}}\biggr{\\}}$
$=\frac{\Gamma_{0}}{4\pi}\biggl{\\{}\pi+a\int_{-1}^{1}\frac{1}{x-a}\,\frac{dx}{\sqrt{1-x^{2}}}\biggr{\\}}.$
(3.1)

It is clear that $w(-a)=w(a)$ so that it is sufficient to consider $0<a<1$.

The integral appearing in (3.1) must be treated as a Cauchy principal value;
that is, we evaluate the integral up to $\epsilon$ either side of $x=a$ and
then let $\epsilon\to 0$. Thus we have, with
$R(w):=(1-a^{2}-2aw-w^{2})^{1/2}$,

$\int_{-1}^{1}\frac{1}{x-a}\,\frac{dx}{\sqrt{1-x^{2}}}=\int_{\epsilon}^{1-a}\frac{dw}{wR(w)}+\int_{1+a}^{\epsilon}\frac{dw}{wR(-w)}$
$=\frac{1}{\sqrt{1-a^{2}}}\biggl{\\{}\biggl{[}\log\,\frac{w}{2(1\\!-\\!a^{2}\\!-\\!aw\\!+\\!\sqrt{1-a^{2}}R(w))}\biggr{]}^{1-a}_{\epsilon}+\biggl{[}\log\,\frac{w}{2(1\\!-\\!a^{2}\\!+\\!aw\\!+\\!\sqrt{1-a^{2}}R(-w))}\biggr{]}^{\epsilon}_{1+a}\biggr{\\}}$
$=-\frac{2a\epsilon}{(1-a^{2})^{3/2}}+O(\epsilon^{2})$

after some routine algebra. Thus the integral vanishes as $\epsilon\to 0$ with
the result that

$\int_{-1}^{1}\frac{x}{x-a}\,\frac{dx}{\sqrt{1-x^{2}}}=\pi$ (3.2)

independent of $a$ [3, p. 202]. Hence $w(a)=\Gamma_{0}/4$, which means that in
the case of an elliptical distribution of the circulation the induced downward
velocity is constant across the wing.

The situation when $d\Gamma/dx=x^{2n+1}\Gamma_{0}/\sqrt{1-x^{2}}$,
$n=0,1,2,\ldots\ $, corresponds to a flattening of the basic elliptic profile
($n=0$). For example, when $n=1,2$ we have the profiles

$\Gamma=\frac{\Gamma_{0}}{2}(2+x^{2})\sqrt{1-x^{2}}\ \ \
(n=1),\qquad\Gamma=\frac{\Gamma_{0}}{15}(8+4x^{2}+3x^{4})\sqrt{1-x^{2}}\ \ \
(n=2);$

as $n$ increases the profile becomes progressively flatter in the central
portion of the wing. Then the extension of the integral (3.1) we consider is

$J_{n}(a;\mu):=\int_{-1}^{1}\frac{x^{2n+1}}{x-a}\,\frac{dx}{(1-x^{2})^{\mu}}\qquad(\mu<1,\
n=0,1,2,\ldots).$ (3.3)

which reduces to (3.1) when $n=0$ and $\mu=\mbox{${\textstyle\frac{1}{2}}$}$.
It is easily seen that $J_{n}(-a;\mu)=J_{n}(a;\mu)$, so that again it is
sufficient to consider $0<a<1$.

3.1 The evaluation of the integral $J_{n}(a;\mu)$

We write $J_{n}(a;\mu)$ in the form

$J_{n}(a;\mu)=\int_{-1}^{1}\frac{x^{2n+1}-a^{2n+1}}{x-a}\,\frac{dx}{(1-x^{2})^{\mu}}+a^{2n+1}\int_{-1}^{1}\frac{1}{x-a}\,\frac{dx}{(1-x^{2})^{\mu}}\equiv
I_{1}+a^{2n+1}I_{2}.$

In the first integral we employ the expansion

$\frac{x^{m+1}-a^{m+1}}{x-a}=\sum_{r=0}^{m}a^{m-r}x^{r}\qquad(m=0,1,2,\ldots)$
(3.4)

to find

$I_{1}=\sum_{r=0}^{n}a^{2n-2r}\int_{-1}^{1}\frac{x^{2r}}{(1-x^{2})^{\mu}}\,dx=\frac{a^{2n}\sqrt{\pi}\Gamma(1-\mu)}{\Gamma(\mbox{${\textstyle\frac{3}{2}}$}-\mu)}\sum_{r=0}^{n}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{r}a^{-2r}}{(\mbox{${\textstyle\frac{3}{2}}$}-\mu)_{r}}$
$=\frac{a^{2n}\sqrt{\pi}\Gamma(1-\mu)}{\Gamma(\mbox{${\textstyle\frac{3}{2}}$}-\mu)}\,{}_{2}F_{1}(\mbox{${\textstyle\frac{1}{2}}$},1;\mbox{${\textstyle\frac{3}{2}}$}-\mu;a^{-2})|_{n},$

where the subscript $n$ on the Gauss hypergeometric function ${}_{2}F_{1}$
denotes that only the first $n+1$ terms are to be taken.

The second integral is

$I_{2}=\sum_{k=0}^{\infty}\frac{(\mu)_{k}}{k!}\int_{-1}^{1}\frac{x^{2k}}{x-a}\,dx=\sum_{k=0}^{\infty}\frac{(\mu)_{k}}{k!}\biggl{\\{}a^{2k}\int_{-1}^{1}\frac{dx}{x-a}+\int_{-1}^{1}\frac{x^{2k}-a^{2k}}{x-a}\,dx\biggr{\\}}$
$=\frac{1}{(1-a^{2})^{\mu}}\log\biggl{(}\frac{1-a}{1+a}\biggr{)}+\sum_{k=1}^{\infty}\frac{(\mu)_{k}}{k!}\int_{-1}^{1}\frac{x^{2k}-a^{2k}}{x-a}\,dx.$

Use of (3.4) in the above integral shows that

$S:=\sum_{k=1}^{\infty}\frac{(\mu)_{k}}{k!}\int_{-1}^{1}\frac{x^{2k}-a^{2k}}{x-a}\,dx=\frac{1}{a}\sum_{k=1}^{\infty}\frac{(\mu)_{k}}{k!}\sum_{r=0}^{k-1}a^{2k-2r}\int_{-1}^{1}x^{2r}dx$
$=2a\sum_{k=0}^{\infty}\frac{(\mu)_{k+1}a^{2k}}{(k+1)!}\sum_{r=0}^{k}\frac{a^{-2r}}{2r+1}=2a\sum_{r=0}^{\infty}\frac{a^{-2r}}{2r+1}\sum_{k=r}^{\infty}\frac{(\mu)_{k+1}a^{2k}}{(k+1)!}$
$=2a\sum_{r=0}^{\infty}\frac{1}{2r+1}\sum_{k=0}^{\infty}\frac{(\mu)_{k+r+1}a^{2k}}{(k+r+1)!},$

where we have put $k\to k+r$ in the last sum. From the identity
$(\mu)_{k+r+1}=(\mu)_{r+1}(\mu+r+1)_{k}$, we then obtain

$\displaystyle S$ $\displaystyle=$ $\displaystyle
2a\sum_{r=0}^{\infty}\frac{(\mu)_{r+1}}{(2r+1)(r+1)!}\sum_{k=0}^{\infty}\frac{(\mu+r+1)_{k}a^{2k}}{(r+2)_{k}}$
$\displaystyle=$ $\displaystyle
2a\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)r!}\,{}_{2}F_{1}(1,\mu+r;r+1;a^{2})$
$\displaystyle=$
$\displaystyle\frac{2a}{1-a^{2}}\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)r!}\,{}_{2}F_{1}\biggl{(}1,1-\mu;r+1;-\frac{a^{2}}{1-a^{2}}\biggr{)}$

upon use of the well-known Euler transformation [2, (15.8.1)] for the
hypergeometric function.

Then we have the result

$J_{n}(a;\mu)=\frac{a^{2n}\sqrt{\pi}\Gamma(1-\mu)}{\Gamma(\mbox{${\textstyle\frac{3}{2}}$}-\mu)}\,{}_{2}F_{1}(\mbox{${\textstyle\frac{1}{2}}$},1;\mbox{${\textstyle\frac{3}{2}}$}-\mu;a^{-2})|_{n}+\frac{a^{2n+1}}{(1-a^{2})^{\mu}}\,\log\biggl{(}\frac{1-a}{1+a}\biggr{)}$
$+\frac{2a^{2n+2}}{1-a^{2}}\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)r!}\,{}_{2}F_{1}\biggl{(}1,1-\mu;r+1;-\frac{a^{2}}{1-a^{2}}\biggr{)}$
(3.5)

valid for $\mu<1$ and non-negative integer $n$. We remark that the convergence
of the infinite sum in (3.5) is slow, being controlled by $r^{\mu-2}$ as
$r\to\infty$. A simple means of accelerating the convergence is given in the
sub-section below.

3.2 Improved convergence of (3.5)

We can write the hypergeometric function appearing in the infinite sum in
(3.5) in the form

${}_{2}F_{1}(1,1-\mu;r+1;-X)=1-\frac{(1-\mu)X}{r+1}+\frac{(1-\mu)_{2}X^{2}}{(r+1)_{2}}{}_{2}F_{1}(1,3-\mu;r+3;-X),$

where, for convenience, we have put $X:=a^{2}/(1-a^{2})$. Then

$S=\frac{2a}{1-a^{2}}\biggl{(}\sigma_{0}-(1-\mu)X\sigma_{1}+(1-\mu)_{2}X^{2}\sigma_{2}\biggr{)}$
$+\frac{2a}{1-a^{2}}(1-\mu)_{2}X^{2}\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)(r+2)!}\,\biggl{\\{}{}_{2}F_{1}(1,3-\mu;r+3;-X)-1\biggr{\\}},$
(3.6)

where

$\displaystyle\sigma_{0}$ $\displaystyle=$
$\displaystyle\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)r!}=1-\frac{\sqrt{\pi}\Gamma(1-\mu)}{\Gamma(\mbox{${\textstyle\frac{1}{2}}$}-\mu)},$
$\displaystyle\sigma_{1}$ $\displaystyle=$
$\displaystyle\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)(r+1)!}=\frac{2-3\mu}{3(1-\mu)}-\frac{2\sqrt{\pi}\Gamma(1-\mu)}{3\Gamma(\mbox{${\textstyle\frac{1}{2}}$}-\mu)},$
$\displaystyle\sigma_{2}$ $\displaystyle=$
$\displaystyle\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)(r+2)!}=\frac{16-25\mu+15\mu^{2}}{30(1-\mu)(2-\mu)}-\frac{4\sqrt{\pi}\Gamma(1-\mu)}{15\Gamma(\mbox{${\textstyle\frac{1}{2}}$}-\mu)}.$

The convergence of the final sum in (3.6) is now controlled by terms of
O$(r^{\mu-5})$ as $r\to\infty$, which is an improvement on that in (3.5). This
procedure can be continued to produce higher rates of convergence, but at the
expense of the evaluation of higher-order sums $\sigma_{m}$, $m\geq 2$.

3.3 The special case $n=0$, $\mu=\mbox{${\textstyle\frac{1}{2}}$}$

We demonstrate that the result in (3.5) reduces to the value given in (3.2)
when $n=0$ and $\mu=\mbox{${\textstyle\frac{1}{2}}$}$. We have from (3.5) with
$X=a^{2}/(1-a^{2})$

$J_{0}(a;\mbox{${\textstyle\frac{1}{2}}$})=\pi+\frac{a}{\sqrt{1-a^{2}}}\,\log\biggl{(}\frac{1-a}{1+a}\biggr{)}+2X\sum_{r=1}^{\infty}\frac{(\mu)_{r}}{(2r-1)r!}\,{}_{2}F_{1}(\mbox{${\textstyle\frac{1}{2}}$},1;r+1;-X).$

The sum can be written in the form

$T:=X\sum_{r=0}^{\infty}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{r+1}}{(r+\mbox{${\textstyle\frac{1}{2}}$})(r+1)!}\sum_{k=0}^{\infty}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{k}(-X)^{k}}{(r+2)_{k}}=X\sum_{k=0}^{\infty}(\mbox{${\textstyle\frac{1}{2}}$})_{k}(-X)^{k}\sum_{r=0}^{\infty}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{r+1}}{(r+\mbox{${\textstyle\frac{1}{2}}$})(r+k+1)!}$
$=X\sum_{k=0}^{\infty}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{k}(-X)^{k}}{(k+1)!}\,{}_{2}F_{1}(\mbox{${\textstyle\frac{1}{2}}$},1;k+2;1).$

Application of Gauss’ summation theorem series [2, (15.4.20)] shows that the
${}_{2}F_{1}(1)$ series has the value
$(k+1)\Gamma(k+\mbox{${\textstyle\frac{1}{2}}$})/\Gamma(k+\mbox{${\textstyle\frac{3}{2}}$})$,
whence

$T=2X\sum_{k=0}^{\infty}\frac{(\mbox{${\textstyle\frac{1}{2}}$})_{k}(\mbox{${\textstyle\frac{1}{2}}$})_{k}}{(\mbox{${\textstyle\frac{3}{2}}$})_{k}k!}\,(-X)^{k}=2X\,{}_{2}F_{1}(\mbox{${\textstyle\frac{1}{2}}$},\mbox{${\textstyle\frac{1}{2}}$};\mbox{${\textstyle\frac{3}{2}}$};-X)$
$=\frac{2X}{\sqrt{1+X}}\,{}_{2}F_{1}\biggl{(}\mbox{${\textstyle\frac{1}{2}}$},1;\mbox{${\textstyle\frac{3}{2}}$};\frac{X}{1+X}\biggr{)}=\frac{a}{\sqrt{1-a^{2}}}\,\log\biggl{(}\frac{1+a}{1-a}\biggr{)}$

upon use of Euler’s transformation followed by evaluation of the resulting
hypergeometric function by [2, (15.4.2)]. Hence, we recover the result

$J_{0}(a;\mbox{${\textstyle\frac{1}{2}}$})=\pi$

as given in (3.2).

Appendix: The inversion of (2.5)

We require the inversion of the expression

$\mbox{erfc}\,x=\frac{2u}{s}\biggl{(}1-\frac{u}{2s}+\cdots\biggr{)}\qquad(u\to
0).$

Since the limit $u\to 0$ corresponds to $x\to+\infty$, we employ the
asymptotic expansion [2, (7.12.1)]

$\mbox{erfc}\,x=\frac{e^{-x^{2}}}{\sqrt{\pi}x}\biggl{\\{}1-\frac{1}{2x^{2}}+\frac{3}{4x^{4}}+O(x^{-6})\biggr{\\}}\qquad(x\to+\infty)$

to yield

$xe^{x^{2}}=\frac{s(1-\mbox{${\textstyle\frac{1}{2}}$}x^{-2}+\mbox{${\textstyle\frac{3}{4}}$}x^{-4}-\cdots)}{au(1-u/(2s)+\cdots)},\qquad
a=2\sqrt{\pi}.$

Taking logarithms, we obtain

$x^{2}=\log\,s\biggl{(}1-\frac{\log\,au}{L}-\frac{\log\,x}{L}-\frac{1}{2x^{2}L}+\frac{5}{8x^{4}L}+\cdots+O((sL)^{-1})\biggr{)},$
(A.1)

where we set $L=\log\,s$, $L_{1}=\log\sqrt{\log\,s}$.

Following the iterative procedure described in [1, pp. 25–26] we have as a
first approximation $x^{2}=\log\,s$, whence

$x^{2}=\log\,s\biggl{(}1+\frac{A(u)}{L}\biggr{)},\qquad A(u)=-\log\,au-L_{1}.$

Then with

$\frac{\log\,x}{L}=\frac{L_{1}}{L}+\frac{A(u)}{2L^{2}}+\ldots,$

we find a second approximation given by

$x^{2}=\log\,s\biggl{(}1+\frac{A(u)}{L}+\frac{B(u)}{2L^{2}}+\cdots\biggr{)},\qquad
B(u)=-A(u)-1.$

This last result produces

$\frac{\log\,x}{L}=\frac{L_{1}}{L}+\frac{A(u)}{2L^{2}}+\frac{B(u)-A^{2}(u)}{4L^{3}}+\cdots,$

so that we obtain the third approximation

$x^{2}=\log\,s\biggl{(}1+\frac{A(u)}{L}+\frac{B(u)}{2L^{2}}+\frac{C(u)}{4L^{3}}+O(L^{-4},(sL)^{-1})\biggr{)},$
(A.2)

where

$C(u)=A^{2}(u)+3A(u)+\mbox{${\textstyle\frac{7}{2}}$}.$

Acknowledgement: The author wishes to acknowledge F. Mainardi for bringing to
his attention the 1933 paper by Tricomi.

## References

  * [1] N.G. De Bruijn, Asymptotic Methods in Analysis, Dover, New York, 1981.
  * [2] F.W.J. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, Cambridge, 2010.
  * [3] L. Prandtl and O.G. Tietjens, Applied Hydro- and Aeromechanics, Dover Publications, New York, 1957.
  * [4] F.G. Tricomi, Determinazione del valore asintotico di un certo integrale, Atti della Reale Accad. Nazionale dei Lincei, Rend. Classe Sci. Fis. Mat. Nat., 17 (1933) 116-119; Collected Works, vol. 1, pp. 399–402, Unione Matematica Italiana, 2011.