\section{QMC methods}
\subsection{Finite temperature auxiliary field QMC method}
As mentioned in the main text, we investigate an half-filling $\pi$-flux extended Hubbard model on square lattice, its Hamiltonian is $H=H_0+H_U$ with non-interacting part $H_0=-\sum_{\langle i j\rangle, \alpha} t_{i j}\left(c_{i \alpha}^{\dagger} c_{j \alpha}+\text { H.c. }\right)$ and interacting part $H_U=U \sum_{\square}\left(n_{\square}-2\right)^{2}$.
Since $H_0$ and $H_U$ do not commute, we use Trotter decomposition to separate $H_0$ and $H_U$ in the imaginary time propagation
\begin{equation}
	Z=\operatorname{Tr}\left[\mathrm{e}^{-\beta(H-\mu N)}\right]=\operatorname{Tr}\left[\left(\mathrm{e}^{-\Delta_{\tau} H_{U}} \mathrm{e}^{-\Delta_{\tau} H_{t}}\right)^{M}\right]+\mathcal{O}\left(\Delta_{\tau}^{2}\right),
\end{equation}
where $Z$ is the partition function, $\beta=M\Delta\tau$ is the inverse temperature, $\mu$ is chemical potential and equals $0$ at half-filling.
In QMC method, we can only deal with quadratic fermionic operator, while $H_U$ contains the quartic term, thus we should employ a Hubbard-Stratonovich (HS) decomposition as
\begin{equation}
e^{-\Delta\tau U(n_{\square}-2)^{2}}=\frac{1}{4}\sum_{\{s_{\square,\tau}\}}\gamma(s_{\square,\tau})e^{-2\alpha\eta(s_{\square,\tau})}e^{\alpha\eta(s_{\square,\tau})n_{\square}}
\label{eq:decompo}
\end{equation}
with $\alpha=\sqrt{-\Delta\tau U}$, $\gamma(\pm1)=1+\sqrt{6}/3$,
$\gamma(\pm2)=1-\sqrt{6}/3$, $\eta(\pm1)=\pm\sqrt{2(3-\sqrt{6})}$,
$\eta(\pm2)=\pm\sqrt{2(3+\sqrt{6})}$ and the sum symbol is taken over the auxiliary fields $s_{\square,\tau}$ on each $\tau$-th time slice square plaquette. Now, the interacting part is transformed into quadratic term but coupled with an auxiliary field.
Following simulations are based on the single-particle basis $\boldsymbol{c} = \{c_1, c_2 \cdots c_{N} \}$, $N=L\times L$ is the system size, so we can use the matrix notation $K$ and $V$ to represent $H_0$ and $H_U$ operators.
We define the imaginary time propagators
\begin{equation}
\begin{aligned}
&U_{s_{\square}}\left(\tau_{2}, \tau_{1}\right)=\prod_{m=m_{1}+1}^{m_{2}} e^{\boldsymbol{c}^{\dagger} V\left(s_{\square,m\Delta\tau}\right) \boldsymbol{c}} e^{-\Delta_{\tau} \boldsymbol{c}^{\dagger} K \boldsymbol{c}}, \\
&B_{s_{\square}}\left(\tau_{2}, \tau_{1}\right)=\prod_{m=m_{1}+1}^{m_{2}} e^{V\left(s_{\square,m\Delta\tau}\right)} e^{-\Delta_{\tau} K},
\end{aligned}
\end{equation}
where $m_1 \Delta\tau=\tau_1$ and $m_2 \Delta\tau=\tau_2$.
Then partition funciton $Z$ can be rewritten as
\begin{equation}
Z=\sum_{\{s_{\square,\tau} \}} \operatorname{Tr}\left[U_{s_{\square}}(\beta, 0)\right] \prod_{m=1}^{M} \gamma(s_{\square,m\Delta\tau})e^{-2\alpha\eta(s_{\square,m\Delta\tau})} = \sum_{\{s_{\square,\tau} \}} \operatorname{det}\left[1+B_{s_{\square}}(\beta, 0)\right] \prod_{m=1}^{M} \gamma(s_{\square,m\Delta\tau})e^{-2\alpha\eta(s_{\square,m\Delta\tau})}.
\end{equation}
Physical observables are measured according to
$\langle O \rangle = \frac{\operatorname{Tr}\left[\mathrm{e}^{-\beta H} O\right]}{\operatorname{Tr}\left[\mathrm{e}^{-\beta H}\right]}$.
The equal-time single-particle Green function is given by
$
G_{i,j}(\tau,\tau)=\left\langle c_{i,\tau} c_{j,\tau}^{\dagger}\right\rangle=\left(1+B_{s_{\square}}(\tau, 0) B_{s_{\square}}(\beta, \tau)\right)_{i, j}^{-1}
$, and the dynamical single-particle Green function is given by
$G_{i,j}(\tau_1,\tau_2)=\left\langle c_{i,\tau_1} c_{j,\tau_2}^{\dagger}\right\rangle=-\left[\left(\mathbf{1}-\mathbf{G}\left(\tau_{1}, \tau_{1}\right)\right) B_{s_{\square}}^{-1}\left(\tau_{2}, \tau_{1}\right)\right]_{i, j}$.
Other physical observables can be calculated from single-particle Green function through Wick's theorem. More technical details of the finite-temperature QMC algorithms can be found in the reference book~\cite{assaadWorld-line2008}.

\subsection{Projection QMC method}
Since we also want to investigate the ground state properties, the projection quantum Monte Carlo (PQMC) method is a good choice.
We can calculate the ground state wave function $\vert \Psi_0 \rangle$ through the projection of a trial wave function $\vert \Psi_T \rangle$ as $\vert \Psi_0 \rangle = \lim\limits_{\Theta \to \infty} e^{-\frac{\Theta}{2} \mathbf{H}} \vert \Psi_T \rangle$, where $\Theta$ is the projection length.
Actually, in PQMC method, we use $\langle\Psi_{T}|e^{-\Theta H}|\Psi_{T}\rangle$ to replace the role of partition function in finite-temperature version.
Physical observables $O$ can be measured according to
\begin{equation}
\label{eq:observablepqmc}
\langle O \rangle = \frac{\langle \Psi_0 \vert O \vert \Psi_0 \rangle}{\langle \Psi_0 \vert \Psi_0 \rangle}
						= \lim\limits_{\Theta \to \infty} \frac{\langle \Psi_T \vert  e^{-\frac{\Theta}{2} \mathbf{H}} O  e^{-\frac{\Theta}{2} \mathbf{H}} \vert \Psi_T \rangle}{\langle \Psi_T \vert  e^{-\Theta \mathbf{H}} \vert \Psi_T \rangle} .
\end{equation}

After the Trotter and HS decomposition, we have
\begin{equation}
\langle\Psi_{T}|e^{-\Theta H}|\Psi_{T}\rangle=\sum_{\{s_{\square,\tau}\}}\left[\left(\prod_{\tau}\prod_{\square}\gamma(s_{\square,\tau})e^{-2\alpha\eta(s_{\square,\tau})}\right)\det\left[P^{\dagger}B(\Theta,0)P\right]\right]
\label{eq:mcweight}
\end{equation}
here $P$ is the coefficient matrix of trial wave function $|\Psi_T\rangle$; $B$ matrix is defined as
$B(\tau+1,\tau)=e^{V_\tau}e^{-\Delta_\tau K}$
and has a property $B(\tau_3,\tau_1)=B(\tau_3,\tau_2)B(\tau_2,\tau_1)$.
In the practice,  we choose the ground state wavefunction of the half-filled non-interacting parts of Hamiltonian as the trial wave function.
The Monte Carlo sampling of auxiliary fields are further performed based on the weight defined in the sum of  Eq.~\eqref{eq:mcweight}. Single particle observables are measured by Green's function directly and many body correlation functions are measured from the products of single-particle Green's function based on their corresponding form after Wick-decomposition.
The equal time Green's function are calculated as
\begin{equation}
G(\tau,\tau)=1-R(\tau)\left(L(\tau)R(\tau)\right)^{-1}L(\tau)
\end{equation}
with $R(\tau)=B(\tau,0)P$, $L(\tau)=P^{\dagger}B(\Theta,\tau)$.
The imaginary-time displaced Green's function
$G(\tau,0) \equiv \left\langle \mathbf{c} \left(\frac{\tau}{2}\right) \mathbf{c}^{\dagger}\left(-\frac{\tau}{2}\right)\right\rangle$
are calculated as
\begin{equation}
\left\langle \mathbf{c} \left(\frac{\tau}{2}\right) \boldsymbol{c}^{\dagger}\left(-\frac{\tau}{2}\right)\right\rangle=\frac{\left\langle\Psi_{T}\left|e^{-\left(\frac{\Theta}{2}+\frac{\tau}{2}\right) H} \boldsymbol{c} e^{-\tau H} \boldsymbol{c}^{\dagger} e^{-\left(\frac{\Theta}{2}-\frac{\tau}{2}\right) H}\right| \Psi_{T}\right\rangle}{\left\langle\Psi_{T}\left|e^{-\Theta H}\right| \Psi_{T}\right\rangle}
\end{equation}
More technique details of PQMC method, please also refer to Refs~\cite{assaadWorld-line2008}.
\subsection{Infinite $U$ QMC algorithm}
For infinite $U$ calculation, we use following formula to perform the $SU(N)$ infinite $U$ projection~\cite{ouyangProjection2021}
\begin{equation}
 \left.e^{-\Delta_\tau U (n_{\square}- 2)^2}\right|_{U\rightarrow+\infty}
=  \frac{1}{M}\sum_{s_{\square,\tau}=1}^M e^{\frac{\text{i} 8\pi s_{\square,\tau}}{M} (n_{\square}- \frac{N}{2})},
\end{equation}
where we can set $M=2N+1$.
As the infinite $U$ term has now been replaced by fermion bilinears coupled to auxiliary fields, we can further use above finite temperature and projection QMC scheme to perform the calculation.

\begin{figure}[htb!]
\includegraphics[width=1.0\columnwidth]{figS1}
\caption{(a) $1/L$ extrapolation of single-particle gap. (b) $1/L$ extrapolation of VBS structure factor. The dash lines are obtained by cubic polynomial fitting with least squares method. (c) The AF structure factors as function of interaction strength. Error bars are smaller than the symbols.}
\label{fig:AF-order}
\end{figure}

\section{Other QMC data}
\subsection{The single-particle gap and VBS order parameter near QCP of SM-VBS }
The Dirac SM to VBS phase transition is a transition from massless Dirac fermion to insulator, thus single-particle gap $\Delta_{\text{sp}}$ will open at Dirac cones, i.e. at the momentum point $\mathbf{K}_0$.
We can extract $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L)$ from a fit to the asymptotic long imaginary time behavior of the single-particle Green's function
$G(\mathbf{k}, \tau, L)=\left(1/L^{4}\right) \sum_{i, j, \sigma} e^{i \mathbf{k} \cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}\left\langle c_{i, \sigma}^{\dagger}(\tau) c_{j, \sigma}(0)\right\rangle \propto e^{-\Delta_{\mathrm{sp}}(\mathbf{k},L) \tau}$, then extrapolate it to the thermodynamic limit (TDL).
As shown in Fig.~\ref{fig:AF-order}~(a), it is clear that $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) \rightarrow 0$ at $U/t=32$ in the SM phase, $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) > 0$ at $U/t=38$ in the VBS phase, and  $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) \sim 0$ at $U_c/t=36$.
The $1/L$ extrapolation of $C_{\text{VBS}}$ has the similar behavior as $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L)$ near $U_c$, as shown in Fig.~\ref{fig:AF-order}~(b).

\subsection{The absence of AF order}
In Hubbard-like model, the AF order usually dominate when VBS order become weak in strong coupling region~\cite{xuKekule2018,liaoGross2022}.
Here, we demonstrate that the AF order is absent in our model.
We define the structure factor of AF order as
$
C_\text{AF} (\mathbf{k}, L) \equiv \frac{1}{L^4}\sum_{\mathbf{r}} e^{-i\mathbf{q}\cdot\mathbf{r}}\sum_{\alpha,\beta}\langle S^{\alpha}_{\beta}(i) S^{\beta}_{\alpha}(j)\rangle
$, $C_\text{AF}$ will peak at momentum point $\mathbf{M}$ for AF order.
We can perform $1/L$ extrapolation of $C_\text{AF}(\mathbf{M},L)$ to get the AF structure factor at TDL, markded as $C_\text{AF}(\mathbf{M},\infty)$.
As shown in Fig.~\ref{fig:AF-order} (c), $C_\text{AF}(\mathbf{M},\infty)$ vanish at whole interaction strength range $U/t \in [0,160]$, which mean the AF order is absent in our model.

\section{Dynamical spin spectrum of free $\pi$-flux model}
We have a free $SU(4)$ $\pi$-flux model on square lattice with Hamiltonian
\begin{equation}
	H=- \sum_{\langle i j\rangle, \alpha} t_{i j}\left(c_{i \alpha}^{\dagger} c_{j \alpha}+\text { H.c. }\right).
\end{equation}
, and consider two-site unit cell, with inner cell coordinates $\mathbf{v}_{1}=(0,0)$,
$\mathbf{v}_{2}=(1,0)$.
We transform the Hamiltonian to the momentum space,
\begin{equation}
H=-t\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\left(\begin{array}{cc}
2\cos k_y & 1+e^{-2ik_{x}}\\
1+e^{2ik_{x}} & -2\cos k_y
\end{array}\right)c_{\mathbf{k}}
\end{equation}
where $c_{\mathbf{k}}=(c_{1,\mathbf{k}},c_{2,\mathbf{k}})^{T}$. We further write
$H$ in the diagonal basis
\begin{align}
H & =\sum_{\mathbf{k},a,b,m,n}c_{a,\mathbf{k}}^{\dagger}(U_{\mathbf{k}})_{am}(D_{\mathbf{k}})_{mn}(U_{\mathbf{k}}^{-1})_{nb}c_{b,\mathbf{k}}\nonumber \\
 & \equiv\sum_{\mathbf{k},a}\epsilon_{a,\mathbf{k}}f_{a,\mathbf{k}}^{\dagger}f_{a,\mathbf{k}}
\end{align}
where we use $a,b,m,n\in [1,2] $ to denote sublattice index, $(D_{\mathbf{k}})_{mn}=\delta_{mn}\epsilon_{m,\mathbf{k}}$, $\epsilon_{m,\mathbf{k}}=\pm 2 \sqrt{\cos^2 k_x + \cos^2 k_y}$ is the dispersion,
and we define $f_{a,\mathbf{k}}\equiv U_{\mathbf{k}}^{-1}c_{a,\mathbf{k}}$. Here the dispersion
is flavor degenerate, and we have omitted the flavor index, . The sum
is over $\mathbf{k}$ in the small (rectangle) BZ, i.e. $k_{x}\in(-\frac{\pi}{2},\frac{\pi}{2})$ and $k_{y}\in(-\pi,\pi$).

To calculate spin spectrum, we can use fluctuation-dissipation theorem,
which relates the spin spectrum to the spin susceptibility. The spin
susceptibility has following form

\begin{equation}
\chi_{ab}^{R}(\mathbf{r}-\mathbf{r}',t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r},a,t),S_{\mu}^{\nu}(\mathbf{r}',b,t')]\rangle .
\end{equation}
Fourier transform
to momentum space, we get
\begin{equation}
\chi_{ab}^{R}(\mathbf{q},t-t')=-i\theta(t-t')\frac{1}{V}\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{q},a,t),S_{\mu}^{\nu}(-\mathbf{q},b,t')]\rangle
\end{equation}
, where $S^{\alpha}_{\beta}(\mathbf{k},a,t)=c^{\dagger}_{\mathbf{k},\alpha}c_{\mathbf{k},\beta}-\frac{\delta_{\alpha,\beta}}{4}\sum_{\gamma}c^{\dagger}_{\mathbf{k},\gamma}c_{\mathbf{k},\gamma}$ are the $SU(4)$ spin full operators with $\alpha,\beta,\gamma \in [1,4]$.

For non-interacting $SU(4)$ case, the time dependence of spin operator is
given by

\begin{align}
S_{\nu}^{\mu}(\mathbf{q},a,t) & =\sum_{\mathbf{k}}e^{iHt}(c_{a,\mathbf{k},\mu}^{\dagger}c_{a,\mathbf{k}+\mathbf{q},\nu}-\frac{\delta_{\mu\nu}}{N}\sum_{\lambda}c_{a,\mathbf{k},\lambda}^{\dagger}c_{a,\mathbf{k}+\mathbf{q},\lambda})e^{-iHt}\nonumber \\
 & =\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(f_{m,\mathbf{k},\mu}^{\dagger}f_{n,\mathbf{k}+\mathbf{q},\nu}-\frac{\delta_{\mu\nu}}{N}\sum_{\rho}f_{m,\mathbf{k},\rho}^{\dagger}f_{n,\mathbf{k}+\mathbf{q},\rho})e^{i(\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}})t}
\end{align}
Therefore
\begin{align}
\chi_{ab}^{R}(\mathbf{q},t-t')= & -i\theta(t-t')\frac{N^{2}-1}{V}\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(U_{\mathbf{k}+\mathbf{q}}^{-1})_{nb}(U_{\mathbf{k}})_{bm} [n_{F}(\epsilon_{m,\mathbf{k}})-n_{F}(\epsilon_{n,\mathbf{k}+\mathbf{q}})]e^{i(\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}})(t-t')}
\end{align}
Perform Fourier transformation to frequency space, we get
\begin{align}
\chi_{ab}^{R}(\mathbf{q},\omega)
= \frac{N^2-1}{V}\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(U_{\mathbf{k}+\mathbf{q}}^{-1})_{nb}(U_{\mathbf{k}})_{bm}\frac{n_{F}(\epsilon_{m,\mathbf{k}})-n_{F}(\epsilon_{n,\mathbf{k}+\mathbf{q}})}{\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}}+\omega+i\eta}
\end{align}

To get the full spin spectrum, we should keep in mind that above derivation
is in small BZ. Now we consider the spin correlation on a square lattice
(one site unit cell), and relate it to above formula.

\begin{equation}
\chi^{R}(\mathbf{r}_{i}-\mathbf{r}_{j},t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r}_{i},t),S_{\mu}^{\nu}(\mathbf{r}_{j},t')]\rangle
\end{equation}
denote $\mathbf{r}_{i}=\mathbf{r}+\mathbf{v}_{a}$, $\mathbf{r}_{j}=\mathbf{r}'+\mathbf{v}_{b}$, we have

\begin{equation}
\chi^{R}(\mathbf{r}+\mathbf{v}_{a}-\mathbf{r}'-\mathbf{v}_{b},t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r},a,t),S_{\mu}^{\nu}(\mathbf{r}',b,t')]\rangle\equiv\chi_{ab}^{R}(\mathbf{r}-\mathbf{r}',t-t')
\end{equation}
Perform fourier transformation, we get following relation
\begin{equation}
\chi^{R}(\mathbf{q},\omega)=\sum_{ab}\chi_{ab}^{R}(\mathbf{q},\omega)e^{-i\mathbf{q}(\mathbf{v}_{a}-\mathbf{v}_{b})},
\end{equation}
and the full spin spectrum is
\begin{equation}
S(\mathbf{q},\omega)=-2\text{Im}\chi^{R}(\mathbf{q},\omega)
\end{equation}

\end{document}


\section*{Funding}
Russian Science Foundation (RSF), project №19-72-10069 (experimental part), Ministry of Science and Higher Education of the Russian Federation, program №FFNN-2022-0016 for the Valiev Institute of Physics and Technology, Russian Academy of Sciences and the Foundation for the Advancement of Theoretical Physics and Mathematics BASIS (project nos. 20-1- 1-34-1, 21-1-3-40-1) (theoretical part).

\section*{Acknowledgments}
This research was performed according to the Development program of the Interdisciplinary Scientific and Educational School of Lomonosov Moscow State University "Photonic and Quantum technologies. Digital medicine". N. N. Skryabin and I. V. Dyakonov acknowledge support from the Russian Roadmap for Quantum computing (Contract No. 868-1.3-15/15-2021 dated October 5, 2021 and Contract No. P2154 dated November 24, 2021)

\section*{Supplemental document}
See Supplement for supporting content.

\bibliographystyle{unsrt}

\section{Derivation of the first-order autocorrelation function of  quasi-thermal light}

Let us derive the correlation function of the quasi-thermal field obtained using a disk rotating with an angular velocity $\Omega$ (Fig.~1 in the main text). To describe the field scattered at an angle to the disk surface normal, we follow the paper~\cite{Crosignani1971}. While in that work a plane wave fallen on the disk was bounded by a rectangular window, for our case we consider Gaussian beam, which amplitude is

\begin{equation}
    E(x,y)=E_0 \exp \left(-\frac{(x-L)^2+y^2}{2 w^2}\right).
\end{equation}
Here,  $w$ is a waist of the incident laser beam, which is focusing on the disk surface, $E_0$ corresponds to the field in the center of the beam waist. Also for simplicity we consider the case when the center of the beam waist is at point $\{x=L, y=0, z=0\}$.

Taking into account the intermediate formulas in \cite{Crosignani1971} and equation (13), the dependence of the first-order correlation function on time delay can be written as
\begin{multline}
        g^{(1)}(\tau) \varpropto \exp[\text{i}\omega \tau] \int\limits_{-\infty}^{+\infty}dx''dy''\exp\left[-\text{i}k_y y''(\cos\left(\Omega\tau\right)-1) \right]\exp\left[\text{i}k_y x''\sin\left(\Omega\tau\right)\right]\times\\
    \times E(x'',y''-\Omega\tau L)E(x'', y''),
\end{multline}
where $\omega$ is the frequency of the laser radiation.
After integrating and using series expansion of the trigonometric functions, one has:

\begin{equation}\label{eq:2_3}
    g^{(1)}(\tau) \varpropto \exp\left[-\cfrac{L^2\Omega^2\tau^2 }{4 w^2}\right]\exp \left[\text{i}(\omega+k_y L \Omega)\tau \right] =\exp\left[-\cfrac{V^2\tau^2 }{4 w^2}\right]\exp \left[\text{i}(\omega+\Delta)\tau \right] .
\end{equation}
Here, $V$ is speed of the disk point at the center of the waist, the first factor is the standard Gaussian envelope corresponding to the correlation properties of the thermal state, with the correlation time $w/V$, which corresponds to the time of passage of the disk point through the beam waist, and the second factor is sinusoidal oscillation, where $\omega$ is the frequency of the radiation incident on the disk, and $\Delta = k_y V$ corresponds to the Doppler shift associated with a non-zero projection of the disk velocity on the direction of light propagation.

\section{Intensity correlation function of two interfering quasi-thermal fields}

The Gaussian (thermal) distribution is remarkable by its higher-order moments can be expressed in terms of first-order moments.

Let's put the intensities at inputs $1$ and $j$ as $I_1$ and $I_j$, and intensities at outputs $1$ and $i$ as $I_1^{\text{(out)}}$ and $I_i^{\text{(out)}}$. In these terms the output second-order cross-correlation function equals:

\begin{multline}\label{eq:124}
    \left\langle {I_1^{\text{\text{(out)}}}I_i^{\text{\text{(out)}}}(\tau )} \right\rangle  = \left\langle {E_1^{\text{(out)}*} E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}{E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle  = \\
    = {\left\langle {E_1^{\text{(out)}*} {E_1^{\text{\text{(out)}}}}} \right\rangle } {\left\langle {E_i^{\text{(out)}*} (\tau ){E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle } + \left\langle {E_1^{\text{(out)}*} {E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle \left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}} \right\rangle \\ = \left\langle I_1^{\text{\text{(out)}}} \right\rangle \left\langle I_i^{\text{\text{(out)}}} \right\rangle +\left\langle {E_1^{\text{(out)}*} {E_i^{\text{\text{(out)}}}}(\tau )} \right\rangle \left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{\text{(out)}}}}} \right\rangle
\end{multline}

We are interested in the last term. Let's describe it in more detail:
\begin{multline}
    \left\langle {E_1^{\text{(out)}*} {E_i^{\text{(out)}}}(\tau )} \right\rangle  = {\left\langle {E_i^{\text{(out)}*} (\tau ){E_1^{\text{(out)}}}} \right\rangle ^ * } = \\ = \left\langle {\left( {M_{11}E_1^*  + M_{1j}E_j^* } \right)\left( {M_{i1}{E_1}(\tau ) + M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}{E_j}(\tau )} \right)} \right\rangle  = \\
 = M_{11}M_{i1}\underbrace {\left\langle {E_1^* {E_1}(\tau )} \right\rangle }_{\left\langle {I_1} \right\rangle g_1^{(1)}(\tau )} + M_{11}M_{ij}{{\mathop{\rm e}\nolimits} ^{\text{i}\varphi }}\underbrace {\left\langle {E_1^* {E_j}(\tau )} \right\rangle }_0 +\\+ M_{1j}M_{i1}\underbrace {\left\langle {E_j^* {E_1}(\tau )} \right\rangle }_0 +M_{1j}M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}\underbrace {\left\langle {E_j^* {E_j}(\tau )} \right\rangle }_{\left\langle {I_j} \right\rangle g_j^{(1)}(\tau )},
\end{multline}

\begin{multline}\label{eq:2_4}
    { {\left\langle {E_1^{\text{(out)}*} {E_i^{\text{(out)}}}(\tau )} \right\rangle } ^2} = { {M_{11}M_{i1}\left\langle {I_1} \right\rangle g_1^{(1)}(\tau ) + M_{1j}M_{ij}{{\mathop{\rm e}\nolimits} ^{i\varphi }}\left\langle {I_j} \right\rangle g_j^{(1)}(\tau )}^2} = \\
 \left\langle {I_1} \right\rangle \left\langle {I_j} \right\rangle M_{11}  M_{1j}  M_{ij} M_i1e^{-\text{i} \varphi } g_1^{(1)}(\tau ) g_j^{(1)}(\tau )^*+
 \left\langle {I_1} \right\rangle \left\langle {I_j} \right\rangle M_{11}  M_{1j}  M_{ij} M_i1 e^{\text{i} \varphi } g_j^{(1)}(\tau ) g_1^{(1)}(\tau )^*+\\
 \left\langle {I_1} \right\rangle^2 M_{11}^2 M_{i1}^2 g_1^{(1)}(\tau ) g_1^{(1)*}(\tau )+\left\langle {I_j} \right\rangle^2 M_{1j} ^2 M_{ij}^2 g_j^{(1)}(\tau ) g_j^{(1)*}(\tau )
 .
\end{multline}

Now we replace the first-order correlation functions here with \eqref{eq:2_3}. Let them differ in Doppler shift for fields at different inputs. This can be achieved by collecting radiation scattered on the same RGGD at different angles. Then assuming

\begin{gather}
    g_1^{(1)}(\tau ) = \exp \left(  - \frac{V^2\tau ^2}{4w^2} \right)\exp \left[ {  \text{i}\left( {{\omega} + {\Delta _1}} \right)\tau } \right],  \\
    g_j^{(1)}(\tau ) = \exp \left(  - \frac{V^2\tau ^2}{4w^2} \right)\exp \left[ {  \text{i}\left( {{\omega} + {\Delta _j}} \right)\tau } \right],\\
\text{and } d = {\Delta _j} - {\Delta _1}
\end{gather}
for the normalized cross-correlation function we have:
\begin{multline}\label{eq:g2t}
    g^{(2)}(\tau) = 1 + \exp{\left(-\frac{\tau ^2 V^2}{2 w^2}\right)} \frac{1}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)} \times \\
    \left[\left\langle I_1 \right\rangle^2 M_{11}^2 M_{i1}^2 +\left\langle I_j \right\rangle^2 M_{1j}^2 M_{ij}^2 +2 \left\langle I_1 \right\rangle \left\langle I_j \right\rangle M_{11} M_{i1} M_{1j} M_{ij} \cos \left(d \tau +\varphi \right)\right] .
\end{multline}
In the special case of the ideal symmetric beam splitter ($M_{11} = M_{12} = M_{21} = M_{22} = 1/\sqrt{2}$ and the only MEP equals $\pi$) \eqref{eq:g2t} is reduced to Eq.~(4) of the main article:
\begin{equation*}\label{eq:g2_ideal_BS}
       g^{(2)}(\tau) = 1+ \exp\left(-\frac{\tau ^2}{2 \sigma^2}\right)\times
       \left[\frac{ \left< I_1 \right>^2+\left< I_2 \right>^2}{\left(\left< I_1 \right>+\left< I_2 \right>\right)^2}-\frac{2 \left< I_1 \right> \left< I_2 \right> }{\left[ \left< I_1 \right>+\left< I_2 \right>\right] ^2}\cos \left(d \tau \right)\right].
\end{equation*}

In general case, \eqref{eq:g2t} can be represented as Eq.~(5) from the main article:
\begin{equation*}\label{eq:g2_approx}
    g^{(2)}\left( \tau  \right) = C + \exp \left( { - {\frac{\tau ^2}{2\sigma^2}}} \right)\left[ {A + B\cos \left( {d \tau  + \varphi_{\text{thermal}} } \right)} \right],
\end{equation*}

where parameters:

\begin{equation}
\begin{aligned}
A &= \frac{\left\langle I_1 \right\rangle^2 M_{11}^2 M_{i1}^2 +\left\langle I_j \right\rangle^2 M_{1j}^2 M_{ij}^2}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)};\\
B &= \frac{2 \left\langle I_1 \right\rangle \left\langle I_j \right\rangle M_{11} M_{1j} M_{ij} M_{i1}}{\left(\left\langle I_1 \right\rangle M_{11}^2+\left\langle I_j \right\rangle M_{1j}^2\right) \left(\left\langle I_1 \right\rangle M_{i1}^2+\left\langle I_j \right\rangle M_{ij}^2\right)};\\
C &=1;\\
\sigma &= \frac{w}{V}.
\end{aligned}
\end{equation}

\section{Implementation of the beam splitter with adjustable matrix element phase}

The arrangement of the polarization elements (Fig.~2 in the main text) in the order PBS, HWP at -22.5$^\circ$, tilted Glan-Thompson prism, and quarter-wave plate QWP at 45$^\circ$ and second PBS allows us to change the MEP of the transfer matrix from $0$ to $\pi$. Swapping the outputs of the circuit, we get the same matrix element phase (MEP) with the opposite sign. Thus, we can change the MEP from $-\pi$ to $\pi$. The  Fig.~\ref{fig:phase_eta} shows the dependence of the MEP and minimum intensity transmission coefficient $T_{min}$ on the prism transmission coefficient of horizontally polarized beam $\eta$. It can be seen that the MEP varies from $0$ to $\pi$, without zeroing out the transmission anywhere.

When calculating this dependence, we used Jones calculus:
\begin{equation}
   \vec{H}=\left(
\begin{array}{c}
 1 \\
 0 \\
\end{array}
\right),
\Vec{V}=\left(
\begin{array}{c}
 0 \\
 1 \\
\end{array}
\right)
\end{equation}

Glan-Thompson prism matrix can be represented as

\begin{equation}
    \hat{P}=\left(
\begin{array}{cc}
 \sqrt{\eta} & 0\\
 0 & 1
\end{array}
\right),
\end{equation}
where $\eta$ is transmission coefficient for muted polarization, which can be changed by tilting the prism around the vertical axis. In our case, in the initial position of the prism, when $\eta=0$, the horizontal polarization is not transmitted.

Matrices of wave plates are represented as
\begin{equation}
    \hat{HWP} = \left(
\begin{array}{cc}
 \cos ^2(\alpha)-\sin ^2(\alpha) & 2 \sin (\alpha) \cos (\alpha) \\
 2 \sin (\alpha) \cos (\alpha) & \sin ^2(\alpha)-\cos ^2(\alpha) \\
\end{array}
\right),
\end{equation}

\begin{equation}
    \hat{QWP} = \left(
\begin{array}{cc}
 \cos ^2(\beta)+\text{i} \sin ^2(\beta) & (1-\text{i}) \sin (\beta) \cos (\beta) \\
 (1-\text{i}) \sin (\beta) \cos (\beta) & \sin ^2(\beta)+\text{i} \cos ^2(\beta) \\
\end{array}
\right),
\end{equation}
where $\alpha$ and $\beta$ are angles of fast axes.

The resulting matrix $\hat{U}$ will be

\begin{equation}
    \hat{U} = \hat{HWP} \times \hat{P} \times \hat{QWP}.
\end{equation}
Then we can find an equivalent matrix with zeroed first row and column and determine the MEP of the element $\{2,2\}$ and the minimal transmission coefficient $T_{min}$.

\begin{figure}[htbp]
\centering\includegraphics[width=0.7\linewidth]{Supplement_phase_realization.pdf}
\caption{MEP value $\varphi/\pi$ and minimum intensity transmission coefficient  $T_{min}$ vs Glan prism transmission coefficient $\eta$ for horizontal polarization, which can be tuned by rotation of the prism. }\label{fig:phase_eta}
\end{figure}

\bibliographystyle{unsrt}


\section{Introduction}

The interplay between disorder and Coulomb interactions has been the subject
of intense research, mostly devoted to degenerate Fermi systems \cite{1,2,3,4
. A common feature observed in disordered-interacting condensed-matter systems
is a local depression in their single-particle density-of-states $\rho
$($\varepsilon$). This feature, referred to as a zero-bias anomaly (ZBA), is
anchored to the chemical potential of the system and appears in both, the
diffusive \cite{5,6}, and insulating regimes \cite{7,8}. A zero-bias anomaly
may be observed when inserting a particle into the many-body system faster
than the time it takes existing particles to relax to its presence. This may
be accomplished in a tunneling or photoemission measurement of disordered
metals and doped semiconductors. When the system is quantum-coherent, the
process is closely related to the Anderson orthogonality catastrophe \cite{9}.

An indirect way to monitor $\rho$($\varepsilon$), becomes possible in systems
where the disorder exceeds the critical value for Anderson-localization; the
competition between disorder and the unscreened Coulomb interaction slows down
the medium relaxation which, in turn, makes the system $\rho$($\varepsilon$)
observable in field-effect experiments \cite{10,11,12,13,14,15,16,17,18}. In
this case a ZBA appears as a modulation of the conductance versus
gate-voltage, G(V$_{\text{g}}$) \cite{14}. This feature, called a memory-dip
(MD), has been observed in several heavily-doped semiconductors where the
disorder necessary to render them Anderson-insulators is large enough to
reduce electronic relaxation rates many decades below the transition-times
associated with their conductivity \cite{19}. Anderson insulators in this
group exhibit glassy dynamics and are referred to as `electron-glasses'. These
should be distinguished from lightly-disordered systems (sometimes called
Coulomb-glasses) that do not exhibit MD in field-effect experiments \cite{19}.
The relation between the single-particle density-of-states and the MD has been
elucidated in \cite{12,14,17,18}.

Theoretical models that incorporate disorder and interactions are usually
concerned with the low-energy part of the single-particle density-of-states
The energy location of the states were expelled from the ZBA region has
received less attention. In particular, the Anderson-localized regime, where
$\rho$($\varepsilon$) has traditionally been derived via a classical
Coulomb-gap approach, the missing states are usually depicted as being evenly
spread outside the depleted region \cite{20,21,22,23}. Given that disorder in
Anderson insulators is in general significantly stronger than the Coulomb
interaction \cite{24}, this picture seems plausible and is consistent with
results of tunneling \cite{25,26}, and field-effect experiments on strongly
disordered systems \cite{27}.

Recently however, a strikingly different behavior was detected in a particular
version of amorphous indium-oxide In$_{\text{x}}$O films; as shown in figure 1

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4785in,
width=3.3399in

{Fig_1.eps
\caption{Conductance G vs. gate-voltage V$_{\text{g}}$ for a 20nm thick
In$_{\text{x}}$O film with carrier-concentration \textit{N}$\approx
$10$^{\text{19}}$cm$^{\text{-3}}$. The sample is separated from a degenerate
Si:B gate by 2$\mu$m layer of SiO$_{\text{2}}$. The figure shows a memory-dip
centered at V$_{\text{g}}$=0 where the sample was equilibrated for an hour
after being cooled from room-temperature. The dashed line delineates the
energy dependence of $\partial$n/$\partial$$\mu$, the material thermodynamic
density-of-states..
\end{figure}
The MD exhibited by specimen of this version reveal structure that resembles
that of a superconductor $\rho$($\varepsilon$) with its characteristic
`coherence-peaks' at the edges of the ZBA. The version of In$_{\text{x}}$O
used in this work is \textit{not} superconducting (at least down to $\approx
$0.28K) even when the system is in the diffusive transport regime \cite{28}.
Moreover, the indium-rich version of In$_{\text{x}}$O, which \textit{is}
superconducting when its disorder is sufficiently small, does show a
memory-dip when strongly-localized but without these side-shoulders \cite{27}.
Different versions of In$_{\text{x}}$O are distinguished by their In/O
composition that determines their carrier-concentration \textit{N} \cite{28}.
In terms of microstructure they are quite similar. Electron-diffraction
patterns of the superconducting version of the compound and that of the low
carrier-concentration sample used in this study are shown in Fig.2. A scan of
the radial intensity distribution of these patterns shown in Fig.3 is required
to be able to detect any difference between the two versions. Detailed study
of the In$_{\text{x}}$O versions including x-ray diffraction and
interferometry, and Raman-spectrometry showed that in terms of structural
properties, different versions only differ by small quantitative aspects like
the position of the boson-peak \cite{29}. However, in terms of disorder
perceived by the charge carriers, the difference may be substantial when
comparing In$_{\text{x}}$O with similar resistivities but different
carrier-concentrations. The system with the lower carrier-concentration has a
lower disorder (given the same resistivity), and therefore it may also be more homogeneous

\begin{figure}[ptb
\centering
\includegraphics[
height=1.8585in,
width=3.3399in

{Fig_2.eps
\caption{Electron diffraction-patterns of two versions of In$_{\text{x}}$O
films 20nm thick. The left pattern (marked High D) has carrier-concentration
\textit{N}$\approx$8x10$^{\text{20}}$cm$^{\text{-3 }}$ the right pattern
(marked Low D) is characteristic of the batch used in the present study with
\textit{N}$\approx$10$^{\text{19}}$cm$^{\text{-3 }}$
\end{figure}

\begin{figure}[ptb
\centering
\includegraphics[
height=2.3644in,
width=3.3399in

{Fig_3.eps
\caption{Intensity profile scan of the two diffraction-patterns shown in
Fig.2. The first ring diameters (in units of pixels) are shown for comparison.
Note the similarity and the small difference in the nearest-neighbor distance
of the two versions.
\end{figure}

Theoretically, piling-up of states at the edges of the ZBA is consistent with
the sum-rule for $\rho$($\varepsilon$). A non-monotonic $\rho$($\varepsilon$)
above the gap was demonstrated in numerical studies of disordered and
interacting quantum-systems, explicitly in strongly-correlated compounds
\cite{30,31,32,33,34} and in quantum-dots \cite{35}. A "shouldered" $\rho
$($\varepsilon$), quite similar in shape to the MD exhibited in Fig.1, was
obtained in Monte-Carlo simulations based on a random-displacement version of
the Coulomb-glass \cite{36,37,38}. The shoulders in this scenario were related
to a charge-ordering process in the spirit of formation of a "disordered
Wigner crystal or, perhaps more appropriately, a "Wigner-glass"
\cite{39,40,41,42}. The visibility of shoulders in \cite{36} is exponentially
reduced with the disorder-strength W. Consequently, this feature has been
resolved in the simulations only for W \textit{much} smaller than the Coulomb
interaction amplitude E$_{\text{C}}$ \cite{38,43}. By contrast, the ratio
W/E$_{\text{C}}$ in Anderson-insulators with Gaussian disorder is typically
larger than unity \cite{24}. In the sample shown in Fig.1 for example, this
ratio is estimated to be $\approx$6, which, according to numerical simulations
\cite{36,38}, seems unfavorable for a charge-ordering mechanism to be
effective. To account for a "shouldered-MD" (SMD), it may perhaps be necessary
to take additional factors into account than included in the considerations
used for the classical Coulomb-glass.

We present in this work extensive data pertaining to this new phenomenon and
offer a heuristic interpretation for the main effects that characterize it.
Some puzzling issues that need further experimental and theoretical
elucidation are pointed out.

\section{Experimental}

\subsection{Sample preparation}

The samples used in this study were amorphous indium-oxide (In$_{\text{x}}$O)
films made by e-gun evaporation of 99.999\% pure In$_{\text{2}}
O$_{\text{3-x}}$ onto room-temperature Si-wafers in a partial pressure of
1.3x10$^{\text{-4}}$mBar of O$_{\text{2}}$ and a rate of 0.3$\pm$0.1\AA /s.
Under these conditions the carrier-concentration \textit{N} of the samples,
measured by the Hall-Effect at room-temperatures, was \textit{N}$\approx
$(1$\pm$0.1)x10$^{\text{19}}$cm$^{\text{-3}}$. Using free-electron formula,
this carrier-concentration is associated with $\partial$n/$\partial\mu\approx
$10$^{\text{32}}$erg$^{\text{-1}}$cm$^{\text{-3}}$. The Si wafers (boron-doped
with bulk resistivity $\rho\leq$2x10$^{\text{-3}}\Omega$cm) were employed as
the gate-electrode in the field-effect experiments. A thermally-grown
SiO$_{\text{2}}$ layer, 2$\mu$m thick, served as the spacer between the sample
and the conducting Si:B substrate. The screening length of the material
$\lambda$ $\approx$($\pi$e$^{\text{2}}\partial$n/$\partial\mu$)$^{\text{-1/2
}$ is $\approx$2nm and therefore the voltage that actually affects the sample
is $\approx$10$^{\text{-3}}$ of V$_{\text{g}}$, the voltage applied between
the sample and the gate. At room-temperature the sheet resistance of the
samples used here ranged between 22.5k$\Omega$ to 25k$\Omega.$

Films thickness was measured in-situ by a quartz-crystal monitor calibrated
against X-ray reflectometry. Samples geometry was defined by the use of a
stainless-steel mask during deposition into rectangular strips 0.9$\pm$0.1mm
long and 1$\pm$0.1mm wide. Four different batches were made in this study and
here we report detailed results of three different samples from one of these
batches with a thickness of 20$\pm$5nm. The dimensionless parameter
k$_{\text{F}}\ell$ for these samples is 0.27-0.29, just on the insulating side
of the critical value for this material (0.31$\pm$0.1 \cite{28}).

\subsection{Measurement techniques}

Conductivity of the samples was measured using a two-terminal ac technique
employing a 1211 ITHACO current preamplifier and a PAR 124A lock-in amplifier
using frequencies of 64-71Hz depending on the RC of the sample-gate structure.
R is the source-drain resistance and C is the capacitance between the sample
and the gate. C in our samples was typically $\cong$10$^{\text{{\small -10}}
$F and R for the samples studied in this work ranged between 1.5-2.5M$\Omega$
$at$ 4.11K. Except when otherwise noted, the ac voltage bias in conductivity
measurements was small enough to ensure near-ohmic conditions. Except where
otherwise noted, measurements were performed with the samples immersed in
liquid helium at T$\approx$4.11K held by a 100 liters storage-dewar. This
allowed up to two months measurements on a given sample while keeping it cold
and in the dark. These conditions are essential for measurements where
extended times of relaxation processes are required at a constant temperature.
Fuller details of the field-effect measurements setup, samples configuration
and characterization are described in \cite{29}. An up to date list of the
studied Anderson-insulators that exhibit a memory-dip along with a description
of their systematic dependence on the carrier-concentration of the material
are given in \cite{19}.

\section{Results and Discussion}

\subsection{The basic facts associated with the SMD}

The SMD state is sensitive to the conditions under which it is measured in the
field-effect experiment. In particular, sweeping V$_{\text{g}}$ over a range
$\vert
$\Delta$V$_{\text{g}}
$\vert$
that exceeds the typical width $\Gamma_{\text{S}}$ of the shoulder results in
a shoulderless MD when a V$_{\text{g}}$-sweep is taken again along the
original interval and polarity. The protocol used throughout this work to
achieve this state (labeled as "reference") involves sweeping V$_{\text{g}}$
from 0V to +70V and back to 0V then the sample is relaxed at this position for
20 minutes allowing the ZBA to reform before taking the next G(V$_{\text{g}}$)
scan. An example is illustrated in Fig.4

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4785in,
width=3.3399in

{Fig_4.eps
\caption{Field-effect scans comparing the reference, shoulderless state (see
text) with the G(V$_{\text{g}}$) of the sample 20 minutes after being
quench-cooled from room temperature to the bath-temperature. The two plots
were taken with the same scan-rate $\partial$V$_{\text{g}}$/$\partial$t of
1.2V/s.
\end{figure}

However, this shoulderless G(V$_{\text{g}}$) turns out to be metastable; the
SMD state of the sample reappears if one let the sample equilibrates under the
initial V$_{\text{g}}$ for long enough time before taking a new sweep. The
recovery of the SMD is illustrated in Fig.5 and Fig.6 for the right-hand and
left-hand side G(V$_{\text{g}}$) scans respectively

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4215in,
width=3.3399in

{Fig_5.eps
\caption{The recovery of the SMD-state with relaxation time t$_{\text{rel}}$
starting from a "reference" state at t=0.08 hours. Data are shown here for the
right-hand side of the field-effect. Before taking each plot, a new reference
state was prepared by sweeping the gate voltage to +70 V and back to 0V. The
plots are displaced along the ordinate for clarity and each was taken with the
same sweep-rate of dV$_{\text{g}}$/dt=1.2V/s.
\end{figure}

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4128in,
width=3.3399in

{Fig_6.eps
\caption{The recovery of the SMD-state as function of the relaxation time
t$_{\text{rel}}$ starting from a reference state at t=0.25 hours. Data are
shown for the left-hand side of the field-effect. Before taking each plot, a
new reference state was prepared by sweeping the gate voltage to +70 V and
back to 0V. The plots are displaced along the ordinate for clarity. All plots
were taken with the same sweep-rate of dV$_{\text{g}}$/dt=1.2V/s. The dashed
line and the red arrow delineate the way we define the magnitude of the
shoulder at a given t$_{\text{rel}}$ (see text and Fig.7).
\end{figure}
These time-dependent data imply that the SMD state has a lower energy than the
reference state.

The dynamics associated with the SMD `rejuvenation' process described by the
data in Figs. 5 and 6 differs from that exhibited by that of the
electron-glass relaxation-process monitored through the conductance G(t)
which, characteristically, is logarithmic \cite{27}. The two types of
relaxations are compared in Fig.7 for one of the samples

\begin{figure}[ptb
\centering
\includegraphics[
height=3.3062in,
width=3.3399in

{Fig_7.eps
\caption{(a): The recovery of the shoulder with time characterized by $\Delta
$G defined by the construction shown in Fig.6 as a red arrow. The dashed line
in the figure is a fit to the stretched-exponential function with $\beta$=0.78
and $\tau$=23,400s. (b): The relaxation of this sample conductance monitored
by recording G(t) at V$_{\text{g}}$=45V after changing the gate-voltage from
the equilibrium state at 0V to 45V. The dashed line delineates the logarithmic
relaxation law. Note the two decades difference in time scales between the two
plots.
\end{figure}

Note first that the SMD evolution persists and is followed for two more
decades in time than the conductance relaxation. The main part of the change
in G(t) is over while that of $\Delta$G is still building up. Secondly,
although the two processes are probably related (and begin and end essentially
at the same times), the functional time dependence of $\Delta$G(t) being
stretched exponential with exponent $\beta$=0.78, suggests a process dominated
by a relatively narrow rate-distribution \cite{44}. This point will be
clarified below.

The disparity in the dynamics of the two types of slow relaxations is also
manifested in the dependence of the SMD magnitude on the sweep-rate. This is
shown in Fig.8 for another sample in the batch studied

\begin{figure}[ptb
\centering
\includegraphics[
height=3.8493in,
width=3.3399in

{Fig_8.eps
\caption{(a): Dependence of the SMD plots on the gate-voltage sweep-rate
dV$_{\text{g}}$/dt. (b): The magnitude of the shoulder relative to that of the
MD, defined as shown in the inset, as function of the sweep-rate
dV$_{\text{g}}$/dt.
\end{figure}
Note that $\Delta$G, the depth of the ZBA in these plots is changing
substantially with the sweep-rate (Fig.8a) owing to the relative fast
relaxation dynamics of these samples. The logarithmic dependence of $\Delta$G
on the sweep-rate is a common feature in electron-glasses \cite{45}. This
dependence carries-over to the ratio $\delta$G/$\Delta$G ( where $\delta$G is
the shoulder height, see: Fig.8b) because of the much weaker change of
$\delta$G with the sweep-rate.

\subsection{The shoulder origin - a heuristic interpretation}

Based on the data presented above, it is conjectured that the shoulders
reflect an increase in the density of states at the edges of the MD in the
same vein that the depression of G(V$_{\text{g}}$) signifies a depletion of
charge at the chemical potential set by v$_{\text{g}}^{\text{eq.}}$, the
equilibrium gate-voltage. The MD has been related to the soft gap at $\rho
$($\varepsilon$) resulting from the long-range Coulomb interaction
\cite{12,13,14,17,18,19,45}.

On its own, the interaction favors a correlated spatial distribution of the
states expelled from the low energy of $\rho$($\varepsilon$). Therefore, a
lower energy configuration should be attained when the distribution of the
nearest-neighbor distances is peaked at a value of the order \textit{N
$^{\text{-1/3}}$ where \textit{N} is the system carrier-concentration. This in
turn would show up as a peak in G(V$_{\text{g}}$) at E$_{\text{C}}\approx
$e$^{\text{2}}\mathit{N}^{\text{1/3}}$/$\kappa$ where $\kappa$ is the
effective dielectric constant of the medium. The shoulder in our samples is
peaked at E$_{\text{C}}\approx$7~meV and it tapers-off towards $\approx
$40-50~meV (Fig.1). For the batch used here, \textit{N}$\approx
10$^{\text{19}}$cm$^{\text{-3}}$, the observed peak of the shoulder is
consistent with E$_{\text{C}}\approx$e$^{\text{2}}\mathit{N}^{\text{1/3}
$/$\kappa$ if $\kappa\approx$45. The bare dielectric constant of the material
is $\kappa\approx$10 but the polarization of localized states may enhance it
to be consistent with this scenario.

\ To further asses the plausibility of this picture, we need to consider the
role of disorder. The precondition for the MD to be observable in the
field-effect experiment is a strong enough disorder \cite{19}. Explicitly, the
disorder has to be strong enough to Anderson-localize the system \textit{and}
to render its relaxation slower \cite{19} than the sweep-rate of V$_{\text{g
}$. On the other hand, too strong a disorder would defeat the
electron-electron tendency to form an ordered structure and randomize the
states spatial-distribution thus suppressing the shoulder. Intuitively, out of
the systems that exhibit MD, the best candidate to show shoulders is where the
disorder to interaction ratio is small. This favors systems with low
carrier-concentration \textit{N} because the disorder necessary to
Anderson-localize a system increases with its Fermi-energy E$_{\text{F
}\propto$\textit{N}$^{\text{\textit{2/3}}}$ while the interaction scales with
\textit{N}$^{\text{\textit{1/3}}}$. Note that the SMD discussed here has the
lowest \textit{N} among all the electron-glasses systems\textit{ }yet reported
\cite{19,46}. This may be in line with our charge-ordering scenario in the
sense that the present system has more favorable parameters than previously
studied electron-glasses to show shoulders in field-effect measurements.

But how favorable are the actual system parameters? The magnitude of the
interaction identified above with the position of the shoulder was
E$_{\text{C}}\approx$7~meV. The relevant \cite{47} disorder W may be estimated
from the width of the shoulder (Fig.1) yielding W/E$_{\text{C}}\approx$6. The
disorder to interaction ratio typical for a degenerate semiconductor ranges
between $\approx$3 to $\approx$100 \cite{24}. As anticipated, the current
system is at the "favorable" limit in terms of parameters. Yet, this ratio is
far from the W/E$_{\text{C}}\ll$1 limit where shoulders are observed in
classical simulations \cite{36,37,38}. Hopefully, the results reported here
may assist in identifying the missing ingredients from these simulations.

The presence of disorder is also essential for understanding the
time-dependent processes depicted in Figs.5 and 6; Sweeping the gate voltage
by a
$\vert
$\Delta$V$_{\text{g}}
$\vert
$>
$\Gamma_{\text{S}}$ (where $\Gamma_{\text{S}}$ is the shoulder-width) results
in a shoulderless state (Fig.4). However, given time, the shoulder is
sluggishly recovered and G(V$_{\text{g}}$) eventually exhibits the SMD
behavior. Presumably, the metastable MD-state is made up of charges trapped by
disorder over the voltage-range swept by the gate. Some of the charges that
were driven by the gate to occupy higher energy-states are trapped by deep
wells of the potential while the gate travels back to v$_{\text{g
}^{\text{eq.}}$. Charges occupying states with shallow wells of the potential
will escape these states during the time the gate is set back in v$_{\text{g
}^{\text{eq.}}$ to rebuild the MD for the next G(V$_{\text{g}}$) run. This
explains both, the shifting of charge away from the shoulder towards higher
energies, and the relative absence of the fast relaxation rates in the
relaxation associated with the shoulder recovery.

The proposed picture for creating the shoulderless state is analogous to the
process by which debris is washed ashore by sea-waves and is held on the beach
slope by friction when the wave pulls back seaward.

Despite the apparent difference in their dynamics, it is clear that the
shoulder and the memory-dip are two parts of the same phenomenon; when the
equilibrium voltage is moved, so is the position of the shoulder. An example,
illustrating the `two-dip-experiment' \cite{46} for the current batch of
In$_{\text{x}}$O, is shown in Fig.9

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4137in,
width=3.3399in

{Fig_9.eps
\caption{G(V$_{\text{g}}$) plots taken consecutively; First, V$_{\text{g}}$
was swept from the original equilibrium v$_{\text{g}}^{\text{eq.}}$=0V to +65V
(triangles). This is followed by sweeping V$_{\text{g}}$ to +20V and parked
there for 20 minutes. Next, V$_{\text{g}}$ was swept to +80V (squares).
V$_{\text{g}}$ was then moved to +20V and parked there for 20 minutes to
relax. Finally, V$_{\text{g}}$ was swept to -65V (circles). Note that the
latter G(V$_{\text{g}}$) plot was taken less than an hour after the first
sweep and the magnitude of the MD at v$_{\text{g}}$=0V has shrank by a factor
of $\approx$60 relative to the original value. All plots were taken with the
same dV$_{\text{g}}$/dt=0.5V/s.
\end{figure}

The feature that strikes the eye in Fig.9 is the barely discernible MD
produced in the third sweep at the original v$_{\text{g}}^{\text{eq.}}$=0V.
The magnitude of this MD is smaller by almost two orders of magnitude relative
to either of the two "fresh" MD's. This is due to the relatively fast dynamics
of this weak-disorder batch \cite{19,46,48}. A closer look at the figure
reveals a new aspect of the phenomenon; establishing a new v$_{\text{g
}^{\text{eq.}}$(at +20V) yields at the edge of the MD a prominent shoulder
despite being in voltage-range covered by V$_{\text{g}}$ in a previous sweep.
In a way, this is similar to that the shoulders appear along with the MD upon
the first cooldown from room-temperatures. Both the MD and the shoulder appear
in this case after a brief relaxation period (typically 20 minutes), no
days-long waiting-time is required. In fact, as will be shown next, once a new
equilibrium conditions are set, the full SMD appears in the G(V$_{\text{g}}$) plot.

\subsection{Shoulder destruction and rejuvenation}

As shown above, taking a G(V$_{\text{g}}$) measurement over a V$_{\text{g}
$-range larger than $\Gamma_{\text{S}}$, eliminates the shoulder from
showing-up in the next scan. A prolonged relaxation of the sample is required
for its recovery. However, there are several ways to restore the SMD without a
long time-delay. Two of these were mentioned in the paragraph above. Three
other schemes effective in rejuvenating the shoulder are described next.

The SMD's resulting from applying the first two are illustrated in Fig.10

\begin{figure}[ptb
\centering
\includegraphics[
height=2.4422in,
width=3.3399in

{Fig_10.eps
\caption{Field-effect plots illustrating two examples of the shoulder
reappearance following a quench from a nonequilibrium, higher-energy state:
Stressing the sample with a non-ohmic field F=10$^{\text{4}}$V/m held on for
20 minutes (circles), and exposing the sample to 1mA of the IR source for 3s
(triangles). In both cases the shown G(V$_{\text{g}}$) plots are taken after
letting the sample relax for 20 minutes since the respective quench. These are
contrasted with the reference plot produced after a scan to +70V. All three
plots were measured in the linear-response and with the same dV$_{\text{g}
$/dt=1V/s.
\end{figure}

Both schemes produce shoulders that exhibit higher visibility than would be
expected by allowing the system to relax at v$_{\text{g}}^{\text{eq.}}$ for
the time required for completing the scheme (compare with Fig. 5).

The third scheme is the easiest to implement and it turns out to be the most
effective in pulling back the shoulder from the reference state. In this case,
V$_{\text{g}}$ is first wept to V$_{\text{g}}$=
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert
$_{\text{1}}$ staying there for $\approx$2 seconds and then set back at
V$_{\text{g}}$=0V for 20 minutes relaxation before recording a G(V$_{\text{g
}$) plot from 0V to +60V. As noted before, this trip erases the shoulder. This
"v$_{\text{g}}^{\text{bias}}$" protocol is repeated for different (negative)
bias-values thus generating plots as function of
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert
. In each such plot the shoulder is "rejuvenated" by pre-biasing the sample
with a finite
$\vert
v$_{\text{g}}^{\text{bias}}$. These bias values were taken at random order to
verify that `history' does not play a role. A series of such plots with
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert$
values ranging from 5V to 80V is shown in Fig.11a

\begin{figure}[ptb
\centering
\includegraphics[
height=4.1796in,
width=3.3399in

{Fig_11.eps
\caption{The reappearance of the shoulder in the G(V$_{\text{g}}$) plots
recorded by the "v$_{\text{g}}^{\text{bias}}$ protocol" described in the text.
These are labeled by their pre-scan (negative) bias (a). The plots were taken
consecutively using the same sweep rate dV$_{\text{g}}$/dt=1V/s. (b) The
magnitude of the shoulder vs. the pre-bias voltage. This magnitude is taken as
the distance from the peak of the shoulder at 6.8V to the value of the
reference-curve at this voltage{\protect\small .}
\end{figure}

Note (Fig.11b) that the relative magnitude of the shoulder $\delta$G increases
with
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert$
and tends to saturate for v$_{\text{g}}^{\text{bias}}\simeq$-50V. The range
over which the gate-voltage bias affects $\delta$G is similar to the typical
width of the shoulder $\Gamma_{\text{S}}$ that, in turn, is related to the
disorder associated with the phenomenon W. It is emphasized that these data
are generated using random values for
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert$
in the consecutive measurements. Therefore the correlation between
$\vert
v$_{\text{g}}^{\text{bias}}
$\vert$
and $\delta$G is meaningful and suggestive. Actually, these results lead us to
recognize the common element in the four different schemes that produces a SMD
from a "reference" state: \textit{They all involve a quench from a state where
the charge-carriers are endowed with excess energy }E$_{\text{excess}
$\textit{ relative to the equilibrium state}. These schemes differ by how
E$_{\text{excess}}$ is structured (distributed in energy) relative to W, which
presumably determines its efficiency in pulling-up the shoulder. When the
excess-energy of charge-carriers exceeds W, they are no longer trapped by
potential wells of the disorder. Rather, their spatial distribution is
controlled by the Coulomb repulsion forcing them apart. In this stage the
system is in a fluid, spatially-correlated state that, following a quench and
brief relaxation, is well primed to form SMD.

The efficiency of a scheme to produce SMD is not necessarily reflected in the
conductance enhancement accompanying the process. For example, stressing the
sample with a large longitudinal field enhances the sample conductance by a
factor of $\approx$2 while the shoulder that results from this scheme is
visibly weaker (Fig.10) than that obtained by the brief IR illumination that
had a much smaller impact on the conductance (see Fig.12a,b)

\begin{figure}[ptb
\centering
\includegraphics[
height=4.2298in,
width=3.3399in

{Fig_12.eps
\caption{The effects produced on the sample by a brief IR illumination and a
large applied electric-field. (a) The conductance increase and subsequent
relaxation in response to a 3s exposure to the 1mA IR source. This is compared
to the similar effect produced by a sudden change of the gate-voltage from 0V
to 45V (a 5 seconds trip-time). (b) The sample resistance as function of the
applied electric field.
\end{figure}

\subsubsection{Temperature dependence}

Temperature has a marked effect on the electron-glass features. Both the width
of the memory-dip and its relative magnitude are systematically affected by
raising the temperature \cite{49}. Here it is shown that the shoulder follows
the same trend; its peak is slightly shifted towards a higher V$_{\text{g}}$
upon increasing the bath temperature, and its magnitude is exponentially
suppressed (Fig.13b). This temperature dependence was obtained before in the
electron-glass phase of In$_{\text{2}}$O$_{\text{3-x}}$ films \cite{49}

\begin{figure}[ptb
\centering
\includegraphics[
height=4.2021in,
width=3.3399in

{Fig_13.eps
\caption{(a): The dependence of the G(V$_{\text{g}}$) plot on temperature for
a typical SMD sample. Each plot was taken after first sweeping the gate to
-70V and back to 0V to allow shoulder rejuvenation (see text) and relaxation.
Then the bath temperature was changed and the same protocol repeated with
V$_{\text{g}}$/dt=0.5V/s. (b): The relative magnitude of the ZBA and the
shoulder extracted from the data in (a). $\Delta$G/G for the ZBA is the value
of G(V$_{\text{g}}$) at V$_{\text{g}}$=7V. $\Delta$G/G values for the shoulder
are taken as in Fig. 6 relative to the dashed lines in (a). Dotted-lines fit
exp[-T/1.4] (see text). Note that both axes in plot (b) use logarithmic
scale.
\end{figure}

Apart from reaffirming that the shoulder is an integral part of the
memory-dip, it illustrates how sensitive is its visibility to the addition of
high-energy components.

\subsubsection{Nonequilibrium steady-state under IR illumination}

Even more striking is the effect of a continuous IR illumination which
suppresses the shoulder despite the extremely weak power-load on the sample (Fig.14)

\begin{figure}[ptb
\centering
\includegraphics[
height=3.5976in,
width=3.3399in

{Fig_14.eps
\caption{ Sample conductance versus time while being exposed it to a IR
radiation with a constant power of 2\textperiodcentered10$^{\text{-7}}$Watt.
Note that G(t) tends to reach a saturated value for t$\gtrsim$1 hour (b): The
sample G(V$_{\text{g}}$) plot in the dark (squares) compared with the
G(V$_{\text{g}}$) plot taken under the constant IR illumination after the
conductance reached a "saturated" conditions (explicitly, the change with time
for the duration of the G(V$_{\text{g}}$) plot is negligible). Both curves
were taken after first sweeping V$_{\text{g}}$ to -70V, then relaxing the
sample for 20 minutes, and using the same sweep-rate dV$_{\text{g}
$/dt=0.5V/s.
\end{figure}
Under a continuous IR illumination the sample SMD magnitude was reduced to
half its dark-value while the conductance increased by only $\simeq$2.8\%
(see, Fig.14a). Based on the conductance versus temperature of this sample
shown in Fig.15 this is the equivalent of raising the temperature by $\approx
$20mK. To achieve a similar reduction of the SMD magnitude of the same sample,
the bath temperature must be raised to T$\simeq$5.5K (see Fig.13) which
increases the conductance by $\simeq$300\%.

The difference between these two protocols is in how the energy added to the
system is distributed among the available degrees of freedom. Under the IR
illumination, the system is populated with an excess of high-energy phonons.
These are generated by a cascade process: Electrons are excited to
high-energy, then relax by optical-phonons emission \cite{50}. This
nonequilibrium steady-state has much fewer low-energy phonons than when the
sample is heated-up to $\approx$5.5K. The density of high-energy phonons is
also increased by raising the bath temperature, but most of the absorbed
energy is spent on enhancing the conductance. The indication is that the
SMD-state is sensitive to the presence of high-energy phonons. This may hint
on an underlying mechanism that counteracts the disorder in the competition
with interactions. High-energy phonons is a potent source of dephasing which
raises the possibility that quantum effects play a role. Effects that hinge on
quantum-coherence and may be relevant to the problem at hand are those that
promote spatially-extended wavefunctions like resonances \cite{51,52,53} and
other wavefunction forms resulting from hybridization. The existence of
extended wavefunctions naturally mimics a weaker disorder. There is already
evidence for quantum coherence in this system \cite{50, 55} so there is reason
to further explore this direction

\begin{figure}[ptb
\centering
\includegraphics[
height=2.9758in,
width=3.3399in

{Fig_15.eps
\caption{The temperature dependence (on a double-log scale) of the conductance
for the sample used in figures 8-14 inclusive
\end{figure}

\subsection{Discussion}

There are several issues related to the appearance of the shoulders that need
to be further clarified. First is the absence of this feature in other
versions of In$_{\text{x}}$O. In contrast with the samples reported here, all
studied Anderson-insulating In$_{\text{x}}$O samples with \textit{N }ranging
from\textit{ }$\approx$4x10$^{\text{19}}$cm$^{\text{-3}}$ to $\approx
$8x10$^{\text{21}}$cm$^{\text{-3}}$ exhibit electron-glass characteristics
with MD magnitude and width varying systematically with \textit{N } but
without shoulders \cite{19,27,48}. It seems therefore that, for In$_{\text{x
}$O, as \textit{N }increases from $\approx$1x10$^{\text{19}}$cm$^{\text{-3}}$
to $\approx$4x10$^{\text{19}}$cm$^{\text{-3}}$, the structure manifested by
the shoulders melts-away. In terms of In/O ratio, which is the way \textit{N}
is determined in the compound, this range is merely $\approx$2\% wide
\cite{28}. It would be interesting to fine-tune the system \textit{N} and
track the evolution of the phenomenon across this region. This is a much
harder undertaking than it may seem and will require establishing a new level
of material deposition and post-treatment control. Nonetheless, on the basis
of what is already established, this SMD$\rightarrow$MD transition is
surprisingly sharp.

It appears that one has to decide between two possibilities: Either the
transition is a real phase-change, or it is a crossover to a regime with ultra
slow dynamics. Explicitly, as \textit{N} increases, the shoulders-dynamics
becomes too slow to be resolved on a realistic time scale. This dilemma is
similar to the quandary that has haunted the glass community \cite{54}, and
perhaps it may take as long to reach a conclusion. In the present case the
possibility of a diverging relaxation time has empirical support; The
relaxation time of In$_{\text{x}}$O electron-glasses with \textit{N}$\approx
$4x10$^{\text{19}}$cm$^{\text{-3}}$ to \textit{N}$\approx$10$^{\text{20}
$cm$^{\text{-3}}$ was shown to increase dramatically with \textit{N}
\cite{45,46}, and as seen here. the dynamics associated with the shoulders may
actually be slower.

Yet, the possibility that the SMD to MD changeover is a true phase-transition
cannot be ruled-out either. If the MD$\rightarrow$SMD transition is indeed a
change between two distinct phases what are their nature? The SMD state,
characterized by a shoulder in the G(V$_{\text{g}}$) plot has been interpreted
here as a preponderance of charge at energy of the order E$_{\text{C}}$
relative to the chemical-potential. It is illuminating to contrast this
picture with the spatial-order of the amorphous structure manifested in the
electron-diffraction patterns (Fig.2, and Fig.3 above). The shoulders and the
diffraction-pattern represent two amorphous structures; one composed of ions,
and the other of electronic-charges. The first ring in the former is related
to the nearest-neighbor distance between ions which is the analog of the
shoulder position in energy. Both the diffraction-ring and the shoulder have a
width that is a measure of the disorder associated with the respective
phenomenon. The diffraction-ring is relatively much sharper than the
shoulder-width. This is due to the difference in the two-particle potential
that constrain spatial-freedom of the constituents; Lennard-Jones versus
Coulomb-interaction for the ions and charge-carriers respectively. For
energies outside the Coulomb gap, this description depicts the SMD-state as a
"frozen-liquid" while the MD-state, lacking `medium-range-order', is a
"frozen-gas". The experiments that were performed up to and including this
work indicate that there is a transition between these states controlled by
the carrier-concentration \textit{N}. This was ascribed to the competition
between disorder and interaction that is tilted in favor of the former when
\textit{N} increases. To further test this picture one needs to change
disorder, or the interaction, independently which is a challenging task.

A related issue is the use of lower carrier-concentrations than \textit{N
$\approx$10$^{\text{19}}$cm$^{\text{-3}}.$ As mentioned above, in the
Anderson-localized regime (a precondition to show electron-glass features
\cite{19}), both disorder W and the Coulomb-interaction E$_{\text{C}}$
increase with \textit{N} but W$\propto$\textit{N}$^{\text{2/3}}$ while
E$_{\text{C}}\propto$\textit{N}$^{\text{1/3}}$ and therefore a lower
disorder-interaction ratio is achievable by reducing \textit{N}. This however
has a limited scope; the ratio is a weak function of \textit{N} while the
dynamics, mainly controlled by W \cite{19,46}, depends on it exponentially.
When W decreases the relaxation will quickly become too fast for the
experimental $\partial$V$_{\text{g}}$/$\partial$t to expose a memory-dip in
the field-effect experiments. To our knowledge, there is no Anderson-insulator
with \textit{N}$\lesssim$10$^{\text{18}}$cm$^{\text{-3 }}$that exhibits either
MD or SMD.

There are other factors that may be detrimental for charge-ordering besides
too strong disorder. Obviously, sample inhomogeneity beyond that inherent in
the sample disorder (estimated `globally', usually on the basis of measuring
the resistance) will smear-out the shoulders even when all other conditions
are favorable. In particular, an excess of randomly distributed deep
potential-wells, should be avoided. Generally, the detailed structure of the
background potential, not just its average amplitude, is relevant for the
charge-ordering scenario.

A subtle question that is yet to be answered is the effect of a superimposed
periodic potential. The effect of underlying crystallinity on the
disorder-interaction balance is a pertinent question that may be tested by
numerical simulations. However, numerical studies usually use a lattice-model
with half-filling, which probably introduces a spatial constraint that likely
is incommensurate with the optimal charge arrangement. It is not clear that it
may suffice to employ just a very small filling-factor in the simulation to
capture the physics of charge-ordering.

Attempts to answer some of these questions may benefit from further
experimental work. Systems that exhibit electron-glass features, such as GeTe,
Tl$_{\text{x}}$O, In$_{\text{2}}$O$_{\text{3-x}}$, GeSb$_{\text{2}
$Te$_{\text{5}}$, and Be \cite{19} do not show shoulders in their
G(V$_{\text{g}}$) plots. It would be useful to know whether this is due to
unfavorable W/E$_{\text{C}}$ or because of their crystallinity. Filling factor
in these compounds is naturally small but their relatively high
carrier-concentration may be unfavorable. It is worth trying to tweak this by
varying their composition. These systems should also be tested for
quantum-coherence as this may play a role as alluded to above. A methodical
study of these aspects may establish a new platform for the interplay between
disorder and interaction in Anderson insulators.

\begin{acknowledgments}
Discussions with Oded Agam, and Ady Vaknin are gratefully acknowledged.
\end{acknowledgments}