Probing weak dipole-dipole interaction using phase-modulated non-linear spectroscopy

TL;DR

This paper explores how phase-modulated non-linear spectroscopy can detect weak dipole-dipole interactions by analyzing peak intensities and phase shifts in the demodulated signals, offering a new approach to measure interaction strength.

Contribution

It introduces a theoretical framework combining numerical solutions and analytical perturbation theory to identify dipole-dipole interactions through spectral peak intensities and phase differences.

Findings

Dipole-dipole interactions can be detected via peak intensities in the spectrum.

The phase between harmonic signals depends on the sign of the interaction.

Parameter dependence of signal intensities provides insights into interaction strength.

Abstract

Phase-modulated non-linear spectroscopy with higher harmonic demodulation has recently been suggested to provide information on many-body excitations. In the present work we theoretically investigate the application of this method to infer the interaction strength between two particles that interact via weak dipole-dipole interaction. To this end we use full numerical solution of the Schr\"odinger equation with time-dependent pulses. For interpretation purpose we also derive analytical expressions in perturbation theory. We find one can detect dipole-dipole interaction via peak intensities (in contrast to line-shifts which typically are used in conventional spectroscopy). We provide a detailed study on the dependence of these intensities on the parameters of the laser pulse and the dipole-dipole interaction strength. Interestingly, we find that there is a phase between the first and second harmonic demodulated signal, whose value depends on the sign of the dipole-dipole interaction.

Tables

Table 1. Table 1: Parameters that are fixed in all numerical calculations shown. All energis are expressed in units of ω 0 subscript 𝜔 0 \omega_{0} . Time is in units of 1 / ω 0 1 subscript 𝜔 0 1/\omega_{0} .

Parameters	${\omega }_{\mathrm{fe}}$	${\mu }_e$	${\mu }_f$	${\omega }_M$	${\varphi }_{21}^{\left(0\right)}$	$m$	${\Omega }_{21}$	$T_{\mathrm{Rep}}$
Values	$0.05$	$0.75$	$1.054$	$0$	$0$	$1000$	$2\pi \times {10}^{-3}$	$1$

Equations

$$H$$

$$H_n$$

$$V_{n{n}^{\prime }}$$

$$V_{n{n}^{\prime }}^{\alpha \beta }$$

$$H_n^{\mathrm{int}}(t)$$

$$E(t)$$

$$S_{\mathrm{Fluor}}(t_{21},{\tau }_m)$$

$$P_{\mathrm{Fluor}}=\sum _n(|\mathrm{e}{⟩}_n⟨\mathrm{e}|+|\mathrm{f}{⟩}_n⟨\mathrm{f}|).$$

$${\tilde{S}}_{\mathrm{Fluor}}(t_{21},\kappa )$$

$$S_{\mathrm{ref}}(t_{21},{\tau }_m,\kappa )$$

$$S(\omega ,\kappa )=\frac{1}{\sqrt{2\pi }}{\int }_0^{\infty }d{t}_{21}{e}^{-i\omega t_{21}}{\tilde{S}}_{\mathrm{Fluor}}(t_{21},\kappa )$$

$$A_j(t-t_j-{\tau }_m)=E_j^{(0)}{e}^{-(t-t_j-{\tau }_m{)}^2/2{\sigma }^2}.$$

$${\mathcal{A}}_j(\omega )$$

$${\mathcal{A}}_{jj}(\omega ,\tilde{\omega })$$

$$A_j^{(\mathrm{g}\alpha ,1)}={\mathcal{A}}_j({\omega }_{\alpha \mathrm{g}}+V_{\alpha \alpha }).$$

$$A_{jj}^{(\alpha \alpha )}$$

$$A_{jj}^{(\mathrm{ef},\alpha )}$$

$$A_{211}^{(\mathrm{ge},2)}$$

$$A_{211}^{(\mathrm{gf},\bar{2})}$$

$$A_{211}^{(\mathrm{ge},\bar{2})}$$

$$A_{211}^{(\mathrm{gf},2)}$$

$$A_{211}^{(\mathrm{ge},3)}$$

$$A_{211}^{(\mathrm{gf},3)}$$

$$S^{(1)}(t_{21})=S_{\mathrm{e}}(t_{21})+S_{\mathrm{f}}(t_{21})$$

$$S_{\alpha }(t_{21})$$

$${\Lambda }_{\alpha }$$

$${\Phi }_{\alpha }$$

$$S^{(2)}(t_{21})=S_{\mathrm{ee}}(t_{21})+S_{\mathrm{ef}}(t_{21})+S_{\mathrm{ff}}(t_{21})$$

$$S_{\alpha \alpha }$$

$$S_{\mathrm{ef}}$$

$${\Lambda }_{\alpha \alpha }$$

$${\Phi }_{\alpha \alpha }$$

$${\Lambda }_{\mathrm{ef}}$$

$${\Phi }_{\mathrm{ef}}$$

$$A_{j}(t-t_{j}-\tau_{m})=\left\{\begin{array}[]{ccc}E_{0}&;&-\frac{\Delta}{2}<t<\frac{\Delta}{2}\\ &&\\ 0&;&{\rm{otherwise}}\end{array}\right.$$

$${\mathcal{A}}_j(\omega )$$

Full text

I Introduction

The study of many-body effects in weakly interacting ensembles of particles poses challenges to experiments. In many cases strong single-particle signals or broadening effects mask the weak collective signals. Recently, coherent time-resolved non-linear spectroscopy introduced a very sensitive approach to probe weak many-body effects with double-quantum two-dimensional (2D) spectroscopy StoneCundiffNelson09Science ; KaraiskajCundiff10PRL ; TurnerNelson10Nature ; DaiCundiff12PRL ; GaoCundiffLi16 where the two-body response is background-free detected, hence the mere presence of a double-quantum signal indicates the existence of inter-particle interactions in the system. Double quantum coherence measurements in 2D spectroscopy are very useful to detect weak interactions. However, most 2D schemes do not provide sufficient sensitivity to investigate dilute samples.

A multidimensional spectroscopy approach suitable for highly dilute samples is based on phase modulation (PM) TekavecMarcus07JCP . This approach has been demonstrated in 2D electronic spectroscopy TekavecMarcus07JCP but also in electronic wave packet interferometry TekavecMarcus06JCP . In the latter scheme, one records the signal from an observable (such as fluorescence intensity, photocurrent or photoionization TekavecMarcus06JCP ; NardinCundiff13 ; BruderStienkemeier2015PCCP ) emitted by the particles subject to two short laser pulses with a specific time-delay and a slowly modulated relative phase. Demodulation with respect to this relative phase and a subsequent Fourier transformation with respect to the time-delay yield a complex-valued spectrum. By using a variant of this technique, in Ref. BruderStienkemeier15PRA, the demodulation has been performed by using higher harmonics of the reference signal.

In the present work we aim at understanding the effect of long-range dipole-dipole interaction on the spectra resulting from higher order demodulation. To this end we investigate a minimal model consisting of two particles which can interact via transition dipole-dipole interaction. The model that we use is in particular related to interacting dye molecules which have a parallel arrangement. To infer the dipole-dipole interaction in such systems is of considerable interest (see e.g. Refs. WeSt03_125201_ ; FiSeEn08_12858_ ; Ei07_321_ ; GuZuCh09_154302_ ; WoMo09_15747_ ; DuNaPr15_072002_ ). The model of our measurement results effectively in a one-body measurement operator (number of detected photons). Our work will thus provide reference for the predictions from such a single-particle operator.

In our theoretical investigations we will on the one hand use a non-perturbative fully numerical approach, on the other hand we will derive analytical formulas valid if a perturbation treatment with respect to the light-matter interaction is applicable. In our numerical calculations we aim at following the experimental procedure as close as possible.

The structure of this paper is as follows. In Sec. II, we introduce a model of multiple particles that interact with laser fields and also a theoretical background of the phase demodulation technique used to reveal collective behaviors of particles. Taking two particles as a typical example, in Sec. III we present the fully numerical simulation results for Gaussian pulses. In Sec. IV, basing on perturbation theory, we develop analytical formulas of phase-(de)modulated fluorescence signal responsible for that measured in the experiment and present analytical results of signals for rectangular pulses in the appendix. Discussions and conclusions are finally given in Sec. VI. In a supporting information we provide details about our numerical and analytical procedures. Throughout the paper we set $\hbar=1$.

II The model

II.1 System and interaction with laser pulses

We consider an ensemble of particles that do not move and that are so far apart that only long-range interaction of the dipole-dipole type has to be taken into account. The Hamiltonian for $N$ particles coupled to a pulsed laser field reads

[TABLE]

The Hamiltonian $H_{n}$ describes the individual particles, $V_{nn^{\prime}}$ is the dipole-dipole interaction between particle $n$ and $n^{\prime}$, and $H_{n}^{\rm int}$ is the interaction of particle $n$ with the electromagnetic field. In the following subsections we discuss the various parts of the Hamiltonian above. Throughout the work we focus on the case of two particles (i.e. $N=2$). Our basic setup is sketched in Fig. 1.

II.1.1 Individual particles $H_{n}$

For each particle we consider three energy eigenstates $\{|{\rm g}\rangle_{n},|{\rm e}\rangle_{n},|{\rm f}\rangle_{n}\}$ and write the Hamiltonian of particle $n$ as

[TABLE]

Here $|{\rm g}\rangle$, $|{\rm e}\rangle$ and $|{\rm f}\rangle$ represent the ground, first and second excited states, respectively. We assume for simplicity and to elucidate the basic effects most clearly that there are dipole allowed transitions only between the ground state and each of the excited states, i.e. between $|{\rm g}\rangle$ and $|{\rm e}\rangle$, and between $|{\rm g}\rangle$ and $|{\rm f}\rangle$. The corresponding dipoles we denote by $\vec{\mu}_{{\rm ge}}^{(n)}$ and $\vec{\mu}_{{\rm gf}}^{(n)}$, respectively, and take them to be parallel (and identical for the two particles).

II.1.2 Dipole-dipole interaction $V_{nn^{\prime}}$

The basis states of two particles $n$ and $n^{\prime}$ can be constructed as product states of single particles, namely $\{|{\rm gg}\rangle,|{\rm ge}\rangle,|{\rm eg}\rangle,|{\rm gf}\rangle,|{\rm fg}\rangle,|{\rm ee}\rangle,|{\rm ef}\rangle,|{\rm fe}\rangle,|{\rm ff}\rangle\}$, where $|\alpha\beta\rangle=|\alpha\rangle_{i}|\beta\rangle_{j}$ ($\alpha,\beta={\rm g},{\rm e},{\rm f}$). The dipole-dipole interaction between two particles is then written in terms of these basis states as May-Kuehn

[TABLE]

For parallel transition dipoles (see Fig. 1) the interaction strengths are given by

[TABLE]

with $\vec{R}_{nn^{\prime}}$ being the vector between the two point-dipoles ($R_{nn^{\prime}}=|\vec{R}_{nn^{\prime}}|$), $\theta_{nn^{\prime}}$ is the angel between $\vec{R}_{nn^{\prime}}$ and the dipole orientation. Note that depending on the angle $\theta_{nn^{\prime}}$ the sign of the interaction can change from positive to negative. To facilitate the analytical considerations we neglect off-resonant contributions coupling the two excited states via dipole-dipole interaction $V_{nn^{\prime}}^{\rm ef}$. This approximation is good when the dipole-dipole interaction is much smaller than the energy separation between the two states (which is the case in which we are primarily interested in).

Note that from the definition Eq. (4) the interaction $V$ depends on the distance between the two particles $R$ and on the angle $\theta$. That means that for a given value of $V$ there exist many possible configurations $R$, $\theta$. Note also, that for certain angles $\theta$ the interaction $V$ can be negative.

Two particles: Eigenenergies and eigenstates: For two particles the Hamiltonian Eq. (1) becomes $H_{1}+H_{2}+V_{12}$ where $H_{n}$ and $V_{12}$ are given by Eqs. (2) and (3) respectively, where $V_{12}^{\rm ee}$ and $V_{12}^{\rm ff}$ will be abbreviated by $V_{\rm ee}$ and $V_{\rm ff}$ in the following. For the interpretation of the results it is convenient to use the eigenstates and eigenenergies. The energy levels without and with the dipole-dipole interaction are shown in Figs. 2(a) and (b) respectively.

The symmetric states, i.e. $|\psi_{{\rm gg}}\rangle$, $|\psi_{{\rm ge}}^{+}\rangle$, $|\psi_{{\rm gf}}^{+}\rangle$, $|\psi_{{\rm ee}}\rangle$, $|\psi_{{\rm ef}}^{+}\rangle$, and $|\psi_{{\rm ff}}\rangle$ span one subspace of the Hilbert space, while the other subspace is spanned by antisymmetric states, i.e. $|\psi_{{\rm ge}}^{-}\rangle$, $|\psi_{{\rm gf}}^{-}\rangle$, and $|\psi_{{\rm ef}}^{-}\rangle$. For the definition of these states see Fig. 2.

II.1.3 Particle-field interaction $H_{n}^{{\rm int}}$

The interaction between each particle and an electromagnetic field is described by

[TABLE]

with time-dependent coupling strengths $\Omega_{n}^{{\rm e}}(t)=-\vec{\mu}_{{\rm ge}}^{(n)}\vec{E}(t)$ and $\Omega_{n}^{{\rm f}}(t)=-\vec{\mu}_{{\rm gf}}^{(n)}\vec{E}(t)$ with $\vec{E}(t)$ the laser pulse field. Again for simplicity we take $\vec{E}(t)$ to be identical for the two particles and parallel to the transition dipoles of the particles (see Fig. 1(a)). For identical particles in our consideration, i.e. $\mu_{{\rm ge}}^{(n)}=\mu_{\rm e}$ and $\mu_{{\rm gf}}^{(n)}=\mu_{\rm f}$, the coupling strengths then become $\Omega_{n}^{{\rm e}}(t)=\Omega^{{\rm e}}(t)$ and $\Omega_{n}^{{\rm f}}(t)=\Omega^{{\rm f}}(t)$.

II.2 Modelling the measurement procedure

In the considered PM pump-probe experiment, two phase-modulated femtosecond (fs) laser pulses excite the system. Fluorescence intensity is recorded via time-integrated detection in a small angle perpendicular to the excitation laser propagation direction. Thereby it is not distinguished between fluorescence emitted from a many-particle or single-particle population. The many-body information is encoded in the time evolution between the two pulses resulting in a time-dependence of the fluorescence yield.

Each pulse is tagged with an individual phase signature stemming from the phase modulation, allowing for filtering of certain contributions of the final signal, as described below. This phase signature is manifested in a linear phase sweep $\phi_{j}(t)=\Omega_{j}t+\phi_{j}^{(0)}$ imparted by acousto-optic devices. The laser field consisting of a train of pulse pairs and illustrated in Fig. 1(b) is hence expressed as

[TABLE]

where $A_{j}(t-t_{j}-\tau_{m})$ is the pulse envelope function. Here $t_{j}$ is the arrival time of the $j$th pulse in a pulse pair and $\tau_{m}=mT_{\mathrm{Rep}}$ with $m$ being the number of pulse pairs and $T_{\mathrm{Rep}}$ the pulse repetition time. In Eq. (6), $\omega_{j}$ is the laser central frequency, $\phi_{j}^{(0)}$ is the phase offset, and the term $\Omega_{j}t$ describes the linear phase sweep due to the pulses passing individual acousto-optic modulators.

The PM $\phi_{j}(t)$ results in a slow intensity modulation of the fluorescence which allows one to use lock-in techniques after recording the fluorescence intensity with a photodetector. Lock-in detection can even be used to filter the fluorescence signal with respect to harmonics of the modulation frequency TianWarren2002OptLett ; KaKrMa16_015504_ ; BruderStienkemeier15PRA . In this way individual many-body responses are separated into different signal channels. Fourier-transforming the final demodulated signals yields complex-valued frequency spectra.

We assume that the signal produced by the photodetector is proportional to the number of photons emitted by the system which we take to be proportional to the number of excited particles after each pulse pair, ignoring emission during the short time-scale of the pulses. We take $T_{\mathrm{Rep}}$ much larger than the radiation life time of the particles. Therefore, we can consider each pulse pair individually. We use the state after this pulse pair (at time $T_{\rm F}$ sufficiently large) to calculate the number of photons that will be emitted after the pulse pair:

[TABLE]

where $\rho(t;t_{21},\tau_{m})$ is the density matrix at time t, obtained for the pulse delay $t_{21}=t_{2}-t_{1}$ and the modulation time $\tau_{m}$ and

[TABLE]

II.2.1 Signal demodulation

The signal Eq. (7) is modulated as a function of $\tau_{m}$ by the frequency difference $\Omega_{21}=\Omega_{2}-\Omega_{1}$. In the experiment harmonic demodulation with a lock-in amplifier as described in Ref. BruderStienkemeier15PRA, is applied. Basically, the lock-in amplifier multiplies the signal with a known reference waveform and applies subsequently a low-pass filter, thus the demodulated output signal can be described by:

[TABLE]

where $\tau_{\rm LI}$ is the lock-in time, $\kappa$ is the demodulation order, and $S_{\mathrm{ref}}(t_{21},\tau_{m},\kappa)$ is a complex-valued reference signal given by

[TABLE]

where $\phi_{21}^{(0)}=\phi_{2}^{(0)}-\phi_{1}^{(0)}$. In this work we focus on the first ($\kappa=1$) and second ($\kappa=2$) harmonic demodulation. We will abbreviate them as 1HD and 2HD, respectively. Note that for sufficiently large $\tau_{\rm LI}$ the demodulation Eq. (9) essentially results in extracting the frequency components with $\kappa\Omega_{21}\tau_{m}$ out of the signal and shifting the frequencies of the resulting function by $\kappa\omega_{\rm M}$. In our theoretical considerations we will use $\omega_{\rm M}=0$.

For visualization it is convenient to perform a Fourier-transformation with respect to $t_{21}$. The resulting (complex-valued) quantity

[TABLE]

will be denoted as the spectrum in the $\kappa$th order of demodulation.

III Fully numerical simulation

In our fully numerical calculations we propagate the Hamiltonian specified in section II (details on the numerical implementation can be found in the Supporting Information, section I.A).

In the following we consider two identical pulses, i.e. $\omega_{1}=\omega_{2}=\omega_{\rm L}$ and $E_{1}^{(0)}=E_{2}^{(0)}=E_{0}$ in a pulse pair. The pulse envelope was taken to be of Gaussian shape, i.e.

[TABLE]

To understand the dependence on the parameters we varied the pulse width $\sigma$, the strength of the dipole-dipole interaction $V$, and the field strength $E_{0}$ of the pulses. All parameters that were kept fixed for all calculations are given in Table 1. The unit of energy we denote by $\omega_{0}$ (possible choices are for example bandwidth of the pulses or even $\omega_{0}=1/T_{\rm Rep}$).

For each parameter set we perform propagations for 4500 values of $t_{21}$ equally spaced between $t_{21}=0$ and $t_{21}=2250$ and 1000 values of $\tau_{\rm m}$ with $T_{\rm Rep}=1$, resulting in the populations after the second pulse $P(t_{21},\tau_{\rm m}$). Performing a numerical demodulation (described in the Supporting Information I.B), for each $t_{21}$ we have obtained a signal $\tilde{S}_{\rm F1uor}(t_{21})$. This signal is then numerically Fourier transformed after multiplication with a Gaussian window function $e^{-{t^{2}_{21}}/{2\bar{\sigma}^{2}}}$ with $\bar{\sigma}=500$ to avoid Fourier transform artefacts. The final time and spacing of the $t_{21}$ values are chosen such that the peaks can be well resolved. Larger times $t_{21}$ are preferable to reduce pulse overlap effects. The temporal FWHM of the window function is chosen to be $2\sqrt{2\ln 2}\bar{\sigma}\sim 1177.41$ which results in a FWHM of $2\sqrt{2\ln 2}\frac{1}{\bar{\sigma}}\sim 0.0047$ in frequency domain. The resulting spectra exhibit sharp resonances which sit on a broad background. We identified this background to be caused by pulse overlap effects. In all results presented in the following, we have subtracted this background (details about the numerical procedure can be found in the Supporting Information, I.C).

Figs. 3 and 4 show example-spectra. Each row shows for a set of parameters the real and imaginary parts of the $1$st harmonic demodulation (1HD) (left) and $2$nd harmonic demodulation (2HD) (right) spectrum. From top to bottom the absolute value of the dipole-dipole interaction is increased. In Fig. 3 all dipole-dipole interactions are positive and in Fig. 4 they are negative. The other parameters are specified in the figure caption. The vertical dashed lines indicate the peak positions obtained from the analytical treatment of section IV.

We will refer to the two peaks in the 1HD spectra as ${\rm S_{e}}$ and ${\rm S_{f}}$, respectively (in Figs. 3 and 4 the ${\rm S_{e}}$ peak is located at a smaller frequency than that of the ${\rm S_{f}}$ peak). The three peaks in the 2HD spectra are denoted as ${\rm S_{ee}}$, ${\rm S_{ef}}$, and ${\rm S_{ff}}$ (in the order of increasing frequency).

Before we present results for a wide range of parameters, let us first briefly discuss some features that can be observed in Fig. 3. First note that the peak positions in the 1HD are shifted with respect to the position of the non-interacting particles by the energy of the dipole-dipole interaction. In contrast, in the 2HD they are always at the constant frequencies $2\omega_{\rm eg}$, $2\omega_{\rm fg}$ and $\omega_{\rm eg}+\omega_{\rm fg}$. While in the 1HD peaks appear in the real part, they appear in the imaginary part of the 2HD and have a negative sign. To obtain a similar behaviour as for the 1HD one has to multiply the 2HD by $e^{i\pi/2}$ for positive dipole-dipole coupling and by $e^{-i\pi/2}$ for negative dipole-dipole interaction. The amplitude of each peak in 1HD decreases slightly ($\sim 10$ %) upon increasing the dipole-dipole interaction. Remarkably, in the 2HD the peak amplitude increases by an order of magnitude. The last point is highlighted in Fig. 5 in more detail. In this figure peak amplitudes are plotted as a function of the strength of the dipole-dipole interaction, now also including negative values. The top panel shows 1HD and the bottom panel shows 2HD. The filled circles, stars and squares are from numerical calculations, and the solid, dashed and dot-dahsed lines from the analytical results of section IV.

The peaks at $2\omega_{\rm eg}$ and $2\omega_{\rm fg}$ have the same absolute value of the peak amplitude for the same absolute value of the dipole-dipole interaction strength. However, the peak at $\omega_{\rm eg}+\omega_{\rm fg}$ has a different absolute value of the intensity when changing the sign of the dipole-dipole interaction (as indicated by the thin dashed blue line and the unfilled triangles).

The effect of pulse width on the peak amplitude is shown in Fig. 6. Here we kept the ‘pulse-area’ $I=\int A_{j}(t){\rm d}t$ constant when changing the pulse width. For the used Gaussian pulses this implies a scaling of $E_{0}$ according to $E_{0}\sim 1/\sigma$. As expected the 1HD spectrum has a very weak dependence on the pulse width. Essentially it reflects strength of the pulse-spectrum at the field strength of the respective transition. The decrease at large pulse width comes from the fact that the spectral width of the pulse (which scales as $1/\sigma$) becomes comparable to the energetic separation between the transition frequency and the central frequency of the laser. It is remarkable that the peak amplitude in the 2HD has a very strong dependence on the pulse width. In particular, for short pulse-width (i.e. large bandwidth of the pulse) the peak-amplitudes are decreasing for decreasing pulse-width. This effect is related to the curvature (and magnitude) of the pulse-spectrum at the frequency of the relevant transitions. In addition it also depends in a non-trivial way on the difference between the pulse-spectrum at the ‘two-photon transition’ (i.e. $|{\rm gg}\rangle\rightarrow|{\rm ef}\rangle$) and the pulse-spectrum at the ‘one-photon transition’. In the Appendix B analytic calculations for rectangular pulses are performed in the perturbative regime with respect to the electric field (times transition dipole). From these results (in particular Eqs. (B.0.1)–(40)) one sees that for the case when the detuning of the laser with respect to the single particle transitions is large compared to the dipole-dipole interaction, then the amplitude of the peaks consists of the product of two parts: one part leads to a ‘trivial’ scaling (electric field $\times$ pulse width)n, where $n=2$ for 1HD and $n=4$ for 2HD. The normalization of the electric field strength used in Fig. 6 keeps this factor constant. The second factor does for the 1HD signal in leading order not depend on the dipole-dipole interaction and approaches one for decreasing pulse width. For the 2HD signal however, the leading term depends linearly on the dipole-dipole interaction (as we also found numerically in Fig. 5). The slope has non-trivial dependence on the pulse-width, and approaches zero for decreasing pulse-width. Let us note that often simulations of non-linear spectroscopy signals are performed using delta pulses. From the foregoing considerations one sees that one has to be careful not to miss important features.

Finally, in Fig. 7 we show the dependence on the electric field strength. The upper panel shows the 1HD peak intensities as a function of $E_{0}^{2}$ and the lower panel shows 2HD peak intensities as a function of $E_{0}^{4}$. The linear slope indicates that we are in the perturbative regime for field strength smaller than about $E_{0}\approx 0.01$. For larger values deviation from the linear behaviour becomes apparent.

IV Perturbation theory with respect to the light-matter interaction

In the following subsections we firstly present our general formulas for harmonic demodulated signals that are based on perturbation theory and applicable to arbitrary pulses. For Gaussian pulses considered in our fully numerical simulation, our formulas together with a numerical integration of time-ordered double integrals show consistent results with the full numerical calculations and are displayed together with the numerical results. In Appendix B we apply our formulas to rectangular pulses for which fully analytical expressions can be obtained. These results are not directly comparable to the results of our numerical simulations but they help to develop physical understanding about the dependence of the signal on the various parameters.

IV.1 General formulas

In section II of the Supporting Information we provide information on our perturbative treatment and present intermediate results. Corresponding double-sided Feynman diagrams are provided in Appendix C. Here, in the main text, we will only provide final results after demodulation.

Note that in the derivation of the formulas presented below, we have neglected overlap between the two pulses.

IV.1.1 Some definitions

It is convenient to define for the pulses $j=1,2$ a single-interaction spectral amplitude

[TABLE]

and a double-interaction spectral amplitude

[TABLE]

Furthermore it is convenient to define the following ‘amplitudes’:

[TABLE]

[TABLE]

IV.1.2 The first harmonic demodulated signal

The first harmonic demodulated signal is

[TABLE]

with

[TABLE]

where

[TABLE]

From the expression Eq. (24) one sees that the 1HD spectrum consists of two contributions. From Eq. (25) together with Eq. (27) one sees that each contribution results in a single peak located at the transition frequencies $\omega_{\alpha\mathrm{g}}+V_{\alpha\alpha}$ of the interacting particles ($\alpha={\rm e},\,{\rm f}$). The intensity of these peaks is determined by Eq. (26), which is a real function.

IV.1.3 The second harmonic demodulated signal

The second harmonic demodulated signal is

[TABLE]

with

[TABLE]

The coefficients and phases are given by

[TABLE]

and

[TABLE]

From the expression Eq. (28) one sees that the 2HD spectrum consists of three contributions. These contributions correspond to peaks at two times the transition frequencies of the uncoupled particles (see Eq. (29) and Eq. (32)) and to one peak that has a frequency of the sum of these transition frequencies (see Eq. (34)). The ‘intensity’ of these peaks is determined by Eq. (31) and Eq. (33), which are in general complex functions.

IV.1.4 Spectra in frequency domain

The ‘spectrum’ can be obtained by applying the Fourier transformation, namely, $\mathcal{F}\{u(t)\cos(\omega_{0}t)\}=\frac{1}{2}\sqrt{\frac{\pi}{2}}[\delta(\omega+\omega_{0})+\delta(\omega-\omega_{0})]-\frac{1}{2}\frac{i}{\sqrt{2\pi}}(\frac{1}{\omega+\omega_{0}}+\frac{1}{\omega-\omega_{0}})$ and $\mathcal{F}\{u(t)\sin(\omega_{0}t)\}=\frac{i}{2}\sqrt{\frac{\pi}{2}}[\delta(\omega+\omega_{0})-\delta(\omega-\omega_{0})]+\frac{1}{2}\frac{1}{\sqrt{2\pi}}(\frac{1}{\omega+\omega_{0}}-\frac{1}{\omega-\omega_{0}})$ where $u(t)$ is the unit step function ($u(t)=1$ for $t\geq 0$ and [math] for $t<0$) and $\mathcal{F}\{x(t)\}=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}x(t)e^{-i\omega t}dt$. The factor $\omega_{0}$ has to be chosen according to the phases $\Phi_{\alpha}$, $\Phi_{\alpha\alpha}$ and $\Phi_{\rm ef}$ defined in Eqs. (27), (32) and (34).

In order to compare with results given by the numerical simulation, the signals are further multiplied by the same Gaussian window function $e^{-{t_{21}^{2}}/{2\bar{\sigma}^{2}}}$ 111By using the convolution theorem of Fourier transformation, we have $\Re[\mathcal{F}\{u(t)\cos(\omega_{0}t)w(t)\}]=\frac{\bar{\sigma}}{4}[e^{-\frac{(\omega+\omega_{0})^{2}\bar{\sigma}^{2}}{2}}+e^{-\frac{(\omega-\omega_{0})^{2}\bar{\sigma}^{2}}{2}}]$ and $\Im[\mathcal{F}\{u(t)\sin(\omega_{0}t)w(t)\}]=\frac{\bar{\sigma}}{4}[e^{-\frac{(\omega+\omega_{0})^{2}\bar{\sigma}^{2}}{2}}-e^{-\frac{(\omega-\omega_{0})^{2}\bar{\sigma}^{2}}{2}}]$. Here we used $w(t_{21})=e^{-\frac{t_{21}^{2}}{2\bar{\sigma}^{2}}}$ and $\mathcal{F}\{w(t_{21})\}=\bar{\sigma}e^{-\frac{\omega_{0}^{2}\bar{\sigma}^{2}}{2}}$..

IV.2 Gaussian pulses – comparison to numerics

The result for Gaussian pulses can be given by our formulas above together with a numerical integration of time-ordered double integrals in Eq. (14). The results are shown in Figs. 5, 6, and 7 as solid, dashed and dot-dashed lines. One sees a very good agreement with our fully numerical simulation. The reason for this good agreement is that the electric field strength of the pulses in the numerical simulation has been chosen to be in the perturbative regime. This can be seen in Fig. 7 where we present numerical calculations as a function of the field strength. The linear dependence of the 1HD with $E_{0}^{2}$ and the 2HD with $E_{0}^{4}$ indicates the validity of the perturbative treatment. The solid, dashed and dot-dashed lines are again results from the analytical perturbative treatment. One sees that only for the last two numerical points deviations from the perturbative calculations become visible.

Deviations between analytical and numerical calculation are caused by the long pulse length used in the numerical simulations, where overlap of the first and second pulse is fully included (while it is ignored in the analytical treatment). We have also performed numerical calculations with shorter pulse length and found much better agreement. One example is shown in Appendix A.

V Applicability to molecular and atomic systems

Although in the present study our focus is not on a particular system, we will in the following comment on two different systems. The first one is the case of organic molecules in solution or on a substrate, the second one are atoms in a gas, similar to the situation of the experiments of Bruder et al. BruderStienkemeier15PRA .

V.1 Organic dye molecules

One application of the present method is the investigation of dipole-dipole interacting organic dye molecules. For such molecules, the electronic transition dipoles are typically in the order of several Debye. Dimers with specific interaction strength could be formed for example by linking the organic molecules by spacers ZhLiZh15_20108_ ; C5SC04956C , by self-assembly in helium-nanodroplets WeSt03_125201_ ; RoEiDv11_054907_ or by depositing the molecules on a surface with a small coverage MPaMa13_064703_ ; MIzVl15_121408_ (using single molecule spectroscopy techniques one then can select close-by molecular pairs). The distances between the molecules can range from a few Ångström to several nanometer. For close molecules the strength of the dipole-dipole interaction can be of the order of $\pm$100  cm*-1*. At a distance of 10  nm one has interactions in the order of $\pm$0.1  cm*-1*. While interactions of the order of $\pm$100  cm*-1* lead to line shifts that can usually be well resolved by linear absorption spectroscopy, for interactions around $|V|\approx 1\;{\rm cm}^{-1}$ this becomes difficult because of line-broadening. Let us briefly relate these numbers to the present calculations. Taking as unit of energy $\omega_{0}=100\,{\rm cm}^{-1}$, then in Figs. 3, 4 and 5 the used pulse width $\sigma=10.4$ corresponds to $\approx 3\,{\rm ps}$ and the detuning $\omega_{\rm eg}-\omega_{\rm L}=-2.5\,{\rm cm}^{-1}$. For most molecules the electronically excited states are energetically well separated. Then the states $|\rm e\rangle$ and $|\rm f\rangle$ would correspond to vibronic states. If their energetic separation is much larger than the spectrum of the pulse, then essentially only one state contributes. However, for energetically very close energies (which is in particular the case for highly excited vibrational states) then two (or even more) states can contribute. The above example shows that one might be able to obtain information on the dipole-dipole coupling strength in the case when this interaction is weak compared to the energetic disorder of the individual molecules.

V.2 Dilute potassium gas

In the experiments of Ref. BruderStienkemeier15PRA, a dilute gas ($10^{10}\,$cm*-3*) of potassium atoms has been investigated. In the resulting spectra at higher harmonic demodulation, peaks are visible at the energies that correspond to multiple excited particles. However, it is not clear how to model the resulting spectra (e.g. understand the intensities of the various peaks). It has been suggested that the applied PM scheme reduces the experimental observable to an effective many-particle operator Mukamel2016JCP . As a consequence, collective resonances should be observed independent of the presence of interactions among the atoms in the system. However, it is not obvious why the experimental procedure should result in a many-body measurement operator since the fluorescence is detected via a photodiode in the proportional range.

Although our model may neglect some aspects relevant in an atomic gas of many particles, we will nevertheless briefly relate the present results to the above mentioned experiment. For a potassium atom, we make the following assignments: $|{\rm g}\rangle=|4{}^{2}\!S_{1/2}\rangle$, $|{\rm e}\rangle=|4{}^{2}\!P_{1/2}\rangle$, and $|{\rm f}\rangle=|4{}^{2}\!P_{3/2}\rangle$. In the 1HD one sees two peaks corresponding to the transition ${\rm g}\rightarrow{\rm e}$ and ${\rm g}\rightarrow{\rm f}$ at frequencies $\omega_{\rm eg}$ and $\omega_{\rm fg}$, respectively (for potassium one has $\frac{\mu_{\rm f}^{2}}{\mu_{\rm e}^{2}}\sim 1.96$). The 2HD spectrum has peaks at $2\omega_{\rm eg}$, $2\omega_{\rm fg}$ and $\omega_{\rm eg}+\omega_{\rm fg}$. In the experiment, the magnitude of these peaks is approximately one order of magnitude smaller as the one of the 1HD peaks.

To obtain a rough prediction from our theoretical model we use experimental parameters for pulse intensities and width. The dipole-dipole interaction we estimate from the mean particle separation in the gas. Using these numbers we find that the 2HD is six orders of magnitude smaller than the 1HD. This discrepancy between our model and the experiment could have several reasons: In the PM experiments Ref. BruderStienkemeier15PRA, high repetition rate lasers have been used. The repetition rate is faster than the radiative lifetime of a single potassium atom. Thus the system may not be relaxed to the ground state before a consecutive pulse train perturbs the sample again. Such an effect is discussed in Ref. OsShHa16_053845_, for the case of a two-photon transition. In the experiment not two atoms at a well-defined distance are addressed but the laser interacts with a large number of atoms in a gas. In our estimations we have assumed that the dominant contributions come from nearest neighbours and simply considered the average nearest neighbor distance. In a more detailed calculation one has to average over distances and relative orientations (also the assumption that the two atoms see the same field at a certain time is no longer correct, since the distance between the atoms is much larger than the wavelength). Note that depending on the relative orientation one can have positive and negative values of the dipole-dipole interaction. Regarding the orientation of the transition dipoles some additional care is necessary. For atoms the transition dipoles are not simple vectors. And in particular for the dipole-dipole interaction one has to take the detailed level substructure into account (see e.g. the discussions in Refs. RoHeTo04_042703_, and MoeWueAt11_184011_, ). In our theoretical modelling we have ignored correlation effects in the emission (like stimulated emission or superradiance).

VI Discussion and conclusion

In this work, we have investigated systematically the collective excitations of particles probed via phase-modulated non-linear spectroscopy within the framework of a one-body measurement operator. We have focused on two particles that can exchange interaction via transition dipole-dipole interaction. The laser-pulses had frequencies overlapping with the transition frequencies of the individual particles. We developed analytical formulas for the (de)modulated signals in terms of a single- and double-interaction pulse spectral amplitudes. We obtained analytical results for rectangular pulses and presented fully numerical simulation for Gaussian pulses. Below we summarize our main findings:

The peak positions in the 1HD spectra are shifted with respect to the position of the non-interacting particles by the energy of the dipole-dipole interaction. In contrast, in the 2HD they are always at the constant frequencies $2\omega_{\rm eg}$, $2\omega_{\rm fg}$ and $\omega_{\rm eg}+\omega_{\rm fg}$.

The complex 1HD spectrum has the typical features of absorption and dispersion: Peaks appear in the real part (of the complex spectrum) and the corresponding dispersion is in the imaginary part. However, for the 2HD spectrum there is an important difference: Peaks appear in the imaginary part and dispersion-like features are in the real part. If one wants to have a similar functional form as the 1HD feature (to interpret the signal as absorption and dispersion) one has to multiply the 2HD by a phase-factor. This phase factor depends on the sign of the interaction; it is $e^{i\pi/2}$ for positive dipole-dipole coupling and by $e^{-i\pi/2}$ for negative dipole-dipole interaction.

For the amplitudes of the 2HD spectrum we found a linear dependence with respect to the strength of the dipole-dipole interaction.

We found a strong effect of the pulse-shape on the peak amplitude in the 2HD. In particular we found that it has a non-monotonic behaviour as a function of pulse width. For delta-pulses the 2HD signal vanishes.

Let us briefly comment on some of our approximations. In the present work we have neglected all decoherence effects. We expect that this is a good approximation, as long as the linewidth of the transitions of the individual atoms is small compared to the bandwidth of the pulses. In the present treatment all correlations in the field have been neglected. For example stimulated emission or Dicke-superradiance could have an effect on the observed fluorescence DorfmanMukamel13PhysRevA .

In the present work we have tried to follow the experimental procedure of Ref. BruderStienkemeier15PRA, as closely as possible. However, instead of sampling with $\tau_{m}$ and performing a demodulation one can also employ a phase-cycling scheme, experimentally as well as numerically. In particular in the numerics this will give a considerable speed up and one can still treat non-perturbative effects and arbitrary pulse-shapes. However, it will be difficult to treat effects specific to the demodulation procedure. For example in the present numerical approach one can model the lock-in amplifier in more detail or consider the effect of fast pulse repetition compared to the decay time of the particles.

Finally, let us note that our results indicate that phase modulation approach with higher order demodulation offers a possibility to detect weak interactions between particles when there is strong inhomogeneous broadening present. This is based on the fact that in the 2HD spectra dipole-dipole interaction manifests itself not in a line-shift (which is difficult to extract from a broadened spectrum), but from a peak amplitude.

The coupling strength could be obtained from the non-trivial dependence of the peak amplitude of the 2HD spectrum on the pulse width, which is shown in Fig. 6 and also illustrated by analytical results e.g. Eqs. (43) and (44) for rectangle pulses and small pulse width. The pulse width could be experimentally readily varied with a pulse shaper incorporated in the experimental setup.

Appendix A Pulse length effect for the difference between analytics and numerics

In this appendix, we provide an example that shows how the agreement between our perturbation treatment and the full numerics depends on the length of the laser pulses. In Fig. 5, it has been shown that there is deviation between analytical and numerical calculations. In this appendix, we perform calculations with shorter pulse length $\sigma=6$. A much better agreement shown in Fig. 8 supports our claim that the deviation in Fig. 5 is caused by the long pulse length $\sigma=10.4$ used in the numerical simulations.

Appendix B Fully analytical solutions for rectangular pulses

In order to have a better understanding of physics behind and gain some insights from analytical results one would like to evaluate the various $\mathcal{A}$ of Eqs. (13) and (14). Due to time-ordered double integrals in Eq. (14), analytical solutions can only be obtained in some special cases. Here we consider such a case where both pulses have a rectangle shape, namely,

[TABLE]

with $\Delta$ being pulse width. The single-interaction spectral amplitude defined in Eq. (13) becomes

[TABLE]

where the sinc function is defined as $\mathrm{sinc}[x]=\frac{\sin x}{x}$. The double-interaction spectral amplitude given by Eq. (14) becomes

[TABLE]

We want to emphasize that these two amplitudes, i.e. Eqs. (36) and (LABEL:eq:double_pulse_spectrum_rectangle) or their original definitions in Eqs. (13) and  (14) are key building blocks of harmonic signals of collective excitations.

It is important to note that the $\mathcal{A}_{jj}$ is in general complex-valued, while $\mathcal{A}_{j}$ is always real.

The expression Eq. (LABEL:eq:double_pulse_spectrum_rectangle) is still quite complicated. In the following, we perform series expansions to consider a regime of weak dipole-dipole interaction and also investigate the effect of finite duration of pulses.

B.0.1 Dipole-dipole interactions small compared to detuning

Here we consider the case $|\frac{V_{\rm ee}}{\delta_{\rm eg}}|\ll 1$ and $|\frac{V_{\rm ff}}{\delta_{\rm fg}}|\ll 1$ where the detunings with respect to the laser central frequency are defined as $\delta_{\rm eg}=\omega_{\rm eg}-\omega_{\rm L}$ or $\delta_{\rm fg}=\omega_{\rm fg}-\omega_{\rm L}$.

Then Eq. (26) for the first harmonic signals becomes

[TABLE]

The second harmonic ones, Eqs. (31) and (33) become

[TABLE]

with

[TABLE]

The intensity of collective resonances scales linearly with the dipole-dipole interaction. Collective signals become zero when the dipole-dipole interaction is absent.

Note that the lowest order term of $\Lambda_{\rm e}$ and $\Lambda_{\rm f}$ is real and positive, while for $\Lambda_{\rm ee}$, $\Lambda_{\rm ff}$, and $\Lambda_{\rm ef}$, they are positive and imaginary. After inserting $\Lambda_{\rm ee}$, $\Lambda_{\rm ff}$, and $\Lambda_{\rm ef}$ into $S_{\rm ee}$, $S_{\rm ff}$, and $S_{\rm ef}$ in Eqs. (29) and (30) respectively, the second harmonic demodulated signals given by the demodulation of $S_{\rm ee}$, $S_{\rm ff}$, and $S_{\rm ef}$ have a phase shift $\frac{3\pi}{2}$ or $\frac{\pi}{2}$ for positive and negative dipole-dipole interaction respectively, from the first harmonic demodulated signal (i.e. phase demodulation of Eq. (25)). This phase shift that does not appear when considering delta pulses, has also been observed in our simulation in Figs. 3, 4, and 5.

B.0.2 Dipole-dipole interaction small compared to bandwidth of pulse

Here we consider $|V_{\rm ee}\Delta|\ll 1$ or $|V_{\rm ff}\Delta|\ll 1$ . Up to the third order, we have

[TABLE]

where the detunings with respect to the laser central frequency are again defined as $\delta_{\rm eg}=\omega_{\rm eg}-\omega_{\rm L}$ or $\delta_{\rm fg}=\omega_{\rm fg}-\omega_{\rm L}$ and are also taken to be small. The expression for $\mathcal{B}$ is

[TABLE]

As in the previous subsection, the intensity of collective resonances scales in lowest order linearly with the dipole-dipole interaction and collective signals become zero when the dipole-dipole interaction is absent. Also the second harmonic demodulated signals given by the demodulation of $S_{\rm ee}$, $S_{\rm ff}$, and $S_{\rm ef}$ have a phase shift $\frac{3\pi}{2}$ or $\frac{\pi}{2}$ for positive and negative dipole-dipole interaction respectively, from the first harmonic demodulated signal. This disappears for delta-pulses.

It is interesting to note that the imaginary part of $\Lambda_{\alpha\alpha}$, which gives rise to the absorptive features in the spectrum, is an odd function of $V_{\alpha\alpha}$. Therefor the intensities for positive and negative dipole-dipole interaction $|V_{\alpha\alpha}|$ have the same magnitude but different signs. This is not the case for $\Lambda_{\rm ef}$. This is because of the third order term containing $\mathcal{B}$. One finds that the absolute value of signal for $V_{\alpha\alpha}<0$ becomes larger than that for $V_{\alpha\alpha}>0$. This is exactly the phenomenon that we have observed in simulation for Gaussian pulses.

Appendix C Relation to double-sided Feynman diagrams

The origin of the theoretical formulas of section IV can be easily understood by considering the relevant double-sided Feynman diagrams using the general rules provided e.g. in Ref. Mukamel95, . In Fig. 9 all diagrams are displayed that result in a phase $\Omega_{21}\tau_{m}$ (a and b) and $2\Omega_{21}\tau_{m}$, and therefore contribute to the 1HD and 2HD signal, respectively. Below each diagram the resulting amplitude is provided. Note that the appearance of a factor two in front of a diagram reflects the fact that in the final state two particles are excited and consequently will emit two photons.