Time resolved multi-photon effects in the fluorescence spectra of two-level systems at rest and in motion

We study the time-resolved fluorescence spectrum in two-level systems interacting with an incident coherent field, both in the weak and intermediate coupling regimes. For a single two-level system in the intermediate coupling case, as time flows, the spectrum develops distinct features, that are not captured by a semi-classical treatment of the incident field. Specifically, for a field on resonance with the atomic transition energy, the usual Mollow spectrum is replaced by a four peak structure, and for a frequency that is half of the atomic transition energy, the time-dependent spectrum develops a second harmonic peak with a superimposed Mollow triplet. In the long-time limit, our description recovers results previously found in the literature. After analyzing why a different behavior is observed in the quantum and classical dynamics, the reason for the occurrence of a second harmonic signal in a two-level system is explained via a symmetry analysis of the total (electron and photon) system, and in terms of a three level system operating in limiting regimes. We find an increased second harmonic signal in an array of two-level systems, suggesting a superradiance-like enhancement for multiple two-level systems in cavity setups. Finally, initial explorative results are presented for two-level model atoms entering and exiting a cavity, which hint at an interesting interplay between cavity-photon screening and atomic dynamics effects.


I. INTRODUCTION
Fluorescence, a type of luminescence [1][2][3], is a hallmark of quantum mechanics at work: A system that has absorbed electromagnetic radiation re-emits it at a later time, while the spins of the electrons involved in the deexcitation process conform to specific selection rules.
In addition to being an operational mechanism in several biological systems [4], fluorescence serves in many technologies of different complexity, ranging from simple indoor lighting to in-depth spectroscopic characterization at the atomic scale. Investigations of fluorescence started before the advent of quantum mechanics, but it would be the latter that provided the conceptual framework for a microscopic description [1][2][3].
Spectroscopic methods have a rich history as a means of investigating the internal structure of matter. In particular, with the development of ever more sophisticated laser systems, the electronic dynamics in atomic and solid state systems can now be mapped out in real time while maintaining a high frequency resolution [5,6]. This allows to study in a precise manner the basic processes of light-matter interaction, and even to characterize properties of light itself, as for example emitted via fluorescence.
A minimal-complexity model to study fluorescence emission and fluorescence spectroscopy is a two-level system [7][8][9] (for example, a spinless electron that can be in either of two nondegenerate quantum states, or in any linear superposition thereof) interacting with a single radiation mode via dipolar coupling. If the two levels are thought to be selected from an atom, the model is also referred to as a "two-level atom" (with the additional option that the center of mass of the atom can be either at rest or in motion).
Two-level systems came into prominence with Rabi et al.'s work for a magnetic moment exposed to a classical circularly polarized field [10]. In a subsequent study by Bloch and Siegert the linearly polarized case was then considered [11] (in this situation, the solution is more complicated compared to Rabi's original case [12]). The next important development took place when the radiation mode was also treated quantum-mechanically, and the so-called rotating wave approximation (RWA) was introduced [13][14][15][16][17]. Designed for weak-coupling and nearresonance regimes, the RWA permits an explicit treatment of the time dependence [18], and provides a convenient route to the so-called dressed-level (or -atom) approach [9,19,20], where the levels of the system split and are renormalized (shifted) by the radiation field. In turn, within this approach a clear picture emerges [21] of the Mollow spectrum [22], a three-peak structure due to the fluorescent response of a two-level system to a resonant or quasi resonant radiation mode.
Nevertheless, a number of interesting physical situations are outside the reach of RWA, as e.g. the intense pulsed regime, where field monochromaticity is absent and off-resonant coupling cannot be avoided. The need under some circumstances to go beyond the RWA has in fact been recognized in several contexts (also by comparing exact and RWA solutions [23,24]). For example, when discussing modifications of the shape of the three peaks in Mollow spectra [25][26][27], spontaneous emission in three-level systems [28], or when center of mass dynamics is included [29] (for a recent review, see e.g. [30]).
As these few, incomplete remarks suggest, two-level systems coupled to radiation in different "flavors" remain of capital relevance to this day to probe and redefine the knowledge boundaries in (quantum) optics [31][32][33]. This can occur via generalization of the basic model(s) together with deeper mathematical analysis (see e.g. [34][35][36]), to address unexplored coupling regimes [37,38], or novel areas of applications. For example, cavity quantum optics and the Unruh effect [29], quantum mechanical interference [39], two-photon relaxation [40], quantum phase transition [41], interaction of a photons matter qubits [42] and Mollow spectra in ultracold atoms [43].
The quantum nature of light manifests in a clear way at low photon number and for large light-matter coupling. These two aspects contribute distinctly. This is different from the strong-field regime, where a semiclassical treatment becomes appropriate and where the effective coupling parameter is the product of field strength and coupling strength [12]. Concerning the few-photon limit, this can be e.g. reached in high quality-factor cavities [44][45][46]. On the other hand, to attain the strong coupling (also denoted polaritonic) regime, a possibility is offered by the insertion of a quantum well into a distributed Bragg reflector cavity [47], i.e. by coupling the photon field and an optical inter-band transition (a Wannier exciton).
Scope of this work.-In the present study we consider two-level, one-electron systems interacting with two optical modes (the coherent-pump and a de-excitation field). Specifically, we address the (so far largely unexplored) multi-photon effects in fluorescence spectra, that depend separately on field intensity and light-matter coupling strength. This will be done in situations of progressive complexity: a single two-level atom at rest, an array of two-level atoms at rest, and finally a single two-level atom moving through an optical cavity. To this end, it is necessary to employ a theoretical framework suitable for both non-linear effects (in relation to certain experimental setups [48][49][50]) and an explicitly time-dependent light-matter coupling.
Several years ago, three of the present authors introduced an exact solution method [51,52] for a large class of multi-photon spectroscopy models. Their method is based on a recursion technique (see e.g. Ref. [53]) in the frequency domain, i.e. for the stationary limit of fluorescence. They also pointed out the difference between the exact solutions and the RWA. The aim was to address systems where interactions with the environment are as weak as possible. That is, where both inhomogeneous (such as due to static and dynamical disorder, mode leakage, etc. ) and homogeneous (atomic collisions, but also non-radiative decay, etc.) decoherence factors plays a minor role. Put differently, the focus was on a regime where both energy-dependent broadening (not considered in the rest of this paper, because assumed to be controllable), and energy-independent one (denoted by Γ and retained in the paper) are as small as possible. As discussed above, current experimental capabilities provide practical and close-to-ideal realisations of these premises with optical cavity setups, and make time-resolved studies possible.
Schematic of the systems considered in the paper. In panel (a), a single two-level system of transition energy interacting with a coherent field of frequency ωa (with coupling strength ga) and a fluorescent field of frequency ω b (with coupling strength g b ). In panel (b), an array of N two-level systems interacting with a coherent field and fluorescent field. In panel (c), a two-level atom of momentum p passing through a cavity of length L where it interacts with a coherent cavity field and emits fluorescent photons.
In this way, it is possible to investigate the actual development of the fluorescence signal before the steady state signal sets in. With these considerations in mind, here we take a different methodological route from that in Refs. [51,52], by working in a real-time (and again free of RWA) picture, and computing the exact time-dependent fluorescence response. To establish the effect of multi-photon contributions and counter-rotating terms, we specialize to the Mollow regime (on-resonance situation) and to secondharmonic generation (SHG) (off-resonance situation). In the literature, Mollow spectra are discussed in terms of two level systems [22], and thus our work conforms to previous treatments. In contrast, SHG is commonly discussed in terms of three level systems [54]; however, we will show that a genuine two-level system admits SHG. Overall, our study here can thus be summarized as an exploration of multi-photon effects in Mollow and SHG fluorescence spectra across three different two-level-system setups (a single system at rest, many systems at rest, a single system in motion).
A remark about the units used in this work: unless otherwise stated at specific points in the paper, the energy unit is = 2 − 1 , the distance between the two levels in the system, and the time unit is / [55].
Organisation of the paper.-The rest of this article is organized in three parts: In Sec. II we study the timedependent fluorescence spectrum of an isolated two-level system (see Fig. 1a). In Sec. III, we consider an array of two-level systems interacting with a common coherent field (see Fig. 1b). In this case, the fluorescence signal shows an enhancement compared to the single two-level system case, consistently with a superradiance-like mechanism. Finally, in Sec. IV we consider a two-level model atom passing through an optical cavity (see Fig. 1c). We explicitly treat the quantum motion of the atomic cen-ter of mass, and show that this results in a fluorescence spectrum which differs from that obtained from a semiclassical treatment of the atomic motion. Conclusive remarks and outlook are in Sec. V.

II. A SINGLE TWO-LEVEL SYSTEM
We start with the simplest of the three situations, i.e. a single two-level system interacting with two optical modes. After presenting model and method of solution, we investigate the time evolution and the long time limit of the system's fluorescence spectrum in the Mollow regime. We consider the case when the Rabi frequency g becomes a moderate fraction of the level spacing of the two-level system (about 10 %). The results show that for a field in resonance with the atomic transition energy, the usual three peak Mollow spectrum is replaced by a four peak structure. Instead, when the frequency is half of the atomic transition energy, we obtain an SHG spectrum with a superimposed Mollow splitting. The emergence of an SHG signal in a two-level system is anticipated by the analysis of the spectrum of a three-level system, and further validated by a symmetry analysis of the coupled electron-photon states.

A. Model and method
Our two-level system interacting with an incident and a fluorescent light-field [51,52] is described by the following Hamiltonian: We assume that the atom is occupied by a single spinless electron, so that Hereĉ i destroys an electron in the orbital |i with energy i ,σ z is the z-component Pauli operator, and = 2 − 1 .
The free radiation modes are described by the Hamilto-nianĤ whereâ annihilates a photon of the incident field with frequency ω a , which we assume to be in a coherent state defined byâ|α = α|α . Similarlyb annihilates a photon of the fluorescent field with frequency ω b . The lightmatter interaction Hamiltonian is given bŷ where g a (t) and g b (t) are the (time-dependent) couplings of the electron to the incident and fluorescent fields respectively, and can have any time-dependence. Also,σ x is the x-component Pauli operator. In the following we consider the case g b (t) = g b e −Γt , which introduces a frequency independent, phenomenological damping of rate Γ. The use of Γ takes into account in a qualitative way effects left out, e.g. non-radiative transitions and/or mode leakages in a cavity geometry (see also Appendix A 1). Concerning the role of spontaneous decay in our description, we note that in the stationary regime (e.g. due to a steady photon pump) it is often legitimate to overlook this type of decay with respect to the stimulated one. Away from stationarity, other factors come into play, depending on the situation: i) in the single-atom case (and with Einstein's description of radiation-matter interaction as conceptual reference), spontaneous decay does not induce a thermal bath; here, the coherence of the overall optical response is not altered by neglecting such decay (see e.g. [56]) ii) in the the many-atom case (e.g the Dicke's regime as discussed in Sect. III), thermalbath effects and the significance of spontaneous decay are hindered by the collective effect of super-radiance [57].
To describe the dynamics according to Eq. (1), we use the exact configuration interaction method. In this way, the full wavefunction of the coupled atom-light system is represented in the basis |i, n, m ≡ |i |n |m , with |i the state of the atom (where |1 is the ground state and |2 the excited state), |n a number state of the incident field and |m a number state of the fluorescent field. We start from the initial state |ψ 0 = |1, α, 0 (here, the number state |n has been replaced by the coherent state |α of the pump field), and time-evolve it with the full Hamiltonian H using the short iterated Lanczos technique [58,59] (see also Appendix A 2).
We note here that the recursion method originally employed to study fluorescence in the stationary limit [51,52] can be used for more general setups, e.g. to generate exact solutions with several electronic levels and several bosonic modes [53]. However, its applicability is not immediate for genuinely time-dependent Hamiltonians, and to address the transient response of a system. Hence, the need to proceed here with a real-time approach.
It has been shown [60] that to make contact with experimental time-resolved light spectra, the transition probability needs to be convolved with the resolution function of a Fabry-Perot spectrometer. Also, depending on the experiment performed on a quantum system, a detection of N -photon correlations can be used [61].
However, here we employ a definition of the spectrum different from the one considered in [60,61]. Namely, we consider the probability P that at least one photon (at a given frequency) is emitted. Such a spectrum would not be obtained by an interferometer, Fabry-Perot or otherwise, but rather by an apparatus including e.g. a prism and a photomultiplier. In other words, this we would correspond to measure particle-like photons, rather than waves, and probing simultaneously the atomic state.
Accordingly, our observable of main interest is the probability to find m photons in the fluorescent field of frequency ω b at time t, given by where the dependence on ω b in RHS is implicitly contained inĤ(t). Since g b (t) has an exponential decay, P will be independent of t in the long time limit. In the following we will focus on the quantity P = m>0 P m , giving the probability that at least one fluorescent photon has been emitted. The semi-classical limit of our model is obtained by taking α → ∞ and g → 0, while keeping gα constant. In this limitĤ r = ωb †b and the interaction Hamiltonian iŝ The t → ∞ limit of the model has been studied in previous works [51,52]. In that case, the Hamiltonian was split asĤ(t) =Ĥ 0 (t) +Ĥ (t) withĤ (t) = g b (t)(ĉ † 1ĉ 2 +ĉ † 2ĉ 1 )(b † +b), and the probability P 1 was evaluated under the assumption thatĤ (t) acts only once during the time evolution. Physically this corresponds to a first-order treatment of the fluorescent field, in contrast to the exact numerical solution of the present work, which retains the effects of the interaction at all orders.
B. Two-level system and Mollow spectra In the following we let the transition energy = 1 define our unit of energy, and further fix the parameters g b = 0.01 and Γ = 0.02. We start by studying the fluorescence spectrum for ω a = 1 as a function of α and g a , keeping the product g a α fixed. The results are displayed in Fig. 2, and show both the asymptotic spectrum (as t → ∞) and the explicit time-evolution. For α = 5 and g a = 0.02 we see the well-known Mollow spectrum, that is qualitatively reproduced by the semi-classical approximation. Keeping the product g a α fixed and taking α = 1 and g a = 0.1 the semi-classical result is unchanged, while the full quantum treatment gives a fluorescence spectrum with four peaks and additional substructure.
To understand the spectra we consider the states |i, n with the atom in state i and n photons in the incident field. For ω a = the levels |1, n and |2, n − 1 are degenerate, and mixed by the light-matter interaction. Diagonalizing the Hamiltonian in this subspace (i.e. neglecting the counter-rotating terms), we find the energies For a coherent state with large α the Poisson distribution is sharply peaked around α 2 , so that n ∼ α 2 and the energies differences between successive n are In this limit the energy splittings are independent of n, and transitions between successive levels give a threepeaked structure. For small α, this is no longer the case, since the splitting +,n − −,n = 2g a √ n. The nonuniformity in the level spacing is a clear sign of the quantum nature of the light (the low photon number limit), and is e.g. responsible for photon blockade effect [62,63]. In the present context, it gives rise to the additional features observed in the fluorescence spectrum of Fig. 2.
C. Prelude to SHG in a two level system: frequency doubling in an ordinary three-level system Before addressing SHG in a two-level system, we make a detour into the more familiar theory of SHG in threelevel systems. Let the Hamiltonian beĤ(t) =Ĥ e +Ĥ r + H i (t), withĤ r as in previous sections, and the electronic Hamiltonian bê For the light-matter interaction we consider two scenarios, illustrated in Fig. 3. In the first the incident field couples to the transitions |1 ↔ |2 and |2 ↔ |3 , and the fluorescent field to the transition |1 ↔ |3 . The interaction Hamiltonian is then In the second case the incident field also couples to the transition |1 ↔ |3 , with a strength g(t), and the Hamiltonian iŝ If the levels |1 and |3 are of definite and different symmetry, while the level |2 is assumed to be of mixed symmetry, both these models allow SHG using a perturbative treatment. However, in the case of a parity-invariant electronic Hamiltonian, where all electronic states have definite parity, both models forbid SHG in a perturbative approximation. The reason is that there is no way to arrange the parities such that both the excitation and emission steps are allowed: the parities π 1 and π 3 of levels |1 and |3 need to be different (for the fluorescent transition to be allowed), while the parity π 2 of level |2 needs to be different both from π 1 and π 3 (in order for the exciting transitions to be allowed).
We now study the fluorescence spectrum in the nonperturbative limit. We take 2 − 1 = 0.5 and 3 − 1 = 1, and let the incident field be resonant with the transition energies (ω a = 0.5). We choose the couplings as f (t) = g a (t) = 0.1, take α = 1, and as before let g b = 0.01 and Γ = 0.02. The results are reported in Fig. 3

. For the HamiltonianĤ
(1) i we see a broadened SHG peak centered around ω b = 1, while forĤ (2) i the spectrum has two contributions corresponding to Rayleigh scattering and SHG centered at ω b = 0.5 and ω b = 1 respectively. As expected, both models predict a non-zero SHG signal.
A connection with a simpler two-level system can be made by taking the limit 2 → ∞, illustrated here by considering 2 − 1 = 1, 1.5 and 2 (see Fig. 3). With the interactionĤ (1) i the fluorescence signal narrows around ω = 1 as 2 − 1 is increased, since the coupling between degenerate states causing the broadening decreases as the levels are energetically separated (cf. the discussion of the Mollow spectrum). For even larger values of 2 − 1 (not shown), the peak tends to zero since excitation of the atom becomes increasingly unlikely. With the interaction H (2) i both the Rayleigh and SHG peaks narrow as 2 − 1 is made large. However, due to the presence of the coupling between the incident field and the transition |1 ↔ |3 , a finite fluorescence signal remains even for 2 − 1 → ∞, i.e when effectively the HamiltonianĤ (2) i collapses onto that of the two-level system. The specific shape of the SHG profile is discussed in the next Section.
As a final consideration, we note that in quantum optics it is sometimes useful to perform an adiabatic elimination of intermediate levels, to end up with a reducedspace effective Hamiltonian [64,65]. In appendix B, such a reduction is performed to map the three-level system Fluorescence spectra for a three-level system interacting with a coherent field with a frequency ωa = 1 and an average number of photons α 2 = 1. The coupling between the system and the coherent field is given by f = ga = 0.1, while the coupling to the fluorescent field is g b = 0.01. The energy levels of the atom are 1 = 0 and 3 = 1, with 2 taking the values 0.5, 1, 1.5 and 2. In panel (a) the coherent field is coupled to the transitions 1 ↔ 2 and 2 ↔ 3, while in panel (b) the coherent field in addition couples to transition 1 ↔ 3. The couplings in the system are indicated schematically in the panels to the left. For ease of reading, the z-axis has been multiplied by 100. The unit of energy is given by = 3 − 1.
of Fig. 3 onto a two-level one [66], and to see if SHG can be exactly described in a two-level system in the limit pling the semi-classical calculation is seen to be in good agreement with the quantum treatment, while for large coupling the results differ by the presence of a superimposed Mollow spectrum on the second harmonic peak. This additional feature can be understood similarly to the Mollow spectrum of Fig. 2, as due to an energy level splitting depending on the coupling strength g a and photon number α 2 . The energy levels are shown in Fig. 5, and we see that for small g a and large α the splittings are small and uniform (between successive levels), while for large g a and small α the splittings are large and nonuniform. We note that for small g a , in contrast to the Mollow spectrum discussed above, the level splitting is expected to be ∼ g 2 a to leading order, since the coupling between the levels is of second order in the interaction. To summarize, i) the ability of the quantum treatment to discriminate between photon number and coupling strength (they always appear via their product in the semiclassical treatment) and ii) the photon fluctuations at low photon number, are likely the reasons for the Mollow structure in the quantum SHG signal. This is in line with the evidence provided by the three-level system results in Fig. 3 where, depending on the closure of the "three-level triangle" via the pump field, the Mollow structure in the SHG peak is observed or not.
For the results just presented the photon energy ω a of the coherent field is commensurate with the atomic transition energy, i.e. ω a = /n. Although our approach is completely general and we can in principle study any frequency, a systematic scan of the parameters is outside the scope of the of the present work. However, we mention that according to additional calculations (not shown) the spectra for frequencies ω a / ≈ 0.3 − 0.8 resemble the second harmonic spectrum above, although with the position of the Rayleigh peak displaced to ω b ≈ ω a . Instead, for frequencies ω a / ≈ 0.9 − 1.1 the Rayleigh and harmonic peaks blend together to form the four peak Mollow structure discussed above.
Parity conservation.-Although borne by the exact numerical results of Sect. II D (and also supported by the connection to three-level physics as discussed in Sect. II C), the occurrence of a SHG in a two-level system might remain at some extent counter-intuitive. For second harmonic generation (SHG) to occur, two photons are needed to excite the atom. This requires atomic levels of equal parity, since in the low intensity limit where the absorption process is well described by perturbation theory, the two-photon absorption probability otherwise vanishes. However, the emission of a double frequency photon requires the atomic levels to have opposite parity. This apparent paradox however vanishes at stronger light-matter coupling, where a more appropriate description of the system is in terms of dressed atomic levels. In this regime the electronic states are mixed by the coupling to the light field, and no longer have definite parity.
To understand how SHG can happen in a two-level system, we look at the parity of the eigenstates ofĤ. A parity operator can be defined through [68] wheren a andn b are the photon number operators of the incident and fluorescent fields respectively. It is straightforward to show that [Ĥ,Π] = 0, from which it follows that the eigenstates ofĤ can be classified according to the eigenvalues ofΠ. The general structure of the eigenstates with the electron and the incident field coupled is [68] |Ψ k e = n c k 2n |1, 2n + n c k 2n+1 |2, 2n + 1 (13) where e and o denote even-and odd-parity states respectively. We see that an eigenstate with a well-defined parity in the coupled system has in general an undefined parity in the electronic and photonic subspaces. Therefore an argument against SHG in the two-level system, which relies on the conservation of only electronic parity is not applicable. Consider now the evolution of a system starting from the initial state |Ψ = |1, α . Since the coherent state is a superposition of number states, it explicitly breaks parity symmetry. We therefore expect the time-evolution to induce transitions between the initial state and all eigenstates allowed by energy conservation, which in particular means that an electron may be excited to the upper atomic level. However, even for an initial state of definite parity (as for example |Ψ = |1, n ) the time-evolution will mix eigenstates, but now in a definite parity sector. Thus, in stark contrast with the predictions of perturbation theory, an SHG signal also occurs for parity conserving dynamics. In addition, and differently from a three-level system, the closure of the multi-photonic triangle (which seems needed for the appearance of the three-peaked Mollow structure), always occurs in a two level system, as an emerging symmetry mixing behavior. Furthermore, Mollow-like overtones can also be present in higher-order harmonics [51,52], a feature that could become relevant in the time domain (e.g. in the ultrafast regime). These summary notions provide rigorous foundation and motivation to investigate SHG in more complicated two-level system setups and out of the stationary limit, i.e. when the effects discussed so far will appear within characteristic timescales (for example, when investigating an atom moving through a cavity of finite length). Accordingly, time-resolved multi-photon fluorescence is the main theme of the rest of the paper.

III. AN ARRAY OF TWO-LEVEL SYSTEMS
We now consider N two-level systems interacting with one incident and one fluorescent light mode. This brings in the possibility that the different two-level systems interact cooperatively with the radiation, with an enhancement of the radiation field that goes under the name of superradiance [69] (more precisely, depending on the nature of the initial state, one can speak of superradiance or supercoherence [70]). The standard Dicke model [69], which describes N twolevel systems interacting with a single optical mode, plays an archetypal role in the study of superradiance and has been the focus of extensive investigations (see e.g. [71][72][73][74]). An exact solution within the RWA of the Dicke model at resonance has been known for quite some time (for this reason, the model is also referred to as the Tavis-Cummings model [75]) but some aspects related to this model are still being debated [76]. A notable case concerns the existence of a no-go theorem for superradiance, also in connection to the role of dipolar couplings between electrons and photons and the interactions among the different two-level systems. Furthermore, consideration is needed for the superradiant behavior when N → ∞, i.e. the periodic array limit. Mathematically, the intensity diverges in that situation, but in fact a continuous electronic band structure emerges in the thermodynamical limit, and the bulk polaritonic regime applies, with a finite optical response [77].
These different aspects are not addressed here. Rather, our simple analysis is aimed to gain qualitative insight into how the behavior of the single two-level system discussed above changes in a more complex setup. As before, we position ourselves in the fully-resonant and the SHG regimes for the incident field. We will consider the problem from two different but complementary perspectives: (i) We first study the equilibrium properties of the system in the large-N limit using a generalization of the Dicke model. (ii) Secondly, we investigate the nonequilibrium properties of the system for finite N (similar to what done earlier for N = 1) using an exact numerical treatment. We find that, compared to the single two-level system case, the fluorescence signal shows an enhancement compatible with a superradiance-like mechanism, both in the resonant and SHG regimes.
As an obvious extension of the N = 1 case considered earlier, the total Hamiltonian becomes Hereĉ 2i−1 refers to the ground state andĉ 2i to the excited state of two-level system i. Due to the specific way the fields couple to the two-level systems, we can rewrite the Hamiltonian in a more compact form aŝ where ω i = 2i − 2i−1 . When all excitation energies ω i = ω s , the fields only couple to the total spin operatorsŜ z = σ z,i andŜ x = σ x,i . The Hamiltonian then closely resembles the Dicke Hamiltonian [69], and is written A. Ground state We start by discussing the ground state properties. For an initial state with all two-level systems in their ground state, we can further simplify the Hamiltonian using the Holstein-Primakoff transformation [78]. Writing the spin operators in terms of bosonic creation and annihilation operators,Ŝ z = −(N/2)+ŝ †ŝ andŜ + = s † √ N −ŝ †ŝ , it is straightforward to check that the commutation relations defining the spin algebra are preserved. Taking the limit where N ŝ †ŝ , theŜ + operator can be expanded aŝ S + ≈ s † N + O(N −1 ), and the Hamiltonian becomeŝ (21) and the system maps into three coupled oscillators. For g b = 0, the Hamiltonian of Eq. (21) reduces to the usual Dicke Hamiltonian, where a transition to a superradiant state takes place at a critical coupling g a ≈ (2N ) −1 , with the field intensity of the radiation becoming proportional to N 2 [76]. To see this, we follow closely and reproduce here the discussion given in Ref. [76], starting by writing the Hamiltonian in terms of canonical coordinates: where, to keep the notation light, we have introduced the coupling parameter λ = N g a . We also assume that the field is resonant with the two-level systems, i.e. ω a = ω s = 1. To identify the superradiant transition, we look for the point λ c where the number of photons in the lowest normal mode ofĤ diverges. The lowest mode is identified by moving to normal canonical coordinates and momenta ofĤ, given byx ± = (x s ±x a )/ √ 2 andp ± = (p s ±p a )/ √ 2, respectively, and with normal frequencies Ω 2 ± = 1 ± 2λ. In the ground state, the expectation value ofx 2 ± is x 2 ± = (2Ω ± ) −1 , which diverges at coupling λ c , whilst p 2 ± remains well behaved. We then go back to the original coordinateŝ x a/b ,p a/b , expressing them in terms ofx ± ,p ± , and note that n a = â †â = ( p 2 a + x 2 a − 1)/2. Then, when λ → λ c , the diverging contribution comes from x 2 ± , i.e.
which diverges for |λ| = λ c = 1/2. Thus, for λ > 0 we can identify the ground state superradiant transition as the point where the lowest eigenvalue Ω − vanishes. To extend the above discussion from Ref. [76] to the original Hamiltonian Eq. (21) with the fluorescent field included, we writê where again λ a = N g a and λ b = N g b . At this point, if we take ω s = ω a = ω b = 1, the lowest normal mode eigenvalue Ω 0 of H is given by which is the generalization to the case of two fields both in resonance with the N two-level systems. If instead, to address the SHG regime, we assume that ω s = ω b = 1 and ω a = 1/2, we obtain where A = α − α 2 − 4β 3 , α = 1/32 + 9λ 2 a − 18λ 2 b , and β = 1/16 + 12λ 2 a + 12λ 2 b . As before, we search for signatures of a superradiant transition by looking at the points where Ω 0 vanishes. In the resonant case this is easily done, and one obtains the semi-circle solution set λ a = 1/4 − λ 2 b . We note that for either λ a = 0 or λ b = 0, we recover the critical coupling of the standard Dicke model discussed above. In the non-resonant case, we solve instead Eq. (26) numerically to find the result in Fig. 6a. Interestingly, we find that the signature of a superradiant transition occurs for smaller values of the coupling when one of the fields is non-resonant with the atomic transitions.
The above argument indicates that the lowest eigenstate of the Hamiltonian undergoes a superradiant transition. However, we have not yet determined how this state is related to the original photon fields, and therefore at this point it is not clear how these fields behave at the transition. In analogy with the Dicke model, the average number of photons in each of the original light fields (in the ground state) is proportional to |c a/b | 2 (2Ω 0 ) −1 , where |c i | 2 = | i|Ω 0 | 2 are the projections of the normal mode |Ω 0 onto the original oscillators. The fields should therefore undergo a transition to a superradiant state for values of λ a and λ b such that (i) Ω 0 = 0, and (ii) |c a/b | 2 are finite. The coefficients are found by numerical diagonalization of the Hamiltonian in Eq. (24), and shown in Fig. 6b for values of λ a and λ b at the superradiant transition. It is apparent that for λ a/b > 0 we always have |c a/b | 2 > 0, so that condition (ii) above is always satisfied. We thus find that at the transition points found above, both the incident and the coherent fields behave as if a superradiant state is attained.

B. Real-time simulations for finite N
Having discussed some ground state features of a twolevel system array coupled to radiation, we now return to explore the time evolution of the system defined by Eq. (19). We start from an initial state with all atoms in their ground state, and therefore we do not expect to see a superradiant emission burst. However, we are interested in exploring how the fluorescence spectrum changes as we approach coupling strengths close to the (equilibrium) superradiant transition. We consider the parameters ω s = ω b = 1, ω a = 0.5, and choose |α| 2 = M = 9 for the average number of photons in the cavity. With M of the same of order as the number of atoms N , the energy of the field should be enough to simultaneously excite all the atoms. This means that for a two-level array with bare couplings given by e.g. g a = 0.03 and g b = 0.01 (to be used below in the actual simulations), the minimal value of N needed to get the critical values of λ a and λ b indicated in Fig. 6 is given by N ≈ 14. This is done by checking for which value of N that λ a = N g a and λ b = N g b cross the blue line in Fig. 6a.
Since the arguments of the previous section are only strictly valid in the ground state and for large N (due to the use of the Holstein-Primakoff transformation), there is no guarantee that they would hold in real time. Thus, our estimate for N just provides a hint of the order magnitude of N where we can expect superradiant effects to appear. In addition, the non-equilibrium signatures of superradiance are typically expressed through the scaling of the duration and intensity of the superradiant burst with the number of two-level systems N , given respectively by 1/N and N 2 . We therefore consider below the fluorescence spectrum for N ∈ {1, 2, . . . , 10}, and look for signatures consistent with these scaling laws.
In Fig. 7 we show the fluorescence spectra for an array of N = 2 and N = 10 two-level systems. We find that with increasing N , the time it takes for the peaks to develop is reduced, as indicated by the vertical lines in the figure. For N = 2 the resonant Rayleigh peak develops before the second harmonic peak, while for N = 10 the order is opposite. In addition, for N = 10 the SHG peak transiently exceeds the Rayleigh peak also in magnitude.
To quantify these observations, we define T as the time it takes a peak to reach half its maximum value. As shown in Fig. 8, we find that T as a function of N shows a crossover around N = 2, from a regime where the Rayleigh peak develops first to a regime where the SHG peak comes first. Further, we find for the Rayleigh peak that the dependence of T on N is approximately linear, while for the SHG peak is behaves as 1/N . For the latter case the scaling is consistent with a superradiant behavior, where the duration of the superradiant burst decreases as 1/N . In Fig. 8 we also show the number of photons P(ω b = 1) in the fluorescent field as a function of N , and we find it increases as N 2 for both the Rayleigh and SHG peaks. Again this is consistent with a superradiant mechanism.
Taken together, the results presented here indicate that, compared to a single two-level system, the SHG signal can be enhanced by considering an array of N of two-level systems. Furthermore, the dependence on N of both the emission time and intensity of the field are consistent with a superradiance behavior.

IV. MOTION OF A TWO-LEVEL SYSTEM IN A CAVITY
The two-level systems considered in the previous sections were at fixed positions in space, like e.g. for a given pair of levels in a quantum dot. However, if a two-level system is meant to model a pair of atomic orbitals, then the atomic center-of-mass motion can be an important, if not crucial, element to take into account.
On the experimental side, compelling evidence comes for example from studies of quantum control, where atoms moving across an optical cavity provide information about cavity photons [79], or laser beams across an ion trap provide information about the internal state of the ions [80]. On the theoretical side, the role of atom dynamics has been extensively considered [9,21,29,32,[81][82][83][84][85][86][87][88][89][90], often in terms of a generalized Jaynes-Cummings model where the standard two-level, one-mode Hamiltonian is augmented by a kinetic energy operator (for the center-of-mass motion). Furthermore, the light-matter coupling can become position dependent [9], for example when the atom is moving inside/outside a optical cavity.
The solution of the generalized Jaynes-Cummings model has been approached in many different ways [9,21,29,32,[81][82][83][84][85][86][87][88][89][90]. For example, with or without the RWA, with the center-of-mass motion described classically or quantum mechanically, using density-matrix techniques or resorting to a direct solution of the generalized Rabi equations in wavefunction space. Here, we again consider the exact numerical time-evolution of the full system wavefunction, thus avoiding the RWA. Since we are interested in how the atom dynamics affects fluorescence spectra, we consider a generalized Jaynes-Cummings model with two (pump and fluorescence) modes, and with a center-of-mass that moves longitudinally across a cavity of finite length. Transverse motion is not considered (i.e., space-wise, our system is strictly one-dimensional). We treat the atom dynamics quantum mechanically, but for comparison we also consider a classical description via the Ehrenfest approximation. Our approach includes all these element on equal-footing, and in a single coherent description. This can offer an advantage: for example, the cavity boundaries can have nontrivial effects on the spectra which depend also on the level of description.
Since we intend to look only in a preliminary and explorative way at fluorescence spectra in this setup, we already here anticipate that in our calculations the "atomic" mass value is taken rather small, but not so small that it is necessary to take into account spatial dispersion effects in the radiation-matter interaction. The purpose of this choice is twofold: on the one hand, the atom moves "faster", which alleviates the costs of the numerical time evolution to reach the long-time limit. On the other hand the role of quantum effects in the nuclear motion is enhanced, since on increasing the value of the atomic mass a classical description becomes increasingly appropriate.
It is worth to mention that excitons in solid-state systems (e.g. heterostructures [91]) can also be used for two-level atom optics in quantized light fields, with the exciton dynamics manipulated by optical means. This option has the merit that the value of the exciton electron/hole effective masses (and thus of the total mass) can be tailored by manipulating the band-edge curvatures. Furthermore, using a mass-scaling transformation as described in Appendix A 3, the calculations and results to follow can qualitatively relate to microwave transitions of atoms with realistic masses.
Out of this discussion, the Hamiltonian to consider iŝ wherep andx are the momentum and position operators of the atomic center of mass, M is the atomic mass. As mentioned above, Eq. (27) satisfies a useful scaling property (see Appendix A 3.) For the definition of the other quantities, we refer to Eqs. (1)(2)(3)(4). As before, we assume that the atom is occupied by a single spinless electron, that the cavity field is of frequency ω a and in a coherent state defined byâ|α = α|α . Further, the fluorescent field is of frequency ω b and initially in the vacuum statê b|0 = 0.
The interaction between the light and the atom is given by the couplings g a (x, t) and g b (x, t) respectively, with the spatial dependence of the coupling g a coming from the spatial profile of the cavity mode. We assume that the atomic motion happens only along the cavity axis, and restrict the length of the coordinate axis to the set X = [0, L] of length L. We further divide the coordinate axis into two parts X in and X out , corresponding respectively to inside and outside the cavity, where the set X in = [x 1 , x 2 ] is of length l = x 2 − x 1 and X out = X \ X in is of length L − l. For consistency we need to take 0 < x 1 < x 2 < L. The cavity electric field E is assumed to be in the lowest mode, with a spatial profile given by E(x) = sin((x − x 1 )π/l) for x ∈ X in and E(x) = 0 otherwise. Using the characteristic function χ I , which is unity on the interval I and zero otherwise, the lightmatter couplings are then given by For the coupling to the fluorescent field we have assumed a constant coupling g 1 (g 2 ) inside (outside) the cavity, with a phenomenological decay Γ 1 (Γ 2 ) taking into account collision effects and an effective coupling to additional radiation modes in the continuum.

A. Time dependent fluorescence spectrum
Taking into account the quantum motion of the center of mass significantly increases the numerical effort necessary to obtain the fluorescence spectrum. To simplify the calculations we therefore consider the fluorescence response in the one-photon limit, where the coupling between the fluorescent field and the atom only acts once during the time evolution. For the initial state, it is assumed that i) the atom is prepared with a nuclear wavefunction the form of which is and with the electron in level 1 (with energy 1 ) ii) the cavity field is in a coherent state α, and there are zero photons of the fluorescent field. Thus the system's initial state is denoted by |1, φ, α (the label for fluorescent state being omitted, because of the zero-photon assumption). The spectrum is defined as the probability that at time t there is one photon in the fluorescent field: Here, the state |i, x, n also contains zero photons of the fluorescent field (note theb operator immediately to the right of i, x, n|). In the final state, i denotes the electronic level (1 or 2), x the atomic position, n the number of photons in the cavity field, and the trace over i, x, n. The perturbative limit is obtained by assuming that the time-evolution operator can be written as In Appendix C we show that in this limit the probability is given by where λ are the eigenenergies ofĤ 0 . The coefficients S 1/2 λλ are given by where λ, λ label complete sets of states ofĤ 0 (but again with zero photons in the fluorescent field), and X in and X out have been defined earlier. This expression, used in the next section to calculate the fluorescence spectrum, is valid for all times t and has the practical advantage of requiring only a single diagonalization of the Hamiltonian H 0 . However, in actual calculations the numerical grid for the nuclear coordinate x is confined to the interval (0 ≤ x ≤ L), and, for a given L, the maximum useful t is limited by the need to avoid wavepacket reflection at the interval boundaries.

B. Classical versus quantum nuclear motion
Before studying the full dynamics of the system, we consider the classical limit of the nuclear dynamics as given by the Ehrenfest approximation. Assuming thatx andp in Eq. (27) are replaced by classical variables x and p evolving under the force F (t) = − ∂ xĤt ) , the spatial dependence of g a turns into a time-dependence through g a (t) = g a (x(t)) = g a sin((x(t) − x 1 )π/L). For the coupling to the fluorescent field we assume (temporarily) that there is no decoherence in the cavity, so that with t 0 the time at which the atom has passed through the cavity. The quantum evolution is obtained by solving the time-dependent Schrödinger equation on a grid x n , as described in more detail later on. In the following we work in atomic units (a.u.) and take M = 10 a.u., = 2 − 1 = 0.043 a.u., ω a = , α = 1, g a = 0.1 , g 1 = 0.1 , g 2 = 0.01 and Γ = 0.02 . Except for the coupling g 1 of the fluorescent field to the atom inside the cavity, these values as the same as earlier in the paper with the identification = 0.043 a.u.. However, the value of g 1 was increased to enhance the emission into the fluorescent field, with b †b 1 still applying. A physical notion of the chosen parameters can be gathered by noting that for a cavity of length l c = λ/2 and = ω a /κ (for example, κ = 0.5 for SHG), we have = ( cπ)/(κl c ). Choosing l c = 10 4 a.u., one obtains ≈ 0.043/κ a.u., which can be a reasonable value for example for excitons, and fairly consistent with the value of M = 10 a.u. introduced above [92].
Within the given units, the spatial simulation interval (cavity and outside) of our calculations is L = 10l c (i.e. about 5µm), and the cavity boundaries are set at x 1 = 4l c and x 2 = 5l c . As initial conditions we take p 0 = 0.5 a.u. and x 0 = 4l c for the classical simulations, while in the quantum simulations the initial wave packet is given by the expression in Eq. (30) with p 0 = 0.5 a.u., x 0 = 3.5l c and σ = 3l c . This momentum corresponds to a velocity v ≈ 10 4 cm/s. In Fig. 9 we compare the results obtained with the Ehrenfest approximation with the results of the full quantum evolution. To characterize the atomic motion, we look in the classical case at the functions x(t) and p(t), and in the quantum case at the nuclear probability den- We see in Fig. 9 that in the classical case the atom moves through the cavity with little resistance, and note that the effects oscillations of p is not visible in x due to the scale of the figure. In contrast, the quantum results show a splitting of the atomic wave packet. This is precisely the regime where the Ehrenfest approximation fails [93], since in a classical description the atom must be either reflected or transmitted. However, at t = 6400 a.u. the quantum particle has both a (larger) reflected and (smaller) transmitted contribution, while the classical particle is out of the cavity already at t = 3000 a.u. Thus the atom has a reduced amplitude in the barrier and a reduced coupling to the field. Compared to the classical case, treating the atomic motion at the quantum level also has a large impact on the fluorescence spectrum. This is addressed in the bottom panels for different time snapshots (spectra at different times are rather similar to each other, and only the one at the latest time is fully visible). The Ehrenfest result resembles at great extent that of a stationary atom (cf. Fig. 2), while the spectrum corresponding to the quantum motion looks qualitatively different: It contains two main peaks instead of four, and is asymmetric with respect to the central frequency.
To understand this dissimilarity in behavior, we note that in the classical approach the nuclear wavepacket is perfectly localized both in position and momentum. By contrast, the quantum amplitude gets smaller in the repulsive barrier region. Further, the classical atom sees only a single resonant frequency (Doppler shifted due to the motion) and coupling to the cavity field at each given time of its travel through the cavity. Since the field is strongest in the center of the cavity, where it takes the same value as in the stationary case discussed above (α = 1 and g a = 0.1 ), the main contribution to the fluorescence signal comes from when the atom is in this region. However, compared to the stationary case there is an enhancement of the spectrum for frequencies ω b ≈ , which most likely comes from fluorescent photons emitted in the regions where g a < 0.1 and the fact that the classical atom couples to the photons more strongly.
In contrast, the quantum atom simultaneously experiences a range of resonance frequencies and field-atom couplings, there is a lot of structure in the corresponding spectrum, and its wave function gets low within the barrier region. From the shape of the nuclear wavepacket we expect the dominant contribution to the fluorescence signal to come from when the atom is in the initial and final part of the cavity, the probability distribution being mainly localized to these regions (see Fig. 9). Consequently the spectrum is closer to what could be expected for a stationary atom weakly interacting with a light field (g a < 0.1 ), leading to a smaller splitting between the Mollow peaks. However, a detailed explanation of the asymmetric form of the spectrum is difficult to give, but plausibly related to the varying Doppler shifts associated with the different parts of the atomic wavepacket.

C. Fluorescence and quantum motion
We now consider the fluorescence spectra resulting from a quantum evolution of the coupled atom-photon system. We take L = 10l c , and, as before, the cavity is placed again between x 1 = 4l c and x 2 = 5l c . To solve the Schrödinger equation we consider a grid x n for the atomic position, with 500 points in the interval [0, L]. The fluorescence spectrum is computed from Eq. (33), and to get the atomic probability density N (x, t) = in | i, x, n|Ψ(t) | 2 we solved the Schrödinger equation without the fluorescent field. We have verified that the atomic dynamics is highly insensitive to the presence of the fluorescent field, by explicitly solving the Schrödinger equation with the complete Hamiltonian for a number of values of the fluorescence frequency. This insensitivity is due to the weak atom-field coupling and the absence of the spatial dispersion effects. The weak coupling also guarantees that the first order fluorescence spectrum is a good approximation of the exact one. In the following we let M = 10 a.u. and ω a = 0.043 a.u. be fixed, and take = ω a or = 2ω a for the Mollow or SHG regimes respectively. The remaining parameters are α = 1, g a = 0.1 , g 1 = 0.1 , g 2 = 0.05 , Γ 1 = 0.01 and Γ 2 = 0.02 .
In Fig. 10 we show N (x, t) and P (ω, t) for = ω a . We see that for a stationary atom placed in the center of the cavity, x 0 = 4.5l c and p 0 = 0, the spectrum resembles the Mollow spectrum in Fig. 2. We also find that although the atomic wave packet is initially contained in the cavity, parts of the probability distribution are ejected as time progresses. For higher initial momentum p 0 , we see that an atom initially outside the cavity (at x 0 = 3.5l c ) is either split (for p 0 = 0.5 a.u.) or travels through the cavity (for p 0 = 2 a.u.). This is in agreement with the expectation based on an atom moving in the presence of a dipole force [94], where the force on the particle is proportional to the negative of the detuning and the gradient of the light intensity F ∼ −[ω a −( 2 − 1 )]∂ x I(x). For a field on resonance, an atomic motion in the positive direction leads to a positive detuning via a Doppler shift, so that atom is expelled from regions of higher intensity. This is why a minimal non-zero momentum is needed to pass through the cavity.
In Fig. 11 we show N (x, t) and P (ω, t) for = 2ω a . We note that the results for N (x, t) are rather similar to those in Fig. 10, presumably because the coupling to the radiation is too weak to make a larger difference. As for the Mollow regime above, we find that for an atom at rest with l x 0 = 4.5l c and p 0 = 0, the spectrum resembles the stationary SHG spectrum in Fig. 4. For this initial state the atomic probability distribution is trapped in the cavity, consistent with motion under a dipole force as discussed above. For non-zero initial momentum p 0 , we find that the SHG signal is strongly suppressed, and that the elastic scattering peak is broadened. Following considerations similar to those for Fig. 10, this effect is likely ascribable to the finite extent of the atomic wavepacket.
Interestingly, and as for Fig. 10, on increasing p 0 the intensity of the SHG spectrum exhibits a non monotonic behavior. Further, due to the emission of a fluorescent photon being delayed with respect to the atomic excitation, the resonance frequency of the atom has time to change slightly between the two events. In fact, being a second order process, SHG is expected to be more sensitive to this type of detuning than the resonant scattering. Thus the combination of detuning and the decoherence induced by emission from different parts of the atomic wavepacket is likely the cause of the suppression of the SHG signal.

V. CONCLUSIONS
In optics and in photonics, the two-level system plays the role of a Rosetta-stone for light-matter interactions, at the interface of quantum with classical and linear with nonlinear behavior. In this work, we used this paradigmatic system to address basic aspects of multi-photon fluorescence in the time-dependent and stationary regimes. The fluorescence response was considered in three cases of increasing complexity, namely in a single two-level system interacting with an incident coherent field, as well as in an array of two-level systems, and a finally for a twolevel atom moving across a cavity.
By solving the Schrödinger equation through exact numerical time-propagation, we showed that, depending on the system and field parameters, the time-dependent fluorescence spectrum develops distinct features that in some cases are not captured by a semi-classical treatment of the incident field. Some of these features offer direct evidence that the usual selection rules of perturbative optics, which consider photons and electrons separately, do not apply in the strong coupling regime.
As clear-cut example, we showed that a second harmonic signal (SHG) can occur in a two-level system. This result was analyzed in terms of the parity of the coupled electron-photon states, and we argued that this nonlinear process is allowed even for parity invariant Hamiltonians. We also studied the SHG process in a three-level system, and showed that in an appropriate limit the three-level system reproduces the results of the two-level system.
The SHG signal gets enhanced in a setup with N twolevel systems (compared to the N = 1 case), with trends suggestive of superradiant behavior. This conclusion was gathered by looking analytically at possible signatures of a phase transition in the large N limit, and by analyzing the onset of the SHG response in exact numerical calculations for N 10. However, compared to the case of one atom at rest, both Mollow and SHG signals are suppressed by atomic motion. This effect is especially strong in the quantum case: For quantum atomic motion, the spread in space of the traveling nuclear wavepacket is greatly increased by the presence of the barrier represented by the optical cavity, giving rise to a position dependent Doppler shift. This in turn results in a large frequency dispersion in the fluorescence signal. It is important to specify here that the noted quantum effects were enhanced by choosing artificially small nuclear masses. However, although we focussed here on model atomic systems, it is a fair assumption that many of the results obtained should carry over to solid state two-level systems, where suitable (excitonic) masses could be engineered.
The aforementioned effects are weak, and they manifest at low intensities; this is confirmed by the fact that, for the chosen parameters, a linear and an all-orders treatment of the fluorescence field provide an identical scenario. This also means that our situation does not correspond to standard heterodyne setups, and no intensity renormalization is needed. Even so, said effects should be of some conceptual (if not practical) interest. It can in fact be argued that the found superimposed Mollow-like structure to the SHG signal is specifically distinctive of the genuine "two-level" character of the material system, as also gathered by looking at three-levels results. We add that a similar (albeit weaker) superimposed structure also occurs for higher-order harmonics.
Concerning dissipation effects, we expect on purely speculative grounds that the peculiar four-peak structure of the Mollow spectrum could be dimmed, while the SHG signal would considerably change but still survive for not too strong dissipation.
In conclusion, we have addressed general features of multi-photon fluorescence, but only in very simple model systems: The inclusion of additional radiation modes, more general time-dependent couplings, a careful inclusion of bath effects, and more realistic atomic and solidstate setups are possible directions for future work, and to confirm in a broader sense the robustness of present results. Ultimately, true validation comes from experiment, and we hope that our work will stimulate investigations in that direction.

ACKNOWLEDGMENTS
We wish to thank W. P. Schleich for useful discussions. E. V. B was supported by Crafoordska Stiftelsen. C. V. was supported by the Swedish Research Council.

Appendix A
We provide here some additional motivation and detail about the model and the method of solution.

About dissipative effects
For the first two typologies of systems considered, i.e. one or many two-level system(s) at rest, we assume that (via e.g. a cavity-geometry or an high-optical quality sample in ultra-high vacuum and at helium temperature), the inhomogeneous broadening has been made as negligible as possible. In this way, the focus is solely on the homogeneous broadening. This is due to both radiative and non-radiative components, that in our treatment are accounted for by a total phenomenological Lorentzian damping.
For the the third typology of system (i.e. the "atom" in motion, and where the dissipative environment can be considerably different) we have followed a common practice in the literature, considering only the moving material system and the relevant modes of interest. Ultimately, this choice is also dictated by computational convenience, since the full quantum treatment of nuclear motion adds considerable complexity to the numerics.
Looking ahead, a possible way to include dissipative effects is via Lindblad-type master equations or, alternatively, via non-equilibrium Green's functions (NEGF). NEGF permit to include memory effects and the dispersive contribution of the environment in a very direct and systematic way, and it would be rather interesting to perform a comparison between NEGF and master equation results. These calculations and comparisons are deferred to future work.

Computational details
The short iterated Lanczos method is an efficient algorithm to approximate the time-evolution operator U . This is done by constructing U in a small optimized subspace (the Krylov space), which allows to maintain unitarity of U (in contrast to a straightforward Taylor expansion) while being numerically efficient. We used this algorithm to propagate the many-particle Schrödinger equation, and additional details can be found in Ref. [58].
Regarding the choice of basis, we use two basis states |1 and |2 two describe the "atomic" electron states, and the number states |n a and |n b for the coherent and fluorescent fields respectively. In the last part of the manuscript, where we study the motion through a cavity, we use the position basis |x n on an equidistant grid to describe the "atomic" center of mass.
Finally, the choice of cut-off number(s) for the radiation modes is determined by the convergence (i.e. by increasing the number of states until the results are converged within machine accuracy). For the fluorescent field we found it was sufficient in all cases considered to use a maximum n b ≈ 10. For the coherent field the numerical cut-off depends on α, since the coherent field follows a Poissonian distribution when written in terms of number states. For α = 1 and α = 5 respectively, we found that a cut-off at n a ≈ 30 and n a ≈ 150 is enough.

A scaling property
The Hamiltonian in Eq. (27) satisfies a scaling property relating the full quantum dynamics of systems with different masses. Specifically, we start by considering a scaling parameter Z, and the Schrödinger equation i∂ t ψ(t) =Ĥ(t)ψ(t). Dividing by Z, and setting t = Zt, we get i∂ t φ(t ) = Z −1Ĥ (t /Z)φ(t ), where φ(t ) = ψ(t /Z). By relabeling the time variable, t → t, we then have i∂ t φ(t) =H(t)φ(t), whereH(t) = Z −1Ĥ (t/Z). According to this scaling prescription, a given numerical calculation represents in fact a entire one-parameter set of numerical simulations, where the integration interval, the time dependence inĤ and the fermion-boson interactions are changed.

Appendix B
We consider here the adiabatic elimination approximation (AEA) for the three-level Hamiltonian given by the sum of Eqs. (9),(10) and the free field part. Most often, the AEA is done in connection with the rotating wave approximation (RWA) (see e.g. [64][65][66]95]), and choosing the order in which AEA and RWA are performed can be important [95]. Since our considerations here aim to be qualitative and general in character, we use for sim-plicity the more common protocol where the RWA is introduced before the AEA [65], and before the pump field undergoes a transformation to a coherent photon picture. Proceeding in this way, we obtain At this point the coherent state picture could be introduced. However, already at this stage, the AEA two-level model of Eq. (B2) seems rather different from the original two-level model of Eqs. (1)-(4). Therefore, the results of the main text for the SHG and Mollow spectra in a two-level system should not be ascribed to an adiabatic suppression of the virtual level.

Appendix C
We want to calculate the probability defined in Eq. (31) and repeated here for convenience: (C1) In the perturbative limit, the time-evolution operator becomes: and to further simplify the analysis, we define the probability amplitude SinceĤ 0 is independent of time, the probability amplitude can be found through a straightforward expansion in the eigenstates ofĤ 0 . In the expression for the probability P above, we trace over a complete set of final states |i, x, n , but since any complete set is allowed we can instead choose to trace over the eigenstates ofĤ 0 . In the following we therefore consider the probability amplitude A λ (t, ω), and by inserting a set of complete states we find λ|be −iĤ0(t−t )Ĥ (t )e −iĤ0t |1, φ, α = λ λ λ|b|λ e −i( λ +ω)(t−t ) H λ λ (t )e −i λ t λ |1, φ, α = λ e −i( λ +ω)(t−t )−i λ t H λλ (t ) λ |1, φ, α .
Now the matrix elements H λλ (t ) can be broken into two parts according to Integrating over t we find the probability amplitude to be and inserting this into the expression for the probability we find If necessary, it is possible to further simplify Eq. (C8) by going at long times (i.e. where the exponentials e −Γ k t tend to zero) provided that L is correspondingly taken large enough to avoid atom reflection at the ends of the x-coordinate domain. Arguing that the cross terms vanish for t → ∞, the asymptotic limit becomes This latter result makes contact with the long time limit of the static-atom case of Sec. II. However, to calculate the time-dependent fluorescence spectrum for a moving atom we will go back to the full expression for P (t, ω) in Eq. (C8) here, or Eq. (33) in the main text.