Resonance fluorescence of noisy systems

Light scattering from resonantly or nearly resonantly excited systems, known as resonance fluorescence, has been gaining importance as a versatile tool for investigating quantum states of matter and readout of quantum information, recently including also the inherently noisy solid state systems. In this work we develop a general theory of resonance fluorescence in the low excitation limit on systems in which the transition energy is subject to noise for two important classes of noise processes: white noise fluctuations that lead to phase diffusion and an arbitrary stationary Markovian noise process on a finite set of states. We apply the latter to the case of random telegraph noise and a sum of an arbitrary number of identical random telegraph noise contributions. We show that different classes of noise influence the RF spectrum in a characteristic way. Hence, the RF spectrum carries information on the characteristics of noise present in the physical system.


Introduction
Light-matter interaction is one of the main tools for studying various properties of physical systems. In particular, resonant or nearly resonant light scattering, known as resonance fluorescence (RF) [1,2] has been used for a long time to characterize systems of various kinds [3,4,5,6,7]. More recently, the RF technique has found a variety of applications in condensed matter systems, both in the physical investigation of quantum emitters, as well as in manipulating and reading out the quantum information encoded in solid-state qubits [8,9]. It has been used to observe spin dynamics in semiconductor quantum dots (QDs) [10,11], to interface [12] and entangle [13] QD spins with single photons, to generate indistinguishable photons [14], to read out spin states in QDs [15,16,12,11,10] and silicon defects [17] as well as to demonstrate quantum-optical effects in macroscopic superconducting qubits [18,19]. The recent observation of RF from a waveguide-coupled solid-state emitter [20] opens a perspective of on-chip device integration.
The resonance fluorescence from a single unperturbed quantum emitter shows different properties under weak and strong excitation. In the former case, energy conservation for each single-photon scattering event leads to a single line with the broadening limited by the laser line width [1]. Under strong excitation, the modulation of the system state due to Rabi rotations gives rise to a triplet of broadened lines separated by the Rabi frequency, referred to as the Mollow triplet [21,1].
No physical system is completely isolated from its environment. The environmental impact is of particular importance in solid-state systems, where the optical properties are to a large extent influenced by the coupling to charge [22,23] or spin fluctuations (both nuclear [24] and electronic [25]), as well as to lattice vibrations. The latter can be induced in a coherent way, leading to controllable modification of the scattering spectrum [26,27], but in most cases is a source of noise [28,29]. When treated as a classical background noise, the environmental fluctuations are often modeled using Gaussian distributions [22,23]. However, it is generally believed that the underlying physics involves discrete state changes of nearby physical objects, like nuclear or dopant spin flips, or charging and discharging of defects, and such dynamics is indeed observed in certain experiments [25]. Therefore, a more fundamental model for the description of noise needs to be based on telegraph-noise dynamics.
While numerous studies considered fluctuations in the phase [30,31,32,33,34], amplitude [30,35,36], and frequency [30,32,37] of the laser beam, much less attention has been devoted to the effect of the environmental noise on RF. The existing studies include the particular case of interaction with phonons [38,39] and the dynamics of inversion [40], coherence [41], entanglement [42], as well as non-linear wave-mixing response [43] in two-level systems subject to environmental fluctuations. Apart from studying the detrimental effects of noise on the dynamics of quantum systems, the latter can also be used as noise sensors in order to determine the properties of the noise itself. Such noise characterization is crucial for the robustness of quantum systems, hence considerable effort has recently been invested in the development of noise spectroscopy techniques [44,45,46,47].
In this work, we generalize the recently proposed description of RF from a deterministically modulated two-level system [26,27] to systems subject to random fluctuations. We formulate a theory of low-excitation RF (coherent Rayleigh scattering) spectra for a system with the transition energy subject to white noise or an arbitrary Markovian random process that shifts the transition energy between a number of discrete spectral positions. In particular, we consider single-source symmetric and asymmetric telegraph noise and multi-source telegraph noise.
By relating the scattering spectrum to the formal characteristics of the underlying noise process, we are able to show that the spectrum bears clear fingerprints of the properties of noise: In the white noise (phase diffusion) case, the noise leads to Lorentzian line broadening. In contrast, a slow discrete process on a small set of states results in a multiple of Lorentzian lines that merge into a broad Gaussian feature when the number of process states increases so that they become dense on the energy axis. This picture changes when the process is fast. In this case a motional narrowing effect leads to the appearance of a single Lorentzian line.
The paper is organized as follows. In Sec. 2 we describe the system and define its general model. Sec. 3 contains the essential definitions and presents the general framework of our theoretical description; here we also present the theory for the simplest case of white noise. Sec. 4 contains the central formal result of the paper: the theory of the resonance fluorescence spectrum for a Markovian noise process. Sec. 5 presents the results obtained by applying the theory to selected noise processes. The paper is concluded by Sec. 6.

System and Model
We consider light scattering on a two-level system which is subject to environmental fluctuations that randomly shift the energy of the excited state. As in the standard model of resonance fluorescence [1], the system is driven by a resonant or nearly resonant monochromatic laser light and undergoes spontaneous emission.
We denote the laser frequency by ω L and the system states by |0⟩ and |1⟩. Let ℏω 0 (t) be the time-dependent (fluctuating) energy difference between these states. The system is then described by the Hamiltonian where is the laser field (treated classically) with the amplitude E 0 and d is the dipole moment operator. We assume ⟨0|d|0⟩ = ⟨1|d|1⟩ = 0. The system relaxation due to spontaneous emission is accounted for by the Lindblad dissipator where ρ is the density matrix of the system, {A, B} = AB + BA, and σ + = σ † − = |1⟩⟨0|. The system state in the rotating frame is defined bỹ The Hamiltonian in this rotating frame and in the rotating wave approximation is where we define the detuning ∆(t) = ω L − ω 0 (t) and the Rabi frequency Ω = (E * 0 /ℏ) · ⟨0|d|1⟩ that we assume real. The Master Equation describing the evolution of the system state has the form The randomly changing detuning ∆(t), which reflects the environmental fluctuations, is the central feature of our model. Formally, it is described by a stochastic process, the properties of which may depend on the physical situation. Here we assume that ∆(t) is a stationary Markov process with a unique stationary probability distribution p (st) reached asymptotically in the long time limit. Two particular examples will be discussed in the following.

The resonance fluorescence spectrum of a noisy system
In this section we define the formal quantities relevant to the RF spectrum, present the general framework of the theoretical model and discuss the simplest case of phase diffusion due to white noise.
In the Markov approximation, the detected spectrum of scattered light can be related to the system autocorrelation function G(t, t + τ ) = ⟨σ + (t)σ − (t + τ )⟩ by [1] Here σ ± (t) are operators in the Heisenberg picture relative to the rotating-frame evolution as defined by Eq. (2) and Γ is the (Lorentzian) instrumental broadening accounting for the finite resolution of the detection. Let us formally denote the solution of Eq. (2) byρ(t) = L t 0 ,t [ρ(t 0 )], where L t 1 ,t 2 is the evolution superoperator. The Lax quantum regression theorem then yields the autocorrelation function in the form [2,27] Here t 0 is the initial moment of the evolution, while t − t 0 is a sufficiently long time for the system to reach its steady state. The total scattering intensity is I tot = ∞ −∞ F (ω)dω. In the noise-free limit, the RF spectrum consists of a Dirac delta (broadened by the instrumental resolution) corresponding to elastic light scattering, which survives to various extent in the noisy case. Its intensity will be denoted by I el . The remaining part of the spectrum is due to inelastic scattering induced by the random fluctuations. Its intensity is I inel .
The equations of motion for the elements ρ 01 , ρ 10 , ρ 11 of the density matrix, following from Eq. (2), have the formρ jl = a jl ρ jl + iΩ mn b jl,mn ρ mn , where a 11 = −γ, a 01 = a * 10 = i∆ − γ/2, and b 11,10 = −b 11,01 = b 10,11 = −b 01,11 = 1/2. The same holds for an arbitrary matrix, not necessarily a density matrix. Since Eq. (2) is trace-preserving, one has ρ 00 = c 0 − ρ 11 , where c 0 is a constant determined by the initial values (c 0 = 1 for a density matrix). In the absence of the laser field (Ω = 0) the equation of motion can be solved trivially to yield the zeroth-order propagation In the weak excitation regime, one can then solve Eq. (2) iteratively in the subsequent orders r > 0 in Ω, These equations fully define the perturbative expansion of the evolution superoperator L t 1 ,t 2 in powers of Ω. Substituting this evolution into Eq. (4) one finds, in the leading order of Ω 2 , the autocorrelation function for an arbitrary time-dependent energy shift ∆(t) in the form Here we changed the variables according to s = t + u, set t − t 0 → ∞ (steady-state regime), averaged over the realizations of the noise (denoted by a line above the averaged quantities) and used the fact that the noise is stationary, hence e iΦ(t d +s,tc+s)±iΦ(t b +s,ta+s) = e iΦ(t d ,tc)±iΦ(t b ,ta) .
The detailed derivation along with a graphical interpretation of the evolution paths contributing to the RF signal can be found in the Supplementary Material to Ref. [27]. Before developing the theory for an arbitrary Markovian noise process on a discrete state space, we find the explicit form of the correlation function in the case of simple phase diffusion. We assume that where ∆ W N (t) is a stationary white noise process with ⟨∆ W N (t)⟩ = 0 and ⟨∆ W N (t)∆ W N (t + τ )⟩ = Dδ(τ ). Then, the phase shift given in Eq. (6) has a normal distribution with the mean value ∆(t b −t a ) and variance D(t b −t a ), hence D is the phase diffusion coefficient. With this Gaussian distribution, Eq. (5) can easily be evaluated by using the statistical independence of phase shifts over non-overlapping periods of time. This yields The RF spectrum obtained by substituting Eq. (8) to Eq. (3) is a sum of elastic and inelastic peaks that are, respectively, given by and Thus, the elastic scattering line is located at the laser frequency, while the inelastic line appears at the average system transition frequency. By integrating we find the total intensity as well as the intensities of the elastic and inelastic components Here and in the following we relate the scattering intensities to the standard resonance fluorescence intensity in weak excitation limit [1], I 0 = πΩ 2 /γ 2 .

Theory of the RF spectrum for arbitrary Markovian noise
In this section we develop the theory of the resonance fluorescence spectrum for a system subject to noise that can be described by a continuous-time stationary Markov process on a finite set of states {∆ i } k i=1 . This model can also be used as an approximation to more general Markov processes, based on a physically motivated truncation and discretization of the state space of the noise process.
The process is characterized by the transition probabilities P m,n (τ ) = P[∆(t + τ ) = ∆ m |∆(t) = ∆ n ] (with P denoting the conditional probability) forming the transition matrix P (τ ) = exp(Cτ ), with the generator C = dP (τ )/dτ | τ =0 . As shown in Appendix A, for such a process where B = C + i diag(∆ 1 , ∆ 2 , . . . , ∆ k ), P ∞ = lim τ →∞ P τ = p (st) (1, . . . , 1), with p (st) representing the column vector of stationary probabilities, and the conjugation in the last term refers to the '−' sign on the left-hand side. Substituting Eq. (13) to Eq. (5) one gets Upon substituting to Eq. (3), the first term in Eq. (14) can be factorized by reordering the integrals with respect to τ and u ′ and then changing the variable according to τ = u ′ + s. The second term can be evaluated directly. As a result one gets While this closed analytical form may be convenient for evaluating the spectrum in the case of a small state space of the process, much more insight is gained by relating the RF spectrum to the spectral properties of the generator C. To this end, we transform C to the Jordan form (over the field of complex numbers) by the similarity transformation where C j are Jordan blocks belonging to the respective eigenvalues λ j of algebraic multiplicity d j . We then apply the Jordan- Let Π j be the projector on the subspace supporting C j . Since In the case of a diagonalizable matrix C, d j = 1 for all j, the above procedure reduces to simple diagonalization, and λ j become eigenvalues of C in the most common sense. Using this result in the first term of Eq. (14), substituting to Eq. (3), and performing the integrations one finds The essential factors that describe the form of the spectrum are those depending on the frequency ω.
The first term defines Lorentzian features and more general line shapes in the case of degenerate eigenvalues (d j > 1), as well as the corresponding dispersive contributions, with central frequencies at ω L + Im λ j and broadened by − Re λ j (here additionally broadened by the instrumental resolution Γ). The existence of a stationary state implies that one of the eigenvalues (say, λ 0 ) is zero and P 00 = P ∞ . This contribution leads to an unbroadened (apart from the instrumental broadening) elastic scattering peak centered at ω = ω L . For λ 0 there is no dispersive function. Each real non-zero eigenvalue leads to a broadened spectral feature at ω = ω L .
The second term is a sum of simple Lorentzians and the corresponding dispersive contributions with positions and widths determined by the spectrum of B and further broadened due to spontaneous emission.
The total and elastic scattering intensities, found by integrating Eq. (16), are respectively and As an application of this formalism, we study in detail the special case of N identical noise sources, each generating telegraph noise (TN).
A single noise source has two states that contribute ∆ = ± ∆ 0 2 √ N to the system energy shift. The switching rates between the two states of the noise are β ↑ and β ↓ , leading to stationary probabilities β ↓ / (β ↑ + β ↓ ) and β ↑ / (β ↑ + β ↓ ) for the two noise states. For N identical and additively contributing noise sources the space of possible noise states is composed of N + 1 values of the total system detuning, corresponding to j sources in the "upper state". Here ∆ is the mean detuning (the laser detuning from the noise-free transition energy) and the renormalization by a factor 1/ √ N assures convergence in the limit of N → ∞. For this case the generator C is a N + 1-dimensional tridiagonal matrix, where the non-zero elements are The stationary probability follows the binomial distribution

Results
In this section we present the results for the RF spectrum based on the theory developed in Sec. 3 and Sec. 4 for three noise models. For the white noise calculations we set Γ/γ = 0.25. For the discrete process we choose γ = 4Γ = 0.02∆ 0 . As a natural unit for presenting and comparing the RF spectra we will use the maximum value of the spectrum for an unperturbed system under weak resonant excitation [1], F 0 = Ω 2 /(γ 2 Γ), which can be obtained from Eq. (9) with D = ∆ = 0 and ω = ω L . All the spectra presented in the following will be related to this quantity.

White noise
We start the discussion with the case of a system affected by white noise. Fig. 1 shows the results for this case, calculated from Eq. (9) and Eq. (10), under resonant and detuned excitation in the left and right columns, respectively. In both cases, the RF spectrum is composed of two Lorentzians. In the resonant case they overlap, while for a detuned excitation they are split. As expected, the overall intensity is also lower in the latter case.
The total spectrum, shown in Fig. 1(a,b), is decomposed into the elastic and inelastic contributions in Fig. 1(c,d) and Fig. 1(e,f), respectively. While the positions of the spectral features do not change, the evolution of their intensities and of the width of the non-elastic line are clearly visible. The width of the inelastic contribution grows with D. The intensity of the elastic contribution starts to decrease when D ∼ γ, while the intensity of the inelastic one changes non-monotonically with D, reaching a maximum for D ∼ γ. The decrease of the inelastic contribution for weak noise is an obvious consequence of restoring the noise-free limit of purely elastic scattering. For strong noise both components decrease because the increasing spread of the transition energy reduces the effective overlap with the excitation frequency, which affects the excitation intensity. The appearance of an additional broadened spectral line at a spectral position bound to the transition energy in addition to the elastic line at the laser position can be understood by invoking the model of a classical charged harmonic oscillator driven by a periodic force. This analogy is formally validated by the fact that in the leading order in Ω, the whole emitted light is coherent, that is, originating from the transition dipole induced coherently by the laser field. In its steady state, the classical system oscillates periodically with the laser frequency, which gives rise to a sharp line at this spectral position. However, any perturbation of the steady-state evolution leads to the appearance of a damped transient at the eigenfrequency of the system renormalized by damping. Here the noise serves as a perturbation that permanently excites the transient oscillations and simultaneously damps the coherence due to phase diffusion, leading to a broadened line.
A systematic study of the intensities as a function of the noise strength D is presented in Fig. 1(g,h). As can also be seen directly from Eq. (12), the dependence has asymptotically a power-law character. In the case of resonant excitation, ∆ = 0, the absence of noise (D → 0) leads to permanent resonance condition, maximizing the total scattering intensity. For a sufficiently large detuning, ∆ > γ/2, the interplay of the decreasing elastic scattering and non-monotonic inelastic one leads to the appearance of a maximum of the total scattering intensity at D = 2∆ − γ. In this case one observes a noise-induced enhancement of scattering: At the maximum, the total intensity is larger than in the noise-free case by a factor of (γ 2 + 4∆ 2 )/(4γ∆).

Single-source telegraph noise (TN)
In this section we discuss the results for the scattering spectra and line intensities for a single source of random telegraph noise, depending on the characteristics of the noise dynamics (the switching rates β ↑,↓ ) and laser detuning ∆ from the average transition energy. In this case, the detuning takes randomly two values, hence the set of states of the stochastic process is reduced to ∆ = ∆ ± ∆ 0 /2. The generator in Eq. (20) reduces to Its eigenvalues are λ 0 = 0 and λ 1 = − (β ↑ + β ↓ ). Both eigenvalues are nondegenerate, hence the spectrum is composed of four simple Lorentzians and, possibly, the corresponding dispersive functions. For the sake of presentation we set β ↑ = β(1−x) and We begin with the case of symmetric noise (x = 0, that is, β ↑ = β ↓ = β) and then describe the effects of noise asymmetry.

Symmetric switching
We start our analysis by discussing the total scattering intensity I tot , depending on the spectral position of the laser. The full RF intensity as a function of the laser detuning and noise switching rate in the case of symmetric noise is shown in Fig. 2, where we plot the total scattering intensity obtained by numerical evaluation of Eq. (17). At low switching rates (slow noise), the scattering intensity is the largest when the laser is tuned to one of the two randomly alternating spectral positions of the optical transition, while at high switching rates (fast noise) the strongest scattering occurs at the average value of the transition energy.
The slow-noise case is easily understood as the quasi-static limit of the random dynamics: over times much longer than the characteristic time scale of the system evolution, 1/∆ 0 , the system transition energy remains constant at one of the two spectral Thus, in the limit of slow noise the scattering intensity reaches half of the noise-free resonant scattering intensity at each of the two resonant spectral positions.
In the opposite limit of fast noise the random switching takes place many times during the characteristic time 1/∆ 0 , hence the accumulated dynamical phase slip with respect to the laser light depends only on the averaged transition energy, leading to the resonance condition at ∆ = 0. By directly evaluating Eq. (17) in the limit of β → ∞ one finds in this case hence the full standard intensity is recovered at the average spectral position. For moderate values of β one can neglect γ in Eq. (17), hence the noise speed at which the transition between the slow and fast regimes takes place can only depend on ∆ 0 . Indeed, the scattering intensities at ∆ = ∆ 0 /2 and ∆ = 0 become equal for β/∆ 0 = 1/4 + O(γ/∆ 0 ). The particular form of the RF spectrum at a given spectral position of the laser is determined by the poles of Eq. (16), which we denote by ω j , j = 0, . . . , 3, with Correspondingly, the Lorentzian features located at the spectral positions Re ω j will be labeled by L 0 , . . ., L 3 , respectively. Both ω 0 − ω L and ω 1 − ω L are purely imaginary, corresponding to resonant peaks. The positions of the spectral features (Re ω j ) as a function of the switching rate β, calculated from fEqs. (25a)-(25c), are shown in Fig. 3(a) for ∆ = 0. Currently, we focus on symmetric switching, x = 0, represented by solid lines. Apart from the resonant (central) features L 0 and L 1 the system in the slow switching limit shows two side peaks L 2 and L 3 . Their positions evolve from ω L ±∆ 0 /2 (the locations of the transition energy) in the quasi-static limit towards ω L , undergoing a qualitative transition at β = ∆ 0 /2, where the characteristic frequencies collapse to a single value of ω = ω L . Hence, at β = ∆ 0 /2 the RF spectrum changes its form from three lines to a single line. This resembles the properties of a damped harmonic system. Here, however, the transition is driven by the switching rate of the noise instead of the damping magnitude, with the cases of slow and fast noise corresponding to the underdamped and overdamped regimes, respectively. According to Eqs. (25a)-(25c), changing the laser detuning does not affect the spectral positions of the lines L 2 and L 3 (with respect to the fixed transition energy), while the lines L 0 and L 1 follow the frequency of the laser. The broadening of the spectral features, | Im ω j |, is presented in in Fig. 3(b), where we set the instrumental broadening Γ = 0, keep the color coding from Fig. 3(a), and omit the L 0 peak that has zero width. It follows from Eqs. (25a)-(25c) that the widths are independent of the mean detuning ∆. As follows directly from Eq. (25b), the resonant peak is broadened by 2β and becomes unbroadened (corresponding to purely elastic scattering) in the quasi-static limit. The side peaks are symmetric in the slow switching regime, with the broadening decreasing to γ/2 in the quasi-static limit. In the fast regime, when they take the same spectral location, one of them is narrowing asymptotically as γ/2 + O[(∆ 0 /β) 2 ], while the other one is broadening asymptotically as 2β + γ/2.
i Obviously, the peak positions and widths do not provide the complete information about the spectrum, as long as the intensities are not known. In order to fully analyze the spectra we evaluate F (ω) from Eq. (15) and plot the result in Fig. 4. We start the discussion with the excitation frequency located symmetrically, mid-way between the two positions of the fluctuating transition energy, i.e., ∆ = 0. The spectra for this case, for a few selected values of the switching rate β, are shown in Fig. 4(a) and Fig. 4(b) for slow and fast switching, respectively. gAs the switching rate grows, the form of the spectrum evolves from a single narrow line, via a triplet of broadened lines that eventually merge into a single line that subsequently narrows down. The quasi-static limit of β → 0 ( Fig. 4(a), black line) corresponds to the standard result for low-excitation resonance fluorescence from a two-level system, where the scattering spectrum consists exclusively of one narrow line at the spectral position of the laser [1]. Since the laser is detuned from both positions of the transition energy, the overall intensity is very weak. The subsequent evolution of the spectrum is consistent with the structure of the poles discussed above. The appearance of a single line in the fast switching regime can be interpreted again in terms of the averaging of the transition energy on time scales The intensity of the particular components of the RF spectra for the two excitation conditions, respectively, split into the slow (c,g) and fast (d,h) switching case. In the intensity plots, we group certain lines and show only their total intensity in some cases, as explained in the text.
The limiting values of the intensities of the RF spectrum components in the static and ultrafast limit, for two spectral positions of the exciting laser.
shorter than 1/∆ 0 and is therefore the counterpart of the effect observed in Fig. 2. As the averaged energy level is resonant with the laser, the scattering intensity considerably grows. The line width reduction when speeding up the noise dynamics follows from self-averaging of the fluctuations and is an example of the motional narrowing effect, by analogy to the narrowing of the nuclear magnetic resonance line for a particle that travels very fast through regions of spatially inhomogeneous magnetic field [48,49]. A quantitative understanding of the spectra is possible by combining the information on peak positions and widths, presented in Fig. 3, with the peak intensities. The latter are extracted directly by evaluating the pre-factors of the Lorentzian terms in Eq. (16) and are shown, for ∆ = 0, as a function of the switching rate in Fig. 4(c,d), where we split the result into the slow and fast noise regimes (β < ∆ 0 /2 and β > ∆ 0 /2). Although the spectrum is always positive, its decomposition into individual peaks is to some extent artificial and some of the components defined in this way may have negative amplitudes if the peaks overlap. Therefore, in some cases we group a few lines that have the same position or the same physical nature and show only the sum of their intensities, so that the presented intensities are positive. For β ≪ ∆ 0 the total scattering intensity is dominated by the nominally broadened contribution L 1 . However, as discussed above, in the limit of β → 0 its width decreases to zero, hence the fully elastic scattering is recovered in the static limit. On the other hand, for β ≫ ∆ 0 total intensity reaches the value characteristic of resonant scattering (Fig. 4(d), purple line) and is dominated by the elastic contribution L 0 (blue line), which is consistent with the resonance with the averaged transition energy, leading again to the situation known from a two level system at resonance [1]. The exact limiting values of the intensities of all the spectral components are collected in Tab. 1.
We now turn to the case when the laser is tuned to one of the two possible transition frequencies, ∆ = ∆ 0 /2. Fig. 4(e) and Fig. 4(f) show the RF spectrum in this case for β < ∆ 0 /2 and β ≥ ∆ 0 /2, respectively. In the quasi-static regime we again observe a single sharp line at laser frequency, but now the intensity is much larger than for ∆ = 0 (Fig. 4(e), black line). As the switching rate grows, this line is accompanied by two broadened lines, that initially appear around the transition energies (one of which now coincides with the laser frequency) and then merge around the central spectral position to disappear again for β ≫ ∆ 0 (Fig. 4(f), blue line) The position of the broadened features is the same as in the previous case ( Fig. 4(a,b)) with respect to the transition energies and the position of the sharp peak follows the laser frequency, while the overall intensity now decreases as the switching rate grows. In this case the laser is tuned to resonance with one of the transition energies, leading to strong scattering in the quasistatic case, which again reproduces the known result for resonant light scattering [1]. On the other hand, in the fast-switching regime, the averaged transition energy is detuned from the laser, hence in this motionally narrowed limit the spectrum corresponds to resonance fluorescence with strongly detuned excitation, showing a weak narrow line at the laser frequency (Fig. 4(f), blue line).
A quantitative analysis of the intensity of the spectral features contributing to the RF spectrum for ∆ = ∆ 0 /2 is shown in Fig. 4(g) (β < ∆ 0 /2) and Fig. 4(h) (β ≥ ∆ 0 /2). In the slow switching regime, only the spectral lines at the laser position contribute, with L 2 of negligible intensity (see Tab. 1) and L 1 becoming narrow, as discussed previously. Hence, the elastic scattering fully dominates, as expected for the quasi-static limit. However, the total RF intensity is now lower than the standard resonant scattering intensity I 0 (roughly by half), because the probability that the laser is resonant to the transition is now only 50%. For fast switching the intensities are consistent with the concept of detuned averaged transition energy, with elastic light scattering (L 0 contribution) dominating (Fig. 4h, blue line) and low total intensity.

Asymmetric switching
In this section we extend our considerations to the case of asymmetric switching, that is, β ↑ ̸ = β ↓ , or x ̸ = 0. At the beginning we discuss the total scattering intensity I tot , obtained numerically from Eq. (17) and now depending on the spectral position of the laser and on the degree of noise asymmetry, Fig. 5, analogous to Fig. 2, shows the impact of the asymmetry of the noise. As the preference for the upper position of the transition energy grows with increasing asymmetry, the spectrum gradually evolves into a single line at this spectral position. For slow noise, this happens via transferring the intensity to the right line, without changing the line positions. For fast noise, the position of the line shifts to the right without changing the intensity. As a complementary view on the same parameter dependence, Fig. 6 presents I tot as a function of ∆ and x for several values of β. At low switching rates ( Fig. 6(a,b)), the areas of high RF intensity extend around ∆ = ±∆ 0 /2, i.e., when the laser is tuned to one of the two randomly alternating spectral positions of the optical transition. As β increases (Fig. 6(c-e)), high intensity areas merge, forming finally one diagonal line (Fig. 6(f)). The intensity in the slow switching regime is a consequence of the quasi-static dynamics, with the two spectral positions of the transition occurring with the probabilities p ± = (1 ± x)/2. Indeed, Eq. (23) is generalized in this case to In the opposite limit of fast noise, the resonance appears at the averaged transition energy, where the average is now weighted by the probabilities p ± , leading to the averaged energy level of x∆ 0 /2. The total intensity in this limit is given by with I tot reaching its maximum for ∆ = x∆ 0 /2. We next analyze how the positions and widths of the spectral features change with noise asymmetry, parametrized by the parameter x. As follows from Eq. (25a) and Eq. (25b), the positions and widths of the two peaks L 0 and L 1 , located at the laser frequency, are independent of the asymmetry. The other two peaks are represented in Fig. 3 with dashed and dotted lines for two values of the asymmetry parameter x. The spectral positions of these lines (Fig. 3(a)) are again bound to the actual spectral positions of the transition at slow switching and converge towards the average frequency as the switching rate grows. However, contrary to the case of symmetric switching, they do not overlap completely but remain separated by a splitting proportional to the asymmetry parameter x. Indeed, from Eq. (25c) one finds for β ≫ ∆ 0 the peak positions ω 2,3 = ω L − ∆ ± |x|∆ 0 /2. The widths of the peaks L 2 and L 3 are shown with dashed and dotted lines in Fig. 3(b) with Γ = 0. The asymptotics of the width of the peaks L 2 and L 3 for very slow and very fast noise is the same as in the symmetric case but the behavior at intermediate switching rates is different. For small asymmetry of the noise (dashed lines) the widths evolve with β in a way similar to the symmetric case. As x grows (dotted lines), the picture changes considerably and one of these lines attains the broadening close to 2β, while the other remains narrow in the whole range of β. Fig. 7 presents RF spectra for various values of the laser detuning ∆, noise switching rate β, and noise asymmetry x. Each panel corresponds to a certain choice of ∆ and β, and compares the spectrum obtained in the presence of symmetric noise (black lines) with spectra at asymmetric noise (green and blue lines), showing how the intensities and positions of the spectral features evolve with asymmetry. In general, as the asymmetry grows to the limiting values of x = ±1, the static limit is achieved irrespective of the switching rate β, so that the spectrum evolves towards a single narrow line at the spectral position of the laser. Figs. 7(a,b) show how this happens for the central spectral position of the laser (∆ = 0). In this case the spectra are mirror-symmetric under the change of the sign of x, so only the results for x > 0 are shown. The initially symmetric spectrum first develops an asymmetry in favor of the most frequently visited spectral position, followed by the decay of the side peaks. As discussed previously, for slow noise, Fig. 7(a), the overall scattering intensity in the symmetric case is low, as the excitation is detuned from the two spectral positions of the laser, while for faster noise, Fig. 7(b), the intensity is larger as the role of the averaged spectral position increases. However, in the limit of x = ±1 the noise rate becomes irrelevant and the spectra must converge to the same limit. Hence, the intensity of the central elastic line increases in the former case and decreases in the latter. In Fig. 7(c) we show the spectra for the excitation tuned to the upper spectral position of the transition (∆ = ∆ 0 /2) and for slow noise (β/∆ 0 = 0.02). The spectrum is dominated by the spectral line at the position of the laser (composed of the lines L 0 , L 1 and L 2 ) that gains considerable intensity as x evolves towards +1, which means that the excitation becomes resonant with an increasing probability. The other spectral feature (line L 3 ) is enabled dynamically and is always very weak when the noise dynamics is slow (here we magnify it by a factor of 100). It has to vanish at x → ±1 and reaches a maximum intensity at x near 0. As can be deduced from Eq. (25c), the position of this line very weakly depends on x when β is small. The scattering spectrum for the same spectral position of the laser but faster noise (β/∆ 0 = 0.5) is shown in Fig. 7(d). Here not only the intensity but also the position of the off-resonant peak changes, in accordance with Eq. (25c). As follows from Fig. 6(e), the overall intensity in this case gradually increases as x changes from −1 to 1, which is reflected in the spectra, both for the resonant and off-resonant peaks. Ultimately, in the static limit of x → 1 (not shown), the intensity of the resonant lines increases by many orders of magnitude and the spectra in Figs. 7(c,d) reach the same form of a single, narrow, strong resonant  line.
In Fig. 8 we analyze quantitatively the intensities of the individual spectral features as a function of the noise asymmetry parameter x ∈ [−1, 1]. We restrict our discussion to the excitation tuned to the upper spectral position of the transition, i.e., ∆ = ∆ 0 /2. For the sake of clarity of the presentation, we plot the results for x < 0 and x > 0 in a logarithmic scale in Fig. 8(a) and Fig. 8(b), respectively, which reveals power-law dependence as a function of 1 − |x| as x approaches its limiting values. In the quasistatic limit of x → ±1, elastic light scattering (L 0 spectral line) dominates, as discussed above (blue line in Fig. 8). Obviously the intensity of scattering differs by orders of magnitude in these two limits, as they correspond to strongly detuned and resonant excitation, respectively. In particular, I tot → I 0 when x → 1 (Fig. 8(b)). In a wide range of intermediate values of noise asymmetry, the noise-induced inelastic scattering dominates (red and green lines in Fig. 8). The inelastic side line (L 3 , red line in Fig. 8) has its maximum for a slightly asymmetric noise.

N-source telegraph noise
In this section we present the results for a system subject to noise originating from N identical additive sources, restricting the discussion to symmetric switching. The total scattering intensity as a function of the spectral position of the laser is shown in Fig. 9 for N = 5 and N = 100. These results were calculated by numerical evaluation of Eq. (17). In the slow-switching regime, when scanning the laser frequency, we observe N + 1 resonant frequencies (see Fig. 9(a)). For large N , these resonances form a broad spectral feature, with the maximum intensity for the laser tuned centrally (∆ = 0) ( Fig. 9(b)). As the switching rate grows, the resonances merge into a single motionally narrowed line. The appearance of multiple resonances in the slow-noise regime is obviously related to the N + 1 positions of the transition energies in this case. In the quasi-static limit (β → 0), the matrix B becomes diagonal and Eq. (17) trivially yields a series of Lorentzian features weighted by the probabilities p (st) that follow the binomial distribution according to Eq. (21). As a result, the envelope of these resonances forms an approximately Gaussian line (by virtue of the standard Gaussian approximation of the binomial distribution) with a width of ∆ 0 , which is a consequence of our choice to renormalize the noise amplitudes by √ N in Eq. (19). For β ≫ ∆ 0 the resonances merge, like in the previously discussed case of N = 1, forming a single narrow Lorentzian resonance at the average transition energy.
The RF spectrum for different number of noise sources, calculated numerically from Eq. (15), is shown in Fig. 10. For slow switching and small N , the RF spectrum has N + 1 visible side peaks and the central peak ( Fig. 10(a), black and red line). As N increases, the side peaks start to overlap and form a broad feature centered at the laser frequency (blue line in Fig. 10(a)). At N = 50 the spectrum has reached its asymptotic form and does not change when the number of sources is increased further (which is again due to the normalization of noise amplitudes assumed here). In the fast noise regime, the side peaks are merged into a single feature, as in the previously discussed case of a single source, and there is no visible dependence on N ( Fig. 10(b)). The central feature corresponds to the first term in Eq. (16) and is composed of N + 1 Lorentzians localized at the laser frequency with widths 2nβ, where n = 0, 1, 2, . . . , N . In the slowswitching regime, the remaining part of Eq. (16) yields N +1 Lorentzian side peaks with an approximately Gaussian envelope at small β. As in the single-source case, they can be related to the spectral positions of the transition energy resulting from the states of the noise sources. For β ≈ ∆ 0 the spectrum is restructured and all peaks are localized at the laser frequency.

Conclusions
In this paper we have studied resonance fluorescence from a two-level system subject to classical external noise that leads to fluctuations of the transition energy. We have formulated the general theory of the resonance fluorescence spectrum in the low-excitation regime in the presence of noise that can be described as a stationary Markovian random process on a finite state space, which can approximate an even wider class of Markovian processes. We have also described the resonance fluorescence spectrum under uncorrelated noise leading to phase diffusion. Our theory relates the light scattering spectrum to the formal characteristic of the stochastic noise process.
We have applied our theory to the cases of a single two-state noise source (random telegraph noise) and an arbitrary number of identical sources, where many characteristics can be extracted in an analytical form. Our results show essential differences not only between the phase diffusion and random-telegraph-like processes but also between the regimes of slow and fast dynamics of the random telegraph noise. Most remarkably, the resonance fluorescence spectrum changes its form from multiple spectral features or a broad Gaussian feature (depending on the number of noise sources) to a single motionally narrowed line as the noise dynamics gets faster. In this way, we have demonstrated that light scattering on a two-level system in a noisy environment can yield information on the character of the noise processes experienced by the system. These findings may be useful in particular for interpreting experiments on the inherently noisy solid-state systems, where resonance fluorescence finds a constantly growing range of applications. In these systems, the typical lifetimes γ −1 are in the nanosecond range, setting the limit on the field amplitudes for which our lowexcitation theory is valid (Ω ≪ γ). The noise induced by electrical or spin environment [24,25,22], is typically slow compared to the dynamical time scales of the system. On the other hand, a carefully designed optical experiment [23] shows the coexistence of slow (nanosecond time scale) noise with a fast noise component, on picosecond or shorter time scales, which may be due to lattice vibrations. This might open the path to direct verification of our theory. One must note, however, that the transition between the slow and fast regimes in our theory is controlled by the ratio of the noise dynamical rate β and noise amplitude ∆ 0 . The former is an inherent feature of the noisy environment, why the latter may only be modified by engineering the coupling between the emitter and the environment. It may therefore turn out that the most straightforward way to validate the theory would be to use artificially generated mechanical noise, taking advantage of the high flexibility of mechanical signal generation and controllability of the acoustic coupling to solid-state emitters [27,26].

Acknowledgments
This work was supported by the Polish National Science Centre (NCN) under Grant No. 2016/23/G/ST3/04324.

Data availability
Any data that support the findings of this study are included within the article.

Appendix A. Averaging over the random process
In this Appendix we present the technical details of the averaging in Eq. (13). We start from the basic formula for averaging an arbitrary function of the process state at a finite set of time instants, f (∆(t 1 ), ..., ∆(t n )) = k j 1 =1 ... k jn=1 p jn,j n−1 (t n −t n−1 )p j 2 ,j 1 (t 2 −t 1 )p Here we take advantage of the fact that the process is Markovian and stationary, hence the joint probability can be written as a chain of conditional (transition) probabilities with the initial probability distribution assumed to be the stationary probability distribution of the system and the transition probabilities depending only on the time difference.
To find expectation values of exponential terms, e iΦ(t d ,tc)+iηΦ(t b ,ta) , η = ±1 we divide the time intervals (t a , t b ) and (t c , t d ) into N and N ′ pieces, respectively. Then where N δs = t b − t a , N ′ δs ′ = t d − t c , t j = t a + j − 1 2 δs, and t ′ l = t c + l − 1 2 δs ′ . Upon using Eq. (A.1), one gets e iΦ(t d ,tc)+iηΦ(t b ,ta) = lim where q (η) l,j (∆t) = e iηx l ∆t p l,j (∆t). This formula can be rewritten in terms of a product of transition matrices P τ and k × k matrices Q(τ ) with matrix elements q where the conjugation in the last term corresponds to the '−' sign on the left-hand side.