\section{Introduction}

\IEEEPARstart
{N}{on} intrusive load monitoring (NILM)\cite{Hart1992Nonintrusive} represents a pure computation
technique which identifies and extracts the power consumption of
individual appliances from measurements of the aggregate power usage for the entire home.
Comparing to using appliance-level energy monitors, NILM minimizes maintenance and installation costs of sensors.
A significant application of NILM is to produce the individual electricity bills using readings from a single smart meter, which can deepen energy feedback leading to more efficient use of appliances and reduce the energy consumption.

In the last decade, with the smart energy meters
have been being deployed in many countries, there is a growing number of datasets developed specifically for this research field. Thus a wide
variety of artificial intelligence and machine learning techniques applied to this problem, such as Hidden Markov Model (HMM) and its
variants\cite{Parson2012Non}, \cite{Zoha2012Non}, \cite{Kim2012Unsupervised}, \cite{Egarter2014PALDi}, decision tree\cite{Jing2014Non}, graph signal processing\cite{He2018Non}, support vector machines\cite{Lin2010Applying}, and non-negative tensor factorization\cite{Figueiredo2014Electrical}.

More recently, with the current breakthrough of deep neural networks in image classification\cite{krizhevsky2012imagenet}, speech recognition\cite{graves2014towards}, machine translation\cite{sutskever2014sequence},
DNN have regained their interests in addressing the disaggregation problem. Mauch and Yang\cite{mauch2015new} exploited a generic two-layer bidirectional
Recurrent Neural Network architecture featuring Long Short Term Memory (LSTM)\cite{hochreiter1997long} units in extracting single appliance profiles. In a latter work, Mauch
and Yang\cite{mauch2016novel} continued using a combination of discriminative and generative models in a two-stage eventless extraction of appliance profiles.
Kelly and Knottenbelt\cite{kelly2015neural} evaluated and compared three neural network architectures, a RNN architecture with LSTM units similar to \cite{mauch2015new}, a de-noising Auto-Encoder(dAE)\cite{vincent2010stacked}, and a regression-based disaggregator which estimates the main key points of an activation cycle of the target appliance. Nascimento \cite{do2016applications} applied three deep neural network architectures, a basic convolutional dAE, a RNN, and a ResNet-based model\cite{he2016deep}, by introducing several improvements such as redefining the loss function, exploiting batch normalization \cite{ioffe2015batch}, and applying residual connections\cite{he2016deep}. All these
papers use synthetic data by summing all sub-meters, which limits the amount of noise as appliances not sub-metered would be excluded.
He and Chai\cite{he2016empirical} applied two architectures, a convolutional dAE and an RNN, with different kernel sizes applied to parallel convolutional layers.  Zhang et al.\cite{zhang2018sequence} simplified the objective of the dAE architecture in\cite{kelly2015neural} to predict a single time instance of the target appliance profile for a given window of the aggregated power. Murray et al.\cite{murray2018transferability} applied two architectures, a Convolutional Neural Network (CNN) and a bidirectional RNN featuring Gated Recurrent Units\cite{chung2014empirical}, to both estimates the state and the average consumption of targeted appliances. However, among these
works, each disaggregation window length (and consequently the width of subsequent layers) depends on the specified appliance being monitored, which
is not feasible in practical application. Said and Yang\cite{barsim2018feasibility} applied a deep fully convolutional neural networks to estimate a variety of load categories, but the model contains large number (i.e., 44) layers which is costly computational.
Additionally, Lange et al.\cite{lange2016bolt} adopted a deep neural network with constrained binary and linear activation units
in the last two layers that estimate the on-off activation vector of each load, however it was applied on very high frequency current measurements.

Review all above works, each deep neural network is designed to disaggregate a window of aggregate power and to estimate the on-off state or the power consumption of each target appliance.
In the training process, each neural network needs to use large fraction or even entire power data of several houses, and the test usually carried only in one house (and some with data preprocessing), which is lack of persuasiveness. In fact, except estimate the state or the power consumption of appliance, there is an emerging field of NILM in estimating the operational state change of each appliance. In NILM, operational state change is defined as substantial statistical change in the aggregate power measurement occurs that indicates that one or more appliances have been switched on or off, or change their operational state. After such operational state change are identified, it will be classified into the predefined appliance categories. Through estimating the operational sate change of appliance, the state of appliance at each time instance can be determined by the sign of the power variation. Recently, the operational state change estimation in NILM gains some breakthroughs by using  GSP\cite{He2018Non}, but there is few research that applied deep neural network to estimate the operational state change of appliance. In this work, we apply a DNN to estimate the operational state change of each appliance with single aggregate data. And show competitive performance through experiments comparing to other two benchmark methods and a RNN architecture.
The main contributions of this paper are:

(1) We show the neural network can perform well in estimating the operational change of appliance with a simple architecture.

(2) Our proposed solution represents a significant reduction in complexity compared to previous works\cite{mauch2015new}, \cite{mauch2016novel}, \cite{kelly2015neural}, \cite{do2016applications}, \cite{he2016empirical}, \cite{zhang2018sequence}, \cite{murray2018transferability}, \cite{barsim2018feasibility}, \cite{lange2016bolt} in estimating the state or the
power consumption of appliance.

(3) Unlike previous deep architecture which need to input a window of aggregate power, our proposed solution can estimate the results at each time instance, which is more meaningful and feasibility across a variety of appliances.

(4) The training process only use small percentage of data, and the test
is performed on the rest of raw dataset without any processing, which is
occurs in many previous works, like balanced test data\cite{murray2018transferability} and sythesis data\cite{mauch2015new}, \cite{mauch2016novel}, \cite{kelly2015neural}, \cite{do2016applications}.

The rest of the paper is organized as follows. In section \ref{section 2}, we describe the NILM task and the operational state change. In section \ref{section 3}, we describe the architecture of the proposed neural network. Section \ref{section 4} shows the experiment. Conclusion and future work are discussed in section \ref{section 5}.

\section{Load Disaggregation Based on Operational State Change}
\label{section 2}
In this section, we will introduce the load disaggregation task and the operational state change of each appliance.

\subsection{Load Disaggregation}
Let $\mathcal{A}$ be the set of all known appliances in a house and $P(t_i)$ be the
aggregate power of the entire house measured at time $t_i$. Without
loss of generality, in the following, we denote $P(t_i)$ as $P(t_i) = P_i \geq 0$.
Let $P^a_j \geq 0$ be the power load of appliance $a \in \mathcal{A}$ at time instance
$t_j$.  Then,
\begin{equation}
P_i = \sum_{a=1}^{|A|} P^a_i + n_i,
\end{equation}
where $n_j$  is the measurement noise that not submetered. The disaggregation task is for $i=1,...,T$   and  $a \epsilon \mathcal{A}$, given the $P_i$  to estimate the  $P^a_i$.

\subsection{Operational State Change}
Naturally, $\Delta P_i = P_{i+1} - P_i, i=1,...,T$ and $\Delta P^a_i = P^a_{i+1} - P^a_i, i=1,...,T$
respectively correspond to the variation of the aggregate power and appliance $a$ power
measured at time $t_i$. The classification labels of each appliance $a$ at time $t_i$ is denoted
as $s^a_i$. The $s^a_i$ is defined
as following:
\begin{equation}
 s^a_i =
 \begin{cases}
  1, & \text{for $|\Delta p^a_i| \geq Thr_a$} \\
  0, & \text{for $|\Delta p^a_i| < Thr_a$},
 \end{cases}
\end{equation}
where $Thr_a \geq 0$ is a power threshold for appliance $a$. Like in previous literature\cite{He2018Non}, $Thr_a$ is set to half of difference between mean values of appliance $a$'s adjacent states, which is observed through training process. When $s^a_i$ is set to 1, it means appliance $a$ changed its operational state
 (e.g., switched on/off) at time $t_i$; When $s^a_i$ is set to 0, it means appliance $a$ didn't change its operational state at time $t_i$. $N_{pos}$ and $N_{neg}$ denote the number of the positive samples ($s^a_i = 1$) and negative samples ($s^a_i = 0$) respectively.

\subsection{Detect Operational State Change}
As can be seen in Fig. \ref{Figure 1}, at each time sample $t_i$, given the $\Delta P_i$ (i.e., the
rising edge or the falling edge or the horizontal line), we estimate the $s^a_i$ which show whether appliance $a$ changed its operation state at time $t_i$.
\begin{figure}
\centering
\includegraphics[scale=0.5]{powervariation3.png}
\caption{An example of the operational state change.}
\label{Figure 1}
\end{figure}

\section{Proposed Network Architecture}
\label{section 3}
In the previous subsection, we have framed the operational state change detection task, now we describe how we adapted neural network to solve it.

Fig. \ref{Figure 2} show the detail architecture of the neural network. We propose a deep neural network (DNN) to calculate the classification labels for each appliance $a$.  The model consists of five layers reaching 1200 trainable parameters. Specifically, given the  aggregate power variation of the entire house $\Delta P_i$ at time $t_i$, firstly we initialize the input,
\begin{gather}
 x^{(0)} = |\Delta P_i|,
\end{gather}
by using the absolute value of the aggregate power variation, then input it through the neural network. Each layer includes a sequence of elementary operations shown in the figure and briefly introduced in the sequel.

\textbf{Fully Connected:} the linear operation of each layer defined as a weighted multiplication and adding a bias:
\begin{equation}
f(x) = W^{(d)}x^{(d-1)} + b^{(d)},
\end{equation}
where $d$ denotes the layer number, $1 \leq d \leq D$, $W^{d}$ and $b^{d}$ denote the weight and bias of this layer respectively.

\textbf{Tanh:} is the non-linear activation function defined as:
\begin{equation}
f(x) = \frac{sinh(x)}{consh(x)} = \frac{e^x - e^{-x}}{e^x + e^{-x}},
\end{equation}
the non-linear activation is the component which convert the linear operation to non-linear.

\textbf{Batch Normalization:} is a composition of two affine transformation applied to the output of each layer based on mini-batch statistic:
\begin{equation}
f(x) = \gamma \widehat{x} + \beta = \gamma \frac{x - \mu_{\beta}}{\sigma_\beta} + \beta,
\end{equation}
where $x$  is the original output of a unit, $\mu_\beta$  and $\sigma_\beta$  are the sample mean and standard variation of all outputs of this neuron over the mini-batch  $B$,  $\gamma$ and $\beta$  are two learnable parameters. The batch normalization is very important in neural network. When we apply it to the neural network, the performance is improved a lot.

\textbf{Softmax:} is an activation function applied to the output of last layer in the model:
\begin{equation}
f(x)_i = \frac{e^{x_i}}{\sum^K_{k=1} e^{x_k}} \text{ for } i = 1, ..., K,
\end{equation}
which calculate the normalized probability distribution of a output vector, and $K$  denotes the number of the class (in this paper,  $K=2$).

 After the features are extracted through the neural network, it will be passed to the:
\begin{gather}
\mathit{ p_i}^a = softmax(x^{(D)}),
\end{gather}
to calculate a binary likelihood probabilities $\mathit{p_i}^a =(p_0, p_1)$, which represent the probabilities of whether the appliance changed it's state or not respectively.

\begin{figure}
\centering
\includegraphics[scale=0.5]{model_architecture.png}
\caption{The detail architecture of the proposed model.}
\label{Figure 2}
\end{figure}

\subsection{Loss Function}
The cost for our neural network training is the negative log-likelihood:
\begin{gather}
J({\mathbf{p}^a}, \mathbf{{s}^a}) = \sum_{i} -\log \mathit{p}^a_i[s^a_i],
\end{gather}
where $\left[ \cdot \right]$ is the index operate.

\section{Experiment}
\label{section 4}

\subsection{Dataset}
We evaluate the proposed model on the REDD dataset\cite{kolter2011redd} downsampled to 1 min resolution as in \cite{He2018Non}.

\textbf{REDD} dataset\cite{kolter2011redd} is a dataset for
energy disaggregation.  The dataset contains about half month power consumption from real
homes in US, for the whole house as well as for each individual circuit in
the house (labeled by the main type of appliance on that circuit). The
main types of the appliances are: Dishwasher (DW), Refrigerator (REFR),
Microwave (MW), Kitchen outlet (KO), Stove (ST), air-conditioning (AC), Electronics (EL), Wash Dryer (WD). In this experiment, we use three
houses energy data: House 1, House 2, House 6. A week of of the REDD House 1 data can be seen in Fig. \ref{Figure 3}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{a_week_of_redd_house1.png}
\caption{A week power data of REDD House1.}
\label{Figure 3}
\end{figure}

\subsection{Training}
We train the proposed architecture using a small part data of House1, House2, House6 in REDD. The total samples of House1, House2 and House6
are 25946, 19856, 17605 respectively. The training detail of three house are respectively shown in Table \ref{Table 1}, Table \ref{Table 2}, Table \ref{Table 3}
\begin{table}[h]
\caption{THE TRAINING DETIALS OF HOUSE 1}
 \label{Table 1}
\begin{center}
 \begin{tabular}{|c|ccccc|}
  \hline
  Appliance & REFR & MW & DW & KO & WD \\
  \hline
  Training Samples & 2000 & 2000 & 5000 & 2000 & 8000 \\
  Thresold(Watts) & 150 & 750 & 210 & 550 & 1300 \\
  \hline
 \end{tabular}
\end{center}
\end{table}

\begin{table}[h]
 \caption{THE TRAINING DETAILS OF HOUSE 2}
 \label{Table 2}
 \begin{center}
  \begin{tabular}{|c|cccc|}
   \hline
   Appliance & REFR & MW & KO & ST \\
   \hline
   Training Samples  & 2000 & 2000 & 2000 & 4920\\
   Thresold(Watts) & 85.5 & 920 & 528 & 204\\
   \hline
  \end{tabular}
 \end{center}
\end{table}

\begin{table}[h]
 \caption{THE TRAINING DETAILS OF HOUSE 6}
  \label{Table 3}
 \begin{center}
  \begin{tabular}{|c|ccccc|}
   \hline
  Appliance & REFR & AC & EL & KO & ST \\
  \hline
  Training Samples & 2000 & 2000 & 4640 & 3000 & 3445 \\
  Thresold(Watts) & 74.5 & 862  & 225 & 660 & 1700 \\
  \hline
  \end{tabular}
 \end{center}
\end{table}

\subsection{Augment positive samples}
In real life, many appliances are not commonly used, such as stove, wash dryer etc., which means the positive samples
is much less than the negative samples in the training set, and that will lead the model output a unbalance result (classify all the samples to
negative).

Let us denote $N_{pos}$ and $N_{neg}$ as the number of the positive samples ($s^a_i=1$) and negative samples ($s^a_i = 0$) in the
training set respectively.
In order to solve the unbalance issue, during training, we augment the positive samples
in the following way:
\begin{enumerate}
\item calculate the ratio between the negative samples and positive samples in the
training dataset, $\eta = N_{neg} / N_{pos}$;
\item determine a ratio of positive to negative $\alpha$, then calculate the positive
scaling factor $\sigma = \eta * \alpha$;
\item duplicate the positive samples by $\sigma$ times, and randomly insert them into the
original training dataset.
\end{enumerate}
The $\alpha$ in this paper is set to 1:8.

\subsection{Evaluation metrics}
The evaluation metrics used are precision (PR), recall (RE) and F-Measure ($F_M$)\cite{Delen2008Advanced} defined as:
\begin{gather}
 PR = TP / (TP + FP)\\
 RE = TP / (TP + FN)\\
 F_M = 2 * (PR * RE)/(PR + RE),
\end{gather}
where true positive (TP) is recorded when the state of the detected appliance was actually changed, false positive (FP) is recorded when the
state of the detected appliance was not changed, and false negative (FN) indicates that the state changed appliance was not detected.
Precision captures the correctness of detection, and the high Recall implies a higher percentage of appliance state changes are detected
correctly.

\subsection{Setup}
We use PyTorch\cite{Paszke2017Automatic} to develop our model and optimize the model by  Adam\cite{Kingma2014Adam} optimizer with a base learning rate 1e-4, momentum 0.99.  The data
is processed by the Numpy and Pandas.

\subsection{Comparison with benchmarks}
We compare our model with two benchmark methods, HMM-based approach\cite{Parson2012Non} and GSP-based approach\cite{He2018Non}. The results are
shown in Table \ref{Table 4}, Table \ref{Table 5}, Table \ref{Table 6} for Houses 1, 2, 6
respectively.

\begin{table}[h]
 \caption{COMPARISION RESULTS OF THREE METHODS IN HOUSE 1}
 \label{Table 4}
\begin{center}
 \begin{tabular}{|c|ccccc|}
  \hline
  Appliance & REFR & MW & DW & KO & WD \\
  \hline
  $F_{M_{NN}}$ & 0.88 & 0.76 & 0.47 & 0.64 & 0.88 \\
  $F_{M_{GSP}}$ & 0.88 & 0.70 & 0.57 & 0.39 & 0.89 \\
  $F_{M_{HMM}}$ & 0.97 & 0.50 & 0.13 & 0 & 0 \\
  \hline
 \end{tabular}
\end{center}
\end{table}

\begin{table}[h]
 \caption{COMPARISION RESULTS OF THREE METHODS IN HOUSE 2}
 \label{Table 5}
 \begin{center}
  \begin{tabular}{|c|cccc|}
   \hline
   Appliance & REFR & MW & KO & ST \\
   \hline
   $F_{M_{NN}}$ & 0.85 & 0.97 & 0.91 & 0.83\\
   $F_{M_{GSP}}$ & 0.84 & 0.93 & 0.88 & 0.86\\
   $F_{M_{HMM}}$ & 0.90 & 0.47 & 0.68 & 0.21\\
   \hline
  \end{tabular}
 \end{center}
\end{table}

\begin{table}[h]
 \caption{COMPARISION RESULTS OF THREE METHODS IN HOUSE 6}
 \label{Table 6}
 \begin{center}
  \begin{tabular}{|c|ccccc|}
   \hline
  Appliance & REFR & AC & EL & KO & ST \\
  \hline
  $F_{M_{NN}}$ & 0.80 & 0.89 & 0.70 & 1 & 0.90 \\
  $F_{M_{GSP}}$ & 0.77 & 0.88 & 0.66 & 0.88 & 0.92 \\
  $F_{M_{HMM}}$ & 0.99 & 0.12 & 0.03 & 0 & 0 \\
  \hline
  \end{tabular}
 \end{center}
\end{table}

As can be seen from Tables \ref{Table 4}, Table \ref{Table 5} and Table \ref{Table 6}, our proposed solution
outperforms the other two methods in many cases, which shows
the superiority of our method in estimating the operational state change of appliance.
Specifically, our method significantly outperforms the HMM-based method in all appliances except the refrigerator.
This is mainly due to continuous and sole operation of refrigerator, hence there is large
available data for learning and improving initial HMM model. The poor performance of HMM for other appliances can be attributed to
the short training period.
The proposed model shows better or similar performance to the GSP-based method, especially Kitchen Outlets in three houses. This is mainly
due to the large fluctuations during operation, thus the GSP-based method cannot accurately capture the appliance operation.

All the results for multi-state appliances (dishwasher in House 1) are generally worse for all three methods. This is due to the similarity between the
refrigerator load and low-state of the dishwasher, so they are often `hidden' in the baseload and noise. And multi-state appliances are not
used frequently, thus it is more difficult to extract during the training phase.

Additionally, the training details in Table \ref{Table 1}, Table \ref{Table 2} and Table \ref{Table 3} show the training condition differences between different kinds of appliances. Stove usually
needs more training data, which is due to the fact that Stove normally has short operation time
and relative high power, thus the neural network needs more data to generate the probabilistic
models to capture the appliance operational state change. Electronics and
Wash Dryer  are not used often in common houses, so the model also needs more data to
learn the statistically pattern of these two appliances.

\subsection{Comparison with RNN}
We also compare our solution with a GRU-based RNN architecture. The GRU is a variant of the LSTM unit, especially designed for time series
data. Comparing to LSTM, GRUs have fewer parameter and are more suited to online learning.

We adapt the GRU-based RNN to estimate two time-step and three time-step window of aggregate data, and the results are shown in Table \ref{Table 7}, \ref{Table 8}, and \ref{Table 9} for REDD House 1, 2, 6 respectively. The architecture details of RNN can be seen in Fig. \ref{Figure 4}. For consistency, the number of layers and hidden features in RNN are same as those in DNN. And the training and test setup are also same as DNN setup.

\begin{figure}
\centering
\includegraphics[scale=0.5]{rnn_architecture.png}
\caption{Architecture of compared RNN model.}
\label{Figure 4}
\end{figure}

\begin{table}[h]
 \caption{COMPARISION RESULTS WITH RNN IN HOUSE 1}
 \label{Table 7}
\begin{center}
 \begin{tabular}{|c|ccccc|}
  \hline
  Appliance & REFR & MW & DW & KO & WD \\
  \hline
  $F_{M_{NN}}$ & 0.88 & 0.76 & 0.47 & 0.64 & 0.88 \\
  $F_{M_{RNN_2}}$ & 0.88 & 0.65 & 0.44 & 0.39 & 0.82 \\
  $F_{M_{RNN_3}}$ & 0.88 & 0.59 & 0.21 & 0.2 & 0.75 \\
  \hline
 \end{tabular}
\end{center}
\end{table}

\begin{table}[h]
 \caption{COMPARISION RESULTS WITH RNN IN HOUSE 2}
 \label{Table 8}
 \begin{center}
  \begin{tabular}{|c|cccc|}
   \hline
   Appliance & REFR & MW & KO & ST \\
   \hline
   $F_{M_{NN}}$ & 0.85 & 0.97 & 0.91 & 0.83\\
   $F_{M_{RNN_2}}$ & 0.87 & 0.84 & 0.70 & 0.4\\
   $F_{M_{RNN_3}}$ & 0.86 & 0.70 & 0.71 & 0.21\\
   \hline
  \end{tabular}
 \end{center}
\end{table}

\begin{table}[h]
 \caption{COMPARISION RESULTS WITH RNN IN HOUSE 6}
 \label{Table 9}
 \begin{center}
  \begin{tabular}{|c|ccccc|}
   \hline
  Appliance & REFR & AC & EL & KO & ST \\
  \hline
  $F_{M_{NN}}$ & 0.80 & 0.89 & 0.70 & 1 & 0.90 \\
  $F_{M_{RNN_2}}$ & 0.80 & 0.24 & 0.13 & 0.61 & 0.47 \\
  $F_{M_{RNN_3}}$ & 0.79 & 0.37 & 0.22 & 0.66 & 0.20 \\
  \hline
  \end{tabular}
 \end{center}
\end{table}
As can be seen in Tables \ref{Table 7}, \ref{Table 8}, \ref{Table 9}, without any other preprocessing (like synthesizing data or balancing data), the results of the GRU-based RNN architecture in estimating small window aggregate data are much worse than DNN in estimating a single aggregate data in most cases. This is mainly because of the
low correlation between the aggregate sequence in operational state change. Intuitively, it's hard to estimate whether the appliance changes
its operational state according to the situation in surrounding time. And in many types of appliance, we have many more negative samples than positive samples, which also limits the performance of RNN in estimating the operational state change.

\section{Conclusion}
\label{section 5}
This paper adapts a simple DNN to detect the operational state change of appliance in the NILM task.
Comparing to previous deep learning architecture which estimates the on-off state or power consumption, the proposed
solution is low-complexity and more feasibility across variety kinds of appliance.
Through the simulated
experiments from three real houses in the REDD dataset comparing to other two benchmarks and a RNN architecture, we prove the
neural network's competitive performance in estimating the operational state change of appliance.
Unlike previous deep learning architectures which are trained on large part or even entire power data of several houses, our proposed method only needs a small percentage of one house data to achieve
a competitive performance. So the results also indicate the statistically
regularity of appliance operational state change in one house, which reveal the potential of
neural network in detecting the
operational state change in NILM task.

In the future work, we will try to use the time information in a more reasonable way to improve the performance on more appliances,  and we will also try to disaggregate several or all appliances at once.

\section*{Acknowledgment}
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 734331 and the Fundamental Research Funds for the Central Universities no. 0800219369.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\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}