\subsection{Phase diagram}
\vspace{0.3cm}
\noindent \textbf{Phase diagram.} As mentioned earlier, depending on the magnitude and angle of the stress applied, the system can be BG, RG, or slides continuously. When the system is subjected to a shear stress $\pmb{\tau}=\mathbf{F}/A$, where $\mathbf{F}$ is the applied force and $A$ the area of the flake, it can be studied using the bilayer enthalpy $H(\mathbf{r}) = V(\mathbf{r}) - \pmb{\tau} \cdot \mathbf{r} A$. As we will now see, the number of minima of $H$ determines the phase the system is in, giving a stress-angle phase diagram (Fig.~\ref{fig:potential_2d}g) for multilayer graphene.

In the pairwise model, if a stress $\pmb{\tau}$ is applied to the upper layer in a quasistatic transformation, then the layer below exerts all the remaining stress $-\pmb{\tau}$. The same applies to subsequent layers. Thus, for each pair of layers, their enthalpy $H$ is the same. Depending on the angle and magnitude of $\pmb{\tau}$, there are three possible situations, which correspond to the 3 colored regions of Fig.~\ref{fig:potential_2d}g:\\
(i) Blue region: 2 minima. If no stress is applied, $H=V$, and there are 2 minima $m_1$ and $m_2$ (see Fig.~\ref{fig:potential_2d}h). If shear stress is sufficiently low, $H$ still has 2 minima, and the system has $2^{N-1}$ minima. For each pair of layers, they do not escape the local minimum they are currently in. Since this holds for all layers, the full system does not escape the local minimum it is currently in. In particular, if BG is chosen as the starting structure, the system remains BG in the blue region (in the pairwise model, all stacking sequences have the same energy, but in reality BG is the most stable structure).\\
(ii) Green region: 1 minimum. As stress increases, $m_1$ disappears and $m_2$ remains close to its $\pmb{\tau}=0$ position. Barrier $b_1$ disappears, while $b_2$ remains. This occurs because the stress is more aligned with the direction in which the saddle point $b_1$ has a maximum (that is, the direction in which $b_1$ acts as a barrier) than with the corresponding direction of $b_2$ (except at $\theta=\ang{30}$, where both barriers are affected in the same way, and the system transitions directly from 2 minima to 0 minima). Since there is only one minimum, all layers are in the same position relative to the lower layer, and the resulting stacking is rhombohedral. This is a key observation of our work. From the convention in Fig.~\ref{fig:mechanical_model}a, $m_2$ corresponds to configuration AC, and layers are ACB-stacked. If stress is applied in the opposite direction, $m_1$ is the only remaining minimum, and layers are ABC-stacked.\\
(iii) Orange region: 0 minima. For larger stresses, there are 0 minima. Since there are no local minima, layers keep on sliding without reaching a stable configuration.
If the stress is eventually removed, the system will not necessarily be RG.
However, transformations at different angles as in Fig.~\ref{fig:force_6l} suggest that the system will still have a high degree of rhombohedral order.

In particular, for small angles, the structure remains fully RG. Also, $\theta=\ang{0}$ is the angle with the largest range of stress that results in RG, of about 0.3 GPa (from 0.2 to 0.5 GPa). Thus, the armchair direction is the most robust direction to obtain RG.

Fig.~\ref{fig:phase_diagram}g also shows with plus signs '+' the critical stress values. For each angle, the small critical stress corresponds to the lower value, and the big critical stress to the upper value. They coincide with the blue-green border, and green-orange border, respectively. This excellent agreement shows that the number of minima of $H$ does indeed define the stress-angle phase diagram.

It is worth pointing out that when $m_1$ becomes shallow enough, it might be possible for thermal fluctuations to excite layers from $m_1$ from $m_2$. This will depend on experimental conditions, like duration of the experiment, temperature, and size of the flakes. Thus, the curve that separates multiple minima from 1 minimum is actually an upper bound. Also, BG is more stable than RG and is presumably located in a deeper (local) minima (in an analogous fashion to the two minima of the blue curve of Supplementary Fig.~\ref{fig:6layer}). So for $\theta$ close to 30$^\circ$, it might occur that the system transitions directly from BG to continuous sliding.

\vspace{0.3cm}
\noindent \textbf{Pressure.} Other hydrostatic pressures were also considered.
Supplementary Fig.~\ref{fig:force_pressure} shows the phase boundaries at P=0 (blue-green and green-orange borders of Fig.~\ref{fig:phase_diagram}g, or lower and upper borders) and P=2 GPa. The values of the lower and upper boundary at $\theta=0$
increase approximately linearly with pressure, at about 0.07 and 0.18 GPa of shear stress per GPa of hydrostatic pressure, respectively. Thus, the amount of stress needed to obtain RG increases, but also the range of allowed values to obtain RG (which might increase the robustness of an experiment).


\vspace{0.3cm}
\noindent \textbf{Shear in previous works and superlubricity.} The values of stress to produce RG suggested by our calculations  are very similar to values already published in experimental reports.
The configurations we have described,
where layers are commensurate with each other, are referred to as \textit{lock-in} states in graphene literature related to friction or shear. In this type of systems, values of shear stress of the order of 0.1 GPa were measured \cite{Dienwiebel2004}, and observed to be in good agreement with previous calculations\cite{Verhoeven2004}.
In another experiment\cite{Liu2012b}, a microtip of a micromanipulator was used to apply a shear force on a graphene flake to ``unlock'' it (remove it from the minima), and based on the deformation of the tip, a value of 0.14 GPa was reported, also lower than 0.20 GPa .
On the other hand, when the layers are incommensurate with each other, the values of friction are 2 or 3 orders of magnitude lower. This phenomenon is referred to as superlubricity and has sparked a lot of interest.
Optimal conditions for superlubricity include big flakes, low temperatures and low loads\cite{Wijn2010}. Since layers have to be moved out of the local minimum in the mechanism we have proposed, depending on the experimental conditions, care might need to be taken to avoid the upper layer to rotate into a superlubricant state.

To conclude, we have described a mechanism to transform multilayer graphene into RG through shear stress, which implies that applying sufficient shear strain to graphite results in long-range RG. Also, existing experimental values of shear stress are similar to the ones suggested by our results.
Our model suggests a compelling method
for experimental groups trying to obtain multiple layers of rhombohedral-stacked graphene.





\clearpage
\small
\noindent \textbf{\large Methods}

\noindent \textbf{DFT calculations.}
Calculations were performed in Quantum Espresso (QE) \cite{QE} using an LDA functional. The bilayer potential and the six layer first-principles transformations were obtained with a cutoff of 80 Ry, a $k$-grid of $56\times 56\times 1$ and an electronic temperature (Fermi-Dirac smearing) of 284 K. The energies considered in this work are always in meV per interface atom. This means that the total energy of the system is divided by 2, whether the number of layers is 2 or 6 (which permits a direct comparison between energies of both systems). In each system, the energy is set to 0 at the relative position where the energy is the lowest. To converge the energy differences of Supplementary Table 1 with a precision below 0.01 meV/atom, we used an energy cutoff of 120 Ry and an electronic grid of $80 \times 80 \times 1$. In the rvv10 calculations of Supplementary Fig.~\ref{fig:rvv10}, 100 Ry were used.

The most important parameter in the calculations is basically the small barrier $V_\textrm{SP}$ of Fig.~\ref{fig:mechanical_model}a between the two minima, since it determines the amount of stress necessary to go from one local minimum to another. The value of $V_\textrm{SP}$ can be expected to be somewhat accurate if the curvature around the minima is accurate. The curvature is directly proportional to the frequency squared of the shear mode LO' at $\Gamma$ \cite{Mounet2005} (this is the mode in which layers move in the plane in opposite directions). Experimental values of the shear frequency vary between 42 cm$^{-1}$ and 45 cm$^{-1}$(ref. \cite{Lebedeva2011}). This is in good agreement with the value that results from fitting a parabola close to $m_1$ in Fig.~\ref{fig:potential_1d}a, 42 cm$^{-1}$. The curvature of other functionals like rvv10 (ref. \cite{Sabatini2013}), which includes van der Waals, is actually lower in our calculations, so they agree less well with experiments (see Supplementary Fig.~\ref{fig:rvv10}). In previous works, rvv10 also gives a lower frequency than LDA\cite{Grande2019}.

To be confident that it is sufficient to consider the 2 layer potential in our analysis, we considered a similar calculation, but with 6 layers (dotted-blue line in Supplementary Fig.~\ref{fig:potential_1d}a) instead of 2. In each calculation, the $x_\textrm{CM}$ and $y_\textrm{CM}$ of each of the three upper layers was moved with respect to 
the three lower layers, and the system was relaxed. This corresponds to a generalized stacking fault energy \cite{Wang2015}. The curves are very similar, and the energy difference between the minima is about 6\% of $V_\textrm{SP}$, so the transformations could be analyzed in terms of the bilayer potential. Indeed the curves in Fig.~\ref{fig:force_6l} are almost identical. Also, the value obtained for the stacking fault energy, 1.58 meV/atom, compares well with 1.53 meV/atom obtained with RPA \cite{Zhou2015}. A detailed comparison is included in Supplementary Table 1.

\vspace{0.3cm}
\noindent \textbf{Mechanical model.} Let us label L$i$ the $i$th layer from the bottom to the top, and let us consider that a force is applied to the right on the upper layer L6, as in the top part of Fig.~\ref{fig:mechanical_model}c.
First, L6 moves from B to C without resistance, where it locks with L5. Then, it pushes L5 from A to B, which in turn pushes L4 from B to C. The lower layers have not moved yet. In the next step, L2 and L3 move as well. The transformation is (labeling always from bottom to top): ABABAB$\rightarrow$ABABAC$\rightarrow$ABACBA$\rightarrow$ACBACB. If the direction is reversed, the steps of the transformation are: ABABAB$\rightarrow$ABABCA$\rightarrow$ABCABC. In both cases, RG is obtained.

\vspace{0.3cm}
\noindent \textbf{Transformation calculations.}
In Fig.~\ref{fig:force_6l}, in addition to moving the upper layer in steps of $d/12$ along the direction $\theta$ (and fixing it, while relaxing in the perpendicular direction), where $d$ is the bond length, after each step layer $i=2,..,5$ is actually moved as well in $[(i-1)/5]d/12$ (other layers are fully relaxed), giving a configuration closer to the equilibrium position. This gives a more accurate value of the position at which the abrupt changes in the stress take place, and which particular layers slide relative to each other in $y$ or in the perpendicular direction $x$. 

At $\theta=\ang{0}$ there is a bifurcation of behavior when the big jump occurs (with the upper layer jumping in $-a/2$ for $\theta <0$, and $a/2$ for $\theta>0$). So the result of Fig.~\ref{fig:force_6l}a actually corresponds to $\theta= \ang{0.1}$, to avoid an artificial behavior at $\theta=\ang{0}$.


\vspace{0.3cm}
\noindent \textbf{Fourier interpolation.}
\label{sec:Fourier}
To obtain the bilayer potential at an arbitrary $\mathbf{r}$ point as in Fig.~\ref{fig:potential_2d}b (see coordinate system there), and to then obtain the phase diagram, the potential was first calculated at a set of $N$ discrete points $\mathbf{r}_i=(x_{\textrm{CM},i},y_{\textrm{CM},i})$ of the upper layer within a primitive cell (separated by $a$/15 along the lattice vectors, with $a=2.46$ \AA\hspace{1mm} the length of the lattice vectors, totaling $N=225$). As mentioned earlier, the interlayer distance and relative coordinates were relaxed. Let $\mathbf{G}$ be a set of $N$ reciprocal lattice vectors around the origin with the symmetry of the crystal. Then, we can write

\begin{equation}
\begin{split}
V(\mathbf{r}) = \frac{1}{N} \sum_\mathbf{G} V_\mathbf{G} e^{i \mathbf{G}.\mathbf{r}} \\
\textrm{with } V_\mathbf{G} = \sum_\mathbf{r_{i}} V(\mathbf{r_i}) e^{-i \mathbf{G}.\mathbf{r_i}} \\
\end{split}
\end{equation}

\noindent where $\mathbf{r}$ can take any value.


\vspace{0.3cm}
\noindent \textbf{Hydrostatic pressure.}
In order to obtain a phase diagram as in Fig.~\ref{fig:phase_diagram}g when there is an external pressure $P$, we have to add an additional term to the enthalpy $H$ defined earlier, resulting in the new enthalpy

\begin{equation}
H(\mathbf{r}) = V(\mathbf{r}) - \pmb{\tau} \cdot \mathbf{r} A + P z A,
\end{equation}

\noindent where $z$ is the interlayer distance of the center of mass at pressure $P$.

\vspace{0.3cm}
\noindent \textbf{Acknowledgements}

\noindent The authors would like to thank Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the high-performance SeaWulf computing system, which was made possible by a \$1.4M National Science Foundation grant (\#1531492).
We also acknowledge support from the Graphene Flagship
Core 2 grant number 785219 and Agence Nationale de
la Recherche under references ANR-17-CE24-0030. This
work was performed using HPC resources from GENCI,
TGCC, CINES, IDRIS (Grant 2019-A0050901202) and PRACE.

\vspace{0.3cm}
\noindent \textbf{Author Contributions}

\noindent The project was conceived by all authors. JP.N. performed the numerical calculations, prepared the figures and wrote the paper with inputs from M.C. and F.M. All authors contributed to the editing of the manuscript.

\vspace{0.3cm}
\noindent \textbf{Competing Interests}

\noindent The authors declare no competing interests.


\setcounter{figure}{0} 

\begin{figure*}[!ht]
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{coord_6l_xy_both.pdf}}
\caption{
\textbf{a}, Change of $y_\textrm{CM}$	and $x_\textrm{CM}$ of each layer, as a function of the center of mass position $y_\textrm{CM}$ of the upper layer (L6), for the 6 layer pairwise model calculation (see Fig.~\ref{fig:force_6l}a). Starting from ABABAB, the upper layer is ``pushed'' along the bond direction: 1/12 of the bond length $d$ is added to $y_{\textrm{CM}}$ of the upper layer in each step, and then the structure is relaxed. The relaxed coordinates are used in the subsequent calculation. In the first sliding steps, layers only move in $y$. After RG is formed,
stress increases until the upper layers jumps in $x$ in $a/2$. \textbf{b}, Analogous plot, with $\theta=\ang{15}$. In this case, the structure does not remain fully rhombohedral after the jump in $x$, but is soon recovered after a sliding step in $y$.}
\label{fig:coord_6l}
\end{figCaption}
\end{figure*}

\begin{figure*}[h!]
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{potential_1d_6l.pdf}}
\caption{Full-black: Energy of the bilayer system along the armchair direction. Same as in Fig.~\ref{fig:mechanical_model}a.
Dotted-blue: Similar calculation with 6 as
opposed to 2 layers. The center of mass of the lower layers is fixed in the configuration ABA in all calculations, while the upper
3 layers move (the center of mass of each layer is fixed in each calculation, while other coordinates are relaxed). Labels $m_1$
and $m_2$ correspond to ABABAB and ABACBC, respectively. The energy difference corresponds to a stacking fault energy. M
corresponds to ABAACA. The similarity between the curves indicates that the 2 layer potential is a good approximation to
analyze how shear affects the stacking order. See Supplementary Table I for a more detailed comparison.}
\label{fig:6layer}
\end{figCaption}
\end{figure*}

\begin{figure*}[t]
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{LDA_fix_rvv10.pdf}}
\caption{Energy of the bilayer system along the armchair direction. Full-black: LDA relaxed calculation (same as Fig.~\ref{fig:potential_1d}a). Dashed-red: Fixed z (at the relaxed value of the AB configuration). Dotted-blue: Relaxed rvv10. Although in the AB-AC region the agreement is good, the disagreement between both curves increases as the configuration approaches AA. This is consistent with Supplementary Fig.~\ref{fig:dist_pressure}, which shows the variation of the interlayer distance $z_\textrm{CM}$.
$z_\textrm{CM}$ at the small barrier is 3.33 \AA, while at AA  is 3.60 \AA. Thus, it is important to relax the configuration. Otherwise, the potential differences are artificially high.}
\label{fig:rvv10}
\end{figCaption}
\end{figure*}

\begin{figure*}[t]
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{zpr.pdf}}
\caption{Interlayer distance of the bilayer system along the armchair direction, for various hydrostatic pressures. Full-black: $P=0$. Dotted-blue: $P=$ 1 GPa. Dashed-red: $P=$ 2 GPa. }
\label{fig:dist_pressure}
\end{figCaption}
\end{figure*}

\begin{figure*
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{energies_pr.pdf}}
\caption{Plots analogous to those of Fig.~\ref{fig:mechanical_model}a at varying pressures. Full-black: $P=0$ (same as Fig.~\ref{fig:potential_1d}a). Dotted-blue: $P=$ 1 GPa. Dashed-red: $P=$ 2 GPa. Supplementary Fig. 3 shows that fixed $z$ close to AA corresponds to artificial pressures of over 1 GPa.}
\label{fig:energy_pressure}
\end{figCaption}
\end{figure*}

\begin{figure*
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{phase_diagram_pressure.pdf}}
\caption{Phase boundaries at different hydrostatic pressures. For each color, the lower line corresponds to the border between the blue and green region of the phase diagram (Fig.~\ref{fig:phase_diagram}g), and the upper curve to the border between the green and orange regions. Blue: $P=0$ (same as in Fig.~\ref{fig:phase_diagram}f). Red: $P$= 2 GPa. In this range of pressures, the values of the curves at $\theta=\ang{0}$ increase linearly with pressure.
}
\label{fig:force_pressure}
\end{figCaption}
\end{figure*}


\begin{figure*}[!ht]
\begin{figCaption}
\centering
	{\includegraphics[width=0.5\textwidth]{jumps_vs_Nl.pdf}}
\caption{Shear stress on the upper layer as a function of $y_\textrm{CM}$ of the upper layer. The jump of stress, as sliding steps take place, decreases with the number of layers. The jump is also smaller after RG is formed.
This can be useful to apply a less varying and more controllable stress. If for example a cantilever is used to move the upper layer, the number of layers of the sample would have to be taken into account to decide on its optimal elastic constant. It can also be seen how the position of the first jump increases linearly with the number of layers.}
\end{figCaption}
\label{fig:jump}
\end{figure*}


\newpage

\begin{table*}[h]
\begin{tabCaption}
\caption{Comparison between the 2 layer potential and the 6 layer calculation of Fig.~\ref{fig:potential_1d}, and previous works. Ref.~\cite{Wang2015,Zhou2015} use the adiabatic-connection fluctuation-dissipation theorem within the random phase approximation (ACFDT-RPA). Ref.~\cite{Telling2003} uses anisotropic elasticity theory, revised values for the elastic constants and the experimental data of ref.~\cite{Baker1961,Amelinckx1965}, to obtain an average stacking fault energy of
0.14 meV per interface atom, the same value as the RPA calculation. The difference between the minima (stacking fault energy) is small compared to the barrier $V_\textrm{SP}$ separating them, which implies the 2 layer potential provides a good approximation to study the transformation through shear.
{\begin{tabular}{cccc}
\hline \hline
\multicolumn{1}{l}{}  & \multicolumn{3}{c}{Energy (meV/atom)} \\
\multicolumn{1}{l}{}  & \begin{tabular}[c]{@{}c@{}}2 layer\\ potential \end{tabular}  & \begin{tabular}[c]{@{}c@{}}6 layer\\ potential\end{tabular}  & \begin{tabular}[c]{@{}c@{}} Previous \\works \end{tabular} \\
\hline
\begin{tabular}[c]{@{}c@{}} Stacking fault\\ Bernal \end{tabular} & 0                 & 0.10   &  \begin{tabular}[c]{@{}c@{}} 0.14 (RPA) \cite{Wang2015}           \\ 0.14 (Exp.)\cite{Telling2003} \end{tabular} \\
\hline \\
Small barrier $V_\textrm{SP}$    & 1.58              & 1.59   & 1.53 (RPA) \cite{Zhou2015}           \\
\hline \\
Large barrier $V_\textrm{AA}$    & 9.7	               & 9.5  & \begin{tabular}[c]{@{}c@{}}8.8 (RPA) \cite{Zhou2015} \\ 12.4 (QMC) \cite{Mostaani2015} \end{tabular}             \\
\hline \hline \\
\end{tabular}} \\
\end{tabCaption}
\label{tab:energies}
\end{table*}

\clearpage

