Absorption Spectrum of a Two-Level System Subjected to a Periodic Pulse Sequence
H. F. Fotso V. V. Dobrovitski





Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA
Department of Physics, University at Albany (SUNY),  Albany, New York 12222, USA

Ames Laboratory US DOE, Ames, Iowa, 50011, USA
QuTech and Kavli Institute of Nanoscience, TU Delft, Lorentzweg 1, 2628 CJ Delft, the Netherlands


We investigate how the quantum control of a two-level system (TLS) coupled to photons can modify and tune the TLS's photon absorption spectrum. Tuning and controlling the emission and the absorption is of much interest e.g. for the development of efficient interfaces between stationary and flying qubits in modern architectures for quantum computation and quantum communication. We consider the periodic pulse control, where the TLS is subjected to a periodic sequence of the near-resonant Rabi driving pulses, each pulse implementing a 180^∘ rotation. For small inter-pulse delays, the absorption spectrum features a pronounced peak of stimulated emission at the pulse frequency, as well as equidistant satellite peaks with smaller spectral weights. As long as the detuning between the carrier frequency of the driving and the TLS transition frequency remains moderate, this spectral shape shows little change. Therefore, the quantum control allows shifting the absorption peak to a desired position, and locks the absorption peak to the carrier frequency of the driving pulses. Detailed description of the spectrum, and its evolution as a function time, the inter-pulse spacing and the detuning, is presented. 




Absorption Spectrum of a Two-Level System Subjected to a Periodic Pulse Sequence.
    V. V. Dobrovitski
    December 30, 2023
=================================================================================





§ INTRODUCTION

An interface between stationary and flying qubits, that enables a long-range entanglement between different quantum network nodes, is essential for quantum information processing <cit.>. It is of particular importance for the solid state qubits, such as quantum dots or color centers <cit.>, which can be efficiently coupled to each other via photons and thus employed for quantum communications and distributed quantum information processing. However, the slow fluctuations in the environment of the solid-state qubits (e.g. the local strain and/or the local electric fields) constitute a lingering challenge, because they unpredictably shift the optical transition frequency of the qubits <cit.>. This slow drift of the transition frequency (spectral diffusion) makes it difficult to achieve the precise matching between the photons originating from different qubits that is required for efficient entanglement. To mitigate the spectral diffusion problem, various methods have been proposed and successfully used <cit.>, focusing primarily on the tuning of the emission spectrum and on improving the indistinguishability of the photons emitted from different qubits. In particular, it has been recently suggested <cit.> that the application of a periodic sequence of the optical control pulses to a quantum emitter (a two-level system coupled to the electromagnetic radiation bath) can re-direct most of the emission into a peak located at a preset target frequency (determined by the carrier frequency of the pulse driving field), and therefore greatly improve the  indistinguishability of the photons coming from different emitters.

At the same time, there is a growing interest, accompanied by impressive progress <cit.> in the long-range entanglement schemes based on the photon absorption, and the theoretical developments which allow control and tuning of the absorption spectra have become timely and interesting. Correspondingly, a question arises whether the absorption-based entanglement can also be improved using the pulse control of the emitters, i.e. whether the absorption spectrum of a two-level system (TLS) coupled to the radiation bath can be modified and tuned by the control pulses. Besides, the studies of absorption of a TLS subjected to an external control are of fundamental interest due to the intimate connection between emission and absorption <cit.>. For instance, if the TLS is continuously driven by a strong coherent laser field then the TLS emission spectrum has an interesting three-peak structure, with two additional side peaks located at the frequencies ±Ω_R (where Ω_R is the laser Rabi driving frequency), and the absorption spectrum of the same system also acquires additional structure, displaying the regions of gain, corresponding to an amplification of the probing weak field instead of attenuation <cit.>. 



The emission spectrum of the pulse-controlled TLS exhibits similarities with the continuously driven TLS emission: it has a central peak at the carrier frequency of the pulses ω_0, as well as the  satellite peaks at ω_0 ±π/τ,  ± 2π/τ,  ⋯, where τ is the inter-pulse distance. Thus, it is reasonable to expect that absorption also can be controlled with the periodic pulses, and that the resulting absorption spectrum also has non-trivial features. In this work we study the absorption spectrum of a TLS driven by a periodic sequence of optical π-pulses, and examine its dependence on the pulse sequence period and the detuning of the emitter with respect to the pulse frequency (Fig.<ref>). We show that both expectations above are correct, and therefore the pulse control indeed can be a useful tool for controlling and tuning the absorption spectrum of a TLS. The absorption spectrum has a pronounced peak of stimulated emission at the carrier frequency of the pulses, and equidistant satellite peaks with smaller spectral weights. The qualitative features of this absorption spectrum do not change much as long as the detuning between the carrier frequency of the driving pulses and the TLS transition frequency remains moderate. Therefore, we show that the optical control enables creation of pairs of quantum nodes (one node working as an emitter and the other as an absorber) with precisely matching frequencies, and therefore greatly increased entanglement efficiency. This approach can also be used to improve the coupling of the emitters and the absorbers to the optical cavities, since the laser pulses can tune both the emission and the absorption lines of the respective quantum nodes, bringing them in the resonance with the respective cavities, and stabilizing both the emission and the absorption peaks at the desired location.

The rest of the paper is organized as follows. In Sec. <ref> we describe the model of the two-level system coupled to the photon bath and controlled by the pulses, the master equations governing the system dynamics, and the two methods, analytical and numerical, used for calculating the absorption spectrum. In Sec. <ref> we present analytical and numerical results demonstrating the control and tunability of the absorption spectrum. In Sec. <ref> we present conclusions.



§ MODEL OF THE TWO-LEVEL SYSTEM COUPLED TO THE ELECTROMAGNETIC RADIATION BATH



We model the quantum emitter as a TLS with the ground state |g⟩ and the excited state |e⟩, separated in energy by E_e - E_g = ħω_1; below we set ħ=1. Initially, at time t=0, the excited state is occupied and the ground state is empty. The TLS is coupled to a photon bath, and is periodically driven by pulses of the laser field with the Rabi frequency Ω. Within the rotating-wave approximation (RWA) <cit.>, in the reference frame rotating at frequency ω_0, the system in question is described by the Hamiltonian 

    H = ∑_kω_k a^†_ka_k   +   Δ/2σ_z - i ∑_k g_k( a^†_kσ_- - a_kσ_+ ) 
       +   Ω_x(t)/2(σ_+ + σ_-),

where Δ=ω_1-ω_0 is the detuning of the TLS's transition frequency from the carrier frequency of the pulses; here we introduced the standard pseudo-spin Pauli operators for the TLS, namely σ_z = |e⟩⟨ e| - |g⟩⟨ g|, σ_+ = |e⟩⟨ g| and σ_- = |g⟩⟨ e| = (σ_+)^†. Furthermore, a^†_k and a_k are respectively the creation and the annihilation operator for a photon of mode k with the frequency ω_k, and g_k is the strength of coupling to the TLS. Note that in the rotating frame all frequencies are measured from the pulse carrier frequency ω_0, so that the zero frequency in the rotating frame corresponds to ω_0 in the lab frame; we take it as the target frequency for our TLS. 

The time-dependent driving Ω_x(t) in Eq. (<ref>) represents the control pulses; here we consider the simple situation of the square-shaped pulses, with Ω_x(t)=Ω during the pulses and zero otherwise. In fact, below we assume that the pulses are almost instantaneous, i.e. that Ω is much larger than all other relevant energy scales, and that each pulse performs an almost instantaneous 180^∘ rotation of the TLS around the x-axis, interchanging |e⟩ and |g⟩; this assumption will be discussed further below. In the absence of control (Ω_x(t)≡ 0), the system exhibits spontaneous decay, and the corresponding emission rate is Γ = 2π∫ g_k^2  δ(ω_k - Δ)   dk; we normalize our energy and time units so that Γ = 2, and the corresponding spontaneous emission line has a simple Lorentzian shape 1/(ω^2 + 1), with the half-width equal to 1. The absorption spectrum is defined here in a standard way, as the energy absorbed by the TLS from a weak probing field of frequency ω. The probing field is assumed to be weak enough that it does not significantly affect the population of each state <cit.>; our goal is to calculate the absorption as a function of frequency and time.

To understand the dynamics of the system, we analyze the time evolution of the density matrix of the emitter, which is written as

    ρ(t)    =   ρ_ee(t) |e⟩⟨ e| + ρ_eg(t) |e⟩⟨ g| 
       +   ρ_ge(t) |g ⟩⟨ e|  + ρ_gg(t) |g⟩⟨ g|   ,

with ρ_ge^* = ρ_eg. For the TLS described by the above Hamiltonian (<ref>), within the Markovian approximation, the density matrix operator is governed by the master  equations <cit.> in the rotating-wave approximation:

    ρ̇_ee   = i Ω_x(t)/2(ρ_eg - ρ_ge) - Γρ_ee  ,   
    
    ρ̇_gg   = -i Ω_x(t)/2(ρ_eg - ρ_ge) + Γρ_ee  ,  
    
    ρ̇_ge   = (iΔ -Γ/2) ρ_ge -i Ω_x(t)/2(ρ_ee - ρ_gg)   , 
    
    ρ̇_eg   = (-iΔ -Γ/2) ρ_eg +i Ω_x(t)/2(ρ_ee - ρ_gg) .

Since the pulse driving is assumed to be strong and short (Ω≫Δ, Γ), the pulses can be considered as instantaneous. Each of them inverts the populations of the excited and ground state and swaps the values of ρ_eg and ρ_ge, i.e. 

    ρ(nτ + 0) = σ_x ρ(nτ - 0) σ_x

where ρ(nτ - 0) and ρ(nτ + 0) are the density matrices immediately before and after the pulse, correspondingly, with n being an integer and τ the period of the pulse sequence; in other words, the pulses interchange ρ_ee with ρ_gg, and ρ_eg with ρ_ge. 

We want to determine the energy absorbed from a weak probing field by the TLS subjected to the periodic sequence of the π-pulses. Since the effect of the probing field is small, the absorption spectrum can be calculated within the linear response theory, so that at long times T the absorbed energy Q(ω) is given by <cit.>

    Q(ω)    =    2 A^2 
       ×   Re{∫_0^T  dt ∫_0^T-t dθ ⟨[ σ_-(t) , σ_+(t+θ) ] ⟩e^-i ωθ},

where [O_1, O_2] is the commutator of the operators O_1 and O_2, and the angled brackets represent the expectation values evaluated in the absence of the probing field. σ_-(t) and  σ_+(t+θ) are the time-dependent operators in the Heisenberg representation, and the expectation values are taken with respect to the initial state of the two-level system (in our case, fully occupied excited state and empty ground state). The constant A is independent of the pulse parameters, and does not affect the spectral shape, determining only the absolute scale of the absorption. The expression (<ref>) can be rewritten as

    Q(ω)    =    2 A^2 Re{𝒫_2(ω) - 𝒫_1(ω) }
       =    P_2(ω) - P_1(ω)

where

    𝒫_2(ω) =  ∫_0^T   dt ∫_0^T-t dθ ⟨σ_-(t) σ_+(t+θ) ⟩e^-i ωθ

and

    𝒫_1(ω) =  ∫_0^T   dt ∫_0^T-t dθ ⟨σ_+(t+θ) σ_-(t) ⟩e^-i ωθ

The term P_1(ω)=2 A^2 Re{𝒫_1(ω) } can be viewed as the direct emission of the two-level system and P_2(ω)=2 A^2 Re{𝒫_2(ω) } as the direct absorption so that the difference yields the net absorption <cit.>. We evaluate the terms P_1(ω) and P_2(ω) separately, and obtain the total absorption spectrum Q(ω) by taking the difference.



To evaluate the emission spectrum, it is convenient to re-express the two-time correlation function ⟨σ_+(t+θ) σ_-(t) ⟩ as a single-time expectation value <cit.>, according to

    ⟨σ_+(t+θ) σ_-(t) ⟩ = 
             =Tr[ ρ(0) U^-1(0,t+θ) σ_+ U(0,t+θ)U^-1(0,t) σ_- U(0,t)] 
             = Tr[ σ_- ρ(t) U^-1(t,t+θ) σ_+ U(t,t+θ) ] 
             = Tr[ ρ'(t,t + θ) σ_+ ]





where σ_+ and σ_- are the time-independent Pauli operators in the Schrödinger representation, and U(t_1,t_2) is the evolution operator of the emitter from time t_1 to time t_2, as  determined by the master equations  (<ref>). The calculations are simplified by introducing the matrix ρ'(t,s); its initial value at s=t is ρ'(t,t)=σ_-ρ(t), and its further evolution from s=t to s=t+θ is governed by the emitter's evolution operator U(t,t+θ), so that ρ'(t,t+θ)=U(t,t+θ)ρ'(t,t) U^-1(t,t+θ). In this way the evaluation of the two-time correlators becomes rather straightforward (although lengthy, see Appendix for details), and the function P_1(ω) can be obtained analytically and/or numerically. In order to calculate the function P_2(ω), we use the same procedure, simplifying the two-time correlation function as

    ⟨σ_-(t) σ_+(t+θ) ⟩ = Tr[ σ_+ ρ”(t,t + θ) ],

by introducing the matrix ρ”(t) = ρ(t)σ_-, whose time evolution is also governed by U(t,t+θ), i.e. ρ”(t,t+θ)=U(t,t+θ) ρ”(t,t) U^-1(t,t+θ). Note that ρ' and ρ” are not density matrices, and the symmetries of the proper density matrix ρ(t) (such as ρ_gg=1-ρ_ee and/or ρ_ge^* =ρ_eg) are not applicable to ρ' and ρ”.



In the absence of the pulses, the absorption spectrum has a Lorentzian-shaped profile centered at the emitter's frequency that equals to the detuning Δ (Fig.<ref>). In the presence of the pulses, we calculated the absorption spectrum both analytically and numerically by iteratively evolving the density matrix operator between successive pulses on a discrete time grid, using the equations of motion (<ref>), with the initial conditions ρ_ee = 1, ρ_eg=ρ_ge=ρ_gg=0, and then making use of (<ref>) and (<ref>) to calculate the two-time correlation functions.





 §.§ Numerical solution



To find the solution numerically, we divide the time axis in the intervals of length τ (equal to the inter-pulse separation), and each interval between the pulses is further discretized into smaller steps of length Δ t. The goal is to find the two-time correlators ⟨σ_-(t) σ_+(t+θ) ⟩ and ⟨σ_-(t) σ_+(t+θ) ⟩ for each value of t and θ on this time grid, and use Fourier transform to find P_1(ω) and P_2(ω), whose difference gives the absorption spectrum Q(ω). 

We start at t=0 with the known initial conditions for ρ(t), and use Eqs. <ref> to evolve all elements of the density matrix ρ(t) from time t to t+Δ t, and repeat this integration up to t = τ. Then the π-pulse is applied to the system, transforming the density matrix in accordance with Eq. (<ref>), and the iterative integration is resumed to propagate the density matrix from t=τ to t=2τ, until another pulse is applied at 2τ. The process is repeated until time T = N_pτ is reached, where N_p is the total number of pulses. 
In this way we can obtain the elements of ρ'(t,t) and ρ”(t,t) for every t∈[0,N_pτ]. Then, for each time t we propagate the matrices ρ' and ρ” from time t to time t+θ by solving the master equations (<ref>); the values ρ'(t,t) and ρ”(t,t) serve as initial conditions. As a result, we obtain ρ'(t,t+θ) and ρ”(t,t+θ) for all values of θ∈ [0,T-t]. This procedure produces the two-time correlators ⟨σ_+(t+θ) σ_-(t)⟩ and ⟨σ_-(t) σ_+(t+θ)⟩, see Eqs. (<ref>) and (<ref>). Finally, Fourier transform with respect to θ and integration over t give us 𝒫_1(ω) and 𝒫_2(ω), thus determining the absorption spectrum Q(ω).




 §.§ Analytical solution



The analytical solution for the density matrix evolution between the pulses can be obtained directly from Eqs. <ref>, and combined with the analytically calculated transformation of the density matrix by pulses as described by Eq. (<ref>), thus providing a fully analytical solution to the problem. The corresponding calculation is quite lengthy, and is presented in detail in the Appendix. In the limit of long T (i.e. large number of pulses N_p), the resulting expression for 𝒫_1(ω) is


    𝒫_1(ω)    =   1/(1 + e^-Γτ)γ_0[ ( 1-e^-Γτ/Γ - e^-γ_0 τe^γ_2 τ-1/γ_2
     + e^γ_2 τ-1/γ_2 1-e^-γ_0 τ/e^2 γ_1 τ-1) 
    
    (N_p + e^-Γτ/1 + e^-Γτ) 
       -   e^γ_2 τ-1/γ_2 1-e^-γ_0 τ/e^2 γ_1 τ-1 (  
    2 e^-N_p γ_1 τ - 1 /e^-2 γ_1 τ - 1 + (e^-Γτ - e^ -2 Γτ)e^-N_p γ_1 τ/e^-2 γ_1 τ - e^-2 Γτ) 
    
    
    
    ]


We have also performed the similar calculation for 𝒫_2(ω), expressing it as  𝒫_2(ω)=𝒫_3(ω)-𝒫_1(ω), and the resulting expression for 𝒫_3(ω) in the long-time limit is

    𝒫_3 =  N_pτ/γ_0 -  N_p/γ_0^2(1 - e^-γ_0 τ)
    + e^γ_0 τ+e^-γ_0 τ-2/γ_0^2(e^2γ_1τ-1) [N_p - 2/1-e^-2γ_1τ]


where γ_0 = i(ω -Δ ) + Γ/2, γ_1 = iω  + Γ/2, and γ_2 = i(ω -Δ ) - Γ/2. 



§ RESULTS



Fig. <ref> shows the absorption spectrum obtained using both analytical and numerical approaches for a two-level system with Δ = 3.0 and a pulse sequence with τ = 0.2 after N_p=8 pulses. The panels <ref>(a) and <ref>(b) show P_1(ω) and P_2(ω), respectively, while the panel <ref>(c) shows the absorption spectrum obtained by taking their difference according to Eq. <ref>. In the presence of the pulse control we see the main peak in the absorption spectrum at the carrier frequency of the pulses (ω = 0 in the rotating frame), and the satellite peaks at the multiples of ±π/τ. A good agreement is clearly seen, despite the large number of pulse limit used in the analytical result and the numerical Fourier transform of the finite time step data in the numerical result. The agreement is further improved by considering the spectrum at longer times, see Appendix. These results provide clear validation of the tools used in these studies. Moreover, note that the absorption spectrum, as given by Eq.(<ref>), is the difference of two terms of comparable magnitude in a broad frequency range. As a result, it is influenced by numerical errors, but the small discrepancy between the analytical and the numerical results shows that this kind of errors is not critical. Thus, the numerical solution can be used in the future studies of more complex driving protocols, which may not be amenable to an analytical solution. 

 

Fig.<ref> shows the analytic results for the time evolution of the absorption spectrum for Δ = 3.0 and τ = 0.2. The snapshots of the spectrum are presented after N_p=8  (black), N_p=12 (red), N_p=16 (green), and N_p=20 (blue) pulses. The absorption spectra feature a positive part and a negative part. The latter corresponds to the stimulated emission and the former to the ”true absorption” <cit.>. The central peak corresponds to the stimulated emission at the pulse frequency and satellite peaks at multiples of ±π/τ with amplitudes that decrease away from the central frequency and are strongly suppressed at large frequencies. The lineshape is established early, and the amplitude of the peaks increases with time.



In Fig.<ref> we present the dependence of the absorption spectrum for Δ = 3.0 on the period τ of the pulse sequence. The absorption spectrum is shown after 8 pulses for τ = 0.2, τ = 0.3, τ = 0.4, and τ = 0.5. The satellite peaks move closer to the central peak and their relative amplitude increases as τ becomes longer. Also note the increase in the positive fraction of spectral weight with increasing τ.

It is also interesting to study the dependence of the absorption spectrum on the detuning Δ. The corresponding results are shown in Fig.<ref> which presents the spectrum under a pulse sequence of period τ = 0.2 for the detuning values of Δ = 3.0(black), Δ = 4.0(red), Δ = 5.0(green), and Δ = 6.0(blue) after 8 pulses. The lineshape remains almost the same for all four values of the detuning parameter when τ is kept constant. In fact, we observe that the lineshape of the absorption spectrum shows little dependence on Δ as long as Δ·τ≲ 1.  The figure shows that the fraction of the spectral weight contained in the positive-frequency satellites (with ω>ω_0) slightly increases with Δ, while the spectral weight of the negative-frequency satellites (ω<ω_0) correspondingly decreases.



§ CONCLUSIONS



We have studied the absorption spectrum of a two-level system driven by a periodic sequence of π-pulses. This absorption spectrum is determined by the energy absorbed from a probing field weak enough to not significantly affect the population of the excited and ground state. We have solved the problem by integrating the master equation analytically and numerically and obtained from both methods results that are in excellent agreement. Our results show that for moderate values of Δ·τ, the absorption spectrum has a lineshape with little dependence on Δ. It has a pronounced peak of stimulated emission at a the pulse frequency along with satellite peaks at multiples of ±π/τ away from this frequency. The weights of these satellite peaks are strongly suppressed away from the central peak. 

By using the optical control considered in this work (with, possibly, more complex pulse protocols), it is possible to create pairs of quantum nodes, with one node working as an emitter and the other as an absorber, with precisely matching frequencies, and therefore greatly increased entanglement efficiency. In a similar manner, one can use it to improve the coupling of the emitters and absorbers to the optical cavities, using the laser pulses to tune both the emission and the absorption lines of the respective quantum nodes, bringing them in the resonance with the respective cavities, and stabilizing both the emission and the absorption peaks at the desired location.




§ ACKNOWLEDGMENTS
 

We thank D. D. Awschalom and M. E. Flatté for helpful discussions. Work at the Ames Laboratory was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Division of Materials Sciences and Engineering. The Ames Laboratory is operated for the US Department of Energy by Iowa State University under Contract No. DE-AC02-07CH11358. This work was partially supported by AFOSR MURI program. 


99



JeffKimble_qtmInternet H. J. Kimble, Nature 453, 1023 (2008).







Imamoglu_Awschalom_QDOT_PRL_99 A. Imamoǧlu, D. D. Awschalom, G. Burkard, D. P. DiVincenzo, D. Loss, M. Sherwin and A. Small, Phys. Rev. Lett. 83, 4204 (1999). 

ChildressRepeater L. Childress, J. M. Taylor, A. S. Sørensen, and M. D. Lukin, Phys. Rev. Lett. 96, 070504 (2006).

HansonAwschalom_QIP_ss_08 R. Hanson and D. D. Awschalom, Nature 453 1043 (2008).














Bernien_Hanson_Nature2013  H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen,
L. Childress and R. Hanson, Nature 497, 86 (2013).

Pfaff_Hanson_Science2014 W. Pfaff, B. J. Hensen, H. Bernien, S. B. van Dam, M. S. Blok, T. H. Taminiau, M. J. Tiggelman, R. N. Schouten,
 M. Markham, D. J. Twitchen, R. Hanson, Science 345 6196, 532 (2014).

Gao_Imamoglu_NatComm2013 W. B. Gao, P. Fallahi, E. Togan, A. Delteil, Y.S. Chin, J. Miguel-Sanchez and A. Imamoglu, Nature Com. 4:2744 (2013) .

Hanson_loopholeFree_Nature2015 B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb,	M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau and R. Hanson, Nature 526, 682 (2015).

Basset_Awschalom_PRL2011 L. C. Basset, F. J. Heremans, C. G. Yale, B. B. Buckley, and D. D. Awschalom, Phys. Rev. Lett. 107, 266403 (2011).
FaraonEtalNatPhot A. Faraon, P. E. Barclay, C. Santori, K.-M. C. Fu, and R. G. Beausoleil, Nature Photonics 5, 301 (2011)

NV_Review_PhysRep2013 M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, L. C. L. Hollenberg, Physics Reports 528, 1 (2013).


Sipahigil_SiV A. Sipahigil, K. D. Jahnke, L. J. Rogers, T. Teraji, J. Isoya, A. S. Zibrov, F. Jelezko, and M. D. Lukin, 
Phys. Rev. Lett. 113, 113602 (2014).

Rogers_SiV L. J. Rogers, K. D. Jahnke, M. H. Metsch, A. Sipahigil, J. M. Binder, T. Teraji, H. Sumiya, J. Isoya, M. D. Lukin, P. Hemmer, and F. Jelezko, 
Phys. Rev. Lett. 113, 263602 (2014).




SantoriVuckovicYamamoto C. Santori, D. Fattal, J. Vučković, G. S. Solomon, and Y. Yamamoto, Nature 419, 594 (2002).
CarterGammonQDcavity S. G. Carter, T. M. Sweeney, M. Kim, C. S. Kim, D. Solenov, S. E. Economou, T. L. Reinecke, L. Yang, A. S. Bracker, and D. Gammon, Nat. Photonics 7, 329 (2013).

Fu_Beausoleil_PRL2009 K.-M. C. Fu, C. Santori, P. E. Barclay, L. J. Rogers, N. B. Manson, and R. G. Beausoleil, Phys. Rev. Lett. 103, 256404 (2009).

AmbroseMoerner_spectralDiffusion_Nature1991 W. P. Ambrose and W. E. Moerner, Nature 349, 225 (1991).

Hansom_Atature_APL2014 J. Hansom, C. H. H. Schulte, C. Matthiesen, M. J. Stanley, and M. Atatüre,  Appl. Phys. Lett. 105, 172107 (2014).

Acosta_Beausoleil_PRL2012 V. M. Acosta, C. Santori, A. Faraon, Z. Huang, K.-M. C. Fu, A. Stacey, D. A. Simpson, K. Ganesan, S. Tomljenovic-Hanic, A. D. Greentree, S. Prawer, and R. G. Beausoleil, Phys. Rev. Lett. 108, 206401 (2012).

Kuhlmann_Warburton_NatPhys2013 A. V. Kuhlmann, J. Houel,	A. Ludwig, L. Greuter, D. Reuter, A. D. Wieck, M. Poggio, and  R. J. Warburton, Nat. Phys. 9, 570 (2013).

Crooker_Bayer_PRL2010 S. A. Crooker, J. Brandt, C. Sandfort, A. Greilich, D. R. Yakovlev, D. Reuter, A. D. Wieck, and M. Bayer, Phys. Rev. Lett. 104, 036601 (2010).

Matthiesen_Atature_SciRep2014 C. Matthiesen, M. J. Stanley, M. Hugues, E. Clarke, and M. Atatüre, Sci. Rep. 4, 4911 (2014).

FotsoEtal_PRL2016 H. F. Fotso, A. E. Feiguin, D. D. Awschalom and V. V. Dobrovitski, Phys. Rev. Lett. 116, 033603 (2016).

LenhardEtal_PRA2015 A. Lenhard, M. Bock, C. Becher, S. Kucera, J. Brito, P. Eich, P. Müller, and J. Eschner, Phys. Rev. A 92, 063827 (2015).
TrautmannAlber_PRA2015 N. Trautmann and G. Alber, Phys. Rev. A 93, 053807 (2016).
YangWrachtrupEtal_NatPhot2016 S. Yang, Y. Wang, D. D. Bhaktavatsala Rao, T. H. Tran, A. S. Momenzadeh, M. Markham, D. J. Twitchen, P. Wang, W. Yang, R. Stöhr, P. Neumann, H. Kosaka, and J.  Wrachtrup, Nat. Photonics 10, 507 (2016).

einstein_phys_Z_1917 A. Einstein, Phys. Z. 18, 121 (1917).

Mollow_PhysRevA_5_1972 B. R. Mollow, Phys. Rev. A 3, 2217 (1972).

Wu_et_al_Mollow_PRL1977 F. Y. Wu, S. Ezekiel, M. Ducloy and B. R. Mollow, Phys. Rev. Lett. 38, 1077 (1977). 

RF_Mollow_PhysRev1969 B. R. Mollow, Phys. Rev. 188, 1969 (1969).

Autler_Townes_PhysRev1955 S. H. Autler and C. H. Townes, Phys. Rev. 100, 703 (1955).

KnightMilonni_PhysRep1980 P. L. Knight and P. W. Milonni, Phys. Rep. 66, 21 (1980).

RF_Heitler_Book1960 W. Heitler, The quantum Theory of Radiation (Oxford University Press, London, 1960, Third Ed.).

Cohen_Tannoudji_Book1992 C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Atom-Photon Interactions, Basic Processes and Applications (John Wiley & Sons, Inc., New York, 1992).

Mollow_PhysRevA_3_1972 B.R. Mollow, Phys. Rev. A 5, 1522 (1972).

BloomMargenau_PhysRev1953 S. Bloom and H. Margenau, Phys. Rev. 90, 791 (1953).



Scully_Zubairy_book1997 M.O. Scully and M.S. Zubairy, Quantum Optics (Cambridge University Press, New York 1997).

Loudon_book1983 R. Loudon The Quantum Theory of Light (Clarendon Press, Oxford, 1983).








§ APPENDIX

Here we present the details of the analytical calculation of the absorption spectrum for a two-level system subjected to a periodic sequence of control pulses. As shown in Refs. Mollow_PhysRevA_5_1972,BloomMargenau_PhysRev1953, the absorption spectrum can be determined from the two-time correlation functions of the TLS:

    Q(ω)    =    2A^2 ×
           Re{∫_0^T  dt ∫_0^T-t dθ ⟨[ σ_-(t) , σ_+(t+θ) ] ⟩e^-i ωθ},

where [ , ] is the commutator of the two enclosed operators, and the angled brackets represent the expectation values evaluated in the absence of the probing field. A is a proportionality constant. This can be rewritten as:

    Q(ω)    =    2 A^2 Re{𝒫_2(ω) - 𝒫_1(ω) }
       =    P_2(ω) - P_1(ω),

with

    𝒫_2(ω) =  ∫_0^T   dt ∫_0^T-t dθ ⟨σ_-(t) σ_+(t+θ) ⟩e^-i ωθ

and

    𝒫_1(ω) =  ∫_0^T   dt ∫_0^T-t dθ ⟨σ_+(t+θ) σ_-(t) ⟩e^-i ωθ.

The terms P_1(ω)=2 A^2 Re{𝒫_1(ω) } and P_2(ω)=2 A^2 Re{𝒫_2(ω) } can be evaluated separately and the absorption spectrum obtained by taking the difference. To find 𝒫_2(ω), we express the correlation function as


    ⟨σ_-(t) σ_+(t+θ) ⟩   =   Tr[ ρ(0) U^-1(0, t) σ_- U(0, t) U^-1(0, t+θ) σ_+ U(0, t+θ)] 
       =   Tr[ σ_+  U(t, t+θ) U(0, t) ρ(0) U^-1(0, t) σ_- U(0, t)U^-1(0, t) U^†(t, t+θ) ] 
       =   Tr[ σ_+ U(t, t+θ) ρ(t)σ_-  U^-1(t, t+θ)] 
       =   Tr[ σ_+ U(t, t+θ) ρ”(t,t) U^-1(t, t+θ)] 
       =   Tr[ σ_+ ρ”(t,t+θ)]


where σ_+ and σ_- are the Pauli operators, and U(t_1,t_2) is the operator of the emitter's evolution from t_1 to t_2, as determined by the master equations (<ref>). The subsequent calculations are facilitated by introducing the matrix ρ”(t,s); its initial value at s=t is defined as ρ”(t,t)=ρ(t)σ_-, and its further evolution from s=t to s=t+θ is governed by the emitter's evolution operator U(t,t+θ), so that ρ”(t,t+θ)=U(t,t+θ)ρ”(t,t) U^-1(t,t+θ). 

It is informative to write ρ”(t,s) explicitly as

    ρ”(t,s) = [ ρ”_ee(t,s) ρ”_eg(t,s); ρ”_ge(t,s) ρ”_gg(t,s) ]

so that

    σ_+ ρ^''(t + θ) = [   ρ”_ge(t,t+ θ) ρ^''_gg(t,t+ θ);               0               0 ],

and the corresponding two-time correlation function is obtained directly as

    ⟨σ_-(t)σ_+(t+θ)⟩   =   Tr[ σ_+ ρ”(t,t+θ)] 
       =   ρ”_ge(t,t+ θ).

The initial condition for ρ”, corresponding to s=t, has a form

    ρ”(t,t) = ρ(t)σ_- = [ ρ_eg(t)       0; ρ_gg(t)       0 ],

being determined by the elements of the “true” density matrix ρ_eg(t)=⟨ e|ρ(t)|g⟩ and ρ_gg(t)=⟨ g|ρ(t)|g⟩, see Eq. <ref>. 
Similarly, for ρ'(t,s) the initial condition at s=t are 

    ρ'_ee(t,t)   =   ρ'_eg(t,t)=0,
    ρ'_gg(t,t)   =   ρ_eg(t),  ρ'_ge(t,t)=ρ_ee(t),

and the corresponding two-time correlator is 

    ⟨σ_+(t+θ)σ_-(t)⟩   =   Tr[ ρ'(t,t+θ) σ_+ ] 
       =   ρ'_ge(t,t+θ).

Therefore, our task is reduced to to determining ρ”_ge(t,t+θ) and ρ'_ge(t,t+θ).

The master equations characterizing the time evolution of the TLS density matrix are given by Eqs. <ref> and Eq. <ref>; the time development of the matrices ρ”(t,s) and ρ'(t,s) also obeys these equations of motion as s increases from t to t+θ. Specifically, when s corresponds to the time interval between the pulses, we have

    d/ds ρ”_ee(t,s)    =    -Γ ρ”_ee(t,s) 
    d/ds ρ”_gg(t,s)    =   Γ ρ”_ee(t,s) 
    d/ds ρ”_ge(t,s)    =   (iΔ - Γ/2) ρ”_ge(t,s) 
    d/ds ρ”_eg(t,s)    =   (-iΔ - Γ/2) ρ”_eg(t,s)

for any value of the parameter t; the same equations govern the dynamics of ρ'. The effect of the pulses on ρ” and ρ' is also easily derived from Eq. <ref>: when t+s coincides with the time of the pulse application, i.e. when s=nτ for some integer n, the matrix transforms as

    ρ”(t,nτ + 0) = σ_x ρ”(t,nτ - 0) σ_x

where ρ”(t,nτ-0) and ρ”(t,nτ+0) are the matrices immediately before and after the pulse, correspondingly; in other words, each pulse interchanges ρ”_ee with ρ”_gg, and ρ”_eg with ρ”_ge; the transformation of ρ' is the same.

Let us start with establishing the initial condition for ρ”(t,s) at s=t, which is determined by ρ_gg(t) and ρ_eg(t), see Eq. <ref>. First, we note that ρ_eg(t)≡ 0. Indeed, the initial condition at t=0 for the density matrix ρ are

    ρ_ee(0) = 1,   ρ_gg(0) = ρ_ge(0) = ρ_eg(0) = 0.

As the master equations (<ref>) show, both quantities ρ_eg and ρ_ge remain zero before the first pulse (when Ω_x(t)≡ 0). The effect of the pulse is to interchange these two values, i.e. they both remain zero after the pulse. The same considerations can be applied for the second, third, etc. pulse, showing that ρ_eg(t)=ρ_ge(t)=0 for all t. 
Thus, to determine ρ”(t,t) we only need to find ρ_gg(t). We assume that the time instant t is between the M-th and the (M+1)-th pulse, i.e. t = M τ + (τ - τ_1) for some τ_1∈ [0,τ], as shown in Fig. <ref>. Immediately before the first pulse, at the time moment τ-0, we have:

    ρ_ee(τ)   =   e^-Γτ
    ρ_gg(τ)   =    1- e^-Γτ;

then at time 2τ -0 we have:

    ρ_ee(2τ)   =    (1 - e^-Γτ) e^-Γτ
    ρ_gg(2τ)   =    1- e^-Γτ +e^-2 Γτ,

at time 3τ -0: 

    ρ_ee(3τ)   =   e^-Γτ- e^-2Γτ +e^-3 Γτ
    ρ_gg(3τ)   =    1- e^-Γτ +e^-2 Γτ-e^-3 Γτ,



 ⋮



so that eventually, right before the M-th pulse, at time Mτ-0

    ρ_ee(Mτ-0) =  ∑_k=1^M (-1)^k-1e^- kΓτ = -∑_k=1^M (-1)^ke^- kΓτ ,

and right after the M-th pulse, which interchanges ρ_ee and ρ_gg, 

    ρ_ee(Mτ+0) =  1- ρ_ee(Mτ-0) = 1 + ∑_k=1^M (-1)^ke^- kΓτ.

Thus, at the time instant t=Mτ+(τ-τ_1), we have

    ρ_ee(t)   =   ( 1+ ∑_k=1^M (-1)^ke^- kΓτ) e^-Γ (τ -τ_1)
       =   ( ∑_k=0^M (-1)^ke^- kΓτ) e^-Γ (τ -τ_1)
       =   1-(-1)^M+1e^- (M+1) Γτ/1 + e^ - Γτe^-Γ (τ -τ_1),

and, since ρ_gg(t)=1-ρ_ee(t), we obtain

    ρ_gg(t) = 1-1-(-1)^M+1e^- (M+1) Γτ/1 + e^ - Γτe^-Γ (τ -τ_1).









Having established the explicit initial value of ρ”(t,t), now we can proceed evaluating the value of ρ”_ge(t,t+θ). Between the pulses both ρ^''_ge and ρ^''_eg evolve according to Eqs. <ref> and <ref>.

 

Thus, if t and t+θ belong to the same inter-pulse interval (i.e. when Mτ <t+θ < (M+1)τ), we have

    ρ^''_ge(t,t+θ) = e^(iΔ - Γ/2)θρ_gg(t)    and   ρ^''_eg = 0.

With increasing θ, at some point it will become equal to τ_1, and then the instant t+θ will coincide with the instant when a pulse is applied: t+θ=(M+1)τ+0. At this point the time instants t and t+θ will become separated by one pulse, and the value of ρ”_eg will be interchanged with ρ”_ge, i.e. when θ=τ_1+0 we will have already 

    ρ^''_ge(t,t+θ)    =    0 
    ρ^''_eg(t,t+θ)    =   e^(iΔ - Γ/2)τ_1ρ_gg(t)  .

At this point, the accumulation rate of the phase in ρ”_ge and ρ”_eg changes sign: note that the factors on the right-hand sides of Eqs. <ref> and <ref> have opposite imaginary parts, iΔ and -iΔ, respectively. Thus, right before the next pulse, when t+θ=(M+2)τ-0 (i.e. when θ=τ+τ_1-0), we will have 

    ρ^''_ge(t,t+θ)    =    0 
    ρ^''_eg(t,t+θ)    =   e^(iΔ - Γ/2)τ_1e^(-iΔ - Γ/2)τρ_gg(t)

and right after the pulse, when t+θ=(M+2)τ + 0 (i.e. when θ=τ+τ_1 + 0), the values will be interchanged again:

    ρ^''_ge(t,t+θ)    =   e^(iΔ - Γ/2)τ_1e^(-iΔ - Γ/2)τρ_gg(t) 
    ρ^''_eg(t,t+θ)    =    0.







Proceeding further in this way, right before the next pulse, at θ=τ_1+2τ - 0, we get

    ρ^''_ge(t,t+θ)    =   e^(iΔ - Γ/2)τ_1e^- Γτρ_gg(t) 
    ρ^''_eg(t,t+θ)    =    0.

Note that the phase of ρ”_ge still equals to iΔτ_1, because after each pulse the phase accumulation rate changes sign. Further, at θ=τ_1+3τ - 0, 

    ρ^''_ge(t,t+θ)    =    0 
    ρ^''_eg(t,t+θ)    =   e^(iΔ - Γ/2)τ_1 -iΔτ - 3Γ/2τρ_gg(t)




Thus we obtain that for θ=τ_1+(m-1)τ - 0 with m even, as shown in Fig. <ref>,

    ρ^''_ge(t,t+θ)    =    0 
    ρ^''_eg(t,t+θ)    =   e^(iΔ - Γ/2)τ_1 -iΔτ - (m-1)Γ/2τρ_gg(t),

and for θ=τ_1+(m-1)τ + τ_2 with m even and τ_2 < τ,

    ρ^''_ge(t,t+θ)    =   e^(iΔ - Γ/2)τ_1 -iΔτ - (m-1)Γ/2τe^iΔτ_2 -Γ/2τ_2 ρ_gg(t) 
    
    		        =   e^-Γθ /2e^iΔ ( τ_1 + τ_2 - τ)ρ_gg(t) 
    ρ^''_eg(t,t+θ)    =    0.

Altogether we can write

    ρ^''_ge(t,t+θ) = f(t,θ) ρ_gg(t)

with ρ_gg(t) given above by Eq. <ref>, and





 
  * for t and t+θ in the same pulse interval,
 
    f(t, θ) = e^(i Δ - Γ/2)θ


  * for t and t+θ separated by an odd number of pulses,

    f(t, θ) = 0


  * for t and t+θ separated by an even number m of pulses, i.e. when θ=τ_1+(m-1)τ+τ_2 with even m and τ_2<τ (Fig. <ref>), 

    f(t, θ) = e^-Γθ /2e^iΔ ( τ_1 + τ_2 - τ)=e^-Γθ /2e^iΔ (θ - mτ).


Note that, due to the pulses, the phase of the function f(t,θ) does not grow linearly with θ, being confined to the interval [-τΔ,τΔ] at all values of t and θ. This is the reason why, for small inter-pulse delay τ≪Δ^-1, both emission and absorption are concentrated in the vicinity of ω=0 instead of ω=Δ.

With this result, we can now rewrite the direct absorption integral 𝒫_2(ω)



in the form

    𝒫_2(ω) = ∫_0^T   dt  ρ_gg(t) ∫_0^T-t dθ f(t, θ) e^-i ωθ


First let us evaluate the inner integral, that we will denote as I_θ, using the explicit form of f(t,θ) above:


    I_θ   =   ∫_0^τ_1  dθ e^-i ωθe^(iΔ-Γ/2)θ + ∫_τ_1+τ^τ_1 + 2τ  dθ e^-i ωθe^-i 2 Δτ + ( i Δ-Γ/2) θ
       +   ∫_τ_1+3τ^τ_1 + 4τ  dθ e^-i ωθe^-i 4 Δτ + ( i Δ-Γ/2) θ + ⋯ +
    ∫_τ_1+(m-1)τ even m^τ_1 + mτ dθ e^-i ωθe^-imΔτ + (i Δ-Γ/2) θ + ⋯
       =   e^[i(Δ -ω) - Γ/2]τ_1  -1/ i(Δ -ω) - Γ/2 + ∑_m=2 even m^m_max∫_τ_1 + (m-1)τ^τ_1 + mτ dθ e^-i ωθe^[ -i mΔτ + (i Δ-Γ/2) θ]

where the summation is over even values of m, and m_max is the maximum value of m; since it has to be even, its specific value depends on whether M is odd or even, see below for details. 
Defining γ_0 = i(ω -Δ ) + Γ/2, we can write

    I_θ   =   1 - e^-γ_0 τ_1/γ_0 + ∑_m=2 even m^m_maxe^ -imΔτe^-γ_0(τ_1+(m-1)τ)-e^-γ_0(τ_1+mτ)/γ_0
    
    
                  =   1-e^-γ_0τ_1/γ_0 + e^-γ_0τ_1/γ_0(e^γ_0 τ-1) 1-e^-m_maxγ_1τ/e^2γ_1 τ-1 ,


where we have introduced γ_1 = Γ/2 + iω. Note that this result correctly reproduces the situation of m_max<2, i.e. when m_max=0; this happens when t belongs to the last inter-pulse interval of the sequence, and θ varies only from zero to τ_1. Then the value of I_θ is given by the first integral in Eq. <ref>, while the remaining sum over m is zero. Thus, we do not need to worry about this special case in the calculations below. 

Now we need to evaluate the outer integral:

    𝒫_2(ω) =  ∫_0^T   dt  ρ_gg(t) I_θ


with the quantity ρ_gg calculated earlier,

    ρ_gg(t)    =    1 - ρ_0(M) e^-Γ (τ -τ_1)
       =    1 -  1- (-1)^M+1e^-(M+1)Γτ/ 1 + e^-Γτe^-Γ (τ -τ_1),

where we introduced the shorthand notation ρ_0(M) for the awkward fraction appearing on the second line. In this way, we represent 𝒫_2 as

    𝒫_2    =   ∫_0^T dt   I_θ - ∫_0^τ dt  ρ_0(0) e^-Γ(τ-τ_1) I_θ
       -   ∫_τ^2τ dt  ρ_0(1) e^-Γ(τ-τ_1) I_θ
       -   ⋯ - ∫_(N_p-1)τ^N_pτ dt  ρ_0(N_p-1) e^-Γ(τ-τ_1) I_θ
       =   ∫_0^T dt   I_θ
       -   ∑_M=0^N_p-1 ρ_0(M) ∫_0^τ dt_1  e^-Γ t_1(1/γ_0 + e^-γ_0 τ_1 I_θ^(1))

where we have defined t_1 =t-Mτ= τ - τ_1 and

    I_θ   =   1/γ_0 + e^-γ_0 τ_1 I_θ^(1)
    
     I_θ^(1)   =    -1/γ_0 + e^γ_0 τ -1/γ_0 1-e^-m_maxγ_1τ/e^2γ_1τ-1
    
        =    -1/γ_0 + I_θ^(2)

where we introduced the shorthand notation I_θ^(2) for another awkward fraction, the second summand on the second line above.

Now we have

    𝒫_2 =    -   ∑_M=0^N_p-1ρ_0(M) ∫_0^τ dt_1  e^-Γ t_1(1/γ_0 + e^-γ_0 (τ-t_1) I_θ^(1)) 
       +   ∫_0^T dt   I_θ
    
    =    -   ∑_M=0^N_p-1ρ_0(M) [1 - e^-Γτ/γ_0 Γ  + I_θ^(1)e^-γ_0 τe^γ_2 τ -1/γ_2] 
       +   𝒫_3

with 𝒫_3 = ∫_0^T dt   I_θ and γ_2 = γ_0 - Γ = i(ω - Δ) - Γ/2. 

Below we will show that the first sum gives exactly the contribution from the stimulated emission 𝒫_1(ω). It is convenient to calculate the simpler term 𝒫_3 first. We can rewrite 𝒫_3 as:

    𝒫_3    =   ∫_0^T dt/γ_0 +  ∫_0^τ dt  e^-γ_0 τ_1 I_θ^(1)
       +   ⋯ + ∫_(N_p-1)τ^N_pτ dt  e^-γ_0 τ_1 I_θ^(1)
       =   N_pτ/γ_0 + e^-γ_0 τ∑_M=0^N_p-1 I_θ^(1)∫_0^τ dt_1  e^γ_0 t_1
    
    
    
       =   N_pτ/γ_0 + ∑_M=0^N_p-1 I_θ^(1)1- e^-γ_0 τ/γ_0



With the explicit form of I_θ^(1) given above, we have 

    𝒫_3 
       =   N_pτ/γ_0 - ∑_M=0^N_p-11 - e^-γ_0 τ/γ_0^2 + ∑_M=0^N_p-11 - e^-γ_0 τ/γ_0^2 I_θ^(2)
       =   N_pτ/γ_0 - N_p/γ_0^2(1 - e^-γ_0 τ) + 𝒫_4

where 

    𝒫_4    =   1-e^-γ_0τ/γ_0 e^γ_0τ-1/γ_0∑_M=0^N_p-11-e^-γ_1τ m_max/e^2γ_1 τ-1
    
    
    
       =   e^γ_0 τ+e^-γ_0 τ-2/γ_0^2 (e^2γ_1τ-1)[ N_p -∑_M=0^N_p-1e^-γ_1τ m_max]


In order to calculate the last sum in the equation above, we need to determine m_max. To do this let us consider the case of N_p=2K, i.e. when the number K of the full cycles of the sequence has been applied to the TLS. Let us recall that we represent t=Mτ+(τ-τ_1), i.e. M is the number of pulses between zero and t. The number of pulses separating t and t+θ is m, and the maximum value of θ is θ_max=T-t, which limits the maximum value of m; however, f(t,θ) is zero if m is odd, so that m_max should be even.
Therefore, starting from larger values of t, we obtain:

 
  * if t∈[T-τ,T] then θ_max=(T-t) ∈[0,τ], so that if M=2K-1 then m_max=0,
 
  * if t∈[T-2τ,T-τ] then θ_max=(T-t) ∈[τ,2τ], so that if M=2K-2 then m_max=0 because m_max should be even,
 
  * if t∈[T-3τ,T-2τ] then θ_max=(T-t) ∈[2τ,3τ], so that if M=2K-3 then m_max=2,
 
  * if t∈[T-4τ,T-3τ] then θ_max=(T-t) ∈[3τ,4τ], so that if M=2K-4 then m_max=2 (should be even),
 
  * if t∈[T-5τ,T-4τ] then θ_max=(T-t) ∈[4τ,5τ], so that if M=2K-5 then m_max=4, 
 
  * if t∈[T-6τ,T-5τ] then θ_max=(T-t) ∈[5τ,6τ], so that if M=2K-6 then m_max=4 (should be even),

  * ⋯
 
  * if t∈[τ,2τ] then θ_max=(T-t) ∈[T-2τ,T-τ], so that if M=1 then m_max=2K-2,
 
  * if t∈[0,τ] then θ_max=(T-t) ∈[T-τ,T], so that if M=0 then m_max=2K-2 (should be even).

To summarize, if we parametrize M=2n for even M and M=2n+1 for odd M, where n varies from 0 to K-1, then m_max=2(K-n-1) for both M=2n and M=2n+1.

Thus, the last sum in Eq. <ref> is calculated as

    ∑_M=0^N_p-1e^-γ_1τ m_max   =    2∑_n=0^K-1e^-2γ_1τ (K-n-1)
       =    21-e^-γ_1τ N_p/1-e^-2γ_1τ,

where the factor 2 appears because m_max is the same for both M=2n and M=2n+1, so that the sums over odd M and even M are combined. Putting all terms together, we obtain


    𝒫_3 =  N_pτ/γ_0 -  N_p/γ_0^2(1 - e^-γ_0 τ)
    + e^γ_0 τ+e^-γ_0 τ-2/γ_0^2(e^2γ_1τ-1) [N_p - 21-e^-γ_1τ N_p/1-e^-2γ_1τ]




Now, the calculation of the emission term 𝒫_1 can be simplified if we notice that ρ'(t,s) obeys the same equations of motion as ρ”(t,s), and is transformed by the pulses in exactly the same way. Therefore, the quantity ρ'_ge(t,s) (that determines 𝒫_1) evolves in exactly the same way as ρ”_ge(t,s), and the difference between them is only in the initial condition: at s=t we have ρ'_ge(t,t)=ρ_ee(t), while ρ”_ge(t,t)=ρ_gg(t)=1-ρ_ee(t). Thus, the reasoning that was used in deriving 
Eqs. <ref>–<ref> can be directly applied to ρ'_ge(t,s) if ρ_gg(t) is substituted by ρ_ee(t), due to the linearity of the master equations. As a result, we immediately see that ρ'_ge(t,t+θ) has the form 

    ρ'_ge(t,t+θ) = f(t,θ) ρ_ee(t)

with the same function f(t,θ). Thus, the integral I_θ can be used without modifications in the calculation of 𝒫_1, and, since ρ_gg(t)=1-ρ_ee(t), we immediately obtain

    𝒫_1(ω) = ∫_0^T  ρ_ee(t) I_θ  dt  = ∫_0^T   [1-ρ_gg(t)] I_θ  dt = 𝒫_3 - 𝒫_2.

Comparing this expression with Eq. <ref> above, we obtain an explicit expression

    𝒫_1    =   ∑_M=0^N_p-1ρ_0(M) ∫_0^τ dt_1  e^-Γ t_1(1/γ_0 + e^-γ_0 (τ-t_1) I_θ^(1)) 
       =   ∑_M=0^N_p-1ρ_0(M) [1 - e^-Γτ/γ_0 Γ  + I_θ^(1)e^-γ_0 τe^γ_2 τ -1/γ_2]











Now let us evaluate the sums appearing in this expression. First, we need the sum

    ∑_M=0^N_p-1   ρ_0(M)    = 1/1 + e^-Γτ∑_M=0^N_p-1[ 1+e^-Γτ (-e^-Γτ)^M   ] 
       =   N_p/1+e^-Γτ + e^-Γτ/1+e^-Γτ1 -(-e^-Γτ)^N_p/1 + e^-Γτ

and for sufficiently large N_p, when the exponentially small terms can be omitted, this yields


    ∑_M=0^N_p-1ρ_0(M) ≈N_p/1+ e^-Γτ + e^-Γτ/(1 + e^-Γτ)^2,


The second required sum is

    ∑_M=0^N_p-1ρ_0(M) e^-γ_1 τ m_max,

and in order to evaluate it we use the same parametrization as above, M=2n for even M and M=2n+1 for odd M, with n=0,…,K-1. We pair the neighboring terms, i.e.
    ∑_M=0^N_p-1   ρ_0   (M) e^-γ_1 τ m_max
       =    [∑_ even M +∑_ odd M] ρ_0(M) e^-γ_1 τ m_max
       =    ∑_n=0^K-1[ ρ_0(2n) + ρ_0(2n+1)] e^-2γ_1 τ (K-n-1).

Since 

    ρ_0(2n) + ρ_0(2n+1) = 2 + e^-Γτ (2n+1) - e^-Γτ (2n+2) / 1 + e^-Γτ,

we obtain

    ∑_M=0^N_p-1   ρ_0   (M) e^-γ_1 τ m_max
       =   ∑_n=0^K-12 - e^-Γτ (2n+1) + e^-Γτ (2n+2) / 1 + e^-Γτe^-2γ_1 (K-n-1) τ
    
       =   e^-2γ_1 (K-1)τ/ 1 + e^-Γτ[2e^2Kγ_1τ-1/e^2γ_1τ-1. 
       +   . e^-Γτ(1- e^-Γτ) e^2K(γ_1-Γ)τ-1/e^2(γ_1-Γ)τ-1]


















Substituting these results into Eq. <ref> for 𝒫_1, we obtain in the limit of large N_p




    𝒫_1(ω)    =   1/(1 + e^-Γτ)γ_0[ ( 1-e^-Γτ/Γ - e^-γ_0 τe^γ_2 τ-1/γ_2
     + e^γ_2 τ-1/γ_2 1-e^-γ_0 τ/e^2 γ_1 τ-1) 
    
    (N_p + e^-Γτ/1 + e^-Γτ) 
       -   e^γ_2 τ-1/γ_2 1-e^-γ_0 τ/e^2 γ_1 τ-1 (  
    2 e^-N_p γ_1 τ - 1 /e^-2 γ_1 τ - 1 + (e^-Γτ - e^ - 2 Γτ)e^-N_p γ_1 τ/e^-2 γ_1 τ - e^-2 Γτ) 
    
    
    
    ]


and the net absorption spectrum is obtained as

    Q(ω)    =    2 A^2 Re{𝒫_2(ω) - 𝒫_1(ω) }
       =    2 A^2 Re{𝒫_3(ω) -2𝒫_1(ω) }

with the explicit analytical expressions for 𝒫_1 and 𝒫_3 given above.

In Fig.<ref> we show a comparison of the numerical result and the analytical result described above for the absorption spectrum of a two-level system with detuning Δ = 3.0 driven by a periodic pulse sequence of period τ = 0.2 after 20 pulses. The comparison reveals a very good agreement between the solutions, and the agreement indeed improves as the number of pulses increases.

