Accumulation effects in modulation spectroscopy with high repetition rate pulses: recursive solution of optical Bloch equations

TL;DR

This paper investigates how accumulation effects influence modulation spectroscopy with high repetition rate pulses, providing a recursive solution to optical Bloch equations to accurately interpret experimental data.

Contribution

It introduces a recursive method to account for accumulation effects in quantum two-level systems driven by phase-modulated pulses, enhancing analysis accuracy.

Findings

Accumulation effects significantly alter harmonic amplitudes in fluorescence signals.

Corrections due to accumulation are essential for precise quantum measurements.

The recursive solution improves the interpretation of modulation spectroscopy data.

Abstract

Application of the phase modulated pulsed light for advance spectroscopic measurements is the area of growing interest. The phase modulation of the light causes modulation of the signal. Separation of the spectral components of the modulations allows to distinguish the contributions of various interaction pathways. The lasers with high repetition rate used in such experiments can lead to appearance of the accumulation effects, which become especially pronounced in systems with long-living excited states. Recently it was shown, that such accumulation effects can be used to evaluate parameters of the dynamical processes in the material. In this work we demonstrate that the accumulation effects are also important in the quantum characteristics measurements provided by modulation spectroscopy. In particular, we consider a model of quantum two-level system driven by a train of phase-modulated light pulses, organised in analogy with the 2D spectroscopy experiments. We evaluate the harmonics' amplitudes in the fluorescent signal and calculate corrections appearing from the accumulation effects. We show that the corrections can be significant and have to be taken into account at analysis of experimental data.

Figures

Equations

$$\mathcal{F}\left[P(t)\right](\nu )=\sum _{n=-\infty }^{\infty }{e}^{\mathrm{i}\nu t_{0}n}{\int }_0^{t_0}P(t-t_{0}n|P_n)e^{-\mathrm{i}\nu t}dt.$$

$$P_n=\kappa (P_{n-1},I_n;K).$$

$$\Omega (t)=\frac{V}{\Delta }\sum _{n=-\infty }^{\infty }\sum _{i=1,2,3,4}{\Theta }_{\Delta }(t-t_i-t_{0}n)\cos ({\omega }_{i}t).$$

$$\frac{dR}{dt}=MR-\frac{\Gamma}{2}\left(\begin{array}[]{c}0\\ 0\\ 1\end{array}\right),\quad M=\left(\begin{array}[]{ccc}-\gamma&-\omega_{0}&0\\ \omega_{0}&-\gamma&2\Omega\\ 0&-2\Omega&-\Gamma\end{array}\right).$$

$${\gamma }^{-1}\sim {\tau }_1\sim T\sim {\tau }_2\ll {\Gamma }^{-1}\sim t_1\sim t_0\ll \frac{2\pi }{{\varphi }_0}\sim \frac{2\pi }{{\varphi }_1}\sim \frac{2\pi }{{\varphi }_2}\sim \frac{2\pi }{{\varphi }_3}.$$

$$R(t)=U(t,t_{\mathrm{init}})\left[R(t_{\mathrm{init}})-\frac{\Gamma }{2}{\int }_{t_{\mathrm{init}}}^t{U}^{-1}(t^{\prime },t_{\mathrm{init}})\Xi d{t}^{\prime }\right],$$

$$R_1^{(n)}$$

$$R_2^{(n)}$$

$$R_3^{(n)}$$

$$R_4^{(n)}$$

$$PL(\nu)=\int_{-\infty}^{+\infty}\left(Z(t)+\frac{1}{2}\right)e^{-\mathrm{i}\nu t}dt\\ =\sum_{n=-\infty}^{+\infty}e^{-\mathrm{i}\nu t_{0}n}\int_{0}^{t_{0}}\left(Z\big{(}t^{\prime}+t_{0}n\big{)}+\frac{1}{2}\right)e^{-\mathrm{i}\nu t^{\prime}}dt^{\prime}$$

$$\begin{gathered} PL(\nu )\approx \sum _{n=-\infty }^{+\infty }{e}^{-\mathrm{i}\nu t_{0}n}{\int }_0^{2t_1\approx t_0}\left(Z_4^{(n)}+\frac{1}{2}\right)e^{-\Gamma t-\mathrm{i}\nu t}dt \\ \approx \frac{1}{\Gamma +\mathrm{i}\nu }\sum _{n=-\infty }^{+\infty }{e}^{-\mathrm{i}\nu t_{0}n} \\ \times \left(\left(Z_4^{(n)}+\frac{1}{2}\right)-\left(Z_4^{(n-1)}+\frac{1}{2}\right)e^{-2\Gamma t_1}\right) \end{gathered}$$

$$Z_{4\mathrm{steady}}^{(n)}=-\frac{1}{2}+Z_{\mathrm{pr}}^{(n)}{e}^{-\Gamma t_1}\left(1-Z_{\mathrm{nu}}^{(n)}{e}^{-2\Gamma t_1}\right)+\mathcal{O}(e^{-5\Gamma t_1}).$$

$$\begin{gathered} PL(\nu )\approx \frac{1-e^{-2\Gamma t_1}}{\Gamma +\mathrm{i}\nu }\sum _n{e}^{-\mathrm{i}\nu t_{0}n}{Z}_{pr}^{(n)}\left(1-Z_{nu}^{(n)}{e}^{-2\Gamma t_1}\right) \\ =\sum _{\mathbf{m}}\delta (\mathbf{m}\cdot \mathbf{\varphi }-\nu )A_{\mathbf{m}}, \end{gathered}$$

$$A_{\mathbf{m}}=\frac{(1-e^{-2\Gamma t_1}){\sin }^{2}V}{4(\Gamma +\mathrm{i}\nu )t_0}\left({\mathcal{B}}_{\mathbf{m}}+e^{-2\Gamma t_1}{\sin }^{2}V\delta {\mathcal{B}}_{\mathbf{m}}\right).$$

$$A_{1,-1,-1,1}$$

$$A_{1,1,-1,-1}$$

$$U_{0}(t)=\left(\begin{array}[]{ccc}e^{-\gamma t}\cos\omega_{0}t&-e^{-\gamma t}\sin\omega_{0}t&0\\ e^{-\gamma t}\sin\omega_{0}t&e^{-\gamma t}\cos\omega_{0}t&0\\ 0&0&e^{-\Gamma t}\end{array}\right).$$

$$\frac{dU_{I}(t)}{dt}=2\Omega(t)U_{0}^{-1}(t)\left(\begin{array}[]{ccc}0&0&0\\ 0&0&1\\ 0&-1&0\end{array}\right)U_{0}(t)U_{I}(t).$$

$$U^{(n)}_{i}=U_{0}(\tilde{t}_{i}-\tilde{t}_{i-1})\\ \times\left[\mathds{1}+\left(\begin{array}[]{ccc}-2\sin^{2}f_{i}\sin^{2}\frac{V}{2}&-\sin 2f_{i}\sin^{2}\frac{V}{2}&e^{-(\Gamma-\gamma)(t_{i}-\tilde{t}_{i-1})}\sin f_{i}\sin V\\ -\sin 2f_{i}\sin^{2}\frac{V}{2}&-2\cos^{2}f_{i}\sin^{2}\frac{V}{2}&e^{-(\Gamma-\gamma)(t_{i}-\tilde{t}_{i-1})}\cos f_{i}\sin V\\ -e^{(\Gamma-\gamma)(t_{i}-\tilde{t}_{i-1})}\sin f_{i}\sin V&-e^{(\Gamma-\gamma)(t_{i}-\tilde{t}_{i-1})}\cos f_{i}\sin V&-2\sin^{2}\frac{V}{2}\end{array}\right)\right];$$

$$F^{(n)}_{i}=\frac{1}{2}\left(\begin{array}[]{c}e^{-\gamma(\tilde{t}_{i}-2\tilde{t}_{i-1}+t_{i})}(1-e^{-\Gamma(t_{i}-\tilde{t}_{i-1})})\sin V\sin(f_{i}-\omega_{0}(\tilde{t}_{i}-\tilde{t}_{i-1}))\\ e^{-\gamma(\tilde{t}_{i}-2\tilde{t}_{i-1}+t_{i})}(1-e^{-\Gamma(t_{i}-\tilde{t}_{i-1})})\sin V\cos(f_{i}-\omega_{0}(\tilde{t}_{i}-\tilde{t}_{i-1}))\\ e^{-\Gamma(\tilde{t}_{i}-2\tilde{t}_{i-1}+t_{i})}-1+(1-e^{\Gamma(\tilde{t}_{i}-2\tilde{t}_{i-1}+t_{i})})\cos V\end{array}\right),$$

$$Z_{\mathrm{pr}}^{(n)}$$

$$Z_{\mathrm{nu}}^{(n)}$$

$$\delta {\mathcal{B}}_{1,-1,-1,1}$$

$$\delta {\mathcal{B}}_{1,1,-1,-1}$$

Full text

Introduction. High repetition rate, tens of MHz, laser systems produce trains of short pulses at relatively low energy. Still, because the pulse can be very short ($<$ 10 fs), the field strength during the pulses is sufficient for inducing nonlinear optical effects BoydBook . In addition, if the studied material does not reach the ground state equilibrium during the time interval between the pulses ($\sim$ 10 ns) accumulation effects can occur. In recent years, application of the phase modulated pulsed light for the purposes of advanced spectroscopic measurements became an area of growing interest our ; KWSMLPM2014 ; LBSE2017 ; KKMP2016 ; BBS2015 ; TDM2003 ; TLM2007 ; BBS2017 ; M2016 . In this technique the phase modulation of light pulses generates harmonics in the measured signal allowing extraction of information about the various light-induced coherent and dissipative dynamics. For the systems with long characteristic life-time of the excited states the accumulation effects may become essential. Accounting of such effect as well as the question of new information, which can be obtained from it, is an important theoretical and practical problem.

The accumulation effects occurring due to irradiation of matter by the impulsive lasers can have several manifestations, such as heat accumulation TPV2009 or particle shielding NJSH2014 . In this article we discuss somewhat more delicate experiments addressing the spectroscopic measurements of matter characteristics. In the latter context the accumulation effects were discussed in Ref. our , where dynamics of the system was modelled by a (classical) rate equations. There the effects were described theoretically and verified experimentally. In the experiment the light pulses are organised in two collinear beams and arrive to the sample at the same time. The laser frequency is chosen in such a way that the system is excited by the two-photon absorption. The base optical frequencies of the two beams are modulated by acoustic frequencies $\phi_{1}$ and $\phi_{2}$, respectively. Due to the non-linear absorption, the effective intensity, $I_{n}$, of the light absorbed from the $n$th pair of pulses oscillates with $n$ on the frequencies $\phi=\left|\phi_{1}-\phi_{2}\right|$ and $2\phi$, it is $I_{n}\propto\cos(\phi t_{0}n)+\frac{1}{4}\cos(2\phi t_{0}n)$, where $t_{0}$ is the time-interval between the pulses (in practice $t_{0}\sim$10ns, which is around $10^{3}$ times less then the modulation period, $2\pi/\phi$). As the result, the population of the system excited state oscillates in time as well and generates, in turn, an oscillating signal. It can be seen from the following consideration. The signal (fluorescence) is proportional to the population of the excited state $P(t)$, i.e. fraction of the total number of molecules that have been excited. As it has been shown in our the Fourier transform of the signal is given by a series of integrals. Each integral is the Fourier transform of the population between the neighbouring pulses,

[TABLE]

Behaviour of the population $P(t-t_{0}n|P_{n})$ after the $n$th pulse parametrically depends on the population $P_{n}$, which is the population taken at the instance of time right before the $n$th pulse. One can show, that $P_{n}$ satisfies a non-linear recurrence

[TABLE]

The function $\kappa$ is structurally determined by the type of relaxation kinetics $K$ and includes the effective intensity, $I_{n}$ absorbed by the system from the $n$th pulse. Indeed, as soon as the population of the excited state taken right before the $n$th pulse is $P_{n-1}$, the light pulse can excite $(1-P_{n-1})$ fraction of molecules only, so that the initial condition for solution of kinetic equations, denoted by letter $K$, is $P_{n-1}$, while $I_{n}$ enters as a parameter. Being solved in the time-interval starting right before the $n$th pulse and ending right before the $(n+1)$th pulse the kinetic equations give the value $P_{n}$. The dependence of $P_{n}$ on $K$, $P_{n-1}$ and $I_{n}$ is encoded in the formula (2). For large $n$ the solution of the above recurrence (2) does not depend on the way the system was prepared, this solution determines the ”dynamical” steady-state of the system. Obviously, in case of fast decay of the excited state, meaning that $P_{n-1}=0$, the population as well as the generated fluorescence oscillates according to $I_{n}$. Therefore, the harmonics appearing in the signal are determined by $I_{n}$ only. Detection of the harmonics and their relative amplitudes can give information on the light-matter interaction pathways, while periodic repetition of the process allows to improve the signal to noise ratio for the harmonics’ amplitudes, i.e. amplitudes of the peaks in the frequency domain $\nu$ of the signal.

Below, we call the excited state generated by the laser pulses with a given $n$ a primary state until the next $(n+1)$th pulse comes. The part of the primary excited state, which remains even after the $(n+1)$th pulse, we call accumulated state. The modulated frequencies of the signal component due to the primary and accumulated excited states follow different peaks. It was proposed in Ref. our , that the cumulative effects can be taken into account, but also allow to extract information of the dynamical processes in the material. Indeed, when the time interval between pulses is short and the system does not relax to its equilibrium completely, the population $P_{n-1}\neq 0$, the non-linearity of function $\kappa$ causes generation of new harmonics and bring changes into the relative amplitudes of the original harmonics (harmonics which are contained in $I_{n}$). Thus, measuring of the harmonics’ relative amplitudes in the signal gives us information on $K$, as far as other experimental parameters are known. In particular, it was shown, that comparison of the linear kinetics and the one with quadratic term, which are causing the spontaneous decay in a molecular system ($\dot{P}(t)=-\Gamma P(t)$, where $P(t)$ is the population of the excited state) and the band-to-band annihilation of hole-electron pairs in a semicoductor ($\dot{P}(t)=-gP^{2}(t)$, where $P(t)$ is the concentration of free carriers) respectively, have a well defined distinct signatures. Thus, the idea to utilise the frequency modulated laser pulses for optical measurements of dynamical processes in various media seems promising, while accounting of the emerging cumulative effects requires revision of our theoretical models.

In the present work we extend the theoretical approach developed in Ref. our for the analysis of the experimental scheme in fluorescence detected 2D spectroscopy. In the conventional photon echo 2D spectroscopy 3 noncollinear laser pulses generate coherent signal in phase matched directions. The signal is mixed with the 4th one, so called local oscillator pulse, which allows phase sensitive field detection J2003 ; HZ2011 . Analysis of the signal coming out with a given combination of $k$’s for various choices of the time delays between the pulses gives information of quantum dynamics and dissipation of the system. In the modulation spectroscopy approach, all four beams are collinear, but the base frequencies are modulated by the acoustic frequencies $\phi_{1}$, $\phi_{2}$, $\phi_{3}$, $\phi_{4}$. As the result, the signal (fluorescence or photocurrent) is expected to contain harmonics at frequencies ${\bm{m}}\cdot{\bm{\phi}}=m_{1}\phi_{1}+m_{2}\phi_{2}+m_{3}\phi_{3}+m_{4}\phi_{4}$, ($m_{i}$ are integers, positive or negative, subject to the constrain $m_{1}+m_{2}+m_{3}+m_{4}=0$). The corresponding relative signal amplitudes should give information of the internal parameters of the system. The spatial separation of signals in the original 2D scheme is replaced by filtering out the proper Fourier component of the signal. A special attention in our research is paid to the description of the accumulation effects, i.e. the case of a system with a long life-time excited state.

Our analysis is based on solution of a system of Bloch equations for an atom driven by a train of phase-modulated light pulses. Contrary to Ref. our , where the two-photon absorption process has been involved, we consider a near-resonance processes. As we demonstrate below, the theoretical estimation of fluorescence signal contains the sought harmonics. We estimate their amplitudes and corrections to the amplitudes appearing due to the accumulation effects in the case of long-living excitation and analyse nontrivial features in the structure of the obtained harmonics.

Model of two-level system driven by frequency modulated pulsed field. The Hamiltonian of two-level system interacting with light can be written in terms of Pauli matrices Pauli , it has the form $\mathcal{H}=\frac{\hbar}{2}\omega_{0}\sigma_{3}-\hbar\Omega(t)\sigma_{1}$. The $\omega_{0}$ is the transition energy and $\Omega(t)$ describes interaction with the time-dependant external field,

[TABLE]

Here $n$ counts the quad-pulses (we use the notion quad-pulse as a single word for the neighbouring four pulses, see figure 1), while summation over $i$ runs over four single pulses within the single quad-pulse, $\omega_{i}=\omega+\phi_{i}$ and $V$ is the coupling of the system with the field. The function $\Theta_{\Delta}(t)$ describes the envelope of the single pulse with the pulse duration $\Delta$.

The standard theoretical approach for prediction of 2D spectra is based on the time-dependant perturbation theory Mbook developed for the solution of density matrix evolution equation and is a useful tool for bookkeeping of diagrams representing various quantum pathways of light-matter interaction. This theory, while being very efficient in for calculation of the response functions in spectroscopy, is hard to be used for derivation of the solution describing the accumulation effects. To this end, the language of Bloch equation is more convenient. In this language various diagrams are counted automatically by multiplication of the evolution matrices. Nevertheless, the language of ladder diagrams is used on the figure 4.

Evolution of the Bloch vector $R$ with components $X$, $Y$ and $Z$ obeys the matrix equation AEbook

[TABLE]

Below for the constant vector in (4) we use the notation $\Xi$, i.e. $\Xi=\left(0,0,1\right)^{T}$. Remind, that the components of the vector $R=(X,Y,Z)^{T}$ are connected with the components of the density matrix, $\rho$, such that Pauli $X=\frac{1}{2}\textrm{Tr}\;\rho\sigma_{1}$, $Y=\frac{1}{2}\textrm{Tr}\;\rho\sigma_{2}$, and $Z=\frac{1}{2}\textrm{Tr}\;\rho\sigma_{3}$. The component $Z$ is equal to the difference of the diagonal elements, $(\rho_{ee}-\rho_{gg})/2$, and, since $\textrm{Tr}\;\rho=1$, the value $Z=-1/2$ corresponds to the system in the ground state ($\rho_{gg}=1$), while the $Z+\frac{1}{2}$ is the probability to find system in the excited state. The constants $\gamma$ and $\Gamma$ describe decay rate of coherence (the off-diagonal elements of the density matrix) and of the excited state due to fluorescence and the non-radiative processes, respectively. In the physically relevant picture $\gamma$ is large, while the other parameter, $\Gamma$, is small. The relation between these two parameters means that the fluorescence is generated mainly between the quad-pulses and is proportional to the exponent $e^{-2\Gamma t_{1}}$, while the decoherence effects are essential on the time scales comparable with the time intervals between the pulses, see figure 1 top.

In our model we employ the usual approximations for the 2D spectra analysis. We assume the time ordering of pulses and semi-impulsive limit Mbook . It means that the pulses are well separated and short compared to any time scale (including the ones corresponding to the detuning $\epsilon=\omega_{0}-\omega$ and the coupling $V$ of the system with the field), but long compared to the oscillation period of the light field, $\omega^{-1},\omega_{0}^{-1}\ll\Delta$. Notations for the time-intervals used in the text are shown on the figure 1 bottom. The time-scales are ordered in the following way:

[TABLE]

In experiments often the time interval $T$ is much longer then the intervals $\tau_{1}$, $\tau_{2}$, but this fact do not influence on our further consideration, as soon as we assume that the fluorescence during interval $T$ is negligible.

The solution of the equation (4) at time $t$, for the given initial condition $R(t_{\mathrm{init}})$ taken at time $t_{\mathrm{init}}$, can be generically represented in the form PY1990

[TABLE]

where $U(t,t_{\mathrm{init}})$ is an orthogonal 3$\times$3 matrix, which satisfies the homogeneous matrix equation $\frac{d}{dt}U=MU$, subject to some initial condition $U(t_{\mathrm{init}},t_{\mathrm{init}})$. In case of zero field, $\Omega(t)=0$, the matrix $U(t,t_{\mathrm{init}})$ is a matrix of free evolution of the system, $U_{0}(t-t_{\mathrm{init}})$, see equation (17) in the appendix. Since the pulses are well separated, the electric field of light between the pulses is zero, $\Omega(t_{\mathrm{init}})=0$, so that $U(t_{\mathrm{init}},t_{\mathrm{init}})=\mathds{1}$. Therefore, one can seek for solution around each single pulse and glue them up. The whole calculation scheme can be presented in the following way

[TABLE]

The vectors $R_{i}^{(n)}$ are the vectors obtained after evolution around a single pulse, R_{i}^{(n)}=R\big{(}(n-1)t_{0}+\tilde{t}_{i}\big{)} ($\tilde{t}_{i}$ is the time after the $i$th pulse where $\Omega(\tilde{t}_{i})=0$, i.e. $\tilde{t}_{1}=t_{1}+\tau_{1}/2$, $\tilde{t}_{2}=t_{2}+T/2$, $\tilde{t}_{3}=t_{3}+\tau_{2}/2$, and $\tilde{t}_{4}=t_{0}$), which, in turn, serves as initial vector for the next stage of evolution. The evolution matrix, $U((n-1)t_{0}+t^{\prime},(n-1)t_{0}+\tilde{t}_{i-1})$ calculated in the local time frame, $\tilde{t}_{i-1}<t^{\prime}<\tilde{t}_{i}$, generate matrices $U^{(n)}_{i}\equiv U((n-1)t_{0}+\tilde{t}_{i},(n-1)t_{0}+\tilde{t}_{i-1})$. Note that the $n$-dependence of $U^{(n)}_{i}$ arises due to the global time dependence of the excitation field phase in $\Omega(t)$, eq. (3). The vector $F^{(n)}_{i}$ originates from the second term in the solution (5). Within the above assumptions the problem can be solved analytically. Explicit forms of $U^{(n)}_{i}$ and $F^{(n)}_{i}$, and comparison with numerics are given in appendix, see equations (19) and (20) and the text around.

For the purposes of the present investigation we are interested in certain Fourier components of the fluorescent signal. Such analysis can be done by analogy with the one provided in Ref. our . Fluorescence is proportional to the population of the excited state, $Z(t)+1/2$. As soon as we are not interested in the way the system was prepared the Fourier transform of the fluorescent signal, $PL(\nu)$, can be presented as a series of Fourier integrals calculated in the vicinity of each quad-pulse, such that the initial instant is in $-\infty$:

[TABLE]

The main contribution to the fluorescence comes form the longest time interval $2t_{1}$, where the population decays exponentially with the rate $\Gamma$, namely the time behaviour of $Z(t)$ in this time interval is given by $\left(Z_{4}^{(n)}+\frac{1}{2}\right)e^{-\Gamma t}$, where $Z_{4}^{(n)}$ is $Z$-component of the vector $R_{4}^{(n)}$ calculated from equation (9). The function $PL(\nu)$ after this assumption can be approximated as

[TABLE]

Therefore, $Z^{(n)}_{4}$ is the only component of the solution we are interested in. The form of expression (11) for the fluorescent signal reproduces the one obtained in ref. our for the case of linear kinetics. Thus, the result (11) has a universal nature and reflects the exponential decay of the excited state.

Calculation of the fluorescent signal in case of long-living excited state. In our analysis we assume that (i) $t_{1}\gg 1/\gamma$ and the quantum coherence does not contribute to the formation of the accumulated state; (ii) the oscillations attributed by the frequencies $\phi_{i}$ are slow, $\phi_{i}t_{0}\ll 1$, so that the nearest quad-pulses generate identical excitations of the system; (iii) the excited state life-time $1/\Gamma$ is comparable with $t_{0}\simeq 2t_{1}$, but less then $2t_{0}$, so that $e^{-\Gamma t_{1}}$ is a small parameter of the model. Therefore, to find the accumulated state solution it is enough to require equivalence of $Z$ components of the Bloch vector before and after the quad-pulse, $Z_{4}^{(n-1)}=Z_{4}^{(n)}$. Indeed, $X$ and $Y$ components decay faster then the time interval between the quad-pulses, $t_{1}$, so in the initial condition for the second quad-pulse they can be put to zero and the initial vector $R_{4}^{(n-1)}$ is proportional to $\Xi$. Population of the accumulated state is the same before and after the quad-pulse, see figure 1.

From the set of vector equations (6)-(9), which becomes closed by setting $R_{4}^{(n-1)}=Z_{4}^{(n)}\Xi$, one can find the solution for $Z_{4\,\mathrm{steady}}^{(n)}$. To do this we solve the system (6)-(9) perturbatively. The solution is an expansion over the small factor $e^{-\Gamma t_{1}}$:

[TABLE]

The primary contribution $Z_{\mathrm{pr}}^{(n)}$ is the one obtained from interaction of the quad-pulse with two-level system, which initially is in the ground state. Formally, $-\frac{1}{2}+e^{-\Gamma t_{1}}Z_{\mathrm{pr}}^{(n)}$ is $Z$ component of $R_{4}^{(n)}$ obtained from (6)-(9) at $R_{4}^{(n-1)}$ set to $-\frac{1}{2}\Xi$. The second term in (12) describes ”interference“ of the newly generated excitation with the primary state. The part $Z_{\mathrm{nu}}^{(n)}$ can be found from the solution of the problem (6)-(9) with $R_{4}^{(n-1)}=z\,\Xi$. Then $Z$ component of $R_{4}^{(n)}$ becomes $Z_{4}^{(n)}=-\frac{1}{2}+Z_{\mathrm{pr}}^{(n)}e^{-\Gamma t_{1}}+\left(z+\frac{1}{2}\right)Z_{\mathrm{nu}}^{(n)}e^{-2\Gamma t_{1}}$. Setting $z=-\frac{1}{2}$ we return back to the primary contribution. To find the accumulated state one has to take $z$ to be again equal to $Z_{4}^{(n)}$ and solve the resulting linear equation. Keeping the lower order terms, up to the $e^{-3\Gamma t_{1}}$, we arrive to the formula (12). Note, that the higher order terms, which are denoted by $\mathcal{O}(e^{-5\Gamma t_{1}})$ come from accounting of three and more quad-pulses. The explicit forms of $Z_{\mathrm{pr}}^{(n)}$ and $Z_{\mathrm{nu}}^{(n)}$ were derived by method of functional programming, their expressions can be obtained by the formulas (21) and (22) and the tables given in the appendix.

Having calculated the $Z$ component of the accumulated state one can find the fluorescence by the formula (11). The laser pulse leads the system through the coherence to the fully excited state, in which the fluorescent signal is generated. After interaction with the $i$th pulse the components of Bloch vector get a term with the phase factor $e^{\mathrm{i}m_{i}\phi_{i}t_{0}(n-1)}$, since the $i$th carrier frequency from the $n$th quad-pulse has a phase shift $\phi_{i}t_{0}(n-1)$. The value of $m_{i}$ is determined by the effective number of interactions with the $i$th pulse in the primary contribution and by the sum of the effective numbers of interactions from the nearest quad-pulses for the accumulation part. Eventually $Z$ component contains a sum of terms distinguished by the phases $e^{\mathrm{i}(\bm{m}\cdot\bm{\phi})t_{0}(n-1)}$. Summation in the formula (11) can be easily implemented using the formula $\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}e^{\mathrm{i}n\alpha}=\lim_{N\to\infty}\frac{\sin(N+1/2)\alpha}{2\pi\sin\alpha/2}=\delta(\alpha)$. The resulting expression for the fluorescence after substitution of (12) into (11) has the form

[TABLE]

with the amplitudes

[TABLE]

The sum in (13) runs over the integer-valued vectors $\bm{m}=(m_{1},m_{2},m_{3},m_{4})$ subject to the constraint $m_{1}+m_{2}+m_{3}+m_{4}=0$.The correction $\delta\mathcal{B}_{\bm{m}}$ can be obtained by the formula $\sin^{2}V\delta\mathcal{B}_{\bm{m}}=\sum_{\tilde{\bm{m}},\tilde{\bm{m}}^{\prime}}\mathcal{B}_{\tilde{\bm{m}}}\left(\mathcal{B}_{\tilde{\bm{m}}^{\prime}}+\mathcal{C}_{\tilde{\bm{m}}^{\prime}}\right)\delta_{\tilde{\bm{m}}+\tilde{\bm{m}}^{\prime}-\bm{m}},$ here $\delta_{\bm{m}}$ is the Kronecker symbol. Explicit formulas for the complex-valued coefficients $\mathcal{B}_{\bm{m}}$ and $\mathcal{C}_{\bm{m}}$ are given in the appendix.

Positions of the peaks, amplitudes and corrections. Each $\delta$-function in (13) is responsible for one peak on the Fourier transform plot, figure 2. Within the assumptions made in the model, there are 13 primary peaks (orange lines on the figure) generated by each quad-pulse and 49 secondary peaks generated by the accumulated steady-state. This numbers do not include the peaks appearing due to the symmetry of the coefficients $\mathcal{B}_{-\bm{m}}=\mathcal{B}^{*}_{\bm{m}}$, $\delta\mathcal{B}_{-\bm{m}}=\delta\mathcal{B}^{*}_{\bm{m}}$. These peaks are positioned on the negative half-line of the $\nu$ axis symmetrically to the ones plotted on the picture (figure 2) and have complex conjugated amplitudes. One more peak appears at the origin, $\nu=0$, and is not shown.

Positions of the peaks in the frequency domain $\nu$ (see eq. (13)) form an unstructured set, see figure 2. There is no regular ordering of the peaks. To find the indexing vector $\bm{m}$ one has to calculate the values $(\bm{m}\cdot\bm{\phi})$ and find the appropriate value of $\nu$. The strength of amplitudes is mainly determined by the order of interaction with the quad-pulse, $V^{2}$, $V^{4}$ or $V^{6}$. The minimal order of interaction can be calculated as the sum $\left|m_{1}\right|+\left|m_{2}\right|+\left|m_{3}\right|+\left|m_{4}\right|$. The peak amplitudes also include exponential factors $e^{-\gamma(\tilde{t}_{i}-\tilde{t}_{i-1})}$ and $e^{-\Gamma(\tilde{t}_{i}-\tilde{t}_{i-1})}$, i.e. they depend on the time intervals between the pulses.

The set of the peaks becomes structured in another representation, which is obtained after taking the Fourier transformation of the expression (13) for the fluorescence over the parameters $\tau_{1}$ and $\tau_{2}$ (or $T$). In coordinates $(\omega_{\tau_{1}},\omega_{\tau_{2}})$ (or $(\omega_{\tau_{1}},\omega_{T})$) it forms some regular lattice-like structure, see figure 3. The peaks are positioned around the multiples of the detuning frequency $\epsilon$ and arranged with respect to their indices $\bm{m}$. The shift from the exact values of the detuning is due to the factor $\bm{m}\cdot\bm{\phi}$ coming additively to the $\epsilon$. The primary peaks (orange dots on the plot, fig 2) are placed around six points $\omega_{\tau_{1}},\left|\omega_{\tau_{2}}\right|,\left|\omega_{T}\right|\lesssim\epsilon$, while the secondary peaks can have positions around $\pm 2\epsilon$. The latter peaks is a results of effective multiplication of similar or identical interaction processes taken place in the neighbouring quad-pulses.

Two important remarks has to be done here. Presentation in coordinates $(\omega_{\tau_{1}},\omega_{\tau_{2}})$ allows us to group the signals systematically. While the detuning frequency $\epsilon=\omega-\omega_{0}$ has a clear definition, in case of short spectrally broad laser pulses, it cannot be experimentally determined. Therefore the plots on the figure 3 is a convenient way to represent all peaks at once, but not the real picture that can be observed in experiment. The second remark, is that keeping of the constant phase differences in the long series of quad-pulses is important for observation of harmonics in the fluorescent signal. Randomisation of the phase would lead to destruction of the obtained picture. Indeed, the summation described in the text above the equation (13) with any randomisation of phases would result in the sum \frac{1}{2\pi}\sum_{n=-N}^{N}e^{\mathrm{i}\big{(}(\bm{m}\cdot\bm{\phi})-\nu\big{)}t_{0}n+\mathrm{i}\alpha_{n}} with $\alpha_{n}$ randomly chosen for each $n$. In the limiting case of statistically independent random $\alpha_{n}$ distributed uniformly from $-\pi$ to $\pi$ the Fourier transform would result in one peak at $\nu=0$ as $N$ goes to infinity, while in the intermediate cases for distribution of $\alpha_{n}$ one expects AHN2002 some distribution of peak positions around their canonical values given by $(\bm{m}\cdot\bm{\phi})$, which can be sensitive to the particular choice of $N$.

Compare now the amplitudes of primary peaks and corrections to them from accumulated state. As an example, we consider the most interesting peak ${\bm{m}}=(1,-1,-1,1)$ and the newly formed peak ${\bm{m}}=(1,1,-1,-1)$. The exact expressions for the corrections are given in the appendix, eqs. (23), (24). It is instructive to compare the peak amplitude in the leading orders over $V$ and $e^{-2\Gamma t_{1}}$. From the eqs. (23), (24) we obtain

[TABLE]

As one can see the correction to $A_{1,-1,-1,1}$, the last term in the brackets in eq. (15), as well as the amplitude $A_{1,1,-1,-1}$ has the same order of magnitude as the amplitude of the primary peak and only the factor $e^{-\Gamma t_{1}}$ can significantly suppress the above corrections. The ratio of the amplitudes $A_{1,1,-1,-1}/A_{1,-1,-1,1}\approx 2e^{-2\Gamma t_{1}-2\gamma T}$ gives us an access to the parameter $\gamma$. Note, that the combination of phases $\phi_{1}+\phi_{2}-\phi_{3}-\phi_{4}$ should correspond to the signal attributed in the literature as the double quantum coherence signal MOY2007 , which can appear in the three or more level atom due to formation of the coherence between the ground state and the upper excited level. In present theory the signal is the result of accumulation. Namely, the signal amplitude is a product of two amplitudes with the indexes $(1,0,0,-1)$ and $(0,1,-1,0)$. Since they both have the same order of $V^{2}$ this explains the strength of the total signal proportional to $V^{4}$. Some more explanations regarding the structure of eqs. (15) and (16) are given in the caption of the figure 4.

Conclusion. In this work we studied the model of two-level system driven by light coming in quad-pulses with modulated frequencies. Our analysis was based on analytic solution of the system of the optical Bloch equations under standard assumptions used in analysis of spectroscopic problems. It was shown that the Fourier transform of the fluorescent signal contains a number of oscillating modes with the frequencies ${\bm{m}}\cdot{\bm{\phi}}$. Amplitudes of the primary peaks, originated from each single quad-pulse, get corrections when one takes into account the accumulation effects. Moreover, these effects lead to formation of new (secondary) peaks. The corrections as well as amplitudes of the secondary peaks are of the same order of magnitude as the amplitudes of the primary peaks, the only suppressing factor is the exponential decay of the excited state from one quad-pulse to the other. Therefore, the cumulative effects have to be taken into account in the analysis of phase-modulated harmonic light spectroscopy data. Note, also that the ratio of the amplitudes of primary and secondary peaks can be used for revealing of the internal parameters of the system, see for instance eqs. (15) and (16).

There are other effects that can change the peaks’ amplitudes, which have not been considered above. First, our consideration is based on certain physical assumptions, such as semi-impulsive limit and perfect time ordering of pulses. Accounting of effects beyond these simplifications can bring additional corrections into the amplitude values. More important, however, is the following observation. The peaks obtained in our derivation, based on Bloch equation, should be considered as a two-level system response on light excitation averaged over ensemble of many identical systems. In the work Ref. M2016 it is noticed that generation of higher harmonics in the signal can be caused by many-particle effects. In fact, the same many-particle arguments lay in the basis of our classical consideration of the modulation spectroscopy in Ref. our . Therefore, the corrections to the primary peaks can have various sources. This part of the problem certainly requires further theoretical efforts to formulate a complete physical picture of the interplay between the accumulation and many-particle phenomena.

Acknowledgements. Authors thanks Khadga J. Karki and Shaul Mukamel for useful discussions.

Appendix A

Explicit forms of matrices used in the paper. Matrix of free evolution of the system has the form

[TABLE]

It is solution of the equation $\frac{dU_{0}(t)}{dt}=\tilde{M}U_{0}(t)$, where $\tilde{M}$ is matrix $M$ from eq. (4) with $\Omega=0$. To find matrix $U(t)$ we make the substitution $U(t)=U_{0}(t)U_{I}(t)$ and consider the problem in interaction picture PY1990

[TABLE]

Under assumptions made in the paragraph around the equation (5) the matrix equation (18) can be solved approximately in the vicinity of each pulse to obtain expressions for $U^{(n)}_{i}$ and $F^{(n)}_{i}$. They are

[TABLE]

[TABLE]

where we use the notations $f_{i}=\omega_{0}(t_{i}-\tilde{t}_{i-1})-\omega_{i}t_{i}-\omega_{i}t_{0}(n-1)$.

To check the obtained analytic approximation we compared it with a numerical solution for a single pulse, see figures below.

a)    b)

On the figure a) we compared $X$ and $Z$ components of the numeric solution and of the one obtained within our approximations for Gaussian choice of the pulse envelop, $\Theta_{\Delta=2\sqrt{\pi}\sigma}(t)=e^{-\frac{(t-2)^{2}}{4\sigma^{2}}}$ and the frequency $\omega_{0}\simeq 285\Delta$, such that there are around six periods of the wave lie within the pulse duration. The detuning $\epsilon=2.82\Delta$, meaning that system is in resonance. Time dependence of the analytic solution is obtained by replacing $V$ in equations (19) and (20) with $\frac{V}{\Delta}\int_{-\infty}^{t}\Theta_{\Delta}(t^{\prime})dt^{\prime}$. On the figure b) the value of relative error $\left|\frac{Z-Znum}{Z}\right|$ is plotted. The discrepancy reaches not more then 2%. Note that the further one from the resonance the larger becomes the error.

The non-zero values of coefficients $\mathcal{B}_{\bm{m}}$ in eq. (13) are given in the table below,

[TABLE]

Expression for the solution (12) can be obtained from the above table by using the formulas (c.c. means complex conjugation)

[TABLE]

supplemented by the table of non-zero coefficients $\mathcal{C}_{\bm{m}}$:

[TABLE]

Explicit form of corrections to the coefficients $\mathcal{B}_{1,-1,-1,1}$ and $\mathcal{B}_{1,1,-1,-1}$ from the non-trivial steady state

[TABLE]