Discriminating Cascading Processes in Nonlinear Optics; A QED Analysis Based on Their Molecular and Geometric Origin

TL;DR

This paper uses a quantum electrodynamic approach to distinguish direct and cascading nonlinear optical signals based on their phase, considering molecular and geometric factors, to improve interpretation of complex multidimensional signals.

Contribution

It introduces a QED-based method to predict the phase difference between direct and cascading signals, aiding experimental differentiation and understanding of many-body effects in nonlinear optics.

Findings

Predicted phase dependence on microscopic dynamics and geometry.

Identified conditions to distinguish direct and cascading signals.

Discussed higher-order multi-photon exchange processes.

Abstract

The nonlinear optical response of a system of molecules often contains contributions whereby the products of lower-order processes in two separate molecules give signals that appear on top of a genuine direct higher-order process with a single molecule. These many-body contributions are known as cascading and complicate the interpretation of multidimensional stimulated Raman and other nonlinear signals. In a quantum electrodynamic (QED) treatment, these cascading processes arise from second-order expansion in the molecular coupling to vacuum modes of the radiation field, i.e., single-photon exchange between molecules, which also gives rise to other collective effects. We predict the relative phase of the direct and cascading nonlinear signals and its dependence on the microscopic dynamics as well as the sample geometry. This phase may be used to identify experimental conditions for distinguishing the direct and cascading signals by their phase. Higher-order cascading processes involving the exchange of several photons between more than two molecules are discussed as well.

Figures

Equations

$$\left(\frac{{\partial }^2}{\partial t^2}-c^2{\nabla }^2\right)\mathbf{E}(\mathbf{r},t)=-\frac{{\partial }^2}{\partial t^2}\mathcal{P}(\mathbf{r},t)$$

$$\mathbf{E}(\mathbf{k},\omega )=\frac{{\omega }^2}{{\omega }^2-k^2{c}^2+i\eta }\mathcal{P}(\mathbf{k},\omega )$$

$$\mathbf{E}(\mathbf{k},\omega )=i\pi \delta ({\omega }^2-k^2{c}^2){\omega }^2\mathcal{P}(\mathbf{k},\omega ).$$

$$\hat{\mathbf{E}}(\mathbf{r},t)\equiv \sum _j{\mathbf{E}}_j(\mathbf{r},t)+{\hat{\mathbf{E}}}_{\text{v}}(\mathbf{r},t)$$

$${\hat{H}}_{\text{int}}=-\int d\mathbf{r}\hat{\mathbf{E}}(\mathbf{r},t)\cdot \hat{\mathcal{V}}(\mathbf{r},t)=-\int d\mathbf{r}\sum _j{\mathbf{E}}_j(\mathbf{r},t)\cdot \hat{\mathcal{V}}(\mathbf{r},t)-\int d\mathbf{r}{\hat{\mathbf{E}}}_{\text{v}}(\mathbf{r},t)\cdot \hat{\mathcal{V}}(\mathbf{r},t)={\hat{H}}_{\text{LM}}+{\hat{H}}_{\text{vM}}$$

$$\displaystyle S(\bm{k}_{s},\Lambda)=\Im\bigg{[}\int d\bm{r}dt\mathbf{E}^{*}_{s}(\bm{r},t)\cdot\langle\hat{\bm{\mathcal{V}}}(\bm{r},t)\rangle\bigg{]}=\Im\left[\mathbf{E}_{s}^{*}(\omega_{s})\cdot\bm{\mathcal{P}}(\bm{k}_{s},\omega_{s})\right]$$

$$\mathcal{P}(\mathbf{r},t)=\text{Tr}\left[\hat{\mathcal{V}}(\mathbf{r},t)\mathcal{T}{e}^{-i{\int }_{-\infty }^{t}d\tau {\hat{H}}_{-,\text{LM}}(\tau )}{e}^{-i{\int }_{-\infty }^{t}d\tau {\hat{H}}_{-,\text{vM}}(\tau )}{\rho }_{\text{eq}}\right]$$

$${\mathcal{P}}^{\nu }({\mathbf{k}}_s,\omega )=P^{\nu }({\mathbf{k}}_a,\omega )f({\mathbf{k}}_a-{\mathbf{k}}_s)+\int d{\omega }_b{\tilde{P}}^{\nu {\nu }_{\text{v}}}({\mathbf{k}}_a,{\omega }_s;-{\omega }_b)P^{{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_b,{\omega }_b)G^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_s-{\mathbf{k}}_a,-{\mathbf{k}}_b,{\omega }_b).$$

$$f(\mathbf{k})\equiv \sum _a{e}^{i\mathbf{k}\cdot {\mathbf{r}}_a}.$$

$$G(\mathbf{k},{\mathbf{k}}^{\prime },\omega )=\sum _{ab}{e}^{-i(\mathbf{k}\cdot {\mathbf{r}}_a+{\mathbf{k}}^{\prime }\cdot {\mathbf{r}}_b)}G({\mathbf{r}}_a,{\mathbf{r}}_b,\omega )=\sum _{ab}{e}^{-i(\mathbf{k}\cdot {\mathbf{r}}_a+{\mathbf{k}}^{\prime }\cdot {\mathbf{r}}_b)}\left(-{\nabla }^2{\delta }_{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}+{\nabla }_{{\nu }_{\text{v}}}{\nabla }_{{\nu }_{\text{v}}^{\prime }}\right)\frac{e^{i\frac{{\omega }_b}{c}{r}_{ab}}}{2\pi r_{ab}}$$

$$G^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_s-{\mathbf{k}}_a,-{\mathbf{k}}_b,{\omega }_b)\to f({\mathbf{k}}_a+{\mathbf{k}}_b-{\mathbf{k}}_s)G^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}(-{\mathbf{k}}_b,{\omega }_b)\to f({\mathbf{k}}_a+{\mathbf{k}}_b-{\mathbf{k}}_s)\frac{{\omega }_b^2}{{\omega }_b^2-k_b^2{c}^2+i\eta }\frac{-4{\delta }_{{\nu }_{\text{v}},{\nu }_{\text{v}}^{\prime }}}{3}$$

$${\mathcal{P}}^{\nu }({\mathbf{k}}_s,\omega )=P^{\nu }({\mathbf{k}}_a,{\omega }_s)f({\mathbf{k}}_a-{\mathbf{k}}_s)-i\frac{4\pi n{k}_{b}c}{3}{\tilde{P}}^{\nu {\nu }_{\text{v}}}({\mathbf{k}}_a,{\omega }_s;\mp k_{b}c)P^{{\nu }_{\text{v}}}({\mathbf{k}}_b,\pm k_{b}c)f({\mathbf{k}}_a+{\mathbf{k}}_b-{\mathbf{k}}_s).$$

$$G^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_s-{\mathbf{k}}_a,-{\mathbf{k}}_b,{\omega }_b)=\sum _a{e}^{i({\mathbf{k}}_a+{\mathbf{k}}_b-{\mathbf{k}}_s)\cdot {\mathbf{r}}_a}{I}^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)$$

$$I^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)=\sum _b{e}^{i{\mathbf{k}}_b\cdot ({\mathbf{r}}_a-{\mathbf{r}}_b)}{G}^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{r}}_b,{\omega }_b),$$

$$I^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)\to I^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_b,{\omega }_b)=G^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{k}}_b,{\omega }_b)$$

$$H_{\mathrm{int}}=-\int d\mathbf{r}\phi (\mathbf{r})\rho (\mathbf{r}),$$

$$G(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=\frac{e^{i\frac{\omega }{c}|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$I({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)={\int }_{V}d{\mathbf{r}}_b{e}^{i{\mathbf{k}}_b\cdot ({\mathbf{r}}_b-{\mathbf{r}}_a)}G({\mathbf{r}}_a,{\mathbf{r}}_b;{\omega }_b),$$

$$I={\int }_0^R{r}_b^{2}d{r}_b{\int }_0^{\pi }\sin \theta d\theta {\int }_0^{2\pi }d\phi {\xi }^3(\theta ,\phi )(\xi (\theta ,\phi )r_b{)}^{-1}{e}^{i\xi (\theta ,\phi )(\frac{{\omega }_b}{c}+k_b\cos \theta )r_b},$$

$$I=-{\int }_0^{2\pi }d\phi {\int }_0^{\pi }\sin \theta d\theta \frac{1}{(\frac{{\omega }_b}{c}+k_b\cos \theta {)}^2}+{\int }_{S^2}d\mathbf{n}\left(\frac{1}{(\frac{{\omega }_b}{c}+{\mathbf{k}}_b\cdot \mathbf{n}{)}^2}-\frac{iR\xi (\mathbf{n})}{\frac{{\omega }_b}{c}+{\mathbf{k}}_b\cdot \mathbf{n}}\right)e^{iR\xi (\mathbf{n})(\frac{{\omega }_b}{c}+{\mathbf{k}}_b\cdot \mathbf{n})}.$$

$$I=-\frac{4\pi }{{\left(\frac{{\omega }_b}{c}\right)}^2-k_b^2}.$$

$$I={\int }_0^R{r}_{b}d{r}_b{\int }_0^{2\pi }d\phi {\xi }^2(\phi )(\xi (\phi )r_b{)}^{-1}{e}^{i\xi (\phi )(\frac{{\omega }_b}{c}+k_b\cos \phi )r_b},$$

$$I=-i{\int }_0^{2\pi }d\phi \frac{1}{\frac{{\omega }_b}{c}+k_b\cos \phi }+i{\int }_0^{2\pi }d\phi \frac{1}{\frac{{\omega }_b}{c}+k_b\cos \phi }{e}^{i(\frac{{\omega }_b}{c}+k_b\cos \phi )a\xi (\phi )}=I_0+I_1.$$

$$I=I_0=-i{\int }_C\frac{dz}{2\pi iz}\frac{1}{\frac{{\omega }_b}{c}+k_b(z+z^{-1})/2}=-i{\int }_C\frac{dz}{2\pi i}\frac{1}{k_b{z}^2+\frac{{\omega }_b}{c}z+k_b}=\frac{-i}{\sqrt{{\left(\frac{{\omega }_b}{c}\right)}^2-k_b^2}}.$$

$$I({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)=2\pi i\frac{c}{{\omega }_b}{z}_a{e}^{i{k}_b{z}_a}.$$

$$G^{\nu {\nu }^{\prime }}({\mathbf{r}}_a,{\mathbf{r}}_b,{\omega }_b)=\left(-{\nabla }^2{\delta }_{\nu {\nu }^{\prime }}+{\nabla }_{\nu }{\nabla }_{{\nu }^{\prime }}\right)\frac{e^{i\frac{{\omega }_b}{c}r}}{2\pi r}=\frac{-e^{i\frac{{\omega }_b}{c}r}}{2\pi r^3}\left[({\delta }_{\nu {\nu }^{\prime }}-3{\hat{r}}_{\nu }{\hat{r}}_{{\nu }^{\prime }})(1-i\frac{{\omega }_b}{c}r)-({\delta }_{\nu {\nu }^{\prime }}-{\hat{r}}_{\nu }{\hat{r}}_{{\nu }^{\prime }}){\left(\frac{{\omega }_b}{c}\right)}^2{r}^2\right],$$

$$G({\mathbf{k}}_b,{\omega }_b)=\frac{-4}{3}\frac{{\omega }_b^2}{{\omega }_b^2-k_b^2{c}^2+i\eta }$$

$$S_{\text{c}}({\mathbf{k}}_s,\Lambda )=-n^2\Im \left[E_s^{\ast {\nu }_s}({\omega }_s)\int d{\omega }_b{P}^{{\nu }_s{\nu }_v}({\omega }_s;{\omega }_b)P^{\perp {\nu }_v}({\mathbf{k}}_b,{\omega }_b)\delta ({\mathbf{k}}_a+{\mathbf{k}}_b-{\mathbf{k}}_s)\frac{k_b^2{c}^2}{{\omega }_b^2-k_b^2{c}^2+i\eta }\right],$$

$$I^{\nu {\nu }^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)=\frac{-{\delta }_{\nu {\nu }^{\prime }}{e}^{i(\frac{{\omega }_b}{c}r-{\mathbf{k}}_b\cdot \mathbf{r})}}{2\pi r^3}\left[(1-3{\delta }_{\nu {\nu }^{\prime }z})(1-i\frac{{\omega }_b}{c}r)-(1-{\delta }_{\nu {\nu }^{\prime }z}){\left(\frac{{\omega }_b}{c}\right)}^2{r}^2\right].$$

$${\mathcal{P}}_{3-\text{body}}^{\nu }({\mathbf{k}}_s,{\omega }_s)$$

$$\times \sum _{abc}{e}^{i({\mathbf{k}}_a+{\mathbf{k}}_b+{\mathbf{k}}_c-{\mathbf{k}}_s)\cdot {\mathbf{r}}_a}{I}^{{\nu }_{\text{v}}{\nu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_b,{\omega }_b)I^{{\mu }_{\text{v}}{\mu }_{\text{v}}^{\prime }}({\mathbf{r}}_a,{\mathbf{k}}_c,{\omega }_c)$$

Full text

I introduction

Ultrafast nonlinear optical signals in the condensed phase usually use heterodyne detection, whereby the signal interferes with a reference local oscillator beam. The signal then scales linearly with the number of molecules $N$ and is distinguished by its power dependence in the incoming fileds and its direction (phase matching) Mukamel (1995); Boyd (1992). Matter information is contained in multipoint correlation functions of the dipole operator. We will refer to such signals as direct. While linear spectroscopy reveals the one-photon excitation spectrum of the sample, nonlinear spectroscopies contain multi-photon transitions and can probe dynamics in a controlled way. For example, in a pump-probe stimulated Raman experiment, a sample is excited by two interactions with an initial light pulse to a vibrational coherence which is then probed by a second pulse following a time delay $T_{1}$. This technique probes the vibrational structure of the molecule and is a 3rd-order process in the electric field (and is therefore related to $\chi^{(3)}$, the 3rd-order susceptibility). An extension of this technique is a 5th-order Raman process in which, following the initial delay $T_{1}$, a second pair of interactions transfers the sample from one vibrational coherence to another. The system then evovles freely during a second interpulse delay $T_{2}$ before being probed. This more general 5th-order ($\chi^{(5)}$) technique offers the possibility of probing the coupling between participating vibrational modes.

Another type of signal is also generated by a sequential process: some fields interact with one molecule to generate a field which then acts on another molecule, together with the external fields, to finally produce a signal. Such signals scale as $N^{2}$ and are known as cascading. The cascading signals have the same power dependence in the incoming fields and the same phase matching and are thus hard to distinguish from the direct signal Blank et al. (1999); Golonzka et al. (2000); Astinov et al. (2000, 2000); Kubarych et al. (2003); Li et al. (2008); Mehlenbacher et al. (2009); Bennett and Mukamel (2014). In neat liquids, cascading dominates the direct signal due to the quadratic vs. linear scaling in molecular density. It is then hard to separate the direct signal, which carries a higher level of molecular information. An important example is the direct 5th-order Raman $\chi^{(5)}$ process which is accompanied by a product of two 3rd-order $\chi^{(3)}$ signals. Separating the two had drawn considerable attention and took several years to recognize Tanimura and Mukamel (1993); Tominaga and Yoshihara (1996a, b); Steffen and Duppen (1997a, b); Tokmakoff and Fleming (1997); Blank et al. (1999); Golonzka et al. (2000); Kubarych et al. (2003); Li et al. (2008). This has been the main obstacle for multidimensional Raman spectroscopy (5th-and-higher orders) Mehlenbacher et al. (2009).

A number of ideas have been pursued to separate these undesired cascading signals. The $N$ vs. $N^{2}$ scaling is an obvious way to separate the cascading and direct contributions but in many cases (e.g., 5th-order Raman in neat molecular liquids), it is not possible to vary the molecular density over a sufficiently-large range in order to seperate the $N$- and $N^{2}$-scaling contributions. Another idea comes from the fact that, in an infinite homogeneous sample, it follows from Maxwell’s equations that the electric field emitted by a polarization is $\frac{\pi}{2}$ phase-shifted relative to the polarization itself. From a macroscopic perspective, the electric field $\mathbf{E}$ created by a polarization $\bm{\mathcal{P}}$ is given by the Maxwell equation

[TABLE]

which in the frequency domain reads

[TABLE]

where $\eta$ is a positive infinitessimal selecting the appropriate boundary conditions of an outgoing wave. In the limit of perfect phase matching $\omega^{2}\to k^{2}c^{2}$, we have

[TABLE]

This gives a $\pi/2$ phase shift between the polarization and the resulting field. We will derive a more general, geometry-dependent result and show under what conditions it reduces to this simple macroscopic relation. The extra emission event (from the source molecule) in a cascading process thus renders the cascading signal out of phase relative to the direct nonlinear signal ($E\propto P-iPP$) Mukamel (1995); Li et al. (2008); Frostig et al. (2015). This property has been parlayed into a number of successfully-implemented techniques to select for the direct 5th-order Raman signal and eliminate the cascading contributions Golonzka et al. (2000); Blank et al. (1999); Zhao et al. (2011); Frostig et al. (2015) (techniques based on polarization-sensitive measurements have also been pursued Gelin and Domcke (2013)).

In this paper, we extend our earlier treatment of cascading to study this phase shift on a microscopic level Bennett and Mukamel (2014). We show that the $\pi/2$ phase-shift described above occurs in the limit of macroscopic homogeneity only for certain geometries and does not hold for samples consisting of only a few particles. The essence of this geometry-related phase comes from 2D integrations as was described by Feynman for the emitted field by an infinite sheet of dipoles Feynman et al. (2013). After tracing over quantum field degrees of freedom, general expressions for cascading signals are given that are valid in the non-perturbative regime with respect to arbitrary external fields. The result can then be expressed solely as a product of two single-molecule dipole correlation functions and a geometric factor describing the relative positions of the molecules in the sample. This presentation allows the calculation of cascading signals using only the single-molecule calculations familiar from standard treatments of nonlinear spectroscopy and a sum of c-numbers associated with the geometry. This paper extends our previous line of reasoning Bennett and Mukamel (2014) which framed cascading signals as a subset of a more general class of vacuum-mediated interactions which includes local-field effects and other higher-order processes. We then obtained microscopic expressions for the first few orders (in external fields) for these vacuum-mediated effects. The focus of this paper is the study of the macroscopic continuum limit and a comparison with a microscopic collection of (few) molecules.

II Direct vs. Cascading Signals

In the semiclassical approach to nonlinear optical spectroscopy, a quantum material system interacts with only a few classical modes of the radiation field and the remaining infinite number of modes, which are in their vacuum states, are neglected. However, taking the vacuum modes into account in a fully quantum electrodynamical (QED) treatement can result in several many-body effects caused by photon exchange, i.e., the emission of a photon by one molecule and its absorption by another. These effects include the scrambling of the time ordering of incoming short pulses where the free-induction decay produced by one molecule can be long-lived and interacts with another molecule, a local-field $\chi^{(1)}$ effect Mukamel (1995); Lozovoy et al. (2001); Cundiff (2002); Bennett and Mukamel (2014). Additional local-field effects include corrections to the transmission/reflection of a thin film Benedict et al. (1991) and the Rabi oscillations of a quantum dot Paspalakis et al. (2006). Nonlinearities are also induced in ensembles of noninteracting harmonic oscillators which are otherwise linear Glenn et al. (2015). Dipole-dipole coupling, responsible for, e.g., Forster resonant energy transfer (FRET) and spontaneous quantum synchronization, is also induced by the exchange of photons Thirunamachandran (1980); Zhu et al. (2015). Superradiance, a cooperative spontaneous emission process, is another well-studied effect that finds its origin in the quantum nature of the radiation field Dicke (1954). Quantum-field effects due to quantum fluctuations of the laser field rather than the unoccupied modes of the electromagnetic field have also been noted in stimulated signals León-Montiel et al. (2016). This paper analyzes the aforementioned cascading contribution to heterodyne-detected nonlinear spectroscopic signals by treating them as the lowest-order many-body correction due to photon exchange.

Cascading shares a common origin with other collective effects observed in spectroscopic experiments. For example, superradiant emission also arises from collective coupling of the constituent molecules of the sample to a photon mode and scales quadratically in the number of molecules Gross and Haroche (1982). However, in superradiance the spontaneously emitted photons (or the associated excited-state population decay) are ordinarily the object of detection while in cascading the exchanged photons are virtual and the relevant experiment is a heterodyne detection with respect to the externally-applied signal (laser) field. The Lamb shift and collective analogues are also due to virtual-photon exchange but these effects are observed as alterations of a material resonance rather than an entirely new signal as in cascading Scully (2009). Forster transfer (discussed more thoroughly in the conclusion) can also be derived perturbatively as a vacuum-mediated interaction. When the field-matter interaction is treated in the dipole approximation, all of the above depend on the dipole-dipole coupling tensor, giving the derivations a similar flavor Thirunamachandran (1980); Salam (2010a).

We consider a sample made of identical, noninteracting molecules with non-overlapping charge distributions for which the coupling to the radiation field can be treated in the dipole approximation. The total system dipole operator then takes the form of a sum over molecular dipoles $\hat{\bm{\mathcal{V}}}(\bm{r},t)=\sum_{a}\hat{\mathbf{V}}(t)\delta(\bm{r}-\bm{r}_{a})$. The material is subjected to a set of classical laser modes $\mathbf{E}_{j}$ so that the total electric field is

[TABLE]

where $\hat{\mathbf{E}}_{\text{v}}$ is an electric field operator representing the infinitely-many vacuum modes. The material system is coupled to this total electric field via the interaction Hamiltonian

[TABLE]

where the laser modes and the quantum vacuum interact with the material via Hamiltonians $\hat{H}_{\text{LM}}$ and $\hat{H}_{\text{vM}}$ respectively.

Heterodyne-detected nonlinear spectroscopic signals are given by the rate-of-change of photon number in some signal mode $s$, taken to be in a coherent state. Commuting the photon number operator with Eq. (5) and integrating over time to get the total photon number change gives the signal

[TABLE]

where $\mathbf{E}_{s}$ is the signal field, which we have assumed to have a precisely-defined propagation direction $\hat{\bm{k}}_{s}$, and $\Lambda$ stands for the set of parameters defining the classical fields which must be specified to simulate particular signals. The total polarization is given by the expectation value of the total dipole operator $\langle\hat{\bm{\mathcal{V}}}(\bm{r},t)\rangle=\bm{\mathcal{P}}(\bm{r},t)$, and can be written in terms of time-ordered exponentials

[TABLE]

where $\rho_{\text{eq}}$ is the equilibrium field+matter density matrix and “$-$” subscript stands for the commutator $\hat{H}_{-}\rho=\hat{H}\rho-\rho\hat{H}$. Expanding Eq. (7) order-by-order in $\hat{H}_{\text{LM}}$ but to zeroth order in $\hat{H}_{\text{vM}}$, i.e., neglecting the vacuum modes altogether, yields the standard semiclassical nonlinear optical signals, which we term the direct signal $S_{\text{d}}$. All field-matter interactions in the perturbative expansion then occur on a given molecule and the signal is obtained by summing over molecules. These signals thus scale linearly with the number of molecules in the sample $N$, are proportional to the single-molecule signal suitably averaged, and contain no cooperative many-body contributions.

The presence of the quantum vacuum modes of the radiation field creates corrections to this picture. We can systematically generate such corrections by expanding the second exponent in Eq. (7) (the vacuum-mode terms). Since $\text{Tr}[\hat{a}^{(\dagger)}|0\rangle\langle 0|]=0$, the 2nd-order expansion is the lowest non-vanishing correction. This represents successive emission and re-absorption of a photon by the material. When the same molecule both emits and absorbs this photon, the result is the lowest-order radiative correction, or one-loop correction, to the energy (Lamb shift with radiative decay Cohen-Tannoudji et al. (1992)). The terms in which the absorber and emitter are different molecules lead to a transfer of a coherent excitation between the two molecules which is the origin of cascading. This is depicted diagrammatically in Fig. 1.

A derivation given in Appendix A gives the total system polarization $\mathcal{P}$ in terms of the individual molecular polarizations $P$ as

[TABLE]

where the $\nu$’s denote cartesian coordinates (with summation implicit), $P^{\nu}(\bm{k},\omega)$ is the polarization of a single molecule and $\tilde{P}^{\nu\nu^{\prime}}(\bm{k},\omega;\omega^{\prime})$ is the polarization of a molecule resulting from a single perturbative interaction with the polarization of another molecule in the sample and arbitrarily-many interactions with the classical fields (defined formally as a dipole correlation function in Eq. (35)). In Eq. (8), $\bm{k}_{a}$ and $\bm{k}_{b}$ stand for any linear combination of the set of incoming classical field modes and represent the set of laser modes that interact with molecules $a$ and $b$ respectively in a perturbative expansion. Since we work in the dipole approximation, the $\bm{k}$-dependence of the polarizations comes only as $\delta$-functions, originating in the spatial phase factor, that specify the linear combinations. In practice, one must sum over the possible subsets that generate different choices of the $\bm{k}_{a}$, $\bm{k}_{b}$ (we have omitted integration over $d\bm{k}_{a}d\bm{k}_{b}$ for brevity). Finally, $G$ is the photon Green’s function defined in $(\bm{r},t)$-space in Eq. (53).

The first term in Eq. (8) is the direct nonlinear signal and comes proportional to the sample’s form factor

[TABLE]

This factor carries information on the position of the molecules in the sample and, in the continuum limit, goes over to the delta function $f(\bm{k})\to(2\pi)^{3}n\delta(\bm{k})$, where $n=\frac{N}{V}$ is the molecular concentration. This corresponds to momentum conservation and yields the phase matching condition. The second term in Eq. (8) is the cascading signal in which the polarization of molecule $b$ serves as a source, along with the external fields, for the polarization of molecule $a$. Note that when this is expanded to order $m$ in the classical modes, we will have $\sum_{p+q=m-1}P^{(p)}P^{(q)}$. The traditional nomenclature is to refer to those terms with $p,q=1$ as local-field corrections and the remainder as cascading (the most familiar being $p=q=3$ but 3,5 cascading has also been discussed Cho et al. (2000)). Even though Eq. (8) represents a more general class of 2nd-order, vacuum-mediated interactions, we will refer to them collectively as cascading. The key quantity in Eq. (8) that connects the polarization emitted by molecule $b$ with the effective field “felt” by molecule $a$ is the $(\bm{k},\omega)$-space photon Green’s function

[TABLE]

which contains all information about the sample geometry and determines the relative phase of the cascading vs. direct contributions to the signal. Moreover, Eq. (8) reveals the important point that, in the off-resonant limit where $P(\omega)$ only has the phase imprinted by the external fields, the relative phase of cascading versus direct processes will be controlled entirely by $G(\bm{k},\bm{k}^{\prime},\omega)$. In the limit of a large, homogeneous sample, the spatial dependence of the photon Green’s function is reduced to the single variable $\bm{r}_{a}-\bm{r}_{b}$, the difference vector, and we have

[TABLE]

where $G(\bm{k},\omega)$ directly controls the phase (i.e., if this quantity is real then the cascading process is in phase with the direct process and if it is purely imaginary the two processes are out of phase) and the second relation follows on taking the integral in the continuum limit. In this regime, the total polarization is thus

[TABLE]

where we have assumed perfect phase-matching $k_{b}c=\omega_{b}$ and both choices of $\pm$ must be summed over.

Equation (12) is in the form of the simple macroscopic expression discussed in the introduction that is in common use to understand the cascading contribution to nonlinear signals. It reveals that the cascading contribution carries an additional factor of molecular density $n$ compared to the direct signal. Recalling that the form factor also scales linearly with $n$ in the continuum limit, the direct signal is linear while the cascading is quadratic. The crucial factor of $i$ responsible for the phase shift that is used to filter out cascades is found to originate in the phase-matching condition. Finally, it is worthwhile to note the linear dependence on $k_{b}c$ in the cascading contribution, indicating that it will be stronger at higher, such as x-ray, frequencies.

In making the assumption that $k_{b}c=\omega_{b}$, we have neglected the principal value of the denominator in Eq. (11) which generates a cascading contribution that is in phase with the direct signal. Additionally, the situation can be expected to be different in the case of few-molecule samples or oddly-shaped macroscopic systems. In the next section, we examine these effects in greater detail.

III Sample Geometry and the Phase of the Cascading Signal

From the previous section, it is clear that, for the off-resonant response, the $\bm{k}$-space photon Green’s function (Eq. (11) determines the phase of the cascading signal. In the semiclassical description, the effective electric field emitted by molecule $b$ is given by a phase-shift of the polarization $E\propto iP$. The extra factor of $i$, which in this semiclassical picture comes from the additional emision event inherent in cascading processes, renders the cascading terms out of phase with the direct nonlinear terms. This effect has been used to discriminate between the two Golonzka et al. (2000); Blank et al. (1999); Zhao et al. (2011); Frostig et al. (2015). In this section, we evaluate the photon Green’s function for different geometries and show under what conditions the cascaded signal acquires this well-defined phase shift.

To simplify the analysis, we will separate the summations (integrations) over the positions over the two molecules. We thus re-write

[TABLE]

where the second summation differs from the discrete Fourier transform of $G$ with respect to $\bm{r}_{b}$ only by the spatial phase factor $e^{i\bm{k}_{b}\cdot\bm{r}_{a}}$. In the macroscopic homogenous limit, $I$ does not depend on the choice of $\bm{r}_{a}$ and we have

[TABLE]

in line with Eq. (11). Thus, we will primarily be concerned with the evaluation of Eq. (14).

III.1 A Simplified Treatment for a Scalar Field

The presence of the dipole coupling tensor $\left(-\nabla^{2}\delta_{\nu\nu^{\prime}}+\nabla_{\nu}\nabla_{\nu^{\prime}}\right)$ in Eq. (10) complicates matters somewhat. For clarity, we first treat the problem by replacing the electromagnetic field, which is a vector gauge field, with a scalar field. This simplified model still maintains all qualitative features of the original model. The effects will be extended for the original electromagnetic field setting later and we will show that the scalar model describes the relation between the cascaded and direct signals absolutely adequately.

We assume the following interaction Hamiltonian

[TABLE]

where $\rho(\bm{r})$ is the scalar polarization, and $\varphi(\bm{r})$ is the scalar field that replaces the full electromagnetic counterpart, so that the Maxwell equations are replaced with the Helmholtz equation ($(\nabla^{2}+k^{2})\phi(\bm{r})=0$), whose Green’s function (that replaces the Green’s function of the Maxwell equations, also known as the Green’s function of the electromagnetic field) has the form

[TABLE]

We will further switch to the uniform continuum limit wherein all properties of the cascaded signal, compared to its direct counterpart are fully contained in the integral

[TABLE]

where $V$ denotes the integration region, occupied by the sample. We will focus on three special cases. In cases (i) and (ii) there is a poor phase matching, whereas in case (iii) phase matching is good.

Case (i): the region $V$ is a convex $3D$ region with the size large, compared to the wavelength. Introducing an inhomogeneous spherical coordinate system, associated with the convex region $V$, centered at $r_{0}$, with the $z$-axis oriented in the direction of $\bm{k}_{0}$ (see Appendix E) we recast the integral $I$ in the form

[TABLE]

where $R$ and $\xi(\theta,\varphi)$, are determined by the shape and size of the region $V$, as well as the position $\bm{r}_{a}$. Assuming a natural condition $k_{b}<\omega_{b}$, integration over $r_{a}$ can be performed for any values of $\theta$ and $\varphi$, resulting in

[TABLE]

The second integral over the unit sphere represents the integral over the boundary of the sample. If the point $\bm{r}_{a}$ is not within the wavelength region from the boundary (which is the typical case), the integrand is a fast oscillating function and the integration can be performed using the saddle point method. In Appendix F we estimate this contribution and demonstrate that the second integral is negligible compared to the first one. Performing the integration in the first term in the r.h.s. of Eq. (20) explicitly, we arrive at

[TABLE]

This result is natural and just means that, when $\bm{r}_{a}$ is not within a wavelength distance from the boundary, the result is the same as in an infinite medium. Since the contribution of the thin layer in the vicinity of the boundary is minor, the cascaded signal can be computed using the simple expression for $I$, given by Eq. (21). Note that we have implicitly assumed poor phase matching, i.e., Eq. (21) only applies in a principal value sense and neglects the perfect phase matching point $\omega_{b}=k_{b}c$. In the absense of such phase matching, $I$ is clearly real and the cascading signal is in phase with the direct signal. We will see later that the singular point gives an imaginary, out-of-phase contribution which is dominant under good phase matching.

Case (ii): $V$ is a convex $2D$ region with the size large, compared to the wavelength, such as a molecular monolayer or thin film with thickness smaller than the wavelength. We can use an inhomogeneous polar coordinate system, resulting in the following expression for the integral $I$

[TABLE]

where in this $2D$ case $\bm{k}_{b}$ naturally denotes the projection of the cascading field wave vector onto to the monolayer/film plane. Performing integration over $r_{b}$ we obtain

[TABLE]

Neglecting the second contribution coming from a fast-oscillating integral as in the $3D$ case, we obtain

[TABLE]

This means that, in the off-resonant case where the susceptibilities are real, the direct and cascaded signals are in phase for $3D$, whereas in the $2D$ case there is a $\pi/2$ phase shift between the two signals.

Case (iii): the region $V$ is a cylinder of generically irregular shape located on a reflecting plane. Let $\bm{r}=(x,y,z)$ describing Cartesian coordinates. We assume the region $V$ to be a cylinder of thickness $l$, so that $0\leq z\leq l$, where $z=0$ identifies a reflecting plane. The shape of the cross section of $V$ at $0\leq z\leq l$ is generally irregular but convex and may depend on $z$. We also assume $\bm{k}_{b}=(0,0,k_{b})$ to be directed along the $z$-axis, i.e., the cascading wavevector is normal to the plane. The derivation of this signal follows a similar form to case (ii) but is followed by an integration over $z$. It is given in detail in Appendix C and results in

[TABLE]

where we have taken the limit of good phase matching so that $|\omega_{b}-k_{b}|l\ll 1$ and $z_{a}$ is not within the wavelength scale from the borders.

To summarize, we note that the $i$ factor in the r.h.s. of Eq. (72) appeared due to the $2D$-integration over a cross section (Eq. 65)), exactly along the lines of the Feynman’s argument Feynman et al. (2013). The third contribution in Eq. (70), which is due to the reflecting surface, vanishes in the good phase matching limit. This tells us that in this limit the effect of the reflecting boundary is negligible, even if it is not perfectly reflecting. Finally, we note that, for a co-linear beam geometry, we have $\bm{k}_{b}c=\omega_{b}$ and this phase-matching condition is automatically satisfied. In this case, as long as a sufficient number of molecules are involved so as to permit the continuous $2D$-integration over a cross-section of the interaction region, we will obtain a cascading signal exactly out-of-phase with respect to the direct signal.

III.2 The Full Vector Guage Field

The previous section demonstrated how to perform the sorts of integrals needed to evaluate $G(\bm{r}_{a},\bm{k}_{b},\omega_{b})$ for different geometries utilizing a simplified scalar field. Restoring the vector nature of the coupling gives for the real space Green’s function,

[TABLE]

where we have used $r=|\bm{r}_{a}-\bm{r}_{b}|$ for brevity. Equation (26) has terms proportional to $r^{-3}$, $r^{-2}$, and $r^{-1}$. The last of these terms is purely transverse and, aside from the $\omega_{b}^{2}$ factor, matches the Helmholtz Green’s function used as an example in the previous section. In the limit of an isotropic sample, we have $\hat{r}_{\nu_{v}}\hat{r}_{\nu^{\prime}_{v}}\to\frac{1}{3}\delta_{\nu_{v}\nu^{\prime}_{v}}$ and this $r^{-1}$ term is the only to survive yielding

[TABLE]

where we have omitted the cartesian components since the result is proportional to $\delta_{\nu_{v}\nu^{\prime}_{v}}$. This is the equivalent to Eq. (21) for the vector coupling field and holds under the same geometric assumptions as that equation. In fact, in the isotropic limit, we may directly obtain the result for the full vector Green’s function by simply multiplying the corresponding result obtained above for a scalar field by the factor $\frac{\omega_{b}^{2}}{3\pi c^{2}}$. More generally, the far-field term is dominant when $\frac{\omega_{b}r}{c}\gg 1$, corresponding to seperations of greater than $2\pi\lambda$ or roughly 6 times the wavelength of the cascading light. Thus, for intermolecular seperations much less than this, the static dipole-dipole coupling is dominant while for much larger distances, the long-range $\frac{1}{r}$ term is dominant. As demonstrated above, it is this far-field term that leads to a factor of $i$ in the case of a thin sample under good phase matching $|k_{b}-\frac{\omega_{b}}{c}|l\ll 1$. Note that the $r^{-2}$ term comes with an $i$ factor and hence, when the sample is anisotropic, we may expect that the cascading signal has both in-phase and out-of-phase contributions relative to the direct nonlinear signal.

The limit of an infinite, homogeneous sample can be obtained in a somewhat simpler fashion by integrating over space before handling any vacuum mode summations. The result is a $\delta$-function selecting the participating vacuum mode and we obtain (see Appendix A.1)

[TABLE]

where we omit factors of $2\pi$ for brevity. It is thus apparent that the cascading signal consists of a real, principal value part and an imaginary part with a $\delta$-function selecting $\omega_{b}^{2}c^{2}=k_{b}^{2}$. Note that $\bm{k}_{b}$ is a linear combination of classical modes interacting with molecule $b$ and $\omega_{b}$ is the corresponding linear combination of frequencies. Both parameters are therefore externally controlled (up to permutation since which particular classical modes interact with molecule $b$ versus $a$ must clearly be summed over) and we obtain a $\pi/2$ phase shift (factor of $i$) only under perfect phase matching $\omega_{b}^{2}=k_{b}^{2}c^{2}$. As demonstrated in the previous section, this phase-shifted component is dominant for finite samples when the phase-matching is such that $|\omega_{b}-k_{b}c|l\ll 1$ where $l$ is the optical path length.

Equation (28) strongly resembles the commonly invoked marcroscopic relation that the cascading field comes as a product of two polarizations $E_{\text{c}}\propto PP$. Our microscopic derivation reveals the precise sense in which they are related, i.e., through integration over the $(\bm{k},\omega)$-space photon Green’s function.

The case of two-molecules is also easily treated and instructive. Arranging our coordinate system such that the two molecules both lie on the $z$-axis, we obtain the simplified form of the integrated Green’s function

[TABLE]

We immediately notice that the Green’s function is now a diagonal tensor. In this two-molecule case, the phase of the cascading signal is sensitively dependent on the distance between molecules and the angle between the cascading beam and the intermolecular axis and comes via the factor $e^{i(\frac{\omega_{b}}{c}r-\bm{k}_{b}\cdot\bm{r})}$. Considering, for example, a co-linear beam geometry perpendicular to the intermolecular axis, we have $\bm{k}_{b}\cdot\bm{r}=0$ and the phase $e^{i\frac{\omega_{b}}{c}r}$ is a sensitive probe of the intermolecular distance. Similarly, with a co-linear beam geometry parallel to the intermolecular axis, we obtain a phase of $e^{i(\frac{\omega_{b}}{c}r-k_{b}r)}=0$ and the long-range cascading signal is in phase with the direct signal. Finally, we note that the long-range term which usually generates cascading does not contribute to the $G^{zz}$ component. This means that, if the dipoles are aligned along the axis connecting the two molecules, no long-range cascading takes place.

IV Conclusions

Cascading is a vacuum-mediated exchange of coherent polarization between two molecules in a sample. Being a coherent process, it generates terms with the same phase-matching and scaling with external field amplitudes as direct signals that are ordinarily of more interest since they reveal higher nonlinearities. In this paper, we have provided a microscopic QED derivation of cascading processes to arbitrary order in the classical modes and connected this to the common macroscopic result obtained via Maxwell’s equations. In particular, we have demonstrated that the $i$ factor used to discriminate cascading from direct signals originates in phase-matching and geometric concerns. In few-molecule samples for example, the cascading signal will generally produce both in-phase and out-of-phase terms of varying dependence on the intermolecular separation vector $\bm{r}$ and cascading wavevector $\bm{k}_{b}$. The uniform integration over the interaction volume cross-section was found to be responsible for this phase-shift, along the same lines as explicated by Feynman for a sheet of dipoles Feynman et al. (2013). A phase-matching condition $|\frac{\omega_{b}}{c}-k_{b}|l\ll 1$, where $l$ is the thickness of the sample, was identified as necessary for the in-phase component of the cascaded signal to vanish. This condition can be achieved in thin films or guaranteed by a co-linear beam geometry. Because they scale quadratically in molecular concentration, the cascading terms often dominate direct nonlinear signals. Additionally, the cascading signal scales linearly with the frequency $\omega_{b}$ of the emitting polarization. This implies that cascading processes will become even more dominant at higher frequencies, such as in x-ray Raman experiments, necessitating a thorough understanding of cascading processes for data analysis purposes Mukamel et al. (2013).

We are now in a position to appreciate the difference between cascading and fluorescence resonant energy transfer (FRET). FRET processes are evaluated by taking the square of the Hilbert space amplitude for molecule $b$ to emit a photon then aborbed by molecule $a$ Craig and Thirunamachandran (1984); Salam (2010b). But this square transition amplitude involves 4 orders in the vacuum mode, and is thus 2 orders higher than cascading. In fact, in a FRET process the emitting molecule populates a photon mode while in cascading the only intermediate photon state is a coherence between the 0- and 1-photon states. Despite these differences, the derivations in Appendix A have much the same flavor as QED derivations of FRET Craig and Thirunamachandran (1984); Andrews (1989); Salam (2010b); Scholes and Andrews (2005), with the principal difference being that cascading depends directly on the dipole coupling tensor rather than its square. The derivation of Appendix A can thus be extended to account for FRET and other higher-order processes. At fourth-order, the same as FRET, we also encounter a 3-body cascading contribution that comes as

[TABLE]

and will, under phase matching, provide terms that scale cubically in the molecular concentration $n$. Higher-order corrections generate $n$-body cascading terms that follow similarly. From the perspective of Feynman diagrams, cascading is like a vertex insertion (with four free branches corresponding to the two in-states and two out-states of the participating molecules) while the Lamb shift is the corresponding self-energy insertion. Both come from vacuum interactions but only the former can be incorporated into a Dyson equation. The $n$-body cascading terms similarly behave like vertex insertions with $2n$ free branches.

An interesting future extension of this work would be to consider manipulating the cascading signal from a system of molecules embedded in an optical cavity, systems which have drawn recent interest Hutchison et al. (2012); Shalabney et al. (2015); Simpkins et al. (2015); Kowalewski et al. (2016a, b); Herrera and Spano (2016). Optical cavities alter the density of electromagnetic field modes from its free-space value, suppressing cascading in all but the cavity mode. In particular, we note that the field mode participating in a cascading process is determined only for infinite samples (when the mode summation collapses to a $\delta$-functino as in Appendix A.1) while, for few-molecule samples, all vacuum modes participate. By strongly coupling a few-molecule sample to an optical cavity, we can effectively force the cascading to occur with a particular field mode (the cavity mode). On the other hand, by suppressing the density of field-modes at the cascading wave-vector for a macroscopic system, the cascading signal can be suppressed for arbitrary sample-sizes and beam geometries. Moreover, molecular coupling to the cavity mode can be tuned by the cavity volume and made much larger than the coupling to the vacuum field. This would allow for enhancement and control of the cavity-mode cascading.

Acknowledgements.

The support of the National Science Foundation (grant CHE-1361516) as well as from the Chemical Sciences, Geosciences, and Biosciences division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy through award No. DE- FG02-04ER15571 is gratefully acknowledged. Support for K.B. was provided by DOE.

Appendix A Microscopic Derivation of Cascading Signals

In this section, we will derive an expression for the cascading signal from a fully microscopic QED perspective while keeping the perturbative order in the external fields completely arbitrary. This will allow us to make very general conclusions without considering particulars of the laser fields.

Expanding equation (6) to 2nd order in the vacuum modes of the electric field

[TABLE]

results in

[TABLE]

where the factor of $i^{2}$ comes from the two orders of expansion in the interaction with the vacuum modes. the $\{\nu\}$ stand for cartesian coordinates coming from the dot products and are implicitly sumed over, and $\langle\dots\rangle_{0}\equiv\text{Tr}[\dots|0\rangle\langle 0|]$ stands for an expectation value taken over the vacuum state. The $-(+)$ subscripts stand for the commutator (anti-commutator) as usual and we have selected the only appropriate contributing terms when performing the initial spatial integrations. In particular, 3-body terms (in which each of the three dipole operators occur at a different molecule) and that in which the action on molecule $b$ preceedes that on $a$ all vanish since $\text{Tr}[\hat{O}_{-}\rho]=0$ for all operators $\hat{O}$ and density matrices $\rho$ (this argument is explained in more detail in section II of Bennett and Mukamel (2014)). Note however, that 3-body cascading does occur when expanding to 4th order in the vacuum modes and cascading 3rd-order processes have been considered in 7th-order nonlinear techniques Cho et al. (2000). Equation (32) is represented diagramatically in Fig. 1.

It is important to note that the dipole expectation values in Eq. (32) includes propagation with respect to the full Hamiltonian (excepting only the vacuum modes). That is, the externally applied semiclassical modes are still included in the propagators and the expression incorporates all orders of interaction between the matter and laser modes. Explicitly, we have

[TABLE]

where the time-dependence of the dipole operators is through $\hat{H}_{0}$, the free material Hamiltonian, and the new interaction Hamiltonian $\hat{H}_{\text{LM}}$ only includes interactions with the laser fields. The dependence on molecule position is notated outside the expectation values since all interactions within each expectation value occur on a particular molecule. Within the dipole approximation, this spatial dependence comes as a simple exponential while the full multi-polar polarization operator has a more general dependence on position.

Equation (32) can be simplified by introducing the polarization of a molecule

[TABLE]

as well as the quantity

[TABLE]

which represents the polarization of a molecule in the laser fields due to a perturbative interaction with the polarization of another molecule. Finally, we identify the photon Green’s function

[TABLE]

which comes as a time-ordered, vacuum expectation value of electric field operators. In terms of these quantities, we may write the cascading signal as

[TABLE]

In appendix B, we simplify the photon Green’s function. Remaining in the time domain for now, we substitute Eq. (61) to obtain

[TABLE]

where the spatial derivatives (the $\nabla$’s) act on $r_{ab}$, the distance between molecules. This expression is very intuitive; the polarization of molecule $b$, evaluated at the retarded time and adjusted by the action of the dipole coupling tensor and factor of $\frac{1}{r}$, is in place of an external field interaction. We may make this identification explicit, writing the effective electric field that is felt by molecule $a$ and caused by the polarization emitted from molecule $b$

[TABLE]

in terms of which the cascading signal is

[TABLE]

Equations (37)-(40) have straightforward physical interpretations. In particular, it is clear that the phase of the cascading signal will depend on the phase of the effective electric field (Eq. (39)) which in turn depends on the geometry of the sample via the summation over positions of molecule $b$. In physical experiments, the material is subjected to a set of lasers with well-defined propagation vectors $\bm{k}_{j}$ and, within the dipole approximation, the real-space dependence of the polarizations will always come as a spatial phase factor $e^{i\sum_{j}^{n}\bm{k}_{j}\cdot\bm{r}}$ where $n$ is the order to which the polarization is expanded with respect to the laser-matter interaction. Thus, the $\bm{k}$-space polarization will have a delta function setting $\bm{k}$ to some linear combination of the incoming wavevectors. In contrast, the spatial structure of the sample can be quite complicated. We thus transform to $\bm{k}$-space and, in the interests of brevity, we omit the integrations over $d\bm{k}_{a}$ $d\bm{k}_{b}$ with the understanding that they will collapse to sums over different choices of $\bm{k}_{a}$, $\bm{k}_{b}$. To completely put all geometric dependence on a single term, all else that is required is to change to the frequency domain with respect to the retarded polarization. We obtain

[TABLE]

where we have used the fact that the detected “signal” field is $E^{*\nu_{s}}_{s}(\bm{r}_{a},t)\to E^{*\nu_{s}}_{s}(\omega_{s})e^{-i(\bm{k}_{s}\cdot\bm{r}_{a}-\omega_{s}t)}$ representing ideal frequency resolution. We explicitly highlight the signal’s dependence on this detected mode rather than continuing to include it implicitly in the set of field parameters $\Lambda$. Additionally, we have substituted the spatiotemporal Fourier transforms of Eqs. (34)-(35). In terms of the discrete Fourier transform of the photon Green’s function

[TABLE]

we have

[TABLE]

or

[TABLE]

where we have substituted the quantity

[TABLE]

that forms the basis of our discussion of the geometric differences between cascading and direct signals in section III. Formally, Eq. (38) in the time domain and Eqs. (43) or (44) in the frequency domain ar our general results. They give the cascading contribution to a heterodyne-detected signal in terms of the molecular polarization and the photon Green’s function, which encodes all geometric information that shapes the cascading signal.

As an aside, we note that one may wish to remain in the time-domain with respect to molecule $a$, while still handling the geometry more completely than in Eq. (38). This can be accomplished by writing

[TABLE]

[TABLE]

where the $\bm{r}_{a}$-dependence of $E_{\text{v}}$ goes away in the homogeneous limit.

In summary, we first derived an intuitive formula (Eq. (38)) for the cascading signal based on the polarization of molecule $b$ at the retarded time acting as a source to interact with molecule $a$. We then obtained a convenient, compact expression (Eq. (44)) for the cascading terms in the heterodyne signal based on the photon Green’s function. It can be readily expanded to any order to obtain cascading corrections to particular nonlinear signals and can also be recast as interaction between an effective electric field and the molecule (Eqs. (39)-(40) and (46)-(47)). The effect of the sample geometry is contained in the exponential phase factor, which approaches the molecule number under good phase matching, and the integral over the photon Green’s function $I(\bm{r}_{a},\bm{k}_{b},\omega_{b})$. Under off-resonant excitation, the phase of the cascading signal is determined solely by this Green’s function. Finally, we note that everything is currently written in terms of discrete summations over molecular positions but we may take the continuum limit by replacing summation by integration $\sum_{\bm{r}}\to\int d\bm{r}n(\bm{r})$ where $n$ denotes molecular concentration.

A.1 Alternative Derivation of the Macroscopic Homogeneous Limit

In Appendix A, we performed all possible simplifications, such as vacuum mode summations, bundling all geometric dependence into the photon Green’s function and a phase-matching exponential factor. Such an approach allows full generality in treating different geometries and provides the necessary ingredients for a simulation of cascading processes in few-molecule samples. The resulting expressions can then be applied directly to microscopic geometries or, after converting summations over molecular locations to spatial integrations, various macroscopic geometries and smoothly taken to the infinite limit by integrating over all space. Some additional insight into this may be obtained by pursuing an alternative derivation in which this spatial integration is performed first, generating $\delta$-functions that determine the participating vacuum mode and collapse the mode sum.

We begin with

[TABLE]

obtained from inserting Eq. (56) into (32) which we then use Eqs. (57) and (9) to rewrite as

[TABLE]

where we have substituted the spatial dependence of the dipole expectation values and the signal field $E_{s}$ as before. Switching to frequency domain, this can be written as

[TABLE]

where we have used the definitions of the frequency-domain polarizations of the previous section to simplify the expression.

In the infinite homogeneous limit, we have $f(\bm{k})\to(2\pi)^{3}n\delta(\bm{k})$

[TABLE]

where we have identified the transverse part of the polarization

[TABLE]

Finally, we can combine the two terms above to obtain Eq. (28)

Appendix B Field Vacuum Expectation Value

In this section, we evaluate the releveant photon Green’s function for this problem defined as a time-ordered, vacuum expectation value of electric field operators

[TABLE]

Ignoring the prefactors for the moment and expanding first the $\hat{\mathcal{E}}_{-}$ operator gives

[TABLE]

where $\rho_{\text{v}}\equiv|0\rangle\langle 0|$ is the vacuum density matrix and we have kept only non-vanishing terms. Due to the cyclic invariance of the trace, we may consider only action from the left with respect to the $\hat{\mathbf{E}}_{+}$ operator that remains. This results in

[TABLE]

The basic commutation relation then gives $\text{Tr}[\hat{a}_{\mathbf{k}^{\prime}_{v},\lambda^{\prime}}\hat{a}^{\dagger}_{\mathbf{k}_{v},\lambda}\rho_{v}]=\text{Tr}[\hat{a}^{\dagger}_{\mathbf{k}^{\prime}_{v},\lambda^{\prime}}\rho_{v}\hat{a}_{\mathbf{k}_{v},\lambda}]=\delta_{\mathbf{k}^{\prime}_{v}\mathbf{k}_{v}}\delta_{\lambda^{\prime}\lambda}$ and we obtain the expression

[TABLE]

To simplify, we perform the polarization sum and change the vacuum mode summation to an integration (see Craig and Thirunamachandran (1984); Salam (2010b))

[TABLE]

after which the integration over solid angle can be carried out

[TABLE]

which then allows us to easily perform the $\omega_{v}$ integration

[TABLE]

Combining this with Eq. (53) yields

[TABLE]

with temporal Fourier transform

[TABLE]

where acting the differential operators results in Eq. (26).

Appendix C Cylindrical Geometry

In this section, we describe the derivation of $I$ for the case of a cylindrical geometry as discussed in Section III, case (iii). The field Green function in this case has a form

[TABLE]

where $G_{M}$ is the field Green function in the complete space $\mathbb{R}^{3}$, and $G_{S}$ is the correction associated with the surface $S$, determined by the condition $z=0$, so that $G$ is the complete Green function in the half-space $z\geq 0$, with the reflecting boundary conditions. The representation of Eq. (63) represents the integral $I$ of Eq. (18)

[TABLE]

We start with computing $I_{M}$; the second integral $I_{S}$ is computed in a similar fashion. Fixing $z$ we compute the integral over the corresponding cross section

[TABLE]

by using inhomogeneous polar coordinates with the center in $(x_{a},y_{a})$ followed by integration over the inhomogeneous radius, resulting in

[TABLE]

The integration is performed in a similar way to how it has been carried out for the $2D$ case, using the identity

[TABLE]

and the fact that the projection of $\bm{k}_{b}$ onto the plane is [math]. For the surface contribution we obtain in a similar way

[TABLE]

By a similar to case (ii) argument, which involves saddle-point approximation for computing a fast oscillating integral, the second contribution in the r.h.s. of Eq. (66) and (68) can be neglected resulting in

[TABLE]

Finally performing integration over $z$ we obtain

[TABLE]

where we have expressed the answer in terms of the sinc function

[TABLE]

with the properties ${\rm sinc}(x)\approx 1$, for $|x|\ll 1$, and $|{\rm sinc}(x)|\leq|x|^{-1}$. In the case $|\omega_{b}-k_{b}|l\ll 1$ of a good phase matching and if $z_{a}$ is not within the wavelength scale from the borders, the first contribution is dominating, so that we reproduce a well-known result

[TABLE]

Appendix D Integration Measure for Inhomogeneous Polar Coordinates

In this section we compute the integration measure for inhomogeneous polar coordinates. Inhomogeneous in this context mans that the radius depends on the angle. Formally, let $\xi(\varphi)$ be a positively-defined function of angle that has a property that the region $\{\bm{r}\in\mathbb{R}^{2}\;|\;r(\bm{r})\leq\xi(\varphi\bm{r})\}$ is convex, where $r(\bm{r})$ and $\varphi(\bm{r})$ are the standard polar coordinates of $\bm{r}$. The inhomogeneous polar coordinates are defined by

[TABLE]

It is easy to demonstrate that the inhomogeneous polar coordinates is a well-defined coordinate system for the region defined by the condition $r<r_{0}$ for any $r_{0}$, provided the convexity condition is satisfied. A reverse statement is also true, i.e., for any connected convex region and any reference point in it one can (and very easily) identify $r_{0}$ and $\xi(\varphi)$ so that the region is defined by the condition $r<r_{0}$ in the inhomogeneous polar coordinate system, centered in the reference point.

To compute the integration measure we first differentiate

[TABLE]

and then apply the concept of the wedge (i.e., antisymmetric) product to obtain

[TABLE]

where in the derivation we have used the antisymmetric properties of the wedge product, namely $dr\wedge dr=d\varphi\wedge d\varphi=0$ and $d\varphi\wedge dr=-dr\wedge d\varphi$. This provides a very simple result for the integration measure

[TABLE]

in particular the integration measure does not contain derivative of the function $\xi(\varphi)$ that describes the inhomogeneity of the coordinate system. We will refer to the above property as the homogeneous property of the inhomogeneous coordinate system.

Appendix E Integration Measure for Inhomogeneous Spherical Coordinates

An inhomogeneous spherical coordinate system is introduced in a similar way:

[TABLE]

Obviously, the inhomogeneous spherical system satisfies the same convexity properties as the inhomogeneous polar system. A much more tedious computation, similar to the one for the polar system, presented in Eq. (75), provides a very simple result

[TABLE]

which shows that the measure is also homogeneous. This means that the homogeneity of an inhomogeneous measure takes place in all dimensions, and there is a general geometrical argument that derives Eq. (78) for all dimensions that bypasses tedious computations.

Appendix F Estimates for the Boundary Contributions Using the Saddle-Point Approximation

We start with the $2D$ case Denoting in Eq. (23)

[TABLE]

we obtain in the saddle-point approximation

[TABLE]

Estimating $\lambda_{0}\sim a\omega_{0}\gg 1$ we see that $I_{1}$ is small compared to $I_{0}$ by a factor $\sim(\sqrt{a\omega_{0}})^{-1}$.