\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.


\section{Introduction}
Because erroneous or delayed responses in self-driving applications can compromise road safety, equipment, and/or the lives of the passengers themselves, Autonomous Driving Systems (ADS) are required to achieve outstanding performances on core driving tasks, such as perception and localization. Consequently, ADS platforms are designed today to run highly-sophisticated algorithms on intricate hardware architectures to realize such desired levels of performance while being robust to any adverse driving contexts. For that, modern ADS platforms have adopted a multi-modal processing approach for constructing ensemble perspectives of driving scenes using a diverse set of sensory inputs, as in how the Tesla Autopilot systems possess 8 cameras and 12 ultrasonic sensors \cite{tesla}.

Such a multi-sensor approach leads to the generation of an enormous volume of high-dimensional data that requires tremendous resources for real-time processing, further adding to the power demands of the entire system. Addressing this, a heterogeneous collection of hardware components, as in Application-Specific Integrated Circuits (ASICs) and GPUs, are commonly integrated onto ADS platforms to balance performance demands and power efficiency \cite{lin2018architectural}. Still, hardware advancements are met with growing algorithmic complexity and the requirement for supporting new features, leading the power footprint to remain relatively high. For instance, if we compare two generations of ADS platforms: the Nvidia Drive PX2, which was used by Tesla and Audi Q7 for their autopilot programs \cite{px2tesla, px2audi}, against its successor, the Nvidia Drive AGX Orin \cite{Abuelsamid2020orin}, we find that performance efficiency aside, the baseline power demands increased from 250 W to 800 W, which in theory can have adverse effects on both the thermal comfort of the passengers and the vehicle's driving range \cite{lin2018architectural}.
\vspace{-1.2ex}
\begin{figure}[!h]
\centering
{\includegraphics[,width = 0.46\textwidth]{Figures/edge_computing.pdf}}
\vspace{-1.5ex}
\caption{Vehicular Edge Computing (VEC) Architecture}
\label{fig:hierarchy}
\vspace{-2ex}
\end{figure}

Given how the bulk of processing in ADS modules is largely dominated by deep neural networks (DNNs), compression techniques, e.g., quantization and pruning, have been considered to reduce the modules' complexity, and in turn, their resource requirements  \cite{li2018auto-tuning}. However, experienced performance degradation poses a concern with regards to adopting such techniques' for this class of critical applications.
Alternatively, recent research efforts have targeted exploiting the emerging edge computing paradigm
for autonomous driving and other vehicular services, in which cloud computing capabilities are brought to the edge of the network through edge servers deployed close to the edge devices, enabling offloading of cumbersome processing burdens to these edge servers for better resource management \cite{liu2019edge}. In the context of vehicular applications, the edge devices are known as connected vehicles (CVs) and the paradigm is further specified as vehicular edge computing (VEC), or vehicular fog computing (VFC) in other cases \cite{baidya2020vehicular}.

Figure \ref{fig:hierarchy} illustrates the hierarchical architecture of VEC and its reliance on the wireless infrastructure, where the higher bandwidths and ultra-low latencies promised by the forthcoming 5G and Dedicated Short Range Communications (DSRC) technologies are to be instrumental in VEC's wide-scale adoption \cite{liu2019edge}. In this regard, VEC computing servers are expected to be deployed either at Road Side Units (RSUs) -- as part of the V2X paradigm -- or cellular base stations, where recent works have proposed to optimize the offloading process to minimize the overall latency and the energy consumed by the ADS \cite{malawade2021sage, baidya2020vehicular}. Still, we find that current approaches in the literature are lacking in the following departments:

 \begin{itemize}
    \item The driving context -- crucial to robustness -- is not factored in the offloading decision. Given how a scene's complexity directly correlates with the risk level, and delayed responses could lead to consequences with different levels of severity
    \item Adapting the offloading load according to the runtime conditions while accounting for the underlying ADS composition (e.g., concurrent sensor processing pipelines) is overlooked
\end{itemize}
\vspace{-2ex}

\subsection{Motivational Example}
In Figure \ref{fig:motivation}, we show two frames from the Radiate dataset \cite{sheeny2021radiate}, and compare their mean average Precision (mAP) scores on the object detection task. As shown, the left frame instantiates a complex scene with numerous objects of diverse classes, some of them superimposed or obstructed from view leading to relatively low mAP scores. Still, we observe that through fusing the outputs from all sensory pipelines (2 cameras, lidar, and radar), an mAP score of 17.6\% is realized, which surpasses the highest score achieved by a standalone sensory pipeline output -- 11.7\% from the right stereo camera. This alludes to the power of sensor fusion as each sensor can capture its own unique set of features that complement those from other sensors to provide more comprehensive views of the driving scenes. On the flip side, the right frame contains a mere single vehicle that is easily detectable by the standalone camera, achieving a 100\% mAP score. From here, we can contemplate the desired behavior when VEC is supported with regards to tuning the operating mode. Specifically, an ADS experiencing plain driving scenes can opt for \textit{offloading} processing loads from a subset of sensory pipelines for resource efficiency, because although delayed server responses could cause some partial outputs to be absent by the execution deadline, partial fusion of the available local outputs would suffice for this time window due to the relative simplicity of the scene. Contrarily, complicated scenes should have all sensor outputs available for fusion to stimulate robustness, which is only achieved during \textit{local execution} mode as the uncertainty of the wireless networks is avoided. This behavior would be learned by our proposed solution as will be detailed in the following sections.

\subsection{Novel Contributions}

To address the above limitations, we present a methodology for \underline{R}obust Task \underline{O}ffloading in Modular \underline{M}ulti-Sensor \underline{A}uto\underline{n}omo\underline{us} Driving Systems, namely \textsc{Romanus}. From here, we can summarize the main contributions of this paper as follows:
\begin{itemize}
    \item We present \textsc{Romanus}, a methodology to support efficient and robust offloading for modular ADS platform comprising multiple sensory pipelines with support for sensor fusion.
    \item As far as our knowledge goes, we are the first to factor the driving context in the offloading decision for the robustness of autonomous driving.
    \item We integrate optimal offloading points within each sensor processing model to realize a dynamic decision space for the runtime operating modes of the ADS.
    \item We implement a Deep Reinforcement Learaning (DRL) based runtime solution that leverages contextual and temporal correlations in the data to optimize the offloading process for latency, energy, and robustness given concurrent pipelines.
    \item Experiments on the object detection use-case using a real-world driving dataset and an industry-grade ADS indicate that our approach is 14.99\% more energy-efficient than local execution while achieving a 77.06\% reduction in risky behavior form a robust-agnostic baseline.
\end{itemize}

 \begin{figure}[!tbp]
 \centering
{\includegraphics[,width = 0.46\textwidth]{Figures/Difficulty_motivation.pdf}}
\vspace{-2.5ex}
\caption{Two frames of different complexities showing single camera and late fusion mAP scores and the selected operational modes by our learning-based solution. The bounding boxes indicate the ground truths from the dataset.}
\label{fig:motivation}
\vspace{-4.5ex}
\end{figure}

\vspace{-1.5ex}

\section{Related Works} \label{sec:related-works}

\textit{\textbf{Mutli-Sensor Perception: }}To maximize information extraction from a driving scene, data is collected from a diverse set of sensors, e.g., cameras, lidar, and radar, to promote perception robustness. Mainly, There are two primary schemes for processing these multi-sensory inputs: \textit{early fusion} \cite{yoo20203d,shahian2019real} and \textit{late fusion} \cite{xu2018pointfusion}. The former combines all sensory features to a single feature at an early point in the ADS pipeline, but is susceptible to sensing noise.
Conversely, the latter offers more resilience at the expense of more redundancy across the sensor pipelines. Recent works \cite{malawade2022hydrafusion, malawade2022ecofusion} have also explored the potential of \textit{hybrid fusion} approaches to leverage the best of both worlds, albeit with added implementation complexities. Here, we concentrate our analysis on the standard \textit{late fusion} approach as it is more challenging and understudied from an offloading perspective.

\textit{\textbf{Vehicular Edge Computing (VEC): }} Numerous research efforts have targeted system-wide resource optimization for VEC through optimal task offloading and scheduling strategies given a variety of servers, vehicles, and tasks \cite{wu2019delay, zhang2016energy}. Typically, such strategies are complemented with runtime solutions that can tune the operation according to variations in the deployment environment, such as the network connectivity conditions \cite{cui2020offloading}. Nonetheless, delayed responses from edge servers are not tolerated in autonomous driving application as the safety of the road, vehicles, and passengers \cite{baidya2020vehicular} can be compromised. Hence, \cite{wang2019auto} proposed a customized communication protocol for a stable and fast offloading of autonomous driving tasks. Even more so, the authors in \cite{malawade2021sage} proposed a fail-safe routine to re-invoke local computation if responses are delayed beyond a certain threshold to account for the uncertainty of wireless links. Such schemes would be even more convoluted when offloading from multiple concurrent pipelines is considered.

\textit{\textbf{DNN Split Computing: }}To identify optimal offloading points within DNN architectures, \cite{neurosurgeon, odema2021lens} analyzed the expected computation and communication costs for each potential offloading layer. For a considerable number of architectures, either direct raw inputs offloading or pure local execution represented the most efficient option. Therefore, works in \cite{eshratifar2019bottlenet, matsubara2019distilled, matsubara2020head} proposed to modify a DNN's structure to include an early optimal offloading layer that shrinks the size of transmissible data, minimizing the costs of both computation and communication. This split-computing concept was applied for end-to-end control in autonomous vehicles \cite{malawade2021sage}, and here, we extend its applicability to multi-sensor modular ADS platforms.

\section{System and Problem Overview}

\subsection{Autonomous Driving System Composition}

For perspective, we briefly describe the primary modules that compose a state-of-the-art ADS shown in Figure \ref{fig:av-pipeline} as follows:

 \begin{figure}[!tbp]
 \centering
{\includegraphics[,width = 0.46\textwidth]{Figures/ADS.pdf}}
\vspace{-1ex}
\caption{Modular industry-grade ADS. The fusion point following the object detection module is this work's use-case.}
\label{fig:av-pipeline}
\vspace{-3ex}
\end{figure}
\textbf{\emph{Perception:}} As the main receptor of the raw sensory data, the perception module is responsible for processing the data over two successive computing blocks. The first is an \textit{object detector} to identify and classify objects of interest, e.g., pedestrians and vehicles, that surround the ego vehicle. A \textit{tracking} module ensues to receive identified objects and associate them with their past movements so as to predict the current movement trajectory.

\textbf{\emph{Localization:}} Another module taking in the raw inputs is the localization module, whose task is to pinpoint the position of the vehicle at high precision using SLAM/GPS modules.

\textbf{\emph{Planning:}} Outputs from the perception and localization blocks are fused together onto the same 3D co-ordinate space for further processing by the behavioral and motion planning block. From here, a series of sequential path information can be generated from starting position until the endpoint.

\textbf{\emph{Control:} }The final block tasked with mapping the information generated from the planning block onto control instructions for the actuators (e.g., driving wheel, brakes, accelerator).

\subsection{Problem Formulation}

In a modular ADS pipeline, the perception block is the dominant entity affecting the end-to-end performance the most \cite{lin2018architectural}, and thus, directing offloading optimizations at this module can lead to substantial efficiency gains across the entire system. Still, sub-optimal operating points can be reached if the following two aspects are not considered properly: (\emph{i}) the nominal safety considerations of the autonomous driving application, and/or (\emph{ii}) the structural composition of the ADS modules themselves.
For the former, an ADS is required to conclude end-to-end processing under stringent execution time limits to maintain road safety -- a 100 ms deadline at the worst \cite{lin2018architectural, baidya2020vehicular}. Hence, when VEC is supported, expected additional delays due to wireless channel impairments should be considered as part of the overall end-to-end latencies to determine the best offloading decision. Still, additional abrupt delays could threaten the integrity of the self-driving application considering the tightness of the execution windows.
Whereas for the latter, understanding the underlying structure of an ADS is key to determine the optimal placement of an offloading point that effectively balances the inherent trade-off between communication and computation. For instance, offloading prior to the fusion point can incur a substantial transmission overhead, as opposed to offloading after it which could incur a sizeable computational overhead due to prolonged periods of local processing.

Formally, a module employing late fusion comprises $N$ processing models $\{f_1, f_2, ..., f_N\}$ for every supported sensor. Thus, for an input vector $X:=x_{1:N}$, the fusion block output can be given by:
\begin{equation}
    {y} = \mathcal{H}(f_1(x_1), f_2(x_2), ...., f_N(x_N))
    \label{eqn:late_fusion}
\end{equation}
where $\mathcal{H}$ is the fusion algorithm whose inputs are the $N$ outputs $f_i(x_i) \forall i \in N$. When offloading is supported, the goal is to avoid excessive computational overheads. Thus, each model $f_i$ would incorporate an offloading point to be further defined as:
\begin{equation}
f_i(x_i) = f_i^T(f_i^H(x_i)
\end{equation}
where $f_i^H$ and $f_i^T$ are the head and tail parts of the $i_{th}$ model placed prior to and after the offloading point, respectively. The former sub-model is to be deployed locally while the latter is to be replicated across the local and edge server platforms. As server responses could peak due to the wireless channel uncertainty, some model outputs may \textit{not be} \textit{available} for fusion given the strict execution deadlines. Thus, we obtain instead partial fusion outputs given by:
\begin{equation}
    {\hat{y}} = \mathcal{H}(\mathcal{I}_1\times f_1(x_1), \mathcal{I}_2\times f_2(x_2), ...., \mathcal{I}_N \times f_N(x_N))
\end{equation}
where the random variable $\mathcal{I}_i \in \{0,1\}$ indicates whether $f_i(x_i)$ is available for fusion. Naturally, the lesser number of inputs available the more robustness is compromised. Therefore, given $M$ operating offloading modes, the objective is to identify the mode satisfying:
\begin{equation}
    \min\limits_{m\in M} E(X|m), \; s.t. \;\; L(X|m) \leq L_{th}, \;\Delta(y, \hat{y}) \leq C_{th}
    \label{eqn:objective}
\end{equation}
where $E(\cdot)$ and $L(\cdot)$ are the respective end-to-end energy consumption and execution latency for processing the input vector $X$ given processing mode $m$. This formulation is regulated by a latency constraint $L_{th}$ for nominal safety, and a robustness constraint $C_{th}$ based on the difference in prediction quality between $y$ and $\hat{y}$. In the following sections, we present our methodology for solving the optimization objective in (\ref{eqn:objective}), which will entail applying DNN structural optimizations and a runtime learning-based approach. We demonstrate our analysis for the use-case of late fusion following the object detection module as illustrated in Figure \ref{fig:av-pipeline}.

\begin{figure}[!tbp]
{\includegraphics[,width = 0.46\textwidth]{Figures/Object_detection.pdf}}
\vspace{-2ex}
\caption{Object Detection with Late Fusion and offloading support. Blue blocks/variables are passed to/from the DRL domain in section \ref{sec:DRL}. Transparent blocks are inactive.}
\label{fig:object-detection}
\vspace{-3ex}
\end{figure}

\section{Concurrent Pipelines Design \label{sec:design}}

Figure \ref{fig:object-detection} depicts the proposed processing domain for object detection with our applied modifications given concurrent DNNs and a late fusion scheme detailed as follows.

\subsection{Object Detection with Late Fusion}

This scheme entails processing each sensory input separately before aggregating the outputs together through fusion. Specifically, there are two primary computational tasks:

\textit{\textbf{Object Detection:}} For each sensor, an object detection pipeline is implemented to identify and classify objects in a scene. Initially, each model entails a feature extractor based on a Convolutional Neural Network (CNN), e.g., ResNet-18 here \cite{he2016deep}, responsible for abstracting raw sensory data into smaller-sized features for the following detection model, e.g., Faster R-CNN network \cite{ren2015faster}, which consists of a regional proposal network (RPN) to suggest regions of interest where objects may exist, a classification stage to categorize the objects within each proposal, and a final post processing stage to convert classified proposals into bounding box predictions.

\textit{\textbf{Fusion:}} As the outputs from each pipeline are bounding boxes, we can directly fuse them together using Non-Maximum Suppression (NMS) \cite{ren2015faster} to calculate the intersection over union (IoU) and obtain an estimate on the degree of overlapping between each pair of bounding boxes from the overall set of predictions. If an IoU for a pair of boxes exceeds a predetermined threshold, the bounding box with less confidence score is discarded, and this operation repeats until all possible pairing combinations are covered.

\subsection{Implementing DNNs to support Offloading}
To avoid the overhead of offloading raw inputs, we scale the optimal offloading point injection technique in \cite{malawade2021sage, matsubara2020head} to each concurrent pipeline without compromising the overall utility as follows:

\textit{\textbf{Structural Alterations:}} We alter the structure of the feature extractors (ResNet-18 here) within each pipeline to minimize local computation overhead prior to the offloading point and downsize the transmissible data. Specifically, we substitute a considerable portion from the earliest parts of a DNN with an encoder-decoder like structure of two functional components: (\emph{i}) an encoder, $\mathcal{M}_E$, which offers an efficient offloading option at its output through shrinking the input data into a lower-dimensional representation that retains the most relevant of features, i.e., small output sizes translate to low communication overheads., and (\emph{ii}) a decoder, $\mathcal{M}_D$, to receive outputs from the encoder and cast them back to higher dimensional representations of dimensions compatible with the remainder of the network. Here, we replace the first two residual blocks from each ResNet-18 with an encoder-decoder structure.

\begin{figure}[!tbp]
\centering
{\includegraphics[,width = 0.33\textwidth]{Figures/distill.pdf}}
\vspace{-2ex}
\caption{Our ResNet-18 feature extractors undergo Knowledge Distillation to train $\mathcal{M}_E\cdot\mathcal{M}_D$ using the first 2 blocks}
\label{fig:KD}
\vspace{-3.5ex}
\end{figure}

\textit{\textbf{Knowledge Distillation:}} Next, modified architectures need to be retrained to maintain utility. We apply knowledge distillation \cite{matsubara2020head} to train $\mathcal{M}_E$ and $\mathcal{M}_D$ through minimizing a loss function, e.g., mean squared error ($L_{MSE}$), between $\mathcal{M}_D$ outputs and those from the original parts, $\mathcal{M}_{orig}$. Figure \ref{fig:KD} illustrates this for our ResNet-18 with the loss component for a single input $x$ given by:
\begin{equation}
    L_{SE} = || \mathcal{M}_{orig}(x) - \mathcal{M}_D (\mathcal{M}_E(x)) ||_2^{2} \label{eqn:kd}
\end{equation}

Hence, unaltered DNN components can retain their weight values with only the parameters of the new structure trained to produce the same output values as the originals.

\textit{\textbf{Deployment for Inference:}} After retraining, the modified architectures would be deployed for each concurrent detection pipeline on the local ADS. Furthermore, each $\mathcal{M}_D$ and its succeeding blocks would be replicated across the VEC servers to enable online autonomous driving services. Thus during runtime, servers can receive outputs from $\mathcal{M}_E$ components, process them, and return predictions, e.g., bounding boxes coordinates, to the ADS platforms. On the ADS, available local and received predictions are fused to provide the final outputs for the following blocks.

\section{Reinforcement Learning Control} \label{sec:DRL}

VEC operation is reliant on the surrounding conditions with regards to the wireless channel state and the server load. Hence, we propose a learning-based approach based on deep reinforcement learning (DRL) to adapt the mode of operation so as to maximize performance efficiency while maintaining robustness -- which we account for in the offloading decision through leveraging the abstract feature representations already computed within the processing pipelines.

\subsection{Hierarchical Agent}
As shown in Figure \ref{fig:DRL}, our DRL solution constitutes a hierarchical agent whose main components are as follows:

\subsubsection{Contextual Encoder} In order to estimate the complexity of the corresponding scene, we leverage the computed feature set, $\mathcal{F}_{t}$, at time window $t$ from the main sensor processing pipelines to guide the decision for the following window $t+1$, given as $\mathcal{F}_{t} = \{(\mathcal{F}_1)_{t}, (\mathcal{F}_2)_{t}, .., (\mathcal{F}_N)_{t}$\}. The rationale behind using the feature set of the preceding time window is twofold: (\emph{i}) features do not need to be computed from scratch as they have already been generated within the primary processing pipelines (see the global pooling blocks in Figure \ref{fig:object-detection}), and, (\emph{ii}) the small window size for autonomous driving ($\leq$ 100 ms) means that successive frames share similar driving contexts due to the high spatio-temporal scenic correlations.

Given how $\mathcal{F}_{t}$ can outweigh other DRL inputs due to its relatively larger size, $\mathcal{F}_{t}$ needs to be initially encoded into a further lower-dimensional representation. Hence, $\mathcal{F}_{t}$ is processed through a \emph{contextual encoder} comprising a sequence of fully-connected layers to obtain the final abstraction $\mathcal{F}^{*}_{t}$. In our experiments, $\mathcal{F}^{*}_{t}$ was of 256$\times$ smaller in size than $\mathcal{F}_{t}$.

\begin{figure}[!tbp]
{\includegraphics[,width = 0.49\textwidth]{Figures/DRL.pdf}}
\vspace{-6ex}
\caption{Our hierarchical agent for runtime mode selection}
\label{fig:DRL}
\vspace{-3.5ex}
\end{figure}

\subsubsection{State Encoder} The next component is the state encoder whose input is the final state representation $s_{t} = \{\mathcal{F}^{*}_{t}, \phi_{t}, q_t\}$ formed from aggregating the contextual encoder outputs, $\mathcal{F}^{*}_{t}$, the channel capacity $\phi_{t}$, and server queuing delays $q_t$. Practically, the latter two metrics can be estimated by probing the edge server.

\subsubsection{Action Space} Represented by the final fully-connected layer in the \emph{state encoder}, the action space covers the set of all possible modes of operation that can be selected by the DRL at runtime. We define it as $A = \{offload_{0}, offload_{1}, offload_{2}, ... offload_{N-1}\}$, where an action $offload_{i}$ is for choosing the offloading option for $i$ sensory pipelines, with $offload_{0}$ being pure local execution. In the case that the same DNN structure is shared across all pipelines, only the number of offloading pipelines matter. We do not consider $offload_{N}$ as a viable option so that the ADS is always guaranteed a new output every time window since at least one pipeline is always processed locally. This way, even under a worst-case scenario when tasks from $N-1$ pipelines are offloaded and results are not received within the time limit, the vehicle can still operate in a safe manner. In practice, we merely need a subset of actions $A^* \subseteq A$, with \{$offload_{0}, offload_{N-1}\} \subseteq A^*$, where $A^*$ can contain the actions that exhibit notable variability in performance. At runtime, action vector, $a_{t}$, is mapped onto the control of each processing pipeline.

\subsection{DRL Environment}

We detail the emulated DRL training environment for learning a policy $\pi$ that makes offloading decisions based on the current state.

\subsubsection{Training and Reward}

Reinforcement learning approaches rely on having a Q function to provide value estimates for each state-action pair so as to select the optimal action $\hat{a}$ = $\argmax_{a \in A} Q_{\pi}(\hat{s}, a)$ for each $\hat{s}$ under a learnt policy $\pi$. However, estimating state-action pair values in continuous state spaces is challenging, and DRL offers to approximate $Q_{\pi}$ by a policy network trained to maximize a reward. With no loss in generality, our DRL employs a Double Deep Q-Network \cite{van2016deep} with a compounded reward function $\mathcal{R}$ as:
\begin{equation}
    \mathcal{R} =
    \begin{cases}
        \mathcal{A}, \;\;\;\;\;\;\;\;\;\;\; \text{if mAP($y$) $<$ mAP$_{th}$} \\

        \mathcal{B}, \;\;\;\;\;\;\;\;\;\;\; otherwise
    \end{cases}
\end{equation}
which evaluates to different functions based on a measure of robustness, which we associate here with the degree of uncertainty in the final fused predictions $y$ in (\ref{eqn:late_fusion}), determined by the mean Average Precision (mAP) scores for object detectors as in \cite{rahman2021per}. In brief, our goal is for the agent to realize a policy that deters from offloading actions when prediction confidence is low, which we achieve here through leveraging the contextual information in $\mathcal{F}_{t-1}$ to assess the scene's complexity, and make offloading decisions accordingly with the goal of minimizing prediction uncertainty. Thus, if $mAP_{th}$ is not met, $\mathcal{R}$ evaluates to $\mathcal{A}$ defined as:
\begin{equation}
    \mathcal{A} =
    \begin{cases}
        0, \;\;\;\;\;\;\;\;\;\; \text{if $\hat{a} == offload_0$} \\
        \frac{P}{N-i},\;\;\;\;\; \text{if $\hat{a} == offload_i$}
    \end{cases}
    s.t., i\neq 0, i<N
\end{equation}
for penalizing the agent whenever an offloading action is selected, with the penalty value being proportionate to the number of offloading pipelines, $i$, out of $N$ total, with a maximum \textit{negative} penalty of $P$. Recall that $offload_N \not\in A$ as one pipeline always executes locally to ensure  at least one output is available irrespective of the wireless network conditions. On the flip side, when mAP$_{th}$ is satisfied, $\mathcal{R}$ evaluates to $\mathcal{B}$ as follows:
\begin{equation}
    \mathcal{B} =
    \begin{cases}
        P, \;\;\;\;\;\;\;\;\;\;\; \text{if $L(X|\hat{a}) > L_{th}$} \\
        \mathcal{C}, \;\;\;\;\;\;\;\;\;\;\; otherwise
    \end{cases}
\end{equation}
which penalizes the agent by $P$ when its selected action $\hat{a}$ causes the overall execution latency for inputs $X$, $L(X|\hat{a})$, to exceed the critical execution latency constraint, $L_{th}$. In other words, this means that the agent is penalized when not all partial outputs are available in time for late fusion.
In reality, state-of-the-art ADS platforms are designed to meet the application latency demands, and hence, we set the value of $L_{th}$ to that of \textit{local execution}.  Contrarily, when $L_{th}$ is satisfied, $\mathcal{R}$ finally evaluates to $\mathcal{C}$ given by:
\begin{equation*}
    \mathcal{C} =
    \begin{cases}
        0, \;\;\text{if $E(X|\hat{a}) == min(E(X|a)|L(X|a) \leq L_{th})$ } \\
        P, \;\; otherwise
    \end{cases}
\end{equation*}
\begin{equation}
    \forall a \in A^*, A^*\subseteq A
\end{equation}
penalizing the agent by $P$ if the energy consumption footprint $E(X|\hat{a})$ from selecting action $\hat{a}$ is not the minimal from amongst those of all viable actions $a\in A^*$ that are projected to meet $L_{th}$.

\subsubsection{Latency and Energy Estimation \label{subsec:latency_energy}}

In order to evaluate $\mathcal{R}$ for each selected $\hat{a}$, the end-to-end estimates for energy and latency can be approximated every time window as follows:
\begin{align}
    L &= L_{local} + L_{Tx} + L_{server} + L_{Rx} \\
    E &= E_{local} + E_{Tx} + E_{idle} + E_{Rx}
\end{align}
where the latency $L$ can be broken down into the respective local, transmission, server, and receiving latencies. Similarly, energy consumption constitutes the same components except for incorporating idling energy as we are only concerned about the ADS energy footprint. From here, the local components are given by:
\begin{align}
    L_{local} =
        N \times L_{\mathcal{M}_E} + (N-i) \times L_{tail}\mid \text{$\hat{a} == offload_{i}$}
\end{align}
in which $L_{\mathcal{M}_E}$ and $L_{tail}$ are the respective latencies for executing the encoder $\mathcal{M}_E$ and the remaining tail parts of the model, respectively. When the selected action is to offload processing from $i$ processing pipelines (i.e. $\hat{a}$ $== offload_{i}$), the total local execution latency accounts for processing across the $N$ encoders and the $N-i$ tail models. This additive form represents the most direct approach for modeling local execution. However in reality, the concurrency of pipelines can speed up local execution depending on the available hardware resources at the expense of a larger power consumption footprint, $P_{local}$. We approximate this trade-off through considering energy for performance evaluation, defining $E_{local}$ as:
\begin{equation}
   {E}_{local} = L_{local} \times P_{local} \label{eqn:Elocal}
\end{equation}

\subsubsection{Channel Estimation}To estimate the communication overheads, we first fit a Rayleigh distribution curve with scale $\sigma$ to throughput traces $\Phi$ collected from the real-world for different wireless technologies, i.e., $\Phi \sim \text{Rayleigh}(\sigma)$. Then, we use the constructed distribution to sample independent and identically distributed (i.i.d.) random variables as the channel capacity $\phi$ to be used for the training and evaluation processes of the DRL agent where data transmission parameters can be evaluated as:
\begin{equation}
    L_{Tx} =
        \frac{i \times {b}}{\phi}\mid\text{$\hat{a} == offload_i$};\ E_{Tx} = L_{Tx} \times P_{Tx} \label{eqn:tx_latency}
\end{equation}
where $b$ is the transmissible data size from one sensory pipeline while $P_{Tx}$ is the transmission power incurred by the ADS. Similarly, the formulation for the receiving parameters, $L_{Rx}$ and $E_{Rx}$, can be provided given corresponding estimates for channel capacity and data sizes in the downlink.

\subsubsection{Server Queuing\label{eqn:serverQing}}
Lastly, we model the server latency $L_{server}$ using queuing delays where we have:
\begin{equation}
    q_{c} = \frac{(1-\rho)(\rho)^{c}}{1-\rho^{C+1}}
\end{equation}
representing the probability that the offloaded task would encounter $c$ other tasks before it in the server's processing queue, with
$\rho$ being the average server load, and $C$ being the queue size. From here, we are able to generate a probability density function (pdf) for values within 0-C from which we can sample queuing positions, and consequently approximate $L_{server}$.

\section{Experiments and Results}

\subsection{Experimental Setup}

\subsubsection{Dataset}
We use the RADIATE multimodal perception dataset \cite{sheeny2021radiate} for its diverse driving scenarios and adverse weather conditions such as snow, fog, and rain. The dataset covers 8 object classes with annotations from a Navtech CTS350-X radar, a Velodyne HDL-32e LiDAR, and a ZED stereo camera. The variety of scenes provides a varying degree of difficulty for ADS and enables the robustness assessment. For instance, cameras obstructed by snow offer poor visibility indicating higher difficulty that can cause sub-optimal object detection performance. Here, we implemented 4 object detection DNN pipelines: 2 stereo cameras, radar, and lidar. All inputs are mapped onto the forward-facing perspective for late fusion.

\subsubsection{Training and Metrics}
As was mentioned in Section \ref{sec:design}, the original processing pipelines for each sensing modality comprise a ResNet-18 followed by a Faster R-CNN. These models were trained using a batch size of 1, learning rate of 0.005, and the multi-task loss function in \cite{ren2015faster} which combines both classification and box regression losses. For the NMS fusion, we use a fusion IoU threshold of 0.4. We employ mAP as our evaluation metric with boxes IoU $\geq$ 0.5 since it is widely adopted for object detection tasks \cite{everingham2010pascal} where the average precision is estimated using the precision and recall values. More details about evaluating these values are in \cite{everingham2010pascal,ren2015faster}.

\subsubsection{Hardware and Performance Evaluation} We use the industry-grade Nvidia Drive PX2 Autochauffer as our ADS hardware. The concurrent DNN models are compiled using the TensorRT library becoming inference engines. The local execution power $P_{local}$ is estimated as the difference in the ADS power measurements when processing and idling. For the transmission power $P_{Tx}$, we follow \cite{malawade2021sage} and evaluate it using the data transfer power models in \cite{close}.

\subsubsection{Encoder-Decoder Structure} The input frame's resolution for each of the sensory pipelines is 672 $\times$ 376 ($\approx$ 740.25 kB). The encoder, $\mathcal{M}_E$, comprises 3 layers (2 convolutional and 1 pooling), each with a stride of 2 with only 3 channels at the output. Therefore, when the outputs from $\mathcal{M}_E$ are quantized to 8 bits for offloading \cite{malawade2021sage}, the transmissible data size $b$ in equation \ref{eqn:tx_latency} becomes $\approx$ 11.57 kB (64$\times$ less than the input's). The decoder $\mathcal{M}_D$ mimics the structure presented in \cite{matsubara2020head} to have its output of the same dimensions as that from the original second ResNet-18 block.

\subsubsection{DRL Settings} For safety, we always execute the radar pipeline locally \cite{liu2017computer} and define $A^* = \{offload_0, offload_2, offload_3\}$. We set $P=-2$, $C=4000$, $\rho = 0.9$, and mAP$_{th}$ = 0.68 unless otherwise stated. We set $L_{th}=68.12$ based on pure local execution latency.

\begin{table}[ht]
    \centering
    \caption{Loss and mAP (\%) before (orig) and after (dist) integrating $\mathcal{M}_E\cdot \mathcal{M}_D$ across various late fusion combinations.}
    \vspace{-2ex}
    \begin{tabular}{l | c c | c c}
    \hline
    Sensor & Loss (orig) & Loss (dist) & mAP (orig) & mAP (dist) \\
    \hline
    2 Cameras & 0.15 & 0.17 & 67.14 & 67.14 \\

    Radar+Lidar& 0.10 & 0.11 & 67.14 & 67.14\\

    Full Fusion & 0.13 & 0.15 & 71.24 & 70.38 \\

    \hline
    \end{tabular}
    \label{tab:map}
\end{table}
\vspace{-3ex}

\subsection{Object Detection and Performance}
We first assess how the inclusion of $\mathcal{M_E}$ and $\mathcal{M_D}$ impacts the loss and prediction accuracy of object detection. Table \ref{tab:map} shows the changes in these metrics across different late fusion combinations on the Radiate evaluation dataset. As seen, full sensor fusion has the best performance in mAP, asserting how prediction robustness relates to the number of fused outputs. It is also observed that the new DNN structures maintain the same level of performance as their original counterparts, with the highest degradation in mAP from 71.24\% to 70.38\% experienced by the full fusion case, but still offering a better score than that of the simpler sensor combinations.

\begin{table}[ht]
    \centering
    \caption{Hardware Measurements on the Nvidia Drive PX2}
    \vspace{-2ex}
    \begin{tabular}{l c c c c}
    \hline
    DNN & $L_{local}$ (ms) & $E_{local}$ (J) & Memory (MB) \\
    \hline
    Encoder  & 3.78 & 0.03 & 0.025\\

    1 pipeline & 17.03 & 0.12 & 80.3\\

    4 pipelines & 68.12 & 0.48 & 321.2 \\

    DRL Agent & 0.66 & 0.005 & 5.4\\

    \hline
    \end{tabular}
    \label{tab:px2}
\end{table}

Table \ref{tab:px2} displays the processing overheads for different DNN components deployed on the PX2 hardware. The encoder $\mathcal{M_E}$ and DRL agent take 3.78 and 0.66 ms, respectively, emphasizing how the decision $a_t$ is obtained before the generation of any transmissible outputs. Moreover, the execution latency for 4 pipelines on the PX2 can add up to 68.12 ms given the same power $P_{local}$.

\subsection{Channel Capacity and Queuing Analysis}

In this experiment, we analyze the influence of the experienced channel capacity, $\phi$, and queuing delay, $q_t$, on the optimal action choice when optimizing for energy consumption under the latency constraint $L_{th}$. To elaborate, we illustrate in Figure \ref{fig:channel_sweep} parametric sweeps with respect to $\phi$ given $q_t=$ 15 ms. As shown, offloading options are consistently more energy efficient than the pure local option (\textit{offload$_0$}), but the $L_{th}$ constraint dictates which action should be chosen considering how poor values of $\phi$ could disqualify some offloading choices. When $\phi>4$ Mbps, the latency overhead for \textit{offload$_2$} does not exceed $L_{th}$ making it the optimal action until $\phi>7$ Mbps, at which the most energy-efficient option, \textit{offload$_3$}, becomes valid. Similarly, this analysis is repeated in Figure \ref{fig:queue_sweep} when sweeping across $q_t$ given $\phi=8$ Mbps. Naturally, the latency overhead is linearly proportional to $q_t$ under fixed network conditions, demonstrating the influence of server load over the optimal offloading decision. From here, the key takeaway is that based on the wireless infrastructure and VEC server capabilities, the maximum number of concurrent offloading pipelines that meet $L_{th}$ can be determined, and used accordingly to construct the decision space.

\begin{figure}[!tbp]
\centering
{\includegraphics[,width = 0.46\textwidth]{Figures/channel_sweep.pdf}}
\vspace{-2ex}
\caption{Variation of Latency and Energy Analysis w.r.t. Channel Capacity. Energy as bar charts, Latency as plot lines}
\label{fig:channel_sweep}
\vspace{0.5ex}
\end{figure}

\begin{figure}[!tbp]
\centering
{\includegraphics[,width = 0.46\textwidth]{Figures/queue_sweep.pdf}}
\vspace{-2ex}
\caption{Variation of Latency and Energy Analysis w.r.t. Queuing Delay. Energy as bar charts, Latency as plot lines}
\label{fig:queue_sweep}
\vspace{-3.5ex}
\end{figure}

\begin{table}[ht]
    \centering
    \caption{Robustness analysis at $\rho$=0.97}
    \vspace{-2ex}
    \begin{tabular}{l c c c c }
    \hline
    \multirow{2}{*}{Policy} & \multirow{2}{*}{mAP$_{th}$} & \multicolumn{3}{c}{Average mAP (\%)} \\ \cline{3-5} & &  $offload_0$ & $offload_2$ & $offload_3$ \\ \hline

    R-agnostic & N/A & 64.93 & 64.50 & 64.85 \\ \hline
    \multirow{3}{*}{DRL} &
    0.50 & 60.39 & 67.69 & 68.58 \\ \ &
    0.68 & 60.68 & 70.93 & 72.55 \\  &
    0.98 & 61.82 & 72.49 & 73.14 \\ \hline
    \multirow{3}{*}{Oracle} &
    0.50 & 49.58 & 85.34 & 85.26 \\ \ &
    0.68 & 50.58 & 93.68 & 93.62 \\ \ &
    0.98 & 54.84 & 99.51 & 99.49 \\ \hline
    \end{tabular}
    \label{tab:robustness_mAP}
\vspace{-3ex}
\end{table}

\subsection{Robustness Analysis}
To analyze the DRLs capacity to make robust decisions, we define 2 baseline policies for comparison: (\emph{i}) a robustness-agnostic or \textit{R-agnostic} policy that is aware of $\phi$ and $q_t$ to optimize for energy so long as $L_{th}$ is satisfied, and (\emph{ii}) an \textit{Oracle} resembling an optimal strategy which in addition to the information available to the \textit{R-agnostic} policy, also possesses the true per-frame mAP estimate apriori, and consequently, the optimal sequence of decisions satisfying the robustness constraint mAP$_{th}$. All of the mentioned strategies are more energy-efficient than pure local execution, and the mAP$_{th}$ values are set in the experiments to 0.5, 0.68, and 0.98, estimated based on the cumulative distribution of the evaluation dataset such that 30\%, 50\%, and 70\% of the evaluation mAP scores fall under the corresponding thresholds.

To evaluate robustness across each policy, we compute the average experienced mAP per action (AMAP) given the action selection frequencies. Mainly, a robust behavior cause frames of high uncertainty (mAP $\leq$ mAP$_{th}$) to be processed locally, implying how the AMAP experienced locally should be low compared to those from the offloading actions. We illustrate this concept in Figure \ref{fig:action_map} across the 3 policies for mAP$_{th}$ = 0.68 and $\rho$ = 0.97. As seen, the \textit{R-agnostic} policy only considers performance efficiency for its action selection, and subsequently, its AMAP across \textit{offload$_0$}, \textit{offload$_2$}, and \textit{offload$_3$}, are equivalent with values of 64.39\%, 65.01\%, and 65.25\%, respectively. Conversely, the $Oracle$ resembles the ideal embodiment of robustness, assigning high uncertainty frames to \textit{offload$_0$}, despite performance gains from offloading. In contrast to the \textit{R-agnostic} policy, \textit{66.94\%} of the \textit{Oracle} policy decisions are \textit{offload$_0$} with an AMAP of 50.58\%, and an AMAP as high as \textit{93.68\%} for the remaining offloading decisions. From here, our proposed \emph{DRL} approach strives to learn the \textit{Oracle}'s behavior, through the observed action selection breakdown, with 63.65\% of actions belonging to \textit{offload$_0$}. Moreover, the AMAP for \textit{offload$_2$} and \textit{offload$_3$} are 70.93\% and 72.55\%, respectively, which despite outperforming the \textit{R-agnostic} policy, are far from that of the \textit{Oracle}. This is expected considering the \textit{Oracle} policy is the unrealistic ideal behavior with apriori mAP knowledge. We extend this analysis to other thresholds values in Table \ref{tab:robustness_mAP}, where we observe that as the robustness constraint becomes smaller, the \textit{DRL} exhibits a behavior closer to the \textit{R-agnostic} and farther from the \textit{Oracle} and vice versa, indicating the \textit{DRL}'s capacity to adapt to various robustness requirements.
\vspace{-2ex}

\begin{figure}[!tbp]
\centering
{\includegraphics[,width = 0.46\textwidth]{Figures/robustness.pdf}}
\vspace{-2ex}
\caption{Action selection frequencies (\%) breakdown across the 3 policies at mAP$_{th}$=0.68 and $\rho$=0.97. Numbers next to the bars indicate the average experienced mAP (AMAP) for the evaluation dataset inputs mapped to each action.}
\label{fig:action_map}
\vspace{-3.5ex}
\end{figure}

\begin{table}[ht]
    \centering
    \caption{Action frequency analysis at mAP$_{th}$ = 0.68}
    \vspace{-2ex}
    \begin{tabular}{l c c c c}
    \hline
    \multirow{2}{*}{Policy} & \multirow{2}{*}{$\rho$} & \multicolumn{3}{c}{Action Frequency (\%)} \\ \cline{3-5}
    & & $offload_0$ & $offload_2$ & $offload_3$ \\ \hline
    \multirow{3}{*}{R-agnostic}
    & 0.90 & 11.17 & 14.73 & 74.10 \\
    & 0.97 & 37.18 & 17.17 & 45.65 \\
    & 0.99 & 70.34 & 9.25 & 20.41  \\ \hline
    \multirow{3}{*}{DRL}
    & 0.90 & 51.73 & 8.83 & 39.44   \\
    & 0.97 & 63.65 & 10.95 & 25.39  \\
    & 0.99 & 85.14 & 5.46 & 9.4     \\ \hline
    \multirow{3}{*}{Oracle}
    & 0.90 & 53.19 & 8.00 & 38.82 \\
    & 0.97 & 66.94 & 9.17 & 23.89 \\
    & 0.99 & 84.37 & 5.04 & 10.60 \\ \hline
    \end{tabular}
    \label{tab:serverload}
\end{table}

Furthermore, we vary the server load $\rho$ in Table \ref{tab:serverload} and show how the action frequency varies for each policy. As $\rho$ increases, the selection of local processing becomes more frequent across all policies, irrespective of the energy or robustness due to the $L_{th}$ constraint. Such behavior is learned by our $DRL$ solution given how the action selection frequencies closely imitate that of the $Oracle$.
\vspace{-1.5ex}

\begin{table}[ht]
    \centering
    \caption{Energy analysis relative to local at mAP$_{th}$ = 0.98}
    \vspace{-2ex}
    \begin{tabular}{l c c c c}
    \hline
    Metric & Local & R-agnostic & DRL & Oracle \\
    \hline
    Risky Actions (\%) & 0 & 63.37 & \textbf{14.54} & 0 \\
    Robust Actions (\%) & 100 & 36.63 & \textbf{85.46} & 100 \\
    Total Energy (kJ) & 2.916 & \textbf{1.729} & 2.479 & 2.487 \\
    Total Energy Red. (\%) & 0 & \textbf{40.7} & 14.99 & 14.72 \\

    \hline
    \end{tabular}
    \label{tab:energy0.98}
\end{table}

\vspace{-3ex}

\begin{figure}[!tbp]
\centering
{\includegraphics[,width = 0.49\textwidth]{Figures/risk.pdf}}
\vspace{-5ex}
\caption{Comparing Energy and Risky Actions}
\label{fig:energy_risk}
\vspace{-4ex}
\end{figure}

\subsection{Energy Reduction vs Risky Actions}
We also compare the energy savings relative to the pure local execution, \textit{offload$_0$}, in addition to their risky behaviors. We first define \textit{Risky Actions} as the fraction of offloading actions whose respective mAP scores fall below mAP$_{th}$, and \textit{Robust Actions} as the fraction whose scores exceed the mAP$_{th}$. We compare the performance of each policy in Table \ref{tab:energy0.98}, where although \textit{R-agnostic} offers the highest energy reduction of 40.7\% compared to \textit{DRL}'s 14.99\%, 63.37\% of \textit{R-agnostic}’s energy savings are \textit{Risky Actions}, unlike \textit{DRL} whose \textit{Risky Actions} constitute 14.54\% of the offloading decisions. Through extending this analysis further to entail multiple mAP$_{th}$ values, i.e., a higher threshold means a stricter offloading constraint, we observe in Figure \ref{fig:energy_risk} that the robustness-aware \textit{DRL} at higher mAP$_{th}$ substantially reduces the percentage of risky offloads compared to the \textit{R-agnostic} policy, with reductions of 24.50\% and 77.06\%, at mAP$_{th}=0.50$ and mAP$_{th}=0.98$ respectively.

\section{Conclusion}
In this work, we presented \textsc{Romanus}, a methodology for robust and efficient task offloading for multi-sensor autonomous driving systems (ADS). We first showed how to integrate optimal offloading points along the processing pipelines in a multi-sensor object detection module with late fusion, and then implemented a DRL-based runtime offloading that achieves 14.99\% energy efficiency over pure local execution with up to 77.06\% decrease in risky offloading actions from a robust-agnostic solution. This methodology can be generalized to a variety of sensors and fusion strategies depending on the underlying system structure and how robustness is characterized with regards to the application's primary task, which would form the basis for future research works along this direction.

\bibliographystyle{ACM-Reference-Format}