
\section{More on the one dimensional Hamiltonians} \label{app: 1d models}

In this section, we show several properties of the 1d Hamiltonian mentioned in the main text. We will show that (i) the only ground states of these Hamiltonians are the ones given in the main text, (ii) all these Hamiltonians are gapped, (iii) the ground states are nontrivial symmetry-protected topological (SPT) states under the remaining $\mathbb{Z}_n^x\times\mathbb{Z}_n^z$ symmetry, and (iv) under the remaining $\mathbb{Z}_n^x\times\mathbb{Z}_n^z\times\sigma$ symmetry they are beyond the conventional group-cohomology-based classification developed in Refs. \cite{Jiang2016, Thorngren2016}. Note that the first three properties apply to the Hamiltonians defined for any $n\geqslant 2$, while the last one applies to the case with $n\geqslant 3$. In the main text, we have shown that the relevant $G_0\rightarrow G$ symmetry breaking orders are entanglement-enabled when $n\geqslant 3$. For $n=2$, the $G_0\rightarrow G$ symmetry breaking order can be realized by local product states, {i.e., } $\prod_{j}(|00\rangle_j+(-1)^j|11\rangle_j)$ with $j$ the site index, so it is not entanglement-enabled.

Here we repeat the Hamiltonian for reader's convenience:
\begin{equation}\label{eq:1ddimerH_SM}
    H_L=\sum_{j=1}^{L-1} (I-P_{j,j+1}) + r(I-P_{L,1})~,
\end{equation}
where $P_{j,j+1}$ is the projector onto the linear space 
$\mathcal{P}_{j,j+1}=\text{span}\{ |D^{(\alpha)}_{a,b}\rangle ,|D^{(\beta)}_{a,b}\rangle;  a,b=1 \dots n \}$ and
\begin{eqnarray}
|D^{(\alpha)}_{a,b}\rangle &=&\sum_{d=1}^n|\alpha_j\!=\!d,\beta_j\!=\!a,\alpha_{j+1}\!=\!d,\beta_{j+1}\!=\!b\rangle \notag \\
|D^{(\beta)}_{a,b}\rangle &=&\sum_{d=1}^n|\alpha_j\!=\!a,\beta_j\!=\!d,\alpha_{j+1}\!=\!b,\beta_{j+1}\!=\!d \rangle \notag
\end{eqnarray}
are the maximum entangled states on $\alpha$ or $\beta$ degrees of freedom. We may refer to the space spanned by these states the ``dimer" subspace of the $\alpha$ or $\beta$ degrees of freedom, respectively.
In particular, we will consider the Hamiltonian with a open boundary condition ($r=0$) or a periodic boundary condition ($r=1$).

It is sometimes useful to express the projector in terms of the operators~\cite{Jiang2019}. 
Consider the projectors projecting the states into the dimer subspace on $\alpha$ or $\beta$ degrees of freedom $P_{j,j+1}^{(\alpha/\beta)}$, we have $P_{j,j+1}^{(\alpha/\beta)}=P^x_{(\alpha/\beta)}P^z_{(\alpha/\beta)}$, where
\begin{eqnarray}
    P^x_{(\alpha)}=\frac{1}{n} \sum_{d=0}^{n-1} (\mu_j^x \mu_{j+1}^x)^d \quad \text{and} \quad  P^z_{(\alpha)}=\frac{1}{n} \sum_{d=0}^{n-1} (\mu_j^z)^{d} (\mu_{j+1}^z)^{-d}~,
\end{eqnarray}
while 
\begin{eqnarray}
    P^x_{(\beta)}=\frac{1}{n} \sum_{d=0}^{n-1} (\nu_j^x \nu_{j+1}^x)^d \quad \text{and} \quad  P^z_{(\beta)}=\frac{1}{n} \sum_{d=0}^{n-1} (\nu_j^z)^{d} (\nu_{j+1}^z)^{-d}~.
\end{eqnarray}
We then have $P_{j,j+1}=P^{(\alpha)}_{j,j+1}+P^{(\beta)}_{j,j+1}-P^{(\alpha)}_{j,j+1}P^{(\beta)}_{j,j+1}$.


\subsection{Uniqueness of the ground states}

First we show the uniqueness of the ground states. 
Consider the following states depicted in Fig~, 
\begin{eqnarray} \label{eq: ground states}
|\psi^{(A)}_{a,b}\rangle &=& \sum_{\{\alpha_j,\beta_j\}}\psi^{(A)}_{\{\alpha_j\beta_j\}}(a,b)|\alpha_1\beta_1\dots \alpha_L\beta_L\rangle ~,\notag \\
|\psi^{(B)}_{a,b}\rangle &=& \sum_{\{\alpha_j,\beta_j\}}\psi^{(B)}_{\{\alpha_j\beta_j\}}(a,b)|\alpha_1\beta_1\dots \alpha_L\beta_L\rangle ~, \quad a,b=1 \dots n
\end{eqnarray}
where 
\begin{eqnarray}
    \psi^{(A)}_{\{\alpha_j,\beta_j\}}(a,b) &=& \delta_{\alpha_1a}\left(\prod_{j=1}^{L/2-1}\delta_{\beta_{2j-1}\beta_{2j} } \delta_{\alpha_{2j}\alpha_{2j+1} }\right)\delta_{\beta_{L-1}\beta_L}\delta_{\alpha_Lb}~, \notag \\
    \psi^{(B)}_{\{\alpha_j,\beta_j\}}(a,b) &=& \delta_{\beta_1a}\left(\prod_{j=1}^{L/2-1}\delta_{\alpha_{2j-1}\alpha_{2j} }\delta_{\beta_{2j}\beta_{2j+1} }\right)\delta_{\alpha_{L-1}\alpha_L} \delta_{\beta_L b} ~,\quad \text{if $L$ is even},
\end{eqnarray}
 and
\begin{eqnarray}
    \psi^{(A)}_{\{\alpha_j,\beta_j\}}(a,b) &=& \delta_{\alpha_1a}\left(\prod_{j=1}^{(L-1)/2}\delta_{\beta_{2j-1}\beta_{2j} } \delta_{\alpha_{2j}\alpha_{2j+1} }\right)\delta_{\beta_Lb}~, \notag \\
    \psi^{(B)}_{\{\alpha_j,\beta_j\}}(a,b) &=& \delta_{\beta_1a}\left(\prod_{j=1}^{(L-1)/2}\delta_{\alpha_{2j-1}\alpha_{2j} } \delta_{\beta_{2j}\beta_{2j+1} }\right)\delta_{\alpha_L b} ~, \quad \text{if $L$ is odd}.
\end{eqnarray}

Since the Hamiltonian is of the form of a sum over projectors, the eigenvalues of $H_L$ is bounded below by $0$.
For the open boundary condition, it is easy to verify $H_L|\psi^{(A)}_{a,b}\rangle=H_L|\psi^{(B)}_{a,b}\rangle=0$ for all $a,b=1\dots n$.
We therefore would like to show that $\mathcal{G}_L \equiv \text{span}\{ |\psi^{(A)}_{a,b}\rangle, |\psi^{(B)}_{a,b}\rangle, a,b=1 \dots n \}=\text{Ker}{H_L}$.
Note that these $2n^2$ states are linearly independent (but not orthogonal) when $L \geqslant 3$, so $\text{dim}(\mathcal{G}_L)=2n^2$ when $L\geqslant 3$. If $L=2$, $\text{dim}(\mathcal{G}_L)=2n^2-1$.


To show $\text{Ker}{H_L}=\mathcal{G}_L$ for any $L \geqslant 2$, we use mathematical induction.
Assume that for a 1d lattice with $L$ sites (onsite Hilbert space dimension $n^2$) and the corresponding Hilbert space $\mathcal{H}_L$, the ground state space $\text{Ker}(H_L) = \mathcal{G}_L \subset \mathcal{H}_L$. 
Now we would like to add one more site and to find $\text{Ker}(H_{L+1}) = \text{Ker}(H_L)\cap \text{Ker}(I-P_{L,L+1})$, where both $\text{Ker}(H_L)$ and $\text{Ker}(I-P_{L,L+1})$ should be understood as a linear subspace in $\mathcal{H}_{L+1}$. 
(This means $\text{dim}(\text{Ker}(H_L))=2n^2 \times n^2$ as a linear subspace of $\mathcal{H}_{L+1}$.)

To find the intersection of the two linear subspaces, assume $|\psi \rangle \in \text{Ker}(H_L) \subset \mathcal{H}_{L+1}$.
Let us first consider the case where $L$ is odd.
We can write
\begin{eqnarray}
    |\psi\rangle&= \sum_{\{\alpha_j\beta_j\}}\left( A_{\alpha_1\beta_L\alpha_{L+1}\beta_{L+1}}\chi^{(A)}_{\beta_1\alpha_2 \dots \alpha_L}+ B_{\beta_1 \alpha_{L} \alpha_{L+1}\beta_{L+1}}\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}\right)|\alpha_1\beta_1 \dots \alpha_{L+1}\beta_{L+1}\rangle~,
\end{eqnarray}
where $\chi^{(A)}_{\beta_1\alpha_2 \dots \alpha_L}=\prod_{j=1}^{(L-1)/2}\delta_{\beta_{2j-1}\beta_{2j} } \delta_{\alpha_{2j}\alpha_{2j+1} }$ and $\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}=\prod_{j=1}^{(L-1)/2}\delta_{\alpha_{2j-1}\alpha_{2j} } \delta_{\beta_{2j}\beta_{2j+1} }$ (the dimer part of the wavefunction), and there are $2n^2 \times n^2$ coefficients $A_{\alpha_1\beta_L\alpha_{L+1}\beta_{L+1}}$ and $B_{\beta_1 \alpha_{L} \alpha_{L+1}\beta_{L+1}}$.
Requiring $|\psi\rangle \in \text{Ker}(I-P_{L,L+1})$ as well, we can also express
\begin{eqnarray}
    |\psi\rangle&= \sum_{\{\alpha_j\beta_j\}}\left( X_{\alpha_1\beta_1\dots\alpha_L\alpha_{L+1}}\delta_{\beta_L\beta_{L+1}} +Y_{\alpha_1\beta_1\dots\beta_L\beta_{L+1}}\delta_{\alpha_L\alpha_{L+1}}\right)|\alpha_1\beta_1 \dots \alpha_{L+1}\beta_{L+1}\rangle~,
\end{eqnarray}
where $X_{\alpha_1\beta_1\dots\alpha_L\alpha_{L+1}}$ and $Y_{\alpha_1\beta_1\dots\beta_L\beta_{L+1}}$ are the coefficients of the linear combination of the vectors in $\text{Ker}(I-P_{L,L+1})$.

Equating the above two expressions, we obtain $n^{2(L+1)}$ linear equations ($\alpha_j,\beta_j =1 \dots n$)
\begin{eqnarray}
A_{\alpha_1\beta_L\alpha_{L+1}\beta_{L+1}}\chi^{(A)}_{\beta_1\alpha_2 \dots \alpha_L}+ B_{\beta_1 \alpha_{L} \alpha_{L+1}\beta_{L+1}}\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}=X_{\alpha_1\beta_1\dots\alpha_L\alpha_{L+1}}\delta_{\beta_L\beta_{L+1}} + Y_{\alpha_1\beta_1\dots\beta_L\beta_{L+1}}\delta_{\alpha_L\alpha_{L+1}}~, \notag
\end{eqnarray} 
which pose constraints on $A_{\alpha_1\beta_L\alpha_{L+1}\beta_{L+1}}$ and $B_{\beta_1 \alpha_{L} \alpha_{L+1}\beta_{L+1}}$ (or $X_{\alpha_1\beta_1\dots\alpha_L\alpha_{L+1}}$ and $Y_{\alpha_1\beta_1\dots\beta_L\beta_{L+1}}$ conversely).
For convenience, we will use the notation $\bar{\alpha}_j$ to denote the numbers in $\{\bar{\alpha}_j = 1 \dots n, \bar{\alpha}_j \neq \alpha_j\}$ and similarly for $\bar{\beta}_j$.


First, considering the subset of the equations where $\beta_{2j}=\beta_{2j-1}$, $\alpha_{2j+1}=\alpha_{2j}$ for $j=1 \dots (L-1)/2$ (so that $\chi^{(A)}=1$), $\beta_{L} \neq \beta_{L+1}$, $\alpha_{L} \neq \alpha_{L+1}$ (so that the right-hand side is zero) and $\beta_{L-1} \neq \beta_L$ (so taht $\chi^{(B)}=0$), we have
\begin{eqnarray}
A_{\alpha_1\beta_L\alpha_{L+1}\bar{\beta}_L}=0~. \notag
\end{eqnarray} 
Similarly, the subset of the equations where $\alpha_{2j}=\alpha_{2j-1}$, $\beta_{2j+1}=\beta_{2j}$ for $j=1 \dots (L-1)/2$ (so that $\chi^{(B)}=1$), $\beta_{L} \neq \beta_{L+1}$, $\alpha_{L} \neq \alpha_{L+1}$ (so that the right-hand side is zero) and $\alpha_{L-1} \neq \alpha_L$ (so that $\chi^{(A)}=0$) give us 
\begin{eqnarray}
B_{\beta_1\alpha_L\bar{\alpha}_{L}\beta_{L+1}}=0~. \notag
\end{eqnarray}

Next, we consider again the equations where $\beta_{2j}=\beta_{2j-1}$, $\alpha_{2j+1}=\alpha_{2j}$ for $j=1 \dots (L-1)/2$ (so that $\chi^{(A)}=1$), but now with $\beta_{L} = \beta_{L+1}$ and $\alpha_{L} \neq \alpha_{L+1}$ (so that $B=0$), giving us
\begin{eqnarray}
A_{\alpha_1\beta_L\alpha_{L+1}\beta_L}= X_{\alpha_1\beta_1\dots \bar{\alpha}_{L+1}\alpha_{L+1}}~. \notag
\end{eqnarray} 
Note that on the right-hand side, it is independent of $\beta_L$ (but depending on $\alpha_1$ and $\alpha_{L+1}$). 
We therefore conclude that $A_{\alpha_1\beta_L\alpha_{L+1}\beta_{L+1}} = a_{\alpha_1 \alpha_{L+1}} \delta_{\beta_L\beta_{L+1}}$.
Similarly, we consider the equations where $\alpha_{2j}=\alpha_{2j-1}$, $\beta_{2j+1}=\beta_{2j}$ for $j=1 \dots (L-1)/2$ (so that $\chi^{(B)}=1$), but now with $\beta_{L} \neq \beta_{L+1}$ (so that $A=0$) and $\alpha_{L} = \alpha_{L+1}$.
We have
\begin{eqnarray}
B_{\beta_1\alpha_L\alpha_L\beta_{L+1}}= Y_{\alpha_1\beta_1\dots \bar{\beta}_{L+1}\beta_{L+1}}~. \notag
\end{eqnarray} 
Again, since it is independent of $\alpha_L$ on the right-hand side, we conclude that $B_{\beta_1\alpha_L\alpha_{L+1}\beta_{L+1}}=b_{\beta_1\beta_{L+1}} \delta_{\alpha_L\alpha_{L+1}}$.
This shows that $\mathcal{G}_{L+1} \supseteq \text{Ker}(H_{L+1})$.
Since we already know $\mathcal{G}_{L+1} \subseteq \text{Ker}(H_{L+1})$, we indeed have $\mathcal{G}_{L+1} = \text{Ker}(H_{L+1})$.





If $L$ is even, a similar argument can be made with some modification. 
In particular, the set of linear equations constraining the coefficients $A$, $B$, $X$ and $Y$ are now
\begin{eqnarray}
A_{\alpha_1\alpha_L\alpha_{L+1}\beta_{L+1}}\chi^{(A)}_{\beta_1\alpha_2 \dots \beta_{L-1}\beta_L}+ B_{\beta_1 \beta_{L} \alpha_{L+1}\beta_{L+1}}\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \alpha_{L}\beta_L}=X_{\alpha_1\beta_1\dots\alpha_L\alpha_{L+1}}\delta_{\beta_L\beta_{L+1}} + Y_{\alpha_1\beta_1\dots\beta_L\beta_{L+1}}\delta_{\alpha_L\alpha_{L+1}}~, \notag
\end{eqnarray} 
where
$\chi^{(A)}_{\beta_1\alpha_2 \dots \beta_{L-1}\beta_L}=\prod_{j=1}^{L/2-1}\delta_{\beta_{2j-1}\beta_{2j} } \delta_{\alpha_{2j}\alpha_{2j+1} }\delta_{\beta_{L-1}\beta_L}$ and $\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \alpha_{L}\beta_L}=\prod_{j=1}^{L/2-1}\delta_{\alpha_{2j-1}\alpha_{2j} }\delta_{\beta_{2j}\beta_{2j+1} }\delta_{\alpha_{L-1}\alpha_L}$.

Again, considering the subset of equations where $\beta_{2j-1}=\beta_{2j}$ for $j=1 \dots L/2$,  $\alpha_{2j}=\alpha_{2j+1}$  for $j=1 \dots L/2-1$, $\beta_L\neq\beta_{L+1}$, $\alpha_L\neq\alpha_{L+1}$ and $\alpha_{L-1}\neq\alpha_{L}$, we have $A_{\alpha_1\alpha_L\bar{\alpha}_{L}\beta_{L+1}}=0$;
while considering the subset of equations where $\alpha_{2j-1}=\alpha_{2j}$ for $j=1 \dots L/2$,  $\beta_{2j}=\beta_{2j+1}$  for $j=1 \dots L/2-1$, $\alpha_L\neq\alpha_{L+1}$, $\beta_L\neq\beta_{L+1}$ and $\beta_{L-1}\neq\beta_{L}$, we have $B_{\beta_1\beta_L\alpha_{L}\bar{\beta}_{L+1}}=0$.

Now the set of equations where $\beta_{2j-1}=\beta_{2j}$ for $j=1 \dots L/2$,  $\alpha_{2j}=\alpha_{2j+1}$  for $j=1 \dots L/2-1$, $\beta_L\neq\beta_{L+1}$, $\alpha_L=\alpha_{L+1}$ give us
$A_{\alpha_1\alpha_L\alpha_{L}\beta_{L+1}}=Y_{\alpha_1\dots \bar{\beta}_{L+1}\beta_{L+1}}$, resulting in $A_{\alpha_1\alpha_L\alpha_{L+1}\beta_{L+1}}=a_{\alpha_1\beta_{L+1}}\delta_{\alpha_L\alpha_{L+1}}$.
Similarly, the set of equations where $\alpha_{2j-1}=\alpha_{2j}$ for $j=1 \dots L/2$,  $\beta_{2j}=\beta_{2j+1}$  for $j=1 \dots L/2-1$, $\beta_L=\beta_{L+1}$, $\alpha_L\neq\alpha_{L+1}$ give us
$B_{\beta_1\beta_L\alpha_{L+1}\beta_{L}}=X_{\alpha_1\dots \bar{\alpha}_{L+1}\alpha_{L+1}}$, resulting in $B_{\beta_1\beta_L\alpha_{L+1}\beta_{L+1}}=b_{\beta_1\alpha_{L+1}}\delta_{\beta_L\beta_{L+1}}$.
We therefore can conclude $\mathcal{G}_{L+1} = \text{Ker}(H_{L+1})$ if $L$ is even.
The mathematical induction is therefore establish since $\mathcal{G}_{2} = \text{Ker}(H_2)=\text{Ker}(I-P_{1,2})$.

Finally, we show that if $L \geq 4$ is even and $r=1$ (periodic boundary condition), the ground state space is spanned by the states depicted in the main text.
We would like to find $\text{Ker}[H_L(r=1)]=\text{Ker}[H_L(r=0)] \cap \text{Ker}(I-P_{L,1}) \subset \mathcal{H}_L$.

By expressing $|\psi \rangle$ as 
\begin{eqnarray}
    |\psi\rangle &=& \sum_{\{\alpha_j\beta_j\}}\left( A_{\alpha_1\alpha_{L}}\chi^{(A)}_{\beta_1\alpha_2 \dots \alpha_L}+ B_{\beta_1 \beta_{L}}\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}\right)|\alpha_1\beta_1 \dots \alpha_{L+1}\beta_{L+1}\rangle \notag \\
    &=&\sum_{\{\alpha_j\beta_j\}}\left( X_{\alpha_1\alpha_2\beta_2\dots\alpha_L}\delta_{\beta_L\beta_{1}} +Y_{\beta_1\alpha_2\dots\beta_L}\delta_{\alpha_L\alpha_{1}}\right)|\alpha_1\beta_1 \dots \alpha_{L+1}\beta_{L+1}\rangle~,
\end{eqnarray}
where $\chi^{(A)}_{\beta_1\alpha_2 \dots \beta_{L-1}\beta_L}=\prod_{j=1}^{L/2-1}\delta_{\beta_{2j-1}\beta_{2j} } \delta_{\alpha_{2j}\alpha_{2j+1} }\delta_{\beta_{L-1}\beta_L}$ and $\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \alpha_{L}\beta_L}=\prod_{j=1}^{L/2-1}\delta_{\alpha_{2j-1}\alpha_{2j} }\delta_{\beta_{2j}\beta_{2j+1} }\delta_{\alpha_{L-1}\alpha_L}$,
we have the linear equations
\begin{eqnarray}
     A_{\alpha_1\alpha_{L}}\chi^{(A)}_{\beta_1\alpha_2 \dots \alpha_L}+ B_{\beta_1 \beta_{L}}\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}= X_{\alpha_1\alpha_2\beta_2\dots\alpha_L}\delta_{\beta_L\beta_{1}} +Y_{\beta_1\alpha_2\dots\beta_L}\delta_{\alpha_L\alpha_{1}}~.
\end{eqnarray}
The equations such that $\chi^{(A)}_{\beta_1\alpha_2 \dots \beta_{L-1}\beta_L}=1$, $\alpha_{L-1}\neq \alpha_{L}$ and $\beta_{1}\neq \beta_{L}$ give us $A_{\alpha_1\alpha_L} = Y_{\dots}\delta_{\alpha_1\alpha_L} \equiv a \delta_{\alpha_1\alpha_L}$; while the equations such that $\chi^{(B)}_{\alpha_1\alpha_2\beta_2 \dots \beta_{L-1}\beta_L}=1$, $\beta_{L-1}\neq \beta_{L}$ and $\alpha_{1}\neq \alpha_{L}$ give us 
$B_{\beta_1\beta_L} = X_{\dots}\delta_{\beta_1\beta_L} \equiv b \delta_{\beta_1\beta_L}$.






\subsection{Existence of the spectral gap in the thermodynamic limit}
To show the existence of the spectral gap in the thermodynamic limit, we will use Theorem 2 in Ref.~\cite{nachtergaeleSpectral1996}. Our ground states in a open chain are indeed a type of the \textit{generalized valence bond solid} (GVBS) states defined in Ref.~\cite{nachtergaeleSpectral1996}.
Consider an interval $[M,N]$ in an infinite 1d lattice, we can define our Hamiltonian in this subregion $H_{[M,N]} = \sum_{j=M}^{N-1}h_{j}$, where $h_j = (1-P_{j,j+1})$. 
Now in the previous section, we have shown that the ground state subspace is indeed the GVBS states (condition $\mathcal{F}_{\omega}=\mathcal{F}_{h}$ in Theorem 2 of Ref.~\cite{nachtergaeleSpectral1996}). Due to the projector construction of the Hamiltonian, we also have the condition $\text{Ker}(H)_{[M,N]} = \mathcal{G}_{M,N}$, where $\mathcal{G}_{M,N}$ is the span of the eigenvectors obtained from the reduced density matrix of $\omega$ on the region $[M,N]$.
We therefore conclude that our model and ground states satisfy the conditions of Theorem 2 in Ref.~\cite{nachtergaeleSpectral1996}, which shows the existence of the spectral gap in the thermodynamic limit.


\subsection{Symmetry-protected edge states under open boundary condition}

As pointed out in the main text, the ground states of Eq.~(\ref{eq:1ddimerH_SM}) have dangling edge degrees of freedom, and here we show that they are indeed protected by the symmetry $\mathbb{Z}^x_n \times \mathbb{Z}^z_n$, by showing that these edge degrees of freedom transform in projective representations under $\mathbb{Z}^x_n \times \mathbb{Z}^z_n$.

The ground states under open boundary condition are given by Eq. \eqref{eq: ground states}.
Recall the symmetry $\mathbb{Z}^x_n$ is implemented by $\prod_j \mu^x_j \nu^x_j$ and $\mathbb{Z}^z_n$ by $\prod_{j \in \text{even}} (\mu^z_j \nu^z_j)\prod_{j \in \text{odd}}(\mu^z_j \nu^z_j)^{-1}$.
To determine the effective symmetry actions on the edges, it suffices to examine the action of these symmetries on the ``dimer" $\sum_{d=1}^n|dd\rangle$, which is the building block of the bulk. One can easily verify that $\mu^x_1 \mu^x_2\sum_{d=1}^n|dd\rangle=\sum_{d=1}^n|dd\rangle$, and similarly for $\nu^x$; similarly, $\mu^z_1 (\mu^z_2)^{-1}\sum_{d=1}^n|dd\rangle=\sum_{d=1}^n|dd\rangle$ and likewise for $\nu^z_1 (\nu^z_2)^{-1}$.
The bulk of the wavefunction is therefore invariant under the symmetry action, and the symmetry action only transforms the edge degrees of freedom. 
Moreover, the ground states transform among the same $A$-type $|\psi^{(A)}_{a,b}\rangle$ or same $B$-type of states $|\psi^{(B)}_{a,b}\rangle$.
It is then easy to identify the edge symmetry operator. 
For example, if $L$ is even, within the $A$-type ground state space $\text{span}\{|\psi^{(A)}_{a,b}\rangle\}$,
\begin{eqnarray}
\begin{split}
    &W^x_L = \mu^x_{1}~, \quad &W^x_R = \mu^x_{L}~, \notag \\
    &W^z_L = (\mu^z_{1})^{-1}~, \quad &W^z_R = \mu^z_{L}~;
\end{split}
\end{eqnarray}
while 
\begin{eqnarray}
\begin{split}
    &W^x_L = \nu^x_{1}~, \quad &W^x_R = \nu^x_{L}~, \notag \\
    &W^z_L = (\nu^z_{1})^{-1}~, \quad &W^z_R = \nu^z_{L}~,
\end{split}
\end{eqnarray}
within the $B$-type ground state space $\text{span}\{|\psi^{(B)}_{a,b}\rangle\}$. The case for an odd $L$ is similar.


It is straightforward to check that the actions of $\mathbb{Z}_n^x$ and $\mathbb{Z}_n^z$ no longer commute on the edge DOF, which means that these edge DOF transform in projective representations under $\mathbb{Z}_n^x\times\mathbb{Z}_n^z$. Projective representations of $\mathbb{Z}_n^x\times\mathbb{Z}_n^z$ are classified by $H^2(\mathbb{Z}_n^x\times\mathbb{Z}_n^z)=\mathbb{Z}_n$, and they can be characterized by an integer $\eta$ that is defined modulo $n$, which is given by $W^x_{L, R}W^z_{L, R}=e^{i\frac{2\pi}{n}\eta_{L, R}}W^z_{L, R}W^x_{L, R}$. Using the above symmetry actions, we see that $\eta_L=-1$ and $\eta_R=1$, for both the $A$-type and $B$-type ground states. Since the degeneracy between the $A$-type and the $B$-type subspace is protected by the spontaneous breaking of the reflection symmetry $T\sigma$, we conclude that all these edge modes are protected by the unbroken symmetries and spontaneous symmetry breaking.






\subsection{Beyond the group-cohomology-based classification}

In the above section, we have shown that the ground states we obtain are nontrivial SPTs under the remaining $\mathbb{Z}_n^x\times\mathbb{Z}_n^z$ symmetry. Since the remaining symmetry in fact contains a $\mathbb{Z}_n^x\times\mathbb{Z}_n^z\times\sigma$ symmetry, we can ask what kind of SPT our ground states are under this symmetry. We will see that our ground states are actually beyond the conventional group-cohomology-based classification developed in Refs.~\cite{Jiang2016, Thorngren2016}.

The argument is a generalization of that in Ref.~\cite{Jiang2019}, which only discusses the case with $n=4$ and has no spontaneous symmetry breaking. The group-cohomology-based classification of $1+1$ dimensional $\mathbb{Z}_n^x\times\mathbb{Z}_n^z\times\sigma$ SPTs is
\begin{eqnarray}
H^2(\mathbb{Z}_n^x\times\mathbb{Z}_n^z\times\sigma, U(1)_\sigma)
=
\left\{
\begin{array}{lr}
\mathbb{Z}_2, & {\rm odd\ }n\\
\mathbb{Z}_2^4, & {\rm even\ }n
\end{array}
\right.
\end{eqnarray}
where the subscript $\sigma$ in $U(1)$ means that the $U(1)$ phase factor should be complex conjugated when acted by $\sigma$. 
For odd $n$, the nontrivial SPT is protected only by $\sigma$, which has no protected edge state. Since our states do have protected edge states, they are beyond this classification. For even $n$, the only possible edge states protected by the $\mathbb{Z}_n^x\times\mathbb{Z}_n^z$ symmetry from this classification have $\eta=n/2$. In our case, $\eta=\pm 1$, so our states are also beyond the group-cohomology-based classification if $n>2$. As noted in Ref. \cite{Jiang2019}, here going beyond the group-cohomology-based classification is possible because the states at each site transform in projective representations of the on-site symmetries. To capture such states, one needs a more refined classification that takes into account the possible projective representations of the states at each site, such as the ones proposed in Refs.~\cite{Else2019, Jiang2019}.




\section{The $G_0=SO(3)\times\mathbb{Z}^d\rightarrow G=\mathbb{Z}_2\times\mathbb{Z}^d$ symmetry-breaking orders} \label{app: maintext example with continuous SSB}

In this section, we discuss in more detail the EESBO in spin-$S$ lattice system with $G_0=SO(3)\times\mathbb{Z}^d$ and $G=\mathbb{Z}_2\times\mathbb{Z}^d$. We will first explain why the symmetry-breaking order is ``classical" in the case with $S\geqslant 1$. Then we will show that the remaining symmetry of the wavefunction presented in the main text is indeed $G$. Next, we will show that this wavefunction represents a $G$-symmetric short-range entangled state that is beyond the conventional group-cohomology-based classification of Refs. \cite{Jiang2016, Thorngren2016}, and in some sense even beyond the more refined classification of Refs. \cite{Jiang2019, Else2019}. Finally, we will show that for all $S\geqslant 1/2$, this symmetry-breaking order generically has two gapless Goldstone modes, where one of them has a linear dispersion and the other has a quadratic dispersion.


\subsection{Classical symmetry-breaking orders for $S\geqslant 1$}

To start, first recall that any spin-$S$ state can uniquely be represented $2S$ points on the Bloch sphere \cite{Majorana1932, bloch1945, Makel2009}; this representation is known as the Majorana representation, and it is a generalization of the familiar statement that any spin-$1/2$ state can be uniquely represented by a single point on the Bloch sphere. The validity of this representation is based on the fact that all spin-$S$ states can be constructed by symmetrizing $2S$ spin-1/2 states. The advantage of the Majorana representation is that the $SO(3)$ spin rotation acts as a rigid body rotation of these $2S$ points.

Now it is easy to see that for all $S\geqslant 1$, there are single-site spin-$S$ states with a $\mathbb{Z}_2$ symmetry generated by the $\pi$ rotation with respect to $S_z$, and one such example is to put 2 of the $2S$ points at spherical coordinates with polar angle anything but not $\pi/2$, and azimuthal angles $0$ and $\pi$ for each of them. The other $2S-2$ points are put on the north pole. It is easy to check that such a state has no symmetry other than the $\mathbb{Z}_2$. Therefore, the $G_0\rightarrow G$ symmetry-breaking order can be represented by an ultralocal product state for $S\geqslant 1$, and a representative wavefunction is simply taking the spin at each site to be the $\mathbb{Z}_2$ symmetric single-site state discussed here.


\subsection{EESBO for $S=1/2$}

Next, we will show that the one dimensional matrix product state (MPS) represented by the 
matrices
\begin{eqnarray} \label{eq: MPS app}
A^{[S_z=\frac{1}{2}]}=
\left(
\begin{array}{cc}
1+a & 0 \\
0 & 1-a
\end{array}
\right),
\quad
A^{[S_z=-\frac{1}{2}]}=
\left(
\begin{array}{cc}
0 & b \\
c & 0
\end{array}
\right)
\end{eqnarray}
is short-range entangled, and the only subgroup of $SO(3)$ that is still a symmetry of this state is a $\mathbb{Z}_2$ symmetry generated by the $\pi$ rotation with respect to $S_z$, for generic nonzero real-valued $a, b, c$. This implies that the stack of chains described in the main text is also a short-range entangled state, whose symmetry is $G=\mathbb{Z}_2\times\mathbb{Z}^d$.

Showing that this state is short-range entangled amounts to show that this MPS is injective. To this end, we first construct the tensor
\begin{eqnarray} \label{eq: double tensor app}
T_{i, j, k, l}=\sum_{S_z=\pm\frac{1}{2}}A^{[S_z]}_{ij}\left(A^{[S_z]}_{kl}\right)^*~,
\end{eqnarray}
from which we construct the $4\times 4$ transfer matrix
\begin{eqnarray} \label{eq: transfer app}
\tilde T_{(ik), (jl)}=T_{i, j, k, l}~.
\end{eqnarray}
Here the relation between the indices of the tensor $T$ and those of the transfer matrix $\tilde T$ is that $(11)\leftrightarrow 1$, $(12)\leftrightarrow 2$, $(21)\leftrightarrow 3$ and $(22)\leftrightarrow 4$. Using Eqs. \eqref{eq: MPS app} and \eqref{eq: double tensor app}, we get
\begin{eqnarray}
\tilde T=
\left(
\begin{array}{cccc}
 (a+1)^2 & 0 & 0 & b^2 \\
 0 & 1-a^2 & b c & 0 \\
 0 & b c & 1-a^2 & 0 \\
 c^2 & 0 & 0 & (a-1)^2 \\
\end{array}
\right)
\end{eqnarray}
which has eigenvalues $\{1+a^2\pm\sqrt{4a^2+b^2c^2}, 1-a^2\pm bc\}$. The largest eigenvalue is non-degenerate for all $a\neq 0$. That is, this MPS is injective and represents a short-range entangled state as long as $a\neq 0$.

Below we show that $\mathbb{Z}_2$ (generated by the $\pi$ rotation with respect to $S_z$) is the only subgroup of the $SO(3)$ spin rotation symmetry that is still a symmetry of this MPS. Our strategy is to first consider an operation whose action on the physical states is given by $U(\theta)=\exp\left(iS_z\theta\right)$, where $0\leqslant\theta<2\pi$. These are just all spin rotating with respect to $S_z$.
We will see that the MPS given by Eq.~\eqref{eq: MPS app} is invariant only if $\theta=0, \pi$. This means that the $\mathbb{Z}_2$ generated by the $\pi$ rotation with respect to $S_z$ is indeed a symmetry of this MPS. We then show that this is the only symmetry of that MPS, by showing that $\langle S_x\rangle=\langle S_y\rangle=0$ while $\langle S_z\rangle\neq 0$ for this MPS.

$U(\theta)=\exp\left(iS_z\theta\right)$ is a symmetry of the MPS if and only if \cite{Perez-Garcia2008, Sanz2009}
\begin{eqnarray}
U(\theta)_{S_zS_z'}A^{[S_z']}=e^{i\alpha}VA^{[S_z]}V^\dag
\end{eqnarray}
with $\alpha\in\mathbb{R}$ and $V\in SU(2)$. Here $\alpha$ and $V$ can depend on $\theta$. In the eigenbasis of $S_z$, $U(\theta)$ is diagonal, and this equation simply becomes
\begin{eqnarray}
e^{i\frac{\theta}{2}}A^{[S_z=\frac{1}{2}]}=e^{i\alpha} VA^{[S_z=\frac{1}{2}]}V^\dag, \quad
e^{-i\frac{\theta}{2}}A^{[S_z=-\frac{1}{2}]}=e^{i\alpha} VA^{[S_z=-\frac{1}{2}]}V^\dag~,
\end{eqnarray}
Now we use the fact that any $2\times 2$ matrix has a unique expansion in terms of the identity matrix and the Pauli matrices, $\sigma_{x, y, z}$. In our case, we can rewrite $A^{[S_z=\frac{1}{2}]}=1+a\sigma_z$ and $A^{[S_z=-\frac{1}{2}]}=\frac{b+c}{2}\sigma_x+\frac{b-c}{2}i\sigma_y$, so the above equation becomes
\begin{eqnarray}
e^{i\frac{\theta}{2}}(1+a\sigma_z)=e^{i\alpha} V(1+a\sigma_z)V^\dag, \quad
e^{-i\frac{\theta}{2}}\left(\frac{b+c}{2}\sigma_x+\frac{b-c}{2}i\sigma_y\right)=e^{i\alpha} V\left(\frac{b+c}{2}\sigma_x+\frac{b-c}{2}i\sigma_y\right)V^\dag~.
\end{eqnarray}
Under the conjugation by $V\in SU(2)$, the identity matrix is invariant while the Pauli matrices transform as a vector under $SO(3)$. So, when $a\neq 0$, for the first equation to hold, we must have
$e^{i\alpha}=e^{i\frac{\theta}{2}}$ and $V=\exp(i\frac{\sigma_z}{2}\phi)$ for some $\phi\in\mathbb{R}$. Then the second equation becomes
\begin{eqnarray}
\begin{split}
e^{-i\theta}\left(\frac{b+c}{2}\sigma_x+\frac{b-c}{2}i\sigma_y\right)
&=\frac{b+c}{2}(\cos\phi\sigma_x+\sin\phi\sigma_y)+i\frac{b-c}{2}(-\sin\phi\sigma_x+\cos\phi\sigma_y)\\
&=\frac{be^{-i\phi}+ce^{i\phi}}{2}\sigma_x+\frac{be^{-i\phi}-ce^{i\phi}}{2}i\sigma_y~,
\end{split}
\end{eqnarray}
which implies that
\begin{eqnarray}
e^{-i\theta}(b+c)=be^{-i\phi}+ce^{i\phi},
\quad
e^{-i\theta}(b-c)=be^{-i\phi}-ce^{i\phi}~.
\end{eqnarray}
If $bc\neq 0$, then $e^{i\theta}=e^{i\phi}=e^{-i\phi}$, giving us $e^{i\theta}=\pm 1$, {i.e., } $\theta=0$ or $\theta=\pi$. Therefore, within the group of $U(1)$ spin rotations with respect to $S_z$, the $\mathbb{Z}_2$ generated by $\pi$ rotation is the only symmetry of the MPS in Eq. \eqref{eq: MPS app}. Note that if any of $b$ and $c$ vanishes, any $\theta=\phi$ solves these equations and the MPS has a $U(1)$ symmetry corresponding to the $S_z$ rotation.

The above $\mathbb{Z}_2$ symmetry implies that $\langle S_x\rangle=\langle S_y\rangle=0$ for the MPS in Eq. \eqref{eq: MPS app}. As long as we can show that $\langle S_z\rangle\neq 0$, we know that the only subgroup of $SO(3)$ spin rotation which can possibly be a symmetry of Eq. \eqref{eq: MPS app} must be a subgroup of $U(1)$ spin rotations with respect to $S_z$. Then combined with the previous result, we conclude that the only subgroup of $SO(3)$ spin rotation that is still a symmetry of Eq.~\eqref{eq: MPS app} is the $\mathbb{Z}_2$ generated by $\pi$ rotation with respect to $S_z$.

In the following, we use the standard transfer matrix method to calculate $\langle S_z\rangle$ for a spin-1/2 described by the MPS in Eq.~\eqref{eq: MPS app}, which contains $L$ sites and has periodic boundary condition. To this end, we first consider the tensor
\begin{eqnarray}
M_{i, j, k, l}=\sum_{m, m'=\pm\frac{1}{2}}A_{ij}^{[m]}\hat S_{z, mm'}\left(A_{kl}^{[m']}\right)^*
\end{eqnarray}
and convert it into a matrix in a way similar to Eq.~\eqref{eq: transfer app}:
\begin{eqnarray}
\tilde M_{(ik), (jl)}=M_{i, j, k, l}
\end{eqnarray}
where $\hat S_z={\rm diag}(1/2, -1/2)$. Explicitly, we have
\begin{eqnarray}
\tilde M=
\left(
\begin{array}{cccc}
 (a+1)^2 & 0 & 0 & -b^2 \\
 0 & 1-a^2 & -b c & 0 \\
 0 & -b c & 1-a^2 & 0 \\
 -c^2 & 0 & 0 & (a-1)^2 \\
\end{array}
\right)
\end{eqnarray}
Then
\begin{eqnarray} \label{eq: Sz formula}
\langle S_z\rangle=\frac{{\rm Tr}(\tilde M\tilde T^{L-1})}{{\rm Tr}(\tilde T^{L})}
\end{eqnarray}

To evaluate the above expression, we perform Jordan decomposition for $\tilde T$ and $\tilde M$ to obtain $\tilde T=S_TJ_TS_T^{-1}$ and $\tilde M=S_MJ_MS_M^{-1}$,
where
\begin{eqnarray}
S_T=
\left(
\begin{array}{cccc}
 0 & 0 & \frac{2 a-\sqrt{4 a^2+b^2 c^2}}{c^2} & \frac{\sqrt{4 a^2+b^2 c^2}+2 a}{c^2} \\
 -1 & 1 & 0 & 0 \\
 1 & 1 & 0 & 0 \\
 0 & 0 & 1 & 1 \\
\end{array}
\right),
\quad
S_M=
\left(
\begin{array}{cccc}
 0 & 0 & \frac{\sqrt{4 a^2+b^2 c^2}-2 a}{c^2} & \frac{-\sqrt{4 a^2+b^2 c^2}-2 a}{c^2} \\
 1 & -1 & 0 & 0 \\
 1 & 1 & 0 & 0 \\
 0 & 0 & 1 & 1 \\
\end{array}
\right)
\end{eqnarray}
and $J_T=J_M={\rm diag}(1-a^2-bc, 1-a^2+bc, 1+a^2-\sqrt{4a^2+b^2c^2}, 1+a^2+\sqrt{4a^2+b^2c^2})$. Plugging these into Eq. \eqref{eq: Sz formula} and taking the limit $L\rightarrow\infty$ yield
\begin{eqnarray}
\langle S_z\rangle=\frac{1+a^2+\frac{4 a^2-b^2 c^2}{\sqrt{4 a^2+b^2 c^2}}}{1+a^2+\sqrt{4 a^2+b^2 c^2}}~.
\end{eqnarray}
Clearly $\langle S_z\rangle\neq 0$ for generic nonzero $a, b, c\in\mathbb{R}$.

Therefore, we have completed the proof that, for generic nonzero $a, b, c\in\mathbb{R}$, the only subgroup of $SO(3)$ that is still a symmetry of the MPS in Eq.~\eqref{eq: MPS app} is the $\mathbb{Z}_2$ generated by the $\pi$ rotation with respect to $S_z$, which further implies that the remaining symmetry of the stack of chains described in the main text is $G=\mathbb{Z}_2\times\mathbb{Z}^d$.


\subsection{$G$-symmetric short-range entangled state beyond the group-cohomology-based classification}

In the previous section, we have already shown that the stack of spin-1/2 chains described in the main text is a $G$-symmetric short-range entangled state. It is then natural to ask how this state fits into the classification of $d+1$ dimensional $G$ symmetry-protected topological states ($G$-SPTs). In this subsection, we show that this stack of spin-1/2 chains is beyond the classification of $d+1$ dimensional $G$-SPTs develeped in Refs. \cite{Jiang2016, Thorngren2016}.

According to Refs. \cite{Jiang2016, Thorngren2016}, the classification of $d+1$ dimensional $G$-SPTs is $H^{d+1}(G, U(1))$, with $G=\mathbb{Z}_2\times\mathbb{Z}^d$. This group cohomology has been calculated in Ref. \cite{Cheng2015}, and it is found that 
\begin{eqnarray}
H^{d+1}(G, U(1))=\prod_{k=1}^{d+1}H^{k}\left(\mathbb{Z}_2, U(1)\right)^{\left(\begin{array}{c} d \\ k-1 \end{array}\right)}
\end{eqnarray}
This classification has a simple interpretation in terms of dimensional reduction: to construct a $G$-SPT in $d$ spatial dimensions, we can first pick up $k-1$ spatial dimensions, build up $(k-1)+1$ spacetime dimensional $\mathbb{Z}_2$-SPTs along these dimensions, and stack these lower dimensional $\mathbb{Z}_2$-SPTs in a way that preserves the $\mathbb{Z}^d$ translation symmetry. Using the fact that $H^k(\mathbb{Z}_2, U(1))=\mathbb{Z}_2$ for odd $k$ and $H^k(\mathbb{Z}_2, U(1))=\mathbb{Z}_1$ for even $k>0$, we can get the classification of $G$-SPTs in various dimensions. When $d=1, 2$ and 3, the classification is $\mathbb{Z}_2$, $\mathbb{Z}_2^2$ and $\mathbb{Z}_2^4$, respectively. The analysis can be extended to higher $d$, in which case it is still true that our stack of spin-1/2 chains is beyond this group-cohomology-based classification, but below we will focus on the case with $d\leqslant 3$.

Let us start with the case where $d=1$ and the classification is $\mathbb{Z}_2$. The dimensional reduction picture implies that the nontrivial $1+1$ dimensional $G$-SPT, roughly speaking, has each of its translation unit cells hosting a $\mathbb{Z}_2$ odd state.
More precisely, it means that when the length of the chain increases by 1, the $\mathbb{Z}_2$ eigenvalue of the ground state changes. In other words, the two different $1+1$ dimensional $G$-SPTs are expected to be characterized by $\lambda\equiv\lim_{L\rightarrow\infty}\frac{\lambda_{L+1}}{\lambda_L}$, where $\lambda_L$ is the eigenvalue of the ground state with length $L$, $|\psi\rangle_L$, under the $\mathbb{Z}_2$ symmetry action $X$, {i.e., } $X|\psi\rangle_L=\lambda_L|\psi\rangle_L$. The trivial (nontrivial) $G$-SPT should have $\lambda=1$ ($\lambda=-1$). A simple example of a trivial (nontrivial) $1+1$ dimensional $G$-SPT is a many-qubit state $\prod_j|0\rangle_j$ ($\prod_j|1\rangle_j$), where the on-site $\mathbb{Z}_2$ symmetry action $X_j$ satisfies $X_j|0\rangle_j=|0\rangle_j$ and $X_j|1\rangle_j=-|1\rangle_j$. Note that both of these SPTs are ultralocal product states.

For our MPS state described by Eq. \eqref{eq: MPS app}, the results from the previous subsection imply that $\frac{\lambda_{L+1}}{\lambda_L}=e^{i\alpha}=i$, which is neither $1$ nor $-1$. This means that this simple MPS is a $G$-symmetric short-range entangled state beyond the conventional group-cohomology classification, which is possible because in our case the $\mathbb{Z}_2$ on-site symmetry acts on states at each site as a $\mathbb{Z}_4$ symmetry. Notice that this is different from having projective representations under the on-site symmetry ($\mathbb{Z}_2$ actually has no projective representation), so in some sense such a state is not even captured by the more refined classification proposed in Refs. \cite{Else2019, Jiang2019}. Instead, it is more appropriate to view our state as being described by a torsor over $H^2(\mathbb{Z}_2\times\mathbb{Z}, U(1))$, rather than an element of $H^2(\mathbb{Z}_2\times\mathbb{Z}, U(1))$.

Next, we move to the case with $d=2$, where the classification is $\mathbb{Z}_2^2$. It is easy to see that one of the two root SPTs is simply the nontrivial SPT protected only by $\mathbb{Z}_2$, where the translation symmetry is unimportant, and the other root SPT is simply the higher dimensional generalization of the $1+1$ dimensional $G$-SPT discussed above. Roughly speaking, each translation unit cell in this SPT hosts a $\mathbb{Z}_2$ odd state. Clearly, our stack of spin chains is neither of these two states, which means it is beyond the group-cohomology-based classification.

Finally, we turn to the case with $d=3$, where the classification is $\mathbb{Z}_2^4$. One of the four root states is simply again the generalization of the $1+1$ dimensional state, where each translation unit cell hosts a $\mathbb{Z}_2$ odd state. The other three root states can be viewed as stacks of $2+1$ dimensional $\mathbb{Z}_2$-SPTs that have $\mathbb{Z}^3$ translation symmetry. Again, it is clear that our stack of spin chains is none of these states, and it is therefore beyond the group-cohomology-based classification.
Instead, such a state should be thought of as a torsor over $H^2(\mathbb{Z}_2\times\mathbb{Z}, U(1))$, rather than an element in $H^2(\mathbb{Z}_2\times\mathbb{Z}, U(1))$.


\subsection{Properties of the Goldstone modes}

In this subsection, we discuss the properties of the Goldstone modes associated with the $G_0\rightarrow G$ symmetry-breaking orders. This discussion applies to all $S\geqslant 1/2$.

We will use the results in Ref.~\cite{Watanabe2014}, which connect the number of broken generators of a continuous symmetry and the expectation values of the charge density in the ground states to the number and dispersion of the Goldstone modes. The expectation values of the charge density enter through the anti-symmetric matrix $\rho$ defined as
\begin{eqnarray}
\rho_{ab}\equiv -i\langle[\hat S_a, \hat S_b]\rangle/V~,
\end{eqnarray}
where $\hat S_{a, b}$ are the broken generators, and $V$ is the volume of the system. Denote the number of the broken generators by $n_{BG}$, Ref.~\cite{Watanabe2014} shows that the total number of Goldstond modes is $n_{BG}-\frac{1}{2}{\rm rank}(\rho)$, where $n_{BG}-{\rm rank}(\rho)$ of them generically have linear dispersion, and the other $\frac{1}{2}{\rm rank}(\rho)$ of them generically have quadratic dispersion.

In our case, the number of broken generators is $n_{BG}=3$. Because the remaining symmetry of $G$ does not require $\langle\vec S\rangle=0$, the matrix $\rho$ generically has rank $2$. So there will be $n_{BG}-\frac{1}{2}{\rm rank}(\rho)=3-2/2=2$ gapless Goldstone modes, where one of them generically has linear dispersion and the other has quadratic dispersion.



\section{More examples of EESBOs with spontaneously broken continuous symmetries} \label{app: more EESBOs with continuous SSB}

In this section, we discuss two more examples of EESBOs with spontaneously broken continuous symmetries.

The first example is a $d$-dimensional spin-1 lattice system with one site per unit cell, where $d>1$. Suppose its Hamiltonian has $PSU(3)$ and translation symmetries, {i.e., } $G_0=PSU(3)\times \mathbb{Z}^d$, where each site of the lattice hosts a spin-1 moment that transforms in the fundamental representation of $SU(3)$, which is a projective representation of $PSU(3)$. Consider the case $G=SO(3)\times \mathbb{Z}^d$, {i.e., } the translation is intact while $PSU(3)$ is spontaneously broken to $SO(3)$. Here the translation symmetry ensures that all local product states must be ultralocal, because if any two sites were entangled, translation symmetry requires that extensively many sites are entangled. But ultralocal product states are clearly not invariant under the remaining on-site $SO(3)$ symmetry. So this symmetry-breaking order is entanglement-enabled. Its representative wave function can be simply taken as a stack of many copies of the Affleck-Kennedy-Lieb-Tasaki (AKLT) chains arranged to have the $\mathbb{Z}^d$ translation symmetry \cite{AKLT}. If desired, one can also construct a cat state that preserves the full $G_0$ symmetry by summing over the symmetry-breaking orbits. Because each AKLT is short-range entangled, the entire stack is also short-range entangled. In other words, there is no LSM constraint associated with the remaining symmetry $G$. Also, this representative wavefunction shows that this EESBO is a within-group-cohomology state under the remaining symmetry $G$.

The second example is a honeycomb lattice spin-1/2 system with $G_0= SO(3)\times\mathbb{Z}_2^T\times p6m$, where the spins live on the sites of a honeycomb lattice with a $p6m$ lattice symmetry, and they carry spin-1/2 under the $SO(3)$ spin rotation symmetry and transform as a Kramers doublet under the $\mathbb{Z}_2^T$ time reversal symmetry. Such a system is free of LSM constraints \cite{Kim2015}. Now consider a spin-nematic order where the remaining symmetry $G=\mathbb{Z}_2\times\mathbb{Z}_2\times\mathbb{Z}_2^T\times p6m$, {i.e., } the lattice symmetry and time reversal symmetry are intact, while the $SO(3)$ symmetry is broken to $\mathbb{Z}_2\times\mathbb{Z}_2$ generated by $\pi$ spin rotations around three orthogonal axes. It is straightforward to check that the lattice symmetry again requires a local product state to be ultralocal, which must then break the remaining $\mathbb{Z}_2\times\mathbb{Z}_2\times\mathbb{Z}_2^T$ on-site symmetry, because each spin transforms in its projective representation. So such a spin-nematic order is entanglement-enabled. To obtain the wave function of this state, one can start from the $G_0$-symmetric wavefunction in Ref.~\cite{Kim2015} and perturb it so that $G_0$ is broken to $G$. As a spin-nematic order, there are gapless Goldstone modes carrying spin-2. Furthermore, Ref.~\cite{Ye2021} predicted that the vicinity of this EESBO may realize an exotic quantum spin liquid that cannot be described by the usual parton approach.



