Emission spectra and intrinsic optical bistability in a two-level medium
Scattering of resonant radiation in a dense two-level medium is studied theoretically with account for local field effects and renormalization of the resonance frequency. Intrinsic optical bistability is viewed as switching between different spectral patterns of fluorescent light controlled by the incident field strength. Response spectra are calculated analytically for the entire hysteresis loop of atomic excitation. The equations to describe the non-linear interaction of an atomic ensemble with light are derived from the Bogolubov-Born-Green-Kirkwood-Yvon hierarchy for reduced single particle density matrices of atoms and quantized field modes and their correlation operators. The spectral power of scattered light with separated coherent and incoherent constituents is obtained straightforwardly within the hierarchy. The formula obtained for emission spectra can be used to distinguish between possible mechanisms suggested to produce intrinsic bistability in experiments.
pacs:42.65.Pc, 32.50.+d, 42.50.Ct
It has long been known that optical properties of dense atomic ensembles or complex materials may be very different from those exhibited by independent atoms. Some of these systems can produce two different output signals responding to the same input intensity from a driving laser. For the entire range of applied laser intensities the output function would form a hysteresis loop. When no external feedback is required for this phenomenon to occur it is recognized as intrinsic optical bistability (IOB). Experimentally IOB has been observed in the optical response from rare earth ions and ion pairs in glasses and crystallines [1, 2, 3, 4, 5, 6, 7, 8]. However, its physical nature is still subject to debate. Most theoretical models suggest IOB originates from the fact that the transition frequency of a light emitter in a driven system is somehow renormalized and takes the effective value . The latter is a linear function of inversion where the renormalization constant follows from the coupling mechanism between the light emitter and its surrounding. Since is dependent, it provides the system with an intrinsic feedback and nonlinearity and, therefore, produces bistability in the steady state solution. Most commonly such renormalization is introduced through near dipole-dipole interactions (NDD) between emitters [9, 10, 11]. Alternatively, it can be achieved by considering model interactions for emitter pairs [12, 13, 14], coupling between a single emitter and the local vibration mode in its host medium [15, 16] or in an interacting charge-transfer exciton system . But, certainly, the IOB community has always been open to other ideas.
Originally IOB was predicted by Bowden with co–workers for a dense ensemble of two-level atoms [18, 19, 20, 21]. Currently, in a simple semiclassical view to the problem, is a result of the famous Lorentz local-field correction (LFC). It states that in a dense homogeneous ensemble exposed to a laser light each atom is actually driven by an effective or local field strength
where is the macroscopic field in the medium commonly approximated as the laser field and is the macroscopic polarization. If one substitutes this condition to the optical Bloch equations this corrected field can be manipulated to give the effective for which is the Lorentz frequency or NDD interaction parameter. It depends on the number density of atoms and the matrix element of the atomic transition dipole moment . In the literature is often referred to as Bowden-Lorentz redshift . Thus it must be widely understood that the Bowden type excitation-dependent renormalization of the resonance frequency follows straightforwardly from the locally effective Rabi frequency . Additionally, interpretations or enhancements of IOB involve cooperative processes [1, 2, 22] and laser heating . In particular, the work  pioneered the long series of papers studying local field effects in thin films. Various numerical models using nonlinear energy transfer processes in rare-earth-doped crystals have also been reported to show IOB behavior [23, 24, 25].
Most IOB theories were built to describe either atomic excitation or radiation intensity hysteresis. However, not much attention has been paid to studying spectral properties of the IOB response. In the steady state limit the “classic” resonance fluorescence of a still two-level atom in vacuum is known to be a function of atomic excitation or inversion . The line shape of the inelastic component of fluorescent light is determined by three frequencies: the spontaneous decay rate , the detuning between the laser and the transition frequencies , and the Rabi frequency [26, 27, 28]. It is apparent that as soon as the atom is placed in a host medium the frequencies take effective values , , and . Such metamorphose could be detected by an observer. We hypothesize that for an intrinsically bistable system with hysteresis the spectral patterns would be evidence of which frequency is “responsible”. This is potentially a good technique to identify the adequate IOB model.
Our work studies IOB as switching between fluorescence spectra and changing of spectral patterns with alternative mechanisms of bistability occurrence. As the first step we will study the original model of IOB based on the Lorentz local-field correction. We suggest a new method for calculations of emission spectra as a convenient and consistent addition to the conventional techniques . The equations to describe the non-linear interaction of a two-level atom in the ensemble with the laser field and the properties of the scattered light are derived from the Bogolubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy for reduced single particle density matrices of atoms and quantized field modes and their correlation operators. This method is a significant improvement to the atom-photon density operator formalism described in books [26, 29] and numerous papers.
2 BBGKY and Master Equation
We consider an ensemble of atomic particles and quantized field modes which we denote by subscript indices and respectively. The atoms are assumed to be “motionless” so that there is no interaction between them rather than via the field modes. The general structure of the Hamiltonian for such a system is given by
where and are the energies of free particles and describes their binary interactions. In order to introduce our approach in a general manner we do not specify the explicit terms until we finalize the model. The evolution and properties of the system are found from the many-particle total density operator which obeys the von Neumann equation:
However, in our study we need to find the solutions for the reduced single-particle operators and . The initial conditions for them can be defined using the fact that if at no interaction is observed can be factorized to a product of all single-atom and single-photon density operators. Further development of the single particle operators as constituents of the “unfactorable” system and the steady state mode are to be revealed. This goal can be approached if we substitute the von Neumann equation with the BBGKY hierarchy for reduced density and correlation operators [30, 31, 32]. Some issues of building such an hierarchy for atom-photon systems and proper decoupling of equations were discussed in [33, 34]. In this paper we will skip the cumbersome derivations and only point out the main steps to the analyzable system of equations. First of all, let us note that in order to treat the problem properly one needs to consider the actual number of particles and avoid using the thermodynamic limit where the problem is treated in terms of constant densities of material particles. We formally assume all atoms to be individuals and turn to their bulk concentration only when performing spatial integration. Correspondingly, the reduced matrices are obtained from by taking a partial trace, i.e.
From this definition it follows that in order to meet (3) the trace or even a partial trace of the correlation operator should be strictly zero:
With all these starting points the von Neumann equation could be equivalently replaced with an infinite series of coupled kinetic equations. The first two types of equations are precisely
Here, the l.h.s. terms constitute the von Neumann equations for single-particle operators coupled on the r.h.s. to each other and the remaining system through the atom-field correlations found from
This equation is written in the limit of the generalized second Born approximation [31, 34] for which the contributions from other two-particle correlations are neglected and only partial account for the three-particle correlations is taken:
The primes at the particle-reference indices are used to distinguish different particles of the same sort. The higher order correlations are beyond the scope of the current analysis. Equations (6)-(2) are closed and energy conservative. The single particles here are also coupled to the remaining system via the semi-averaged interaction operators or Hartree (mean field) terms in the effective Hamiltonians:
The averages must be understood as . Note that the interaction operator (2) corrected with the Hartree terms is the one that enables the correlation operators be strictly zero when traced over any correlating particle as stated in (5).
Now we can specify the light emitters as two-level atoms and for free particles write
The atom-photon interactions will be treated in the electric-dipole approximation
for which the dipole moment operator of an atom localized at some point is
and a single mode of the electric field at is given by
Here, are the population inversion, raising and lowering operators for atoms with transition frequency and is its dipole transition matrix element. Then and are the creation and annihilation operators that refer to a field mode with frequency , wave vector , , and polarization . denotes the quantization volume. The action of the classic external field is reflected by setting the initial condition for a number of selected modes:
where indicates a coherent state of the field. We will further use the space–time commutation relation
for the electric field
for which we introduced the field correlation tensor
where denote normal ordering and are the retarded and advanced Green’s tensors respectively. The important averages we use below are
Equation (6) can be shown to give the master equation if one eliminates its field operators using the standard techniques of quantum optics. First, it is beneficial to transform to the interaction picture for which
In this picture the remaining Hamiltonian operators (13) and (14) describe interactions and acquire the explicit exponential time dependencies and . Now we can substitute the formal integrals of (7) and (2) over time into (6). Then we assume . Retaining the terms which contain and only and making the use of (16)-(18) with we come to
The first commutator on the r.h.s. of (20) describes interaction of atom with the local field (21) in which is the external field. It originated properly in accordance with the condition defined in (15). The integrals in (21) describe contributions from the other atoms in the ensemble to the resulting effective field acting at position . Here we used the continuous medium approximation and replaced the summation over atoms with a proper integration over the sample volume . The prime reflects the property of BBGKY which excludes atom from the mean field (2) and eliminates its contribution to (21). We will apply the conventional approach to such integration where one should mind an “empty” sphere around atom .
The next two integrals in (20) describe atomic relaxation due to spontaneous emission. The last term makes account for atomic transitions induced by incoherent photons. However the latter processes are to be neglected in this paper. In order to make a comparison with the “classic” works on IOB we restrict our model to the case when there is no reabsorption of scattered light.
The equation for the local field does not make an account for the complete Green’s tensor which has the delta-function peculiarity at zero point [33, 36, 37] because this region has been removed from the integration volume. However, if we rewrite (21) with the full tensor and separate integration over volumes and the local field becomes easier to handle with for homogenous and isotropic media:
where we introduced macroscopic polarization of the sample . Equation (22) is exactly the well-known Lorentz local field condition (1) which describes the difference between the macroscopic field in a medium and the field acting on each atom therein. This is the first result of this paper which demonstrates that the master equation making an account for the local-field correction can be self-consistently obtained from BBGKY.
3 Bistability and emission spectra
The master equation obtained in section 2 can now be rewritten in a more transparent form. Without its last term (20) becomes independent and may be easily transformed to the IOB equation using the Born-Markov and the rotating wave approximations. Changing back to matrix in (19) and operating in the frame rotating with the laser frequency we get
where is the detuning, is the effective Rabi frequency and is the spontaneous decay rate. Since all operators and frequencies are associated with a single arbitrary atom we may now drop the unnecessary subscripts to simplify the notation. Given the model described above these parameters are exactly the following: , , and . The contribution of the Lamb shift is neglected.
This form of the master equation has been very well studied in the steady state limit and is known to describe IOB (see section 1). When (23) is written as the set of equations for matrix elements and resolved for the population difference function it produces the third-order algebraic equation:
where the polynomial coefficients are:
Coefficient in our case is precisely the NDD parameter that follows from (22). Consequently, we are to write , where and . Three real roots of (24) found for various laser field strengths () form the hysteresis loop of atomic excitation as shown in figure 1. Similarly, it is also the function which describes IOB in terms of the total or frequency-integrated power of fluorescent light as . The coexisting steady state solutions are notated as and and refer to the lower and the upper curves respectively.
Equation (24) can as well be derived from (23) in which and . As discussed earlier it could be, for example, either the case suggested in [15, 16] or in . In this paper we will not reproduce those models using BBGKY formalism as it would require chains of cumbersome transformations similar to those outlined in section 2. With no new features revealed it could be reasonable to introduce the effective detuning phenomenologically. In this approach we assume , where is an arbitrary factor. In any case with either effective detuning or Rabi frequency we use (24) with to find the source of fluorescent light.
Calculation of fluorescent spectra requires the use of (7) and (2). In [33, 34] the spectrum was defined as irreversible loss of radiation in a nonconservative system. It was shown that such treatment was necessary to calculate the time-dependent or transient mode spectra correctly. The spectra relevant to the steady-state hysteresis do not demand this kind of special care and can be evaluated within the set of equations (6)-(2) that conserve energy. However in this paper we will follow the approach of [33, 34] to demonstrate succession of the research techniques applied for related problems.
Let us note first that elimination of three-particle correlation in (2) produces operators similar to those found in r.h.s. of (20). Dropping the terms dependent on the number of photons we get the atomic contribution to the correlation function associated with the Lorentz local field and spontaneous emission. Besides, we will supply (7) and (2) with the field damping operator:
where is the loss rate at zero temperature. This damping could be introduced accurately within BBGKY if the system is completed with the field reservoir. However, since it is a well known result we simply write it in. Later on in our calculations the damping rate is considered to be .
Meanwhile we have a complete set of equations to study spectral properties of resonance fluorescence. The power density of a field mode is a function of
The time derivative of is defined by (7) with attachment of (25). Since we are to count photons leaving the system we need to consider the nonconservative term only. The latter is nothing but (25). It is also necessary to write the opposite sign to become an “observer”. In such a way we come to the function that describes the power of the emission signal absorbed by an idealized detector:
It must be understood that the photon counter is characterized by a continuity of . Now our task is to find in the steady state. The proper equation for the number of photons follows straightforwardly from (7):
where . The atom-field coupling is represented here by the coupling constant . In the steady state with evaluation of becomes obvious. Functions and are eliminated naturally using (7) and (25) transformed to the form
where we introduced the scale frequency and applied . In derivation of these equations we always assumed that all the approximations used to obtain the master equation (23) were also valid. In the steady state for light emitters with no account for the spatial part we get
where the off-diagonal matrix elements and are yet to be found from (23) and (30) respectively with the time derivatives put to zero. is a power density dimension factor. It is important to note that (30) is linear and can be solved explicitly as shown further. When the final expression becomes the famous formula for the fluorescence spectrum with coherent and incoherent parts. However, unlike the method described in [26, 29] we performed no artificial steps to separate the two parts. Therefore, the technique we present is convenient for studying spectral properties of systems exhibiting nonlinear response to external light for which analytical solution of (23) is hard or impossible to find.
4 Comparison of spectral patterns
As we have already pointed out the equations obtained in the previous section are suitable to describe two mechanisms of bistability, i.e., IOB produced by either or . The latter model will give (30) effective detuning and . In this section we will show the differences in spectral patterns of inelastic fluorescent emission expected for and based IOB models given and excitation curves are identical as plotted in figure 1. We can perform a formal generalization of (30) by writing and and deciding which is to represent an effective value later in the final expression. From (32) it follows that we need the function of which can be found from (30). It can be expanded to the system of four equations for four matrix elements . However, we can reduce the number of equations to three using the general property of a correlation operator (5) which in this case gives . Consequently, the set of equations can be represented in the matrix form as follows:
where M is the matrix:
and and are vectors. The first represents the correlation function in terms of its matrix elements
Vector has two constituents where represents the source of spontaneously emitted photons
and vector eliminates the Rayleigh singularity
Solution of (33) gives
In order to obtain the desired expression for let us multiply the numerator and denominator by the complex conjugate term and take the real part of the numerator:
The corresponding polynomial coefficients are:
Equation (39) with proper substitutes for either or is the main results of our paper. It determines the shapes of the inelastic emission spectra in a wide range of parameters. Below we present several plots of the spectrum profiles corresponding to four important points on the hysteresis curve in figure 1. Points and give the atomic excitations around the threshold value of that refers to the switch from the lower to the upper branch of the hysteresis loop. For a single two-level atom in the steady state one cannot expect the excitation to be as low as indicated by point in the limit of a strong field. The threshold point for the lower curve is sufficiently far in the saturation region. Consequently, the spectral function (39) for is notably different from conventional triplets [27, 28]. Figure 2 demonstrates the plots for and variants of IOB. One can see that the line splittings and peak values are sharply “outlined” by the respective models. The distance between the satellites and the central peak is over times greater in terms of units for IOB response conditioned by effective as compared to the Lorentz model. Moreover at this point the linear response spectrum is sharply distinctive as it is nothing but the “classic” Mollow’s profile.
For the upper branch at point 2 we have and the emission spectra plots are shown in figure 3. Since here , the contribution from the effective detuning is weak, i.e. , and . Obviously, the respective triplet coincides with the classic shape. In its turn the Lorentz IOB spectrum still reflects some difference between and . However, it will gradually vanish with further increase of . This can be very well seen from the strict relation
As we reduce the strength of excitation we move to point on the upper steady state curve. It corresponds to the weak threshold field , high excitation and small . For these parameters the emission spectra are shown in figure 4. Here we find smaller splitting between spectral lines as compared to the single atom picture. Opposite to point 1 the bistability triplet has a narrower range in comparison with IOB.
Further reduction of the pump strength will make the system transit to point at which . Since is proportional to the signal intensity should be very small and disadvantageous for registering the line shapes.
Let us consider the denominator in equation (39) for the emission spectrum:
It can be written as a product of the peak roots and two positive terms in the following form:
From (42) the triplet structure of fluorescence becomes clear. The central peak at is symmetrically surrounded by the satellites with maximums at . It follows straightforwardly from the fact that if function (39) reaches its maximum value. In addition to the conventional analysis of fluorescence these expressions contain the effective or renormalized frequencies. Equation (43) gives the exact position of the side peaks for both models of IOB. One may write the proper substitute for or to introduce the respective mechanism and obtain a picture similar to the one shown in figure 5. It is demonstrated here that the frequency shifts of the satellites in a strong fluorescence signal may vividly characterize which mechanism is likely to produce a saturation curve with an IOB hysteresis loop. Moreover, (43) is also suitable for studying joint contribution of effective parameters. However, this issue as well as calculation of to find a complex IOB picture is out of discussion in this paper. Finally, we may note that the Bowden-Lorentz type of IOB is best detected at point or, alternatively, on the lower branch as narrowing of the spectral range. At the same time the real renormalization of the atomic transition frequency could be very well seen when uprising to the threshold point .
We have performed a purely analytical study of IOB as switching between the patterns of fluorescence spectra. The local field induced hysteresis was shown to originate naturally for the model of a dense collection of like emitters when described thoroughly within the BBGKY hierarchy. The steady state spectra were found self-consistently within the developed formalism using the atom-field correlation matrix. We demonstrated that in terms of spectral analysis the contribution from the Lorentz local field cannot be described via introduction of an effective excitation depended detuning. In other words no evidence of actual renormalization of the atomic transition frequency can follow from observing light scattering in a monocomponent system. Within this concept the incoherent Mollow triplets would reflect abnormally small effective pumping of light emitters along the lower branch of the hysteresis loop. At the same time, a significant increase of the field induced splitting must be expected for the high excitation state. The latter though should gradually disappear in the saturation limit.
In order to make a comparison with the models suggesting renormalization of the transition frequency due to interactions of the atom with its complex environment we introduced this property phenomenologically. Since the IOB equation for atomic excitation is valid for both concepts the system parameters were chosen to give identical hysteresis loops. Analysis of the spectral profiles for this case revealed substantial differences in the character of line splitting at weak and moderate excitation. The results of our comparison are supported by analytical equations (39)-(40) and (43).
This analysis was carried out to emphasize the importance of spectroscopic data in studying optical response of complex materials. The problem could be described theoretically with the BBGKY formalism adjusted for atom-photon interactions. It is very advantageous to use this approach because the theory could be created from the first principles avoiding common phenomenology. The BBGKY method applied in this paper gave us a self-consistent set of equations with remarkable features which revealed some optical properties of the atomic system and the scattered field. It was due to the BBGKY concept and the general property of correlation function that we succeeded in reducing the number of equations to obtain the fluorescence spectrum function. The Rayleigh elastic constituent of the scattered light was automatically separated from the broad inelastic components which simplified our calculations greatly as compared to other spectroscopic techniques using the density matrix. The formula obtained for emission spectra can be used to distinguish among possible contributions from effective atomic frequencies in a medium showing IOB in experiments.
-  M. P. Hehlen, H. U. Güdel, Q. Shu, J. Rai, S. Rai, and S. C. Rand. Cooperative bistability in dense, excited atomic systems. Phys. Rev. Lett., 73(8):1103–1106, Aug 1994.
-  Markus P. Hehlen, Hans U. Güdel, Qize Shu, and Stephen C. Rand. Cooperative optical bistability in the dimer system Cs3Y2Br9:10Yb3+. The Journal of Chemical Physics, 104(4):1232–1244, 1996.
-  S. R. Lüthi, M. P. Hehlen, T. Riedener, and H. U. Güdel. Excited-state dynamics and optical bistability in the dimer system Cs3Lu2Br9:Yb3+. Journal of Luminescence, 76-77:447 – 450, 1998. Proceedings of the Eleventh International Conference on Dynamical Processes in Excited States of Solids.
-  Markus P. Hehlen, Amos Kuditcher, Stephen C. Rand, and Stefan R. Lüthi. Site-selective, intrinsically bistable luminescence of ion pairs in . Phys. Rev. Lett., 82(15):3050–3053, Apr 1999.
-  A. Kuditcher, M. P. Hehlen, C. M. Florea, K. W. Winick, and S. C. Rand. Intrinsic bistability of luminescence and stimulated emission in yb- and tm-doped glass. Phys. Rev. Lett., 84(9):1898–1901, Feb 2000.
-  Daniel R. Gamelin, Stefan R. Lüthi, and Hans U. Güdel. The role of laser heating in the intrinsic optical bistability of yb-doped bromide lattices. The Journal of Physical Chemistry B, 104(47):11045–11057, 2000.
-  S. M. Redmond and S. C. Rand. Intrinsic chromatic switching of visible luminescence in yb3,er3:cscdbr3. Opt. Lett., 28(3):173–175, 2003.
-  Jonathan M. Ward, Danny G. O’Shea, Brian J. Shortt, and Síle Nic Chormaic. Optical bistability in er-yb codoped phosphate glass microspheres at room temperature. Journal of Applied Physics, 102(2):023104, 2007.
-  Charles M. Bowden and Michael E. Crenshaw. Cooperativities in two-level systems. Optics Communications, 179(1-6).
-  Michael E. Crenshaw. Comparison of quantum and classical local-field effects on two-level atoms in a dielectric. Phys. Rev. A, 78(5):053827, Nov 2008.
-  Denis V. Novitsky. Local field effect as a function of pulse duration. Phys. Rev. A, 82(1):015802, Jul 2010.
-  O. Guillot-Noël, L. Binet, and D. Gourier. Towards intrinsic optical bistability of rare earth ion pairs in solids. Chemical Physics Letters, 344(5-6):612 – 618, 2001.
-  O. Guillot-Noël, L. Binet, and D. Gourier. General conditions for intrinsic optical bistability at the atomic and molecular scale: An effective spin-hamiltonian approach. Phys. Rev. B, 65(24):245101, May 2002.
-  O. Guillot-Noël, Ph. Goldner, and D. Gourier. Dynamics of intrinsic optical bistability in two weakly interacting quantum systems. Phys. Rev. A, 66(6):063813, Dec 2002.
-  F. Ciccarello, A. Napoli, A. Messina, and S. R. Lüthi. A new monomeric interpretation of intrinsic optical bistability observed in yb3+-doped bromide materials. Chemical Physics Letters, 381(1-2):163 – 167, 2003.
-  F. Ciccarello, A. Napoli, A. Messina, and S. R. Lüthi. A microscopic monomeric mechanism for interpreting intrinsic optical bistability observed in yb 3+ -doped bromide materials. Journal of Optics B: Quantum and Semiclassical Optics, 6(3):S118, 2004.
-  Ka-Di Zhu and Wai-Sang Li. The nonlinear optical response of interacting charge-transfer excitons with a weak probe in a strong pump field. Journal of Physics B: Atomic, Molecular and Optical Physics, 33(20):4255, 2000.
-  C. M. Bowden and C. C. Sung. First- and second-order phase transitions in the dicke model: Relation to optical bistability. Phys. Rev. A, 19(6):2392–2401, Jun 1979.
-  F. A. Hopf, C. M. Bowden, and W. H. Louisell. Mirrorless optical bistability with the use of the local-field correction. Phys. Rev. A, 29(5):2591–2596, May 1984.
-  Y. Ben-Aryeh, C. M. Bowden, and J. C. Englund. Intrinsic optical bistability in collections of spatially distributed two-level atoms. Phys. Rev. A, 34(5):3917–3926, Nov 1986.
-  Yacob Ben-Aryeh and Charles M. Bowden. Mirrorless optical bistability in a spacially distributed collection of two-level systems. Optics Communications, 59(3):224 – 228, 1986.
-  M. G. Benedict, V. A. Malyshev, E. D. Trifonov, and A. I. Zaitsev. Reflection and transmission of ultrashort light pulses through a thin resonant medium: Local-field effects. Phys. Rev. A, 43(7):3845–3853, Apr 1991.
-  M. A. Noginov, M. Vondrova, and D. Casimir. Analysis of intrinsic optical bistability in tm-doped laser-related crystals. Phys. Rev. B, 68(19):195119, Nov 2003.
-  Li Li, Xinlu Zhang, and Lixue Chen. Numerical analysis of intrinsic bistability and chromatic switching in tm 3+ single-doped systems under photon avalanche pumping scheme. Journal of Physics D: Applied Physics, 41(19):195105, 2008.
-  Li Li, Xinlu Zhang, Jinhui Cui, and Lixue Chen. A theoretical study of intrinsic optical bistability dynamics in tm 3+ /yb 3+ codoped systems with an upconversion avalanche mechanism. Journal of Optics A: Pure and Applied Optics, 11(10):105203, 2009.
-  P.A. Apanasevich. Foundations of the Theory of Interaction of Light (in Russian). Nauka i Tekhnika, Minsk, 1977.
-  B.R. Mollow. In E. Wolf, editor, Progress in optics XIX. Norht-Holland, Amsterdam, 1981.
-  Emil Wolf Leonard Mandel. Optical coherence and Quantum Optics. Cambridge Univerisity Press, New York, 1995.
-  S. Stenholm. Foundations of Laser Spectroscopy. John Wiley, Wiley-Interscience, New York, 1984.
-  A.I. Akhiezer and S.V Peletminsky. Methods of Statistical Physics. Pergamon, Oxford, 1981.
-  M. Bonitz. Quantum Kinetic Theory. B.G.Teubner, Stuttgart-Leipzig, 1998.
-  I.A. Kvasnikov. Statistical Physics. Editoral URSS, Moscow, 2002.
-  M. Gladush, A. Panteleev, and Vl. Roerich. Local field effects and stimulated multimode scattering of resonance radiation in a two-level medium. Journal of Experimental and Theoretical Physics, 103:206–217, 2006. 10.1134/S1063776106080048.
-  M.G. Gladush, A.A. Panteleev, and Vl.K. Roerich. Local field effects and stimulated multimode scattering of resonance radiation in a two-level medium. Acta Phys. Hung., 26:19, 2006.
-  A.I. Akhiezer and V.B. Berestetskiy. Quantum Electrodynamics. Nauka, Moscow, 1981.
-  Pedro de Vries, David V. van Coevorden, and Ad Lagendijk. Point scatterers for classical waves. Rev. Mod. Phys., 70(2):447–466, Apr 1998.
-  Michael Fleischhauer and Susanne F. Yelin. Radiative atom-atom interactions in optically dense media: Quantum corrections to the lorentz-lorenz formula. Phys. Rev. A, 59(3):2427–2441, Mar 1999.