\section{\label{secSm:action}Renormalization group and effective action}
In this Section, we present renormalization group (RG) equations for the cosine perturbations appearing in the main text and we provide details for the calculation of the effective action for the Hamiltonian $H_{\rm DSG}$ defined in Eq.~(7) of the main text.

In order to derive the RG equations, one has to specify the form of the small-momentum interaction matrix $N$ appearing in the Hamiltonian $H_0$,
\begin{equation}
H_0 = \sum_{\sigma,\sigma',r,r'}\int dx \hspace{1mm}n_{r \sigma}(x)N_{r\sigma r'\sigma'}n_{r' \sigma'}(x),\label{eq:H}
\end{equation}
where $n_{R\sigma}(x)=R_{\sigma}^{\dagger}(x)R_{\sigma}(x)$ and $n_{L\sigma}(x)=L_{\sigma}^{\dagger}(x)L_{\sigma}(x)$ are the electronic densities for each channel. The $4\times 4$-matrix $N$ is given by
\begin{equation}
N = \left(
\begin{matrix}
v + g_{\parallel} & g_{\perp} & g_{\parallel} & g_{\perp}\\ g_{\perp} & v + g_{\parallel} & g_{\perp} & g_{\parallel} \\g_{\parallel} & g_{\perp} & v + g_{\parallel} & g_{\perp}\\ g_{\perp} &  g_{\parallel} & g_{\perp} & v + g_{\parallel} 
\end{matrix}\right),
\end{equation}
where $v=2\mu_{1/2}/k_F$ is the velocity and $g_{\parallel}$ ($g_{\perp}$)  the coupling constants for density-density interaction between electrons with the same (opposite) spins, respectively. It is useful to express the Hamiltonian in the basis in which $H_0$ is diagonal, which is given by the following bosonic fields:
\begin{equation}
\left(\begin{matrix}
\theta_{\rho}(x) \\ \phi_{\rho}(x)\\\theta_{\sigma}(x)\\\phi_{\sigma}(x) 
\end{matrix}\right) = \frac{1}{2\sqrt{2}}\left(\begin{matrix}
1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & -1 & 1
\end{matrix}\right) \left(\begin{matrix}
\phi_{L\uparrow}(x) \\ \phi_{R\uparrow}(x)\\\phi_{L\downarrow}(x)\\\phi_{R\downarrow}(x) 
\end{matrix}\right).
\end{equation}
In order to reproduce the correct commutation relations among the two perturbations in the original fermionic representation [see Eq.~(2) in the main text], the bosonic fields have to satisfy the following commutation relations:
\begin{align}
&\left[\phi_{\rho}(x),\theta_{\rho}(x')\right]=i \pi\ \text{sign}(x-x'),\\
&\left[\phi_{\sigma}(x),\phi_{\sigma}(x')\right]=-i \pi\ \text{sign}(x-x'),\\
&\left[\theta_{\sigma}(x),\theta_{\sigma}(x')\right]=i \pi\ \text{sign}(x-x').
\end{align} 
The commutators between  right/left mover fields are given by $\left[\phi_{r\sigma}(x),\phi_{r'\sigma'}(x')\right] = i \pi M_{r\sigma ,r' \sigma'} \text{sign}(x-x')$. Here, the $4\times 4$-matrix $M$ can be expressed in the basis $\left(\phi_{R\uparrow},\phi_{R\downarrow},\phi_{L\uparrow},\phi_{L\downarrow}\right)$ as
\begin{equation}
M = \left(\begin{matrix}
1& 1& 1 & -1 \\ 1 & 1 & -1 & 1 \\ 1 & -1 & -1 & -1 \\ -1 & 1 & -1 & -1 
\end{matrix}\right).\label{eq:boso_comm_matrix_sol}
\end{equation}
The Hamiltonian becomes 
\begin{align}
H &= H_0 + \frac{\Lambda_1}{4\pi^2a^2} \int dx ~ \cos\left[\sqrt{2}\left(2\phi_{\rho} - \theta_{\sigma} + \phi_{\sigma}\right)\right] + \frac{\Lambda_2}{4\pi^2a^2} \int dx ~ \cos\left[\sqrt{2}\left(2\phi_{\rho} - \theta_{\sigma} - \phi_{\sigma}\right)\right],
\end{align}
where
\begin{equation}
H_0 = \frac{1}{2\pi} \sum_{\mu=\rho,\sigma}\, \int dx \left[ u_{\mu} K_{\mu} ~\left[\partial_x\theta_{\mu}(x)\right]^2 +  \frac{u_{\mu}}{K_{\mu}}  ~\left[\partial_x\phi_{\mu}(x)\right]^2 \right].
\end{equation}
The RG equations for the two cosine perturbations are given by
\begin{align}
\frac{d}{d l}\tilde{\Lambda}\left(l\right) &= \left[2-2 \tilde{K}_{\rho}(l)-\frac{1}{2}\left(\tilde{K}_{\sigma}(l)+\tilde{K}^{-1}_{\sigma}(l)\right)\right]\tilde{\Lambda}(l)\label{eq:RG1}\\\frac{d}{d l}\tilde{K}_{\rho}\left(l\right) &=- 4\tilde{K}^2_{\rho}(l)\tilde{\Lambda}^2(l),\\\frac{d}{d l}\tilde{K}_{\sigma}\left(l\right) &= \left[1-\tilde{K}^2_{\sigma}(l)\right]\tilde{\Lambda}^2(l).
\end{align}
Here, we introduce $\Lambda\equiv\Lambda_1=\Lambda_2$ and use the rescaled amplitude $\tilde{\Lambda}= \Lambda/v$.
These equations have been used to derive the phase diagram shown in Fig.~2 of the main text. Importantly, we note that the Luttinger liquid parameters
$K_{\rho}$ and $K_{\sigma}$
are flowing under the RG and, according to the above equations, they flow to the values 
\begin{align}
K_{\rho}&\rightarrow 0,\\K_{\sigma}&\rightarrow 1.
\end{align}
In order to derive the effective action, we change the basis as follows:
\begin{equation}
\left(\begin{matrix}\eta_1\\ \eta_2\\\eta_3\\\eta_4\end{matrix}\right)  =  \frac{1}{2\sqrt{2}}\left(\begin{matrix}
0 & -2  & 1& -1 \\ 0 & -2 & 1 & 1 \\ 0 & 1 & 0 & 0 \\ \frac{1}{2}& 0 &1 & 0 
\end{matrix}\right)\left(\begin{matrix}\theta_{\rho}\\\phi_{\rho}\\ \theta_{\sigma}\\ \phi_{\sigma} \end{matrix}\right).
\end{equation}
The Hamiltonian becomes
\begin{equation}
H =\frac{1}{2\pi} \int dx ~\partial_x\vec{\eta}^T(x) \tilde{N}~ \partial_x\vec{\eta}(x)+ \frac{\Lambda_1}{4\pi^2a^2} \int dx ~ \cos\left[4\eta_1(x)\right] + \frac{\Lambda_2}{4\pi^2a^2} \int dx ~ \cos\left[4\eta_2(x)\right],
\end{equation}
where
\begin{equation}
\tilde{N} =\left(
\begin{array}{cccc}
v_1 & g_{12} & g_{13} &g_{14} \\
g_{12} &v_2 & g_{23} & g_{24} \\
g_{13} & g_{14} & v_3 &g_{34} \\
g_{14} & g_{24} & g_{34} & v_4 \\
\end{array}
\right)
\end{equation}
and
\begin{align}
v_1 &= v_2 = \frac{2K_{\sigma }^2 u_{\sigma}+8 u_{\rho}  K_{\rho } K_{\sigma }+2u_{\sigma}}{ K_{\sigma }},\\
g_{12}& =  \frac{2K_{\sigma }^2 u_{\sigma}+8 u_{\rho}  K_{\rho } K_{\sigma }-2u_{\sigma}}{ K_{\sigma }},\\
g_{13} &= g_{23} =8 K_{\sigma}u_{\sigma}+32 K_{\rho
}u_{\rho},\label{eq:param1}\\
g_{14} &= g_{24} =  -16K_{\rho }u_{\rho},\\
v_3 & =128
K_{\rho }u_{\rho}+32K_{\sigma }u_{\sigma}+8\frac{u_{\rho}}{K_{\rho }},\\
v_4 &=32 K_{\rho}u_{\rho},\\
g_{34} &=   -64K_{\rho }u_{\rho}.\label{eq:param2}
\end{align}
The commutators between the bosonic fields are given by
\begin{align}
\left[\eta_1(x),\eta_2(x')\right] &= i \frac{\pi}{4} \text{sign}(x-x'),\\
\left[\eta_3(x),\eta_4(x')\right] &=  i  \frac{\pi}{16}\text{sign}(x-x'),\\
\left[\eta_4(x),\eta_4(x')\right] &=i  \frac{\pi}{8} \text{sign}(x-x'),\\
\left[\eta_1(x),\eta_1(x')\right] &= \left[\eta_2(x),\eta_2(x')\right] = \left[\eta_1(x),\eta_3(x')\right] = \left[\eta_1(x),\eta_4(x')\right] = \left[\eta_2(x),\eta_3(x')\right] = \left[\eta_2(x),\eta_4(x')\right] = 0.
\end{align}

The total action can be divided into four contributions,
\begin{equation}
\mathcal{S}=\mathcal{S}_{12}+\mathcal{S}_{34}+\mathcal{S}_{\rm int}+\mathcal{S}_{\rm cos},
\end{equation}
where 
\begin{align}
\mathcal{S}_{12}&=\frac{1}{2\pi}\int dq\ d\omega ~ \left(\eta_1(q,\omega),\eta_2(q,\omega)\right)\left(\begin{matrix}
-v_1 q^2 & g_{12}q^2 + i \frac{q}{4} \omega\\ g_{12}q^2 + i \frac{q}{4} \omega & -v_2 q^2
\end{matrix}\right)\left(\eta_1(-q,\omega),\eta_2(-q,\omega)\right)^T,\\
\mathcal{S}_{34}&=\frac{1}{2\pi}\int dq\ d\omega ~ \left(\eta_3(q,\omega),\eta_4(q,\omega)\right)\left(\begin{matrix}
-v_3 q^2 & g_{34}q^2 + i \frac{q}{16} \omega\\ g_{34}q^2 + i \frac{q}{16} \omega & -v_4 q^2+ i \frac{q}{4} \omega
\end{matrix}\right)\left(\eta_3(-q,\omega),\eta_4(-q,\omega)\right)^T,\\
\mathcal{S}_{\rm int}&=\frac{1}{2\pi}\int dq\ d\omega ~ q^2\sum_{i,k}\frac{g_{ik}}{2}\eta_i(q,\omega) \eta_k(-q,\omega),\\
\mathcal{S}_{\rm cos}&=\frac{1}{4\pi^2 a^2}\int dx\int d\tau \left[\Lambda_1 \cos(4 \eta_1(x,\tau)) + \Lambda_2 \cos(4 \eta_2(x,\tau))\right].
\end{align}
Then, by using the following relation,
\begin{equation}
\prod_{k}\left(\int~\frac{d u_k d u_k^*}{2\pi i}\right)e^{-\sum_{ij} u_i^*A_{ij}u_j+\sum_{i}h_i^* u_i+\sum_{i}u_i^* h_i}=\frac{e^{\sum_{ij}h_i^* (A^{-1})_{ij}h_j}}{\text{Det}A},
\end{equation}
we can integrate out the bosonic fields $\eta_3$ and $\eta_4$, thus obtaining the following effective action
\begin{equation}
\mathcal{S}_{\rm eff} = \mathcal{S}_{12} + \tilde{\mathcal{S}} + \mathcal{S}_{\rm cos},
\end{equation}
where
\begin{align}
\tilde{\mathcal{S}}&=\frac{1}{2\pi}\int dq\ d\omega ~ \left(g_{13}\eta_1(q,\omega)+g_{23}\eta_2(q,\omega),g_{14}\eta_1(q,\omega)+g_{24}\eta_2(q,\omega)\right)q^2\nonumber\\&\times\left(
\begin{array}{cc}
\frac{64 (4 q v_4-i \omega )}{256 q g_{34}^2+16 i \left(q^2+1\right) \omega 
	g_{34}-q \left(256 v_3 v_4 q^2-64 i v_3 \omega  q+\omega
	^2\right)} &\frac{256 g_{34} q+16 i \omega }{256 q g_{34}^2+16 i \left(q^2+1\right) \omega 
	g_{34}-q \left(256 v_3 v_4 q^2-64 i v_3 \omega  q+\omega
	^2\right)}\\
\frac{256 g_{34} q+16 i \omega }{256 q g_{34}^2+16 i \left(q^2+1\right) \omega 
	g_{34}-q \left(256 v_3 v_4 q^2-64 i v_3 \omega  q+\omega
	^2\right)} & \frac{256 q v_3}{256 q g_{34}^2+16 i \left(q^2+1\right) \omega
	g_{34}-q \left(256 v_3 v_4 q^2-64 i v_3 \omega  q+\omega
	^2\right)} \\
\end{array}
\right) \nonumber\\&\times
\left(g_{13}\eta_1(-q,\omega)+g_{23}\eta_2(-q,\omega),g_{14}\eta_1(-q,\omega)+g_{24}\eta_2(-q,\omega)\right)^T.
\end{align}
We note that, since the Luttinger liquid parameters flow as $K_{\rho}\rightarrow 0$, $K_{\sigma}\rightarrow 1$, the coefficients in Eqs.~\eqref{eq:param1}-\eqref{eq:param2} become such that $\tilde{\mathcal{S}}\rightarrow 0$. Moreover, one also finds that $g_{12}\rightarrow 0$ and $v_1 = v_2  = 4u_{\sigma} \equiv v$. As a result, the effective action becomes
\begin{equation}
\mathcal{S}_{\rm eff} = \mathcal{S}_{12}+\mathcal{S}_{\rm cos},
\end{equation}
with 
\begin{align}
\mathcal{S}_{12}&=\frac{1}{2\pi}\int dq\ d\omega ~ \left(\eta_1(q,\omega),\eta_2(q,\omega)\right)\left(\begin{matrix}
-v q^2 & + i \frac{q}{4} \omega\\ + i \frac{q}{4} \omega & -v q^2 \end{matrix}\right)\left(\eta_1(-q,\omega),\eta_2(-q,\omega)\right)^T,
\end{align}
which corresponds to the Hamiltonian $H_{\rm DSG}$ defined in Eq.~(7) of the main text.

\section{\label{secSm:para}Low-energy limit of the $\mathbb{Z}_k$ clock model }

In this Section, the bosonized forms of the operators $\sigma_{p,i}$ ($p=1,2$), given in Eq~(11)  of the main text, are used to prove that the DSG Hamiltonian $H_{\rm DSG}$, defined  in Eq.~(7) of the main text, is also describing the low-energy limit of the $\mathbb{Z}_4$ parafermion clock chain model. Since this low-energy correspondence is valid for the general case of $\mathbb{Z}_k$ symmetry, we provide the mapping for an arbitary value of $k>1$. The complete mapping proceeds in two  steps. First, we remind the reader of the well-known mapping from the $\mathbb{Z}_k$ clock model to the $\mathbb{Z}_k$ parafermion chain~\cite{Fendley}. We emphasize that, since in this step no assumption is necessary, these two models, the $\mathbb{Z}_k$ parafermion chain and the $\mathbb{Z}_k$ clock model,  are entirely equivalent: for this reason, we denote both of the corresponding Hamiltonians with the same symbol $H_{\rm cl}$. In a second step, we introduce a  representation of $\sigma_{1,i}$ and $\sigma_{2,i}$ in terms of the bosonic fields $\eta_1$ and $\eta_2$ introduced in the main text. We prove that this bosonic representation implements the correct commutation relations for $\sigma_{1,i}$ and $\sigma_{2,i}$. Then, we exploit them to prove that the $\mathbb{Z}_k$ self-dual sine-Gordon model, $H_{\rm DSG}$, is the low-energy limit of the clock model and, therefore, of its parafermion representation.\\
\subsection{From $\mathbb{Z}_k$ clock model to $\mathbb{Z}_k$ parafermion chain}
The one-dimensional lattice Hamiltonian for the $\mathbb{Z}_k$ clock model reads~\cite{Fendley,Fateev} 
\begin{equation}
H_{\rm cl}=- \sum_{\alpha=1}^{k-1}\left[\sum_{i=1}^{n-1}J_{1,\alpha}\left(\sigma_{i}^{\dagger}\sigma_{i+1}\right)^{\alpha} +\sum_{i=1}^{n} J^*_{2,\alpha}\left(\tau_i\right)^{\alpha}\right],\label{eq:pottsSM}
\end{equation}
where $n=L/a$, with $a$ being the lattice constant, and $k>1$. The operators $\sigma_{i}$ and $\tau_i$ satisfy the following set of relations:
\begin{align}
\sigma_{i}^k = 1, \hspace{3mm} &\sigma_{i}^{k-1} = \sigma_{i}^{\dagger},\hspace{3mm}\tau_{i}^k = 1, \hspace{3mm} \tau_{i}^{k-1} = \tau_{i}^{\dagger}\label{eq:commsigma1}\\
&\sigma_{i}\tau_{i} = \omega\tau_i\sigma_{i},\label{eq:commsigma2}
\end{align}
where $\omega= e^{ \frac{2\pi i}{k}}$ and the local operators $\sigma_{i}$ and $\tau_j$ commute at different sites $i\neq j$.
In order for the Hamiltonian $H_{\rm cl}$ in Eq.~(\ref{eq:pottsSM}) to be hermitian, the coefficients must satisfy  the following relations: $J_{1,\alpha}^* = J_{1,k-\alpha}$ and $J_{2,\alpha}^* = J_{2,k-\alpha}$. 
Next, $H_{\rm cl}$ expressed in terms of the clock operators $\sigma_{i}$ and $\tau_{i}$ can be mapped onto a parafermion  representation with the help of the  operators $\chi_{p,i}$ defined as 
\begin{equation}
\chi_{1,i}=\sigma_{i}\prod_{l<i} \tau_l, ~~~ \chi_{2,i} =\omega^{(k-1)/2}  \chi_{1,i}\tau_i .
\end{equation}
The operators $\chi_{p,i}$ (for $p=1,2$) obey  $\mathbb{Z}_k$ parafermion statistics,
\begin{align}
&\chi_{p,i}^k = 1, ~\chi_{p,i}^{k-1} = \chi_{p,i}^{\dagger},\label{eq:para_stat1_SM}\\
&\chi_{p,i} \chi_{p',m} = \omega \chi_{p',m} \chi_{p,i},\,\, i< m,\,\, {\rm and}\,\,\, \chi_{1,i} \chi_{2,i} = \omega \chi_{2,i} \chi_{1,i}
\label{eq:para_stat2_SM}.
\end{align} 
Using these relations we can map the clock model defined in Eq.~(\ref{eq:pottsSM}) onto the parafermion chain~\cite{Fradkin} 
\begin{equation}
H_{\rm cl}=-\sum_{\alpha=1}^{k-1}\left[\sum_{i=1}^{n-1}J_{1,\alpha}\left(\omega^{-(k-1)/2}\chi_{2,i}^{\dagger}\chi_{1,i+1}\right)^{\alpha} +\sum_{i=1}^{n} J^*_{2,\alpha}\left(\omega^{-(k-1)/2}\chi_{1,i}^{\dagger}\chi_{2,i}\right)^{\alpha}\right].\label{seq:parafermions}
\end{equation}
This Hamiltonian can be rewritten as
\begin{align}
H_{\rm cl}=-\sum_{\alpha=1}^{k-1}\left[\sum_{j=1}^{n-1}J_{1,\alpha}\left(\omega^{-(k-1)/2}\chi_{2,j}^{\dagger}\chi_{1,j+1}\right)^{\alpha} +\sum_{j=1}^{n} J_{2,\alpha}\left(\omega^{(k-1)/2}\chi_{2,j}^{\dagger}\chi_{1,j}\right)^{\alpha}\right]\nonumber\\=-\sum_{\alpha=1}^{k-1}(-1)^{\alpha}
\left[\sum_{j=1}^{n-1}J_{1,\alpha}\left(e^{i\pi /k}\chi_{2,j}^{\dagger}\chi_{1,j+1}\right)^{\alpha}+\sum_{j=1}^{n}J_{2,\alpha}\left(e^{-i\pi /k} \chi_{2,j}^{\dagger}\chi_{1,j}\right)^{\alpha}\right],
\label{eq:para2}
\end{align}
where [with respect to Eq.~\eqref{seq:parafermions}] we rewrote the second term  using the hermitian conjugate of the same term in Eq.~\eqref{seq:parafermions}, thus obtaining the form given in the main text [see Eq.~(8)].

\subsection{$\mathbb{Z}_k$ double sine-Gordon model (DSGM)}
As preparation for the mapping of the clock model Eq.~\eqref{eq:pottsSM} onto the DSGM in the  continuum limit $a/L \ll 1$, we 
recall some of  the essential properties of the DSGM given by 
\begin{align}
H_{\rm DSG} &=\int_{0}^{L} dx \left[\frac{u K}{2\pi}   \left(\partial_x \theta(x)\right)^2  + \frac{u}{2\pi K}   \left(\partial_x \phi(x)\right)^2+ \frac{\Lambda_1}{4\pi^2a^{2}} \cos\left(\sqrt{2k}\phi(x)\right) + \frac{\Lambda_2}{4\pi^2a^{2}} \cos\left( \sqrt{2k}\theta(x)\right)\right] ,
\label{eq:DSGphitheta}
\end{align}
where we fix $K\le 1$ and where the dual bosonic fields $\phi$ and $\theta$ satisfy the commutation relation
\begin{equation}
\label{eq:thetacomm}
\left[\phi(x),\theta(x')\right] = \frac{i \pi}{2}\text{sign}(x-x')
\end{equation}
and are related to the $\eta_p$ fields as
\begin{align}
\eta_1(x)&=\sqrt{\frac{2}{k}}\phi(x),\\\eta_2(x)&=\sqrt{\frac{2}{k}}\theta(x).
\end{align}
The $\mathbb{Z}_k$ DSGM in terms of $\eta_p$ reads
\begin{align}
H_{\rm DSG} &= \int_{0}^{L} dx \left[\frac{v K}{2\pi}   \left(\partial_x \eta_2(x)\right)^2  + \frac{v}{2\pi K}   \left(\partial_x \eta_1(x)\right)^2+ \frac{\Lambda_1}{4\pi^2a^{2}} \cos\left(k \eta_1(x)\right) + \frac{\Lambda_2}{4\pi^2a^{2}} \cos\left(k \eta_2(x)\right)\right] ,
\label{eq:DSGeta}
\end{align}
where $v = 2u /k$.

Next, we list some useful relations for bosonic operators needed in the following derivation~\cite{Giamarchi,Manisha}. If  $A$ and $B$ are bosonic operators that are linear functions of bosonic creation and annihilation operators, we have
\begin{equation}
e^Ae^B = :e^{A+B}: e^{\frac{1}{2}\left\langle A^2 + 2AB + B^2\right\rangle }\, , \label{eq:expord}
\end{equation}
where $:\dots:$ stands for normal ordering and $\left\langle...\right\rangle$ denotes the bosonic ground state expectation value~\cite{Giamarchi}. 
Further we will make use of the relations~\cite{Manisha,Giamarchi}:
\begin{align}
\left\langle\left[\phi(x,\tau)- \phi(0,0)\right]^2\right\rangle &= K\ln \frac{\sqrt{x^2 + u^2\tau^2}}{a}, ~~~\left\langle\left[\theta(x,\tau)- \theta(0,0)\right]^2\right\rangle = \frac{1}{K}\ln \frac{\sqrt{x^2 + u^2\tau^2}}{a},\\
\left\langle \phi^2(0,0)\right\rangle &= -\frac{K}{2}\ln \frac{2\pi a}{L}, ~~~\left\langle \theta^2(0,0)\right\rangle = -\frac{1}{2 K}\ln  \frac{2\pi a}{L},
\end{align} 
where, again, 
$a$ is the lattice spacing and the expectation values and time evolution (in imaginary time $\tau$), $\phi(x,\tau)$ and $\theta(x,\tau)$, are governed by the kinetic term of Eq.~\eqref{eq:DSGphitheta}. These relations can be expressed in terms of the bosonic fields $\eta_p$ as
\begin{align}
\left\langle\left[\eta_1(x,\tau)- \eta_1(0,0)\right]^2\right\rangle &= \frac{2K}{k}\ln \frac{\sqrt{x^2 + u^2\tau^2}}{a}, ~~~\left\langle\left[\eta_2(x,\tau)- \eta_2(0,0)\right]^2\right\rangle = \frac{2}{Kk}\ln \frac{\sqrt{x^2 + u^2\tau^2}}{a},\\
\left\langle \eta^2_1(0,0)\right\rangle &= -\frac{K}{k}\ln \frac{2\pi a}{L},~~~\left\langle \eta^2_2(0,0)\right\rangle = -\frac{1}{kK}\ln \frac{2\pi a}{L}.
\end{align} 
In the following, we will be interested in the case $\tau = 0$ and suppress the $\tau$-argument. Note that above relations are valid for
translationally invariant systems, e.g. satisfied for periodic boundary conditions. Thus, our mapping is strictly speaking restricted to this case. However,
it is straightforward to describe boundary effects in the continuum theory by allowing for domain walls [see main text and section \ref{secSm:para}].
By using these relations, we obtain
\begin{align}
\label{eq:loga}
&e^{-\frac{1}{2}c^2\langle \left[\eta_{1}(0,0)\right]^2\rangle} = \left(\frac{2\pi a}{L}\right)^{K c^2 / 2k},\\
&e^{-\frac{1}{2}c^2\langle \left[\eta_{2}(0,0)\right]^2\rangle}= \left(\frac{2\pi a}{L}\right)^ {c^2 / \left(2k K\right)},
\end{align}
where $c$ is some real constant.
Using Eq.~\eqref{eq:expord} for $A=i c_1 \eta_p(x_j)$ and $B=i c_2 \eta_p(x_j)$, we find
\begin{equation}
e^{i c_1 \eta_1(x_j)}e^{i c_2 \eta_1(x_j)} =\, :e^{i (c_1+c_2) \eta_1(x_j)}: \left(\frac{2\pi a}{L}\right)^{K (c_1+c_2)^2/2k}, ~~~e^{i c_1 \eta_2(x_j)}e^{i c_2 \eta_2(x_j)} = \, :e^{i (c_1+c_2) \eta_2(x_j)}: \left(\frac{2\pi a}{L}\right)^{ (c_1+c_2)^2/\left(2k K\right)}, 
\end{equation}
while for $A=-i c \eta_p(x_j)$ and $B=i c \eta_p(x_{j+1})$ we get 
\begin{align}
\label{eq:normalordereta}
e^{-i c \eta_1(x_j)}e^{i c \eta_1(x_{j+1})} &=\, :e^{-i c \left[\eta_1(x_j)-\eta_1(x_{j+1})\right]}: e^{-c^2\frac{K}{k}\ln \frac{\sqrt{a^2}}{a}} = :e^{i a c \partial_j\eta_1(x_j)}: + \dots \\ 
~~~e^{-i c \eta_2(x_j)}e^{i c \eta_2(x_{j+1})} &=\, :e^{-i c \left[\eta_2(x_j)-\eta_2(x_{j+1})\right]}: e^{-c^2\frac{1}{kK}\ln \frac{\sqrt{a^2}}{a}}= :e^{i a c \partial_j\eta_2(x_j)}: + \dots
\end{align}
where $\partial_j \equiv \partial_{x_j}$. The last step, where the derivative is introduced, is valid only in the continuum limit $\frac{a}{L}\rightarrow 0$; the dots $\dots$ stand for the subleading terms that we drop by taking this limit.

For the following calculations, it is important to comment about the $\mathbb{Z}_k$ DSGM when the cosine terms are normal ordered. In this case, the 
Hamiltonian is rewritten as
\begin{align}
H_{\rm DSG} &= \int_{0}^{L} dx \left[\frac{v K}{2\pi}   \left(\partial_x \eta_2(x)\right)^2 + \frac{v}{2\pi K}   \left(\partial_x \eta_1(x)\right)^2  
+ \frac{\Lambda_1}{4\pi^2L^{2}}\left(\frac{2\pi a}{L}\right)^{\frac{k K}{2}-2} :\cos\left(k \eta_1(x)\right): + \frac{\Lambda_2}{4\pi^2L^{2}}\left(\frac{2\pi a}{L}\right)^{\frac{k}{2 K}-2} :\cos\left(k \eta_2(x)\right):\right].
\label{eq:DSGeta_normal}
\end{align}
In this representation, the cosine terms seem to be of higher-order  in $\frac{a}{L}$ compared to the kinetic term and one might argue of dropping them. However, this argument corresponds only to a tree-level RG analysis. For large enough values of $\Lambda_i$, one has to consider also the effect of the higher-order corrections to the RG flow. It has been shown that, for the $\mathbb{Z}_k$ DSGM, according to the third-order RG equations, these cosine terms are relevant and that the system flows to  multicritical fixed points, separating the $\Lambda_1$-dominated and the $\Lambda_2$-dominated phases~\cite{Katharina,Boyanovsky,Chiral_sup}. We also note that, when $K \ne 1$, the system  flows to the phase dominated by that cosine term with the largest scaling dimension (i.e. only a single cosine term remains in the Hamiltonian). 
In conclusion, in order to show consistenly that the continuum limit of the $\mathbb{Z}_k$ clock model is indeed given by the $\mathbb{Z}_k$ DSGM, one has to keep the lowest order cosine terms in $\frac{a}{L}$ in the expansion of the clock model Hamiltonian, even though they might be of higher order compared to other terms in this expansion.


\subsection{Introducing new operators}
For later purposes, it is convenient to rewrite Eq.~\eqref{eq:pottsSM} as
\begin{equation}
H_{\rm cl}=- \sum_{\alpha=1}^{k-1} \left[\sum_{i=1}^{n-1}J_{1,\alpha}\left(\sigma_{1,i}^{\dagger}\sigma_{1,i+1}\right)^{\alpha} +\sum_{i=1}^{n}J_{2,\alpha} \left(\sigma_{2,i}^{\dagger}\sigma_{2,i+1}\right)^{\alpha}\right]\label{eq:pottsSM2},
\end{equation}
where we define $\sigma_{1,i}\equiv \sigma_i$ and $\sigma_{2,i}$ is related to $\tau_i$ as $\tau_i=\sigma_{2,i}^{\dagger}\sigma_{2,i+1}$ and is assumed to possess the following  properties:
\begin{align}
\label{eq:sigma2a}
&\sigma_{2,j}^k = 1, ~~~ \sigma_{2,j}^{k-1} = \sigma_{2,j}^{\dagger}\\
&\sigma_{1,m}\sigma_{2,j}^{\dagger} =e^{i\pi \text{sign}\left(m-j+\epsilon\right)/k} \sigma_{2,j}^{\dagger}\sigma_{1,m}, \label{eq:sigma2b}\\
&\sigma_{1,m}\sigma_{2,j} =e^{-i\pi \text{sign}\left(m-j+\epsilon\right)/k} \sigma_{2,j}\sigma_{1,m},\ \label{eq:sigma2c}
\end{align}
where this expression is assumed in the limit $\epsilon \rightarrow 0^+$. With these properties, the commutation relations between $\sigma_{p,i}$, $p=1,2$, indeed implements the correct commutation relations between $\sigma_{1,m}$ and $\tau_j$, as we can easily verify:
\begin{align}
\label{eq:tausigma}
\sigma_{1,m}\tau_j &= \sigma_{1,m}\sigma_{2,j}^{\dagger}\sigma_{2,j+1}
=e^{i\pi \text{sign}\left(m-j+\epsilon\right)/k} \sigma_{2,j}^{\dagger}\sigma_{1,m}\sigma_{2,j+1} 
\nonumber \\
&= e^{i\pi \text{sign}\left(m-j+\epsilon\right)/k}e^{-i\pi \text{sign}\left(m-j-1+\epsilon\right)/k} \sigma_{2,j}^{\dagger}\sigma_{2,j+1} \sigma_{1,m} \nonumber\\&= e^{i\pi\left[ \text{sign}\left(m-j+\epsilon\right)- \text{sign}\left(m-j-1+\epsilon\right)\right]/k}\tau_j\sigma_{1,m} 
\nonumber \\
&=e^{i2\pi \delta_{m,j}/k} \tau_j\sigma_{1,m},
\end{align}
where we used the fact that, when $m$ and $j$ are integers, the following relation holds  true: 
\begin{equation}
\lim\limits_{\epsilon\rightarrow 0^+}\left[\text{sign}\left(m-j+\epsilon\right)- \text{sign}\left(m-j-1+\epsilon\right)\right]= 2\delta_{m,j}.
\end{equation}
By using Eqs.~\eqref{eq:sigma2a}-\eqref{eq:sigma2c}, we see that the desired properties of $\tau_i$ [see Eqs.~\eqref{eq:commsigma1} and \eqref{eq:commsigma2}] are also satisfied: 
\begin{equation}
\tau_j^k = \left(\sigma^{\dagger}_{2,j}\right)^k\sigma_{2,j+1}^k = 1,
\end{equation}
and
\begin{equation}
\tau_j^{k-1} = \left(\sigma^{\dagger}_{2,j}\right)^{k-1}\sigma_{2,j+1}^{k-1} = \sigma_{2,j}\sigma_{2,j+1}^{\dagger} = \sigma_{2,j+1}^{\dagger}\sigma_{2,j} = \tau^{\dagger}_j.
\end{equation}
It is instructive to note that the parafermion operators $\chi_{p,j}$ can be expressed in a simplified form in terms of the operators $\sigma_{p,j}$ by using the property that $\sigma^{\dagger}_{p,j}\sigma_{p,j}=1$. We find then
\begin{align}
\chi_{1,j} = \sigma_{1,j}\prod_{l<j}\sigma_{2,l}^{\dagger}\sigma_{2,l+1} = \sigma_{1,j}\sigma_{2,1}^{\dagger}\sigma_{2,2}\sigma^{\dagger}_{2,2}\sigma_{2,3}\dots\sigma_{2,j-2}^{\dagger}\sigma_{2,j-1}\sigma_{2,j-1}^{\dagger}\sigma_{2,j} = \sigma_{1,j}\sigma_{2,1}^{\dagger}\sigma_{2,j},
\end{align}
and
\begin{align}
\chi_{2,j} =\omega^{(k-1)/2} \chi_{1,j}\sigma_{2,j}^{\dagger}\sigma_{2,j+1} =\omega^{(k-1)/2} \sigma_{1,j}\sigma_{2,1}^{\dagger}\sigma_{2,j}\sigma_{2,j}^{\dagger}\sigma_{2,j+1} = \omega^{(k-1)/2}  \sigma_{1,j}\sigma_{2,1}^{\dagger}\sigma_{2,j+1}.
\end{align}

\subsection{Bosonic mapping and continuum limit}
Next, we introduce a bosonic representation for $\sigma_{p,j}$ in terms of $\eta_p(x_j)$:
\begin{equation}
\sigma_{p,j} = A_ke^{i \eta_p(x_j)}+B_ke^{-i (k-1) \eta_p(x_j)},\label{eq:sigma_p}\,\,\,\,\,\,  p=1,2, 
\end{equation}
where $A_k$ and $B_k$ are real coefficients which depend on $k$ (but not on $a$); below we will specify the constraints they must satisfy. The fields $\eta_p$ satisfy the following commutation relations [see Eq.~\eqref{eq:thetacomm}]:
\begin{align}
\left[\eta_p(x_m),\eta_p(x_j)\right] &= 0,\\
\left[\eta_1(x_m),\eta_2(x_j)\right] &= i\frac{\pi}{k}\text{sign}(x_m-x_j).\label{eq:commSM}
\end{align} 

We recall that the $\mathbb{Z}_k$ DSGM obeys the global symmetry $\mathbb{Z}_k$:$ ~ \eta_1\rightarrow \eta_1 + \frac{2\pi}{k}$,  $\eta_2\rightarrow \eta_2$ and its dual symmetry $\mathbb{Z}^{\rm dual}_k$:$ ~ \eta_1\rightarrow \eta_1$, $ \eta_2\rightarrow \eta_2+ \frac{2\pi}{k}$. 
Under these transformations, it follows from Eq. \eqref{eq:sigma_p} that the $\sigma_{p,j}$ fields behave as
\begin{align}
&\mathbb{Z}_k: ~ \sigma_{1,j} \rightarrow \omega\sigma_{1,j}, ~ \sigma_{2,i}\rightarrow \sigma_{2,j}\\
&\mathbb{Z}^{\rm dual}_k: ~ \sigma_{1,j}\rightarrow \sigma_{1,j} , ~ \sigma_{2,j} \rightarrow \omega \sigma_{2,j},
\end{align} 
as expected for clock model operators~\cite{Fendley}. 

Next we prove that $\sigma_{p,j}$ in the bosonic representation, Eq. \eqref{eq:sigma_p},
indeed satisfies the commutation relations given in Eqs. \eqref{eq:commsigma1}, \eqref{eq:commsigma2}, \eqref{eq:sigma2a}, \eqref{eq:sigma2b}, and \eqref{eq:sigma2c}. 
First we consider the property
\begin{align}
\sigma_{p,j}^{\dagger}=\sigma_{p,j}^{-1}, 
\end{align}
or equivalently 
\begin{align}
\sigma_{p,j}^{\dagger}\sigma_{p,j}=1.
\end{align}
By using the bosonic representation of these operators and making use of Eqs. \eqref{eq:loga} - \eqref{eq:normalordereta}, we find
\begin{align}
\sigma_{1,j}^{\dagger}\sigma_{1,j} &=\left[A_ke^{-i \eta_1(x_j)}+B_k e^{i (k-1) \eta_1(x_j)}\right]\left[A_k e^{i \eta_1(x_j)}+B_k e^{-i (k-1) \eta_1(x_j)}\right]\nonumber\\&= A_k^2+ B_k^2 + A_k B_k \left(\frac{2\pi a}{L}\right)^{Kk/2} \left[:e^{-i k \eta_1(x_j)}: + :e^{i k \eta_1(x_j)}:\right]= A_k^2+ B_k^2 + \dots,\\
\sigma_{2,j}^{\dagger}\sigma_{2,j} &=\left[A_ke^{-i \eta_2(x_j)}+B_k e^{i (k-1) \eta_2(x_j)}\right]\left[A_k e^{i \eta_2(x_j)}+B_k e^{-i (k-1) \eta_2(x_j)}\right]\nonumber\\&= A_k^2+ B_k^2 + A_k B_k \left(\frac{2\pi a}{L}\right)^{k/2K}\left[:e^{-i k \eta_2(x_j)}: + :e^{i k \eta_2(x_j)}:\right]= A_k^2+ B_k^2 + \dots,
\end{align}
where the dots stand for terms proportional to positive powers of $\frac{a}{L}$, which we will drop in the continuum limit $\frac{a}{L}\to 0$. 
By imposing that $A_k^2+ B_k^2 = 1$, we arrive in the continuum limit at 
\begin{align}
\sigma_{p,j}^{\dagger}\sigma_{p,j} = 1,
\end{align}
which indeed proves $\sigma_{p,j}^{\dagger}=\sigma_{p,j}^{-1}$. As a second relation, we show that
\begin{align}
\sigma_{p,j}^k& = 1.
\end{align}
In this case, we find similarly
\begin{align}
\sigma_{1,j}^k = k A_k^{k-1}B_k +A_k^k \left(\frac{2\pi a}{L}\right)^{kK/2}:e^{i k \eta_1(x_j)}:+ \sum_{q=2}^{k} \left(\begin{matrix} k\\q\end{matrix}\right) A_k^{k-q}B_k^q :e^{i k \left(1-q\right) \eta_1(x_j)}: \left(\frac{2\pi a}{L}\right)^{\left(1-q\right)^2 kK/2}= k A_k^{k-1}B_k +  \dots,\\
\sigma_{2,j}^k = k A_k^{k-1}B_k +A_k^k \left(\frac{2\pi a}{L}\right)^{k/2K}:e^{i k \eta_2(x_j)}:+ \sum_{q=2}^{k} \left(\begin{matrix} k\\q\end{matrix}\right) A_k^{k-q}B_k^q :e^{i k \left(1-q\right) \eta_2(x_j)}: \left(\frac{2\pi a}{L}\right)^{\left(1-q\right)^2 k/2K}= k A_k^{k-1}B_k +  \dots
\end{align}
Imposing a further condition $k A_k^{k-1}B_k = 1$, we see that the relation $\sigma_{p,j}^k = 1$ is also satisfied. By putting these results together, we have shown that
\begin{align}
\sigma_{p,j}^{\dagger}&=\sigma_{p,j}^{-1},\label{eq:comm1}\\
\sigma_{p,j}^k& = 1,\label{eq:comm2}
\end{align}
provided the coefficients $A_k$ and $B_k$ satisfy both conditions $A_k^2+ B_k^2 = 1$ and $k A_k^{k-1}B_k = 1$. This eventually proves the validity of the relations given in
Eqs. \eqref{eq:commsigma1} in the bosonic representation and in the continuum limit.

The two relations in Eqs.~\eqref{eq:comm1} and \eqref{eq:comm2} together give
\begin{equation}
\sigma_{p,j}^{k-1}=\sigma_{p,j}^\dagger.\label{eq:comm3}
\end{equation}

As a final step, we show that, in the bosonic representation, the correct commutation relations between $\sigma_{1,j}$ and $\tau_j$, Eq. \eqref{eq:tausigma}, are satisfied. We have
\begin{align}
\tau_j &=\sigma_{2,j}^{\dagger}\sigma_{2,j+1} =\left[A_ke^{-i \eta_2(x_j)}+B_k e^{i (k-1) \eta_2(x_j)}\right]\left[A_k e^{i \eta_2(x_{j+1})}+B_k e^{-i (k-1) \eta_2(x_{j+1})}\right]\nonumber\\
&= A_k^2 :e^{-i\left[\eta_2(x_j)-\eta_2(x_{j+1})\right]}: +B_k^2 :e^{i\left[(k-1)\eta_2(x_j)-(k-1)\eta_2(x_{j+1})\right]}: \nonumber \\&+A_k B_k e^{-i\left[\eta_2(x_j)-k\eta_2(x_{j})\right]}e^{i\eta_2(x_{j+1})}+A_k B_k e^{i\left[(k-1)\eta_2(x_j)-k\eta_2(x_{j})\right]}e^{-i(k-1)\eta_2(x_{j+1})} \nonumber\\& = 
A_k^2 :e^{ia\partial_j\eta_2(x_j)}: + B_k^2 :e^{-ia(k-1)\partial_j\eta_2(x_j)}: \nonumber\ \\
&+\left(\frac{2\pi a}{L}\right)^{k/2K}\left[A_k B_k :e^{ia\partial_j\eta_2(x_j)}: :e^{ik \eta_2(x_j)}:+A_k B_k :e^{-ia(k-1)\partial_j\eta_2(x_j)}:  :e^{-ik \eta_2(x_{j})}:\right]\nonumber \\&=A_k^2 :e^{ia\partial_j\eta_2(x_j)}: + B_k^2 :e^{-ia(k-1)\partial_j\eta_2(x_j)}: + \dots \label{eq:esp_tau}
\end{align}
Next, we make use of the Baker-Hausdorff-Campbell relation for two operators $A$ and $B$ whose commutator
is a c-number:
\begin{equation}
\label{eq:BCH}
e^A e^B = e^B e^A e^{\left[A,B\right]}.
\end{equation} 
Using the commutator
\begin{equation}
\left[\eta_1(x_m),\partial_j \eta_2(x_j)\right] = -\frac{i 2\pi}{k a}\delta_{m,j},
\end{equation}
which is obtained by taking the derivative with respect to $x_j$ of Eq.~\eqref{eq:commSM},
we then obtain with Eq.~\eqref{eq:BCH}
\begin{equation}
e^{i c_1 \eta_1(x_m)}e^{i c_2 a \partial_j\eta_2(x_j)} = e^{i c_2 a \partial_j\eta_2(x_j)}e^{i c_1 \eta_1(x_m)} e^{i c_1 c_2 2\pi \delta_{m,j}/k}\label{eq:exp_comm}.
\end{equation}
For the expression $\sigma_{1,m} \tau_j$, one has three different possible values of $c_1 c_2$, which gives the same exponential factor
\begin{align}
c_1 c_2 &= 1 \rightarrow e^{i 2\pi \delta_{m,j}/k},\label{eq:c1c2_1}\\
c_1 c_2 &= -(k-1) \rightarrow e^{-i(k-1) 2\pi \delta_{m,j}/k} = e^{-i2\pi \delta_{m,j}} e^{i 2\pi \delta_{m,j}/k} = e^{i 2\pi \delta_{m,j}/k},\label{eq:c1c2_k-1}\\
c_1 c_2 &= (k-1)^2 \rightarrow e^{i(k-1)^2 2\pi \delta_{m,j}/k} = e^{i2\pi(k-2) \delta_{m,j}} e^{i 2\pi \delta_{m,j}/k} = e^{i 2\pi \delta_{m,j}/k}.\label{eq:c1c2_(k-1)^2}
\end{align}
By using the definition of $\sigma_{1,m}$ (see Eq.~\eqref{eq:sigma_p} for $p=1$) and the expansion in Eq.~\eqref{eq:esp_tau} for $\tau_i$, one obtains
\begin{align}
\sigma_{1,m} \tau_j &= \left[A_k e^{i \eta_1(x_m)}+B_k e^{-i (k-1) \eta_1(x_m)}\right]\left[A_k^2 :e^{ia\partial_j\eta_2(x_j)}: + B_k^2 :e^{-ia(k-1)\partial_j\eta_2(x_j)}: + \dots\right] \nonumber\\=&A^3_k e^{i \eta_1(x_m)}:e^{ia\partial_j\eta_2(x_j)}: + A_k B^2_k e^{i \eta_1(x_m)}:e^{-ia(k-1)\partial_j\eta_2(x_j)}: \nonumber\\+& B_k A_k^2 e^{-i (k-1) \eta_1(x_m)}:e^{ia\partial_j\eta_2(x_j)}:  + B^3_k e^{-i (k-1) \eta_1(x_m)}:e^{-ia(k-1)\partial_j\eta_2(x_j)}:+\dots
\end{align}
Then, in the second line of the previous equation, we can apply Eq.~\eqref{eq:c1c2_1} to the first term, Eq.~\eqref{eq:c1c2_k-1} to the second and third terms, and Eq.~\eqref{eq:c1c2_(k-1)^2} to the last term to commute the exponentials with $\eta_1$ past the exponentials with $\eta_2$, thus obtaining
\begin{align}
\sigma_{1,m} \tau_j =&A^3_k e^{i2\pi \delta_{m,j}/k} :e^{ia\partial_j\eta_2(x_j)}:e^{i \eta_1(x_j)} + A_k B^2_k  e^{i2\pi \delta_{m,j}/k} :e^{-ia(k-1)\partial_j\eta_2(x_j)}:e^{i \eta_1(x_j)} \nonumber\\+& B_k A_k^2  e^{i2\pi \delta_{m,j}/k}:e^{ia\partial_j\eta_2(x_j)}: e^{-i (k-1) \eta_1(x_j)} + B^3_k  e^{i2\pi \delta_{m,j}/k}:e^{-ia(k-1)\partial_j\eta_2(x_j)}: e^{-i (k-1) \eta_1(x_j)}+\dots  \nonumber
\end{align}
From this, we find the following commutation relation for the bosonic representations of $\sigma_{1,m}$ and $\tau_j$ in the limit $a/L \rightarrow 0$
\begin{equation}
\sigma_{1,m}\tau_j  = e^{i2\pi \delta_{m,j}/k} \tau_j  \sigma_{1,m}+\dots \label{sigmatau}
\end{equation}
This eventually proves that Eq. \eqref{eq:tausigma} is satisfied in the bosonic representation in the continuum limit.

Next, by using these operators, we derive now the DSGM as the low-energy limit of the clock Hamiltonian. For this, we  expand $\sigma_{p,j}^{\dagger}\sigma_{p,j+1}$  in powers of $\frac{a}{L}$, 
\begin{align}
\sigma_{p,j}^{\dagger}\sigma_{p,j+1} &= \left[A_ke^{-i \eta_p(x_j)}+B_k e^{i (k-1) \eta_p(x_{j})}\right]\left[A_k e^{i \eta_p(x_{j+1})}+B_k e^{-i (k-1) \eta_p(x_{j+1})}\right]\nonumber\\&=A_k^2 :e^{ia\partial_j\eta_p(x_j)}: + B_k^2 :e^{-ia(k-1)\partial_j\eta_p(x_j)}: + 2 A_k B_k \cos\left[k \eta_p(x_j)\right]+ \dots\nonumber\\&= 1 + 2 A_k B_k \cos\left[k \eta_p(x_j)\right] + \left[A_k^2 - (k-1) B_k^2\right]  i a \partial_j \eta_p(x_j) - \left[A_k^2 + (k-1)^2 B_k^2\right]a^2 :\left[\partial_j\eta_p(x_j)\right]^2:  + \dots\nonumber\\&= 1 + 2 A_k B_k \cos\left[k \eta_p(x_j)\right] + \left[A_k^2 - (k-1) B_k^2\right]  i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j) - \left[A_k^2 + (k-1)^2 B_k^2\right]\left(\frac{a}{L}\right)^2 :\left[\partial_{\bar{j}}\eta_p(\bar{x}_j)\right]^2:  + \dots,
\end{align}
where we introduced the dimensionless variable $\bar{x}_j = \frac{x_j}{L}$ and the corresponding derivative $\partial_{\bar{j}} = L \partial_{j}$ and, in contrast with the similar derivation in Eq.~\eqref{eq:esp_tau}, we kept the lowest order cosine term in $\frac{a}{L}$, in accordance with the discussion below Eq.~\eqref{eq:DSGphitheta}. This rescaling is possible because the variable $x_j$ appears in the fields $\eta_p$ as $\frac{x_j}{L}$.
Since for $k>3$, the Hamiltonian in Eq.~\eqref{eq:pottsSM2} includes also powers of $\sigma_{p,j}^{\dagger}\sigma_{p,j+1}$, the low-energy mapping is complete only when also the powers $(\sigma_{p,j}^{\dagger}\sigma_{p,j+1})^\alpha$ are provided. In order to obtain these expansions, we use the following expression 
\begin{align}
&\left\{1 + \frac{a}{L} C_1 + \left(\frac{a}{L}\right)^2 C_2 + 2 A_k B_k \cos\left[k \eta_p(x_j)\right] + \dots \right\}^{\alpha} = \alpha  \frac{a}{L} C_1 +  \frac{\alpha}{2}\left(\frac{a}{L}\right)^2\left[C_1^2\left(\alpha-1\right) + 2 C_2\right] \nonumber\\&+  \left\{1 + 2 A_k B_k \cos\left[k \eta_p(x_j)\right]\right\}^{\alpha} + \dots,
\end{align}
where $C_1$ and $C_2$ stands for the operators contributing at first order and second order in $\frac{a}{L}$ in $\sigma_{p,j}^{\dagger}\sigma_{p,j+1}$ and where we kept the lowest order cosine term in $\frac{a}{L}$, in accordance with the discussion below Eq.~\eqref{eq:DSGphitheta}. Let us consider the cosine term for $p=1$:
\begin{align}
&\left\{1 + 2A_k B_k \cos\left[k \eta_1(x_j)\right]\right\}^{\alpha} = 1+  \sum_{n=1}^{\alpha}\left(\frac{2A_k B_k}{2}\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\left(e^{ik \eta_1(x_j)}+e^{-ik \eta_1(x_j)}\right)^{n}\nonumber\\& =1 + \sum_{n=1}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\sum_{\beta=0}^{n}\left(\begin{matrix} n\\ \beta\end{matrix}\right) :e^{i \left(n-2\beta\right) k \eta_{1}(x_j)}: \left(\frac{2\pi a}{L}\right)^{\frac{K k}{2}\left(n-2\beta\right)^2} \nonumber\\& =1 + \sum_{n=1}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\sum_{\beta=0}^{n}\left(\begin{matrix} n\\ \beta\end{matrix}\right) \frac{1}{2}\left[:e^{i \left(n-2\beta\right) k \eta_{1}(x_j)}: + :e^{-i \left(n-2\beta\right) k \eta_{1}(x_j)}:\right] \left(\frac{2\pi a}{L}\right)^{\frac{K k}{2}\left(n-2\beta\right)^2} \nonumber\\& =1 + \frac{\left(\frac{2\pi a}{L}\right)^{\frac{K k}{2}} }{2}\left[:e^{i  k \eta_{1}(x_j)}: + :e^{-i  k \eta_{1}(x_j)}:\right]\sum_{n=1, \text{odd } n}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\left[\left(\begin{matrix} n\\ 
\frac{n-1}{2}\end{matrix}\right) + \left(\begin{matrix} n\\ 
\frac{n+1}{2}\end{matrix}\right)\right]  + \dots \nonumber\\ &= 1 + 2\cos\left[k \eta_1(x_j)\right] \sum_{n=1, \text{odd } n}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\left(\begin{matrix} n\\ 
\frac{n+1}{2}\end{matrix}\right)  + \dots, 
\end{align}
where the sum over $n$ is performed only for odd integers and we kept the lowest order cosine term in $\frac{a}{L}$. In the case $p=2$, one has to replace $K\rightarrow \frac{1}{K}$ in the intermediate steps, but the final result is unchanged. 

In conclusion, we have 
\begin{align}
\left(\sigma_{p,j}^{\dagger}\sigma_{p,j+1}\right)^{\alpha} & = \left\{1 + 2 A_k B_k \cos\left[k \eta_p(x_j)\right] + \left[A_k^2 - (k-1) B_k^2\right]  i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j) - \left[A_k^2 + (k-1)^2 B_k^2\right]\left(\frac{a}{L}\right)^2 :\left[\partial_{\bar{j}}\eta_p(\bar{x}_j)\right]^2: + \dots \right\}^{\alpha}  \nonumber\\& = 2\cos\left[k \eta_p(x_j)\right] \sum_{n=1, \text{odd } n}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\left(\begin{matrix} n\\ 
\frac{n+1}{2}\end{matrix}\right)   \nonumber\\  &+\alpha \left[A_k^2 - (k-1) B_k^2\right]  i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j) + \frac{\alpha(\alpha-1)}{2} \left\{\left[A_k^2 - (k-1) B_k^2\right]  i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j)\right\}^2\nonumber\\&-\alpha\left[A_k^2 + (k-1)^2 B_k^2\right]\left(\frac{a}{L}\right)^2 :\left[\partial_{\bar{j}}\eta_p(\bar{x}_j)\right]^2:  + \dots,\label{eq:expansion_sigma}
\end{align}
where the dots $\dots$ stand for higher powers of $\frac{a}{L}$ and constant terms. One can write the expansion in Eq.~\eqref{eq:expansion_sigma} in a more compact form as
\begin{align}
\left(\sigma_{p,j}^{\dagger}\sigma_{p,j+1}\right)^{\alpha} = -f_{k}^{(\alpha)}  \cos\left[k \eta_p(x_j)\right]  +\alpha \left[A_k^2 - (k-1) B_k^2\right]  i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j) -d_{k}^{(\alpha)} \left(\frac{a}{L}\right)^2 :\left[\partial_{\bar{j}}\eta_p(\bar{x}_j)\right]^2:  + \dots,
\end{align}
where the term $ \left\{ i \frac{a}{L} \partial_{\bar{j}} \eta_p(\bar{x}_j)\right\}^2$ has been written in normal ordered form as $-\left(\frac{a}{L}\right)^2 :\left[\partial_{\bar{j}}\eta_p(\bar{x}_j)\right]^2:$ by adding some unimportant constants absorbed in the remainder $\dots$ Here,
\begin{align}
d_{k}^{(\alpha)} &= \frac{\alpha}{2}\left[2\left(A_k^2 + (k-1)^2 B_k^2\right) + \left(A_k^2-(k-1) B_k^2\right)^2(\alpha-1)\right],\\
f_{k}^{(\alpha)} &=-2 \sum_{n=1, \text{odd } n}^{\alpha}\left(A_k B_k\right)^n\left(\begin{matrix}\alpha \\ n\end{matrix}\right)\left(\begin{matrix} n\\ 
\frac{n+1}{2}\end{matrix}\right) ,
\end{align}
with $d_{k}^{(\alpha)}  = d_{k}^{(k-\alpha)} $ and $f_{k}^{(\alpha)}  = f_{k}^{(k-\alpha)}  $, since, due to the properties $\sigma^{k-1}_{p,i} = \sigma^{\dagger}_{p,i}$, one has that   $(\sigma_{p,j}^{\dagger}\sigma_{p,j+1})^\alpha = (\sigma_{p,j}^{\dagger}\sigma_{p,j+1})^{k-\alpha}$. Then,  taking the continuum limit of Eq.~\eqref{eq:pottsSM} (with $\sum_{j} \rightarrow \int \frac{dx}{a}$), we find 
\begin{align}
H_k = &\int_{0}^{L} dx ~ \left\{\frac{v_1}{2\pi}:\left[\partial_x\eta_1(x)\right]^2: + \frac{v_2}{2\pi}:\left[\partial_x\eta_2(x)\right]^2:\right\} + \frac{\Lambda_1}{4\pi^2a^{2}} \int_{0}^{L} dx ~ \cos\left[k\eta_1(x)\right] +  \frac{\Lambda_2}{4\pi^2a^{2}} \int_{0}^{L} dx ~ \cos\left[k\eta_2(x)\right].
\end{align}
We note that the integral $\int_{0}^{L} dx$ gives zero for the term proportional to $\partial_{\bar{j}} \eta_p(\bar{x}_j)$. Here,
\begin{align}
{v_1} / {(2\pi a)}&=\sum_{\alpha=1}^{k-1} d^{(\alpha)}_k J_{1,\alpha},~~~
{v_2} / {(2\pi a)}=\sum_{\alpha=1}^{k-1} d^{(\alpha)}_k J_{2,\alpha},\\
{\Lambda}_1 / {(4\pi^2 a)}&=\sum_{\alpha=1}^{k-1} f^{(\alpha)}_k J_{1,\alpha},~~~
{\Lambda_2} / {(4\pi^2 a)}=\sum_{\alpha=1}^{k-1} f^{(\alpha)}_k J_{2,\alpha}.
\end{align} 
Starting from $v_1$ and $v_2$ one can also find $K$ and $v$ as $K=\sqrt{v_2/v_1}$ and $v=\sqrt{v_1 v_2}$. We note that in our case we have $K=1$ and $\Lambda_1 = \Lambda_2$.

It is instructive to write down a few special cases for $k$. For $k=2$, we find
\begin{align}
v_1 / (2\pi a)&=  J_{1,1}, ~~~ v_2 / (2\pi a)= J_{2,1},\\{\Lambda}_1 / {(4\pi^2 a)} &= - J_{1,1}, ~~~ {\Lambda}_2 / {(4\pi^2 a)} = - J_{2,1}.
\end{align}

For $k=3$ we find two solutions for $A_k, B_k$. First, we get $A_3 = \sqrt{\frac{1}{3} \left(1+2 \cos \left(\frac{2 \pi }{9}\right)\right)}$ and $B_3 = \frac{\sqrt{2 \sin \left(\frac{\pi }{9}\right)}}{\sqrt[4]{3}}$, and
\begin{align}
v_1 / (2\pi a) &=  \frac{2}{3} \left(1+8 \sqrt{3} \sin \left(\frac{\pi }{9}\right)+2 \cos \left(\frac{2 \pi
}{9}\right)\right) \text{Re}\left[J_{1,1}\right], ~~~ v_2 / (2\pi a) = \frac{2}{3} \left(1+8 \sqrt{3} \sin \left(\frac{\pi }{9}\right)+2 \cos \left(\frac{2 \pi
}{9}\right)\right)\text{Re}\left[J_{2,1}\right],\\{\Lambda}_1 / {(4\pi^2 a)} &= - 2 \frac{\sqrt{2 \sin \left(\frac{\pi }{9}\right) \left(1+2 \cos \left(\frac{2 \pi
		}{9}\right)\right)}}{3^{3/4}}\text{Re}\left[J_{1,1}\right], ~~~ {\Lambda}_2 / {(4\pi^2 a)} = - 2\frac{\sqrt{2 \sin \left(\frac{\pi }{9}\right) \left(1+2 \cos \left(\frac{2 \pi
		}{9}\right)\right)}}{3^{3/4}}\text{Re}\left[J_{2,1}\right].
\end{align}
Second, we get
$A_3 =\frac{1}{2} \sqrt{\frac{1}{3} \left(4+4 \sqrt{3} \sin \left(\frac{2 \pi }{9}\right)-4
	\cos \left(\frac{2 \pi }{9}\right)\right)}$ and $B_3 = \sqrt{1+\frac{1}{12} \left(-4-4 \sqrt{3} \sin \left(\frac{2 \pi }{9}\right)+4 \cos
	\left(\frac{2 \pi }{9}\right)\right)}$,
\begin{align}
v_1 / (2\pi a) &=  \left[3-\sqrt{3} \sin \left(\frac{2 \pi }{9}\right)+\cos \left(\frac{2 \pi }{9}\right) \right]\text{Re}\left[J_{1,1}\right], ~~~ v_2 / (2\pi a) = \left[3-\sqrt{3} \sin \left(\frac{2 \pi }{9}\right)+\cos \left(\frac{2 \pi }{9}\right)\right]\text{Re}\left[J_{2,1}\right],\\{\Lambda}_1 / {(4\pi^2 a)} &=- \frac{\sqrt{-\sin \left(\frac{\pi }{9}\right)+\sin \left(\frac{2 \pi }{9}\right)+\cos
		\left(\frac{\pi }{18}\right)}}{3^{3/4}}\text{Re}\left[J_{1,1}\right], ~~~ {\Lambda}_2 / {(4\pi^2 a)} = \frac{\sqrt{-\sin \left(\frac{\pi }{9}\right)+\sin \left(\frac{2 \pi }{9}\right)+\cos
		\left(\frac{\pi }{18}\right)}}{3^{3/4}}\text{Re}\left[J_{2,1}\right].
\end{align}

For $k=4$, we find again two solutions. First, $A_4=B_4 = \frac{1}{\sqrt{2}}$, and
\begin{align}
v_1 / (2\pi a) &=  10\text{Re}\left[J_{1,1}\right] + 11 J_{1,2}, ~~~ v_2 / a = 10\text{Re}\left[J_{2,1}\right] + 11 J_{2,2},\\{\Lambda}_1 / {(4\pi^2 a)} &= - \text{Re}\left[J_{1,1}\right] -2 J_{1,2}, ~~~ {\Lambda}_2 / {(4\pi^2 a)} = - \text{Re}\left[J_{2,1}\right] - 2 J_{2,2}.
\end{align}
Second, we get $A_4 = \frac{1}{\sqrt{\frac{6}{1+\sqrt[3]{19-3 \sqrt{33}}+\sqrt[3]{19+3 \sqrt{33}}}}}$ and $B_4 = \sqrt{1+\frac{1}{6} \left(-1-\sqrt[3]{19-3 \sqrt{33}}-\sqrt[3]{19+3 \sqrt{33}}\right)}$, and
\begin{align}
v_1 / (2\pi a) &=  \left[\frac{1}{3} \left(23-4 \sqrt[3]{19-3 \sqrt{33}}-4 \sqrt[3]{19+3 \sqrt{33}}\right)\right]\text{Re}\left[J_{1,1}\right] + \nonumber\\& + \left[\frac{1}{9} \left(219-52 \sqrt[3]{19-3 \sqrt{33}}+4 \left(19-3 \sqrt{33}\right)^{2/3}-52
\sqrt[3]{19+3 \sqrt{33}}+4 \left(19+3 \sqrt{33}\right)^{2/3}\right)\right]  J_{1,2},\\ v_2 / (2\pi a) &=  \left[\frac{1}{3} \left(23-4 \sqrt[3]{19-3 \sqrt{33}}-4 \sqrt[3]{19+3 \sqrt{33}}\right)\right]\text{Re}\left[J_{2,1}\right]\nonumber\\& + \left[\frac{1}{9} \left(219-52 \sqrt[3]{19-3 \sqrt{33}}+4 \left(19-3 \sqrt{33}\right)^{2/3}-52
\sqrt[3]{19+3 \sqrt{33}}+4 \left(19+3 \sqrt{33}\right)^{2/3}\right)\right] J_{2,2},\\{\Lambda}_1 / {(4\pi^2 a)} &= -\frac{1}{3} \sqrt{-3+4 \sqrt[3]{19-3 \sqrt{33}}-\left(19-3 \sqrt{33}\right)^{2/3}+4
	\sqrt[3]{19+3 \sqrt{33}}-\left(19+3 \sqrt{33}\right)^{2/3}} \left\{\text{Re}\left[J_{1,1}\right] + 2J_{1,2}\right\},\\ {\Lambda}_2 / {(4\pi^2 a)} &= - \frac{1}{3} \sqrt{-3+4 \sqrt[3]{19-3 \sqrt{33}}-\left(19-3 \sqrt{33}\right)^{2/3}+4
	\sqrt[3]{19+3 \sqrt{33}}-\left(19+3 \sqrt{33}\right)^{2/3}}\left\{\text{Re}\left[J_{2,1}\right] + 2J_{2,2}\right\}.
\end{align}


