Witnessing entanglement by proxy

As mentioned in the introduction, the mean energy of a system can be used to witness entanglement in the corresponding quantum state. Namely any state with mean energy less than

$$\begin{eqnarray}&&{E}_{\mathrm{min},\mathrm{sep}}=\underset{\rho \in \mathrm{sep}}{\mathrm{min}}\mathrm{Tr}\rho H\end{eqnarray} \tag{ 1 }$$

is entangled. By convexity the minimum is attained in a pure state. $\mathrm{sep}$ can be chosen to be the set of fully separable states or the set of k-separable states. For k = 2, genuinely multipartite entanglement is detected. If the system is in thermal equilibrium, it is possible to derive an analogous criterion for the temperature.

The goal of this chapter is to find a condition that is able to detect entanglement at higher energies than ${E}_{\mathrm{min},\mathrm{sep}}$. The idea is to add additional constraints to (1). For example, a lower bound on the (von Neumann) entropy of ρ:

$$\begin{eqnarray}&&{E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}=\underset{\rho \in \mathrm{sep},\ {\rm{S}}(\rho )\geqslant {{\rm{S}}}^{\prime }}{\mathrm{min}}\mathrm{Tr}\rho H.\end{eqnarray} \tag{ 2 }$$

If $E(\rho )\lt {E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}$ and $S(\rho )\geqslant {S}^{\prime }$, ρ is entangled. ${S}^{\prime }$ can be varied between 0 and $\mathrm{ln}d$. Clearly, it holds ${E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}\geqslant {E}_{\mathrm{min},\mathrm{sep}}$. But the relevant question is whether there are entropies for which the strict inequality holds. While (1) is minimised by a pure state, the minimum does not have to be unique. Hence a mixture of minimisers could also have mean energy ${E}_{\mathrm{min},\mathrm{sep}}$, but at non-zero entropy. The exact behaviour depends, of course, on the Hamiltonian. However, it is possible to show that if there exists an ${S}^{\prime }$, for which the strict inequality holds, it will hold for any larger entropy. This follows from

Lemma 1.  ${E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}$ as function of ${S}^{\prime }$ is convex.

Proof. Let $0\leqslant {S}_{\mathrm{1,2}}\leqslant \mathrm{ln}d$ and $0\leqslant p\leqslant 1$. Let ${\rho }_{\mathrm{1,2}}$ be separable states with entropies ${S}_{1,2}^{\prime }$ respectively, that minimise (2). Then

$$\begin{eqnarray}&&{E}_{\mathrm{min},\mathrm{sep},{{\rm{S}}}^{\prime }}({{pS}}_{1}^{\prime }+(1-p){S}_{2}^{\prime })\leqslant \mathrm{Tr}H(p{\rho }_{1}+(1-p){\rho }_{2})={{pE}}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}({S}_{1}^{\prime })+(1-p){E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}({S}_{2}),\end{eqnarray} \tag{ 3 }$$

where the inequality is due to the fact that $p{\rho }_{1}+(1-p){\rho }_{2}$ is a feasible point of (2), which follows from the concavity of the entropy.

Since ${E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}$ is convex, it is strictly monotonically increasing as soon as it exceeds ${E}_{\mathrm{min},\mathrm{sep}}$. Let us call the smallest entropy where the constraint hits ${S}_{\mathrm{min}}$. See also figure 1.

Zoom In Zoom Out Reset image size

Since for ${S}^{\prime }\geqslant {S}_{\mathrm{min}}$, ${E}_{\mathrm{min},\mathrm{sep},{S}^{\prime }}$ is strictly monotonically increasing, we can obtain the same curve, i.e. detect the same entanglement, by computing

$$\begin{eqnarray}&&{S}_{\mathrm{max},\mathrm{sep},E}=\underset{\rho \in \mathrm{sep},\ \mathrm{Tr}H\rho =E}{\mathrm{max}}S(\rho )\end{eqnarray} \tag{ 4 }$$

and varying ${E}_{\mathrm{min},\mathrm{sep}}\leqslant E\leqslant \frac{\mathrm{Tr}H}{d}$. Equivalently we could also demand $\mathrm{Tr}H\rho \leqslant E$. Note that if we remove the separability constraint (4) will become the Gibbs state entropy

$$\begin{eqnarray}&&{S}_{\mathrm{Gibbs}}=\underset{\rho \mathrm{state},\ \mathrm{Tr}H\rho =E}{\mathrm{max}}S(\rho ),\end{eqnarray} \tag{ 5 }$$

attained by the Gibbs state. Since the optimisation is of the same form as the one yielding the regular Gibbs ensemble we will call the resulting state a separable Gibbs ensemble. In order to compare its entropy with the Gibbs state entropy, let us define the entropy gap ${\rm{\Delta }}S={S}_{\mathrm{Gibbs}}-{S}_{\mathrm{max},\mathrm{sep},E}$.

Zoom In Zoom Out Reset image size

Theorem 2. For mean energy E and corresponding temperature T, any state ρ with ${S}_{\mathrm{max},\mathrm{sep},E}\lt S(\rho )\leqslant {S}_{\mathrm{Gibbs}}$ is entangled. In particular the Gibbs state is entangled if ${\rm{\Delta }}S\gt 0$.

Note that for any non-zero gap between ${S}_{\mathrm{max},\mathrm{sep},E}$ and ${S}_{\mathrm{Gibbs}}$ there always exist states different from the Gibbs state which we detect, for example there is always a state with entropy $\lambda {S}_{\mathrm{max},\mathrm{sep},E}+(1-\lambda ){S}_{\mathrm{Gibbs}}$, where $0\leqslant \lambda \leqslant 1$.

The optimisation in (4) is difficult to deal with because it contains the separability constraint. The main idea is to relax this constraint using sets of states which remain positive semidefinite after the application of positive (yet not completely positive) maps Λ. A prominent example would be the partial transposition. While these are of course supersets of the separable states, it is clear that for every entangled state in principle there exists a map Λ and thus a semidefinite relaxation that will still yield optimal results for the constrained optimisation.

Let us use the fact that every set of Λ-positive states forms a convex set that can be approximated by a suitable set of entanglement witnesses W_i with suitable weights ${\nu }_{i}\geqslant 0$ [24]. As every optimal entanglement witness for the set of Λ-positive states for partition A can be written as ${{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[| \psi \rangle \langle \psi | ]$ one can find suitable entanglement witnesses for the system in question. While this of course comes at the expense of finding suitable entanglement witnesses for specific systems, it also opens the possibility to constrain the entropy beyond just states which are separable under fixed bi-partitions. I.e. it enables us to find also genuine multipartite entanglement or any other non-partially separable set, by choosing corresponding witnesses. In order to remain fully general we include both an arbitrary set of entanglement witnesses and positive maps in the following considerations.

The relaxed problem is thus given by

$$\begin{eqnarray}{S}_{\mathrm{max},\mathrm{sep},E}\leqslant {S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E} & = & \mathrm{max}S(\rho )\\ & & {\rm{s}}.{\rm{t}}.\rho \mathrm{state},\mathrm{Tr}(\rho H)=E,\\ & & {{\rm{\Lambda }}}_{A}\otimes {{\mathbb{1}}}_{\bar{A}}[\rho ]\geqslant 0,\mathrm{Tr}(\rho {W}_{i})\geqslant 0.\end{eqnarray} \tag{ 6 }$$

${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}$ will be an upper bound on ${S}_{\mathrm{max},\mathrm{sep},E}$, it is tightness depending on the choice of witnesses or maps. Since the relaxation just provides an upper bound for the separable entropy it trivially follows that according to theorem 2:

Corollary 3. For mean energy E and corresponding temperature T, any state ρ with ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\lt S(\rho )\leqslant {S}_{\mathrm{Gibbs}}$ is entangled. In particular the Gibbs state is entangled if ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\lt {S}_{\mathrm{Gibbs}}$. For an illustration see figure 2.

While the relaxation weakens the detection criterion, we are now dealing with a convex optimisation problem with only linear and semidefinite constraints (since maximising a concave function is of course equivalent to minimising a convex one). Such programs can be solved efficiently numerically thanks to so-called interior point methods [25]. In addition, they have a duality theory which can be used to give certified upper bounds on ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}$, as we will discuss in the next section. One might ask why instead of applying corollary 3, one cannot simply apply the witnesses or positive maps involved directly. While this is possible theoretically, let us note again that witnesses or positive maps are in general not easily accessible in experiments, while macroscopic variables such as the mean energy and entropy are. Since corollary 3 makes use of witnesses without the need to measuring them directly, but instead requires the measurement of mean energy and entropy, we call those two quantities proxy witnesses.

Let us conclude this section by noting that apart from the von Neumann entropy, in principle any concave function can serve as a proxy witness. An example is the so called linear entropy, also known as impurity of a quantum state, which is defined by ${S}_{L}(\rho )=1-\mathrm{Tr}{\rho }^{2}$ [26]. The linear entropy is upper bounded by the von Neumann entropy and can be seen as a measure of mixedness of a quantum state, as well. The optimisation of the linear entropy is easier to deal with than of the von Neumann entropy, as it is only quadratic in ρ.

For convex optimisation problems, such as (6), it is possible to derive a dual optimisation problem. To do so, one defines the Lagrangian, a function incorporating both objective, i.e. the function to be maximised, and the constraints. The constraints are added by means of Lagrange multipliers. The Lagrange multipliers are referred to as dual variables, whereas the variables of the original problem are referred to as primal. Maximisation of the Lagrangian over all admissible primal points yields the dual objective. The dual problem is then given by the minimisation of the dual objective with respect to the dual variables. It can be shown that any dual feasible point provides an upper bound on the original problem, which is referred to as weak duality. For details, please refer to [25]. For (6),the Lagrangian is given by

$$\begin{eqnarray}L(\rho ,\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}) & = & S(\rho )+\lambda (\mathrm{Tr}\rho -1)+\mu (\mathrm{Tr}\rho H-E)+\displaystyle \sum _{i}{\nu }_{i}\mathrm{Tr}\rho {W}_{i}\\ & & +\mathrm{Tr}\rho {X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}\mathrm{Tr}{{\rm{\Lambda }}}_{A}\otimes {{\mathbb{1}}}_{\bar{A}}[\rho ]{X}_{A,{\rm{\Lambda }}},\end{eqnarray} \tag{ 7 }$$

where $\lambda ,\mu$ and the ${\nu }_{i}$ are real Lagrange dual variables, corresponding to the trace, mean energy and witness constraints, respectively. As the witness constraint is given by an inequality, we can restrict to ${\nu }_{i}\geqslant 0$. X₀ and the ${X}_{A,{\rm{\Lambda }}}$, the Lagrange dual variables corresponding to the positivity and positive maps constraints respectively, are positive semidefinite matrices.

Here it becomes clear that witnesses are easier to deal with numerically than positive maps, as they only require scalar variables. The dual objective function reads as follows

$$\begin{eqnarray}&&{\ell }(\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}):= \underset{\rho }{\mathrm{max}}L(\rho ,\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}),\end{eqnarray} \tag{ 8 }$$

where the maximisation is over all admissible ρ, i.e. those for which the right-hand side is defined, in particular not necessarily obeying the primal constraints in (6). Note that the dual problem is also a convex optimisation problem: the constraints are indeed linear and semidefinite, while ${\ell }$ is convex, owed to the linearity of L in the dual variables.

Clearly, for every primal feasible ρ and dual feasible $\lambda ,\mu ,{\nu }_{i},{X}_{0}$ and ${X}_{A,{\rm{\Lambda }}}$ it holds

$$\begin{eqnarray}&&S(\rho )\leqslant L(\rho ,\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}})\leqslant {\ell }(\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}).\end{eqnarray} \tag{ 9 }$$

In particular we have what is referred to as weak duality:

$$\begin{eqnarray}&&{S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\leqslant \mathrm{min}{\ell }(\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}),\end{eqnarray} \tag{ 10 }$$

where the minimisation is over the dual feasible region, i.e. $\lambda ,\mu \in {\mathbb{R}},{\nu }_{i}\geqslant 0,{X}_{0}\geqslant 0,{X}_{A,{\rm{\Lambda }}}\geqslant 0$. The right-hand side is referred to as the Lagrange dual problem. Note that any dual feasible bound yields an upper bound on ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}$. This allows us to obtain analytical upper bounds from numerics. One simply has to numerically optimise the dual problem and check if the so obtained optimisers are dual feasible. If they are they can be inserted into ${\ell }$ yielding the bound.

Let us now evaluate the dual objective function (8) for the von Neumann entropy, which is well-defined on any positive semidefinite matrix. Since the entropy is concave, so is the Lagrangian. Hence it is sufficient to find a critical point, that is positive semidefinite, of L as function of ρ. Using ${{\rm{\nabla }}}_{\rho }S(\rho )=-{\mathbb{1}}-\mathrm{ln}\rho$ and ${{\rm{\nabla }}}_{\rho }\mathrm{Tr}\rho M=M$ with respect to the trace inner product on Hermitian matrices, as well as the fact that $\mathrm{Tr}({{\rm{\Lambda }}}_{A}\otimes {{\mathbb{1}}}_{\bar{A}}[\rho ]{X}_{A,{\rm{\Lambda }}})=\mathrm{Tr}(\rho {{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])$ one can obtain

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\rho }L=-{\mathbb{1}}-\mathrm{ln}\rho +\lambda {\mathbb{1}}+\mu H+\displaystyle \sum _{i}{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}],\end{eqnarray} \tag{ 11 }$$

which vanishes for

$$\begin{eqnarray}&&{\rho }^{\mathrm{crit}}=\mathrm{exp}((\lambda -1){\mathbb{1}}+\mu H+\displaystyle \sum _{i}{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}]).\end{eqnarray} \tag{ 12 }$$

Hence the dual objective is given by

$$\begin{eqnarray}&&{\ell }(\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}})=L({\rho }^{\mathrm{crit}},\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}}),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&=\mathrm{Trexp}((\lambda -1){\mathbb{1}}+\mu H+\displaystyle \sum _{i}{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])-\lambda -\mu E,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&={e}^{(\lambda -1)}\mathrm{Trexp}(\mu H+\displaystyle \sum _{i}{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])-\lambda -\mu E.\end{eqnarray} \tag{ 15 }$$

We will have to minimise ${\ell }$ with respect the dual variables. Since ${X}_{0}\geqslant 0$, the minimum will be attained at ${X}_{0}=0$. ${\ell }$ can be easily minimised for λ. This is due to the fact that the eigenvalues of a hermitian matrix cannot decrease if a positive semi-definite matrix is added, which follows from theorem 4.3.1 of [27]. By convexity the minimum is attained where the derivatives vanishes, yielding a function $\tilde{{\ell }}$ just of $\mu ,{\nu }_{i},{X}_{A,{\rm{\Lambda }}}$

$$\begin{eqnarray}&&\tilde{{\ell }}(\mu ,{\nu }_{i},{X}_{A,{\rm{\Lambda }}})=\mathrm{ln}\;\mathrm{Tr}\;\mathrm{exp}(\mu H+\displaystyle \sum _{i}{\nu }_{i}{W}_{i}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])-\mu E.\end{eqnarray} \tag{ 16 }$$

By weak duality, it holds

$$\begin{eqnarray}&&{S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\leqslant \underset{\mu \in {\mathbb{R}},{\nu }_{i}\geqslant 0,{X}_{A,{\rm{\Lambda }}}\geqslant 0}{\mathrm{min}}\tilde{{\ell }}(\mu ,{\nu }_{i},{X}_{A,{\rm{\Lambda }}}).\end{eqnarray} \tag{ 17 }$$

As ${\ell }$ is convex, $\tilde{{\ell }}$ is convex, as well [25]. If only witnesses are used, $\tilde{{\ell }}$ only has scalar variables, which results in a greatly enhanced numerical performance compared to the primal problem.

It is an interesting observation that $\tilde{{\ell }}$ is reminiscent of a grand canonical ensemble and

$$\begin{eqnarray}&&\tilde{{\ell }}=\mathrm{ln}{Z}^{\prime }+\beta E=\mathrm{ln}\;\mathrm{Tr}\;\mathrm{exp}(-\beta {H}^{\prime })+\beta E,\end{eqnarray} \tag{ 18 }$$

where ${H}^{\prime }=H-\displaystyle \frac{1}{\beta }({\displaystyle \sum }_{i}{\nu }_{i}{W}_{i}+{\displaystyle \sum }_{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])$. Instead of particle numbers the constraints stem from the specific witnesses or maps used. We will sometimes refer to equation (18) as a witness canonical ensemble entropy (WCEE).

In other words, when choosing $\mu =-\beta$, the operator ${\displaystyle \sum }_{i}{\nu }_{i}{W}_{i}+{\displaystyle \sum }_{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}]$ can be seen as some sort of grand canonical ensemble with additional 'chemical' or rather witness potentials. In particular the entropy gap can be lower bounded as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}S\geqslant \mathrm{ln}\displaystyle \frac{\mathrm{Tr}\;\mathrm{exp}(-\beta H)}{\mathrm{Trexp}(-\beta H+{\displaystyle \sum }_{i}{\nu }_{i}{W}_{i}+{\displaystyle \sum }_{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])}\end{eqnarray} \tag{ 19 }$$

for all ${\nu }_{i}\geqslant 0$ and ${X}_{A,{\rm{\Lambda }}}\geqslant 0$.

It is also possible to compute ${\ell }$ for the linear entropy. After rewriting the Lagrangian as

$$\begin{eqnarray}&&{L}_{\mathrm{LIN}}=1-\lambda -\mu E+\mathrm{Tr}\rho (\lambda {\mathbb{1}}+\mu H+{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])-\mathrm{Tr}{\rho }^{2},\end{eqnarray} \tag{ 20 }$$

it is easy to show that the optimiser is given by

$$\begin{eqnarray}&&{\rho }^{\mathrm{crit}}=\displaystyle \frac{1}{2}(\lambda {\mathbb{1}}+\mu H+{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}]).\end{eqnarray} \tag{ 21 }$$

Hence

$$\begin{eqnarray}&&{{\ell }}_{\mathrm{LIN}}(\lambda ,\mu ,{\nu }_{i},{X}_{0},{X}_{A,{\rm{\Lambda }}})=\displaystyle \frac{1}{4}\mathrm{Tr}{(\lambda {\mathbb{1}}+\mu H+{\nu }_{i}{W}_{i}+{X}_{0}+\displaystyle \sum _{A,{\rm{\Lambda }}}{{\rm{\Lambda }}}_{A}^{*}\otimes {{\mathbb{1}}}_{\bar{A}}[{X}_{A,{\rm{\Lambda }}}])}^{2}+1-\lambda -\mu E.\end{eqnarray} \tag{ 22 }$$

Minimisation of ${{\ell }}_{\mathrm{LIN}}$ yields an upper bound on the maximal linear entropy achievable by separable states. Again, this minimisation is easier to deal with as in the von Neumann case, because it is only a quadratic function.

Let us now briefly discuss how the primal and dual optimisation problems introduced in the preceding sections can be implemented numerically. As mentioned before, (6) and (17) have concave and convex objectives, respectively, as well as linear and semidefinite constraints. Note that for our purposes it is sufficient to only compute the dual problems. It can however be instructive to also compute the primal problem in order to obtain the optimiser and check if strong duality holds. While interior-point methods can in principle solve such problems [25], readily available solvers such as Sedumi [28] or SDPT3 [29] can only handle linear and quadratic objectives. This is sufficient to solve (6) and minimise (22) for the Linear entropy.

For the von Neumann entropy, however, there is a way to obtain an approximate solution: it is a well known fact that the von Neumann entropy of a state is equal to the Shannon entropy of its eigenvalues. As for $\tilde{{\ell }}$, note that for a hermitian n × n matrix M with eigenvalues ${\lambda }_{i}$ it holds $\mathrm{Trexp}(M)={\displaystyle \sum }_{i}\mathrm{exp}({\lambda }_{i})$. Hence $\tilde{{\ell }}$ is a function of the eigenvalues of the exponent. Both S and $\tilde{{\ell }}$ are invariant under permutation of the eigenvalues. This allows us to reformulate (6) as

$$\begin{eqnarray}{S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E} & = & \underset{{\bf{v}},\rho }{\mathrm{max}}H({\bf{v}}),\\ & & {\rm{s}}.{\rm{t}}.{\bf{v}}=\mathrm{eig}(\rho ),\ \rho \mathrm{state},\ \mathrm{Tr}(\rho H)=E,\\ & & {{\rm{\Lambda }}}_{A}\otimes {{\mathbb{1}}}_{\bar{A}}[\rho ]\geqslant 0,\mathrm{Tr}(\rho {W}_{i})\geqslant 0,\end{eqnarray} \tag{ 23 }$$

where $\mathrm{eig}(\rho )$ denotes the vector of eigenvalues of ρ. Note that this is not a semidefinite constraint. There is, however, a trick to include the eigenvalues into a semidefinite programme[30]: while the eigenvalues are generally not SDP-representable, the sum of the k largest eigenvalues of a matrix is [31]. We can now replace the eigenvalue constraint in (23) by the constraint that ${\bf{v}}$ has to majorise $\mathrm{eig}(\rho )$, i.e.

$$\begin{eqnarray*} & & \mathrm{Tr}(\rho )=\displaystyle \sum _{i}{v}_{i}\\ & & {s}_{1}(\rho )\leqslant {v}_{1}\\ & & {s}_{2}(\rho )\leqslant {v}_{1}+{v}_{2}\\ & & ...\\ & & {s}_{n-1}(\rho )\leqslant \displaystyle \sum _{i=1}^{n-1}{v}_{i},\end{eqnarray*}$$

where s_k(M) denotes the sum of the k largest eigenvalues of a matrix M. Then the optimising ${\bf{v}}$ will be equal to the eigenvalues of the optimising ρ. To see this, recall that if ${\bf{v}}\succ {\bf{w}}$, it holds $H({\bf{v}})\leqslant H({\bf{w}})$ [32]. Let us now assume that the (unique) optimising ${\bf{v}}$ of the Shannon entropy majorises but is not equal to the eigenvalues of the optimal ρ. Then $S(\rho )$ would be greater or equal to $H({\bf{v}})$, which is a contradiction. Since ${\bf{v}}\succ {\bf{w}}$ implies $f({\bf{v}})\geqslant f({\bf{w}})$ for any convex function [32], the same argument can be applied to (17).

So far, we have replaced the von Neumann by the Shannon entropy and, in the dual problem, the matrix exponential by scalar exponential functions and transformed the arising eigenvalue constraints into semidefinite constraints. In order to apply Sedumi or SDPT3, all that is left to do is to approximate the objectives by piecewise linear functions. While this will only give us approximate solutions, let us note that the optimisers found in this way, can be easily checked to be feasible and inserted in the original objective.

When using PPT or other positive maps the dimension of the matrix variables increases exponentially with the number of qubits. As we will show in the next section, this makes it difficult to go beyond five qubits. The computation of (16) greatly simplifies when only witnesses are used. In this case there will only be scalar variables $\mu$ and ${\nu }_{i}$ and no semidefinite constraints. This allows for application of a standard nonlinear solver, such as FMINCON [33]. Using FMINCON we were able to obtain results for up to 13 qubits, as we will present in the next section.

If we want to show if the Gibbs state is entangled, we will also need to compute the Gibbs state entropy $S({\rho }_{\mathrm{Gibbs}})=\beta E+\mathrm{ln}Z$. To do so numerically it is sufficient to compute the eigenvalues of the Hamiltonian, which allows for computation of the partition function $Z={\mathrm{Tre}}^{-\beta H}$ and the mean energy $E=\mathrm{Tr}{\rho }_{\mathrm{Gibbs}}H=\displaystyle \frac{1}{Z}{\displaystyle \sum }_{i}{{\rm{e}}}^{-\beta {E}_{i}}{E}_{i}$.

In order to test our method, we have implemented it for Heisenberg model, which was introduced in order to simplify the analysis of systems of spins, such as ferro- or antiferromagnets. It only takes into account the nearest neighbour exchange interaction between the spins as well an external magnetic field. In the one-dimensional case, i.e. a chain of spins, the Heisenberg model is described by the following Hamiltonian

$$\begin{eqnarray}&&H=-\displaystyle \sum _{i=1}^{N}({J}_{x}{\sigma }_{i}^{x}{\sigma }_{i+1}^{x}+{J}_{y}{\sigma }_{i}^{y}{\sigma }_{i+1}^{y}+{J}_{z}{\sigma }_{i}^{z}{\sigma }_{i+1}^{z})+B\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{z},\end{eqnarray} \tag{ 24 }$$

where ${J}_{x},{J}_{y},{J}_{z}$ are the coupling constants for the $x,y,z$-components of the spins, N the number of spins and B the external magnetic field in zdirection. ${\sigma }_{i}^{x},{\sigma }_{i}^{y},{\sigma }_{i}^{z}$ denote the Pauli operators for the ith spin. Let us assume periodic boundary conditions, i.e. a ring of spins. If ${J}_{x}={J}_{y}={J}_{z}=:J$, we are talking about an isotropic XXX Heisenberg model. $J\gt 0$ and $J\lt 0$ correspond to ferromagnetic and antiferromagnetic systems, respectively. If ${J}_{x}={J}_{y}\geqslant {J}_{z}$, the system is called an XXZ system and so on. If only one component of the spin is considered, i.e. only one ${J}_{i}\ne 0$, the Heisenberg model reduces to the Ising model.

The numerical results presented below have been obtained using either Sedumi or SDPT3 as well as Yalmip [34]. Using PPT constraints we have applied the method for up to five qubits, for witnesses for up to 13 qubits.

2.4.1. Antiferromagnetic Heisenberg model

For the one-dimensional antiferromagnetic Heisenberg model (XXX with $J=-1$), it is possible to detect entanglement at higher energies than [12–14], both using the von Neumann and the linear entropy as proxies. The computations have been performed using both the partial transpose with respect to the partition ${A}_{i\mathrm{even}}| {A}_{j\mathrm{odd}}$ and all possible partitions.

Let us start with the results for the even versus odd partition. In table 1 the energy ranges where entanglement is detected are shown for the von Neumann entropy and for the linear entropy, respectively. Here N denotes the number of qubits, E₀ the ground state energy and ${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}$ denotes the smallest mean energy allowing for PPT w.r.t. the even versus odd partition . Below that mean energy any state is guaranteed to be entangled. Note that ${E}_{\mathrm{min},\mathrm{sep}}\geqslant {E}_{\mathrm{min},\mathrm{PPT}-\mathrm{all}}\geqslant {E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}$. Between ${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}$ and ${E}_{\mathrm{max},\mathrm{gap}}$ all states falling into the entropy gap are entangled. I.e. ${E}_{\mathrm{max},\mathrm{gap}}$ is the largest energy for which the proxy method can work. This includes the Gibbs state in the von Neumann case. Note that the Gibbs state is not necessarily the maximiser for the linear entropy. Hence, instead of the Gibbs state we have computed the state with maximal linear entropy

$$\begin{eqnarray}&&{S}_{L,\mathrm{max}}=\underset{\rho \mathrm{state},\ \mathrm{Tr}H\rho =E}{\mathrm{max}}{S}_{L}(\rho )\end{eqnarray} \tag{ 25 }$$

and compared it to the constraint linear entropy. ${\mathcal{E}}$ is defined as the fraction of the energy range where entanglement is detected

$$\begin{eqnarray}&&{\mathcal{E}}=\displaystyle \frac{{E}_{\mathrm{max},\mathrm{gap}}-{E}_{0}}{{E}_{\mathrm{max}}-{E}_{0}},\end{eqnarray} \tag{ 26 }$$

where ${E}_{\mathrm{max}}$ denotes the largest energy eigenvalue.

Table 1.  Energy ranges in which the von Neumann entropy can be successfully used a proxy for detecting entanglement in the antiferromagnetic Heisenberg model ($J=-1$).

N	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
3	−1.000	−1.000	−0.610	0.195
4	−2.000	−1.000	−0.660	0.447
5	−1.494	−1.008	−0.695	0.320

Table 2.  Energy ranges in which the linear entropy can be successfully used a proxy for detecting entanglement in the antiferromagnetic Heisenberg model ($J=-1$).

N	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
3	−1.000	−1.000	−0.600	0.200
4	−2.000	−1.000	−0.370	0.543
5	−1.494	−1.008	−0.302	0.478

As can be seen in tables 1 and 2, entanglement can be detected in a large area of the energy spectrum. The maximum mean energy where entanglement can be detected is also substantially higher than $-N$, which is the maximum mean energy where the methods of [12–14] work. Let us also note that the ground state energies given in [12–14] are only correct in the limit of large N, as has been noted in [14]. The fact that the linear entropy detects more than the von Neumann entropy could be a result of the piecewise linear approximation of the von Neumann entropy in the maximisation, which results in weaker bounds on ${S}_{\mathrm{max},\mathrm{sep}}$.

The entropy gap can be seen in figures 3(a) and (b) for the von Neumann entropy and for the linear entropy, respectively. Here the Gibbs state entropy or ${S}_{{\rm{L}},\mathrm{max}}$, as well as ${S}_{\mathrm{max},\mathrm{PPT}-\mathrm{even}-\mathrm{dd}}$ are plotted versus the mean energy.

Zoom In Zoom Out Reset image size

Entanglement can also be detected in the presence of a magnetic field. Tables 3 and 4 show the detected energy ranges in a three qubit system for the von Neumann and linear entropies, respectively. The entropies are plotted in figure 4 for B = 3.

Zoom In Zoom Out Reset image size

Table 3.  Energy ranges for three qubits in the antiferromagnetic Heisenberg model ($J=-1$) in a magnetic field B = 3 using von Neumann entropy.

B	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
0	−1.000	−1.000	−0.610	0.195
1	−1.333	−1.333	−0.720	0.184
2	−1.667	−1.667	−1.033	0.136
3	−2.000	−2.000	−1.540	0.077

Table 4.  Energy ranges for three qubits in the antiferromagnetic Heisenberg model ($J=-1$) with a magnetic field B = 3 using linear entropy.

B	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
0	−1.000	−1.000	−0.600	0.200
1	−1.333	−1.333	−0.560	0.232
2	−1.667	−1.667	−0.717	0.204
3	−2.000	−2.000	−0.940	0.177

Going from PPT w.r.t. the even–uneven partition to all possible partitions greatly increases the computation time as the number of partitions grows exponentially with the number of qubit. Since the constrained entropy increases monotonically in E, it holds ${E}_{\mathrm{max},\mathrm{gap},\mathrm{PPT}-\mathrm{all}}\geqslant {E}_{\mathrm{max},\mathrm{gap},\mathrm{PPT}-\mathrm{even}-\mathrm{odd}}$, possibly increasing the energy range, where entanglement can be detected. In the example considered here, however, only a small increase can be obtained. See tables 5 and 6, as well as figure 5. For three qubits it holds ${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{all}}/N=-0.6\gt -1$, implying that ${E}_{\mathrm{min},\mathrm{sep}}/N\geqslant -0.6\gt -1$. This shows that the result of [12, 14] is suboptimal for an odd number of qubits, as mentioned in [14]. The reason is that they use a partition into two sub-lattices such that every neighbouring sites belong to different sub-lattices, which is not possible for an odd number of qubits.

Zoom In Zoom Out Reset image size

Table 5.  Results for the von Neumann entropy in the antiferromagnetic Heisenberg model ($J=-1$) with a magnetic field B = 3 using all possible partitions.

N	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{all}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
3	−1.000	−0.600	−0.600	0.200
4	−2.000	−1.000	−0.660	0.447
5	−1.494	−0.809	−0.809	0.275

Table 6.  Results for the linear entropy in the antiferromagnetic Heisenberg model ($J=-1$) with a magnetic field B = 3 using all possible partitions.

N	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{PPT}-\mathrm{all}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${\mathcal{E}}$
3	−1.000	−0.600	−0.594	0.203
4	−2.000	−1.000	−0.360	0.547

2.4.2. Dicke states

Dicke states were first considered in the theory of coherent spontaneous light emission [35]. They are defined by

$$\begin{eqnarray}&&| {D}_{m}^{n}\rangle ={\displaystyle (\displaystyle \genfrac{}{}{0em}{}{n}{m})}^{-\frac{1}{2}}\displaystyle \sum _{\alpha }| {d}_{\alpha }\rangle ,\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&| {d}_{\alpha }\rangle =\underset{i\notin \alpha }{\displaystyle \otimes }{| 0\rangle }_{i}\underset{i\in \alpha }{\displaystyle \otimes }{| 1\rangle }_{i},\end{eqnarray} \tag{ 28 }$$

m is the number of excitations, $\alpha$ denote sets of indices of excited subsystems and the sum is taken over all inequivalent sets of m indices.

More recently Dicke states turned out to be a useful resource for quantum information processing task as they are LOCC transferable to GHZ or W states [36]. Several experiments have successfully created Dicke states, e.g. [36–38].

In [39], it has been shown using perturbation theory that the approximate ground states of anisotropic ferromagnetic XXZ Heisenberg Hamiltonians (i.e. ${J}_{x}={J}_{y}\geqslant {J}_{z}\gt 0$) are Dicke states. For B = 0 and ${J}_{x}={J}_{z}\gt 0$ the ground states are given by $\{| {D}_{m}^{n}\rangle \}{}_{m=0}^{n}$ with $n+1$-fold degeneracy, which makes entanglement at low energies unlikely. This is because the system will soon go into a mixture of the Dicke states, which is separable. For ${J}_{x}\gt {J}_{z}$ and $B\ne 0$, however, the degeneracy vanishes. In first order perturbation theory $| {D}_{m}^{n}\rangle$ is the ground state for

$$\begin{eqnarray}&&-\displaystyle \frac{n-2m+1}{n-1}{\rm{\Delta }}J\lt B\lt -\displaystyle \frac{n-2m-1}{n-1}{\rm{\Delta }}J,\end{eqnarray} \tag{ 29 }$$

where ${\rm{\Delta }}J={J}_{x}-{J}_{z}$ is the anisotropy parameter.

A number of methods have been developed to detect genuinely multipartite entanglement in Dicke states [40–46], of which [44] is the most generally applicable one. It consists of a nonlinear witness, i.e. an nonlinear inequality that has to hold for every biseprable state. In order apply the witness only a polynomial (in the number of qubits) number of local measurements is necessary.

Combining [44] with the use of proxy witnesses, entanglement can be detected even more easily experimentally. To do so, let us introduce a (weaker) linear version of the nonlinear witness given in [44]. It is given by

$$\begin{eqnarray}&&{W}_{m}^{n}=\displaystyle \frac{1}{2}\displaystyle \sum _{(\alpha ,\beta )\in \gamma }(-| {d}_{\beta }\rangle \langle {d}_{\alpha }| -| {d}_{\alpha }\rangle \langle {d}_{\beta }| +| {d}_{\alpha \cap \beta }\rangle \langle {d}_{\alpha \cap \beta }| +| {d}_{\alpha \cup \beta }\rangle \langle {d}_{\alpha \cup \beta }| ),\end{eqnarray} \tag{ 30 }$$

where $\gamma =\{(\alpha ,\beta ):| \alpha \cap \beta | =m-1\}$. Wⁿ_m detects $| {D}_{m}^{n}\rangle$ in a maximal way. Using these witnesses as a constraint in (6), it is possible to detect entanglement for different choices of J_x and J_z and B satisfying (29). Results for 11 and 13 qubits are presented in table 7 and figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7.  Detection energies and entropies for the XXZ Hamiltonian with J_z = 1 using (30) as witness.

N	m	B	${\rm{\Delta }}J$	${E}_{0}/N$	${E}_{\mathrm{min},\mathrm{gap}}/N$	${E}_{\mathrm{max},\mathrm{gap}}/N$	${T}_{\mathrm{min},\mathrm{gap}}$	${T}_{\mathrm{max},\mathrm{gap}}$	${\mathcal{E}}$
11	5	−1	10	−13.753	−13.599	−12.890	3.970	7.930	0.031
11	5	−3	20	−26.707	−26.299	−24.698	8.460	16.420	0.038
11	5	−1	20	−26.525	−26.344	−24.802	6.940	15.850	0.033
13	6	−1	12	−16.276	−16.149	−15.577	3.970	7.930	0.021

In order to extend the above results to the thermodynamic limit, let us now focus on witnesses, which are diagonal in the energy eigenbasis. A general example of such a witness is given as $W=\alpha {\mathbb{1}}-| {E}_{0}\rangle \langle {E}_{0}|$, where $| {E}_{0}\rangle$ is the ground state of H and $\alpha ={\mathrm{max}}_{| \phi \rangle ,| \varphi \rangle }{| \langle {E}_{0}| | \phi \rangle | \varphi \rangle | }^{2}$ [47]. Instead of $| {E}_{0}\rangle \langle {E}_{0}|$, one could also use any projector onto one or several energy eigenstates. The only obvious requirement being that the states are entangled themselves, otherwise the resulting operators would not be witnesses. We can then either maximise over a fixed bipartition, or over all possible ones. In the former case it is possible to detect bipartite entanglement w.r.t. the partition chosen, in the latter case we can detect genuinely multipartite entanglement. Inserting the witness into (16), we see that the exponent becomes diagonal

$$\begin{eqnarray}&&\tilde{{\ell }}(\mu ,\nu )=\mathrm{ln}\;\mathrm{Tr}\;\mathrm{exp}(\displaystyle \sum _{i}(\mu {E}_{i}+\nu (\alpha -{\delta }_{i0}))| {E}_{i}\rangle \langle {E}_{i}| )-\mu E,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&=\mathrm{ln}\displaystyle \sum _{i}\mathrm{exp}(\mu {E}_{i}+\nu (\alpha -{\delta }_{i0}))-\mu E,\end{eqnarray} \tag{ 32 }$$

where E_i and $| {E}_{i}\rangle$ are the energy eigenvalues and eigenvectors, respectively and ${\delta }_{i0}$ is the Kronecker delta. Setting the gradient equal to zero, we obtain

$$\begin{eqnarray*}\displaystyle \frac{\partial }{\partial \mu }\tilde{{\ell }}(\mu ,\nu ) & = & 0\Rightarrow ({\rm{i}})\ \displaystyle \sum _{i}\mathrm{exp}(\mu {E}_{i}+\nu (\alpha -{\delta }_{i0}))({E}_{i}-E)=0,\\ \displaystyle \frac{\partial }{\partial \nu }\tilde{{\ell }}(\mu ,\nu ) & = & 0\Rightarrow (\mathrm{ii})\ \displaystyle \sum _{i}\mathrm{exp}(\mu {E}_{i}+\nu (\alpha -{\delta }_{i0}))(\alpha -{\delta }_{i0})=0,\end{eqnarray*}$$

which, by convexity of $\tilde{{\ell }}$, are sufficient conditions for a minimum. Solving those transcendent equations w.r.t. the Lagrange variables is only possible numerically for small N. Still it is possible to obtain a result for the thermodynamic limit, as we will now demonstrate.

Let us first consider the case where $\alpha =1$, i.e. H has a separable ground state. Since the first term in the sum of (ii) vanishes, (ii) cannot be fulfilled for any $\nu \geqslant 0$. Since $\tilde{{\ell }}$ increases with $\nu$, the minimum is attained at $\nu =0$. At $\nu =0$, however, (16) is an upper bound on the unconstrained problem, namely

$$\begin{eqnarray*}&&\underset{\mu }{\mathrm{min}}\mathrm{ln}\displaystyle \sum _{i}\mathrm{exp}(\mu {E}_{i})-\mu E\geqslant S({\rho }_{\mathrm{Gibbs}}).\end{eqnarray*}$$

Since $S({\rho }_{\mathrm{Gibbs}})=\mathrm{ln}Z+\beta E$, the minimum is attained at $\mu =-\beta$. In this case no entanglement is detected.

Let us now move on to the case where $\alpha \lt \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{Z}$. Assuming that the minimum is attained at $\nu =0$, i.e. that (ii) is fulfilled at $\nu \leqslant 0$, (ii) tells us that

$$\begin{eqnarray*}&&\alpha =\displaystyle \frac{\mathrm{exp}(-\beta {E}_{0}+\nu (\alpha -1))}{{\displaystyle \sum }_{i}\mathrm{exp}(-\beta {E}_{i}+\nu (\alpha -{\delta }_{i0}))}=\displaystyle \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{{\displaystyle \sum }_{i}{{\rm{e}}}^{-\beta {E}_{i}}{{\rm{e}}}^{\nu (1-{\delta }_{i0})}}=\displaystyle \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{{{\rm{e}}}^{-\beta {E}_{0}}+{\displaystyle \sum }_{i\gt 0}{{\rm{e}}}^{-\beta {E}_{i}}{{\rm{e}}}^{\nu }}\geqslant \displaystyle \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{Z},\end{eqnarray*}$$

with equality for $\nu =0$, which is a contradiction. Hence the minimum is attained at $\nu \gt 0$. This implies that ${\mathrm{min}}_{\nu \geqslant 0,\mu }\tilde{{\ell }}\lt {\mathrm{min}}_{\mu }\tilde{{\ell }}(0,\mu )=S({\rho }_{\mathrm{Gibbs}})$, hence ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\lt S({\rho }_{\mathrm{Gibbs}})$.

Theorem 4. Let $\mathrm{SEP}$ be the set of separable pure states w.r.t. some partitions and $\alpha ={\mathrm{max}}_{| \phi \rangle \in \mathrm{SEP}}{| \langle {E}_{0}| | \phi \rangle | }^{2}$. Then, if $\alpha \lt \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{Z}$, the Gibbs state as well as any state with entropy ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\lt S\lt S({\rho }_{\mathrm{Gibbs}})$ will be inseparable w.r.t. those partitions.

Note that while for $\alpha \lt \displaystyle \frac{{{\rm{e}}}^{-\beta {E}_{0}}}{Z}$ entanglement of the Gibbs state can also be detected by applying the witness directly, the same is not true for states with entropy ${S}_{\mathrm{max},{{\rm{\Lambda }}}_{A},{W}_{i},E}\lt S\lt S({\rho }_{\mathrm{Gibbs}})$. Since such states always exist, this shows that our method generically works. Theorem 4 can be applied in the thermodynamic limit. Since E₀ generally scales linear in N while Z does so exponentially, theorem 4 can yield non-trivial results where $\alpha ={\mathcal{O}}({e}^{-N})$.

We consider a special case of the Heisenberg Hamiltonian, given by

$$\begin{eqnarray*}&&H=-\displaystyle \sum _{i=1}^{N}(\displaystyle \frac{1+r}{2}{\sigma }_{i}^{x}{\sigma }_{i+1}^{x}+\displaystyle \frac{1-r}{2}{\sigma }_{i}^{y}{\sigma }_{i+1}^{y}+h{\sigma }_{i}^{z}).\end{eqnarray*}$$

This is known as the XY model in a transverse magnetic field. If r = 1, we talk about the Ising model, if r = 0, the XX model. The XY model undergoes a phase transition at h = 1 [2]. The ground state energy as well as the partition function have been computed in the thermodynamic limit by [48]. According to a conjecture numerically tested in [49], $\alpha$ w.r.t. full separability is given by

$$\begin{eqnarray*}&&\underset{N\to \infty }{\mathrm{lim}}\displaystyle \frac{\mathrm{ln}\alpha }{N}=2\underset{\xi }{\mathrm{max}}{\displaystyle \int }_{0}^{\frac{1}{2}}{\rm{d}}\mu \mathrm{ln}| \mathrm{cos}\theta {\mathrm{cos}}^{2}\displaystyle \frac{\xi }{2}+\mathrm{sin}\theta {\mathrm{sin}}^{2}\displaystyle \frac{\xi }{2}\mathrm{cot}\pi \mu | ,\end{eqnarray*}$$

where $\mathrm{tan}2\theta =\frac{r\mathrm{sin}2\pi \mu }{h-\mathrm{cos}2\pi \mu }$ and $-\frac{\pi }{2}\leqslant \theta \leqslant \frac{\pi }{2}$. Using this, we can numerically show entanglement for a wide range of T and h. See figure 8. In particular we can show entanglement in vicinity of the phase transition, in accordance with [2].

Zoom In Zoom Out Reset image size

Unfortunately the entropy gap closes in the limit when using the regular ground state witnesses, a fact which is related to the lack of robustness of the ground state witness itself. It would be most interesting to find other, asymptotically robust witnesses, which could result in a macroscopic entropy gap. Let us now present an artificial example of a Hamiltonian which, with the right choice of witness, allows for an arbitrarily big entropy gap in the thermodynamic limit, thus proving the feasibility of this endeavour in principle.

The Hamiltonian is given by

$$\begin{eqnarray}&&H=\displaystyle \sum _{k=0}^{{d}^{2}-1}k| {{\rm{\Psi }}}_{k}\rangle \langle {{\rm{\Psi }}}_{k}| ,\end{eqnarray} \tag{ 33 }$$

where $d={2}^{\displaystyle \frac{n}{2}}$ and $\{| {{\rm{\Psi }}}_{k}\rangle \}{}_{k=0}^{{d}^{2}-1}$ are a Bell state basis of an n-qubit Hilbert space. The basis is ordered in such a way that $| {{\rm{\Psi }}}_{k}\rangle =\frac{1}{\sqrt{d}}{\displaystyle \sum }_{i=0}^{d-1}| i\rangle | i+k\rangle$ for $k=0...d-1$. Note that these first d basis elements have orthogonal support. As witness, we choose

$$\begin{eqnarray}&&W=\alpha {\mathbb{1}}-\displaystyle \sum _{k=0}^{d-1}| {{\rm{\Psi }}}_{k}\rangle \langle {{\rm{\Psi }}}_{k}| ,\end{eqnarray} \tag{ 34 }$$

which, again, is diagonal in the energy eigenbasis. α is the defined as the maximal overlap of the projector ${\displaystyle \sum }_{k=0}^{d-1}| {{\rm{\Psi }}}_{k}\rangle \langle {{\rm{\Psi }}}_{k}|$ with a separable state. As all the $| {{\rm{\Psi }}}_{k}\rangle$ have orthogonal support, α is the maximal overlap of any of the $| {{\rm{\Psi }}}_{k}\rangle$ with a separable state, hence $\alpha =\displaystyle \frac{1}{d}$ because the $| {{\rm{\Psi }}}_{k}\rangle$ are maximally entangled. Inserting this into (19), we obtain

$$\begin{eqnarray}&&{\rm{\Delta }}S\geqslant \mathrm{ln}\displaystyle \frac{{\displaystyle \sum }_{k=0}^{{d}^{2}-1}{{\rm{e}}}^{-\beta k}}{{\displaystyle \sum }_{k=0}^{d-1}{{\rm{e}}}^{-\beta k}{e}^{-\nu (1-\frac{1}{d})}+{\displaystyle \sum }_{k=d}^{{d}^{2}-1}{{\rm{e}}}^{-\beta k}{e}^{\frac{\nu }{d}}}.\end{eqnarray} \tag{ 35 }$$

Letting n, hence d, go to infinity, we obtain

$$\begin{eqnarray}&&\underset{n\to \infty }{\mathrm{lim}}{\rm{\Delta }}S\geqslant \nu ,\end{eqnarray} \tag{ 36 }$$

where we have used the fact that ${\mathrm{lim}}_{d\to \infty }{\displaystyle \sum }_{k=0}^{{d}^{2}-1}{{\rm{e}}}^{-\beta k}={\mathrm{lim}}_{d\to \infty }{\displaystyle \sum }_{k=0}^{d-1}{{\rm{e}}}^{-\beta k}=\frac{{{\rm{e}}}^{\beta }}{{{\rm{e}}}^{\beta }-1}$. Note that ν can be chosen arbitrarily large.

In macroscopic systems assuming the exact form of the Hamiltonian is always an idealisation. While it is fair from a physical point of view one might wonder about the impact of mis-characterised Hamiltonians on our entropy gap ${\rm{\Delta }}S$. If we assume that the real Hamiltonian is given as $\tilde{H}=H+\epsilon P$ we can that the dual entropy of the actual Hamiltonian for the chosen parameters

$$\begin{eqnarray}&&{\tilde{S}}_{\mathrm{max},W,E}(\mu ,\nu )=\mathrm{log}\;\mathrm{Tr}\;\mathrm{exp}(-\mu H-\mu \epsilon P+\nu W)+\mu E+\mu \epsilon \mathrm{Tr}(\rho P),\end{eqnarray} \tag{ 37 }$$

can be bounded from above using the Golden–Thompson $\mathrm{Tr}({e}^{A+B})\leqslant \mathrm{Tr}({e}^{A})\mathrm{Tr}({e}^{B})$ and Hölder's ${| | {AB}| | }_{1}\leqslant {| | A| | }_{\infty }{| | B| | }_{1}$ inequalities to yield

$$\begin{eqnarray}&&{\tilde{S}}_{\mathrm{max},W,E}(\mu ,\nu )\leqslant {S}_{\mathrm{max},W,E}(\mu ,\nu )+\mu \epsilon \mathrm{Tr}(\rho P)-\mu \epsilon {P}_{o}.\end{eqnarray} \tag{ 38 }$$

This shows that small perturbations or inaccuracies in the description of the Hamiltonian will only have a correspondingly small impact on the validity of the entanglement certification by proxy witnesses.

While the mean energy and the entropy are both macroscopic properties of quantum states, measurements of the latter are not possible directly as they do not correspond to a quantum observable. There are however various ways that the global entropy can be determined, and all we need is a lower bound on the entropy. The most straightforward way would of course consist of equilibrating the system with a thermal bath at temperature T. This is naturally the case, as any system found in nature at ambient temperature T, will most accurately be described by its corresponding Gibbs state. Thus, once the system is equilibrated we know that its corresponding entropy should correspond to the Gibbs entropy $S({\rho }_{\mathrm{th}}(\beta (E)))$. Starting from equilibrated systems one can introduce global quenches of system parameters, such that the entropy will still be bounded from below by the initial Gibbs entropy due to the second law of thermodynamics, despite the system being far out of equilibrium. In that way one directly receives a lower bound on the system entropy for a wide range of out of equilibrium systems and experimental preparations. There are of course other methods, such as reasonable assumptions about symmetry in the state, that can be used to estimate linear entropy from macroscopic spin observables (as would be possible e.g. in [50]). Another straightforward way to obtain lower bounds would be access to a d × d sub-matrix of the global state, whose entropy will always yield a lower bound to the global entropy. The exact ways of experimentally estimating entropies is of course highly dependent on the experimental setup, access and reasonable assumptions about system properties. As any non-zero amount of entropy already provides an advantage in entanglement detection it is however fairly straightforward to infer correspondingly useful entropies.

This highlights the main advantage of our approach: while entanglement of Gibbs states could be directly inferred from its description, it would strictly work only if we know that the state is in fact in thermal equilibrium. This would require strict assumptions about the state of the system and exact characterisation of its Hamiltonian. Using our entropy gap there is no need whatsoever to assume any particular form of the system's state: measurements of the mean energy and lower bounds on the global entropy are completely sufficient to prove that the underlying system is entangled in a robust way. We believe that future work will uncover witnesses robust enough to be amenable to the proxy method, yielding a sizeable entropy gap in the thermodynamic limit.