\section{Supplemental material}

\maketitle


\subsection{Projective QMC approach}

We used the  projective QMC algorithm of the ALF-library \cite{ALF17}.
This canonical  algorithm  filters out the ground  state, $|\psi_0\rangle $,
from a trial wave function,  $|\psi_T  \rangle$,  that is required to be non-orthogonal to the ground state:
\begin{equation}
 \frac{  \langle \psi_0 |  \hat{O} | \psi_0  \rangle   }{\langle \psi_0 |   \psi_0  \rangle } = \lim_{\Theta \rightarrow \infty}
  \frac{ \langle \psi_T |  e^{- \Theta \hat{H}}  \hat{O}  e^{- \Theta \hat{H}}  | \psi_T \rangle }{
   \langle \psi_T |  e^{ - 2\Theta \hat{H} } | \psi_T \rangle   }.
\end{equation}

The trial wave function $ | \psi_T \rangle  $ is chosen to be  a Slater determinant  with
$N_\text{p}$ particles ( $ \hat{N} |\psi_T \rangle = N_\text{p} | \psi_T \rangle  $ ).    In particular,

\begin{equation}
|\psi_T \rangle =  |\psi_T^{\uparrow} \rangle   \otimes  |\psi_T^{\downarrow} \rangle
\end{equation}
with
\begin{equation}
  |\psi_T^{\sigma} \rangle \equiv   \prod_{n=1}^{N_\text{p}/2} \left(  \sum_{\bm{i}} \hat{c}^{\dagger}_{\bm{i},\sigma}  U_{ \bm{i}, n} \right)   | 0 \rangle.
\end{equation}
 $ U_{ \bm{i}, n } $ is the $n^{\text{th}}$  single-particle  eigenstate, ordered in ascending   energy eigenvalues,  of   the spinless fermion Hamiltonian
\begin{equation}
  \hat{H}  = - t  \sum_{ \langle \bm{i}, \bm {j} \rangle } (\hat{c}^{\dagger}_{\bm{i} } \hat{c}^{\phantom\dagger}_{\bm{j}}  + H.c.)
  +   \sum_{ \langle \bm{i}, \bm {j} \rangle }   \xi_{ \bm{i}, \bm{j} }
   (\hat{c}^{\dagger}_{\bm{i} } \hat{c}^{\phantom\dagger}_{\bm{j}}  + H.c.).
\end{equation}
The first term corresponds to the tight-binding Hamiltonian on the honeycomb lattice.
We require  the perturbing  hopping  matrix elements  $ | \xi_{ \bm{i}, \bm{j}} |  \ll  t $  and  $\Im {\xi_{\bm{i}, \bm{j} } } = 0 $.
The sign and modulus of  $\xi_{\bm{i}, \bm{j} }$ are chosen randomly so that all energy eigenvalues of the spinless Hamiltonian are non-degenerate.
Our trial wave function hence breaks  lattice and point group  symmetries.
Crucially, however,  time-reversal symmetry is present.   Since $\lambda > 0$    (see Eq.~(\ref{Eq:Ham})), we can decouple  the  interaction with a real Hubbard-Stratonovich  transformation  such that both the  imaginary time propagation and the  trial wave function  are invariant under  time reversal:
\begin{equation}
  T \alpha  \begin{pmatrix} \hat{c}_{\ve{i}, \uparrow} \\  \hat{c}_{\ve{i}, \downarrow}  \end{pmatrix}  T^{-1}
 =  \Bar{\alpha}    \begin{pmatrix} \hat{c}_{\ve{i}, \downarrow} \\ - \hat{c}_{\ve{i}, \uparrow}  \end{pmatrix}.
\end{equation}
Hence, the     eigenvalues of the fermion matrix  come in   complex conjugate pairs  and no  negative sign problem occurs.

A projection length $\Theta =L $ was found to be sufficient to converge to the finite-size ground state
for all of our system sizes.
We have used an imaginary time  step   $\Delta_{\tau}=0.2$
and  a  symmetric  Trotter  decomposition  to
guarantee the Hermiticity of the imaginary time propagator.  All calculations were carried out at  $\lambda= 0.026$   in units where $t=1$. 


\subsection{Equal-time structure factor}

In Fig.~\ref{fig:S_K}, we show the  momentum dependence of the  equal-time  QSH and SSC structure factors  at $\delta =0$ and  at
$\delta = 1/36$.   Upon doping, the QSH structure factor does not develop
incommensurate   features.   At $\delta = 1/36$,  the QSH data
(Fig.~\ref{fig:S_K}(b)) are consistent with the absence of long-range order,
whereas the SSC  structure factor (Fig.~\ref{fig:S_K}(d)) shows a marked
increase as a function of system size.

The onset of long-range order as well as a measure for the correlation length
can be obtained by considering $1/S(\bm{Q}=0)$ as function of $\delta$ (see
Fig.~\ref{fig:S_delta}).  The SSC ordering appears immediately at $\delta>0$,
characterized by the quick decay of $1/S_{\text{SSC}}$ as function of system
size. In particular, $1/S_{\text{SSC}}$ shows no saturation as a function of
system size.  On the other hand, $1/S_{\text{QSH}}$ shows a clear saturation
at \textit{large} doping.  For a given doping, the lattice
size at which this quantity converges is a measure of the correlation
length. Upon inspection, one will see that
larger lattice sizes are required to achieve convergence upon approaching half-filling.  In particular,
following the envelope of these curves again suggests that the correlation
length of the QSH fluctuations grows continuously and diverges as
$\delta \rightarrow 0$.  This is consistent with the data shown in the main
text.

\begin{figure}[h]
\centering
\includegraphics[width=0.48\textwidth]{S_K.pdf}
\caption{\label{fig:S_K}
Momentum dependence of the equal-time structure factor for ((a),(b)) QSH and
((c),(d)) SSC operators for ((a),(c)) $\delta=0$ and
and ((b),(d))  $\delta=\frac{1}{36}$.
}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.48\textwidth]{S_delta.pdf}
\caption{\label{fig:S_delta}
$\delta$ dependence of $1/S(\bm{Q}=0)$ for SSC (a) and QSH (b), as a function of $\delta$, for $L=9$, 12, 15, 18, 21 and $24$}
\end{figure}

\subsection{Consistency check of the pairing gap}

\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{Pair_Gap_N.pdf}
\caption{\label{fig:Gap_Energy}
Three ways of pairing gap evaluation for (a) $L=9$, (b) $12$, (c) $15$, and (d) $18$, for doped particle number
$N_\text{p} - 2L^2 = 0, -2, -4, -6$ and $-8$.
}
\end{figure}

We check the consistency of our evaluation of  the ground-state energy difference between different even particle-number
sectors:
\begin{equation}
  E(N_\text{p}) - E(N_\text{p} -2) = \Delta_{ \eta^{-} }(N_\text{p})  = \Delta_{ \eta^{+} }(N_\text{p} -2)
\end{equation}
where $E(N_\text{p})$ is the ground state energy measured within the PQMC in the  $N_\text{p}$ particle number sector;
$ \Delta_{ \eta^{-} }(N_\text{p})$ is the s-wave pairing ($\eta^{-}$) gap extrapolated from
the time-displaced correlation function.

The imaginary-time domain $\beta$, in which we measure the time-displaced
correlation function, is set to $\beta=L$ for $L=9$, 12, 15, and $18$, and to $\beta=10$ for $L=21$.
To extrapolate the pairing gap, we use sequential fits
\begin{equation}
 \langle \eta^{+}( m \tau_0 + \tau ) \eta^{-} (m \tau_0)  \rangle  \propto  e^{ -\Delta_m \tau }   \  \  \   m = 0,1,2,3...
\end{equation}
where $\tau_0=1.0$ and $\tau \in [0, \tau_0)$.  The gap is extrapolated as
\begin{equation}
 \Delta_m - \Delta( m \rightarrow \infty )  \propto  e^{-a m}
\end{equation}
where $a$ is optimized for the best fit.

In  Fig.~\ref{fig:Gap_Energy}, we show that   the three different  ways of evaluating the  gap give  consistent results  for
$L=9$, 12, 15, and 18 for several particle-number sectors near half-filling.
In particular, one can compute the  ground-state energy and take the
difference or  measure time-displaced correlation functions of the pair
adding or removal operator.   Using the energy difference  generically
produces  bigger error bars. Here, we  carry out two independent simulations
and thereby  have to add the errors on two extensive  quantities (total
energies) to estimate the error on an intensive one, the total energy
difference.   Hence to keep the error bar on the total energy
difference, we have to scale the error on the energy per site as $1/L^2$.  Even
taking into account self-averaging on large system sizes, this proves to  be
numerically expensive.



\subsection{Gap extrapolation at half filling}

In this section  we discuss the calculations carried out to determine the  phase boundary of Fig.~\ref{fig:Phase_diagram}.  
As explained in the main text, the critical chemical potential is obtained from half  the pairing gap $\Delta_{\eta}$ at half filled case.   

Fig.~\ref{fig:Gap_Pair_Half_Filling} shows the extrapolation of  the  pairing gap $ \Delta_{\eta} $ at $\mu = 0$, for four different values of $\lambda$  
inside the  QSH phase.  An exponential finite size behavior is assumed for the extrapolation:  
\begin{equation}
    \Delta_{\eta} (L)  =   \Delta_{\eta} ( L \longrightarrow \infty )  + a e^{ -L / \xi } 
\end{equation}
For $\lambda=0.021, 0.023$ and $0.028$, an acceptable $\chi^2$ is obtained for a collective fit using sizes $L=6,9,...,21$.  
For $\lambda=0.031$,  an acceptable $\chi^2$ is obtained when the  $L=6$  data is omitted.  


Additionally we show the single particle gap at $\mu=0$  as extracted from  the imaginary time displaced  Green function:    
\begin{equation}
\begin{aligned}
    \langle  c_{\bm{k} } (\tau)  c^{\dagger}_{ \bm{k}, \alpha }(0) \rangle  \propto  e^{ -\Delta_{ \text{sp} } \tau }. 
\end{aligned}
\end{equation}
Due to the shift of the minimal gap   in  momentum space  at   large $\lambda$,     
 we choose $ \bm{k}  \equiv ( \frac{4\pi}{3}, 0 ) $ for $ \lambda = 0.021, 0.023$ and $0.026$,  and $ \bm{k}  \equiv ( \pi, \frac{\pi}{ \sqrt{3}} ) $ for 
$ \lambda = 0.028$ and $0.031$.  

To extrapolate to the thermodynamic limit, we have again used an exponential fit,
\begin{equation}
    \Delta_{ \text{sp} } (L)  =   \Delta_{ \text{sp} } ( L \longrightarrow \infty )  + a e^{ -L / \xi },
\end{equation}
that is supported by the data. 

Before proceeding, we would like to comment on the exponential fit, since it is unexpected. In principle,  a doped hole (or pair)  will couple to the Goldstone modes  
such that  one expects the finite size effects of the gap to pick up the finite size  behavior of the Goldstone mode.  The latter follow a polynomial law  in $1/L$.  
This behavior is explicitly seen in the context of antiferromagnetic Mott insulators \cite{Assaad96a}.   If on the other hand the doped  charge carriers do not couple 
strongly to the low-lying  Goldstone modes,  then we expect dominant exponential finite size effects.  This point of view will be confirmed  in the  next section. 



\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{Pair_Diagram.pdf}
\caption{\label{fig:Gap_Pair_Half_Filling} 
Pairing  and single particle gap extrapolation at half-filling,  for $\lambda=0.021$, $0.023$, $0.028$ and $0.031$.  
}
\end{figure}




\subsection{Single-particle spectrum at finite doping}

We consider the single-particle spectral function at finite doping.    Away
from half-filling,  particle-hole symmetry is broken  and we have to
separately calculate the spectra for electron addition and removal,
\begin{equation}
\begin{aligned}\label{Eq:Spectrum}
A(\bm{k}, \omega) & = \frac{1}{Z}  \sum_n
  ( | \langle  n | c_{\bm{k}} | 0 \rangle |^2  \delta(E_n - E_0 - \omega)  )  \\
& +  \frac{1}{Z}  \sum_m ( | \langle  m | c^{\dagger}_{\bm{k}} | 0 \rangle |^2  \delta(E_m - E_0 + \omega)  )
\end{aligned}
\end{equation}
via the independent analytical continuations
\begin{equation}
\begin{aligned}
  &  \langle  c_{\bm{k} } (\tau)  c^{\dagger}_{ \bm{k}, \alpha }(0) \rangle  =  \int d\omega  e^{ -\tau \omega }  A_{+}(\omega )      \\
  &  \langle  c^{\dagger}_{\bm{k} } (\tau)  c_{ \bm{k}, \alpha }(0) \rangle  =  \int d\omega  e^{ -\tau \omega }  A_{-}(\omega )
\end{aligned}
\end{equation}
with $A(\omega)=A_{+}(\omega) + A_{-}(-\omega)$.  Here, $|0 \rangle$ in
Eq.~(\ref{Eq:Spectrum}) is the ground state at finite doping and $\langle n |$
is an eigenstate of the Hamiltonian with energy $E_n$ and an additional particle
(hole) relative to the ground state.  In
Fig.~\ref{fig:Single_particle_spectrum}, we plot the spectral functions for $L=21$ and
$\delta=0, \frac{1}{441}, \frac{3}{441}$ and $\frac{4}{441}$
($2L^2-N_\text{p} =0,2,6$ and $8$).  The dominant feature follows the
mean-field BCS form $E(\ve{k}) = \pm \sqrt{ (\epsilon(\ve{k}) -\mu )^2 + |\Delta|^2} $, where
$\pm \epsilon(k)$ denotes the Dirac dispersion of the honeycomb lattice.
This result shows that the Goldstone modes do not strongly couple to single-particle excitations.

In fact, the Green's function at the $\Gamma$ point has a special property,
due to a commutation rule between the fermion operator and the interaction
term of Hamiltonian (we use the notation $ \ve{c}^{\dagger}_{\ve{i}} =  ( \hat{c}^{\dagger}_{\bm{i}, \uparrow },\hat{c}^{\dagger}_{\bm{i}, \downarrow } )$)
\begin{equation}\label{Eq:Conserve_Int}
\begin{aligned}
\left[ \sum_{ \bm{i} }  \hat{c}_{\bm{i}, \alpha }  ,  \sum_{\varhexagon}  \left( \sum_{\langle \langle \bm{i} \bm{j} \rangle \rangle
\in  \varhexagon }  i \nu_{ \bm{i} \bm{j} }
   \hat{\ve{c}}^{\dagger}_{\bm{i}} \bm{\sigma} \hat{\ve{c}}^{\phantom\dagger}_{\bm{j}}  + H.c. \right)^2    \right]      =  0.
\end{aligned}
\end{equation}
The above relation follows directly from
\begin{equation}\label{Eq:Conserve_Hexagon}
\begin{aligned}
\left[ \sum_{ \bm{i} }  \hat{c}_{\bm{i}, \alpha }  ,    \sum_{\langle \langle \bm{i} \bm{j} \rangle \rangle  \in  \varhexagon  }  i \nu_{ \bm{i} \bm{j} }
   \hat{\ve{c}}^{\dagger}_{\bm{i}} \bm{\sigma} \hat{\ve{c}}^{\phantom\dagger}_{\bm{j}}  + H.c.      \right]      =  0
\end{aligned}
\end{equation}
which holds for the summation of spin-orbit  operators inside each hexagon
and for any vector $\ve{\sigma}$ in the space of Pauli matrices $\sigma_x$, $\sigma_y$, $\sigma_z$.
Hence, the Green's function  $\langle  \hat{\ve{c}}^{\dagger}_{\bm{k}} (\tau)  \hat{\ve{c}}^{}_{ \bm{k} }(0) \rangle$  at the  $\bm{k} = \Gamma$ point is identical to that
of the non-interacting Hamiltonian.

\begin{figure}
\centering
\includegraphics[width=0.50\textwidth]{Green_dop_L21.pdf}
\caption{\label{fig:Single_particle_spectrum}
Single-particle spectrum at dopings  (a) $\delta=0$, (b)
$\frac{1}{441}$, (c) $\frac{3}{441}$, and (d) $\frac{4}{441}$.
The green dotted line is the  chemical potential $ \mu $ evaluated from
Eq.~(\ref{Eq:cp-pqmc}) of the main text.
}
\end{figure}



\subsection{Finite-temperature calculation}

\begin{figure}
\centering
\includegraphics[width=0.50\textwidth]{Dop_FT.pdf}
\caption{\label{fig:Dop_FT}
Doping factor $ \delta $ as a function of chemical potential $\mu$ from FTQMC
simulations, with sizes $L=9$, $12$ and $15$, and $\beta=\frac{1}{3}L^2$.
}
\end{figure}

In contrast to the  projective  approach, the finite-temperature auxiliary
field QMC (FTQMC) algorithm is formulated in the grand-canonical ensemble.
Adding a chemical potential  term, $-\mu   \sum_{ \bm{i} } \hat{\ve{c}}^{\dagger}_{\bm{i} } \hat{\ve{c}}^{\phantom\dagger}_{\bm{i}} $,  to  the Hamiltonian, see  Eq.~(\ref{Eq:Ham}) of the
main text,  allows us to compute the doping  away from half-filling ($\mu = 0$):
 \begin{equation}
  \delta \equiv  \frac{ \langle \sum_{\bm{i}}  \hat{n}_{\bm{i}}  \rangle }{2 L^2} - 1.
 \end{equation}
Figure~\ref{fig:Dop_FT} shows the corresponding result for the case where the
inverse temperature $\beta$ and the system size $L$  scale  as $\beta=L^{2} /3$.
This implicitly  makes the assumption that $z=2$. Overall, our limited data
are consistent with the more efficient PQMC calculation.  At large values of
$\delta$, results for $L=12$ and $L=15$ are consistent with a linear dependence of
$\delta$ on $\mu$ that overshoots the critical chemical potential and suggest $z > 2$.

From the numerical point of view, the   FTQMC    is not as  efficient as the
PQMC.    The numerical cost to reach the low-temperature limit  scales as
$V^{3} \beta^{z} $.    We have also noticed  long warm-up  and autocorrelation
times  to equilibrate the particle number in the vicinity of $\mu = \mu_c $.

\subsection{Mean-field calculation}

In this section, we summarize the details of our mean field calculation.
Expanding interacting part of Eq.~(\ref{Eq:Ham}) of the main text as
\begin{equation}
\begin{aligned}\label{Eq:MF_decouple}
 H_V = &  -\lambda \sum_{\varhexagon}  \left(\sum_{\langle \langle \bm{i} \bm{j} \rangle \rangle }  i \nu_{ \bm{i} \bm{j} }
   \hat{c}^{\dagger}_{\bm{i}} \bm{\sigma} \hat{c}_{\bm{j}}  + H.c. \right)^2    \\
   = &  -\lambda  \sum_{\varhexagon}  \sum_{ \langle \langle \bm{i} \bm{j} \rangle \rangle }
   \sum_{ \langle \langle \bm{i'} \bm{j'} \rangle \rangle \neq  \langle  \langle \bm{i} \bm{j} \rangle \rangle  }
  \bm{ \hat{J}_{ \langle \langle i, j \rangle \rangle }  \cdot  \hat{J}_{ \langle \langle i', j' \rangle \rangle }  }    \\
   & - \lambda \sum_{\varhexagon}  \sum_{ \langle \langle \bm{i} \bm{j} \rangle \rangle }
   [ + 6 \hat{\eta}^{\dagger}_{\bm{i}}  \hat{\eta}_{\bm{j}} + h.c  -4 \bm{ \hat{S}_i } \cdot \bm{ \hat{S}_j }   \\
   &  -5 \hat{n}_{\bm{i}} \hat{n}_{\bm{j}}  +  5 ( \hat{n}_{\bm{i}} + \hat{n}_{\bm{j}} )    ]
\end{aligned}
\end{equation}
\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{Half_fill.pdf}
\caption{\label{fig:MF_Half}
Mean-field solution as a function of $\lambda$ at half-filling. (a) QSH and SSC order parameters. (b) Fermionic single-particle gap.
}
\end{figure}
where
\begin{equation}
\begin{aligned}
 &  \bm{ \hat{J} }_{ \langle \langle i, j \rangle \rangle }  \equiv    i \nu_{ \bm{i} \bm{j} }
  \hat{c}^{\dagger}_{\bm{i}} \bm{\sigma} \hat{c}_{\bm{j}}  + H.c.,         \\
 &  \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},    \\
 &  \bm{ \hat{S}_i }  \equiv  \frac{1}{2}\hat{c}^{\dagger}_{\bm{i}}  \bm{\sigma} \hat{c}_{\bm{i}}.
\end{aligned}
\end{equation}
The self-consistent calculation is  based on selecting a polarization direction
for the three  (two) components of the QSH (SSC) order parameter.
The calculation is done by numerically minimizing the free energy in the space of the two order parameters
\begin{equation}
\begin{aligned}\label{Eq:MF_F}
  f(\beta)_{\phi}  =  \frac{-1}{\beta V} \ln{ \mathop{\mathrm{Tr}}{ e^{-\beta (H_T + H_V) - 15 \beta  V \lambda    {\phi_{\text{QSH}}}^2   -  36 \beta V \lambda    {\phi_{\text{SSC}}}^2    }  } }
\end{aligned}
\end{equation}
where
\begin{equation}
\begin{aligned}
 H_T = & - t  \sum_{ \langle \bm {i}, \bm {j} \rangle } (\hat{\ve{c}}^{\dagger}_{\bm{i} } \hat{\ve{c}}^{\phantom\dagger}_{\bm{j}}  + H.c.  )
   + \mu \sum_{\bm{i}} \hat{\ve{c}}^{\dagger}_{\bm i} \hat{\ve{c}}^{\phantom\dagger}_{\bm i}     \\
 H_V
   = &  - 5  \lambda  \sum_{\varhexagon}  \sum_{ \langle \langle \bm{i} \bm{j} \rangle \rangle }
     \phi_{\text{QSH}}
    \cdot  \hat{J^z}_{ \langle \langle i, j \rangle \rangle }       \\
   & -  36 \lambda   \sum_{ \bm{i} }
      \phi_{\text{SSC}}  \hat{\eta}^{x}_{\bm{i}}.
\end{aligned}
\end{equation}

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{lamb_01.pdf}
\caption{\label{fig:MF_l0.1}
Mean-field solution as a function of chemical potential $\mu$ at $\lambda=0.1$. (a) QSH and SSC order parameters. (b) Doping factor $\delta$.
}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.42\textwidth]{lamb_02.pdf}
\caption{\label{fig:MF_l0.2}
Mean-field solution as a function of chemical potential $\mu$ at $\lambda=0.2$. (a) QSH and SSC order parameters. (b) Doping factor $\delta$.
}
\end{figure}

We  consider a   paramagnetic saddle point  with   $ \langle \hat{\ve{S}}_{\ve{i}}  \rangle = 0 $  and $\langle \hat{n}_{\bm{i}} \rangle = 1 - \delta $.
Thus, for any local minimum of Eq.~(\ref{Eq:MF_F}) with $ \frac{\partial f}{
  \partial \phi_{\text{QSH}}} = 0 $ and  $ \frac{\partial f}{ \partial
  \phi_{\text{SSC}}} = 0 $,
\begin{equation}
\begin{aligned}
 & \phi_{\text{QSH}}  =  \frac{1}{6} \sum_{ \langle \langle \bm{i} \bm{j} \rangle \rangle }  \langle   \hat{J^z}_{ \langle \langle i, j \rangle \rangle }   \rangle  \\
 & \phi_{\text{SSC}}  =  \frac{1}{2} \langle  \hat{\eta}^{x}_{\bm{i},A} + \hat{\eta}^x_{\bm{i},B}    \rangle
\end{aligned}
\end{equation}
which holds locally due to  translational  symmetry.  We numerically integrated over the Brillouin zone of an $L=120$ lattice and took the zero-temperature limit $\beta \rightarrow \infty$.

The two order parameters as a function of $\lambda$ in the  half-filled case
are shown in Fig.~\ref{fig:MF_Half}(a).  We  observe a semimetal
($\phi_{\text{QSH}} = \phi_{\text{SSC}}=0$ ), a pure QSH  state
($\phi_{\text{QSH}} \neq 0, \phi_{\text{SSC}}=0$) as well as a coexistence
(QSH+SSC) state ($\phi_{\text{QSH}} \neq 0, \phi_{\text{SSC}} \neq 0$).
Since charge conservation is a  protecting symmetry of  the QSH  insulator,
the transition  between the QSH and QSH+SSC  states can be continuous without
a closing of the single-particle gap (see Fig.~\ref{fig:MF_Half}(b)).

On the other hand, upon doping the pure QSH state, the phase diagram  exhibits
two distinct mean-field scenarios.  Two representative examples at
$\lambda=0.1$ and $\lambda=0.2$ are shown in Figs.~\ref{fig:MF_l0.1} and
\ref{fig:MF_l0.2}, respectively.  In the case of $\lambda=0.1$
(Fig.~\ref{fig:MF_l0.1}), which is close to the Gross-Neveu transition, a
clear first-order transition between the QSH and SSC phases is observed.
Doping at $\lambda=0.2$ (Fig.~\ref{fig:MF_l0.2}) leads to two phase
transitions: (i) a $z=2$ transition from the pure QSH state to the coexistence
state at $\mu \approx 0.4$, characterized by a linear growth of $\delta$ and
(ii) a first-order phase transition to an SSC state at $\mu \approx 0.8$.  Such
first-order transitions are characterized by a level crossing corresponding
to two local minima in the free-energy density in Eq.~(\ref{Eq:MF_F}).

\end{document}


