\section{Introduction}
Compressed sensing is a framework that enables an $N$-dimensional sparse signal $\bm{s}=(s_i)$ to be recovered from $M(<N)$ linear measurements of its elements, $\bm{y}=\bm{Fs}$, by exploiting the prior knowledge that $\bm{s}$ contains many zero elements \cite{intro_CS}. A simple consideration guarantees that $\ell_0$-recovery,
\begin{eqnarray}
\hat{\bm{s}}=\mathop{\rm argmin}_{\bm{x}} || \bm{x} ||_0 \quad
\mbox{subj. to} \quad\bm{y}=\bm{F} \bm{x},
\label{l0}
\end{eqnarray}
where $||\bm{x}||_0$ denotes the number of non-zero elements in $\bm{x}$, is theoretically optimal in terms of minimizing the number of measurements $M$ necessary for perfectly recovering any original signal $\bm{s}$. However, carrying out $\ell_0$-recovery for a general measurement matrix $\bm{F}$ is NP-hard. To avoid such computational difficulties, an alternative approach, $\ell_1$-recovery
\begin{eqnarray}
\hat{\bm{s}}=\mathop{\rm argmin}_{\bm{x}} || \bm{x} ||_1 \quad
\mbox{subj. to} \quad \bm{y}=\bm{F} \bm{x},
\label{l1}
\end{eqnarray}
where $||\bm{x}||_1=\sum_{i=1}^N |x_i|$, is widely employed, as (\ref{l1}) is generally converted into a linear programming problem, and therefore, signal recovery is mathematically guaranteed in an $O(N^3)$ computational time through the use of the interior point method. Nevertheless, the $O(N^3)$ cost of computation can still be unacceptably high in many practical situations, and much effort is being put into finding more computationally feasible and accurate recovery schemes \cite{OMP,IRLS,Bregman,IHT,AMP,Krzakala}.

Among such efforts, the recovery scheme recently proposed by Krzakala et al.~\cite{Krzakala} is worth special attention. Their scheme basically follows the Bayesian approach. Namely, the signal recovery problem is formulated as one of statistical inference from the posterior distribution,
\begin{eqnarray}
P(\bm{x}|\bm{F},\bm{y})
=\frac{\delta(\bm{F}\bm{x}-\bm{y}) P(\bm{x})}{Z(\bm{F},\bm{y})},
\label{posterior}
\end{eqnarray}
where $Z(\bm{F},\bm{y})$ is a normalization factor imposing the condition $\int d\bm{x}P(\bm{x}|\bm{F},\bm{y})=1$ and a component-wise prior distribution $P(\bm{x})=\prod_{i=1}^N [(1-\rho) \delta(x_i)+ \rho \phi(x_i)]$ is assumed. $\rho$ and $\phi(\bm{x})$ represent the density of the non-zero signal elements and a Gaussian distribution, respectively. Exactly inferring $\bm{s}$ from (\ref{posterior}) is NP-hard, similarly to (\ref{l0}). However, by employing the belief propagation (BP) in conjunction with the expectation-maximization (EM) algorithm for estimating $\rho$ and the parameters of $\phi(x)$, they developed an approximation algorithm, termed EM-BP, which has better recovery performance than (\ref{l1}) with a computational cost of only $O(N^2)$. Furthermore, they showed that, by employing a peculiar type of ``seeded'' matrix $\bm{F}$, the threshold of the compression rate $\alpha=M/N$ of EM-BP, above which the original signal is typically recovered successfully, can approach very close to that of $\ell_0$-recovery, $\alpha_{s-EMBP}=\rho_0$, where $\rho_0$ is the actual signal density of $\bm{s}$ and $s-EMBP$ stands for `seeded EM-BP'. The seeded matrix is composed of blocks along the diagonal densely filled with Gaussian random variables. It is important to remember that this result is achieved for the first time with an approach different from $\ell_0$-recovery, being the threshold for $\ell_1$-recovery that is much higher than the optimal one: $\alpha_{\ell_1}>\rho_0$. The optimality of the $\ell_0$-recovery is guaranteed for EM-BP, which means that this scheme can practically achieve the theoretically optimal threshold of signal recovery with an $O(N^2)$ computational cost.
This remarkable property was recently proved in a mathematically rigorous manner
in the case that the matrix entries satisfy certain conditions concerning their statistics~\cite{Montanari}.
However, it is still unclear whether their scheme is optimal in terms of the computational complexity; there might be a certain design of the measurement matrix $\bm{F}$ that makes it possible to further reduce the computational cost while keeping the same signal recovery threshold.

The purpose of the present study is to explore such a possibility. For this, we focus on a class of matrices that are characterized by the following properties:
\begin{itemize}
\item {\bf sparsity:} The matrix $\bm{F}$ has only $O(1)$ non-zero elements per row and column. This implies that the measurements can be performed in a time linear in the signal length. This situation is highly preferred for the sake of practicality, given that such an operation typically needs to be done in real time, during data acquisition.
\item {\bf integer values:} The matrix elements are not real valued, but take on small integer values. This means that an optimized code for the measurements can work with bitwise operations, thus achieving much better performance without any loss of precision.
\item {\bf no block structure:} The block structure used in \cite{Krzakala} may not be necessary for reaching the optimal threshold. As an alternative possibility, we study a structure made of a square matrix in the upper left corner (the seed) plus a stripe along the diagonal. This structure is much more amenable for analytic computations, since it corresponds to a one-dimensional model homogeneous in space.
\end{itemize}

The use of sparse matrices for compressed sensing has already been suggested in several earlier studies \cite{Berinde,Gilbert,Akcakaya,KabashimaWadayama,Baron}. Our approach is particularly similar to that of \cite{Baron} in the sense that both sides are based on the Bayesian framework and use integer-valued sparse measurement matrices. Nevertheless, these two approaches differ considerably in the following two points. Firstly, we adopt EM-BP, which updates only a few variables per node for the signal recovery, while the recovery algorithm of \cite{Baron} involves functional updates and needs significantly more computational time than ours. Secondly, we thoroughly explored a simple design of $\bm{F}$ that achieves nearly optimal recovery performance. In contrast, the problem of the matrix design is not fully examined in \cite{Baron}. By carrying out extensive numerical experiments in conjunction with an analysis based on density evolution \cite{RichardsonUrbanke}, we show that a threshold close to the theoretical limit $\alpha=\rho_0$ can be achieved by using matrices with the above properties with an almost {\em linear} computational cost in the measurement and recovery stages.

This paper is organized as follows. In Section \ref{sec:EM-BP}, the EM-BP algorithm is briefly explained and the results for dense matrices are summarized. The algorithm is applied to homogeneously sparse matrices in Section \ref{sec:sparse} and to structured sparse block matrices in Section \ref{sec:block}. A new type of ``striped'' sparse matrix without blocks is introduced in Section \ref{sec:stripe}. The last Section summarizes our work, focusing on its importance for practical use, and touches on future issues.

\section{Expectation Maximization Belief Propagation\label{sec:EM-BP}}
The new algorithm based on BP in conjunction with the EM proposed in Ref. \cite{Krzakala} starts from Eq.~(\ref{posterior}).
A similar idea was also proposed in Ref. \cite{Schniter}.
In order to solve it with BP, $O(MN)$ messages for the probability distributions of the variables $x_i$ are constructed in the following way:
\begin{align*}
m_{\mu\rightarrow i}(x_i)&=\frac{1}{Z^{\mu\rightarrow i}}\int\prod_{j\neq i}dx_j m_{j\rightarrow \mu}(x_i) \delta\Big(y_{\mu}-\sum_k F_{\mu k}x_k\Big)\\
m_{i\rightarrow\mu}(x_i)&=\frac{1}{Z^{i\rightarrow\mu}}\Big[(1-\rho) \delta(x_i)+ \rho \phi(x_i)\Big] \prod_{\gamma\neq \mu} m_{\gamma\rightarrow i}(x_i)
\end{align*}
where $Z^{i\rightarrow\mu}$ and $Z^{\mu\rightarrow i}$ are normalization factors. This EM-BP equations are very complicated because the messages are distribution functions. In order to make them simpler, the messages can be approximated by assuming that they are Gaussian, thus obtaining the equations for the mean $a_{i\rightarrow\mu}$ and the variance $v_{i\rightarrow\mu}$ of $m_{i\rightarrow \mu}(x_i)$. This approximation was introduced for sparse matrices in Refs.~\cite{Kabashima2003,Guo,MontanariTse,Rangan}, and it becomes asymptotically exact if $\bm{F}$ is dense.
In fact, it is derived from an expansion in small $F_{\mu i}$, and in the dense case $F_{\mu i}=O(1/\sqrt{N})$. Supposing that the elements of the original signal follow a Bernoulli-Gaussian distribution with parameters $\rho_0$, $\overline{x}_0$ and $\sigma_0$, the update rules for the messages are the following:
\begin{align}
a_{i\rightarrow\mu}&=f_a\left(\sum_{\gamma\neq\mu}A_{\gamma\rightarrow i},\sum_{\gamma\neq\mu}B_{\gamma\rightarrow i}\right) \nonumber \\ \nonumber
a_{i}&=f_a\left(\sum_{\gamma}A_{\gamma\rightarrow i},\sum_{\gamma}B_{\gamma\rightarrow i}\right)\\ \nonumber
v_{i\rightarrow\mu}&=f_c\left(\sum_{\gamma\neq\mu}A_{\gamma\rightarrow i},\sum_{\gamma\neq\mu}B_{\gamma\rightarrow i}\right) \\ \nonumber
v_{i}&=f_c\left(\sum_{\gamma}A_{\gamma\rightarrow i},\sum_{\gamma}B_{\gamma\rightarrow i}\right)\\ \nonumber
A_{\mu\rightarrow i}&=\frac{F_{\mu i}^2}{\sum_{j\neq i}F_{\mu j}^2v_{j\rightarrow\mu}}\\
B_{\mu\rightarrow i}&=\frac{F_{\mu i}\left(y_{\mu}-\sum_{j\neq i}F_{\mu j}a_{j\rightarrow\mu} \right)}{\sum_{j\neq i}F_{\mu j}^2v_{j\rightarrow\mu}}
\label{Eq:updates}
\end{align}
where $f_a$ and $f_c$ are some analytical functions depending on the parameters $\rho$, $\overline{x}$ and $\sigma$. For details, see Ref.~\cite{Krzakala}.

In general, the original $\rho_0$, $\overline{x}_0$ and $\sigma_0$ are not known, but one can use EM to derive the update rules for them, using the property that the partition function
\begin{equation*}
Z(\rho,\overline{x},\sigma)=\int d\bm{x} P(\bm{x})\delta(\bm{y}-\bm{Fx})
\end{equation*}
is the likelihood of the parameters $(\rho,\overline{x},\sigma)$ and is maximized by the true parameters $\rho_0$, $\overline{x}_0$, and $\sigma_0$. Thus, after the update of all the messages, the inferred parameters of the original distribution are updated following these rules:
\begin{align*}
\overline{x}&\leftarrow\frac{1}{\rho N}\sum_i a_i\;, \qquad
\sigma^2\leftarrow\frac{1}{\rho N}\sum_i (v_i+a_i^2)-\overline{x}^2\;,\\
\rho&\leftarrow\frac{\sum_i \frac{1/\sigma^2+U_i}{V_i+\overline{x}/\sigma^2}a_i}{\sum_i\left[1-\rho+\frac{\rho}{\sigma(1/\sigma^2+U_i)^{\frac{1}{2}}}e^{\frac{(V_i+\overline{x}/\sigma^2)^2}{2(1/\sigma^2+U_i)}-\frac{\overline{x}^2}{2\sigma^2}}\right]^{-1}}\;,
\end{align*}
with $U_i=\sum_{\gamma}A_{\gamma\rightarrow i}$ and $V_i=\sum_{\gamma}B_{\gamma\rightarrow i}$.
If the algorithm converges to the correct solution, $a_i=s_i$ and $v_i=0$.

To reduce the number of messages from $O(NM)$ to $O(N)$, one can see that in the large $N$ limit, the messages $a_{i\rightarrow\mu}$ and $v_{i\rightarrow\mu}$ are nearly independent of $\mu$. Thus, we can derive the equations involving only a variable per each measurement node and a variable per signal node, if we are careful to keep the correcting Onsager reaction term as in the TAP equations of statistical physics \cite{TAP}. This method was introduced in the context of compressed sensing in Ref. \cite{AMP} and is called approximated message passing (AMP).

In general, the correct distribution of the original signal is unknown. However, in Ref. \cite{Krzakala}, it is demonstrated that if $\alpha>\rho_0$, the most probable configuration of $\bm{x}$ with respect to $P(\bm{x})=\prod_{i=1}^N [(1-\rho) \delta(x_i)+ \rho \phi(x_i)]$ with $\rho<1$, restricted to the subspace $\bm{y}=\bm{Fx}$, is the original signal $\bm{s}$, even if the signal is not distributed according to $P(\bm{x})$. So our choice of a Gaussian distribution for $\phi(x)$ should be perfectly fine even if the original signal has a different distribution.

The free entropy $\Phi(D)$ at a fixed mean square error $D=(1/N)\sum_{i=1}^N(x_i-s_i)^2$ can be computed if a dense matrix is used. For $\alpha>\rho_0$, the global maximum of the function $\Phi(D)$ is at $D=0$, that corresponds to the correct solution. However, below a certain threshold $\alpha<\alpha_{BP}$ that depends on the distribution $P(\bm{s})$, the free entropy develops a secondary, local maximum at $D\neq0$. As a consequence, the EM-BP algorithm can not converge to the correct solution for $\rho_0<\alpha<\alpha_{BP}$, because a dynamical transition occurs. Nonetheless, the threshold $\alpha_{BP}$ is lower than $\alpha_{\ell_1}$.

\section{EM-BP with a sparse matrix\label{sec:sparse}}

\begin{figure}[t]
\centering
\includegraphics[width=0.8\columnwidth]{sparse.eps}
\includegraphics[width=0.8\columnwidth]{DE_sparse.eps}
\caption{\label{sparse}Top: Probability of perfect recovery versus the signal sparsity $\rho_0$ using sparse matrices with $\alpha=0.5$ and $K=20$. The threshold is the same as with dense matrices. Bottom: Probability of perfect recovery computed with density evolution has the same threshold. Here, $N$ is the population size.}
\end{figure}

First of all, we want to verify if the use of a sparse matrix can reach the same results as the use of a dense one. For a sparse random matrix, the AMP equations can not be used; thus, we can use the update rules in Eq.~(\ref{Eq:updates}) for inferring the original signal $\bm{s}$. In particular, we choose the matrix $\bm{F}$ to have only $K=O(1)$ elements different from zero in each row and $H=\alpha K=O(1)$ elements in each column, extracting them from the distribution,
\begin{eqnarray}
P(F_{\mu i})=\frac{1}{2}\delta(F_{\mu i}-J)+\frac{1}{2}\delta(F_{\mu i}+J)
\label{Eq:PF}
\end{eqnarray}
with $J=1$. The use of the messages $a_{i\rightarrow\mu}$ and $v_{i\rightarrow\mu}$ instead of the AMP equations does not involve an extra cost in memory, because the number of the messages is $O(N)$ from the sparsity of the matrix. In principle, the messages $m_{i\rightarrow\mu}(x)$ are not Gaussian if the matrix is sparse, so the use of only the two parameters $a_{i\rightarrow\mu}$ and $v_{i\rightarrow\mu}$ is not exact. However, the convolution of $K$ messages (with $K=20$ in a typical matrix we use) is not far from a Gaussian, and indeed, we can verify a posteriori that this approximation is valid, because it gives good results.

In all our numerical simulations, we use a Bernoulli-Gaussian distributed signal and a compression rate $\alpha=0.5$.

Figure~\ref{sparse} (top) shows the probability of perfect recovery as a function of the sparsity of the signal $\rho_0$ for different sizes, by applying the EM-BP algorithm using a sparse matrix with $K=20$. The threshold for perfect recovery in the thermodynamic limit ($N\to\infty$) is $\rho_{BP}\simeq 0.315$, which is the same as the one obtained in Ref.~\cite{Krzakala} with a dense matrix. We can not analytically compute the free entropy $\Phi(D)$ as in \cite{Krzakala}, because we use sparse matrices and cannot use methods such as the saddle point one.
However, we performed a numerical density evolution analysis, as shown in Fig.~\ref{sparse} (bottom), and found that the threshold is almost the same as the one computed with the matrices $\bm{F}$.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{local_stability_sparse.eps}
\caption{\label{stability}Stability check of the solution determined by the EM-BP message passing algorithm. Starting the recovery process with a sparse matrix from an initial condition differing less than $\Delta$ from the correct solution, the latter is recovered as long as $\rho_0 < \alpha_{stab}(\Delta)$. In the limit $\Delta\to 0$, the stability limit $\alpha_{stab}(\Delta)$ tends to the theoretical bound $\alpha$ (which is 0.5).}
\end{figure}

Next, we will verify that the correct solution is always the global maximum of $\Phi(D)$ and it is locally stable up to $\alpha \simeq \rho_0$ when using EM-BP with a sparse matrix. Since we can not analytically compute the free entropy, we must resort to a numerical method. We start EM-BP with an initial condition very close to the correct solution: $a_{i\rightarrow\mu}^0=s_i+\delta_{i\rightarrow\mu}$, with $\delta_{i\rightarrow\mu}$ a random number uniformly distributed in $[-\Delta,\Delta]$. In this way, we have verified (see Fig.~\ref{stability}) that if $\Delta$ is sufficiently small, the correct solution can be found up to $\alpha \simeq \rho_0$, as in the case of a dense matrix.

For the algorithms based on the $\ell_1$ minimization, it is known that the threshold with a sparse matrix is lower than that with a dense one. However, these algorithms are not optimal, because the correct solution disappears below the threshold $\alpha_{\ell_1}$. In this sense, the EM-BP algorithm is optimal, because the global maximum of the free entropy is always on the correct solution. Thus, one can expect that, if the rank of the sparse matrix is the same as that of the dense one, a similar threshold can be reached (as we have demonstrated numerically).

In summary, we can say that the EM-BP algorithm of Ref. \cite{Krzakala} seems to reach the same threshold $\alpha_{BP}$, either using a dense Gaussian matrix or a sparse binary one
if the numbers of non-zero elements per row/column are $O(1)$ but sufficiently large.
However, the use of a sparse matrix is computationally much faster than the use of a dense one. Moreover, the use of binary elements, instead of Gaussian real values, allows for better code optimization and eventually for hard-wired encoding of the compression process.

\section{Block-structured sparse matrices\label{sec:block}}

\begin{figure}[t]
\includegraphics[width=\columnwidth]{structured.eps}
\caption{\label{structured}$\rho_c(N)$ with a sparse, structured matrix with blocks, for different sizes $N$ and values of $L$ (see text). The thresholds for the $\ell_1$-recovery and for EM-BP without any structure are also drawn. For comparison with the data in Fig.~\ref{unid}, we used the exponent $a \simeq 0.18$ best fitting those data.}
\end{figure}

To avoid the secondary maximum of the free entropy, in Ref. \cite{Krzakala}, the authors use a structured block matrix that helps to nucleate the correct solution. The idea is that the correct solution is found for the first variables, and then it propagates to the whole signal. This idea is similar to the so-called spatial coupling that is very useful for solving many different problems \cite{spat-coup}. With this trick, the authors of Ref.~\cite{Krzakala} reach perfect recovery for almost any $\alpha > \rho_0$ in the large $N$ limit while Ref.~\cite{Kudekar} reports that the gain is quite small when different recovery algorithms are used. Here, we try to use a matrix with the same block structure, but sparsely filled. We divide the $N$ variables into $L$ groups of size $N/L$ and $M$ measurements into $L$ groups of size $M_p=\alpha_pN/L$ in such a way that $M=\sum_{p=1}^L M_p=\alpha N$ and $1/L\sum_{p=1}^L\alpha_p=\alpha$.

In this way, the matrix ${\bm F}$ is divided into $L^2$ blocks, labeled with indices $(p,q)$. Each block is a sparse binary matrix with $k$ elements different from zero for each row and $h_p=\alpha_p k$ elements for each column, distributed according to Eq.~(\ref{Eq:PF}), with $J=J_{p,q}$. As in Ref. \cite{Krzakala}, we choose $J_{p,p-1}=J_1$, $J_{p,p}=1$, $J_{p,p+1}=J_2$, and $J_{p,q}=0$ otherwise. The important ingredient to nucleate the correct solution is that in the first block $\alpha_1=(M_1/N)L>\alpha_{BP}$ holds. For simplicity, we can choose $\alpha_1=1$ and $\alpha_p=(L\alpha-1)/(L-1)$ for $p\neq1$.
The recovery strongly depends on the parameters $J_1$ and $J_2$, and the best results for $\alpha=0.5$ are obtained around $J_1=4$ and $J_2=1$. Moreover, we used these two values in the experiments described below because we wanted to work with matrices with elements having small integer values.

Similarly to the dense case, the use of a sparse structured matrix with blocks allows to overcome the dynamical transition at $\alpha_{BP}$ and to nucleate the correct solution until $\alpha$ is very close to $\rho_0$. Figure~\ref{structured} shows the mean critical threshold $\rho_c(N)$ for different signal lengths at a fixed compression rate $\alpha=0.5$. The $x$ axis uses the same scaling variable as in Fig.~\ref{unid}, and the best parameter $a$ obtained from the fit of data in Fig.~\ref{unid} also interpolates the data in Fig.~\ref{structured} quite well. In the thermodynamic limit, $\rho_c$ extrapolates to a value compatible with the optimal one, $\alpha$, and it is certainly much higher than the thresholds for $\ell_1$-recovery and for EM-BP without any structure. We have also done a density evolution analysis that confirms this result.

For each value of $N$ and $L$, the mean critical threshold $\rho_c$ is computed as follows. We randomly generate a block structured matrix ${\bm F}$ with the given $N$ and $L$. We start with a sufficiently sparse original signal $\bm{s}$, which has been recovered by the algorithm; we then add non-zero entries to the signal and check whether the new signal can be recovered by the algorithm; we go on adding non-zero elements to the signal until a failure in a perfect recovery occurs. The previous to the last value for $\rho_0$ is the critical threshold for the matrix ${\bm F}$. The mean critical threshold is obtained by averaging over many different random matrices and signals, with the same values of $N$ and $L$. The number of such random extractions goes from $10^3$ for the largest $N$ value up to $10^4$ for the smallest $N$ value.

The values of $(N,L,k)$ used for the simulations shown in Fig.~\ref{structured} are the following: $(2250,10,9)$, $(19000,20,19)$, $(31200,40,39)$, and $(49000,50,49)$. We need to increase both $N$ and $L$ if we want to obtain good results in the thermodynamic limit. However, if we change $L$, we must change $k$ too. Indeed, in order to have the same number of elements per row and column for each of the $L^2$ blocks, we must satisfy the conditions: $(N/L)/M_p=k/h_p$ with $k$ and $h_p$ integer valued. Here, we have used the smaller possible value for $k$, that is $k=L-1$. The fact that it is impossible to keep $k$ constant while increasing $L$ implies that these kinds of block-structured matrices always become dense in the thermodynamic limit. This is a limitation of the block structure that we want to eliminate with the matrix proposed in the following Section.

\section{And without blocks?\label{sec:stripe}}
The matrix proposed in Ref. \cite{Krzakala} is not the only one that allows the optimal threshold to be reached. Reference \cite{Krzakala2} analyzes the use of other good dense, block-structured matrices. However, the block-structure is not so simple to handle if one wants to do analytical calculations in the continuum limit. Moreover, in making these block-structured matrices sparse, one has to be careful to find the right values of $L, M, N, M_p, \alpha_p, k$. For these reasons, we want to know if the block structure is crucial, and, if not, we want to eliminate it.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{stripeMatr.eps}
\caption{\label{matrixunid}A nearly one-dimensional sparse matrix with a square block of size $L \times L$ at the top left and non-zero elements in the stripes around the diagonal can achieve compression and perfect recovery close to the theoretical bound in linear time.}
\end{figure}

We tried a different structured sparse matrix (see Fig.~\ref{matrixunid}), that we called a striped matrix. It has one sparse square block of size $L$ on the top left of the matrix with $K=O(1)$ elements for each row and column extracted from (\ref{Eq:PF}) with $J=1$. This arrangement is fundamental for nucleating the correct solution. Apart from this first block, the residual compression rate is $\alpha'=\frac{M-L}{N-L}$. Then, we construct a one-dimensional structure around the diagonal of the remaining matrix. For each column $c>L$, we randomly place $2K\alpha'$ non-zero elements, again extracted from (\ref{Eq:PF}), in the interval of the width $2L\alpha'$ around the diagonal. One element with $J=1$ is always placed on the diagonal (actually on the position closest to the diagonal). For the remaining elements, we use the following rules. If the element is at a distance $d\leq L\alpha'/3$ from the diagonal, we use $J= 1$. Otherwise, if its distance is $d>L\alpha'/3$, we use $J=J_1$ below the diagonal and $J=J_2$ above the diagonal. In this way, the number of elements per column is constant, while the number of elements per row is a truncated Poisson random variable with mean $2K$: indeed, there are no empty rows, thanks to the rule of placing the first element of each column closest to the diagonal. When constructing the matrix, we apply exactly the same rule to each column, but in the last $L$ columns it may happen that a non-zero element has a row index larger than $M$: these elements are then moved below the first square matrix by changing the row and column indices as follows: $r \leftarrow r - (M-L)$ and $c \leftarrow c - (N-L)$.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{unid.eps}
\caption{\label{unid}$\rho_c(N)$ with a sparse, striped matrix, as described in Section \ref{sec:stripe} for different sizes (from $N=2000$ to $N=40000$). The thresholds for $\ell_1$-recovery and for EM-BP without structure are drawn for comparison. The best fitting parameter is $a \simeq 0.18$, and it leads to an extrapolation of $\rho_c(N)$ in the thermodynamic limit compatible with $\alpha$.}
\end{figure}

In this way, we have some kind of continuous one-dimensional version of the block-structured matrix discussed in the previous Section. Within this striped matrix ensemble, the thermodynamic limit at a fixed matrix sparsity can be calculated without any problem, by sending $N,L \to \infty$ at a fixed $L/N$ and fixed $K=O(1)$. In Fig.~\ref{unid}, we show the mean critical threshold reached by using striped matrices with a fixed ratio $L/N=1/50$ (the same used in the plot of Fig.~\ref{matrixunid}) and different signal lengths. Perfect decoding up to $\rho_c$ is again achieved by using the EM-BP algorithm. We extrapolated the $\rho_c(N)$ data to the thermodynamic limit by assuming the following behavior in the large $N$ limit:
\begin{equation}
\rho_c(N) = \rho_c(\infty) - b N^{-a}
\end{equation}
The data in Fig.~\ref{unid} are plotted with the best fitting parameter $a \simeq 0.18$, and the extrapolated value $\rho_c(\infty)$ is perfectly compatible with the theoretical bound $\alpha$.

Hence, we can conclude that the important ingredient to reach optimality is not the block structure, but the nearly one-dimensional structure, associated with the initial block with $\alpha_1 > \alpha_{BP}$ to nucleate the correct solution.

It is worth noticing that the corresponding statistical mechanics model for these striped random matrices is a one-dimensional disordered model with an interaction range growing with the signal length, as in a Kac construction. Models of this kind are analytically solvable and have shown very interesting results \cite{SilvioAndrea}.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{tempi.eps}
\caption{\label{Fig:time}Actual time (in seconds) for recovery of a signal with dense and sparse matrices for different data lengths $N$. The data are fitted respectively by a quadratic and a linear function.}
\end{figure}

The use of our striped sparse matrices allows for a great reduction in computational complexity. Indeed, the measurement and recovery times grow linearly with the size of the signal if sparse matrices are used, while they grow quadratically if dense matrices are used. Figure \ref{Fig:time} shows the measurement and recovery times of a signal for different signal lengths $N$. For this test, we used dense block-structured matrices and sparse striped matrices. The number of EM-BP iterations to reach the solution is roughly constant for different $N$. A quadratic fit for the dense case and a linear fit for the sparse one perfectly interpolate the data.

\section{Conclusions and future development\label{sec:conclusions}}
We introduced an ensemble of sparse random matrices ${\bm F}$ that, thanks to their particular structure (see Figure \ref{matrixunid}), allow us to perform the following operations in linear time:\\
(i) measurement of a $\rho_0$-sparse vector ${\bm s}$ of length $N$ by using a linear transformation, ${\bm y} = \bm{Fs}$, to a vector ${\bm y}$ of length $\alpha N$ \\
(ii) perfect recovery of the original vector ${\bm s}$ by using a message passing algorithm (Expectation Maximization Belief Propagation) for almost any parameter satisfying the theoretical bound $\rho_0 < \alpha$.

These striped sparse matrices that have such good performance because there is a `seeding' sparse square matrix in the upper left corner that nucleates a seed for the right solution and the one-dimensional structure along the diagonal propagates the initial seed to the complete right solution. Both seeding and a one-dimensional structure have been used in the past \cite{spat-coup,Krzakala}, but in our new ensemble, the matrices are sparse, and this permits us to perform all the operations in a time linear in the signal length.

We also checked that sparse matrices perform as well as dense ones in the case of block-structured matrices and for matrices with no structure at all.

Apart from the compressed sensing case, several other applications require a sparse matrix or equivalently linear time complexity \cite{Gilbert}

In data streaming computing, one is typically interested in doing very quick measurements in constant time. For example, if the task is to measure the number of packets $s_i$ with destination $i$ passing through a network router, it is not possible to keep a vector $\bm{s}$ because it is generally too long. Instead, a much shorter sketch of it, $\bm{y}=\bm{Fs}$, is measured in such  a way that a very sparse vector $\bm{s}$ can be recovered from $\bm{y}$. The matrix $\bm{F}$ must be sparse in order to be able to update the sketch $\bm{y}$ in a constant time for each new packet passing through the router.

Another interesting application is the problem of group testing, where a very sparse vector $\bm{s}\in\{0,1\}^N$ is given and one is interested in performing the fewest linear measurements, $\bm{y}=\bm{Fs}$, that allow for detection of the defective elements ($s_i=1$). In this case, the experimental constraints require a sparse matrix $\bm{F}$: only if the tested compound $y_\mu = \sum_i F_{\mu i} s_i$ is made of a very few elements of $\bm{s}$, the linear response holds and non-linear effects can be ignored.

However, in the more general case, one does not directly observe the sparse signal $\bm{s}$ but rather a linear transformation of it, $\bm{x}=\bm{Ds}$, made with a dictionary matrix $\bm{D}$ (which is typically a Fourier or wavelet transformation, and thus is a dense matrix). In this more difficult case, one would like to design a sparse measurement matrix $\bm{A}$ such that the measurement/compression operation, $\bm{y}=\bm{Ax}$, is fast, and the resulting observed data $\bm{y}$ is short, thanks to the sparseness of $\bm{s}$. The conflicting requirement is to have a fast recovery scheme, because, now, to recover the original signal one should solve $\hat{\bm{s}}= {\rm argmin}\,||\bm{s}||_0$ subject to $\bm{y}=(\bm{AD})\bm{s}$, where $\bm{AD}$ is typically dense (e.g.\ in case of Fourier and wavelet transformations). So a very interesting future development of the present approach is to extend it to this more complex case.

\section*{Acknowledgments}
FR-T is grateful for his useful discussions with F.~Krzakala, M.~M\'ezard and L.~Zdeborova and financial support from the Italian Research Minister through the FIRB Project No. RBFR086NN1 on ``Inference and Optimization in Complex Systems: From the Thermodynamics of Spin Glasses to Message Passing Algorithms''.

YK is supported by grants from the Japan Society for the Promotion of Science (KAKENHI No. 22300003) and the Mitsubishi foundation.

\IEEEtriggeratref{16}


\section{Introduction.}
Low-temperature transport measurements have revealed level crossings in Coulomb blockaded quantum dots and single-molecule transistors, induced either by tuning  gate voltage\cite{Kogan03, Hauptmann08, Roch08a, Osorio10} or external magnetic field\cite{Nygard00, Fasth07, Katsaros10, Jespersen11}. These show up in a bias spectroscopy as inelastic cotunneling thresholds crossing zero, when tuning the control parameter towards the level crossing. If the two different ground states have the same elastic cotunneling amplitudes, the linear conductance shows hardly any trace of the crossing, except perhaps for a peak right at the crossing if the degeneracy gives rise to Kondo-effect~\cite{Hauptmann08, Jespersen11}. In the case of a transition from a singlet to a triplet ground state, however, the elastic cotunneling amplitude can easily be different on the two sides, and the triplet state can even give rise to Kondo effect~\cite{Pustilnik01prb, Kogan03, Roch08a, Osorio10}. The difference in elastic cotunneling leads to a sudden change in the linear conductance at the crossing, which may be smeared by a small mixing of the two competing ground states. In the case of a singlet-triplet crossing, a mixing due to spin-orbit coupling leads to an avoided crossing, which can be observed as a lower bound on the inelastic cotunneling threshold~\cite{Fasth07, Katsaros10, Jespersen11}.

Two recent experiments on widely different single-molecule transistors have shown gate-induced singlet-triplet crossings. These experiments revealed an extra cotunneling line in the differential conductance pinned to the crossing at zero bias and extending up to bias voltages much larger than the other excited spin states, thus demarcating the boundary between two different regions~\cite{Roch08a, Osorio10}. In (the supplementary material to) Ref.~\onlinecite{Osorio10}, see Figs. 6S and 7S, reproduced here in Fig.~\ref{fig:FigExpCuts}, it was pointed out that such a feature could easily arise if the device has sufficiently asymmetric coupling to source and drain. In this  case, even at finite bias voltage the molecule is largely in equilibrium with the strongest coupled electrode. Like at zero bias, a difference in the elastic cotunneling amplitude may therefore show up as a step edge extending into the finite bias region.
\begin{figure}[t]
\includegraphics[width=0.35\textwidth]{LevelCross.pdf}
\caption{\label{fig:LevCross}Sketch of differential cotunneling conductance as a function of bias voltage, $V$, and control parameter, $x$. Thick blue (solid) lines indicate inelastic cotunneling thresholds with an avoided crossing near zero bias, caused by a finite mixing, $\lambda$, of the two involved states. The red (dashed) line indicates the line along which the levels cross (possibly resulting in a conductance peak), including the possibility of a slight bias dependence of the level splitting. The two different sides of this line may have a difference in $dI/dV$ (indicated by different shades of blue) due to a difference in elastic cotunneling amplitudes for the two states.}
\end{figure}
This scenario is sketched in Fig.~\ref{fig:LevCross}, with differential conductance as a function of bias voltage $V$ and control parameter $x$, which in these two experiments is gate voltage. In both experiments this line has a finite slope, implying that the bias voltage somehow influences the tuning of the crossing as part of the total effective gating of the relevant states. In this case, the step edge is accompanied by a peak deriving from differentiating a slanted step in current. In Fig.~\ref{fig:LevCross} the difference in conductance is reflected by two different shades of blue, and a conductance peak may be observed along the border marked by a red (dashed) line. This behavior can indeed be seen from Fig. S6a in Ref.~\onlinecite{Roch08a}, but the corresponding feature observed in Ref.~\onlinecite{Osorio10} revealed a surprising dip-peak Fano-like structure, which remained unexplained.

In this paper, we show that the dip-peak lineshape observed in Ref.~\onlinecite{Osorio10} can arise from a competition between nonequilibrium effects and spin-orbit mixing. In general, a mixing, $\lambda$, of the two states can lead to an avoided crossing as indicated in Fig.~\ref{fig:LevCross}. However, in the experiment \cite{Osorio10} the spin-orbit coupling is weak enough that no avoided crossing is observed in the inelastic cotunneling, i.e., it is smaller than the tunnel-induced broadening of the cotunneling lines. Since the line is only weakly bias dependent (nearly vertical), its detailed line shape still provides otherwise inaccessible information about the spin-orbit coupling. We shall restrict our attention to the singlet-triplet crossing observed in Ref.~\onlinecite{Osorio10}, but the basic mechanism is of general validity and this phenomenon could be observed in a variety of cotunnel junctions.

The paper is organized as follows. In Sect.~\ref{sec:model} we introduce the double-dot model which will be considered throughout the paper. In Sect.~\ref{sec:rateeq} we present the results of a perturbative generalized master equation calculation of transport as a function of gate and bias voltage. This calculation reveals the singlet-triplet boundary in finite-bias cotunneling, and shows that it should connect to kinks in the Coulomb diamond edges arising from sequential tunneling. In Sect.~\ref{sec:kondo}, we focus on the interior of the Coulomb diamond and recalculate the inelastic cotunneling spectrum within an effective exchange cotunneling model. The effective cotunneling model is first established from a Schrieffer-Wolff transformation of the original double-dot model, and the features of Sect.~\ref{sec:rateeq} are reproduced. Finally, we include a singlet-triplet mixing in the form of a spin-orbit induced Dzyaloshinskii-Moriya term, and demonstrate how this causes a nonequilibrium population inversion which in turn gives rise to the peak-dip structure observed experimentally in Ref.~\onlinecite{Osorio10}.

\section{Double dot model}\label{sec:model}

As a simple model which can undergo a transition between singlet and triplet groundstates, we consider an interacting double quantum dot, where the two dots are connected in parallel to a left ($L$) and a right ($R$) electrode and capacitively coupled to a gate electrode, see Fig.~\ref{fig:DDmodel}. The Hamiltonian is given by $H = H_{DD} + H_R + H_T$, where
\begin{align}
\label{eq:HD}
	H_{DD} &= \sum_{i \sigma} \epsilon_{i \sigma} n_{i \sigma} + \sum_{i} U_i n_{i \uparrow} n_{i \downarrow} + U_{12} n_1 n_2 \nonumber \\
	    &+ J_{12} \mathbf{S}_1 \cdot \mathbf{S}_2 + B \left( S_1^z + S_2^z \right), \\
\label{eq:HR}
	H_R &= \sum_{r k \sigma} \epsilon_{r k} n_{r k \sigma}, \\
\label{eq:HT}
	H_T &= \sum_{r k i \sigma} t_{r i} c_{r k \sigma} d_{i \sigma}^\dagger + h.c.
\end{align}
The double dot system is described by $H_{DD}$, where $\epsilon_{i \sigma}$ is the energy of the single-particle orbital on dot $i$, which has number operator $n_i = \sum_\sigma n_{i \sigma} = \sum_{\sigma} c_{i \sigma}^\dagger c_{i \sigma}$ for spin projection $\sigma = \uparrow,\downarrow$. $U_i$ and $U_{12}$ are the onsite and inter-dot Coulomb charging energies, respectively. $\mathbf{S}_i$ is the spin on dot $i$ (where $\mathbf{S}_i = \mathbf{0}$ if dot $i$ is empty or doubly occupied). There is an interaction between the spins on the different dots, which can be dominated by either superexchange ($J_{12} > 0$) or Hund's rule coupling ($J_{12} < 0$), and arises as an effective low-energy description of a more complicated multi-orbital double dot model, where orbitals on different dots are partially overlapping~\cite{Osorio10}. $B$ is the magnetic field, which we assume to be applied in the $z$-direction and therefore coupled to the $z$-projection of the dot spins, $S_i^z$.

\begin{figure}
\includegraphics[width=0.4\textwidth]{DDmodel.pdf}
\caption{\label{fig:DDmodel}(Color online) Sketch of double dot model. Electrons on different dots interact via their charge (Coulomb repulsion, $U_{12}$) as well as via their spin (exchange, $J_{12}$). Tunneling between electrode $r$ and  dot $i$ takes place with amplitude $t_{ir}$. The gate electrode is not shown.}
\end{figure}
$H_R$ describes non-interacting electrons in reservoir $r = L,R$, with energy $\epsilon_{r k}$ and number operator $n_{r k \sigma} = c_{r k \sigma}^\dagger c_{r k \sigma}$. Tunneling between the dots and the electrodes is described by $H_T$. For simplicity we assume that the tunnel amplitudes, $t_{r i}$, are independent of $k$ and that the density of states in the electrodes, $\rho_r$, is energy independent. Then also the tunnel rates, $\Gamma_{r i} = 2\pi\rho_r |t_{r i}|^2$, which set the inverse time-scale for resonant single electron tunneling (SET), are energy-independent.

\section{Full gate and bias spectroscopy}\label{sec:rateeq}

We calculate the nonequilibrium reduced density matrix of the double quantum dot, and the current flowing between the electrodes, $I$, using a generalized master equation (GME) approach~\cite{Leijnse08a, Koller10}, based on the real-time transport theory~\cite{Koenig96a}. This approach treats all the local interactions ($U_i, U_{12}, J_{12}$) non-perturbatively. The tunneling between the dots and electrodes is treated perturbatively up to fourth order in $H_T$, which includes all coherent one- and two-electron tunnel processes, such as SET, elastic and inelastic cotunneling~\cite{DeFranceschi01}, pairtunneling~\cite{Koch05c, Leijnse09a}, and tunnel-induced level broadening and shift. The Kondo effect is not included at this level of perturbation theory, which is therefore applicable only at temperatures much larger than the Kondo temperature. In fact, for the theory to be applicable also in the SET regime, the thermal energy scale must dominate over the tunnel broadening, $T \gg \Gamma_{ri}$.

The results of the GME calculations are shown in Fig.~\ref{fig:diamonds}, represented as two-dimensional conductance maps, with the absolute value of the differential conductance, $|dI/dV|$, plotted on a color scale as a function of bias voltage, $V$, and gate voltage, $V_g$.
Because of capacitive effects, the energy difference between states with $N = \sum_i n_i$ and $N-1$ electrons on the dot is $\propto - N V_g$, where we have set the proportionality constant (gate coupling) to one for simplicity. We assume equal capacitive couplings to left and right electrodes and apply the bias voltage symmetrically with chemical potentials $\mu_{L,R} = \pm V/2$, in which case the bias voltage does not result in direct capacitive energy-shifts. However, as was shown in Ref.~\onlinecite{Osorio10}, $J_{12}$ may depend on both $V$ and $V_g$ (for simplicity we assume a linear dependence).
Therefore the $N=2$ groundstate will change from singlet to triplet along a line in the $(V,V_g)$ plane.

Figure~\ref{fig:diamonds}(a) shows the conductance within a voltage range where $N = 2$ in equilibrium.
\begin{figure}
\includegraphics[width=0.5\textwidth]{Fig2ab.pdf}
\caption{\label{fig:diamonds}(Color online) Results of GME calculation presented as conductance maps, $|dI/dV|$ plotted on a logarithmic color scale as a function of $V$ and $V_g$ (the absolute value is plotted since negative values, occuring in certain regions because of the tunnel-rate asymmetry, are not well represented on a logarithmic scale). The red arrows point out the singlet-triplet line. Here the tunnel amplitudes are $t_{1R} = 30, t_{1L} = 0.7, t_{2R} = 0.8, t_{2L} = 0.5$ (arbitrary units), the Coulomb charging energies are $U_{11} = U_{22} = 2000 T$, $U_{12} = 400 T$, and we choose $\epsilon_2 = \epsilon_1 + 2T$. $B = 0$ in (a) and $B = 50 T$ in (b). }
\end{figure}
Already the rather simple model Hamiltonian~(\ref{eq:HD}) has a complicated many-body eigenspectrum and gives rise to a rich conductance map, but we focus here on the features originating from the transition between singlet and triplet ground states. The strong (yellow) conductance lines with a large gate dependence mark the onset of SET. The crossing of such lines at zero bias voltage indicates a charge-degeneracy point.
The lines emmanating from the zero-bias crossing points at the left (right) indicates the onset of SET involving the $N=1$ ($N=3$) charge states.
Inside the central region (Coulomb diamond), the charge is fixed to $N=2$, and the current is dominated by coherent two-electron cotunneling processes, where an electron is effectively transferred from one electrode to the other, involving only virtual occupation of either the $N=1$ or $N=3$ charge states. A cotunneling process which does not transfer any net energy to or from the double dot is called elastic and gives rise to a constant background conductance. In an inelastic cotunneling process~\cite{DeFranceschi01}, in contrast, energy is transferred to the double dot, which becomes possible when that energy can be supplied by the bias voltage. This gives rise to conductance features whenever the bias voltage equals an excitation energy of the double dot system. In first approximation the conductance feature is a step, but the joint effects of nonequilibrium population of excited states and Kondo correlations may give rise to a more complicated lineshape~\cite{Paaske06}.

The inelastic cotunneling step corresponding to excitations from the singlet ground state to the triplet excited states emerges from the left diamond edges in Fig.~\ref{fig:diamonds}(a) (seen as an abrupt color change from white to blue). Because of the gate dependence of the singlet-triplet splitting (due to the gate dependence of $J_{12}$), the inelastic cotunneling step acquires a finite slope and crosses zero bias close to the center of the Coulomb diamond. At this point the double dot ground state switches from singlet to triplet. To the right of this crossing a cotunneling step again emerges, corresponding to the opposite excitation (from the triplet, which is now the ground state, to the singlet). This step ends when it intersects the right diamond edge. Also emerging from the zero-bias crossing is a weak conductance peak (for $V>0$) or conductance dip (for $V<0$), marked with red arrows, which signals the finite-bias position of the singlet-triplet transition (the slope is a result of the combined $V$ and $V_g$ dependencies of $J_{12}$). This "singlet-triplet line", which was experimentally observed in Refs.~\onlinecite{Osorio10, Roch08a}, is thus a signature of the ground state transition as a function of the applied voltages and appears since the elastic cotunneling conductance of the triplet is somewhat larger than that of the singlet (the lineshape, and the conditions under which it is visible, are discussed in detail in Sect.~\ref{sec:kondo}).
At the points where the singlet-triplet line intersects the edge of the Coulomb diamond, there is a discontinuous change in the slope of the diamond edges. The reason is that when crossing this line, the molecular ground state is changed into one with a different voltage dependence, which directly determines the diamond slopes.

Figure~\ref{fig:diamonds}(b) shows the same conductance map, but under the influence of a magnetic field. Those SET resonances which correspond to excitations with a finite total spin are now split by the Zeeman effect. In addition, the Zeeman splitting of the triplet state results in a splitting of the corresponding inelastic cotunneling steps. When the ground state is a singlet, excitations to the three triplet states lead to three distinct gate-dependent cotunneling steps. To the right of the ground state crossing point, where the ground state is the $T_{-1}$ triplet, excitations to the singlet give rise to a gate-dependent cotunneling step and excitations to the $T_{0}$ triplet result in a gate-independent step (excitations from the $T_{-1}$ to the $T_{+1}$ triplet are suppressed by spin-selection rules for cotunneling, $\Delta S_{z} = 0, \pm 1$). In an experiment, the Zeeman splitting of the inelastic cotunneling steps provides an important check of the spin assignment of the excitations.

In both Figs.~\ref{fig:diamonds}(a) and (b), weak conductance peaks with the same gate dependence as the Coulomb diamonds are seen inside the Coulomb blockade region, close to the diamond edges. These peaks result from SET transitions starting from excited states~\cite{Schleser05}, which can be populated in a nonequilibrium situation due to inelastic cotunneling processes.

\section{Lineshape of the singlet-triplet transition}\label{sec:kondo}

In Fig.~\ref{fig:diamonds}, the singlet-triplet line appeared as a peak or a dip in the conductance, for respectively positive or negative bias voltage. This was also the result of the calculations presented in the supplementary information to Ref.~\onlinecite{Osorio10}, but as pointed out already there, the experimental data from Ref.~\onlinecite{Osorio10}, reproduced here in Fig.~\ref{fig:FigExpCuts}, do not adhere to this behavior. Instead, the experiment shows a conductance peak, for both positive and negative bias voltage, which then evolves into a peak-dip, for $V>0$, or dip-peak, for $V<0$, roughly when the magnitude of the bias voltage becomes larger than all the relevant inelastic cotunneling thresholds. In this Sect., we will show that these Fano-like lineshapes along the singlet-triplet boundary can be reproduced by introducing a small mixing of the competing singlet and triplet states. In contrast to the scenario in Fig.~\ref{fig:diamonds} and the one discussed in the supplementary information to Ref.~\onlinecite{Osorio10}, this effect relies entirely on nonequilibrium pumping of excited states and therefore shows up only in devices with not too different couplings to source and drain electrodes.

\subsection{Effective exchange cotunneling model}

\begin{figure}
\begin{center}
\subfigure[]
{
\centering
\label{fig:FigExpCuts}
\includegraphics[width=1.0\columnwidth]{FigExpCuts.pdf}
}
\subfigure[]
{
\centering
\label{fig:DensPlot}
\includegraphics[width=1.0\columnwidth]{DensPlot62.pdf}
}
\caption{(a) The lower density plot shows $d^{3}I/dV^{3}$ vs. bias voltage, $V$, and gate voltage, $V_{g}$, for the Mn complex measured in Ref.~\onlinecite{Osorio10}. An external magnetic field of 10 Tesla clearly reveals a singlet groundstate with an excited triplet towards the left, and a triplet groundstate with an excited singlet towards the right. The zoom-in shows the same data now plotted as $dI/dV$. A steep faint line is seen to be pinned to the groundstate crossing. Cuts at various bias voltages, marked by dashed lines, reveal a distinct lineshape as a function of $V_g$ developing at large applied bias: {\it dip-peak} at $V<0$ and {\it peak-dip} at $V>0$. The two nearly vertical dashed lines indicate the range in $V_g$ used for each of the cuts. The vertical range is the same in all the cuts, and has been centered at the mean value within each cut. (b) Density plot of differential conductance as calculated below in Sect.~\ref{sec:STmix}, with cuts at bias voltages $V/T=\pm269,\pm76,0$ (parameters are chosen as in Fig.~\ref{fig:symcutsV}). The peak-dip (dip-peak) structure from panel (a) is clearly reproduced. All density plots are on a grey scale where darker indicates higher values.}
\end{center}
\end{figure}

We now focus our attention to the region deep inside the Coulomb diamond with $N=2$, where real charge fluctuations (SET) are suppressed, and second order processes (cotunneling), involving virtual fluctuations to the nearby $N=3$ and $N=1$ charge states, dominate. If $U_i \gg U_{12}$, there is always one electron in each dot when $N = 2$. In this limit, the system described by the Hamiltonian~(\ref{eq:HD}) is effectively a two-spin system, and it is natural to work in the singlet-triplet basis $\{| T_{+1}\rangle,| T_{0}\rangle,| T_{-1}\rangle,| S\rangle\}$. In this basis, the {\it double-dot spin} Hamiltonian reads:
\begin{equation}\label{eq:HD2}
 H_{DDS}=\sum_{\eta}|\eta\rangle\varepsilon_{\eta}\langle\eta|,
\end{equation}
with eigenenergies
\begin{align}
\varepsilon_{S}=&-\frac{3}{4}J_{12},\\
\varepsilon_{T_{-1}}=&-B+\frac{1}{4}J_{12},\\
\varepsilon_{T_{0}}=&\frac{1}{4}J_{12},\\
\varepsilon_{T_{1}}=&B+\frac{1}{4}J_{12},
\end{align}
where we have omitted the constant term $U_{12} + \epsilon_1 + \epsilon_2$. As discussed above, $J_{12}$ is assumed to depend linearly on $V$ and $V_g$, and taking the singlet-triplet crossing point to be the zero of
our gate voltage axis, we thus assume an exchange coupling of the form:
\begin{align}\label{eq:Jfct}
	J_{12}=-\alpha V_g - \beta V,
\end{align}
in terms of two dimensionless numbers, $\alpha$ and $\beta$.

The kinetic energy of the electrode-electrons is still described by $H_R$ and, applying a Shrieffer-Wolff transformation~\cite{schriefferWolff66} to effectively eliminate the charge fluctuations, we arrive at the following exchange cotunneling Hamiltonian:
\begin{align}\label{eq:Hcot}
H_{cot}&=\!\!\sum_{\substack{r,r'=R,L\\i=1,2\\\sigma\sigma'=
\uparrow,\downarrow\\k,k'}}
\!\left(
J_{i,rr'}\bm{S}_{i}\cdot\bm{\tau}_{\sigma\sigma'}
+\frac{W_{i,rr'}}{2}\delta_{\sigma \sigma'}\right)
c^{\dag}_{r k\sigma}c_{r' k'\sigma'}.
\end{align}
where $\bm{\tau}$ is the vector of Pauli matrices. Written in the singlet-triplet basis, the local spin-operators,
\begin{align}
\bm{S}_{i} =\left((S_{i}^+ + S_{i}^-)/2, -i (S_{i}^+ - S_{i}^-)/2, S_{i}^z\right),
\end{align}
take the following form:
\begin{align}
S^{+}_{i}&=\frac{1}{\sqrt{2}}\left(
 \mid T_{+1}\rangle\langle T_{0}\mid
+\mid T_{0}\rangle\langle T_{-1}\mid
\right.\label{eq:Sop1}\\&\left.\hspace*{14mm}
\pm\mid T_{+1}\rangle\langle S\mid
\mp\mid S\rangle\langle T_{-1}\mid
\right),\nonumber\\
S^{-}_{i}&=\frac{1}{\sqrt{2}}\left(
\mid T_{0}\rangle\langle T_{+1}\mid
+\mid T_{-1}\rangle\langle T_{0}\mid
\right.\\&\left.\hspace*{14mm}
\pm\mid S\rangle\langle T_{+1}\mid
\mp\mid T_{-1}\rangle\langle S\mid
\right),\nonumber\\
S^{z}_{i}&=\frac{1}{2}\left(
\mid T_{+1}\rangle\langle T_{+1}\mid
-\mid T_{-1}\rangle\langle T_{-1}\mid
\right.\label{eq:Sop3}\\&\left.\hspace*{14mm}
\mp\mid S\rangle\langle T_{0}\mid
\mp\mid T_{0}\rangle\langle S\mid
\right).\nonumber
\end{align}
The spin exchange, and the potential-scattering amplitudes are related to the tunneling amplitudes appearing in~\eqref{eq:HT} as:
\begin{align}
 J_{i,rr'} &=t^{*}_{ir}t_{ir'}\left(\frac{1}{E_2-E_3}+\frac{1}{E_2-E_1}\right),\\
 W_{i,rr'} &=t^{*}_{ir}t_{ir'}\left(\frac{1}{E_2-E_3}-\frac{1}{E_2-E_1}\right),
\end{align}
where the energy denominators keep track of the energy cost for the virtual charge fluctuations to respectively the $N=3$ and the $N=1$ states. Assuming that we are deep inside the Coulomb diamond, we have omitted the slight differences in these denominators for the different spin-states. The potential scattering term vanishes close to the center of the diamond, where $E_3-E_2=E_1-E_2$. Since this gate voltage dependence has little influence on the line shape of the singlet-triplet crossing studied below, we shall assume that $|J_{i,rr'}| = |W_{i,rr'}|$. This has been done in the results shown in Figs.~\ref{fig:asymcuts} and~\ref{fig:symcutsV}.

The cotunneling current can now be found as~\cite{Bruus04book}
\begin{align}\label{eq:curr}
I=e\sum_{\eta'\eta}(\Gamma^{RL}_{\eta'\eta}
-\Gamma^{LR}_{\eta'\eta})P(\eta).
\end{align}
The non-equilibrium occupation numbers for the four singlet-triplet basis states, $P(\eta)$, are found by solving the rate equations
\begin{align}
0=\sum_{\eta'}\left(
\Gamma_{\eta'\eta}P(\eta)-\Gamma_{\eta\eta'}P(\eta')\right),
\end{align}
with the constraint that
\begin{align}
1=\sum_{\eta}P(\eta),
\end{align}
and with rates
\begin{align}
\Gamma_{\eta'\eta}=\sum_{r'r}
\Gamma^{r'r}_{\eta'\eta},
\end{align}
determined from Fermi's golden rule as
\begin{align}
\Gamma^{r'r}_{\eta'\eta}=
\gamma^{r'r}_{\eta'\eta}
\Delta\varepsilon^{r'r}_{\eta'\eta}
n_{B}(\Delta\varepsilon^{r'r}_{\eta'\eta}).
\end{align}
Here $n_{B}$ denotes the Bose function, and we have introduced the energy differences
\begin{align}
\Delta\varepsilon^{r'r}_{\eta'\eta}=
\varepsilon_{\eta'}-\varepsilon_{\eta}+\mu_{r}-\mu_{r'},
\end{align}
and the tunneling probabilities
\begin{align}
\gamma^{r'r}_{\eta'\eta}=&
\pi\rho^{2}\sum_{i=1,2}\left(\sum_{j=x,y,z}
|J_{i,r'r}|^{2}|T_{i,\eta'\eta}^{j}|^{2}
+|W_{i,r'r}|^{2}\delta_{\eta'\eta}\right),
\end{align}
where the matrix elements $T_{i,\eta'\eta}^{j}=\langle\eta'|S^{j}_{i}|\eta\rangle$ can be read off from Eqs. (\ref{eq:Sop1})--(\ref{eq:Sop3}).

If the device is coupled much stronger to one electrode than the other, the nonequilibrium occupations of excited states are strongly suppressed, and the finite-bias conductance becomes a simple succession of steps at bias voltages corresponding to energies of the excited states. Already from the elastic cotunneling amplitude, $\varGamma^{RL}_{\eta\eta}$, one can see that the conductance is larger for a triplet than for a singlet ground state
\begin{align}
\varGamma^{RL}_{SS}=&4\pi\rho^{2} V n_{B}(V)|W_{RL}|^{2},\\
\varGamma^{RL}_{T_{-1}T_{-1}}=&4\pi\rho^{2} V n_{B}(V)\left(|J_{RL}|^{2}+|W_{RL}|^{2}\right),
\end{align}
which amounts to a factor of two in difference between conductance on the singlet, or the triplet side if $|J_{i,rr'}|=|W_{i,rr'}|$. If the singlet-triplet boundary is completely vertical in the $(V_g,V)$-plane ($J_{12}$ depends only on $V_g$), this difference in conductance only shows up as a change in contrast in the conductance map as indicated in Fig.~\ref{fig:LevCross}. If, on the other hand, the singlet-triplet crossing has a slight bias-dependence, this already explains why Fig.~\ref{fig:diamonds} shows a peak for positive, and a dip for negative bias voltage. This is seen clearly in Fig.~\ref{fig:asymcuts}(a), which shows $dI/dV$ as a function of $V$ for two different values of $V_g$, calculated from formula \eqref{eq:curr}. The singlet-triplet crossing gives rise to a peak at positive bias (blue line) and a dip at negative bias (red line). Figure~\ref{fig:asymcuts}(b), shows the corresponding nonequilibrium occupation numbers, and since we have chosen a much stronger coupling to the right electrode, nonequilibrium effects are weak and one observes only slight changes in a small region about the crossing.
\begin{figure}
\includegraphics[width=0.9\columnwidth]{asymcuts2.pdf}
\caption{\label{fig:asymcuts}(a) Differential conductance $dI/dV$ vs. $V/T$ at $V_{g}/T=-300$ (blue upper line) and at $V_g/T=200$ (red lower line). The curves are offset for clarity. The tunneling amplitudes are very asymmetric, $t_{1R} = 30, t_{1L} = 0.7, t_{2R} = 0.8, t_{2L} = 0.5$ (arbitrary units). We have furthermore chosen
$B/T = 50$, $\alpha = 0.5$, and $\beta = 0.1$.
(b) Occupation probabilities of the different spin states at fixed bias $V/T=250$ with parameters as in (a). The state $|+\rangle$ refers to the lowest lying state, i.e. $|S\rangle$ on the l.h.s. and $|T_{-1}\rangle$ on the r.h.s.. $|-\rangle$ refers to the next lowest, i.e. $|T_{-1}\rangle$ on the l.h.s. and $|S\rangle$ on the r.h.s.}
\end{figure}

\begin{figure}
\includegraphics[width=0.9\columnwidth]{symcuts3.pdf}
\caption{\label{fig:symcutsV}(a) Differential conductance $dI/dV$ vs. $V/T$ at $V_{g}/T=-300$ (blue upper line) and at $V_g/T = 200$ (red lower line). The curves are offset for clarity. The tunnel amplitudes are now chosen more symmetric than in Fig.~\ref{fig:asymcuts}, with $t_{1R} = 0.3$,
and a finite spin-orbit coupling, $D/T=2\sqrt{2}$ has been included. All other parameters are the same as in Fig.~\ref{fig:asymcuts}.
(b) Corresponding occupation probabilities of the different spin states at fixed bias $V/T=250$. Notice the occupation inversion between the two states $|-\rangle$ and $|+\rangle$ at their would-be degeneracy point at $V_{g}/T=-150$.}
\end{figure}

\subsection{Including singlet-triplet mixing}\label{sec:STmix}

The cotunneling spectroscopy of the Mn complex in Ref.~\onlinecite{Osorio10} reveals a transition between the lowest lying and the highest lying triplet states (cf. r.h.s. of lower panel of Fig.~\ref{fig:FigExpCuts}). This implies a violation of the normal spin-selection rules for cotunneling, from which we infer that the relevant spin-states are not really true eigenstates of $S^{z}$, but instead a mixture of these inflicted by spin-orbit coupling in the molecule. We now add this information to our model by including a simple Dzyaloshinskii-Moriya term~\cite{Dzyaloshinsky58a,Moriya60a,Herzog10}, $H_{DM}=({\bf S}_{1}\times{\bf S}_{2})\cdot{\bf D}$, which reads
\begin{align}
H_{DM}=\frac{D}{2\sqrt{2}}\left[(|T_{+1}\rangle +|T_{-1}\rangle)\langle S|+|S\rangle(\langle T_{+1}| + \langle T_{-1}|)\right],
\end{align}
for ${\bf D}\perp{\bf B}$. This has the important effect of mixing the competing singlet and triplet states and will in general give rise to an avoided crossing as illustrated in Fig.~\ref{fig:LevCross}, with $\lambda\sim D$. If, however, the coupling strength $D$ is smaller than the bias width of the cotunneling peaks, this avoided crossing need not be observed. Interestingly, even such a small spin-orbit mixing may still have an observable effect on the line shape of the rather steep finite-bias singlet-triplet boundary. For simplicity we only include the part of $H_{DM}$ which mixes the singlet and lowest triplet ($T_{-1}$) states, which is the only relevant term close to the singlet-triplet line.

Including $H_{DM}$, we may determine the new eigenstates at a given bias voltage, $V$, which still enters via $J_{12}(V_g,V)$. Using the same procedure as outlined above, the cotunneling current is readily recalculated. The result obtained for $D/T=2\sqrt{2}$ is shown in Fig.~\ref{fig:symcutsV}(a), revealing a pronounced peak-dip, and dip-peak structure at respectively positive and negative bias, in qualitative correspondence with the experimental observation. Whereas the peak (dip) in Fig.~\ref{fig:asymcuts}, with no spin-orbit coupling included, could only be obtained with a substantial difference in the couplings to respectively source and drain electrodes, the peak-dip (dip-peak) structure in Fig.~\ref{fig:symcutsV} only arises when the couplings to source and drain are not too different. The observed peak-dip structure is thus a nonequilibrium effect arising from the bias dependent occupation numbers of the various spin-states, shown in Fig.~\ref{fig:symcutsV}(b) to exhibit a population inversion of the two competing ground states near the crossing, around $V_g/T\sim -150$ in this plot. Apart from the requirement of nearly equal tunnel couplings, this population inversion and the concomitant peak-dip structure does not rely on any fine-tuning of parameters. Fig.~\ref{fig:DensPlot} shows this result in a full density plot of the nonlinear conductance with a number of gate voltage traces taken at different bias voltages, all plotted in the same style as the experimental data in Fig.~\ref{fig:FigExpCuts}. The characteristic Fano-like lineshape is seen both when the conductance is plotted as a function of $V$ and of $V_g$.

Reading off the slope of the nearly vertical singlet-triplet boundary and the slope of the inelastic singlet-triplet cotunneling threshold in zero magnetic field, one can now readily determine the constants $\alpha$ and $\beta$ from formula \eqref{eq:Jfct}. Doing this, we find that $\alpha\approx 0.002$ and $\beta\approx 0.020$. In Fig.~\ref{fig:FigExpCuts}, the width in gate voltage of the peak-dip structure is roughly $0.05$ V, and using $\alpha$ this now translates into a corresponding change in $J_{12}$ of the order of 0.1 meV, which must be an indirect measure of the effective singlet-triplet mixing, $D$, in this Mn-complex.

\section{Conclusions}

Inspired by the experimental results for the conductance map of a Mn complex reported in Ref.~\onlinecite{Osorio10}, we have studied the signatures in the nonlinear conductance near a finite-bias singlet-triplet transition. For asymmetrically coupled devices the nonequilibrium effects are weak. A mere difference in elastic cotunneling amplitudes for respectively singlet and triplet states can give rise to current plateaus of different heights and hence a conductance peak along the transition line in the $(V,V_{g})$-plane. For more symmetrically coupled devices, this effect will be washed out at sufficiently large bias voltages. In this case, however, a new possibility arises, if the two competing states are mixed by e.g. spin-orbit coupling. In this case, the nonequilibrium occupations exhibit a momentary population inversion near the transition line, which in turn gives rise to a peak-dip anomaly in the conductance. This provides an explanation for the observations in Ref.~\onlinecite{Osorio10}, and it provides an indirect method of gauging the size of the singlet-triplet mixing.

This type of cotunneling signature should be rather general, and the considerations about the singlet-triplet transition studied here could readily be exported to entirely different Coulomb blockaded systems with level crossings induced by gate-voltage. If the competing states are sensitive to gate voltage, it is only natural to expect also a slight sensitivity to bias voltage and thus a non-vertical transition line as depicted in Fig.~\ref{fig:LevCross}. As was the case for the Mn complex discussed here, the steepness of this line can help bring out otherwise hidden spectroscopic details.

\section{acknowledgements}
The research leading to these results has received funding from the European Union Seventh Framework Programme
(FP7/2007-2013) under agreement no 270369 (“ELFOS”).

\bibliographystyle{apsrev}