\section{Supplementary Materials}
\begin{appendix}
\input{SM_raw.tex}
\end{appendix}

\bibliographystyle{apsrev4-1}

\end{document}
\section{Lattice formulation, duality and simulation details}

Using the lattice gauge action from Eq.~(2) of the letter and the standard discretization of bosonic matter, i.e., the Higgs fields couple 
to $U(1)$-valued link variables
$U_\mu(x) = e^{iA_{x,\mu}}$, the lattice-discretized partition function for the $U(1)$ gauge-Higgs model we consider reads
\begin{equation}
Z \; = \; \int \! D[A] \, D[\Phi] \; B_G[A;\theta] \, e^{\, - \, S_m[\Phi,A]} \; ,
\label{Zdef}
\end{equation}
where $B_G[A;\theta]$ is the gauge field Boltzmann factor
\begin{equation}
B_G[A;\theta] \; = \; \prod_{x\in\Lambda} \sum_{n_x\in\mathbb{Z}} 
e^{\, - \, \frac{\beta}{2}\left(F_x+2\pi n_x\right)^2 \, - \, i \theta n_x} \; ,
\label{B_G}
\end{equation}
with $F_x = A_{x+\hat{1},2} - A_{x,2} - A_{x+\hat{2},1} + A_{x,1}$. The matter action $S_M[\Phi,A]$ is given by
\begin{eqnarray}
S_H[\Phi,A] & \; = \; & \sum_{x}\Big[ M \Phi_x^\dagger \Phi_x  + 
\lambda \big(\Phi_x^\dagger \Phi_x\big)^2 
\nonumber \\
& & \quad - \sum_{\mu=1}^{2} 
\big(\Phi_x^\dagger e^{iA_{x,\mu}} \Phi_{x + \hat \mu} \, + \, c.c.\big) \Big] \; ,
\label{S_Higgs}
\end{eqnarray}
with the mass parameter $M$ defined as $M = 4 + m^2$, where $m$ is the bare mass. The variable $x$ runs over all sites
of an $L_1 \times L_2$ square lattice where all fields obey periodic boundary conditions. The measures in the partition function 
(\ref{Zdef}) are the usual product measures over all degrees of freedom on the sites and links of the lattice.  

With the help of Poisson resummation discussed in Eq.~(3) of the letter we can linearize the dependence of the Boltzmann factor (\ref{B_G}) 
on $F_x$ and after an expansion of the nearest neighbor term Boltzmann factors of the Higgs field integrate out the gauge fields
$A_{x,\mu}$ and the matter degrees of freedom $\Phi_x$ (see \cite{Gattringer:2018dlw,Sulejmanpasic:2019ytl} for more details). 

The result is an exact rewriting of the partition sum (1) in terms of flux variables $j_{x,\mu}, k_{x,\mu} \in \mathds{Z}$ 
which describe the two matter field components and the plaquette based 
variables $m_x \in \mathds{Z}$ for the gauge degrees of freedom. In this dual form the partition sum is a sum over the 
configurations of the new variables,
\begin{eqnarray}
\hspace*{-0mm} && Z  =   \sum_{\{m\}}  W_G[m]  \; \sum_{\{j,k\}} W_H[j,k] \;
\prod_x  \delta\!\left(\!\vec{\nabla}\cdot \vec{j}_x\! \right) \delta\!\left(\!\vec{\nabla}\cdot \vec{k}_x\! \right) 
\label{dualZ} \\
\hspace*{-0mm} &&
\times \, \prod_x \!  
\delta\!\left(j_{x,1}\!+\!k_{x,1}\!+\!p_{x}\!-\!m_{x-\hat{2}}\right)  \delta\!\left(j_{x,2}\!+\!k_{x,2}\!-\!m_{x}\!+\!m_{x-\hat{1}}\right)\! .
\nonumber
\end{eqnarray}
In (\ref{dualZ}) the sums over the configurations of the new variables are defined as 
$\sum_{\{m\}} \equiv  \prod_x \sum_{m_x \in \mathds{Z}}$ and 
$\sum_{\{j,k\}} \equiv \prod_{x,\mu} \sum_{j_{x,\mu}\in\mathds {Z}} \sum_{k_{x,\mu}\in\mathds{Z}}$. 
The weight factor $W_G[m]$ for the plaquette occupation numbers $m_x \in \mathds{Z}$ is given by 
($V = L_1 L_2$ denotes the number of sites) 
\begin{equation}
W_G[m] \; = \; \left( 2\pi \beta \right)^{-\frac{V}{2}} 
 e^{\, - \frac{1}{2 \beta} \sum_x\big(m_x + \frac{\theta}{2\pi}\big)^2} \; .
\label{W_G}
\end{equation}
The weight factor $W_H[j,k]$ for the scalar fields is a sum over link-based auxiliary variables 
$a_{x,\mu}, b_{x,\mu} \in \mathds{N}_0$ with  $\sum_{\{a,b\}} \equiv 
\prod_{x,\mu} \sum_{a_{x,\mu}\in\mathds {N}_0} \sum_{b_{x,\mu}\in\mathds{N}_0}$,
\begin{eqnarray}
&& W_H[j,k] =  \sum_{\{a,b\}} \, \prod_x \! I\left(f_x,g_x\right) 
\nonumber \\
&& \hspace{10mm} \times \, \prod_{x,\mu} \, 
\frac{1}{\left(|j_{x,\mu}|\!+\!a_{x,\mu}\right)! \ a_{x,\mu}!} \, 
\frac{1}{\left(|k_{x,\mu}|\!+\!b_{x,\mu}\right)! \ b_{x,\mu}!} \; ,
\nonumber \\
&& 
I\!\left(f_x,g_x\right)  =  \int_{0}^{\infty} \!\!\!\!\! dr \!\!\!  \int_{0}^{\infty} \!\!\!\!\! ds \, r^{\,f_x+1} \, s^{\,g_x+1} \, 
e^{\, - \, M (r^2 + s^2) \, - \, \lambda (r^2 +  s^2)^2 } ,
\nonumber \\
&& 
f_x = \sum_{\mu}
\left[ |j_{x,\mu}| + |j_{x-\hat{\mu},\mu}| + 2\left( a_{x,\mu} + a_{x-\hat{\mu},\mu} \right) \right] \; ,
\nonumber \\
&& 
g_x =\ \sum_{\mu}
\left[ |k_{x,\mu}| + |k_{x-\hat{\mu},\mu}| + 2\left( b_{x,\mu} + b_{x-\hat{\mu},\mu} \right) \right] \; .
\label{W_H}
\end{eqnarray}
Obviously the weight factors (\ref{W_G}) and (\ref{W_H}) are real and positive for all  
$\theta$, such that the sign problem is solved. 

The dual variables $j_{x,\mu}, k_{x,\mu}$ and $m_x$ obey constraints that are written as products of
Kronecker deltas in (\ref{dualZ}), where we use the notation $\delta(n) = \delta_{n,0}$. These constraints 
enforce vanishing divergence $\vec{\nabla}\cdot \vec{j}_x = \sum_\mu [ j_{x + \hat \mu,\mu} - j_{x,\mu}]$ at all sites $x$
for both, the $j$ and the $k$ variables, which implies that they must form closed loops of flux. 
The second set of constraints involves all three types of the dual variables and enforces that for all links the combined oriented flux 
of $j$, $k$ and the plaquette variables $m_x$ on the plaquettes that contain the link must add up to 0. 
  
A suitable Monte Carlo update of the dual form (\ref{dualZ}) must take into account the constraints. In our simulation this is 
implemented by a mix of updates that ensure ergodicity and a reasonably fast decorrelation of the configurations of the dual variables.
More specifically we change plaquette occupation numbers $m_x$ by $\pm 1$ and simultaneously change the flux of $j$ or $k$
around that plaquette. This is combined with worms \cite{Prokofev:2001ddj} for doubly occupied loops where we jointly update  
$j$- and $k$- fluxes of opposite sign, as well as surface worms \cite{Mercado:2013yta} that jointly update flux and plaquette variables. 
Furthermore we include a global change of the plaquette occupation numbers by proposing to
 change all $m_x$ by the same value $\pm 1$, and finally also make use of the $\mathds{Z}_2$ symmetries of our system by performing 
 charge conjugation and flavor swapping updates on our configurations.

We consider various lattice sizes $V = L \times L$ with $L$ ranging from $L = 32$ up to $L = 112$. The couplings
$\lambda$ and $\beta$  are kept fixed at $\lambda = 0.5$  and $\beta = 3.0$ 
for this work and we study the system as a function of the mass parameter 
$M = 4 + m^2$. Typically we use statistical sample sizes of $10^4$ to $10^7$ configurations, which are separated by 
10 to 100 decorrelation sweeps. The initial equilibration is performed with $5 \times 10^5$ sweeps. The error bars are estimated by Jackknife, combined with binning to account for autocorrelations.

\section{Stiffness in compact scalar theory}

The compact scalar Lagrangian in Euclidean space is given by
\begin{equation}
\mathcal L=\frac{J}{2}(\partial_\mu\phi)^2\;,
\end{equation}
where compactness implies the identification $\phi\sim \phi+2\pi$. For computing the spin-stiffness 
we couple the model to a constant gauge field  by replacing $\partial_1\phi \rightarrow \partial_1\phi+A_1$. Expanding the action we find
\begin{equation}
\mathcal S=\int d^2x\;\Big[\frac{J}{2}\left(\partial_\mu\phi\right)^2+J\partial_1\phi A_1+\frac{J}{2}A_1^2\Big]\;.
\end{equation}
We are interested in the theory on a torus with $L_1,L_2$ being the lengths of the two cycles.
Now we write 
\begin{equation}
\phi(x_1,x_2)=\frac{2\pi Q x_1}{L_1}+\tilde\phi(x_1,x_2)\;,
\end{equation}
where $\tilde \phi(x_1,x_2)$ is periodic, and $Q=\int dx^1\partial_1\phi$ is the winding number along the cycle $1$. The action turns into
\begin{equation}
\mathcal S=\frac{J}{2}\int d^2x\; (\partial_\mu\tilde\phi)^2+\frac{JL_2}{2L_1}\left({2\pi Q}+A_1L_1\right)^2 .
\end{equation}
Obviously the  field $\tilde\phi$ decouples from the rest, such that we can consider the partition function
\begin{equation}\label{eq:Z}
Z=\sum_{Q\in \mathbb Z}e^{-\frac{JL_2}{2L_1}\left({2\pi Q}+A_1L_1\right)^2}\;.
\end{equation}
Note that the sum above can be expressed in terms of the Jacobi theta function $\vartheta_3$, but we will not need this form for what follows.

The stiffness is given by
\begin{equation}
\rho_s=-\frac{1}{L_1L_2}\frac{\partial^2\log Z}{\partial A_1^2}\Big|_{A_1=0}=-(2\pi)^2J^2 \frac{L_2}{L_1}\avg{Q_1^2}+J\;,
\end{equation}
where the average is taken with the partition function \eqref{eq:Z} at $A_1=0$. 

If the combination $\frac{JL_2}{L_1}$ is sufficiently large, all the values of $Q$ away from zero are exponentially suppressed in (\ref{eq:Z}), 
and we immediately find $\rho_s=J$. This is the usual result, which is valid to significant accuracy, for the usual BKT transition for $L_1=L_2$, that happens at $J=\frac{2}{\pi}$ \cite{prokof2000two}. Here, however, we are interested in the point $J=\frac{1}{2\pi}$, which is not sufficiently large to make the above approximation. The plot in Fig.~\ref{fig:fig1} shows the comparison of the stiffness result from the full partition sum (\ref{eq:Z}) at 
$L_1=L_2$ with the function $J$ to illustrate this.

On the other hand the partition function \eqref{eq:Z} is periodic with respect to the magnetic flux $A_1L_1$ with a period $2\pi$, as can be seen by inspection. The Fourier expansion thus is given by
\begin{equation}
Z=\sum_{w\in \mathbb Z} \sqrt{\frac{L_1}{2\pi JL_2}}e^{-\frac{w^2L_1}{2JL_2}+i w A_1L_1}\;.
\end{equation}
The stiffness therefore reads
\begin{equation}\label{eq:Zdual}
\rho_s=\frac{L_1}{L_2}\avg{w^2}\;,
\end{equation}
where again the average is taken at $A_1=0$. Now if $\frac{JL_2}{L_1}$ is sufficiently small, we find that $\rho_s=O(e^{-\frac{L_1}{2JL_2}})$. 


We are interested in the case $L_1=L_2=L$ and $J=\frac{1}{2\pi}$, which is between the two extremes of small and large $\frac{JL_2}{L_1}$. However, in this case it is easy to see that $\avg{w^2}=\avg{Q^2}$. Now identifying \eqref{eq:Z} and \eqref{eq:Zdual}, we find
$\avg {Q^2}=\avg{w^2}=\frac{1}{4\pi}$, and obtain
\begin{equation}
\rho_s=\frac{1}{4\pi}\;.
\end{equation}



\begin{figure}[tb]
   \centering
   \includegraphics[width=.44\textwidth]{SM_fig1.pdf} 
   \caption{Plot of the stiffness from the full partition sum (\ref{eq:Z}), and comparison with the result 
   $\rho_s\approx J$ valid for large $\frac{JL_2}{L_1}$.}
   \label{fig:fig1}
\end{figure}



