\section{Appendix}\label{suppmaterial}

In this section we describe how the Hamiltonian is diagonalized, subsequently how one obtains the fermionic correlation functions and from them the spin correlations functions that are analyzed in the main manuscript.


\subsection{Diagonalization of the Hamiltonian}

Despite the presence of three-body interactions the Hamiltonian~(\ref{hamiltonian}) can be mapped into a delocalized fermionic form
via a generalization of the approach described in Ref.~\cite{Lieb1961}. To obtain this map we introduce the Jordan-Wigner transformation~\cite{Jordan1928} that connects the spin operators to the fermionic
operator via
\begin{equation}
 \label{JW_transformation}
 c_j=\left(\bigotimes_{i=-l-1}^{j-1} 2 S_i^z \right) S_j^- \;, \, \, \, \, \, \,
 c_j^\dagger=\left( \bigotimes_{i=-l-1}^{j-1} 2 S_i^z \right) S_j^+\, ,
\end{equation}
where $S_i^\pm=S_i^x\pm i S_i^y$. Applying this transformation to the Hamiltonian~(\ref{hamiltonian}) we obtain
\begin{eqnarray}
 \label{hamiltonian_2}
 H&=& \frac{1-\alpha}{8} \sum_{i=-l}^{l} (c_{i-1}^\dagger c_{i+1}^\dagger + c_{i-1}^\dagger c_{i+1} + h.c.)\nonumber \\
 & +&\frac{\alpha}{4} \sum_{i=-l-1}^{l} (c_{i}^\dagger c_{i+1}^\dagger +c_{i}^\dagger c_{i+1} + h.c)\;.
\end{eqnarray}
This  Hamiltonian corresponds to a so-called free fermion case in which non-interacting fermions can jump between nearest and next to nearest neighbors of the chain via open boundary conditions, i.e. $N=2l+3$.
This Hamiltonian can be brought in a diagonal form
\begin{equation}
 \label{hamiltonian_3}
 H= \sum_k \lambda_k\; \eta_k^\dagger \eta_k + const.
\end{equation}
using the following linear transformation
\begin{eqnarray}
\label{eta}
\eta_k&=& \sum_{i=-l-1}^{l+1} \frac{\phi_{k,i} + \psi_{k,i}}{2} c_i+  \frac{\phi_{k,i} - \psi_{k,i}}{2} c_i^\dagger \nonumber \\
\eta_k^\dagger&=& \sum_{i=-l-1}^{l+1} \frac{\phi_{k,i} + \psi_{k,i}}{2} c_i^\dagger+  \frac{\phi_{k,i} - \psi_{k,i}}{2} c_i \, ,
\end{eqnarray}
where both vectors $\phi_k$ and $\psi_k$ are a set of $N$ vectors with $N$ components that are solutions of the coupled equations
\begin{eqnarray}
\label{psipsh1}
\lambda_k \psi_k&=& \phi_k  (\mathbf{A}+\mathbf{B})\nonumber \\
\lambda_k \phi_k&=& \psi_k (\mathbf{A}-\mathbf{B})  \,.
\end{eqnarray}
Here the matrix $\mathbf{A}$ (symmetric) and $\mathbf{B}$ (anti-symmetric) are, respectively,
\begin{equation}
\label{matA}
\mathbf{A}\!=\!\frac{1}{8}\!\left(\!
\begin{array}{cccccccccc}
 0\!\!&2 \alpha\!\!& \gamma \!\!&0\!\!&0\!\!&0\!\!&0\!\!&0\!\!&0\!\!&0 \\
 2 \alpha\!\!& 0 \!\!& 2 \alpha\!& \gamma\!&0\!&0\!&0\!&0\!&0\!&0 \\
 \gamma\!\!& 2\alpha \!& 0 \!& 2 \alpha\!& \gamma\!&0\!&0 \!&0\!&0\!&0 \\
 \!\!& \ddots\!& \ddots\!& \ddots\!& \ddots\!& \ddots \!& \\
\!\!& \!& \ddots\!& \ddots\!& \ddots\!& \ddots\!& \ddots \!& \\
\!\!& \!& \!& \ddots\!& \ddots\!& \ddots\!& \ddots\!& \ddots \!& \\
\!\!& \!& \!& \!& \ddots\!& \ddots\!& \ddots\!& \ddots\!& \ddots \!& \\
 0\!\!&0 \!&0\!&0\!&0\!&\gamma\!& 2\alpha \!& 0 \!& 2 \alpha\!& \gamma\\
 0\!\!& 0\!&0 \!&0\!&0\!&0\!&\gamma\!& 2\alpha \!& 0 \!& 2 \alpha \\
 0\!\!&0\!& 0\!&0 \!&0\!&0\!&0\!&\gamma\!& 2\alpha \!& 0
\end{array}
\!\right)
\end{equation}
and
\begin{equation}
\label{matB}
\mathbf{B}\!=\!\frac{1}{8}\!\left(\!\!
\begin{array}{cccccccccc}
\!\! 0\!&2 \alpha\!&\! \gamma \!&\!0\!&\!0\!&\!0\!&\!0\!&\!0\!&\!0\!&\!0 \\
\!\! -2 \alpha\!&\! 0 \!&\! 2 \alpha\!&\! \gamma\!&\!0\!&\!0\!&\!0\!&\!0\!&\!0\!&\!0 \\
\!\! -\gamma\!&\! -2\alpha \!&\! 0 \!&\! 2 \alpha\!&\! \gamma\!&\!0\!&\!0 \!&\!0\!&\!0\!&\!0 \\
\!\! \!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots \!&\! \\
\!\!\!&\! \!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots \!&\! \\
\!\!\!&\! \!&\! \!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots \!&\! \\
\!\!\!&\! \!&\! \!&\! \!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots\!&\! \ddots \!&\! \\
\!\! 0\!&\!0 \!&\!0\!&\!0\!&\!0\!&\!-\gamma\!&\! -2\alpha \!&\! 0 \!&\! 2 \alpha\!&\! \gamma\\
\!\! 0\!&\! 0\!&\!0 \!&\!0\!&\!0\!&\!0\!&\!-\gamma\!&\! -2\alpha \!&\! 0 \!&\! 2 \alpha \\
\!\! 0\!&\!0\!&\! 0\!&\!0 \!&\!0\!&\!0\!&\!0\!&\!-\gamma\!&\! -2\alpha \!&\! 0
\end{array}
\!\!\right)
\end{equation}
and $\gamma=1-\alpha$.
Either $\psi_k$ or $\phi_k$ can be eliminated from the equations~(\ref{psipsh1}) resulting in
\begin{eqnarray}
\label{psipsh2}
\lambda_k^2 \psi_k&=& \psi_k (\mathbf{A}+\mathbf{B})(\mathbf{A}-\mathbf{B})\nonumber \\
\lambda_k^2 \phi_k&=& \phi_k (\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})  \,.
\end{eqnarray}
For $\lambda_k \neq 0$ we can choose one of these two equations that gives us one set of vectors and with the help of the euations.~(\ref{psipsh1}) one obtains the second set of vectors. For $\lambda_k = 0$ one obtains $\psi_k$ and $\phi_k$ from the equations~(\ref{psipsh1}).

\subsection{Fermionic correlation functions}

The set of vectors $\psi_k$ and $\phi_k$ plays an important role for the evaluation of the fermionic correlations functions.
To see that let us introduce two new fermionic operators
\begin{equation}
\label{otherfermionic}
 \nu _i= c^\dagger_i + c_i \,; \,\,\,\,\,\,\,\,\,\,\,\,\,\, \mu_i = c^\dagger_i - c_i \, .
\end{equation}
The advantage to use this set of operators, $\mu_i$ and
$\nu_i$, with respect to the set of operators $c_i$ and $c_i^\dagger$ comes from the fact that, taking into account the Eqs.~(\ref{eta}),
we obtain particular simplifications computing the expectation values of the product of two operators, in particular:
\begin{eqnarray}
\label{correlations}
 \langle \mu_i \mu_k\rangle = & \sum_l \phi_{l,i} \phi_{l,j} = & \delta_{ik} \nonumber \\
 \langle \nu_i \nu_k\rangle = & - \sum_l \psi_{l,i} \psi_{l,j} = & -\delta_{ik}  \\
 \langle \nu_i \mu_k\rangle = -\langle \mu_i \nu_k\rangle = & - \sum_l \psi_{l,i} \phi_{l,j} = & G_{i,j}(\alpha)\;. \nonumber
\end{eqnarray}
In the finite size case the $G_{i,j}(\alpha)$ depends on the choice of the site $i$ and $j$ and, therefore, also the spin correlation functions
depend on the position of the spins in the chain. However, if we evaluate the correlation functions between three central spins of
a chains with an odd $N$, the mirror symmetry plays the role of the invariance under spatial translation, imposing that
\begin{equation}
\langle S_{-1}^\alpha S_{0}^\beta \rangle=\langle S_{1}^\alpha S_{0}^\beta \rangle \quad \forall\; \alpha,\beta=x,y,z\;.
\end{equation}

On the contrary, in the thermodynamic limit the invariance under spatial translation is restored by the fact that the end of the chain is
infinitely far away. In this case, imposing the invariance under spatial translation, we can solve Eqs.~(\ref{psipsh1}) considering that
each components of $\phi_k$ and  $\psi_k$ vectors obey to the following rule
\begin{equation}
\phi_{k,l}=\frac{q_{k,l}}{\sqrt{N}}e^{-i \frac{2 \pi}{N} k l} \, \, ; \, \, \, \, \, \psi_{k,l}=\frac{q'_{k,l}}{\sqrt{N}}e^{-i \frac{2 \pi}{N} k l} \, .
\end{equation}
With a long but straightforward evaluation we obtain the fermionic correlation functions in the thermodynamic limit given in Eq.~(\ref{Gfunction}) and Eq.~(\ref{Lambda}), respectively.

\subsection{Spin correlation functions}

Having introduced a method to evaluate the fermionic correlation functions we can compute the spin correlation functions from which we obtain any
reduced density matrix.
We can divide the spin correlation functions in two families. The first is consists of correlation functions from operators that do not break the
symmetry of the Hamiltonian. In this case transforming the spin operators into the fermionic operators $\mu_i$ and $\nu_i$ one always obtains
an operator constructed via an equal number of $\mu_i$ and $\nu_i$. Applying the Wick's theorem we then obtain the expression of the correlation
functions in terms of the $G(n,\alpha)$. Let us also point out that the dependence of the spin correlation functions on the fermionic
correlation functions are independent of the model. A list of these spin correlation functions given in terms of fermionic correlation functions
can be founded in Ref.~\cite{Hofmann2013}.

The situation is much more involved in the case the symmetry of the
Hamiltonian is broken. In this second case, if we try to apply the same strategy as described above, we immediately observe that
the elements due to parity immediately vanish. In fact, re-writing the operators in terms of $\mu_i$ and $\nu_i$, we obtain an
expression with different numbers of fermionic operators, each one referring to distinct sites, and, in agreement with Eq.~(\ref{correlations}),
this implies that the correlation function must be zero. This result holds for all ground states that preserve the symmetry of the
Hamiltonian, however, in the thermodynamic limit and in the anti-ferromagnetic phase one has also ground states that break the $Z_2$
symmetry of the Hamiltonian. For these ground states some of these correlation functions are expected to be different from zero.

\begin{figure}[t]
{\includegraphics[width = 7cm]{supp1.pdf}}\\
{\includegraphics[width = 7cm]{supp2.pdf}}\\
{\includegraphics[width = 7cm]{supp3.pdf}}
\caption{(Color online) The figures show the values of $\langle \sigma_i^y \rangle_B $ (red points) and
$\langle \sigma_{i-1}^y \sigma_{i}^y \sigma_{i+1}^y \rangle_B$ (black squares) for different distances $r$ and for
different values of $\alpha$ in the anti-ferromagnetic phase.}
\label{plotsupp1}
\end{figure}

To evaluate these correlation functions we adopt the approach of Barouch and McCoy~\cite{BarouchMcCoy1971}. Given an operator
$O_{\{k\}}$ defined on a set of spins $\{k\}$ that does not commute with the parity operator along the $Z$ axis one can define the correlation function
$\langle O_{\{k\}}\rangle_B$ via
\begin{equation}
 \label{corr-break}
 \langle O_{\{k\}} \rangle_B=\lim_{r \rightarrow \infty} \sqrt{\langle O_{\{k\}} O_{\{k\}+r}\rangle}\;.
\end{equation}

For analyzing the entanglement properties of three adjacent spins one
has to compute all the one, two and three body correlation functions that do not commute with the parity operators along the $Z$ axis. All
the correlation functions of this type are vanishing in the limit of $r \rightarrow \infty$ except $\langle \sigma_i^y \rangle_B $ and
$\langle \sigma_{i-1}^y \sigma_{i}^y \sigma_{i+1}^y \rangle_B$. In Fig.~\ref{plotsupp1} we have reported, for some values of $\alpha$, the values
of  $\langle \sigma_i^y \rangle_B $ and $\langle \sigma_{i-1}^y \sigma_{i}^y \sigma_{i+1}^y \rangle_B$. We observe that the farer we are away from the quantum critical point the faster is the convergence of the correlation functions to the asymptotic values. Moreover, for all $\alpha$ we observe that $\langle \sigma_i^y \rangle_B = \langle \sigma_{i-1}^y \sigma_{i}^y \sigma_{i+1}^y \rangle_B$ for large enough $r$
(we have tested this relations for more that $1000$ different value of $\alpha>\alpha_c$). The fact that this equality holds impose that the
genuine tripartite entanglement between three adjacent spins does not depends on the particular ground state selected as emphasized in the main body of the manuscript.

\end{document}

