Measurement of linear response functions in Nuclear Magnetic Resonance

We measure multi-time correlation functions of a set of Pauli operators on a two-level system, which can be used to retrieve its associated linear response functions. The two-level system is an effective spin constructed from the nuclear spins of 1H atoms in a solution of 13C-labeled chloroform. Response functions characterize the linear response of the system to a family of perturbations, allowing us to compute physical quantities such as the magnetic susceptibility of the effective spin. We use techniques exported from quantum information to measure time correlations on the two-level system. This approach requires the use of an ancillary qubit encoded in the nuclear spins of the 13C atoms and a sequence of controlled operations. Moreover, we demonstrate the ability of such a quantum platform to compute time-correlation functions of arbitrary order, which relate to higher-order corrections of perturbative methods. Particularly, we show three-time correlation functions for arbitrary times, and we also measure time correlation functions at fixed times up to tenth order.

In nature, closed quantum systems exist only as a convenient approximation. When systems are subjected to perturbations or have strong interactions with their environment, open models yield a more reliable description. A complete statistical characterization of an open quantum system unavoidably involves knowledge on the expectation value of multi-time correlations of observables, which are related to measurable quantities 1 . Time-correlation functions are at the core of optical coherence theory 2 , and can also be used for the quantum simulation of Lindbladian dynamics 3 . A plethora of physical magnitudes, such as susceptibilities and transport coefficients, can be microscopically derived in terms of time correlation functions 4,5 . In a statistical approach 4 , linear response functions represent a powerful tool to compute the susceptibility of an observable to a perturbation on the system. Such functions are constructed in terms of time-correlation functions of unperturbed observables.
Despite the ubiquity of time correlations in physics, their measurement on a quantum mechanical system is not straightforward. This difficulty lies in the fact that in quantum mechanics the measurement process disturbs the system, leaving it unreliable for a later correlated observation. Statistical descriptions typically involve an averaging of the time-correlation functions over an ensemble of particles. In such a case, it is possible to measure operator A at time t 1 over a reduced number of particles of the ensemble, and operator B at time t 2 over particles that were not perturbed by the first measurement. However, it is not always possible to perform measurements discriminating a subset of particles out of an ensemble. Moreover, nowadays, single quantum systems offer a high degree of controllability, which legitimates the interest in measuring time-correlation functions on single quantum systems. A solution to this puzzle can be found in algorithms for quantum computation. It is known that introducing an ancillary two-level system and performing a reduced set of controlled operations, time-correlation functions of a system can be reconstructed from single-time observables of the ancilla [6][7][8] . In this article, we measure n-time correlation functions for pure states, up to = n 10, in a highly-controllable quantum platform as is the case of nuclear magnetic resonance (NMR). Moreover, we frame these correlation functions in the context of linear response theory to compute physical magnitudes including the susceptibility of the system to perturbations. Finally, the scalability of the approach is shown to be efficient for multi-time correlation functions.
The measurement of n-time correlation functions plays a significant role in the linear response theory. For instance, we can microscopically derive useful quantities such as the conductivity and the susceptibility of a system, with the knowledge of 2-time correlation functions. As an illustrative example, we study the case of a spin-1/2 particle in a uniform magnetic field of strength B along the z-axis, which has a natural Hamiltonian  γ σ = − B z 0 , where γ is the gyromagnetic ratio of the particle. We assume now that a magnetic field with a sinusoidal time dependence ′ ω − B e i t 0 and arbitrary direction α perturbs the system. The Hamiltonian representation of such a situation is given by The magnetic susceptibility of the system is the deviation of the magnetic moment from its thermal expectation value as a consequence of such a perturbation. For instance, the corrected expression for the magnetic moment in , where χ α β ω , is the frequency-dependent susceptibility. From linear response theory, we learn that the susceptibility can be retrieved integrating the linear response function as ) . Moreover, the latter can be given in terms of time-correlation functions of the measured and perturbed observables, 0 0 , and the averaging is made over a thermal equilibrium ensemble. Notice that for a two level system, the thermal average can easily be reconstructed from the expectation values of the ground and excited states. So far, the response function can be retrieved by measuring the 2-time correlation functions of the unperturbed system σ σ 〈 〉 α β t ( ) and σ σ 〈 〉 β α t ( ) . It is noteworthy to mention that when α β , and it is enough to measure one of them. All in all, measuring two-time correlation functions from an ensemble of two level systems is not merely a computational result, but an actual measurement of the susceptibility of the system to arbitrary perturbations. Therefore, it gives us insights about the behavior of the system, and helps us characterize it. In a similar fashion, further corrections to the expectation values of the observables of the system will be given in terms of higher-order correlation functions. In this experiment, we will not only measure two-time correlation functions that will allow us to extract the susceptibility of the system, but we will also show that higher-order correlation functions can be obtained.

Results
Theoretical protocol. We will follow the algorithm introduced in ref. 7 to extract n-time correlation func- from a two-level quantum system, with the assistance of one ancillary qubit. Here, φ | 〉 is the quantum state of the two-level quantum system and σ α t ( ) is a time-dependent Pauli operator in the Heisenberg picture, defined as , , , and U t t ( ; ) j i is the evolution operator from time t i to t j . The considered algorithm of ref. 7 is depicted in Fig. 1, for the case where nuclear spins of 13 C and 1 H respectively encode the ancillary qubit and the two-level quantum system, and consists of the following steps: i ( ) The input state of the probe-system qubits is prepared in ( ) It follows a unitary evolution of the system qubit from t k to time + t k 1 , , which needs not be known to the experimenter. In our setup, we engineer this dynamics by decoupling qubit 13 C and 1 H, such that only the system qubit evolves under its free-energy Hamiltonian. If we were to measure the time-correlation functions for the system following a different dynamics, the corresponding Hamiltonian should be imposed on the system at this stage of the protocol, while the system and the ancilla qubits are decoupled. Then, steps ii ( ) and iii ( ) will be iterated n times, taking k from 0 to − n 1 and avoiding step iii ( ) in the last iteration. With this, all n Pauli operators will be interspersed between evolution operators with the time intervals of interest The final state of the probe-system qubits can be written as Figure 1. Two-qubit quantum circuit for measuring general n-time correlation functions. The first line is the ancilla (held by the nuclear spin of 13 C), and second line is the system qubit (held by the nuclear spin of 1 H). The blue zone between the different controlled gates α U k on the line of qubit A represents the decoupling of the 13 C nucleus from the nuclear spin of 1 H, while the latter evolves according to . The measurement of the quantities σ 〈 〉 x and σ 〈 〉 y of the ancillary qubit at the end of the circuit will directly provide the real and imaginary values of the n-time correlation function for the initial state ρ ϕ ϕ The time correlation function is then extracted as a non-diagonal operator of the ancilla, Tr( x y , such that its measurement corresponds to which is in general a complex magnitude, and where integers r and l are the occurrence numbers of Pauli operators σ y and σ z in 1 . Notice that even if between the controlled operations the system can undergo dynamics that are unknown to the experimenter, the system still needs to be controllable in order for us to be able to perform the controlled-operations.
Experimental procedures and results. We will measure n-time correlation functions of a two-level quantum system with the assistance of one ancillary qubit by implementing the quantum circuit shown in Fig. 1. Experiments are carried out using NMR [9][10][11] , where the sample used is 13 C-labeled chloroform. Nuclear spins of 13 C and 1 H encode the ancillary qubit and the two-level quantum system, respectively. With the weak coupling approximation, the internal Hamiltonian of 13 is the chemical shift, and J 12 is the J-coupling strength as illustrated in the methods, while ν 1 0 and ν 2 0 are reference frequencies of 13 C and 1 H, respectively. We set ν ν = such that the natural Hamiltonian of the system qubit is . The detuning frequency ν ∆ is chosen as hundreds of Hz to assure the selective excitation of different nuclei via hard pulses. All experiments are carried out on a Bruker AVANCE 400 MHz spectrometer at room temperature.
The spatial averaging technique can be used to transform the system into the so-called pseudo-pure state (PPS) | 〉〈 | 00 00 from the thermal equilibrium state of an NMR ensemble, which is a highly-mixed state [12][13][14] . In Step ii ( ), all controlled quantum gates α U k are chosen from the set of gates 15-17 . The notation − C U means operator U will be applied on the system qubit only if the ancilla qubit is in state represents a single-qubit rotation on qubit j along the n-axis, with the rotation angle θ. The corresponding pulse sequence can be found in the methods.
We will now apply the described algorithm to a collection of situations of physical interest. These include two-time correlation functions of a system evolving under time-independent and time-dependent Hamiltonians, as well as three-time correlation functions. In Fig. 2 we give the detailed NMR sequences employed for the measurement of the time-correlation functions in each of these cases. More especifically, Fig. 2(a) shows the experimental sequence for measuring σ σ x z , can be measured in a similar fashion by replacing the corresponding controlled quantum gates. Figure 2(b) describes the NMR sequence for measuring the three-time correlation function σ σ σ for different values of t 1 and t 2 . Here we used a pair of π pulses, which change the sign of the Hamiltonian H 0 , if t 1 is greater than t 2 . Finally, Fig. 3(c) illustrates the NMR sequence corresponding to the measurement of the time-correlation function σ σ 〈 〉 t ( ) . The dynamics corresponding to this Hamiltonian are generated by a time-dependent radio-frequency pulse applied on the resonance of the nuclear spin of 1 H, that is to say on the system qubit.
In Fig. 3, we show the measured two-time correlation functions σ σ 〈 〉 α β t ( ) (0) for a collection of α and β, and different initial states. In this experiment the two-level system was evolving under the Hamiltonian . The observed oscillations correspond to the rotation of the two-level system along the z-axis of its Bloch sphere, as dictated by the evolution Hamiltonian. Consistently, the bottom plot of Fig. 3(d) shows no oscillations, as the time-dependent operator in this case is σ t ( ) z , which is aligned with the oscillation axis. The plotted times correspond to the time-scales of the implemented dynamics. The chosen millisecond time-range is especially convenient, as the decoherence effects become significant only at longer times. In the experiment, 0  is realized by setting ν ν = Similarly, a π rotation on the second qubit is additionally needed to prepare ρ = |+〉〈+| ⊗ | 〉〈 | 1 1 in CH as the input state of the ancilla-system compound, or alternatively a π R ( /2) y 2 rotation to generate the initial state ρ = |+〉〈+| ⊗ |+〉〈+| in CH . These correlation functions are enough to retrieve the response function for a number of physical situations corresponding to different magnetic moments and applied fields. On the other hand, extracting correlation functions for initial states | 〉 0 and | 〉 1 will allow us to reconstruct such correlation functions for a thermal state of arbitrary temperature.
In Fig. 4, we show the measured time-correlation function σ σ 〈 〉 t ( ) are applied with a time interval t. A decoupling sequence Waltz-4 18-20 is used to cancel the interaction between the 13 C and 1 H nuclei during the evolution between the controlled operations. During the decoupling period, a time-dependent radio-frequency pulse is applied on the resonance of the system qubit 1 H to create the Hamiltonian, as explained above. From a physical point of view, this kind of correlations would be descriptive of a situation where the system is in a magnetic field with an intensity that is decaying exponentially in time, that is, the unperturbed system Hamiltonian turns now into a time-dependent  γ σ = − B e at y 0 . A degradation in agreement between experiment and theory is observed at the upper end of times in Fig. 4. This is due to the cumulative effects of decoherence mechanisms and the power attenuation of the employed radio-frequency pulses at long times, which results in a weak NMR response of the nuclei and as a consequence in more imprecise spectroscopic results.

Third and higher-order time correlations.
When the perturbation is not weak enough, for instance when the radiation field applied to a material is of high intensity, the response of the system might not be linear. In such situations, higher-order response functions, which depend in higher-order time-correlation functions, will be needed to account for the non-linear corrections 4,21 . For example, the second order correction to an observable B when the system suffers a perturbation of the type = + H t H AF t ( ) ( ) 0 would be given by . The method to decouple the interaction between 13 C and 1 H nuclei is Waltz-4 sequence.

B B t A t A t F t F t dt dt
in Eq. (3). The J-coupling term of Eq. (3) will be canceled by using a refocusing pulse in the circuit. Three controlled quantum gates α U 0 , β U 1 and γ U 2 should be chosen as C R ( ) . The free evolution of the 1 H nuclei between α U 0 and β U 1 is given by the evolution operator  − e i t 0 1 . Accordingly, the free evolution of the 1 H nuclei between β U 1 and γ U 2 is given by  − − e i t t ( ) 0 2 1 . However, when t 2 < t 1 , we perform the evolu-  by inverting the phase of the Hamiltonian 0  , which is realized by using a pair of π pulses at the beginning and at the end of the evolution 15 .
For testing scalability, we also measure higher-order time-correlation functions, up to = n 10. In this case, we consider the free Hamiltonian πσ = −100 z 0  and the input state φ | 〉 = . π | 〉 R (1 41 /2) 0 x , and we measure the following high-order correlation functions as a function of the correlation order n, with time intervals

 
Here, the superscripts and subscripts of  are the order of the correlation and the involved Pauli operators, respectively. Index m runs from 1 to − n ( 1)/2 for odd n, and to n/2 for even n. The quantum circuit used to measure  xx n and  xy n is based in the gradient ascent pulse engineering (GRAPE) technique 22,23 , which is designed to be robust to the static field distributions ( ⁎ T 2 process) and RF inhomogeneities. In Fig. 6, we show the measured results for high-order time-correlation functions, demonstrating the scalability of the technique and its high accuracy even for a 10-time correlation function. With this we demonstrate that high-order correlation functions are efficiently accessible in NMR via our algorithm. In general, high-order time-correlation functions correspond to high-order corrections in perturbation theories. In our example, the jagged pattern of the measured data with the order of the time-correlation function can be explained in terms of each order corresponding to measurements in different axes of the Bloch sphere.

Discussion
For systems of bigger size and complex dynamics our technique should be equally valid, and would be useful for computational purposes, when the dynamics of the system is not reproducible by classical means. In this case, the algorithm would also work with a single ancillary qubit, however, the controlled gates would pass from two-qubit gates to multi-qubit gates, which have been little studied in NMR. Nevertheless, multi-qubit gates like the Mø lmer-Sø rensen gate can always be efficiently decomposed into a circuit of c-NOT gates 7 , allowing for the measurement of multi-qubit time correlations in NMR.
For all cases here discussed, experimental data shows a high degree of agreement with the theoretical predictions. Error bars are not shown, as they are always smaller than the used dots themselves. The dephasing times T 2 for t 1 and t 2 going from 0.5 ms to 5 ms with 0.5 ms time step, showing the agreement of experimental results with theoretical predictions. The quantum circuit for measuring zyy 3  includes three controlled quantum gates U z 0 , U y 1 and U y 2 , which are experimentally implemented by hard pulses.
ScIenTIfIc RePoRTs | 7: 12797 | DOI:10.1038/s41598-017-13037-4 of our spin-qubits are of the order of seconds, while the experimental time of a whole sequence is at most of 10 ms, allowing us to ignore the effect of dephasing effects during the experiment. In our setup, the sources of errors are related to the initialization of the PPS, and imperfections in the width of the employed hard pulses. Moreover, the latter effect is cumulative and can result in a snowball effect. Additionally, factors such as RF inhomogeneities, and measurement errors, bring in a signal loss.
In summary, we have shown that the measurement of time-correlation functions of arbitrary order in NMR is an efficient task, and that it can be used to obtain the linear response function of the system. Although the linear response function could be calculated indirectly with a precise determination of the Hamiltonian parameters of the system, this experiment can be considered its first direct measurement in NMR. For systems of bigger size and complex dynamics, the indirect estimation of this magnitude would become intractable, as an analytical or numerical solution of the dynamics is always required. However, a direct detection would still be possible following the ideas demonstrated in this experiment. Not only that, in this work, we have demonstrated that such magnitudes can be experimentally retrieved with high accuracy. This will be of interest for physicist and engineers, either to characterize systems that follow computationally intractable dynamics, or to use them for computation purposes, opening the door to the quantum simulation of physical models where time correlations play a central role. It is generally accepted, that NMR platforms scale poorly, and there is no indication that this will change in the foreseeable future. However, the central ideas demonstrated in this experiment do not rely on any property which is exclusive of NMR platforms. Therefore, it is our believe that other more scalable quantum platforms may extend the protocol demonstrated here to systems of arbitrary size, where a single ancillary qubit will always suffice.
Details of Experimental implementation. Experiments are carried out using nuclear magnetic resonance (NMR), where the sample 13 C-labeled Chloroform is used as our two-qubit quantum computing processor. 13 C and 1 H in the Chloroform act as the ancillary qubit and system qubit, respectively. Figure 7 shows the molecular structure and properties of the sample. The top plot in Fig. 8 presents some experimental spectra, such as the spectra of the thermal equilibrium and pseudo pure states. The bottom plot in Fig. 8 shows NMR spectra which is created after we measure the 2-time correlation function xy n  ( = n 2). t 0 3 ms. Refocusing pulses are used to decouple the interaction between the nuclei of 13 C and 1 H, during the time intervals ∆t between gates. For this experiment, we use the GRAPE pulsed technique to implement the quantum circuit. Using hard pulses like in the previous experiments would result in poor quality of the measured data due to the high number of pulses required and the cumulative effect of their imperfections.  Here, ( ) . Moreover, any z-rotation θ R ( ) z can be decomposed in terms of rotations around the x and y axes, . On the other hand, the decoupling of the interaction between 13 C and 1 H nuclei can be realized by using refocusing pulses or the Waltz-4 sequence. Figure 8. Experimental spectra of 13 C nuclei. The blue line of the top plot shows the observed spectrum after a π/2 pulse is applied on 13 C nuclei in the thermal equilibrium state. The signal measured after applying a π/2 pulse following the preparation of the PPS is shown by the red line of the top plot. The bottom plot shows the spectrum when we measure xy n  ( = n 2). The red, blue and black lines represent the experimental spectra, fitting results and corresponding simulations, respectively.