# Steady-state photoluminescence and nanoscopy of two near-identical emitters
with dipole-dipole coupling

Natalia A. Lozing  Moscow State Pedagogical University , 1-1 M Pirogovskaya
St, Moscow 119991, Russia  Ekaterina A. Smirnova  Moscow State Pedagogical
University , 1-1 M Pirogovskaya St, Moscow 119991, Russia  Vladimir K. Roerich
Microchip Technology Inc, 3870 N First St, San Jose CA, United States 95134
Maxim G. Gladush  Moscow State Pedagogical University , 1-1 M Pirogovskaya St,
Moscow 119991, Russia

###### Abstract

We report progress in the theory of photoluminescence and light scattering by
two closely spaced particles. This study is based on our original method to
derive the master equation for a system of coupled quantum emitters driven by
a cw laser. We emphasise on the case when the emitters are two-level systems
but notably different in their transition frequencies and transition moments.
The master equation is shown to describe the dipole-dipole coupling and the
cooperative entanglement between the particles naturally. It is provided by
the use of the Bogolyubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchies for
reduced density matrices and correlation operators of the material and
photonic subsystems. Tackling the hierarchies has also provided us with an
elegant way to calculate the photoluminescence excitation spectra subject to
the arrangement of the excitation and the position of the detector against the
emitters. Our findings may be important for developing novel possibilities for
the all-optical nanoscopy and may help construction of entangled systems to be
implemented in quantum technologies.

††preprint: APS/123-QED

##  I Introduction

Light emission by photoluminescence from systems of interacting atom- or
molecule-like particles has been the subject of extensive research in quantum
optics and related areas for a few decades. It is now widely understood that
ensembles of light emitters may be organized in a certain fine manner to show
cooperativity and entanglement. This type of collective behaviour is a result
of the interparticle interactions. The interactions produce collective
excitations alternative to a sum of excited single emitters. In turn, the
excitations are generally expected to damp out through some tricky scenario
[1, 2] and give birth to light with a few distinctive properties. The research
interests are mostly associated with the super- and sub-radiance phenomena
[1]. Both follow from formation of the entangled symmetric and antisymmetric
eigenstates which are the fast and slow decaying collective excitations
respectively. The entanglement of emitters and the dipole-dipole coupling
between them are known to give rise to new resonances in the emission and
absorption spectra [3, 4, 5, 6, 7, 8]. All these properties are currently of
special attention in several research communities. For instance, construction
of the entangled states has become a stand-alone topic in the contemporary
studies for quantum computing and the quantum information science [9, 10, 11].
The sub-radiative states are being studied for different configurations of
emitters and checked as candidates for the elements of optical and quantum
memory devices [12].

The recent surge of interest in cooperative emission has been fueled by the
advances in single molecule and single emitter spectroscopy. It allows
detection and characterization of two and more fluorescing particles with
subwavelength separation. Since the pioneering papers [13, 14] a number of
experimental techniques have now evolved into several routines which help
interdisciplinary studies. These include spectroscopy of single organic
molecules [15, 16], quantum dots [17], ions [18], and light-harvesting
pigment-protein complexes [19]. Observation of single quantum emitters reached
remarkable efficiencies and precisions shortly after it had become possible to
detect zero-phonon lines (see [20] for a review). In most cases it works best
at low temperatures so the vibrational excitations are “frozen” and the
emitters can be thought of as quasi two-level systems with bright
fluorescence. These spectroscopies have been so far the only successful
methods to demonstrate the cooperativity of only two particles located at a
small (less than the radiation wavelength) distance from one another. A
demonstration of coherent optical dipole-dipole coupling of two individual
molecules separated by $12$ nanometers in an organic crystal was reported in
[8]. At the same time, more novel techniques to measure small separations with
a fine accuracy are under way. For example, a new experimental method for
direct measurements of the separation between two semiconductor colloidal
quantum dots (down to $20-10$ nm) from fluorescent images was recently
described in [21]. The growing ability to observe the relationship between the
distance and changes in absorption and emission requires a better theoretical
description of small cooperative ensembles. Besides, recent progress in
preparation of the samples makes it possible to manage and arrange particles
in one-, two-, and three-dimensional structures at micro- and nano-scales.
Many advanced procedures are currently practiced to successfully trap and
confine small groups of emitters within a desired geometry, e.g., ions in
magneto-optical traps [22, 23], artificial atoms in a transmission line [24],
quantum dots in liquids [25], etc.

A system of two coupled light emitters has always been a special case of
cooperativity. It has the least number of parameters which makes it a perfect
object for developing a theory or interpreting an experiment. There have been
numerous studies of the dynamics, entanglement, and fluorescence from paired
emitters (for a review see, e.g., [26]). However, despite the long history of
research, many theories had limitations. Most of the research focused on
separate specific cases, which either allowed obtaining analytical expressions
or simulating a particular situation. At the same time, since the early papers
[27] and [28] the largest number of publications has been devoted to the study
of the cooperative spontaneous decay. Because the problem is geometry-
dependent, several authors have succeeded in proposing ingenious indirect ways
to observe cooperativity and determine the location of emitters with respect
to external conditions like in [29]. In [30] the time evolution of singly and
doubly excited states of a system of two spatially separated qubits (two-level
atoms) coupled to a one-dimensional waveguide was studied. The procedure of
subwavelength optical microscopy of closely spaces particles was reported in,
for example, [31, 32, 33] and [34]. The exchange between symmetric and
antisymmetric states due to the spatial variation of the applied laser
excitation was studied in works [35, 36]. In the work [37] the system of two
non-identical atoms showed the possibility of population trapping effect.
Another direction of the research has been devoted to studying cooperative
couples in presence of a continuous excitation [38, 39, 40]. In works [41, 7]
the total fluorescence intensity, fluorescence spectra and correlation
functions were calculated in cases of running and standing waves, but
conventionally only a case of two identical atoms was considered, and only a
limited geometry of interaction was investigated.

To date pairs of organic molecules and semiconductor quantum dots are
considered being among most favored candidates for creating the smallest, best
theoretically studied and most stable cooperative systems. Both types of
particles can be easily imbedded into a transparent host material and fixed in
unchangeable positions for unlimited time. The problem of controlling the
distance between particles can be solved by preparing a variety of pairs with
stochastically different separations. Controlling the distance in one pair is
also a solvable problem with few theoretical restrictions – it just requires
the use of rather tricky manipulations. For such cooperative couples, two
inevitable circumstances need to be considered. The first is that both types
of emitters are known for their wide inhomogeneous broadening [20]. The second
is that a sufficiently strong influence of the host matrix is possible [42].
Regardless of how these circumstances manifest in a particular case the
emitters in a pair cannot be considered completely identical. In this paper,
we study the steady-state cooperative emission of light and how it depends on
the parameters of the incident pump field. We consider a system of two near-
identical two-level emitters that are driven by a cw laser beam. We derive the
complete cooperative master equation for the atomic density matrix which is
adjustible to the excitation geometry and appears from the conventional
electric-dipole interaction Hamiltonian. The master equation is a part of the
decoupled Bogolyubov-Born-Green-Kirkwood-Yvon hierarchy (BBGKY) which also
provides the equations for the total photoluminescence intensity. Finally we
show how a rotation of the excitation polarization may help restoring the
orientation of the couple and the individual parameters of the emitters.

This paper is organized as follows. In Sec. II we make a review of the
theoretical formalism, state the problem and formulate the basic assumptions.
In Sec. III we derive the master equation in the differential-integral from
and the field correlation tensor. Sec. IV describes reduction of the master
equation to a particular solvable form. In Sec. V we derive the equation for
the total emission intensity and calculate the photoluminescence excitation
spectra. We show three resonances typical for pure two-level emitter and how
they change with different geometries of excitation and observation.

##  II Basic Equations and Assumptions

###  II.1 BBGKY hierarchies for emitters and photons

The evolution of the system which includes the emitters and the field is
described by a many-particle density matrix $\rho(t)$, which obeys the von
Neumann equation:

$i\hbar\frac{d{\rho}}{dt}=[\widehat{H},{\rho}],$ (1)

where $\widehat{H}$ is the full Hamiltonian. We use the conventional notations
for the unit imaginary number, reduced Plank constant, and commutator
brackets. The Hamiltonian $\widehat{H}$ is composed of two parts
$\widehat{H}=\widehat{H}_{0}+\widehat{V}$, where $\widehat{H}_{0}$ is the
energy of free particles
$\widehat{H}_{0}=\sum_{i}\widehat{H}_{i}+\sum_{f}\widehat{H}_{f}$, where $i$
numbers all material particles and $f$ numbers field modes, while
$\widehat{V}$ is the interaction operator.

The first step in our transformations is to replace Eq. (1) with the
equivalent hierarchy of equations for reduced density matrices of the emitting
system $\rho_{\\{a\\}}$, the field modes $\rho_{f}$, and the higher order
many-particle density matrices. This can be performed using the fundamental
reduction property
$\rho_{\left\\{S\right\\}}=Tr_{\left\\{S^{\prime}\right\\}}\rho$, where
$\\{S\\}$ and $\\{S^{\prime}\\}$ \- arbitrary sets of elements of the system
such that $\\{S+S^{\prime}\\}$ is a complete set of photons and material
particles. This way we obtain an infinite system of coupled equations:

$i\hbar\frac{d}{dt}\rho_{\\{a\\}}-\Big{[}\sum_{i}\widehat{H}_{i},\rho_{\small\\{a\\}}\Big{]}=\sum_{f}Tr_{f}\Big{[}\sum_{i}\widehat{V}_{if},\rho_{\\{a\\}f}\Big{]},$
$i\hbar\frac{d}{dt}\rho_{f}-\left[\widehat{H}_{f},\rho_{f}\right]=Tr_{\\{a\\}}\Big{[}\sum_{i}\widehat{V}_{if},\rho_{\\{a\\}f}\Big{]},$
(2)
$i\hbar\frac{d}{dt}\rho_{\\{a\\}f}-\Big{[}\sum_{i}\widehat{H}_{i}+\widehat{H}_{f}+\sum_{i}\widehat{V}_{if},\rho_{\\{a\\}f}\Big{]}\\\
=\sum_{f^{\prime}\neq
f}Tr_{f^{\prime}}\Big{[}\sum_{i}\widehat{V}_{if^{\prime}},\rho_{\\{a\\}ff^{\prime}}\Big{]},$
$i\hbar\frac{d}{dt}\rho_{ff^{\prime}}-\left[\widehat{H}_{f}+\widehat{H}_{f^{\prime}},\rho_{ff^{\prime}}\right]\\\
=Tr_{\\{a\\}}\Big{[}\sum_{i}\widehat{V}_{if}+\sum_{i}\widehat{V}_{if^{\prime}},\rho_{\\{a\\}ff^{\prime}}]\\\
\text{etc.}$

It is advantageous to replace the many-particle density matrices with
corresponding correlation operators using the cluster expansion [43]. Each
expansion is a superposition of the product of single-particle density
matrices and correlation operators, satisfying
$Tr_{\\{y\\}}g_{\\{x\\}}~{}=~{}0,\;\;\\{y\\}~{}\in~{}\\{x\\}$, as follows:

$\displaystyle\rho_{ff^{\prime}}=\rho_{f}\rho_{f^{\prime}}+g_{ff^{\prime}},\;\;\;\rho_{\\{a\\}f}=\rho_{\\{a\\}}\rho_{f}+g_{\\{a\\}f},$
$\displaystyle\rho_{\\{a\\}ff^{\prime}}=\rho_{\\{a\\}}\rho_{f}\rho_{f^{\prime}}+\rho_{\\{a\\}}g_{ff^{\prime}}+\rho_{f}g_{\\{a\\}f^{\prime}}$
$\displaystyle+\rho_{f^{\prime}}g_{\\{a\\}f}+g_{\\{a\\}ff^{\prime}}$ etc.

Substituting the expansions into Eq. (2) produces the hierarchy with one
emitter equation

$i\hbar\frac{d}{dt}\rho_{\\{a\\}}-\Big{[}\sum_{i}\overline{{H}}_{i},\rho_{\\{a\\}}\Big{]}=\sum_{f}Tr_{f}\Big{[}\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\Big{]}$
(3)

and an infinite number of equations for the field modes

$i\hbar\frac{d}{dt}\rho_{f}-\left[\overline{H}_{f},\rho_{f}\right]=Tr_{\\{a\\}}\Big{[}\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\Big{]}.$
(4)

The left hand sides of Eqs. (3) and (4) contain the effective energy
operators. They are the sums of the energy of free particles and Hartree
potentials, which are partial quantum-mechanical averages
$\langle\widehat{O}\rangle_{S}=Tr_{S}\widehat{O}\rho_{S}$ with respect to
particular species:

$\displaystyle\overline{H}_{i}=\widehat{H}_{i}+\sum_{f}\langle\widehat{V}_{if}\rangle_{f},$
$\displaystyle\overline{H}_{f}=\widehat{H}_{f}+\sum_{i}\langle\widehat{V}_{if}\rangle_{i}.$

These Hartree potentials provide the first-order coupling between the material
particles and the field. In Eq. (3) one finds the energy of interaction
between each emitter and the mean electric field while in equation (4) this
potential refers to interaction of each field mode with the mean induced
dipole. The second order coupling arises from the right hand sides where both
equations contain the emitters-field correlation $g_{\\{a\\}f}$. The equation
for $g_{\\{a\\}f}$ follows from the similar algebra applied to Eq. (2) and the
cluster expansion. We write it in a form convenient for further
transformations:

$i\hbar\frac{d}{dt}g_{\\{a\\}f}-\left[\sum_{i}\overline{H}_{i}+\overline{H}_{f},g_{\\{a\\}f}\right]+\widehat{S}=\widehat{I}+\widehat{L}+\widehat{\Pi}+\widehat{C}.$
(5)

The commutator in the left hand side contains the same energy operators as in
Eqs. (3) and (4). Operator $\widehat{S}$ combines corrections to the Hartree
terms which are three commutators:

$\displaystyle\begin{split}\widehat{S}=&\left[\sum_{i}\langle\widehat{V}_{if}\rangle_{i}+\langle\widehat{V}_{if}\rangle_{f},g_{\\{a\\}f}\right]+\rho_{f}Tr_{f}\left[\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\right]\\\
&+\rho_{\\{a\\}}Tr_{\\{a\\}}\left[\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\right].\end{split}$

The four right hand side operators represent the hierarchical structure of how
the emitter-field correlation is coupled to the density matrices and the
higher order correlations:

$\displaystyle\begin{split}\widehat{I}=&\left[\sum_{i}\overline{V}_{if},\rho_{\\{a\\}}\rho_{f}\right],\\\
\widehat{L}=&\left[\sum_{i}{V}_{if},g_{\\{a\\}f}\right],\\\
\widehat{\Pi}=&\sum_{f^{\prime}\neq
f}Tr_{f^{\prime}}\left[\sum_{i}\widehat{V}_{if},\rho_{\\{a\\}}g_{ff^{\prime}}\right],\\\
\widehat{C}=&\sum_{f^{\prime}\neq
f}Tr_{f^{\prime}}\left[\sum_{i}\widehat{V}_{if},g_{\\{a\\}ff^{\prime}}\right].\end{split}$

This hierarchy of couplings corresponds to the formal representation of the
BBGKY described in [43].

The first operator $\widehat{I}$ is a pure inhomogeneity formed by the
operators of the orders lower than $g_{\\{a\\}f}$, i.e., the product of the
reduced density matrices. It contains the interaction operators $V_{if}$
corrected with the corresponding components of the Hartree potentials:

${\overline{V}}_{if}=\widehat{V}_{if}-\langle\widehat{V}_{if}\rangle_{i}-\langle\widehat{V}_{if}\rangle_{f}.$

Operator $\widehat{L}$ is homogeneous with respect to $g_{\\{a\\}f}$ and is
known to collect so-called ladder terms. These are the commutators between the
quantum mechanical interactions of the emitters with the field mode of the
attention and the correlation operator. Obviously this term is of the higher
order as compared to the left hand side of Eq. (5). The other operators are
inhomogeneous terms which provide the links to the other types of
correlations. Operator $\widehat{\Pi}$ is a collection of the terms that
describe contributions from the correlations between two different field
modes. The term $\widehat{C}$ is the highest order inhomogeneity, i.e.,
contains the three particle correlation which is the bridge to the rest of the
BBGKY.

Equation (5) suggests the photon-photon correlations must be included. The
corresponding equations follow the above structure and are the following:

$i\hbar\frac{d}{dt}g_{ff^{\prime}}-\left[\overline{H}_{f}+\overline{H}_{f^{\prime}},g_{ff^{\prime}}\right]=\widehat{\Pi}^{\prime}+\widehat{C}^{\prime},$
(6)

where

$\displaystyle\widehat{\Pi}^{\prime}=Tr_{\\{a\\}}$
$\displaystyle\left[\sum_{i}\widehat{V}_{if},\rho_{f}g_{\\{a\\}f^{\prime}}\right]$
$\displaystyle+Tr_{\\{a\\}}\left[\sum_{i}\widehat{V}_{if^{\prime}},\rho_{f^{\prime}}g_{\\{a\\}f}\right],$
$\widehat{C}^{\prime}=Tr_{\\{a\\}}\bigg{[}\sum_{i}\widehat{V}_{if}+\sum_{i}\widehat{V}_{if^{\prime}},g_{\\{a\\}ff^{\prime}}\bigg{]}.$

For the purposes of this study it is sufficient to remain within the second
order operators in the inhomogeneous terms and in the low order homogeneous
terms, that is, we set $\widehat{L}=\widehat{C}=\widehat{C}^{\prime}=0$ and
disregard the impact of contributions of all higher order correlations between
the subsystems.

It means we go further beyond the Born approximation and take the advantage of
building the hierarchy either in polarization or ladder approximation. In
either case we have a closed system of equations for reduced density and
correlation matrices, and, consequently, can proceed with solving a self-
consistent problem by stating the necessary initial conditions.

Further it is beneficial to transform to the interaction picture in which the
total density matrix is

$\widetilde{\rho}=exp\left\\{\frac{i}{\hbar}\widehat{H}_{0}t\right\\}\rho\>exp\left\\{-\frac{i}{\hbar}\widehat{H}_{0}t\right\\}.$

Hereinafter, we will write the operators without the tildes. Now we have come
to a general closed system of equations for density matrices and correlation
operators which we will use in our analysis:

$i\hbar\frac{d}{dt}\rho_{\\{a\\}}-\sum_{f}\left[\sum_{i}\langle\widehat{V}_{if}\rangle_{f},\rho_{\\{a\\}}\right]\\\
=\sum_{f}Tr_{f}\left[\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\right],$ (7)
$i\hbar\frac{d}{dt}\rho_{f}-\left[\sum_{i}\langle\widehat{V}_{if}\rangle_{i},\rho_{f}\right]=Tr_{\\{a\\}}\left[\sum_{i}\widehat{V}_{if},g_{\\{a\\}f}\right],$
(8)
$i\hbar\frac{d}{dt}g_{\\{a\\}f}=\left[\sum_{i}\overline{V}_{if},\rho_{\\{a\\}}\rho_{f}\right]\\\
+\sum_{f^{\prime}\neq
f}Tr_{f^{\prime}}\left[\sum_{i}\widehat{V}_{if^{\prime}},\rho_{\\{a\\}}g_{ff^{\prime}}\right],$
(9)
$i\hbar\frac{d}{dt}g_{ff^{\prime}}-\left[\sum_{i}\langle\widehat{V}_{if}+\widehat{V}_{if^{\prime}}\rangle_{i},g_{ff^{\prime}}\right]\\\
=Tr_{\\{a\\}}\left[\sum_{i}\widehat{V}_{if},\rho_{f}g_{\\{a\\}f^{\prime}}\right]+Tr_{\\{a\\}}\left[\sum_{i}\widehat{V}_{if^{\prime}},\rho_{f^{\prime}}g_{\\{a\\}f}\right].$
(10)

In this system Eq. (7) is unique while Eqs. (8)-(10) are multiple as they
refer to the photonic modes. Our further developments will be introducing
formal solutions of Eqs. (8)-(10) into (7).

###  II.2 Statement of the problem for two-particle systems

In this work we consider a pair of quantum emitters excited by a cw laser
light. Each emitter is approximated as a two-level atom with a transition
frequency $\omega_{i}$ and a transition dipole moment ${\bf d}_{i}$, where
$i=1,2$, such that $\omega_{1}\neq\omega_{2}$, $|{\bf d}_{1}|\neq|{\bf
d}_{2}|$ and $|\omega_{1}-\omega_{2}|\ll\omega_{1},\omega_{2}$. The excitation
light beam may be, in principle, of arbitrary composition but in this paper as
a convention assumption for a cw laser it is treated as a linearly polarized
plane wave with the wave vector ${\bf k}_{L}$ and the electric field strength
amplitude ${\bf E}_{L}$. The atoms are viewed as attached to static positions
separated by a distance represented by a radius-vector ${\bf r}_{12}$, that is
a parameter to change arbitrary. The photoluminescence of the atomic pair is
to be calculated for arbitrary directions ${\bf r}_{D}$ from the midpoint of
the system. Figure 1 illustrates this statement of the problem and possible
arrangements for the excitation and detection means. In this diagram we set
the Cartesian system such the emitters are placed along the $Z$-axis at equal
distances from the origin. The spatial arrangement (relative orientation) for
the excitation system ${\bf k}_{L}$ and ${\bf E}_{L}$, vectors ${\bf r}_{12}$
and ${\bf r}_{D}$ may be fully defined by setting independent values to four
angles: $\xi$, $\theta$, $\varphi$, and $\psi$. The excitation scheme and the
directions of induced dipole moments are determined by two polar angles $\xi$
and $\theta$ that refer to ${\bf k}_{L}$ and ${\bf E}_{L}$ orientations
respectively. The azimuthal position of ${\bf d}_{i}$ (not marked on the
figure) is a function $f(\xi,\theta)$ which is calculated given that ${\bf
d}_{i}$ and ${\bf E}_{L}$ are collinear and ${\bf k}_{L}$ lies in the
$ZY$-plane. The inclination of ${\bf d}_{i}$ is dependent on $\xi$ and ranges
within $|\pi/2-\xi|\leq\theta\leq\pi-|\pi/2-\xi|$. Direction to the observer
${\bf r}_{D}$ may be varied by hand defining the values for polar and
azimuthal angles $\varphi$ and $\psi$. Consequently, we have sufficient data
to calculate the angle $\beta(\xi,\theta,\varphi,\psi)$ between vectors ${\bf
r}_{D}$ and ${\bf d}_{i}$ that we refer to as the ‘observation parameter’.

Figure 1: The excitation/detection diagram. Two emitters are fastened at a
distance $r_{12}=|{\bf r}_{12}|$ apart. The excitation is described by vectors
${\bf k}_{L}$ and ${\bf E}_{L}$ (${\bf E}_{L}\perp{\bf k}_{L}$) and is
adjustable by varying polar angles $\xi$ and $\theta$. The dipole moments
${\bf d}_{i}$ are defined by polar angle $\theta$. Detector position is
defined by ${\bf r}_{D}$ with polar $\varphi$ and azimuthal $\psi$ angled. The
angle between ${\bf r}_{D}$ and ${\bf d}_{i}$ is $\beta$.

The free Hamiltonian of the system is:

$\widehat{H}_{0}=\sum_{i}\widehat{H}_{i}+\sum_{f}\widehat{H}_{f}=\sum_{i}\hbar\omega_{i}\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{i}+\sum_{f}\hbar\omega_{f}\hat{a}^{\dagger}_{f}\hat{a}_{f},$

where we have two types of particle species, i.e., the emitters and the field
photons. The first sum represents the energies of free material particles with
raising $\hat{\sigma}_{i}^{+}=|e_{i}\rangle\langle g_{i}|$ and lowering
$\hat{\sigma}_{i}^{-}=|g_{i}\rangle\langle e_{i}|$ operators, where
$|e_{i}\rangle$ and $|g_{i}\rangle$ are the eigenvectors for the excited and
the ground states, index ’$i$’ numbers all the emitters ($i\in\\{{1,2}\\}$).
The second sum represents the energy of a free quantized electromagnetic
field, $\hat{a}_{f}^{\dagger}$ and $\hat{a}_{f}^{-}$ are photon creation and
annihilation operators respectively. Index $f=\\{{\bf k},\alpha\\}$ determines
the mode frequency and direction (${\bf k}_{f}$ \- wave vector, $|{\bf
k}_{f}|=\omega_{f}/c$) and polarization $(\alpha\in\\{1,2\\})$.

Operator $\widehat{V}$ describes the interaction of each emitter with
electromagnetic field modes in the electric dipole approximation:

$\widehat{V}=\sum_{f}\sum_{i}\widehat{V}_{if},\;\;\;\;\
\widehat{V}_{if}=-\widehat{\bf d}_{i}\cdot\widehat{\bf E}_{f},$

where $\widehat{\bf d}_{i}$ – quantum dipole moment operator:

$\widehat{\bf d}_{i}=\hat{\sigma}^{-}_{i}{\bf d}_{i}+\hat{\sigma}^{+}_{i}{\bf
d}^{*}_{i},$

and $\widehat{\bf E}_{f}$ \- quantized electric field mode operator:

$\widehat{\bf E}_{f}\left({\bf r}\right)=i\lambda_{f}{\bf
e}_{f}\left(\hat{a}_{f}e^{i{\bf k}_{f}\cdot{\bf
r}}-\hat{a}_{f}^{\dagger}e^{-i{\bf k}_{f}\cdot{\bf r}}\right),$

where ${\bf e}_{f}$ – unit polarization vector,
$\lambda_{f}=\sqrt{2\pi\hbar\omega_{f}/V}$ – coupling constant in the
quantization volume $V$.

In the interaction representation the dipole moment operator and the electric
field mode operator acquire an explicit time dependencies:

$\widehat{\bf d}_{i}(t)=\hat{\sigma}^{-}_{i}{\bf
d}_{i}e^{-i\omega_{i}t}+\hat{\sigma}^{+}_{i}{\bf d}^{*}_{i}e^{i\omega_{i}t},$
$\widehat{\bf E}_{f}({\bf r},t)=i\lambda_{f}{\bf
e}_{f}\left(\hat{a}_{f}e^{-i\omega_{f}t+i{\bf k}_{f}\cdot{\bf
r}}-\hat{a}^{\dagger}_{f}e^{i\omega_{f}t-i{\bf k}_{f}\cdot{\bf r}}\right).$

##  III Equations for emitters and the field-tensor

Equation (7) is the one of our attention and will be transformed into a master
equation. To do this, we first rearrange the terms and write it in the form:

$\frac{d}{dt}\rho_{\\{a\\}}=\frac{i}{\hbar}\sum_{i}\left[\widehat{\bf
d}_{i}(t)\cdot\bm{\mathcal{E}}({\bf r}_{i},t),\rho_{\\{a\\}}(t)\right]\\\
+\frac{i}{\hbar}\sum_{i}\left[\widehat{\bf
d}_{i}(t)\cdot\widehat{\bf\Phi}({\bf r}_{i},t)\right],$ (11)

where the commutator of the vector operators reads as $[\hat{\bf
a}\cdot\hat{\bf b}]=(\hat{\bf a}\cdot\hat{\bf b})-(\hat{\bf b}\cdot\hat{\bf
a})$. The traces in Eq. (7) have produced the quantum mechanical averages with
respect to the field operators which we denote as $\bm{\mathcal{E}}({\bf
r},t)$ and $\widehat{\bf\Phi}({\bf r},t)$:

$\bm{\mathcal{E}}({\bf r},t)=\sum_{f}Tr_{f}\\{\widehat{\bf E}_{f}({\bf
r},t)\rho_{f}(t)\\},$ (12) $\widehat{\bf\Phi}({\bf
r},t)=\sum_{f}Tr_{f}\\{\widehat{\bf E}_{f}({\bf r},t)g_{\\{a\\}f}(t)\\}.$ (13)

Once evaluated at the positions of the emitters ${\bf r}_{i}$ they represent
respectively the action of the effective mean field and weaker contributions
from the emitter-field correlations. Equations (12) and (13) should be
expanded through substitution of the formal solutions of Eqs. (8)-(10). This
will produce the commutators of the the field operators:

$\sum_{f}\left[\widehat{\bf E}_{f}({\bf r},t)\otimes\widehat{\bf E}_{f}({\bf
r}^{\prime},t^{\prime})\right]=-i\hbar\overleftrightarrow{\bf G}({\bf r}-{\bf
r}^{\prime},t-t^{\prime}),$ (14)

where $\overleftrightarrow{\bf G}({\bf r},t)$ is the Green’s tensor of
electromagnetic field propagation in vacuum [44]. Note that the commutation in
this case must be treated as $\big{[}\hat{\bf a}\otimes\hat{\bf
b}\big{]}=\hat{\bf a}\otimes\hat{\bf b}-(\hat{\bf b}\otimes\hat{\bf a})^{T}$
and give a dyadic product.

For the positive and negative frequency components of the electric field
$\widehat{\bf E}_{f}=\widehat{\bf E}_{f}^{+}+\widehat{\bf E}_{f}^{-}$ we use
the following commutation relations:

$\sum_{f}\left[\widehat{\bf E}_{f}^{\pm}({\bf r},t)\otimes\widehat{\bf
E}^{\mp}_{f}({\bf r}_{i},t^{\prime})\right]=-i\hbar\overleftrightarrow{\bf
G}^{\pm}({\bf r}-{\bf r}_{i},t-t^{\prime}),$ (15)

where $\overleftrightarrow{\bf G}^{\pm}$ \- advanced and retarded parts of the
Green’s tensor.

Consequently, for equation (12) we find:

$\bm{\mathcal{E}}({\bf r},t)=\bm{\mathcal{E}}_{0}({\bf
r},t)+\int_{0}^{t}dt^{\prime}\sum_{i}\overleftrightarrow{\bf G}({\bf r}-{\bf
r}_{i},t-t^{\prime}){\bf p}_{i}(t^{\prime}),$ (16)

where $\bm{\mathcal{E}}_{0}({\bf r},t)=\sum_{f}Tr_{f}\\{\widehat{\bf
E}_{f}({\bf r},t)\rho_{f}^{0}\\}$ is determined by the initial value of
$\rho_{f}^{0}=\rho_{f}(t=0)$. This part of the field describes all possible
incident perturbations acting on the emitters as well as noise background. We
may assume a situation when at $t=0$ the quantized field is in the vacuum
state except for one or a few selected modes in coherent states. This
assumption will set up the presence of the incident laser beam. Terms ${\bf
p}_{i}(t)=\langle\widehat{\bf d}_{i}(t)\rangle_{i}$ are namely the induced
dipole moments, so the second part in (16) represents the reaction field.

Using a formal solution of the equation (9) and commutation relations of field
operators (15), we obtain the expression for $\widehat{\bf\Phi}$:

$\widehat{\bf\Phi}({\bf
r},t)=\sum_{i}\int_{0}^{t}dt^{\prime}\bigg{\\{}\overleftrightarrow{\bf
G}^{+}({\bf r}-{\bf r}_{i},t-t^{\prime})\overline{\bf
d}_{i}(t^{\prime})\rho_{\\{a\\}}(t^{\prime})\\\ +\overleftrightarrow{\bf
G}^{-}({\bf r}-{\bf
r}_{i},t-t^{\prime})\rho_{\\{a\\}}(t^{\prime})\overline{\bf
d}_{i}(t^{\prime})\bigg{\\}}\\\
+\frac{i}{\hbar}\sum_{i}\int_{0}^{t}dt^{\prime}\overleftrightarrow{\bf W}({\bf
r},{\bf r}_{i},t,t^{\prime})\left[\widehat{\bf
d}_{i}(t^{\prime}),\rho_{\\{a\\}}(t^{\prime})\right].$ (17)

The first sum in this equation contains retarded Green’s tensors and
$\overline{\bf d}_{i}=\widehat{\bf d}_{i}-{\bf p}_{i}$. The terms in this sum
are evaluated with previous states of the emitters, i.e.,
$\rho_{\\{a\\}}(t^{\prime}),t^{\prime}\leq t$. The second sum contains tensor
we have denoted as $\overleftrightarrow{\bf W}$ and further refer to as the
scattered field correlation tensor:

$\overleftrightarrow{\bf W}({\bf r},{\bf
r}^{\prime},t,t^{\prime})=\sum_{f,f^{\prime}\neq
f}\Big{\\{}\langle:\widehat{\bf E}_{f}({\bf r},t)\otimes\widehat{\bf
E}_{f}({\bf r}^{\prime},t^{\prime}):\rangle_{f}\\\ -\langle\widehat{\bf
E}_{f}({\bf r},t)\rangle_{f}\otimes\langle\widehat{\bf E}_{f}({\bf
r}^{\prime},t^{\prime})\rangle_{f}\\\ +Tr_{ff^{\prime}}\\{\widehat{\bf
E}_{f}({\bf r},t)\otimes\widehat{\bf E}_{f^{\prime}}({\bf
r}^{\prime},t^{\prime})g_{ff^{\prime}}(t^{\prime})\\}\Big{\\}},$ (18)

where : : means the normal ordering of the operators. This tensor describes
the correlation of the components of nonclassical incoherent electromagnetic
field and the dynamic intermode redistribution of photons. Note that for a
field described by a set of coherent states, this tensor will be strictly
equal zero. It must be outlined that this tensor is dependent on the previous
performance of the scatterers (the two emitters in our case). Indeed, in (18)
one can see that $\rho_{f}$ and $g_{ff^{\prime}}$ are linked with
$\rho_{\\{a\\}}$ as it follows from the equations (8)-(10). It means these
terms in the emitter-field correlation provides the feedback loop via the
scattered field.

If the time required for the light signal to transverse the system is small in
the comparison to the time $\Delta t$ required for appreciable changes in
population of the atomic levels, it is possible to use Markov approximation,
i.e., to put $\rho(t^{\prime})~{}\approx~{}\rho(t)$. For the interatomic
distances of the order of the resonant wavelength this condition is satisfied.

Now we need to rewrite the material equation using the above results. Thus,
introducing the formal solution of (8) and (9) into (7) and using Markov
approximation equation (7) transforms into:

$\frac{d}{dt}\rho_{\\{a\\}}=\frac{i}{\hbar}\sum_{i}\left[\widehat{\bf
d}_{i}(t)\bm{\mathcal{E}_{0}}({\bf r}_{i},t),\rho_{\\{a\\}}\right]\\\
+\frac{i}{\hbar}\sum_{i,j}\Bigg{[}\widehat{\bf
d}_{i}(t),\int_{0}^{t}dt^{\prime}\bigg{\\{}\overleftrightarrow{\bf G}^{+}({\bf
r}_{i}-{\bf r}_{j},t-t^{\prime}){\bf d}_{j}(t^{\prime})\rho_{\\{a\\}}\\\
+\overleftrightarrow{\bf G}^{-}({\bf r}_{i}-{\bf
r}_{j},t-t^{\prime})\rho_{\\{a\\}}{\bf d}_{j}(t^{\prime})\bigg{\\}}\Bigg{]}\\\
-\frac{1}{\hbar^{2}}\sum_{i,j}\left[\widehat{\bf
d}_{i}(t),\int_{0}^{t}dt^{\prime}\overleftrightarrow{\bf W}({\bf r}_{i},{\bf
r}_{j},t,t^{\prime})\left[\widehat{\bf
d}_{j}(t^{\prime}),\rho_{\\{a\\}}\right]\right].$ (19)

The first group of commutators in the RHS contains the field
$\bm{\mathcal{E}}_{0}$ from equation (16), which in a particular case may
represent the incident laser. One can immediately recognize the appearance of
the relevant Rabi frequencies, that fulfil the coherent part of the
interaction. Note that the time integral in (16) cancels with the terms with
${\bf p}_{i}$ from (17) under Markov approximation. That means a natural
exclusion of the self-action of the emitters provided by BBGKY.

The second sum collects the commutators which are responsible for various
incoherent interactions (energy exchange) between the emitters and the field.
Because of the properties of the Green’s tensor these terms can be separated
with respect to its real and imaginary parts. Besides, each part will have the
commutations in which the Green’s tensor is taken at the poles (when $i=j$)
and the commutations in which it is evaluated for the distances between the
emitters ($i\neq j$). Thus, one may anticipate these term will describe the
radiative decay and transition frequency shifts respectfully.

The last group of terms contains the field correlation tensor (18). The
operator part suggests that they will produce multiple cross-contributions to
the kinetics of the populations and coherencies in the emitting pair. These
contributions are presumably characterized by the rates determined by the flux
density of incoherent photons. That includes multiple scattering of light and
excitation exchange between the emitters.

##  IV The cooperative master equation

To get the solvable equation for $\rho_{\\{a\\}}$ we use the rotating wave
approximation (RWA) and go back to Schroedinger’s picture into the frame
rotating with the laser frequency. As driving field we consider the plane
monochromatic wave $\bm{\mathcal{E}}_{0}({\bf
r},t)=Re(\bm{\mathcal{E}}e^{i\omega_{L}t-i{\bf k}_{L}\cdot{\bf r}})$ with
complex amplitude $\bm{\mathcal{E}}$, frequency $\omega_{L}$ and wave vector
${\bf k}_{L}$. Thus, we perform the following unitary transformation:

$\widetilde{\widetilde{\rho}}=exp\left\\{i\sum_{i}\Delta_{i}\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{i}t\right\\}\rho\>exp\left\\{-i\sum_{i}\Delta_{i}\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{i}t\right\\}.$

where $\Delta_{i}=\omega_{i}-\omega_{L}$ is the detuning of the laser
frequency from the corresponding transition frequencies of the emitters.

The time integral in the second sum in (19) are computable as the a one-way
Fourier transform of the Green’s function:

$\displaystyle\int_{0}^{t}dt^{\prime}\overleftrightarrow{\bf G}^{\pm}({\bf
r},t-t^{\prime})e^{i\omega(t-t^{\prime})}=\overleftrightarrow{\bf
G}^{\pm}({\bf r},\omega).$

It satisfy the following relations:

$\displaystyle\overleftrightarrow{\bf G}^{-}({\bf
r},\omega)=0,\;\;\;\overleftrightarrow{\bf G}^{+}({\bf r},-\omega)=0,$
$\displaystyle\overleftrightarrow{\bf G}^{+}({\bf
r},\omega)=\overleftrightarrow{\bf G}({\bf r},\omega)$
$\displaystyle\overleftrightarrow{\bf G}^{-}({\bf
r},-\omega)=\overleftrightarrow{\bf G}({\bf
r},-\omega)=\overleftrightarrow{\bf G}^{*}({\bf r},\omega).$

Figure 2: The cooperative relaxation rate (upper) and the parameter of dipole-
dipole interaction (lower) as functions of angle $\theta$ for constant
separation $r_{12}=0.08\lambda$ ($\lambda=\lambda_{1}$).

In considered picture after integration this sum consists of the terms we
denote as $Z_{ij}$, which take the following form:

$Z_{ij}=D_{ij}\left[\hat{\sigma}_{i}^{+},\hat{\sigma}_{j}^{-}\rho_{\\{a\\}}\right]+D_{ij}^{*}\left[\hat{\sigma}_{i}^{-},\rho_{\\{a\\}}\hat{\sigma}_{j}^{+}\right],$

where

$D_{ij}=\frac{1}{\hbar}{\bf d}_{i}^{*}\overleftrightarrow{\bf G}({\bf
r}_{i}-{\bf r}_{j},\omega_{0}){\bf d}_{j}.$

Here, $\omega_{0}=(\omega_{1}+\omega_{2})/2$ is the emitters average
frequency.

The imaginary part of $Z_{ij}$ is responsible for the relaxation of the
density matrix:

$\mathrm{Im}Z_{ij}=\gamma_{ij}\left(\left[\hat{\sigma}_{i}^{+},\hat{\sigma}_{j}^{-}\rho_{\\{a\\}}\right]+\left[\rho_{\\{a\\}}\hat{\sigma}_{i}^{+},\hat{\sigma}_{j}^{-}\right]\right),$

where $\gamma_{ij}=\mathrm{Im}(D_{ij})$. For ${\bf d}_{i}^{*}={\bf d}_{i}$, we
get:

$\gamma_{ij}=\frac{1}{\hbar}{\bf
d}_{i}\mathrm{Im}\left(\overleftrightarrow{\bf G}({\bf r}_{i}-{\bf
r}_{j},\omega_{0})\right){\bf d}_{j}.$

If $i=j$, then $\gamma_{ii}=\gamma_{i}=2|{\bf
d}_{i}|^{2}\omega_{i}^{3}/(3\hbar c^{3})$ is the radiation relaxation rate of
a single emitters; for $i\neq j$, the value $\gamma_{ij}$ determines the
collective damping rate.

The real part of $Z_{ij}$ produces the radiative shifts for $i=j$ which refer
to the transitions of individual emitters while for $i\neq j$ we get the
shifts of the energy levels of the collective ensemble:

$\mathrm{Re}Z_{ij}=\Omega_{ij}\left[\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-},\rho_{\\{a\\}}\right],$

where

$\Omega_{ij}=\mathrm{Re}(D_{ij})=\frac{1}{\hbar}{\bf
d}_{i}\mathrm{Re}\left(\overleftrightarrow{\bf G}({\bf r}_{i}-{\bf
r}_{j},\omega_{0})\right){\bf d}_{j}$

is the individual radiative shift for $i=j$ and the dipole-dipole interaction
parameter for $i\neq j$. In the current work we will neglect the individual
radiative shifts. Taking these terms into account leads to small changes in
the transition frequencies ($\omega_{i}\rightarrow\omega_{i}+\Omega_{ii}$).
Obviously, we may rely on the convention to assume these shifts as included to
$\Delta_{i}$.

To calculate the collective relaxation rate and the dipole – dipole
interaction parameter, we use the explicit form of the Green’s tensor of the
electromagnetic field propagation:

$\displaystyle\overleftrightarrow{\bf G}({\bf
r}_{ij},\omega)=k^{2}\frac{e^{ikr_{ij}}}{r_{ij}}\left(R(ikr_{ij})\overleftrightarrow{\bf
I}+Q(ikr_{ij}){\bf n}_{ij}\otimes{\bf n}_{ij}\right),$ $\displaystyle{\bf
r}_{ij}={\bf r}_{i}-{\bf r}_{j},\;\;\;{r}_{ij}=|{\bf
r}_{ij}|,\;\;\;k=\frac{\omega}{c},\;\;\;{\bf n}_{ij}\otimes{\bf
n}_{ij}=\frac{{\bf r}_{ij}\otimes{\bf r}_{ij}}{{\bf r}_{ij}^{2}},$
$\displaystyle
P(z)=1-\frac{1}{z}+\frac{1}{z^{2}},\;\;\;\;Q(z)=-1+\frac{3}{z}-\frac{3}{z^{2}}.\;\;\;\;\;\;\;\;$

Thus, the relaxation rate and dipole-dipole interaction parameters take the
known form [27, 37]:

$\gamma_{12}=\gamma_{21}=\frac{3}{2}\frac{\sqrt{\gamma_{1}\gamma_{2}}}{k_{0}r_{12}}\Bigg{(}\cos(k_{0}r_{12})\left[\frac{1}{k_{0}r_{12}}-\frac{3\cos^{2}\theta}{k_{0}r_{12}}\right]\\\
+\sin(k_{0}r_{12})\left[1-\cos^{2}\theta+\frac{3\cos^{2}\theta-1}{k_{0}^{2}r_{12}^{2}}\right]\Bigg{)},$

and

$\Omega_{12}=\Omega_{21}=\frac{3}{2}\frac{\sqrt{\gamma_{1}\gamma_{2}}}{k_{0}r_{12}}\Bigg{(}\sin(k_{0}r_{12})\left[\frac{3\cos^{2}\theta}{k_{0}r_{12}}-\frac{1}{k_{0}r_{12}}\right]\\\
+\cos(k_{0}r_{12})\left[1-\cos^{2}\theta+\frac{3\cos^{2}\theta-1}{k_{0}^{2}r_{12}^{2}}\right]\Bigg{)},$

Here, $\theta$ is the angle between the direction of the dipole moments ${\bf
d}_{i}$ and the vector connecting two emitters ${\bf r}_{12}$, and
$k_{0}=(k_{1}+k_{2})/2$. The relaxation rate and dipole-dipole interaction
parameters describe the collective damping and the collective shift of energy
levels and determine the collective properties of the two-atom system. The
interaction strongly depends on interatomic separation $r_{12}$ and the
spatial orientation of dipole moments. Figure 2 illustrates the dependence of
the collective parameters $\gamma_{12}$ and $\Omega_{12}$ (hereinafter all
parameters are scaled to $\gamma_{1}$ and, for simplicity, we take
$\gamma_{1}=1$) on the angle $\theta$, which reflects the incident field
polarization (see Fig.1).

Finally, we obtain the master equation describing the time evolution of the
atomic system:

$\frac{d}{dt}\rho_{\\{a\\}}=-i\sum_{i}\Delta_{i}\left[\hat{\sigma}^{+}_{i}\hat{\sigma}^{-}_{i},\rho_{\\{a\\}}\right]+i\sum_{i}\left[\Omega_{i}({\bf
r}_{i})\hat{\sigma}^{-}_{i}+\Omega^{*}_{i}({\bf
r}_{i})\hat{\sigma}^{+}_{i},\rho_{\\{a\\}}\right]+i\sum_{i\neq
j}\Omega_{ij}\left[\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-},\rho_{\\{a\\}}\right]\\\
-\sum_{i,j}\gamma_{ij}\left(\left[\hat{\sigma}_{i}^{+},\hat{\sigma}_{j}^{-}\rho_{\\{a\\}}\right]+\left[\rho_{\\{a\\}}\hat{\sigma}_{i}^{+},\hat{\sigma}_{j}^{-}\right]\right)-\frac{1}{\hbar^{2}}\sum_{i,j}\Big{(}\left[\hat{\sigma}_{i}^{+},\left[\hat{\sigma}_{j}^{-},\rho_{\\{a\\}}\right]\right]\int_{0}^{t}dt^{\prime}{\bf
d}_{i}\overleftrightarrow{\bf W}({\bf r}_{i},{\bf r}_{j},t,t^{\prime}){\bf
d}_{j}e^{i\omega_{j}(t-t^{\prime})}\\\
+\left[\hat{\sigma}_{i}^{-},\left[\hat{\sigma}_{j}^{+},\rho_{\\{a\\}}\right]\right]\int_{0}^{t}dt^{\prime}{\bf
d}_{i}\overleftrightarrow{\bf W}({\bf r}_{i},{\bf r}_{j},t,t^{\prime}){\bf
d}_{j}e^{-i\omega_{j}(t-t^{\prime})}\Big{)}.$ (20)

Here $\Omega_{i}({\bf r})=\Omega_{i}e^{-i{\bf k}_{L}{\bf\cdot r}}$,
$\Omega_{i}={\bf d}_{i}\cdot\bm{\mathcal{E}}/2\hbar$ is the Rabi frequency.

The last sum in this equation is responsible for the various exchange channels
between the populations and coherences induced by multiple scattering of
photons. This follows from the operator part, i.e., the commutations of the
projectors with density matrix. The time integral originates from the term
commented on under equation (18). In this paper we will assume that all
photoluminescence leaves the emitter system immediately with no chances to
return. In other words, we state the conditions when we can ignore the
multiple scattering. The choice of conditions that would allow approximate
evaluation of time integrals deserves special attention. Solving the master
equation with the feedback loop is a subject of our separate study and will be
reported in another paper.

Thus, we have prepared a new master equation that describes the evolution of
two near-identical two-level emitters driven by an incident field. This master
equation is purely universal in terms of the incident plane wave geometry and
positions of the emitters. For large interatomic separations, i.e., for
$kr_{12}\gg 1$ and $\gamma_{12}=\Omega_{12}=0$ in accordance with their
dependence on $r_{12}$, the cooperative coupling is very weak or can be
assumed zero and the master equation is reduced to the Lindblad master
equations for single atoms.

Note that we obtained the master equation for cooperative system after we
started with the conventional electric dipole Hamiltonian and performed the
other conventional steps, i.e., Born-Markov and the rotating wave
approximations. No by-hand terms were needed as dipole-dipole coupling and
collective relaxation are organized by the BBGKY from the electric dipole
potential for the interaction between the emitters and the quantized field.

##  V Emission intensity and excitation spectra

In this section we are to calculate the total photoluminescence intensity from
the coupled emitters as function of the pump frequency. We will use this
dependence to plot photoluminescence excitation spectra for different
arrangements of the excitation and detection facilities. From the literature
we know that the emission intensity has been conventionally calculated upon
the following reasoning: it is determined from the one-time first-order
correlation function of the electromagnetic field at point ${\bf r}$ in the
far-field zone of the radiation emitted by the atomic system [39]. However we
have at hand a new tool which may help us make a shortcut to the equation for
the total photoluminescence intensity. It was shown in the previous section
that BBGKY gives birth to field correlation tensor $\overleftrightarrow{\bf
W}({\bf r},{\bf r}^{\prime},t,t^{\prime})$. Being a natural internal part of
our system it may have different physical interpretations along with different
values of its arguments. If taken at one point $({\bf r},t)$, where ${\bf r}$
is the position of the ”observer” and $t$ is the time of registration, the
tensor would determine the amount of photoluminescence signal captured by the
detector. Thus, the electromagnetic field correlation tensor defined in (18)
taken at the point $({\bf r},t)$ allows us to calculate the total radiation
intensity. One may combine the definition in equation (18) and the Poynting
vector to show that ${\bf I}({\bf r},t)\sim(c/4\pi)\overleftrightarrow{\bf
W}({\bf r},{\bf r},t,t)$. Equation 18 can be evaluated if we use the formal
solutions for $\rho_{f}$ (8) and $g_{ff^{\prime}}$ (10) from the equations 8
and 10 respectfully. By completing two more steps, which are the summation
over the photonic modes and applying the Markov approximation, i.e.
$\rho_{\\{a\\}}(t^{\prime})=\rho_{\\{a\\}}(t)$, we come to the expression for
the two-point tensor:

$\overleftrightarrow{\bf W}({\bf r},{\bf r}^{\prime},t,t^{\prime})\\\
=\int_{0}^{t^{\prime}}d\tau
Tr_{\\{a\\}}\Bigg{\\{}\sum_{i}\overleftrightarrow{\bf G}({\bf r}-{\bf
r}_{i},t-\tau)\widehat{\bf d}_{i}(\tau)\otimes\widehat{\bf\Phi}({\bf
r}^{\prime},t^{\prime})\\\ +\widehat{\bf\Phi}({\bf
r},t)\otimes\sum_{i}\widehat{\bf d}_{i}(\tau)\overleftrightarrow{\bf G}({\bf
r}^{\prime}-{\bf r}_{i},t^{\prime}-\tau)\Bigg{\\}}.$

Here one can see that the solution for $\overleftrightarrow{\bf W}({\bf
r},{\bf r}^{\prime},t,t^{\prime})$ produces $\widehat{\bf\Phi}({\bf r},t)$,
which contains $\overleftrightarrow{\bf W}({\bf r},{\bf
r}^{\prime},t,t^{\prime})$ itself according to (17). It means that certain
contributions will be reoccurring on a loop principle producing higher order
terms. On the other hand, our current objective is to extract and study pure
photoluminesce process. By looking again at (17) and (7) one can clearly see
that there is only a single term in $\widehat{\bf\Phi}({\bf r},t)$ that
produces the radiative relaxation operator in the master equation (7).
Consequently, it is sufficient to keep only this part of the complete
$\widehat{\bf\Phi}({\bf r},t)$ when substituting into the tensor
$\overleftrightarrow{\bf W}({\bf r},{\bf r}^{\prime},t,t^{\prime})$ to
describe spontaneous emission:

$\widehat{\bf\Phi}_{\mathcal{L}}({\bf
r},t)=\sum_{j}\int_{0}^{t}dt^{\prime}\bigg{\\{}\overleftrightarrow{\bf
G}^{+}({\bf r}-{\bf r}_{j},t-t^{\prime})\widehat{\bf
d}_{j}(t^{\prime})\rho_{\\{a\\}}\\\ +\overleftrightarrow{\bf G}^{-}({\bf
r}-{\bf r}_{j},t-t^{\prime})\rho_{\\{a\\}}\widehat{\bf
d}_{j}(t^{\prime})\bigg{\\}},$

Substitution of the solution for the function
$\hat{\bf\Phi}_{\mathcal{L}}({\bf r},t)$ and taking $r^{\prime}\to
r,t^{\prime}\to t$ must be followed with applying the rotating wave
approximation. This brings us to a new tensor that we refer to as intensity
tensor $\overleftrightarrow{\bf I}({\bf r},t)$:

$\overleftrightarrow{\bf I}({\bf r},t)=\frac{c}{4\pi}\overleftrightarrow{\bf
W}({\bf r},{\bf
r},t,t)=\sum_{i,j}\langle\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-}\rangle\overleftrightarrow{\bf
I}_{ij},$

which is represented in the form of components $\overleftrightarrow{\bf
I}_{ij}$ with factors that are the corresponding elements of the atomic
density matrix $\langle\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-}\rangle\equiv
Tr_{\\{a\\}}\\{\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-}\rho_{\\{a\\}}\\}$. The
tensor components can now be written in the form where the corresponding time
integrations have been performed using the same steps as in previous section.
It produces the Fourier transforms of the Green’s tensors:

$\overleftrightarrow{\bf I}_{ij}=\frac{c}{\pi}{\bf
d}_{i}^{*}\overleftrightarrow{\bf G}^{*}({\bf r}-{\bf
r}_{i},\omega_{i})\otimes\overleftrightarrow{\bf G}({\bf r}-{\bf
r}_{j},\omega_{j}){\bf d}_{j}.$ (21)

Here, the spatial arguments describe the distances from each emitter to the
observer. It is beneficial for further calculations to modify the Green’s
tensor as follows:

$\overleftrightarrow{\bf G}({\bf r-r}_{i},\omega_{i})=\overleftrightarrow{\bf
G}({\bf r},\omega_{i})e^{-ik_{i}r_{ij}({\bf n}\cdot{\bf n}_{ij})/2},$

where $r_{ij}=|{\bf r}_{i}-{\bf r}_{j}|=r_{ji}$, ${\bf n}_{ij}={\bf
r}_{ij}/r_{ij}=-{\bf n}_{ji}$ (note that ${\bf n}_{ii}=0$),
$k_{i}=\omega_{i}/c$. Equation (21) is general as it allows a solution for
arbitrary orientations of the emitting dipoles. However, for the problem
stated in section IV we may assume that the dipole moments are collinear and
represent the vector of the dipole moment in the form ${\bf d}_{i}={d}_{i}{\bf
n}_{d}$. Thus, we get:

$\overleftrightarrow{\bf I}({\bf
r},t)=\sum_{i,j}\langle\hat{\sigma}_{i}^{+}\hat{\sigma}_{j}^{-}\rangle\overleftrightarrow{\bf
I}_{ij}e^{i\phi_{ij}+i(k_{i}-k_{j})r},$ (22)

where $\phi_{ij}=r_{ij}k_{0}({\bf n}\cdot{\bf
n}_{ij})=r_{ij}k_{0}\cos{\varphi}$ is the observation parameter associated
with the angle between vectors ${\bf n}_{ij}$ and ${\bf n}={\bf r}/r$ and

$\displaystyle\overleftrightarrow{\bf I}_{ij}=$
$\displaystyle\frac{\omega_{i}^{2}\omega_{j}^{2}{d}_{i}{d}_{j}}{\pi
r^{2}c^{3}}\Big{(}({\bf n}_{d}\otimes{\bf n}_{d})$ $\displaystyle-({\bf
n}_{d}\otimes{\bf n}+{\bf n}\otimes{\bf n}_{d})({\bf n}_{d}\cdot{\bf n})+({\bf
n}\otimes{\bf n})({\bf n}_{d}\cdot{\bf n})^{2}\Big{)}.$

In order to get a scalar expression for the intensity in the desired direction
with respect of the orientation of the emitting dipole moments it is necessary
to evaluate the values ${I}_{ij}={\bf n}_{x_{1}}\overleftrightarrow{\bf
I}_{ij}{\bf n}_{x_{1}}+{\bf n}_{y_{1}}\overleftrightarrow{\bf I}_{ij}{\bf
n}_{y_{1}}$. Here vectors $n_{x_{1}}$ and $n_{y_{1}}$ are perpendicular to the
vector directed at the observer and orthogonal to each other. Thus, we get a
scalar value:

$\displaystyle{I}_{ij}=\frac{\omega_{i}^{2}\omega_{j}^{2}{d}_{i}{d}_{j}}{\pi
r^{2}c^{3}}\sin^{2}{\beta},$

where we introduce $\beta$ as the angle between observation direction ${\bf
r}$ and the direction of the dipole moments ${\bf d}_{i}$.

For a final evaluation of the total intensity we need the solutions for
density matrix elements. The equations for these follow from the master
equation (20). The master equation gives 15 kinetic equations describing the
evolution of the atomic variables (components of density matrix). Our goal is
to find steady-state intensity so we put $\dot{\rho}_{\\{a\\}}=0$. It is
favourable to rewrite the elements of the atomic density matrix using the
collective states representation [1]:

$\displaystyle|G\rangle=|g_{1}g_{2}\rangle$
$\displaystyle|S\rangle=1/\sqrt{2}(|e_{1}g_{2}\rangle+|g_{1}e_{2}\rangle)$
$\displaystyle|A\rangle=1/\sqrt{2}(|e_{1}g_{2}\rangle-|g_{1}e_{2}\rangle)$
$\displaystyle|U\rangle=|e_{1}e_{2}\rangle$

In the collective states basis, the two-atom system behaves as a single four-
level system with the ground state $|G\rangle$, the upper state $|U\rangle$,
and two intermediate states: the symmetric state $|S\rangle$ and the
antisymmetric state $|A\rangle$. The energies of the intermediate states
depend on the dipole-dipole interaction.

In this way, the total emission intensity from a two-particle cooperative
system is written as

${\bf I}({\bf
r})=\frac{I_{11}+I_{22}}{2}(2\rho_{UU}+\rho_{SS}+\rho_{AA})-2I_{12}\sin{\alpha}\mathrm{Im}\rho_{SA}\\\
+I_{12}\cos{\alpha}(\rho_{SS}-\rho_{AA})+(I_{11}-I_{22})\mathrm{Re}\rho_{SA},$
(23)

where $\alpha=ir_{ij}k_{0}\cos{\varphi}+i(k_{1}-k_{2})r$.

Figure 3: PLE spectra of two identical atoms (a) with $\Delta_{\omega}=0$
$d_{1}=d_{2}$, and non-identical atoms (b) and (c) with $\Delta_{\omega}=20$,
$d_{2}=1.5d_{1}$, for different distances between them. Other parameters:
$\Omega_{1}=2$, $\Omega_{2}=\Omega_{1}d_{2}/d_{1}$, $\xi=\pi$, $\theta=\pi/2$,
$\varphi=0$.

Our objective is to investigate the photoluminescence excitation (PLE)
spectra, i.e., the steady-state emission intensity as a function of the
detuning $\Delta=(\omega_{0}-\omega_{L})$, where
$\omega_{0}=(\omega_{1}+\omega_{2})/2$. Registration of PLE spectra is a
common experimental technique implemented with tunable lasers. First, we will
compare the entanglements occurring between two identical emitters with the
entanglements possible for two near-identical emitters, i.e., with individual
transition frequencies and moments. Since entanglement is enforced by the
dipole-dipole coupling and collective decay it is of special interest to see
how the PLE spectra are changing with the distance between the emitters. The
incident beam parameters are also important in this regard. We will use the
concept of observing photoluminescence from emitters with fixed positions when
excitation may be subject to choice. Figures 3(a)-(c) show the PLE spectra
depending on the separation between the particles (in terms of resonance
wavelength of one of the emitters $\lambda_{1}=\lambda$). The rear spectral
profiles in Figs.3(a)-(b) demonstrate weak entanglement when the separation is
large enough to let the emitters act independently ($\Omega_{12}$ becomes
small). For identical atoms in Fig.3(a) one can see a single aggregated peak
corresponding to the resonance frequency of both emitters. For non-identical
atoms in Fig.3(b) one finds two distinct peaks at the frequencies
corresponding to each resonance. Figure 3(c) illustrates the situation when
the emitters are at the distances when dipole-dipole coupling and collective
decay are already negligible small. In this case one can clearly see the two
peaks which are largely individual PLEs. However when the excitation is tuned
between the transitions we see the fluctuations of the resulting interference
between individual emissions. This occurs due to non-vanishing matrix elements
in (23). The aparency of this interference is determined by angles $\xi$ and
$\varphi$.

For smaller distances the spectra profiles become notably different (the front
stage profiles in Figs.3(a)-(b)). Convergence of the atoms split the picture
into three peaks. The three peaks demonstrate the well-known configuration of
the cooperative resonances, observed experimentally in [8]. As was shown, the
central peak corresponds to the simultaneous absorption of two photons and
subsequent radiation of the pair of correlated photons. It is evidenced by
photon-bunching in the autocorrelation function, which was also measured in
the experiment. The side peaks correspond to the collective states $|S\rangle$
and $|A\rangle$, where $|S\rangle$ is superrariative and $|A\rangle$ is
subradiative. The only so-far successfully observation of the cooperative
triplet [8] indeed has brought the evidence of the physical picture we
simulate in out study. They measured the PLE spectra from two terrylene
molecules embedded in a para-terphenyl crystal. Terrylene molecules are very
well studied and known to demonstrate wide inhomogeneous broadening. The
emitters in [8] were claimed to have different transition frequencies and be
separated by 12 nanometers. Besides, the conditions of the experiment were
such that the authors suggested the two-level atom model was a good
approximation for terrylene with regard to dipole-dipole coupling. Figure 4
shows how the spectra changes under different excitation intensities. It is
seen that the central peak corresponding to two-photon process arises at
certain incident power and increases with the increasing intensity. This
behaviour is in a good qualitative agreement with the PLE plots in [8]. The
spectral positions of the cooperative resonances, the ratios between the peak
maximums and the widths are generally the same in both the experiment and the
theory. However we have prepared the theory to show in detail how the PLE
shape is dependent on the geometry of the problem. Moreover, for a
quantitative comparison with the measurements it may be also necessary to take
care of some special issues like the following: 1) the choice of the
appropriate approximation for the emitters; 2) the role of the conic
collection of photoluminescence photons; and 3) possible effects of the data
processing.

Figure 4: PLE spectra for different excitation strengths with
$\Delta_{\omega}=100$, $d_{1}=1.5d_{2}$, $\xi=\pi/2$, $\varphi=\pi/2$,
$r_{12}=0.04\lambda$.  Figure 5: PLE spectra for identical atoms (solid line),
and non-identical atoms with $\Delta_{\omega}=10$, $d_{1}=d_{2}$ (dashed line)
and $\Delta_{\omega}=20$, $d_{1}=d_{2}$ (dash-dotted line). Other parameters:
$\Omega_{1}=\Omega_{2}=2$, $\xi=\pi/2$, $\theta=\pi/2$, $\varphi=\pi/2$,
$r_{12}=0.08\lambda$.

Now let us look at how geometry of the spatial arrangement of the emitters
against the excitation and detection modify the PLE profiles. First, we
consider the case when the excitation is perpendicular to the interatomic
axis, i.e., ${\bf k}_{L}\perp{\bf r}_{12}$. It is interesting in terms of
comparison with the results known in the literature for two identical
emitters. That variant ($\omega_{1}=\omega_{2}$, ${\bf d}_{1}={\bf d}_{2}$)
was shown to demonstrate photoluminesce intensity with only two distinct peaks
at $\Delta=0$ and $\Delta=-\Omega_{12}$ (solid line in Fig.5). It means that
at this geometry the peak corresponding to the antisymmetric state $|A\rangle$
is not visible. At the same time, for different atoms the spectra exhibits
three peaks at $\Delta=0$ and
$\Delta\simeq\pm\sqrt{\Omega_{12}^{2}+\Delta_{\omega}^{2}}$,
$\Delta_{\omega}=(\omega_{1}-\omega_{2})/2$ (dashed and dash-dotted lines in
Fig.5).

The distinctions between the emitters may as well be vividly demonstrated at
another variant of geometry. The key feature is to observe photolimescence
from different angles (Fig.6). For ${\bf r}\perp{\bf r}_{12}$
($\varphi=\pi/2$) the fluorescence intensity for identical emitters exhibits
two peaks centered at $\Delta=0$ and $\Delta=-\Omega_{12}$ and a dip at
$\Delta=\Omega_{12}$ (solid line in Fig. 6). For $\varphi<\pi/2$ this dip
changes to a peak (dashed line in Fig. 6). The situation is different when
atoms are non-identical. The fluorescence intensity exhibits three peaks for
all observation parameters with small changes in peak height(inset in Fig. 6).

Figure 6: PLE spectra for identical atoms for different directions of
observation: $\varphi=\pi/2$ (solid line), $\varphi=0$ (dashed line). The
inset show the same picture for non-identical atoms with $\Delta_{\omega}=10$,
$d_{1}=d_{2}$. Other parameters: $\Omega_{1}=\Omega_{2}=2$, $\xi=\pi$,
$\theta=\pi/2$, $r_{12}=0.09\lambda$. Figure 7: Migration of the PLE spectra
upon rotation of the polarization of an external field for
$\Delta_{\omega}=10$, $d_{1}=d_{2}$, $\xi=\pi/2$, $\varphi=\pi/2$,
$r_{12}=0.08\lambda$, $\Omega_{1}=\Omega_{2}=2$.  Figure 8: Migration of the
PLE spectra upon rotation of the polarization of an external field for
$\Delta_{\omega}=10$, $d_{1}=d_{2}$, $\xi=\pi/2$, $\varphi=0$,
$r_{12}=0.08\lambda$, $\Omega_{1}=\Omega_{2}=2$.

From the analysis above it appears that two different emitters may happen to
be entangled in the way that provides better conditions for efficient
realisation of the subradiative state. Experimentally one may hope to observe
the asymmetric state emission because it has a greater spectral separation
from the other cooperative resonances and notably stronger intensity.

Figure 9: Position of the maximum fluorescence peak relative to the central
peak as a function of polarization ($\theta$) for $\Delta_{\omega}=0$ (solid
line), $\Delta_{\omega}=10$ (dashed line), $\Delta_{\omega}=25$ (dotted line),
$\Delta_{\omega}=40$ (dashed-dotted line) for $r_{12}=0.08\lambda$.

Now we focus on the particular parameters of the incident field. It is clear
that dipole-dipole coupling parameter is a strong function of angle $\theta$
as it is shown in Fig.2 to change from zero to considerable absolute values.
Therefore is can be straightforwardly anticipated that the form of the PLE
spectra strongly depends on the incident field polarisation, i.e., on the
angle $\theta$ (see Fig.1 and the corresponding comments in section II.2). The
actual relation between the configuration of the beam polariser and $\theta$
is to be calculated for a particular experimental setup. For an illustration
of the effect it is sufficient to choose the incidence such that $\theta$ is
the actual rotation angle of ${\bf E}_{L}$ around ${\bf k}_{L}$. The
transformation of PLE pattern with polarization rotation were calculated for a
fixed distance between emitters as shown in Fig.7. This Figure shows how the
spectral distance between the peaks appears to be shrinking with growing
$\theta$ until the pattern becomes a two-peak PLE. That is because at certain
angle $\theta$ we have $\Omega_{12}=0$ and the distance between peaks is
determined by the difference in the transition frequencies of the emitters.
With further increase of $\theta$ the PLE is again a cooperative triplet in
which the symmetric and antisymmetric peaks have exchanged their positions
against the origin. From the quarter-turn of $\theta$ the PLE pattern starts
to unwind back to the initial picture. Obviously, with a half-turn the
positions of the peaks restore completely. Figure 8 shows similar behaviour as
in Fig.7 but with finer discreteness., while the observation angle was chosen
to be such that no emission is observable at both margins, i.e., $\sin\beta=0$
at $\theta=0,\pi$. In this case one would see how a half-turn of the
excitation polarization makes a photoluminescence signal grow from null, reach
maxima and disappear. For each excitation frequency the emission intensity
will change in some particular manner. In other words, it may be an individual
function of $\theta$. The general picture appears as a gorge formed by the
symmetric and asymmetric summits. The two narrow passes characterize the
differences of the transition frequencies $\Delta_{\omega}$, while the width
of the gorge at the quarter-turn, $\theta=\pi/2$, indicates the distance
between the emitters. Note that the small central chine is peculiar for weak
or moderate excitation powers. Another interpretation of the changed brought
by turning the polarization of the excitation beam is the ”spectral jump” of
the symmetric (or asymmetric) maximum. Figure 9 demonstrates the spectral
migration of the corresponding maximum. For identical emitters this migration
produces the indissoluble trail. On the contrary, the emitters with different
transition frequencies produce broken trails shown in Fig.9, where the gaps
are regulated by $\Delta_{\omega}$. These properties suggest that measuring
the amplitude of spectral migration and the length of spectral jumps can help
detection and characterization of paired emitters.

##  VI Conclusion

In conclusion, we have studied the steady-state photoluminescence from two
near-identical emitters with dipole-dipole coupling. We considered two two-
level particles with individual transition frequencies and transition moments.
The excitation of the system was assumed as an incident continuous
monochromatic beam of light with linear polarization. The geometrical
arrangements of the excitation beam, positions of the emitters and the
detector were included into the solutions. The calculations were performed
using the formalism based on the BBGKY hierarchy of equations, adapted to
solving light-matter interactions and spectroscopic problems. The von Neumann
equation for the density matrix of a many-particle system was written with the
conventional electric-dipole Hamiltonian for two emitters and photons. The
excitation signal was set via the initial conditions for the photonic modes.
The full von Neumann equation was expanded into the BBGKY chain which was
later truncated into a closed system of equations for reduced density matrices
and correlation operators. The system of equations described the two-particle
emitter, the photonic modes, and the emitter-field and intrafield
correlations. It was treated with the use of the conventional approximations
to give the cooperative master equation and the equation for the total
photoluminescence intensity. The cooperative entanglement and radiative decay,
as well as the dipole-dipole interaction appeared naturally in the master
equation as the result of the decoupling procedures applied to the BBGKY.

The master equation and the intensity equation were used to perform computer
simulations of PLE spectra, i.e., the total emission against the excitation
frequency. We showed the properties of the spectral shapes which can reveal
the differences between the emitters and demonstrate the importance of the
excitation-detection geometry relative to the emitting pair. The spectra were
proved to fully agree with the previous thoeries and the experiment. Two non-
identical emitters exhibit a pronounced cooperative PLE triplet in which the
central peak corresponds to the excitation equally detuned from both emitters
transitions while two side peaks correspond accordingly to the super- and sub-
radiative states of the pair or the cooperative resonances. At the same time
we have found new properties of photoluminescence from the near-identical
emitters which are the following:

1\. Non-identical emitters versus identical emitters need smaller separation
to start showing a pronounced cooperative PLE triplet. However, the non-
identical emitters are more favourable for experimental investigations because
the spectral distance between the peaks in PLE is significantly greater while
the antisymmetric peak (sub-radiative state) is substantially inflated.
Besides, there are several configurations of the excitation-detection geometry
for which the symmetric and antisymmetric peaks become comparable in both the
intensity and the width.

2\. At fixed incidence of the excitation beam and fixed registration direction
the rotation of the excitation beam polarisation may become a handy tool in
either nanoscopy or recovering of the diversity of the emitters. We showed
that there is spectral migration of the brightest peak in PLE that follows the
changing polarisation. This can be a good indication of the presence and the
strength of the dipole-dipole interaction as the migration occurs within a
certain range. Moreover, we found the existence of the spectral pseudo-jump of
the brightest peak which may serve as the measure of the difference in
transition frequencies. In this case a PLE trail would have a point
corresponding to a certain polarisation at which the dipole-dipole interaction
vanishes and the spectral distance between symmetric and antisymmetric
resonances is determined solely by the difference of the transition
frequencies.

Note that our results were obtained with the least possible order of the BBGKY
decoupling approximation. The use of the higher order terms provided in the
BBGKY would surely reveal fine details of the cooperative photoluminescence.
In this paper we have paved the direction to further important improvements of
the theory. A strong contribution may be expected from the field correlation
tensor which is present in the full master equation. This contribution is
generally non-vanishing because it appears in the same order of the BBGKY
decoupling. It provides additional terms in kinetics of populations and
coherences due to multiple scattering of photons. It may be of special
importance for particular conditions which provide confinement or trapping of
the emitting pair. Our findings may be important for heading towards all-
optical nanoscopy and construction of entangled systems to be used in quantum
technologies.

###### Acknowledgements.

This work was carried out among the State Contract of the Moscow Pedagogical
State University (MPGU) “Physics of the perspective materials and
nanostructures: basic researches and applications in material sciences,
nanotechnologies and photonics” supported by the Ministry of Education of the
Russian Federation (AAAA-A20-120061890084-9) in collaboration with the Centre
of collective usage “Structural diagnostics of materials” of the Federal
Research Center RAS “Crystallography and photonics”.

## References

  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. 93, 99 (1954).
  * Kuraptsev and Sokolov [2019] A. S. Kuraptsev and I. M. Sokolov, Many-body cooperative effects in an ensemble of pointlike impurity centers near a charged conductive surface, Phys. Rev. A 100, 063836 (2019).
  * Kuś and Wódkiewicz [1981] M. Kuś and K. Wódkiewicz, Two-atom resonance fluorescence, Phys. Rev. A 23, 853 (1981).
  * Agarwal  _et al._ [1980] G. S. Agarwal, R. Saxena, L. M. Narducci, D. H. Feng, and R. Gilmore, Analytical solution for the spectrum of resonance fluorescence of a cooperative system of two atoms and the existence of additional sidebands, Phys. Rev. A 21, 257 (1980).
  * Richter [1983] T. Richter, Cooperative resonance fluorescence from two atoms experiencing different driving fields, Optica Acta: International Journal of Optics 30, 1769 (1983), https://doi.org/10.1080/713821122 .
  * Ficek and Sanders [1990] Z. Ficek and B. C. Sanders, Two-atom resonance fluorescence spectrum in a squeezed vacuum including the dipole-dipole interaction, Quantum Optics: Journal of the European Optical Society Part B 2, 269 (1990).
  * Rudolph  _et al._ [1995] T. G. Rudolph, Z. Ficek, and B. J. Dalton, Two-atom resonance fluorescence in running- and standing-wave laser fields, Phys. Rev. A 52, 636 (1995).
  * Hettich  _et al._ [2002] C. Hettich, C. Schmitt, J. Zitzmann, S. Kühn, I. Gerhardt, and V. Sandoghdar, Nanometer resolution and coherent optical dipole coupling of two individual molecules, Science 298, 385 (2002).
  * Streltsov  _et al._ [2017] A. Streltsov, G. Adesso, and M. B. Plenio, Colloquium: Quantum coherence as a resource, Rev. Mod. Phys. 89, 041003 (2017).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * Pastori  _et al._ [2019] L. Pastori, M. Heyl, and J. C. Budich, Disentangling sources of quantum entanglement in quench dynamics, Phys. Rev. Research 1, 012007 (2019).
  * Kornovan  _et al._ [2019] D. F. Kornovan, N. V. Corzo, J. Laurat, and A. S. Sheremet, Extremely subradiant states in a periodic one-dimensional atomic array, Phys. Rev. A 100, 063832 (2019).
  * Moerner and Kador [1989] W. E. Moerner and L. Kador, Optical detection and spectroscopy of single molecules in a solid, Phys. Rev. Lett. 62, 2535 (1989).
  * Orrit and Bernard [1990] M. Orrit and J. Bernard, Single pentacene molecules detected by fluorescence excitation in a p-terphenyl crystal, Phys. Rev. Lett. 65, 2716 (1990).
  * Basché  _et al._ [1992] T. Basché, W. E. Moerner, M. Orrit, and H. Talon, Photon antibunching in the fluorescence of a single dye molecule trapped in a solid, Phys. Rev. Lett. 69, 1516 (1992).
  * Piliarik and Sandoghdar [2014] M. Piliarik and V. Sandoghdar, Direct optical sensing of single unlabelled proteins and super-resolution imaging of their binding sites, Nature Communications 5, 4495 (2014).
  * Nirmal  _et al._ [1996] M. Nirmal, B. O. Dabbousi, M. G. Bawendi, J. J. Macklin, J. K. Trautman, T. D. Harris, and L. E. Brus, Fluorescence intermittency in single cadmium selenide nanocrystals, Nature 383, 802 (1996).
  * Eichhammer  _et al._ [2015] E. Eichhammer, T. Utikal, S. Götzinger, and V. Sandoghdar, Spectroscopic detection of single pr3+ ions on the 3h4-1d2 transition, New Journal of Physics 17, 083018 (2015).
  * van Oijen  _et al._ [1999] A. M. van Oijen, M. Ketelaars, J. Köhler, T. J. Aartsma, and J. Schmidt, Unraveling the electronic structure of individual photosynthetic pigment-protein complexes, Science 285, 400 (1999).
  * Naumov [2013] A. V. Naumov, Low-temperature spectroscopy of organic molecules in solid matrices: from the shpol’skii effect to laser luminescent spectromicroscopy for all effectively emitting single molecules, Phys. Usp. 183, 633 (2013).
  * Eremchev  _et al._ [2018] I. Y. Eremchev, N. A. Lozing, A. A. Baev, A. O. Tarasevich, M. G. Gladush, A. A. Rozhentsov, and A. V. Naumov, Luminescence microscopy of single quantum dot pairs with nanometer spatial resolution, JETP Letters 108, 30 (2018).
  * Sauter  _et al._ [1986] T. Sauter, R. Blatt, W. Neuhauser, and P. Toschek, ‘quantum jumps’ observed in the fluorescence of a single ion, Optics Communications 60, 287 (1986).
  * Brewer  _et al._ [1992] R. G. Brewer, R. G. DeVoe, and R. Kallenbach, Planar ion microtraps, Phys. Rev. A 46, R6781 (1992).
  * van Loo  _et al._ [2013] A. F. van Loo, A. Fedorov, K. Lalumière, B. C. Sanders, A. Blais, and A. Wallraff, Photon-mediated interactions between distant artificial atoms, Science 342, 1494 (2013).
  * Ropp  _et al._ [2010] C. Ropp, R. Probst, Z. Cummins, R. Kumar, A. J. Berglund, S. R. Raghavan, E. Waks, and B. Shapiro, Manipulating Quantum Dots to Nanometer Precision by Control of Flow, Nano Letters 10, 2525 (2010).
  * Ficek and Tanaś [2002] Z. Ficek and R. Tanaś, Entangled states and collective nonclassical effects in two-atom systems, Physics Reports 372, 369 (2002).
  * Milonn`i and Knight [1975] P. W. Milonn`i and P. L. Knight, Retarded interaction of two nonidentical atoms, Phys. Rev. A 11, 1090 (1975).
  * Ficek  _et al._ [1987] Z. Ficek, R. Tanaś, and S. Kielich, Quantum beats and superradiant effects in the spontaneous emission from two nonidentical atoms, Physica A: Statistical Mechanics and its Applications 146, 452 (1987).
  * Mokhov and Makarov [2019] A. I. Mokhov and A. A. Makarov, Quantum beat spectroscopy in a system of n identical, closely spaced atoms. i. n = 2, Optics and Spectroscopy 127, 7 (2019).
  * Redchenko and Yudson [2014] E. S. Redchenko and V. I. Yudson, Decay of metastable excited states of two qubits in a waveguide, Phys. Rev. A 90, 063829 (2014).
  * Chang  _et al._ [2006] J.-T. Chang, J. Evers, M. O. Scully, and M. Suhail Zubairy, Measurement of the separation between atoms beyond diffraction limit, Phys. Rev. A 73, 031803 (2006).
  * Evers and Qamar [2011] J. Evers and S. Qamar, Subwavelength position measurements with running-wave driving fields, Phys. Rev. A 84, 023802 (2011).
  * Redchenko  _et al._ [2018] E. S. Redchenko, A. A. Makarov, and V. I. Yudson, Nanoscopy of pairs of atoms by fluorescence in a magnetic field, Phys. Rev. A 97, 043812 (2018).
  * Sun  _et al._ [2011] Q. Sun, M. Al-Amri, M. O. Scully, and M. S. Zubairy, Subwavelength optical microscopy in the far field, Phys. Rev. A 83, 063818 (2011).
  * Das  _et al._ [2008] S. Das, G. S. Agarwal, and M. O. Scully, Quantum interferences in cooperative dicke emission from spatial variation of the laser phase, Phys. Rev. Lett. 101, 153601 (2008).
  * Sete and Das [2011] E. A. Sete and S. Das, Quantum interference in timed dicke basis and its effect on bipartite entanglement, Phys. Rev. A 83, 042301 (2011).
  * Akram  _et al._ [2000] U. Akram, Z. Ficek, and S. Swain, Decoherence and coherent population transfer between two coupled systems, Phys. Rev. A 62, 013413 (2000).
  * Lehmberg [1970] R. H. Lehmberg, Radiation from an $n$-atom system. ii. spontaneous emission from a pair of atoms, Phys. Rev. A 2, 889 (1970).
  * Agarwal [1974] G. S. Agarwal, Quantum statistical theories of spontaneous emission and their relation to other approaches, in  _Quantum Optics_ , edited by G. Höhler (Springer Berlin Heidelberg, Berlin, Heidelberg, 1974) pp. 1–128.
  * Allen and Eberly [1978] L. Allen and J. Eberly, Optical Resonance and Two-Level Atoms (1978) p. 222.
  * Richter [1982] T. Richter, Power broadening and two-photon correlations in two-atom resonance fluorescence, Optica Acta: International Journal of Optics 29, 265 (1982).
  * Naumov  _et al._ [2018] A. V. Naumov, A. A. Gorshelev, M. G. Gladush, T. A. Anikushina, A. V. Golovanova, J. Köhler, and L. Kador, Micro-refractometry and local-field mapping with single molecules, Nano Letters 18, 6129 (2018).
  * Bonitz [2016] M. Bonitz, Quantum kinetic theory (Springer International Publishing, 2016) 2nd ed.
  * Mandel and Wolf [1995] L. Mandel and E. Wolf, Optical Coherence and Quantum Optics, in  _American Journal of Physics_ (1995) p. 1192.
  * Eschner  _et al._ [2001] J. Eschner, C. Raab, F. Schmidt-Kaler, and R. Blatt, Light interference from single atoms and their mirror images, Nature 413, 495—498 (2001).
  * Kuznetsov  _et al._ [2011] D. V. Kuznetsov, V. K. Roerich, and M. G. Gladush, Using bbgky hierarchies to study the effect of the local field on the rate of radiative relaxation of quantum systems in a dielectric medium, Theoretical and Mathematical Physics 168, 1078 (2011).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, 1997).
  * Gladush  _et al._ [2011] M. G. Gladush, D. V. Kuznetsov, and V. K. Roerich, Emission spectra and intrinsic optical bistability in a two-level medium, The European Physical Journal D 64, 511 (2011).
  * Rephaeli  _et al._ [2011] E. Rephaeli, i. m. c. E. Kocabaş, and S. Fan, Few-photon transport in a waveguide coupled to a pair of colocated two-level atoms, Phys. Rev. A 84, 063832 (2011).
  * Fang  _et al._ [2014] Y.-L. L. Fang, H. Zheng, and H. U. Baranger, One-dimensional waveguide coupled to multiple qubits: photon-photon correlations, EPJ Quantum Technology 1, 3 (2014).
  * Roerich  _et al._ [2015] V. C. Roerich, A. N. Starostin, and A. A. Panteleev, Many-body approaches to time-dependent resonance fluorescence, Progress in Nonequilibrium Green’s Functions II , 103 (2015).
  * Trubilko and Basharov [2019] A. I. Trubilko and A. M. Basharov, On the Decay of an “Isolated” Oscillator Nonlinearly Coupled to a Damped Oscillator, JETP Letters 110, 517 (2019).