From entanglement certification with quench dynamics to multipartite entanglement of interacting fermions

Multipartite entanglement, such as witnessed through the quantum Fisher information (QFI), is a crucial resource for quantum technologies, but its experimental certification is highly challenging. Here, we propose an experimentally friendly protocol to measure the QFI. It relies on recording the short-time dynamics of simple observables after a quench from a thermal state, works for spins, bosons, and fermions, and can be implemented in standard cold-atom experiments and other platforms with temporal control over the system Hamiltonian. To showcase the protocol, we simulate it for the one-dimensional Fermi--Hubbard model. Further, we establish a family of bounds connecting the QFI to multipartite mode entanglement for fermionic systems, which enable the detection of multipartite entanglement at sizable temperatures. Our work paves a way to experimentally accessing entanglement for quantum enhanced metrology.


Introduction.
A central question for quantum many-body physics is to understand the structure of entanglement and how it translates into observable features. Besides its potential to explain certain salient many-body phenomena [1][2][3][4][5][6], it may take a decisive role as a resource in upcoming quantum technologies. Hence, as these technologies mature, scalable protocols for detecting entanglement become increasingly necessary [7,8]. This demand is already a reality for quantum metrology [9] where the quantum Fisher information (QFI) [10], a witness for multipartite entanglement, determines the metrological quantum enhancement [11][12][13][14]. Although lower bounds of the QFI have been obtained in recent groundbreaking experiments [15][16][17][18], general and efficient procedures to directly extract its precise value in many-body systems are lacking.
To tackle this challenge, we develop an experimentally accessible technique for measuring the QFI for states in thermal equilibrium. In contrast to a previous proposal relying on frequency-dependent dynamic susceptibilities [19], our protocol only requires measuring the short-time dynamics of mean expectation values after a quench. This straightforward procedure is ideally suited, e.g., for standard experiments on ultra-cold atoms [20][21][22][23][24]. This measurement protocol for the QFI is our first main result.
Motivated by this, we derive bounds that relate multipartite fermionic mode entanglement to the QFI, by generalizing the concept of k-producibility to fermionic systems. This framework for fermionic multipartite en- tanglement is our second main result. We illustrate these bounds as well as our quench-based measurement protocol for the QFI at a paradigmatic example, the Fermi-Hubbard model in one dimension (1D). As shown in Fig. 1 and discussed further below, we certify the presence of multipartite mode entanglement for a broad region of the parameter space.
The article is organized as follows: First, we review some basic notions regarding the QFI. We proceed to derive our quench protocol. Afterwards, we rigorously define multipartite mode entanglement for fermions and determine the correct fermionic entanglement bounds for the QFI. Subsequently, we discuss the results shown in Fig. 1 in detail and provide a guideline to experiments aiming at certifying entanglement in the Fermi-Hubbard model. We conclude the article with a brief outlook.
Background on the QFI. The quantum Fisher information, F Q [ρ, O], is a central concept in quantum metrology. It quantifies the metrological sensitivity obtained from a given quantum state ρ in a phase estimation setup, in which a unitary generated by an operator O rotates ρ by an angle θ, ρ(θ) = e −iθO ρe iθO [32]. The aim in this scenario is to precisely estimate the parameter θ, whose variance after m measurements is bounded through the Cramr-Rao bound, Var θ ≥ 1/ (mF Q [ρ, O]) [10].
Moreover, the QFI witnesses multipartite entanglement. Specifically, we show this below for fermionic states ρ defined by |M | = dk + r fermionic modes. If ρ satisfies it must be, at least, (k + 1)-partite mode entangled. This result complements existing, analogous bounds for spin systems [11][12][13] so it can be calculated efficiently. However, the formula for an arbitrary density matrix ρ = λ ρ λ |λ λ|, requires diagonalizing the state, which is a challenging undertaking for quantum many-body systems, both theoretically and experimentally. In what follows, we show how to circumvent this difficulty by extracting the QFI for thermal states from expectation values using a quench.
Derivation of the quench protocol. In a previous work [19], a connection between the QFI and linear response theory was found that enables one to compute the QFI for systems in equilibrium at temperature T , This formula requires knowledge of χ"(ω, T ) = (χ(ω, T )), the imaginary part of the Fourier transform of the response function This function characterizes the linear response of an observable O(t) to a time-dependent perturbation from where H 0 is the Hamiltonian with respect to which the system was at thermal equilibrium. For deviations ∆O(t) from the equilibrium value, the Kubo formula gives [33] By transforming the integral from frequency to time domain, an equation analogous to Eq. (3) follows which allows the QFI to be obtained directly from the Kubo response function. The time domain expression has computational advantages compared to Eq. (3) and has been used for computing the QFI [19,34]. Conceptually, these expressions represent a significant advance as they explicitly relate the QFI to correlations encoded in the response functions. However, their application still presents practical problems as measurements of unequal-time correlation functions, such as [O(t), O] , are often challenging.
We overcome such limitations by introducing a protocol that solely relies on measurements of expectation values O(t) . To realize such a simplified protocol only requires a weak, abrupt quench, as can be conveniently implemented, e.g., in cold-atom experiments [21][22][23][24]. In this scenario, the drive function is simply f (τ ) = q θ(τ ), so, from Eq. (5), the dynamics are governed by Here, q denotes the quench amplitude and we introduced ξ(t, T ) = ∆O(t) quench /q . Using χ(t, T ) = dξ(t, T ) / dt in Eq. (6), we arrive at (d) Due to the functional form of κ(t, T ), the convergence is exponentially fast with a decay constant set by the tempera- protocol applies to arbitrary quench operators and quantum many-body systems, including fermionic, bosonic, and spin systems. Moreover, it has a series of advantageous properties. For example, it simplifies the requirements for extracting the QFI in many situations, as no time-time correlations are required, and the exponential decrease of κ(t, T ) = 4πT 2 [sinh (πtT ) tanh (πtT )] −1 with time implies only short measurement times are required.
Entanglement bounds for fermionic systems. Based on the concept of k-producibility, for spin systems bounds on the QFI have been derived that only states with multipartite entanglement can overcome [11][12][13]. However, for fermionic systems, such bounds do not exist.
To remedy this situation, we first need to adapt the notion of k-producibility. To see why this is necessary, we recall the condition for a state of N spins to be kproducible, |ψ spin k-prod. = |ψ 1 ⊗ |ψ 2 ⊗ · · · ⊗ |ψ P −1 ⊗ |ψ P , (9) where |ψ j is a state of N j ≤ k spins and j N j = N . Such a decomposition is not meaningful in the fermionic case due to the antisymmetric structure of the wave function. Fortunately, this also suggests what is the correct criteria, which we now introduce.
Consider a set of fermionic modes M , with associated creation and annihilation operators c † m and c m , labeled by m ∈ M . A k-partition of the system is defined as a partition M = M 1 ∪ M 2 ∪ · · · ∪ M P subject to |M j | ≤ k. Now, we introduce the following definition: a pure fermionic state |ψ is k-producible if there is a k-partition such that where the operator C j is restricted to act on M j . The C j can be written as linear combinations of products of creation operators c † m acting within M j , Here, without loss of generality, we fix some order for applying the creation operators. The summands are labeled by numbers η j (m) ∈ {0, 1}, which one can envision as the possible occupations of the modes, with associated amplitudes φ j (η j ) ∈ C. An explicit connection with the spin definition is possible if we introduce |ψ j = C j | and notice that Eq. (10) can be written as |ψ ∼ |ψ 1 ∧ |ψ 2 ∧ · · · ∧ |ψ P , with the exterior product ∧ acting as an antisymmetric analogue of the tensor product. The 1-producible decomposition with the exterior product has been used before to study mode entanglement in fermionic systems [35]. Nonetheless, for our purposes the operator language as in Eq. (11) is more convenient. The same formulation can be adapted to bosonic and spin systems, where it reproduces the usual definition of k-producible states. For 2-partite entanglement, our definition is equivalent to the one through the Slater number [36,37].
The extension of these concepts to mixed states ρ is standard [11][12][13]: a mixed state ρ k-sep. is k-separable if it can be written as a convex hull of k-producible states |λ k-prod. . This formulation introduces a hierarchy for mixed states that defines multipartite entanglement of fermionic modes: a state is (k + 1)partite mode entangled if it is not k-separable.
Using this notion, we can now establish bounds on multipartite mode entanglement. To connect to the QFI, we focus on operators of the form with w(m) ∈ R weighting the occupation of different modes. Given a k-producible state ρ = |ψ k-prod. ψ| k-prod. , one can define a probability distributions p j (η j ) = |φ j (η j )| 2 for the η j and associated ran- [38] to bound Var w j , it follows that Additional knowledge about the state |ψ k-prod. leads to restrictions on the allowed occupations η j and permits the derivation of tighter bounds for Eq. (14). In particular, if |ψ k-prod. has a fixed occupation number where we used the decomposition |M | = dk + r. See the supplementary material for a detailed discussion and tighter bounds for the case where the occupation number is known. Equation Eq. (15) and related bounds immediately extend to k-separable mixed states due to the convexity of F Q [ρ, O] . As a consequence, any state that overcomes this bound cannot be k-separable and must be, at least, (k + 1)-partite mode entangled.
Results for Fermi-Hubbard chain. We illustrate our main results on the 1D Fermi-Hubbard model, a paradigmatic model for an interacting, fermionic manybody system. Its Hamiltonian reads The fermions live on lattice sites x = 1, 2, . . . L and have two internal states, σ =↑, ↓. J governs hopping between neighboring sites and U controls on-site interactions. The Hamiltonian H 0 commutes with total occupation and magnetization, and we choose to work on the magnetization-free subspace at half-filling.
To evaluate the QFI via Eq. (8), we consider quenches using the staggered magnetization, O + , and density, O − , This choice is motivated by limit cases: at U → +∞ and half-filling, the fermions form a Nel state with homogeneous density and alternating internal state. Here, O + differentiates between the two degenerate ground states describing two possible alternating orders. For U → −∞, the fermions pair up to form a charge-density wave with homogeneous magnetization. Here, O − distinguishes two possibilities of alternating large and low density. These limiting situations can be described analytically by an effective antiferromagnetic theory [31]. Based on intuition from previous work [19], we expect O ± to give a large QFI as one goes from the free theory at U/J → 0 to the antiferromagnetic limit.
To simulate the quench protocol, we extract ξ(t, T ) from exact diagonalization and use it to calculate A quench with O± amounts to abruptly modifying the chemical potential in a staggered fashion, which can be simply implemented through superlattices, which are spin-dependent for O+, without the need for quantum gas microscopes. The relevant observable O±(t) can be measured through sitedependent imaging [39]. theory describes the system, multipartite entanglement is detected at temperatures as large as T /J = 0.4. The system-size dependence suggests the entanglement to be especially robust in this strongly interacting region, making it a prime candidate to search for experimental signatures of multipartite entanglement. Figure 4 illustrates how one can straightforwardly realize the quenches with O − and O + using optical superlattices. Ultracold atoms are now reaching stronglycorrelated many-body states of the Fermi-Hubbard model at temperatures as low as T /J = 0.25 [40][41][42], well within the region where multipartite entanglement can be detected (see Fig. 1). Moreover, as shown in Fig. 2c,d, at such temperatures the QFI converges within few hopping events (Jt 8), i.e., on time scales faster than typical decoherence rates [42]. Thus, our quench protocol enables the detection of multipartite entanglement within existing experimental setups.
Conclusion. Though discussed in the context of ultracold fermionic gases, the simplicity and generality of our protocol make it readily applicable across different platforms. It is also straightforward to replace the simple quench we have chosen by other time-dependent functions f (t), which just requires modifying the kernel function κ(t, T )(See the supplementary material for details). Recent works have studied the dynamical behavior of the quantum Fisher information after a quantum quench [43,44]. Here, turning things on their head, we have demonstrated the power of induced dynamics to extract the quantum Fisher information. Beyond the setup developed here, there is the possibility of applying our protocol to different thermodynamical ensembles [45] and even extend it outside the realm of thermodynamical states [46].

Supplementary Material
In this supplementary material, we give details on the mathematical derivations of the quench protocol, including the general expression for arbitrary driving functions f (t) as well as numerical studies of the convergence with the quench strength q . Further, we provide additional details on the derivations of the entanglement bounds as well as refined expressions for a fixed particle number.

Quench Protocol
In this section, we give further details on the derivation of the detection scheme, in particular considering for arbitrary quench protocols. Moreover, we discuss experimental issues, such as a noise, finite ramping times, and non-infinitesimal quench amplitudes.

Details on the derivation of the detection scheme for arbitrary quench protocol
To derive Eq. (6), one first notices that, in the canonical ensemble at temperature T , where ω λλ = λ − λ denotes the energy splitting between the two energy levels |λ and |λ . Combining this expression with Eq. (2) and inserting a Heaviside step function θ (t) yields All that remains to be shown is that the term in parenthesis is χ(t, T ). To see this, it is sufficient to expand Eq. (4) in the eigenbasis of the equilibrium Hamiltonian H 0 as Once Eq. (6) is available, it is possible to obtain expressions such as Eq. (8) by applying a deconvolution procedure to the Kubo formula to extract χ(t, T ) from ∆O(t) and f (t). Importantly, the procedure works for arbitrary wellbehaved quench protocols f (t). Analytically, the deconvolution can be carried out in frequency space where the Kubo formula reads ∆ O(ω) = f (ω) χ(ω, T ), which can be formally rewritten as χ(ω, T ) = f (ω) −1 ∆ O(ω). By applying an inverse Fourier transform, we obtain and their Fourier transforms should be understood as generalized functions. The inverse relation that defines f (ω) −1 encodes that φ * f * v f = φ must hold for any well behaved test function φ. Combining Eq. (6) and Eq. (S17) with the fact that, due to causality, ∆O(t) = 0 for t < 0, we obtain

Examples of different quench protocols
The simplest example is that of a delta-pulse, f (t) = q δ(t). In that case, f (ω) −1 = 1 q , so v delta-pulse (t) = δ(t) q and which just comes from the fact that χ(t, T ) = q −1 ∆O delta-pulse (t). The case presented in the main text uses a step-type quench f (t) = q θ(t) so and we recover Eq. (8), It is straightforward to check the validity of Eq. (S19) in frequency space as ωδ(ω) = 0 or, alternatively, in real time as Importantly, the deconvolution procedure discussed in the previous section is only possible if f (t) is well behaved. In particular, f (ω) must only have isolated zeros in the support of χ(ω, T ) since f (ω) −1 has to be well defined as a generalized function. A simple drive that does not fulfill this condition is f (t) = sin (ω 0 t) as f (ω) = iπ (δ(ω + ω 0 ) − δ(ω − ω 0 )) cannot be inverted. From a physical standpoint, this simply highlights that the time dependent perturbation under consideration must probe all frequencies of the Kubo response function χ(ω, T ).

Experimental considerations
In an experimental setting, where the drive function f (t) might not have a simple functional expression and both f (t) and ∆O(t) will contain some noise, one can employ a direct deconvolution procedure such as the Wiener deconvolution [47]. This yields an approximation for v f (ω) given by where S χ(n) (ω) denotes the mean power spectral density of χ (the noise).
Regarding the quench scenario, there are two potential concerns that we would like to discuss as they are relevant for experimental implementation. First, in a real experiment one does not have an ideal quench f (t) = q θ(t) but some ramp f (t) = q r(t) with a smooth function r(t). This is of no concern as long as the timescales where the ramp reaches r(t) ≈ 1 are much smaller than the relevant timescales for the system dynamics. Even when that is not the case, it is possible to account for the ramp rigorously by deriving the correct κ f for the specific ramp profile. For instance, for a well defined ramp r t0 (t) one has whereḡ(t) is a filter defined by its Fourier transformḡ(ωt 0 ) = e −iωt0 + iωt 0 1 0 dy e −iyωt0 g(y) and κ(t, T ) = 4πT 2 / [q sinh (πtT ) tanh (πtT )] is the kernel of the ideal instantaneous quench with t 0 = 0 as in Eq. (S20). For the example of a linear ramp, g(t/t 0 ) = t/t 0 , we haveḡ(ωt 0 ) = iωt 0 / e −iωt0 − 1 .
The second potential concern is that, in principle, ∆O(t) contains higher-order corrections, while we are only interested in the part that is described by linear response theory. Formally, one has while the actual term that goes into Eq. (8) is ξ(t, T ) = ∆O (1) quench (t). Fortunately, the linear part dominates for short times, so the exponential decay of the kernel function κ(t, T ) mitigates any errors coming from non-linear effects. Even more, it is possible to obtain better estimates on the linear contribution by using different values of the quench parameter q and combining the results through a polynomial fit. The simplest application relies on performing the quench with some small q and with −q ; the two measurements can then be combined to yield which removes the quadratic contributions and enables one to get accurate results over larger timescales. One can also apply this principle directly to the values of F Q [ρ, O] as any deviations from the correct value, due to higher-order terms, will also depend algebraically on q . In Fig. S5, this is shown through the convergence of the deviations to zero as q becomes infinitesimal.

Entanglement Bounds
In this section, we derive entanglement bounds for k-producible states of fixed particle number, as given in Eq. (15). For this, we start by computing the pure-state variance of the relevant operators, and find upper bounds for it through Popoviciu's inequality. By convexity of the QFI, these yield bounds also for the case of mixed quantum states. Afterwards, we show how tighter bounds can be achieved in the case of fixed fermion number, as is relevant for cold-gas experiments performed at fixed atom number.  Analogously, and we get as the crossed terms j = j cancel out. Hence, we conclude that the QFI of a k-producible pure state is given by F Q [ρ, O] = j 4Var w j as claimed in the main text.
To obtain a useful, general bound on the QFI, it is necessary to find a bound for j Var w j that depends neither on the probability distributions p j nor on the specific partitions M j , as these are state dependent, but solely on the k-producibility and the values of the w(m). To do so, we first find a bound that assumes a given partition and then optimize the bound to find the worst-case scenario. Let and it follows from Popoviciu's inequality that where the sum over M + j corresponds to the maximum of w(η j ) over all η j and the sum over M − j to the minimum. The right hand side of Eq. (S29) is an upper bound for the Eq. (14) in the main text.
The k-partitionM j that maximizes the right hand side of Eq. (S29) can be constructed by concentrating the modes with the highest weight, in absolute value, in the same partition. More explicitly, if we enumerate the modes m 1 , m 2 , . . . m |M | such that |w(m 1 )| ≥ |w(m 2 )| ≥ · · · ≥ |w(m |M | )|, then where |M | = dk + r. Hence, we can calculate specific values for a given choice of the weights to obtain the bounds, for any k-producible pure state and it follows from the convexity of the QFI [10,13] that the same holds for any k-separable mixed state.
Tighter bounds at fixed fermion number The lower and upper bounds for w(η j ), shown in Eq. (S28), cannot always be reached if there are restrictions on the occupation numbers. For instance, if we assume that |ψ ∼ |ψ 1 ∧ |ψ 2 ∧ · · · ∧ |ψ P has a fixed total occupation number N , i.e., N = m∈M c † m c m is a conserved quantity, then the states of individual partitions, |ψ j = C j | , must also have fixed occupation numbers N j such that N j = N . If that is the case, then the number of modes inM + j andM − j has to be the same and equal to N j if the η j 's should be capable of achieving the limits of Eq. (S28). While the bound Eq. (S30) still holds, it is possible to derive tighter bounds by exploiting these facts.
To obtain the improved bounds, let us again consider an arbitrary partition and divide each M j into M j = M l j ∪ M i j ∪ M u j where the lower portion M l j contains the N j modes with the lowest weights, the upper portion M u j contains the N j modes with the highest weights, and the intermediary portion M i j contains the rest of the modes. It follows that, for any η j that respects the constraint of having occupation number N j , The task is now to find the k-partitionM j that optimizes the right hand side of Eq. (S32). This can be done by concentrating as much as possible the modes with the highest weights into the same M u j 's and those with the lowest weights into the same M l j 's in a similar fashion as above. If we enumerate the modes m 1 , m 2 , . . . m |M | such that w(m 1 ) ≥ w(m 2 ) ≥ · · · ≥ w(m |M | ), then Here, for simplicity we assume k is even, but it is straightforward to correct the indexes when k is odd. The above partition already takes into account the fact that we also have to optimize over the numbers N j as the bound can only depend on the total number N . The algorithm below describes how to allocate the optimal choiceN j . It essentially consists in keeping as many partitions at half-filling as possible, with priority given to the initial ones as they contribute more to Eq. (S32):