Detecting initial correlations via correlated spectroscopy in hybrid quantum systems

Generic mesoscopic quantum systems that interact with their environment tend to display appreciable correlations with environment that often play an important role in the physical properties of the system. However, the experimental methods needed to characterize such systems either ignore the role of initial correlations or scale unfavourably with system dimensions. Here, we present a technique that is agnostic to system–environment correlations and can be potentially implemented experimentally. Under a specific set of constraints, we demonstrate the ability to detect and measure specific correlations. We apply the technique to two cases related to Nitrogen Vacancy Centers (NV). Firstly, we use the technique on an NV coupled to a P1 defect centre in the environment to demonstrate the ability to detect dark spins. Secondly, we implement the technique on a hybrid quantum system of NV coupled to an optical cavity with initial correlations. We extract the interaction strength and effective number of interacting NVs from the initial correlations using our technique.

Quantifying complex dynamics of a correlated system and environment is a necessary step for realizing practical quantum technologies. Traditionally, majority of the methods describing the dynamical evolution of an open-quantum system assume factorization approximation 1 , which states that the initial system-environmental correlation is negligible. The advantage in assuming absence of correlations between system and environment is the guarantee that the evolution of the reduced system is given by completely positive trace preserving maps (CPTP) between the reduced system states 2 . A typical dynamical equation that appears for describing complex quantum systems derived under factorization approximation alongside Markov approximation and rotating wave approximation is the GKLS (Gorini-Kossakowski-Lindblad-Sudarshan) master equation 3 . However, these descriptions fail in the presence of initial correlations, which are present in strongly correlated nanoscale system-environment dynamics [4][5][6] .
For such initially correlated system-environment dynamics, the reduced dynamical map describing the system evolution is described by "not completely positive" (NCP) maps. NCP-maps capture the effect of initial correlations on the subsequent dynamics, since in the presence of correlations, the marginal states and the correlation matrix are not independent. Such NCP-maps have been discussed in the literature from a mathematical point of view [7][8][9] , and several theoretical results have been established relating to the violation of laws of physics as a consequence of NCP-maps 10,11 . Process tensor tomography does not assume CP maps and is often used to characterize such maps. While state tomography of a d × d density matrix scales as O (d 2 ) , process tomography scales as O (d 4 ) and process tensor tomography scales as O (d 6 ) 12 , making it impractical to characterize even modest systems. Thus experimental efforts till date for characterization have been majorly capable of only witnessing the presence of NCP maps by detecting initial correlations. Most of these techniques [13][14][15] rely on measurement of the change in trace distance in time between two states which initially have same system marginal states. Another method 16 to witness initial correlation is by measuring the conditional past-future correlations.
The fact that characterizing NCP dynamics is cumbersome has direct consequences on several forms of spectroscopic techniques used. The experimental inconvenience has forced these techniques to be based on local master equations (like GKLS) that are based on CP maps formalism which in turn rely on Hilbert space dimensionality of the system and the Born-Markov assumption. Hence, it is essential to develop spectroscopic techniques that neither rely strongly on number of dimensions of the Hilbert space of the system nor on whether these degrees of freedom are interacting with an environment.
In this work, we propose a new method we call Prepare-Probe-Spectroscopy (PrePSy). Using PrePSy, we show selective detection of particular system-environment correlations and furthermore quantify them under some assumptions. Our spectroscopic method extracts information about the system-environmental initial correlations by utilizing the role measurements play on correlated systems. Such measurement based methods were previously used to characterize NCP maps 17 , modify information theoretic bounds 18 , characterize correlations in the quantum baths 19 and witness initial correlations 20 .
We demonstrate PrePSy by detecting the initial system-environment correlations between NV and a P1 defect centre to predict the presence of a dark spin in the environment. Here P1 is a dark spin which refers to defects in the diamond lattice that does not yield fluorescence under optical excitation. NV centres can be engineered to control these dark spins [21][22][23] and further employ these dark spin as a resource for quantum sensing 24 and quantum information 25 . The spectroscopic technique applied in PrePSy additionally facilitates deriving information related to system-environment dynamics embedded in the initial correlations. This functionality of PrePSy is shown using a spin-photon hybrid quantum system-environment represented by color centers present in diamonds 26 or SiC 27,28 placed in an optical resonator. By implementing PrePSy on this setup, information such as the photon-mediated interaction strength between two color centers and decoherence led effects are derived. Such hybrid QED setups have been used to implement quantum information processing applications 29,30 , quantum key distribution 31,32 and entanglement distribution [33][34][35] . Among various attempts for constructing such solid-state cavity QED systems, NVs coupled to a photonic crystal cavity 36,37 , microsphere resonator 38,39 , or micro-toroids 40,41 have emerged as a promising candidate. The prevalence of color centers for quantum information applications is because of their spin degree of freedom, which has advantages such as a long spin decoherence time at room temperature [42][43][44] and convenient optical readout of spin state 45 in addition to the photonic degree of freedom.
The organization of the manuscript is as follows-in the section "Prepare Probe Spectroscopy (PrePSy)", we describe the three steps-Prepare, Evolve, and Probe, that constitute PrePSy. In section "Toy model", the PrePSy technique is demonstrated using a toy model to detect and measure initial correlations. Section "Application: Dark spins coupled to NV" describes the results of applying PrePSy on an NV centre to detect the presence of a P1 defect centre in its environment. In section "Application: Cavity QED with NVs", the capability of PrePSy to decode information related to dynamics of system-environment from the initial correlations is displayed for multiple NVs positioned within a photonic cavity signifying a spin-photon hybrid system. Information such as the coupling constant and effective number of NVs interacting given the decoherence is extracted.

Prepare Probe Spectroscopy (PrePSy)
To study systems which are initially correlated with their environment, it is well known that operation of an entanglement breaking channel 11,17,46 on the system can be used to generate environmental states conditioned on the outcome of the system states. In addition to entanglement breaking channel, we use a two dimensional spectroscopic technique to probe dynamics of the system interacting with a conditioned state of the environment. We demonstrate how the two techniques in tandem can detect and measure initial or intermediate correlations in quantum systems alongside detecting the hidden Hilbert space dimensionality of the environment.
Step 1: Conditional preparation. The first step entails preparing the system and environment state. In this step, a projective measurement is performed on the system which acts as an entanglement breaking channel. This is followed by a unitary transformation of the resultant system marginal state to a standard state.
For a mathematical description of this step, consider a general initial state of the system-environment written as where R is the total state of the system-environment, ρ and τ are the marginal state of the system and environment respectively, and χ is the correlation matrix. The projective measurement E m projects the system state onto |m� , which leads to Here σ (m) = τ + Tr s (|m��m| ⊗ I · χ) is the post measurement state of the environment where Tr s (·) is the partial trace over the system dimensions. Thus in presence of initial system-environment correlations, projective measurement performed on the system generates environment state conditioned on the measurement. The choice of projective measurement operator is crucial as it selectively allows parts of the initial correlation matrix χ , that are along the projective operator as shown in section 1 of Supplementary. Next, the system state is prepared by a unitary transformation to a fiducial state |0� . The total state after preparation is as follows, where U m transforms system state |m� to |0� . Hence, after preparation any difference between two states which were projected onto different outcomes ( |m� , |n� ), is only between the environmental states ( σ (m) and σ (n) ). This difference between the environmental marginal states is caused only due to the presence of initial correlations. We perform step 1 multiple times with variations in the states chosen for projective measurement so that pairwise differencing as shown later in Eq. (4) can be applied. We choose orthogonal states for the set of variations in projective measurements as it facilitates generation of independent data sets. The number of maximum possible iterations scales favourably with the dimensions of the system only. www.nature.com/scientificreports/ Step 2: Two dimensional spectroscopy. The next step is inspired from 2D Phase coherent spectroscopy which is routinely deployed for nuclear magnetic resonance (NMR) 47,48 . An elementary example of 2D spectroscopy is shown in Fig. 1a. The rotation caused by the first control-pulse is for preparation of spins. The delay t 1 and t 2 is individually varied and separated by a similar rotation generating control pulse. Consequently the signal is recorded by measuring an appropriate property for e.g. magnetization is measured for nuclear spin ensemble. The signal will be a function of t 1 and t 2 , and on Fourier transformation will generate a 2D spectrum common to spectroscopy. Though one-dimensional spectroscopy suffices the need for characterizing initial correlations, we choose two-dimensional spectroscopy. The reason is that with two-dimensional spectroscopy in addition to characterizing correlations we can also predict system-environment dynamics based on the correlations measured. Thus the added benefit of 2D spectroscopy is measuring coherent population transfer channels i.e., valid transitions, in a multipartite system required for calculating the system-environment dynamics.
Step 1-"Conditional Preparation" is integrated with the two-dimensional spectroscopy by substituting the step 1 in place of the first pulse of the pulse sequence for two dimensional spectroscopy as shown in Fig. 1a. The substitution generates a procedure as shown in the first and second part of Fig. 1b. Thus after the conditional preparation, the system and environment evolve for time t 1 followed by a π/2 pulse generated by suitable rotation operator ( Â ). We use the NMR definition of π/2 pulse and loosely extend it to any system Hilbert space to denote a π/2 rotation in a plane defined by the operator Â acting on the prepared state. Here, π 2 pulse implies traversing half the angle between two states, preferably orthonormal to each other. The angle subtended between them is θ = 2 * cos −1 (|�ψ|γ �|) where |ψ� and |γ � are normalized states. This is then followed by free evolution for time t 2 .
Step 3: Measurement. The final step is probing/measuring the system with a suitably chosen system observable. This system observable is kept constant for all the runs of variations of the projective measurement performed in step 1. From the set of these variations, pairs are created for their respective measured signal to be differenced as follows where M i (t 1 , t 2 ) is the measured signal where in step 1 the projective measurement was performed along |i� . Presence of initial correlation leads to difference in the environmental marginal states after step 1. This type of difference is now also seen in the system marginal states due to spectroscopy performed in step 2. Thus pairwise differencing at the end implies that effectively we are observing the result of spectroscopy performed on only the terms representing the presence of initial correlations. Each of these pairs is then Fourier transformed to generate the spectrum. The maximum number of pairs possible is d C 2 , however all pairs are not required. Thus PrePSy scales as O (d 2 ) with the system dimensions. Figure 1b shows the complete procedure for Prepare Probe Spectroscopy (PrePSy).
For zero initial correlations, the total state after step 1 (conditional preparation) will be the same for any projective measurement outcome. Thus the total state will also be the same after step 2 (spectroscopy) for any projection operator chosen in step 1 (conditional preparation). This is because the same process is performed for all states. Hence for any two pairs of different projection operators chosen in step 1 (conditional preparation), the pairwise difference in step 3 (measurement) will give zero if there is no initial correlation. Thus if PrePSy detects some signal, it implies the presence of correlations. An appropriate choice of the projective measurement, rotation operator of the π/2 pulse and final measurement operators will ensure that the initial correlations are always almost captured by PrePSy as shown in sections 1 and 2 of Supplementary. This will reduce the need to measure all d C 2 pairs.
Though the scaling of PrePSy as O (d 2 ) is less than process tensor tomography, which scales as O (d 6 ) , this restricted scaling is sufficient to detect the system-environment correlations outside a set of pairs of projective measurement choices, which lead to zero detection as shown in Section 1 of Supplementary. This independence from the environmental dimension provides a huge benefit experimentally. Moreover if the entities in the environment interacting with the system increase then on an average the correlations with the system will mostly decrease. This effect is due to the fact that a small number of environment spins approximates a thermal bath under generic interactions. For typical Hamilonians the correlations always decrease with the increase in dimension of the environment generating thermalization 49 . This renders PrePSy and other system-environment characterizing tools irrelevant in this regime of large environment dimensions. www.nature.com/scientificreports/ PrePSy scales not only with the system dimensions but also with sampling being done in t 1 and t 2 . For a sampling size "S" for t 1 & t 2 and "A" number of averages per sample, there needs to be (A × S) 2 measurements. Hence the mean number of runs required such that all the data points for a particular choice of preparation is (1/P(a)) · (AS) 2 where P(a) is the probability of state "a". Thus the mean number of runs required for a particular choice of preparation depends directly on the sampling size and indirectly on the system dimensions (P(a) depends on system dimensions). Thus PrePSy scales exponentially with the sampling size.

Toy model
In this section, we devise a toy model of a system and environment with initial correlation and characterize it using PrePSy. A generic Hamiltonian for the system and the environment can be written as a sum of two parts. The first part is a diagonal Hamiltonian that defines the joint energy levels of the system and environment. The other part is an interaction Hamiltonian, which may be off-diagonal to promote transitions between different levels. We can write a generic density matrix (R) of such a system-environment in the Fano representation 50 as ..D 2 −1 } are elements of a Hermitian traceless matrix basis, d(D) is the Hilbert space dimension of the system(environment) and r, s and T ij are real coefficients such that R is a well-defined density operator. Using these guidelines, we construct our toy model.
Here we consider a toy model of a spin-qubit system coupled to a spin-qubit environment as an example and demonstrate the effectiveness of PrePSy by measuring initial correlations. Though our constructions are simple, the derived conclusions hold true for higher dimensions.

Total Hamiltonian and initial state. A general two-qubit Hamiltonian can be written as
i ) is the spin−1/2 angular momentum of the system (environment), ω (s) ( ω (e) ) is the energy difference between the two states of the system (environment) qubit and (ij) is the coupling between the two spins.
For the initial state, the general density matrix of two qubits can be written as where � σ = (σ x , σ y , σ z ) is a Pauli vector. For R to be a well defined density matrix, � u, � v ∈ R 3 and T j,k ∈ R are appropriately chosen. Using singular value decomposition, matrix T = {T jk } can be written as T = O a diag{c 2 , c 2 , c 3 }O b where O a and O b are orthogonal matrices. One can see that T has local unitary equivalence to where a, b, c are real coefficients chosen such that R is a valid density matrix. Thus a general two qubit state can always be reduced, up to local unitary equivalence, to a state in the above form, Since focus is on the correlations in the bipartite system, we set � a = � b = 0 , thus considering only the maximally mixed marginals. The system-environment state is given as Here, {c j } are real parameters consistent with R being a well-defined density operator. For the numerical example, arbitrarily but consistently chosen c j and {ω (s) , ω (e) , (ij) } are used. PrePSy requires minimum of two variations in the projection operator in Step 1 for pairwise differencing, and the number of orthogonal variations can go up to the size of the system. Here as the system is itself two dimensional, only two runs with projection operators |x��x| and |−x��−x| are performed. The π/2 pulse is chosen to be in a direction perpendicular to both the projective operators i.e. about the y-axis for an optimal rotation. The final measurement operator chosen is the population measurement of |x� . In section 2 of Supplementary, we present a calculation to understand the relation between the various matrix elements discussed in producing the signal detected. The simulation performed assumes a closed unitary dynamics of two qubit. PrePSy on toy model. For a maximally mixed initial state with non-zero correlation, applying PrePSy gives the peaks shown in Fig. 2a. For the simulation, we chose 150 sample point in t 1 and t 2 . Thus a total of 150 2 = 22,500 experiments were run for each outcome of the conditional preparation step. The probability of each outcome is 1/2 since the initial state is chosen to be a maximally mixed marginal state. Thus roughly for There is an extra 2 to account for the two variations required to form a pair. Hence we can see PrePSy scales as (1/P(a)) · (AS) 2 where A is the number of averages, S is the number of samples and P(a) is the probability of one of the outcomes. Thus obtaining non-zero data in the form of peaks implies that it does detect correlation correctly as discussed previously. The position of the peaks in the 2D diagram in Fig. 2a describes the frequency corresponding to the energy gap of two states between which a population transfer occurred during t 1 and t 2 . The peak positions depends on the interaction Hamiltonian which dictates the energy gap. The intensity of the peak represents the quantity of the population transferred and hence depends on the initial density matrix and thus the initial correlations. Upon measuring the variation of intensity of any peak with the correlation strength, a linear dependency is revealed. This is shown in section 2 of Supplementary, where Eq. (7) in Supplementary implies that derivative of the signal measured with initial correlation is independent of the initial correlation proving the aforementioned statement. Linearity is not surprising because the master equation is linear in the density matrix which contains the correlation. Plot of total signal intensity (2D sum of the signal over all frequency) versus correlation is shown in Fig. 2b.
If the equation of this line is known a priori, then by applying PrePSy and measuring the total intensity of the total signal generated by PrePSy, initial correlation can be calculated from the equation of the line. Thus to measure correlations the equation of line must be known.
To know the equation of a line two data points are required. Out of which first is trivial i.e. the line passes through (0, 0). One more data point, if acquired, should be sufficient. However, knowing the system-environment correlation a priori is nontrivial. If the Hamiltonian of the system-environment is entirely known then a simple method for calculating the other data point is using the Gibb's thermal state where initial correlation can be calculated. Performing PrePSy on the thermal state of the system and environment will give an additional data point. If the Hamiltonian is unknown, then the equation of the line defining the correlation for the measured signal cannot be estimated. In this case, for a small enough system-environment Hilbert space, brute force method can be used to learn Hamiltonian parameters from the PrePSy data.
As mentioned earlier one-dimensional spectroscopy suffices for characterizing initial correlations. The results of PrePSy, if 1D spectroscopy is performed on the toy model instead of 2D spectroscopy is shown in section 3 of Supplementary. Unlike 2D spectroscopy, for PrePSy with 1D spectroscopy, Hamiltonian parameters cannot be learned from the measured signal as the signal generated by a particular Hamiltonian is no longer unique. The loss of uniqueness is briefly discussed in the conclusion section.

Application: Dark spins coupled to NV
The environment of an NV centre is a spin-bath that generates fluctuating magnetic field leading to the decoherence of the NV centre spin. Various decoupling schemes presented 51,52 work to increase the coherence of the NV centre spin. A category of the spins in the environment is dark spins i.e., electronic and nuclear spins that cannot be initialized or detected optically 53 .
Here we consider substitutional nitrogen (P1) defect centres which are dark electron spins among a diverse set of dark spins. Coupling with proximal P1 centres is possible by tuning a static magnetic field to a level anticrossing 21 , but this technique fails if the dark defect is unknown. We present PrePSy as a technique to detect dark spins in the environment irrespective of its gyromagnetic ratio or hyperfine coupling.
The interaction between NV center and P1 defect is described by the magnetic dipole-dipole interaction and can be written as  Fig. 3. For PrePSy, the initial state was chosen to be a thermal state at room temperature of NV centre and P1 centre combined. The presence of a signal generated by PrePSy indicates the presence of dark spins in the vicinity of the NV centre. Further characterization of the defect by analyzing the PrePSy data is feasible by fitting Hamiltonian models. Few techniques exist that employ 2D spectroscopy similar to PrePSy to characterize nearby spins 54,55 . There are even techniques 56,57 that can detect and calculate the position and type of the dark spins by manipulation of the dark spins is not required by PrePSy.

Application: Cavity QED with NVs
Developing reliable interfaces between photons and the NVs 58 are particularly tricky due to a significantly small fraction of fluorescence contributed by zero phonon line and low coupling strength of NVs interacting with cavity photon. Moreover, these cavities have dimensions in the nanoscale regime which generates low yield and in addition to that positioning of quantum emitters is not deterministic 59 . Thus scaling to multiple quantum emitters coupled with photons scales the difficulty.
Additionally, the interactions of individual NV center with cavity photon is measured with the Purcell enhancement of the spontaneous emission 60 and hybridization of the electronic state 61 . However, for multiple emitters interacting with photons, although dipole-dipole interaction between the emitters is negligible, emitter-emitter interaction will exist through photon mediation. This will give rise to entanglement 62,63 . Thus for many quantum emitters interacting with cavity photons, characterization using previous methods becomes nontrivial. For such type of setups, we demonstrate that PrePSy can quickly provide valuable insights.
As most of the system-environment properties are encoded in the correlations, we present PrePSy as an auxiliary characterizing technique for such systems where coupling between entities is sufficiently stronger than collective decoherence effects. With PrePSy, we demonstrate for multiple NVs trapped within a cavity the ability to read off parameters like coupling constant, the quantity of NVs coupled to the cavity. We choose a single NV to be the system and the cavity, along with other NVs to be the environment, as shown in Fig. 4a. Hamiltonian dynamics. The NV has an S = 1 electronic ground state labeled as | 3 A 2 � = |E 0 � ⊗ |m s = 0, ±1� where |E 0 � is the orbital angular momentum state and |m s = 0, ±1� are the spin angular momentum states. The optical transition between the ground and excited state manifold is spin preserving but changes the orbital angular momentum. The spin-photon coupling can be reduced to an effective pairwise photon mediated Jaynes-Cummings model with the help of laser-induced Raman transition between two centers via the exchange of virtual cavity photons 64 .
To show this consider that the cavity mode in optical regime dispersively couples to the jth NV center's transition between the excited state ( |e� = |A 2 �:=1/ √ 2(|E − �|m s = 1� + |E + �|m s = −1�) and a states in the ground state manifold ( |0�:=|E 0 �|m s = −1� ), with coupling constant g j and detuning as shown in Fig. 4b. The selectivity of the states is because of the spin selectivity nature of the optical transition and this particular -type www.nature.com/scientificreports/ transition was recently used for spin-photon interaction 61 . Under the condition ≫ g j , the cavity operators can be adiabatically eliminated. To perform adiabatic elimination, time-averaging of the Hamiltonian is performed given that the cavity is detuned 65 . The corresponding Hamiltonian is reduced to In addition to cavity interaction, a largely detuned σ + laser is coupled to the same transition with Rabi frequency ′ j and detuning ′ ( ≫ ′ j ). The purpose of this laser transition is to eliminate the stark shift term of the state |0� generated by the vacuum state of the cavity.
Another adiabatic elimination of the excited state manifold of NVs is possible if a σ − laser is coupled to |e� ↔ |1� (where |1� = |E 0 �|m s = +1� and Rabi frequency is j ) to create a two photon Raman process Fig. 4C. Thus the effective Hamiltonian in the subspace of {|10�, (|e0� + |0e�)/ √ 2, (|e0� − |0e�)/ √ 2, |01�} is where ξ i,j = � i � j /� ij and � ij = g i g j /� . The adiabatic elimination is valid only when the effective NV-NV coupling is much larger than the laser coupling ( ij ≫ i , j ). The derivation of the effective Hamiltonian as shown in Eq. (13) which is twice adiabatically approximated is derived in section 4 of Supplementary.
If effective pairwise coupling is uniform i.e. ξ i,j = ξ ∀i, j then Hamiltonian for the group of two level systems can be written in terms of collective angular momentum operators ( J x , J y , J z ) as where N is the number of NVs in the cavity. The Hamiltonian derived is a type of non-linear rotator of spins and is called a one-axis twisting Hamiltonian. We simulate the density matrix evolution given in Eq. (15) where Ĥ is given by Eq. (12).
The Hilbert space spanned by ρ includes the system's (single NV) and the environment's (other NVs) Hilbert spaces. The system-environment is weakly coupled to a super environment (bath). The decoherence processes due to this bath considered for the Lindblad master equation are spontaneous emission and spin-lattice deexcitation. Here, σ αβ = |α i ��β i | , κ is the cavity decay rate, γ i 10 is the spin-lattice de-excitation rate and γ i ej is the spontaneous emission from the excited state of NV. www.nature.com/scientificreports/ PrePSy on cavity QED with NVs. The NV-cavity setup weakly coupled to a super environment is initialized to Gibb's state at room temperature where the initial correlation between system NV and other NVs is primarily dependent on the interaction strength and decoherence rate. For the NV-cavity system, the coupling strength of the interaction Hamiltonian can be few orders higher than the cavity decay rate 66,67 . Hence the thermal state will have significant correlations between the system NV and other NVs in the cavity generated by photon-mediated coupling. The choice of the initial thermal state also simplifies state preparation. For applying PrePSy, conditional preparation is performed with for two variations of the projection operator |+x��+x| and |−x��−x| defined in the Hilbert subspace of |0�, |1� and the system is then prepared to |+x� . The π/2 rotation is generated by σ y in the same Hilbert subspace. The final measurement is the population of the system NV in the state |+x� . Results for PrePSy's simulation for an open quantum system dynamics of 6 NV centers trapped in an optical cavity is shown in Fig. 5 below where a single NV is considered as the system and the rest of the NVs are defined as the environment.
The result of PrePSy displays peaks placed uniformly in the frequency domain. Similar to standard 2D coherent spectroscopy, the position of the peak corresponds to the gap between the energy levels. Consider Hamiltonian for six spins given by Eq. (14), where the energy separation between two eigen-states ( |j h , m h � , |j k , m k � ) in term of quantum numbers (j, m) is For all possible combination of (h, k) a peak will be observed in the 2D map. A peak at (x, y) corresponds to energy transition x during t 1 and energy transition y during t 2 respectively. For example the smallest energy gap according to Eq. (16) is between |0, 0� and |1, 1� which is 4ξ = 0.004 and the peak closest to the origin in the figure at (0.004, 0.004). Thus from PrePSy pairwise coupling for simple systems can be easily calculated. At the positions predicted by the Eq. (16), multiple smaller peaks are visible rather than a single peak because of the adiabatic approximation. The simulation is performed for Eq. (12), which does not assume adiabatic elimination of the excited state |e� ; however, the Eq. (16) does.
Since there is decoherence in the form of spontaneous emission and cavity decay, higher energy levels eventually are depopulated, thus all peaks will not appear in the image and the higher frequency peaks are generally less visible.
The total number of peaks visible corresponds total number of transitions possible. For 6 spins there are 10 distinct eigenenergy levels and 12 distinct energies of the transition. However Fig. 5 displays approximately 6 distinct energies which indicates the presence of decoherence and is equivalent to 4 spins interacting without decoherence. From this point of view, in the presence of noisy environment 6 NVs in the cavity are effectively just 4 NVs in the cavity without any decoherence.

Conclusions and discussion
PrePSy is capable of witnessing specific correlations, as demonstrated with the help of a toy model in Section "Toy Model". In addition to this, we show that PrePSy can also measure specific correlations because of the linear dependence of initial correlations on the signal measured. In Section "Application: Dark spins coupled to NV", by predicting dark spins in the environment, we demonstrate the capability of PrePSy to detect correlation of NV with electron spins near it. In section "Application: Cavity QED with NVs", we illustrate the potential of PrePSy to derive information related to the system-environment dynamics from the initial correlation using a cavity (16) E h,k = 4ξ j h − j k j h + j k + 1 + m 2 h − m 2 k Figure 5. PrePSy applied on a single NV placed in a cavity with an environment of five other NVs. Parameters are chosen as = 10g and = 0.01g , ξ = 0.001g. www.nature.com/scientificreports/ QED setup. More specifically, PrePSy measures the coupling constant of a spin-photon interaction between NVs and photons of an optical cavity and estimates the effective number of NVs interacting with the cavity in the presence of decoherence. Experimentally deploying PrePSy is feasible for a vast array of practically relevant physical systems since PrePSy assumes little about the system-environment state and is motivated for generic quantum systems. For color center-cavity QED systems, PrePSy is achievable with current experimental techniques 36,68,69 where experimentally spin-photon coupling has been demonstrated. The only requirement for applying PrePSy is a clear demarcation of system and environment. The system should be such that it is amenable to measurement, projection and rotation. Moreover, pulse sequences applied to the system should not affect the environment beyond producing conditional post-measurement states. This demarcation can restrict the physical systems PrePSy is applicable to. An example in NMR is the inability to separate out a single molecule as the system and the rest as environment because it works on ensemble statistics.

Scientific Reports
When experimentally performing PrePSy, 2D spectroscopy is preferred as it aids in quantifying the type of interaction Hamiltonian. As population transfer channels vary with different interaction Hamiltonian, so will the 2D image obtained. However, for 1D spectroscopy, two different Hamiltonians can give results that are not differentiable. For example, if 1D spectroscopy is applied to 3 NVs with pairwise interactions and a closed chain of 4 NVs, the same number of peaks will be visible for both the cases. Hence 2D spectroscopy differentiates between these two cases. The result of PrePSy for both these cases, as shown in Fig. 6, is strikingly different.
As our method scales favorably with the system size, it can be conveniently applied on NVs with any other form of ancilla [70][71][72] . We hope this will find applications in characterizing complex quantum systems that have initial or intermediate correlations.
where |a 5 � = |a 0 � for a = 0, 1 . The number of major peaks visible with 1D spectroscopy is 6, shown above the 2D image (ignoring the low frequency peaks). In 2D spectroscopy approximately 7 major peaks in the 2nd and 4th quadrant each. (b) PrePSy applied on a single NV which is in the cavity with two other NVs interacting in a pairwise manner. The Hamiltonian is H eff = 3 i,j;i� =j |1 i 0 j ��0 i 1 j | + H.C. . The number of major peaks visible with 1D spectroscopy is 6, shown above the 2D image (ignoring the low frequency peaks). In 2D spectroscopy approximately 9 major peaks in all quadrants.