Structure of Quantum Entanglement at a Finite Temperature Critical Point

TL;DR

This paper investigates the nature of quantum entanglement near finite temperature critical points, revealing that certain long-range correlations are classical and introducing a quantum correlation length that remains finite at the transition.

Contribution

It introduces a scheme using tripartite entanglement negativity to characterize quantum entanglement at finite temperature critical points and demonstrates its classical nature in a mean-field Ising model.

Findings

Tripartite negativity shows no singularity across the transition.

Long-distance critical fluctuations are purely classical.

Mixed bosonic states with certain symmetries have zero entanglement.

Abstract

We propose a scheme to characterize long-range quantum entanglement close to a finite temperature critical point using tripartite entanglement negativity. As an application, we study a model with mean-field Ising critical exponents and find that the tripartite negativity does not exhibit any singularity in the thermodynamic limit across the transition. This indicates that the long-distance critical fluctuations are completely classical, allowing one to define a `quantum correlation length' that remains finite at the transition despite a divergent physical correlation length. Motivated by our model, we also study mixed state entanglement in tight-binding models of bosons with $U(1)$ and time-reversal symmetry. By employing Glauber-Sudarshan `P-representation', we find a surprising result that such states have zero entanglement.

Equations

$$E_N=E_{N,\text{local}}+E_{N,\text{non-local}},$$

$$E_N(A)=\log \text{norm}{{\rho }^{T_A}}_1=\log \text{norm}{({\rho }^{T_{\bar{A}}}{)}^T}_1=\log \text{norm}{{\rho }^{T_{\bar{A}}}}_1,$$

$$H=\frac{1}{2}\sum _{\mathbf{r}}\left({\pi }_{\mathbf{r}}^2+m^2{\varphi }_{\mathbf{r}}^2\right)+\frac{1}{2}\sum _{\text{expectationvalue}\mathbf{r},{\mathbf{r}}^{\prime }}K{\left({\varphi }_{\mathbf{r}}-{\varphi }_{{\mathbf{r}}^{\prime }}\right)}^2,$$

$$m(g)=\{\begin{array}{l}\sqrt{g-g_c}\text{for}g>g_c\\ \sqrt{2(g_c-g)}\text{for}g<g_c.\end{array}$$

$$E_{N,\text{local}}=\left[\alpha +b_{\pm }|t|\right]L^{d-1}$$

$${\Delta }_3{E}_N=E_{N,A}+E_{N,B}+E_{N,C}-E_{N,AB}-E_{N,BC}-E_{N,CA}+E_{N,ABC}$$

$${\Delta }_3{E}_N\sim e^{-L/{\xi }_Q}$$

$$H_{1D,k_y}=\frac{1}{2}\sum _x\left({\pi }_x^2+m_{k_y}^2{\varphi }_x^2\right)+\frac{1}{2}\sum _{\text{expectationvalue}x,x^{\prime }}K{\left({\varphi }_x-{\varphi }_{x^{\prime }}\right)}^2,$$

$$E_{N,k_y=0}=\frac{1}{2}\ln \left[\frac{\beta }{\pi a}\sinh \left(\frac{\pi L_A}{\beta }\right)\right]-\frac{\pi L_A}{2\beta }+f(e^{-2\pi L_A/\beta })+2\ln c_{1/2},$$

$$=\left[\frac{1}{2}\log \left(\frac{\beta }{2\pi a}\right)+c_{1/2}\right]+\left[\log \left(1-e^{-2\pi L_A/\beta }\right)+f(e^{-2\pi L_A/\beta })\right].$$

$$H=-\sum _{ij}{t}_{ij}{a}_i^{†}{a}_j,$$

$$\rho ={\int }_{{\mathbb{C}}^{\mathbb{N}}}\prod _{i=1}^N{d}^2{\alpha }_{i}P(\mathbf{\alpha })|\mathbf{\alpha }⟩⟨\mathbf{\alpha }|,$$

$$E_{N,\text{local, singular}}\sim |t{|}^{1-\alpha }{L}^{d-1}$$

$$H=\frac{1}{2}\int d^{d}x{m}^2\varphi (x)+{\left(\nabla \varphi (x)\right)}^2,$$

$$Z=\sqrt{\frac{{\left(2\pi T\right)}^V}{{\prod }_k{\lambda }_k}},$$

$$F=-T\log Z=-\frac{1}{2}VT\log (2\pi T)+\frac{1}{2}T\log \left[\prod _k{\lambda }_k\right],$$

$$S=\frac{V}{2}\left[1+\log \left(2\pi T\right)\right]-\frac{1}{2}\log \left[\prod _k{\lambda }_k\right].$$

$$\prod _{n=0}^{L-1}\left[a^2+4{\sin }^2(n\frac{\pi }{L})\right]=g(a),$$

$$\prod _k{\lambda }_k=\prod _{\begin{array}{c}{k}_1,\cdots ,k_d\end{array}}{\lambda }_k=\prod _{\begin{array}{c}{k}_1,\cdots ,k_{d-1}\end{array}}g\left(\sqrt{m^2+4\sum _{n=1}^{d-1}{\sin }^2(\frac{1}{2}{k}_i)}\right).$$

$$\prod _k{\lambda }_k=g(m)\prod _{(k_1,\cdots ,k_{d-1})\ne 0}g\left(\sqrt{m^2+4\sum _{n=1}^{d-1}{\sin }^2(\frac{1}{2}{k}_i)}\right).$$

$$\log \left[\prod _k{\lambda }_k\right]=2\log (mL)+\sum _{(k_1,\cdots ,k_{d-1})\ne 0}\log \left[g\left(2\tilde{\omega }\right)\right],$$

$$\log \left[\prod _k{\lambda }_k\right]\to 2\log (mL)+2L^d\int \frac{d^{d-1}k}{(2\pi {)}^{d-1}}\log \left[\sqrt{{\tilde{\omega }}^2+1}+\tilde{\omega }\right].$$

$$H=\frac{1}{2}\sum _{\mathbf{r}}\left({\pi }_{\mathbf{r}}^2+m^2{\varphi }_{\mathbf{r}}^2\right)+\frac{1}{2}\sum _{\text{expectationvalue}\mathbf{r},{\mathbf{r}}^{\prime }}K{\left({\varphi }_{\mathbf{r}}-{\varphi }_{{\mathbf{r}}^{\prime }}\right)}^2,$$

$$\begin{array}{rl}& {\gamma }_{\varphi }(\mathbf{r},{\mathbf{r}}^{\prime })=2\text{expectationvalue}{\varphi }_{\mathbf{r}}{\varphi }_{{\mathbf{r}}^{\prime }}=\frac{1}{N}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)}\frac{1}{{\omega }_{\mathbf{k}}}\coth (\frac{1}{2}\beta {\omega }_{\mathbf{k}})\\ & {\gamma }_{\pi }(\mathbf{r},{\mathbf{r}}^{\prime })=2\text{expectationvalue}{\pi }_{\mathbf{r}}{\pi }_{{\mathbf{r}}^{\prime }}=\frac{1}{N}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)}{\omega }_{\mathbf{k}}\coth (\frac{1}{2}\beta {\omega }_{\mathbf{k}}),\end{array}$$

$$P(\mathbf{r},{\mathbf{r}}^{\prime })=\{\begin{array}{l}{\delta }_{\mathbf{r},{\mathbf{r}}^{\prime }}\text{for}\mathbf{r}\in A\\ -{\delta }_{\mathbf{r},{\mathbf{r}}^{\prime }}\text{for}\mathbf{r}\in B,\end{array}$$

$$E_N=-\frac{1}{2}\sum _{i=1}^N\min \{0,\log {\lambda }_i\},$$

$$\psi (\mathbf{r})=e^{i{\sum }_{i=2}^d{k}_i{r}_i}{\psi }_1(r_1).$$

$$E_N=\sum _{k_2,\cdots k_d}{E}_N^{(1)}(k_2,\cdots ,k_d),$$

$$m^2\to m^2+4K\sum _{i=2}^d\sin [2](\frac{1}{2}{k}_i).$$

$$\begin{array}{rl}& a_{\pm }=\frac{1}{2}\left[\frac{1}{m}\coth (\frac{1}{2}\beta m)\pm \frac{1}{\sqrt{m^2+\frac{1}{4}}}\coth (\frac{1}{2}\beta \sqrt{m^2+4})\right]\\ & b_{\pm }=\frac{1}{2}\left[m\coth (\frac{1}{2}\beta m)\pm \sqrt{m^2+\frac{1}{4}}\coth (\frac{1}{2}\beta \sqrt{m^2+4})\right].\end{array}$$

Full text

I Introduction

Qualitatively, there are two distinct classes of phase transitions in quantum mechanical Hamiltonians: those that occur at the absolute zero temperature, and those at a finite (i.e. non-zero) temperature. Heuristically, the zero temperature phase transitions result due to quantum fluctuations while the finite temperature ones typically result from thermal fluctuations. For example, consider the transverse field Ising model on a $d$-dimensional hypercubic lattice, $H_{\textrm{TFI}}=-\sum_{<i,j>}Z_{i}Z_{j}-h\sum_{i}X_{i}$. This Hamiltonian supports two phases: a ferromagnetic phase and a paramagnetic phase. The critical exponents associated with the zero temperature transition belong to the $d+1$-dimensional Ising universality while those for the finite temperature transition belong to the $d$-dimensional Ising universality i.e. at finite temperature, one may as well set $h=0$ to obtain the critical exponents Suzuki (1976); Sachdev (2011). This is consistent with the conventional wisdom that quantum mechanics does not play any role in the long-distance equilibrium physics of finite temperature phase transitions.

However, there also exist models such as the four-dimensional toric code Dennis et al. (2002) which host finite temperature ‘quantum memory’, that is, one can encode a qubit non-locally in this model at finite temperature such that it is well protected for an infinite time even when coupled to a heat bath. This model also exhibits a finite temperature phase transition across which the quantum memory is destroyed. Therefore, one suspects that in this model, both at and below the critical temperature there exist intrinsically quantum effects even at long-distances. This raises the question: is there a quantity that sharply distinguishes the finite temperature transition in a transverse field Ising model from that in the 4D Toric code?

A related question is: how easy is it to prepare thermal (mixed) states on a classical computer? It has been argued that a thermal state can be prepared efficiently if the system does not possess a finite temperature quantum memory, and has a short correlation length Swingle and McGreevy (2016a, b); Brandão and Kastoryano (2019). Similar results have been argued for the preparation of thermofield double state which corresponds to a purification of a thermal density matrix Wu and Hsieh (2018); Martyn and Swingle (2018); Cottrell et al. (2019); Maldacena and Qi (2018); Chapman et al. (2019). However, if the long-range correlations in a system arise purely due to classical effects (e.g., consider $H_{\textrm{TFI}}$ at the finite temperature critical point for $|h|\ll 1$), one might wonder if the corresponding state can again be prepared efficiently?

The above discussion motivates us to ask: How does one separate quantum-mechanical correlations from classical correlations at a finite temperature, and in particular, in the vicinity of a phase transition? At zero temperature, a system can typically be described by a pure state and correspondingly, the von Neumann entropy of a reduced density matrix corresponding to a subsystem is a faithful measure of long-range entanglement in the critical ground state Holzhey et al. (1994); Calabrese and Cardy (2004); Ryu and Takayanagi (2006); Casini and Huerta (2007); Metlitski et al. (2009). In contrast, at a finite temperature $T$, the system is described by a thermal (i.e. Gibbs) state $\propto e^{-H/T}$, which is a mixed density matrix. To probe an intrinsic quantum correlation at finite temperature, one must therefore resort to an entanglement measure for mixed states Werner (1989); Horodecki et al. (2009). To this end, here we will employ ‘entanglement negativity’ (henceforth just ‘negativity’ for brevity) which has the property that it is an entanglement monotone and unlike most other mixed-state measures, does not require optimizing a function over all possible quantum states Eisert and Plenio (1999); Vidal and Werner (2002); Plenio (2005); Horodecki et al. (2009). As shown in Ref. Sherman et al. (2016), one very interesting property of negativity is that for thermal states of local Hamiltonians it satisfies an ‘area-law’, akin to the von Neumann entanglement entropy of (pure) ground states of gapped Hamiltonians Bombelli et al. (1986); Srednicki (1993); Eisert et al. (2010), and in strong contrast to the volume law for pure finite energy density eigenstates Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Santos et al. (2012); Deutsch (2010); D’Alessio et al. (2016); Garrison and Grover (2018).

A recent work, Ref.Lu and Grover (2019), constructed a class of exactly solvable models which host finite-$T$ order-disorder transition and for which negativity can be calculated analytically. It was found that in all models considered, whenever negativity is non-zero in the vicinity of the transition, it is a singular function of the tuning parameter driving the transition, despite the fact that these phase transitions are driven purely by thermal fluctuations and do not host finite-$T$ quantum memory. As argued in Ref.Lu and Grover (2019), the area law coefficient of negativity receives contribution from the expectation value of local operators close to the entangling boundary, and since the expectation value of local operators is singular across the transition, the area-law coefficient is singular as well.

The aforementioned singularity in area-law coefficient across a finite-$T$ transition leaves open the question of extracting purely quantum correlations that are sensitive only to long-distance physics, unlike the area-law coefficient of negativity which certainly depends on short-distance physics despite containing information about critical exponents. This is the topic of this paper. Inspired by the methods used to extract universal entanglement encoded in the ground states of gapped Hamiltonians Kitaev and Preskill (2006); Levin and Wen (2006), we propose a tripartite negativity to probe long-distance, universal quantum correlations at finite temperature. We study the tripartite negativity, denoted as $\Delta_{3}E_{N}$ below, for a simple model that exhibits singularity in correlation functions, von Neumann entropy as well as the area-law coefficient of negativity across a finite-$T$ transition. We find that $\Delta_{3}E_{N}$ completely cancels out the aforementioned singularity associated with the transition, and is exponentially small in the system size $\Delta_{3}E_{N}\sim e^{-L/\xi_{Q}}$ where $\xi_{Q}$ defines a ‘quantum correlation length’ which, in contrast to the physical correlation length, does not diverge at the finite-$T$ transition (Fig.1). As $T\rightarrow 0$, $\xi_{Q}$ diverges resulting in a non-zero $\Delta_{3}E_{N}$ which corresponds to the universal non-zero subleading term for Renyi entropy $S_{1/2}$ at the quantum phase transition. Note that at $T=0$, for a gapped, topological ordered phase $\Delta_{3}E_{N}$ also equals the topological entanglement entropy Lee and Vidal (2013); Castelnovo (2013).

Partly inspired by our model, we also study mixed state entanglement in tight-binding models of free bosons with time-reversal and $U(1)$ symmetry. We show that the corresponding thermal state is separable, and therefore, any measure of mixed state entanglement for such a state, such as entanglement of negativity or entanglement of formation, is zero.

The paper is organized as follows: In Sec. II.1, we discuss the general structure of negativity for local Hamiltonians. In Sec. II.2, we introduce our scheme for calculating the universal part of negativity, and implement it for a model that shows a finite temperature transition. In Sec. II.3 we discuss cross-over towards the zero temperature quantum phase transition, and in Sec.II.4 we discuss the eigenvalues and eigenfunctions of correlation matrix that determines negativity. In Sec. III, we discuss our aforementioned result on separability of bosonic Gaussian states with $U(1)$ and time-reversal symmetry. In Sec. IV, we conclude with a summary and possible implications of our results.

II Negativity across finite-T critical points

II.1 General structure of negativity: local versus non-local contributions

Given a state $\rho$ acting on the Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{\bar{A}}$, the negativity is defined as $E_{N}(A)=\log\norm{\rho^{T_{A}}}_{1}$. Due to the area-law of negativity for thermal states of local Hamiltonians (Ref.Sherman et al. (2016)), the problem to characterize the universal part of their negativity is somewhat analogous to the characterization of long-distance entanglement in ground states of gapped Hamiltonians. Following Ref.Grover et al. (2011), we consider a coarse-grained, continuum description, and write

[TABLE]

where $E_{N,\textrm{local}}$ is expressible as a sum of local terms along the entangling surface: $E_{N,\textrm{local}}=\int_{\partial A}F(\{\kappa,\partial_{i}\kappa,...,\})$ where $\kappa$ is the local curvature along the entangling surface. Similar to the von Neumann entropy of pure states, negativity of mixed states satisfies $E_{N}(A)=E_{N}(\bar{A})$. This is because

[TABLE]

where the last equality results from invariance of eigenspectrum under matrix transpose operation. Since under the exchange $A\leftrightarrow\bar{A}$, the curvature $\kappa\leftrightarrow-\kappa$, $F$ must be an even functional of the curvature $\kappa$, and consequently $E_{N,\textrm{local}}=\alpha_{d-1}L_{A}^{d-1}+\alpha_{d-3}L_{A}^{d-3}+...$, i.e., only alternate terms in the expansion in terms of $L_{A}$ are allowed, again similar to the discussion of gapped ground states (Ref.Grover et al. (2011)). One implication of this is that in two-dimensions, a non-zero constant term $\gamma$ in $E_{N}\sim L_{A}-\gamma$ necessarily implies a non-zero $E_{N,\textrm{non-local}}$.

Again motivated by the theory of gapped ground states Kitaev and Preskill (2006); Levin and Wen (2006), below we use a subtraction scheme to cancel out $E_{N,\textrm{local}}$ to understand the behavior of $E_{N,\textrm{non-local}}$, and specifically whether it has a non-zero value at a finite-$T$ phase transition in the thermodynamic limit. It’s worth emphasizing that the coefficients $\alpha_{i}$ that enter $E_{N,\textrm{local}}$ will generically be singular functions of the tuning parameter driving the transition. This is because these coefficients will depend on the expectation value of local operators, such as energy density, which themselves are a singular function of the tuning parameter. This leads to the singularity in the area-law coefficient of negativity as discussed in Ref.Lu and Grover (2019), and which we will again encounter below.

It is also important to note that for pure states which exhibit power-law correlations, such as a ground state corresponding to a conformal field theory, or a Fermi surface, the non-local part of entanglement is necessarily non-zero Callan and Wilczek (1994); Holzhey et al. (1994); Calabrese and Cardy (2004); Gioev and Klich (2006); Wolf (2006); Swingle (2010). In contrast, for thermal (mixed) states of quantum systems that display power-law correlations, such as those corresponding to a finite-$T$ phase transition, the entanglement structure close to the transition remains completely unexplored. To that end, we now turn to studying negativity in a specific model that displays a finite-$T$ phase transition.

The aforementioned curvature expansion for negativity relies on a coarse-grained continuum description. For such a description to be valid, one requires that all length scales involved are much larger than the short-distance lattice-cutoff $a$. Close to a finite temperature critical point, the physical correlation length $\xi$ of course satisfies $\xi\gg a$, but as hinted above, we will also encounter a length scale $\xi_{Q}$ that does not diverge. As we will see below, $\xi_{Q}\gg a$ at low-temperature, and therefore, the argument is valid over a range of temperatures. This is similar to the discussion in the context of gapped ground states Grover et al. (2011) where the correlation length is assumed to stay large compared to the lattice cut-off.

II.2 Universal negativity of a model with finite-T phase transition

We consider a $d$ dimensional cubic lattice of $N=L^{d}$ sites, where a site at $\mathbf{r}$ is associated with a degree of freedom described by a canonically conjugate pair $(\phi_{\mathbf{r}},\pi_{\mathbf{r}})$. The Hamiltonian reads

[TABLE]

where the physical mass obeys

[TABLE]

$g_{c}$ sets the critical point where correlation length diverges due to the vanishing mass term. The motivation to study this specific model comes from the transverse-field Ising model: $H_{\textrm{TFI}}=-\sum_{<\mathbf{r},\mathbf{r}^{\prime}>}Z_{\mathbf{r}}Z_{\mathbf{r}^{\prime}}-h\sum_{\mathbf{r}}X_{\mathbf{r}}$. Both $H$ and $H_{TFI}$ have $\mathbb{Z}_{2}$ symmetry, and $H$ can be thought of as a mean-field approximation to $H_{TFI}$ where $\pi_{\mathbf{r}}^{2}/2$ plays the role of the transverse field $X_{\mathbf{r}}$, and $\phi_{\mathbf{r}}$ plays the role of $Z_{\mathbf{r}}$. The function $m(g)$ corresponds to the effective mass (= inverse correlation length) within the mean-field theory. Close to $g_{c}$, at finite temperature, the system exhibits long-range correlation functions. For example, in 2d, at $g_{c}$ $\langle\phi_{\mathbf{r}}\,\phi_{\mathbf{r}^{\prime}}\rangle\sim\log(|\mathbf{r}-\mathbf{r}^{\prime}|)$. Another signature of the divergent correlation length is that the von Neumann entropy $S=-\tr\left(\rho\log\rho\right)$ shows a logarithmic divergence as approaching to $g_{c}$: $S\sim-\log(mL)$Gelfand and Fisher (1988) as shown in Fig.2(b) ( see Appendix A for derivation, and footnote foo regarding how the massless limit is taken). Note that, instead of tuning $g$, one can also tune the temperature to drive the finite-$T$ transition (see Fig.2(a)). All of our results below are unchanged if one simply replaces $g-g_{c}$ by $(T-T_{c})/T_{c}$.

Next, we analyze the structure of quantum correlations close to the transition using entanglement negativity. Since $H$ is quadratic, a thermal state $\rho$ at inverse temperature $\beta$, i.e. $\rho\sim e^{-\beta H}$, is a Gaussian state. This allows an efficient calculation of negativity Audenaert et al. (2002); Ferraro et al. (2008); Anders (2008); Cavalcanti et al. (2008); Marcovitch et al. (2009) because one can perform partial transposition on $\rho$ using the covariance matrix technique. For concreteness, we consider a two dimensional lattice ($d=2$), and set $K=1$. We first plot the negativity of a subregion $A$ and its complement as a function of $(g-g_{c})$ in Fig.3(a). In agreement with the result of Ref.Lu and Grover (2019), we find that the area law coefficient of the negativity is singular at $g=g_{c}$. As mentioned above, this singularity originates from the singular behavior of the expectation value of local operators close to the entangling boundary Lu and Grover (2019). What is the precise nature of this singularity? We numerically find that in the limit $a<<\beta<<1/m$, at the lowest order, the singular part of area-law coefficient is proportional to $m^{2}$. Since $m^{2}\sim|t|$ where $t$ is the thermal tuning parameter (either $T-T_{c}$ or $g-g_{c}$),

[TABLE]

where $\alpha$ is an analytic function of underlying parameters, and $b_{+}~{}(b_{-})$ is coefficient of $|t|$ for $g>g_{c}$ ($g<g_{c}$) with $b_{-}/b_{+}=2$ due to the singular dependence of mass close to the transition (Eq. 4). As shown in Appendix.C, one can gain some intuition for the $m^{2}$ dependence by analytically studying the negativity between two sites. We comment on the relation of this singularity to critical exponents in Sec.IV.

To isolate long-distance quantum correlations, we now define a quantity analogous to ‘topological entanglement entropy’ in the context of ground-states of gapped Hamiltonians Kitaev and Preskill (2006); Levin and Wen (2006) that cancels out the $E_{N,\textrm{local}}$ component of the negativity. Consider three subregions $A,B$ and $C$, and define a tripartite negativity,

[TABLE]

Here $E_{N,A}$ denotes the negativity between the region $A$ and its complement, and similarly for all the other quantities present in $\Delta_{3}E_{N}$. On a square lattice of size $L$ with the periodic boundary condition, we consider the partition shown in the inset of Fig. 3(b), where $A,B$ and $A\bigcup B\bigcup C$ are squares of size $2/5L,2/5L$ and $4/5L$ respectively. Fig. 3(b) shows the dependence of $\Delta_{3}E_{N}$ on $g-g_{c}$ for various system sizes. Despite the fact that each of the seven individual terms that enters the definition of $\Delta_{3}E_{N}$ (Eq.6) is singular (see Fig.3(a)), one finds that $\Delta_{3}E_{N}$ itself does not exhibit any singularity across $g_{c}$, upto terms that are exponentially small in the total system size $L$. In fact, $\Delta_{3}E_{N}$ itself vanishes exponentially with $L$ at all non-zero temperatures (see Fig.3(c)):

[TABLE]

which defines a ‘quantum correlation length’ $\xi_{Q}$ that remains finite even at the critical point. The peak in $\Delta_{3}E_{N}$ at the critical point (Fig.3(b)) indicates that $\xi_{Q}$ is largest at the critical point. We discuss the detailed behavior of $\xi_{Q}$ below for the case of a straight bipartition without any corners. It is worth emphasizing that the singularity in quantities that are sensitive both to classical and quantum correlations survives in the thermodynamic limit after an analogous subtraction scheme. For example, an analogously defined tripartite von Neumann entropy $\Delta_{3}S$ continues to show singular behavior identical to Fig.2(b).

A partial analytical understanding of finiteness of $\xi_{Q}$ at the critical point is provided by considering a bipartition without any curvature or corners by dividing a torus of size $L\times L$ into two strips of equal size $L/2\times L$. Based on the discussion above, we expect that for such bipartition, $E_{N,\textrm{local}}\propto L^{d-1}_{A}$, i.e., it is strictly an area-law without any subleading corrections. The universal part of the negativity can now be extracted by studying $\Delta_{2}E_{N}=E_{N}(2L)-2E_{N}(L)$ which cancels out the aforementioned $E_{N,\textrm{local}}$ contribution, leaving only $E_{N,\textrm{non-local}}$, the term of interest. The subtraction scheme based on $\Delta_{2}E_{N}$ is conceptually quite similar to that based on $\Delta_{3}E_{N}$, the former is more suited towards a bipartition without any curvature, while the latter is more general. Note that setting $m=0$ explicitly leads to numerical instability in the diagonalization of the covariance matrix, which we regularize by setting $m=10^{-5}$, and confirm that a further decrease of the mass (while keeping it non-zero) does not change the numerical value of negativity. We find that in this massless limit, $\Delta_{2}E_{N}$ continues to decay exponentially with system size $L$ for all finite temperatures (Fig.4(a)), similar to the behavior of $\Delta_{3}E_{N}$ studied above. Furthermore, we find that the quantum correlation length scale $\xi_{Q}$ is roughly proportional to the inverse temperature $\beta$ as shown in the inset. Using conservation of momentum along the $\hat{y}$ direction, one can express the negativity for this bipartition as the sum of negativities corresponding to the following 1D Hamiltonians $H_{1D,k_{y}}$ labeled by $k_{y}$, the momentum along the $\hat{y}$ direction (see Appendix.B for derivation):

[TABLE]

where $m^{2}_{k_{y}}=4K\sin[2](\frac{1}{2}k_{y})$. One expects that any non-local contribution to negativity can arise only when the mass $m_{k_{y}}=0$, i.e., the contribution of $k_{y}=0$ mode. Curiously, the contribution of $k_{y}=0$ mode is identical to negativity corresponding to the thermal state of a central charge $c=1$ 1+1-D CFT studied in Ref.Calabrese et al. (2015). Using results from Ref.Calabrese et al. (2015), one finds that

[TABLE]

where $f$ is a universal scaling function which tends to a constant when $L_{A}\gg\beta$, $c_{1/2}$ is not universal and $a$ is the lattice constant. When $L_{A}\gg\beta$, the expression can be written as

[TABLE]

This expression implies that when $L_{A}/\beta\rightarrow\infty$, $E_{N,k_{y}=0}$ approaches a constant value over a characteristic length-scale $\beta$ in line with our numerical results for $\Delta_{2}E_{N}$.

We also study the behavior of $\xi_{Q}$ in the vicinity of the transition, and find that it is maximum at the transition and exhibits a cusp singularity (Fig.4(b)). This is expected since $\xi_{Q}$ is a function of the mass $m$ which itself is singular across the transition.

II.3 Approach to Quantum Critical Point

The quantum correlation length diverges as $T\rightarrow 0$, and in the (pure) ground state the negativity equals $S_{1/2}$ where $S_{n}=-\frac{1}{n-1}\log\tr\rho_{A}^{n}$ is the $n$’th Renyi entropy. Since a massless scalar has long-range entanglement in its ground state which is reflected in Renyi entropies as well, one expects that the non-local part of negativity will be non-zero in the ground state. Fig.4(c) shows how the non-local negativity interpolates between its exponentially small value at any non-zero $T$ to a non-zero, universal $O(1)$ value at $T=0$. The scaling collapse of the data when plotted as a function of $e^{-L/\beta}$ again indicates that the quantum correlation length $\xi_{Q}\sim\beta$. We verified that the O(1) constant contribution to negativity at zero temperature agrees with the known result for a massless scalar Witczak-Krempa et al. (2017); Chen et al. (2017).

II.4 Eigenvalue structure of the partial transposed correlation matrix

We find that the eigenvalues and eigenvectors of the correlation matrix that determines negativity have a very specific pattern which reveals more information about the mixed state entanglement. For a Gaussian density matrix such as ours, the negativity is determined by the eigenvalues of the matrix $\sqrt{\gamma_{\phi}P\gamma_{\pi}P}$ which we will denote as the ‘partial transposed correlation matrix’. Here $\gamma_{\phi}(\mathbf{r},\mathbf{r}^{\prime})=2\expectationvalue{\phi_{\mathbf{r}}\phi_{\mathbf{r}^{\prime}}}$, $\gamma_{\pi}(\mathbf{r},\mathbf{r}^{\prime})=2\expectationvalue{\pi_{\mathbf{r}}\pi_{\mathbf{r}^{\prime}}}$, and $P(\mathbf{r},\mathbf{r}^{\prime})=\delta_{\mathbf{r},\mathbf{r}^{\prime}}\quad\text{for}\quad\mathbf{r}\in A$ and $-\delta_{\mathbf{r},\mathbf{r}^{\prime}}\quad\text{for}\quad\mathbf{r}\in\bar{A}$. Specifically, $E_{N}=-\sum_{i=1}^{N}\min\{0,\log\nu_{i}\}$ where $\{\nu_{i}\}$ are the eigenvalues of $\sqrt{\gamma_{\phi}P\gamma_{\pi}P}$. Therefore only eigenvalues less than unity contribute to negativity.

Consider, for instance, the Hamiltonian in Eq.3 in $d=1$ (as discussed above, for a bipartition without any corners, the eigenvalues in $d>1$ can be determined in terms of eigenvalues for the $d=1$ problem). We find that the spectrum consists of a ‘bulk’ continuum part and and two discrete eigenvalues isolated from the continuum spectrum (see Fig.5(a)). We numerically find that the bulk continuum part in fact matches with the eigenvalues of $\sqrt{\gamma_{\phi}\gamma_{\pi}}$, the correlation matrix without any partial transpose operation. Quite strinkingly, as one notes from Fig.5(a), only the isolated discrete eigenvalue less than unity contributes to the negativity.

Given the distinctive nature of eigenvalue spectrum, it is instructive to contrast the eigenfunctions corresponding to the continuum eigenvalues with that for the isolated eigenvalue that contributes to negativity. We find that while the eigenfunctions corresponding to the bulk continuum spectrum essentially behave as a plane wave, the eigenfunction corresponding to the discrete eigenvalue is localized at the bipartition boundary (see Fig.5(b)), signifying the fact that quantum entanglement is localized close to the boundary. Furthermore, the localized eigenfunction decays exponentially away from the entanglement cut even in the massless limit (Fig.5(c)). In particular, the characteristic decay length is proportional to the inverse temperature $\beta$ as indicated by the scaling collapse analysis (see Fig.5(c) inset), quite similar to the aforementioned behavior of the ‘quantum correlation length’ determined using the decay of the non-local negativity.

For plots in Fig.5, we impose open boundary condition. We have checked that these observations apply to periodic boundary condition as well, the only difference being that now the discrete eigenvalues will be two-fold degenerate due to the presence of two entanglement cuts.

III Separability of bosonic Gaussian states with $U(1)$ and Time-reversal symmetries

Above, we studied negativity for a Gaussian density matrix with Ising symmetry. In this section, we report a somewhat surprising observation on mixed state entanglement for a closely related problem. We consider a bosonic tight-binding model in arbitrary spatial dimension:

[TABLE]

where $i$ index labels the lattice sites and $a_{i},a_{i}^{\dagger}$ are the corresponding annihilation and creation bosonic operator. Apart from the $U(1)$ symmetry corresponding to the particle number conservation, we also impose time-reversal symmetry so the hopping amplitude $t_{ij}$ is real. Crucially, we work in grand-canonical ensemble, i.e., we do not fix the number of bosons exactly, but only on average via a chemical potential (which will correspond to diagonal elements of the hopping matrix $t_{ij}$ in $H$). This model is of interest since it exhibits Bose-Einstein condensation (BEC) when cooled below a critical temperature. Therefore, one may wonder how quantum correlations behave across such a finite-$T$ phase transition. Surprisingly, we find that the Gibbs state, i.e. $\rho=\exp(-\beta H)/Z$, for this model can be written as a convex combination of product states and thus is separable. Therefore, all measures of mixed state entanglement, including negativity, are identically zero at all temperatures. It is important to note that negativity can be zero even for non-separable states Horodecki et al. (1998, 2009), and therefore this result is much stronger than just showing that negativity is zero for this system.

The central idea in our proof is to employ the Glauber-Sudarshan ‘$P$ representation’ Sudarshan (1963); Glauber (1963) for the density matrix $\rho$:

[TABLE]

where $\ket{\bm{\alpha}}=\otimes_{i}\ket{\alpha_{i}}$ is a tensor product of coherent states at all sites. $P(\bm{\alpha})$ is a quasiprobability distribution since it can be negative in general. We find that when $\rho$ is a Gibbs state corresponding to the aforementioned tight-binding model, $P(\bm{\alpha})$ is a proper probability distribution function, and thus $\rho$ is separable for all inverse temperature $\beta$. See appendix D for details.

At a first glance, this result seems puzzling since as $T\to 0$, one might expect that $\rho$ will correspond to a pure ground state of $H$, which can be entangled, and is contrary to our finding. This tension is resolved by noticing that $\rho$ is not pure even at $T=0$. To see this, consider the thermal state $\rho\propto e^{-\beta H}$ for a single mode Hamiltonian $H=\epsilon a^{\dagger}a$. A simple calculation shows purity $\tr\rho^{2}=1/(2\expectationvalue{a^{\dagger}a}_{\beta}+1)$, where $\expectationvalue{a^{\dagger}a}_{\beta}$ is the expectation value of $a^{\dagger}a$ in the thermal state, and thus the state is never pure for any non-zero number of bosons. Alternatively, the tension can be traced back to the difference between the canonical ensemble and the grand canonical ensemble. As mentioned above, the Gibbs state $\rho\propto\exp{-\beta H}$ is treated within grand canonical ensemble, where the particle number is fixed only on average by a chemical potential. On the other hand, in the canonical ensemble, the Gibbs state is restricted to a fixed particle number sector $\rho\propto\exp{-\beta H}\delta(N-\sum_{i}a_{i}^{\dagger}a_{i})$. Due to the delta function constraint, all the bosons are enforced to occupy the lowest single particle state, which is a pure state that may very well be entangled.

IV Summary and Discussion

In this work we set out to reconcile the tension between the following two observations: (a) The universal long-distance correlations for typical finite-$T$ transitions are described by a low-energy effective theory that is fully classical. Sachdev (2011) (b) The mixed state quantum entanglement, as quantified by entanglement negativity, is singular across several such transitions Lu and Grover (2019). We studied a specific model that hosts a finite-$T$ transition in a mean-field universality class to understand and resolve this tension. Conceptually, our basic idea is to separate the negativity into a local term, i.e., a term which can be written as sum of local terms along the entangling boundary, and contributes to the leading area-law behavior, and a non-local term, which cannot be written in this way and therefore encodes long-distance quantum correlations. We found that in the thermodynamic limit, the singularity of negativity originates only from the local term, and can be fully canceled out by a subtraction scheme that leaves only the non-local term, which vanishes exponentially with the total system size $L$. Therefore, in the model studied, the long-distance quantum correlations are non-singular across the transition in the thermodynamic limit. We defined a length scale $\xi_{Q}$ over which quantum correlations exist, and showed that at non-zero $T$, $\xi_{Q}$ remains finite even when the physical correlation length $\xi$ diverges. Therefore quantum mechanically, the system continues to be short-range correlated, despite a diverging physical correlation length. This provides a sharp distinction between a ‘quantum phase transition’ and a ‘classical phase transition’: the quantum correlation length diverges only at a quantum phase transition.

Our discussion was focused on a Gaussian theory that exhibits mean-field critical exponents. Interactions at a finite $T$-transition in a quantum Ising model modify the critical exponents, but the critical field theory again belongs to the universality class of a classical Ising model Sachdev (2011). Therefore, our expectation is that the tripartite negativity $\Delta_{3}E_{N}$ will continue to decay exponentially with system size even at the Wilson-Fisher fixed point. Note that it is already rather non-trivial to find a quantity, namely $\Delta_{3}E_{N}$, that decays exponentially in the Gaussian critical theory. All correlation functions of local operators decay as power-law, simply because the physical correlation length is infinite. As an analogy, consider the quantum phase transition in the quantum Ising model at $T=0$. Quantum entanglement at the interacting fixed-point has the same general structure as that at the Gaussian critical theory at $T=0$, the only difference being that the value of the universal $O(1)$ subleading term in entanglement is modified Metlitski et al. (2009); Whitsitt et al. (2017); Inglis and Melko (2013); Casini and Huerta (2007). It will of course be very interesting to do an actual calculation of negativity at the finite-$T$ transition using field-theoretic techniques, to check this intuition.

We also studied the singularity associated with the area law coefficient in detail. Although all our calculations are restricted to the Hamiltonian in Eq.3, we expect these conclusions to generalize to other finite-$T$ transitions. Based on results in Ref.Lu and Grover (2019), we conjecture that the leading singular part of area-law coefficient originates from the expectation value of a local operator that is invariant under the symmetries of the Hamiltonian, and has the lowest scaling dimension. In the case of Ising model, this operator corresponds to the energy density, and therefore, we expect

[TABLE]

where $t$ is the thermal tuning parameter i.e. $t=(T-T_{c})/T_{c}$ or $t=g-g_{c}$, and $\alpha$ is the critical exponent that defines the divergence of specific heat $C\sim|t|^{-\alpha}$ (and consequently, the energy density has a singular part that scales as $|t|^{1-\alpha}$). This scaling matches with the results for our mean-field model: within mean-field theory, $\alpha=0$, and therefore $E_{N,\textrm{local, singular}}\sim|t|\sim m^{2}$.

As shown in Ref.Vidal and Werner (2002), negativity upper bounds ‘entanglement of distillation’, which intuitively corresponds to best rate at which one can extract near-perfect EPR singlets from multiple copies of a state using local operations and classical communications (LOCC). The absence of long-range negativity in our model even at the finite-$T$ critical point indicates that the distilled EPR pairs originate only from correlations close to the entangling boundary, as suggested in Fig.1. Relatedly, one expects that operators that contribute to the violation of Bell’s inequality between regions $A$ and $\bar{A}$ are located close to the entangling boundary.

Motivated by the calculation of negativity in our model that has a $\mathbb{Z}_{2}$ symmetry, we also studied negativity in closely related models that instead have an $U(1)$ symmetry. We found a surprising result that the thermal density matrices corresponding to bosonic tight-binding models with $U(1)$ and time-reversal symmetry are separable, and therefore, any mixed state measure of entanglement, including negativity, vanishes for such states. One consequence of this result is that a convex sum of density matrices $\rho=\sum_{i}p_{i}\rho_{i}$, where each $\rho_{i}$ is Gaussian and has $U(1)$ and time-reversal symmetry, is also separable. Note that $\rho$ itself is not Gaussian and corresponds to an interacting, albeit generically non-local, Hamiltonian.

Our results raise an intriguing question. It has been argued that a thermal state can be efficiently prepared if a system does not possess finite-$T$ topological order, and if it is above any finite-$T$ phase transition, so that the correlation functions of local operators are short-ranged Swingle and McGreevy (2016a, b); Brandão and Kastoryano (2019); Martyn and Swingle (2018); Wu and Hsieh (2018); Cottrell et al. (2019); Maldacena and Qi (2018); Chapman et al. (2019). Our results indicate that quantum correlations can be short-ranged even when the correlation functions of local operators are long-ranged. Although our calculations were specific to a rather simple Hamiltonian (Eq.3), we suspect that this is true more generally as long as finite-$T$ topological order is absent Hastings (2011). This raises the possibility that even density matrices that have infinite correlation length, such as quantum systems below or even at a symmetry breaking finite temperature transition, might be efficiently preparable. A starting point could be to consider a purely classical density matrix below or at the critical temperature $T_{c}$, and variationally apply a finite-depth quantum channel on it, so as to minimize the trace distance between the resulting density matrix and the actual density matrix $\rho\propto e^{-\beta H}$.

It’s worth comparing our protocol with other measures introduced previously to detect quantum coherence at finite temperature. In particular, Ref.Malpetti and Roscilde (2016) introduced a measure called ‘quantum correlation function’ (QCF) that takes the form $\langle\delta O_{A}\,\delta O_{B}\rangle_{Q}=\langle\delta O_{A}\,\delta O_{B}\rangle-\frac{1}{\beta}\int_{0}^{\beta}d\tau\langle\delta O_{A}(\tau)\,\delta O_{B}(0)\rangle$ where $\delta O=O-\langle O\rangle$, and all averages are with respect to the thermal density matrix. It was argued in Ref.Malpetti and Roscilde (2016) that this quantity is smooth across finite-$T$ transitions, and decays exponentially, allowing one to define a ‘quantum coherence length’. One advantage of QCF is that it is relatively simple to calculate compared to entanglement based measures such as negativity. On the other hand, due to its definition in terms of local operators, QCF is not suitable to capture non-local many-body entanglement. As a concrete illustration, consider a gapped, topologically ordered phase at zero temperature. Since $\Delta_{3}E_{N}$ equals the topological entanglement entropy, the ‘quantum correlation length’ introduced in our paper is infinite throughout the gapped phase. In strong contrast, the ‘quantum coherence length’ based on QCF is finite and just equals the correlation length defined via local operators.

Let us mention a few other future directions. Since our proposed method is designed to isolate the non-local part of negativity, it will be worthwhile to apply it to characterize finite-$T$ topological order in models such as 4D Toric code Dennis et al. (2002); Hastings (2011). It will also be interesting to study Renyi versions of negativity for interacting models using quantum Monte Carlo method Chung et al. (2014); Hastings et al. (2010); Isakov et al. (2012); Wu et al. (2019), and implement the subtraction scheme introduced here numerically. Another direction is to study negativity for interacting fermion systems that exhibit finite temperature phase transition such as the model of two SYK systems Sachdev and Ye (1993); Kitaev (2015) coupled to each other Maldacena and Qi (2018). We note that for fermions, two different definitions of partial transpose have been introduced Eisler and Zimborás (2015); Shapourian et al. (2017); Shapourian and Ryu (2018) and it will be interesting to understand the qualitative differences between the two in interacting theories.

Acknowledgements.

Acknowledgments: TG is supported by an Alfred P. Sloan Research Fellowship. We thank Tim Hsieh and Sung-Sik Lee for helpful discussions, and especially John McGreevy for helpful discussions and comments on the manuscript. We thank Perimeter Institute for Theoretical Physics for their hospitality where part of this manuscript was completed.

Appendix A von Neumann entropy in the Gaussian theory

Consider a $d$ dimensional system of size $L$ described by the classical Hamiltonian

[TABLE]

we prove that as $mL\to 0$, $L\to\infty$, the von Neumann entropy of the Gibbs state $\rho\sim\exp{-\beta H}$ contains a subleading term $\Delta S=-\log mL$ with the periodic boundary condition imposed on all spatial directions.

To proceed, we first discretize the continuum theory in a finite box of volume $V=L^{d}$ with lattice cutoff being $1$. A standard calculation of the Gaussian theory gives the thermal partition function

[TABLE]

where $\lambda_{k}=m^{2}+4\sum_{n=1}^{d}\sin^{2}(\frac{1}{2}k_{i})$, and $k_{i}\in 2\pi n/L$ with $n=0,\cdots,L-1$. It follows that the thermal free energy $F$ is

[TABLE]

and the von Neumann entropy is given by $S=-\partial F/\partial T$:

[TABLE]

Since the first term contributes to the volume law part of $S$, we will only focus on the second term. To proceed, we employ the identityBrankov and Danchev (1993)

[TABLE]

where $g(a)=\left\{2^{-L}\left[\sqrt{a^{2}+4}+a\right]^{L}-2^{L}\left[\sqrt{a^{2}+4}+a\right]^{-L}\right\}^{2}$. Thus one can first evaluate the product along the $d$’th spatial direction to obtain

[TABLE]

By singling out the contribution of the zero mode, one finds

[TABLE]

As $L\to\infty$ with $mL\to 0$, one finds $g(m)\to(mL)^{2}$, and thus

[TABLE]

where $\tilde{\omega}=\sqrt{\sum_{n=1}^{d-1}\sin^{2}(\frac{1}{2}k_{i})}$. By noticing that as $\tilde{\omega}=O(1)$, $L\to\infty$, $\log g(2\tilde{\omega})=2L\log\left[\left(\sqrt{\tilde{\omega}^{2}+1}+\tilde{\omega}\right)\right]$. We find

[TABLE]

Note that the above result is not exact since $\tilde{\omega}$ can be of order $1/L$, in which case $\log g(2\tilde{\omega})\sim O(1)$ instead of $O(L)$. Nevertheless it does not affect the logarithmic divergence caused by the first term. By plugging this result into Eq.17, we find the logarithmic divergence for the entropy $S\sim-\log(mL)$ for $mL\to 0$ at arbitrary dimension. In this calculation, we start from a purely classical Gaussian theory for deriving the logarithmic divergence in the massless limit. However, this result is applicable to the quantum Gaussian model studied in the main text as well by mapping a $d$ dimensional quantum problem to a classical problem with one extra dimension of size $\beta$.

Appendix B Mapping of negativity of a d-dimensional problem to a 1-dimensional problem

Consider a $d$ dimensional lattice of $N$ sites with the Hamiltonian

[TABLE]

Imposing the periodic boundary condition for all spatial direction, a standard calculation for two-point functions gives

[TABLE]

where $\mathbf{k}=\left(k_{1},k_{2},\cdots,k_{d}\right)=\frac{2\pi}{L}\left(n_{1},n_{2},\cdots,n_{d}\right)$ for $n_{i}=0,1,\cdots L-1$, and $\omega_{\mathbf{k}}=\sqrt{m^{2}+4K\sum_{i=1}^{d}\sin^{2}(\frac{1}{2}k_{i})}$. To calculate the negativity between the subregion $A$ and its complement $B$, we follow the correlation matrix technique introduced in Ref.Audenaert et al. (2002). One defines a matrix $P$ with the matrix element

[TABLE]

and the negativity is given by

[TABLE]

where $\lambda_{i}$ is the eigenvalue of $\gamma_{\phi}P\gamma_{\pi}P$.

Suppose that the boundary between $A$ and $B$ only cut through one of the spatial direction, say direction labeled by $1$, the matrix $\gamma_{\phi}P\gamma_{\pi}P$ preserves the translational invariance along the other $d-1$ directions, and therefore its eigenfunction $\psi(\mathbf{r})$ takes the form

[TABLE]

Given this ansatz, one can reduce the $d$-dimensional problem to a $1$-dimensional problem. In particular, the negativity is

[TABLE]

where $E_{N}^{(1)}(k_{2},\cdots,k_{d})$ is the negativity of the $1$-dimensional theory with a modified mass term :

[TABLE]

Note that this implies the leading contribution in $E_{N}$ is an area law term if $E_{N}^{(1)}$ in $1$-dimension follows an area law.

Appendix C Negativity between two sites

Consider a system of two sites, the dispersion relation reads $\omega_{k}=\sqrt{m^{2}+4\sin^{2}(k/2)}$ with $k=0,\pi$, and the covariance matrices are $\gamma_{\phi}=\begin{pmatrix}a_{+}&a_{-}\\ a_{-}&a_{+}\end{pmatrix}$, $\gamma_{\pi}=\begin{pmatrix}b_{+}&b_{-}\\ b_{-}&b_{+}\end{pmatrix}$ where

[TABLE]

The eigenvalues of $\gamma_{\phi}P\gamma_{\pi}P$ are

[TABLE]

While $\lambda_{1}\geq 1$ for arbitrary $\beta$ and $m$, $\lambda_{2}$ can be less than $1$ for some parameters. Thus the negativity is $E_{N}=-\frac{1}{2}\min\{0,\log\lambda_{2}\}$. In particular, while $\lambda_{1}$ diverges at the massless limit $m\to 0$, the eigenvalue $\lambda_{2}$, which contributes to the negativity, remains finite and is perturbative in $m^{2}$. Specifically, $-\log\lambda_{2}=-\log(\frac{\coth\beta}{\beta})+\left(\frac{1}{8}-\frac{\beta^{2}}{12}+\frac{\beta}{8}\coth\beta-\frac{\beta}{8}\tanh\beta\right)m^{2}+O(m^{4})$. Thus the leading order correction for the negativity in the massless limit is of order $O(m^{2})$.

Appendix D Separability of the thermal density matrix corresponding to the bosonic tight-binding model

We define the $N$-sites coherent state: $\ket{\bm{\alpha}}=\ket{\alpha_{1},\cdots,\alpha_{N}}$ such that $a_{i}\ket{\bm{\alpha}}=\alpha_{i}\ket{\bm{\alpha}}$. Note that the coherent state $\ket{\bm{\alpha}}$ can be generated by acting the displacement operator $\hat{D}(\bm{\alpha})\equiv\exp{\sum_{i=1}^{N}\left[\alpha_{i}a_{i}^{\dagger}-\alpha_{i}^{*}a_{i}\right]}$ on a vacuum state: $\ket{\bm{\alpha}}=D(\bm{\alpha})\ket{0}$.

The central ingredient we employ is the Glauber-Sudarshan P representationSudarshan (1963); Glauber (1963), which can cast any bosonic state into a diagonal matrix in the basis of coherent states:

[TABLE]

where $d^{2}\alpha_{i}\equiv d\text{Re}(\alpha_{i})d\text{Im}(\alpha_{i})$. To define $P(\bm{\alpha})$, we first introduce the characteristic function $\chi(\bm{\alpha})$Serafini (2017):

[TABLE]

$\chi(\alpha)$ serves as a generating function in the sense that its derivatives give the expectation value of the observables:

[TABLE]

Then $P(\bm{\alpha})$ is given by the complex Fourier transform of $\chi$:

[TABLE]

Up until this point, the discussion applies to any bosonic state and generically, $P(\bm{\alpha})$ is not necessarily positive. Now we specialize to Gaussian states, such as the tight-binding model of bosons (Eq.11), the subject of our focus. By defining $\alpha_{i}=\frac{1}{\sqrt{2}}\left(x_{i}+ip_{i}\right)$, one can show that for a Gaussian state, $P$ can be expressed as Serafini (2017):

[TABLE]

Here $\mathbf{r}^{T}=(x_{1},p_{1},\cdots,x_{N},p_{N})$ with $\overline{\mathbf{r}}$ being the first moment ($\overline{r}_{i}=\tr(\rho r_{i})$) and $\sigma$ being the $2N\crossproduct 2N$ covariance matrix ($\sigma_{ij}=\tr(\rho\{r_{i},r_{j}\})$). $\Omega$ denotes the symplectic matrix: $\Omega=\bigoplus_{i=1}^{N}\begin{pmatrix}0&1\\ -1&0\end{pmatrix}$. The state $\rho$ will be separable if $1-\Omega^{T}\sigma\Omega\leq 0$ since then $P(\bm{r})$ is a well-defined probability function. Below we show that a bosonic tight-binding model with $U(1)$ and time-reversal symmetry (Eq.11) indeed satisfies this condition, which proves the separability in the form of Eq.32.

D.1 Proof that $\Omega^{T}\sigma\Omega\geq 1$

We divide the proof into two-parts, we first show that $\sigma\geq 1$, i.e., all eigenvalues of $\sigma$ are greater than unity, and next we show that $\Omega^{T}\sigma\Omega$ and $\sigma$ have identical eigenspectrum, and thus $\Omega^{T}\sigma\Omega\geq 1$.

D.1.1 Proof that $\sigma\geq 1$

Given $H=-\sum_{i,j}t_{ij}a^{\dagger}_{i}a_{j}$ on a $N$-site lattice in arbitrary spatial dimension, we can define the conjugate variables: $x_{i}=\frac{1}{\sqrt{2}}\left(a_{i}+a_{i}^{\dagger}\right)$ and $p_{i}=\frac{1}{i\sqrt{2}}\left(a_{i}-a_{i}^{\dagger}\right)$. By introducing $\bm{r}^{T}=\left(x_{1},\cdots,x_{N},p_{1},\cdots,p_{N}\right)$, our goal is to calculate the covariance matrix $\sigma_{ij}=\expectationvalue{\{r_{i}-\overline{r}_{i},r_{j}-\overline{r}_{j}\}}$, where $\{A,B\}=AB+BA$, and $\overline{r_{i}}=\expectationvalue{r_{i}}$ with the expectation values taken with respect to a thermal state $\rho=e^{-\beta H}/Z$. To proceed, we write $a_{i}=\sum_{k}u_{ik}b_{k}$ to diagonalize the Hamiltonian: $H=\sum_{k}E_{k}b_{k}^{\dagger}b_{k}$. Given this result, we find

[TABLE]

and $\overline{r_{i}}=0$. Note that $u$ can be chosen as an orthogonal matrix since $t_{ij}$ is real (due to time reversal symmetry of $H$). Therefore, $\sigma=\gamma_{x}\bigoplus\gamma_{p}$, with

[TABLE]

or equivalently, $\gamma_{x}=u\Lambda u^{\dagger}$, where \Lambda=\text{ diag }\left\{1+2\expectationvalue{b_{k}^{\dagger}b_{k}}\big{|}\text{for }k=1,2,\cdots,N\right\} storing the eigenvalues of $\gamma_{x}$. This completes the proof that $\sigma=\gamma_{x}\bigoplus\gamma_{p}\geq 1$.

D.1.2 Proof that $\sigma$ and $\Omega^{T}\sigma\Omega$ have identical spectrum

To prove this, first notice that an orthogonal matrix $O$ can diagonalize $\sigma$ since it is symmetric: $\sigma=u\Lambda u^{T}$, where $\Lambda$ is a diagonal matrix. Then,

[TABLE]

Recall that $\Omega^{T}\Omega=\Omega\Omega^{T}=-1$, we then have

[TABLE]

where $O^{\prime}\equiv\Omega^{T}O\Omega$ and $\Lambda^{\prime}=\Omega^{T}\Lambda\Omega$. A straightforward calculation shows that $O^{\prime}$ is an orthogonal matrix, and $\Lambda^{\prime}$ is a diagonal matrix. In particular, the set of the diagonal entries of $\Lambda^{\prime}$ is exactly the same as that of $\Lambda$, and hence $\Omega^{T}\sigma\Omega$ and $\sigma$ have the same eigenspectrum. Explicitly, if

[TABLE]

then

[TABLE]

D.2 An aside: alternative proof for zero negativity for bosonic tight-binding model

Recall that for a separable state, all measures of entanglement, including negativity, are identically zero, but a state with zero negativity can still be non-separable Horodecki et al. (1998, 2009). Therefore, the above proof already establishes that negativity corresponding to a tight-binding Hamiltonian of bosons with $U(1)$ and time-reversal is identically zero. Still, as an aside, we provide an alternate proof of zero negativity for this system.

From the fact that the eigenvalues of the correlation matrix $\gamma_{x}$ are greater than or equal to $1$, i.e. $\gamma_{x}\geq 1$, we prove that $\gamma_{x}P\geq 1$, and thus all eigenvalues of $\gamma_{x}P\gamma_{x}P$ are greater than or equal to 1. It follows that the negativity is zero from Eq.26.

Proof for $\gamma_{x}P\geq 1$:

Given a $N\crossproduct N$ symmetric matrix $\gamma_{x}$ with all eigenvalues $\lambda_{m}\geq 1$ and a diagonal matrix $P$ with $1$ or $-1$ diagonal entries. First we perform an eigendecomposition of $\gamma_{x}$:

[TABLE]

where $O$ is an orthogonal matrix and $\Lambda=\text{diag}\left(\lambda_{1},\cdots,\lambda_{N}\right)$. Next the eigenequation of $\gamma_{x}P$ reads

[TABLE]

where $\tilde{\lambda}_{m}$ and $v_{m}$ are the eigenvalue and eigenvector respectively. Taking the transpose of the above equation gives

[TABLE]

From the product of the above two equations and employing Eq.43, we obtain

[TABLE]

where $\tilde{v}_{m}=OPv_{m}$. Since $\tilde{v}_{m}^{T}\tilde{v}_{m}=\sum_{i=1}^{N}\tilde{v}^{2}_{m,i}=1$, $\tilde{\lambda}^{2}_{m}\geq\sum_{i=1}^{N}\tilde{v}^{2}_{m,i}\lambda_{\text{min}}^{2}\geq 1$, where $\lambda_{\text{min}}$ is the minimum eigenvalue of $\gamma_{x}$. Thus, all eigenvalues $\tilde{\lambda}_{m}$ of $\gamma_{x}P$ satisfy $\tilde{\lambda}_{m}\geq 1$.