1 Introduction↩︎

One of the key challenges in achieving quantum advantage lies in understanding how to properly measure the complexity of quantum states, a factor that fundamentally shapes the power of quantum computers and simulators. Historically, entanglement has been regarded as the most prominent manifestation of this complexity and it has long been considered the primary quantum resource [1]–[4]. However, there are many aspects beyond entanglement that contribute to the intrinsic complexity of a quantum state, and additional notions are required to get a clear assessment of quantum resources.

A central concept in this search is non-stabilizerness, or magic, defined as the distance from the convex hull of stabilizer states (stabilizer polytope). Since these latter states can be simulated efficiently on classical computers within the Clifford formalism [5]–[8], magic quantifies the difficulty to simulate classically a state. Furthermore, stabilizer states can be implemented fault-tolerantly in several quantum architectures [9]–[12], and the presence of non-Clifford resources has been shown to be essential for universal quantum computation [13], and even connected to direct measures of complexity in concrete quantum algorithms, like Shor’s factorization [14]. Hence non-stabilizerness serves as a meaningful measure of the intrinsic quantum complexity of a state.

Non-stabilizerness has gained attention in the field of quantum many-body systems [15]–[29], where it provides an additional lens for studying quantum dynamics and phase transitions, complementing the already extensive use of entanglement measures as an investigation tool.

In this work, we focus on the quantum magic of interacting fermionic systems. This investigation inevitably faces with extra layers of complexity. In particular, already non-interacting (Gaussian) fermionic states display entanglement even if they cannot generate universal quantum complexity [30]–[37]. As a consequence, fermonic Gaussian circuits describing their evolution – also known as matchgates – can be efficiently simulated with classical algorithms [38]–[40].

Despite this non-trivial character, a Gaussian state describes in principle non-interacting fermions. For interacting fermions, an approximate effective Gaussian description is obtained via the Hartree-Fock scheme, which is however expected to break down when the correlations between fermions dominate. This strongly correlated regime leads to a rich and celebrated phenomenology, ranging from high-temperature superconductivity to a variety of unconventional quantum states. In this context, non-Gaussianity, also called non-freeness [41]–[46], has emerged as a fundamental tool to quantify deviations from the set of free-fermion states, and has proven to be necessary to achieve relevant tasks in quantum information processing [36], [37], [47]–[53]. Moreover, in close analogy with magic-state injection in stabilizer circuits, suitable non-Gaussian resources are sufficient to recover universal features in random matchgate circuits [54], [55]. Finally, most previous studies of non-stabilizerness for fermions are limited to pure states, even if a proper treatment of mixed states is however necessary to connect with experimental realizations, both in cold-atom and in solid-state platforms.

In this work we address these issues by focusing on a simple, yet non-trivial, strongly interacting fermionic model, namely a two-site Hubbard model, or Hubbard dimer, where we can perform simple exact calculations of different estimates of magic and compare them with non-Gaussianity for both pure and mixed states.

The Hubbard model is a cornerstone of condensed matter physics, that provides a conceptual framework for strong correlation physics and plays a central role in the study of high-\(T_\mathrm{c}\) superconductivity and the Mott metal–insulator transition. Within Hubbard-like models, considerable effort has been devoted to characterizing entanglement [56]–[63], given its well-established relevance for quantum information processing tasks. At the same time, it has been recognized that entanglement between particles can be entirely absent at the local level, as single-site non-Gaussianity is fully determined by classical inter-flavour correlations [64], [65]. These findings indicate that entanglement alone does not exhaust the structure of correlations and complexity in interacting fermionic systems. Yet, a systematic investigation of non-stabilizerness in Hubbard systems is still missing.

The half-filled Hubbard dimer provides an ideal playground for exact calculations at both zero and finite temperature, where mixed states naturally arise. Despite the small size, two-site models feature a non-trivial behavior that mirrors that of large systems and it is often used as a testbed of approximate theories for the Hubbard [66], [67] and related models [68]–[70]. The simplicity of the model allows us to address also the non-equilibrium dynamics following a quantum quench, explicitly demonstrating how properly accounting for mixed states becomes crucial in the presence of decoherence.

By comparing different estimators of quantum magic, we uncover features that are not reflected by entanglement and non-Gaussianity, establishing non-stabilizerness as a distinct resource for characterizing fermionic mixed states. Our results establish a direct connection between strong nonlocal magic [71]–[73] and the onset of localized fermionic states, and reveal that nonlocal magic can disappear when temperature and decoherence are taken into account.

Overall, we identify the robustness of magic [12] as the most trusted quantity to quantify non-stabilizerness in mixed fermionic states, while the widely used stabilizer Rényi entropy [74] turns out to provide unfaithful results for mixed states.

Our work demonstrates that quantum magic can be successfully studied in strongly interacting fermionic systems, and it offers a different quantum resource with respect to non-Gaussianity. Our simple analysis paves the way for investigations beyond exactly solvable models, including applications in accurate numerical approaches and for richer modeling of actual correlated materials.

This article is structured as follows. After reviewing the fundamentals of non-stabilizerness in 2, and briefly describing the model and its exact solution in 3, we show our results at zero temperature, in 4, as well as at finite temperature in 5. In section 6 we focus on the dynamics of non-stabilizerness after a quantum quench in the Coulomb repulsion, and finally we devote 7 for final discussions and outlooks.

2 Preliminaries↩︎

In this section we briefly review the non-stabilizerness measures that we use throughout the text, starting from a short introduction to the Majorana-Clifford group and stabilizer states, in 2.1, which are the building blocks for the theory of magic. We then introduce the robustness of magic in 2.2, a magic monotone well suited for mixed states, and finally present the Stabilizer Renyi entropies in 2.3, highlighting their advantages and limitations when applied to mixed states.

2.1 Majorana-Clifford group and stabilizer states↩︎

We consider a system of \(N\) fermionic spin-orbitals (qubits) \(\hat{c}_i\), \(\hat{c}^\dagger_i\), which satisfy anti-commutation relations \(\{\hat{c}_i, \hat{c}^\dagger_j\}=\delta_{ij}\) and span a Hilbert space of dimension \(d=2^N\). Notice that in this context the label \(i\) includes all possible indexes, such as lattice site, spin and atomic orbital. We denote the local Pauli operators by \(\{ \mathbb{\hat{1}}_i, \hat{X}_i, \hat{Y}_i, \hat{Z}_i \}\) and define Majorana operators by the standard Jordan-Wigner mapping, as follows: \[\begin{align} \hat{\gamma}_{2i-1} & = \hat{Z}_1 \otimes ... \otimes \hat{Z}_{i-1} \otimes \hat{X}_i \otimes \mathbb{\hat{1}}_{i+1} \otimes ... \otimes \mathbb{\hat{1}}_N \\ \hat{\gamma}_{2i} & = \hat{Z}_1 \otimes ... \otimes \hat{Z}_{i-1} \otimes \hat{Y}_i \otimes \mathbb{\hat{1}}_{i+1} \otimes ... \otimes \mathbb{\hat{1}}_N, \end{align}\] where \(i=1,...,N\). Notice that these \(2N\) operators are hermitian, and satisfy fermionic anti-commutation relations \(\{\hat{\gamma}_i, \hat{\gamma}_j\}=2\delta_{ij}\). Furthermore, they can be rewritten in terms of fermionic creation and annihilation operators as \(\hat{\gamma}_{2i-1}=\hat{c}_i+\hat{c}_i^\dagger\) and \(\hat{\gamma}_{2i}=i(\hat{c}_i-\hat{c}_i^\dagger)\). Through these operators we can construct the set of \(2^{2N}\) Majorana strings: \[\hat{M}_\boldsymbol{v} = i^{\boldsymbol{v}^T\omega_L \boldsymbol{v}}\hat{\gamma}_1^{v_1}\hat{\gamma}_2^{v_2} ... \hat{\gamma}_{2N-1}^{v_{2N-1}}\hat{\gamma}_{2N}^{v_{2N}},\] where v\(\in (\mathbb{Z}_2)^{2N}\) is a binary vector of length \(2N\) whose components indicate whether the correspondent Majorana fermion is present or not, while the phase factor \(i^{\boldsymbol{v}^T\omega_L \boldsymbol{v}}\) is needed to ensure hermitianicity. The matrix \(\omega_L\) has elements equal to one in the lower triangle and zero everywhere else. It was shown in [75] that Majorana strings are in one-to-one correspondence with Pauli strings, which we denote \(P_j\) (\(1\leq j\leq4^N\)) and can be obtained as generic tensor products of local Pauli operators, namely \(P_j \in \{ \mathbb{\hat{1}}_i, \hat{X}_i, \hat{Y}_i, \hat{Z}_i \}^{\otimes N}\). For this reason, they form a complete orthogonal basis for the space of operators, since \(\text{Tr}\left[\hat{M}_\boldsymbol{v}\hat{M}_{\boldsymbol{v}\prime}\right]=d\,\delta_{\boldsymbol{v}\boldsymbol{v}\prime}\). Furthermore, we can define the Majorana-Clifford group \(\mathcal{C}_N\) as the group of unitary operators that map Majorana strings to Majorana strings, as follows \[\mathcal{C}_N=\{U: U^\dagger \hat{M}_\boldsymbol{v} U = \hat{M}_{\boldsymbol{v}\prime}\}.\]

The states that can be obtained by means of Clifford operations starting from the vacuum state \(|0\rangle^{\otimes N}\) are called stabilizer states [5]. We denote as \(\mathcal{S}_N\) the set of all pure N-qubit stabilizer states, whose number of elements is [76], [77] \[|\mathcal{S}_N| = 2^N\prod^N_{k=1}(2^k+1),\] which is super-exponential in the number of qubits. We further define the convex hull of stabilizer states, as the classical mixture of pure stabilizer states \[\text{STAB}_N = \left\{ \;\sum_{i=1}^{|\mathcal{S}_N|} p_i\sigma_i \;\Big| \;\sigma_i \in \mathcal{S}_N, \; p_i \geq 0, \; \sum_i p_i = 1 \;\right\},\] which constitutes the stabilizer polytope illustrated in 1, and contains all mixed stabilizer states [78], [79].

Clifford operations can be generated using only Hadamard, \(\pi/4\) phase, and controlled-not gates [76]. Although circuits composed exclusively of Clifford operations can produce arbitrarily large amounts of entanglement, they can nonetheless be efficiently simulated on classical computers via the Gottesman–Knill theorem [7]. As a consequence, Clifford circuits alone are not sufficient for universal quantum computation and cannot provide any quantum computational advantage [13]. The same limitation applies to quantum states belonging to the stabilizer polytope STAB\(_N\). Indeed, if a state \(\rho\) admits a convex decomposition of the form \[\rho = \sum_i p_i \sigma_i\] then one may efficiently mimic a quantum computer by classically sampling the stabilizer \(\sigma_i\) with probability \(p_i\) and simulate its subsequent evolution using the Gottesman–Knill protocol [77], [78]. Hence, mixed stabilizer states are equally useless as computational resources.

It therefore becomes essential to quantify how far a given quantum state lies from the stabilizer polytope. In this context, non-stabilizerness—often referred to as magic—emerges as an additional resource required to achieve quantum computational advantage, as it measures the distance between a given state and STAB\(_N\). Motivated by this perspective, in the following sections we focus on two measures that have been proposed to quantify magic, with particular emphasis on whether or not they are well defined for mixed quantum states.

2.2 Robustness of magic (\(L\mathcal{R}\))↩︎

The resource theory of magic can be developed in close analogy to the well-established concept of robustness of entanglement, that was put forward in Ref. [80]. This quantity generally quantifies the endurance of entanglement against noise, by accurately quantifying the minimal mixing with a separable state that is required to completely wash out entanglement. In the context of non-stabilizerness, one can proceed analogously by considering the convex hull of stabilizer states as the set of free states. Since pure stabilizer states \(\sigma_i \in \mathcal{S}_N\) form an overcomplete basis for the set of d-dimensional matrices, we can write any density matrix as an affine combination of pure stabilizer states \(\rho = \sum_{i=1}^{|\mathcal{S}_N|} x_i\sigma_i\). In this expression, the coefficients \(x_i\) of the decomposition form a quasi-probability distribution, as they satisfy \(\sum_i x_i = 1\) but may be negative. Furthermore, the vector \(\boldsymbol{x}\) is not unique. Then, the robustness of magic \(\mathcal{R}\) is defined as the minimal \(l_1\)-norm \(||\boldsymbol{x}||_1 = \sum_i|x_i|\) over all possible decompositions [12], [15], [78] \[\label{eq:32RoM} \mathcal{R}(\rho) = \min_{\boldsymbol{x}} \left\{\, ||\boldsymbol{x}||_1 \, \; \Big| \; \rho = \sum_{i=1}^{|\mathcal{S}_N|} x_i\sigma_i, \; \sigma_i \in \mathcal{S}_N \, \right\},\tag{1}\] and quantifies the minimum overlap between the density matrix \(\rho\) and the stabilizer polytope. The \(l_1\)-norm \[||\boldsymbol{x}||_1 = \sum_{i=1}^{|\mathcal{S}_N|}|x_i| = 1 + 2\sum_{i; \, x_i<0}|x_i|\] measures the amount of negativity in the affine decomposition, which has been related to the simulation runtime in the context of quantum computation. Indeed, in Monte Carlo simulation, the number of samples needed to achieve a particular accuracy scales as \(\mathcal{O}\left(\mathcal{R}(\rho)^2\right)\) [12], [81].

By collecting terms with the same sign in the above expression, one can express the density matrix as a combination of two (mixed) stabilizer states and a real positive number, as \[\rho = (1+s)\sigma_+-s\sigma_-,\] as shown in 1. The robustness of magic is then equivalently defined as \[\mathcal{R}(\rho) = \min_{\sigma_\pm\in \text{STAB}_N} \Bigg\{ \, 2s+1 \quad \Big| \quad \rho = (1+s)\sigma_+-s\sigma_-, \, s\geq0 \, \Bigg\},\] which has a clean geometrical interpretation (see 1) and represents the minimum weight of the combination of stabilizers that can reproduce the state \(\rho\).

The problem of 1 can be equivalently reformulated as the minimization over the solutions of the system of linear equations \[\mathcal{R}(\rho) = \min_{\boldsymbol{x}} \left\{\, ||\boldsymbol{x}||_1 \, \; \Big| \; \boldsymbol{A}_N \boldsymbol{x} = \boldsymbol{b} \right\},\] where we have used the unique decomposition of the quantum state onto the basis of N-qubit Pauli strings, which leads to the definitions \(b_j = \text{Tr}(\rho P_j)\) and \((\boldsymbol{A}_N)_{ij} = \text{Tr}(\sigma_i P_j)\), where \(P_j\) is the j-th Pauli string (\(1\leq j\leq4^N\)). In this way, one can relate the optimization over the large set of stabilizer states to a standard problem of linear algebra [82].

The robustness of magic is a good magic measure, in the sense that it satisfies all the required properties. Indeed, it is faithful, meaning that \(\mathcal{R}(\rho)\geq 1\), and it is equal to 1 if and only if \(\rho\) belongs to the convex hull of stabilizers. As a magic monotone, it is non-increasing under all trace preserving stabilizer channels \(\mathcal{E}\), namely \(\mathcal{R}(\mathcal{E}(\rho))\leq \mathcal{R}(\rho)\). Furthermore, it is submultiplicative, i.e. \(\mathcal{R}(\rho_1 \otimes \rho_2)\leq \mathcal{R}(\rho_1) \mathcal{R}(\rho_2)\) and convex, \(\mathcal{R}(\sum_k p_k\rho_k)\leq \sum_k |p_k|\mathcal{R}(\rho_k)\). Together, these properties ensure that \(\mathcal{R}\) is well-behaved for arbitrary quantum states.

By taking the logarithm of \(\mathcal{R}\) one obtains the log-free robustness of magic \[\label{eq:32LRoM} L\mathcal{R}(\rho) = \log_2\left(\mathcal{R}(\rho)\right)\tag{2}\] where we used the base-2 logarithm so that entropies are expressed in units of bits. Inheriting the properties of \(\mathcal{R}\), the log-free robustness of magic is also a valid measure in the context of magic resource theory. In particular, it is subadditive, meaning \(L\mathcal{R}(\rho_1 \otimes \rho_2)\leq L\mathcal{R}(\rho_1) + L\mathcal{R}(\rho_2)\), and satisfies \(L\mathcal{R}(\rho)\geq 0\), where the equality holds only for stabilizer states.

2.3 Stabilizer Renyi entropy (SRE)↩︎

Despite being a valuable measure of non-stabilizerness for mixed states, \(L\mathcal{R}(\rho)\) involves computationally hard optimizations over the large set of stabilizers, limiting its applicability to systems of size \(N\lesssim 5\) qubits. To overcome this, the stabilizer Rényi entropy (SRE) has been introduced [74], and has been proven to be significantly more tractable in large many-body systems, where it can be estimated using sampling-based techniques developed for Majorana and Pauli strings [23], [27], [83].

For a pure quantum state \(\rho\), the SRE is defined as the classical Rényi entropy of order \(\alpha\) associated with the probability distribution \(\pi_\rho(\boldsymbol{v}) = \text{Tr}^2(\rho \hat{M}_\boldsymbol{v}) / d\) [74] \[\label{eq:32SRE} \mathcal{M}_\alpha(\rho) = \frac{1}{1-\alpha} \log_2 \left[ \sum_{\boldsymbol{v} \in (\mathbb{Z}_2)^n} \pi^\alpha_\rho(\boldsymbol{v}) \right] - N,\tag{3}\] up to a constant shift. This entropy captures how broadly the quantum state \(\rho\) is distributed over the complete basis of Majorana strings \(\hat{M}_\boldsymbol{v}\).

Due to its efficiency and scalability, the SRE has been successfully applied to various many-body systems, including spin chains [16], [18]–[21], [84], fermionic states [27], [28], [85], lattice gauge theories [22], [23], neural quantum states [86], [87], matrix product states [24]–[26] and even experimental realizations [88]–[92]. Despite this broad applicability, the SRE is well defined only in the context of pure states.

While for the special case \(\alpha = 2\) the SRE has been extended to a subset of mixed states [74]: \[\label{eq:32mixedSRE} \tilde{\mathcal{M}}_2(\rho) \equiv \mathcal{M}_2(\rho) - \mathcal{S}_2(\rho) = -\log_2 \left[ \frac{ \sum_\boldsymbol{v} \text{Tr}^4\left(\rho \hat{M}_\boldsymbol{v}\right) }{ \sum_\boldsymbol{v} \text{Tr}^2\left(\rho \hat{M}_\boldsymbol{v}\right) } \right],\tag{4}\] where \(\mathcal{S}_2(\rho) = -\log_2 \text{Tr}(\rho^2)\) is the 2-Rényi entropy of \(\rho\), this correction term only partially accounts for the state’s mixedness, as we detail below.

Indeed, the expression of 4 only quantifies deviations from stabilizer states of the form [74] \[\label{eq:32stabSRE} \rho = \frac{1}{d} \left( \mathbb{1} + \sum_{\hat{M}_\boldsymbol{v} \in G} \phi_\boldsymbol{v} \hat{M}_\boldsymbol{v} \right),\tag{5}\] where \(\phi_\boldsymbol{v} = {\pm 1}\) and \(G\) is a subgroup of the full set of Majorana strings (\(|G| < d - 1\)). However, 5 fails to capture the full set of mixed stabilizers, which corresponds to the convex hull of pure stabilizer states. Consequently, 4 is a well defined magic measure only for a restricted subset of mixed states, and generally yields an overestimation when applied more broadly. In this sense, it cannot be considered a faithful measure for mixed states—a serious limitation, whose consequences will be explored in the upcoming sections.

3 Model↩︎

We consider the two-site Hubbard model at half-filling, defined by the following Hamiltonian \[\label{eq:32HDimer} \mathcal{H} = -t\sum_\sigma \left( \hat{c}^\dagger_{1,\sigma}\hat{c}_{2,\sigma} + \hat{c}^\dagger_{2,\sigma}\hat{c}_{1,\sigma} \right) + U \sum_{i=1,2}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow},\tag{6}\] often referred to as the Hubbard dimer. The operator \(\hat{c}_{i\sigma}\) (\(\hat{c}^\dagger_{i\sigma}\)) annihilates (creates) an electron with spin \(\sigma\) at the site \(i\) of the lattice, \(\hat{n}_{i\sigma} = \hat{c}^\dagger_{i\sigma} \hat{c}_{i\sigma}\) is the local spin-resolved density, \(t\) is the nearest-neighbor hopping amplitude, and \(U\) is the local Coulomb repulsion.
When extended to an infinite lattice, the Hamiltonian 6 hosts the celebrated Mott metal-to-insulator transition as the ratio \(U/t\) increases. Since in our case we only consider two sites, one cannot really speak of a phase transition, but nonetheless we can identify two distinct regimes: for \(U\ll t\) electrons will be well delocalized on the whole dimer, while for \(U\gg t\) they strongly localize on either of the two single sites.

Even though an analytical solution for the Hubbard model on an infinite lattice is currently not available – except for the limits of one and infinite dimensions, for the Hubbard dimer we can derive an exact, analytical expression for the eigenstates restricted to the sector with \(N_\uparrow=N_\downarrow=1\), since the Hamiltonian conserves the total number of particles and magnetization. In the basis \(| n_{1,\uparrow} n_{1,\downarrow} n_{2,\uparrow} n_{2,\downarrow} \rangle\) the unique ground state reads \[\label{eq:32gs} \left|\psi_-\right\rangle = \frac{1}{\mathcal{N}_+}\left( \left|\uparrow\downarrow,\circ\right\rangle + \Delta_+ \left|\uparrow,\downarrow\right\rangle - \Delta_+ \left|\downarrow,\uparrow\right\rangle + \left|\circ,\uparrow\downarrow\right\rangle \right),\tag{7}\] while the excited states are \[\label{eq:32excited95states} \begin{align} \left|\psi_+\right\rangle &= \frac{1}{\mathcal{N}_-}\left( \left|\uparrow\downarrow,\circ\right\rangle + \Delta_- \left|\uparrow,\downarrow\right\rangle - \Delta_- \left|\downarrow,\uparrow\right\rangle + \left|\circ,\uparrow\downarrow\right\rangle \right) \\ \left|D\right\rangle &= \frac{1}{\sqrt{2}}\left( \left|\uparrow\downarrow, \circ\right\rangle - \left|\circ, \uparrow\downarrow\right\rangle \right) \\ \left|t_0\right\rangle &= \frac{1}{\sqrt{2}}\left( \left|\uparrow,\downarrow\right\rangle + \left|\downarrow,\uparrow\right\rangle \right) \end{align}\tag{8}\] where we have defined \(\Delta_\pm = U / 4t \pm \sqrt{1+(U/4t)^2}\), while \(\mathcal{N_\pm}=\sqrt{2(1+\Delta_\pm^2)}\) is a suitable normalization factor. The corresponding eigenenergies are \(\{E_-, E_+, U, 0\}\), where we have defined \(E_\pm = U/2 \pm 2t\sqrt{1+(U/4t)^2}\). Finally, we point out that all eigenstates can be written in the occupation basis of the local spin-orbitals \(| n_{i,\uparrow}\rangle \otimes | n_{i,\downarrow} \rangle\), via the identification \(\left|\circ\right\rangle = \left|0\right\rangle\otimes\left|0\right\rangle\), \(\left|\uparrow\right\rangle = \left|1\right\rangle\otimes\left|0\right\rangle\), \(\left|\downarrow\right\rangle = \left|0\right\rangle\otimes\left|1\right\rangle\), \(\left|\uparrow\downarrow\right\rangle = \left|1\right\rangle\otimes\left|1\right\rangle\).

4 Zero temperature↩︎

We start by considering the system at zero temperature. Clearly, in this case all physical properties are solely determined by the pure ground state \(|\psi_-\rangle\) of 7 . For instance, the double occupancy per site, i.e. the probability for a site to be occupied simultaneously by two electrons with opposite spin, can be computed as \[\langle d \rangle = \frac{1}{2}\left \langle \psi_-|\hat{n}_{1\uparrow}\hat{n}_{1\downarrow} + \hat{n}_{2\uparrow}\hat{n}_{2\downarrow}|\psi_-\right\rangle= \frac{1}{\mathcal{N}_+^2}.\] This quantity, which is crucial to asses the mobility of electrons within the dimer, is shown in panel (d) of 2: as the Coulomb repulsion grows, double occupancies are greatly suppressed, signaling the onset of the aforementioned charge localization on the dimer.

Owing to the fact that the system under study can be considered as a collection of four qubits, we can thoroughly analyze its quantum magic for a wide range of values of \(U/t\). We display the \(L\mathcal{R}\) and the 2-SRE and 1-SRE computed on the pure ground state \(|\psi_-\rangle\) in the top panel of 2. Remarkably, all quantities exhibit a pronounced peak at intermediate values of \(U/t\), and vanish in both the weakly interacting and strongly interacting limits. This behavior aligns with physical intuition: in the non-interacting limit (\(U = 0\)), the wavefunction is completely delocalized over the dimer, and the system is well described by single-particle states. In the opposite limit (\(U \gg t\)), charge fluctuations are entirely suppressed, and the system reduces to a spin-singlet, which is a stabilizer, despite being maximally entangled. Both extremes do not support any significant degree of quantum complexity. It is therefore natural that the non-stabilizerness reaches its maximum at intermediate interaction strengths, where the system deviates most strongly from both limits. We also observe that, while the fully localized state in the limit \(U\rightarrow\infty\) does not retain any quantum magic, the maximum of non-stabilizerness occurs in conjunction with a strong depletion of doubly occupied states, indicating that the maximum non-stabilizerness is tightly bound to the onset of localization on the dimer. To provide a more complete assessment of quantum resources in the model, we compute inter-site entanglement and non-Gaussianity as functions of the Coulomb repulsion \(U\) (panels (b) and (c) of 2, respectively). While for this two-site pure ground-state one could readily evaluate the inter-site entanglement by computing the local von Neumann entropy [93]–[96], if one wants to single out the portion of the two-site entanglement that constitute a genuine quantum resource, appropriate superselection rules must be applied [56]–[60], [97]–[99]. While for mixed states this minimization has no general analytical solution, in this case, as one starts from a pure state and the reduced states are single orbitals, one obtains closed formulas for the parity and charge superselected entanglement between the sites of the Hubbard dimer [56], which we denote \(E^\text{P-SSR}\) and \(E^\text{N-SSR}\) respectively: \[\begin{align} E^\text{N-SSR} &= \left(1-2\langle d \rangle\right) \log{2} , \tag{9}\\ E^\text{P-SSR} &= \log{2}. \tag{10} \end{align}\] Analogously, non-Gaussianity, rigorously quantifies how much a state deviates from the set of non-interacting states, and is formally defined as the relative entropy between the given state and the closest Gaussian state. Remarkably, this minimization over the set of Gaussian states can be carried out analytically, so that the non-Gaussianity of an arbitrary state \(\rho\) can be expressed as [41]–[43] \[\mathcal{N}(\rho) = s(\gamma_\rho) + s(\mathbb{I}-\gamma_\rho) - s(\rho), \label{eq:32Nonfreeness}\tag{11}\] where \((\gamma_\rho)_{ij} = \langle c_j^\dagger c_i\rangle_\rho\) is the one-body density matrix of the system, and \(s(\rho)\) is the von Neumann entropy. Non-Gaussianity quantifies in a well-defined way correlations between electrons [41]–[45], constitutes a genuine resource of state-complexity [47], [100], [101] and has been connected to the information theory of many-body systems [46], [64], [65].

Our results highlight the striking difference between these quantities and magic. Both entanglement and non-Gaussianity, indeed, grow monotonically with the Coulomb repulsion, and saturate to a constant value in the limit of large \(U\). In particular, the non-Gaussianity per site, displayed in 2 (c) is zero only in the non-interacting case \(U=0\), while the entanglement (2 (b)) is everywhere larger than zero. Interestingly, both entanglement and non-Gaussianity become maximal in the large-\(U\) limit, where the ground state of the model is a stabilizer. This proves that magic, entanglement and non-Gaussianity probe complementary aspects of complexity.

4.1 Absence of local magic↩︎

We now turn our attention to an analysis of local resources in the Hubbard dimer, which provides a first insight on the unfaithfulness of the SRE. Starting from the ground state density matrix \(|\psi_-\rangle\langle\psi_-|\) of the Hubbard dimer, we can consider the density matrix of one site (its local component) by tracing out the other one, yielding the definition of the local reduced density matrix (LRDM): \[\rho_1 = \text{Tr}_2 \left(|\psi_-\rangle\langle\psi_-|\right).\] The LRDM represents the density matrix of a single fermionic orbital, and can therefore be expressed in the occupation basis of its spin-orbitals \(| n_{1,\uparrow}\rangle \otimes | n_{1,\downarrow} \rangle\), as described in 3. In this basis, \(\rho_1\) takes the form \[\begin{align} \rho_1 & = \frac{1}{\mathcal{N}_+^2} \left(|\circ\rangle\langle\circ\right| +\Delta_+^2\left|\uparrow\rangle\langle\uparrow\right| +\Delta_+^2\left|\downarrow\rangle\langle\downarrow\right| +\left|\uparrow\downarrow\rangle\left\langle\uparrow\downarrow\right|\right) \nonumber \\ & = \langle d \rangle \left(\left|\circ\rangle\langle\circ\right| + \left|\uparrow\downarrow\rangle\langle\uparrow\downarrow\right| \right) + \nonumber \\ & \qquad \qquad \qquad \qquad + \left(\frac{1}{2}-\langle d \rangle \right) \left(\left|\uparrow\rangle\langle\uparrow\right|+\left|\downarrow\rangle\langle\downarrow\right|\right), \nonumber \\ & = \langle d \rangle \left(\left|0,0\rangle\langle0,0\right| + \left|1,1\rangle\langle1,1\right| \right) + \nonumber \\ & \qquad \quad + \left(\frac{1}{2}-\langle d \rangle \right) \left(\left|1,0\rangle\langle1,0\right|+\left|0,1\rangle\langle0,1\right|\right), \label{eq:32localRDM} \end{align}\tag{12}\] where in the second line we have recasted the expression in terms of the double occupancy, using the fact that \(\langle d\rangle = 1/\mathcal{N}^2_+\) and \(1/2-\langle d\rangle = \Delta^2_+/\mathcal{N}^2_+\), and in the third line we have explicitly written it in the spin-orbital occupation basis. The form of 12 explicitly shows that the LRDM is diagonal and uniquely composed of two-qubit stabilizers, since all states that appear in \(\rho_1\) can be obtained as tensor products of \(\left|0\rangle\langle0\right|\) and \(\left|1\rangle\langle1\right|\), which are single-qubit stabilizer states. It is therefore clear that for each value of \(U/t\) we are in the presence of a classical mixture of stabilizer states, which does not host any quantum magic. The log-free robustness of magic captures this by vanishing exactly, unambiguously marking the absence of local quantum magic, as shown in 3 (a).

We stress here that the diagonal form of the LRDM is a completely general feature of Hubbard-like Hamiltonians with U\((1)\) and SU\((N)\) symmetries (\(N\) being the number of fermionic flavours, in the most popular single band Hubbard model \(N=2\)), even on an infinite lattice, and is valid also away from half filling with minor modifications [64], [65], [93]. Hence, this result generally proves the absence of local magic for Hubbard-like models, suggesting that nonlocal magic [71]–[73] characterizes these physical systems.

The vanishing of local magic aligns with the fact that entanglement is also absent at the local level. Indeed, 12 represents a mixture of separable states, therefore ruling out the possibility of intra-orbital entanglement. It follows that all correlations between spin-orbitals are due to the classical mixture of states, and are entirely captured by the local non-Gaussianity [64], [65], which grows monotonically with the Coulomb repulsion and saturates to \(\log(2)\) as shown in the bottom panel of 3. This is the maximum value for classical correlations in a 2-qubit system, in stark contrast with the \(2\log(2)\) value of the nonlocal non-Gaussianity per site in 2, which is instead characteristic of maximal quantum non-Gaussian correlations. Once again, a clear signature of the classical nature of local Hubbard states.

On the other hand, the mixed SRE \(\mathcal{\tilde{M}}_2\) behaves very differently from the \(L\mathcal{R}\). Indeed, it is everywhere greater than zero and reaches its maximum at intermediate values of \(U/t\), as shown in 3 (a). The failure of \(\mathcal{\tilde{M}}_2\) in capturing the stabilizer nature of \(\rho_1\) arises from the fact that the SRE vanishes only if the density matrix can be written in the form \(\rho_1=\tfrac{1}{4}(\hat{\mathbb{1}}+\sum_\boldsymbol{v}\phi_\boldsymbol{v}M_\boldsymbol{v})\) with \(\phi_\boldsymbol{v}=\pm 1\) (cf. 5 ), which however is not the most general mixed stabilizer state. Indeed, for any value of \(U\) the local reduced density matrix takes the form \[\rho_1 = \frac{1}{4} \left( \hat{\mathbb{1}} - \delta \hat{P} \right),\] where \(\delta = \frac{2(\Delta^2_+ - 1)}{\mathcal{N}^2_+}=1-4\langle d\rangle\) is a real number that depends on \(U/t\), not necessarily \(\pm 1\), and \(\hat{P}\) is the parity operator, i.e. the Majorana string corresponding to the vector \(\boldsymbol{v}=(1,1,1,1)\). This leads to \[\tilde{\mathcal{M}}_2(\rho_1) = -\log_2 \left( \frac{1+\delta^4}{1+\delta^2} \right),\] which vanishes only for \(\delta=0\) (at \(U=0\)) and \(\delta=1\) (as \(U \to \infty\)). Due to this clear limitation of the 2-SRE in discriminating between genuine mixed non-stabilizers and states that are convex mixtures of stabilizers, we refrain from computing it in the next section, where we focus on mixed thermal states.

5 Finite temperature↩︎

We now turn our attention to the case of a finite temperature, for which the system is described by the density matrix \[\rho = \frac{e^{-\mathcal{H}/T}}{Z}\] where \(T\) is the temperature in units of \(k_B\) and \(Z=\sum_i e^{-E_i/T}\) is the partition function. Expanding on the basis of energy eigenstates, the density matrix for fixed \(U\) and \(T\) reads \[\label{eq:rho95T} \begin{align} \rho(U,T) = \frac{1}{Z}\Big(& |\psi_-\rangle\langle\psi_-| + e^{-\frac{E_+ - E_-}{T}} |\psi_+\rangle\langle\psi_+| \\ &+ e^{-\frac{U - E_-}{T}} |D\rangle\langle D| + e^{\frac{E_-}{T}} |t_0\rangle\langle t_0| \Big), \end{align}\tag{13}\] where we used the eigenstates and eigenvalues defined in 7 8 .

In this situation, the \(L\mathcal{R}\) reveals a richer structure in the joint dependence on temperature and Coulomb repulsion, giving rise to a nontrivial phase diagram and to a thermal death of magic. For a fixed ratio \(U/t\), increasing the temperature induces a thermal mixing between the non-stabilizer states \(|\psi_-\rangle\) and \(|\psi_+\rangle\) and the stabilizer states \(|D\rangle\) and \(|0\rangle\). At sufficiently high temperatures, this mixing drives the thermal state inside the stabilizer polytope, causing the \(L\mathcal{R}\) to vanish. At first sight, one might expect this vanishing to be caused by mixing with the stabilizer states. However, this is not the case: admixture with stabilizer states can only bring the non-stabilizer ground state arbitrarily close to the stabilizer polytope, but never inside it. Instead, the crucial mechanism is the thermal mixing between the two non-stabilizer states \(|\psi_-\rangle\) and \(|\psi_+\rangle\). As we discuss more thoroughly in Appendix 8, when the thermal weight of \(|\psi_+\rangle\) becomes sufficiently large, its mixture with the ground state \(|\psi_-\rangle\) leads to a (mixed) stabilizer state, thereby eliminating magic. As a consequence, the \(L\mathcal{R}\) vanishes above a critical temperature \(T_\mathrm{c}(U/t)\), which depends on the Coulomb repulsion, as highlighted in 4 (b). Moreover, 4 (d) shows that at finite temperature magic is absent for small values of \(U/t\) and emerges only when the Coulomb repulsion exceeds a critical value \(U_\mathrm{c}\), which can be determined by inverting the relation \(T_\mathrm{c}(U)\). Physically, this behaviour reflects the fact that, at small vales of \(U\), the contribution of the excited non-stabilizer \(|\psi_+\rangle\)—which drives the thermal state inside the stabilizer polytope—remains significant, suppressing non-stabilizerness. Only for sufficiently large values of the Coulomb repulsion is the Boltzmann weight of \(|\psi_+\rangle\) effectively reduced, allowing magic to develop.

To illustrate how double occupancies shape the \(U\) vs \(T\) phase diagram, we display in 4 (c) and (e) their evolution as functions of temperature and Coulomb repulsion. The behavior of the double occupancy at high temperatures is straightforward to understand: as \(T\) increases, doubly occupied configurations become more thermally populated, causing \(\langle d \rangle\) to increase and eventually approach its noninteracting value of 0.25. As a consequence, as shown in 4 (e), double occupancies are suppressed only at low temperatures, while they are generally enhanced at higher temperatures.

At intermediate temperatures, however, the system exhibits a tendency toward localization as the temperature is raised, as shown in 4 (c). The thermal non-monotonicity of the double occupancy has been studied extensively in the half-filled Hubbard model and can be traced back to the fact that localization leads to higher entropy, and is thus thermodynamically favored at intermediate temperatures, until the thermal population of doubly occupied states becomes dominant at higher temperatures [102], [103]. This non-monotonic behavior underlies the temperature-dependent shift of the maximum of non-stabilizerness, indicated by the dashed azure line in the phase diagram and visible in 4 (d). Indeed, the zero-temperature analysis of 4 shows that the maximum robustness of magic coincides with strong localization on the dimer. For \(T\lesssim 1\), localization on the dimer is thermodynamically favored and sets in already at lower interaction strengths. As a consequence, the position of the maximum in the \(L\mathcal{R}\) shifts towards smaller values of \(U\) as the temperature is increased. Conversely, at higher temperatures, a strong depletion of double occupancies—and hence the development of large non-stabilizerness—requires progressively larger values of the Coulomb repulsion in order to overcome thermal fluctuations.

6 Dynamics after a quantum quench↩︎

Finally, we study how the quantum dynamics of the Hubbard dimer affect its non-stabilizer content. To this end, we consider a sudden quench in the Coulomb repulsion \(U\), from an initial value \(U=U_\mathrm{i}\) at times \(t<0\) to a final value \(U=U_\mathrm{f}\) at times \(t>0\). Owing to an increasing level of control in cold-atom platforms—especially in optically lattices, where the ratio \(U/t\) can be accurately tuned via the laser intensity—this protocol represents a situation of significant physical relevance. For \(t<0\) the Hamiltonian is \(\mathcal{H}_0 = \mathcal{H}(U_\mathrm{i})\) and the system is in the ground state \(|\psi_0\rangle = |\psi_-(U_\mathrm{i})\rangle\). The quench \(U_\mathrm{i} \rightarrow U_\mathrm{f}\) takes place at time \(t=0\), and for subsequent times the system evolves as \[|\psi(t)\rangle = e^{-\frac{i}{\hbar} \mathcal{H} \, t}|\psi_0\rangle,\] under the new Hamiltonian \(\mathcal{H} = \mathcal{H}(U_\mathrm{f})\). To compute the time evolution for \(t>0\) we rewrite the state \(|\psi_0\rangle\) on the basis of eigenstates of \(\mathcal{H}\). Since the states \(|D\rangle\) and \(|t_0\rangle\) do not depend on \(U\), they are simultaneously eigenstates of \(\mathcal{H}(U_\mathrm{f})\) and \(\mathcal{H}(U_\mathrm{i})\), and hence they are always orthogonal to \(|\psi_0\rangle\). Then, one is left with \[|\psi_0\rangle = \alpha |\psi_-(U_\mathrm{f})\rangle + \beta |\psi_+(U_\mathrm{f})\rangle,\] with \[\begin{align} \alpha &= \langle \psi_-(U_\mathrm{f})|\psi_0\rangle = \frac{2}{\mathcal{N}_+^0} \left( \frac{1+\Delta^0_+\Delta_+}{\mathcal{N}_+} \right), \\ \beta &= \langle \psi_+(U_\mathrm{f})|\psi_0\rangle = \frac{2}{\mathcal{N}_+^0} \left( \frac{1+\Delta^0_+\Delta_-}{\mathcal{N}_-} \right), \end{align}\] where we introduced the following shorthands: \(\mathcal{N}_\pm = \mathcal{N}_\pm(U_\mathrm{f})\), \(\mathcal{N}^0_\pm = \mathcal{N}_\pm(U_\mathrm{i})\), \(\Delta_\pm = \Delta_\pm(U_\mathrm{f})\) and \(\Delta^0_\pm = \Delta_\pm(U_\mathrm{i})\). Finally, the wavefunction evolved at time \(t\) reads \[|\psi(t)\rangle = \alpha \, e^{-\frac{i}{\hbar} E_- \, t}|\psi_-\rangle + \beta \, e^{-\frac{i}{\hbar} E_+ \, t}|\psi_+\rangle.\]

Since the system is closed, the evolution is unitary and the time dynamics exhibit coherent oscillations, which are clearly visible in 5 (a-d). Interestingly, the quench induces periodic oscillations of the double occupancies, with frequency \((E_+-E_-)/\hbar\). These oscillations, in turn, drive corresponding oscillations in the non-Gaussianity, in the N-SSR inter-site entanglement, as well as in the non-stabilizerness. Moreover, the resulting magic oscillations are twofold in the sense that within each period of the double-occupancy oscillations there are two oscillations of the non-stabilizer content, due the fact that the magic increases when \(\langle d \rangle\) deviates from an intermediate value, as shown in 5 (a).

In real experiments, however, interactions with the environment cause the system to progressively lose coherence, and the oscillatory behaviour discussed above is expected to be damped as the dephasing sets in. To account for this effect, we model the dissipative evolution in the density matrix formalism and apply the following trace-preserving map \[\begin{align} |\psi_-\rangle\langle\psi_+| &\rightarrow e^{-\Gamma t} |\psi_-\rangle\langle\psi_+| \\ |\psi_+\rangle\langle\psi_-| &\rightarrow e^{-\Gamma t} |\psi_+\rangle\langle\psi_-|, \end{align}\] which induces decoherence during the time evolution when \(\Gamma>0\). We observe that \(\Gamma^{-1}\) is the corresponding decoherence time.

For any finite value of \(\Gamma\), the amplitude of the oscillations is progressively reduced, as shown in 5 (e-h), and for \(t \gg \Gamma^{-1}\) the non-stabilizer content eventually reaches a saturation value.

A primary effect of decoherence is that the system is described by a mixed density matrix \(\rho(t)\), which at long times approaches \[\rho(t \gg \Gamma^{-1}) \simeq|\alpha|^2|\psi_-\rangle\langle\psi_-| + |\beta|^2|\psi_+\rangle\langle\psi_+|,\] corresponding to a time-independent mixture of two non-stabilizer states. As a result, all observables—including measures of non-stabilizerness—saturate to constant values. Moreover, the mixedness of the state implies that the 2-SRE is no longer a faithful measure of magic, as highlighted by the qualitative difference with the \(L\mathcal{R}\), shown in 5 (i) and (j). Interestingly, while the saturation value of the 2-SRE is always positive, this is not generally the case for the \(L\mathcal{R}\), which can vanish for specific values of the quench parameter \(U_\mathrm{f}\). As we detail in Appendix 8, this behavior originates from the fact that particular mixtures of \(|\psi_-\rangle\) and \(|\psi_+\rangle\) lie inside the stabilizer polytope, and therefore correspond to stabilizer states. Moreover, it is important to notice that unlike magic, entanglement and non-Gaussianity never saturate to zero in presence of dephasing, once again underscoring the difference between these markers of state complexity. We note in passing that the formulas for the inter-site entanglement under charge and parity SSRs given in 9 10 are only valid when the system is described by a pure state, and thus were not used for the case \(\Gamma>0\). The correct formulas in presence of mixed states have been derived in [57], and are not reported here for brevity.

One can further investigate the saturation value of the \(L\mathcal{R}\) as a function of the quench parameters \(U_\mathrm{i}\) and \(U_\mathrm{f}\), corresponding to the initial and final values of the Coulomb repulsion. This quantity is also experimentally relevant, as it can be interpreted as the long-time average of the non-stabilizer content, \[\overline{L\mathcal{R}} = \lim_{T\rightarrow\infty} \frac{1}{T}\int_0^T L\mathcal{R}(t) dt .\] The resulting phase diagram is shown in 6. Interestingly, magic vanishes over a wide range of parameters, indicating that for these values of \(U_\mathrm{i}\) and \(U_\mathrm{f}\) the quantum evolution, in the presence of decoherence, drives the system to a (mixed) stabilizer state.

7 Conclusions and outlooks↩︎

In this work we have studied the non-stabilizerness (magic) of strongly correlated fermions, using the Hubbard dimer as a paradigmatic model for interacting electrons. This simple, yet non trivial model is solved exactly, providing us with a full access to different measures of quantum magic and other quantum resources in a variety of scenarios, including zero and finite temperature, as well as the non-equilibrium dynamics after a quench of the interaction. Thermal and non-equilibrium protocols call for a proper treatment of mixed states.

In particular, we have quantified non-stabilizerness using the robustness of magic \(L\mathcal{R}\) and the Stabilizer Renyi Entropy SRE. Our results prove that the latter is not a reliable estimate of non-stabilizerness for mixed states.

Relying on symmetries of the Hubbard Hamiltonian, we proved the completely general result that local magic is absent in these systems, since the local reduced density matrix of the model is always a classical mixture of stabilizer states [64]. This property is missed by the SRE, that severely overestimates magic in the local reduced density matrix, for any eigenstate of the model.

By means of \(L\mathcal{R}\), we have shown that the nonlocal magic on the dimer grows to a maximal value in correspondence to the onset of electron localization on the dimer at zero temperature, while it becomes zero in the non-interacting and strong-interacting limits.

A finite temperature analysis, and the dynamics after a quantum quench in the presence of dephasing allow us to solidify and enrich the physical picture. In both cases, we have shown that non-stabilizerness can vanish whenever the mixture of the non-stabilizer energy eigenstates lies within the stabilizer polytope. In the former case, this leads to the thermal death of magic at large temperatures and small Coulomb repulsion while, in the latter, the dynamics can lead the state into the stabilizer polytope at large times, depending on the starting and ending strength of the interaction. Also in this case, the main features are completely missed by the SRE, as they are inherent to the mixed character of the state.

In order to fully assess the role of quantum magic as a marker for complexity in the Hubbard dimer, we have compared it with more traditional tools of quantum information theory, particularly entanglement and non-Gaussianity, which have been intensively studied in the context of strongly correlated fermions and quantum chemistry in the recent years [44]–[46], [57]–[65]. Our analysis clearly highlights that non-stabilizerness is a distinct and complementary measure of quantum complexity, necessary in order to provide a complete assessment of available quantum resources in systems of correlated fermions.

While our statement on the absence of local magic holds for the case of the single-band Hubbard model and its SU(\(N\)) version, the extension to the most general multiorbital setup is expected to introduce richer physics, which might lead to the development of non-stabilizerness even at the local level. Moreover, within non-local extensions of dynamical mean-field theory [104], it is possible to characterize the quantum complexity of finite clusters [60], [61], enabling the characterization of magic for subsystems embedded in the lattice.

Finally, recent analyses have established the notion of nonlocal magic [71]–[73], namely the amount of non-stabilizerness that cannot be removed by local basis changes, as the central character in understanding the interplay between magic and entanglement. While remarkable advances have been made towards the understanding of nonlocal magic in Gaussian fermionic systems [72], [73], the same cannot be said about interacting fermionic systems. In this framework, our results suggest that magic in the two-site Hubbard model is decidedly nonlocal, and extending this analysis to larger systems would be important to assess whether this behavior persists in more generic interacting fermionic models. All in all, our work paves the way to a more complete understanding of quantum resources in the ground and thermal states of fermionic non-Gaussian systems.

8 Mixing and the vanishing of magic↩︎

In this appendix, we highlight more explicitly the mechanisms behind both the thermal and the dynamical death of magic, shown respectively in [sec: FiniteT] [sec: Dynamics]. As temperature increases for a fixed interaction strength \(U/t\), thermal mixing occurs between the energy eigenstates of the system, namely the non-stabilizer states \(|\psi_-\rangle\) and \(|\psi_+\rangle\) and the stabilizer states \(|D\rangle\) and \(|0\rangle\). This mixing alters the structure of the density matrix, and at sufficiently high temperatures the thermal state enters the stabilizer polytope, causing the \(L\mathcal{R}\) to vanish. Interestingly, this vanishing is not primarily due to the mixing with stabilizer states, which can only bring the non-stabilizer ground state \(|\psi_-\rangle\) arbitrarily close to, but not inside, the stabilizer polytope. Rather, the crucial mechanism is the thermal population of the excited non-stabilizer state \(|\psi_+\rangle\): when its Boltzmann weight becomes sufficiently large, the resulting mixture of \(|\psi_-\rangle\) and \(|\psi_+\rangle\) produces a density matrix that lies entirely within the stabilizer polytope, thus eliminating magic.

To further illustrate this effect, we consider two families of linearly interpolated density matrices, and we compute their non-stabilizerness. In the first case, the mixture \(\rho = \lambda\,|\psi_+\rangle\langle\psi_+| + (1-\lambda)\,|\psi_-\rangle\langle\psi_-|\) shows that \(L\mathcal{R}(\rho)\) vanishes over a finite range of \(p\), indicating that the linear combination of the two non-stabilizer states lies within the stabilizer polytope for intermediate mixing ratios. In contrast, if we mix the ground state with a stabilizer state, \(\rho = \lambda\,|D\rangle\langle D| + (1-\lambda)\,|\psi_-\rangle\langle\psi_-|\), the robustness only vanishes continuously when \(\lambda=1\), meaning that the state gradually approaches the polytope but only enters it when the density matrix is fully composed of stabilizers. Our analysis is summarized in 7. This comparison highlights that it is specifically the combination of non-stabilizer states, rather than their admixture with stabilizers, that can drive the system into a stabilizer state. Consequently, the evolution of magic with temperature is closely tied to the redistribution of thermal weight among non-stabilizer eigenstates, which ultimately determines the critical temperature above which magic disappears.

References↩︎