\section*{Introduction}

Recent years have witnessed remarkable interest in application of tools of quantum information theory in chemistry \cite{Legeza-2003b, Legeza-2004b, Huang2005, rissler_2006, Pipek-2009, Barcza-2011, McKemmish-2011, boguslawski_2012, boguslawski_2012b, Boguslawski-2013, Kurashige-2013, Fertitta-2014, Duperrouzel-2014, murg_2014, Knecht-2014, Boguslawski-2015, Szalay-2015a, Boguslawski-2015, Freitag-2015, Zhao-2015, Szilvasi-2015, Molina2015, krumnow_2016, stein_2016, stein_2017, kovyrshin_2017, Szalay-2017, Stemmle2018}.
As a prominent example, the performance of state-of-the-art tensor product methods for electronic structure calculations \cite{Kurashige-2009, Murg-2010a, Nakatani-2013, Szalay-2015a, Chan-2015, Reiher-MPO, Wouters-2014e, Gunst-2018} heavily relies on proper manipulation of entanglement \cite{Legeza-2003b, rissler_2006, Barcza-2011, Fertitta-2014, murg_2014, Szalay-2015a, krumnow_2016}. These include
density matrix renormalization group (DMRG) method \cite{White-1992b, white_1993}, which variationally optimizes wave functions in the form of matrix product states (MPS).\cite{schollwock_2011} 

Other important examples represent characterization of electron correlation into its static (strong) and dynamic contributions \cite{boguslawski_2012b}, automatic (black-box) selection of the active spaces \cite{Legeza-2003b, Barcza-2011, Szalay-2015a, stein_2016, stein_2017, Faulstich2018}, or the self-adaptive tensor network states with multi-site correlators \cite{kovyrshin_2017},
all of which harness single- and two-orbital entanglement entropies.
Last but not least, correlation measures based on the single- and two-orbital entanglement entropies have also been employed for the purposes of bond ana\-ly\-sis \cite{Boguslawski-2013, Szilvasi-2015}.

In the preceding work \cite{Szalay-2017}, we have presented the very general \emph{correlation theory of the chemical bond} based on multiorbital correlation measures which goes beyond the scope of two-orbital picture. It is able to properly describe multiorbital bonds, and we have demonstrated its performance on a representative set of organic molecules (aliphatic as well as aromatic).

In the present article, we apply this theory to systems with electron-deficient bonds, i.e., to compounds which have too few valence electrons for the connections between atoms to be described as covalent bonds, and which have always fascinated chemists. First we apply the theory to the notoriously known textbook example of the diborane(6)\footnote{The number in parentheses denotes the number of hydrogen atoms.}
molecule (B$_2$H$_6$) with two-electron three-center bridge bonds and then also to recently characterized diborane(4) \cite{Chou2015} (B$_2$H$_4$) and zero-valent complexes of beryllium \cite{Arrowsmith-2016, Brabec2018}. The neutral form of the latter compound exhibits surprising stability, which was attributed to a strong three-center two-electron $\pi$ bond stretching across the C-Be-C core \cite{Arrowsmith-2016}. Unlike in the previous study \cite{Szalay-2017}, here we work in the bigger detail in a sense that we also employ eigenstates of multiorbital reduced density matrices, which give us additional insights into the character of bonding.

The article is organized as follows: in Sec. II, we briefly present the studied systems, in Sec. III we review the main concepts of the theory of multiorbital correlations, Sec. IV presents the computational details and Sec.~V the results of our calculations which are followed by their discussion, the final Section closes with conclusions.


\section*{Studied systems}

\subsection*{Diboranes}
In its ground state, diborane(6) (Figure \ref{d6}) adopts its most stable conformation with two 
bridging B-H-B bonds, and four terminal B-H bonds. Its structure was first 
correctly measured in 1943 from infrared spectra of gaseous samples 
by an undergraduate student, Longuet-Higgins.\cite{LonguetHiggins1946,LonguetHiggins1943}
Subsequent measurements with electron diffraction confirmed his conclusions\cite{Eberhardt1954}, and X-ray diffraction detected 
further systems with bridging hydrogen bond.\cite{Lammertsma1996}
The B-H-B bridging was considered an atypical electron-deficient covalent chemical bond.\cite{Lipscomb1973}
Diborane(6) is a prominent example of a molecule with three-center two-electron bonds.\cite{Neeve2016}  
As it is a well studied system, we use its B-H-B linkage as a reference to compare with 
bond strengths and properties of more complex systems featuring three-center two-electron bonds.

According to quantum-chemical calculations \cite{Vincent1981,Mohr1986,Curtiss1989,Curtiss1989b,Demachy1994,Alkorta2011}, different species of diborane 
with less than six hydrogens should exist, also featuring the bridging B-H-B bonds. 
However, all candidates are short-lived reaction intermediates, difficult to prepare and to identify.
Hence, no neutral species has been identified experimentally until 2015, when Chou irradiated 
diborane(6) dispersed in neon at $3\;\text{K}$ with far-ultraviolet light, detecting diborane(4), B$_{2}$H$_{4}$ (Figure \ref{d4}).\cite{Chou2015} 
This new species with two terminal hydrogen atoms possesses two bridging hydrogen atoms,
 and so it became the simplest neutral boron hydride identified with such a structural feature.

\begin{figure}[!ht]
  \subfloat[diborane(6)\label{d6}]{%
    \includegraphics[width=4cm]{b2h6.png}
  }
  \hfill
  \subfloat[diborane(4)\label{d4}]{%
    \includegraphics[width=4cm]{b2h4.png}
  }
  \caption{Structures of diborane(6) and diborane(4).}
\end{figure}

\subsection*{Beryllium complexes}
Complexes of metal atoms of the s-block of the periodic table are often found 
in their zero oxidation state due to their exceptional electron donation. For 
their interesting reactivities, these became frequent synthetic targets, competing 
with traditional transition metal complexes.\cite{Power2010,Power2011,Giffin2011} 
We follow the recent experimental work of Arrowsmith, who isolated, for the first time, 
neutral compounds with zero-valent s-block metal, beryllium. \cite{Arrowsmith-2016}
These brightly coloured molecules have very short Be-C bonds and beryllium in 
linear coordination geometries.\cite{Niemeyer1997,Naglav2015,Arnold2015,Lerner2003} 
This indicates strong multiple Be–C bonding. 
According to the theoretical and spectroscopic results, the molecules adopt 
a closed-shell singlet configuration with a Be(0) metal centre. \cite{Arrowsmith-2016}
The complexes are surprisingly stable, and this was ascribed to an unusually 
strong three-center two-electron $\pi$ bond stretching across the C–Be–C unit.
Two bonding mechanisms depicted in Figure \ref{bonding_scheme} are taking place, namely $\sigma$
donation from the carbon doubly occupied sp$^2$ hybrid orbital to empty s and p$_\text{x}$
orbitals on the central Be atom and $\pi$ back donation from the beryllium p$_\text{z}$ orbital
to p$_\text{z}$ orbitals located on C atoms.

We studied two of the proposed systems. First the Be(CAC)$_{2}$ complex (Figure \ref{str}), where CAC 
corresponds to cyclic amino carbene donors, which stabilize the compound 
due to their $\pi$-acidity. \cite{Mondal2013,Li2013} We performed a multireference calculation, 
in order to verify the proposed singlet configuration with a Be(0) metal centre 
and to provide a deeper insight into the bonding scheme.
Next we studied dication [Be(CAC)$_{2}$]$^{2+}$ (Figure \ref{str2+}), in which the removal of two electrons 
disrupts the bridging C–Be–C bond. This allowed us to compare with the former system 
and to determine the stabilization effect of the bridging bond.


\begin{figure}
  \begin{center}
    \includegraphics[width=8cm]{./bonding_scheme.png}
    \caption{Schematic representation of the C-Be-C bonding mechanisms \cite{Arrowsmith-2016}.}
    \label{bonding_scheme}
  \end{center}
\end{figure}

\section*{Methodology}

Recently, the correlation theory of the chemical bond \cite{Szalay-2017} was developed, characterizing bonds based on the correlations among orbitals localized on individual atoms. Simply put, if we think of a simple covalent bond and localize the bonding and antibonding molecular orbitals into their atomic contributions, these localized orbitals will be highly correlated. Therefore, standard two-orbital bonds can be characterized by pairs of strongly correlated localized orbitals, and the strength of the correlation characterizes the strength of the bond from the quantum information theoretical point of view.

The correlation theory of the chemical bond can also be used for the characterization of bonds more involved than the covalent bonds.
The concept in general is to find the finest possible correlation based clustering of the localized orbitals into clusters, so that the clusters are weakly correlated with each other, and the orbitals inside the clusters are strongly correlated.\cite{Szalay-2017} These clusters then form independent bonds of a Lewis structure of a given molecule 
and the strength of the correlation with respect to this clustering refers to the validity of such a representation. The weaker the correlation is, the better the Lewis structure represents bonding.

In order to review the correlation measures\cite{Szalay-2015b}, which will be used in our ana\-ly\-sis,
let us denote the set of (the labels of) localized orbitals with $M$.
We aim at investigating the correlations in an $L\subseteq M$ set of orbitals (cluster).
The state of the full electronic system of the cluster $L$
is given by the density operator $\varrho_L$,
while
the reduced state of a (sub)cluster $X\subseteq L$
is given by the reduced density operator $\varrho_X$
in general.\cite{Ohya-1993,Araki-2003b,Wilde-2013}
If the cluster of orbitals $L$ can be given by a state vector $\ket{\psi_L}$
(for example, when a given eigenstate of the whole molecule is considered),
then its density operator is of rank one, $\varrho_L=\ket{\psi_L}\bra{\psi_L}$, called a pure state.
Its reduced density operator is usually mixed (not of rank one),
which is the manifestation of entanglement\cite{Horodecki-2009} 
between (sub)cluster $X$ and the rest of the cluster $L\setminus X$.
In general, a density operator $\varrho_X$ can be decomposed in infinitely many ways
into state vectors $\ket{\psi_i}$ with mixing weights $p_i\geq0$ as
$\varrho_X=\sum p_i \ket{\psi_i}\bra{\psi_i}$.
The spectral decomposition (where the weights are the $\lambda_i$ eigenvalues,
and the $\ket{\psi_i}$-s are eigenvectors, being \emph{orthogonal})
is a special one, in the sense that its weights are the least mixed.\cite{Schrodinger-1936,Hughston-1993}
Each eigenvector $\ket{\psi_i}$ can be expanded in the occupation number basis,
the square of the absolute value of the coefficients are the weights of the given occupations
in that given eigenvector $\ket{\psi_i}$ \mbox{of weight $\lambda_i$}.


On the first level,
the correlation is defined with respect to a partition\cite{Davey-2002} of the $L$ set of the orbitals,\cite{Szalay-2012,Szalay-2015b,Szalay-2017,Szalay2018cc}
denoted with $\xi =\{X_1,X_2,\dots,X_{\abs{\xi}}\} \equiv X_1|X_2|\dots|X_{\abs{\xi}}$,
where the clusters $X\in\xi$, called parts, are disjoint subsets of the cluster $L$,
and $\cup_{X\in\xi} X = L$.
The measure of correlation among the parts $X\in\xi$ is the
$\xi$-correlation,\cite{,Szalay-2015b,Szalay-2017}
\begin{equation}
\label{eqr:xiCorr}
	C(\xi) := \sum_{X\in\xi} S(X) - S(L).
\end{equation}
Here $S(X)=-\tr(\varrho_X\log_4\varrho_X)$ is the von Neumann entropy.\cite{Ohya-1993,Wilde-2013}
(Note that we use the logarithm to the base $4$, 
which is the dimension of the Hilbert space of an orbital. The resulting numerical values are then the same as of the original measures with natural logarithm\cite{Szalay-2017} given in the units of $\ln4$. Note that $S(X)\leq |X|$, where $|X|$ is the number of orbitals in cluster $X$.)
As a special case, the correlation of two single orbitals,
\begin{equation}
	\label{eqr:twoorbitMutInf}
	C(i|j) = S(i) + S(j) - S(i,j) = I(i|j),
\end{equation}
is the well-known (two-orbital) mutual information,\cite{Ohya-1993,Wilde-2013,Adesso-2016}
which has already been considered in chemistry.%
\cite{Legeza-2006a,Rissler-2006,Barcza-2011,Boguslawski-2013,Kurashige-2013,Mottet-2014,Fertitta-2014,Duperrouzel-2014,Boguslawski-2015,Freitag-2015,Szilvasi-2015,Murg-2015,Barcza-2015,Zhao-2015,Szalay-2017}
(For convenience, we omit the curly brackets $\{\;\}$ in the cases 
when this does not cause confusion.)
For a general partition $\xi$, we have the bound\cite{Szalay-2017} 
\begin{subequations}
	\begin{equation}
		C(X|Y)\leq 2 (|L| -  \max_{X\in\xi} |X|),
	\end{equation}
	which for a bipartition $\xi=X|Y$ reduces to 
	\begin{equation}
		C(X|Y)\leq 2 \min\{|X|,|Y|\}.
	\end{equation}
\end{subequations}

Note that $C(\xi)$ 
is zero for the trivial split $\xi=\top=\{L\}$,
and it takes its maximum, $C(\bot)$, for the finest split $\xi = \bot = \{ \{i\} \;:\; i\in L \}$.
The latter quantity
is also called total correlation,\cite{Lindblad-1973,Horodecki-1994,Legeza-2004b,Legeza-2006b,Herbut-2004}
\begin{equation}
	\label{eqr:totalCorr}
	C_\text{tot}(L) := C(\bot)=\sum_{i\in L}S(i)-S(L).
\end{equation}
(Note that if cluster $L$ is described by a pure state then $S(L)=0$,
and the correlation is entirely quantum entanglement.\cite{Horodecki-2009,Modi-2010,Szalay-2015b}
Moreover, the correlation in a pure state with respect to a bipartition $\xi=X|(L\setminus X)$ 
is just two times the usual entanglement entropy\cite{Bennett-1996,Nielsen-2000,Wilde-2013}
\begin{equation}
        \label{eqr:entEnt}
	C(X|(L\setminus X)) = 2S(X) = 2S(L\setminus X),
\end{equation}
because of the Schmidt decomposition of pure states.\cite{Schmidt-1907,Nielsen-2000,Wilde-2013})


On the second level,
the correlations can be defined
in an overall sense, that is, without respect to a given partition.\cite{Szalay-2017,Szalay-2015b,Szalay2018cc}
The $k$-partitionability correlation and
the $k$-producibility correlation, are\cite{Szalay-2017}
\begin{subequations}
\label{eqr:kppCorr}
\begin{align}
	C_\text{$k$-part}(L) &:= \min_{\xi:\; \abs{\xi}\geq k}C(\xi),\\
	C_\text{$k$-prod}(L) &:= \min_{\xi:\; \forall X\in\xi,\: \abs{X}\leq k }C(\xi),
\end{align}
\end{subequations}
for $1\leq k\leq\abs{L}$.
These characterise the strength of two different (one-parameter-) notions of multiorbital correlations;
those which cannot be restricted inside at least $k$ parts,
and
those which cannot be restricted inside parts of size at most $k$, respectively.\cite{Szalay-2017}


For the cluster $L$,
as special cases,
$C_\text{$\abs{L}$-part}=C_\text{$1$-prod}=C(\bot)$ grabs all the correlations,
it is zero if and only if there is no correlation at all in the cluster $L$.
On the other hand,
$C_\text{$2$-part}=C_\text{$(\abs{L}-1)$-prod}$ is sensitive only for the strongest correlations,
it is nonzero if and only if the cluster $L$ is globally correlated.
Note also that $C_\text{$1$-part}=C_\text{$\abs{L}$-prod}=C(\top) = 0$, by definition.
Beyond these,
there are no such coincidences among the partitionability and producibility correlations
for other values of $k$,
however, the relation $C_\text{$k$-part}\geq C_\text{$(\abs{L}-k+1)$-prod}$ holds.\cite{Szalay-2017}
Also, the following (non stricht) bounds hold \cite{Szalay-2017}
\begin{subequations}
\label{eqr:bounds}
\begin{align}
	C_\text{$k$-part} &\leq 2(k-1), \\
	C_\text{$k$-prod} &\leq 2(\abs{L}-k).
\end{align}
\end{subequations}

\section*{Computational details}

In case of diborane molecules, 
the ground state geometries were optimized with the B3LYP/cc-pVDZ method. For the multiorbital correlation studies, the Pipek-Mezey \cite{Pipek-1989} localized HF/STO-3G molecular orbitals (MOs) were employed \cite{Szilvasi-2015, Szalay-2017} and they were manually hybridized (rotated) to better reflect the chemical environment. 
The quantum chemical (QC-) DMRG method was applied to study the multiorbital correlations in the full orbital following the procedure
outlined in Ref.~\citen{Szalay-2017}.

The ground state geometries of both forms of the beryllium complex, namely Be(CAC)$_2$ and [Be(CAC)$_2$]$^{2+}$ were taken from Ref.~\citen{Arrowsmith-2016} and they correspond to the BP86/def2-TZVPP level of theory. Due to the size of the problem, multiorbital correlation studies by means of QC-DMRG are not feasible in the full orbital space. We have rather chosen a different strategy. Since we were interested only in the bonding of the C-Be-C atomic core, we have selected the complete active space (CAS) of relevant orbitals participating or influencing these bonds. In particular, 2p$_{\text{x}}$ orbitals on both C and Be atoms and 2s orbital on Be, all of them contributing to the $\sigma$ bonds and 2p$_{\text{z}}$ orbitals on both C and neighbouring N atoms and Be, forming or directly influencing the $\pi$ bonds. The CAS orbitals were optimized by means of the CASSCF(10,9)/cc-pVDZ method in case of the neutral complex and CASSCF(8,9)/cc-pVDZ method in case of the dication and again localized using the Pipek-Mezey \cite{Pipek-1989} procedure. They were not hybridized in order to directly compare with the previous work \cite{Arrowsmith-2016} making conclusions about atomic-like orbitals.

All the quantum chemistry calculations except the QC-DMRG ones were performed with the MOLPRO package \cite{MOLPRO}.
The QC-DMRG calculations were carried out using the Budapest QC-DMRG code \cite{budapest_qcdmrg}. Molecular orbitals were visualized with Charmol \cite{charmol}.


\section*{Results}

The results on diborane(6) are summarized in Figure \ref{fig:B2H6} and Table \ref{table:B2H6}, whereas results on diborane(4) are presented in Figure \ref{fig:B2H4} and Table \ref{table:B2H4}. All the Figures depict mutual information of pairs of localized orbitals, defined in \ref{eqr:twoorbitMutInf}, while the Tables contain numeric values of measures of the relevant kinds of multiorbital correlations,  which are discussed below.

In a similar fashion, the results on beryllium complexes are presented in Figures \ref{fig:CBeCv2}, \ref{fig:dication}, and \ref{fig:cbecrotatedsmall} and Table \ref{table:cbec01}.

In the Figures, 
individual localized orbitals are 
represented as black dots and dashed blue lines encircles orbitals belonging to one atom.
The mutual information is plotted as grayscaled edges between the orbitals. Black lines correspond to the strongest
correlations, while light gray lines connect the weakly correlated orbitals.
Based on the mutual information structure, the orbitals are grouped into strongly correlated clusters, which in our examples correspond to either core orbitals or chemical bonds.
The clusters are encircled by red borders.

\section*{Discussion}

\subsection*{Diborane(6)}

\begin{figure}[!ht]
\centering
\includegraphics{fig_B2H6}
\caption{Schematic view of diborane(6) with mutual information: 
	each dot represents a localized orbital, 
	dashed blue line encircles individual atoms,
	edges correspond to mutual information (plot shaded by a 
	logarithmic scale depending on strength) and 
red circles show how the orbitals group into clusters, i.e. independent bonds. Sorted values of the two-obital mutual information are plotted in Appendix \ref{appendix}.}
\label{fig:B2H6}
\end{figure}

\begin{table}[!ht]
\centering
\begin{tabular}{l c c}
 \hline
 \hline
 correlation & abs. value & rel. value \\ [0.5ex]
 \hline
 $C(X_1 \;|\; \text{rest})$ & 0.515 & 8.6\% \\
 $C(X_1 \cup X_2 \;|\; \text{rest})$ & 0.412 & 3.4\% \\
 $C(\orbh{B}{1}{2}{1}, \orb{H}{5}{1}{s} \;|\; \text{rest})$ & 1.852 &  46\% \\
 $C(\orb{H}{5}{1}{s},\orbh{B}{2}{2}{1} \;|\; \text{rest})$ & 1.852 & 46\% \\
 $C(\orbh{B}{1}{2}{2},\orbh{B}{2}{2}{1} \;|\; \text{rest})$ & 2.114 & 53\% \\
 $C(\orbh{B}{1}{2}{1} \;|\; \orb{H}{5}{1}{s})$ & 0.894 & 45\% \\
 $C(\orb{H}{5}{1}{s} \;|\; \orbh{B}{2}{2}{1})$ & 0.894 & 45\% \\
 $C(\orbh{B}{1}{2}{1} \;|\; \orbh{B}{2}{2}{1})$ & 0.605 & 30\% \\
 $C_{\text{$2$-part}}(X_1)$ & 1.500 & 75\% \\
 $C_{\text{$3$-part}}(X_1)$ & 2.394 & 60\% \\
 $C(X_1 \;|\; X_2)$ & 0.309 & 5.2\% \\
 $C(\orb{H}{5}{1}{s} \;|\; \orb{H}{6}{1}{s})$ & 0.042 & 2.1\% \\

 \hline
 \hline
\end{tabular}
\caption{Correlation measures for diborane(6). Relative values are related to the upper bounds. Labeling of localized orbitals corresponds to Figure \ref{fig:B2H6}.}
\label{table:B2H6}
\end{table}

Figure \ref{fig:B2H6} shows that our results fit well the established bonding picture of diborane(6) with two
bridging \mbox{B-H-B} bonds. Let us now discuss in detail how the ana\-ly\-sis of correlations leads to the clustering and to the bonding picture presented in Figure \ref{fig:B2H6}.
We will discuss only the bridging bonds, the core orbitals as well as terminal B-H bonds are well separated, i.e. not correlated with the rest. (This is confirmed by the weak correlation \mbox{$C(X_1 \cup X_2 \;|\; \text{rest})$} in Table \ref{table:B2H6}.)

First we consider the cluster $X_1$ containing sp$^3$ hybrid orbitals on B atoms and 1s orbital on the bridging H atom.
Because of the point group symmetry, the same results hold for cluster $X_2$. As can be seen in Table \ref{table:B2H6}, the correlation (entanglement) of $X_1$ with the remaining orbitals is very weak, only 8.6\% of the maximum value, which indicates that $X_1$ forms an independent three-center bond. 

However, to confirm this conclusion, we have to show that $X_1$ cannot be split further.
If we take separate pairs of orbitals from $X_1$
and measure the correlation with the remaining orbitals, we obtain significantly higher values, in particular 46\% and 53\% of the maximum values (see Table \ref{table:B2H6}), which justifies existence of the three-center bond.

The mutual information of pairs of orbitals within $X_1$ reach rather 
small relative values (45\%, 30\%, see Table \ref{table:B2H6}), however according to our numerical experience,
even strong multicenter bonds typically yield low percentage, never
approaching near the theoretical maxima.\cite{Szalay-2017}
Intuitive perspective suggests that the correlation of one orbital with the others can be thought of as a resource shared among the orbitals.
In other words, all the pairs inside $X_1$ cannot reach the maximum simultaneously,
bounded by entanglement monogamy.\cite{Osborne2006,Coffman2000} 
The formulation of an inequality bounding the mutual 
correlations inside orbital clusters still remains an open problem, to our best knowledge. 
The smaller value of the mutual information between sp$^3$ hybrid orbitals on B atoms reflects their larger internuclear distance.

We employ $k$-partitionability in order to 
quantify and benchmark the strength of the diborane(6)
three-center bonds (in terms of correlation).
As can be seen in Table \ref{table:B2H6},
$C_{\text{$2$-part}}(X_1)$ reaches 75\% of the upper bound and $C_{\text{$3$-part}}(X_1)$ 60\%, which points at a strong bond in $X_1$.

The very weak correlation (entanglement) of $X_1$ with the remaining orbitals  
also indicates that the state of the cluster $X_1$ is close to a pure state. 
Indeed, the eigenstate ana\-ly\-sis of the reduced density operator $\rho_{X_1}$ shows that there is the following two-electron eigenstate 
with a corresponding eigenvalue (probability) of $0.94$
\begin{align*}
	\vspace{-0.2cm}
	&\ket{\psi_{X_1}} = \\
		&+0.2146\ket{--\uparrow\mathrel{\mspace{-1mu}}\downarrow}
		+0.4313\ket{-\uparrow\mathrel{\mspace{-1mu}}\downarrow -}
		+0.2146\ket{\uparrow\mathrel{\mspace{-1mu}}\downarrow--}\\
		&+0.3787\ket{-\downarrow\;\uparrow}
		-0.3787\ket{-\uparrow\;\downarrow}
		+0.2721\ket{\downarrow - \uparrow}\\
		&+0.3787\ket{\uparrow \; \downarrow -}
		-0.3787\ket{\downarrow \; \uparrow -}
		-0.2721\ket{\uparrow - \downarrow},
\end{align*}

\noindent
where the ordering of orbitals in a ket corresponds to \mbox{B$^1$sp$^3$}, \mbox{H$^5$1s}, and \mbox{B$^2$sp$^3$}.
Other eigenstates have probabilities below $0.01$. 
The principal two-electron eigenstate together with the 
above discussion on correlations imply that the three orbitals of $X_1$ form a three-center two-electron bond.
Note that the electron pair exhibits a preferred occupation on H atom, which is 
due to its higher electronegativity when compared to B, as we can see from the principal eigenstate. 
It is in agreement with the expectation values of particle-number-operators
($\bra{\psi} \hat{n}^{(i)}_{\text{el}} \ket{\psi}$)
which for \mbox{B$^1$sp$^3$}, \mbox{H$^5$1s} and \mbox{B$^2$sp$^3$} equal 0.53, 0,95 and 0.53, respectively.

As one can observe in Table \ref{table:B2H6}, 
the main source of correlation between $X_1$ and the remaining orbitals is the correlation with the other three-center bond, $X_2$. Specifically,
the correlation between two bridging H atoms is the strongest, which is caused by higher electron density on these atoms.


\subsection*{Diborane(4)}

\begin{figure}[!ht]
\centering
\includegraphics{fig_B2H4}
\caption{Schematic view of diborane(4) with mutual information: 
	each dot represents a localized orbital, 
	dashed blue line encircles individual atoms,
	edges correspond to mutual information (plot shaded by a 
	logarithmic scale depending on strength) and 
red circles show how the orbitals group into clusters, i.e. independent bonds. Sorted values of the two-obital mutual information are plotted in Appendix \ref{appendix}.}
\label{fig:B2H4}
\end{figure}

\begin{table}[!ht]
\centering
\begin{tabular}{l c c}
 \hline
 \hline
 correlation & abs. value & rel. value \\ [0.5ex]
 \hline
 $C(X_1 \;|\; \text{rest})$ & 1.328 & 22\% \\
 $C(\orbh{B}{1}{2}{1} \;|\; \orb{H}{3}{1}{s})$ & 0.647 & 32\% \\
 $C(\orbh{B}{1}{2}{1} \;|\; \orbh{B}{2}{2}{1})$ & 0.701 & 35\% \\
 $C_{\text{$2$-part}}(X_1)$ & 1.388 & 69\% \\
 $C_{\text{$3$-part}}(X_1)$ & 2.089 & 52\% \\
 $C(X_3 \;|\; \text{rest})$ & 1.438 & 36\% \\
 $C(\orbh{B}{1}{2}{3} \;|\; \orbh{B}{2}{2}{3})$ & 1.245 & 62\% \\
 $C(\orbh{B}{1}{2}{3} \;|\; \orb{H}{3}{1}{s})$ & 0.130 & 6.5\% \\
 $C(X_1\cup X_2\cup X_3 \;|\; \text{rest})$ & 0.535 & 6.7\% \\
 $C(X_1 \;|\; X_2)$ & 0.066 & 1.1\% \\
 $C(X_1 \;|\; X_3)$ & 0.639 & 16\% \\ [1ex]
 \hline
 \hline
\end{tabular}
\caption{Correlation measures for diborane(4). Relative values are related to the upper bounds. Labeling of localized orbitals corresponds to Figure \ref{fig:B2H4}.}
\label{table:B2H4}
\end{table}

In case of diborane(4), the two terminal H atoms are missing and instead a direct covalent bond connecting both B atoms is present \cite{Chou2015}. This is also the picture resulting from our ana\-ly\-sis and depicted in Figure \ref{fig:B2H4}. In comparison to diborane(6), we have the similar three-orbital clusters $X_1$ and $X_2$, but also the two-orbital cluster $X_3$ containing sp$^3$ hybrid orbitals on B atoms and corresponding to the aforementioned \mbox{B-B} bond.

Considering the cluster $X_1$, one can observe in Table \ref{table:B2H4} that
it is more correlated with the remaining orbitals than in case of diborane(6). The value of $C(X_1 \;|\; \text{rest})$ is more than 
two times 
larger, 
but the picture of $X_1$ as a standalone chemical bond is still justifiable. 
Consequently, 
the reduced density operator $\rho_{X_1}$ is more mixed, with 
the principal eigenvalue $0.7803$. The remaining eigenstates share low
probabilities (below $0.081$), and therefore, the picture of $X_1$ as a standalone chemical bond is still a reasonable qualitative description. 
The principal eigenstate is again two-electron, i.e. electron-deficient,  and it has the following form
\begin{align*}
	\vspace{-0.2cm}
&\ket{\psi_{X_1}} =\\ 
	  &-0.2287\ket{--\uparrow\mathrel{\mspace{-1mu}}\downarrow}
		-0.3671\ket{-\uparrow\mathrel{\mspace{-1mu}}\downarrow -}
		-0.2287\ket{\uparrow\mathrel{\mspace{-1mu}}\downarrow--}\\
		& +0.3823\ket{-\downarrow\;\uparrow}
		-0.3823\ket{-\uparrow\;\downarrow}
		-0.2967\ket{\downarrow - \uparrow}\\
		& +0.3823\ket{\uparrow \; \downarrow -}
		-0.3823\ket{\downarrow \; \uparrow -}
		+0.2967\ket{\uparrow - \downarrow}.
\end{align*}

\noindent
Similarly to diborane(6), higher electron density is on the bridging H atom, which is due to its higher electronegativity.

Comparing the two-orbital correlations inside $X_1$ with diborane(6) (Tables \ref{table:B2H6} and \ref{table:B2H4}), one can see a weaker correlation between sp$^3$ hybrid orbitals on B atoms and \mbox{H 1s} orbital, but a slightly stronger correlation between both B-atom-orbitals. This stronger correlation can be certainly assigned to a shorter distance of B atoms
(1.477\AA~vs. 1.784\AA).
Based on the values of $C_{\text{$2$-part}}(X_1)$ and $C_{\text{$3$-part}}(X_1)$, 
the covalent bond corresponding to the cluster $X_1$ is slightly weaker (in terms of correlation) than the same bond in diborane(6).

For the cluster $X_3 = \{\orbh{B}{1}{2}{3},\,\orbh{B}{2}{2}{3}\}$, correlation
with the remaining orbitals 
is stronger than for $X_1$ and $X_2$, which in turn weakens the internal two-orbital correlation.
The major contribution to the correlation of $X_3$ with the remaining orbitals 
originates from 
$C(\orbh{B}{1}{2}{3} \;|\; \orb{H}{3}{1}{s})$, which is still very weak compared to other correlations in the molecule (see Table \ref{table:B2H4}).

The overall correlation of the three bonding clusters $X_1$, $X_2$, and $X_3$ with the rest of the system is similarly 
weak as in diborane(6) 
so the considered bonding can be described independently of the rest of the molecule.


\subsection*{Beryllium complexes}

\begin{figure*}
\centering
\includegraphics{fig_CBeCneutral}
\caption{Schematic view of Be(CAC)$_2$ with mutual information. In order to be consistent with Figures \ref{fig:B2H6} and \ref{fig:B2H4}, all the atoms are depicted (dashed blue line circles), even though only a subset of orbitals (black dots) formed the complete active space. Note that one ring is artificially flipped for clearer correlation picture.
	Sorted values of the two-obital mutual information are plotted in Appendix \ref{appendix}.}
\label{fig:CBeCv2}
\end{figure*}

In order to check how well our bonding picture of Be(CAC)$_2$ fits the one proposed by Arrowsmith \textit{et al.} \cite{Arrowsmith-2016}, we consider the clusters $X_1$ and $X_2$ from Figure \ref{fig:CBeCv2}. The cluster $X_1$ contains p$_\text{z}$ orbitals on C and Be atoms and corresponds to the suggested three-center two-electron $\pi$ bond, whereas the cluster $X_2$ contains \mbox{C p$_\text{x}$} orbitals and \mbox{Be s and p$_\text{x}$} orbitals and corresponds to the $\sigma$ bonds.

Let us start with $X_1$. In Table \ref{table:cbec01}, one can see that the correlation of $X_1$ with the remaining orbitals is higher than 30\% of the maximum value, which means that the picture of the three-orbital C-Be-C $\pi$ bond might be good as a qualitative description, but for a quantitatively adequate description, we might seek to include further orbitals into X$_1$, as shown below. This is also demonstrated by weaker pairwise correlations within $X_1$, especially $C(\orb{C}{1}{2}{p_z} \;|\; \orb{C}{2}{2}{p_z})$, than in the three-orbital bonds discussed above. The correlations of the internal pairs in $X_1$ with the remaining orbitals are large (61\% and 86\%) and clearly cannot be considered as standalone bonds. 

The inaccuracy of the picture of a standalone three-orbital bond is also demonstrated by the reduced density operator $\rho_{X_1}$ (see in Appendix \ref{appendix}), which is much more mixed than in previous cases. It has three dominant eigenvalues, instead of just one. The most significant, nevertheless, corresponds to the two-electron state, which is in agreement with the overall picture of the three-orbital two-electron bond.

\begin{table}[h!]
\centering
\begin{tabular}{l c c} 
 \hline
 \hline
 correlation & abs. value & rel. value \\ [0.5ex] 
 \hline
 $C(\orb{C}{1}{2}{p_z},\orb{C}{2}{2}{p_z} \;|\; \text{rest})$ & 3.424 & 86\% \\
  $C( \orb{C}{1}{2}{p_z},\orb{Be\,}{}{2}{p_z} \;|\; \text{rest})$ & 2.432 & 61\% \\
  $C(X_1 \;|\; \text{rest})$ & 2.032 & 34\% \\
  $C(\orb{C}{1}{2}{p_z} \;|\; \orb{C}{2}{2}{p_z})$ & 0.194 & 9.7\% \\
  $C(\orb{C}{1}{2}{p_z} \;|\; \orb{Be\,}{}{2}{p_z})$ & 0.681 & 34\% \\
  $C_{2-part}(X_1)$ & 1.153 & 58\% \\
  $C_{3-part}(X_1)$ & 1.834 & 46\% \\
  $C(X_1 \;|\; \orb{N}{1}{2}{p_z})$ & 0.915 & 46\% \\
  $C(X_1' \;|\; \text{rest})$ & 0.209 & 2.6\% \\
  $C(\orb{N}{1}{2}{p_z}, \orb{C}{1}{2}{p_z} \; | \; \mathrm{rest\; of}\;X_1')$ & 1.737  & 43\% \\
  $C(\orb{N}{1}{2}{p_z} \;|\; \orb{C}{1}{2}{p_z})$ & 0.560 & 28\% \\
  $C(X_2 \;|\; \text{rest})$ & 0.209 & 2.6\% \\
  $C(\orb{C}{1}{2}{p_x} \;|\; \orb{Be\,}{}{2}{s})$ & 0.765 & 38\% \\
  $C(\orb{C}{1}{2}{p_x} \;|\; \orb{Be\,}{}{2}{p_x}) $ & 0.771 & 39\% \\
  $C_{\text{$2$-part}}(X_2)$ & 1.936 & 97\% \\
  $C(\orb{C}{1}{2}{p_x} \;|\; \orb{Be\,}{}{}{sp_2})$ & 1.912 & 96\% \\
 \hline
 \hline
\end{tabular}
\caption{Correlation measures for Be(CAC)$_2$. Relative values are related to the upper bounds. Labeling of localized orbitals corresponds to Figure \ref{fig:CBeCv2}.}
\label{table:cbec01}
\end{table}

As can be seen in Figure \ref{fig:CBeCv2}, the strongest external correlation of $X_1$ is $C(X_1 \;|\; \orb{N}{1}{2}{p_z})$. 
It results from a conjugation of p$_\text{z}$ orbitals and has a stabilization effect. Notice that the N$^1$-C$^1$-Be-C$^2$-N$^2$ group of atoms form perfectly planar structure (the dihedral angle $\alpha_\text{N-C-C-N} = 179.97^{\circ}$) enabling an efficient overlap of all p$_{\text{z}}$ orbitals, which is necessary for the aforementioned conjugation.

The more accurate bonding picture can thus be obtained by considering the enlarged cluster $X_1^{\prime}$
\begin{equation*}
	X_1 \longmapsto X_1'\equiv X_1\cup\{\orb{N}{1}{2}{p_z}, \orb{N}{2}{2}{p_z}\}.
\end{equation*}
\noindent
The corresponding structure of Be(CAC)$_2$ is depicted in Figure \ref{str}. The cluster $X_1^{\prime}$ is independent of the rest of the molecule. This follows from the negligible correlation of $X_1^{\prime}$ with the remaining orbitals (see Table \ref{table:cbec01}). Employing the standard notation, the $\pi$ electron bond can be denoted as $\Pi^{6}_{5}$, i.e. six-electron (two electrons from the Be atom and two from each N atom lone pair) five-center bond, which is confirmed by the particle number expectation value of $6.003$. 

\begin{figure}[!ht]
  \subfloat[Be(CAC)$_2$\label{str}]{%
    \includegraphics[width=4cm]{str.png}
  }
  \hfill
  \subfloat[Be(CAC)$_2^{2+}$\label{str2+}]{%
    \includegraphics[width=4cm]{str2+.png}
  }
  \caption{Structures of studied beryllium complexes suggested by our correlation ana\-ly\-sis.}
\end{figure}

In order to verify the suggested $\pi$ back-donation mechanism \cite{Arrowsmith-2016}, or in other words probe the local electronic configuration of the Be atom, we have also performed the correlation ana\-ly\-sis for the dication [Be(CAC)$_2$]$^{2+}$. As can be seen in Figure \ref{fig:dication}, the difference between the correlation picture of Be(CAC)$_2$ and [Be(CAC)$_2$]$^{2+}$ are almost missing correlations inside the cluster $X_1$, which is for example demonstrated by the negligible value of $C_{2-part}(X_1) = 0.093$. Also the reduced density operator $\rho_{X_1}$ is highly mixed and without the dominating two-electron eigenstate (see in Appendix \ref{appendix}). On the other hand, the correlations in the cluster $X_2$ remain practically unchanged. 

\begin{figure*}
\centering
\includegraphics{fig_CBeCplus}
\caption{Schematic view of [Be(CAC)$_2$]$^{2+}$ with mutual information. Note that one ring is artificially flipped for clearer correlation picture.}
\label{fig:dication}
\end{figure*}

This is in agreement with the picture of Be atom having originally two electrons in the p$_z$ orbital. When the C-Be-C bond is formed, they are shared with C-atom p$_z$ orbitals through the back donation mechanism, as is depicted in Figure \ref{bonding_scheme}. These two $\pi$ electrons are missing in case of the dication and the aforementioned $\pi$ bond is clearly not formed.

Another feature of the correlation picture from Figure \ref{fig:dication} is that there are considerably stronger pairwise correlations between C and N-atom p$_z$ orbitals \mbox{[$C(\orb{N}{1}{2}{p_z} \;|\; \orb{C}{1}{2}{p_z}) = 1.691$]}. They are indeed of the strength of donor-acceptor bonds \cite{szilvasi_2015}. We thus assign double bonds between N and C atoms to the dication, as is depicted in Figure \ref{str2+}. The $\pi$ bonds are formed from the originally doubly filled N p$_\text{z}$ orbitals and empty C p$_\text{z}$ orbitals. The existence of these $\pi$ bonds is also confirmed by almost perfectly planar environment with the dihedral angle $\alpha_\text{H-N-C-Be} = 1.26^{\circ}$. Note that in Figure \ref{fig:dication}, we can see only the part of the double bond corresponding to the $\pi$ bond - the $\sigma$ bonds along the rings are excluded from the active space.

Let us now turn to the $X_2$ cluster of Be(CAC)$_2$, i.e. to the $\sigma$ bonding. The correlation of $X_2$ with the remaining orbitals is insignificant, however splitting of the four-orbital cluster into two $\sigma$ bonds is not possible in this basis. It may therefore seem that the correct $\sigma$ bond is also multiorbital. It is, however, only the artefact of the atomic-like basis, which was used in order to directly compare with Ref.~\citen{Arrowsmith-2016}. By simple rotation of Be s and p$_x$ orbitals (forming the sp hybrids as in case of diboranes), one can form two independent $\sigma$ bonds, essentially without influencing the rest of the system, which can be seen in Figure \ref{fig:cbecrotatedsmall}. Note that in the correlation theory of the chemical bond, superposing orbitals is allowed if this does not affect their locality too much. So superposing orbitals on different atoms is usually forbidden, while doing the same on a given atom is allowed.

\begin{figure}
\centering
\vspace{-0.3cm}
\includegraphics{fig_CBeCneutralrotsmall}
\caption{Mutual information of the Be(CAC)$_2$ $\sigma$ bonding channel after the rotation of $\orb{Be\,}{}{2}{s}$ and $\orb{Be\,}{}{2}{p_x}$ orbitals.} 
\label{fig:cbecrotatedsmall}
\vspace{-0.1cm}
\end{figure}

Last but not least, we would like to compare the strength of both contributions to the bonding of the \mbox{C-Be-C core,} namely the $\pi$ bond (clusters $X_1$ and $X_1^{\prime}$) and the $\sigma$ bond (cluster $X_2$). In the previous study \cite{Arrowsmith-2016}, it was shown by means of the 
energy decomposition ana\-ly\-sis combined with natural orbitals for chemical valence (EDA-NOCV) \cite{Krapp2007} that the $\pi$ bonding channel is considerably more energy stabilizing than the $\sigma$ channel. Let us now check the strength of both contributions by means of the correlations.

When using the more rough (three-orbital two-electron) description of the $\pi$ bond (cluster $X_1$) and by looking at Table \ref{table:cbec01}, one can see on the values of $C_{2-part}$ that $\sigma$ bonds are considerably stronger than the $\pi$ bond. 
The situation, however, dramatically changes  when we use more accurate conjugated description (cluster $X_1^{\prime}$).
Note that since we are not interested in a split dissecting N and C, 
(because the N-C bond on the ring is stabilized by another $\sigma$ bond, not visible on the plots),
we have to consider the p$_z$ orbitals together on N and C atoms of the same ring,
and calculate \mbox{$C_{2-part}( \orb{N}{1}{2}{p_z}, \orb{C}{1}{2}{p_z} \; | \; \orb{Be\,}{}{2}{p_z}\; | \;\orb{N}{2}{2}{p_z}, \orb{C}{2}{2}{p_z})$,} which turns out to be
$C(\orb{N}{1}{2}{p_z}, \orb{C}{1}{2}{p_z} \; | \; \mathrm{rest\; of}\;X_1')$.
Also note that in this paragraph we compared the absolute values of the correlation measures, because the clusters are of different sizes. 

In Table \ref{table:cbec01}, one can see that the more accurate description of the $\pi$ bond ($X_1^{\prime}$) makes it of a similar strength as the $\sigma$~bonds \mbox{(1.737 vs 1.936)}. We believe that our results describe the nature of a Be(CAC)$_2$ bonding reliably, especially because we have used the genuine multireference description unlike in Ref.~\citen{Arrowsmith-2016}, where the ana\-ly\-sis was based on the density functional theory (BP86 functional).
We would also like to note that we have studied slightly different system than in Ref.~\citen{Arrowsmith-2016}. In our case all substituents were replaced by hydrogen atoms. This, however, should not influence the electronic structure of the C-Be-C core.
Also note that using the s and p$_{\text{x}}$ orbitals on beryllium in $X_2$ was only for the purpose of comparison with the previous study\cite{Arrowsmith-2016}. For having a more physical picture, we should use the hybridized orbitals \mbox{(Figure \ref{fig:cbecrotatedsmall})}, by which $X_2$ consists of two simple covalent bonds. Then, in order to characterize the strength of the bond, we would wave to consider the Be sp$_1$ and sp$_2$ orbitals together and calculate 
$C_{2-part}( \orb{C}{1}{2}{p_x} \; | \; \text{Be\;sp}_1,\,\text{Be\;sp}_2 \; | \; \orb{C}{2}{2}{p_x}) = 1.935$. 
Nevertheless, since this value is nearly identical to $C_{2-part}(X_2)$, the conclusion is the same.

\section*{Conclusions}

In this article, we have reviewed the recently developed correlation theory of the chemical bond \cite{Szalay-2017} and applied it on molecules with multicenter electron-deficient bonds. We have demonstrated the usefulness of our methodology in characterizing molecular bonding properties by
fingerprints of correlations among individual orbitals which form these types of bonds.

We have verified the computational procedure on a textbook molecule with electron-deficient bonds, namely diborane(6), and further characterized bonding in diborane(4) and zero-valent complexes of beryllium with intricate bonding patterns. In all the cases, our results fit well with known bonding pictures or previous theoretical predictions.
We have therefore proved capabilities of our new method to reliably describe bonding in complex molecular systems.

In case of the Be(CAC)$_2$ molecule, we have also compared both contributions to the C-Be-C bonding ($\sigma$ and $\pi$), finding, in contrast to the previous study \cite{Arrowsmith-2016}, the $\sigma$ and $\pi$ contributions of a similar strength, in the sense of correlational quantities. We believe that our result is reliable and attribute the discrepancy with the previous study to the single reference description employed in Ref.~\citen{Arrowsmith-2016}, which may not be accurate enough 
in this multireference case.

Finally, we would like to note that, despite employing the DMRG method \cite{White-1992b, white_1993} for calculations of subsystem reduced density matrices, the theory presented in this article is general and other correlated methods
can in principle be employed as well \cite{Boguslawski-2015}. Especially in cases of large molecules with the electronic structure dominated by the dynamical correlation for which the DMRG description may be unnecessary and computationally prohibitive.

\section*{Acknowledgements}
This work has been supported by the \textit{Czech Science Foundation} (grants no. 16-12052S, 18-24563S, and 18-18940Y), 
\textit{Czech Ministry of Education, Youth and Sports} (project no. LTAUSA17033),
and the \textit{Hungarian-Czech Joint Research Project MTA/16/05}.
G.B., Sz.Sz. and {\"O}.L.~are supported by the National Research, 
Development and Innovation Fund of Hungary (NRDIFH) within the Researcher-initiated Research Program (project \mbox{Nr: NKFIH-K120569)} 
and the \textit{``Lend\"ulet'' Program} of the Hungarian Academy of Sciences (HAS).
Sz.Sz. and {\"O}.L. are supported by the Quantum Technology National Excellence Program 
(project \mbox{Nr: 2017-1.2.1-NKP-2017-00001)} of NRDIFH.
G.B. and Sz.Sz. are also supported by the \textit{``Bolyai'' Research Scholarship} of HAS.
{\"O}.L. also acknowledges financial support from the Alexander von Humboldt foundation.



\bibliographystyle{apsrev4-1}
