$\displaystyle\left(\begin{array}[]{ccc}0&G^{{}^{\pm}}_{0,j}+4\Re\big{[}\underset{m\in\mathbb{E}}{\sum}g_{m,0}^{*}(h_{m,j}+g_{m,j})\big{]}&G^{{}^{\pm}}_{{}_{\\!0,L+1}}+8\Re\big{[}\underset{m\in\mathbb{E}}{\sum}g_{m,0}^{*}g_{m,L+1}\big{]}\\\
0&G^{{}^{\pm}}_{i,j}+2\Re\big{[}\underset{m\in\mathbb{E}}{\sum}(g_{m,0}^{*}-h^{*}_{m,0})(h_{m,j}+g_{m,j})\big{]}&\quad
G^{{}^{\pm}}_{i,L+1}+4\Re[\underset{m\in\mathbb{E}}{\sum}(g_{m,i}^{*}-h^{*}_{m,i})g_{m,L+1}]\\\
0&0&0\end{array}\right).$ (A.34)

##  Appendix B Reduced density matrix calculations

In this appendix, we are presenting the calculations of (108) in details.
Starting from (104), we can write

$\displaystyle\rho_{\mathbf{1}}^{\beta}\\!(\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=$
$\displaystyle\int\prod_{l\in\mathbf{2}}\\!\text{d}\bar{\xi_{l}}\text{d}\xi_{l}\;e^{-\underset{n\in\mathbf{2}}{\sum}\bar{\xi}_{n}\xi_{n}}\left\langle\xi_{1},\cdots,\xi_{k},-\xi_{k+1},\cdots,-\xi_{L}\right|\rho^{\beta}\left|\xi_{1}^{\prime},\cdots,\xi_{k}^{\prime},\xi_{k+1},\cdots,\xi_{L}\right\rangle$
(B.1) $\displaystyle=$
$\displaystyle|C^{\beta}|^{2}e^{\tfrac{1}{2}(\mathbf{R}_{11})_{ij}\bar{\xi_{i}}\bar{\xi}_{j}-\tfrac{1}{2}(\mathbf{R}_{11}^{*})_{ji}\xi_{j}^{\prime}\xi_{i}^{\prime}}\int\prod_{l\in\mathbf{2}}\\!\text{d}\bar{\xi_{l}}\text{d}\xi_{l}\;e^{-\bar{\xi}_{n}\xi_{n}}\mathcal{F}(\\{\xi_{n}\\},\\{\xi_{i}^{\prime}\\},\\{\bar{\xi}_{j}\\},\\{\bar{\xi}_{m}\\})$
$\displaystyle\times\big{[}e^{-\tfrac{1}{2}(\mathbf{R}_{12})_{in}\bar{\xi}_{i}\bar{\xi}_{n}-\tfrac{1}{2}(\mathbf{R}_{21})_{ni}\bar{\xi}_{n}\bar{\xi}_{i}}e^{\tfrac{1}{2}(\mathbf{R}_{22})_{mn}\bar{\xi}_{m}\bar{\xi}_{n}}e^{-\tfrac{1}{2}(\mathbf{R}_{22}^{*})_{nm}\xi_{n}\xi_{m}+\tfrac{1}{2}(\mathbf{R}_{12}^{*})_{jm}\xi^{\prime}_{j}\xi_{m}+\tfrac{1}{2}(\mathbf{R}_{21}^{*})_{mj}\xi_{m}\xi_{j}^{\prime}}\big{]},$

where in the above Einstein summation convention is used. To clarify the
notation above, the indices $i,j$ belong to subsystem $\mathbf{1}$ and indices
$n,m$ to subsystem $\mathbf{2}$. Eventually, the function $\mathcal{F}$ is
given by

$\displaystyle\mathcal{F}(\\{\xi_{n}\\},\\{\xi_{i}^{\prime}\\}\\!\\!\\!\\!\\!\\!$
$\displaystyle\\!\\!\\!\\!,\\{\bar{\xi}_{j}\\},\\{\bar{\xi}_{m}\\})=1+\beta(\boldsymbol{\mathfrak{M}}_{1})_{0j}\bar{\xi}_{j}-\beta(\boldsymbol{\mathfrak{M}}_{2})_{0m}\bar{\xi}_{m}+\beta^{*}(\boldsymbol{\mathfrak{M}}_{1}^{*})_{0i}\xi_{i}^{\prime}$
(B.2)
$\displaystyle+\beta^{*}(\boldsymbol{\mathfrak{M}}_{2}^{*})_{0n}\xi_{n}+|\beta|^{2}(\boldsymbol{\mathfrak{M}}_{1})_{0j}(\boldsymbol{\mathfrak{M}}_{1}^{*})_{0i}\bar{\xi}_{j}\xi_{i}^{\prime}+|\beta|^{2}(\boldsymbol{\mathfrak{M}}_{1})_{0j}(\boldsymbol{\mathfrak{M}}_{2}^{*})_{0n}\bar{\xi}_{j}\xi_{n}$
$\displaystyle-|\beta|^{2}(\boldsymbol{\mathfrak{M}}_{2})_{0m}(\boldsymbol{\mathfrak{M}}_{1}^{*})_{0i}\bar{\xi}_{m}\xi_{i}^{\prime}-|\beta|^{2}(\boldsymbol{\mathfrak{M}}_{2})_{0m}(\boldsymbol{\mathfrak{M}}_{2}^{*})_{0n}\bar{\xi}_{m}\xi_{n}.$

In (B.1) we have divided $\mathbf{R}$ into four submatrices $\mathbf{R_{11}}$,
$\mathbf{R_{12}}$, $\mathbf{R_{21}}=-\mathbf{R_{12}}^{T}$ and
$\mathbf{R_{22}}$, according to the part we are tracing out (or not). Although
these submatrices do not need to have same size, $\mathbf{R_{11}}$ and
$\mathbf{R_{22}}$ should be square matrices. We write the (B.1) in the compact
form:

$\rho_{\mathbf{1}}(\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=|C^{\beta}|^{2}e^{\tfrac{1}{2}\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\\!\boldsymbol{\xi}^{\prime}\\\
\end{matrix}\right)\left(\begin{matrix}\mathbf{R}_{11}&0\\\
0&-\mathbf{R}_{11}^{*}\end{matrix}\right)\left(\begin{matrix}\boldsymbol{\bar{\xi}}\\\
\boldsymbol{\xi}^{\prime}\end{matrix}\right)}\int\text{{D}}\boldsymbol{\eta}\mathcal{F}\\!(\boldsymbol{\xi_{i}^{\prime}},\boldsymbol{\bar{\xi}_{j}},\boldsymbol{\eta})e^{\tfrac{1}{2}\boldsymbol{\eta}^{T}\boldsymbol{\mathcal{A}}\boldsymbol{\eta}+\boldsymbol{\lambda}^{T}\\!\boldsymbol{\eta}}.$
(B.3)

where $\textbf{D}\boldsymbol{\eta}=\text{d}\bar{\xi}\text{d}\xi$ and

$\displaystyle\boldsymbol{\eta}^{T}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\boldsymbol{\xi}\\\
\end{matrix}\right),\qquad\boldsymbol{\lambda}^{T}=\left(\begin{matrix}-\boldsymbol{\bar{\xi}}\mathbf{R}_{12}&-\boldsymbol{\xi}^{\prime}\mathbf{R}_{12}^{*}\\\
\end{matrix}\right),$ (B.4) $\displaystyle\boldsymbol{\mathcal{A}}$
$\displaystyle=$
$\displaystyle\left(\begin{matrix}\mathbf{R}_{22}&-\mathbf{I}\\\
\mathbf{I}&-\mathbf{R}_{22}^{*}\\\ \end{matrix}\right).$ (B.5)

It is much easier to solve the integration in (B.3) with this new variables.
In the expression of $\mathcal{F}$, the first, second, forth and sixth terms
do not depend on the variables of integration. Therefore, for those terms we
can write:

$\rho_{\mathbf{1}}^{\beta}(1^{\text{st}},2^{\text{nd}},4^{\text{th}},6^{\text{th}},\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=|C^{\beta}|(\text{det}\big{[}\boldsymbol{\mathcal{A}}\big{]})^{\tfrac{1}{2}}\big{(}1+\beta\boldsymbol{\mathfrak{M}}_{1}\boldsymbol{\bar{\xi}}\big{)}\big{(}1+\beta^{*}\boldsymbol{\mathfrak{M}}_{1}^{*}\boldsymbol{\xi^{\prime}}\big{)}e^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)},$ (B.6)

where

$\boldsymbol{\Omega}=\left(\begin{matrix}\mathbf{R}_{11}&0\\\
0&-{\mathbf{R}_{11}}^{*}\\\
\end{matrix}\right)+\left(\begin{matrix}\mathbf{R}_{12}&0\\\
0&{\mathbf{R}_{12}}^{*}\\\
\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\left(\begin{matrix}{\mathbf{R}_{12}}^{T}&0\\\
0&{\mathbf{R}_{12}}^{\dagger}\\\ \end{matrix}\right).$ (B.7)

For the terms with linear Grassmann variables in the integration we first
substitute
$\boldsymbol{\eta}\rightarrow\boldsymbol{\eta}+\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda}$
in (B.3). Using the Berezin integration techniques in presented in [56], we
are left with

$\rho_{\mathbf{1}}^{\beta}(3^{\text{rd}},5^{\text{th}},7^{\text{th}},8^{\text{th}},\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=|C^{\beta}|(\text{det}\big{[}\boldsymbol{\mathcal{A}}\big{]})^{\tfrac{1}{2}}(\text{det}\big{[}\boldsymbol{\mathcal{W}}\big{]})^{\tfrac{1}{2}}F(\boldsymbol{\bar{\xi}}_{\mathbf{1}},\boldsymbol{\xi^{\prime}}_{\mathbf{1}})e^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)},$ (B.8)

with

$\displaystyle
F(\boldsymbol{\bar{\xi}}_{\mathbf{1}},\boldsymbol{\xi^{\prime}}_{\mathbf{1}})$
$\displaystyle=$
$\displaystyle|\beta|^{2}\text{Pf}\big{[}\boldsymbol{\mathcal{W}}\big{]}+\left(\begin{matrix}-\beta\boldsymbol{\mathfrak{M}}_{2}&0\\\
0&\beta^{*}\boldsymbol{\mathfrak{M}}_{2}^{*}\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda}+|\beta|^{2}\boldsymbol{\mathfrak{M}}_{1}\boldsymbol{\bar{\xi}}\times\left(\begin{matrix}0&0\\\
0&\boldsymbol{\mathfrak{M}}_{2}^{*}\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda}$
$\displaystyle+|\beta|^{2}\left(\begin{matrix}-\boldsymbol{\mathfrak{M}}_{2}&0\\\
0&0\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda}\times\boldsymbol{\mathfrak{M}}_{1}^{*}\boldsymbol{\xi^{\prime}}+|\beta|^{2}\left(\begin{matrix}-\boldsymbol{\mathfrak{M}_{2}}&0\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda}\times\left(\begin{matrix}0&\boldsymbol{\mathfrak{M}_{2}}^{*}\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\lambda},$
$\displaystyle\boldsymbol{\mathcal{W}}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\boldsymbol{\mathfrak{M}}_{2}&\mathbf{0}\\\
\mathbf{0}&-\boldsymbol{\mathfrak{M}}_{2}^{*}\\\
\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-T}\left(\begin{matrix}\boldsymbol{\mathfrak{M}}_{2}^{T}&\mathbf{0}\\\
\mathbf{0}&-\boldsymbol{\mathfrak{M}}_{2}^{\dagger}\\\ \end{matrix}\right).$
(B.10)

Putting (B.6) and (B.8) together, after some algebraic manipulations, we get
to (105).

We could have used another trick to solve the Berezin integrations to get RDM.
This method ends up in having two exponentials in the final result. The trick
is to write the linear Grassmann variable in the integration as exponentials,
like:

$(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})_{k_{1}}\cdots(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}}\big{)}_{k_{N}}\\!\big{(}\boldsymbol{\mathfrak{M}}^{*}\boldsymbol{\xi}^{\prime}\big{)}_{k_{1}}\cdots(\boldsymbol{\mathfrak{M}}^{*}\boldsymbol{\xi}^{\prime}\big{)}_{k_{N}}=\int\prod_{i}\text{d}\bar{\theta}_{i}\text{d}\theta_{i}e^{\sum_{k_{j}}\theta_{k_{j}}(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})_{k_{j}}+\bar{\theta}_{k_{j}}(\boldsymbol{\mathfrak{M}}^{*}\boldsymbol{\xi})_{k_{j}}}.$
(B.11)

In the above relation, $\theta$ and $\bar{\theta}$ are Grassmann variables
too.

If we start again from density matrix expression in coherent basis (103), and
rewrite it as:

$\displaystyle\rho^{\beta}\\!(\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=$
$\displaystyle|C^{\beta}|^{2}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}R_{ij}\bar{\xi}_{i}\bar{\xi}_{j}}}(1+\beta\mathfrak{M}_{0k}\bar{\xi}_{k}+\beta^{*}\mathfrak{M}_{0l}^{*}\xi_{l}^{\prime}+|\beta|^{2}\mathfrak{M}_{0k}\mathfrak{M}_{0l}^{*}\bar{\xi}_{k}\xi_{l}^{\prime})\scalebox{1.3}{$e$}^{{}^{-\tfrac{1}{2}R^{*}_{nm}\xi_{n}^{\prime}\xi_{m}^{\prime}}}$
(B.12) $\displaystyle=$
$\displaystyle\rho^{\beta}\\!(1,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})+\varrho^{\beta}\\!(2,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})+\varrho^{\beta}\\!(3,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})+\rho^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime}),$

where the $\rho^{\beta}\\!(1,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})$ and
$\rho^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})$ are the
density matrices corresponding to the vacuum-vacuum state and excited-excited
state, and the
$\varrho^{\beta}\\!(2,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})$ and
$\varrho^{\beta}\\!(3,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})$ are density
cross terms corresponding to the vacuum-excited (and exited-vacuum) state
terms. We are going to partial trace each term separately and then put the
results together afterwards. For the first term in (B.12), we have:

$\rho_{\mathbf{1}}^{\beta}(1,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=\frac{\mathcal{C}^{\beta}}{|\beta|^{2}}e^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}^{\beta}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)},$ (B.13)

where again

$\mathcal{C}^{\beta}=\frac{|\beta|^{2}\sqrt{\det\big{[}\mathbf{I}+{\mathbf{R}_{22}}^{\dagger}\mathbf{R}_{22}\big{]}}}{(1+|\beta|^{2})\sqrt{\det\big{[}\mathbf{I}+{\mathbf{R}}^{\dagger}\mathbf{R}\big{]}}},\qquad\boldsymbol{\mathcal{A}}=\left(\begin{matrix}\mathbf{R}_{22}&-\mathbf{I}\\\
\mathbf{I}&-\mathbf{R}_{22}^{*}\\\ \end{matrix}\right),$ (B.14a)
$\boldsymbol{\Omega}=\left(\begin{matrix}\mathbf{R}_{11}&0\\\
0&-{\mathbf{R}_{11}}^{*}\\\
\end{matrix}\right)+\left(\begin{matrix}\mathbf{R}_{12}&0\\\
0&{\mathbf{R}_{12}}^{*}\\\
\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\left(\begin{matrix}{\mathbf{R}_{12}}^{T}&0\\\
0&{\mathbf{R}_{12}}^{\dagger}\\\ \end{matrix}\right).$ (B.14b)

For the last term in (B.12), we can write:

$\rho^{\beta}(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=|C^{\beta}\beta|^{2}\int\text{d}\bar{\theta}\text{d}\theta
e^{\theta(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})+\bar{\theta}(\boldsymbol{\mathfrak{M}}^{*}\boldsymbol{\xi}^{\prime})}e^{\tfrac{1}{2}R_{ij}\bar{\xi_{i}}\bar{\xi_{j}}}e^{-\tfrac{1}{2}R^{*}_{ij}\xi^{\prime}_{i}\xi^{\prime}_{j}}.$
(B.15)

For the reduced density matrix we get

$\displaystyle\rho_{\mathbf{1}}^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=$
$\displaystyle|C^{\beta}\beta|^{2}\\!\int\text{d}\bar{\theta}\text{d}\theta
e^{\tfrac{1}{2}(\mathbf{R}_{11})_{ij}\bar{\xi_{i}}\bar{\xi}_{j}-\tfrac{1}{2}(\mathbf{R}_{11}^{*})_{ji}\xi_{j}^{\prime}\xi_{i}^{\prime}}e^{\theta(\boldsymbol{\mathfrak{M}_{1}}.\bar{\boldsymbol{\xi}})+\bar{\theta}(\boldsymbol{\mathfrak{M}_{1}^{*}}.\boldsymbol{\xi}^{\prime})}\\!\int\prod_{l=k+1}^{L+1}\\!\text{d}\bar{\xi_{l}}\text{d}\xi_{l}\>e^{-\bar{\xi}_{l}\xi_{l}}$
(B.16) $\displaystyle\times
e^{-\theta(\boldsymbol{\mathfrak{M}_{2}}\bar{\boldsymbol{\xi}})+\bar{\theta}(\boldsymbol{\mathfrak{M}_{2}^{*}}\boldsymbol{\xi})}\big{[}e^{-(\mathbf{R}_{12})_{in}\bar{\xi}_{i}\bar{\xi}_{n}+\tfrac{1}{2}(\mathbf{R}_{22})_{mn}\bar{\xi}_{m}\bar{\xi}_{n}}e^{-\tfrac{1}{2}(\mathbf{R}_{22}^{*})_{nm}\xi_{n}\xi_{m}+(\mathbf{R}_{12}^{*})_{jm}\xi^{\prime}_{j}\xi_{m}}\big{]}$

We denote the terms on the left of second integral as
$\mathcal{F}(\bar{\boldsymbol{\xi}},\boldsymbol{\xi},\theta,\bar{\theta})$. If
we introduce new Grassmann variables
$\boldsymbol{\eta}^{T}=\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\boldsymbol{\xi}\\\
\end{matrix}\right)$ and
$\boldsymbol{\lambda}^{T}=\left(\begin{matrix}-\boldsymbol{\bar{\xi}}\mathbf{R}_{12}-\theta(\boldsymbol{\mathfrak{M}})_{0}&-\boldsymbol{\xi}^{\prime}\mathbf{R}_{12}^{*}+\bar{\theta}(\boldsymbol{\mathfrak{M}^{*}})_{0}\\\
\end{matrix}\right)$, then we can write the integral as

$\rho_{{}_{\mathbf{1}}}^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=\int\text{d}\bar{\theta}\text{d}\theta\mathcal{F}(\bar{\boldsymbol{\xi}},\boldsymbol{\xi}^{\prime},\theta,\bar{\theta})\int\text{{D}}\boldsymbol{\eta}\>e^{\tfrac{1}{2}\boldsymbol{\eta}^{T}\\!\boldsymbol{\mathcal{A}}\boldsymbol{\eta}+\boldsymbol{\lambda}^{T}\\!\boldsymbol{\eta}}$
(B.17)

where $\boldsymbol{\mathcal{A}}$ has been defined previously. Solving the
first integral we get:

$\rho_{{}_{\mathbf{1}}}^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=\mathcal{C}^{\beta}e^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)}\int\text{{D}}\boldsymbol{\Theta}e^{\tfrac{1}{2}\boldsymbol{\Theta}^{T}\boldsymbol{\omega}\boldsymbol{\Theta}+\boldsymbol{\eta^{\prime}}^{T}\boldsymbol{\mathcal{J}}\boldsymbol{\Theta}}$
(B.18)

where $\boldsymbol{\Theta}^{T}=\left(\begin{matrix}\bar{\theta}&\\!\theta\\\
\end{matrix}\right)$,
$\textbf{D}\boldsymbol{\Theta}=\text{d}\bar{\theta}\text{d}\theta$,
$\boldsymbol{\eta^{\prime}}^{T}\\!=\\!\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\\!\boldsymbol{\xi}^{\prime}\\\
\end{matrix}\right)$, and $\boldsymbol{\Omega}$ and $\mathcal{C}^{\beta}$ are
given by (106b) and (106a) respectively. Finally, if we take the last integral
in the equation (B.18), we get

$\rho_{{}_{\mathbf{1}}}^{\beta}\\!(4,\boldsymbol{\xi},\boldsymbol{\xi}^{\prime})=\mathcal{C}^{\beta}\text{Pf}[\boldsymbol{\omega}]e^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}^{\prime}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)}$ (B.19)

where

$\boldsymbol{\Omega}^{\prime}=\boldsymbol{\Omega}+\boldsymbol{\mathcal{J}}\boldsymbol{\omega}^{-1}\boldsymbol{\mathcal{J}}^{T}$
(B.20a)
$\boldsymbol{\omega}=\left(\begin{matrix}0&(\boldsymbol{\mathfrak{M}_{2}^{*}})_{0}\\\
-(\boldsymbol{\mathfrak{M}}_{2})_{0}&0\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\left(\begin{matrix}0&-(\boldsymbol{\mathfrak{M}_{2}})_{0}^{T}\\\
(\boldsymbol{\mathfrak{M}}_{2})_{0}^{\dagger}&0\end{matrix}\right),$ (B.20b)
$\boldsymbol{\mathcal{J}}=\left(\begin{matrix}-\mathbf{R}_{12}&\mathbf{0}\\\
\mathbf{0}&-{\mathbf{R}_{12}}^{*}\\\
\end{matrix}\right)\boldsymbol{\mathcal{A}}^{-1}\left(\begin{matrix}\mathbf{0}&\boldsymbol{\mathfrak{M}}_{2}^{*}\\\
-\boldsymbol{\mathfrak{M}}_{2}&\mathbf{0}\\\
\end{matrix}\right)-\left(\begin{matrix}\mathbf{0}&\boldsymbol{\mathfrak{M}}_{1}^{*}\\\
\boldsymbol{\mathfrak{M}}_{1}&\mathbf{0}\\\ \end{matrix}\right).$ (B.20c)

To do the same calculation for cross term
$\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})$ in equation
(B.12), we have

$\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=|C^{\beta}|^{2}\beta\>\boldsymbol{\mathfrak{M}}_{{}_{0l}}\bar{\xi}_{{}_{l}}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}R_{ij}\bar{\xi}_{i}\bar{\xi}_{j}-\tfrac{1}{2}R_{ij}^{*}\xi_{i}^{\prime}\xi_{j}^{\prime}}}$
(B.21)

If we use the equation (B.11) and write the $\xi$ behind the exponential as
$(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})_{0}=\int\\!\text{d}\theta
e^{\theta(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})_{0}}=\int\\!\text{d}\theta
e^{\theta\boldsymbol{\mathfrak{M}}_{0l}\bar{\xi}_{l}}$, where $\theta$ is a
new Grassmann variables. We can write the density matrix as

$\displaystyle\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=$
$\displaystyle|C^{\beta}|^{2}\beta\int\text{d}\theta
e^{\theta(\boldsymbol{\mathfrak{M}}\bar{\boldsymbol{\xi}})_{0}}e^{\tfrac{1}{2}R_{ij}\bar{\xi_{i}}\bar{\xi_{j}}}e^{-\tfrac{1}{2}R^{*}_{ij}\xi^{\prime}_{i}\xi^{\prime}_{j}}.$
(B.22)

Equivalently, If we introduce new Grassmann variables
$\boldsymbol{\eta}^{T}\\!=\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\boldsymbol{\xi}\\\
\end{matrix}\right)$ then we can write the integral as

$\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=\int\text{d}\theta\mathcal{F}(\bar{\boldsymbol{\xi}},\boldsymbol{\xi}^{\prime},\theta)\int\text{{D}}\boldsymbol{\eta}\>\scalebox{1.2}{$e$}^{\tfrac{1}{2}\boldsymbol{\eta}^{T}\\!\boldsymbol{\mathcal{A}}\boldsymbol{\eta}+\boldsymbol{\lambda}^{T}\\!\boldsymbol{\eta}},$
(B.23)

where we have denoted the terms on the left of second integral by
$\mathcal{F}(\bar{\boldsymbol{\xi}},\boldsymbol{\xi},\theta)$ and

$\boldsymbol{\lambda}^{T}=\left(\begin{matrix}-\boldsymbol{\bar{\xi}}\mathbf{R}_{12}-\theta(\boldsymbol{\mathfrak{M}_{2}})_{0l^{\prime}}&-\boldsymbol{\xi}^{\prime}\mathbf{R}_{12}^{*}\\\
\end{matrix}\right)$ (B.24)

Taking the integration on $\boldsymbol{\eta}$, we have:

$\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=\frac{\mathcal{C}^{\beta}}{\beta^{*}}\scalebox{1.3}{$e$}^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)}\int\text{d}\theta\scalebox{1.2}{$e$}^{\theta\boldsymbol{\zeta}}$
(B.25)

where $\boldsymbol{\Omega}$ is the same as (106b) and $\mathcal{C}^{\beta}$ is
given by $\eqref{eq:BetaCcalConst}$. Then, we can take the integral over the
$\theta$ above and write the final result as

$\varrho^{\beta}(2,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=\frac{\mathcal{C}^{\beta}}{\beta^{*}}\big{(}\boldsymbol{\Upsilon}\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\mathcal{R}}^{T}+\boldsymbol{\Upsilon^{\prime}}\big{)}\boldsymbol{\eta^{\prime}}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}\boldsymbol{\eta^{\prime}}^{T}\boldsymbol{\Omega}\boldsymbol{\eta}^{\prime}}}$
(B.26)

where
$\boldsymbol{\eta^{\prime}}^{T}=\left(\begin{matrix}\boldsymbol{\bar{\xi}}&\\!\boldsymbol{\xi}^{\prime}\\\
\end{matrix}\right)$, and Also we have

$\boldsymbol{\Upsilon}=\begin{bmatrix}\boldsymbol{\mathfrak{M}_{2}}&\mathbf{0}\end{bmatrix}_{{}_{1\times
2L^{\prime}}}\\!\qquad\boldsymbol{\Upsilon}^{\prime}=\begin{bmatrix}\boldsymbol{\mathfrak{M}_{1}}&\mathbf{0}\end{bmatrix}_{{}_{1\times
2l}}\\!\qquad\boldsymbol{\mathcal{R}}=\begin{bmatrix}\mathbf{R}_{12}&0\\\
0&\mathbf{R}_{12}^{*}\end{bmatrix}_{{}_{2l\times 2L^{\prime}}}$ (B.27)

For the other cross density (third term in equation (B.12)) we have:

$\varrho^{\beta}(3,\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=\frac{\mathcal{C}^{\beta}}{\beta}\big{(}-\boldsymbol{\Upsilon}\boldsymbol{\mathcal{A}}^{-1}\boldsymbol{\mathcal{R}}^{T}+\boldsymbol{\Upsilon^{\prime}}\big{)}\boldsymbol{\eta^{\prime}}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}\boldsymbol{\eta^{\prime}}^{T}\boldsymbol{\Omega}^{\beta}\boldsymbol{\eta}^{\prime}}}$
(B.28)

which

$\boldsymbol{\Upsilon}=\begin{bmatrix}\mathbf{0}&\boldsymbol{\mathfrak{M}_{2}}^{*}\end{bmatrix}_{{}_{1\times
2L^{\prime}}}\\!\qquad\boldsymbol{\Upsilon}^{\prime}=\begin{bmatrix}\mathbf{0}&\boldsymbol{\mathfrak{M}_{1}}^{*}\end{bmatrix}\\!\qquad\boldsymbol{\mathcal{R}}=\begin{bmatrix}\mathbf{R}_{12}&0\\\
0&\mathbf{R}_{12}^{*}\end{bmatrix}_{{}_{2l\times 2L^{\prime}}}$ (B.29)

Putting equations (B.13), (B.19), (B.26) and (B.28) together and moving and
reordering some parts we can write the RDM for the state
$\left|\beta\right\rangle$ as

$\displaystyle\rho_{{}_{\mathbf{1}}}^{\beta}(\boldsymbol{\xi},\boldsymbol{\xi^{\prime}})=$
$\displaystyle\;\mathcal{C}_{\beta}\Big{[}\frac{1}{|\beta|^{2}}+\text{Pf}[\boldsymbol{\Omega_{2}}]\scalebox{1.3}{$e$}^{\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\xi^{\prime}\\\
\end{matrix}\right)\boldsymbol{\mathcal{J}}\boldsymbol{\omega}^{-1}\boldsymbol{\mathcal{J}}^{T}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)}+\boldsymbol{\mathcal{K}}^{\beta}\left(\begin{matrix}\bar{\xi}\\\
\xi^{\prime}\end{matrix}\right)\Big{]}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{11})_{ij}\bar{\xi}_{i}\bar{\xi}_{j}}}\scalebox{1.3}{$e$}^{{}^{\mathcal{Y}_{ij}\bar{\xi_{i}}\xi_{j}^{\prime}}}\scalebox{1.3}{$e$}^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{22})_{ij}\xi_{i}^{\prime}\xi_{j}^{\prime}}},$
(B.30)

where
$2\boldsymbol{\mathcal{Y}}=\boldsymbol{\Omega}_{12}-{\boldsymbol{\Omega}_{21}}^{T}$
and

$\boldsymbol{\mathcal{K}}^{\beta}=\left(\begin{matrix}\tfrac{1}{\beta^{*}}\boldsymbol{\mathfrak{M}_{1}}&\tfrac{1}{\beta}\boldsymbol{\mathfrak{M}_{1}}^{*}\\\
\end{matrix}\right)+\left(\begin{matrix}\tfrac{1}{\beta^{*}}\boldsymbol{\mathfrak{M}_{2}}&-\tfrac{1}{\beta}\boldsymbol{\mathfrak{M}_{2}}^{*}\\\
\end{matrix}\right)\scalebox{1.2}{${\boldsymbol{\mathcal{A}}}$}^{-1}\left(\begin{matrix}{\mathbf{R}_{12}}^{T}&0\\\
0&{\mathbf{R}_{12}}^{\dagger}\\\ \end{matrix}\right)$ (B.31)

To write the RDM in operator format, we have to reorder it. Then, the final
result is

$\displaystyle\rho_{\mathbf{1}}^{\beta}(c,c^{\dagger})=$
$\displaystyle\>\mathcal{C}_{\beta}e^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{11})_{ij}c_{i}^{\dagger}c_{j}^{\dagger}}}\Big{[}\frac{1}{|\beta|^{2}}e^{{}^{(\ln\frac{\boldsymbol{\Omega}_{12}-{\boldsymbol{\Omega}_{21}}^{T}}{2})_{ij}c_{i}^{\dagger}c_{j}}}+(\boldsymbol{\mathcal{K}}^{\beta}_{1})_{k}c^{\dagger}_{k}e^{{}^{(\ln\frac{\boldsymbol{\Omega}_{12}-{\boldsymbol{\Omega}_{21}}^{T}}{2})_{ij}c_{i}^{\dagger}c_{j}}}$
(B.32)
$\displaystyle+e^{{}^{(\ln\frac{\boldsymbol{\Omega}_{12}-{\boldsymbol{\Omega}_{21}}^{T}}{2})_{ij}c_{i}^{\dagger}c_{j}}}(\boldsymbol{\mathcal{K}}^{\beta}_{2})_{l}c_{l}\Big{]}e^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{22})_{ij}c_{i}c_{j}}}$
$\displaystyle+\mathcal{C}^{\beta}\text{Pf}[\boldsymbol{\Omega_{2}}]e^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{11}^{\prime})_{ij}c_{i}^{\dagger}c_{j}^{\dagger}}}e^{{}^{(\ln\frac{\boldsymbol{\Omega}_{12}^{\prime}-{\boldsymbol{\Omega}_{21}^{\prime}}^{T}}{2})_{ij}c_{i}^{\dagger}c_{j}}}e^{{}^{\tfrac{1}{2}(\boldsymbol{\Omega}_{22}^{\prime})_{ij}c_{i}c_{j}}}.$

The $\boldsymbol{\mathcal{K}}^{\beta}_{1}$ stands for the first block part of
vector $\boldsymbol{\mathcal{K}}^{\beta}$ (and the same argument for
$\boldsymbol{\mathcal{K}}^{\beta}_{2}$). To find a shorter notation for the
$\rho_{\mathbf{1}}^{\beta}(c,c^{\dagger})$ we use the relation [57]

$e^{{}^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}_{1}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}}\>e^{{}^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}_{2}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}}\>e^{{}^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}_{3}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}}=e^{{}^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}}$ (B.33)

where
$e^{\boldsymbol{\mathcal{M}}}=e^{\boldsymbol{\mathcal{M}}_{1}}\>e^{\boldsymbol{\mathcal{M}_{2}}}\>e^{\boldsymbol{\mathcal{M}_{3}}}$.
Therefore, we can write (B.32) as

$\displaystyle\rho_{\mathbf{1}}^{\beta}(c,c^{\dagger})=$
$\displaystyle\>\mathcal{C}_{\beta}\>\Big{[}\frac{1}{|\beta|^{2}}e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln\boldsymbol{\mathcal{Y}}}}+\text{Pf}[\boldsymbol{\Omega_{2}}]e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}^{\prime}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln(\boldsymbol{\mathcal{Y}}+\boldsymbol{\mathfrak{Y}})}}$
(B.34)
$\displaystyle+(\boldsymbol{\mathcal{K}}^{\beta}_{1})_{k}c^{\dagger}_{k}e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln\boldsymbol{\mathcal{Y}}}}+e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln\boldsymbol{\mathcal{Y}}}}(\boldsymbol{\mathcal{K}}^{\beta}_{2})_{l}c_{l}\Big{]}.$

To move the exponentials to one side (left or right) of the (B.34), we use the
relation (110). Therefore, we can write the (B.34) as

$\displaystyle\rho_{\mathbf{1}}^{\beta}(c,c^{\dagger})=$
$\displaystyle\>\mathcal{C}_{\beta}\>e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln(\tfrac{1}{2}\boldsymbol{\Omega}_{12}-\tfrac{1}{2}\boldsymbol{\Omega}_{21}^{T})}}\Big{[}\frac{1}{|\beta|^{2}}+\left(\begin{matrix}\boldsymbol{\mathcal{L}}^{\beta}_{1}&\boldsymbol{\mathcal{L}}^{\beta}_{2}\end{matrix}\right)\left(\begin{matrix}\mathbf{c}^{\dagger}\\\
\mathbf{c}\end{matrix}\right)\Big{]}$ (B.35)
$\displaystyle\>+\mathcal{C}_{\beta}\>\text{Pf}[\boldsymbol{\Omega_{2}}]e^{\frac{1}{2}\left(\begin{matrix}\mathbf{c}^{\dagger}&\\!\mathbf{c}\end{matrix}\right)\begin{matrix}\boldsymbol{\mathcal{M}^{\prime}}\end{matrix}\left(\begin{matrix}\\!\mathbf{c}\\\
\mathbf{c}^{\dagger}\end{matrix}\right)}e^{{}^{\tfrac{1}{2}\text{tr}\ln(\tfrac{1}{2}\boldsymbol{\Omega}_{12}^{\prime}-\tfrac{1}{2}{\boldsymbol{\Omega}_{21}^{\prime}}^{T})}}$

where we have:

$\boldsymbol{\mathcal{M}}=\ln{\left(\begin{matrix}\tfrac{1}{2}\boldsymbol{\Omega}_{12}-\tfrac{1}{2}{\boldsymbol{\Omega}_{21}}^{T}+2\boldsymbol{\Omega}_{11}(\boldsymbol{\Omega}_{12}^{T}-{\boldsymbol{\Omega}_{21}})^{-1}\boldsymbol{\Omega}_{22}&\;2\boldsymbol{\Omega}_{11}(\boldsymbol{\Omega}_{12}^{T}-{\boldsymbol{\Omega}_{21}})^{-1}\\\
2(\boldsymbol{\Omega}_{12}^{T}-\boldsymbol{\Omega}_{21})^{-1}\boldsymbol{\Omega}_{22}&2(\boldsymbol{\Omega}_{12}^{T}-{\boldsymbol{\Omega}_{21}})^{-1}\end{matrix}\right)},$
(B.36a)
$\boldsymbol{\mathcal{M}}^{\prime}=\ln\left(\begin{matrix}\tfrac{1}{2}\boldsymbol{\Omega}_{12}^{\prime}-\tfrac{1}{2}{\boldsymbol{\Omega}_{21}^{\prime}}^{T}+2\boldsymbol{\Omega}_{11}^{\prime}({\boldsymbol{\Omega}_{12}^{\prime}}^{T}-{\boldsymbol{\Omega}_{21}}^{\prime})^{-1}\boldsymbol{\Omega}_{22}^{\prime}&\;2\boldsymbol{\Omega}_{11}^{\prime}({\boldsymbol{\Omega}_{12}^{\prime}}^{T}-{\boldsymbol{\Omega}_{21}^{\prime}})^{-1}\\\
2({\boldsymbol{\Omega}_{12}^{\prime}}^{T}-\boldsymbol{\Omega}_{21}^{\prime})^{-1}\boldsymbol{\Omega}_{22}^{\prime}&2({\boldsymbol{\Omega}_{12}^{\prime}}^{T}-{\boldsymbol{\Omega}_{21}^{\prime}})^{-1}\end{matrix}\right)$
(B.36b) $\boldsymbol{\Omega}=\left(\begin{matrix}\mathbf{R}_{11}&0\\\
0&-{\mathbf{R}_{11}}^{*}\\\
\end{matrix}\right)+\left(\begin{matrix}\mathbf{R}_{12}&0\\\
0&{\mathbf{R}_{12}}^{*}\\\
\end{matrix}\right){\boldsymbol{\mathcal{A}}}^{-1}\left(\begin{matrix}{\mathbf{R}_{12}}^{T}&0\\\
0&{\mathbf{R}_{12}}^{\dagger}\\\
\end{matrix}\right);\qquad\boldsymbol{\mathcal{A}}=\left(\begin{matrix}\mathbf{R}_{22}&-\mathbf{I}\\\
\mathbf{I}&-\mathbf{R}_{22}^{*}\\\ \end{matrix}\right)$ (B.36c)
$\boldsymbol{\Omega}^{\prime}=\boldsymbol{\Omega}+\boldsymbol{\mathcal{J}}\boldsymbol{\omega}^{-1}\boldsymbol{\mathcal{J}}^{T},\qquad\boldsymbol{\omega}=\left(\begin{matrix}\boldsymbol{\mathfrak{M}_{2}^{*}}(\boldsymbol{\mathcal{A}}^{{-1}})_{22}\boldsymbol{\mathfrak{M}}_{2}^{\dagger}&-\boldsymbol{\mathfrak{M}_{2}^{*}}(\boldsymbol{\mathcal{A}}^{{-1}})_{21}\boldsymbol{\mathfrak{M}}_{2}^{T}\\\
-\boldsymbol{\mathfrak{M}}_{2}(\boldsymbol{\mathcal{A}}^{{-1}})_{12}\boldsymbol{\mathfrak{M}}_{2}^{\dagger}&\boldsymbol{\mathfrak{M}_{2}^{*}}(\boldsymbol{\mathcal{A}}^{{-1}})_{11}\boldsymbol{\mathfrak{M}}_{2}^{T}\end{matrix}\right)$
(B.36d)
$\boldsymbol{\mathcal{J}}=\left(\begin{matrix}-\mathbf{R}_{12}&\mathbf{0}\\\
\mathbf{0}&-{\mathbf{R}_{12}}^{*}\\\
\end{matrix}\right){\boldsymbol{\mathcal{A}}}^{-1}\left(\begin{matrix}\mathbf{0}&\boldsymbol{\mathfrak{M}}_{2}^{*}\\\
-\boldsymbol{\mathfrak{M}}_{2}&\mathbf{0}\\\
\end{matrix}\right)-\left(\begin{matrix}\mathbf{0}&\boldsymbol{\mathfrak{M}}_{1}^{*}\\\
\boldsymbol{\mathfrak{M}}_{1}&\mathbf{0}\\\ \end{matrix}\right),$ (B.36e)
$\boldsymbol{\mathcal{K}}^{\beta}=\left(\begin{matrix}\tfrac{1}{\beta^{*}}\boldsymbol{\mathfrak{M}_{1}}&\tfrac{1}{\beta}\boldsymbol{\mathfrak{M}_{1}}^{*}\\\
\end{matrix}\right)+\left(\begin{matrix}\tfrac{1}{\beta^{*}}\boldsymbol{\mathfrak{M}_{2}}&-\tfrac{1}{\beta}\boldsymbol{\mathfrak{M}_{2}}^{*}\\\
\end{matrix}\right){\boldsymbol{\mathcal{A}}}^{-1}\left(\begin{matrix}{\mathbf{R}_{12}}^{T}&0\\\
0&{\mathbf{R}_{12}}^{\dagger}\\\ \end{matrix}\right)$ (B.36f)
$\boldsymbol{\mathcal{L}}^{\beta}_{1}=2\boldsymbol{\mathcal{K}}^{\beta}_{1}(\boldsymbol{\Omega}_{12}^{T}-{\boldsymbol{\Omega}_{21}})^{-1},\qquad\boldsymbol{\mathcal{L}}_{2}^{\beta}=2\boldsymbol{\mathcal{K}}^{\beta}_{1}(\boldsymbol{\Omega}_{12}^{T}-{\boldsymbol{\Omega}_{21}})^{-1}\boldsymbol{\Omega}_{22}+\boldsymbol{\mathcal{K}}^{\beta}_{2}.$
(B.36g)

$\boldsymbol{\mathcal{K}}^{\beta}_{1}$ stands for the first block part of
vector $\boldsymbol{\mathcal{K}}^{\beta}$. Since, we can think of
$\boldsymbol{\mathcal{K}}^{\beta}$ as a vector made from two block vectors
(and the same argument for $\boldsymbol{\mathcal{K}}^{\beta}_{2}$).

##  Appendix C Entanglement calculations of general $\beta$ parity broken
state

To calculate the $n=2$ Rényi EE, we start by:

$\text{tr}({\rho_{A}^{\beta}}^{2})=\text{tr}(\rho_{A}^{\beta}\mathbf{I}\rho_{A}^{\beta}).$
(C.1)

We can use the Identity resolution and tracing formula of Grassmann variables
to write:

$\mathbf{I}=\int\prod_{l}\text{d}\bar{\xi}_{l}\text{d}\xi_{l}e^{-\boldsymbol{\bar{\xi}.\xi}}\left|\boldsymbol{\xi}\right\rangle\\!\left\langle\boldsymbol{\xi}\right|,$
(C.2a)
$\text{tr}\mathcal{O}=\int\prod_{i}\text{d}\bar{\xi_{i}}\text{d}\xi_{i}e^{-\boldsymbol{\bar{\xi}.\xi}}\left\langle-\boldsymbol{\xi}\right|\mathcal{O}\left|\boldsymbol{\xi}\right\rangle.$
(C.2b)

Therefore, it is possible to calculate the trace of $\rho^{2}$ in terms of
Berezin integrals of Grassmann variables:

$\text{tr}({\rho_{A}^{\beta}}^{2})=\int\footnotesize\prod_{i}\text{d}\bar{\xi_{i}}\text{d}\xi_{i}\prod_{l}\text{d}\bar{\eta}_{l}\text{d}\eta_{l}e^{\small-\boldsymbol{\bar{\xi}.\xi}}e^{-\boldsymbol{\bar{\eta}.\eta}}\left\langle-\boldsymbol{\xi}\right|\rho_{A}^{\beta}\left|\boldsymbol{\eta}\right\rangle\\!\left\langle\boldsymbol{\eta}\right|\rho_{A}^{\beta}\left|\boldsymbol{\xi}\right\rangle.$
(C.3)

We can use the result of (105) to write the integrand in the expression above
as:

$\displaystyle\text{tr}({\rho_{A}^{\beta}}^{2})=$
$\displaystyle{\mathcal{C}^{\beta}}^{2}\\!\int\footnotesize\prod_{i}\\!\text{d}\bar{\xi_{i}}\text{d}\xi_{i}\prod_{l}\\!\text{d}\bar{\eta}_{l}\text{d}\eta_{l}e^{-\boldsymbol{\bar{\xi}.\xi}}e^{-\boldsymbol{\bar{\eta}.\eta}}{\footnotesize\big{[}(\frac{1}{\beta}-\boldsymbol{\mathcal{L}}_{1}.\bar{\boldsymbol{\xi}}+\boldsymbol{\mathcal{L}}_{2}.\boldsymbol{\eta})(\frac{1}{\beta^{*}}-\boldsymbol{\mathcal{L}}_{3}.\bar{\boldsymbol{\xi}}+\boldsymbol{\mathcal{L}}_{4}.\boldsymbol{\eta})-\text{Pf}[\boldsymbol{\mathcal{W}}]\big{]}}$
(C.4) $\displaystyle
e^{\tiny\tfrac{1}{2}\left(\begin{matrix}\bar{\xi}&\\!\eta\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}^{\prime}\end{matrix}\left(\begin{matrix}\bar{\xi}\\\
\eta\end{matrix}\right)}{\footnotesize\big{[}(\frac{1}{\beta}+\boldsymbol{\mathcal{L}}_{1}.\bar{\boldsymbol{\eta}}+\boldsymbol{\mathcal{L}}_{2}.\boldsymbol{\xi})(\frac{1}{\beta^{*}}+\boldsymbol{\mathcal{L}}_{3}.\bar{\boldsymbol{\eta}}+\boldsymbol{\mathcal{L}}_{4}.\boldsymbol{\xi})-\text{Pf}[\boldsymbol{\mathcal{W}}]\big{]}}e^{\tiny\tfrac{1}{2}\left(\begin{matrix}\bar{\eta}&\\!\xi\\\
\end{matrix}\right)\begin{matrix}\boldsymbol{\Omega}\end{matrix}\left(\begin{matrix}\bar{\eta}\\\
\xi\end{matrix}\right)},$

By defining a new Grassmann variable such as:

$\boldsymbol{\theta}=\left(\\!\begin{matrix}\bar{\boldsymbol{\xi}}\\\
\bar{\boldsymbol{\eta}}\\\ \boldsymbol{\eta}\\\
\boldsymbol{\xi}\end{matrix}\\!\right),$ (C.5)

the integral (C.3) can be written in a simpler form as:

$\text{tr}({\rho_{A}^{\beta}}^{2})={\mathcal{C}^{\beta}}^{2}\int\boldsymbol{\mathrm{D}\theta}\mathrm{f}(\boldsymbol{\theta})e^{\small\tfrac{1}{2}\boldsymbol{\theta}^{T}\\!.\boldsymbol{\mathcal{B}}.\boldsymbol{\theta}},$
(C.6)

where in the above
$\boldsymbol{\mathrm{D}\theta}=\prod_{i}\\!\text{d}\bar{\theta_{i}}\text{d}\theta_{i}$,
and also

$\mathrm{f}(\boldsymbol{\theta})=\big{[}(\frac{1}{\beta}+\mathbb{C}.\boldsymbol{\theta}|_{1})(\frac{1}{\beta^{*}}+\mathbb{C}.\boldsymbol{\theta}|_{2})+\text{Pf}[\boldsymbol{\mathcal{W}}]\big{]}\big{[}(\frac{1}{\beta}+\mathbb{C}.\boldsymbol{\theta}|_{3})(\frac{1}{\beta^{*}}+\mathbb{C}.\boldsymbol{\theta}|_{4})+\text{Pf}[\boldsymbol{\mathcal{W}}]\big{]}.$

The new $\boldsymbol{\mathcal{B}}$ and $\boldsymbol{\mathbb{C}}$ matrices are
defined below

$\boldsymbol{\mathcal{B}}=\left(\begin{matrix}\boldsymbol{\Omega}_{{}_{11}}&\boldsymbol{0}&-\boldsymbol{\Omega}_{{}_{12}}&\mathbf{I}\\\
\boldsymbol{0}&\boldsymbol{\Omega}_{{}_{11}}&\mathbf{I}&\boldsymbol{\Omega}_{{}_{12}}\\\
-\boldsymbol{\Omega}_{{}_{21}}&-\mathbf{I}&\boldsymbol{\Omega}_{{}_{22}}&\boldsymbol{0}\\\
-\mathbf{I}&\boldsymbol{\Omega}_{{}_{21}}&\boldsymbol{0}&\boldsymbol{\Omega}_{{}_{21}}\\\
\end{matrix}\right),\qquad\mathbb{C}=\left(\begin{matrix}-\boldsymbol{\mathcal{L}}_{1}&\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{2}&\boldsymbol{0}\\\
-\boldsymbol{\mathcal{L}}_{3}&\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{4}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{1}&\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{2}\\\
\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{3}&\boldsymbol{0}&\boldsymbol{\mathcal{L}}_{4}\\\
\end{matrix}\right).$ (C.7)

The rest of matrices are given by (106). Before proceeding to solve the
integral above, we have to clarify couple of points here. First, the notation
$\mathbb{C}.\boldsymbol{\theta}|_{r}$ stand for the $r^{\text{th}}$ row of the
matrix product $\mathbb{C}.\boldsymbol{\theta}$. Second, the
$\mathrm{f}(\boldsymbol{\theta})$ produce terms with product of odd number of
Grassmann variables such as $\mathbb{C}.\boldsymbol{\theta}|_{r}$ or
$\mathbb{C}.\boldsymbol{\theta}|_{r_{3}}\mathbb{C}.\boldsymbol{\theta}|_{r_{2}}\mathbb{C}.\boldsymbol{\theta}|_{r_{1}}$,
the integration on these terms will be zero automatically. Therefore, The
above integral is straightforward to solve, and the final result of tracing is
given by (133).

In some cases $\boldsymbol{\mathcal{B}}$ does not have an inverse for odd $n$.
However, it is still possible to take the Grassmann integrations.

##  Appendix D Exact diagonalization of A=B=0

In the case where $\mathbf{A}$ and $\mathbf{B}$ are zero then the Hamiltonian
(5) becomes

$\displaystyle H=$
$\displaystyle\sum_{j=1}^{L}\big{[}\alpha_{j}^{0}(c_{0}c_{j}-c_{0}^{\dagger}c_{j})+\alpha_{j}^{L+1}(c_{j}c_{L+1}^{\dagger}+c_{j}c_{L+1})+{\alpha_{j}^{0}}^{*}(c_{j}^{\dagger}c_{0}^{\dagger}-c_{j}^{\dagger}c_{0})+{\alpha_{j}^{L+1}}^{*}(c_{L+1}c_{j}^{\dagger}+c_{L+1}^{\dagger}c_{j}^{\dagger})\big{]}.$
(D.1)

This model can be related to a linear fermionic model such as the model
discussed in [45]

$H=\>\sum_{j=1}^{L}\alpha_{j}c_{j}+\alpha_{j}^{*}c_{j}^{\dagger}.$ (D.2)

In fact, with a specific projection, one can find the eigenstates of (D.2) in
the Hilbert space of (D.1). We can write the new form of the $\mathbf{M}$
matrix (9)

$\mathbf{M}=\left(\begin{array}[]{ccccccc|ccccccc}0&-{\alpha_{1}^{0}}&-{\alpha_{2}^{0}}&\cdots&-{\alpha_{L-1}^{0}}&-{\alpha_{L}^{0}}&0&0&-{\alpha_{1}^{0}}^{*}&-{\alpha_{2}^{0}}^{*}&\cdots&-{\alpha_{L-1}^{0}}^{*}&-{\alpha_{L}^{0}}^{*}&0\\\
-{\alpha_{1}^{0}}^{*}&{0}&{0}&\cdots&{0}&{0}&-{\alpha_{1}^{L+1}}^{*}&{\alpha_{1}^{0}}^{*}&0&0&\cdots&0&0&-{\alpha_{1}^{L+1}}^{*}\\\
-{\alpha_{2}^{0}}^{*}&0&{}\hfil&{}\hfil&{}\hfil&0&-{\alpha_{2}^{L+1}}^{*}&{\alpha_{2}^{0}}^{*}&0&{}\hfil&{}\hfil&{}\hfil&0&-{\alpha_{2}^{L+1}}^{*}\\\
\vdots&\vdots&{}\hfil&{}\hfil&{}\hfil&{}\hfil&\vdots&\vdots&\vdots&{}\hfil&&{}\hfil&\vdots&\vdots\\\
-{\alpha_{L-1}^{0}}^{*}&0&&&&0&-{\alpha_{L-1}^{L+1}}^{*}&{\alpha_{L-1}^{0}}^{*}&0&&&&0&-{\alpha_{L-1}^{L+1}}^{*}\\\
-{\alpha_{L}^{0}}^{*}&0&0&\cdots&0&0&-{\alpha_{L}^{L+1}}^{*}&{\alpha_{L}^{0}}^{*}&0&0&\cdots&0&0&-{\alpha_{L}^{L+1}}^{*}\\\
0&-{\alpha_{1}^{L+1}}&-{\alpha_{2}^{L+1}}&\cdots&-{\alpha_{L-1}^{L+1}}&-{\alpha_{L}^{L+1}}&0&0&{\alpha_{1}^{L+1}}^{*}&{\alpha_{2}^{L+1}}^{*}&\cdots&{\alpha_{L-1}^{L+1}}^{*}&{\alpha_{L}^{L+1}}^{*}&0\\\
\hline\cr
0&{\alpha_{1}^{0}}&{\alpha_{2}^{0}}&\cdots&{\alpha_{L-1}^{0}}&{\alpha_{L}^{0}}&0&0&{\alpha_{1}^{0}}^{*}&{\alpha_{2}^{0}}^{*}&\cdots&{\alpha_{L-1}^{0}}^{*}&{\alpha_{L}^{0}}^{*}&0\\\
-{\alpha_{1}^{0}}&0&0&\cdots&0&0&{\alpha_{1}^{L+1}}&{\alpha_{1}^{0}}&0&0&\cdots&0&0&{\alpha_{1}^{L+1}}\\\
-{\alpha_{2}^{0}}&0&{}\hfil&{}\hfil&{}\hfil&0&{\alpha_{2}^{L+1}}&{\alpha_{2}^{0}}&0&{}\hfil&{}\hfil&{}\hfil&0&{\alpha_{2}^{L+1}}\\\
\vdots&\vdots&{}\hfil&{}\hfil&{}\hfil&\vdots&\vdots&\vdots&\vdots&{}\hfil&&{}\hfil&\vdots&\vdots\\\
-{\alpha_{L-1}^{0}}&0&{}\hfil&{}\hfil&{}\hfil&0&{\alpha_{L-1}^{L+1}}&{\alpha_{L-1}^{0}}&0&{}\hfil&{}\hfil&{}\hfil&0&{\alpha_{L-1}^{L+1}}\\\
-{\alpha_{L}^{0}}&0&0&\cdots&0&0&{\alpha_{L}^{L+1}}&{\alpha_{L}^{0}}&0&0&\cdots&0&0&{\alpha_{L}^{L+1}}\\\
0&-{\alpha_{1}^{L+1}}&-{\alpha_{2}^{L+1}}&\cdots&-{\alpha_{L-1}^{L+1}}&-{\alpha_{L}^{L+1}}&0&0&{\alpha_{1}^{L+1}}^{*}&{\alpha_{2}^{L+1}}^{*}&\cdots&{\alpha_{L-1}^{L+1}}^{*}&{\alpha_{L}^{L+1}}^{*}&0\end{array}\right).$
(D.3)

It can be shown that $\mathbf{\Lambda}$ has only two nonzero diagonal
elements, $\lambda_{1}$ and $\lambda_{2}$ and the rest of eigenvalues are
zero. For nonzero eigenvalues, we have

$\lambda_{1},\lambda_{2}=\sqrt{2}\sqrt{\sum_{j=1}^{L}(|\alpha_{j}^{0}|^{2}+|\alpha_{j}^{L+1}|^{2})\pm\sqrt{\Big{|}\sum_{j=1}^{L}(\alpha_{j}^{0}+\alpha_{j}^{L+1})(\alpha_{j}^{0*}-\alpha_{j}^{L+1*})\Big{|}^{2}}}.$
(D.4)

Then, $\mathbf{M}$ would have only four nonzero eigenvalues and the rest are
zero independent of size. The next step is to find the eigenvectors to
construct the $\mathbf{U}$ matrix. In (12), we have already introduced two
eigenvectors of zero mode. Using the orthogonality condition, and the
condition that if $\left|u\right\rangle$ is an eigenvector with eigenvalue
$\lambda$ then $\mathbf{J}\left|u\right\rangle$ is also an eigenvector with
eigenvalue $-\lambda$, therefore, we can find and present only half of the
eigenvectors. For simplicity, we assume
$\alpha_{j}^{0}=\alpha_{j}^{L+1}=\alpha_{j}$, then eigenvectors of positive
modes. we have:

$\left|u^{+}_{1}\right\rangle=\left(\begin{array}[]{c}u\\\
v\end{array}\right);\qquad\quad\begin{array}[]{l}u_{L+1}=u_{0}=-v_{0}=v_{L+1}=1\\\
u_{i}=-\frac{2\alpha_{i}^{*}}{\sqrt{\sum_{j=1}^{L}|\alpha_{j}|^{2}}}\\\
\text{else}=0\end{array}$ (D.5)
$\left|u^{+}_{2}\right\rangle=\left(\begin{array}[]{c}u\\\
v\end{array}\right);\qquad\quad\begin{array}[]{l}u_{L+1}=-u_{0}=v_{0}=v_{L+1}=\frac{\alpha_{L}}{|\alpha_{L}|}\\\
v_{i}=\frac{2\alpha_{i}}{\sqrt{\sum_{j=1}^{L}|\alpha_{j}|^{2}}}\frac{\alpha_{L}}{|\alpha_{L}|}\\\
\text{else}=0\end{array}$ (D.6)

For eigenvectors of zero modes, we have:

$\left|u^{0}_{1}\right\rangle=\left(\begin{array}[]{c}u\\\
v\end{array}\right);\qquad\quad\begin{array}[]{l}u_{L+1}=u_{0}=v_{0}=-v_{L+1}=1\\\
\text{else}=0\end{array}$ (D.7)
$\qquad\;\;\left|u^{0}_{k}\right\rangle=\left(\begin{array}[]{c}u\\\
v\end{array}\right);\qquad\quad\begin{array}[]{l}u_{i}=-\frac{\alpha_{k}\alpha_{i}^{*}}{\sum_{j=1}^{k-1}|\alpha_{j}|^{2}};\quad
i<k\neq 0\\\ u_{k}=1\\\ \text{else}=0,\end{array}$ (D.8)

where $k=2,3,\cdots,L$. With this expressions for eigenvectors, we can
construct the unitary matrices $\mathbf{U}$. Having the exact $\mathbf{U}$, we
can calculate correlation matrices (see section 4). putting these eigenvectors
together, one can construct the $\mathbf{U}$.

Having the $\mathbf{U}$ matrix, we can directly calculate the $\delta_{+}$
introduced in (42). The calculation of $\delta_{+}$ shows that

$\delta_{+}=\begin{cases}+1&\text{if}\;\>L\;\>\text{even}\\\
-1&\text{if}\;\>L\;\>\text{odd}\\\ \end{cases}$ (D.9)

###  D.1 Vacuum state in configuration basis

Another method to calculate the entanglement of subsystem in a particular
state is to use a more direct way. We could use the equation (28) to calculate
the entanglement. Using the results of subsection 3.4 for
$\mathbf{A},\mathbf{B}=0$ we can find the form of $\mathbf{R}$ matrix as:.

$\mathbf{R}=\left(\begin{array}[]{c|c|c}0&\frac{\alpha_{j}}{\sqrt{|\alpha_{1}|^{2}+\cdots+|\alpha_{L}|^{2}}}&1\\\
\hline\cr\frac{-\alpha_{i}}{\sqrt{|\alpha_{1}|^{2}+\cdots+|\alpha_{L}|^{2}}}&0&\frac{\alpha_{i}}{\sqrt{|\alpha_{1}|^{2}+\cdots+|\alpha_{L}|^{2}}}\\\
\hline\cr-1&\frac{-\alpha_{j}}{\sqrt{|\alpha_{1}|^{2}+\cdots+|\alpha_{L}|^{2}}}&0\\\
\end{array}\right)$ (D.10)

where,

$\left|0\right\rangle_{\eta}=\frac{1}{\sqrt{\mathrm{det}[\mathbf{I}+\mathbf{R}^{\dagger}\mathbf{R}]}}\scalebox{1.5}{$e$}^{\tfrac{1}{2}\sum_{i,j}R_{ij}c^{\dagger}_{i}c_{j}^{\dagger}}\left|0\right\rangle_{c}.$

Numerical investigations suggest that for even sizes we get exactly the same
expression, and for odd sizes we get the minus of the vacuum
($\alpha_{i}\in\mathbb{R}$).

###  D.2 Correlations

For the special case of $\mathbf{A},\mathbf{B}=0$, the correlation matrices
look pretty much simple in terms of the parameters $\alpha_{i}$’s. We are
going to demonstrate the matrix form of correlations for general
$\alpha_{i}\in\mathbb{C}$. For start, in the vacuum state, correlations looks
like below. (In the following expressions indexes $n$ and $m$ are the number
of rows and columns respectively.)

$\mathbf{C}^{0}=\left(\begin{array}[]{c|c|c}\frac{1}{2}&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&-\frac{1}{4}\\\
\hline\cr\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{\alpha_{n}^{*}\alpha_{m}}{2(\sum_{k}|\alpha_{k}|^{2})}&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr-\frac{1}{4}&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{1}{2}\\\
\end{array}\right),\quad\mathbf{F}^{0}=\left(\begin{array}[]{c|c|c}0&\frac{\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{1}{4}\\\
\hline\cr\frac{-\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr-\frac{1}{4}&\frac{-\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\end{array}\right)$ $\mathbf{K}^{0}=\left(\begin{array}[]{c|c|c}1&0&\qquad
0\\\ \hline\cr
0&\;\delta_{nm}-\frac{m-n}{|m-n|}\frac{i\Im[\alpha_{n}^{*}\alpha_{m}]}{\sum_{k}|\alpha_{k}|^{2}}&\frac{-i\Im[\alpha_{n}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr 0&\;\frac{i\Im[\alpha_{m}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&1\\\
\end{array}\right),\quad\mathbf{\bar{K}}^{0}=\left(\begin{array}[]{c|c|c}1&\frac{i\Im[\alpha_{m}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\hline\cr\;\frac{-i\Im[\alpha_{n}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\;\delta_{nm}-\frac{m-n}{|m-n|}\frac{i\Im[\alpha_{n}^{*}\alpha_{m}]}{\sum_{k}|\alpha_{k}|^{2}}&0\\\
\hline\cr 0&0&1\\\ \end{array}\right)$
$\mathbf{G}^{0}=\left(\begin{array}[]{c|c|c}0&\frac{\Re[\alpha_{m}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\qquad
0\\\ \hline\cr
0&\;-\delta_{nm}+\frac{\Re[\alpha_{n}^{*}\alpha_{m}]}{\sum_{k}|\alpha_{k}|^{2}}&\;\frac{\Re[\alpha_{n}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr\qquad-1&0&0\\\ \end{array}\right)$

For the case where the state is ZME state ($\left|\emptyset\right\rangle$)
then the correlation matrices look like below.

$\mathbf{C}^{\emptyset}=\left(\begin{array}[]{c|c|c}\frac{1}{2}&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{1}{4}\\\
\hline\cr\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{\alpha_{n}^{*}\alpha_{m}}{2(\sum_{k}|\alpha_{k}|^{2})}&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr\frac{1}{4}&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{1}{2}\\\
\end{array}\right),\quad\mathbf{F}^{\emptyset}=\left(\begin{array}[]{c|c|c}0&\frac{\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&-\frac{1}{4}\\\
\hline\cr\frac{-\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr\frac{1}{4}&\frac{-\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\end{array}\right)$
$\mathbf{G}^{\emptyset}=\left(\begin{array}[]{c|c|c}0&\frac{\Re[\alpha_{m}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\qquad
0\\\ \hline\cr
0&\;-\delta_{nm}+\frac{\Re[\alpha_{n}^{*}\alpha_{m}]}{\sum_{k}|\alpha_{k}|^{2}}&\;\frac{\Re[\alpha_{n}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr\;1&0&0\\\
\end{array}\right),\qquad\mathbf{K}^{\emptyset}=\mathbf{K}^{0},\qquad\bar{\mathbf{K}}^{\emptyset}=\bar{\mathbf{K}}^{0}.$

Finally, for this case when the state is $\left|G_{\pm}\right\rangle$ then the
correlations look like below.

$\mathbf{C}^{\pm}=\left(\begin{array}[]{c|c|c}\frac{1}{2}&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\hline\cr\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{\alpha_{n}^{*}\alpha_{m}}{2(\sum_{k}|\alpha_{k}|^{2})}&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr 0&\frac{\alpha_{m}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\frac{1}{2}\\\
\end{array}\right),\quad\mathbf{F}^{\pm}=\left(\begin{array}[]{c|c|c}0&\frac{\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\hline\cr\frac{-\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0&\frac{\alpha_{n}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr 0&\frac{-\alpha_{m}^{*}}{4\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&0\\\
\end{array}\right),$
$\mathbf{G}^{\pm}=\left(\begin{array}[]{c|c|c}0&\frac{\Re[\alpha_{m}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}&\qquad
0\\\ \hline\cr
0&\;\>-\delta_{nm}+\frac{\Re[\alpha_{n}^{*}\alpha_{m}]}{\sum_{k}|\alpha_{k}|^{2}}&\;\frac{\Re[\alpha_{n}]}{\sqrt{\sum_{k}|\alpha_{k}|^{2}}}\\\
\hline\cr 0&0&0\\\
\end{array}\right),\qquad\mathbf{K}^{\pm}=\mathbf{K}^{0},\qquad\bar{\mathbf{K}}^{\pm}=\bar{\mathbf{K}}^{0}.$

##  Appendix E M matrix for modified $XY$ chain

In this appendix, the exact forms of the matrices A and $\mathbf{B}$ for
modified $XY$ spin chain is demonstrated. For $|\vec{b}_{1}|=1$ and
$|\vec{b}_{L}|=1$, they can be written as:

$\mathbf{A}=\left(\begin{array}[]{cccccccc}0&-\frac{\sin\theta_{1}e^{i\varphi_{1}}}{2}&0&...\\\
-\frac{\sin\theta_{1}e^{-i\varphi_{1}}}{2}&-h+\cos\theta_{1}&-\frac{J}{2}&0&...\\\
0&-\frac{J}{2}&-h&-\frac{J}{2}&0&...\\\
0&0&-\frac{J}{2}&-h&-\frac{J}{2}&0&...\\\ .&.&.\\\ .&.&.\\\
&&&&.&.&-h+\cos\theta_{L}&-\frac{\sin\theta_{L}e^{-i\varphi_{L}}}{2}\\\
&&&&\cdots&0&-\frac{\sin\theta_{L}e^{i\varphi_{L}}}{2}&0\end{array}\right),$
(E.1)
$\mathbf{B}=\begin{pmatrix}0&-\frac{\sin\theta_{1}e^{-i\varphi_{1}}}{2}&0&\cdots\\\
\frac{\sin\theta_{1}e^{-i\varphi_{1}}}{2}&0&-\frac{J\gamma}{2}&0&\cdots\\\
0&\frac{J\gamma}{2}&0&-\frac{\gamma J}{2}&0&...\\\
\vdots&0&\frac{J\gamma}{2}&0&-\frac{\gamma J}{2}&0&...\\\ &.&.\\\ &.&.\\\
&.&.&.&.&.&\frac{J\gamma}{2}&0&-\frac{\sin\theta_{L}e^{-i\varphi_{L}}}{2}\\\
&.&.&.&.&.&0&\frac{\sin\theta_{L}e^{-i\varphi_{L}}}{2}&0\end{pmatrix}.$ (E.2)

With the help of (8) and (9), one can construct the actual $\mathbf{M}$
matrix.

## References

  * [1] E. Lieb, T. Schultz and D. Mattis,  _Two soluble models of an antiferromagnetic chain_ ,  Annals of Physics 16(3), 407 (1961),  https://doi.org/10.1016/0003-4916(61)90115-4.
  * [2] M.-C. Chung and I. Peschel,  _Density-matrix spectra of solvable fermionic systems_ ,  Phys. Rev. B 64, 064412 (2001),  10.1103/PhysRevB.64.064412.
  * [3] I. Peschel,  _On the reduced density matrix for a chain of free electrons_ ,  Journal of Statistical Mechanics: Theory and Experiment 2004(06), P06004 (2004),  10.1088/1742-5468/2004/06/p06004.
  * [4] G. Vidal, J. I. Latorre, E. Rico and A. Kitaev,  _Entanglement in quantum critical phenomena_ ,  Phys. Rev. Lett. 90, 227902 (2003),  10.1103/PhysRevLett.90.227902.
  * [5] B.-Q. Jin and V. E. Korepin,  _Quantum spin chain, toeplitz determinants and the fisher—hartwig conjecture_ ,  Journal of statistical physics 116(1), 79 (2004).
  * [6] J. P. Keating and F. Mezzadri,  _Random matrix theory and entanglement in quantum spin chains_ ,  Communications in Mathematical Physics 252(1), 543 (2004),  10.1007/s00220-004-1188-2.
  * [7] I. Peschel and V. Eisler,  _Reduced density matrices and entanglement entropy in free lattice models_ ,  Journal of Physics A: Mathematical and Theoretical 42(50), 504003 (2009),  10.1088/1751-8113/42/50/504003.
  * [8] F. Franchini and A. G. Abanov,  _Asymptotics of toeplitz determinants and the emptiness formation probability for the XY spin chain_ ,  Journal of Physics A: Mathematical and General 38(23), 5069 (2005),  10.1088/0305-4470/38/23/002.
  * [9] K. Najafi and M. A. Rajabpour,  _Formation probabilities and shannon information and their time evolution after quantum quench in the transverse-field xy chain_ ,  Phys. Rev. B 93, 125139 (2016),  10.1103/PhysRevB.93.125139.
  * [10] M. N. Najafi and M. A. Rajabpour,  _Formation probabilities and statistics of observables as defect problems in free fermions and quantum spin chains_ ,  Phys. Rev. B 101, 165415 (2020),  10.1103/PhysRevB.101.165415.
  * [11] V. Alba, M. Fagotti and P. Calabrese,  _Entanglement entropy of excited states_ ,  Journal of Statistical Mechanics: Theory and Experiment 2009(10), P10020 (2009),  10.1088/1742-5468/2009/10/p10020.
  * [12] F. C. Alcaraz, M. I. n. Berganza and G. Sierra,  _Entanglement of low-energy excitations in conformal field theory_ ,  Phys. Rev. Lett. 106, 201601 (2011),  10.1103/PhysRevLett.106.201601.
  * [13] M. I. Berganza, F. C. Alcaraz and G. Sierra,  _Entanglement of excited states in critical spin chains_ ,  Journal of Statistical Mechanics: Theory and Experiment 2012(01), P01016 (2012),  10.1088/1742-5468/2012/01/p01016.
  * [14] F. Ares, J. G. Esteve, F. Falceto and E. Sánchez-Burillo,  _Excited state entanglement in homogeneous fermionic chains_ ,  Journal of Physics A: Mathematical and Theoretical 47(24), 245301 (2014),  10.1088/1751-8113/47/24/245301.
  * [15] L. Vidmar, L. Hackl, E. Bianchi and M. Rigol,  _Entanglement entropy of eigenstates of quadratic fermionic hamiltonians_ ,  Phys. Rev. Lett. 119, 020601 (2017),  10.1103/PhysRevLett.119.020601.
  * [16] L. Hackl, L. Vidmar, M. Rigol and E. Bianchi,  _Average eigenstate entanglement entropy of the xy chain in a transverse field and its universality for translationally invariant quadratic fermionic models_ ,  Phys. Rev. B 99, 075123 (2019),  10.1103/PhysRevB.99.075123.
  * [17] A. Jafarizadeh and M. A. Rajabpour,  _Bipartite entanglement entropy of the excited states of free fermions and harmonic oscillators_ ,  Phys. Rev. B 100, 165135 (2019),  10.1103/PhysRevB.100.165135.
  * [18] J. Zhang and M. A. Rajabpour,  _Universal Rényi Entropy of Quasiparticle Excitations_ ,  arXiv e-prints arXiv:2010.13973 (2020),  2010.13973.
  * [19] J. Zhang and M. A. Rajabpour,  _Universal Rényi entropy in quasiparticle excited states of quantum chains_ ,  arXiv e-prints arXiv:2010.16348 (2020),  2010.16348.
  * [20] N. Laflorencie, E. S. Sørensen, M.-S. Chang and I. Affleck,  _Boundary effects in the critical scaling of entanglement entropy in 1d systems_ ,  Phys. Rev. Lett. 96, 100603 (2006),  10.1103/PhysRevLett.96.100603.
  * [21] H.-Q. Zhou, T. Barthel, J. O. Fjærestad and U. Schollwöck,  _Entanglement and boundary critical phenomena_ ,  Phys. Rev. A 74, 050305 (2006),  10.1103/PhysRevA.74.050305.
  * [22] O. Legeza, J. Sólyom, L. Tincani and R. M. Noack,  _Entropic analysis of quantum phase transitions from uniform to spatially inhomogeneous phases_ ,  Phys. Rev. Lett. 99, 087203 (2007),  10.1103/PhysRevLett.99.087203.
  * [23] E. Szirmai, O. Legeza and J. Sólyom,  _Spatially nonuniform phases in the one-dimensional $\mathrm{SU}(n)$ hubbard model for commensurate fillings_,  Phys. Rev. B 77, 045106 (2008),  10.1103/PhysRevB.77.045106.
  * [24] O. A. Castro-Alvaredo and B. Doyon,  _Bi-partite entanglement entropy in massive qft with a boundary: the ising model_ ,  Journal of Statistical Physics 134 (2009),  10.1007/s10955-008-9664-2.
  * [25] I. Affleck, N. Laflorencie and E. S. Sørensen,  _Entanglement entropy in quantum impurity systems and systems with boundaries_ ,  Journal of Physics A: Mathematical and Theoretical 42(50), 504009 (2009),  10.1088/1751-8113/42/50/504009.
  * [26] L. Taddia, J. C. Xavier, F. C. Alcaraz and G. Sierra,  _Entanglement entropies in conformal systems with boundaries_ ,  Phys. Rev. B 88, 075112 (2013),  10.1103/PhysRevB.88.075112.
  * [27] M. Fagotti and P. Calabrese,  _Universal parity effects in the entanglement entropy of xx chains with open boundary conditions_ ,  Journal of Statistical Mechanics: Theory and Experiment 2011(01), P01017 (2011),  10.1088/1742-5468/2011/01/p01017.
  * [28] F. Iglói and I. Peschel,  _On reduced density matrices for disjoint subsystems_ ,  EPL (Europhysics Letters) 89(4), 40001 (2010),  10.1209/0295-5075/89/40001.
  * [29] S. B. Bravyi and A. Y. Kitaev,  _Fermionic quantum computation_ ,  Annals of Physics 298(1), 210 (2002),  https://doi.org/10.1006/aphy.2002.6254.
  * [30] M.-C. Bañuls, J. I. Cirac and M. M. Wolf,  _Entanglement in fermionic systems_ ,  Phys. Rev. A 76, 022311 (2007),  10.1103/PhysRevA.76.022311.
  * [31] H. Moriya,  _On separable states for composite systems of distinguishable fermions_ ,  Journal of Physics A: Mathematical and General 39(14), 3753 (2006).
  * [32] H. J. de Vega and A. G. Ruiz,  _Boundary k-matrices for the six vertex and the n(2n-1)an-1vertex models_ ,  Journal of Physics A: Mathematical and General 26(12), L519 (1993),  10.1088/0305-4470/26/12/007.
  * [33] R. I. Nepomechie,  _Solving the open xxz spin chain with nondiagonal boundary terms at roots of unity_ ,  Nuclear Physics B 622(3), 615 (2002),  https://doi.org/10.1016/S0550-3213(01)00585-5.
  * [34] R. I. Nepomechie,  _Bethe ansatz solution of the open XXZ chain with nondiagonal boundary terms_ ,  Journal of Physics A: Mathematical and General 37(2), 433 (2003),  10.1088/0305-4470/37/2/012.
  * [35] J. Cao, H.-Q. Lin, K.-J. Shi and Y. Wang,  _Exact solution of xxz spin chain with unparallel boundary fields_ ,  Nuclear Physics B 663(3), 487 (2003),  https://doi.org/10.1016/S0550-3213(03)00372-9.
  * [36] J. Cao, W.-L. Yang, K. Shi and Y. Wang,  _Off-diagonal bethe ansatz solution of the xxx spin chain with arbitrary boundary conditions_ ,  Nuclear Physics B 875(1), 152 (2013),  https://doi.org/10.1016/j.nuclphysb.2013.06.022.
  * [37] J. Cao, W.-L. Yang, K. Shi and Y. Wang,  _Off-diagonal bethe ansatz solutions of the anisotropic spin-12 chains with arbitrary boundary fields_ ,  Nuclear Physics B 877(1), 152 (2013),  https://doi.org/10.1016/j.nuclphysb.2013.10.001.
  * [38] G. Niccoli,  _Antiperiodic spin-1/2 xxz quantum chains by separation of variables: Complete spectrum and form factors_ ,  Nuclear Physics B 870(2), 397 (2013),  https://doi.org/10.1016/j.nuclphysb.2013.01.017.
  * [39] S. Belliard, N. Crampé _et al._ ,  _Heisenberg xxx model with general boundaries: eigenvectors from algebraic bethe ansatz_ ,  SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 9, 072 (2013).
  * [40] R. I. Nepomechie and C. Wang,  _Boundary energy of the open XXX chain with a non-diagonal boundary term_ ,  Journal of Physics A: Mathematical and Theoretical 47(3), 032001 (2013),  10.1088/1751-8113/47/3/032001.
  * [41] S. Faldella, N. Kitanine and G. Niccoli,  _The complete spectrum and scalar products for the open spin-1/2 XXZ quantum chains with non-diagonal boundary terms_ ,  Journal of Statistical Mechanics: Theory and Experiment 2014(1), P01011 (2014),  10.1088/1742-5468/2014/01/p01011.
  * [42] N. Kitanine, J. M. Maillet and G. Niccoli,  _Open spin chains with generic integrable boundaries: Baxter equation and bethe ansatz completeness from separation of variables_ ,  Journal of Statistical Mechanics: Theory and Experiment 2014(5), P05015 (2014),  10.1088/1742-5468/2014/05/p05015.
  * [43] Y.-Y. Li, J. Cao, W.-L. Yang, K. Shi and Y. Wang,  _Thermodynamic limit and surface energy of the xxz spin chain with arbitrary boundary fields_ ,  Nuclear Physics B 884, 17 (2014),  https://doi.org/10.1016/j.nuclphysb.2014.04.010.
  * [44] B. Pozsgay and O. Rákos,  _Exact boundary free energy of the open XXZ chain with arbitrary boundary conditions_ ,  Journal of Statistical Mechanics: Theory and Experiment 2018(11), 113102 (2018),  10.1088/1742-5468/aae5a5.
  * [45] J. H. P. Colpa,  _Diagonalisation of the quadratic fermion hamiltonian with a linear part_ ,  Journal of Physics A: Mathematical and General 12(4), 469 (1979),  10.1088/0305-4470/12/4/008.
  * [46] R. Bariev and I. Peschel,  _Non-universal critical behaviour in a two-dimensional ising model with a field_ ,  Physics Letters A 153(4), 166 (1991),  https://doi.org/10.1016/0375-9601(91)90786-8.
  * [47] U. Bilstein and B. Wehefritz,  _The -model with boundaries: Part i. diagonalization of the finite chain_ ,  Journal of Physics A: Mathematical and General 32(2), 191 (1999),  10.1088/0305-4470/32/2/001.
  * [48] M. Campostrini, A. Pelissetto and E. Vicari,  _Quantum ising chains with boundary fields_ ,  Journal of Statistical Mechanics: Theory and Experiment 2015(11), P11015 (2015),  10.1088/1742-5468/2015/11/p11015.
  * [49] J. C. Xavier and M. A. Rajabpour,  _Entanglement and boundary entropy in quantum spin chains with arbitrary direction of the boundary magnetic fields_ ,  Phys. Rev. B 101, 235127 (2020),  10.1103/PhysRevB.101.235127.
  * [50] G. C. Wick,  _The evaluation of the collision matrix_ ,  Phys. Rev. 80, 268 (1950),  10.1103/PhysRev.80.268.
  * [51] G. Vidal, J. I. Latorre, E. Rico and A. Kitaev,  _Entanglement in quantum critical phenomena_ ,  Phys. Rev. Lett. 90, 227902 (2003),  10.1103/PhysRevLett.90.227902.
  * [52] J. I. Latorre, E. Rico and G. Vidal,  _Ground state entanglement in quantum spin chains_ ,  arXiv preprint quant-ph/0304098 (2003).
  * [53] I. Peschel,  _Calculation of reduced density matrices from correlation functions_ ,  Journal of Physics A: Mathematical and General 36(14), L205 (2003),  10.1088/0305-4470/36/14/101.
  * [54] J. L. van Hemmen,  _A note on the diagonalization of quadratic boson and fermion hamiltonians_ ,  Zeitschrift für Physik B Condensed Matter 38(3), 271 (1980),  10.1007/BF01315667.
  * [55] M. Cozzini, P. Giorda and P. Zanardi,  _Quantum phase transitions and quantum fidelity in free fermion graphs_ ,  Phys. Rev. B 75, 014439 (2007),  10.1103/PhysRevB.75.014439.
  * [56] S. Caracciolo, A. D. Sokal and A. Sportiello,  _Algebraic/combinatorial proofs of cayley-type identities for derivatives of determinants and pfaffians_ ,  Advances in Applied Mathematics 50(4), 474 (2013),  https://doi.org/10.1016/j.aam.2012.12.001.
  * [57] R. Balian and E. Brezin,  _Nonunitary bogoliubov transformations and extension of wick’s theorem_ ,  Il Nuovo Cimento B (1965-1970) 64(1), 37 (1969),  10.1007/BF02710281.


# Entropic Corrections to Friedmann Equations

Ahmad Sheykhi 111<EMAIL_ADDRESS>Department of Physics, Shahid Bahonar
University, P.O. Box 76175, Kerman, Iran
Research Institute for Astronomy and Astrophysics of Maragha (RIAAM), Maragha,
Iran

###### Abstract

Recently, Verlinde discussed that gravity can be understood as an entropic
force caused by changes in the information associated with the positions of
material bodies. In the Verlinde’s argument, the area law of the black hole
entropy plays a crucial role. However, the entropy-area relation can be
modified from the inclusion of quantum effects, motivated from the loop
quantum gravity. In this note, by employing this modified entropy-area
relation, we derive corrections to Newton’s law of gravitation as well as
modified Friedman equations by adopting the viewpoint that gravity can be
emerged as an entropic force. Our study further supports the universality of
the log correction and provides a strong consistency check on Verlinde’s
model.
PACS: 04.20.Cv, 04.50.-h, 04.70.Dy.

It was first pointed out by Jacobson Jac  that the hyperbolic second order
partial differential Einstein equation has a predisposition to the first law
of thermodynamics. This profound connection between the first law of
thermodynamics and the gravitational field equations has been extensively
observed in various gravity theories Elin ; Pad . Recently the study on the
connection between thermodynamics and gravity has been generalized to the
cosmological situations Cai2 ; Cai3 ; CaiKim ; Wang ; Cai33 ; Shey0 , where it
was shown that the differential form of the Friedmann equation on the apparent
horizon in the FRW universe can be rewritten in the form of the first law of
thermodynamics. The extension of this connection has also been carried out in
the braneworld cosmology Shey1 ; Shey2 . The deep connection between the
gravitational equation describing the gravity in the bulk and the first law of
thermodynamics on the apparent horizon reflects some deep ideas of holography.
Although Jacobson’s derivation is logically clear and theoretically sound, the
statistical mechanical origin of the thermodynamic nature of general
relativity remains obscure.

Recently, Verlinde Ver  has put forward an idea similar in spirit to
Jacobson’s thermodynamic derivation of the Einstein equations, where it is
argued that Newton’s law of gravitation can be understood as an entropic force
caused by information changes when a material body moves away from the
holographic screen. Quantitatively, when a test particle or excitation moves
apart from the holographic screen, the magnitude of the entropic force on this
body has the form

$F\triangle x=T\triangle S,$ (1)

where $\triangle x$ is the displacement of the particle from the holographic
screen, while $T$ and $\triangle S$ are the temperature and the entropy change
on the screen, respectively. Verlinde’s derivation of Newton’s law of
gravitation at the very least offers a strong analogy with a well understood
statistical mechanism. Therefore, this derivation opens a new window to
understand gravity from the first principles. The study on the entropic force
has raised a lot of enthusiasm recently (see Cai4 ; Other ; Modesto ; Yi  and
references therein).

It is important to note that in Verlinde discussion, the black hole entropy
$S$ plays a crucial role in the derivation of Newton’s law of gravitation.
Indeed, the derivation of Newton’s law of gravity depends on the entropy-area
relationship $S=A/4\ell_{p}^{2}$ of black holes in Einstein s gravity, where
$A=4\pi R^{2}$ represents the area of the horizon and
$\ell_{p}^{2}=G\hbar/c^{3}$ is the Planck length. However, this definition can
be modified from the inclusion of quantum effects, motivated from the loop
quantum gravity (LQG). The quantum corrections provided to the entropy-area
relationship leads to the curvature correction in the Einstein-Hilbert action
and vice versa Zhu ; Suj . The corrected entropy takes the form Zhang

$S=\frac{A}{4\ell_{p}^{2}}-\beta\ln{\frac{A}{4\ell_{p}^{2}}}+\gamma\frac{\ell_{p}^{2}}{A}+\mathrm{const},$
(2)

where $\beta$ and $\gamma$ are dimensionless constants of order unity. The
exact values of these constants are not yet determined and still an open issue
in loop quantum cosmology. These corrections arise in the black hole entropy
in LQG due to thermal equilibrium fluctuations and quantum fluctuations
Rovelli . Taking the corrected entropy-area relation (2) into account, we will
derive the corrections to the Newton’s law of gravitation as well as the
modified Friedman equations. We rewrite Eq. (2) in the following form

$S=\frac{A}{4\ell_{p}^{2}}+{s}(A),$ (3)

where $s(A)$ stands for the correction terms in the entropy expression. We
adopt the viewpoint of Ver . Suppose we have two masses one a test mass and
the other considered as the source with respective masses $m$ and $M$.
Centered around the source mass $M$, is a spherically symmetric surface
$\mathcal{S}$ which will be defined with certain properties that will be made
explicit later. To derive the entropic law, the surface $\mathcal{S}$ is
between the test mass and the source mass, but the test mass is assumed to be
very close to the surface as compared to its reduced Compton wavelength
$\lambda_{m}=\frac{\hbar}{mc}$. When a test mass $m$ is a distance $\triangle
x=\eta\lambda_{m}$ away from the surface $\mathcal{S}$, the entropy of the
surface changes by one fundamental unit $\triangle S$ fixed by the discrete
spectrum of the area of the surface via the relation

$\triangle S=\frac{\partial S}{\partial A}\triangle
A=\left(\frac{1}{4\ell_{p}^{2}}+\frac{\partial{s}(A)}{\partial
A}\right)\triangle A.$ (4)

The energy of the surface $\mathcal{S}$ is identified with the relativistic
rest mass of the source mass:

$E=Mc^{2}.$ (5)

On the surface $\mathcal{S}$, there live a set of “bytes” of information that
scale proportional to the area of the surface so that

$A=QN,$ (6)

where $N$ represents the number of bytes and $Q$ is a fundamental constant
which should be specified later. Assuming the temperature on the surface is
$T$, and then according to the equipartition law of energy Pad1 , the total
energy on the surface is

$E=\frac{1}{2}Nk_{B}T.$ (7)

Finally, we assume that the force on the particle follows from the generic
form of the entropic force governed by the thermodynamic equation of state

$F=T\frac{\triangle S}{\triangle x},$ (8)

where $\triangle S$ is one fundamental unit of entropy when $|\triangle
x|=\eta\lambda_{m}$, and the entropy gradient points radially from the outside
of the surface to inside. Note that $N$ is the number of bytes and thus
$\triangle N=1$, hence from (6) we have $\triangle A=Q$. Now, we are in a
position to derive the entropic-corrected Newton’s law of gravity. Combining
Eqs. (4)- (8), we easily obtain

$F=-\frac{Mm}{R^{2}}\left(\frac{Qc^{3}}{2\pi
k_{B}\hbar\eta}\right)\left[\frac{1}{4\ell_{p}^{2}}+\frac{\partial{s}}{\partial
A}\right]_{A=4\pi R^{2}},$ (9)

This is nothing but the Newton’s law of gravitation to the first order
provided we define $Q=8\pi k_{B}\eta\ell_{p}^{4}$. Thus we reach

$F=-\frac{GMm}{R^{2}}\left[1+4\ell_{p}^{2}\frac{\partial{s}}{\partial
A}\right]_{A=4\pi R^{2}},$ (10)

Finally, using Eq. (2) we obtain the modified Newton’s law of gravitation as

$F=-\frac{GMm}{R^{2}}\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma}{4\pi^{2}}\frac{\ell_{p}^{4}}{R^{4}}\right],$
(11)

Thus, with the correction in the entropy expression, we see that the Newton’s
law will modified accordingly. As we mentioned, these corrections are well
motivated from bottom-up quantum gravity theories. The log correction to the
area-entropy relation appears to have an almost universal status, having been
derived from multiple different approaches to the calculation of entropy from
counting microscopic states in different quantum gravity models. Since the
last two terms in Eq. (11) can be comparable to the first term only when $R$
is very small, the corrections make sense only at the very small distances.
When $R$ becomes large, the entropy-corrected Newton’s law reduces to the
usual Newton’s law of gravitation.

Let us compare our result with that obtained in Modesto . The first correction
term originates from the log correction in Eq. (17) of Modesto  is similar to
ones we obtained in Eq. (11), however, it seems that the second correction
term in Eq. (17) of Modesto  is not reasonable. Physically, the effect of the
correction terms on the quantity should be less than the uncorrected quantity.
Similarly, the contribution of the first correction term in the physical
quantity (force here) should be more from the second term and so on. For all
above reasons, we think the second correction term in Eq. (17) of Modesto  is
not correct, and the corrected form is that presented in our note. The origin
of this difference is due to the fact that the second (volume) correction term
to the entropy expression in Modesto  is not indeed a correction term,
although it was argued by Modesto  that this term is also emerged in a model
for the microscopic degrees comprising the black hole entropy in LQG smolin .
Let us stress here that although in the literature there is doubt about the
second correction term in entropy-corrected relation, however, it is widely
believed Zhang  that the next quantum correction term to black hole entropy
have the form $\ell_{p}^{2}/A$, which leads to the resonable correction terms
to Newton’s law of gravitation (11) as we have shown in the present work and
will lead to corrected modified Friedmann equation as we will see later.

Next, we extend our discussion to the cosmological setup. Assuming the
background spacetime to be spatially homogeneous and isotropic which is given
by the Friedmann-Robertson-Walker (FRW) metric

$ds^{2}={h}_{\mu\nu}dx^{\mu}dx^{\nu}+R^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2}),$ (12)

where $R=a(t)r$, $x^{0}=t,x^{1}=r$, the two dimensional metric
$h_{\mu\nu}$=diag $(-1,a^{2}/(1-kr^{2}))$. Here $k$ denotes the curvature of
space with $k=0,1,-1$ corresponding to open, flat, and closed universes,
respectively. The dynamical apparent horizon, a marginally trapped surface
with vanishing expansion, is determined by the relation
$h^{\mu\nu}\partial_{\mu}R\partial_{\nu}R=0$. A simple calculation gives the
apparent horizon radius for the FRW universe

$R=ar=\frac{1}{\sqrt{H^{2}+k/a^{2}}}.$ (13)

We also assume the matter source in the FRW universe is a perfect fluid with
stress-energy tensor

$T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}.$ (14)

Due to the pressure, the total mass $M=\rho V$ in the region enclosed by the
boundary $\mathcal{S}$ is no longer conserved, the change in the total mass is
equal to the work made by the pressure $dM=-pdV$ , which leads to the well-
known continuity equation

$\dot{\rho}+3H(\rho+p)=0,$ (15)

where $H=\dot{a}/a$ is the Hubble parameter. It is instructive to first derive
the dynamical equation for Newtonian cosmology. Consider a compact spatial
region $V$ with a compact boundary $\mathcal{S}$, which is a sphere with
physical radius $R=a(t)r$. Note that here $r$ is a dimensionless quantity
which remains constant for any cosmological object partaking in free cosmic
expansion. Combining the second law of Newton for the test particle $m$ near
the surface, with gravitational force (11) we get

$F=m\ddot{R}=m\ddot{a}r=-\frac{GMm}{R^{2}}\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma}{4\pi^{2}}\frac{\ell_{p}^{4}}{R^{4}}\right],$
(16)

We also assume $\rho=M/V$ is the energy density of the matter inside the the
volume $V=\frac{4}{3}\pi a^{3}r^{3}$. Thus, Eq. (16) can be rewritten as

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}\rho\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma}{4\pi^{2}}\frac{\ell_{p}^{4}}{R^{4}}\right],$
(17)

This is nothing but the entropy-corrected dynamical equation for Newtonian
cosmology. The main difference between this equation and the standard
dynamical equation for Newtonian cosmology is that the correction terms now
depends explicitly on the radius $R$. However, we can remove this confusion.
Assuming that for Newtonian cosmology the spacetime is Minkowskian with $k=0$,
then we get $R=1/H$, and we can rewrite Eq. (17) in the form

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}\rho\left[1-\frac{\beta\ell_{p}^{2}}{\pi}\left(\frac{\dot{a}}{a}\right)^{2}-\frac{\gamma\ell_{p}^{4}}{4\pi^{4}}\left(\frac{\dot{a}}{a}\right)^{4}\right].$
(18)

It was argued in Cai4  that for deriving the Friedmann equations of FRW
universe in general relativity, the quantity that produces the acceleration is
the active gravitational mass $\mathcal{M}$ Pad2 , rather than the total mass
$M$ in the spatial region $V$. With the entropic corrections terms, the active
gravitational mass $\mathcal{M}$ will also modified as well. On one side, from
Eq. (17) with replacing $M$ with $\mathcal{M}$ we have

$\mathcal{M}=-\frac{\ddot{a}a^{2}}{G}r^{3}\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma}{4\pi^{2}}\frac{\ell_{p}^{4}}{R^{4}}\right]^{-1}.$
(19)

On the other side, the active gravitational mass is defined as Cai4

$\mathcal{M}=2\int_{V}{dV\left(T_{\mu\nu}-\frac{1}{2}Tg_{\mu\nu}\right)u^{\mu}u^{\nu}}.$
(20)

A simple calculation leads

$\mathcal{M}=(\rho+3p)\frac{4\pi}{3}a^{3}r^{3}.$ (21)

Equating Eqs. (19) and (21) we find

$\frac{\ddot{a}}{a}=-\frac{4\pi
G}{3}(\rho+3p)\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma}{4\pi^{2}}\frac{\ell_{p}^{4}}{R^{4}}\right].$
(22)

This is the modified acceleration equation for the dynamical evolution of the
FRW universe. Multiplying $\dot{a}a$ on both sides of Eq. (22), and using the
continuity equation (15), after integrating we find

$H^{2}+\frac{k}{a^{2}}=\frac{8\pi
G}{3}\rho\left[1-\frac{\beta}{\pi}\frac{\ell_{p}^{2}}{\rho
R^{2}}\int{\frac{d(\rho
a^{2})}{a^{2}}}-\frac{\gamma}{4\pi^{2}}\frac{a^{2}\ell_{p}^{4}}{\rho
R^{4}}\int{\frac{d(\rho a^{2})}{a^{4}}}\right].$ (23)

Now, in order to calculate the integrations in the correction terms we need to
find $\rho=\rho(a)$. Assume the equation of state parameter $w=p/\rho$ is a
constant, the continuity equation (15) can be integrated immediately to give

$\rho=\rho_{0}a^{-3(1+w)},$ (24)

where $\rho_{0}$, an integration constant, is the present value of the energy
density. Inserting relation (24) in Eq. (23), after integration, we obtain

$H^{2}+\frac{k}{a^{2}}=\frac{8\pi
G}{3}\rho\left[1-\frac{\beta(1+3w)}{3\pi(1+w)}\frac{\ell_{p}^{2}}{R^{2}}-\frac{\gamma(1+3w)}{4\pi^{2}(5+3w)}\frac{\ell_{p}^{4}}{R^{4}}\right].$
(25)

Using Eq. (13) we can further rewrite the above equation as

$\displaystyle\left(H^{2}+\frac{k}{a^{2}}\right)\left[1-\frac{\beta\ell_{p}^{2}(1+3w)}{3\pi(1+w)}\left(H^{2}+\frac{k}{a^{2}}\right)\right.$
$\displaystyle\left.-\frac{\gamma\ell_{p}^{4}(1+3w)}{4\pi^{2}(5+3w)}\left(H^{2}+\frac{k}{a^{2}}\right)^{2}\right]^{-1}=\frac{8\pi
G}{3}\rho.$ (26)

If $\beta$ and $\gamma$ are viewed as small quantities, then the above
equation can be expanded up to the linear order of $\beta$ and $\gamma$. The
result is

$\displaystyle\left(H^{2}+\frac{k}{a^{2}}\right)+\frac{\beta\ell_{p}^{2}(1+3w)}{3\pi(1+w)}\left(H^{2}+\frac{k}{a^{2}}\right)^{2}$
$\displaystyle+\frac{\gamma\ell_{p}^{4}(1+3w)}{4\pi^{2}(5+3w)}\left(H^{2}+\frac{k}{a^{2}}\right)^{3}=\frac{8\pi
G}{3}\rho,$ (27)

which is in complete agreement with the result of Cai33  (see also Shey0 ). In
this way we derive the entropy-corrected Friedmann equation of FRW universe by
considering gravity as an entropic force caused by changes in the information
associated with the positions of material bodies. In the absence of the
correction terms $(\beta=0=\gamma)$, one recovers the well-known Friedmann
equation in standard cosmology. Since the last two terms in Eq. (25) can be
comparable to the first term only when $a$ is very small, the corrections make
sense only at early stage of the universe where $a\rightarrow 0$. When the
universe becomes large, the entropy-corrected Friedmann equation reduces to
the standard Friedman equation.

In summary, we have shown that with the entropy corrections to the area-
relation of the black hole entropy, the Newton’s law of gravitation and the
Friedmann equations will be modified accordingly. These corrections are
motivated from the LQG which is one of the promising theories of quantum
gravity. We derived the correction terms to the Newton’s law of gravity as
well as modified Friedmann equations of the FRW universe starting from the
holographic principle and the equipartition law of energy by using Verlinde’s
argument that gravity appears as an entropic force. In particular, we have
found that an apparently universal log correction to the area-entropy, yields
deviations from Newton’s law and Friedman equations that are identical in form
to those obtained from perturbative quantum gravity. This at once sheds light
on the reason for the universality of the log correction and provides a strong
consistency check on Verlinde’s model.

###### Acknowledgements.

I am grateful to Prof. B. Wang and Prof. R.G. Cai for helpful discussions and
reading the manuscript. This work has been supported by Research Institute for
Astronomy and Astrophysics of Maragha.

## References

  * (1) T. Jacobson, Phys. Rev. Lett. 75, 1260 (1995).
  * (2) C. Eling, R. Guedens, and T. Jacobson, Phys. Rev. Lett. 96, 121301 (2006).
  * (3) T. Padmanabhan, Phys. Rept. 406 (2005) 49;
T. Padmanabhan, Rept. Prog. Phys. 73 (2010) 046901\.

  * (4) M. Akbar and R. G. Cai, Phys. Rev. D 75, 084003 (2007).
  * (5) R. G. Cai and L. M. Cao, Phys.Rev. D 75, 064008 (2007).
  * (6) R. G. Cai and S. P. Kim, JHEP 0502, 050 (2005).
  * (7) B. Wang, E. Abdalla and R. K. Su, Phys.Lett. B 503, 394 (2001);
B. Wang, E. Abdalla and R. K. Su, Mod. Phys. Lett. A 17, 23 (2002);
R. G. Cai and Y. S. Myung, Phys. Rev. D 67, 124021 (2003);
R. G. Cai and L. M. Cao, Nucl. Phys. B 785 (2007) 135.

  * (8) R. G. Cai, L. M. Cao and Y. P. Hu, JHEP 0808 (2008) 090.
  * (9) S. Nojiri and S. D. Odintsov, hep-th/0506212;
S. Nojiri and S. D. Odintsov, hep-th70506212;
A. Sheykhi, in preparation.

  * (10) A. Sheykhi, B. Wang and R. G. Cai, Nucl. Phys. B 779 (2007)1.
  * (11) A. Sheykhi, B. Wang and R. G. Cai, Phys. Rev. D 76 (2007) 023515;
A. Sheykhi, B. Wang, Phys. Lett. B 678 (2009) 434;
A. Sheykhi, Class. Quantum Grav. 27 (2010) 025007.

  * (12) E. P. Verlinde, arXiv:1001.0785 [hep-th].
  * (13) R.G. Cai, L. M. Cao and N. Ohta, Phys. Rev. D 81, 061501(R) (2010) .
  * (14) R.G. Cai, L. M. Cao and N. Ohta, arXiv:1002.113;
Yi Ling, Jian-Pin Wu, arXiv:1001.5324;
Yun Soo Myung, arXiv:1002.0871;
R. Banerjee, B. R. Majhi, arXiv:1003.2312;
S.W. Wei, Y. X. Liu, Y. Q. Wang, arXiv:1001.5238;
Y. X. Liu, Y. Q. Wang, S. W. Wei, arXiv:1002.1062;
R.A. Konoplya, e-Print: arXiv:1002.2818;

  * (15) L. Modesto, A. Randono, arXiv:1003.1998.
  * (16) Y.F. Cai, J. Liu, H. Li, arXiv:1003.4526.
  * (17) T. Zhu and J.R. Ren, Eur. Phys. J. C 62 (2009) 413;
R.G. Cai et al, Class. Quantum Grav 26 (2009) 155018;
D. A. Easson, arXiv: 1003.1528.

  * (18) M. R. Banerjee, S. K. Modak, JHEP 0905 (2009) 063;
S. K. Modak, Phys. Lett. B 671 (2009) 167;
R. Banerjee, S. Gangopadhyay, S. K. Modak, Phys. Lett. B 686 (2010) 181.

  * (19) J. Zhang, Phys. Lett. B 668 (2008) 353;
R. Banerjee and B. R. Majhi, Phys. Lett. B 662, 62 (2008);
R. Banerjee and B. R. Majhi, JHEP 0806 (2008) 095;
R. Banerjee, B. R. Majhi and S. Samanta, Phys. Rev. D 77, 124035 (2008);
R. Banerjee, B. R. Majhi, Phys. Lett. B 674 (2009) 218;
B. R. Majhi, Phys. Rev. D 79 (2009) 044005;
S. Nojiri and S. D. Odintsov, Int. J. Mod. Phys. A 16 (2001) 3273;
J. Lidsey et. al., Phys. Lett. B 544 (2002) 337;

  * (20) C. Rovelli, Phys. Rev. Lett. 77 (1996) 3288;
A. Ashtekar, J. Baez, A. Corichi, and K. Krasnov, Phys. Rev. Lett. 80 (1998)
904.

  * (21) T. Padmanabhan, Class. Quantum Grav. 21, 4485 (2004);
T. Padmanabhan, arXiv:0912.3165;
T. Padmanabhan, arXiv:1003.5665

  * (22) T. Padmanabhan, arXiv:1001.3380 [gr-qc].
  * (23) L. Smolin, arXiv:1001.3668.