Two-dimensional spectroscopy of Rydberg gases

Two-dimensional (2D) spectroscopy uses multiple electromagnetic pulses to infer the properties of a complex system. A paradigmatic class of target systems are molecular aggregates, for which one can obtain information on the eigenstates, various types of static and dynamic disorder and on relaxation processes. However, two-dimensional spectra can be difficult to interpret without precise knowledge of how the signal components relate to microscopic Hamiltonian parameters and system-bath interactions. Here we show that two-dimensional spectroscopy can be mapped in the microwave domain to highly controllable Rydberg quantum simulators. By porting 2D spectroscopy to Rydberg atoms, we firstly open the possibility of its experimental quantum simulation, in a case where parameters and interactions are very well known. Secondly, the technique may provide additional handles for experimental access to coherences between system states and the ability to discriminate different types of decoherence mechanisms in Rydberg gases. We investigate the requirements for a specific implementation utilizing multiple phase coherent microwave pulses and a phase cycling technique to isolate signal components.


Introduction
Assemblies of molecules can exchange electronic excitation energy via long-range interaction of transition dipoles [1]. This process leads to eigenstates that are coherently delocalised over several molecules, where the details depend on the molecular arrangement and the interaction with internal vibrations and the environment. Prominent examples are J-and H-aggregates of organic dyes [2] or the light harvesting systems in photosynthesis [3].
To understand the resultant eigenstates and dynamics in molecular systems, multidimensional optical spectroscopy has become an indispensable tool [4,5,6,7,8,9,10]. A sequence of ultrashort laser pulses with well defined relative phases results in a signal that is measured as a function of the time-delays between the pulses. Typically, a Fourier-transform with respect to some (or all) time delays yields a spectrum in the frequency domain. The most common technique, referred to as two-dimensional (2D) spectroscopy, employs four pulses and thus three time delays. Usually, the Fourier transformation is made with respect to the first time interval τ (between the first and second pulse) and the last time interval t (between the third and fourth pulse). One obtains a separate 2D spectrum for each time delay between the second and third pulse, called the waiting time T . This sequence of spectra holds a multitude of information about the dynamics of the system, its coherences and relaxation pathways as well as different line broadening mechanisms, see e.g. Ref. [10].
However, these spectra are difficult to interpret and one easily can draw flawed conclusions based on incomplete or inappropriate models of the system under study. Numerical calculations should on the one hand reproduce the experimental data, on the other hand they should provide insight into how the signal depends on the microscopic properties of the system, e.g. by systematically changing the model parameters. Unfortunately due to the theoretical complexity this is not possible for many systems of interest without making severe approximations, especially regarding the treatment of internal vibrations and the coupling to an environment. Although sophisticated numerical methods have been developed to treat the dynamics of molecular systems (e.g. Refs. [11,12,13,14]) in practice one can still handle only a small regime of parameters and model Hamiltonians.
An alternative to numerical simulations might be offered by quantum simulators, well controllable physical systems sharing a common Hamiltonian with the target system that can be used to study quantum dynamics in regimes where numerical simulations fail. There have been several proposals of quantum simulators for assemblies of interacting molecules [15,16,17,18,19] , which could also benefit from sophisticated probing techniques such as multidimensional spectroscopy. Enabling the latter on these platforms, would finally facilitate direct comparisons with real molecular systems. In the present work we propose a quantum simulator based on dipole-coupled assemblies of Rydberg atoms that closely mimics the physics of small molecular aggregates and is compatible with multidimensional spectroscopy techniques.
In the case of molecular assemblies, the transport of electronic excitation energy proceeds on ultrafast timescales by optical dipole couplings between molecular subunits on sub-wavelength distance scales [20,10,2]. In a Rydberg quantum simulator, analogous dipole-dipole interactions in the microwave regime allow the migration of Rydberg excitations through an assembly of atoms on microsecond and micrometer scales [21]. In contrast to the molecular case, the monomers (atoms) in the quantum simulator can be positioned at will [22,23], their transition energies and interactions tuned through choice of Rydberg quantum state [24] and Markovian [25,15] or non-Markovian [26,27,28] environments can be engineered to simulate coupling of an excitation to vibrations. Because transitions between the relevant Rydberg states correspond to microwave frequencies, it is not a-priori clear, whether parameters such as interaction strength, lifetime, dephasing rates and interrogation pulse durations and spacings scale well from the molecular setting to the Rydberg setting such that the technique remains feasible. We show in the following that it is indeed possible to find a suitable parameter regime. An important difference though, is the unavailability of phase matching in this type of microwave spectroscopy. We propose to use phase cycling [29,30,31] and explore the limitations of this requirement for Rydberg experiments.
In this way we establish a clean platform for testing basic features of multidimensional spectroscopy, in which a large variety of disorder and decoherence sources can be engineered [15,25,26,27,28]. Our main focus in this article is to adapt 2D spectroscopy to ultracold Rydberg systems in order to understand and simulate the corresponding situation in optical spectroscopy of molecular systems. In addition 2D spectroscopy would provide an additional diagnostics tool that may complement the already extensive arsenal available to ultracold Rydberg physics [32,33].
While we focus on the specific case of Rydberg atoms, microwave domain 2D spectroscopy can more generally be useful whenever quantum dynamics involves energy differences in the microwave realm, such as in quantum dots [34], Nitrogen-vacancycenters or superconducting circuit arrays [35,16,17], and much of this section should apply to these systems too. The technique is already being used, for example, in rotational spectroscopy of poly-atomic molecules [36] and collisional population transfer [37].
This article is organized as follows: In section 2, we present the general method of 2D microwave spectroscopy involving a series of short microwave pulses interacting with a set of coupled two-level systems. We discuss the pulse sequence, spectroscopic signal and phase cycling method for a measurement of the nonlinear response of the system briefly in that section, with further details in Appendix B, Appendix C and Appendix D. In section 3 we discuss how to port microwave 2D spectroscopy to the field of ultracold Rydberg gases. We show how decoherence in the system can be experimentally controlled in these systems and derive effective equations for its description in Appendix A. In section 4 we present proof-of-principle numerical calculations of 2D spectra for simple models like Rydberg dimer and Rydberg trimer aggregates, which are embedded within an environment of perturbing atoms giving rise to both homogenous and Sketch of Rydberg-dimer (blue) or trimer (blue+grey), to be interrogated by the microwave pulses, while green ultracold ground state atoms provide tuneable disorder and decoherence. (c) For the specific example of a Rydberg dimer, four twobody states play a role. Variables are discussed in the text. (d) The corresponding energy differences show up in a 2D spectrum as diagonal peaks, with off-diagonal peaks indicating coherences between them. Inhomogeneous broadening (∼ ∆E) and homogeneous broadening (∼ 1/Γ) can be distinguished as indicated, since they affect different frequency axes in the 2D spectrum.
inhomogeneous broadening mechanisms. We show that experimental constraints arising from the total spectrum acquisition time might require special consideration, possibly necessitating parallel interrogation of many identical Rydberg aggregates. Finally, we discuss future prospects for Rydberg quantum simulators and the experimental simulation of photosynthetic complexes.

Generic aspects of microwave 2D spectroscopy
In this section, we briefly outline the basic scenario that we propose and provide some general principles of microwave 2D spectroscopy. We keep this discussion quite general. Details specific to Rydberg quantum simulators are discussed in the subsequent section 3.

The system and its environment
We consider a collection of N (long-range) interacting two-level systems. We are in particular interested in a situation where the energy scale of the interaction is small compared to the transition energies between the two states of the two-level systems. An incoming electromagnetic field can drive transitions between these states. For weak fields it is convenient to classify the eigenstates of such a system according to the number of 'excitations', i.e, the number of excited two level systems. We will be mainly be interested in the subspaces with zero, one or two excitations, which are usually denoted as the ground state, the single exciton manifold and the twoexciton manifold, respectively. Typically, there are two effects that complicate the interpretation of measurements on this system: 1.) The measurement is performed on an ensemble of non-identical systems (static disorder). Such an ensemble emerges for example via imperfect preparation, sample inhomogeneity or fluctuations that are slow compared to the interrogation time scale. 2) Coupling to additional degrees of freedom that leads to dephasing and relaxation. Both effects usually result in broadening of absorption lineshapes. One strength of 2D spectroscopy is the ability to disentangle the contributions of these two effects.
Because of the coupling to the environment, the dynamics of the system is typically described by a time-evolution equation for the reduced density matrix of the system ρ S (t), which can be formally written as where U(t) is a time-evolution super operator.

Pulse sequence
This system is irradiated by a sequence of electromagnetic pulses, which are characterized by a time-dependent electric field E(t). The electric field E(t) is assumed to be given by a train of four short pulses where the A(t) are the (slowly varying) pulse-shapes, roughly of duration δt, and ϕ j control the relative phases of pulses, as sketched in figure 1(a), and E 0 their amplitude and polarisation vector. We define time delays as in figure 1, such that τ = t 2 − t 1 , The carrier frequency ω 0 is chosen to be close to the relevant transition energies h n of the two-level systems. For a system with spatial extent much smaller than the wavelength, the interaction of the system with impinging radiation in the dipole approximation takes the form whereμ (n) is the relevant dipole operator of object n and we have assumed a linearly polarized field aligned parallel with the transition dipole moments, using E(t) = |E(t)|. Furthermore, we assume that the wavelength is large compared to the spatial extent of the system. All discussed assumptions will be fulfilled for the Rydberg setup presented in the following sections. The matrix elements of the total system dipole operator µ ≡ nμ (n) are used to determine whether or not transitions between eigenstates are allowed.
The duration of each pulse δt and the delays τ , t are dictated by the details of the spectrum ofĤ sys to be interrogated, since pulse durations and delays have to be chosen to resolve the interactions of interest. We thus require where E and E belong to the same exciton manifold, so that microwave pulses are spectrally wide enough to reach all transitions in the system, see figure 1(c), and π/τ , π/t σ ω , where σ ω is the linewidth of a certain transition, in order for the spectrum to resolve details of resonance peaks.

Spectroscopic signal and phase-cycling
In optical multi-dimensional spectroscopy, the signal of interest generated by the evolution ofρ s (t) can frequently be isolated by choice of a specific geometry of incoming beam directions and outgoing signal direction, exploiting the phase-matching enforced by the propagation of fields through the interrogated medium.
For many of the quantum systems amenable to microwave interrogation that we listed in section 1, there will be no significant propagation through a medium and thus phase-matching is not readily available. On the other hand, these systems often allow additional access to some observables within the system, for example the determination of the total spin projection. Let in the followingF denote any observable, measured immediately after the last pulse, at time t end . The signal depends crucially on the time delays τ , T , t and the three relative phases between subsequent pulses. We write where ξ is an abbreviation for the triplet of relative phases. For weak pulses S ξ (τ, T, t ) is typically dominated by the linear response of the system, while 2D spectroscopy relies on the consequently weaker nonlinear response. Therefore one requires a technique to extract specific higher order contributions. One particular technique is phase cycling. In this approach the signal S ξ (τ, T, t ) is measured for several specific choices of the relative phases ξ. By choosing suitable phase combinations one can isolate a specific signal of interest by a summation of individual signals for the pulses where the c ξ coefficients are determined by the choices of ξ. In optical 2D spectroscopy there has been much interest in the so-called photon-echo signal [4]. Throughout this article, we will use this type of signal as our example to illustrate microwave Rydberg 2D spectroscopy. We briefly introduce phase-cycling in Appendix C, and discuss how the phases and the coefficients c ξ have to be chosen in this case.
As the final step in 2D spectroscopy, we Fourier transform the time-domain signal with respect to τ and t to reach the frequency domain 2D spectrum 3. Rydberg aggregate 2D spectroscopy We now move to the specific application of the 2D spectroscopy protocol to Rydberg aggregates.

Rydberg Aggregate
As a Rydberg aggregate [38] we consider a collection of dipole-dipole coupled atoms, where each atom can be in either of two electronic states | s or | p , with energies E s and E p , respectively. We choose | s = | ν = 43, l = 0, m = 0 and | p = | ν = 43, l = 1, m = 0 , where ν is the principal quantum number and l the angular momentum. Many other pairs of dipole coupled Rydberg states would be equally useful. The Hamiltonian for the aggregate including dipole-dipole interactions is then ( = 1): whereσ (n) ab = [| a b |] n is a transition operator on atom n. In general, dipole dipole interactions would allow transition to additional states with m = 0 [39], which can be suppressed by Zeeman-shifting them out of resonance with a magnetic field, see e.g. [40,41]. Then, describes the strength of dipole-dipole coupling between atom pairs located at r n , r m , with angle θ between the direction of transition dipoles and the inter-atomic axis. The transition dipole between the two states defined above points along the quantisation axis, given by the direction of the applied static magnetic field. Let us now consider additional microwave irradiation of that aggregate, assuming the microwave field to be linearly polarized along the quantisation axis. The central frequency ω 0 of the microwave pulses is chosen to be roughly resonant with the transition energy between these two states ω 0 = ω sp = E p − E s . For the example above ω 0 /(2π) = 49 GHz. Usingμ (n) = µ 0 (σ (n) sp +σ (n) ps )/2 in (3), the Hamiltonian for lightmatter interaction is, after conversion to a rotating frame at the field-frequency and performing the rotating wave approximation: Here, µ 0 is the transition dipole moment between the two states andẼ(t) is the complex electric field envelope of the microwave pulses, By ∆ mw = ω mw − ω 0 we denote the microwave detuning. Spatially, all atoms are assumed to be trapped in the quantum ground-state of very tight optical traps. Thus interaction parameters in (8) are effectively spatially averaged over the ground-state wavefunctions. Recent experiments can position these optical traps almost at will, using digital micromirror devices or spatial light modulators [22,23].
We assume the Rydberg aggregate is initialized in the stateρ s = | s s |, where | s = | ss..s , with all aggregate atoms in | s . It then interacts with the microwavepulse train (2), which will drive transitions to states which contain a non-zero p state population.
Most ultracold Rydberg experiments can routinely measure the total number of atoms in a specific Rydberg state quite accurately. We thus choseF = n [|p p|] n as measurement operator, which counts the number of p-excitations generated in the system by the microwave pulse train.

Controllable environment for decoherence
To assess the effects of decoherence we embed the Rydberg aggregates in a cold atom environment that allows tailoring a variety of different decoherence mechanisms [15,25,26,27].
For this we consider a few additional environment atoms, separately trapped at specific locations x n in the vicinity of the aggregate atoms. The environment atoms are then optically coupled from their ground state to two other states in the scheme of electromagnetically induced transparency (EIT) involving an auxiliary Rydberg state | r , shown in figure 2(b). By populating the | r state the | r → | s, p interactions introduce on-site static energy disorder, which implies that the energy of a certain system state then depends on which atoms are in the | s or | p state, rather than just their total populations. Additionally, when environment atoms reach the decaying state | e , the resultant light absorption allows the discrimination of aggregate states | s , | p and hence causes measurement induced dephasing. For a detailed discussion of these effects, we refer to [15]. The complete many-body Hamiltonian for the assembly of Rydberg aggregate plus environment atoms is described in (A.1)-(A.4). Calculating the corresponding time evolution of this open quantum system while explicitly including the non-trivial environment degrees of freedom is a formidable task. However, for the selected parameter regimes that we focus on below, it is possible to simplify the numerical treatment considerably. These are the cases in which it is possible to adiabatically eliminate the dynamics of environment atoms using techniques of [42] as in [25,15], which requires relaxation that is fast compared to the timescale for dipoledipole interactions, Γ p max nm [W nm ]. In that regime one can obtain the effective master equatioṅρ where the super-operator LÔ[ρ] =ÔρÔ † −(Ô †Ôρ +ρÔ †Ô )/2 accounts for the decoherence in the system. The density matrix is written asρ = kl ρ kl | χ k χ l |, where | ξ k is a many body state for the Rydberg aggregate only, given in the form of a tensor product where T k i ∈ {s, p} and k ∈ {1, 2, ...2 N } as we have 2 N states for N aggregate atoms. In (12),Ĥ agg ,Ĥ mw ,Ĥ eff andL (α) eff are the aggregate Hamiltonian (8), microwave Hamiltonian (10), effective Hamiltonian and Lindblad jump operator, respectively. The index α numbers the environment atoms used to control decoherence. We can writê with details given in Appendix A. Both,Ĥ eff andL (α) eff depend on the probe and coupling Rabi frequencies in the EIT scheme implemented for the environment atoms, Ω p and Ω c respectively, the corresponding detunings ∆ p and ∆ c , as well as the spontaneous decay rate of the intermediate state Γ p . Crucially, they also depend on the van-der-Waals interactions between environmental atoms and aggregate atoms. All of these can be varied in experiments, rendering both disorder and dephasing distributions widely tunable.
Later we shall engineer disorder and dephasing in this setup through a synthetically controlled distribution p x of locations of environmental atoms. These can for example be chosen from a Gaussian random distribution, with a fresh disorder realization in each experiment. Through the resultant distribution of H (k) eff for varying locations of environment atoms, the energy for a given aggregate state k is fluctuating. We call the distribution of these energies P E with width ∆E. Similarly, the L (k,α) eff define a distribution of dephasing rates P γ with width ∆γ for later reference.

Two-dimensional Rydberg spectra
Here we show that with Rydberg atoms and microwave pulses one can simulate the basic features of 2D optical spectroscopy of molecular aggregates. For the example of small aggregates and a specific choice of environmental parameters, we will discuss relevant aspects of the obtained spectra such as peak broadening caused by dephasing and static disorder. We will then discuss aspects specific to an implementation in Rydberg experiments.

Exemplary parameters
For the following demonstrations, we use the parameters summarized in table 1. The different parameter sets are chosen on one hand to demonstrate various aspects of 2D spectroscopy and on the other hand to allow tractable numerical and analytical analysis: the numerically simpler equation (12) is valid, and even diagonal when expressed in the basis of the eigenstates of the effective Hamiltonian, as discussed in Appendix D. The former aids numerical solutions, since the more general equation (1) would be very timeconsuming to solve, and latter allows simple analytical results. However, importantly, 2D spectroscopy of embedded Rydberg aggregates can also be performed for a broad range of parameters well outside the range of validity of (12) that are known to lead to more interesting dynamics [26,27], which is however challenging for numerical or analytical considerations.
For the parameters in the table, the order of magnitude of the most energetic eigenstate will be |U max | ∼max nm [W nm ], with U max ≈ 5 MHz. According to the discussion in section 2.1 we thus need pulse-lengths δt 200ns, which are technically feasible. Additionally, since we do not include the decay of the Rydberg excitation, the longest possible pulse sequence of duration τ tot , must be within the Rydberg system lifetime τ eff ≈ τ 0 /N = 14.6 µs , where τ 0 = 44 µs is the lifetime of a single atom in s-state ‡ and N = 3.

Rydberg dimer spectra
We first consider the simplest example: Two Rydberg atoms separated by a distance d = 10 µm form the aggregate, shown as blue balls in the inset of figure 3a. This aggregate has two single-exciton states | ± = (| ps ± | sp )/ √ 2 (with energies ω sp ± W (d)) and one two-exciton state | pp (with energy 2ω sp ). Each aggregate atom is flanked by a detector atom placed δ = 1.5 µm away, sketched as green balls.  Using phase cycling as described in section 2.3 with details listed in Appendix C, we show the absolute value of the expected 2D spectrum and its complex phase in figure 3. Here and in the following we plot the absolute value using arcsinh(5 · 10 3 |S|) to make small features more visible. We have obtained the spectrum shown in figure 3 by numerically simulating (12) for 64 different time delays in τ and t each. Simulations employ the high-level package XMDS [43,44]. For the pulses we used smoothened step functions with rise-time t rise δt §. We also show in Appendix C how the spectrum looks for a specific single pulse sequence. In these calculations, we have made sure to stay in regime of very low p-excitation probability and thus low signal strength, such that perturbation theory remains valid. In practice, one would want to venture to the limit of the perturbative regime with higher probability of p-excitation, in order to avoid the need for too many repetitions of the experiment, see section 4.5.
For this simple case of a symmetric dimer, we can easily understand the position, intensity and lineshape of peaks in the spectrum from perturbation theory [4]. One can identify and represent the relevant contributions to the signal by so-called double-sided Feynman diagrams, shown in figure 3 and discussed in Appendix B and Appendix D. Of crucial importance are the energy differences and transition dipoles between eigenstates of H agg + H eff , see (12), which are listed in table D1. Using the notation of figure 1(c), we find that from the initial ground-state | ss transitions only proceed via the states | ss ↔ | + ↔ | pp . The corresponding frequencies are ω +,ss = (2π) × 1.56 MHz and ω +,pp = (2π)×1.67 MHz and transition strengths µ +,ss = µ +,pp = √ 2µ. The state | − is not involved due to its anti-symmetry. The first interaction creates a coherence between the ground state | ss and | + . For the parameters used, the environment does not induce coupling between different elements of the reduced density matrix of the system. Therefore, during the time interval τ , the system accumulates a phase ω +,ss τ . The Fourier transformation with respect to τ gives a peak at the position ω τ = ω +,ss . The width of this peak is determined by the corresponding dephasing rate Γ +,ss . In order to provide a contribution to the signal, which is the number of p-excitations, the second pulse must bring the system either to the population of the | + state (diagrams c-e) contributing a single p, or to the doubly p-excited population (diagram f). For the present choice of parameters the environment does not affect populations. For diagram (c) and (d), the system is again in the coherence between ground state | ss and | + during the time-evolution t , after the third pulse. As before, the Fourier transform then gives a peak at ω t = ω +,ss . Thus diagrams (c) and (d) give rise to the diagonal peak at (2π) × (1.56, 1.56) MHz. Its intensity is proportional to |µ +,ss | 4 = 4µ 4 . In contrast, for diagrams (e) and (f), the third pulse brings the system into a coherence between the doubly excited state | pp and | + , during the time-evolution t . Since the frequency ω pp,+ is different from the frequency ω +,ss , see figure 1(c), this results (after Fourier transformation) in an off-diagonal peak at (2π) × (1.56, −1.67) MHz. The complete spectrum for the symmetric dimer can to a good approximation be written analytically as (ω τ − Ω +,ss )(ω t − Ω +,ss ) § The specific pulse shape we use is , with t eff = δt 2 − 3t rise , t rise = 0.01µs and δt = 0.05µs.
with Ω a,b = ω a,b − iΓ a,b . For details see Appendix B-Appendix D.1. Here, in particular, infinitely short pulses have been assumed. For the present choice of parameters and a single environment atom near each site, we find Γ +,ss ≈ (2π) × 0.2155 MHz and Γ pp,+ ≈ (2π) × 0.2150 MHz. There is very good agreement between the full nonperturbative calculation and the analytic result, as can be seen in figure D1. By changing the parameters of lasers acting on the environment atoms, our system allows tuning the widths Γ +,ss and Γ pp,+ . Note, that in parameter regimes where the simple master equation (12) is no longer valid the environment can cause more complicated effects than simple dephasing. If one introduces asymmetry to the system, for example by using a different placement of environment atoms or a different number of them near each of the two aggregate atoms, the | − state is no longer perfectly antisymmetric and additional weak contributions from this transition arise.

Rydberg trimer spectra
The Rydberg trimer is a significant extension compared to the dimer, as it now allows different geometries ranging from ones that lead to only one bright state (as in the previous dimer case) to cases where all transitions are dipole-allowed. In the following we show an example for each, but focus on the latter. As in the previous section for the dimer, we first also take one environment atom per site. These atoms are arranged as shown in the insets of figure 4.
In the first case, of panel (a), the trimer is arranged linearly, but slightly nonequidistant. This gives rise to several dipole-allowed transitions, which can be seen by the emergence of new peaks. These peaks are however very small compared to those involving the dominant symmetric state, which behaves akin to the dimer case in section 4.2.
To obtain a richer spectrum we focus on the L-shaped arrangement shown in the inset of figure 4b. Here, because of the anisotropy of the dipole-dipole interactions, the eigenstates have a very different form compared to the linear arrangement, which results in the emergence of several transitions that have comparably large strength. The spectrum in figure 4(b) thus shows a multitude of diagonal and off-diagonal peaks. As for the dimer case their positions, intensities and lineshapes can be understood from perturbation theory. The peak positions are again determined by the eigenenergy differences of the effective Hamiltonian in (B.1). We will comment here only on the main features. The trimer has three single-exciton states and three double-exciton states. This gives rise to three transitions between ground state | sss and the single excited states and 9 transitions between single-exciton states and double exciton states. The corresponding frequencies are given in table D2. Similar to the discussion in the previous section, the relevant contributions arising during the first time interval τ stem from a

Separating homogeneous and inhomogenous broadening
A strength of higher dimensional spectroscopy is the disentangling of different types of line broadening mechanisms. This can be clearly illustrated and studied in the embedded Rydberg setup, where one can independently tune the two types of broadening  Table 1. Each aggregate atom is flanked by two environment atoms, the trapping location of which is synthetically disordered independently, with width σ = 565 nm along the indicated directions, see text. We average 1000 realizations in each case. (c) Resultant histogram for the distribution of energy disorder P E and (d) dephasing rates P γ , see section 3.2, for the same parameters as in (a). (e,f), the same for parameters of panel (b). [15,25,26,27].
As stated earlier, we assume all atoms to be tightly trapped in the harmonic oscillator ground-state of their respective trap. Detector atoms should be easily trappable against the weak, dressed van-der-Waals Forces exerted on them, while aggregate atoms might have to undergo additional ground-state cooling after each phase set, or sequence, since they experience much stronger bare dipole-dipole interactions. If we assume ground-state trapping, one can envisage an effective Hamiltonian, obtained after integrating out spatial degrees of freedom. This would not contain any disorder without experimental intervention. However disorder can be introduced and controlled synthetically, by moving selected optical trap centres after measuring one complete spectrum, which we assume to happen in the following.
We then explore averaged spectra for a trimer as in the previous section, but with parameters adjusted to enhance diagonal disorder, arising viaĤ eff in (14) compared to dephasing, arising from (15). Synthetical disorder is included by randomly varying the trapping location of two environment atoms per aggregate atom, after each spectrum acquisition. By "spectrum acquisition" we imply gathering the complete dataset: measuring the mean signal strength (p-population) for all sets of phases and all required time-delays. The detector traps in realisation number k are then at x k = x 0 + ση k n, where the central positions x 0 are described in the caption of figure 5, n are normal unit vectors along the axes indicated by black arrows in the insets, η k is drawn from a Gaussian distribution of unit variance, and σ thus sets the width (standard deviation) of this position distribution.
The resultant spectra are shown in figure 5(a,b), together with the corresponding dephasing and disorder distributions. The shift, width and shape of 2D peaks along the diagonal reflects the on-site energy disorder distribution in figure 5(c,d), defined as P E in section 3.2. In contrast, the width in the off-diagonal direction matches the dephasing distributions in figure 5(e,f).

Spectrum acquisition time and ensemble spectroscopy
To experimentally measure any of the spectra shown, one has to gather enough statistics to obtain a sufficiently high signal to noise ratio. Since the signal is of fourth order in the weak electric field, the probability to promote atoms to the p-state is quite small. Now, phase cycling relies on a precise cancellation of an undesired large leading order signal in the sum (6). Hence, we determine the signal to noise ratio required for this cancellation to work. We explicitly explored this aspect in the simulations shown in figure 6 and discuss it in more detail in Appendix E. For the parameters of figure 3, we infer the need for about N rep ≈ 10 5 repetitions of each measurement, where the latter refers to a single sequence of pulses with subsequent measurement of the total p-population.
Based on this one can estimate the duration of acquisition of a single spectrum. The maximum interval between pulses, τ max is linked via τ max = π/∆ω with the desired frequency resolution ∆ω. In order to resolve the shape of peaks, ∆ω ought be smaller than typical dephasing rates γ in the system, which govern the width of peaks. For the example in section 4.2 this yields 2τ max ≈ 10 µs. Let the duration of the two pulses be τ = n∆t and t = m∆t, where the minimal time-step size ∆t = π/ω max , as follows from discrete Fourier theory. In turn the maximum frequency ω max of the target spectrum has to exceed all relevant eigenenergies of the system. The total duration for one sequence of pulses will then be τ tot = 4 N pts /2 n,m (n + m)∆t, where the number of points per frequency axis, N pts , is determined by ∆τ and τ max (and similar for the t axis). The expression gives τ tot = ∆tN 2 pts (1 + N pts ), and when including the number of repetitions discussed above and the number of phases N phases to be cycled over, the total duration is approximately We note that this is a lower limit, since we do not consider any time required for the initialisation of each measurement, such as trap loading. Assuming that a single repetition takes T rep ≈ 0.58 seconds, for the exemplary parameters of Fig. 6 we find T tot ≈ 16 hours. Recall, however, that our parameters are largely chosen to ease theory, and we expect that the total time can be lowered by a significant factor when using different parameters such as higher principal Rydberg quantum numbers up to n ≈ 80. Higher principal quantum numbers increase the dipole-dipole interaction strength C 3 in (9) according to C 3 ∼ n 4 , and consequently the system eigenenergy differences E n − E m . The 2D spectrum of interest then covers a larger frequency range, allowing less frequency resolution and hence much shorter delay-times, if the relative widths of peaks are preserved. For setups where the above duration is nonetheless a hindrance, we propose the simultaneous interrogation of an ensemble of identically prepared replica of the Rydberg aggregate as sketched in figure 6(c) by the same microwave pulse. Replicas should be far enough separated not to affect one another. Such a system can be realised using spatially structured laser fields for trapping the atoms in arbitrary 2D geometries [22,23]. As discussed in section 3.1, our signal is simply the total number of Rydberg atoms in the pstate, hence simultaneous interrogation would automatically average over all aggregates in the system. For the discussion above, one can thus significantly reduce the number of repetitions N rep for the ensemble average, since many repetitions would be provided already by a single experimental run.
In these simulations and those of the preceding section, both detector and aggregate atoms are expected to be at a precisely determined trapping location. We explore in Appendix F what happens if these locations themselves vary during acquisition of a single spectrum, and find that trapping centre fluctuations of spatial widths exceeding σ = 0.2 µm begin to require an increasing number of repetitions.

Conclusions and outlook
We have outlined a protocol for non-linear two-dimensional spectroscopy, which is a powerful tool in material science and physical chemistry, applied to the realm of ultracold Rydberg physics. Through the excellent control over the latter, this would enable a benchmark platform for the further development of algorithms and tweaks for multidimensional spectroscopic techniques, such as phase cycling schemes. Particularly helpful in this regard, is the ability to engineer different decay, disorder and dephasing sources in the Rydberg setting. We demonstrated this with several simple arrangements of a Rydberg dimer and Rydberg trimer embedded in a controllable environment, using realistic parameters.
To illustrate the points above, we did not need to vary the third time delay in our four pulse sequence, the so called waiting time T . By additionally scanning the waiting time, one could not only acquire information about system eigenstates, but also about time-evolution of the system in a non-equilibrium situation [45]. In one further step, by Fourier transforming the third delay, the technique would be augmented to 3D spectroscopy [46]. Finally, the technique can in principle be augmented to an increasingly larger number of pulses, furnishing ND-spectroscopy.
When going beyond the examples discussed above, using 2D spectroscopy on larger, more theoretically challenging Rydberg systems could also add a window on Rydberg dynamics that is complementary to the many existing techniques, providing new ways to look at dipolar-interacting quantum systems, especially those interacting with their environment [47], involving many-body resonant interactions [48], spin relaxation [49] or transport [50].
While we have chosen parameters for demonstrations in this article to yield most tractable theory, performing experiments would of course be even more interesting under conditions that challenge theory, in particular when the environment-aggregate system becomes non-Markovian [26,27]. This can also be achieved for example by allowing motion of the Rydberg aggregate atoms [51,52,53] to be induced by dipole-dipole interactions, which then can mimick for example molecular vibrations [54].

Appendix A. Model for Rydberg aggregate embedded in a decohering environment
The model considered in the present work, shown in figure 2, is discussed in detail in Ref. [15]. However, in that work the focus was on the single excitation manifold, i.e., states of the form | s, · · · , p, · · · , s with a single p-excitation. For 2D spectroscopy one must include also the total ground state and the states of the two exciton manifold. We thus provide an extension of the approximate effective master equation that includes those states in this appendix.
The total Hamiltonian, without the microwave, can be divided into three partŝ WhereĤ agg is the Rydberg aggregate Hamiltonian (8). The HamiltonianĤ EIT of the environment atoms that are driven by two lasers readŝ where | g , | e and | r are ground state, excited state and Rydberg state. All of these are involved in the EIT scheme for the environment atoms [15], which are labelled with the index α. Ω p and Ω c are probe and coupling Rabi frequencies and ∆ p and ∆ c are the respective detunings. The interaction between environment atoms and aggregate atoms is given byĤ The many-body states | χ k are defined in Eq. (13) of the main text andV kα reads, where h k i takes the values 0 for T k i = s and 1 for T k i = p. The interaction potential between the i-th aggregate atom in the state s and the α-th environment atom in the Rydberg state | r is given by V iα rs = C 6,rs |rα−r k | 6 . For the case when the aggregate atom is in the p state we discuss two different types of potential which can be selected by the choice of | r , V iα rp = C 4,rp |rα−r k | 4 , or V iα rp = C 6,rp |rα−r k | 6 , as discussed in Ref. [15]. These different choices are also indicated in table 1.
Finally one has to take into account the spontaneous decay of the intermediate state, on the probe transition from | e back to | e , with rate Γ p . This is accomplished by augmenting the equation of motion to a Lindblad master equation with decay operatorŝ L α = [| g e |] α , as usual [55].

Appendix A.1. Approximate effective master equation
The full many-body Hilbertspace with many environment atoms and several aggregate atoms is too large to handle numerically. For a wide range of laser parameters and couplings it is possible to adiabatically eliminate all excited states | e and | r of the environment atoms to obtain an evolution equation in the space of aggregate atoms alone. Following the method of Ref. [42] we define a projector, P g = k | χ k χ k | ⊗ | g g |, on the relevant state space, in which all environment atoms are in the ground state | g . The first part acts on aggregate atoms and second part on environment atoms (| g = | gg...g ). The complementary operator toP g is given byP e = 1 −P g . As a result of the adiabatic elimination procedure [42], we obtain an effective Hamiltonian and Lindblad operators,Ĥ eff andL (α) eff , which can be written asĤ αL α /2 is a non-Hermitian Hamiltonian, whileĤ g ,Ĥ e ,Ω + andΩ − are given byĤ g =P gĤPg =P gĤaggPg , and Evaluating the expressions above, and using the definitions (14) and (15), we obtain We have defined∆ p = ∆ p + iΓ p /2. As long as | χ k contains only a single excitation, this reduces to the expressions in [15].
Under the conditions discussed in section 3.2, we now can solve (12) using the results of the present section to predict all measurements on the Rydberg aggregate, including 2D spectra. This is computationally much easier than the full many-body evolution (1).

Appendix B. Perturbative calculation of system response
To understand the spectra shown in the present work, it is useful to consider analytical expressions obtained from standard perturbation theory [4] with respect to the aggregate-radiation interaction. For reference we provide here the basic equations needed in the Rydberg context.
We first diagonalizeĤ agg +Ĥ eff to obtain the eigenenergies and eigenstates of the Rydberg aggregate system Both,Ĥ agg andĤ eff , conserve the number of p-excitations, N p , hence the eigenstates can be classified according to N p . States with N p = 1 are denoted as one-exciton states and states with N p = 2 as two exciton states. For zero microwave detuning, (10) becomesĤ mw = E(t) N n=1μ n withμ n = µ 0 (σ (n) sp + σ (n) ps )/2. The matrix elements ofμ ≡ nμ n between eigenstates (B.1) are denoted by χ k |μ|χ k = µ kk and determine whether or not transitions between these states are allowed. The time evolution in (12), including the Lindblad terms from the environment, is written aṡ whereL 0 is the super-operator, indicated by˘, accounting for free evolution of the system without the microwave field. Following [4], we expand the time-evolving ρ to fourth order in the electric field such that the relevant contribution to our signal is withF defined in section 3.1. After the expansion, we find with the Green's matrix in the time domain defined by From (10) we see that we can writeL mw =VE(t), splitting the light-matter interaction into an operator and an electric field part, see (3). Substituting t 1 = τ 2 − τ 1 , t 2 = τ 3 − τ 2 , . . . t 4 = t − τ 4 and making use of the initial time, t 0 → −∞, we obtain, with the response function and electric field product Let us finally decompose the atom-field coupling super-operatorV intoV =V up +V down where the partV up contains the elements µ nσ ps of (10) that cause system excitations, whileV down contains those, µ nσ sp , causing de-excitations. HenceV down =V † up In the original Master-equation (12), the operators inV will act on the density matrix from the left as well as from the right, since these actions originate from a commutator [4]. We can represent this by further splitting the super-operators defined above in a leftand right component [56,57], according toV up =V L +V R with whereV up is the respective operator in Hilbert space. Equation (B.7) is a sum of terms involving four interactions with the electro-magnetic field that involve excitations or de-excitations acting on the left or the right of the density matrix. Between these interactions, the system propagates either in a population or coherence of the density matrix according to (B.5). Each of these terms can thus be represented by a double-sided Feynman diagram as given for example in figure 3. At this stage, the fourth order response (B.6) contains still contributions that are not of much interest, such as such as non-resonant terms or the square of second order responses. These can be removed through phase-matching when the system probed provides a medium for the interrogating fields, or through phase cycling when that is not the case.

Appendix C. Phase cycling technique
After obtaining the simple expression for the fourth order response of our system in (B.6), we now give a brief idea how to isolate contributions of interest using phasecycling. More information is widely available, [31,58,59,60,61,62,63].
We now also perform the rotating wave approximation (RWA), to ultimately reach the form (11) for the fields in the light-matter coupling Hamiltonian. Through the RWA, the complex phase e iφ j of the j'th pulse effectively entersV L , whileV † L contains e −iφ j . In contrast one obtains a factor e −iφ j inV R , whileV † R contains e iφ j . This can be summarized by the rule that all left-pointing arrows in diagrams as in figure 3 contribute the phase −φ j , while right-pointing arrows are contributing with +φ j .
Regardless of these details, it is evident from the preceding discussion and (B.6) and (B.3), is that the dependence of the signal on the phases of a given single pulse sequence will formally take the form where the integers α to δ count how many excitations or de-excitations occurred through a pulse with a given phase. We now recognise (C.1) as discrete Fourier representation of the signal S, hence the coefficientS αβγδ for a specific choice of indices can be extracted using an inverse Fourier transform from the space of phases to that of coefficients. Consider now the photon echo diagrams listed in figure 3. If we inspect the arrow directions and take into account the phase-sign allocation discussed at the beginning of this section, we find that all of them obtain a phase contribution −φ 1 + φ 2 + φ 3 − φ 4 . In terms of the Fourier decomposition (C.1) this corresponds to indices α = −1, β = 1, γ = 1, δ = −1. Since these are also all diagrams with this combination of phases, Figure C1. Two dimensional un-cycled spectrum of a monitored Rydberg dimer, for the same case as shown in figure 3(a). The signal intensity is clearly dominated by independent frequency features for each axis..
we can obtain the signal corresponding to the sum of all photon echo diagrams by the inverse discrete Fourier transform: with the index choice α = −1, β = 1, γ = 1, δ = −1.
The final question is on how many different phase data-points S in (C.1) needs to be sampled. Since each set of phases will take up experimental time, it is typically desirable to minimize the number of sets. Such schemes are discussed for example in Refs. [31,29]. We used a scheme which employs 27 different sets of phases. Since only relative phases between pulses can matter, we set the phase of the last pulse to zero φ 4 = 0, and then cycle the remaining ones in three discrete steps over the complex unit circle, φ 1 = l(2π/3), φ 2 = m(2π/3), φ 1 = n(2π/3), with l, m, n ∈ {0, 1, 2}, yielding a total of 27 combinations. Such a 27 step phase cycling scheme has been experimentally implemented in two dimensional fluorescence spectroscopy [63]. Using this scheme, we thus explicitly obtain the phase cycled signal as where we have now made the remaining dependence on time delays explicit. Using S pe (ω τ , T, ω t ) = 1 2π dτ dt e −iωτ τ e iω t t S pe (τ, T, t ), (C.4) this is then converted to frequency space (ω τ , ω t ).
For the case T = 0 one can reduce the number of phase-sets by usinḡ which requires only 18 phase triples. For a further intuitive idea on the functionality of phase cycling, consider the uncycled spectrum in figure C1, obtained from simulations of (12) for just a single set of phases: φ 1 = φ 2 = φ 3 = φ 4 = 0. This corresponds to the Rydberg dimer of figure 3(a). Clearly, the peak structure of figure 3(a) is not visible. Instead the spectrum is dominated by second order features at ω τ = 0 (or ω t = 0), which imply a signal that is independent of the respective delay τ (t ). This in turn implies that this signal is independent of pulse one (three). Whenever a contribution to S does not depend on one of the first three pulses, it is evident from (C.3) that it drops out in the summation over complex pre-factors, since e.g. 2 n=0 e in( 2π 3 ) = 0. In this manner leading second order contributions are removed and fourth order ones brought to the fore.

Appendix D. Perturbative calculation of photo echo contributions
Let us now continue where we left off in Appendix B, assuming that through phasecycling as discussed in the previous appendix, the system response in (B.7) is reduced to only those contributions corresponding to the photon echo. We show only diagrams where the ket is excited, each diagram has a vertically mirrored partner, which contributes the complex conjugate amplitude, such that the final contribution to the population is real. In the case of interaction on the ket, the first interaction term in (B.7) is V R . To comply with the restrictions discussed in the previous Appendix C, the photon echo component needs an equal number of inward and outward interactions, which reduces the number of terms in (B.7) to four, for which the Feynman diagrams are shown in figure 3. Thus we obtain a response function R(t 1 , t 2 , t 3 , t 4 ) = R 1 +R 2 +R 3 +R 4 with The next step is to convolve the response function R with the electric field pulses as in (B.6), which is simplified by the impulsive limit where the electric field pulse shapes are delta functions, A(t) = δ(t) in (2). We also express the response in the eigen-basis of the effective Hamiltonian, |χ k . Finally taking the Fourier transform of the signal, each response in (D.1) contributes a term Since there are no other time-arguments in the R j , in this impulsive limit for the field, Fourier transforms essentially act directly on the Green's functions.
For all cases considered in this work the Liouville propagator is diagonal in the chosen basis, which is a rather special case and requires the absence of relaxation processes that move population between different eigenstates. However that allows us to obtain simple expression for the 2D spectrum, asS(ω τ , T, ω t ) = 4 i=1S i (ω τ , T, ω t ) with contributions where the indices i, j run over all single exciton eigenstates | β and k runs over all two-exciton states | ζ . The factor 2 in the last equation comes from the fact that the final state in this path-way contains two p-excitations. All expressions use the Green's function in the frequency domain, given bỹ where ω kk = ω k − ω k is the transition frequency between states | k and | k , and is the effective decay rate in this state. In general Γ kk is complex and the imaginary part contributes a shift to the transition energy.
To compare with full numerical calculations, we calculate the total signal, S = 4 i S i , and typically plot the modulus of the resultant complex function. It is thus instructive to consider the absolute value ofG kk (ω), which is a Lorentzian with peak at ω kk − Im[Γ kk ] and a width given by Re[Γ kk ]. To perform a discrete Fourier transform, we additionally multiply the signal with a Gaussian window function prior to transform, where the width which is chosen as σ G = 0.3τ max , where τ max depends on the largest value of τ and t .
We can now apply the results obtained so far to the Rydberg dimer and trimer aggregates discussed in the main article.   Table D1. Parameters entering the analytical calculations (D.17) for the spectrum shown in figure 3, as discussed in the text. It is instructive to compare the analytical result with the numerical nonperturbative calculations, based on the numerical solution of the complete master equation (12). The results are shown in figure D1 along two selected one-dimensional cuts of figure 3. We find good agreement, validating both approaches. Residual small differences can be due to finite pulse duration or non-perturbative effects, which are included in the numerical solution, but not in the analytical one. Note that the impulsive limit is not necessary to obtain analytical results, but it yields simpler expressions, which for our short pulses allow a direct interpretation of the 2D spectra.  and the | ζ n can be obtained from the | β n by swapping all s states with p states. As for the dimer, the interpretation of spectra in the main article will require all transition energies ω kk with the respective transition dipoles as well as Γ kk as provided in table D2. Figure F1. Two dimensional spectrum of Rydberg dimer with position fluctuations during spectrum acquisition. Parameters are as in figure 3, except Ω mw /(2π) = 16 MHz, but the position of aggregate and detector atoms is given a random 2D offset with a Gaussian distribution of width σ = 200 nm. This offset is varied after each individual set of four phases has been simulated. We use the technique of Appendix E to model discrete atoms counts and average over N rep ∼ 10 5 repetitions.

Appendix E. Statistical Ensemble calculation
In order to calculate the number of repetitions of the experiment required, we numerically determine the signal S based on the count of a discrete number of Rydberg atoms in the p-state, as would be the case in experiment, rather than complete knowledge of the time evolved system density matrixρ(t). For this we first determine the reduced density matrix for a single aggregate atom, e.g. for atom 1 in a dimer: ρ 1 = Tr 2 {ρ} = (ρ ss,ss + ρ sp,sp )| s s | + (ρ ps,ps + ρ pp,pp )| p p |. (E.1) From the resultant probabilities, a random number decides whether the atom is measured in s or in p. This is repeated N rep times, to obtain a discretely sampled signalS (ξ) for each set of phases, which after phase cycling and Fourier transform gives rise to the graphs shown in figure 6. We conclude from this simulation, that about N rep ∼ 10 5 repetitions are required to resolve the spectrum for these parameters.

Appendix F. Trap center fluctuations
In this appendix we explore the sensitivity of the phase-cycling procedure, discussed in Appendix C, to position fluctuations of the aggregate atoms or the detector atoms.
Since we assume these to be trapped in the quantum ground state of their respective optical potential, such fluctuations would have to be due to imprecision in the alignment of these. For the simulation, the position of each atom is given a random offset, drawn from a two-dimensional Gaussian distribution with standard-deviation σ. The random offset is newly drawn after each set of phases, in contrast to the typical disorder situation in an optical setting. The results are shown in figure F1 for σ = 200 nm. The simulation employed the technique of the preceding appendix, and modelled the discrete counting of the number of p-excitatons in an experiment, averaging over N rep = 10 5 repetitions. We find, that for larger width σ than shown, the small side-peak is no longer resolved so that N rep would have to be further increased.