Environment mediated multipartite and multidimensional entanglement

TL;DR

This paper explores how an engineered environment can facilitate the formation of multipartite and multidimensional entanglement in quantum systems, challenging the notion that the environment always destroys entanglement.

Contribution

It demonstrates that a reservoir can mediate entanglement among qubits and qudits at low temperatures and intermediate coupling, with a new witness for multipartite entanglement.

Findings

Multipartite entanglement can form via environment mediation.

Entanglement witness based on structure factor is effective.

Possible experimental measurement of entanglement proposed.

Abstract

Quantum entanglement is usually considered a fragile quantity and decoherence through coupling to an external environment, such as a thermal reservoir, can quickly destroy the entanglement resource. This doesn't have to be the case and the environment can be engineered to assist in the formation of entanglement. We investigate a system of qubits and higher dimensional spins interacting only through their mutual coupling to a reservoir. We explore the entanglement of multipartite and multidimensional system as mediated by the bath and show that at low temperatures and intermediate coupling strengths multipartite entanglement may form between qubits and between higher spins, i.e., qudits. We characterise the multipartite entanglement using an entanglement witness based upon the structure factor and demonstrate its validity versus the directly calculated entanglement of formation, suggesting possible experiments for its measure.

Equations

$$H_{\mathrm{S}}$$

$$H_{\mathrm{B}}$$

$$H_{\mathrm{I}}$$

$$P=2\sum _j{s}_z^{(j)}\sum _n\frac{g_n}{{\omega }_n}\left(b_n^{†}-b_n^{}\right).$$

$${\tilde{H}}_{\mathrm{S}}$$

$${\tilde{H}}_{\mathrm{I}}$$

$${\tilde{\rho }}_{\mathrm{S}}=\frac{{\mathrm{Tr}}_{\mathrm{B}}\left(e^{-\beta {\tilde{H}}_{\mathrm{tot}}}\right)}{{\mathrm{Tr}}_{\mathrm{SB}}\left(e^{-\beta {\tilde{H}}_{\mathrm{tot}}}\right)}.$$

$${\tilde{\rho }}_{\mathrm{S},0}$$

$${\tilde{\rho }}_{\mathrm{S},2}$$

$$A=\sum _{\ell }{\int }_0^{\beta }d\tau {\int }_0^{\tau }d{\tau }^{\prime }{C}_{\ell \ell }(\tau -{\tau }^{\prime })e^{-\beta {\tilde{H}}_{\mathrm{S}}}{S}_{\ell }(\tau )S_{\ell }({\tau }^{\prime }),$$

$${\mathfrak{S}}_{\ell {\ell }^{\prime }}(k)=\sum _{i<j}{e}^{\mathrm{i}k(r_j-r_i)}⟨s_{\ell }^{(i)}{s}_{{\ell }^{\prime }}^{(j)}⟩,$$

$$\displaystyle\bar{\Sigma}(k)=\frac{1}{2{N\choose 2}}\Bigl{[}c_{x}{\mathfrak{S}}_{xx}(k)+c_{y}{\mathfrak{S}}_{yy}(k)+{c_{z}}{\mathfrak{S}}_{zz}(k)\Bigr{]}.$$

$$\mathrm{EoF}\ge {\mathrm{EoF}}_{\mathrm{lb}}=-{\log }_2\left(1-\frac{1}{2}{R}^2\right),$$

$$R=\frac{2}{\sqrt{\lvert C\rvert}}\sum_{\begin{subarray}{c}(j,k)\in C\\ j<k\end{subarray}}{\!\!\Bigl{[}\lvert\left\langle j,j\middle|\rho\middle|k,k\right\rangle\rvert-\!\sqrt{\langle j,k|\rho|j,k\rangle\langle k,j|\rho|k,j\rangle}\Bigr{]}},$$

$$H_{\mathrm{S}}$$

$$H_{\mathrm{B}}$$

$$H_{\mathrm{SB}}$$

$$P=\sum _j{\sigma }_z^{(j)}\sum _n\frac{g_n}{{\omega }_n}(b_n^{†}-b_n^{}).$$

$${\tilde{H}}_{\mathrm{S}}$$

$${\tilde{H}}_{\mathrm{B}}$$

$${\tilde{H}}_{\mathrm{SB}}$$

$${\Delta }_{\mathrm{R}}=\Delta \exp \left(-\frac{2}{\pi }{\int }_0^{\infty }\frac{J(\omega )}{{\omega }^2}\coth \left(\frac{\beta \omega }{2}\right)d\omega \right).$$

$$J(\omega )=\gamma {\left(\frac{\omega }{{\omega }_{\mathrm{c}}}\right)}^k{e}^{-\omega /{\omega }_{\mathrm{c}}}$$

$$\ln{\!\left(\frac{\Delta_{\mathrm{R}}}{\Delta}\right)}=-\frac{2}{\pi}\left(\frac{\gamma}{\omega_{c}}\right)\Gamma(k-1)\Bigl{[}2t^{k-1}\zeta(k-1,t)-1\Bigr{]}$$

$$\frac{\Delta_{\mathrm{R}}}{\Delta}\cong\exp{\!\left(-\frac{2}{\pi}\left(\frac{\gamma}{\omega_{\mathrm{c}}}\right)\Bigl{[}1+\frac{\pi^{2}}{3}t^{2}\Bigr{]}\right)}\rightarrow\exp{\!\left(-\frac{2}{\pi}\frac{\gamma}{\omega_{\mathrm{c}}}\right)}.$$

$$\frac{{\Delta }_{\mathrm{R}}}{\Delta }\cong \exp \left(-\frac{4}{\pi }\frac{\gamma }{{\omega }_{\mathrm{c}}}t\right)\to 0.$$

$$V_x$$

$$V_y$$

$$C_{xx}(\tau )$$

$$C_{yy}(\tau )$$

$$\phi(\tau)=\frac{4}{\pi}\int_{0}^{\infty}{\frac{J(\omega)}{\omega^{2}}\frac{\cosh{\bigl{(}\frac{1}{2}\,(\beta-2\tau)\,\omega\bigr{)}}}{\sinh{\bigl{(}\tfrac{1}{2}\beta\omega\bigr{)}}}\,d\omega}.$$

$$J_{\text{Ohmic}}(\omega )$$

$$J_{\text{Lorentz}}(\omega )$$

Full text

Entanglement between systems with few degrees of freedom is often a fragile property as their interaction with a much larger environment drives decoherence. This decoherence often comes from the small system interacting in an incoherent, memory-less manner with the bath. Most theoretical studies of few-body quantum systems make this assumption explicitly or implicitly, i.e., the dynamics is Markovian. Such approximation seems quite appropriate and experimentally valid to high precision in condensed matter systems Gardiner and Zoller (2010). However, it is also a possibility that the environment, if suitably prepared, can mediate the creation of quantum correlation between small systems Braun (2002).

Quantum entanglement for two or even three qubits has been studied extensively Amico et al. (2008); Horodecki et al. (2009) and it is well understood, but less well characterised is multipartite entanglement – entanglement between several systems. For example, for more than four qubits, the creation as well as the characterisation of quantum entanglement becomes exceedingly difficult. Theoretical investigations also face the curse of dimensionality, since the dimension of the corresponding Hilbert-space grows exponentially with the number of subsystems. Despite the hurdles, experimental demonstration of multipartite entanglement has been achieved in ion traps with up to 14 ions Häffner et al. (2005); Leibfried et al. (2005); Monz et al. (2011), or for ten superconducting qubits Vidal et al. (2003); Song et al. (2017) and ten photonic qubits on a linear optical platform Wang et al. (2016).

The main obstacle to achieving significant multipartite entanglement is decoherence and dissipation due to coupling to the environment. The lifetime of entangled states decreases rapidly with the size of the system. Standard approaches to overcome this difficulty involve high fidelity quantum gates and large qubit overheads, which can be more and more challenging as the number of qubits increases. On the other hand, steady state two-qubit entanglement can be generated by engineered dissipation in trapped ions Lin et al. (2013) and superconducting qubits Shankar et al. (2013), and noise assisted quantum transport has also been reported Viciani et al. (2015); Biggerstaff et al. (2016). These experimental results challenge the general assumption that dissipation is always detrimental to quantum information processing, rather they demonstrate that it might be a resource that can be engineered and harnessed.

It has been conjectured that non-Markovian effects in the quantum time evolution, assuming an Ohmic bosonic bath at zero temperature, can leave their fingerprints on the time dependence of coherence Divincenzo and Loss (2005). This conjecture relies on the applicability and accuracy of the Born approximation underpinning the Bloch-Redfield equation. It should be noted that this perturbation approach contains accumulating divergences Slutskin et al. (2011), rendering the Born approximation inappropriate at long times. As a consequence care should be taken when interpreting these theoretical investigations, but the different Markovian and non-Markovian processes which affect relaxation and dephasing do offer opportunities to engineer the effects of the environment upon the system of interest.

Another research direction proposes entangling systems with more numerous degrees of freedom  Mair et al. (2001). In general we may replace a standard two-level system with a $d$-dimensional system – a qudit. Such generalisation is not a simple extension of already existing concepts: the entanglement of two-qubit systems seems somewhat the exception Vollbrecht and Werner (2000) rather than the rule. In addition, there is qualitative difference between bipartite and multipartite entanglement Bengtsson and Zyczkowski (2016) and the latter may provide higher flexibility to achieve more efficient quantum information processing Gröblacher et al. (2006). The main focus of recent experiments has been to generate and characterise multilevel and multipartite entanglement Mair et al. (2001); Martin et al. (2017); Krenn et al. (2013); Molina-Terriza et al. (2007). Our work presented here aims to contribute to these efforts.

Here we build a theoretical model to generate steady state multipartite and multilevel entanglement mediated solely by the environment. First we consider a general open quantum system coupled to a bath and then we examine the equilibrium reduced density matrix of non-interacting qudits through the analytical polaron treatment, complemented by a numerically exact path-integral approach. Next, we demonstrate the use of the structure-factor-based and global entanglement witnesses to compute entanglement in bipartite and multilevel systems. For a non-interacting pair of qubits, we demonstrate that the entanglement is maximised at finite amount of system-bath interaction, and the peak values are temperature dependent.

Our results can be tested experimentally, as scattering experiments have already proved to be useful in quantifying multipartite entanglement Marty et al. (2014). Structure factors are macroscopic quantities measurable in various periodic systems, e.g., Bragg scattering on ultracold atoms in optical lattices Corcovilos et al. (2010) or from ion chains Macchiavello and Morigi (2013) and neutron scattering from solid state systems. Furthermore, the structure factor, particularly the dynamic structure factor, is a central observable quantity in determining single-particle and collective electronic excitations in many-body systems Pitarke et al. (2007); Pines and Nozieres (1966).

Model — We consider the collection of $N$ identical, non-interacting qudits. The energy splittings between levels are assumed to be uniform and denoted by $\epsilon$ while the tunnelling matrix element is $\Delta$. The bath is modelled as an ensemble of harmonic oscillators with frequencies $\omega_{n}$. The strength of coupling between a qudit and an oscillator with angular frequency $\omega_{n}$ is denoted by $g_{n}$. Thus the total system-plus-bath Hamiltonian is formally written as $H_{\mathrm{tot}}=H_{\mathrm{S}}+H_{\mathrm{B}}+H_{{\mathrm{I}}}$, where $H_{\mathrm{S}}$, and $H_{\mathrm{B}}$ govern the internal dynamics of the system and the bath, respectively, while $H_{\mathrm{I}}$ describes the coupling between them. In detail

[TABLE]

where $s_{x}$ and $s_{z}$ are the $x$ and $z$ components of a spin operator of a system with spin $s$, and the $b_{n}$ are the bosonic bath operators. The number of degrees of freedom of the bath is assumed to be much larger than that of the qudits, thus according to the Poincaré recurrence theorem the time-scale on which the bath feeds back the energy to the subsystems is enormous.

This separation of time-scales allows us to smooth the oscillator modes and instead of $\{\omega_{n}\}$ we characterise the bath with a continuous spectral density $J(\omega)$. For analytical calculations below, we assume a super-Ohmic spectrum, $J(\omega)=\gamma(\frac{\omega}{\omega_{\mathrm{c}}})^{3}e^{-\omega/\omega_{\mathrm{c}}}$, of exponential cut-off with a fixed $\omega_{\mathrm{c}}$, and constant system-bath coupling strength, $\gamma$, with the dimension of frequency. The reciprocal of the cut-off frequency, $\tau\propto\tfrac{1}{\omega_{\mathrm{c}}}$ governs the relaxation time of the bath. This form of the spectral density is commonly used in the study of tunnelling effects in solid state systems Leggett et al. (1987). Though the analysis below utilises a super-Ohmic spectrum, numerical results with other spectral densities, e.g., Ohmic and Lorentzian (see supplementary information), show that the conclusions are applicable to a wide range of systems. In order to capture the effects of strong qudit-bath coupling, we employ the polaron transformation defined as $\widetilde{H}_{\mathrm{tot}}=e^{P}H_{\mathrm{tot}}e^{-P}$, where

[TABLE]

The total Hamiltonian preserves its structure under this transformation: $\widetilde{H}_{\mathrm{tot}}=\widetilde{H}_{\mathrm{S}}+\widetilde{H}_{\mathrm{B}}+\widetilde{H}_{\mathrm{I}}$, where $\widetilde{H}_{\mathrm{B}}=H_{\mathrm{B}}$ and

[TABLE]

Several remarks are in order. First, the Hamiltonians are partitioned such that the expectation value $\operatorname{Tr}_{\mathrm{SB}}{(\widetilde{\rho}\widetilde{H}_{\mathrm{I}})}=\langle\widetilde{H}_{\mathrm{I}}\rangle_{\mathrm{SB}}=0$. Second, an effective qudit-qudit interaction term, $\sim s_{z}^{(j)}s_{z}^{(k)}$, appears in $\widetilde{H}_{S}$, induced by the system-bath coupling. This term is the source of the finite equilibrium entanglement in the system since the qudits were originally assumed to be non-interacting. Third, the system-bath coupling in the polaron picture, $\widetilde{H}_{\mathrm{I}}$, assumes a form very different from $H_{\mathrm{I}}$. All spin operators appear in a symmetric way and they are coupled to effective bath operators $V_{x}$, $V_{y}$. The explicit expressions of these bath operators are immaterial for our current analysis but for sake of transparency they are provided in the supplementary material. Finally, the polaron transformation renormalises the tunnelling rate, $\Delta\mapsto\Delta_{\mathrm{R}}$, and it thus picks up a temperature and coupling strength dependence.

The equilibrium reduced density matrix can formally be expressed Weiss (2012) as a partial trace

[TABLE]

Expanding $\tilde{\rho}_{S}$ and taking into account that $\langle\widetilde{H}_{\mathrm{I}}\rangle_{\mathrm{SB}}=0$ guarantees that in perturbative expansion the leading order term for $\widetilde{\rho}_{\mathrm{S}}$ is of second order in the interaction Hamiltonian. The equilibrium state of the system can be approximated Lee et al. (2012) as $\widetilde{\rho}_{\mathrm{S}}\approx\widetilde{\rho}_{\mathrm{S},0}+\widetilde{\rho}_{\mathrm{S},2}$ where

[TABLE]

The leading $\widetilde{\rho}_{\mathrm{S},0}$ term is the expected canonical equilibrium distribution in the polaron picture with $Z_{\mathrm{S},0}=\operatorname{Tr}_{\mathrm{S}}{\!(e^{-\beta\widetilde{H}_{\mathrm{S}}})}$, while the correction term, $\widetilde{\rho}_{\mathrm{S},2}$, picks up some of the correlation between the multi-level subsystems and the bath via

[TABLE]

where $S_{\ell}=\sum_{j}s^{(j)}_{\ell}$ and $Z_{\mathrm{S},2}=\operatorname{Tr}_{\mathrm{S}}{\!(A)}$. Even the zeroth order term captures some correlation through the transformed system Hamiltonian. The exact expressions of $C_{\ell\ell}(\tau)=\langle V_{\ell}(\tau)V_{\ell}(\tau^{\prime})\rangle_{\widetilde{H}_{\mathrm{B}}}$ are in the supplementary information.

It is useful at this point to analyse the behaviour of perturbation theory at strong coupling in the polaron frame. As $\gamma\mapsto\infty$ the system becomes incoherent since the coherent tunnelling element, $\Delta_{\mathrm{R}}$ vanishes. Simultaneously the correlation functions, $C_{\ell\ell}(\mathrm{\tau})$, and hence the second-order correction to the reduced density matrix $\widetilde{\rho}_{\mathrm{S},2}$ vanishes. Therefore the equilibrium density matrix at strong $\gamma$ is dominated by the bath induced coupling term, $\sigma_{z}^{(j)}\sigma_{z}^{(k)}$, e.g., for qubits $\widetilde{\rho}_{\mathrm{S}}\propto\exp(-\beta\widetilde{H}_{\mathrm{S}})$, where $\widetilde{H}_{\mathrm{S}}$ is a diagonal matrix, since it is proportional to $\sigma_{z}$. Therefore $\widetilde{\rho}_{\mathrm{S}}$ is also diagonal. This result, namely that $\widetilde{\rho}_{\mathrm{S}}$ follows the canonical Gibbs’ distribution, at least in the strong coupling regime, is in agreement with Johnson et al. (2014). In the strong coupling (and limit of weak coupling) the density matrix can be factored as a product of single particle density matrices so this qubit system can be considered as a classical simulator in the language of Johnson et al. (2014). It is only in the intermediate coupling regime, in the presence of entanglement, that the system provides a quantum resource.

Entanglement witness — The Hahn–Banach separation theorem guarantees that for any entangled state there is an entanglement witness 111In the physics literature the concept of entanglement witness is predominantly defined as a functional, $W$, of the density matrix, $\rho$, such that $\operatorname{Tr}{\!(W\!\rho)}<0$ if $\rho$ is entangled, and $\operatorname{Tr}{\!(W\!\rho)}>0$ if $\rho$ is separable.. However, currently there is no universal measure of entanglement for an arbitrary number of subsystems each possessing multiple levels, rather there are a plethora of functionals which are useful in specific situations. Unlike qudit and multipartite systems, bipartite qubit systems can be classified by a single measure Vedral et al. (1997). This fact shows that bipartite qubit systems are not typical, rather the exceptions.

In order to quantify bipartite qubit entanglement, we adopt the standard measure of entanglement of formation (EoF) Wootters (1998). However, detecting multipartite entanglement for a mixed state is more complicated Krammer et al. (2009), and we rely on the structure factor as a precursor of an entanglement witness

[TABLE]

where $\ell$, $\ell^{\prime}$ run over the set $\{x,y,z\}$, $k$ is the wavenumber, and $r_{i}$, $r_{j}$ are the positions of spins $i$ and $j$. The distance between two adjacent qubits is normalised to unity 222The physical distance of two qubits is arbitrary. If the qubits are at equidistant locations with distance $a$ then a structure factor ${\mathfrak{S}}_{\ell\ell^{\prime}}(k,1)$ calculated with unit distance and another structure factor ${\mathfrak{S}}_{\ell\ell^{\prime}}(k,a)$ calculated with a different distance, $a$, are related to each other as ${\mathfrak{S}}_{\ell\ell^{\prime}}(k,1)={\mathfrak{S}}_{\ell\ell^{\prime}}(\frac{k}{a},a)$. The entanglement witness, based on ${\mathfrak{S}}_{\ell\ell^{\prime}}(k)$, is defined as $W(k)=\mathbb{1}-\Sigma(k)$, where $\Sigma(k)=\bar{\Sigma}(k)+\bar{\Sigma}(-k)$, and

[TABLE]

It has been proven Krammer et al. (2009) that $W$ detects multipartite entanglement of a mixed state whenever $\left\langle W\right\rangle=\operatorname{Tr}{\!\left(W\!\rho\right)}<0$, or equivalently $\langle\Sigma\rangle>1$. In general the coefficients $c_{\ell}$ can be chosen arbitrarily by the observer provided all $c_{\ell}$ are real and $\lvert c_{\ell}\rvert\leq 1$. Below we use $(c_{x},c_{y},c_{z})=(1,-1,1)$. In the case of qudit systems we use a convenient lower bound Tiranov et al. (2017) of the entanglement of formation

[TABLE]

where

[TABLE]

with $C$ being the set containing all pairs of indices $(j,k)$, while $\lvert C\rvert$ denotes its cardinality.

Results — We first revisit the two-qubit case, i.e., $s_{\ell}\mapsto\tfrac{1}{2}\sigma_{\ell}$, and while this system has been analysed before Braun (2002); McCutcheon et al. (2009) we expand those analyses here. We use an imaginary time path-integral technique for calculating the reduced density matrix for the bare system, while for the dressed particles the previously described polaron transformation and second order approximation is employed since it yields accurate equilibrium density matrices for a wide range of parameters Lee et al. (2012), except for simultaneously low $\omega_{c}$ and temperatures.

In the weak coupling regime, increasing the bath mediated interaction between qubits enhances the induced entanglement. On the other hand, entanglement decreases exponentially as a function of $\gamma$ in the strong coupling regime. The interaction term, $\sim\sigma_{z}^{(j)}\sigma_{z}^{(k)}$ in $\tilde{H}_{\mathrm{S}}$, dominates in the strong coupling limits, and the resulting $\widetilde{\rho}_{\mathrm{S}}$ is localized (i.e., diagonal) in the $\sigma_{z}$ basis with zero entanglement. The scaling limits the entanglement in both weak and strong coupling regimes and gives rise to the maximal entanglement at intermediate $\gamma$. The non-monotonic behavior of EoF as a function of $\gamma$ may be reminiscent of the “stochastic resonance” phenomenon reported by Huelga Huelga and Plenio (2007) and Lee Lee et al. (2011). Unlike the coupled quantum many-body systems considered in these earlier works, the qubits here are not coupled directly, the entanglement is solely induced by the interaction with the common reservoir.

It should be noted that $\widetilde{\rho}_{\mathrm{S}}$ is quite different from $\rho_{\mathrm{S}}$. The transform describes a system “dressed” by the polaron cloud due to the system-bath interaction. However, this effective system is similar to those typically found in solid state quantum dots and BCS pairs. It is interesting to compare the amount of entanglement between two dressed qubits with that between two bare qubits. While the analytical polaron method is capable of exploring the entire range of dissipation strength and temperature, we need to employ the numerically exact imaginary time path-integral technique Moix et al. (2012) to obtain $\rho_{S}$. Given the high computational cost of the path-integral technique, we restricted our analysis to small systems, at moderate $\gamma$ and temperature values.

The EoF between the qubits as a function of dissipation strength, $\gamma$, is shown in Fig. 2. The dependence of entanglement on $\gamma$ between two bare qubits is similar to that of polaron qubits, albeit at a lower value. The difference of EoF between dressed qubits and bare qubits is a consequence of the system-phonon and phonon-phonon entanglement, captured through the polaron transformation but not in the calculation of entanglement between the bare qubits. The polaron transformation dresses the qubit with phonons in the bath and correlations in this phonon gas then contribute to the polaron entanglement. This result does not contradict our assumption of Markovian (memory less) time-evolution, rather it reflects the fact that within the environment equilibrium correlation has built up.

Fig. 2 also illustrates the structure-factor-based entanglement witness, $\langle W\rangle$, for the two-qubit system. The negativity of $\langle W\rangle$ is proportional to the entanglement of formation. Of course $\langle\Sigma\rangle$ could also be used as an approximate measure of the entanglement. The advantage of both $\Sigma$ and $W$, is that they are both directly related to the structure factor, ${\mathfrak{S}}_{\ell\ell^{\prime}}(k)$ and one can thus connect to the wealth of knowledge in designing scattering experiments in order to interrogate the coupled system. Furthermore, adjusting the constants $c_{\ell}$’s one may explore different “pockets” of the high-dimensional Hilbert space Vedral (2014) and can further classify the entangled state.

The inset of Fig. 2 shows the EoF between two polaron qubits over a large range of $\gamma$ and inverse temperature, $\beta$. It is observed that the bath mediated equilibrium entanglement is also sensitive to temperature, in addition to $\gamma$. At high temperature, all density matrices approach a totally mixed state since all eigenstates are equally populated, thus $\widetilde{\rho}_{\mathrm{S}}\rightarrow\mathbb{1}$ and $\widetilde{\rho}_{\mathrm{S},0}\rightarrow\mathbb{1}$ as $A\rightarrow 0$. As a result, entanglement is destroyed at high temperatures, as expected. However, at lower temperatures the entanglement of dressed qubits develops and may become significant ($\sim 0.8$).

Following the illustrative two-qubit case, we next explore the possibility of observing bath induced multipartite entanglement. Let us consider a linear chain of up to six noninteracting qubits ($N=6$) coupled to a common harmonic bath. Fig. 3 demonstrates the behaviour of $\langle W\rangle$ for different system sizes, including $N=2$–6 qubits (red bold labels). As seen, increasing the system size by one additional qubit causes a significant decrease in optimal peak of $\langle W\rangle$ from $-0.926$ to around $-0.4$ with the same noise level. Increasing the system size further, we again observe a lower optimum peak for $\langle W\rangle$. Measuring $\langle W\rangle$ for systems with $5$ and $6$ qubits shows that system size has a direct effect on the correlation and the maximal entanglement moves to lower ranges in the coupling for larger systems.

The inset of the top panel of Fig. 4 displays the entanglement witness, $\langle W\rangle$, as a function of $\gamma$ at different temperatures as indicated by the red, bold numbers which are the dimensionless values of $\beta$. Similar to the bipartite case, it is shown that multipartite entanglement, i.e., $\langle W\rangle<0$, can be generated at intermediate dissipation strength, $\gamma\approx 5\times 10^{-3}$, provided the temperature is low enough. As temperature increases, thermal fluctuations destroy the equilibrium multipartite entanglement as the system reduced density matrix approaches a totally mixed state. Turning to the main graph in the top panel of Fig. 4, a negative $\langle W\rangle$ can be seen at low wavenumbers for low temperatures. This implies the existence of multipartite entanglement occurring around $k\cong 0$ for $\beta>3$.

Finally, the lower heatmap of Fig. 4 displays $\langle W\rangle$ as a function of both $\beta$ and $\gamma$. Dashed black line separates two distinct regimes of multipartite entanglement. It seems worthwhile mentioning the generic u-shape of the solid lines, indicating that multipartite entanglement develops only within a range of $\gamma$ values and below a critical temperature value. While the latter feature, the existence of a critical temperature is expected, the strong dependence on $\gamma$ is –we believe– novel. The steep slope of the solid curve in terms of $\gamma$ suggests that weak coupling to the bath cannot efficiently relay entanglement and the qubits evolve in time more or less independently of each other. On the other hand, strong coupling with the bath connects the qubit with numerous phonon modes, leading to decoherence and the inhibition of any transmission of entanglement between the qubits.

From here onwards we focus on multilevel systems. One may extend the results above to higher spin values by using appropriate spin operators as the formalism has already been presented in this general form. For example, the form of the Hamiltonian for a spin-1 system (qutrits) is similar to that of the spin-$\tfrac{1}{2}$ system, with the main difference being the occurrence of the anisotropic term, $\sim s_{z}^{2}$ in $\widetilde{H}_{\mathrm{S}}$. This term contributes only a constant in the case of $s_{z}=\frac{1}{2}\sigma_{z}$. The polaron transformation can still be defined in this higher dimensional space. In order to apply the transform to $H_{\mathrm{tot}}$, the Baker-Campbell-Hausdorff formula Hall (2003) is needed, but this does not pose any serious limitation on the calculation.

The bottom panel of Fig. 3 compares $\mathrm{EoF}_{\mathrm{lb}}$ for bipartite qudit systems from $s=\tfrac{1}{2}$ up to $s=3$ in steps of $\tfrac{1}{2}$ (red bold labels). Notice the log-log scale of the axes. Increasing the dimensionality of a subsystem manifests itself in more brittle entanglement such that at the same noise-level in the bath lower entanglement can be achieved. Therefore the maximal entanglement shift towards weaker interactions by increasing spin dimensions.

Now looking at Fig. 5, we see that the entanglement of formation exhibits a stochastic resonance behaviour in qutrit systems similar to that in the qubit system, but to a more limited extent.

Moreover, by comparing the curves of qubits and qutrits in the main panel of Fig. 5, we may conclude that observable entanglement in qutrits system can occur only at lower temperatures than in qubit systems. This is expected as the temperature is a proxy for the noise in the phonon-bath. Thus increasing temperature rises noise-levels which destroy correlation. The stochastic resonance thus happens at lower noise levels for the qutrit systems. In summary, the maximal entanglement in both cases depends on temperature and happens at moderate system-bath coupling strengths.

Concluding, we have demonstrated that steady state multipartite and multilevel entanglement can be induced solely through interactions with a common reservoir. The multipartite and multilevel entanglement are experimentally accessible and we propose that an analysis of the system’s structure factors could be used for their detection. Structure factors, both static and dynamic, are commonly used in solid-state physics, and are suited to the capture of correlational structure in few-body systems. They are measurable in scattering experiments, e.g., via neutron scattering in condensed-matter systems, or via light scattering off optical lattices.

We reiterate that the analysis presented here focuses on a bath with a super-Ohmic spectrum for analytical treatment of the polaron perturbation technique. Numerical results with other spectral densities (Ohmic and Lorentzian) show that our conclusions are essentially independent of the precise form of spectral density and applicable to a wide range of systems.

It is clear that the spin-boson models are an important tool in understanding imperfections in quantum gate operation. The results presented above make some first steps in characterising quantitatively steady-state multipartite and multi-dimensional entanglement, thereby helping provide an additional resource for quantum information applications. We envision that entanglement measures based on the structure factor could become a good tool to probe and quantify multipartite qudit entanglement and simultaneously test setups with controllable dissipation, e.g., trapped ions and superconducting qubits.

Fermions (qudits) in optical lattices seem particularly suitable candidates for investigation as these systems are highly controllable, e.g., the potential depth and external magnetic field can control $\epsilon$ and $\Delta$, and how strongly these atoms interact with each others. If a thermal cloud of atoms (the bath) were trapped in a different potential and brought in contact with the strongly trapped atoms in the lattice, their interaction, $\gamma$, could also be tuned. Then one may potentially interrogate the trapped atoms via scattering Marty et al. (2014); Landig et al. (2015), and deduce the structure factor to quantify entanglement.

Acknowledgements — We thank our funding agencies, particularly the NZ Tertiary Education Commission, the National Research Foundation Singapore and the Ministry of Education Singapore. One of us (DAWH) would like to thank Merton College, Oxford, where he is a Visiting Research Fellow, for their generous hospitality.

Appendix A: Details of multipartite qubit system We now discuss in a little more detail the qubit system. The total Hamiltonian for $N$ qubits, $H_{\mathrm{tot}}=H_{\mathrm{S}}+H_{\mathrm{B}}+H_{\mathrm{SB}}$, in the open quantum system formalism is given by

[TABLE]

where $H_{\mathrm{S}}$, $H_{\mathrm{B}}$ and $H_{\mathrm{SB}}$ are the system, bath and the system-bath coupling Hamiltonians, respectively. The system is described by a group of $N$ non-interacting identical two-level systems with energy splitting $\epsilon$ and tunnelling matrix element $\Delta$. The corresponding polaron transformation is

[TABLE]

Up to a constant we obtain $\widetilde{H}_{\mathrm{tot}}=\widetilde{H}_{\mathrm{S}}+\widetilde{H}_{\mathrm{B}}+\widetilde{H}_{\mathrm{SB}}$ with

[TABLE]

For $\widetilde{H}_{\mathrm{S}}$, the tunnelling element in $\widetilde{H}_{\mathrm{S}}$ is renormalized due to the system bath coupling,

[TABLE]

In order to get some insight into the behaviour of the renormalised tunnelling rate, $\Delta_{\mathrm{R}}$, one may evaluate this integral for a generic super-Ohmic spectral density

[TABLE]

where $k$ has to be bigger that 2 to mollify the singularity at $\omega=0$. Substituting this expression into the integral one may arrive at

[TABLE]

where $\zeta(s,a)$ is the generalised Riemann-zeta function and $t=k_{\mathrm{B}}T/(\hbar\omega_{\mathrm{c}})$ is a dimensionless temperature. Let us choose $k=3$ in order to be consistent with the spectral density used in the main body of the work and examine the low- and high-temperature asymptotic of this expression. As $t\rightarrow 0^{+}$ the generalised Riemann-zeta function develops a singularity, however, its pre-factor $t^{k-1}$ cancels that precisely and we obtain

[TABLE]

Thus the renormalised tunnelling term approaches a constant which depends on the coupling strength and cut-off frequency. This result demonstrates that with increasing system-bath coupling the tunnelling term vanishes exponentially. In the opposite, high-temperature limit, the zeta function vanishes as $(k-2)t^{2-k}+2t^{1-k}$, thus in our case ($\gamma>0$)

[TABLE]

As anticipated, for decoupled systems, i.e., $\gamma=0$, the tunnelling term is not affected by the bath in any way, thus $\Delta_{\mathrm{R}}=\Delta$.

After analysing the renormalised tunnelling term in detail, let us return to the formalism used for the two qubit case. As for the general spin case, here the system-bath coupling also induces a qubit-qubit interaction term, $\sigma_{z}^{(j)}\sigma_{z}^{(k)}$, which appears in $\widetilde{H}_{\mathrm{S}}$. The system bath coupling in the polaron picture, $\widetilde{H}_{\mathrm{SB}}$, assumes a form very different from $H_{\mathrm{SB}}$. The bath operators entering into the system-bath coupling are

[TABLE]

with D_{\pm}=\exp{\!\bigl{(}\pm\sum_{n}\frac{g_{n}}{\omega_{n}}\left(b_{n}^{\dagger}-b_{n}^{\phantom{\dagger}}\right)\bigr{)}} and $\langle D_{\pm}^{2}\rangle_{\mathrm{B}}=R$. The bath correlation functions, $C_{\ell\ell^{\prime}}(\tau)=\langle V_{\ell}(\tau)V_{\ell^{\prime}}\rangle_{\mathrm{B}}$ are given by

[TABLE]

where

[TABLE]

The diagonal correlation function, $C_{zz}(\tau)$, together with cross-correlations, e.g., $C_{xz}(\tau)$, vanish identically using the full polaron transformation described in the text.

Appendix B: Ohmic and Lorentzian baths Here we study the equilibrium entanglement induced by a common Ohmic or Lorentzian bath. These baths are represented by the spectral densities

[TABLE]

respectively. Both of these densities behave as $J(\omega)\propto\omega$ for $\omega\ll 1$, thus an essential singularity appears at $\omega=0$ in the expressions for $\Delta_{\mathrm{R}}$ and $\phi(\tau)$. Consequently, the polaron method used in the main text is not applicable for these baths as it suffers from nonphysical divergences Weiss (2012); Zwerger (1983). The integrals in $\Delta_{\mathrm{R}}$ and $\phi(\tau)$ are all divergent for any non-zero coupling strengths, thus the tunnelling element, $\Delta_{\mathrm{R}}$, is always normalised to zero. Therefore, here we only present the numerical results from the imaginary time path integral simulations.

In the spin-boson models, e.g., the one used in this work, where the bath is harmonic, the trace over the bath degrees of freedom for $\widetilde{\rho}_{\mathrm{S}}$ can be performed analytically and –within the path integral formalism– we arrive at the Feynman-Vernon influence functional Feynman and Vernon (1963); Weiss (2012). Employing the Hubbard-Stratonovich transformation and realising that the influence functional can be unravelled by an auxiliary stochastic field, $\xi(\tau)$, leads to a time-dependent Hamiltonian H_{\mathrm{st}}\sim\tfrac{1}{2}\bigl{[}\epsilon\sigma_{z}+\Delta\sigma_{x}\bigr{]}+\xi(\tau)\sigma_{z} governing the imaginary time evolution. All of the effects of the bath are accounted for by the coloured noise term, $\xi(\tau)$. The trace over the bath now corresponds to averaging the imaginary time dynamics over realisations of the noise. The primary benefit of this path integral approach is that it provides the entire reduced density matrix from a single Monte Carlo calculation and arbitrary spectral density, $J(\omega)$, may be used. In our calculations, $10^{8}-10^{11}$ samples achieved convergence. Further details of the numerical implementation can be found in Refs. Lee et al. (2012); Moix et al. (2012).

Due to the high computational cost of this method, we only study the two qubit case ($N=2$) for the weak and intermediate system-bath coupling regimes. In the following, we use $\beta=10$, $\epsilon=0$ and $\Delta=1$. The entanglement of formation and structure factor entanglement witness between two non-interacting qubits in a common Ohmic bath and the same quantities for a Lorentzian bath are shown in Fig. 7. It can be seen that the features in both cases are qualitatively similar to those obtained using a super-Ohmic spectral density in the main text. Therefore, the general observations made in our main text should not be sensitive to the spectral density of the bath.

Appendix C: Scaling peak of entanglement and structure factor Finally, we present an approximate power law scaling relation in the qudit systems. Fig. LABEL:fig:Figure7_Supplementary_ScalingOfPeakEntanglement shows the scaling of peak entanglement for high dimensional bipartite and multipartite systems. The $\mathrm{EoF}_{\mathrm{lb}}$ peaks obey a consistent power law, although some fluctuation is seen in the peak of the structure-factor-based entanglement measure for multipartite system. Relations vary between approximately linear and a square root behaviour. An analytic treatment remains the subject of further work.