Eigenstate phases with finite on-site non-Abelian symmetry

TL;DR

This paper investigates eigenstate phases in disordered spin chains with non-Abelian symmetry, developing a formalism for invariant Hamiltonians and identifying thermal and many-body localized phases with symmetry breaking.

Contribution

It introduces a general group-theoretic framework for constructing symmetric Hamiltonians and applies it to study eigenstate phases with $S_3$ symmetry, revealing distinct thermal and MBL phases.

Findings

Identification of thermal and MBL phases with $S_3$ symmetry

Detection of spontaneous symmetry breaking from $S_3$ to $ ext{Z}_3$

Evidence of an intermediate phase with full $S_3$ symmetry

Abstract

We study the eigenstate phases of disordered spin chains with on-site finite non-Abelian symmetry. We develop a general formalism based on standard group theory to construct local spin Hamiltonians invariant under any on-site symmetry. We then specialize to the case of the simplest non-Abelian group, $S_3$, and numerically study a particular two parameter spin-1 Hamiltonian. We observe a thermal phase and a many-body localized phase with a spontaneous symmetry breaking (SSB) from $S_3$ to $\mathbb{Z}_3$ in our model Hamiltonian. We diagnose these phases using full entanglement distributions and level statistics. We also use a spin-glass diagnostic specialized to detect spontaneous breaking of the $S_3$ symmetry down to $\mathbb{Z}_3$. Our observed phases are consistent with the possibilities outlined by Potter and Vasseur [Phys. Rev. B 94, 224206 (2016)], namely thermal/ ergodic and spin-glass many-body localized (MBL) phases. We also speculate about the nature of an intermediate region between the thermal and MBL+SSB regions where full $S_3$ symmetry exists.

Equations

$$S^x$$

$$S^y$$

$$V(a)=\left(\begin{array}{ccc}\omega & 0& 0\\ 0& 1& 0\\ 0& 0& {\omega }^{\ast }\end{array}\right),V(x)=\left(\begin{array}{ccc}0& 0& 1\\ 0& -1& 0\\ 1& 0& 0\end{array}\right),$$

$$V(a)S^{\pm }V(a{)}^{†}$$

$$V(x)S^{\pm }V(x{)}^{†}$$

$$H(\lambda ,\kappa )$$

$$H_d(\kappa )$$

$$H_t$$

$$H_a$$

$$H_b$$

$$H_c$$

$$U(g)H(\lambda ,\kappa )U(g{)}^{†}=H(\lambda ,\kappa ),U(g)=\underset{i=1}{\overset{L}{⨂}}V(g{)}_i.$$

$$\chi =\frac{1}{N}\sum _{i,j=1}^N\left[⟨s_i{s}_j{⟩}^2\right],$$

$${\chi }^{SG}=\frac{1}{L}\sum _{i,j=1}^L|⟨ϵ|{\sigma }_i^z{\sigma }_j^z|ϵ⟩{|}^2,$$

$${\chi }_{{\mathbb{Z}}_3}^{SG}=\frac{1}{L-1}\sum _{i\ne j=1}^L|⟨ϵ|S_i^z{S}_j^z|ϵ⟩{|}^2.$$

$$r_n=\min ({\delta }_n,{\delta }_{n+1})/\max ({\delta }_n,{\delta }_{n+1}).$$

$$g(\lambda ,L)=L^{a}f((\lambda -{\lambda }_c)L^{\frac{1}{\nu }})$$

$$\begin{gathered} H=\sum _{i=1}^L\sum _{\alpha =1}^{D}E[i{]}_{\alpha }\hat{P}[i{]}_{\alpha } \\ +\sum _{i\ne j=1}^L\sum _{\alpha ,\beta =1}^{D}E[i,j{]}_{\alpha ,\beta }\hat{P}[i{]}_{\alpha }\hat{P}[j{]}_{\beta } \\ +\sum _{i\ne j\ne k=1}^L\sum _{\alpha ,\beta ,\gamma =1}^{D}E[i,j,k{]}_{\alpha ,\beta ,\gamma }\hat{P}[i{]}_{\alpha }\hat{P}[j{]}_{\beta }\hat{P}[k{]}_{\gamma }+\dots \end{gathered}$$

$$g:H\to U(g)H{U}^{†}(g)=H$$

$${\hat{P}}_{\Gamma ,m_{\Gamma };d_{\Gamma }}=|\Gamma ,m_{\Gamma };d_{\Gamma }⟩⟨\Gamma ,m_{\Gamma };d_{\Gamma }|$$

$$g:|\Gamma ,m_{\Gamma };d_{\Gamma }⟩\to \sum _{n_{\Gamma }=1}^{|\Gamma |}\Gamma (g{)}_{n_{\Gamma },m_{\Gamma }}|\Gamma ,n_{\Gamma };d_{\Gamma }⟩$$

$$\begin{gathered} H=\sum _{i=1}^L\sum _{\Gamma =1}^{N_R}\sum _{d_{\Gamma }^i=1}^{D_{\Gamma }^i}E[i{]}_{\Gamma ,d_{\Gamma }^i}\hat{P}[i{]}_{\Gamma ,d_{\Gamma }^i} \\ +\sum _{i\ne j=1}^L\sum _{\Gamma ,{\Gamma }^{\prime }=1}^{N_R}\sum _{d_{\Gamma }^i=1}^{D_{\Gamma }^i}\sum _{d_{{\Gamma }^{\prime }}^j=1}^{D_{{\Gamma }^{\prime }}^j}E[i,j{]}_{\Gamma ,d_{\Gamma }^i,{\Gamma }^{\prime },d_{{\Gamma }^{\prime }}^j}\hat{P}[i{]}_{\Gamma ,d_{\Gamma }^i}\hat{P}[j{]}_{{\Gamma }^{\prime },d_{{\Gamma }^{\prime }}^j} \\ +\dots \end{gathered}$$

$${\hat{P}}_{\Gamma ;d_{\Gamma }}\equiv \sum _{m_{\Gamma }=1}^{|\Gamma |}|\Gamma ,m_{\Gamma };d_{\Gamma }⟩⟨\Gamma ,m_{\Gamma };d_{\Gamma }|.$$

$$H|\{\Gamma ,m_{\Gamma };d_{\Gamma }\}⟩=\mathcal{E}(\{\Gamma ;d_{\Gamma }\})|\{\Gamma ;d_{\Gamma }\}⟩,$$

$$\mathcal{E}(\{\Gamma ;d_{\Gamma }\})=\sum _{i=1}^{L}E[i{]}_{{\Gamma }_i,d_{{\Gamma }_i}^i}+\sum _{i\ne j=1}^{L}E[i,j{]}_{{\Gamma }_i,d_{{\Gamma }_i}^i,{\Gamma }_j^{\prime },d_{{\Gamma }_j^{\prime }}^j}+\dots$$

$$|\{\Gamma ,m_{\Gamma };d_{\Gamma }\}⟩=\underset{i=1}{\overset{L}{⨂}}|{\Gamma }_i,m_{{\Gamma }_i};d_{{\Gamma }_i}^i⟩.$$

$$⟨a,x|a^3=x^2=1,xax=a^{-1}⟩.$$

$$V(g)=\left(\begin{array}{ccc}{\Gamma }^2(g{)}_{11}& 0& {\Gamma }^2(g{)}_{12}\\ 0& {\chi }^{1^{\prime }}(g)& 0\\ {\Gamma }^2(g{)}_{21}& 0& {\Gamma }^2(g{)}_{22}.\end{array}\right)$$

$$V(a)=\left(\begin{array}{ccc}\omega & 0& 0\\ 0& 1& 0\\ 0& 0& {\omega }^{\ast }\end{array}\right),V(x)=\left(\begin{array}{ccc}0& 0& 1\\ 0& -1& 0\\ 1& 0& 0\end{array}\right).$$

$$\hat{\lambda }=\left(\begin{array}{ccc}{\lambda }_{11}& {\lambda }_{12}& {\lambda }_{13}\\ {\lambda }_{21}& {\lambda }_{22}& {\lambda }_{23}\\ {\lambda }_{31}& {\lambda }_{32}& {\lambda }_{33}\end{array}\right)$$

$$V(g)\hat{\lambda }{V}^{†}(g)=\hat{\lambda }$$

$$\begin{gathered} \hat{\lambda }=\left(\begin{array}{ccc}{\lambda }_2& 0& 0\\ 0& {\lambda }_{1^{\prime }}& 0\\ 0& 0& {\lambda }_2\end{array}\right)={\lambda }_{1^{\prime }}|0⟩⟨0|+{\lambda }_2\left(|+⟩⟨+|+|-⟩⟨-|\right) \\ =({\lambda }_2-{\lambda }_{1^{\prime }})(S^z{)}^2+{\lambda }_{1^{\prime }}{1}_3 \end{gathered}$$

$$\hat{J}=\left(\begin{array}{cccc}{J}_{11}& J_{12}& \dots & J_{19}\\ J_{21}& J_{22}& \dots & J_{29}\\ ⋮& & \ddots & ⋮\\ J_{91}& J_{92}& \dots & J_{99}.\end{array}\right)$$

$$V(g)\otimes V(g)\hat{J}{V}^{†}(g)\otimes V^{†}(g)=\hat{J}$$

$$1^{\prime }\otimes 1^{\prime }$$

Full text

I Introduction

Statistical mechanics and thermodynamics are bridges that connect microscopic laws such as Newtonian and quantum mechanics to macroscopic phenomena that we measure in the laboratory. The validity of statistical physics relies on the existence of thermal equilibrium. For an isolated quantum system, the notion of thermalization is understood in the form of the Eigenstate Thermalization hypothesis (ETH) 1991_Deutsch_ETH_PRA ; 1994_Srednicki_ETH_PRE . ETH posits that for a quantum system, an eigenstate embodies an ensemble and thermalization can be diagnosed by monitoring if subsystems are thermal with respect to the rest of the system. Furthermore, if the system thermalizes, all eigenstates are thermal. Integrable systems violate ETH due to the existence of an extensive number of conserved quantities that prevent the system from acting as a bath for itself. However, quantum integrable models are highly fine tuned and one recovers thermalization by any infinitesimal deviation from the integrable point.

Recently, many-body localization (MBL) has emerged as a generic class of interacting and disordered isolated systems which violate ETH. Basko et al. Basko06 showed that all many-body eigenstates remain localized to all orders in perturbation for an effective interacting disordered model. Several numerical works subsequently verified that all many-body eigenstates are localized in one dimensional disordered lattice models with short-range interactions oganesyan2007 ; pal2010 ; bauerNayak2013completeset ; vadim ; bardarson2014 ; huse2015many . Furthermore, there has been a mathematical proof by Imbrie imbrie2014many for the existence of MBL in a particular disordered spin model with short range interactions. The absence of thermalization has been further quantified as a consequence of emergent integrability due to the presence of a complete set of local integrals of motion (LIOMs). The key distinction from fine tuned integrable models is that the LIOMs or ‘lbits’ (for localized bits) in the MBL phase are robust against perturbations. One can use these lbits to construct a phenomenological lbit Hamiltonian that captures the entanglement dynamics serbyn2013 ; chandran2015 ; huseNandkishoreOganesyan2014phenomenology ; ros2015integrals .

Having established the existence of MBL and its violation of ETH in certain models, natural questions that arise are “what is the most robust version of MBL?” and “does it lead to a refined notion of ETH?” To this end, it is worthwhile to consider instabilities to the MBL phase that lead to delocalization and thereby the restoration of thermalization. Recent works have considered instabilities to the MBL phase due to a small bath johri2015many ; hyatt2017many , external drive lazarides2015fate , Griffiths effects and dimensionality de2016absence ; de2017many , topologically protected chiral edge nandkishore2014marginal and a single particle mobility edge li2015many ; modak2015many . Contrary to the common wisdom that these instabilities would lead to the complete restoration of thermalization, preliminary numerical results have indicated that the lack of thermalization tends to survive in some form in all these cases. However, the fate of these exotic phases in the thermodynamic limit is still an open question.

Potter and Vasseur potter2016symmetry have recently added another instability to this list. It was argued that the l-bit Hamiltonian ‘enriched’ with non-Abelian symmetries is unstable to perturbations. The physical intuition of the delocalization and thermalization stems from the resonant splitting of extensive degeneracy in the many-body spectrum associated with the higher-dimensional irreducible representations (irreps) of the non-Abelian group. Any perturbation results in either resonant delocalization of the many-body spectrum and thermalization, a quantum critical glass (QCG)-like phase or a spin-glass (SG) phase by spontaneous symmetry breaking (SSB) to an Abelian subgroup.

Thus symmetry provides a platform to search for exotic violations of ETH beyond MBL in strongly interacting systems. To this end, in this paper, we develop a procedure to construct general Hamiltonians with global symmetries and analyze the thermalization and localization indicators for individual eigenstates. As a particular example, we construct a two-parameter Hamiltonian with on-site $S_{3}$ symmetry. Using numerical exact diagonalization of our model Hamiltonian, we calculate cut-averaged entanglement entropy (CAEE) distributions and level statistics which are indicators of localization and a spin-glass diagnostic which detects symmetry breaking. Within the accuracy of our numerical analysis, we are able to distinctly observe both a thermal phase and an MBL spin-glass (SG) phase with spontaneous symmetry breaking of $S_{3}$ to $\mathbb{Z}_{3}$ symmetry.

We also employ the same diagnostics to quantify an intermediate region between the aforementioned two phases where the full $S_{3}$ symmetry is intact. However, we cannot ascertain the fate of this region in thermodynamic limit, due to the possibility of quantum critical cone like finite-size effects husekhemani2016critical . The paper is organized as follows. In Sec. II we construct a general $S_{3}$ symmetric Hamiltonian using group theory methods. We numerically diagonalize our model and compute indicators of localization and symmetry breaking in  III. We end the paper with discussion of results and conclusion. We provide a review of the conjecture by Potter and Vasseur potter2016symmetry on the incompatibility of MBL with non-Abelian symmetries and other details of our analysis in the Appendices.

II Model $S_{3}$ invariant Hamiltonian

In this section, we analyze a specific spin-1 Hamiltonian that is invariant under the smallest non-Abelian group, $S_{3}$. In terms of the spin angular momentum basis $|S=3,S_{z}=+1,0,-1\rangle$, the spin operators are

[TABLE]

The symmetry group $S_{3}$ contains six elements: $S_{3}=\{1,a,a^{2},x,xa,xa^{2}\}$, and the two generators, $a$ and $x$, satisfy the properties $a^{3}=x^{2}=1$ and $xax=a^{-1}$. Note that in this paper, we refer to the identity element of the group simply as $1$. In the spin basis, they are chosen to have the following representations,

[TABLE]

where $\omega=e^{2\pi i/3}$. It can be verified that the spin operators transform under the generators as follows,

[TABLE]

Using the symmetry arguments detailed in Appendix B, we construct the following Hamiltonian:

[TABLE]

The above Hamiltonian consists of two parts: 1) The disordered part $H_{d}$ with a one body (disordered $h_{i}$ term) and a two-body l-bit term (disordered $J_{i}$ term). The two-body term is designed to drive spontaneous symmetry breaking of non-Abelian $S_{3}$ symmetry down to an Abelian $\mathbb{Z}_{3}$ symmetry. The relative strengths of these two terms are controlled by the $\kappa\in[0,1]$ parameter, where $\kappa=1$ is expected to be the SSB limit. 2) The second term $H_{t}$, the thermalizing term, contains a representative subset of the most general two-body symmetric operators. The intention is to keep $H_{t}$ sufficiently generic while retaining invariance under symmetry action (see Appendix B for details of how this Hamiltonian is constructed and can be generalized to arbitrary symmetry groups). The $\lambda$ parameter controls disorder strength.

Using the transformation of the spin operators listed above, it is straightforward to verify that the Hamiltonian has the desired symmetry, that is, $\forall g\in S_{3}$

[TABLE]

For numerical analysis, the parameters in the Hamiltonian are selected as follows, (a) $h_{i}=w_{h}~{}g_{h}(i)$, $J_{i}=w_{J}~{}g_{J}(i)$ and $g_{h/J}(i)$ are random numbers drawn from a normal distribution with mean 0 and standard deviation 1; (b) $(w_{h},w_{J},\Delta_{t},a,b,c)$ are free real parameters. We arbitrarily fix these to the values $(1.0,0.6,0.17,0.74,0.67,0.85)$ respectively for our numerical study without loss of generality.

III Numerical Results

We perform exact diagonalization of the Hamiltonian in Eq. 4. The local Hilbert space for the $S_{3}$ symmetric Hamiltonian is three dimensional, in contrast to the spin $\frac{1}{2}$ case and this constraints our maximum system size accessible to be $L=10$ sites. We study the properties of the eigenstates pertaining to localization and thermalization of this Hamiltonian for various values of $\lambda\in(0,\infty)$ and $\kappa\in[0,1]$. For clarity of presentation of certain analysis, we use the rescaled variable $\frac{\lambda}{1+\lambda}\in(0,1)$ instead of $\lambda$ (wherever mentioned). We employ periodic boundary conditions for all our analysis. To characterize the phases, we study below several relevant diagnostics that quantify the nature of localization and thermalization of the eigenstates. First, we consider the full entanglement distributions evaluated using the cut averaged entanglement entropy.

III.1 Cut averaged entanglement entropy distributions

Since MBL is a characteristic of a single eigenstate, it is useful to quantify this phase without averaging across different eigenstates. Recent work by Yu-Clark and Pekker clark2016cutavg has proposed cut averaged entanglement entropy (CAEE) to quantify the MBL phase at the level of a single eigenstate. CAEE is obtained by averaging the entanglement entropy for a given subsystem size across all possible locations of the subsystem on the spin chain. The CAEE scaling is then evaluated by repeating this procedure for different subsystem sizes. The key advantage of the CAEE is that strong sub-additivity condition constrains the shape of the entanglement scaling as a function of subsystem size, $i.e.~{}$, CAEE is guaranteed to be a smooth convex function of the subsystem size without any average over disorder or eigenstates clark2016cutavg . This allows us to quantify the entanglement scaling of each eigenstate using the slope of the CAEE (SCAEE), as in Ref clark2016cutavg . We can then construct the full distribution of the slopes across the disorder snapshots and eigenstates. Fig. 1 shows sample CAEE, along with a spline fit for 200 randomly chosen eigenstates for 4 different $\{\lambda,\kappa\}$ from a few disorder realizations of a 10-site Hamiltonian [4].

It can be seen that the eigenstates for small $\lambda$ are mostly volume law, while for large $\lambda$, there exists mixture of area-law and volume-law states. In order to identify the nature of the eigenstate transitions for various parameter regimes, we monitor the full distribution of the slope of the CAEE (SCAEE) evaluated at subsystem size $L/4$ for different values of $\lambda$ and $\kappa$ and disorder realizations. We compute SCAEE as follows: for each eigenstate, we compute the CAEE, $S(d)$ for different subsystem sizes, fit the data to a curve using a spline fit and then evaluate the slope, $S^{\prime}(d)$ for this fit curve at $d=L/4$.

There is however a potential issue because of the non-Abelian nature of the $S_{3}$ symmetry of the Hamiltonian. Generally, at finite-sizes, the eigenstates of a Hamiltonian invariant under the action of a symmetry group $G$ transform as irreps of the same group. For our case, $S_{3}$ has 3 irreps- two 1D irreps ($\mathbf{1}$, $\mathbf{1^{\prime}}$) and one 2D irrep ($\mathbf{2}$). Eigenstates that transform as the $\mathbf{2}$ are two-fold degenerate. We may get different entanglement scaling depending on which precise orthonormal states in this 2D vector space is produced as the eigenstates 111We thank Bela Bauer for pointing this out. To avoid this issue, we diagonalize the Hamiltonian in the 1D irrep sector. In addition to ensuring that we only sample non-degenerate eigenstates which transform as 1D irreps, this also helps us in reaching higher system sizes (See Appendix [C] for more details).

Fig. 2 shows the distribution of $S^{\prime}(L/4)$ for different $\lambda$ and $\kappa$. The values of $\lambda$ are chosen so as to show what the distribution looks like when the system has strong, weak and intermediate disorder strength. For all values of the SSB parameter $\kappa$, for weak disorder $\lambda$, we see that the $S^{\prime}(L/4)$ distribution becomes increasingly narrow with system size with a peak located close to $1.1\approx log(3)$ which is the maximum value possible for $S^{\prime}(d)$ for any state of spin-1 chain. This is also evident from the large number of eigenstates with volume-law entanglement entropy scaling $S(d)\propto d$ in Fig. 1. These properties are consistent with an ergodic/ thermal phase. For high disorder however, we see that the distribution for $\kappa=0$ in Fig. 2 is different from $\kappa=1$. For the former, there is a relatively extended thermal tail which is suppressed for the latter. To gain a better understanding we present the first two moments (mean and variance) of this slope distribution which are indicators of a potential MBL transition.

$\kappa=1$* with $S_{3}\rightarrow\mathbb{Z}_{3}$ symmetry breaking:* The entanglement distribution for the $\kappa=1$ limit Fig. 2 and its moments, displayed in the upper panels of Figs. 3,4 are consistent with the existence of an MBL phase for the large disorder limit and transitions to a fully thermal phase for weak disorder. One way to estimate the transition point is by locating where the mean, $\overline{S^{\prime}(L/4)}$ curves for different system sizes cross on the $\lambda$ axis in Fig. 3. This is roughly at $\lambda/(1+\lambda)\approx 0.72$. Another is by locating the peak of the variance, $\sigma^{2}(S^{\prime}(L/4))$ curve on the $\lambda$ axis which is believed to be close to the point of phase transition in the thermodynamic limit husekhemani2016critical ; husezhang2016floquet . The drift of this point towards the disorder side $i.e.~{}$larger $\lambda$ with increase in system size is considered to be typical for exact diagonalization (ED) studies of MBL husekhemani2016critical . Since our model has a non-Abelian symmetry, the existence of a full MBL phase must accompany SSB to an Abelian subgroup. We confirm SSB ($S_{3}$ to $\mathbb{Z}_{3}$ in this case) by computing a spin-glass diagnostic in III.2.

$\kappa=0$* with full $S_{3}$ symmetry:* The entanglement distribution for the $\kappa=0$ limit, shown in Fig. 2, and its moments, displayed in the lower panels of Figs. 3 and 4 shows an enhanced variance and mean at the $\kappa=0$ for the large disorder limit. The enhanced mean value is an indication of the presence of sub-thermal volume law states and area law states. However, the crossing of the $\overline{S^{\prime}(L/4)}$ curves persists (roughly at $\lambda/(1+\lambda)\approx 0.70$) as does the peak in the $\sigma^{2}(S^{\prime}(L/4))$ plot. How this peak value changes as we approach the thermodynamic limit is an open question and we hope that better tools of numerical analysis like matrix product state methods can shed some light on this issue. We leave this for future work.

III.2 Spontaneous symmetry breaking in excited states

As indicated by the entanglement distributions, the full MBL phase appears only in the $\kappa=1$ limit. This is the limit where the disordered ‘l-bit’ term is dominated by $\sum_{i}J_{i}~{}S^{z}_{i}S^{z}_{i+1}$ which triggers the spontaneous symmetry breaking (SSB). To confirm that SSB has indeed taken place for the many-body excited states, we use a spin-glass diagnostic which we describe below. In the study of classical spin glasses (SG) Edwards_Anderson ; spinglass_binder , one is interested in order parameters sensitive to the spin-glass phase characterized by long-range order in the presence of disorder. One such important quantity of study is the spin-glass susceptibility spinglass_nishimori

[TABLE]

where, $s_{i}$ are classical Ising variables, $\langle*\rangle$ indicates statistical averaging and $\left[*\right]$ indicates disorder averaging. In Ref bardarson2014, , the authors defined a similar quantum mechanical diagnostic to detect spin-glass order arising from SSB of $\mathbb{Z}_{2}\rightarrow~{}trivial~{}group$ in a $\mathbb{Z}_{2}$ invariant Ising- like disordered spin chain:

[TABLE]

where, $|\epsilon\rangle$ is an energy eigenstate and $\sigma^{z}$ is the Pauli-Z operator. In their model, it was shown that the average $\overline{\chi^{SG}}$ scales with system size as $\overline{\chi^{SG}}\sim L$ in the SG phase and approaches a constant value set by normalization for the paramagnetic phase. Similar to eq. 7, for our model, we define the following spin-glass diagnostic that looks for signatures of SG order arising from SSB of $S_{3}\rightarrow\mathbb{Z}_{3}$

[TABLE]

Note that we choose to exclude the $i=j$ term in the summation unlike Eq [7]. We look at the statistics of $\chi^{SG}_{\mathbb{Z}_{3}}$ across randomly sampled eigenstates and disorder realizations. We find signatures for transition to an SG phase as we vary $\lambda$ similar to what was found in Ref bardarson2014, . Fig. 5 shows $\overline{\chi^{SG}_{\mathbb{Z}_{3}}}$ versus $\lambda$ for $\kappa=0,1$.

For $\kappa=1$, we indeed observe that the SG diagnostic increases with system size for large disorder. For $\kappa=0$, we see that the SG diagnostic at best saturates to a constant value independent of system size and at worst reduces with system size but certainly does not increase with it indicating the lack of SSB. To rule out SSB to other subgroups in the $\kappa=0$ regime, we have to use other SG diagnostics that can detect spontaneous breaking of $S_{3}$ down to one of the other subgroups like $\mathbb{Z}_{2}$ and trivial group. In Appendix D, we construct appropriate SG diagnostics and present numerical evidence that the full $S_{3}$ is indeed intact for $\kappa=0$.

III.3 Level statistics

Level statistics is a basis independent diagnostic that indicates localization and thermalization based on the statistics of the adjacent gap ratio ($\delta_{n}=E_{n+1}-E_{n}$) defined as,

[TABLE]

In the presence of symmetry, no interaction term in the Hamiltonian can mix eigenstates that transform as different irreducible representations of the symmetry group. Thus, it is meaningful to compute $r_{n}(\Gamma)$, the level statistics ratio for the eigenstates that transform as each irrep $\Gamma$ of the group separately. In this paper, in Fig. 6, we present $\bar{r}=\frac{\overline{r}(\mathbf{1})+\overline{r}({\mathbf{1^{\prime}}})}{2}$ as a function of $\lambda$ where $\overline{r}(\Gamma)$ is obtained by averaging $r_{n}(\Gamma)$ computed for randomly sampled eigenstates which transform as the $\Gamma$ irrep (see Appendix C for details on how we detect the irrep of each eigenstate), across disorder realizations. On the side with low disorder, for both $\kappa=0$ and $\kappa=1$, we observe $\bar{r}\approx 0.53$. This is typical for a fully thermal phase indicating a Wigner-Dyson (WD) distribution. For a fully localized phase, typically, $\bar{r}\approx 0.39$, corresponding to Poisson distribution. However, the application of this diagnostic to the degenerate spectrum of our $S_{3}$ invariant Hamiltonian can be tricky. In particular, closer to the l-bit point of large $\lambda$, the extensively degenerate spectrum may mimic level clustering and result in $\bar{r}<0.39$ which we do indeed observe. Hence, this diagnostic is only reliable at the thermal side, where we obtain the WD value of $\bar{r}\approx 0.53$. For the strong disorder, the level statistics approaches Poisson value but a clear reading is plagued by the degeneracies for both $\kappa=0,1$.

III.4 Finite size scaling

In order to determine the location and nature of the putative transitions, we perform finite size scaling collapse for MBL and SG diagnosics using the following scaling ansatz used in Ref bardarson2014 .

[TABLE]

Fig. 7 shows the scaling collapse of the diagnostics $\overline{S^{\prime}(L/4)},~{}\sigma^{2}(S^{\prime}(L/4))$ and $\overline{\chi^{SG}_{\mathbb{Z}_{3}}}$ for 8, 9 and 10 site data. For $\kappa=1$, we obtain a good collapse for both the MBL and SG diagnostics at $\lambda_{c}\approx 2.65$ and $\nu\approx 2.5$. For $\kappa=0$, we do not get a scaling collapse for the SG diagnostic $\chi^{SG}_{\mathbb{Z}_{3}}$ for non-zero $\lambda_{c}$ and positive $\nu$ which is consistent with the absence of SSB. On the other hand, we get a reasonable collapse for $\overline{S^{\prime}(L/4)}$ and $\sigma^{2}(S^{\prime}(L/4))$ at $\lambda_{c}\approx 2.35$ and $\nu\approx 2.5$.

Note that the value of the finite size exponent $\nu$ used in Fig. 7 is consistent with the Harris Harris /CCFS CCFS1 ; CCFS2 /CLO CLO criterion of $\nu\geq 2$ for one dimensional spin chains with quenched disorder. However, we emphasize that the our estimate is very rough and even for values of $\nu$ that violates the criterion (eg: $\nu\approx 1.5$ ), we still get a decent collapse. This violation is a common feature of ED studies of small system sizes Louitz_nuviolation ; bardarson2014 . Indeed, recent work by Khemani and Huse khemaniTwoUniversality suggests that this might be an indication of the system not exhibiting true thermodynamic behavior and is expected to undergo a crossover after which we obtain $\nu$ consistent with the Harris/CCFS/CLO bounds. Within the accuracy of our numerical investigation however, we can neither confirm nor rule out the possibility of our system being en-route to such a crossover.

IV Discussion

Based on the above analysis, we are in a position to put together an approximate phase diagram for the model Hamiltonian $H$. Fig. 8 shows a color map of $\overline{S^{\prime}(L/4)}$, $\overline{\chi^{SG}_{\mathbb{Z}_{3}}}$ and $\overline{r}$ plotted in the $\kappa,\lambda$ space. Fig. 9 shows a schematic plot that indicates a thermal phase for the weak disorder limit and a MBL+SSB phase for strong disorder limit. In addition to these two phases, there is a hint of a third regime for strong disorder and no SSB, where there seems to be a coexistence of localized and delocalized states in the many-body spectrum. We now discuss each of this regimes separately.

Thermal phase. For the thermal phase, the distribution of slopes has a mean that saturates at the maximal entropy per unit site. This indicates a substantial presence of volume law scaling eigenstates. The distribution of $\chi^{SG}_{\mathbb{Z}_{3}}$ has mean that does not increase with system size and implies absence of SSB. The level statistics clearly show Wigner-Dyson distribution highlighting the thermal nature of this regime.

MBL+SSB phase. For the MBL+ SSB phase, the distribution of slopes has mean close to 0. This indicates the presence of area-law scaling eigenstates. The distribution of $\chi^{SG}_{\mathbb{Z}_{3}}$ parameter (designed to detect $S_{3}\rightarrow\mathbb{Z}_{3}$ symmetry breaking) increases with system size. This indicates a substantial presence of eigenstates with SSB. The above results are consistent with a spin-glass phase with a residual Abelian symmetry group that supports a full MBL states with strong signatures of SSB at all eigenstates. The level statistics data is more noisy due to the level clustering as a result of degeneracies in this regime but hovers around the Poisson value.

Intermediate phase. This case is shown as dashed region in the schematic Fig 9. The full $S_{3}$ symmetry is intact in this regime and is incompatible with the full MBL phase. At the l-bit point ($\lambda\rightarrow\infty$), the many-body spectrum has extensive number of states with extensive degeneracy. Any infinitesimal thermal perturbation is expected to split this extensive degeneracy down to the minimal values set by the size of the irreps of the group (2 and 1 in our case). This is expected to result in resonant delocalization of an extensive number of states. Thus the most plausible scenario is a fully thermal phase. Our numerical results, however, are not fully consistent with a thermal phase, but do not have the strong signatures of SG or MBL phases, either. A possibility for this regime is a marginal MBL phase that may be separating another MBL+SSB phase that is not explored within our parameter space. This can be uncovered introducing disorder to some of the relevant two-body thermal terms which are non-disordered in our current analysis. We leave this investigation for future work.

Finite-size effects such as the critical-cone region husekhemani2016critical are important in this putative intermediate phase. Indeed, such finite-size effects play an important role at large $\lambda$, and one has to be careful about the order of the thermodynamic and large $\lambda$ limits. At finite-size $L$, large $\lambda$ corresponds to adding a small perturbation of $H_{t}$ on top of the lbit Hamiltonian, and for $1/\lambda$ much smaller than the lbit many-body level spacing $\sim 3^{-L}$ one can apply first-order perturbation theory. One effect is to split the exponentially degenerate lbit states, resulting presumably in volume law eigenstates. However, there are also states in the lbit spectrum which are not highly degenerate, corresponding to sectors where most of the sites sit in the one dimensional irrep, and only a few sit in the two dimensional irrep. These are area law states, and could remain area law upon the addition of such a small perturbation. Thus some vestiges of localization are expected to remain at large $\lambda$ for a finite-size $L$. More numerical work at larger sizes will have to be done to distinguish these effects from a truly thermodynamic intermediate phase.

V Conclusion and outlook

In this paper, we construct and study a spin-1 Hamiltonian invariant under an on-site $S_{3}$ symmetry using various disorder and symmetry diagnostics. We study the eigenstate phases that can arise via ergodicity and symmetry breaking. Within the accuracy of our numerical analysis, we can identify three regions in the two-parameter Hamiltonian space two of which are consistent with thermal, MBL- spin-glass phases and a third whose identity is not established with certainty in the present study. We state our observations about various characteristics of this region and speculate with regard to its identity.

There are other interesting questions that are left for future study. In this work, we are limited to small system sizes by the tools of numerical analysis employed $i.e.$ exact diagonalization. There is much to be gained in designing and employing other numerical techniques to study larger system sizes. In this regard, it would be useful to explore the extension of tensor network techniques, which have been shown to be effective in the case of full-MBL, where eigenstates are expected to have area-law properties, to other settings, in particular to further study the disordered region with unbroken $S_{3}$ symmetry. Another extension of is to repeat our study in a Floquet setting where exotic possibilities have been conjectured in the presence of global symmetries floquet_classification_shivaji ; floquet_classification_vishwanath ; floquet_classification_dominic ; time_crystals_shivaji ; time_crystals_dominic . Here, the arguments for instability of disordered Floquet systems with non-Abelian symmetries to thermalization would be stronger because of energy pumping. However, as we saw in the equilibrium case, such a system might have interesting features worth exploring. Furthermore, we could consider relaxing the setting of indefinitely stable phases, like MBL, to the so-called ‘pre-thermal’ setting prethermal_abanin2015asymptotic ; prethermal_dominic where the system is stable for times that are exponential in system size, and study the role of non-Abelian symmetries.

*Note added: * As we were finishing this draft, we learned of a related paper: “Localization-protected order in spin chains with non-Abelian discrete symmetries” (arXiv:1706.00022) by Aaron J. Friedman, Romain Vasseur, Andrew C. Potter, S. A. Parameswaran.

Acknowledgements

A.P would like to acknowledge the stimulating environment of the Boulder Summer School for Condensed Matter Physics-2016 which is supported by the NSF and the University of Colorado, where this project was initiated. We are grateful to Rahul Nandkishore, Bela Bauer, Sid Parameshwaran and Andrew Potter for helpful discussions. T.-C.W is supported by NSF Grants No. PHY 1314748 and No. 1620252. LF is supported by NSF DMR-1519579 and Sloan FG-2015-65244.

Appendix A Incompatibility of non-Abelian symmetries with full MBL

Here, we review the main hypothesis of Potter and Vasseur [23]. The authors consider the case of a fully MBL system (as opposed to the case of a partially local, partially thermal system with mobility edges), and study the compatibility of MBL in the presence of various global symmetries. The working definition of an MBL system they consider is the existence of a complete set of quasi-local conserved quantities with associated quasi-local projectors in terms of which the Hamiltonian, $H$ can be defined as [44, 13, 45],

[TABLE]

Here, $L$ is the number of spins, $\hat{P}[i]_{\alpha}$ is the projector onto the $\alpha^{th}$ quasi-local conserved quantity at the $i^{th}$ location, $D$ is the number of conserved quantities and the $E$’s are constants that fix the energy eigenvalues. Furthermore, we can apply a finite depth quantum unitary circuit, $W$, that re-expresses the conserved quantities as local degrees of freedom. These local objects are called l-bits in terms of which the original spins (p-bits 222Note that we do not restrict ourselves to a local two-dimensional Hilbert space. This means, we should be talking about p-dits and l-dits instead of p-bits and l-bits. However, to keep with standard terminology, we will use the latter names.) and Hamiltonian can be defined. This is typically called the l-bit Hamiltonian, $H_{lbit}$.

Let us now consider the case when the Hamiltonian $H$ has a global on-site symmetry $i.e$ $H$ commutes with the unitary representation of some group, $G$ of the form $U(g)=\bigotimes_{i=1}^{L}V_{i}(g)$, where, $V_{i}(g)$ acts on each physical spin.

[TABLE]

In this case, representation theory of the group $G$ plays a role in constraining the allowed form of the l-bit Hamiltonian [11]. For spins (p-bits) to allow a well defined on-site group action, the local Hilbert space of each spin must correspond to some faithful representation of that group. This must also be true for a tensor product of the p-bits that constitute an l-bit. In other words, we can write the l-bit- projectors $\hat{P}[i]_{\alpha}$ using a fully reduced basis that can be labeled as $|\Gamma,m_{\Gamma};d_{\Gamma}\rangle$

[TABLE]

We define the different labels below and compare them with the well known case of the representations of the rotation group $SO(3)$:

• •

$\Gamma=1\ldots N_{R}$ labels the irreducible representation (irrep) of the group and is equivalent to the total angular momentum quantum number, $j$. The number of values it can take is equal to the number of irreps of $G$, $N_{R}$.

• •

$m_{\Gamma}=1\ldots|\Gamma|$ is equivalent to the azimuthal quantum number $m_{j}$. The number of values it can take is equal to the dimension of the irrep, $|\Gamma|$.

• •

$d_{\Gamma}=1\ldots D_{\Gamma}$ labels which of the $D_{\Gamma}$ copies of the $\Gamma$ irrep is being considered.

In this basis, the action of the group is

[TABLE]

Demanding the invariance of the Hamiltonian [11] under group action and invoking Schur’s lemma [47] irrep-wise, we get the constrained form of the Hamiltonian compatible with on-site symmetry as

[TABLE]

where,

[TABLE]

We now consider the cases of Abelian and non-Abelian symmetries separately. If $G$ is an Abelian group, all irreps are one dimensional ($|\Gamma|=1$). This means that all projectors $P_{\Gamma;d_{\Gamma}}$ are rank-1 which preserves the form [11]. With sufficient disorder, resulting in sufficiently random $E^{\prime}s$, we can imagine that all degeneracies are lifted and we can obtain an MBL phase stable to perturbations. However, if $G$ is non-Abelian, not all irreps are one dimensional. This means that we invariably have higher-rank local projectors giving us only a partial set of conserved quantities rather than complete which leads to degeneracies that are extensive in system size. This is clearly seen by examining the eigenvalues and eigenvectors of Eq [15].

[TABLE]

Note that none of the eigenvalues, $\mathcal{E}$ have any labels corresponding to the inner multiplicity of the irreps $m_{\Gamma}$. This means that the eigenstate $|\{\Gamma,m_{\Gamma};d_{\Gamma}\}\rangle$ has a degeneracy of $|\Gamma_{1}|\times|\Gamma_{2}|\times\ldots|\Gamma_{L}|$ which is clearly extensive in system size when $G$ is non-Abelian. The authors of Ref 23 state that under the influence of perturbations, such a degeneracy is susceptible to long range resonances which destabilizes MBL. Furthermore, they suggest a possible set of phases for the system to be in depending on the nature of the global symmetry group $G$. Here we list the possibilities for the case of finite groups such as $S_{n\geq 3}$ and $D_{n\geq 3}$, (i) Ergodic/ thermal phase, (ii) The so-called MBL spin-glass (SG) [48] phase characterized by localization with symmetry spontaneously broken (SSB) to an Abelian subgroup, and (iii) The so-called quantum critical glass phase (QCG) [49, 50, 51] characterized by critical scaling of entanglement entropy. Within the accuracy of our numerical analysis, our findings are consistent with the conjecture of Potter and Vasseur.

Appendix B Constructing symmetric Hamiltonians

In this appendix, we give details of how the Hamiltonian used in the main text, Eq [4] was constructed. We also detail a general technique to construct local symmetric operators with which we can build spin Hamiltonians invariant under any on-site symmetry in any dimension. The construction of the 1D $S_{3}$ invariant Hamiltonian of the main text is a specific application of this general technique.

B.1 Building the $S_{3}$ invariant Hamiltonian

B.1.1 Basis properties of $S_{3}$ and its representation used.

We first review some basic properties of the group $S_{3}$ and its representation used in this paper. $S_{3}$, the symmetry group of three objects is the smallest non-Abelian group. It is of order 6 and can be generated using two elements and the following presentation

[TABLE]

It has three irreducible representations, $\mathbf{1},\mathbf{1}^{\prime},\mathbf{2}$ which can be written as

$\chi^{\mathbf{1}}(a)=1$, $\chi^{\mathbf{1}}(x)=1$, 2. 2.

$\chi^{\mathbf{1^{\prime}}}(a)=1$, $\chi^{\mathbf{1^{\prime}}}(x)=-1$, 3. 3.

$\Gamma^{\mathbf{2}}(a)=\begin{pmatrix}\omega&0\\ 0&\omega^{*}\end{pmatrix}$, $\Gamma^{\mathbf{2}}(x)=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$ ,

where $\omega=e^{2\pi i/3}$. The local Hilbert space we have chosen for each spin lives is the three-dimensional reducible representation $\mathbf{2}\oplus\mathbf{1}^{\prime}$. We use the eigenspace of the spin-1 angular momentum operator $S^{z}$ to label the irreps. The $\mathbf{2}$ irrep is encoded in the two-dimensional $\pm 1$ eigenspace of $S^{z}$ (which we will call $|\pm\rangle$) and the $\mathbf{1}^{\prime}$ is encoded in the one-dimensional [math] eigenspace of $S^{z}$ (which we will call $|0\rangle$). The matrix representation of a general group element in this basis looks like

[TABLE]

In particular, the generators have the following representation

[TABLE]

B.1.2 1-spin operator

We first construct a 1-spin invariant operator. We start with a general operator that acts on the space of a single spin

[TABLE]

We now demand that $\hat{\lambda}$ is invariant under conjugation by $V(g)$, $i.e.$, the representation of symmetry on a single spin.

[TABLE]

Schur’s lemma [47] constrains the matrix elements of $\hat{\lambda}$ in the following way:

$\hat{\lambda}$ cannot mix basis states corresponding to different irreps. 2. 2.

$\hat{\lambda}$ must be proportional to the identity operator when acting on basis states corresponding to the internal states of the same irrep. 3. 3.

If there are multiple copies of the same irrep, $\hat{\lambda}$ can mix the basis states corresponding the same internal state of different copies but should still be proportional to the identity operator as an action on the internal states.

The meaning of these constraints should become clearer with the applications that will follow. For a single spin operator, since we have only one copy of each irrep, constraint 3 does not apply. Applying constraints 1 and 2, we get the form of $\hat{\lambda}$:

[TABLE]

From this, we can read off the only non-trivial 1-spin symmetric operator, $(S^{z})^{2}$ which is also Hermitian.

B.1.3 2-spin operator

In order to find a symmetric 2-spin operator, we follow the same logic as that of a 1-spin operator. First, we start with a general operator that acts on the 9 dimensional vector space of 2 spins,

[TABLE]

Next, we demand invariance under conjugation by $V(g)\otimes V(g)$, $i.e.$, the representation of symmetry on two spins.

[TABLE]

We now have to impose the constraints coming from Schur’s lemma. However, in order to do that, we need to find out the irrep content of $V(g)\otimes V(g)$. For this, we first list the Clebsch Gordan (CG) decomposition that gives us the outcomes of fusing different $S_{3}$ irreps. This is the generalization of angular momentum addition of $SU(2)$ irreps. Note that we exclude the trivial case of fusion with the trivial irrep $\mathbf{1}$,

[TABLE]

The irrep content of $V(g)\otimes V(g)$ is obtained from the CG decomposition,

[TABLE]

It is clear that we are bound to have multiple copies of the same irrep in this decomposition for which, using Constraint 3 imposed by Schur’s lemma, unlike the single spin case, we can get off diagonal operators. Let us list and label the different instances of each irrep appearing in the decomposition for convenience.

[TABLE]

The subscripts label the copy of the irrep. We next need the basis of the representation of each irrep in $V(g)\otimes V(g)$. These can be written in terms of the original basis states (labeled by $S^{z}$ eigenvalues) using CG coefficients which we calculate using the technique by Sakata [52] (also see [53]).

[TABLE]

Using this, we have the 2-spin $S_{3}$ symmetric operator constrained by Schur’s lemma

[TABLE]

As in the case of 1-spin invariant operator, we can again read off the several independent symmetric 2-spin operators by simplifying Eq .28. However, since we need the operators to be Hermitian, we take Hermitian combinations of these operators. We finally list the non-trivial independent Hermitian operators expressed in terms of spin-1 operators.

[TABLE]

To construct the Hamiltonian .4, we have used $(S^{z})^{2}$ and $\hat{J}_{1}$ to build the disordered part $H_{d}$ and $\hat{J}_{3}$, $\hat{J}_{4}$, $\hat{J}_{7}$ to build the thermal part of the Hamiltonian, $H_{t}$. Since $\hat{J}_{5}$ and $\hat{J}_{6}$ are mirrored versions of $\hat{J}_{4}$ and $\hat{J}_{7}$ respectively, we can keep our thermal term sufficiently generic even if leave them out. Note that while we have constructed a Hamiltonian for a 1D spin chain, the symmetric operators constructed using this formalism can be used to build Hamiltonians for any spatial dimensions.

B.2 General technique

We now give details of a general procedure that can be applied to obtain n-spin symmetric operators invariant under a representation of any on-site symmetry group. The schematic procedure is as follows

To construct symmetric n-spin operators, we first write down general operators that act on the Hilbert space of n spins. 2. 2.

We then demand the invariance of this operator under symmetry action. 3. 3.

Using Schur’s lemma, we constrain the matrix elements of the n-spin operator and read off the independent operators. 4. 4.

If required, we finally take the hermitian combinations of the independent operators.

If $G$ is the group we are considering, the most general local Hilbert space for the spin compatible with G-action can be an arbitrary number of copies of each irrep of $G$. Like mentioned in Sec A, we can choose the basis as $|\Gamma,m_{\Gamma};d_{\Gamma}\rangle$ where the symbols have the same meaning as before. The matrix representation of the group operators in this basis is block diagonal.

[TABLE]

The “passive” group action on the basis is

[TABLE]

B.2.1 1-spin symmetric operator

Let us now start with symmetric 1-spin operators. The most general 1-spin operator we can write down is

[TABLE]

Note that for notational convenience, here and henceforth, we assume summation over repeated indices. Demanding invariance under conjugation with $V(g)$, we have

[TABLE]

In matrix form, the condition on $\lambda$ becomes

[TABLE]

This means that $[\lambda_{\Gamma^{\prime};d_{\Gamma^{\prime}}}^{\Gamma;d_{\Gamma}}]$ is an intertwiner between the irreps $\Gamma$ and $\Gamma^{\prime}$. Such a matrix is constrained by Schur’s lemma:

[TABLE]

Using this in Eq 32, we get the form of a symmetric 1-spin operator,

[TABLE]

Here, the non-zero matrix elements $\lambda^{\Gamma}_{d_{\Gamma},d^{\prime}_{\Gamma}}$ are free parameters that act on the degenerate subspace associated with the outer multiplicity of each irrep.

B.2.2 2-spin symmetric operator

We now consider 2-spin symmetric operators. The most general operator that acts on the 2-spin Hilbert space is

[TABLE]

For this operator to be symmetric, it needs to be invariant under conjugation by $V_{1}(g)\otimes V_{2}(g)$,

[TABLE]

To use the techniques like we did for the 1-spin operator in the previous subsection, we need to first block-diagonalize $V_{1}(g)\otimes V_{2}(g)$ using a suitable basis change and then use Schur’s lemma. This redefinition of the basis states can be done using Clebsch-Gordan (CG) coefficients. Recall that the irrep content of a direct product of two irreps is schematically given by the CG series:

[TABLE]

$N_{\Delta\Lambda}^{\Gamma}$ denotes the number of copies of the $\Gamma$ irrep that exists in the fusion outcome of $\Delta\otimes\Lambda$. At the level of representations, Eq .39 tells us that there exists a change of basis by a unitary matrix $C_{\Delta\Lambda}$ that fully reduces the direct product of the irreps $\Delta\otimes\Lambda$:

[TABLE]

$C^{\Gamma;\alpha_{\Gamma}}_{\Delta\Lambda}$ are isometries whose matrix elements, $[C^{\Gamma;\alpha_{\Gamma}}_{\Delta\Lambda}]_{m_{\Delta},m_{\Lambda}}^{m_{\Gamma}}$ are the CG coefficients that project $\Delta\times\Lambda$ onto the $\alpha_{\Gamma}^{th}$ copy of the irrep $\Gamma$, where, $\alpha_{\Gamma}=1\ldots N_{\Delta\Lambda}^{\Gamma}$. It is useful to look at the ‘passive’ action of the CG coefficients on the basis kets to remind us of the actual “change-of-basis” action,

[TABLE]

Now, consider the following equivalences

[TABLE]

Here,

[TABLE]

is the number of ‘fusion channels’ of the kind $\Delta\otimes\Lambda\rightarrow\Gamma$ available to produce the irrep $\Gamma$ using the irreps in the 2-spin Hilbert space. In short, Eq [44] tells us that there exists some unitary matrix $W$ which block diagonalizes $V_{1}(g)\otimes V_{2}(g)$

[TABLE]

When viewed as a matrix, $W$ contains operations to both reorder indices appropriately as well as use the CG coefficients to reduce the direct product of irreps form $V_{1}$ and $V_{2}$ block-by-block. If we operate $W$ on both sides of the Eq [38] and call $WJW^{\dagger}=K$, we get

[TABLE]

The matrix elements of $K$ can be written in terms of those of $J$ and CG coefficients.

[TABLE]

Eq [48] is of the same form as Eq [33] and we can use Schur’s lemma again to constrain the elements of $K$,

[TABLE]

Note that $c_{\Gamma}$ is a collective index of compatible $\left(d_{\Delta},d_{\Lambda},\alpha_{\Gamma}\right)$ that runs over the $\mathcal{D}_{\Gamma}$ different fusion channels mentioned above. Also, note that we use the convention $\left[C^{\Gamma;\alpha_{\Gamma}}_{\Delta\Lambda}\right]^{m_{\Gamma}}_{m_{\Delta},m_{\Lambda}}=0$ if $\Gamma\notin\Delta\otimes\Lambda$. Finally, we can undo the transformation of $W$ and get the elements of $\hat{J}$. Since it is important, we expand $c_{\Gamma}$:

[TABLE]

Plugging in Eq [51] into Eq [37], we get the general form of the symmetric 2-spin operator.

[TABLE]

This can be greatly simplified using Eq [43]

[TABLE]

Where, we have once again reintroduced the short hand notation $c_{\Gamma}$ to denote the fusion channels labeled by compatible $(d_{\Delta},d_{\Lambda},\alpha_{\Gamma})$ to produce $\Gamma$. $K^{\Gamma;(d_{\Delta},d_{\Lambda},\alpha_{\Gamma})}_{(d_{\Delta^{\prime}},d_{\Lambda^{\prime}},\beta_{\Gamma})}$ are now the free parameters. In this form, we can clearly see the similarity with the symmetric 1-spin operator Eq [36]. This helps us see the general picture with arbitrary n-local operators. In a fully reduced basis of the n-spin Hilbert space, matrix elements of symmetric operators can only act on the outer multiplicity space of each irrep.

Appendix C Detecting the irrep of the eigenstates

In this appendix, we give details of how we determine the irrep a given eigenstate transforms as. If $U(g)=\bigotimes_{i=1}^{L}V_{i}(g)$ is the many-body representation of the on-site symmetry, we want to find out the irrep $\Gamma$ that an eigenstate $|\epsilon\rangle$ or more generally, a set of degenerate eigenstates $\{|\epsilon\rangle_{a}\}$ transform as

[TABLE]

This is an easy task for the symmetry group $SU(2)$, where, all we need to do is operate the eigenstates with the total angular momentum operator,

[TABLE]

This way, given an eigenstate that transforms as an irrep of $SU(2)$, we can extract the quantum number $j$ which labels the irrep. For general finite groups, we are not aware of an equivalent technique. However, we now present a strategy that works well for the group $S_{3}$ in the form of the following theorem.

Theorem 1**.**

Given a normalized vector $|\epsilon\rangle$, $\langle\epsilon|\epsilon\rangle=1,$ that transforms as some irreducible representation of $S_{3}$, $\Gamma=\mathbf{1},\mathbf{1^{\prime}}$ or $\mathbf{2}$, we can determine $\Gamma$ using the following two real numbers,

[TABLE]

where, $a$ and $x$ are the two generators of $S_{3}$ with presentation $\langle a,x|a^{3}=x^{2}=1,xax=a^{-1}\rangle$. Specifically,

[TABLE]

Proof.

Let us first consider the case when $|\epsilon\rangle$ is some vector in the 2D irrep $\mathbf{2}$. We can expand this vector in the orthonormal eigenbasis of the generator $a$ of the $\mathbf{2}$ representation, $|\omega\rangle,|\omega^{*}\rangle$, with eigenvalues $\omega,\omega^{*}$ respectively, where $\omega=e^{\frac{2\pi i}{3}}$ (see Appendix. B).

[TABLE]

Acting on it by $U(a)$,

[TABLE]

Let us now consider the case when $|\epsilon\rangle$ transforms as either of the 1D irreps. Since the representation of the generator $a$ is simply $1$ for both 1D irreps, we simply have

[TABLE]

Thus, we can see that $\mathcal{A}$ can separate the 2D irrep $\mathbf{2}$ form the 1D irreps. Also, $\mathcal{A}=-0.5$ is necessary and sufficient for $|\epsilon\rangle$ to transform as $\mathbf{2}$. Furthermore, if $\mathcal{A}=1$, we can determine which 1D irrep $|\epsilon\rangle$ transforms as by considering the transformation under the generator $x$ whose representation is $\pm 1$ for $\Gamma$ being $\mathbf{1}$ or $\mathbf{1^{\prime}}$ respectively. Thus, if $\mathcal{A}=1$,

[TABLE]

implies $|\epsilon\rangle$ transforms as $\mathbf{1}$ or $\mathbf{1^{\prime}}$ respectively. This concludes the proof. ∎

In our numerics, we diagonalize our Hamiltonian in the 1D irrep sector by projecting it into the appropriate basis. As discussed above, this means that we need to restrict to the basis states that are left invariant under the action of $U(a)=\bigotimes_{i=1}^{L}V(a)_{i}$. To see how this is done, consider the action of operator $U(a)$ on a many-body basis state labeled by $S^{z}$ eigenvalues on each spin

[TABLE]

We need $\omega^{S^{z}_{tot}}=1$ which means

[TABLE]

Since $3^{L-1}$ of the $3^{L}$ basis states satisfy the condition of Eq .61, this helps us diagonalize an $L$ site Hamiltonian for the price of $L-1$.

Appendix D Spin glass diagnostics for $S_{3}$ subgroups

In this appendix, we give details of how the spin-glass (SG) order parameter used in Sec. III.2 was constructed and also numerical evidence for the assertion that the high disorder region at $\kappa=0$ does not in any form spontaneously break the $S_{3}$ symmetry.

First, let us list the elements of $S_{3}=\{1,a,a^{2},x,xa,xa^{2}\}$ and its five subgroups written in terms of the generators $a,x$ defined in Sec. II and Appendix. B:

$\mathbb{Z}_{3}=\{1,a,a^{2}\}$, 2. 2.

$\mathbb{Z}_{2A}=\{1,x\}$, 3. 3.

$\mathbb{Z}_{2B}=\{1,xa\}$, 4. 4.

$\mathbb{Z}_{2C}=\{1,xa^{2}\}$, 5. 5.

$\{1\}$.

For each subgroup $H\subset G$, we design a SG diagnostic that detects SSB of $G\rightarrow H$ and takes the form

[TABLE]

$\mathcal{O}_{H}$ are local Hermitian order parameters that are chosen to have the following properties under symmetry transformation by $U(g)=\bigotimes_{i=1}^{L}V(g)_{i}$

$\mathcal{O}_{H}$ transforms trivially under $H$:

$U(g)\mathcal{O}_{H}U^{\dagger}(g)$ = $\mathcal{O}_{H}$, $\forall g\in H$. 2. 2.

$\mathcal{O}_{H}$ transforms non-trivially under $G/H$:

$U(g)\mathcal{O}_{H}U^{\dagger}(g)$ $\neq$ $\mathcal{O}_{H}$, $\forall g\notin H$.

Note that $\chi^{SG}_{H}$ is invariant under the redefinition $\mathcal{O}_{H}\rightarrow\mathcal{O}_{H}+\xi\mathbb{1}$. It can be checked that the following operators satisfy the above properties

• •

$\mathcal{O}_{\mathbb{Z}_{3}}$ = $S^{z}$

• •

$\mathcal{O}_{\mathbb{Z}_{2A}}$ = $S^{z}S^{x}+S^{x}S^{z}$ = $\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\ 1&0&-1\\ 0&-1&0\end{pmatrix}$

• •

$\mathcal{O}_{\mathbb{Z}_{2B}}$ = $V(a)\mathcal{O}_{\mathbb{Z}_{2A}}V^{\dagger}(a)$ = $\frac{1}{\sqrt{2}}\begin{pmatrix}0&\omega&0\\ \omega^{*}&0&-\omega\\ 0&-\omega^{*}&0\end{pmatrix}$

• •

$\mathcal{O}_{\mathbb{Z}_{2C}}$ = $V^{\dagger}(a)\mathcal{O}_{\mathbb{Z}_{2A}}V(a)$ = $\frac{1}{\sqrt{2}}\begin{pmatrix}0&\omega^{*}&0\\ \omega&0&-\omega^{*}\\ 0&-\omega&0\end{pmatrix}$

• •

$\mathcal{O}_{\{1\}}$ = $S^{x}$

Note that it is important to make sure that the disconnected part of the two point correlation function, $\langle\epsilon|\mathcal{O}_{H;i}|\epsilon\rangle\langle\epsilon|\mathcal{O}_{H;j}|\epsilon\rangle$ is subtracted when constructing $\chi^{SG}_{H}$. In previous work like Ref [8], the local order parameter $\sigma^{z}$ transformed as a non-trivial irrep of $G/H$ ($\mathbb{Z}_{2}$ in their case). It is then automatically guaranteed that $\langle\epsilon|\sigma^{z}_{i}|\epsilon\rangle=0$. Similarly, the local order parameter used in the main text to detect $S_{3}\rightarrow\mathbb{Z}_{3}$ SSB, $S^{z}_{i}$ transforms as a non-trivial irrep of $S_{3}/\mathbb{Z}_{3}\cong\mathbb{Z}_{2}$ which also ensures $\langle\epsilon|S^{z}_{i}|\epsilon\rangle=0$ and hence we leave out of the definition of $\chi^{SG}_{\mathbb{Z}_{3}}$. This would not be true if we used a different $\mathcal{O}_{\mathbb{Z}_{3}}$ like $S^{z}+\xi(S^{z})^{2}$ or even $S^{z}+\xi\mathbb{1}$ both of which are equally good to detect $S_{3}\rightarrow\mathbb{Z}_{3}$ SSB but would need subtraction of the disconnected part. Similarly, the other SG order parameters we used namely $\mathcal{O}_{\mathbb{Z}_{2A/B/C}}$ and $\mathcal{O}_{\{1\}}$ also need subtraction.

Fig. 10 shows the different SG diagnostics as a function of $\lambda$ for $\kappa=0,1$ averaged over eigenstates across disorder realizations. It can be seen that only $\overline{\chi^{SG}_{\mathbb{Z}_{3}}}$, which detects SSB $S_{3}\rightarrow\mathbb{Z}_{3}$ approaches a value that increases with system size in the region discussed in Sec. III.2. The value of other SG diagnostics becomes increasingly smaller or approaches a constant value with system size for all $\lambda$ and $\kappa$ indicating that SSB to that residual subgroup has not taken place in the eigenstates.