text:
\section*{Long-distance entanglement of purification and reflected entropy\\in conformal field theory -- supplemental material}

\section*{Appendix}

\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\theequation}{S\arabic{equation}}

\mysection{Review of critical Ising model} The Hamiltonian of the transverse Ising model is given by
\begin{equation}\label{eq:TIsingHam}
\hat{H} = - \sum_{k=1}^N \left( 2 J\, \hat{S}^x_k \hat{S}^x_{k+1} + h\, \hat{S}^z_k \right) \ ,
\end{equation}
with spin operators represented by Pauli matrices $\sigma_\alpha$ with $\alpha \in (x,y,z)$ by
\begin{equation}
\hat{S}^\alpha_k \equiv \id^{\otimes (k-1)} \otimes \frac{\sigma_\alpha}{2} \otimes \id^{\otimes (N-k)} \ .
\end{equation}
We also use the identification $\hat{S}^\alpha_{N+1} \equiv \hat{S}^\alpha_1$.
This spin model can be converted to fermions by defining the $2N$ Majorana operators $\operatorname{\gamma}_k$ via
\begin{align}
\operatorname{\gamma}_{2k-1} &= {\sigma_z}^{\otimes (k-1)} \otimes \sigma_x \otimes \id^{\otimes (N-k)} \ , \\
\operatorname{\gamma}_{2k} &= {\sigma_z}^{\otimes (k-1)} \otimes \sigma_y \otimes \id^{\otimes (N-k)} \ .
\end{align}
The Ising Hamiltonian then takes the form
\small
\begin{equation}
\hspace{- 10 pt}\hat{H}_\text{I} = \frac{\mathrm{i}}{2} \left( \operatorname{\gamma}_1 \operatorname{\gamma}_{2N} P + J \sum_{k=1}^{N-1} \operatorname{\gamma}_{2k} \operatorname{\gamma}_{2k+1} + h \sum_{k=1}^{N} \operatorname{\gamma}_{2k-1} \operatorname{\gamma}_{2k} \right) .
\end{equation}
\normalsize
Here $P$ is the total parity operator ${\sigma_z}^{\otimes N}=\prod_k (-\mathrm{i} \operatorname{\gamma}_{2k-1}\operatorname{\gamma}_{2k})$.
At the critical point $J=h$, the Hamiltonian thus simplifies to
\begin{equation}
\hat{H} = \frac{\mathrm{i}}{2} \left( \operatorname{\gamma}_1 \operatorname{\gamma}_{2N} P + \sum_{k=1}^{2N-1} \operatorname{\gamma}_k \operatorname{\gamma}_{k+1} \right)\,,
\end{equation}
which leads for $N\to\infty$ to the lattice model of the $c=\tfrac{1}{2}$ CFT. The critical Ising Hamiltonian as displayed in the main text~\eqref{eq.HIsing}
corresponds to $J = h = 1$.

\mysection{Covariance matrix} For the critical ground state vector $\ket{0}$ which has a positive total parity, all correlations are encoded in the Majorana covariance matrix
\begin{equation}
\Omega_{j,k} = \frac{\mathrm{i}}{2} \bra{0} [\operatorname{\gamma}_j, \operatorname{\gamma}_k] \ket{0} \ ,
\end{equation}
which in the infinite system size limit takes the form
\begin{equation}
\Omega_{j,k} =
\begin{cases}
0 &  k = j  \\
\frac{(-1)^{k-j} - 1}{\pi (k-j)} &  k \neq j
\end{cases} \ .\label{eq:covariance-inf}
\end{equation}
The entropy of a Gaussian mixed state $\rho$ with covariance matrix $\Omega$ is given by
\begin{align}
    S(\rho)=-\sum_{\pm, i}\frac{1\pm\lambda_i}{2}\log\frac{1\pm\lambda_i}{2}\,,\label{eq:Gaussain-SA}
\end{align}
where $\pm\mathrm{i} \lambda_k$ are the eigenvalues of $\Omega$. We will consider mixed states $\rho_{AB}$ or $\rho_{AA'}$, whose mixed state covariance matrices $\Omega_{AB}$ or $\Omega_{AA'}$ result from restricting~\eqref{eq:covariance-inf} to the respective blocks.

\mysection{Fermionic subsystem}
\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{Fig3.pdf}
    \caption{Subsystem setup of our analytical limits for fermions (top) with an inherent ordering and spins (bottom) without one. In both systems, we consider the subsystem $AB$ consisting of two single sites $A$ and $B$ separated by $d/\delta$ sites.
    }
    \label{fig:1+1subsystem}
\end{figure}
We first compute reduced density matrices from the perspective of fermions, {\it i.e.},\  imposing an ordering between modes following from the anti-commuting variables $\gamma_i$ (see Fig.~\ref{fig:1+1subsystem}). A subsystem consisting of $1+1$ sites ($w=\delta$) separated by $d/w=d/\delta$ sites is then fully characterized by the restriction of the covariance matrix in~\eqref{eq:covariance-inf} and explicitly given by
\begin{align}
    \Omega^{\mathrm{fer}}_{AB} =
    \begin{pmatrix}
      & -\frac{2}{\pi } &  & \frac{-2}{(2 d/w+3) \pi } \\
     \frac{2}{\pi } &  & \frac{-2}{(2 d/w+1) \pi } &  \\
      & \frac{2}{(2 d/w+1) \pi } &  & -\frac{2}{\pi } \\
     \frac{2}{(2 d/w+3) \pi } &  & \frac{2}{\pi } &
    \end{pmatrix}\label{eq:omega-fermions}
\end{align}
which corresponds to a lowest-dimension (Majorana) operator with scaling dimension $\Delta=\sfrac{1}{2}$. The associated fermionic density operator is then
\begin{align}\label{eq:rho-fermions}
    \hspace{-2mm}\rho^{\mathrm{fer}}_{AB}\sim\begin{pmatrix}
     D & & &\frac{1}{2\pi}\,\epsilon_{\sfrac{1}{2}}\\
     & E & & \\
     & & E & \\
     \frac{1}{2\pi}\,\epsilon_{\sfrac{1}{2}} & & & F
    \end{pmatrix}
\end{align}
with respect to the basis $(\ket{\downarrow\downarrow}$, $\ket{\uparrow\downarrow},\ket{\downarrow\uparrow},\ket{\uparrow\uparrow})$ and using $D=\frac{1}{4}+\frac{1}{\pi}+\frac{1}{\pi^2}$, $E=\frac{1}{4}-\frac{1}{\pi^2}$, and $F=\frac{1}{4}-\frac{1}{\pi}+\frac{1}{\pi^2}$.
As in the main text, $\epsilon_\Delta \equiv (w/d)^{2\Delta}$ which here becomes $\epsilon_{\sfrac{1}{2}}=w/d$.
If we further restrict to a single site, we find the covariance matrix and density operator
\begin{align}
    \Omega^{\mathrm{fer}}_{A} &=
    \begin{pmatrix}
      & -\frac{2}{\pi }  \\
     \frac{2}{\pi } &
    \end{pmatrix} & \rho_A^{\mathrm{fer}} &=
    \begin{pmatrix}
     \frac{1}{2}-\frac{1}{\pi} & \\
      & \frac{1}{2}+\frac{1}{\pi}
    \end{pmatrix}\,,
\end{align}
where $\rho_A^{\mathrm{fer}}$ is written with respect to the basis $(\ket{\downarrow},\ket{\uparrow})$.

\mysection{Spin subsystem} We can perform a similar calculation in the original Ising spin system whose reduced density matrices can be constructed from the fermionic covariance matrix \cite{Coser:2015mta}. We need not repeat the single interval case, as entanglement entropies of connected regions are equivalent under a Jordan-Wigner transformation. However, we still need the reduced density matrix of a system of $1+1$ sites in the large $d$ limit, which we find to be
\begin{equation}
\label{eq:rhoAB}
    \rho^{\mathrm{spin}}_{AB} \sim
    \begin{pmatrix}
      D &  &  & C\epsilon_{\sfrac{1}{8}} \\
  & E & C\epsilon_{\sfrac{1}{8}} & \\
  & C\epsilon_{\sfrac{1}{8}} & E &  \\
 C\epsilon_{\sfrac{1}{8}} &  &  & F
    \end{pmatrix}
\end{equation}
with $w/\delta=1$ for the setup of $1+1$ sites. As we are considering the spin Ising CFT, the lowest-dimension primary is the ``order field'' $\sigma$ with scaling dimension $\Delta=\sfrac{1}{8}$.
The constant $C$ corresponds to the expectation value of an operator nonlocal in fermions, and can be computed from
\begin{equation}
    C = \lim_{n \to \infty} \left( \frac{2}{\pi} \right)^n \frac{n^{1/4}}{4} \det{M^n} \ ,
\end{equation}
where $M^n$ is defined as the $n \times n$ matrix
\begin{equation}
    M^n_{j,k} =
    \begin{cases}
    \frac{(-1)^{k-j}}{2(k-j)+1} & j \leq k \\
    \frac{(-1)^{j-k+1}}{2(j-k)-1} & j > k
    \end{cases} \ .
\end{equation}
Using this construction, one finds~\cite{McCoy:2010}
\begin{equation}
\label{eq:Cdefinition}
    C =\frac{e^{3\zeta'(-1)}}{2^{23/12}}\approx 0.1612506\,.
\end{equation}

\mysection{MI for fermions}
We now begin computing entanglement measures for the small subsystems whose reduced density matrices we just obtained explicit expressions for.
The continuum limit corresponding to the Ising CFT is obtained by keeping $d/w$ (or, equivalently, $\epsilon_\Delta$) fixed and taking $\delta/w$ to $0$. We will see that taking only a few lattice sites is sufficient to describe the qualitative and approximate quantitative behavior of the continuum limit.
To demonstrate this, we now show that the large distance asymptotics of the MI ~\eqref{eq.IAB} for the case of $1+1$ sites yield results close to the continuum formula~\eqref{eq.MIcftlarged}.
In order to compute the MI, we need to determine the von Neumann entropy of a single site $S_A=S_B$ and of both sites $S_{AB}$. These entropies can be computed from the eigenvalues of the covariance matrix $\Omega$ associated to the respective Gaussian state $\rho$. As an antisymmetric matrix, $\Omega^{\mathrm{fer}}_{AB}$  has pairs of purely imaginary eigenvalues $\pm \mathrm{i}\, \lambda_k$, from which applying~\eqref{eq:Gaussain-SA} leads to
\begin{align}
S_A&=- \tfrac{\pi+2}{2\pi}\log\tfrac{\pi+2}{2\pi} - \tfrac{\pi-2}{2\pi}\log\tfrac{\pi-2}{2\pi} \approx 0.474\\
    S_{AB} &=\sum_{k=1}^n \left( - \tfrac{1+\lambda_k}{2}\log\tfrac{1+\lambda_k}{2} - \tfrac{1-\lambda_k}{2}\log\tfrac{1-\lambda_k}{2} \right)\,,\label{eq:Gaussian-SAB}
\end{align}
where the eigenvalues of $\Omega_{AB}^{\mathrm{fer}}$ are to leading order
\begin{equation}
    \lambda_{1,2} = \frac{1}{\pi} \left( 2 \pm \tfrac{3}{4}\epsilon_{\sfrac{1}{2}}^2  + \dots \right) \ .
\end{equation}
We can similarly expand $S_{1+1} \equiv S_{AB}$ at large $d$, which results in a MI for $w=\delta$ of
\begin{align}\label{eq:MIFerm11}
    I^{\mathrm{fer}}(A:B)&\sim \tfrac{\log\frac{\pi+2}{\pi-2}}{4\pi}\,\epsilon_{\sfrac{1}{2}}^2 \nonumber\\
    &= 0.120 \,\left(\tfrac{w}{d}\right)^2 \ .
\end{align}
This reproduces the correct continuum power law of fermionic MI, but yields a coefficient lower than the continuum value \eqref{eq.MIcftlarged} which also matches the large-distance expansion of earlier results for Dirac fermions~\cite{Casini:2005rm}
\begin{align}\label{eq:MIFerm}
    I(A:B)=\frac{c}{3}\log\frac{(d+w)^2}{d\,(2w+d)} &\sim \frac{1}{6} \left(\frac{w}{d}\right)^2 \ ,
\end{align}
for two blocks of general width $w$.

\mysection{MI for spins}
We compute EE for the spin system directly from the eigenvalue spectrum of the reduced density matrix \eqref{eq:rhoAB}. Its four eigenvalues $\mu_j$ are
\begin{subequations}
\begin{align}\label{eq:EigenSpinRhoAB}
    \mu_{1,2} &= \frac{1}{4} - \frac{1}{\pi} \pm C \,\epsilon_{\sfrac{1}{8}} \ , \\
    \mu_{3,4} &= \frac{1}{4} + \frac{1}{\pi} \pm \sqrt{\frac{1}{\pi^2} + C^2 \, \epsilon_{\sfrac{1}{8}}^2}\ ,
\end{align}
\end{subequations}
where $C$ is given by \eqref{eq:Cdefinition} and from which we can directly compute the EE for $AB$ via
\begin{equation}
\label{eq.Smu}
S = - \sum_j \mu_j \log\mu_j.
\end{equation}
Note that in the following we will denote eigenvalues of any density matrix by $\mu_{j}$.

This analysis leads to the Ising model prediction for spin MI at $w=\delta$ and large separations $d$ of the form
\begin{align}
\label{eq.MIanalytic}
 \hspace{-2mm}   I^\text{spin}(A:B) &\sim
 C^2 \left( \frac{4\pi^2}{\pi^2-4} + \frac{\pi}{2}\log\frac{4+4\pi+\pi^2}{4-4\pi+\pi^2} \right)  \epsilon_{\sfrac{1}{8}}^2 \nonumber\\
 &\sim 0.298\, \sqrt{\frac{w}{d}} \ ,
\end{align}
\normalsize
Comparing this formula with the CFT analytics~\eqref{eq.MIcftlarged} for $\Delta = \sfrac{1}{8}$, we see an exact match in the power-law behavior. Furthermore, the prefactor in~\eqref{eq.MIanalytic} is only $3.6\%$ off from the continuum value $\approx0.309$ predicted by~\eqref{eq.MIcftlarged}.

\mysection{EoP for fermions} Analogous to the MI calculation for free fermions, we now calculate the EoP in the fermionic subsystem of two sites separated by $d/\delta$ sites, expressing all calculations in terms of covariance matrices. We purify $\Omega_{AB}$ in the limit $d/\delta\to\infty$ as
\begin{align}
    \Omega^{(0)}=\left(\begin{array}{cccc|cccc}
      & -G & & &  & L& & \\
    G &    & & & L &  & &\\
      &  & & -G & & & & L\\
     &    & G & & &  & L&\\
     \hline
       & -L & & &  & -G& & \\
    -L &    & & & G & & &\\
      &  & & -L & & &  & -G\\
     &    & -L & & &  & G&\\
    \end{array}\right)
\end{align}
associated to systems $(A,B,A',B')$ with $G=\frac{2}{\pi}$ and $L=\sqrt{1-G^2}$, whose EE $S_{AA'}$ is zero and we thus have $\lim_{d\to\infty}E_P=0$, {\it i.e.},\  the EoP vanishes for large $d/\delta$, as expected.

In order to find the asymptotic behavior of $E_P$, we need to study the variation of the symplectic eigenvalues $\pm\mathrm{i}\lambda_i$ of $\Omega_{AA'}$ when perturbing $\Omega$ according to
\begin{align}
    \Omega\sim \Omega^{(0)}+\epsilon_{\sfrac{1}{2}}\,\Omega^{(1)}+\tfrac{1}{2}\epsilon_{\sfrac{1}{2}}^2\,\Omega^{(2)}\quad\text{as}\quad \epsilon_{\sfrac{1}{2}}\to 0\,.
\end{align}
The requirement of $\Omega$ representing a purification implies $\Omega^2=-\id$, which induces the constraints
\begin{align}
\begin{split}
    \Omega^{(0)}\Omega^{(1)}+\Omega^{(1)}\Omega^{(0)}&=0\,,\\
    2(\Omega^{(1)})^2+\Omega^{(0)}\Omega^{(2)}+\Omega^{(2)}\Omega^{(0)}&=0\,.
\end{split}\label{eq:constraints-J}
\end{align}
We further require that the restrictions $\Omega^{(1)}_{AB}$ and $\Omega^{(2)}_{AB}$ matches the ones of\ \eqref{eq:omega-fermions} expanded in $\epsilon_{\sfrac{1}{2}}$, {\it i.e.},\
\begin{align}
    \Omega^{(1)}_{AB}\!=\!\begin{pmatrix}
    &&& -\frac{1}{\pi}\\
    &&-\frac{1}{\pi}& \\
    &\frac{1}{\pi}&& \\
    \frac{1}{\pi}&&&
    \end{pmatrix},\,\,
    \Omega^{(2)}_{AB}\!=\!\begin{pmatrix}
    &&& \frac{3}{\pi}\\
    &&\frac{1}{\pi}& \\
    &-\frac{1}{\pi}&& \\
    -\frac{3}{\pi}&&&
    \end{pmatrix}\,,\label{eq:constraints2-J}
\end{align}
The equations\ \eqref{eq:constraints-J} and\ \eqref{eq:constraints2-J} can be solved iteratively up to some free variables. We first solve $\Omega^{(1)}$ in terms of $\Omega^{(0)}$ and then $\Omega^{(2)}$ in terms of $\Omega^{(0)}$ and $\Omega^{(1)}$.

In order to find asymptotics of the symplectic eigenvalues $\lambda_i$, we can use the fact that $\operatorname{Tr}(\Omega_{AA'}^2)=-2(\lambda_1^2+\lambda_2^2)$ and $\operatorname{Tr}(\Omega_{AA'}^4)=2(\lambda_1^4+\lambda_2^4)$ to solve for the asymptotics of $\lambda_i$ to be given by
\begin{align}
    \lambda_1=\lambda_2\sim 1-\alpha_{\mathrm{tot}}\,\epsilon_{\sfrac{1}{2}}^2\quad \text{as}\quad \epsilon_{\sfrac{1}{2}}\to 0\,,\label{eq:symp-eig}
\end{align}
where $\alpha_{\mathrm{tot}}$ will depend on some of the free parameters contained in $\Omega^{(1)}$ and $\Omega^{(2)}$. With this trick, one finds
\begin{align}
\begin{split}
    \alpha_{\mathrm{tot}}&=\frac{x_{14}a_{23}-x_{13}x_{24}+\pi^{-2}}{2}+\frac{G(x_{14}-x_{23})\pi^{-1}}{2L}\\
    &\quad +\frac{(x_{14}-x_{23})^2+(x_{13}+x_{24})^2}{4L^2}\,,\label{eq:alpha2}
\end{split}
\end{align}
where the variables $x_{ij}$ represent unconstrained entries in the block $\Omega^{(1)}_{AB,A'B'}$. In order to find the asymptotics of EoP, we need to minimize $\alpha_{\mathrm{tot}}$ over these parameters to find the smallest possible EE $S_{AA'}$. Due to the fact that~\eqref{eq:alpha2} is quadratic in $x_{ij}$, we can calculate this valua analytically as
\begin{align}\label{eq:alphatotEoPFerm}
    \alpha_{\mathrm{tot}}&=\frac{1}{8+2\pi^2}\approx 0.03605\,.
\end{align}
Expanding $S_{AA'}\sim \sum_i(\log{2}-\tfrac{\lambda_i}{2})$ through $\lambda_i$ up to second order in $\epsilon_{\sfrac{1}{2}}$ based on~\eqref{eq:symp-eig} allows us to also find the offset analytically, namely we have
\begin{align}
    S_{AA'}=\epsilon_{\sfrac{1}{2}}^2\left(\alpha_{\mathrm{tot}}\log(\epsilon_{\sfrac{1}{2}}^{-2})+\alpha_{\mathrm{tot}} \log{\tfrac{2e}{\alpha_{\mathrm{tot}}}}\right)\,.
\end{align}
Combining this with the result from~\eqref{eq:alphatotEoPFerm} gives
\begin{align}
\label{eq:eopFerm11}
    E_P^\text{fer}&\sim \left(\frac{1}{8+2\pi^2} \log(\epsilon_{\sfrac{1}{2}}^{-2})+ \frac{\log 2e (8+2\pi^2)}{8+2\pi^2} \right)\, \epsilon_{\sfrac{1}{2}}^2 \nonumber\\
    &\sim \left(0.0361\log\left(\tfrac{d}{w}\right)^2 +0.181\right)\, \left(\tfrac{w}{d}\right)^2  \ ,
\end{align}
which agrees with the form~\eqref{eq.SAAp} in the main text. Note that the simplicity of Gaussian states allowed us to even find the analytical form of the constant offset. The accuracy of this analytical prediction was tested numerically, for which we presented the results in Fig.~\ref{fig:MI_EoP} in the main text.

\mysection{EoP for spins} In the limit of an infinite distance between the two single site subsystems, we purify~\eqref{eq:rhoAB} by the state $\ket{\psi^{(0)}}$ with Schmidt decomposition
\begin{align}
\hspace{-7 pt}    \ket{\psi^{(0)}}=\sqrt{D}\ket{\downarrow\downarrow\downarrow\downarrow}+\sqrt{E}(\ket{\uparrow\downarrow\uparrow\downarrow}+\ket{\downarrow\uparrow\downarrow\uparrow})+\sqrt{F}\ket{\uparrow\uparrow\uparrow\uparrow},\label{eq:schmidt}
\end{align}
\normalsize
where the convention for factors ordering in the purification is $ABA'B'$. Note that in this analysis we assume that a minimal purification from two to four spin degrees of freedom suffices and we will subsequently provide supporting numerical evidence and an additional discussion.

Moving on, we \emph{supplement} this purification with finite distance corrections up to second order in $\epsilon_{\sfrac{1}{8}}$ as
\begin{align}
    \ket{\psi}\sim \ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{8}} \ket{\psi^{(1)}}+\tfrac{1}{2}\epsilon_{\sfrac{1}{8}}^2 \ket{\psi^{(2)}}\,.
\end{align}
\normalsize
We will optimize over $\ket{\psi^{(1)}}$ and $\ket{\psi^{(2)}}$ subject to the normalization constraint $\braket{\psi|\psi}=1$ order by order in $\epsilon_{\sfrac{1}{8}}$. We further require $\rho^{(1)}=\ket{\psi^{(0)}}\bra{\psi^{(1)}}+\ket{\psi^{(1)}}\bra{\psi^{(0)}}$ and $\rho^{(2)}=\ket{\psi^{(0)}}\bra{\psi^{(2)}}+\ket{\psi^{(2)}}\bra{\psi^{(0)}}+2\ket{\psi^{(1)}}\bra{\psi^{(1)}}$ to satisfy
\begin{align}
    \rho^{(1)}_{AB}=\begin{pmatrix}
     &&&C\\
     &&C&\\
     &C&&\\
     C&&&
    \end{pmatrix}\quad\text{and}\quad \rho^{(2)}_{AB}=0\,,\label{eq:nG-constaints2}
\end{align}
which follows from\ \eqref{eq:rhoAB}. We expand the first order perturbation as
\begin{align}
    \ket{\psi^{(1)}}=C\sum^{16}_{i=1} z_i\ket{\phi_i}
\end{align}
where $z_i=x_i+\mathrm{i} y_i$ and $\ket{\phi_i}$ is the basis of $\mathcal{H}_{ABA'B'}$ ordered as $(\ket{\downarrow\downarrow\downarrow\downarrow},\ket{\uparrow\downarrow\downarrow\downarrow},\ket{\downarrow\uparrow\downarrow\downarrow},\ket{\uparrow\uparrow\downarrow\downarrow},\dots,\ket{\uparrow\uparrow\uparrow\uparrow})$. We then need to implement the condition~\eqref{eq:affine-lin} in the main text based on~\eqref{eq:nG-constaints2} together with the normalization constraint~\eqref{eq:constraint1} in the main text. We solve these affine linear constraints by the replacements $x_1=x_6=x_{11}=x_{16}=0$, $z_{5}=-\sqrt{\frac{E}{D}}z_2^*$, $z_{9}=-\sqrt{\frac{E}{D}}z_3^*$, $z_{13}=\frac{1}{\sqrt{D}}-\sqrt{\frac{F}{D}}z_4^*$, $z_{15}=\sqrt{\frac{F}{E}}z_{12}^*$, $z_{10}=\frac{1}{\sqrt{E}}-z_{7}^*$, $z_{14}=-\sqrt{\frac{F}{E}}z_8$.
We can then compute $\alpha_{\mathrm{tot}}$ according to~\eqref{eq:alpha} in the main text as quadratic polynomial in terms of the remaining free variables $z_i$ which leads to the rather involved expression
\begin{widetext}
\begin{align}
\begin{split}

    \hspace{-3mm}\frac{\alpha_{\mathrm{tot}}}{C^2}&
    =\tfrac{(\pi -2)^2 y_1^2}{4 \pi ^2}-\tfrac{2 \sqrt{\pi ^2-4} x_3 x_8}{2+\pi }-\tfrac{2 \sqrt{\pi ^2-4} x_2 x_{12}}{2+\pi }-\tfrac{4 (\pi -2) \pi  x_4}{(2+\pi )^2}-\tfrac{4 \pi  x_7}{\sqrt{\pi ^2-4}}+\tfrac{(\pi -2) \left(-\sqrt{\pi ^2-4} y_6-\sqrt{\pi ^2-4} y_{11}+(2+\pi ) y_{16}\right) y_1}{2 \pi ^2}\\
    &\quad+\tfrac{\left(\pi ^2-4\right) y_6^2}{4 \pi ^2}+\tfrac{\left(\pi ^2-4\right) y_{11}^2}{4 \pi ^2}+\left(\tfrac{1}{4}+\tfrac{1}{\pi ^2}+\tfrac{1}{\pi }\right) y_{16}^2+\tfrac{(2+\pi ) y_6 \left((\pi -2) y_{11}-\sqrt{\pi ^2-4} y_{16}\right)}{2 \pi ^2}-\tfrac{(2+\pi ) \sqrt{\pi ^2-4} y_{11} y_{16}}{2 \pi ^2}\\
    &\quad-\tfrac{2 \sqrt{\pi ^2-4} y_3 y_8}{2+\pi }-\tfrac{2 \sqrt{\pi ^2-4} y_2 y_{12}}{2+\pi }+\tfrac{(\pi -2) |z_2|^2}{2+\pi }+\tfrac{(\pi -2) |z_3|^2}{2+\pi }+\tfrac{2 \left(4+\pi ^2\right) |z_4|^2}{(2+\pi )^2}+2 |z_7|^2+|z_8|^2+|z_{12}|^2+\tfrac{8 \pi ^3}{(\pi -2) (2+\pi )^2}\,.
\end{split}
\end{align}
\end{widetext}
In order to find the EoP, we need to minimize over the $z_i$ to find the smallest possible value $\alpha_{\mathrm{tot}}$, which can be done \emph{analytically} and leads to
\begin{align}
\label{eq.alphatot}
    \alpha_{\mathrm{tot}}=\frac{4\pi^4C^2}{\pi^4-16}\approx 0.12445\,.
\end{align}
The non-vanishing $\alpha_{\mathrm{tot}}$ shows that the resulting EoP obtained from~\eqref{eq.SAAp} again has the form
\begin{align}
\label{EQ_EOP_SPIN_W1}
    E_P^\text{spin} &\sim \left(  \frac{4\pi^4C^2}{\pi^4-16} \log(\epsilon_{\sfrac{1}{8}}^{-2}) + \mathrm{const} \right) \epsilon_{\sfrac{1}{8}}^{2} \nonumber\\
    &\sim \left( 0.124 \log\sqrt{\frac{d}{w}} + 0.440 \right) \sqrt{\frac{w}{d}} \ ,
\end{align}
which, as in the fermion case, exhibits a leading-order long-distance behavior enhanced with respect to that of MI~\eqref{eq.MIanalytic} by a logarithm of the distance.

When it comes to the \emph{subleading} long-distance behavior encapsulated by~$\Big( \sum_{j>0}\alpha_j(1-\log{\alpha_j}) \Big)$, we would need to extract the individual $\alpha_j$ and optimize over the remaining parameters. While it is plausible this can be also done analytically, we determined the value quoted above numerically, as discussed in the main text. Note that for the free fermion case with $w=\delta$ considered above, we determined this term analytically in terms of~$\alpha_{\mathrm{tot}}$.

\mysection{RE for fermions} In the Gaussian case of free fermions, our starting point is the following perturbative expansion of the reduced density matrix $\rho_{AB}$ of a system of $1+1$ fermions in the large $d$ separation, $\rho_{AB}\sim\rho^{(0)}_{A}\otimes\rho^{(0)}_{B}+\epsilon_{1/2}\rho^{(1)}_{AB}$ given by~\eqref{eq:rho-fermions}. We similarly construct the canonical purification of~\eqref{eq:rho-fermions} via $\ket{\sqrt{\rho_{AB}}}=\sum_{i}\sqrt{e_{i}}\ket{e_{i}}\otimes\ket{e_{i}}=\ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{2}}\, \ket{\psi^{(1)}}$ where $\rho_{AB}\ket{e_{i}}=e_{i}\ket{e_{i}}$. Note that in contrast with fermionic MI and EoP, we do not need to phrase our computation of RE in terms of the covariance matrix formalism since we can construct the canonical purification $\ket{\sqrt{\rho_{AB}}}$exactly for the given form of the initial reduced density matrix~$\rho_{AB}$.

In this case, the first-order perturbation $\ket{\psi^{(1)}}$ is simply given by
\begin{align}\label{EQ_PSI1_FERMION}
    \ket{\psi^{(1)}}=\tfrac{1}{2\pi}(\ket{\phi_{4}}+\ket{\phi_{13}}) \, ,
\end{align}
with the same ordering of the basis $\ket{\phi_{i}}$ as in the previous case. From the canonical purification's density matrix $\rho:=\ket{\sqrt{\rho_{AB}}}\bra{\sqrt{\rho_{AB}}}$ we consider a restriction to subsystems $AA'$ given by the reduced density matrix $\rho_{AA'}=\textrm{tr}_{BB'}(\rho)$ which has the perturbative expansion $\rho_{AA'}=\rho^{(0)}_{A}\otimes\rho^{(0)}_{A'}+\epsilon^{2}_{\sfrac{1}{2}}\rho^{(2)}_{AA'}$ explicitly given by
\begin{align}
\label{EQ_RHO_AA_FERMION}
   \rho_{AA'} \sim \rho^{(0)}_{AA'}+\tfrac{1}{2}\epsilon^{2}_{\sfrac{1}{2}}\rho^{(2)}_{AA'}=
   \begin{pmatrix}
  \tilde{G}_{1} &  &  & \tilde{H} \\
  & \tilde{J} &  & \\
  &  &   \tilde{J}&  \\
 \tilde{H} &  &  & \tilde{G}_{2}
    \end{pmatrix} ,
\end{align}
where $\tilde{G}_{1}=\frac{\pi+2}{2\pi}-\frac{\epsilon^{2}_{\sfrac{1}{2}}}{4\pi^{2}}$, $\tilde{G}_{2}=\tilde{G}_{1}-\frac{2}{\pi}$, $\tilde{H}=\frac{\sqrt{\pi^{2}-4}}{2\pi}-\frac{\sqrt{\pi^{2}-4}\epsilon^{2}_{\sfrac{1}{2}}}{4\pi(\pi^{2}-4)}$, and $\tilde{J}=\frac{\epsilon^{2}_{\sfrac{1}{2}}}{4\pi^{2}}$. We once again compute the trace of the square of~\eqref{EQ_RHO_AA_FERMION} according to~\eqref{eq:alpha} in the main text from which we obtain
\begin{align}\label{EQ_APLHA_TOT_RE_FERMION}
    \alpha_{\mathrm{tot}}=\frac{1}{2\pi^{2}}\approx 0.051 \, .
\end{align}
which also shows that the reflected entropy $S_{R}(\rho_{AB})=S_{AA'}(\rho)$ of the fermionic subsystem also exhibits a logarithmic enhancement of the power law decay for $w=\delta$ given by
\begin{align}\label{EQ_RE_FERMION_W1}
    S_{R}^\text{fer}(\rho_{AB}) &\sim
    \left( \frac{1}{2\pi^{2}} \log\epsilon_{\sfrac{1}{2}}^{-2} + \frac{1 + \log(4\pi^2)}{2\pi^2} \right) \epsilon_{\sfrac{1}{2}}^2 \nonumber  \\
    & \sim\left(0.051\log\left(\frac{d}{w}\right)^2+0.237\right) \left(\frac{w}{d}\right)^{2}\, ,
\end{align}
where we also computed the constant term in~\eqref{EQ_RE_FERMION_W1} from the eigenvalues of~\eqref{EQ_RHO_AA_FERMION} according to~\eqref{eq:alpha} in the main text.

\mysection{RE for spins} For the Ising spin case, we now describe the detailed computation of the reflected entropy RE for $w=\delta$ in the large $d$ limit just as for fermions. The reduced density matrix for a spin system of $1+1$ sites in the large $d$ limit can again be computed according to~\eqref{eq.rhogeneral}, {\it i.e.},\  $\rho_{AB} \sim \rho^{(0)}_{A} \otimes \rho^{(0)}_{B}+\epsilon_{\sfrac{1}{8}}\, \rho^{(1)}_{AB}+\ldots$, yielding~\eqref{eq:rhoAB}.

We now construct the canonical purification of~\eqref{eq:rhoAB} via $\ket{\sqrt{\rho_{AB}}}=\sum_{i}\sqrt{e_{i}}\ket{e_{i}}\otimes\ket{e_{i}}=\ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{8}}\ket{\psi^{(1)}}$ for $\rho_{AB}\ket{e_{i}}=e_{i}\ket{e_{i}}$ and where the eigenvalues $e_{i}$ are defined in~\eqref{eq:EigenSpinRhoAB}. In this case, the first order perturbation $\ket{\psi^{(1)}}$ is given by
\begin{align} \label{EQ_PSI1_ISING}
\begin{split}
    \ket{\psi^{(1)}}&= \tfrac{\pi}{\sqrt{\pi^{2}-4}}(\ket{\phi_{7}}+\ket{\phi_{10}})+\ket{\phi_{4}}+\ket{\phi_{13}}\, ,
    \end{split}
\end{align}
where the states $\ket{\phi_{i}}=\ket{\phi_{i}}_{ABA'B'}$ form an orthonormal basis for the purified Hilbert space $\mathcal{H}_{ABA'B'}$ ordered as $(\ket{\downarrow\downarrow\downarrow\downarrow},\ket{\uparrow\downarrow\downarrow\downarrow},\ket{\downarrow\uparrow\downarrow\downarrow},\ket{\uparrow\uparrow\downarrow\downarrow},\dots,\ket{\uparrow\uparrow\uparrow\uparrow})$. From the canonical purification's density matrix $\rho:=\ket{\sqrt{\rho_{AB}}}\bra{\sqrt{\rho_{AB}}}$ we consider a restriction to subsystems $AA'$ given by the reduced density matrix $\rho_{AA'}=\textrm{tr}_{BB'}(\rho)$ which has the perturbative expansion $\rho_{AA'}=\mathrm{tr}_{BB'}(\ket{\psi^{(0)}}\bra{\psi^{(0)}})+\epsilon^{2}_{\sfrac{1}{8}}(2\mathrm{tr}_{BB'}(\ket{\psi^{(1)}}\bra{\psi^{(1)}}))/2$ explicitly given by
\begin{align}
\label{EQ_RHO_AA_SPIN}
   \rho_{AA'} \sim\rho^{(0)}_{AA'}+\tfrac{1}{2}\epsilon^{2}_{\sfrac{1}{8}}\rho^{(2)}_{AA'}=
   \begin{pmatrix}
  \tilde{A}_{1} &  &  & \tilde{F} \\
  & \tilde{B} & \tilde{E} & \\
  & \tilde{E} & \tilde{B} &  \\
 \tilde{F} &  &  & \tilde{A}_{2}
    \end{pmatrix} \, ,
\end{align}
where $\tilde{A}_{1}=\frac{\pi+2}{2\pi}-\frac{2(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{\pi^{2}-4},\tilde{A}_{2}=\tilde{A}_{1}-\frac{2}{\pi},\tilde{B}=\frac{2(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{(\pi^{2}-4)},\tilde{E}=\frac{2\pi C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{\sqrt{\pi^{2}-4}},\tilde{F}=\frac{\sqrt{\pi^{2}-4}}{2\pi}-\frac{2\pi(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{(\pi^{2}-4)^{3/2}}$, where the coefficient $C$ is defined as in~\eqref{eq:Cdefinition}. From here we follow the strategy of the main text and compute the trace of the the square of~\eqref{EQ_RHO_AA_SPIN} according to~\eqref{eq:alpha}.
In this case, we find a value of $\alpha_{\mathrm{tot}}$ computed via~\eqref{eq:alpha} to be
\begin{align}
  \alpha_{\mathrm{tot}}=\frac{4C^{2}(\pi^{2}-2)}{\pi^{2}-4}\approx0.139 \,.\label{eq:alpha-RE-ISING}
\end{align}
As a consequence, the large $d$ leading behaviour of the reflected entropy $S_{R}(\rho_{AB}):=S_{AA'}(\rho)$ exhibits a non trivial logarithmic enhancement of the power law decay according to~\eqref{eq.SAAp} and where the constant contribution can be computed from the eigenvalues of~\eqref{EQ_RHO_AA_SPIN} leading to a reflected entropy $S_{R}$ of the Ising subsystem for $w=\delta$ of
\begin{align}\label{EQ_RE_SPIN_W1}
  S_{R}^\text{spin}(\rho_{AB}) &\sim
  \left( \frac{4C^2\pi^4}{\pi^4 - 16} \log\epsilon_{\sfrac{1}{8}}^{-2} + \mathrm{const} \right)\epsilon_{\sfrac{1}{8}}^2 \nonumber\\
  &\sim \left(0.139\log\sqrt{\frac{d}{w}} + 0.425 \right) \sqrt{\frac{w}{d}}\, .
\end{align}
The constant term is again determined numerically in the main text.

\mysection{Numerical approach and asymmetric purifications}
Our numerical methods are based o
~[29,~70,~71], which outline the construction of an efficient algorithm for local optimization over Gaussian states, based on a gradient descent approach exploiting the natural Lie group parametrization of the state manifolds. Our numerical results are obtained using an adaptation of this algorithm to the non-Gaussian case of interest.

To compute the EoP as given in~\eqref{eq.EoPdef}, we minimise EE $S$ over the manifold $\mathcal{M}$ of purified state density matrices. We first purify our initial mixed density matrix to a $2^N$-dimensional pure $\rho_\id$ via the Schmidt decomposition. Here, $N=\sum_X N_X$ with $N_X$ denoting the physical degrees of freedom in subsystem $X$. We parametrize elements $\rho_U\in\mathcal{M}$ by transformations $U=\id\otimes\widetilde{U}$ with $\widetilde{U}\in\mathrm{U}(2^{N_{A'}+N_{B'}})$, so that $\rho_U=U\rho_\id U^{-1}$. The tensor product signifies that $U$ only acts non-trivially on degrees of freedom in $A'$ and $B'$. We then optimize by performing iterative steps along directions in $\mathcal{M}$ which locally minimize $S_{AA'}$~[29,70]
\begin{align}
    \label{eq.single_step}
    U_{n+1}=U_n\mathrm{e}^{t K_n}\,.
\end{align}
Here, $K_n=\sum_\mu\mathcal{F}^\mu(U_n)\Xi_\mu/||\mathcal{F}||^2$ and $\mathcal{F}^\mu:~\mathcal{M}\to\mathbb{R}$ is the gradient descent vector field
\begin{align}
    \mathcal{F}^\mu(U)=-\frac{\partial}{\partial s}S(U\mathrm{e}^{s\Xi_\mu}\rho_\id\mathrm{e}^{-s\Xi_\mu}U^{-1})|_{s=0}
\end{align}
with $\{\Xi_\mu\}$ as basis of $\mathfrak{u}(2^{N_{A'}+N_{B'}})$. We choose $U_0=\id$ and we pick $0<t<1$ in such a way that the value of $S_{AA'}$ decreases with successive steps.

\begin{table}[t]
    \centering
    \renewcommand{\arraystretch}{1.25}
\begin{tabular}{@{} @{\hskip 7pt}c @{\hskip 7pt}c @{\hskip 7pt}| @{\hskip 7pt} c @{\hskip 7pt}| c @{\hskip 15pt} c c @{\hskip 15pt} c c c @{}
}
\toprule
&& & \multicolumn{6}{c}{$N_{A'}+N_{B'}$} \\
&& & $1+1$ & $1+2$ & $2+1$ & $1+3$ & $2+2$ & $3+1$
\\

\colrule
\parbox[t]{2mm}{\multirow{8}{*}{\rotatebox[origin=c]{90}{$N_A+N_B$}}} & \parbox[t]{2mm}{\multirow{4}{*}{\rotatebox[origin=c]{90}{$1+1$}}} & $d=\delta\;$ & \cellcolor{yellow!35}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$}
\\
&& $d=2\delta$ & \cellcolor{yellow!35}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$}
\\
&& $d=3\delta$ & \cellcolor{yellow!35}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$}
\\
&& $d=2\delta$ & \cellcolor{yellow!35}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$}
\\\hhline{*{1}{~}*{8}{-}}

& \parbox[t]{2mm}{\multirow{4}{*}{\rotatebox[origin=c]{90}{$1+2$}}} & $d=\delta\;$ & \multirow{4}{*}{n.a.} & \cellcolor{yellow!35}{$0.412$} & $0.438$ & \cellcolor{yellow!10}{$0.412$} & \cellcolor{yellow!10}{$0.412$} & $0.440$
\\
&& $d=2\delta$ && \cellcolor{yellow!35}{$0.368$} & $0.412$ & \cellcolor{yellow!10}{$0.368$} & \cellcolor{yellow!10}{$0.368$} & $0.415$
\\
&& $d=3\delta$ && \cellcolor{yellow!35}{$0.345$} & $0.394$ & \cellcolor{yellow!10}{$0.345$} & \cellcolor{yellow!10}{$0.345$} & $0.398$
\\
&& $d=4\delta$ && \cellcolor{yellow!35}{$0.335$} & $0.385$ & \cellcolor{yellow!10}{$0.335$} & \cellcolor{yellow!10}{$0.335$} & $0.389$
\\\botrule

\end{tabular}
    \ccaption{Numerical evidence for optimality of certain minimal purifications}{The table shows the values of the optimization for different choices of the system dimensions and of $d$. The true EoP values (the minimum optimization values) are highlighted in yellow, with the darker shade indicating the lowest-dimensional purification for which the EoP is obtained.}
\label{tab:non-minimal purifications}
\end{table}

The $\{\Xi_\mu\}$ span the tangent space at $U=\id$ and, due to the left-invariance of the Riemannian metric on $\mathcal{M}$, form  orthonormal bases for the tangent spaces at all other points in $\mathcal{M}$, too, where $\Xi_\mu$ is identified with the tangent vector to the curve $\gamma(s)=U\mathrm{e}^{s\Xi_\mu}$ at $\gamma(0)$~[70]
This saves us having to re-evaluate the matrix representation of the metric at each step, as we would have to if we had chosen a coordinate parametrisation of $\mathcal{M}$. While this makes our algorithm more efficient than a naive gradient descent, the numerically accessible range is still highly limited: since $N_{A'}+N_{B'}\geq N_A+N_B$, the dimension of $\mathcal{M}$ is at least $\mathrm{dim}\mathfrak{u}(2^{N_{A'}+N_{B'}})=2^{2N_{A'}+2N_{B'}}-1$ and~\eqref{eq.single_step} requires exponentiation of at least $2^{(N_A+N_B)}\times2^{(N_A+N_B)}$ matrices, with a typical step count of several hundred. This becomes extremely slow on a powerful desktop computer for $N_{A'}+N_{B'}\geq 5$. For the symmetric purifications in the main text this corresponds with $w>2\,\delta$, which explains the regime we were able to explore.

Given this limitation on our numerical capabilities, it is instructive to ask whether an optimization over \textsl{minimal purifications} corresponding with $N_{A'}+N_{B'}=N_A+N_B$ yields the true minimum of EE~\textendash~not least because for large systems this becomes the only numerically viable choice. A natural follow-up question is whether among the choices of minimal purifications, the intuitive choice of $N_{A'}=N_A$ and $N_{B'}=N_B$ suffices to reach the true minimum defined as EoP. More pertinently, we might ask whether it is even possible to reach the true minimum with a minimal purification for which $N_{A'}\neq N_A$ and $N_{B'}\neq N_B$. In~[70],
a combination of numerical and analytical evidence was provided to show that the answer to this question is in affirmative for Gaussian states. While limited by the greater numerical challenge in the non-Gaussian case, we present similar numerical evidence in Table~\ref{tab:non-minimal purifications} to show that the same may be said for our model: the true minimum can only be reached if $N_{A'}\geq N_{A}$ and $N_{B'}\geq N_{B}$, which indicates that the lowest-dimensional purification for which the EoP can be obtained is the minimal purification with $N_{A'}=N_A$ and $N_{B'}=N_B$.

\section*{Long-distance entanglement of purification and reflected entropy\\in conformal field theory -- supplemental material}

\section*{Appendix}

\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\theequation}{S\arabic{equation}}

\mysection{Review of critical Ising model} The Hamiltonian of the transverse Ising model is given by
\begin{equation}\label{eq:TIsingHam}
\hat{H} = - \sum_{k=1}^N \left( 2 J\, \hat{S}^x_k \hat{S}^x_{k+1} + h\, \hat{S}^z_k \right) \ ,
\end{equation}
with spin operators represented by Pauli matrices $\sigma_\alpha$ with $\alpha \in (x,y,z)$ by
\begin{equation}
\hat{S}^\alpha_k \equiv \id^{\otimes (k-1)} \otimes \frac{\sigma_\alpha}{2} \otimes \id^{\otimes (N-k)} \ .
\end{equation}
We also use the identification $\hat{S}^\alpha_{N+1} \equiv \hat{S}^\alpha_1$.
This spin model can be converted to fermions by defining the $2N$ Majorana operators $\operatorname{\gamma}_k$ via
\begin{align}
\operatorname{\gamma}_{2k-1} &= {\sigma_z}^{\otimes (k-1)} \otimes \sigma_x \otimes \id^{\otimes (N-k)} \ , \\
\operatorname{\gamma}_{2k} &= {\sigma_z}^{\otimes (k-1)} \otimes \sigma_y \otimes \id^{\otimes (N-k)} \ .
\end{align}
The Ising Hamiltonian then takes the form
\small
\begin{equation}
\hspace{- 10 pt}\hat{H}_\text{I} = \frac{\mathrm{i}}{2} \left( \operatorname{\gamma}_1 \operatorname{\gamma}_{2N} P + J \sum_{k=1}^{N-1} \operatorname{\gamma}_{2k} \operatorname{\gamma}_{2k+1} + h \sum_{k=1}^{N} \operatorname{\gamma}_{2k-1} \operatorname{\gamma}_{2k} \right) .
\end{equation}
\normalsize
Here $P$ is the total parity operator ${\sigma_z}^{\otimes N}=\prod_k (-\mathrm{i} \operatorname{\gamma}_{2k-1}\operatorname{\gamma}_{2k})$.
At the critical point $J=h$, the Hamiltonian thus simplifies to
\begin{equation}
\hat{H} = \frac{\mathrm{i}}{2} \left( \operatorname{\gamma}_1 \operatorname{\gamma}_{2N} P + \sum_{k=1}^{2N-1} \operatorname{\gamma}_k \operatorname{\gamma}_{k+1} \right)\,,
\end{equation}
which leads for $N\to\infty$ to the lattice model of the $c=\tfrac{1}{2}$ CFT. The critical Ising Hamiltonian as displayed in the main text~\eqref{eq.HIsing}
corresponds to $J = h = 1$.

\mysection{Covariance matrix} For the critical ground state vector $\ket{0}$ which has a positive total parity, all correlations are encoded in the Majorana covariance matrix
\begin{equation}
\Omega_{j,k} = \frac{\mathrm{i}}{2} \bra{0} [\operatorname{\gamma}_j, \operatorname{\gamma}_k] \ket{0} \ ,
\end{equation}
which in the infinite system size limit takes the form
\begin{equation}
\Omega_{j,k} =
\begin{cases}
0 &  k = j  \\
\frac{(-1)^{k-j} - 1}{\pi (k-j)} &  k \neq j
\end{cases} \ .\label{eq:covariance-inf}
\end{equation}
The entropy of a Gaussian mixed state $\rho$ with covariance matrix $\Omega$ is given by
\begin{align}
    S(\rho)=-\sum_{\pm, i}\frac{1\pm\lambda_i}{2}\log\frac{1\pm\lambda_i}{2}\,,\label{eq:Gaussain-SA}
\end{align}
where $\pm\mathrm{i} \lambda_k$ are the eigenvalues of $\Omega$. We will consider mixed states $\rho_{AB}$ or $\rho_{AA'}$, whose mixed state covariance matrices $\Omega_{AB}$ or $\Omega_{AA'}$ result from restricting~\eqref{eq:covariance-inf} to the respective blocks.

\mysection{Fermionic subsystem}
\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{Fig3.pdf}
    \caption{Subsystem setup of our analytical limits for fermions (top) with an inherent ordering and spins (bottom) without one. In both systems, we consider the subsystem $AB$ consisting of two single sites $A$ and $B$ separated by $d/\delta$ sites.
    }
    \label{fig:1+1subsystem}
\end{figure}
We first compute reduced density matrices from the perspective of fermions, {\it i.e.},\  imposing an ordering between modes following from the anti-commuting variables $\gamma_i$ (see Fig.~\ref{fig:1+1subsystem}). A subsystem consisting of $1+1$ sites ($w=\delta$) separated by $d/w=d/\delta$ sites is then fully characterized by the restriction of the covariance matrix in~\eqref{eq:covariance-inf} and explicitly given by
\begin{align}
    \Omega^{\mathrm{fer}}_{AB} =
    \begin{pmatrix}
      & -\frac{2}{\pi } &  & \frac{-2}{(2 d/w+3) \pi } \\
     \frac{2}{\pi } &  & \frac{-2}{(2 d/w+1) \pi } &  \\
      & \frac{2}{(2 d/w+1) \pi } &  & -\frac{2}{\pi } \\
     \frac{2}{(2 d/w+3) \pi } &  & \frac{2}{\pi } &
    \end{pmatrix}\label{eq:omega-fermions}
\end{align}
which corresponds to a lowest-dimension (Majorana) operator with scaling dimension $\Delta=\sfrac{1}{2}$. The associated fermionic density operator is then
\begin{align}\label{eq:rho-fermions}
    \hspace{-2mm}\rho^{\mathrm{fer}}_{AB}\sim\begin{pmatrix}
     D & & &\frac{1}{2\pi}\,\epsilon_{\sfrac{1}{2}}\\
     & E & & \\
     & & E & \\
     \frac{1}{2\pi}\,\epsilon_{\sfrac{1}{2}} & & & F
    \end{pmatrix}
\end{align}
with respect to the basis $(\ket{\downarrow\downarrow}$, $\ket{\uparrow\downarrow},\ket{\downarrow\uparrow},\ket{\uparrow\uparrow})$ and using $D=\frac{1}{4}+\frac{1}{\pi}+\frac{1}{\pi^2}$, $E=\frac{1}{4}-\frac{1}{\pi^2}$, and $F=\frac{1}{4}-\frac{1}{\pi}+\frac{1}{\pi^2}$.
As in the main text, $\epsilon_\Delta \equiv (w/d)^{2\Delta}$ which here becomes $\epsilon_{\sfrac{1}{2}}=w/d$.
If we further restrict to a single site, we find the covariance matrix and density operator
\begin{align}
    \Omega^{\mathrm{fer}}_{A} &=
    \begin{pmatrix}
      & -\frac{2}{\pi }  \\
     \frac{2}{\pi } &
    \end{pmatrix} & \rho_A^{\mathrm{fer}} &=
    \begin{pmatrix}
     \frac{1}{2}-\frac{1}{\pi} & \\
      & \frac{1}{2}+\frac{1}{\pi}
    \end{pmatrix}\,,
\end{align}
where $\rho_A^{\mathrm{fer}}$ is written with respect to the basis $(\ket{\downarrow},\ket{\uparrow})$.

\mysection{Spin subsystem} We can perform a similar calculation in the original Ising spin system whose reduced density matrices can be constructed from the fermionic covariance matrix \cite{Coser:2015mta}. We need not repeat the single interval case, as entanglement entropies of connected regions are equivalent under a Jordan-Wigner transformation. However, we still need the reduced density matrix of a system of $1+1$ sites in the large $d$ limit, which we find to be
\begin{equation}
\label{eq:rhoAB}
    \rho^{\mathrm{spin}}_{AB} \sim
    \begin{pmatrix}
      D &  &  & C\epsilon_{\sfrac{1}{8}} \\
  & E & C\epsilon_{\sfrac{1}{8}} & \\
  & C\epsilon_{\sfrac{1}{8}} & E &  \\
 C\epsilon_{\sfrac{1}{8}} &  &  & F
    \end{pmatrix}
\end{equation}
with $w/\delta=1$ for the setup of $1+1$ sites. As we are considering the spin Ising CFT, the lowest-dimension primary is the ``order field'' $\sigma$ with scaling dimension $\Delta=\sfrac{1}{8}$.
The constant $C$ corresponds to the expectation value of an operator nonlocal in fermions, and can be computed from
\begin{equation}
    C = \lim_{n \to \infty} \left( \frac{2}{\pi} \right)^n \frac{n^{1/4}}{4} \det{M^n} \ ,
\end{equation}
where $M^n$ is defined as the $n \times n$ matrix
\begin{equation}
    M^n_{j,k} =
    \begin{cases}
    \frac{(-1)^{k-j}}{2(k-j)+1} & j \leq k \\
    \frac{(-1)^{j-k+1}}{2(j-k)-1} & j > k
    \end{cases} \ .
\end{equation}
Using this construction, one finds~\cite{McCoy:2010}
\begin{equation}
\label{eq:Cdefinition}
    C =\frac{e^{3\zeta'(-1)}}{2^{23/12}}\approx 0.1612506\,.
\end{equation}

\mysection{MI for fermions}
We now begin computing entanglement measures for the small subsystems whose reduced density matrices we just obtained explicit expressions for.
The continuum limit corresponding to the Ising CFT is obtained by keeping $d/w$ (or, equivalently, $\epsilon_\Delta$) fixed and taking $\delta/w$ to $0$. We will see that taking only a few lattice sites is sufficient to describe the qualitative and approximate quantitative behavior of the continuum limit.
To demonstrate this, we now show that the large distance asymptotics of the MI ~\eqref{eq.IAB} for the case of $1+1$ sites yield results close to the continuum formula~\eqref{eq.MIcftlarged}.
In order to compute the MI, we need to determine the von Neumann entropy of a single site $S_A=S_B$ and of both sites $S_{AB}$. These entropies can be computed from the eigenvalues of the covariance matrix $\Omega$ associated to the respective Gaussian state $\rho$. As an antisymmetric matrix, $\Omega^{\mathrm{fer}}_{AB}$  has pairs of purely imaginary eigenvalues $\pm \mathrm{i}\, \lambda_k$, from which applying~\eqref{eq:Gaussain-SA} leads to
\begin{align}
S_A&=- \tfrac{\pi+2}{2\pi}\log\tfrac{\pi+2}{2\pi} - \tfrac{\pi-2}{2\pi}\log\tfrac{\pi-2}{2\pi} \approx 0.474\\
    S_{AB} &=\sum_{k=1}^n \left( - \tfrac{1+\lambda_k}{2}\log\tfrac{1+\lambda_k}{2} - \tfrac{1-\lambda_k}{2}\log\tfrac{1-\lambda_k}{2} \right)\,,\label{eq:Gaussian-SAB}
\end{align}
where the eigenvalues of $\Omega_{AB}^{\mathrm{fer}}$ are to leading order
\begin{equation}
    \lambda_{1,2} = \frac{1}{\pi} \left( 2 \pm \tfrac{3}{4}\epsilon_{\sfrac{1}{2}}^2  + \dots \right) \ .
\end{equation}
We can similarly expand $S_{1+1} \equiv S_{AB}$ at large $d$, which results in a MI for $w=\delta$ of
\begin{align}\label{eq:MIFerm11}
    I^{\mathrm{fer}}(A:B)&\sim \tfrac{\log\frac{\pi+2}{\pi-2}}{4\pi}\,\epsilon_{\sfrac{1}{2}}^2 \nonumber\\
    &= 0.120 \,\left(\tfrac{w}{d}\right)^2 \ .
\end{align}
This reproduces the correct continuum power law of fermionic MI, but yields a coefficient lower than the continuum value \eqref{eq.MIcftlarged} which also matches the large-distance expansion of earlier results for Dirac fermions~\cite{Casini:2005rm}
\begin{align}\label{eq:MIFerm}
    I(A:B)=\frac{c}{3}\log\frac{(d+w)^2}{d\,(2w+d)} &\sim \frac{1}{6} \left(\frac{w}{d}\right)^2 \ ,
\end{align}
for two blocks of general width $w$.

\mysection{MI for spins}
We compute EE for the spin system directly from the eigenvalue spectrum of the reduced density matrix \eqref{eq:rhoAB}. Its four eigenvalues $\mu_j$ are
\begin{subequations}
\begin{align}\label{eq:EigenSpinRhoAB}
    \mu_{1,2} &= \frac{1}{4} - \frac{1}{\pi} \pm C \,\epsilon_{\sfrac{1}{8}} \ , \\
    \mu_{3,4} &= \frac{1}{4} + \frac{1}{\pi} \pm \sqrt{\frac{1}{\pi^2} + C^2 \, \epsilon_{\sfrac{1}{8}}^2}\ ,
\end{align}
\end{subequations}
where $C$ is given by \eqref{eq:Cdefinition} and from which we can directly compute the EE for $AB$ via
\begin{equation}
\label{eq.Smu}
S = - \sum_j \mu_j \log\mu_j.
\end{equation}
Note that in the following we will denote eigenvalues of any density matrix by $\mu_{j}$.

This analysis leads to the Ising model prediction for spin MI at $w=\delta$ and large separations $d$ of the form
\begin{align}
\label{eq.MIanalytic}
 \hspace{-2mm}   I^\text{spin}(A:B) &\sim
 C^2 \left( \frac{4\pi^2}{\pi^2-4} + \frac{\pi}{2}\log\frac{4+4\pi+\pi^2}{4-4\pi+\pi^2} \right)  \epsilon_{\sfrac{1}{8}}^2 \nonumber\\
 &\sim 0.298\, \sqrt{\frac{w}{d}} \ ,
\end{align}
\normalsize
Comparing this formula with the CFT analytics~\eqref{eq.MIcftlarged} for $\Delta = \sfrac{1}{8}$, we see an exact match in the power-law behavior. Furthermore, the prefactor in~\eqref{eq.MIanalytic} is only $3.6\%$ off from the continuum value $\approx0.309$ predicted by~\eqref{eq.MIcftlarged}.

\mysection{EoP for fermions} Analogous to the MI calculation for free fermions, we now calculate the EoP in the fermionic subsystem of two sites separated by $d/\delta$ sites, expressing all calculations in terms of covariance matrices. We purify $\Omega_{AB}$ in the limit $d/\delta\to\infty$ as
\begin{align}
    \Omega^{(0)}=\left(\begin{array}{cccc|cccc}
      & -G & & &  & L& & \\
    G &    & & & L &  & &\\
      &  & & -G & & & & L\\
     &    & G & & &  & L&\\
     \hline
       & -L & & &  & -G& & \\
    -L &    & & & G & & &\\
      &  & & -L & & &  & -G\\
     &    & -L & & &  & G&\\
    \end{array}\right)
\end{align}
associated to systems $(A,B,A',B')$ with $G=\frac{2}{\pi}$ and $L=\sqrt{1-G^2}$, whose EE $S_{AA'}$ is zero and we thus have $\lim_{d\to\infty}E_P=0$, {\it i.e.},\  the EoP vanishes for large $d/\delta$, as expected.

In order to find the asymptotic behavior of $E_P$, we need to study the variation of the symplectic eigenvalues $\pm\mathrm{i}\lambda_i$ of $\Omega_{AA'}$ when perturbing $\Omega$ according to
\begin{align}
    \Omega\sim \Omega^{(0)}+\epsilon_{\sfrac{1}{2}}\,\Omega^{(1)}+\tfrac{1}{2}\epsilon_{\sfrac{1}{2}}^2\,\Omega^{(2)}\quad\text{as}\quad \epsilon_{\sfrac{1}{2}}\to 0\,.
\end{align}
The requirement of $\Omega$ representing a purification implies $\Omega^2=-\id$, which induces the constraints
\begin{align}
\begin{split}
    \Omega^{(0)}\Omega^{(1)}+\Omega^{(1)}\Omega^{(0)}&=0\,,\\
    2(\Omega^{(1)})^2+\Omega^{(0)}\Omega^{(2)}+\Omega^{(2)}\Omega^{(0)}&=0\,.
\end{split}\label{eq:constraints-J}
\end{align}
We further require that the restrictions $\Omega^{(1)}_{AB}$ and $\Omega^{(2)}_{AB}$ matches the ones of\ \eqref{eq:omega-fermions} expanded in $\epsilon_{\sfrac{1}{2}}$, {\it i.e.},\
\begin{align}
    \Omega^{(1)}_{AB}\!=\!\begin{pmatrix}
    &&& -\frac{1}{\pi}\\
    &&-\frac{1}{\pi}& \\
    &\frac{1}{\pi}&& \\
    \frac{1}{\pi}&&&
    \end{pmatrix},\,\,
    \Omega^{(2)}_{AB}\!=\!\begin{pmatrix}
    &&& \frac{3}{\pi}\\
    &&\frac{1}{\pi}& \\
    &-\frac{1}{\pi}&& \\
    -\frac{3}{\pi}&&&
    \end{pmatrix}\,,\label{eq:constraints2-J}
\end{align}
The equations\ \eqref{eq:constraints-J} and\ \eqref{eq:constraints2-J} can be solved iteratively up to some free variables. We first solve $\Omega^{(1)}$ in terms of $\Omega^{(0)}$ and then $\Omega^{(2)}$ in terms of $\Omega^{(0)}$ and $\Omega^{(1)}$.

In order to find asymptotics of the symplectic eigenvalues $\lambda_i$, we can use the fact that $\operatorname{Tr}(\Omega_{AA'}^2)=-2(\lambda_1^2+\lambda_2^2)$ and $\operatorname{Tr}(\Omega_{AA'}^4)=2(\lambda_1^4+\lambda_2^4)$ to solve for the asymptotics of $\lambda_i$ to be given by
\begin{align}
    \lambda_1=\lambda_2\sim 1-\alpha_{\mathrm{tot}}\,\epsilon_{\sfrac{1}{2}}^2\quad \text{as}\quad \epsilon_{\sfrac{1}{2}}\to 0\,,\label{eq:symp-eig}
\end{align}
where $\alpha_{\mathrm{tot}}$ will depend on some of the free parameters contained in $\Omega^{(1)}$ and $\Omega^{(2)}$. With this trick, one finds
\begin{align}
\begin{split}
    \alpha_{\mathrm{tot}}&=\frac{x_{14}a_{23}-x_{13}x_{24}+\pi^{-2}}{2}+\frac{G(x_{14}-x_{23})\pi^{-1}}{2L}\\
    &\quad +\frac{(x_{14}-x_{23})^2+(x_{13}+x_{24})^2}{4L^2}\,,\label{eq:alpha2}
\end{split}
\end{align}
where the variables $x_{ij}$ represent unconstrained entries in the block $\Omega^{(1)}_{AB,A'B'}$. In order to find the asymptotics of EoP, we need to minimize $\alpha_{\mathrm{tot}}$ over these parameters to find the smallest possible EE $S_{AA'}$. Due to the fact that~\eqref{eq:alpha2} is quadratic in $x_{ij}$, we can calculate this valua analytically as
\begin{align}\label{eq:alphatotEoPFerm}
    \alpha_{\mathrm{tot}}&=\frac{1}{8+2\pi^2}\approx 0.03605\,.
\end{align}
Expanding $S_{AA'}\sim \sum_i(\log{2}-\tfrac{\lambda_i}{2})$ through $\lambda_i$ up to second order in $\epsilon_{\sfrac{1}{2}}$ based on~\eqref{eq:symp-eig} allows us to also find the offset analytically, namely we have
\begin{align}
    S_{AA'}=\epsilon_{\sfrac{1}{2}}^2\left(\alpha_{\mathrm{tot}}\log(\epsilon_{\sfrac{1}{2}}^{-2})+\alpha_{\mathrm{tot}} \log{\tfrac{2e}{\alpha_{\mathrm{tot}}}}\right)\,.
\end{align}
Combining this with the result from~\eqref{eq:alphatotEoPFerm} gives
\begin{align}
\label{eq:eopFerm11}
    E_P^\text{fer}&\sim \left(\frac{1}{8+2\pi^2} \log(\epsilon_{\sfrac{1}{2}}^{-2})+ \frac{\log 2e (8+2\pi^2)}{8+2\pi^2} \right)\, \epsilon_{\sfrac{1}{2}}^2 \nonumber\\
    &\sim \left(0.0361\log\left(\tfrac{d}{w}\right)^2 +0.181\right)\, \left(\tfrac{w}{d}\right)^2  \ ,
\end{align}
which agrees with the form~\eqref{eq.SAAp} in the main text. Note that the simplicity of Gaussian states allowed us to even find the analytical form of the constant offset. The accuracy of this analytical prediction was tested numerically, for which we presented the results in Fig.~\ref{fig:MI_EoP} in the main text.

\mysection{EoP for spins} In the limit of an infinite distance between the two single site subsystems, we purify~\eqref{eq:rhoAB} by the state $\ket{\psi^{(0)}}$ with Schmidt decomposition
\begin{align}
\hspace{-7 pt}    \ket{\psi^{(0)}}=\sqrt{D}\ket{\downarrow\downarrow\downarrow\downarrow}+\sqrt{E}(\ket{\uparrow\downarrow\uparrow\downarrow}+\ket{\downarrow\uparrow\downarrow\uparrow})+\sqrt{F}\ket{\uparrow\uparrow\uparrow\uparrow},\label{eq:schmidt}
\end{align}
\normalsize
where the convention for factors ordering in the purification is $ABA'B'$. Note that in this analysis we assume that a minimal purification from two to four spin degrees of freedom suffices and we will subsequently provide supporting numerical evidence and an additional discussion.

Moving on, we \emph{supplement} this purification with finite distance corrections up to second order in $\epsilon_{\sfrac{1}{8}}$ as
\begin{align}
    \ket{\psi}\sim \ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{8}} \ket{\psi^{(1)}}+\tfrac{1}{2}\epsilon_{\sfrac{1}{8}}^2 \ket{\psi^{(2)}}\,.
\end{align}
\normalsize
We will optimize over $\ket{\psi^{(1)}}$ and $\ket{\psi^{(2)}}$ subject to the normalization constraint $\braket{\psi|\psi}=1$ order by order in $\epsilon_{\sfrac{1}{8}}$. We further require $\rho^{(1)}=\ket{\psi^{(0)}}\bra{\psi^{(1)}}+\ket{\psi^{(1)}}\bra{\psi^{(0)}}$ and $\rho^{(2)}=\ket{\psi^{(0)}}\bra{\psi^{(2)}}+\ket{\psi^{(2)}}\bra{\psi^{(0)}}+2\ket{\psi^{(1)}}\bra{\psi^{(1)}}$ to satisfy
\begin{align}
    \rho^{(1)}_{AB}=\begin{pmatrix}
     &&&C\\
     &&C&\\
     &C&&\\
     C&&&
    \end{pmatrix}\quad\text{and}\quad \rho^{(2)}_{AB}=0\,,\label{eq:nG-constaints2}
\end{align}
which follows from\ \eqref{eq:rhoAB}. We expand the first order perturbation as
\begin{align}
    \ket{\psi^{(1)}}=C\sum^{16}_{i=1} z_i\ket{\phi_i}
\end{align}
where $z_i=x_i+\mathrm{i} y_i$ and $\ket{\phi_i}$ is the basis of $\mathcal{H}_{ABA'B'}$ ordered as $(\ket{\downarrow\downarrow\downarrow\downarrow},\ket{\uparrow\downarrow\downarrow\downarrow},\ket{\downarrow\uparrow\downarrow\downarrow},\ket{\uparrow\uparrow\downarrow\downarrow},\dots,\ket{\uparrow\uparrow\uparrow\uparrow})$. We then need to implement the condition~\eqref{eq:affine-lin} in the main text based on~\eqref{eq:nG-constaints2} together with the normalization constraint~\eqref{eq:constraint1} in the main text. We solve these affine linear constraints by the replacements $x_1=x_6=x_{11}=x_{16}=0$, $z_{5}=-\sqrt{\frac{E}{D}}z_2^*$, $z_{9}=-\sqrt{\frac{E}{D}}z_3^*$, $z_{13}=\frac{1}{\sqrt{D}}-\sqrt{\frac{F}{D}}z_4^*$, $z_{15}=\sqrt{\frac{F}{E}}z_{12}^*$, $z_{10}=\frac{1}{\sqrt{E}}-z_{7}^*$, $z_{14}=-\sqrt{\frac{F}{E}}z_8$.
We can then compute $\alpha_{\mathrm{tot}}$ according to~\eqref{eq:alpha} in the main text as quadratic polynomial in terms of the remaining free variables $z_i$ which leads to the rather involved expression
\begin{widetext}
\begin{align}
\begin{split}

    \hspace{-3mm}\frac{\alpha_{\mathrm{tot}}}{C^2}&
    =\tfrac{(\pi -2)^2 y_1^2}{4 \pi ^2}-\tfrac{2 \sqrt{\pi ^2-4} x_3 x_8}{2+\pi }-\tfrac{2 \sqrt{\pi ^2-4} x_2 x_{12}}{2+\pi }-\tfrac{4 (\pi -2) \pi  x_4}{(2+\pi )^2}-\tfrac{4 \pi  x_7}{\sqrt{\pi ^2-4}}+\tfrac{(\pi -2) \left(-\sqrt{\pi ^2-4} y_6-\sqrt{\pi ^2-4} y_{11}+(2+\pi ) y_{16}\right) y_1}{2 \pi ^2}\\
    &\quad+\tfrac{\left(\pi ^2-4\right) y_6^2}{4 \pi ^2}+\tfrac{\left(\pi ^2-4\right) y_{11}^2}{4 \pi ^2}+\left(\tfrac{1}{4}+\tfrac{1}{\pi ^2}+\tfrac{1}{\pi }\right) y_{16}^2+\tfrac{(2+\pi ) y_6 \left((\pi -2) y_{11}-\sqrt{\pi ^2-4} y_{16}\right)}{2 \pi ^2}-\tfrac{(2+\pi ) \sqrt{\pi ^2-4} y_{11} y_{16}}{2 \pi ^2}\\
    &\quad-\tfrac{2 \sqrt{\pi ^2-4} y_3 y_8}{2+\pi }-\tfrac{2 \sqrt{\pi ^2-4} y_2 y_{12}}{2+\pi }+\tfrac{(\pi -2) |z_2|^2}{2+\pi }+\tfrac{(\pi -2) |z_3|^2}{2+\pi }+\tfrac{2 \left(4+\pi ^2\right) |z_4|^2}{(2+\pi )^2}+2 |z_7|^2+|z_8|^2+|z_{12}|^2+\tfrac{8 \pi ^3}{(\pi -2) (2+\pi )^2}\,.
\end{split}
\end{align}
\end{widetext}
In order to find the EoP, we need to minimize over the $z_i$ to find the smallest possible value $\alpha_{\mathrm{tot}}$, which can be done \emph{analytically} and leads to
\begin{align}
\label{eq.alphatot}
    \alpha_{\mathrm{tot}}=\frac{4\pi^4C^2}{\pi^4-16}\approx 0.12445\,.
\end{align}
The non-vanishing $\alpha_{\mathrm{tot}}$ shows that the resulting EoP obtained from~\eqref{eq.SAAp} again has the form
\begin{align}
\label{EQ_EOP_SPIN_W1}
    E_P^\text{spin} &\sim \left(  \frac{4\pi^4C^2}{\pi^4-16} \log(\epsilon_{\sfrac{1}{8}}^{-2}) + \mathrm{const} \right) \epsilon_{\sfrac{1}{8}}^{2} \nonumber\\
    &\sim \left( 0.124 \log\sqrt{\frac{d}{w}} + 0.440 \right) \sqrt{\frac{w}{d}} \ ,
\end{align}
which, as in the fermion case, exhibits a leading-order long-distance behavior enhanced with respect to that of MI~\eqref{eq.MIanalytic} by a logarithm of the distance.

When it comes to the \emph{subleading} long-distance behavior encapsulated by~$\Big( \sum_{j>0}\alpha_j(1-\log{\alpha_j}) \Big)$, we would need to extract the individual $\alpha_j$ and optimize over the remaining parameters. While it is plausible this can be also done analytically, we determined the value quoted above numerically, as discussed in the main text. Note that for the free fermion case with $w=\delta$ considered above, we determined this term analytically in terms of~$\alpha_{\mathrm{tot}}$.

\mysection{RE for fermions} In the Gaussian case of free fermions, our starting point is the following perturbative expansion of the reduced density matrix $\rho_{AB}$ of a system of $1+1$ fermions in the large $d$ separation, $\rho_{AB}\sim\rho^{(0)}_{A}\otimes\rho^{(0)}_{B}+\epsilon_{1/2}\rho^{(1)}_{AB}$ given by~\eqref{eq:rho-fermions}. We similarly construct the canonical purification of~\eqref{eq:rho-fermions} via $\ket{\sqrt{\rho_{AB}}}=\sum_{i}\sqrt{e_{i}}\ket{e_{i}}\otimes\ket{e_{i}}=\ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{2}}\, \ket{\psi^{(1)}}$ where $\rho_{AB}\ket{e_{i}}=e_{i}\ket{e_{i}}$. Note that in contrast with fermionic MI and EoP, we do not need to phrase our computation of RE in terms of the covariance matrix formalism since we can construct the canonical purification $\ket{\sqrt{\rho_{AB}}}$exactly for the given form of the initial reduced density matrix~$\rho_{AB}$.

In this case, the first-order perturbation $\ket{\psi^{(1)}}$ is simply given by
\begin{align}\label{EQ_PSI1_FERMION}
    \ket{\psi^{(1)}}=\tfrac{1}{2\pi}(\ket{\phi_{4}}+\ket{\phi_{13}}) \, ,
\end{align}
with the same ordering of the basis $\ket{\phi_{i}}$ as in the previous case. From the canonical purification's density matrix $\rho:=\ket{\sqrt{\rho_{AB}}}\bra{\sqrt{\rho_{AB}}}$ we consider a restriction to subsystems $AA'$ given by the reduced density matrix $\rho_{AA'}=\textrm{tr}_{BB'}(\rho)$ which has the perturbative expansion $\rho_{AA'}=\rho^{(0)}_{A}\otimes\rho^{(0)}_{A'}+\epsilon^{2}_{\sfrac{1}{2}}\rho^{(2)}_{AA'}$ explicitly given by
\begin{align}
\label{EQ_RHO_AA_FERMION}
   \rho_{AA'} \sim \rho^{(0)}_{AA'}+\tfrac{1}{2}\epsilon^{2}_{\sfrac{1}{2}}\rho^{(2)}_{AA'}=
   \begin{pmatrix}
  \tilde{G}_{1} &  &  & \tilde{H} \\
  & \tilde{J} &  & \\
  &  &   \tilde{J}&  \\
 \tilde{H} &  &  & \tilde{G}_{2}
    \end{pmatrix} ,
\end{align}
where $\tilde{G}_{1}=\frac{\pi+2}{2\pi}-\frac{\epsilon^{2}_{\sfrac{1}{2}}}{4\pi^{2}}$, $\tilde{G}_{2}=\tilde{G}_{1}-\frac{2}{\pi}$, $\tilde{H}=\frac{\sqrt{\pi^{2}-4}}{2\pi}-\frac{\sqrt{\pi^{2}-4}\epsilon^{2}_{\sfrac{1}{2}}}{4\pi(\pi^{2}-4)}$, and $\tilde{J}=\frac{\epsilon^{2}_{\sfrac{1}{2}}}{4\pi^{2}}$. We once again compute the trace of the square of~\eqref{EQ_RHO_AA_FERMION} according to~\eqref{eq:alpha} in the main text from which we obtain
\begin{align}\label{EQ_APLHA_TOT_RE_FERMION}
    \alpha_{\mathrm{tot}}=\frac{1}{2\pi^{2}}\approx 0.051 \, .
\end{align}
which also shows that the reflected entropy $S_{R}(\rho_{AB})=S_{AA'}(\rho)$ of the fermionic subsystem also exhibits a logarithmic enhancement of the power law decay for $w=\delta$ given by
\begin{align}\label{EQ_RE_FERMION_W1}
    S_{R}^\text{fer}(\rho_{AB}) &\sim
    \left( \frac{1}{2\pi^{2}} \log\epsilon_{\sfrac{1}{2}}^{-2} + \frac{1 + \log(4\pi^2)}{2\pi^2} \right) \epsilon_{\sfrac{1}{2}}^2 \nonumber  \\
    & \sim\left(0.051\log\left(\frac{d}{w}\right)^2+0.237\right) \left(\frac{w}{d}\right)^{2}\, ,
\end{align}
where we also computed the constant term in~\eqref{EQ_RE_FERMION_W1} from the eigenvalues of~\eqref{EQ_RHO_AA_FERMION} according to~\eqref{eq:alpha} in the main text.

\mysection{RE for spins} For the Ising spin case, we now describe the detailed computation of the reflected entropy RE for $w=\delta$ in the large $d$ limit just as for fermions. The reduced density matrix for a spin system of $1+1$ sites in the large $d$ limit can again be computed according to~\eqref{eq.rhogeneral}, {\it i.e.},\  $\rho_{AB} \sim \rho^{(0)}_{A} \otimes \rho^{(0)}_{B}+\epsilon_{\sfrac{1}{8}}\, \rho^{(1)}_{AB}+\ldots$, yielding~\eqref{eq:rhoAB}.

We now construct the canonical purification of~\eqref{eq:rhoAB} via $\ket{\sqrt{\rho_{AB}}}=\sum_{i}\sqrt{e_{i}}\ket{e_{i}}\otimes\ket{e_{i}}=\ket{\psi^{(0)}}+\epsilon_{\sfrac{1}{8}}\ket{\psi^{(1)}}$ for $\rho_{AB}\ket{e_{i}}=e_{i}\ket{e_{i}}$ and where the eigenvalues $e_{i}$ are defined in~\eqref{eq:EigenSpinRhoAB}. In this case, the first order perturbation $\ket{\psi^{(1)}}$ is given by
\begin{align} \label{EQ_PSI1_ISING}
\begin{split}
    \ket{\psi^{(1)}}&= \tfrac{\pi}{\sqrt{\pi^{2}-4}}(\ket{\phi_{7}}+\ket{\phi_{10}})+\ket{\phi_{4}}+\ket{\phi_{13}}\, ,
    \end{split}
\end{align}
where the states $\ket{\phi_{i}}=\ket{\phi_{i}}_{ABA'B'}$ form an orthonormal basis for the purified Hilbert space $\mathcal{H}_{ABA'B'}$ ordered as $(\ket{\downarrow\downarrow\downarrow\downarrow},\ket{\uparrow\downarrow\downarrow\downarrow},\ket{\downarrow\uparrow\downarrow\downarrow},\ket{\uparrow\uparrow\downarrow\downarrow},\dots,\ket{\uparrow\uparrow\uparrow\uparrow})$. From the canonical purification's density matrix $\rho:=\ket{\sqrt{\rho_{AB}}}\bra{\sqrt{\rho_{AB}}}$ we consider a restriction to subsystems $AA'$ given by the reduced density matrix $\rho_{AA'}=\textrm{tr}_{BB'}(\rho)$ which has the perturbative expansion $\rho_{AA'}=\mathrm{tr}_{BB'}(\ket{\psi^{(0)}}\bra{\psi^{(0)}})+\epsilon^{2}_{\sfrac{1}{8}}(2\mathrm{tr}_{BB'}(\ket{\psi^{(1)}}\bra{\psi^{(1)}}))/2$ explicitly given by
\begin{align}
\label{EQ_RHO_AA_SPIN}
   \rho_{AA'} \sim\rho^{(0)}_{AA'}+\tfrac{1}{2}\epsilon^{2}_{\sfrac{1}{8}}\rho^{(2)}_{AA'}=
   \begin{pmatrix}
  \tilde{A}_{1} &  &  & \tilde{F} \\
  & \tilde{B} & \tilde{E} & \\
  & \tilde{E} & \tilde{B} &  \\
 \tilde{F} &  &  & \tilde{A}_{2}
    \end{pmatrix} \, ,
\end{align}
where $\tilde{A}_{1}=\frac{\pi+2}{2\pi}-\frac{2(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{\pi^{2}-4},\tilde{A}_{2}=\tilde{A}_{1}-\frac{2}{\pi},\tilde{B}=\frac{2(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{(\pi^{2}-4)},\tilde{E}=\frac{2\pi C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{\sqrt{\pi^{2}-4}},\tilde{F}=\frac{\sqrt{\pi^{2}-4}}{2\pi}-\frac{2\pi(\pi^{2}-2)C^{2}\epsilon^{2}_{\sfrac{1}{8}}}{(\pi^{2}-4)^{3/2}}$, where the coefficient $C$ is defined as in~\eqref{eq:Cdefinition}. From here we follow the strategy of the main text and compute the trace of the the square of~\eqref{EQ_RHO_AA_SPIN} according to~\eqref{eq:alpha}.
In this case, we find a value of $\alpha_{\mathrm{tot}}$ computed via~\eqref{eq:alpha} to be
\begin{align}
  \alpha_{\mathrm{tot}}=\frac{4C^{2}(\pi^{2}-2)}{\pi^{2}-4}\approx0.139 \,.\label{eq:alpha-RE-ISING}
\end{align}
As a consequence, the large $d$ leading behaviour of the reflected entropy $S_{R}(\rho_{AB}):=S_{AA'}(\rho)$ exhibits a non trivial logarithmic enhancement of the power law decay according to~\eqref{eq.SAAp} and where the constant contribution can be computed from the eigenvalues of~\eqref{EQ_RHO_AA_SPIN} leading to a reflected entropy $S_{R}$ of the Ising subsystem for $w=\delta$ of
\begin{align}\label{EQ_RE_SPIN_W1}
  S_{R}^\text{spin}(\rho_{AB}) &\sim
  \left( \frac{4C^2\pi^4}{\pi^4 - 16} \log\epsilon_{\sfrac{1}{8}}^{-2} + \mathrm{const} \right)\epsilon_{\sfrac{1}{8}}^2 \nonumber\\
  &\sim \left(0.139\log\sqrt{\frac{d}{w}} + 0.425 \right) \sqrt{\frac{w}{d}}\, .
\end{align}
The constant term is again determined numerically in the main text.

\mysection{Numerical approach and asymmetric purifications}
Our numerical methods are based o
~[29,~70,~71], which outline the construction of an efficient algorithm for local optimization over Gaussian states, based on a gradient descent approach exploiting the natural Lie group parametrization of the state manifolds. Our numerical results are obtained using an adaptation of this algorithm to the non-Gaussian case of interest.

To compute the EoP as given in~\eqref{eq.EoPdef}, we minimise EE $S$ over the manifold $\mathcal{M}$ of purified state density matrices. We first purify our initial mixed density matrix to a $2^N$-dimensional pure $\rho_\id$ via the Schmidt decomposition. Here, $N=\sum_X N_X$ with $N_X$ denoting the physical degrees of freedom in subsystem $X$. We parametrize elements $\rho_U\in\mathcal{M}$ by transformations $U=\id\otimes\widetilde{U}$ with $\widetilde{U}\in\mathrm{U}(2^{N_{A'}+N_{B'}})$, so that $\rho_U=U\rho_\id U^{-1}$. The tensor product signifies that $U$ only acts non-trivially on degrees of freedom in $A'$ and $B'$. We then optimize by performing iterative steps along directions in $\mathcal{M}$ which locally minimize $S_{AA'}$~[29,70]
\begin{align}
    \label{eq.single_step}
    U_{n+1}=U_n\mathrm{e}^{t K_n}\,.
\end{align}
Here, $K_n=\sum_\mu\mathcal{F}^\mu(U_n)\Xi_\mu/||\mathcal{F}||^2$ and $\mathcal{F}^\mu:~\mathcal{M}\to\mathbb{R}$ is the gradient descent vector field
\begin{align}
    \mathcal{F}^\mu(U)=-\frac{\partial}{\partial s}S(U\mathrm{e}^{s\Xi_\mu}\rho_\id\mathrm{e}^{-s\Xi_\mu}U^{-1})|_{s=0}
\end{align}
with $\{\Xi_\mu\}$ as basis of $\mathfrak{u}(2^{N_{A'}+N_{B'}})$. We choose $U_0=\id$ and we pick $0<t<1$ in such a way that the value of $S_{AA'}$ decreases with successive steps.

\begin{table}[t]
    \centering
    \renewcommand{\arraystretch}{1.25}
\begin{tabular}{@{} @{\hskip 7pt}c @{\hskip 7pt}c @{\hskip 7pt}| @{\hskip 7pt} c @{\hskip 7pt}| c @{\hskip 15pt} c c @{\hskip 15pt} c c c @{}
}
\toprule
&& & \multicolumn{6}{c}{$N_{A'}+N_{B'}$} \\
&& & $1+1$ & $1+2$ & $2+1$ & $1+3$ & $2+2$ & $3+1$
\\

\colrule
\parbox[t]{2mm}{\multirow{8}{*}{\rotatebox[origin=c]{90}{$N_A+N_B$}}} & \parbox[t]{2mm}{\multirow{4}{*}{\rotatebox[origin=c]{90}{$1+1$}}} & $d=\delta\;$ & \cellcolor{yellow!35}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$} & \cellcolor{yellow!10}{$0.382$}
\\
&& $d=2\delta$ & \cellcolor{yellow!35}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$} & \cellcolor{yellow!10}{$0.333$}
\\
&& $d=3\delta$ & \cellcolor{yellow!35}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$} & \cellcolor{yellow!10}{$0.306$}
\\
&& $d=2\delta$ & \cellcolor{yellow!35}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$} & \cellcolor{yellow!10}{$0.292$}
\\\hhline{*{1}{~}*{8}{-}}

& \parbox[t]{2mm}{\multirow{4}{*}{\rotatebox[origin=c]{90}{$1+2$}}} & $d=\delta\;$ & \multirow{4}{*}{n.a.} & \cellcolor{yellow!35}{$0.412$} & $0.438$ & \cellcolor{yellow!10}{$0.412$} & \cellcolor{yellow!10}{$0.412$} & $0.440$
\\
&& $d=2\delta$ && \cellcolor{yellow!35}{$0.368$} & $0.412$ & \cellcolor{yellow!10}{$0.368$} & \cellcolor{yellow!10}{$0.368$} & $0.415$
\\
&& $d=3\delta$ && \cellcolor{yellow!35}{$0.345$} & $0.394$ & \cellcolor{yellow!10}{$0.345$} & \cellcolor{yellow!10}{$0.345$} & $0.398$
\\
&& $d=4\delta$ && \cellcolor{yellow!35}{$0.335$} & $0.385$ & \cellcolor{yellow!10}{$0.335$} & \cellcolor{yellow!10}{$0.335$} & $0.389$
\\\botrule

\end{tabular}
    \ccaption{Numerical evidence for optimality of certain minimal purifications}{The table shows the values of the optimization for different choices of the system dimensions and of $d$. The true EoP values (the minimum optimization values) are highlighted in yellow, with the darker shade indicating the lowest-dimensional purification for which the EoP is obtained.}
\label{tab:non-minimal purifications}
\end{table}

The $\{\Xi_\mu\}$ span the tangent space at $U=\id$ and, due to the left-invariance of the Riemannian metric on $\mathcal{M}$, form  orthonormal bases for the tangent spaces at all other points in $\mathcal{M}$, too, where $\Xi_\mu$ is identified with the tangent vector to the curve $\gamma(s)=U\mathrm{e}^{s\Xi_\mu}$ at $\gamma(0)$~[70]
This saves us having to re-evaluate the matrix representation of the metric at each step, as we would have to if we had chosen a coordinate parametrisation of $\mathcal{M}$. While this makes our algorithm more efficient than a naive gradient descent, the numerically accessible range is still highly limited: since $N_{A'}+N_{B'}\geq N_A+N_B$, the dimension of $\mathcal{M}$ is at least $\mathrm{dim}\mathfrak{u}(2^{N_{A'}+N_{B'}})=2^{2N_{A'}+2N_{B'}}-1$ and~\eqref{eq.single_step} requires exponentiation of at least $2^{(N_A+N_B)}\times2^{(N_A+N_B)}$ matrices, with a typical step count of several hundred. This becomes extremely slow on a powerful desktop computer for $N_{A'}+N_{B'}\geq 5$. For the symmetric purifications in the main text this corresponds with $w>2\,\delta$, which explains the regime we were able to explore.

Given this limitation on our numerical capabilities, it is instructive to ask whether an optimization over \textsl{minimal purifications} corresponding with $N_{A'}+N_{B'}=N_A+N_B$ yields the true minimum of EE~\textendash~not least because for large systems this becomes the only numerically viable choice. A natural follow-up question is whether among the choices of minimal purifications, the intuitive choice of $N_{A'}=N_A$ and $N_{B'}=N_B$ suffices to reach the true minimum defined as EoP. More pertinently, we might ask whether it is even possible to reach the true minimum with a minimal purification for which $N_{A'}\neq N_A$ and $N_{B'}\neq N_B$. In~[70],
a combination of numerical and analytical evidence was provided to show that the answer to this question is in affirmative for Gaussian states. While limited by the greater numerical challenge in the non-Gaussian case, we present similar numerical evidence in Table~\ref{tab:non-minimal purifications} to show that the same may be said for our model: the true minimum can only be reached if $N_{A'}\geq N_{A}$ and $N_{B'}\geq N_{B}$, which indicates that the lowest-dimensional purification for which the EoP can be obtained is the minimal purification with $N_{A'}=N_A$ and $N_{B'}=N_B$.