On single-frequency asymptotics for the Maxwell-Bloch equations: mixed states

On single-frequenc y asymptotics for the Maxwell–Bloch equations: mix ed states A.I. K omech 1 and E.A. K opylo v a 1 Institute of Mathematics of BOKU University , V ienna, A ustria alexander .komech@boku.ac.at, elena.k opylov a@boku.ac.at Abstract W e consider damped driven Maxwell–Bloch equations which are finite- dimensional approximation of the damped driv en Maxwell–Schr ¨ odinger equa- tions. The equations describe a single-mode Maxwell field coupled to a two-le vel molecule. Our main result is the construction of solutions with single-frequency asymptotics of the Maxwell field in the case of quasiperi- odic pumping. The asymptotics hold for solutions with harmonic initial val- ues which are stationary states of av eraged equations in the interaction pic- ture. W e calculate all harmonic states and analyse their stability . The calcula- tions rely on the Bloch–Feynman gyroscopic representation of von Neumann equation for the density matrix. The asymptotics follo w by application of the av eraging theory of the Bogolyubov type. The ke y role in the application of the av eraging theory is played by a special a priori estimate. MSC classification : 37J40, 58D19, 37J06, 70H33, 34C25, 34C29, 78A40, 78A60. K e ywor ds : Maxwell–Bloch equations; Bloch–Feynman vector equation; Hamilto- nian structure; density matrix; v on Neumann equation; pumping; averaging theory; single-frequency asymptotics; quantum optics; laser . Contents 1 Introduction 2 2 The Bloch-F eynman “gyroscopic” r epresentation 7 3 Dynamics in the interaction picture 8 4 The a veraging 10 1 Supported by Austrian Science Fund (FWF) P A T 3476224 5 Stationary states f or the av eraged dynamics (harmonic states) 11 6 Spectra of linearised equations at the harmonic states 12 7 Attraction to stable submanif old 14 8 Single-frequency asymptotics 16 8.1 KBM vector field . . . . . . . . . . . . . . . . . . . . . . . . . . 16 8.2 The asymptotics in the interaction picture . . . . . . . . . . . . . 17 A The a priori bounds and well-posedness 18 B The Maxwell–Bloch equations for pur e and mixed states 18 2.1 The Maxwell–Bloch equations for pure states . . . . . . . . . . . 19 2.2 The von Neumann equation for mix ed states . . . . . . . . . . . . 20 C On possible treatment of the laser action 20 1 Intr oduction The Maxwell–Bloch equations (MBE) were introduced by Lamb [ 25 ] for the semi- classical description of the laser action [ 1 , 6 , 17 , 31 , 32 , 33 , 34 ]. The equations are the Galerkin approximation of the Maxwell–Schr ¨ odinger system [ 5 , 16 , 19 , 21 , 22 , 27 , 29 ] Our main goal is construction of solutions with single-frequency asymp- totics. The asymptotics seem to correspond to the laser coherent radiation, which remains a k ey mystery of laser action since its disco very around 1960. The damped dri ven MBE for mix ed states read    ˙ A ( t ) = B ( t ) , ˙ B ( t ) = − Ω 2 A ( t ) − γ B ( t ) + c j ( t ) i ¯ h ˙ ρ ( t ) = [ H ( t ) , ρ ( t )]       , j ( t ) = 2 κ Im ρ 21 ( t ) . (1.1) Here A ( t ) , B ( t ) ∈ R , Ω > 0 is the resonance frequency , γ > 0 is the dissipation coef ficient, c is the speed of light, and ¯ h - the Planck constant. Further, H ( t ) is the Hermitian matrix, and ρ ( t ) = ρ is a nonnegati ve Hermitian 2 × 2-density matrix, H ( t ) =   ¯ h ω 1 ia ( t ) − ia ( t ) ¯ h ω 2   ; ρ =   ρ 11 ρ 21 ρ 21 ρ 22   , tr ρ = 1 , ρ ≥ 0 , (1.2) 2 where ¯ h ω 2 > ¯ h ω 1 are the energy levels of activ e molecules, and the function a ( t ) is gi ven by a ( t ) : = κ c [ A ( t ) + A e ( t )] , κ = p ω , ω = ω 2 − ω 1 > 0 , (1.3) where p ∈ R is proportional to the molecular dipole moment. The pumping A e ( t ) we suppose to be a quasiperiodic function: A e ( t ) = Re [ A e e − i Ω t ] + Re N ∑ 1 A e k e − i Ω k t , where A e , A e k ∈ C , Ω k ∈ R and Ω k  = Ω . (1.4) For solutions to the MBE, the conservation tr ρ ( t ) = const holds since the trace of commutators v anish. W e will consider solutions with const = 1; see ( 1.2 ). W e use the Heaviside–Lorentz units and recall the introduction of the equations in Appendix B . In Appendix A , we prove the a priori bound for solutions X ( t ) = ( A ( t ) , B ( t ) , ρ ( t )) to ( 1.1 ) with | p | / γ = r satisfy the following a priori bounds: | X ( t ) | ≤ D r ( | X ( 0 ) | ) , t ≥ 0 , (1.5) which is proved in Appendix A . The bound implies the well-posedness of the MBE. Our main goal is asymptotics of the Maxwell amplitudes A ( t ) and B ( t ) for solutions X ( t ) to the MBE as | p | , γ → 0. Define complex Maxwell amplitudes by M ( t ) = A ( t ) + iB ( t ) / Ω . Then the first tw o equations of the MBE are equiv alent to ˙ M ( t ) = − i Ω M ( t ) − i γ M 2 ( t ) + 2 ic p ω Im ρ 21 / Ω , where M 2 ( t ) = Im M ( t ) . (1.6) Note that the parameters | p | , γ are very small for many types of lasers, see Appendix C . For p = γ = 0, all solutions to ( 1.6 ) are single-frequency: M ( t ) = e − i Ω t M ( 0 ) . For small p and γ , equation ( 1.6 ) implies that M ( t ) = e − i Ω t M ( 0 ) + Z t 0 R ( s ) d s , where sup s ≥ 0 | R ( s ) | = O ( | p | + γ ) by ( 1.5 ). Hence, M ( t ) = e − i Ω t M ( 0 ) + O ( | p | t ) for p / γ = r with arbitrary r  = 0. In particular, for solutions with any fixed initial states M ( 0 ) , max t ∈ [ 0 , | p | − 1 / 2 ] | M ( t ) − e − i Ω t M ( 0 ) | = O ( | p | 1 / 2 ) , p → 0 , p / γ = r . (1.7) 3 Our main results sho w that in the resonance case, when Ω = ω , the time scale p − 1 / 2 in this asymptotics can be extended to | p | − 1 and e ven more in the case of special harmonic states of the MBE. W e will represent the density matrix via the Pauli matrices as ρ ( t ) = 1 2 [ E + S ( t ) · σ ] , where S ( t ) ∈ B = { S ∈ R 3 : | S | ≤ 1 } , σ = ( σ 1 , σ 2 , σ 3 ) , (1.8) and accordingly , solutions X ( t ) to MBE are represented by ( M ( t ) , S ( t )) ∈ X = C × B . The harmonic states are defined via the interaction pictur e of the MBE and the corresponding avera ged equation with the structure ˙ X ( t ) = p F r ( X ( t )) , t ≥ 0; r = p / γ . (1.9) W e define the harmonic states X ∈ X of the MBE as stationary solutions of ( 1.9 ). W e calculate all harmonic states X r = ( M r , S r ) ∈ X of the MBE with p / γ = r and show that the set is a union Z r = Z r 1 ∪ Z r 2 of two smooth 1D manifolds. W e linearise the dynamics ( 1.9 ) at the harmonic states and calculate the spectra of the linearisations. The calculations show that for Ω = ω and c | r | > | A e | , there exist a nonempty submanifold Z r + ⊂ Z r which is attractiv e under the averaged dynamics ( 1.9 ). Denote the matrix V Ω =   0 Ω 0 − Ω 0 0 0 0 0   . (1.10) Our main results are the following asymptotics for solutions X ( t ) = ( M ( t ) , S ( t )) to equations ( 1.1 ) in the representation ( 1.8 ), with a fix ed quotient p / γ = r . Ev- erywhere belo w we consider the case p > 0 only since the extension to p < 0 is obvious. Theorem 1.1. Let Ω = ω and the pumping be quasiperiodic. Then for any r > 0 , the following asymptotics hold. i) Let the initial state X ( 0 ) = ( M , S ) ∈ Z r . Then the corr esponding solutions admit the following adiabatic asymptotics: max t ∈ [ 0 , p − 1 ] h | M ( t ) − e − i Ω t M + | S ( t ) − e V Ω t S | i = O ( p 1 / 2 ) , p → 0 . (1.11) ii) Let cr > | A e | , and D r d denote a suitable subset of the tub ular d -neighborhood of the stable submanifold Z r + with sufficiently small d > 0 . Then for initial states 4 X ( 0 ) ∈ D r d , the corresponding solutions admit the following asymptotics uniformly in initial values X ( 0 ) ∈ D r d : max t ∈ [ 0 , p − 1 ] h | M ( t ) − e − i Ω t M ∗ | + | S ( t ) − e V Ω t S ∗ ( t ) | i = O [ p 1 / 2 + d ] , p → 0 , (1.12) wher e ( M ∗ , S ∗ ( t )) ∈ Z r + and M ∗ = − A e does not depend on r and on X ( 0 ) ∈ D r d . The Lebesgue measur e | D r d | ∼ d 4 , d → 0 . (1.13) Remark 1.2. W ithout the constraint p / γ = r , the a priori bounds ( 1.5 ) and the asymptotics ( 1.11 )–( 1.12 ) do not hold. For example, the bounds must be different for small and big | p | for a gi ven dissipation γ > 0. The asymptotics do not hold without this restriction since the limiting amplitudes M r depend on r . The asymptotics ( 1.11 ) for S ( t ) reads S ( t ) = e V Ω t S + O ( p 1 / 2 ) , t ∈ [ 0 , p − 1 ] (1.14) which is approximately the precession of S ( t ) about axis S 3 . Let us comment on our approach. The asymptotics ( 1.11 ) mean that the MBE admits solutions M ( t ) = e − i Ω t M r ( t ) , S ( t ) = e V Ω t S r ( t ) with slowly varying en velop- ing amplitudes M r ( t ) ∈ C and S r ( t ) ∈ R 3 for small p and γ with p / γ = r . The amplitudes are solutions to the corresponding dynamical system which is the inter - action picture (or “rotating frame representation”) of the MBE. The slo w v ariation of the amplitudes for small p and γ is equi v alent to the fact that the initial state ( M ( 0 ) , S ( 0 )) = ( M ( 0 ) , S ( 0 )) is a harmonic state, i.e., a stationary solution of ( 1.9 ) with p / γ = r . The calculation of harmonic states and of spectra of the linearisations rely on the Bloch–Feynman representation [ 7 , 14 ] of the von Neumann equation via the vector S ( t ) from ( 1.8 ). W e prov e that the harmonic states ( M , S ) with M  = 0 exist only in the resonance case Ω = ω , and only if A e  = 0. In particular , for the single- frequency pumping A e ( t ) = Re [ A p e − i ω p t ] , the harmonic states with M  = 0 exist only in the case A p  = 0 and triple resonance Ω = ω = ω p . (1.15) The asymptotics ( 1.11 ) follo ws for solutions with every initial state ( M ( 0 ) , S ( 0 )) ∈ Z r by the av eraging theory of Bogolyubov type [ 30 , Theorem 4.3.6]. T o prove asymptotics ( 1.12 ), we combine the theorem with the stability of the branch Z r + . 5 Our calculations show that for all harmonic states ( M , S ) ∈ Z r + , the component M = − A e . This is why the limiting amplitude M ∗ in ( 1.12 ) does not depend on r > 0 and on X ( 0 ) ∈ D r d . The ke y role in the application of the averaging theory is played by the special a priori estimate ( 1.5 ). Remark 1.3. Our approach relies on the a v eraging theory which ne glects oscillat- ing terms. So, it gi ves a justification of the “rotating wav e approximation”, which is widely used in Quantum Optics [ 1 , 17 , 31 , 32 , 33 , 34 ]. The asymptotics ( 1.11 ), ( 1.12 ) specify the time scale and the error of such approximations. Let us comment on related results. The problem of existence of time-periodic solutions to the MBE has been discussed since 1960s. The first results in this di- rection were obtained recently in [ 11 ] and [ 36 ] for various versions of the MBE. In [ 36 ], the N-th order time-periodic solutions were constructed by perturbation techniques. For the phenomenological model [ 2 , 3 ], time-periodic solutions were constructed in [ 11 ] in the absence of time-periodic pumping for small interaction constants. The solutions are obtained as the result of a bifurcation relying on ho- motopy in variance of the degree [ 9 ] and developing the a veraging ar guments [ 10 ]. The period is determined by bifurcation. In [ 20 ], we have established the existence of solutions with T -periodic Maxwell amplitude for any T -periodic pumping without smallness conditions. In [ 24 ], we hav e constructed solutions with asymptotics of type ( 1.11 ), ( 1.12 ) for the MBE with pur e states , using the reduction by the symmetry gauge group. In this case, for any r > 0, the set of stationary states is discrete and consists of one or two points. The construction of the asymptotics for system ( 1.1 ) in the present paper , required completely new technique based on the Bloch–Feynman representation. Up to our knowledge, the single-frequency asymptotics for the MBE with mixed states were not constructed till no w . Let us comment on our exposition. In Section 2 , we construct the representa- tion of the Bloch–Feynman type for the v on Neumann equation from MBE. The dynamics on the interaction picture is calculated in Section 3 , its a veraging in Sec- tion 4 , and all stationary states of the a veraged equations are calculated in Section 5 . In Section 6 we analyse the stability of the stationary states, and in Section 8 we prov e the single-frequency asymptotics ( 1.11 )–( 1.12 ). In Appendix A we estab- lish the bound ( 1.5 ). In Appendix B we comment on the introduction of the MBE with mixed states, and in Appendix C we discuss possible treatment of the laser threshold and laser amplification relying on our results. Acknowledgements. The authors thank S. Kuksin, M.I. Petelin, A. Shnirelman and H. Spohn for longterm fruitful discussions, and the Institute of Mathematics 6 of BOKU Uni versity for the support and hospitality . The research is supported by Austrian Science Fund (FWF) P A T 3476224. 2 The Bloch-F eynman “gyr oscopic” r epresentation Here we represent von Neumann equation from MBE as a “gyroscopic equation” ( 2.5 ) applying the Bloch, Fe ynman & al approach [ 6 , 7 , 14 , 15 ]. Hermitian density matrices ρ with tr ρ = 1 admit the representation ρ = 1 2  1 + S 3 S 1 − iS 2 S 1 + iS 2 1 − S 3  , ( S 1 , S 2 , S 3 ) ∈ R 3 . (2.1) Hence, using the P auli matrices and the nonne gati vity det ρ = 1 − S 2 1 − S 2 2 − S 2 3 ≥ 0, we obtain ( 1.8 ): ρ = 1 2  1 0 0 1  + S 1 2  0 1 1 0  + S 2 2  0 − i i 0  + S 3 2  1 0 0 − 1  = 1 2 [ E + S · σ ] . (2.2) In particular , we can e xpand the “Hamiltonian operator” ( 1.2 ) as H ( t ) = ¯ h 2 ( ω 1 + ω 2 ) E + ¯ h 2 ( ω 1 − ω 2 ) σ 3 − a ( t ) σ 2 = ¯ h 2 ( ω 1 + ω 2 ) E + ¯ h 2 ( ω 1 − ω 2 ) e 3 · σ − a ( t ) e 2 · σ , (2.3) where e 3 = ( 0 , 0 , 1 ) and e 2 = ( 0 , 1 , 0 ) . Now the v on Neumann equation in MBE for ρ ( t ) = 1 2 [ E + S ( t ) · σ ] reads as i ¯ h ˙ S ( t ) · σ = 2 i ¯ h 2 ( ω 1 − ω 2 )[ e 3 ∧ S ( t )] · σ − a ( t )[ e 2 ∧ S ( t )] · σ = 2 i h  ¯ h 2 ( ω 1 − ω 2 ) e 3 − a ( t ) e 2  ∧ S ( t ) i · σ , (2.4) where we have used the formula [ a · σ , b · σ ] = 2 i ( a ∧ b ) · σ for a , b ∈ R 3 . This implies the “gyroscopic equation” of type [ 14 , (4)]: ˙ S ( t ) = θ ( t ) ∧ S ( t ) , θ ( t ) = − ω e 3 − 2 a ( t ) ¯ h e 2 , ω = ω 2 − ω 1 . (2.5) The equation can be also written as ˙ S ( t ) = Θ ( t ) S ( t ) , Θ ( t ) =   0 − θ 3 ( t ) θ 2 ( t ) θ 3 ( t ) 0 − θ 1 ( t ) − θ 2 ( t ) θ 1 ( t ) 0   =   0 ω − 2 a ( t ) / ¯ h − ω 0 0 2 a ( t ) / ¯ h 0 0   ∈ so ( 3 ) , 7 (2.6) which corresponds to the cross-product representation of the Lie algebra so ( 3 ) . In particular , the follo wing conserv ation law holds, | S ( t ) | = const , t ≥ 0 . (2.7) The current j ( t ) = κ S 2 ( t ) , so the MBE ( 1.1 ) reduces to (cf. ( 1.6 ))    ˙ M ( t ) = − i Ω M ( t ) − i γ M 2 ( t ) + ic κ S 2 ( t ) / Ω ˙ S ( t ) = Θ ( t ) S ( t )       , M ( t ) = A ( t ) + iB ( t ) / Ω . (2.8) 3 Dynamics in the interaction pictur e By ( 2.6 ), for small p , γ > 0, the system ( 2.8 ) is a small perturbation of the unper- turbed one, ˙ M ( t ) = − i Ω M ( t ) , ˙ S ( t ) = V ω S ( t ) , (3.1) where V ω is the matrix ( 1.10 ) with Ω = ω which coincides with ( 2.6 ) in the case a ( t ) = 0. Solutions to this system are giv en by M ( t ) = e − i Ω t M , S ( t ) = e V ω t S , where M ∈ C , S ∈ B . (3.2) Our goal is construction of similar solutions to the perturbed system ( 2.8 ), M ( t ) = e − i Ω t M ( t ) , S ( t ) = e V ω t S ( t ) (3.3) with slo wly varying en veloping amplitudes: for a wide interval of time [ 0 , T ( p )] sup t ∈ [ 0 , T ( p )]  | M ( t ) − M ( 0 ) | + | S ( t ) − S ( 0 ) |  → 0 , p → 0 , p / γ = r  = 0 . (3.4) Substituting ( 3.3 ) into ( 2.8 ), we obtain the dynamical equations for the en veloping amplitudes (the “interaction picture”):    ˙ M ( t ) = − ie i Ω t  γ Im ( e − i Ω t M ( t )) − ˜ κ ( e V ω t S ( t )) 2  ˙ S ( t ) = e − V ω t ( Θ ( t ) − V ω ) e V ω t S ( t )       , ˜ κ = c κ Ω . (3.5) 8 The equations are called as the interaction picture of ( 2.8 ). By ( 2.5 ), ( 2.6 ), and ( 1.3 ), Θ ( t ) − V ω = 2 a ( t ) ¯ h   0 0 − 1 0 0 0 1 0 0   , a ( t ) : = κ c [ A ( t ) + A e ( t )] , κ = p ω . (3.6) Note that e V ω t is he dynamical group of the gyroscopic equation ( 2.5 ) with the angular velocity θ ( t ) = − ω e 3 . Hence, e V ω t is the rotation about e 3 with the angular velocity ω : e V ω t =   cos ω t sin ω t 0 − sin ω t cos ω t 0 0 0 1   (3.7) Substituting ( 3.6 ) and ( 3.7 ) into ( 3.5 ), we obtain    ˙ M ( t ) = − pie i Ω t  γ 1 ( M 2 ( t ) cos Ω t − M 1 ( t ) sin Ω t ) − κ 1 ( − S 1 ( t ) sin ω t + S 2 ( t ) cos ω t )  , ˙ S ( t ) = − pb [ A ( t ) + A e ( t )]  S 3 ( t )( e 1 cos ω t + e 2 sin ω t ) − e 3 ( S 1 ( t ) cos ω t + S 2 ( t ) sin ω t )  , where γ 1 = γ p = 1 r , κ 1 = ˜ κ p = c ω Ω , b = 2 ω c ¯ h . (3.8) The equations can be written as ˙ M ( t ) = p f r ( M ( t ) , S ( t ) , t ) , ˙ S ( t ) = p g r ( M ( t ) , S ( t ) , t ) , (3.9) where the functions f and g are giv en by    f r ( M , S , t )=  − i cos Ω t + sin Ω t  γ 1 ( M 2 cos Ω t − M 1 sin Ω t ) − κ 1 ( − S 1 sin ω t + S 2 cos ω t )  g r ( M , S , t )= − b  M 1 cos Ω t + M 2 sin Ω t + A e ( t )  S 3 ( e 1 cos ω t + e 2 sin ω t ) − e 3 ( S 1 cos ω t + S 2 sin ω t )        . (3.10) Remark 3.1. It is important that the coef ficients γ 1 , κ 1 , b depend only on Ω , ω , r . Hence, for any fixed Ω , ω > 0 and r  = 0, the asymptotics of solutions to sys- tems( 3.9 ) as p → 0 and p / γ = r can be calculated by methods of the av eraging theory [ 8 , 30 ]. 9 4 The a veraging The av eraged equations ( 3.9 ) read ˙ M ( t ) = p f r ( M ( t ) , S ( t )) , ˙ S ( t ) = p g r ( M ( t ) , S ( t )) , (4.1) where f r ( M , S ) = ⟨ f r ( M , S , · ) ⟩ = lim T → ∞ 1 T Z T 0 f r ( M , S , t ) d t , g r ( M , S ) = ⟨ g r ( M , S , · ) ⟩ . (4.2) Let us calculate the av erages ( 4.2 ). The results differ drastically for the resonance case Ω = ω and non-resonance Ω  = ω . In notation ( 1.4 ), A e = A e 1 + i A e 2 , where A e = 2 ⟨ A e ( t ) e i Ω t ⟩ , A e 1 = 2 ⟨ A e ( t ) cos Ω t ⟩ , A e 2 = 2 ⟨ A e ( t ) sin Ω t ⟩ . (4.3) Resonance case Ω = ω . Using the expressions ( 3.10 ), we obtain the av eraged vector field f r ( M , S ) = − i 2 [ γ 1 M 2 − κ 1 S 2  − 1 2  γ 1 M 1 − κ 1 S 1 ] = − γ 1 2 [ M 1 + i M 2 ] + κ 1 2 [ S 1 + i S 2 ] . (4.4) g r ( M , S ) = − b 2  M 1 [ − e 3 S 1 + S 3 e 1 ] + M 2 [ − e 3 S 2 + S 3 e 2 ] − e 3 ( S 1 A e 1 + S 2 A e 2 ) + S 3 ( e 1 A e 1 + e 2 A e 2 )  = − b 2  S 3 [ e 1 ( M 1 + A e 1 ) + e 2 ( M 2 + A e 2 )] − e 3 [ S 1 ( M 1 + A e 1 ) + S 2 ( M 2 + A e 2 )]  . (4.5) Non-resonance case Ω  = ω . In this case, the calculations ( 4.4 ) simplify to f r ( M , S ) = − γ 1 2 M , (4.6) which implies nonexistence of stationary states with nonzero Maxwell field for the system ( 4.1 ), and also the decay M ( t ) = M ( 0 ) e − γ 2 t . In contrast, the resonance equation ( 4.4 ) includes the interaction term which can prev ent the decay of the Maxwell amplitude, that is expected physically . The a veraged system. By ( 4.4 ) and ( 4.5 ), the averaged equations ( 4.1 ) (or ( 1.9 )) in the resonance case read as    ˙ M ( t ) = − p 2  γ 1 [ M 1 + i M 2 ] − κ 1 [ S 1 + i S 2 ]  ˙ S ( t ) = − pb 2  S 3 [ e 1 ( M 1 + A e 1 ) + e 2 ( M 2 + A e 2 )] − e 3 [ S 1 ( M 1 + A e 1 ) + S 2 ( M 2 + A e 2 )]        . (4.7) 10 5 Stationary states f or the a veraged dynamics (harmonic states) In this section, we calculate all harmonic states ( M , S ) , i.e., stationary states for the av eraged equations ( 4.7 ), in the resonance case Ω = ω . The states satisfy    0 = γ 1 [ M 1 + i M 2 ] − κ 1 [ S 1 + i S 2 ] 0 = S 3 [ e 1 ( M 1 + A e 1 ) + e 2 ( M 2 + A e 2 )] − e 3 [ S 1 ( M 1 + A e 1 ) + S 2 ( M 2 + A e 2 )]       . (5.1) It is important that the stationary equations depend on r but do not depend on p . The first equation of ( 5.1 ) gi ves S 1 = α r M 1 , S 2 = α r M 2 , where α r = γ 1 κ 1 = 1 cr > 0 . (5.2) The second equation of ( 5.1 ) is equi valent to the system    0 = S 3 ( M 1 + A e 1 ) 0 = S 3 ( M 2 + A e 2 ) 0 = S 1 ( M 1 + A e 1 ) + S 2 ( M 2 + A e 2 )       . (5.3) The last equation together with ( 5.2 ) gi ve M 1 ( M 1 + A e 1 ) + M 2 ( M 2 + A e 2 ) = 0 which defines the circle S : | M | 2 + M · A e = 0, or equiv alently , | M + A e / 2 | 2 = | A e | 2 / 4, so M = M ( θ ) = − A e 2 + | A e | 2 e i θ , θ ∈ [ 0 , 2 π ] . (5.4) Recall that | S | ≤ 1 by ( 1.8 ). Hence, the set of all stationary states with a fixed r = p / γ is the union of two real 1D manifolds Z r = Z r 1 ∪ Z r 2 , where Z r 1 = { ( M , S ) ∈ C × B : M ∈ S , S 1 = α r M 1 , S 2 = α r M 2 , S 3 = 0 } . (5.5) The second summand Z r 2 is Z r 2 = { ( M , S ) ∈ C × B : M = − A e , S 1 = α r M 1 , S 2 = α r M 2 } . (5.6) Remark 5.1. i) Z r 1  = / 0 for all r  = 0 since M ( 0 ) = 0. ii) Z r 2  = / 0 if f | α r A e | ≤ 1, i.e. cr ≥ | A e | . iii) The intersection of the manifolds is at most one point: Z r 1 ∩ Z r 2 =  ( − A e , − α r A e , 0 ) , α r | A e | ≤ 1 / 0 , α r | A e | > 1     . (5.7) 11 Remark 5.2. Let ( M ( t ) , S ( t )) be a solution to the av eraged equations ( 4.7 ). Then ( e i θ M ( t ) , e V 1 θ S ( t )) with θ ∈ R is the solution to the same equations with A e re- placed by e i θ A e and rotated vectors e 1 and e 2 . The same correspondence holds for solutions to stationary equations ( 5.1 ) that is particularly obvious for the states ( 5.6 ). 6 Spectra of linearised equations at the harmonic states Here, we calculate i) linearisations of the av eraged system ( 4.7 ) at the harmonic states and ii) their spectra. By ( 4.4 ) and ( 4.5 ), we obtain, omitting the index r :            f i ( M , S ) = − 1 2  γ 1 M i − κ 1 S i  , i = 1 , 2 , g k ( M , S ) = − b 2  S 3 ( M k + A e k )  , k = 1 , 2 , g 3 ( M , S ) = b 2  S 1 ( M 1 + A e 1 ) + S 2 ( M 2 + A e 2 )] Dif ferentiating, we get for i , j , k , l = 1 , 2                        ∂ f i ∂ M j = − γ 1 2 δ i j , ∂ f i ∂ S l = κ 1 2 δ il , ∂ f i ∂ S 3 = 0 ∂ g k ∂ M j = − b 2 S 3 δ k j , ∂ g k ∂ S l = 0 , ∂ g k ∂ S 3 = − b 2 ( M k + A e k ) ∂ g 3 ∂ M j = b 2 S j , ∂ g 3 ∂ S l = b 2 ( M l + A e l ) , ∂ g 3 ∂ S 3 = 0 Thus, we obtain the Jacobian J = J ( M , S ) = 1 2       − γ 1 0 κ 1 0 0 0 − γ 1 0 κ 1 0 − b S 3 0 0 0 − b ( M 1 + A e 1 ) 0 − b S 3 0 0 − b ( M 2 + A e 2 ) b S 1 b S 2 b ( M 1 + A e 1 ) b ( M 2 + A e 2 ) 0       . Let us consider the cases ( M , S ) ∈ Z r 1 and ( M , S ) ∈ Z r 2 separately . I. For ( M , S ) ∈ Z r 1 , the Jacobian simplifies to J = 1 2       − γ 1 0 κ 1 0 0 0 − γ 1 0 κ 1 0 0 0 0 0 − b B 1 0 0 0 0 − b B 2 b S 1 b S 2 b B 1 b B 2 0       , 12 where we denote B k = M k + A e k . Hence, det ( 2 J − λ ) =       − γ 1 − λ 0 κ 1 0 0 0 − γ 1 − λ 0 κ 1 0 0 0 − λ 0 − b B 1 0 0 0 − λ − b B 2 b S 1 b S 2 b B 1 b B 2 − λ       = − ( γ 1 + λ )     − γ 1 − λ 0 κ 1 0 0 − λ 0 − b B 1 0 0 − λ − b B 2 b S 2 b B 1 b B 2 − λ     + κ 1     0 − γ 1 − λ κ 1 0 0 0 0 − b B 1 0 0 − λ − b B 2 b S 1 b S 2 b B 2 − λ     = ( γ 1 + λ ) 2   − λ 0 − b B 1 0 − λ − b B 2 b B 1 b B 2 − λ   − κ 1 ( γ 1 + λ )   0 − λ − b B 1 0 0 − b B 2 b S 2 b B 1 − λ   + κ 1 b S 1   γ 1 + λ κ 1 0 0 0 − b B 1 0 − λ − b B 2   = ( γ 1 + λ ) 2  − λ 3 − λ b 2 B 2 1 − λ b 2 B 2 2  − κ 1 ( γ 1 + λ ) λ b 2 S 2 B 2 − κ 1 ( γ 1 + λ ) λ b 2 S 1 B 1 = − ( γ 1 + λ ) λ  ( γ 1 + λ )  λ 2 + b 2 ( B 2 1 + B 2 2 )  + κ 1 b 2 ( S 1 B 1 + S 2 B 2 )  = − ( γ 1 + λ ) 2 λ  λ 2 + b 2 ( B 2 1 + B 2 2 )  since S 1 B 1 + S 2 B 2 = 0 by ( 5.3 ). Hence, λ 1 = λ 2 = − γ 1 , λ 3 = 0 , λ 4 , 5 = ± ib q B 2 1 + B 2 2 . II. For ( M , S ) ∈ Z r 2 , the Jacobian reads J = 1 2       − γ 1 0 κ 1 0 0 0 − γ 1 0 κ 1 0 − b S 3 0 0 0 0 0 − b S 3 0 0 0 b S 1 b S 2 0 0 0       . (6.1) Hence, det ( 2 J − λ ) =       − γ 1 − λ 0 κ 1 0 0 0 − γ 1 − λ 0 κ 1 0 − b S 3 0 − λ 0 0 0 − b S 3 0 − λ 0 b S 1 b S 2 0 0 − λ       = − λ     − γ 1 − λ 0 κ 1 0 0 − γ 1 − λ 0 κ 1 − b S 3 0 − λ 0 0 − b S 3 0 − λ     = − λ  ( γ 1 + λ )   γ 1 + λ 0 κ 1 0 − λ 0 b S 3 0 − λ   + b S 3   0 κ 1 0 γ 1 + λ 0 κ 1 b S 3 0 − λ    13 = − λ  ( γ 1 + λ )( λ 2 ( γ 1 + λ ) + λ κ 1 b S 3 ) + b S 3 ( κ 2 1 b S 3 + λ κ 1 ( γ 1 + λ )  = − λ  λ 2 ( γ 1 + λ ) 2 + 2 λ ( γ 1 + λ ) b κ 1 S 3 + b 2 κ 2 1 S 2 3  = − λ  λ 2 + λ γ 1 + b κ 1 S 3  2 . Hence, we hav e the follo wing roots: λ 1 , 2 ( S 3 ) = − γ 1 2 + q γ 2 1 − 4 b κ 1 S 3 2 , λ 3 , 4 ( S 3 ) = − γ 1 2 − q γ 2 1 − 4 b κ 1 S 3 2 , λ 5 = 0 . (6.2) By ( 3.8 ), the spectrum of the linearised system ( 4.1 ) at the stationary states ( M , S ) ∈ Z r 2 consists of p λ 1 , 2 ( S 3 ) = − γ 2 + p γ 2 − 4 p 2 b κ 1 S 3 2 , p λ 3 , 4 ( S 3 ) = − γ 2 − p γ 2 − 4 p 2 b κ 1 S 3 2 , λ 5 = 0 . By ( 3.8 ), we hav e b κ 1 = 2 ω ¯ h > 0. Hence, for p  = 0, Re λ k ( S 3 ) < 0 for k = 1 , . . . , 4 iff S 3 > 0 . (6.3) 7 Attraction to stable submanif old Let us assume that r is sufficiently lar ge, so that cr > | A e | . (7.1) Then Z r 2  = / 0 and | α r A e | = | A e | / ( cr ) < 1, and we can rewrite ( 5.6 ) as Z r 2 = { Z r 2 ( S 3 ) = ( − A e , − α r A e 1 , − α r A e 2 , S 3 ) : S 3 ∈ [ − β r , β r ] } , β r : = q 1 − α 2 r | A e | 2 > 0 . (7.2) For 0 ≤ a ≤ b ≤ 1 and d > 0 denote Z r + ( a , b ) = { Z r 2 ( S 3 ) : S 3 ∈ [ a , b ] } , U r d ( a , b ) = { Z r 2 ( S 3 ) + N : N ⊥ T Z r 2 ( S 3 ) , S 3 ∈ [ a , b ] , | N | ≤ d } . (7.3) By ( 6.3 ), for s > 0 and p  = 0, we have max 1 ≤ k ≤ 4 S 3 ∈ [ s , β r ] Re λ k ( S 3 ) < − ν ( s ) < 0 . (7.4) The tubular domains U r d ( a , b ) with 0 < a ≤ b < 1 are attracted to Z r + in the follo wing sense. 14 Lemma 7.1. Let ( 7.1 ) hold, s ∈ ( 0 , β r / 4 ) and ( M ( t ) , S ( t )) be the solution to the avera ged equations ( 4.1 ) with the initial state ( M ( 0 ) , S ( 0 )) ∈ D r d : = U r d ( 2 s , β r − 2 s ) . Then for sufficiently small d , µ > 0 , max t ∈ [ 0 , p − 1 ] dist (( M ( t ) , S ( t )) , Z r + ) = O ( d ) . (7.5) Pr oof. The Jacobian ( 6.1 ) admits a basis of eigen vectors and generalised eigenv ec- tors v k ( S 3 ) with S 3 ∈ ( 0 , 1 ] and k = 1 , . . . , 5 corresponding to eigen v alues ( 6.2 ). The set Z r + is one-dimensional manifold of stationary states of the vector field ( f r , g r ) . Hence, the tangent vectors to Z r + at any point belong to the kernel Ker J . So, the tangent vectors correspond to the eigenv alue λ 5 = 0. Therefore, the eigen- vectors corresponding to the stable eigen v alues ( 6.3 ), are transversal to Z r + . The vectors can be chosen piece-wise continuous in S 3 ∈ ( 0 , β r ] since the multiplicity of the eigen v alues ( 6.2 ) is constant. Define the map R 4 × ( 0 , β r ] → R 5 by ( x 1 , x 2 , x 3 , x 4 , x 5 ) 7→ ( M , S ) = Z r 2 ( x 5 ) + 4 ∑ 1 x k v k ( x 5 ) , x 5 ∈ ( 0 , β r ] . (7.6) Let us check that the map is nondegenerate at e very point ( 0 , 0 , 0 , 0 , x 5 ) with x 5 ∈ ( 0 , β r ] which corresponds to ( M , S ) = Z r 2 ( x 5 ) ∈ Z r + . Indeed, by ( 7.2 ), the last column of the Jacobian ∂ ( M , S ) / ∂ x is the vector ∂ ( M , S ) / ∂ x 5 = ( 0 , 0 , 0 , 0 , 1 ) = v 5 ( x 5 ) . On the other hand for k ≤ 4, we hav e ∂ ( M , S ) / ∂ x k = v k ( x 5 ) . Hence, the columns of the Jacobian are linearly independent. Hence, the coordinates x ( M , S ) can be choosen piece-wise continuous in a neighborhood of Z r 2 ( s , β r − s ) . By ( 7.6 ) and ( 5.6 ), we hav e x 5 = S 3 for x 1 = x 2 = x 3 = x 4 = 0. Let us fix an s ∈ ( 0 , β r / 4 ) and suf ficiently small d > 0 such that the coordi- nates x are well defined in U r d ( s , β r − s ) . Then the distance d = d (( M , S ) , Z r 2 ) for ( M , S ) ∈ U r d ( s , β r − s ) is equi v alent to d = d ( x ) = q ∑ 4 1 x 2 k . For δ > 0 let us denote T r δ ( s ) = { x ∈ R 5 : x 5 ∈ [ s , β r − s ] , d ≤ δ } . Then the coordinates x are well defined in T r δ ( s ) , and the equations ( 4.1 ) become    ˙ x 1 = p [ λ 1 ( x ) x 1 + q 1 ( x )] , ˙ x 2 = p [ λ 2 ( x ) x 2 + ε 1 ( x ) x 1 + q 2 ( x )] ˙ x 3 = p [ λ 3 ( x ) x 3 + q 3 ( x )] , ˙ x 4 = p [ λ 4 ( x ) x 4 + ε 2 ( x ) x 3 + q 4 ( x )] , ˙ x 5 = pq 5 ( x ) ,       , (7.7) where q k ( x ) = O ( d 2 ) and | ε k ( x ) | with k = 1 , 2 are suf ficiently small. By ( 7.4 ), the equations ( 7.7 ) imply that for solutions x ( t ) ∈ T r δ ( s ) and sufficiently small δ > 0, ∂ t d 2 ( x ( t )) ≤ 0 . (7.8) 15 Therefore, for any solution x ( t ) with initial state x ( 0 ) ∈ T r δ ( s ) , d 2 ( x ( t )) ≤ δ 2 until x 5 ( t ) ∈ [ s , β r − s ] . (7.9) This inequality allows us to estimate the exit time t ∗ from T r δ ( s ) of solutions with initial states x ( 0 ) ∈ T r δ ( 2 s ) . Namely , in the case t ∗ < ∞ and sufficiently small δ > 0, the coordinate x 5 ( t ) must pass either the segment [ s , 2 s ] or [ β r − 2 s , β r − s ] . In both cases the last equation of ( 7.7 ) together with ( 7.9 ) imply that s / 2 ≤ Z t ∗ 0 | ˙ x 5 ( t ) | d t ≤ p C Z t ∗ 0 d 2 ( x ( t )) d t ≤ p C δ 2 Z t ∗ 0 d t = p C δ 2 t ∗ . (7.10) Thus, t ∗ > p − 1 s 2 C δ 2 . Finally , for sufficiently small δ > 0 we can choose s = 2 C δ 2 . Then ( 7.9 ) implies ( 7.5 ) for suf ficiently small d > 0. Remark 7.2. By ( 7.4 ), equations ( 7.7 ) with k ≤ 4 imply an exponential approach to the manifold Z r + . Hence, the last equation with k = 5, gives the exponential decay of | ˙ x 5 | . This is why the trajectory remains for a long time in a small neighborhood of Z r + . 8 Single-fr equency asymptotics In this section, we prov e the asymptotics ( 1.11 ), ( 1.12 ) in the resonance case Ω = ω . By the relations ( 3.3 ), the asymptotics are equiv alent to the corresponding asymptotics of solutions ( M ( t ) , S ( t )) to the interaction dynamics ( 3.9 ) which we prov e applying the results of the av eraging theory [ 30 ]. T o justify the application, we are going to check suitable properties of the system ( 3.9 ). 8.1 KBM vector field Let us denote v r = ( f r , g r ) the vector field ( 3.10 ) of the system ( 3.9 ). It is easy to check that in the case of almost periodic pumping A e ( t ) and each r  = 0, for any bounded region D ⊂ X , we hav e sup ( M , S ) ∈ D 1 T    Z T 0 [ v r ( M , S , t ) − v r ( M , S )] d t    → 0 , T → ∞ . (8.1) Moreov er , v r is the Lipschitz continuous vector field in any bounded region D ⊂ X . Hence, v r is a KBM (Krylov–Bogolyubo v–Mitropolsky) vector field in any bounded region D according to [ 30 , Definition 4.2.4]. 16 Furthermore, for the quasiperiodic pumping ( 1.4 ), the formulas ( 3.10 ) implies that δ D ( p ) : = p sup ( M , S ) ∈ D sup T ∈ [ 0 , p − 1 ]    Z T 0 [ v r ( M , S , t ) − v r ( M , S )] d t    = O ( p ) , p → 0 , (8.2) where δ D ( p ) is the corresponding or der function defined in [ 30 , Lemma 4.6.4]. Remark 8.1. For more general almost periodic pumping A e ( t ) , the order function can be dif ferent [ 30 , Section 4.6]. 8.2 The asymptotics in the interaction picture Here we prov e Theorem 1.1 . i) The asymptotics ( 1.11 ) is equi valent to max t ∈ [ 0 , p − 1 ] h | M ( t ) − M | + | S ( t ) − S | i = O ( p 1 / 2 ) , p → 0 , p / γ = r , (8.3) The initial state ( M ( 0 ) , S ( 0 )) = ( M ( 0 ) , S ( 0 )) by ( 3.3 ). Hence, ( M ( 0 ) , S ( 0 )) = ( M , S ) . Ho we ver , ( M , S ) is a stationary state for the av eraged system ( 1.9 ), or equi valently , to ( 4.7 ). Therefore, Theorem 4.3.6 of [ 30 ] implies ( 8.3 ). Indeed, both conditions 1 and 2 of the theorem hold in our case since a) v r is a KBM-vector field in any bounded region D ⊂ X with the order function ( 8.2 ); b) solutions ( M ( t ) , S ( t )) of the system ( 3.5 ) with p / γ = r and small p > 0 are uniformly bounded by the a priori estimate ( A.4 ). Hence, ( 8.3 ) is prov ed. ii) Similarly , the asymptotics ( 1.12 ) is equi valent to max t ∈ [ 0 , p − 1 ] h | M ( t ) − M ∗ | + | S ( t ) − S ∗ ( t ) | i = O [ p 1 / 2 + d ] , p → 0 , p / γ = r . (8.4) Define the domain D r d as in Lemma 7.1 and denote by ( M ( t ) , S ( t )) the solution to the av eraged system ( 4.7 ) with the initial value ( M ( 0 ) , S ( 0 )) = X ( 0 ) ∈ D r d . Then the same Theorem 4.3.6 of [ 30 ] implies that max t ∈ [ 0 , p − 1 ] | ( M ( t ) , S ( t )) − ( M ( t ) , S ( t )) | = O ( p 1 / 2 ) , p → 0 , p / γ = r . (8.5) 17 On the other hand, Lemma 7.1 implies that there exist ( M ∗ ( t ) , S ∗ ( t )) ∈ Z r + such that max t ∈ [ 0 , p − 1 ] h | M ( t ) − M ∗ ( t ) | + | S ( t ) − S ∗ ( t ) | i = O ( d ) . (8.6) No w ( 8.4 ) follows since M ∗ ( t ) = − A e by ( 5.6 ). A The a priori bounds and well-posedness In this section, we pro ve the a priori bounds ( 1.5 ) assuming that A e ( t ) ∈ C [ 0 , ∞ ) . The bounds imply the well-posedness of the MBE in the phase space X = C × R 3 . The density matrices ( 1.8 ) all are bounded by the conserv ation ( 2.7 ) . Hence, it remains to prov e the a priori estimates for the Maxwell amplitudes ( A ( t ) , B ( t )) . The follo wing lemma is prov ed in [ 24 ]. Lemma A.1. Ther e e xists a L yapunov function V ( A , B ) such that a 1 [ A 2 + B 2 ] ≤ V ( A , B ) ≤ a 2 [ A 2 + B 2 ] where a 1 , a 2 > 0 , (A.1) and for solutions to ( 1.1 ), the function V ( t ) = V ( A ( t ) , B ( t )) satisfies the inequality ˙ V ( t ) ≤ − γ b V ( t ) + d p 2 γ , t > 0; b , d > 0 . (A.2) Corollary A.2. Solving the inequality ( A.2 ), we obtain: V ( t ) ≤ V ( 0 ) + d b r 2 , t ≥ 0 , r = p / γ . (A.3) Hence, for solutions to ( 1.1 ) with p / γ = r , the following bounds hold: A 2 ( t ) + B 2 ( t ) ≤ D r ( A 2 ( 0 ) + B 2 ( 0 )) , t ≥ 0 . (A.4) Now ( 1.5 ) follows. B The Maxwell–Bloch equations for pure and mixed states There are various versions of the Maxwell–Bloch equations, see for instance [ 25 ] and [ 2 , 3 ]. In this section we recall introduction of the MBE for pure states and also for mixed states described by density matrix. 18 2.1 The Maxwell–Bloch equations for pur e states In [ 20 ], the MBE for pure states were obtained as the Galerkin approximation of the damped dri ven Maxwell–Schr ¨ odinger system. The approximation consists of a single-mode Maxwell field coupled to two-le vel molecule in a bounded ca vity V ⊂ R 3 : A ( x , t ) = A ( t ) X ( x ) , ψ ( x , t ) = C 1 ( t ) ψ 1 ( x ) + C 2 ( t ) ψ 2 ( x ) , x ∈ V . (B.1) Here A ( x , t ) denotes the vector potential of the Maxwell field, and X ( x ) is a nor- malised eigenfunction of the Laplace operator in V under suitable boundary value conditions with an eigen value − Ω 2 / c 2 . By ψ l we denote some normalised eigen- functions of the Schr ¨ odinger operator H : = − ¯ h 2 2m ∆ + e Φ ( x ) with the corresponding eigen v alues ¯ h ω 1 < ¯ h ω 2 , where Φ ( x ) is the molecular (ion’ s) potential. The MBE read as the Hamiltonian system with a dissipation and an external source: 1 c 2 ˙ A ( t ) = ∂ B H , 1 c 2 ˙ B ( t ) = − ∂ A H − γ c 2 B ; i ¯ h ˙ C l ( t ) = ∂ C l H , l = 1 , 2 . (B.2) Here the Hamiltonian is defined as H ( A , B , C 1 , C 2 , t ) = H ( A X , B X , C 1 ϕ 1 + C 2 ϕ 2 , t ) , where H is the Hamiltonian of the coupled Maxwell–Schr ¨ odinger equations with pumping. Neglecting the spin and scalar potential (which can be easily added), the Hamiltonian H , in the traditional dipole appr oximation , reads as [ 20 , (A.5)]: H ( A , B , C , t ) = 1 2 c 2 [ B 2 + Ω 2 A 2 ] + ¯ h ω 1 | C 1 | 2 + ¯ h ω 2 | C 2 | 2 − 2 κ c [ A + A e ( t )] Im [ C 1 C 2 ] , C = ( C 1 , C 2 ) . (B.3) No w the Hamilton equations ( B.2 ) become    ˙ A ( t ) = B ( t ) , ˙ B ( t ) = − Ω 2 A ( t ) − γ B ( t ) + c j ( t ) i ¯ h ˙ C 1 ( t ) = ¯ h ω 1 C 1 ( t ) + ia ( t ) C 2 ( t ) , i ¯ h ˙ C 2 ( t ) = ¯ h ω 2 C 2 ( t ) − ia ( t ) C 1 ( t )       , j ( t ) = 2 κ Im [ C 1 ( t ) C 2 ( t )] . (B.4) The charge conservation | C 1 ( t ) | 2 + | C 2 ( t ) | 2 = const follo ws by differentiation from the last two equations of ( B.4 ) since the function a ( t ) is real-valued. W e consider solutions with const = 1, cf. ( 1.2 ): | C 1 ( t ) | 2 + | C 2 ( t ) | 2 = 1 , t > 0 . (B.5) 19 2.2 The von Neumann equation f or mixed states The Schr ¨ odinger amplitudes C 1 ( t ) and C 2 ( t ) in ( 1.1 ) can be replaced by density matrix ρ which is a nonne gati ve Hermitian 2 × 2-matrix ( 1.2 ). In particular , for the wa ve function ( B.1 ), the corresponding density matrix reads as ρ ( t ) = | C ( t ) ⟩⟨ C ( t ) | =   | C 1 ( t ) | 2 C 1 ( t ) C 2 ( t ) C 2 ( t ) C 1 ( t ) | C 2 ( t ) | 2   , C ( t ) =  C 1 ( t ) C 2 ( t )  . (B.6) The matrix satisfies the trace condition in ( 1.2 ) by the charge conservation ( B.5 ). The last line of ( B.4 ) in the vector form reads i ¯ h ˙ C ( t ) = H ( t ) C ( t ) , H ( t ) =   ¯ h ω 1 ia ( t ) − ia ( t ) ¯ h ω 2   . (B.7) Accordingly , the density matrix ( B.6 ) satisfies von Neumann equation [ 32 ]: i ¯ h ˙ ρ ( t ) = i ¯ h [ ˙ C ( t ) ⊗ C ( t ) + C ( t ) ⊗ ˙ C ( t )] = [ H ( t ) C ( t )] ⊗ C ( t ) − C ( t ) ⊗ [ H ( t ) C ( t )] = [ H ( t ) , ρ ( t )] . (B.8) No w the Maxwell–Bloch system ( B.4 ) is replaced by ( 1.1 ). General density matrix describes an ensemble of molecules with the pure states C ( k , t ) = ( C 1 ( k , t ) , C 2 ( k , t )) satisfying ( B.5 ): ρ ( t ) = ∑ k p k C ( k , t ) ⊗ C ( k , t ) , p k ≥ 0 , ∑ k p k = 1 , (B.9) where p k are the probabilities of the pure states. C On possible tr eatment of the laser action Here we discuss possible treatment of the laser action relying on the obtained re- sults. On the smallness of the parameters. Note that the dipole moment p and dis- sipation coefficient γ are very small for many types of lasers. In particular , the dissipation coefficient for the Ruby laser is the electrical conduction of corundum which is γ ∼ 10 − 14 in the Heaviside–Lorentz units [ 13 , 33 , 34 ]. For the dipole mo- ment typically | p | ∼ 10 − 18 according to [ 26 ] that agrees with the classical dipole moment ed / 2, where d ∼ 10 − 8 cm is the molecular diameter , and e ∼ 10 − 10 is the elementary charge in the same units. 20 On the laser threshold. The asymptotics ( 1.11 ) hold for solutions to ( 1.1 ) with harmonic initial states X ( 0 ) ∈ Z r , and the Lebesgue measure | Z r | = 0. On the other hand, ( 1.12 ) hold for solutions with initial states from an open domain of attraction. Hence, by ( 1.13 ), the asymptotics ( 1.12 ) appear with a “nonzero prob- ability” in contrast to ( 1.11 ). This fact, provisionally , clarifies the existence of a laser thr eshold to ignite the laser action: the intensity of random pumping must be sufficiently large to bring the solution to the domain of attraction, and then the solution is captured in the domain with the single-frequency asymptotics. On the laser amplification. The equations ( 1.1 ) describe one molecule coupled to the Maxwell field. The limiting amplitudes of the Maxwell field in the asymp- totics ( 1.11 ), ( 1.12 ) do not depend on non-resonance harmonics in the pumping ( 1.4 ) with the frequencies Ω k  = ω . This means that the dynamics ( 1.1 ) acts as a filter , selecting only the resonant harmonics, that itself cannot explain the amplifi- cation of the Maxwell field in laser de vices. The amplification could be explained by a large number of activ e molecules, typically N ∼ 10 20 , under the traditional as- sumption that the molecules interact with the Maxwell field b ut do not interact with each other [ 28 ]. In this case, the amplitude of the total Maxwell field is multiplied by √ N ∼ 10 10 by the Law of Lar ge Numbers. This conclusion essentially depends on the fact that the phases of the ampli- tudes M r and M ∗ in asymptotics ( 1.11 )–( 1.12 ) for all molecules are approximately uniformly distrib uted. The uniformity takes place if it holds for phases of A e . This follo ws from the asymptotics and Remark 5.2 since e i φ e − i Ω t M differs from e − i Ω t M by a shift of time. Self-induced transparency . The harmonic states ( 5.6 ) means that the incident wa ve A e e − i Ω t results in the outgoing wa ve M ( t ) ≈ − A e e − i Ω t with the same ampli- tude and frequency , b ut with the phase jump π . This phenomenon resembles the self-induced transparency [ 33 ]. Refer ences [1] L. Allen, J.H. Eberly , Optical Resonance and T wo-Lev el Atoms, Do ver , Ne w Y ork, 1987. [2] F .T . Arecchi, Chaos and generalized multistability in quantum optics, Phys. Scr . 9 (1985), 85–92. [3] F . Arecchi, R. Bonifacio, Theory of optical maser amplifiers, EEE J ournal of Quantum Electr onics 1 (1965) , no. 4, 169–178. 21 [4] A.V . Babin, M.I. V ishik, Attractors of Evolution Equations, North-Holland Publishing Co., Amsterdam, 1992. [5] I. Bejenaru, D. T ataru, Global wellposedness in the energy space for the Maxwell–Schr ¨ odinger system, Commun. Math. Phys. 288 (2009), 145–198. [6] P .R. Berman, V .S. Malinovsk y , Principles of Laser Spectroscopy and Quan- tum Optics, Princeton Uni versity Press, Princeton, 2011. [7] F . Bloch, Nuclear induction, Phys. Rev . 70 (1946), no. 1, 460–474. [8] N. Bogoliubov , Asymptotic Methods in the Theory of Non-Linear Oscilla- tions, Gordon & Breach, Paris, 1961. [9] M.R. Candido, J. Llibre, D.D. Nov aes, Persistence of periodic solutions for higher order perturbed dif ferential systems via L yapunov–Schmidt reduction, Nonlinearity 30 (2017), 35–60. [10] M.R. Candido, J. Llibre, Stability of periodic orbits in the averaging theory: Applications to Lorenz and Thomas’ differential systems, Int. J. Bifur cat. Chaos 28 (2018), 1830007–14. [11] M. R. Candido, J. Llibre, C. V alls, New symmetric periodic solutions for the Maxwell–Bloch differential system, Math. Phys. Anal. Geom. 20 (2019), no. 16. [12] V .V . Chepyzho v , M.I. V ishik, Attractors for Equations of Mathematical Physics, AMS, Providence, RI, 2000. [13] J. Demaison, W . H ¨ uttner , Dipole Moments and Related Constants of Diamag- netic Molecules, Springer , Berlin, 1992. [14] R.P . Feynman, F .L. V ernon, Jr ., R.W . Hellwarth, Geometrical representation of the Schr ¨ odinger equation for Solving Maser Problems, J. Appl. Phys. 28 (1957), no. 1, 49–52. [15] C.C. Gerry , P .L. Knight, Introductory Quantum Optics, Cambridge University Press, Cambridge, 2004. [16] Y . Guo, K. Nakamitsu, W . Strauss, Global finite-energy solutions of the Maxwell–Schr ¨ odinger system, Comm. Math. Phys. 170 (1995), no. 1, 181– 196. [17] H. Haken, Laser Theory , Springer , Berlin, 1984. 22 [18] A. Haraux, Nonlinear Ev olution Equations – Global Behavior of Solutions, Lecture Notes in Mathematics 841, Springer , Berlin, 1981. [19] A.I. Komech, On quantum jumps and attractors of the Maxwell–Schr ¨ odinger equations, Annales mathematiques du Quebec 46 (2022), 139–159. [20] A. K omech, On periodic solutions for the Maxwell–Bloch equations, Physica D: Nonlinear Phenomena , 475 (2025), 134581. Open access. arXi v:2312.08180 [math-ph] https://www .sciencedirect.com/science/article/pii/S0167278925000600 [21] A. K omech, Quantum Mechanics: Genesis and Achiev ements, Springer , Dor- drecht, 2013. [22] A.I. K omech, Lectures on Quantum Mechanics and Attractors, W orld Scien- tific, Singapore, 2022. [23] A. K omech, E. K opylov a, Attractors of Hamiltonian Nonlinear Partial Differ - ential Equations, Cambridge Uni versity Press, Cambridge, 2021. [24] A. K omech, E. K opylov a, On single-frequency asymptotics for the Maxwell– Bloch equations: pure states, arXiv:2603.17888 [math.AP] [25] W .E. Lamb, jr, Theory of an optical maser , Phys. Rev . 134 (1964), no. 6a, A1429–A1450. [26] MatW eb https://www .matweb .com/search/datasheet.aspx?MatGUID=c8c56ad547 [27] M. Nakamura, T . W ada, Global existence and uniqueness of solutions of the Maxwell–Schr ¨ odinger equations, Comm. Math. Phys. 276 (2007), 315–339. [28] H. Nussenzveig, Introduction to Quantum Optics, Gordon and Breach, Lon- don, 1973. [29] J.J. Sakurai, Advanced Quantum Mechanics, Pearson Education, Incorpo- rated, 2006. [30] J.A. Sanders, F . V erhulst, J. Murdock, A v eraging Methods in Nonlinear Dy- namical Systems, Springer , Ne w Y ork, 2007. [31] M. Sar gent III, M.O. Scully , W .E. Lamb Jr, Laser Physics, Addison W esley , Reading, 1978. [32] M.O. Scully , M.S. Zubairy , Quantum Optics, Cambridge Univ ersity Press, Cambridge, 1997. 23 [33] A.E. Siegman, Lasers, Univ ersity Science Books, 1986. [34] W .T . Silfv ast, Laser Fundamentals, Cambridge Univ ersity Press, Cambridge, 2012. [35] R. T emam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, Springer , New Y ork, 1997. [36] J. W ei, X. W ang, X. Geng, Periodic and rational solutions of the reduced Maxwell–Bloch equations, Comm. Nonl. Sci. & Numer . Simulat. 59 (2018), 1–14. 24