{'fraction_non_alphanumeric': 0.03136922424756253, 'fraction_numerical': 0.011869436201780416, 'mean_word_length': 4.535588580367619, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


Introduction
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 ω A in a driven system is somehow renormalized and takes the effective valueω A = ω A − ζW . The latter is a linear function of inversion W where the renormalization constant ζ follows from the coupling mechanism between the light emitter and its surrounding. Since W isω A 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 [17]. 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,ω A 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 E M is the macroscopic field in the medium commonly approximated as the laser field and P is the macroscopic polarization. If one substitutes this condition to the optical Bloch equations this corrected field can be manipulated to give the effectivē ω A = ω A − ζ L W for which ζ L = (4π/3 )N|d| 2 is the Lorentz frequency or NDD interaction parameter. It depends on the number density of atoms N and the matrix element of the atomic transition dipole moment d. In the literature ζ L W is often referred to as Bowden-Lorentz redshift [10]. 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Ω = d · E L / . Additionally, interpretations or enhancements of IOB involve cooperative processes [1,2,22] and laser heating [6]. In particular, the work [22] 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 W . 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 W 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 [28]. 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.

BBGKY and Master Equation
We consider an ensemble of atomic particles and quantized field modes which we denote by subscript indices A and F 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 H A and H F are the energies of free particles and V AF 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 ρ(t) which obeys the von Neumann equation: However, in our study we need to find the solutions for the reduced single-particle operators ρ A (t) and ρ F (t). The initial conditions for them can be defined using the fact that if at t ≤ 0 no interaction is observed ρ(0) can be factorized to a product of all single-atom ρ A (0) and single-photon ρ F (0) 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.
Here, {S} and {S ′ } denote non-overlapping collections of species such that {S ′ + S} is the system described by ρ. The reduced multi-particle matrices are substituted with their cluster expansions [30,31]:

etc.
From this definition it follows that in order to meet (3) the trace or even a partial trace of the correlation operator g {S} 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)-(9) 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 X S = T r S (Xρ). Note that the interaction operator (10) 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 r A iŝ and a single mode of the electric field at r is given bŷ Here,σ z A ,σ ± A are the population inversion, raising and lowering operators for atoms with transition frequency ω A and d A is its dipole transition matrix element. Thenâ † k,s and a k,s are the creation and annihilation operators that refer to a field mode with frequency ω k , wave vector k, k = ω k /c, and polarizationǫ k,s . V 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 |α k,s indicates a coherent state of the field. We will further use the space-time commutation relation for the electric field is the Green's tensor following from the dyadic product of the field operators [35,28]. At the same time the dyads can be written aŝ for which we introduced the field correlation tensor where :: denote normal ordering and G ± 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 ±iω A t and ±iω k t. Now we can substitute the formal integrals of (7) and (9) over time t ′ < t into (6). Then we assumẽ ρ(t ′ ) ≈ρ(t). Retaining the terms which contain |d A | and |d A | 2 only and making the use of (16)- where The first commutator on the r.h.s. of (20) describes interaction of atom A with the local field (21) in which E 0 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 r A .
Here we used the continuous medium approximation and replaced the summation over atoms with a proper integration over the sample volume V ′ . The prime reflects the property of BBGKY which excludes atom A from the mean field (10) and eliminates its contribution to (21). We will apply the conventional approach to such integration where one should mind an "empty" sphere around atom A. 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 V and δV = V − V ′ the local field becomes easier to handle with for homogenous and isotropic media: where we introduced macroscopic polarization of the sample P = N d A ′ . 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.

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 ω L we get d dt 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: ∆ = ω A − ω L , Ω = d · E L / , and Γ = 2/ d · ImG(0, ω A ) · d = (4/3 )|d| 2 k 3 . 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 ρ A matrix elements and resolved for the population difference function W = −2 σ z 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 Ω = d · E 0 / and ζ = ζ L = (4π/3 )N|d| 2 = πNΓ/k 3 . Three real roots of (24)  Equation (24) can as well be derived from (23) in which ∆ →∆ = ∆ − ζW and Ω = Ω. As discussed earlier it could be, for example, either the case suggested in [15,16] or in [17]. 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 ζ = ζ m = mΓ, where m is an arbitrary factor. In any case with either effective detuning or Rabi frequency we use (24) with ζ = ζ L = ζ m to find the source of fluorescent light.
Calculation of fluorescent spectra requires the use of (7) and (9). 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)-(9) 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 g AF F ′ in (9) 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 (9) with the field damping operator: where η k 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 η k → 0. 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 ρ F 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 η k . Now our task is to find η k â † k,sâ k,s in the steady state. The proper equation for the number of photons follows straightforwardly from (7): where g + k = ξ k (r A )T r F (â k,s g AF ). The atom-field coupling is represented here by the coupling constant ξ k (r) = 2πω k / V(d ·ǫ k,s )e ik·r . In the steady state with d â † k,sâ k,s /dt = 0 evaluation of A k becomes obvious. Functions â k,s and â † k,s are eliminated naturally using (7) and (25) transformed to the form Function g + k and its conjugate are found from the equations obtained from (9) and (25): where we introduced the scale frequency ν k = ω k −ω L and applied η k → 0. 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 N light emitters with no account for the spatial part we get Finally using lim κ→0 1 π κ x 2 + κ 2 = δ(x) to substitute (31)

into (27) the expression for the fluorescence spectrum becomes
where the off-diagonal matrix elements ρ 12 = σ + and (g + k ) 12 = T r A (σ + g + k ) are yet to be found from (23) and (30) respectively with the time derivatives put to zero. F k 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. The solid line represents emission spectrum for ζ L = 50Γ, ζ m = 0, Dotted line -for ζ m = 50Γ, ζ L = 0. The gray shape describes the emission spectrum for an independent two-level atom. A 0 k is maximum value of the emission power function for a single independent atom in the saturation limit.

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∆ = ∆ − ζ m W andΩ = Ω. In this section we will show the differences in spectral patterns of inelastic fluorescent emission expected forΩ and∆ based IOB models given ζ L = ζ m 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 Re (g + k ) 12 which can be found from (30). It can be expanded to the system of four equations for four matrix elements (g + k ) 11 , (g + k ) 12 , (g + k ) 21 , (g + k ) 22 . However, we can reduce the number of equations to three using the general property of a correlation operator (5) which in this case gives (g + k ) 11 + (g + k ) 22 = 0. Consequently, the set of equations can be represented in the matrix form as follows: where M is the matrix: and g and q are vectors. The first represents the correlation function in terms of its matrix elements Vector q = |ξ k | 2 s − |ξ k | 2 ρ 21 ρ A has two constituents where s represents the source of spontaneously emitted photons and vector ρ A eliminates the Rayleigh singularity Solution of (33) gives In order to obtain the desired expression for Re (g + k ) 12 let us multiply the numerator and denominator by the complex conjugate term det M * 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 1 and 2 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 1 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 ρ 22 = ρ low 22 (Ω ↑ ) 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 10 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 ρ 22 = ρ up 22 (Ω ↑ ) and the emission spectra plots are shown in figure 3. Since here W ≈ 0, 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 that follows from (23). As one can see in (39) the line shape depends on |Ω| 2 and the relation to the spectral properties of a single atom is determined by (40). As we reduce the strength of excitation we move to point 3 on the upper steady state curve. It corresponds to the weak threshold field Ω ↓ , high excitation and small W .
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 4 at which W ≈ 1. Since A k is proportional to {ρ low 22 (Ω ↓ )} 2 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: 1 It can be written as a product of the peak roots and two positive terms in the following form: where and From (42) the triplet structure of fluorescence becomes clear. The central peak at ν k = 0 is symmetrically surrounded by the satellites with maximums at ν k = ±ν p . It follows straightforwardly from the fact that if (ν k ± ν p ) 2 = 0 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 |Ω| 2 or∆ 2 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 ρ 22 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 3 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 1.

Summary
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.