The dynamical-decoupling-based spatiotemporal noise spectroscopy

Here we demonstrate how the standard, temporal-only, dynamical-decoupling-based noise spectroscopy method can be extended to also encompass the spatial degree of freedom. This spatiotemporal spectroscopy utilizes a system of multiple qubits arranged in a line that are undergoing pure dephasing due to environmental noise. When the qubits are driven by appropriately coordinated sequences of $\pi$ pulses the multi-qubit register becomes decoupled from all components of the noise, except for those characterized by frequencies and wavelengths specified by the pulse sequences. This allows for employment of the procedure for reconstruction of the two-dimensional spectral density that quantifies the power distribution among spatial and temporal harmonic components of the noise.


I. INTRODUCTION
Single qubits have been successfully used to characterize the frequency spectrum of environmental fluctuations affecting their coherence [1,2]. The basic principle is to subject the qubit to a periodic perturbation-a sequence of unitary operations or measurements-in order to make it sensitive only to certain frequencies of these fluctuations [3][4][5][6][7][8][9]. Currently the most commonly used method of turning a qubit into such a frequencydiscriminating sensor is to drive it with a periodic sequence of pulses that cause an effectively instantaneous π rotations of the qubit's Bloch vector. Such a sequence of rotations acts as a filter of environmental noise [10][11][12]: low frequency noise is strongly suppressed, while at discrete frequencies determined by the pulse spacing the filter has both notches, completely suppressing noise at given frequencies [13], and narrow passbands that single out frequencies of noise that are allowed to have dominant influence on qubit's dephasing [6][7][8][9]. The application of such a sequence leads then to dynamical decoupling [5,[14][15][16][17]] that on one hand extends qubits's coherence time by suppressing the influence of most of environmental fluctuations, and on the other makes the qubit sensitive to noise frequencies determined by the inverse of applied sequence period [5][6][7][8][9]. This method has been successfully implemented with many kinds of qubits: superconducting circuits [7], trapped ions [9], ultracold atoms [18], semiconductor-based quantum dots [19][20][21], donors in silicon [22], and nitrogen-vacancy centers in diamond [23,24]. The information obtained in this way led to multiple new insights into the dynamics of environments affecting various kind of qubits.
Such a wide experimental adoption of single-qubit noise spectroscopy method has motivated recent theoretical research into possibility of using multi-qubit registers to obtain more detailed information about their environment. For example, in [25,26] the authors considered multiple qubits that, although coupled to shared environmental degrees of freedom, were described as being exposed to local-and possibly correlated-noises representing fluctuations at the physical location of each qubit. Then, with an application of properly coordinated sequences of pulses one can perform spectroscopy of these local fluctuations, but more importantly, also of the cross-correlations between noises affecting different qubits. The special cases of perfectly correlated local noises (all cross-correlations equal to local autocorrelations), and of completely independent noises (vanishing cross-correlations), were considered already in first papers devoted to decoherence of multi-qubit states [27,28]. However, quantitative investigations of more general forms of cross-correlations of multiple classical [25,29] and quantum noises [26,[29][30][31] (the term "quantum noise" indicates a case when the interaction between the environment and qubits requires full quantum mechanical treatment as it cannot be approximated with an exposition to effective external stochastic process), have been subjected to closer investigation only quite recently. Such investigations are currently of high importance for ongoing and near-future development of multiqubit quantum registers for various applications. The cross-correlations of local noises are known to have significant influence on quantum metrological applications of systems of multiple entangled qubits [32][33][34][35][36][37][38][39], and quantum error correction [40,41]-a central issue for longterm prospects of quantum computation-that crucially relies on assumptions about correlations in decoherence processes of multiple qubits [42][43][44][45].
However, such a framework for description of noise probed by multi-qubit register has its drawbacks. The main one is that the local noise is defined as temporal fluctuations that are experienced by a concrete qubit placed in a given spatial location. Therefore, when the register is modified-e.g., qubits are relocated, added or removed-the set of local noises defined with respect to the original register configuration ceases to exist and is replaced by a new set that characterizes the decoherence of the modified register. Moreover, when no additional information or assumptions regarding the structure and nature of the environment are brought in form outside (e.g., see Ref. [46]), this framework in itself does not offer any means of relating the properties of local noises affecting qubits in their new configuration with the properties of noises corresponding to the previous, or any other, configuration. Intuitively this looks like a failure of the approach as it is natural to presume that the local noises are ultimately caused by the environment common among all the qubits no matter their current configuration, and thus they should be related in a way that reflects it, and the means for exhibiting those relations should be somehow build-in into the used framework.
Here we propose a subtle but consequential evolution in the description of noise probed by the multi-qubit register that have a potential to satisfy these expectations: instead of a set of local noises attached to qubits we consider a noise field that depends on both time and position. Then, the noise field evaluated at a given spatial argument is identical to local noise attached to a qubit that would be placed at this location, and the cross-correlations between local noises are equivalent to auto-correlations of the field with fixed values of spatial arguments corresponding to positions of qubits that would define these noises. The special cases of independent and perfectly correlated local noises observed in a given configuration of the register are explained by the spatial range of the noise field's auto-correlation, the socalled correlation length; the former case occurs when the range is much shorter than the inter qubit distance, while the latter case is realized when it is much longer than the length of the register. Thus, these two extreme cases are analogous to the white noise and quasi-static noise limits of temporal fluctuations. A straightforward conclusion from this discussion is that the noise field's auto-correlation in its entirety is a much richer source of information about the environmental noise than any finite set of cross-correlations of local noises. The question thus becomes: Is it possible to access this information, and if it is, how can it be accomplished? The main purpose of this paper is to provide a positive answer to this question and to present a setup in which the measurement of multi-qubit coherences is related in possibly most direct way with the field's auto-correlation, or rather with its two-dimensional Fourier transform-the spatiotemporal spectral density of the noise field.
The design of the setup is derived from solutions developed for previously considered multi-qubit spectrometers of local noises and their cross-correlations. The key element is the ability to filter particular noise frequencies together with the wavevectors (or wavelengths) of spatial fluctuations; in our setup it is achieved with coordination of pulse timings between qubits that is specific to given spatial configuration of the register. While this design of spatiotemporal frequency filter is broadly applicable to all noise fields, the detailed scheme for extracting the data on the spectral density from the multi-qubit coherence measurements that we present here works for stationary and spatially uniform fields having Gaussian statistics (i.e., fields completely specified by their average value and the auto-correlation).
The paper is organized in the following way. First, in Section II we give an overview that explains the principles of operation of the method in an intuitive way. Then, in Section III we present all the necessary formal derivations showing that our setup indeed acts as a spatiotemporal frequency filter of the noise field. In Sec. IV we discuss the contribution to decoherence that contains the spectroscopic information about the noise, and we explain how to reliably extract this contribution from the measured data. There we also develop the supporting methods that allow to deal with expected difficulties that would not be an issue for the temporal-only variant of spectroscopy, but here can make or break the whole method. This section is rather technical and can be skipped in the first reading, or if the reader is more interested in the connection between coherence of qubits and spatiotemporal correlations of the noise field, and less in details of how this connection can be practically exploited. In Sec. V we present an example of spectral density reconstruction in a numerically simulated experiment. Finally, in Sec. VI we discuss and summarize the results.

II. OVERVIEW
Consider a qubit localized at a point in space r Q , with energy splittingĤ 1Q = Ωσ z /2, whereσ z is the z-component of qubit's Pauli operator vector. It is coupled to the environmental noise via phase Hamilto-nianV 1Q = ξ(r Q , t)σ z /2, where the stochastic function ξ(r, t) describes the noise field emitted by the environment. During the evolution the qubit is subjected to a control sequence where π-pulses are applied cyclically, with fixed interval τ p between the pulses. The effect of such a driving on qubit's dynamics is most clearly captured using Heisenberg picture of qubit's ladder operator σ + = (σ x + iσ y )/2: where T is the duration of the evolution, and the symbol (. . .) represents the averaging over realizations of the stochastic process ξ. The stochastic phase accumulated for a single trajectory of the noise is given by where f (t) is the time-domain filter function that encapsulate the effects of the periodic pulse sequence described above [2,10,[46][47][48]. This filter function is depicted in Fig. 1: it has a form of a square wave, that switches between +1 and −1 at the moment that coincides with the timing of pulse application. Due to periodicity of the sequence, the resultant filter function oscillates with a well-defined frequency Therefore, f (t) acts as a frequency filter that allows to pass through only the harmonics of the noise that are commensurate with the frequency of the filter's oscillations. The idea is that, by scanning the wide range of settings of pulse sequence periods and measuring the corresponding qubit response, it is possible to reconstruct the frequency distribution of the noise, i.e., to perform the noise spectroscopy [2,[6][7][8].
The single-qubit method described above is, naturally, the temporal type of spectroscopy that can only acquire information concerning the temporal fluctuations of the noise. Here we propose a scheme for performing the full spatiotemporal spectroscopy that provide a more complete characterization of the noise field with a description of both its spatial and temporal fluctuations. The design of this type of spectrometer is derived from the results of Refs. [25,29,46], where it was shown that in order to gain access to non-trivial information of this nature one requires a probe composed of multiple spatially distributed qubits that are driven with properly coordinated pulse sequences, and it consists of three essential elements: 1. The spectrometer is to be constructed out of multiple, non-interacting qubits, labeled with index q = 1, . . . , N , that are arranged in a line-which is the natural geometry for ion trap qubits [49], and which is pursued in recent experiments with quantum dot spin qubits [50][51][52]-i.e., qubits' positions satisfy where n is the unit vector that establishes the spatial orientation of the spectrometer. The set of coordinates x q defines the qubit spatial distribution where δ(z) is the Dirac delta function. For such a geometry of the register the multiqubit phase Hamiltonian is given byV N Q = N q=1 ξ(x q n, t)σ (q) z /2, which commutes with the free HamiltonianĤ 0 = N q=1 Ω qσ 2. The qubits can be addressed with individually tailored pulse sequences that produce the qubit-wise time-domain filter functions f q (t). We require that these filter functions satisfy and that the time shift of qubit's filter is associated with its physical position via linear relation and k p is a real, adjustable parameter, that plays the role analogous to ω p .
3. The linear geometry of qubit register coupled with time-shifted pulse sequences are sufficient to produce spatial and temporal frequency filters for the noise field. However, in order to gain access to the resultant spectroscopic information it is crucial that the qubits' spatial distribution forms a pattern of periodically repeated blocks. Hence, we require that the total number of qubits is N = n s N 0 , where n s is an integer greater than one, N 0 is the number of qubits constituting a single block, and the spatial distribution can be written as Here L 0 is the distribution period, and ρ 0 (x) = N0 q=1 δ(x − x q ) is the spatial distribution of the first block of qubits. For convenience we choose the origin of the coordinate system so that 0 < x 1 < x 2 < . . . < x N0 and L 0 = x N0 + x 1 . The period of distribution can be identified with the length of the qubit block; consequently, the total length of the spectrometer is then given by The simplest example of qubit spatial distribution that adheres to (8) constitutes of N equidistant qubits, which is equivalent to spectrometer composed of N single qubit blocks. However, as we argue in Sec. V, a more complex distributions with blocks containing larger numbers of qubits allow for more precise spectroscopy in most cases.
With the above requirements met, the evolution of multi-qubit ladder operator is given by with the phase acquired for a single noise realization given by where the spatiotemporal filter function is illustrated in Fig. 2 (a).
As it is highlighted in Fig. 2 (b), when one examines the cross section along the spatial dimension of a given temporal harmonic (compare with Fig. 1) the sinusoidal waveform, with wavenumber k p , is revealed. Therefore, it stands to reason that f (x, t) would act as a frequency filter for both temporal and spatial directions; Sec. III is dedicated to proper formal proof that this is the case. This is the main qualitative results of this work.
However, due to the specifics of its physical implementation, the pulse-sequence-generated filters are burdened by a number of limitations. (i) Since a square wave itself can be decomposed into superposition of sine waves (see Fig. 1), the temporal part of the frequency filter possesses infinitely (but countably) many passbands centered at the odd multiples of the filter frequency ω p . (ii) Although the spatial waveform is sinusoidal, one cannot avoid the appearance of side passbands, because of the discrete nature of filter function in x direction. (iii) Due to finite duration of the evolution and the length of the spectrometer, both temporal and spatial passbands have finite width, proportional to T −1 in the former case, and to L −1 in the latter.
In order to overcome these imperfections of the filter, so that an accurate spatiotemporal spectroscopy can be carried out, one has to implement a generalized version of the procedure for data acquisition and analysis, that was originally developed for single-qubit spectrometers in [48]. The details of this procedure are discussed in Sec. IV where we show that the periodicity of qubit spatial distribution and of quibt-wise pulse sequences are indispensable for its successful execution.
The method of spatiotemporal spectroscopy presented in this work is most suited for characterizing stationary and spatially uniform Gaussian noise fields (or when the Gaussian approximation is adequate to describe qubitnoise coupling [2]). In such a case, all statistical proper-ties of ξ are defined by the average value and the autocorrelation function (Gaussianity), that is respectively constant and depends only on the relative position and time (stationarity and uniformity), Then, the final product of the method is reconstruction of the spectral density, that describes the average distribution of power among spatiotemporal harmonic components of the noise, The subscript n indicates that such a spectrum quantifies fluctuation only in one spatial dimension determined by the physical orientation of the spectrometer.

III. COORDINATED PULSE SEQUENCES AS A SPATIOTEMPORAL FREQUENCY FILTER
The goal of this sections is to demonstrate formally that the spatiotemporal filter function f (x, t), introduced in Eq. (12), acts as a two-dimensional passband frequency and wavenumber filter of the spectral density S n (k, ω).

A. The attenuation function
The figure of merit of noise spectroscopy is so-called attenuation function χ(L, T ), defined as (16) where the relation (11) has been invoked to express the stochastic phase in terms of spatiotemporal filter function. The attenuation function is considered as a raw data point obtained from direct measurement of the multi-qubit spectrometer response to noise field [see Eq. (10)], conditioned by the settings of applied pulse sequences. The processing of this data in order to reconstruct the spectral density is what constitutes the core of the method.
For noise field with Gaussian statistics the average in Eq. (16) can be carried out to yield a closed form expression Utilizing the definition of auto-correlation function, and its relation to spectral density, the attenuation function can be transformed into the following form wherẽ is the Fourier transform of spatiotemporal filter function restricted to the duration of the evolution and the length of the spectrometer.

B. Pulse sequence as a frequency filter
Here we formally define the spatiotemporal filter function f (x, t), and provide the result for its two-dimensional Fourier transform (19).
Since the time-domain filter function f (t), produced by the periodic sequence composed of odd 2n t − 1 number of pulses, applied over the duration that encompasses all these periods, does not appear anywhere outside the integral in Eq. (2), it is formally not required to specify it outside the interval [0, T ]; however it is convenient to do it anyway. Here we choose the definition where the filter function is extended beyond the duration of the evolution by repeating its period indefinitely, where Θ(t) is the Heaviside step function that is 1 for t > 0 and 0 otherwise. With definition (21) in hand, the qubit-wise filter functions are straightforward to model mathematically as f q (t) = f (t + ∆τ q ) (with ∆τ q ≥ 0), however the design of pulse sequences that would produce them, require some explanation, see Fig 3. The filter function with zeroshift, f (t), is the result of a sequence where a pulse is applied every τ p , for the total of 2n t − 1 pulses over the duration T = n t T 0 . For the same duration and the shift that satisfies ∆τ q = 2κτ p + δτ , with κ ∈ Integers and 0 < δτ τ p , the required sequence consists of the first pulse timed at τ p − δτ , followed by 2n t − 1 pulses applied every τ p . Finally, when ∆τ q = (2κ + 1)τ p + δτ , the pattern of pulse timings is the same as in the previous case, with an addition of one more pulse applied at t = 0, so that the shifted filter function begins with the negative value.
Now we proceed to find the Fourier transform of restricted spatiotemporal filter (19) (the details of the calculations can be found in Appendix A): (22) where ω p = π/τ p is the filter frequency and k d = 2π/L 0 is the wavenumber associated with the period of qubit spatial distribution. The shape of the passband of each harmonic component of the filter is given by and the corresponding weights equal to They can be identified with the Fourier series coefficients of the time-domain filter function [2,46,48], and of the spatial distribution of qubit block. Note how only the widths of passbands depend on the number of repetitions of qubit blocks n s or the number of pulses n t ; this is a direct consequence of periodicity of qubit spatial distribution and pulse sequences used for construction of the spectrometer. The Eq. (22) readily confirms the wavenumber and frequency filtering properties of f (x, t) described in Sec. II: (i) The frequency-ω-and wavenumber-k-dependent h functions define the temporal and spatial passbands of the filter. Their respective widths are T −1 and L −1 .
(ii) Due to harmonic distribution of f (t)-quantified by Fourier coefficients c mωp -the temporal passbands are located around odd multiples of ω p . (iii) The temporal part of the filter interferes with the spatial part, by shifting the positions of its passbands, so that for each temporal harmonic c mωp , there is a whole set of spatial passbands located around wavenumbers mk p − lk d , each one weighted by the Fourier series coefficient v lk d .

IV. SPATIOTEMPORAL NOISE SPECTROSCOPY
Below we present the detailed overview of the procedure of spectral density reconstruction with an accompanying discussion of its implementation.

A. The spectroscopic formulas
It has been demonstrated in Ref. [48] that the overlap integral of the form with F (ω) being continuous and non-negative function (like any spectral density), can be broken down into combination of terms that are subject to three distinct types of scaling laws in respect to n t . If one substitutes g m = c mωp l v lk d h nsL0 (k − mk p + lk d ), ω m = mω p , and F (ω) = S n (k, ω) in (26), this theorem also applies to the case of attenuation function χ(L, T ) = χ(n s L 0 , n t T 0 ), for which the following decomposition is obtained: where χ S ∝ n t , ∆χ 0 is independent of n t , and ∆χ T ∼ C(0, n t T 0 ). The term that scales linearly with n t , the so-called spectroscopic formula χ S , is given by It has a form of marginal spectral density S ⋆ (mω|L) (i.e., spectral density integrated over spatial degree of freedom) spanned on a frequency comb, the teeth of which are positioned at the central frequencies of temporal passbands. In the case of single-qubit noise spectroscopy where only the temporal degree of freedom is analyzed, so that S ⋆ is replaced with purely temporal spectral density, χ S is the main resource for spectrum reconstruction; it establishes a simple relation between the measured value of attenuation function for a given setting of ω p and the values of spectral density. By combining the information from a set of spectroscopic formulas obtained with properly chosen settings of filter frequencies this relation can be inverted, and thus the spectrum of temporal fluctuations is reconstructed [6,48]. In the case of multi-qubit spectrometer it is only an intermediate step, as S ⋆ can be interrogated further to extract information on spatiotemporal spectral density S n (k, ω). Terms ∆χ T and ∆χ 0 are corrections to spectroscopic formula due to finite widths of the temporal passbands [48]. They diminish accuracy of the spectroscopy, by deviating the attenuation function from the expected form of spectrum-spanned-on-frequency comb required by the reconstruction procedure to recreate spectral density with satisfactory level of fidelity. The n t -independent ∆χ 0 can be effectively eliminated using the method described in the upcoming Sec. IV B. On the other hand, the remaining ∆χ T ∼ C(0, n t T 0 ) cannot be treated in this way. Instead, it can be made negligible (in comparison to e.g., the measurement error) by exploiting the n t -scaling law it adheres to. It is expected on physical grounds that the noise field auto-correlation function has a finite range, both in temporal and spatial dimension, i.e., there exist such a scale of time t c and length x c (called the correlation time and correlation length, respectively) that Therefore, setting the number of pulses n t so that the duration T = n t T 0 is much longer than the correlation time and ∆χ T is as close to zero as C(0, T ), is sufficient to meet the requirement that the n t -dependent correction term is negligible.
Since the form of S ⋆ matches that of (26), the reasoning presented in Ref. [48] that lead to n t -scaling laws of different parts of attenuation function, can be directly applied to the marginal spectral density as well. Therefore, S ⋆ also decomposes into into three terms that correspond to (27), but with the spectrometer length L = n s L 0 replacing the duration T = n t T 0 : where the n s -independent ∆S ⋆ 0 and ∆S ⋆ L ∼ ∞ −∞ e −imωpt C(n s L 0 , t)dt, are the corrections to the spatial version of spectroscopic formula, that enables the full reconstruction of spectral density S n (k, ω). In a manner analogous to (27), meeting the condition n s L 0 ≫ x c is sufficient to treat the n sdependent correction as negligible, while the constant ∆S ⋆ 0 can be taken care of with the method of Sec. IV B, which we now proceed to discuss.
B. The linear fit method for acquiring the spectroscopic formulas In order to perform an accurate noise spectroscopy it is crucial to gain access to parts of the attenuation function χ(n s L 0 , n t T 0 ) that are described by spectroscopic formulas χ S and S ⋆ S [see Eqs. (28) and (31)]. The description ofÁlvarez-Suter method [6] for extracting the values of spectral density from χ S and S ⋆ S by deconvolving them from the frequency combs they are spanned on is postponed to Sec. IV C. Here we discuss an application of a scheme that allows to extract these key parts of χ(n s L 0 , n t T 0 ) by exploiting the differences in the scaling laws of its constituents.
The first step is to generate a dataset composed of attenuation functions gathered in a series of measurements where the number of qubit block repetitions n s and pulses n t are varied in each experiment, starting from the lowest n min s and n min t , and progressively increased to n max s and n max t . On one hand, the required control over the number of applied pulses is both conceptually and practically trivial. On the other hand, an analogous method for extending the length of the spectrometer that involves adding blocks of qubits into the system would be substantially more challenging to implement than any manipulation of the pulse sequences. Therefore, we propose a different, but effectively equivalent approach where the number of qubits is constant, and instead of physically removing or adding qubits to the system one can turn on and off the contribution from specific blocks by focusing on suitable observables. In particular, in the Heisenberg picture the operator (σ (1) evolves only due to contribution from the qubits belonging to the first r blocks. The attenuation function encoded in operator of type (32) is accessed via corresponding set of tomographic measurements performed on spectrometer qubits that were initialized in an adequately prepared state. For such an interrogation to be successful it is necessary that the state, described by the multiqubit density matrix̺, satisfies where Tr q1q2... is a trace over degrees of freedom of qubits q 1 , q 2 , . . ., and | ± z (q) are the eigenstates ofσ z . This requirement is met by, e.g., qubits initialized in a separable state̺ = |x (1) x | ± x (q) = ±| ± x (q) , for which the above matrix element equals 2 −rN0 . Much larger value of 2 −1 is obtained for an entangled state of a form̺ = |Ψ and̺ (rN0+1...N ) is an arbitrary density matrix of qubits rN 0 + 1, . . . , N . These examples show that qubits composing the spectrometer does not necessary have to be correlated, but the correlations-and in particular, the entanglement of the type present in states (34)-can be beneficial, as they increase the magnitude of the measured signal and thus enhance the relative accuracy of tomographic measurements required for extracting the raw spectroscopic data. Such states of multiple qubits were created in ion traps [49,53] and with superconducting qubits [54,55].
The gathered attenuation functions can be conveniently arranged into (n max In the second step of the procedure, for each n s (i.e., for each row of the matrix), the data points {(n min t T 0 , χ ns,n min t ), . . . , (n max t T 0 , χ ns,n max t )} are plotted on the graph. An example of such a plot is presented in Fig. 4. The initial portion of data points exhibit oscillatory behavior with visible decaying envelope-this indicates that the n t -dependent correction term ∆χ T is still active in this region; the time scale of its decay is the correlation time. Interestingly, the observation of finite frequency of oscillations of this term (which reflect the oscillations present in the auto-correlation function) indicates a non-zero characteristic frequency in the spectrum, i.e., the presence of maximum of spectral density (the spectral line) at finite frequency. Note that it would be unreasonable to expect to be able to guess how long n min t T 0 should be to completely avoid the appearance of ∆χ T , since it is not guaranteed that one possesses prior knowledge about spectral density of the noise field, which include the correlation time t c . In fact, it is the analysis performed here that allows to estimate t c (by observing the time scale of the envelope decay), and even the positions of the spectral lines (by analyzing the frequencies present in the oscillatory behavior), prior to complete reconstruction of the spectral density. Proceeding further along the plot, the duration becomes much longer than t c , ∆χ T decays to nearly zero, and data points start to follow a clear linear trend.
The third step is to fit linear function F (T ) = a ns T + b ns to the part of the graph that is free of n t -dependent correction. The intercept of the fit is the n t -independent correction and the slope is the essential part of the spectroscopic formula The intercepts b ns can be safely discarded, while the slopes a ns found for each L = n s L 0 are needed for further processing. The fourth step is analogous to the second step, but instead of attenuation function vs. the duration, the graph is created out of a data set composed of the slopes obtained previously, {(n min s L 0 , a n min s ), . . . , (n max s L 0 , a n max s )}. Since the behavior of S ⋆ as a function of spectrometer length n s L 0 have the same features as that of χ as a function of duration n t T 0 , the data points on their respective graphs would have a similar course, as it can be seen in Fig. 4. Hence, the decay and possible oscillations of n s -dependent correction should be observed, and the linear trend should emerge as n s L 0 becomes much longer than the correlation length x c . The intercept of F (L) = AL + B fitted to this trend is given by the n s -independent corrections, and it can be immediately discarded, while the slope contains the sought after spectroscopic formulas. Let us write it in a form that is fit as an input for spectral density reconstruction procedure: Here, the symmetries of spectral density S n (−k, −ω) = S n (k, ω), and Fourier coefficients c −mωp = c * mωp , v −lk d = v * lk d have been used to limit the sum over m to positive indices only. The arguments of A(k p , ω p ) have been shown explicitly to underline that the slope was generated for a given choice of spatiotemporal filter settings.
The method of linear fit presented here, aside from its main purpose of producing the spectroscopic formulas, also serves as a self-diagnostic tool for the spatiotemporal spectroscopy as a whole. Indeed, if the collected data (35) does not follow the patterns similar to what is shown in Fig. 4, it has to be interpreted as a warning that some of the founding assumptions of spectroscopy were not upheld. For example, if the linear trend does not emerge, even for a very long duration or spectrometer length, it would mean that the noise field is not stationary or is not uniform, and consequently, the spectral density S n (k, ω) does not even exist. Only if the linear fits are possible to be performed with satisfactory precision, one can have a significant degree of confidence that the reconstructed spectrum is an accurate representation of the real thing.

C. TheÁlvarez-Suter method for spectral density reconstruction
The final step of the reconstruction procedure is the deconvolution of the values of spectral density from the frequency and wavenumber combs found in spectroscopic formulas, or more precisely, in the slope (37). This is achieved with theÁlvarez-Suter method [2,6,48] which takes as an input a dataset composed of slopes A(k p , ω p ) generated for a choice of settings of k p and ω p dictated by the shape of temporal and spatial combs. In particular, the required sets of k p and ω p constitute of an arbitrary primary wavenumber k 0 and frequency ω 0 , that are then supplemented with auxiliary settings: ω p = ω 0 , 3ω 0 , . . . , m c ω 0 , and for each multiple of primary frequency mω 0 a subset of wavenumbers, k p = mk 0 , m(k 0 + k d ), . . . , m(k 0 + l c k d ). The slopes produced in such a way can be arranged in a matrix form: m=1,3,...,mc; l=0,...,lc) By substituting the explicit form of Fourier coefficients (24) into Eq. (37), one obtains the following expression for the dataset slopes: in which the linear relation between the spectral density and the input slopes A and δ i,j is the Kronecker delta, that equals 1 if i = j and 0 otherwise. The objective is to invert the relation (39), s (m) = (V (m) ) −1 [(U −1 A (k0,ω0) ) m,l ] (l=0,...,lc) , so that the vector of spectral density values for each m is extracted from the input dataset. Strictly speaking such an operation is impossible to execute, because U and V (m) are rectangular matrices with one finite and one infinite dimension, and as such are not invertible. Therefore, the only way to solve this problem is to adopt the approximation where the infinite sums in (39) are truncated in such a way that they can be rewritten in terms of finite square matricesŨ andṼ (m) , while the reminder is relatively small in comparison to the retained portion of the series. In the case of the temporal comb, for which the Fourier coefficients exhibit the power-law decay |c mωp | 2 ∼ m −2 , the optimal approximation for the sum over m ′′ is also the simplest-the cut-off at index m c , so thatŨ = [U m,m ′ ] (m=1,3...,mc;m ′ =1,3,...,mc) . Due to the choice of auxiliary settings of the frequencies ω p the inverse matrixŨ −1 is guaranteed to exist and it can be calculated analytically [48]. On the other hand, the structure of spatial combs is qualitatively different; in contrast to the temporal comb, the spatial Fourier coefficients |v lk d | 2 are non-monotonic functions of their indices. Indeed, since each coefficient is a combination of finite number of oscillating terms [see Eq. l,l ′ ] (l=0,...,lc;l ′ ∈Lm) . The simplest example is L m = {−l c /2, −l c /2 + 1, . . . , l c /2}, which corresponds to approximating the spatial comb with a segment spanned on an interval [k 0 − l c k d /2, k 0 + l c k d /2]. A more sophisticated example, where more emphasis is placed on the distribution of spatial coefficients |v lk d | 2 , would be to choose a set of indices of l c + 1 largest comb weights. The goal of this strategy would be to prevent the enhancement of the error due to overlap of the comb remainder and parts of spectral density in the case when it is difficult to estimate the positions and widths of spectral lines. WhetherṼ (m) for a given choice of L m are also invertible depends on the Fourier coefficients v lk d , that are in turn determined by the spatial distribution of the qubit block. Assuming that (Ṽ (m) ) −1 does exist, at least for some of the indices m, the reconstruction of spectral density can be finalized: for l ∈ L m = {l

V. EXAMPLE OF SPATIOTEMPORAL SPECTRAL DENSITY RECONSTRUCTION
We illustrate the applicability of the proposed spatiotemporal noise spectroscopy method with the reconstruction of the spectral density following the procedure outlined in Sec. IV. For this purpose we choose a simple model in which the spatial and temporal fluctuations are independent, leading to factorization of the spectral density into S s (k) = ∞ −∞ e −ikx C s (xn)dx and S t (ω) = ∞ −∞ e −iωt C t (t)dt. An example of a type of environ-ment that would be characterized by this kind of spectrum is a collection of uniformly distributed, statistically independent and similar sources of noise that affect the qubits according to coupling law dependent on the relative position µ(r j − x q n), where r j is the position of jth source. Then the noise field is given by ξ(r, t) = j µ(r j − r)η j (t), where the independence and similarity translates to statistical properties of the temporal fluctuations, η j (t)η j ′ (t ′ ) = δ j,j ′ C t (t − t ′ ). The auto-correlation function of such a noise field reads ξ(r, t)ξ(r ′ , t ′ ) = C t (t − t ′ ) j µ(r j − r)µ(r j − r ′ ), with the spatial component determined by the coupling law and the source distribution ρ scr (r ′′ ) = j δ(r ′′ − r j ), (47) where the uniformity of the distribution has been invoked, ρ scr (r ′′ + r) = ρ scr (r ′′ ). Furthermore, we choose for both components of spectral density identical line shapes, given by the Lorentzian distribution S 0 (z) = 2/(1 + z 2 ), Here ±k s and ±ω s are the characteristic wavenumber and frequency of the spectrum (the positions of spectral lines and their mirror images, as required by the symmetry of spectral density), and ν 2 s/t are the maximal intensities of each component of the spectrum. The correlation length and the correlation time of the noise field are denoted by x c and t c , respectively. This choice of spectral line shape is not necessarily representative of actual physical systems. Instead of realism, our approach was to provide an example that is both nontrivial (with the spectra being non-monotonous), and relatively effortless to follow.
Here, the slopes A(k p , ω p ) constituting the dataset matrices [A  (1) is greater than the revival period of Fourier coefficients (i.e., l c + 1 > l 0 ), the lth and (l 0 + l)th rows of the matrix are identical (or are very similar, when (N 0 + 1)/γ is not an integer). In that case, the matrix cannot be inverted because the determinant is zero (if rows are identical) or it is not well conditioned and an attempt at the numerical calculation of its reciprocal would be very inaccurate. To avoid such problems, the positions of qubits in the block used here were chosen at random from Gaussian distributions P (x q<N0 ) ∝ exp{−[x q − qL 0 /(N 0 + 1)] 2 /(2σ 2 )} with σ = 0.2L 0 /(N 0 + 1), and x N0 set in such a way that x 1 + x N0 = L 0 . The random shifts effectively "deregularize" the distribution, thus removing the revivals and periodicity. Instead, the Fourier series coefficients for large l of such an irregular distribution scatter around average value |v 0 | 2 /N 0 = N 0 /L 2 0 in a noisy fashion. We simulated the input matrices A (k0,ω0) for eighteen pairs of primary wavenumber and frequency settings combined from six choices of k 0 : {0.00, 0.05k d , 0.10k d , 0.15k d , 0.20k d , 0.25k d } and three of ω 0 : {0.15 × 2π/t c , 0.20 × 2π/t c , 0.25 × 2π/t c }. The cut-off indexes were chosen to be m c = 3 and L m = {−2, −1, 0, 1, 2}, so that each dataset matrix was 2 × 5 dimensional. This means that, in total, we had to gen- erate 10 × 18 = 180 unique slopes. Then we computed the matricesŨ −1 and (Ṽ (m) ) −1 (for m = 1, 3) according to Eqs. (43) and (44). Finally, we made use of the formula (42) to obtain the values of the reconstructed spectral density, which are presented in Fig. 6. In this example the reconstructed values of spectral density fit well the real course of S n (k, ω), for which we have chosen x c = 1.5L 0 , k s = 0.2k d and ω s = 0.2 × 2π/t c .

VI. CONCLUSIONS
We have demonstrated how a system of many qubits, arranged in a one-dimensional geometry, subjected to pure dephasing due to environmental noise, and controlled with coordinated π pulse sequences, can be used to reconstruct the spatiotemporal spectral density of the noise. The filtering of certain noise frequencies is achieved in a standard way-by using periodic pulse sequences. The key element of the scheme that turns the multi-qubit register into a spectrometer of spatial fluctuations, is the creation of a pattern of relative temporal shifts between π pulse sequences applied to various qubits. When these shifts depend linearly on inter-qubit distance, the slope of the dependence together with the filter frequency of the sequences define the filter wavelength (equivalently a wavevector).
Although the idea is rather simple and appealing, using such a spectrometer to obtain quantitative data on the spatiotemporal spectrum of the noise field requires more caution than when a single qubit is used to reconstruct the spectrum of temporal fluctuations of the local noise affecting it (note that the latter task, while routinely performed in recent years, is not entirely trivial either, especially in case of temporal spectra having peaks at finite frequencies [48,56]). Large part of the paper has been devoted to detailed explanation of methods allowing for reliable extraction of "spectroscopic" information (spectroscopic formulas defined in Sec. IV A) from raw measured data. The procedure of data analysis that we have discussed allows for self-diagnosis of the method, e.g., verifying if the noise is indeed stationary and spatially uniform, so that the spectrum can be in fact defined. It also allows for obtaining qualitative information about the noise-estimating correlation times and lengths, checking if the noise has a characteristic frequency or a wavevector-even before the spectral density is fully reconstructed. Further obstacles that have to be overcome are caused by the the fact that the spectroscopic information obtained from a single set of measurements involves multiple frequencies and wavevectors. This problem is known from single-qubit spectroscopy, in which the temporal spectrum at higher harmonics of the chosen filter frequency of the pulse sequence contributes to the measured signal. In the case considered here this issue becomes more significant due to the character of the wavevector filter generated in our scheme. We have given a thorough discussion of modifications of the single-qubit Alvarez-Suter method necessary to deal with the case of multiple qubits.
The spatiotemporal noise field spectrum reconstructed with methods discussed here has a potential to become a very useful tool for characterizing physical properties of the complex system probed by the register. The main objection that could be raised against this point is the concern regarding the objectivity of the noise field, i.e., to what degree the presence of the register and its inevitable mutual interactions with the environment influences the statistics of the field it is used to probe. If this influence could not be neglected, the spectra reconstructed using variously configured registers would turn out to be manifestly distinct; in other words, different observers would report seeing different noise fields. In such an event, the noise could not be ascribed to the environment alone, and it could not be used as a reliable source of inference of environment's properties. One example of a system in which such a situation might occur are gate-controlled quantum dots hosting spin qubits [57]. While the influence of charge noise on such qubits in well established [58][59][60], the origin of this noise is a subject of discussion [61,62]. Most importantly, it is not clear to what extent the noise is due to processes that depend on the spatial arrangement of metallic gates and voltages applied to them, which determine the positions of the quantum dots (and thus the qubits). In such systems the description based on local noises and their cross-correlations is more suitable, as one usually is more interested in characterizing the decoherence suffered by the given register, which remains in fixed configuration ultimately designed for a different purpose then being a probe of the environment. On the other hand, the existence of objective noise field seems to be generally physically well motivated for systems where the influence of the environment on qubits can be described in the terms of coupling to fluctuating external electric [63][64][65] and magnetic fields [22,[66][67][68][69]. A flagship example of such a system is a nitrogen vacancy center in a diamond acting as a sensor of fluctuating magnetic fields originating from outside of the diamond nanocrystal, e.g., generated by molecules [23,[70][71][72] and magnetic or superconducting materials [69] present in the vicinity of the piece of diamond in which the qubit is localized. What are the quantitative criteria for objectivity understood in this way remains an open question, obtaining the answer to which is one of the most vital issues for further studies on this subject. It is important to note that this question is not a novel one: It is equally applicable to the case of single-qubit spectroscopy, where the objectivity of temporal fluctuations affecting the qubit probe is also of concern. Even though there is no clear resolution to this problem, the success and practical importance of the dynamical-decoupling based spectroscopy method are undeniable. We believe that our multi-qubit extension will also prove to be viable despite these potential objections.
Aside from purely informative qualities of the noise field spectroscopy the reconstructed spectrum S n (k, ω) could be utilized as an input for calculation of decoherence of different multi-qubit system arrangements and geometries, development of realistic quantum error correction protocols that take into account the presence of correlated errors caused by cross-correlations in noises experienced by the qubits, and development of optimized dynamical decoupling protocols, giving maximal enhancement of multi-qubit coherence for given set of practical constraints (number of pulses, minimal time delays between them etc.) [5,[73][74][75].
where the Fourier series coefficients c mωp are given by (24). We shall now use this form to calculate the temporal part of the Fourier transform (19): where the passband filter function h is defined in (23) and the qubit spatial distribution is given by (8).
Since it is defined on a finite interval, the restricted qubit distribution can also be decomposed into Fourier series: Now assume that T = n t 2τ p = n t T 0 and L = n s L 0 where L 0 is the period of ρ(x). The spatial Fourier coefficients calculated in respect to this spectrometer length satisfy the following relation (here k d = 2π/L 0 , so 2π/L = k d /n s ): where v l ′ k d are identical to (25). Therefore, the spatial part of Fourier transform reads