Observation of time-invariant coherence in a room temperature quantum simulator

The ability to live in coherent superpositions is a signature trait of quantum systems and constitutes an irreplaceable resource for quantum-enhanced technologies. However, decoherence effects usually destroy quantum superpositions. It has been recently predicted that, in a composite quantum system exposed to dephasing noise, quantum coherence in a transversal reference basis can stay protected for indefinite time. This can occur for a class of quantum states independently of the measure used to quantify coherence, and requires no control on the system during the dynamics. Here, such an invariant coherence phenomenon is observed experimentally in two different setups based on nuclear magnetic resonance at room temperature, realising an effective quantum simulator of two- and four-qubit spin systems. Our study further reveals a novel interplay between coherence and various forms of correlations, and highlights the natural resilience of quantum effects in complex systems.

The ability to live in coherent superpositions is a signature trait of quantum systems and constitutes an irreplaceable resource for quantum-enhanced technologies. However, decoherence effects usually destroy quantum superpositions. It has been recently predicted that, in a composite quantum system exposed to dephasing noise, quantum coherence in a transversal reference basis can stay protected for indefinite time. This can occur for a class of quantum states independently of the measure used to quantify coherence, and requires no control on the system during the dynamics. Here, such an invariant coherence phenomenon is observed experimentally in two different setups based on nuclear magnetic resonance at room temperature, realising an effective quantum simulator of two-and four-qubit spin systems. Our study further reveals a novel interplay between coherence and various forms of correlations, and highlights the natural resilience of quantum effects in complex systems.
Successfully harnessing genuine nonclassical effects is predicted to herald a new wave of technological devices with a disruptive potential to supersede their conventional counterparts [1]. This prediction is now coming of age, and an international race is on to translate the power of quantum technologies into commercial applications to networked communication, computing, imaging, sensing and simulation [2]. Quantum coherence [3], which incarnates the wavelike nature of matter and the essence of quantum parallelism [4], is the primary ingredient enabling a supraclassical performance in a wide range of such applications. Its key role in quantum algorithms, optics, metrology, condensed matter physics, and nanoscale thermodynamics is actively investigated and widely recognised [5][6][7][8][9][10][11][12][13]. Furthermore, coherent quantum effects have been observed in large molecules [14] and are advocated to play a functional role in even larger biological complexes [15][16][17][18]. However, coherence is an intrinsically fragile property which typically vanishes at macroscopic scales of space, time, and temperature: the disappearance of coherence, i.e. decoherence [19], in quantum systems exposed to environmental noise is one of the major hindrances still threatening the scalability of most quantum machines. Numerous efforts have been thus invested in recent years into devising feasible control schemes to preserve coherence in open quantum systems [20], with notable examples including dynamical decoupling [21,22], quantum feedback control [23] and error correcting codes [24].
In this Letter we demonstrate a fundamentally different mechanism. We observe experimentally that quantum coherence in a composite system, whose subsystems are all affected by decoherence, can remain de facto invariant for arbitrarily long time without any external control. This phenomenon was recently predicted to occur for a particular family of initial states of quantum systems of any (however large) even num-ber of qubits [25], and is here demonstrated in a room temperature liquid-state nuclear magnetic resonance (NMR) quantum simulator [29][30][31][32][33] with two different molecules, encompassing two-qubit and four-qubit spin ensembles. After initialisation into a so-called generalised Bell diagonal state [25], the multiqubit ensemble is left to evolve under naturally occurring phase damping noise. Constant coherence in a reference basis (transversal to the noise direction) is then observed within the experimentally considered timescales up to the order of a second. Coherence is measured according to a variety of recently proposed quantifiers [26], and its permanence is verified to be measure-independent. We also reveal how coherence captures quantitatively a dynamical interplay between classical and general quantum correlations [34], while any entanglement may rapidly disappear [35]. For more general initial states, we prove theoretically that coherence can decay yet remains above a guaranteed threshold at any time, and we observe this experimentally in the two-qubit instance. The present study advances our physical understanding of the resilience of quantum effects against decoherence.
Quantum coherence manifests when a quantum system is in a superposition of multiple states taken from a reference basis. The reference basis can be indicated by the physics of the problem under investigation (e.g. one may focus on the energy eigenbasis when addressing coherence in transport phenomena and thermodynamics) or by a task for which coherence is required (e.g. the estimation of a magnetic field in a certain direction). Here, for an N-qubit system, having in mind a magnetometry setting [36], we can fix the reference basis to be the 'plus/minus' basis {|± ⊗N }, where {|± } are the eigenstates of the σ 1 Pauli operator, which describes the x component of the spin on each qubit [24]. Any state with density matrix δ diagonal in the plus/minus basis will be referred to as incoherent. According to a recently formulated resource states encoded in the 1 H and 13 C nuclear spins of Chloroform; the rf pulse (θ) µ realises a qubit rotation by θ about the spin-µ axis, J is the scalar spin-spin coupling, and time flows from left to right. B: Dynamics of the experimental states ρ(t) (magenta points) in the space of the spin-spin correlation triple c j = σ j ⊗ σ j , j = 1, 2, 3; all BD states (2) fill the light blue tetrahedron, while the subclass of states spanning the inscribed green surface are predicted to have timeinvariant coherence in the plus/minus basis according to any measure of Eq. (1) [25]. C-F: full tomographies of the experimental states as prepared at time t = 0 (C, D) and after t = 0.25 s of free evolution (E, F), recorded in the computational basis (top row) and in the plus/minus basis (bottom row). G: Evolution of (absolute values of) the correlation functions |c j | in the experimental states (points), along with theoretical predictions (solid lines) based on phase damping noise with our measured relaxation times. H: Experimental observation of time-invariant coherence (in the plus/minus basis), measured by relative entropy (red circles) [26], fidelity-based measure (blue triangles) [27], and normalised trace distance (green diamonds) [25], equal to l 1 norm [26] in BD states. The slight negative slope is due to the subdominant effect of amplitude damping. I: Experimental dynamics of coherence and all forms of correlations [28] measured via relative entropies (theoretical curves omitted for graphical clarity). In panels G-I, experimental errors due to small pulse imperfections (0.3% per pulse) result in error bars within the size of the data points.
theory [3,26,27,37], the degree of quantum coherence in the state ρ of a quantum system can be quantified in terms of how distinguishable ρ is from the set of incoherent states, where the distance D is assumed jointly convex and contractive under quantum channels, as detailed in the Supplemental Material [38]. In general, different measures of coherence induce different orderings on the space of quantum states, as it happens e.g. for entanglement or other resources. A consequence of this is that, for states of a single qubit, it is impossible to find a nontrivial noisy dynamics under which coherence is naturally preserved when measured with respect to all possible choices of D in Eq. (1). As predicted in [25], such a counterintuitive situation can occur instead for larger composite systems. Here we observe this phenomenon experimentally. Our NMR setup realises an effective quantum simulator, in which N-qubit states ρ can be prepared by manipulating the deviation matrix from the thermal equilibrium density operator of a spin ensemble [29,31], via application of radiofrequency (rf) pulses and evolution under spin interactions [24,33]. The scalability of the setup relies on availability of suitably large controllable molecules in liquid-state solutions.
We first encoded a two-qubit system in a Chloroform (CHCl 3 ) sample enriched with 13 C, where the 1 H and 13 C spin-1 2 nuclei are associated to the first and second qubit, respectively. This experiment was performed in a Varian 500 MHz liquid-NMR spectrometer at room temperature, according to the plan illustrated in Fig. 1A. The state preparation stage allowed us to initialise the system in any state obtained as a mixture of maximally entangled Bell states, that is, any Bell diagonal (BD) state [47]. These states take the form where {σ j } are the Pauli matrices and I is the identity operator on each qubit; they are completely specified by the spin-spin correlation functions c j = σ j ⊗ σ j for j = 1, 2, 3, and can be conveniently represented in the space spanned by these three parameters as depicted in Fig. 1B. We aimed to prepare specifically a BD state with initial correlation functions c 1 (0) = 1, c 2 (0) = 0.7 and c 3 (0) = −0.7, by first initialising the system in the pseudo-pure state |00 00| as described in Refs. [24,29], and then implementing the sequence of rf pulses shown in Fig. 1A with θ = π and α = arccos(−0.7) ≈ 134 • . After state preparation, the system was allowed to evolve freely during a period of time t, with t increased for each trial in increments of 2/J from 0 to 0.5 s (where J ≈ 215 Hz is the scalar spin-spin coupling constant [38]), in order to obtain the complete dynamics. In the employed setup, the two main sources of decoherence can be modelled as Markovian phase damping and generalised amplitude damping channels acting on each qubit, with characteristic relaxation times T 2 and T 1 , respectively [38]. For our system, the relaxation times were measured as T H 1 = 7.53 s, T H 2 = 0.14 s, T C 1 = 12.46 s, T C 2 = 0.90 s which implies that T H,C 1 T H,C 2 . Therefore, considering also the time domain of the experiment, only the phase damping noise can be seen to have a dominant effect.
The final stage consisted of performing full quantum state tomography for each interval of time t, following the procedure detailed in [38,48]. Instances of the reconstructed experimental states at t = 0 and t = 0.25 s are presented in Fig. 1C-F. The fidelity of the initial state with the ideal target was measured at 99.1%, testifying the high degree of accuracy of our preparation stage. We verified that the evolved state remained of the BD form (2) during the whole dynamics with fidelities above 98.5%: we could then conveniently visualise the dynamics focusing on the evolution of the spin-spin correlation triple {c j (t)}, as indicated by magenta points in Fig. 1B. The time evolution of the triple {c j (t)} is detailed in Fig. 1G.
From the acquired state tomographies during the relaxation progress, we measured the dynamics of quantum coherence in our states adopting all the known geometric coherence monotones proposed in the literature, as shown in Fig. 1H. All quantifiers were found simultaneously constant within the experimental confidence levels, revealing a universal resilience of quantum coherence in the dynamics under investigation. Note that the observed time-invariant coherence is a nontrivial feature which only occurs under particular dynamical conditions. A theoretical analysis [25,38] predicts in fact that, for all BD states evolving such that their spin-spin correlations obey the condition c 2 (t) = −c 1 (t)c 3 (t) (corresponding to the lime green surface in Fig. 1B), any valid measure of coherence as defined in Eq. (1) with respect to the plus/minus basis should remain constant at any time t. As evident from the placement of the data points in Fig. 1B, our setup realised precisely the predicted dynamical conditions for time-invariant coherence, with no further control during the relaxation. Our experiment thus demonstrated a nontrivial spontaneous occurrence of long-lived quantum coherence under Markovian dynamics.
We remark that the observed effect is distinct from the physical mechanism of long-lived singlet states also studied in NMR [49], and from an instance of decoherence-free subspace [50]. In the latter case, an open system dynamics can act effectively as a unitary evolution on a subset of quan-

FIG. 2. A:
Modified preparation stage to engineer non-BD two-qubit states.
We prepared two states ρ 1 and ρ 2 with purity 0.92 and 0.93 respectively, setting phases θ ≈ 0.94 rad, α = π/3 for ρ 1 , and θ ≈ 0.78 rad, α = π/2 for ρ 2 . The evolution and acquisition stages were as in Fig. 1A. Full tomographies of the produced states at t = 0 are presented in: C (ρ 1 , computational basis), D (ρ 1 , plus/minus basis), E (ρ 2 , computational basis), and F (ρ 2 , plus/minus basis). B: Dynamics of the relative entropy of coherence in the prepared states, along with the lower bound inferred from the evolution of their spin-spin correlation functions. The experimental errors are estimated as in Fig. 1. tum states, automatically preserving their entropy and other informational properties. In our case, the states are instead degraded with time, but only their coherence in the considered reference basis remains unaffected. We verified this by measuring other indicators of correlations [28] in our states as a function of time. Fig. 1I shows the dynamics of entanglement, classical, quantum, and total correlations (defined in [38]), as well as coherence. While entanglement is found to undergo a sudden death [35,51] at ≈ 0.21 s, a sharp transition between the decay of classical and quantum correlations is observed at the switch time t = Such a puzzling feature has been reported earlier theoretically [34,52] and experimentally [47,53,54], but here we reveal the prominent role played by coherence in this dynamical picture. Namely, coherence in the plus/minus basis is found equal to quantum correlations before t and to classical ones after t , thus remaining constant at all times. This novel interplay between coherence and correlations, observed in our natural decohering conditions, is expected to manifest for any valid choice of geometric quantifiers used to measure the involved quantities [25]; for instance, in Fig. 1I we picked all measures based on relative entropy.
One might wonder how general the reported phenomena are if the initial states differ from the BD states of Eq. (2). In [38] we prove that, given an arbitrary state ρ with spin-spin correlation functions {c j }, its coherence with respect to any basis is always larger than the coherence of the generalised BD state defined by the same correlation functions. This entails that, even if coherence in arbitrary states may decay under noise, it will stay above a threshold guaranteed by the coherence of corresponding BD states. To demonstrate this, we modified our preparation scheme to engineer more general twoqubit states ( Fig. 2A). We prepared two different pseudo-pure states ρ 1 and ρ 2 , both with matching initial correlation triple c 1 (0) = 0.95, c 2 (0) = 0.62, c 3 (0) = −0.65, within the experimental accuracy. We then measured their coherence dynamics under natural evolution as before. For both of them, we observed a decay of coherence (albeit with different rates) towards a common time-invariant lower bound, which was determined solely by the evolution of the spin-spin correlation functions, regardless of the specifics of the states (Fig. 2B).
Finally, we investigated experimentally the resilience of coherence in a larger system, composed of four logical qubits. To this aim, we performed a more advanced NMR demonstration in a BRUKER AVIII 600 MHz spectrometer equipped with a prototype 6-channel probe head, allowing full and independent control of up to 5 different nuclear spins [30,32]. We used the 13 C O -15 N-diethyl-(dimethylcarbamoyl)fluoromethylphosphonate compound [32], whose coupling topology is shown in Fig. 3A. This molecule contains 5 NMR-active spins ( 1 H, 19 F, 13 C, 31 P and 15 N), therefore we chose to decouple 15 N and encode our four-qubit system in the remaining spins. Each pair of spins k, l = {H, F, C, P} were coupled to each other by suitable scalar constants J kl [38]. We employed an 'Insensitive Nuclei Enhanced by Polarization Transfer' (INEPT)-like procedure [55] to prepare a generalised BD state with initial correlation functions c 1 (0) = 1, c 2 (0) = c 3 (0) = 0.7 [38], as detailed in Fig. 3A. After evolution in a natural phase damping  19 F, 13 C, and 31 P nuclear spins of the 13 C O -15 N-diethyl-(dimethylcarbamoyl)fluoromethyl-phosphonate molecule (whose coupling topology is illustrated as an inset) by an INEPT-like procedure, where d kl = 1/(4J kl ) and J kl is the scalar coupling between spins k and l. Light-gray rectangles denote continuous-wave pulses, used to decouple the 15 N nucleus. The dark grey bar denotes a variable pulse, applied to set the desired correlation triple {c j (0)}. Thicker (red) and thinner (blue) bars denote π and π/2 pulses, respectively; the phases of the striped π/2 pulses were cycled to construct each density matrix element. After the preparation stage, the system was left to decohere in its environment; π pulses were applied in the middle of the evolution to avoid J kl oscillations. The final π/2 pulses served to produce a detectable NMR signal in the 1 H spin channel. B: Evolution of (absolute values of) the correlation functions |c j | in the experimental states (points), along with theoretical predictions (solid lines) based on phase damping noise with an effective relaxation time T 2 ≈ 0.04 s. C: Experimental observation of time-invariant coherence (in the plus/minus basis) in the four-qubit ensemble, measured by relative entropy (red circles) and normalised trace distance (green diamonds), along with theoretical predictions (solid lines). In panels B-C, experimental errors due to pulse imperfections and coupling instabilities result in error bars within the size of the data points. environment as before, the coherence dynamics was measured by a non-tomographic detection method similarly to what was done in [56], reading out the correlation triple (see Fig. 3B) from local spin observables on the 1 H nucleus, whose spectrum exhibited the best signal-to-noise ratio [38]. The results in Fig. 3C demonstrate, albeit with a less spectacular accuracy than the two-qubit case, time-invariant coherence in the plus/minus basis in our four-qubit complex, as measured by normalised trace distance and relative entropy of coherence; the latter quantity also coincides with the global discord, a measure of multipartite quantum correlations [57,58], in generalised BD states.
In conclusion, we demonstrated experimentally in two different room temperature NMR setups that coherence, the quintessential signature of quantum mechanics [3], can resist decoherence under particular dynamical conditions, in principle with no need for external control. While only certain states feature exactly time-invariant coherence in theory, more general states were shown to maintain a guaranteed amount of coherence within the experimental timescales. These phenomena, here observed for two-and four-qubit ensembles, are predicted to occur in larger systems composed by an arbitrary (even) number of qubits [25]. It is intriguing to wonder whether biological systems such as light-harvesting complexes, in which quantum coherence effects persist under exposure to dephasing environments [16][17][18], might have evolved towards exploiting natural mechanisms for coherence protection similar to the ideal one reported here; this is a topic for further investigation [15].
While this Letter realises a proof-of-principle demonstration, our findings can impact on practical applications, specifically on noisy quantum and nanoscale technologies. In par-ticular, in quantum metrology [7], coherence in the plus/minus basis is a resource for precise estimation of frequencies or magnetic fields generated by a Hamiltonian aligned along the spin-x direction. When decoherence with a preferred transversal direction (e.g., phase damping noise) affects the estimation, as in atomic magnetometry [36,59], a quantum enhancement can be achieved by optimising the evolution time [36,60] or using error correcting techniques [61,62]. Here we observed instances in which coherence is basically unaffected by transversal dephasing noise. This suggests that the states prepared here (or others in which similar phenomena occur, such as GHZ states; see also [63]) could be used as metrological probes with sensitivity immune to decoherence. Furthermore, it has been recently shown that the quantum advantage in discriminating phase shifts generated by local spinx Hamiltonians is given exactly by the 'robustness of coherence' [13] in the plus/minus basis, a measure equal to the trace distance of coherence for BD states: this implies that the performance of such an operational task can in principle run unperturbed, if the probes are initialised as in our demonstration, in presence of a natural dephasing environment. We will explore these applications experimentally in future works.
where ω k r f is the frequency of the rf field for nucleus k (on resonance: ω k r f ≈ ω k ), ω k 1 is the nutation frequency and φ k their respective phase. As only the deviation matrix ∆ρ is affected by those unitary transformations, it is convenient to write the resulting state as ρ total = (1 − )(I A ⊗ I B ) + ρ /4, from which we define the logical NMR density matrix as the state ρ ≡ (I A ⊗ I B ) + ∆ρ /4. The state preparation procedure discussed in the main text refers to engineering ρ into any desired two-qubit BD state.
b. Decoherence processes NMR naturally provides well characterised environments, characterised by the Phase Damping (PD) and Generalised Amplitude Damping (GAD) channels acting on each qubit [24]. PD is associated to loss of coherence (in the computational basis) with no energy exchange and is specified by the following Kraus operators, where the q(t) damping function is related to the characteristic relaxation time T 2 by q(t) = (1 − e −t/T 2 ). On the other hand, the GAD channel is associated to energy exchange between system and environment and can be written in Kraus operator form by where u(t) = 1 − e −t/T 1 and p = 1/2 − α with α = ω L /2k B T . As shown in Fig. 1A of the main text, during the evolution period no refocusing pulses were applied. This implies that the phase damping function q(t) decayed naturally according to the characteristic relaxation time T 2 . We note that the effective T 2 for our experiment depends not only on the thermally induced fluctuations of longitudinal fields (standard socalled T * 2 NMR contribution) but also on static field inhomogeneities. This dependence makes the PD decoherence process occur faster, guaranteeing that no GAD effects should be effectively observed during the experiment time domain. The characteristic relaxation times were measured as T H 1 = 7.53 s, T H 2 = 0.14 s, T C 1 = 12.46 s, T C 2 = 0.90 s, which satisfy T 1 T 2 , as desired.

c. Quantum state tomography
The quantum state tomography for this two-qubit system was performed applying the simplified procedure proposed in Ref. [48]. In this case the full matrix reconstruction is obtained after performing the local operations: II, IX, IY, XX on each qubit. Here, I, X and Y, correspond, respectively, to the identity operation, a π/2 rotation around the x-axis, and a π/2 rotation around the y-axis. This set of operations provides a 16 × 16 system of equations, whose solution gives the density matrix elements.

a. NMR setup
The four-qubit experiment was performed in a BRUKER AVIII 600 MHz spectrometer equipped with a prototype 6-channel probe head as described in the main text, see [32] for further details.
The molecule chosen to encode the 4-qubit spin system was the 13 C O -15 Ndiethyl-(dimethylcarbamoyl)fluoromethyl-phosphonate compound, whose coupling topology is shown in Fig. 3A; a detailed description of its synthesis can be found again in [32]. This molecule contains 5 NMR-active spins ( 1 H, 19 F, 13 with I being the set of incoherent states, defines full coherence monotones C D (i.e., satisfying all requirements of the resource theory defined in [26]) if one chooses for example the following distance functionals: relative entropy distance [26] D RE (ρ, τ) = S (ρ||τ), where S (ρ||τ) = Tr[ρ(log(ρ) − log(τ))] is the quantum relative entropy; and fidelity based distance [27] D F (ρ, τ) = 1 − F(ρ, τ), where F(ρ, τ) = Tr √ ρτ √ ρ 2 is the Uhlmann fidelity. It is still currently unknown whether the trace distance D Tr (ρ, τ) = Tr (ρ − τ) 2 induces a full coherence monotone (in particular, it is still unclear whether property C2b of [26] is satisfied by such a measure), even though it is both contractive and jointly convex. However, for BD states the trace distance of coherence is equal to the l 1 norm of coherence (note that we adopt a normalised definition for the trace distance equal to twice the conventional one [24]), which is a full coherence monotone [25,26].
Analogously, one can define faithful measures of correlations such as total correlations, discord-type quantum correlations and entanglement in the following way.
For total correlations, where π = ρ A ⊗ τ B , with ρ A (τ B ) being an arbitrary state of subsystem A (B), form the set of product states P.
For quantum correlations [28], where χ = i j p i j |i A i A | ⊗ | j B j B |, with {p i j } being a joint probability distribution and {|i A } ({| j B }) an orthonormal basis of subsystem A (B), form the set of classical states C.
For entanglement [40], where σ = i p i ρ A i ⊗ τ B i , with {p i } being a probability distribution and ρ A i (τ B i ) arbitrary states of subsystem A (B), form the set of separable states S.
Finally, within this unifying distance-based approach, yet in a quite different way, it is also possible to identify the classical correlations of a state ρ as follows: where the first infimum is taken with respect to the subset of C formed by all the classical states χ ρ which solve the optimisation in Eq. (B3) (i.e., all the classical states which are the closest to ρ according to D), and the second infimum corresponds to the distance between each such χ ρ and the set of product states P [41][42][43].

a. Two qubits
We now outline the general conditions such that constant coherence can be observed for all time [25], constant quantum correlations can be observed up to a switch time t [52], and constant classical correlations can be observed after the switch time t , when considering two-qubit BD states undergoing local nondissipative decoherence.
Consider two qubits A and B initially in a BD state: and such that the following special initial condition is satisfied, where I is the 2 × 2 identity, σ α is the α-th Pauli matrix and {i, j, k} is a fixed chosen permutation of {1, 2, 3}. Let the two (noninteracting) qubits undergo local Markovian flip-type decoherence channels towards the k-th spin direction (i.e. bit-flip noise for k = 1, bit-phase-flip noise for k = 2, and phase-flip noise for k = 3). Their evolved global state at any time t is then represented by where q(t) = 1 − e −γt is the strength of the noise, γ is the decoherence rate, and {i, j, k} is the permutation of {1, 2, 3} fixed in the initial condition, Eq. (B7). One can easily see that the evolved state is still a BD state, whose corresponding correlation function triple is given by Focusing first on coherence, in [25] it was shown that, according to any contractive and jointly convex distance, one of the closest incoherent states δ (α) ρ(t) to the evolved BD state ρ(t), with respect to the basis consisting of tensor products of eigenstates of σ α , is just the Euclidean projection of ρ(t) onto the c α -axis, i.e. it is still a BD state and its triple is given by {δ αβ c α (t)} 3 β=1 , for any α i. Moreover, in [42] it was shown that any contractive distance between the evolved BD state ρ(t) and its Euclidean projection onto the c j -axis must be constant for any t. It immediately follows that any valid distance-based measure of coherence C ( j) D (ρ(t)) of the evolved state, with respect to the product basis consisting of tensor products of eigenstates of σ j , is invariant for any time t.
For quantum correlations, according to any contractive and jointly convex distance, one of the closest classical states χ ρ(t) to the evolved BD state ρ(t) is just the Euclidean projection of ρ(t) onto the closest c-axis, with triple given by {δ αβ c α (t)} 3 β=1 and α set by |c α (t)| = max{|c β (t)|} 3 β=1 [52]. When |c j (0)| > |c k (0)| then α = j until the switch time t = 1 2γ ln c j (0) c k (0) , with α = k afterwards. Combined with the result in [42] that any contractive distance between the evolved BD state ρ(t) and its Euclidean projection onto the c j -axis must be constant for any t, this immediately implies time-invariance of quantum correlations up until the switch time t . No general proof has yet been found for the subsequent time-invariance of classical correlations after the switch time t for any contractive and jointly convex distance, but this has been observed in particular cases (based e.g. on relative entropy, trace and fidelity based distances) [34,42,44]. The time-invariance of classical correlations is related to the finite-time emergence of the pointer basis during the dynamics [45,47].
We note that in our two-qubit experiment we have implemented exactly an instance of the above conditions, specifically in the case i = 2, j = 1 and k = 3; the corresponding phase-flip noise reduces precisely to the PD channel occurring in NMR, as it can be seen by comparing the Kraus operators in Eqs. (A4) and (B8). The decoherence rate in our demonstration was given by γ = With this evolution, the switch time is obtained as t = 1 2γ ln c 1 (0) c 3 (0) . For an initial BD state with c 1 (0) = 1, c 2 (0) = 0.7 , c 3 (0) = −0.7 as we prepared, respecting the constraint in Eq. (B7), the expected switch time was t ≈ 0.043 s, which was found in excellent agreement with the experimental data. Time-invariant coherence in the plus/minus basis (i.e. the eigenbasis of σ 1 ) was observed according to any known valid geometric measure C D (Fig. 1 of the main text). Entanglement and total correlations instead decay monotonically without experiencing any interval of timeinvariance in the considered dynamical conditions. b. Even N qubits BD states are particular instances of a more general class of N-qubit states, that we may call generalised BD states or M 3 N states [25], having all maximally mixed marginals and characterised only by the three correlations functions c α = σ ⊗N α , with j = 1, 2, 3. For any even N ≥ 2, consider N qubits initially in an M 3 N state, and such that the following initial condition is satisfied, where {i, j, k} is a fixed chosen permutation of {1, 2, 3}. Let the N (noninteracting) qubits undergo local Markovian fliptype decoherence channels towards the k-th spin direction as before (notice that such a dynamics is strictly incoherent [26,37] with respect to any product basis {|m } ⊗N , with {|m } being the eigenbasis of any of the three canonical Pauli operators σ m on each qubit). The evolved global state of the N qubits at any time t is then represented by (B13) Once again, the evolved state is still a M 3 N state, whose corresponding correlation function triple is given by Eq. (B10). Then, for any even N ≥ 2, any valid distance-based measure of coherence C ( j) D (ρ(t)) of the evolved state, with respect to the product basis consisting of tensor products of eigenstates of σ j , is invariant for any time t [25].
We note that in our four-qubit experiment (N = 4) we have implemented precisely an instance of the above conditions, specifically in the case i = 2, j = 1 and k = 3, and for an initially prepared M 3 4 state with c 1 (0) = 1, c 2 (0) = 0.7, and c 3 (0) = 0.7, respecting the constraint in Eq. (B12). Timeinvariant coherence in the plus/minus basis (i.e. the eigenbasis of σ 1 ) was observed according to various geometric measures C D (Fig. 3 of the main text).

Coherence lower bound from generalised BD states
We now prove that the coherence of an arbitrary N-qubit state ρ, with correlation functions c j = σ ⊗N j , is lower bounded by the coherence of the M 3 N state defined by the same correlation functions. This holds regardless of the number of qubits N and when considering the basis consisting of tensor products of eigenstates of σ j , for any j = 1, 2, 3 (e.g. the plus/minus basis when j = 1, as demonstrated experimentally in Fig. 2 of the main text for N = 2). Such a proof relies on the following two results.

(B16)
Second, as shown below, the map Θ is an incoherent operation with respect to the basis consisting of tensor products of eigenstates of σ j , for any j = 1, 2, 3. An incoherent operation is a CPTP map that cannot create coherence, i.e. with Kraus operators {K i } satisfying K i IK † i ⊂ I for all i, with I the set of incoherent states with respect to the chosen reference basis. Since any coherence monotone must be non-increasing under incoherent operations [26], it follows that the coherence of Θ(ρ) is less than or equal to the corresponding coherence of ρ.
In what follows we will prove that Θ is an incoherent operation with respect to the basis consisting of tensor products of eigenstates of σ 1 (i.e. the plus/minus basis), although analogous proofs hold when considering the other two Pauli operators. Since the Kraus operators of the map Θ are given by K j = 1 2 N−1 U j , in order for Θ to be an incoherent operation in the plus/minus basis we need that U j δU † j ∈ I for any δ ∈ I and any j ∈ {1, · · · , 2 2(N−1) }, where I is the set of states diagonal in this basis. This obviously holds for j = 1, being U 1 the identity. On the other hand, as it can be seen from Eq. (B15), all the other single-qubit unitaries U j are just products of the single-qubit unitaries U j listed in Eq. (B16), so that we just need to prove that U j δU † j ∈ I for any δ ∈ I and for any j ∈ {1, · · · , 2(N − 1)}. For any j ∈ {1, · · · , N − 1}, U j just leaves any state which is diagonal in the plus/minus basis invariant, it being a tensor product between two σ 1 's acting on two neighbouring qubits and the identity acting on the remaining ones. Otherwise, for any j ∈ {N, · · · , 2(N − 1)}, U j is the tensor product between two σ 2 's acting on two neighbouring qubits and the identity on the rest of the qubits. Consequently, by using σ 2 |± ±|σ 2 = |∓ ∓| and the fact that the general form of a state δ diagonal in the plus/minus basis is δ = j 1 , j 2 ,··· , j N =± p j 1 , j 2 ,··· , j N | j 1 , j 2 , · · · , j N j 1 , j 2 , · · · , j N |, we have that, when e.g. j = N, then U N δU † N = j 1 , j 2 ,··· , j N =± p j 1 , j 2 ,··· , j N U N | j 1 , j 2 , · · ·, j N j 1 , j 2 , · · ·, j N |U † N = j 1 , j 2 ,··· , j N =± p j 1 , j 2 ,··· , j N |π( j 1 ), π( j 2 ), · · ·, j N π( j 1 ), π( j 2 ), · · ·, j N |, where π(±) ≡ ∓, so that U N δU † N ∈ I. Analogously, one can see that all the remaining single-qubit local unitaries U j are such that U j δU j ∈ I, thus completing the proof.