Quantum Theory of a Single Photon in an Arbitrary Medium

1. Introduction

Until now, the quantum theory of single photons has not been sufficiently developed, moreover, the introduction of a wave function for it is often considered a controversial concept, see, for example, [1]. Recall that for the photon there is only the quantum theory of the electromagnetic radiation, which, except for the frequency or wavelength, does not provide any other information about the spacetime structure and other properties of this massless particle [2]. The main reason for the difficulties is that photons are never non-relativistic and can be freely emitted and absorbed. The latter circumstance violates the law of conservation of the number of photons, making their description within the framework of the traditional paradigm of quantum mechanics, if not impossible, then very problematic.

Meanwhile, theoretical attempts to determine the wave functions of the photon have a rich history and date back to the time of the formation of quantum mechanics itself [3,4,5,6]. It should be noted that quantum mechanics, which seeks to describe the most important physical phenomena of nature, must have not only a qualitative but also a quantitative description of the wave function of such a fundamental physical particle as the elementary particle of the electromagnetic field – the photon. A number of recent papers have studied single-photon wave packets that are more or less localized in spacetime using the standard representation of quantum field theory (QFT) [7,8,9]. However, in QFT the photon is associated with a classical solution of the 4-vector potential, which contains singularities that are not physical, since a change in gauge does not result in any change in physical properties. It is obvious that within the framework of the secondary quantization formalism with creation and annihilation operators it is impossible to obtain complete information about a photon and even more so to describe the wave state. In addition, when studying elementary atomic-molecular and nuclear processes with the emission or absorption of photons, it is necessary to be able to describe the wave state of the photon for the correct calculation of the cross sections of the processes under consideration and, ultimately, for their control.

The approaches outlined above make the definition of the photon wave function very vague and difficult to explain, which does not allow us to describe with satisfactory accuracy phenomena, in particular those occurring with the participation of ultra-short light pulses or other spatially complex phenomena with wave functions, the superposition of which leads to interference phenomena. In this connection, Dirac wrote (see [10]): The essential point is the connection between each of the translational states of the photon and one of the wave functions of ordinary wave optics, but he, unfortunately, was unable to formulate this connection. He also put forward another important idea, namely: each photon interferes only with itself, which, in other words, implies the strict stability of the photon and the existence of its wave functions, the superposition of which leads to interference phenomena.

It should be noted that the generally accepted and long-used concept of the corpuscular-wave duality of light in a number of cases does not agree well with the subtle quantum phenomena. The example of the “observer effect" in Young’s experiment on the diffraction of light by two slits, where a photon "recognizes" that it is being observed and changes its "behavior", turning from a wave into a particle, is one of the most difficult to understand in quantum mechanics [10,11]. This raises a key question: can an electromagnetic wave be considered to consist of particles called photons, or should it be considered a field, and the number of photons simply a parameter we attribute to the quantum states of the electromagnetic field [12]?

In any case, experiments studying single photons are currently developing rapidly (see for example [13,14,15]), and a single photon represents a physical reality that has a region of spacetime localization and can be described in terms of quantum mechanics. Moreover, the use of single photons in quantum metrology [16], biology and the fundamentals of quantum physics [17] is today the forefront of scientific and technological development [18,19].

The main goal of this paper is to demonstrate the possibility of constructing a closed, consistent quantum theory of a single photon, allowing one to study the full set of its characteristics under various conditions, taking into account its spin and spacetime structure, which plays an important role in processes involving the photon. Another important goal of the work is to study the probability of the transition of a single photon into an entangled two-photon state as a result of multiple random scattering of a photon on two-level quantum dots with absorption and emission in a nanowaveguide [20], which is of key importance for the organization of quantum communication.

We start our study with the basic first-order partial differential equation for the spin-1 Bose-particle wave function, which has been developed over the last three decades [1,13,21,22,23], and generalize it in various aspects.

In Section 1, we briefly outline the problem of describing a photon using terms of quantum mechanics, i.e., a wave function and a wave equation taking into account its spin.

In Section 2, we present the problem statement and derive a system of second-order PDEs describing the three-component wave function of a single photon in 4D Minkowski spacetime with fabric. Using the components of the vector wave function, the total wave function of an single photon is defined. Taking into account the homogeneity of space along the direction of photon propagation, a certain generalized wave equation of the Klein-Gordon-Fock type is derived for the full wave function of a single photon.

Section 3 examines the quantum motion of a single photon in a nanowaveguide. In particular, it has been proven that under certain conditions the quantum motion of a photon in a nanowaveguide is reduced to the problem of a 2D quantum harmonic oscillator (QHO) with time-dependent frequencies.

In Section 4, we consider the problem of quantum motion of a photon in a nanowaveguide, in which elastic and inelastic processes with absorption and emission of photons continuously occur. The problem is investigated within the framework of a stochastic differential equation (SDE) of the Langevin-Schrödinger (L-S) type. In particular, using a low-dimensional system of stochastic differential equations (SDEs) of the Langevin type, we separate the variables in the original L-Sch equation and obtain the wave function describing orthogonal probabilistic process on the extended space. Using a Langevin-type SDE system, in the limit of statistical equilibrium, we derive a Fokker-Planck-type equation for environmental fields. The geometric and topological features of the arising additional subspace on which the distribution density of the environment fields is specified are analyzed.

In Section 5, the Fokker-Planck measure of the functional space is constructed, which allows, through functional integration, to construct the mathematical expectation (ME) of various random parameters “photon+environment” of the joint system (JS). In particular, using the ME definition, the first four quantized states of a photon in a nanowaveguide were constructed as two-fold integral representations, where the integrand is the solution of the system of two second-order PDEs. In this section, we also construct the MEs of the first four elements of the reduced density matrix (RDM) in the form of double integral representations, where, however, the integrand is the solution of one second-order PDE.

In Section 6, the quantum entropy of a ground-state photon propagating in a nanowaveguide is studied. First, using the element of the reduced density matrix for the ground state, the standard von Neumann quantum entropy is constructed. In this section, we also calculate the quantum entropy of the photon subsystem of the inseparable self-consistent “photon+environment" JS, taking into account the formation of a quantized small environment (SE). Analyzing both expressions for quantum entropies, we are convinced that when two photons are formed as a result of inelastic processes accompanied by the absorption and emission of photons in the medium, they become entangled through the entropies of individual photons.

In Section 7, using the ME of the wave function of a single photon in a nanowaveguide, the elements of the S-matrix of transitions to different Bell states are constructed. or a specific model of the photon wave function, the elements of the S-matrix of transitions are calculated and represented as fourfold integral representations.

In Section 8, we define the Neumann initial-boundary value problem for solving the second-order PDE and for the system of two second-order PDEs needed to construct the ME of the photon wave function scattered by two-level quantum dots randomly located in a nanowaveguide.

In Section 9, a mathematical algorithm for numerical modeling of a system of partial differential equations is developed, which makes it possible to calculate and visualize the probability density of the distribution of electromagnetic field energy in three-dimensional space for various quantum states and thereby prove the existence of structure and the possibility of structuring a photon. Calculations and visualizations of various elements of the reduced density matrix are presented, and, most importantly, it is shown that under certain initial conditions, the developed representation describes the continuous evolution of the quantum state of a photon from the moment of absorption of the particle by the environment and its subsequent transition to a free state, i.e., emission.

In Section 10, we discuss the obtained theoretical and numerical results and outline directions for future research.

Section 11 is an appendix in which we calculate the elements of the transition matrix $\mathrm{S}$ to Bell states using a Gaussian model of the distribution of electromagnetic fields in two single photons.

2. Photon in 4D spacetime with fabric

2.1. Quantum Equations for the Wave Function of a Single Photon

As is known, a photon has electromagnetic energy, which at any given moment in time is localized in space. Based on this, Byalynitsky-Birulya [1,21,22] and Sipe [23] considered it appropriate to describe the distribution of the photon energy density using a wave function. The equation satisfied by this mean-energy density wave function is usually called the photon’s equation of motion — the Byalynitsky-Birulya-Siepe wave equation. As shown, this equation can be derived by parallel reasoning using Dirac’s approach to determining the wave equation for a single electron, a spin-1/2 particle, from Einstein kinematics [24], and it agrees with the equation obtained by the authors [1,21,22,23].

Exactly the same basic three-component vector equation for a massless spin-1 Bose particle can be obtained from the Yang-Mills equation within the $SU\left(2\right)\otimes U\left(1\right)$ gauge symmetry group for Abelian fields [25]:

$${\partial }_t{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)\pm c\left(\mathbf{r}\right)\left(\mathbf{S}·\nabla \right){\mathsf{\Psi }}_{\pm }(\mathbf{r},t)=0,{\partial }_t\equiv \partial /\partial t,\mathbf{r}=\mathbf{r}(x,y,z),$$

(1)

where $c\left(\mathbf{r}\right)$ denotes the speed of light in an in-homogeneous medium, the symbol $\nabla =\left\{\partial /\partial x;\partial /\partial y;\partial /\partial z\right\}$ denotes a vector differential operator, ${\mathsf{\Psi }}_{+}(\mathbf{r},t)$ and ${\mathsf{\Psi }}_{-}(\mathbf{r},t)$ denote the vector wave functions of photons of both helicities, corresponding to left and right polarizations, which can be represented as a three-component wave function:

$${\mathsf{\Psi }}_{\pm }(\mathbf{r},t)=\left(\begin{array}{c}{\psi }_{\pm }^x(\mathbf{r},t)\\ {\psi }_{\pm }^y(\mathbf{r},t)\\ {\psi }_{\pm }^z(\mathbf{r},t)\end{array}\right),$$

(2)

where ${\psi }_{\pm }^x(\mathbf{r},t),{\psi }_{\pm }^y(\mathbf{r},t)$ and ${\psi }_{\pm }^z(\mathbf{r},t)$ are the partial wave functions of the photon. In the future, the equation (1)-(2) can be interpreted as the motion of a quantum particle-photon in four-dimensional Minkowski spacetime with fabric.

The set of matrices $\mathbf{S}=(S_x,S_y,S_z)$ in equation (1) describes the infinitesimal rotation of a massless spin-1 particle, which can be represented as follows:

$$S_x=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\\ 0& i& 0\end{array}\right),S_y=\left(\begin{array}{ccc}0& 0& i\\ 0& 0& 0\\ -i& 0& 0\end{array}\right),S_z=\left(\begin{array}{ccc}0& -i& 0\\ i& 0& 0\\ 0& 0& 0\end{array}\right).$$

Our main goal now will be to derive the wave equation for a single photon from the vector equation (1). It is important to note that when moving from a system of first-order PDEs (1) to a system of second-order PDEs (3), we usually obtain a wider class of solutions, including those that may not satisfy the original system of equations. Nevertheless, this approach is often used in the theory of electromagnetism, quantum mechanics and other areas of physics where the wave properties of matter are studied. This is justified by the fact that the solution of the second-order equation becomes more diverse, and its description requires much more information. The derivation of the second-order Klein-Gordon-Fock equation from the first-order Dirac equation is in this sense a benchmark example (see for example [14]).

Recall that when moving from first-order to second-order equations, information about initial conditions or constraints that are important for first-order equations is often lost, so it is always necessary to check which solutions make physical sense and which do not in the context of the original problem.

Theorem 1. 

Let the photon velocity have dispersion in 3D space, i.e. $c(x,y,z)<c_0$, then its quantum state will be described by a three-component wave function (2) satisfying a system of second-order partial differential equations of the form:

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial t^2}}+\left\{{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\partial {\psi }_{\pm }^y}{\partial y}}+{\displaystyle \frac{\partial {\psi }_{\pm }^z}{\partial z}}\right)-\left({\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial z^2}}\right)\right\}-{\displaystyle \frac{c_x}{c}}\left\{{\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial y}}+{\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial z}}\right\}=0,\\ \hfill {\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial t^2}}+\left\{{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial x}}+{\displaystyle \frac{\partial {\psi }_{\pm }^z}{\partial z}}\right)-\left({\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial z^2}}\right)\right\}-{\displaystyle \frac{c_y}{c}}\left\{{\displaystyle \frac{\partial {\psi }_{\pm }^y}{\partial x}}+{\displaystyle \frac{\partial {\psi }_{\pm }^y}{\partial z}}\right\}=0,\\ \hfill {\displaystyle \frac{1}{c^2}}{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^z}{\partial t^2}}+\left\{{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial x}}+{\displaystyle \frac{\partial {\mathsf{\Psi }}_{\pm }^y}{\partial y}}\right)-\left({\displaystyle \frac{{\partial }^2{\psi }_{\pm }^z}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^z}{\partial y^2}}\right)\right\}-{\displaystyle \frac{c_z}{c}}\left\{{\displaystyle \frac{\partial {\psi }_{\pm }^z}{\partial y}}+{\displaystyle \frac{\partial {\psi }_{\pm }^z}{\partial x}}\right\}=0,\end{array}$$

(3)

where $c_0=const$ is the speed of light in vacuum, and $c_{\eta }={\partial }_{\eta }c(x,y,z)$ denote the first-order derivatives of the speed of light with respect to the coordinates $\eta =x,y,z$, respectively.

Proof. 

Using the following identity:

$$(\mathbf{S}·\nabla )\mathbf{F}=i\nabla \times \mathbf{F},$$

(4)

we can rewrite the equations (1) as follows:

$$i{\partial }_t{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)\mp c\left(\mathbf{r}\right)\nabla \times {\mathsf{\Psi }}_{\pm }(\mathbf{r},t)=0.$$

(5)

Now differentiating the equation (5) with respect to time we find:

$$i{\partial }_t^2{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)\mp c\left(\mathbf{r}\right)\nabla \times {\partial }_t{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)=0.$$

(6)

Substituting the value of the derivative ${\partial }_t{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)$ from (5), we can write equation (6) as:

$${\partial }_t^2{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)+c\left(\mathbf{r}\right)\nabla \times \left[c\left(\mathbf{r}\right)\nabla \times {\mathsf{\Psi }}_{\pm }(\mathbf{r},t)\right]=0.$$

(7)

By opening the vector product of type $a\times [b\times c]$ in equation (7), it can be written as:

$${\partial }_t^2{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)+c\left(\mathbf{r}\right)\left\{\nabla 〈c\left(\mathbf{r}\right)\nabla ,{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)〉-〈\nabla ,c\left(\mathbf{r}\right)\nabla 〉{\mathsf{\Psi }}_{\pm }(\mathbf{r},t)\right\}=0,$$

(8)

where $\langle ,\rangle$ denotes the scalar product.

Finally, by performing simple calculations from equation (8), one can obtain a system of three coupled second-order PDEs (3) that describe the quantum motion of single photons of both left- and right-handed helicity. Considering that when deriving the system of second-order equations (3), only the original system of first-order equations (1) is used, i.e. the identity, it can be considered proven that the equations (1) implies the system of equations (3) and, correspondingly, its solutions are contained in solutions of equations (1). □

Definition 1. The sum of the partial wave functions will be called the total wave function of the photon:

$${\psi }_{\pm }(\mathbf{r},t)={\displaystyle \frac{1}{3}}\left\{{\psi }_{\pm }^x(\mathbf{r},t)+{\psi }_{\pm }^y(\mathbf{r},t)+{\psi }_{\pm }^z(\mathbf{r},t)\right\}.$$

(9)

It is now important to clarify the question of the physical meaning of the individual components of the vector wave function (2). As is easy to see, each component of the triplet (2) describes the quantum state of a photon along a certain coordinate with a changing speed of light.

Since the system of equations (3) is invariant with respect to the change of the sign of helicity, then in what follows it will be convenient to write the photon wave function in the form of the Riemann–Silberstein vector representation [26,27]:

$${\mathcal{A}}_{\pm }(\mathbf{r},t)={\displaystyle \frac{1}{\sqrt{2}}}\left\{{\displaystyle \frac{\mathbf{D}(\mathbf{r},t)}{\sqrt{\epsilon }}}\pm i{\displaystyle \frac{\mathbf{B}(\mathbf{r},t)}{\sqrt{\mu }}}\right\},$$

(10)

where ${\mathcal{A}}_{\pm }(\mathbf{r},t)$ denotes the wave function of a photon, $\mathbf{D}(\mathbf{r},t)$ and $\mathbf{B}(\mathbf{r},t)$ are the vectors denoting the displacement field and the magnetic field, respectively, in addition, $\epsilon$ and $\mu$ determine the magnetic and electrical permeability or refractive indices of the medium.

2.2. Motion of a Photon Under the Condition of Homogeneity of the Medium Along the z-Axis

Now let us assume that the medium along the z axis is homogeneous, that is, the derivative of the photon velocity with respect to the z coordinate is equal to zero $c_z\left(\mathbf{r}\right)=0$. In this case, the z-component of the wave function can be set equal to zero ${\psi }_{\pm }^z=0$. In addition, in order for third equation in (3) to be identically equal to zero, the condition:

$${\partial }_x{\psi }_{\pm }^x(\mathbf{r},t)=-{\partial }_y{\psi }_{\pm }^y(\mathbf{r},t),$$

(11)

must be met.

In other words, in the case under consideration, when the medium is homogeneous along the z axis, the state of the photon is described by a singlet wave function:

$${\mathsf{\Psi }}_{\pm }^{II}(\mathbf{r},t)=\left(\begin{array}{c}{\psi }_{\pm }^x(\mathbf{r},t)\\ {\psi }_{\pm }^y(\mathbf{r},t)\end{array}\right).$$

(12)

Now, assuming that the properties of the medium along the x and y coordinates are the same, taking into account the condition (11), from the system of PDEs (3) we can obtain a system of two second-order PDEs describing the motion of a photon of both helicities:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial t^2}}-c^2\left\{{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^x}{\partial z^2}}\right\}-c{c}_x\left\{{\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial y}}+{\displaystyle \frac{\partial {\psi }_{\pm }^x}{\partial z}}\right\}=0,\\ \hfill {\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial t^2}}-c^2\left\{{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{\psi }_{\pm }^y}{\partial z^2}}\right\}-c{c}_y\left\{{\displaystyle \frac{\partial {\psi }_{\pm }^y}{\partial x}}+{\displaystyle \frac{\partial {\psi }_{\pm }^y}{\partial z}}\right\}=0.\end{array}$$

(13)

It is worth noting that, despite the apparent independence of the PDEs in (13), they are connected by condition (11) and therefore they can only be considered as separate conditionally.

As we can see, the equations in (13) are similar in structure and the only difference is that the first equation describes the motion of a photon with variable speed along the “x” coordinate, while in the second the speed is variable along the “y” coordinate. For a consistent study of the problem, it is important to derive a general equation for the photon, taking into account the dispersion of the photon velocity in both coordinates. In this connection, a natural question arises, namely: what additional condition must be imposed on the partial wave functions ${\psi }_{\pm }^x$ and ${\psi }_{\pm }^y$ in order to obtain from (13) one independent equation for the general wave function of the photon, which can be represented in the form: ${\mathsf{\Phi }}_{\pm }(\mathbf{r},t)={\displaystyle \frac{1}{2}}\left[{\psi }_{\pm }^x(\mathbf{r},t)+{\psi }_{\pm }^y(\mathbf{r},t)\right]$?

Below we will omit the symbol “±" that denotes the polarization of the photon, since we have verified that the equations (13) are invariant with respect to sign reversal.

Theorem 2. 

Let the components ${\psi }^x(\mathbf{r},t)$ and ${\psi }^y(\mathbf{r},t)$ of a singlet wave function ${\psi }_{\pm }^{II}(\mathbf{r},t)$ satisfy the condition:

$$\begin{array}{c}\hfill c_x\left\{{\partial }_y+{\partial }_z\right\}{\psi }^x+c_y\left\{{\partial }_x+{\partial }_z\right\}{\psi }^y=0,\end{array}$$

(14)

then the total wave function of the photon $\mathsf{\Phi }(\mathbf{r},t)$ will be the solution of the following second-order PDE, which will be called the generalized Klein-Gordon-Fock equation:

$$\square \mathsf{\Phi }=\widehat{\mathcal{L}}\mathsf{\Phi },$$

(15)

where $\square =c^{-2}{\partial }_t^2-{\nabla }^2$ is the d’Alembert operator and $\widehat{\mathcal{L}}=c_x{c}^{-1}\left\{{\partial }_y+{\partial }_z\right\}+c_y{c}^{-1}\left\{{\partial }_x+{\partial }_z\right\}$ will be further called the dispersion operator of the speed of light.

Proof. 

If we add the two wave equations of the system (13), we can find the following equation:

$$\square \mathsf{\Phi }-{\displaystyle \frac{c_x}{\sqrt{2}c}}\left\{{\displaystyle \frac{\partial }{\partial y}}+{\displaystyle \frac{\partial }{\partial z}}\right\}{\psi }^x-{\displaystyle \frac{c_y}{\sqrt{2}c}}\left\{{\displaystyle \frac{\partial }{\partial x}}+{\displaystyle \frac{\partial }{\partial z}}\right\}{\psi }^y=0.$$

(16)

Equation (16) can be rewritten as:

$$\begin{array}{c}\hfill \square \mathsf{\Phi }-\left[{\displaystyle \frac{c_x}{c}}\left\{{\displaystyle \frac{\partial }{\partial y}}+{\displaystyle \frac{\partial }{\partial z}}\right\}+{\displaystyle \frac{c_y}{c}}\left\{{\displaystyle \frac{\partial }{\partial x}}+{\displaystyle \frac{\partial }{\partial z}}\right\}\right]\mathsf{\Phi }+\\ \hfill {\displaystyle \frac{c_x}{\sqrt{2}c}}\left\{{\displaystyle \frac{\partial }{\partial y}}+{\displaystyle \frac{\partial }{\partial z}}\right\}{\psi }^x+{\displaystyle \frac{c_y}{\sqrt{2}c}}\left\{{\displaystyle \frac{\partial }{\partial x}}+{\displaystyle \frac{\partial }{\partial z}}\right\}{\psi }^y=0.\end{array}$$

(17)

Now, taking into account the condition (14) relating the partial wave functions, we obtain a closed equation (15) for the general wave function of the photon $\mathsf{\Phi }(\mathbf{r},t)$. □

In the special case when in the $\{x,y\}$ plane the structure of space is such that $\widehat{\mathcal{L}}\equiv 0$, the quantum motion of an individual photon is described by the Klein-Gordon-Fock equation:

$$\square \mathsf{\Phi }=0.$$

(18)

Recall that the Klein-Gordon-Fock equation describes particles with spin 0, however, taking into account the representation of the wave function (10), we can say that the equation also allows us to describe particles with projections of spins 1 and -1.

Finally, we note that when we consider the case $c_x=c_y=0$ or, what is the same, $\widehat{\mathcal{L}}\equiv 0$, this does not mean that the speed of light is constant throughout space.

3. A Photon Moving in a Nanowaveguide

Let the photon propagate freely along the coordinate “z”, then the solution of equation (18) can be represented as:

$$\mathsf{\Phi }(\mathbf{r},t)=e^{i{\omega }_{0}t+i{k}_{0}z}F\left(\mathbf{r}\right),$$

(19)

where ${\omega }_0$ and $k_0$ denote the frequency and momentum of a photon in a free vacuum with the speed of light $c_0$.

Substituting (19) into (18), we can find the following equation:

$$\left\{i{k}_0{\displaystyle \frac{\partial }{\partial z}}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{{\partial }^2}{\partial y^2}}+{\displaystyle \frac{{\partial }^2}{\partial z^2}}\right)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\omega }_0^2}{c^2}}-k_0^2\right)\right\}F=0.$$

(20)

Let us assume that the wave function $F\left(\mathbf{r}\right)$ along the direction of propagation z changes smoothly enough so that:

$$∥{\displaystyle \frac{\partial F}{\partial z}}∥\gg ∥{\displaystyle \frac{{\partial }^{2}F}{\partial z^2}}∥.$$

(21)

Taking into account the estimate (21), equation (20) can be simplified by expanding the term $1/c^2(x,y,z)$ in a Taylor series in the coordinates x and y and retaining the quadratic terms:

$$i{\displaystyle \frac{\partial \mathcal{F}}{\partial \tau }}+{\displaystyle \frac{1}{2}}\left\{{\displaystyle \frac{{\partial }^2\mathcal{F}}{\partial x^2}}+{\displaystyle \frac{{\partial }^2\mathcal{F}}{\partial y^2}}\right\}-{\displaystyle \frac{1}{2}}\left\{A_1\left(\tau \right)x^2+A_2\left(\tau \right)y^2+2A_3\left(\tau \right)xy\right\}\mathcal{F}=0,$$

(22)

where $\tau =z/k_0$ is a new chronologizing parameter and $\mathcal{F}(x,y,\tau )\approx F(x,y,z)$.

Note that in deriving the equation (22) the equality ${({\omega }_0/c_0)}^2-k_0^2=0$ was used and, in addition, the following notations were made:

$$A_1\left(\tau \right)=2{\displaystyle \frac{{\omega }_0^2}{c_0^3}}{c}_{xx}(0,0,z),A_2\left(\tau \right)=2{\displaystyle \frac{{\omega }_0^2}{c_0^3}}{c}_{yy}(0,0,z),A_3\left(\tau \right)=2{\displaystyle \frac{{\omega }_0^2}{c_0^3}}{c}_{xy}(0,0,z),$$

where ${\partial }_{\eta }^{2}c(0,0,z)=c_{\eta \eta }(0,0,z)$, and the coefficients $A_i\left(\tau \right),(i=1,2,3)$ arise when expanding the function $1/c^2\left(\mathbf{r}\right)$ into a Taylor series including terms of the second degree, taking into account that $c_x(z,0,0)=c_y(z,0,0)=0$. If we now assume that the second derivatives of the speed of light are positive, that is:

$$c_{xx}(0,0,z)>0,c_{yy}(0,0,z)>0,c_{xy}(0,0,z)>0,$$

then the functions $\{A_1\left(\tau \right),A_2\left(\tau \right),A_3\left(\tau \right)\}>0$ will also be positive, which will allow us to interpret them as squares of frequencies and reduce the problem of photon evolution in the $\{x,y\}$ plane to the model of a 2D quantum oscillator with a time-dependent frequency.

Thus, the state of a photon in the $\{x,y\}$ plane will be described by a two-dimensional partial differential equation:

$$i{\partial }_{\tau }\mathcal{F}=\widehat{\mathrm{H}}(x_1,x_2,\tau )\mathcal{F},$$

(23)

where the following notations are made $x=x_1$ and $y=x_2$, in addition, $\widehat{H}$ denotes the Hamiltonian of the quantum system:

$$\widehat{\mathrm{H}}(x_1,x_2,\tau )={\displaystyle \frac{1}{2}}\sum _{l=1}^2\left[-{\displaystyle \frac{{\partial }^2}{\partial x_l^2}}+{\mathsf{\Omega }}_l^2\left(\tau \right)x_l^2\right]+\mathsf{\Theta }\left(\tau \right)x_1{x}_2,$$

(24)

in addition, ${\mathsf{\Omega }}_l^2\left(\tau \right)=2{\omega }_0^2{c}_0^{-3}{c}_{x_l{x}_l}(0,0,z)$ and ${\mathsf{\Theta }}_l\left(\tau \right)=2{\omega }_0^2{c}_0^{-3}{c}_{x_1{x}_2}(0,0,z)$.

Further, for the sake of certainty, we will consider the case when the oscillator frequencies along different axes are equal, i.e. ${\mathsf{\Omega }}_1\left(\tau \right)={\mathsf{\Omega }}_2\left(\tau \right)=\mathsf{\Omega }\left(\tau \right)$. The latter circumstance allows us to implement the following coordinate transformations in the equation (23)-(24):

$$q_1={\displaystyle \frac{x_1-x_2}{\sqrt{2}}},q_2={\displaystyle \frac{x_1+x_2}{\sqrt{2}}},$$

(25)

which ultimately reduces the Hamiltonian of the problem to a diagonal form:

$$\widehat{\mathcal{H}}(q_1,q_2,\tau )={\displaystyle \frac{1}{2}}\sum _{l=1}^2\left[-{\displaystyle \frac{{\partial }^2}{\partial q_l^2}}+{\omega }_l^2\left(\tau \right)q_l^2\right].$$

(26)

Recall that ${\omega }_l\left(\tau \right)$ denotes the new effective frequency of the 1D quantum harmonic oscillator (QHO), which is determined by the following formula:

$${\omega }_l^2\left(\tau \right)={\mathsf{\Omega }}^2\left(\tau \right)-{(-1)}^l\mathsf{\Theta }\left(\tau \right),l=1,2.$$

(27)

Note that below we will consider the most general case when:

$$\underset{\tau \to \pm \infty }{lim}{\omega }_l\left(\tau \right)={\omega }_l^{\pm },$$

(28)

where ${\omega }_1^{\pm }$ and ${\omega }_2^{\pm }$ are constants, in addition, in the future we will understand $-\infty$ as the starting point and replace it with 0, i.e. $\tau \to 0$.

The solution of the Schrödinger equation (23) with the Hamiltonian (26) can obviously be represented in factorized form:

$$\mathcal{F}\left(\mathbf{n}\right|q_1,q_2,\tau )=\prod _{l=1}^2{\mathcal{F}}_l\left(n_l\right|q_l,\tau ),\mathbf{n}=(n_1,n_2),$$

(29)

where $n_1$ and $n_2$ denote vibrational quantum numbers, the wave function ${\mathcal{F}}_l\left(n_l\right|q_l,\tau )$ is a solution of $1D$ QHO with an arbitrary non-stationary frequency, which has the following form [28]:

$$\begin{array}{c}\hfill {\mathcal{F}}_l\left(n_l\right|q_l,\tau )=\sqrt{{\displaystyle \frac{g_l}{2^{n_l}{n}_l!{\sigma }_l}}}exp\left\{-i\left(n_l+{\displaystyle \frac{1}{2}}\right){\omega }_l^{-}{\int }_0^{\tau }{\displaystyle \frac{d{\tau }^{\prime }}{{\sigma }_l^2}}+{\displaystyle \frac{1}{2}}\left(i{\displaystyle \frac{\dot{{\sigma }_l}}{{\sigma }_l}}-{\displaystyle \frac{{\omega }_l^{-}}{{\sigma }_l^2}}\right)q_l^2\right\}\times \\ \hfill H_{n_l}\left(\sqrt{{\omega }_l^{-}}{\displaystyle \frac{q_l}{{\sigma }_l}}\right),g_l=\sqrt{{\omega }_l^{-}/\pi },{\omega }_l^{-}=\underset{\tau \to 0}{lim}{\omega }_l\left(\tau \right).\end{array}$$

(30)

Recall that $n_l=0,1...$ denotes the vibrational quantum number and ${\omega }_l^{-}$ is the frequency of the oscillator in the initial state, i.e. at the beginning of the evolution process. In addition, the function ${\sigma }_l\left(\tau \right)$ is a solution to the classical oscillator equation, which we will call the reference equation:

$$\ddot{\xi }+{\omega }_l^2\left(\tau \right)\xi =0,$$

(31)

where

$${\xi }_l\left(\tau \right)={\sigma }_l\left(\tau \right)e^{i{\gamma }_l\left(\tau \right)},{\gamma }_l\left(\tau \right)={\omega }_l^{-}{\int }_0^{\tau }{\displaystyle \frac{d{\tau }^{\prime }}{{\sigma }_l^2\left({\tau }^{\prime }\right)}}.$$

The wave function (30) in the limit $\tau \to -\infty$, i.e. in the channel $\left(in\right)$ goes into the asymptotic state:

$${\mathsf{\Phi }}_{\left(in\right)}\left(\mathbf{n}\right|\mathbf{q},\tau )=\prod _{l=1}^2{e}^{-i(n_l+1/2){\omega }_l^{-}\tau }{\mathsf{\Phi }}_{n_l}\left(q_l\right),\mathbf{q}\equiv (q_1,q_2),$$

(32)

in addition:

$$\mathsf{\Phi }\left(n_l\right|q_l)={\left({\displaystyle \frac{\sqrt{{\omega }_l^{-}/\pi }}{2^{n_l}{n}_l!}}\right)}^{1/2}{e}^{-{\omega }_l^{-}{q}_l^2/2}{H}_{n_l}\left(\sqrt{{\omega }_l^{-}}{q}_l\right),n_l=0,1,...$$

(33)

where ${\mathsf{\Phi }}_{n_l}\left(q_l\right)$ denotes an one particle bosonic Fock state and the function $H_{n_l}\left(q_l\right)$ is a Hermitian polynomial. The set of functions ${\mathsf{\Phi }}_{\left(in\right)\left(\mathbf{n}\right|\mathbf{q},\tau )}$ form orthonormal bases in the Hilbert space.

Finally, it is important to note that the wave function ${\mathcal{F}}_l\left(n_l\right|q_l,\tau )$ in the Hilbert space forms an orthonormal basis:

$${\int }_{-\infty }^{+\infty }{\mathcal{F}}_l\left(n_l\right|q_l,\tau ){\mathcal{F}}_{l^{\prime }}^{*}\left(n_l^{\prime }\right|q_l,\tau )d{q}_l={\delta }_{n_l{n}_l^{\prime }},$$

(34)

where the sign “*” over the wave function denotes the complex conjugate, and ${\delta }_{n_l{n}_l^{\prime }}$ is the Kronecker delta function. It is obvious that the full wave function (29) also satisfies the orthogonality condition.

4. The Wave Function of a Photon at Its Random Motion in a Nanowaveguide

4.1. Description of the Quantum State of a Photon in a Random Environment

Let us now assume that a photon, moving in a nanowaveguide, experiences random elastic and inelastic interactions with the surrounding medium. Recall that by inelastic collisions we mean the possibility of the photon being absorbed by the medium. Then, taking into account the Equations (23) and (24), the motion of the photon will be described by a stochastic differential equation (SDE) of the Langevin-Shcrödinger type:

$$i{\partial }_{\tau }F=\widehat{\mathrm{H}}(x_1,x_2,\tau ;\left\{f\right\})F,-\infty <x_1,x_2<+\infty ,$$

(35)

where $\left\{f\right\}=\{f^{\left(r\right)}\left(\tau \right),f^{\left(i\right)}\left(\tau \right)\}$ denotes a complex random process consisting of real $f^{\left(r\right)}\left(\tau \right)$ and imaginary $f^{\left(i\right)}\left(\tau \right)$ parts, the forms of which will be defined below. Note that $f^{\left(r\right)}\left(\tau \right)$ and $f^{\left(i\right)}\left(\tau \right)$ characterize, respectively, elastic and inelastic scattering of a photon in a medium.

In (35) the evolution operator $\widehat{\mathrm{H}}$ characterizes the Hamiltonian of a photon and its environment and has the following form:

$$\begin{array}{c}\hfill \widehat{\mathrm{H}}(x_1,x_2,\tau ;\left\{f\right\})={\displaystyle \frac{1}{2}}\sum _{l=1}^2\left[-{\displaystyle \frac{{\partial }^2}{\partial x_l^2}}+{\mathsf{\Omega }}^2(\tau ;\left\{f\right\})x_l^2\right]+\mathsf{\Theta }(\tau ;\left\{f\right\})x_1{x}_2,\end{array}$$

(36)

where ${\mathsf{\Omega }}^2(\tau ;\left\{f\right\})$ and $\mathsf{\Theta }(\tau ;\{f\left\}\right)$ are random functions.

By performing similar transformations and calculations as in Section 3, we can find the photon wave function $F(x_1,x_2,\tau ;\left\{f\right\})\to \mathcal{F}(q_1,q_2,\tau ;\left\{f\right\})$, which now represents a random complex probabilistic process (see for example [30]). Note that this wave function can be formally found from the regular photon wave function $\mathcal{F}(q_1,q_2,\tau )$ by replacing ${\sigma }_l\left(\tau \right)\to {\sigma }_l(\tau ;\left\{f\right\})$ in the expression for the wave function of the 1D harmonic quantum oscillator ${\mathcal{F}}_l(q_l,\tau )$ (see Equations (29) and (30)). Recall that the random function ${\sigma }_l(\tau ;\left\{f\right\})=\left|\xi \left(\tau \right)\right|$, where $\xi \left(\tau \right)$ is a complex probabilistic process satisfying the following SDE:

$$\ddot{\xi }+{\omega }_l^2(\tau ;\left\{f\right\})\xi =0.$$

(37)

Note that the reference low-dimensional SDE (i.e. 1D) (37), which is obtained from equation (31) by replacing $\omega \left(\tau \right)\to \omega (\tau ;\{f\left\}\right)$, still needs to be defined, since the canonical SDE must be first order. It is easy to verify that the wave function ${\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f\right\})$ in the $\left(in\right)$ channel goes into an asymptotic state ${\mathsf{\Phi }}_{\left(in\right)}\left(\mathbf{n}\right|\mathbf{q},\tau )$ (see (32)).

Finally, it is important to note that, despite the fact that the wave function of the photon in this case is random complex probabilistic process, it nevertheless constitutes an orthonormal basis in the Hilbert space and satisfies the condition:

$$\int \int \mathcal{F}\left(\mathbf{n}\right|\mathbf{q},\tau ;\left\{f\right\}){\mathcal{F}}^{*}\left(\mathbf{m}\right|\mathbf{q},\tau ;\left\{f\right\})d\mathbf{q}=\prod _{l=1}^2{\delta }_{n_l{m}_l},$$

(38)

where $\mathbf{n}=(n_1,n_2)$ denote vibrational quantum numbers.

Note that the wave function $\mathcal{F}\left(\mathbf{n}\right|\mathbf{q},\tau ;\left\{f\right\})\in L_2\left(\mathsf{\Xi }\right)$, in addition, $\mathsf{\Xi }\cong {\mathbb{R}}^2\otimes {\mathsf{\Xi }}_{\left\{\xi \right(\tau \left)\right\}}$ is an extended space, where ${\mathbb{R}}^2$ and ${\mathsf{\Xi }}_{\left\{\xi \right(\tau \left)\right\}}$ denote 2D Euclidean and uncountable-dimensional functional spaces respectively. Since the functional space ${\mathsf{\Xi }}_{\left\{\xi \right(\tau \left)\right\}}$ will play an important role in subsequent constructions of the theory, it will be very important to construct its exact measure.

4.2. Distribution of Environmental Fields in the Limit of Statistical Equilibrium

To correctly determine the SDE (37), it is first necessary to determine the form of the random function ${\omega }_l(\tau ;\left\{f\right\})$ included in this equation. For simplicity we will represent the effective frequency as the sum of two functions:

$${\omega }_l^2(\tau ;\left\{f\right\})={\omega }_l^2\left(\tau \right)+f_l\left(\tau \right),f_l\left(\tau \right)=f_l^{\left(r\right)}\left(\tau \right)+i{f}_l^{\left(i\right)}\left(\tau \right),$$

where ${\omega }_l\left(\tau \right)>0$ is a regular function, for which we will further use the following analytical representation:

$$\begin{array}{c}\hfill {\omega }_l\left(\tau \right)=a_l+{\displaystyle \frac{1}{{\gamma }_l}}\left(b_l+tanh\left({\nu }_l\tau \right)\right),\end{array}$$

(39)

where $a_l,b_l,{\gamma }_l,{\nu }_l>0$ are some constants, the values of which are selected based on the conditions of the problem. As for $f_l^{\left(\upsilon \right)}\left(\tau \right),(\upsilon =r,i)$, they define the random influence of the environment on a 1D QHO. We will assume that above defined functions satisfy the following conditions:

$$\begin{array}{c}\hfill \underset{\tau \to \mp \infty }{lim}{\omega }_l(\tau ;f_l)={\omega }_l^{\mp }>0,\underset{\tau \to -\infty }{lim}{f}_l\left(\tau \right)=0.\end{array}$$

(40)

We will suppose that $f_l^{\left(\upsilon \right)}\left(\tau \right)$ is an independent Gauss-Markov process with zero mean and delta-shaped correlation function:

$$〈f_l^{\left(\upsilon \right)}〉=0,〈f_l^{\left(\upsilon \right)}\left(\tau \right)f_l^{\left(\upsilon \right)}\left({\tau }^{\prime }\right)〉=2{ϵ}_l^{\left(\upsilon \right)}\delta (\tau -{\tau }^{\prime }),$$

(41)

where ${ϵ}_l^{\left(\upsilon \right)}$ denotes the fluctuations power.

The solution of the SDE (37) can be represented as:

$${\xi }_l\left(\tau \right)=\{\begin{array}{c}exp\left(i{\omega }_l^{-}\tau \right),\tau \le {\tau }_0,\hfill \\ {\xi }_l\left({\tau }_0\right)exp\left\{{\mathsf{\Pi }}_l\left(\tau \right)\right\},\tau >{\tau }_0,\hfill \end{array}$$

(42)

where ${\mathsf{\Pi }}_l\left(\tau \right)={\int }_0^{\tau }{\chi }_l\left({\tau }^{\prime }\right)d{\tau }^{\prime }$ and $\tau ={\tau }_0=0$ denotes the moment of switching on the random environment or the initial moment.

Substituting (42) into (37) we find the following nonlinear complex SDE:

$$\dot{{\chi }_l}+{\chi }_l^2+{\omega }_l^2\left(\tau \right)+f_l\left(\tau \right)=0,$$

(43)

where $\dot{{\chi }_l}=d{\chi }_l/d\tau$ and ${\chi }_l\left(0\right)=i{\omega }_l^{-}$.

Let us represent the complex function ${\chi }_l\left(\tau \right)$ in the form of the sum of a real and imaginary parts:

$${\chi }_l\left(\tau \right)=u_1^{\left(l\right)}\left(\tau \right)+i{u}_2^{\left(l\right)}\left(\tau \right),$$

(44)

where $u_1^{\left(l\right)}\left(\tau \right)\in (-\infty ,+\infty )$ and $u_2^{\left(l\right)}\left(\tau \right)\in (0,+\infty ).$

Taking into account the expressions (43) and (44), we can find the following system of nonlinear SDEs [30]:

$$\{\begin{array}{c}{\dot{u}}_1^{\left(l\right)}={\left[u_2^{\left(l\right)}\right]}^2-{\left[u_1^{\left(l\right)}\right]}^2-{\omega }_l^2\left(\tau \right)-f_l^{\left(r\right)}\left(\tau \right),\hfill \\ {\dot{u}}_2^{\left(l\right)}=-2u_1^{\left(l\right)}{u}_2^{\left(l\right)}+f_l^{\left(i\right)}\left(\tau \right),\hfill \end{array}$$

(45)

where fields satisfy the asymptotic conditions:

$${\dot{u}}_1^{\left(l\right)}\left(0\right)=Re[{\dot{\xi }}_l\left(0\right)/{\xi }_l\left(0\right)]=0,{\dot{u}}_2^{\left(l\right)}\left(0\right)=Im[{\dot{\xi }}_l\left(0\right)/{\xi }_l\left(0\right)]={\omega }_l^{-}.$$

Note that depending on the character of random forces $f_1\left(\tau \right)$ and $f_2\left(\tau \right)$ the probabilistic processes ${\chi }_1\left(\tau \right)$ and ${\chi }_2\left(\tau \right)$ can be realized by way of two different scenarios.

At first we will consider the scenario when random forces are independent. It follows that the distribution of environmental fields under conditions of statistical equilibrium can be represented as follows:

$$P(\mathbf{u},\tau |{\mathbf{u}}_0,{\tau }_0)=\prod _{l=1}^2{P}_l\left(u_1^{\left(l\right)},u_2^{\left(l\right)},\tau \right)=\prod _{l=1}^2〈\prod _{j=1}^2\delta \left(u_j^{\left(l\right)}\left(\tau \right)-u_{0j}^{\left(l\right)}\right)〉,$$

where $u_{0j}^{\left(l\right)}=u_j^{\left(l\right)}\left(0\right)$ is the component of the field in the asymptotic ($in$) state, in addition: $\mathbf{u}\equiv \left\{{\mathbf{u}}^{\left(1\right)}=(u_1^{\left(1\right)},u_2^{\left(1\right)});{\mathbf{u}}^{\left(2\right)}=(u_1^{\left(2\right)},u_2^{\left(2\right)})\right\}$. For the simplification of notations in equations below upper indexes of random fields will be omitted.

Using SDE’s (45) for the probability distribution we can find the Fokker-Planck type equation (see [30,31]):

$${\partial }_{\tau }{P}_l={\widehat{L}}_l\left({\omega }_l^{-},{ϵ}_l^{\left(i\right)},{ϵ}_l^{\left(r\right)}|u_1^{\left(l\right)},u_2^{\left(l\right)},\tau \right)P_l.$$

(46)

Recall that the partial evolution operator ${\widehat{L}}_l$ has the following form:

$${\widehat{L}}_l={ϵ}_l^{\left(r\right)}{\left\{{\displaystyle \frac{\partial }{\partial u_1^{\left(l\right)}}}\right\}}^2+{ϵ}_l^{\left(i\right)}{\left\{{\displaystyle \frac{\partial }{\partial u_2^{\left(l\right)}}}\right\}}^2+k_1{\displaystyle \frac{\partial }{\partial u_1^{\left(l\right)}}}+k_2{\displaystyle \frac{\partial }{\partial u_2^{\left(l\right)}}}+4u_1^{\left(l\right)}$$

(47)

where

$$k_1\left(u_1^{\left(1\right)},u_2^{\left(1\right)},\tau \right)={\left[u_1^{\left(l\right)}\right]}^2-{\left[u_2^{\left(l\right)}\right]}^2+{\omega }_l^2\left(\tau \right),k_2\left(u_1^{\left(1\right)},u_2^{\left(1\right)},\tau \right)=2u_1^{\left(l\right)}{u}_2^{\left(l\right)}.$$

It is quite naturally to assume that the solution of the equation (46) satisfies the initial conditions:

$$\begin{array}{c}\hfill P_l(u_1^{\left(l\right)},u_2^{\left(l\right)},{\tau }_0)=\prod _{j=1}^2\delta (u_j^{\left(l\right)}-u_{0j}^{\left(l\right)}),P_l(u_1^{\left(l\right)},u_2^{\left(l\right)},\tau ){|}_{\parallel {\mathbf{u}}^{\left(l\right)}\parallel \to \infty }\to 0,\end{array}$$

(48)

where $∥{\mathbf{u}}^{\left(l\right)}∥=\sqrt{{\left(u_1^{\left(l\right)}\right)}^2+{\left(u_2^{\left(l\right)}\right)}^2}$.

At realization of the second scenario the source of the random forces is the same one and respectively the distribution of fields obviously is defined in the following form:

$$\mathcal{P}(\mathbf{u},\tau |{\mathbf{u}}_0,{\tau }_0)=〈\prod _{l=1}^2\prod _{j=1}^2\delta \left(u_j^{\left(l\right)}\left(\tau \right)-u_{0j}^{\left(l\right)}\right)〉,$$

(49)

which is the solution of the following equation:

$${\partial }_{\tau }\mathcal{P}=\sum _{l=1}^2{\widehat{L}}_l\left({\omega }_l,{ϵ}_l^{\left(i\right)},{ϵ}_l^{\left(r\right)}|u_1^{\left(l\right)},u_2^{\left(l\right)},\tau \right)\mathcal{P}.$$

(50)

The equation (50) can be solved with the initial conditions type (48). Since there is only one source of fluctuations in this case, it should be assumed that: ${ϵ}_1^{\left(i\right)}={ϵ}_2^{\left(i\right)}={ϵ}^{\left(i\right)}$ and ${ϵ}_1^{\left(r\right)}={ϵ}_2^{\left(r\right)}={ϵ}^{\left(r\right)}.$

If the fluctuation powers are small, i.e. $\{{ϵ}^{\left(r\right)},{ϵ}^{\left(i\right)}\}\ll 1$, then we can assume that the Fokker-Planck equation (46) is defined in an additional subspace ${\mathbb{R}}_{\{{\mathbf{u}}^{\left(l\right)},\tau \}}^2$, which is described by Euclidean geometry, and the coordinates are accordingly specified within the following limits:

$${\mathbf{u}}^{\left(l\right)}=\left\{u_1^{\left(l\right)},u_2^{\left(l\right)}\right\}\in {\mathbb{R}}_{\{{\mathbf{u}}^{\left(l\right)},\tau \}}^2\cong (-\infty ,+\infty )\times (0,+\infty ).$$

In a similar way, we can define the domain of coordinate variation for the equation (50):

$$\mathbf{u}=\left\{u_1^{\left(1\right)},u_2^{\left(1\right)};u_1^{\left(2\right)},u_2^{\left(2\right)}\right\}\in {\mathbb{R}}_{\{\mathbf{u},\tau \}}^4\cong (-\infty ,+\infty )\times (0,+\infty )\times (-\infty ,+\infty )\times (0,+\infty ).$$

In the case where the fluctuation amplitudes are large $\{{ϵ}^{\left(r\right)},{ϵ}^{\left(i\right)}\}\ge 1$, then it is necessary to analyze in detail and determine the properties of the generating additional subspace in which the equation for the distribution of the fields of the environment (46) and (50) is determined.

Finally, it is important to note that equations (46) and (50) describe the distributions of the fields of the environment in the limit of statistical equilibrium. In particular, at realization of the first scenario the environment fields are generated by two independent random forces while at realization of the second scenario these random forces has the common nature. Below we will explore the problem associated with the second scenario.

4.3. Geometric and Topological Features of the Compactified Subspace

To determine the geometric and topological features of the additional subspace formed in the limit of statistical equilibrium, equation (50) can be conveniently represented in tensor form:

$${\partial }_t\mathcal{P}={\nabla }^2\mathcal{P}+K_0(v_1,v_3)\mathcal{P},{\nabla }^2={\displaystyle \frac{1}{\sqrt{\parallel g\parallel }}}\sum _{i,j=1}^4{\displaystyle \frac{\partial }{\partial v^i}}\left(\sqrt{\parallel g\parallel }{g}^{ij}{\displaystyle \frac{\partial }{\partial v^j}}\right),$$

(51)

where the following notations are made: the free term of the equation $K_0(v_1,v_3)=4(v_1+v_3)$, in addition:

$$u_1^{\left(1\right)}=v^1=v_1,u_2^{\left(1\right)}=v^2=v_2,u_1^{\left(2\right)}=v^3=v_3,u_2^{\left(2\right)}=v^4=v_4,$$

in addition, $\mathbf{u}\equiv \mathbf{v}=\{v_1,v_2,v_3,v_4\}$ and $g^{ij}(\mathbf{v},\tau )$ denotes a fourth-rank tensor.

By comparing equations (50) and (51) and equating the corresponding coefficients, one can write the following explicit form for the Laplace-Beltrami operator ${\nabla }^2$:

$$\begin{array}{c}\hfill {\nabla }^2=g^{11}{\displaystyle \frac{{\partial }^2}{\partial v_1^2}}+{\displaystyle \frac{1}{\sqrt{\parallel g\parallel }}}\left[{\displaystyle \frac{\partial }{\partial v_1}}\left(\sqrt{\parallel g\parallel }{g}^{11}\right)+{\displaystyle \frac{\partial }{\partial v_2}}\left(\sqrt{\parallel g\parallel }{g}^{21}\right)\right]{\displaystyle \frac{\partial }{\partial v_1}}+g^{12}{\displaystyle \frac{{\partial }^2}{\partial v_1\partial v_2}}\\ \hfill +g^{22}{\displaystyle \frac{{\partial }^2}{\partial v_2^2}}+{\displaystyle \frac{1}{\sqrt{\parallel g\parallel }}}\left[{\displaystyle \frac{\partial }{\partial v_2}}\left(\sqrt{\parallel g\parallel }{g}^{22}\right)+{\displaystyle \frac{\partial }{\partial v_1}}\left(\sqrt{\parallel g\parallel }{g}^{12}\right)\right]{\displaystyle \frac{\partial }{\partial v_2}}+g^{21}{\displaystyle \frac{{\partial }^2}{\partial v_2\partial v_1}}\\ \hfill +g^{33}{\displaystyle \frac{{\partial }^2}{\partial v_3^2}}+{\displaystyle \frac{1}{\sqrt{\parallel g\parallel }}}\left[{\displaystyle \frac{\partial }{\partial v_3}}\left(\sqrt{\parallel g\parallel }{g}^{33}\right)+{\displaystyle \frac{\partial }{\partial v_4}}\left(\sqrt{\parallel g\parallel }{g}^{43}\right)\right]{\displaystyle \frac{\partial }{\partial v_3}}+g^{34}{\displaystyle \frac{{\partial }^2}{\partial v_3\partial v_4}}\\ \hfill +g^{44}{\displaystyle \frac{{\partial }^2}{\partial v_4^2}}+{\displaystyle \frac{1}{\sqrt{\parallel g\parallel }}}\left[{\displaystyle \frac{\partial }{\partial v_4}}\left(\sqrt{\parallel g\parallel }{g}^{44}\right)+{\displaystyle \frac{\partial }{\partial v_3}}\left(\sqrt{\parallel g\parallel }{g}^{34}\right)\right]{\displaystyle \frac{\partial }{\partial v_4}}+g^{43}{\displaystyle \frac{{\partial }^2}{\partial v_4\partial v_3}},\end{array}$$

(52)

where $g^{31}=g^{41}=g^{32}=g^{42}=g^{13}=g^{14}=g^{23}=g^{24}=0$, in addition;

$$g^{11}=g^{33}={ϵ}^{\left(r\right)},g^{22}=g^{44}={ϵ}^{\left(i\right)},g=(g^{11}{g}^{22}-g^{12}{g}^{21})(g^{33}{g}^{44}-g^{34}{g}^{43}).$$

(53)

Comparing equations (50) and (51) it is easy to see that the following equalities must hold: $g^{12}=-g^{21}$ and $g^{34}=-g^{43}$. The latter is equivalent to the fact that the additional 4D-subspace (manifold) ${\mathsf{\Xi }}_{\{\mathbf{v},\tau \}}^{\left(4\right)}$, arising after the relaxation of the photon+medium system, has the properties of non-commutative geometry.

In the case under consideration, the metric of the space $g^{ij}(\mathbf{u},\tau )$ has the form:

$$g^{ij}(\mathbf{u},\tau )\equiv g^{ij}(\mathbf{v},\tau )=\left(\begin{array}{cccc}{g}^{11}& g^{12}& 0& 0\\ g^{21}& g^{22}& 0& 0\\ 0& 0& g^{33}& g^{34}\\ 0& 0& g^{43}& g^{44}\end{array}\right).$$

(54)

Taking into account the form of the metric tensor $g^{ij}(\mathbf{u},\tau )$, we can represent the solution of equation (50) in factored form:

$$\mathcal{P}(\mathbf{v},\tau )=\mathcal{P}(\mathbf{u},\tau )=\prod _{l=1}^2{P}_l({\mathbf{u}}^{\left(l\right)},\tau ),$$

(55)

where $P_l({\mathbf{u}}^{\left(l\right)},\tau )$ is a solution to the equation (50).

From the above it follows that the additional 4D manifold can be represented as a direct product of two 2D submanifolds:

$${\mathsf{\Xi }}_{\{\mathbf{v},\tau \}}^{\left(4\right)}\equiv {\mathsf{\Xi }}_{\{\mathbf{u},\tau \}}^{\left(4\right)}\cong {\mathsf{\Xi }}_{\{{\mathbf{u}}^{\left(1\right)},\tau \}}^{\left(2\right)}⨂{\mathsf{\Xi }}_{\{{\mathbf{u}}^{\left(2\right)},\tau \}}^{\left(2\right)}.$$

Since the equation (46) for the probability distribution $P_l({\mathbf{u}}^{\left(l\right)},\tau )$ is defined on the manifold ${\mathsf{\Xi }}_{\{{\mathbf{u}}^{\left(l\right)},\tau \}}^{\left(2\right)}$, the main task now will be to determine its geometric and topological properties. In other words, we need to derive an equation to determine the antisymmetric element of the metric tensor $y^{\left(1\right)}=g^{12}({\mathbf{u}}^{\left(1\right)},\tau )\equiv g^{12}(u_1^{\left(1\right)},u_2^{\left(1\right)},\tau )$. Note that the equation for $y^{\left(2\right)}=g^{34}({\mathbf{u}}^{\left(2\right)},\tau )\equiv g^{34}(u_1^{\left(2\right)},u_2^{\left(2\right)},\tau )$ will be derived in a similar manner.

Following the works [32], we can compare the operators (47) and (52) write the following two partial differential equations for definition $y=y^{\left(1\right)}$ (the index in round brackets is omitted everywhere):

$$\begin{array}{c}\hfill {ϵ}^{\left(r\right)}\chi {\displaystyle \frac{\partial y}{\partial u_1}}-\left(1+y\chi \right){\displaystyle \frac{\partial y}{\partial u_2}}=k_1(u_1,u_2,\tau ),\\ \hfill {ϵ}^{\left(i\right)}\chi {\displaystyle \frac{\partial y}{\partial u_2}}+\left(1+y\chi \right){\displaystyle \frac{\partial y}{\partial u_1}}=k_2(u_1,u_2,\tau ),\end{array}$$

(56)

where $\chi =y/(a+y^2)$ and $a={ϵ}^{\left(i\right)}{ϵ}^{\left(r\right)}$.

As follows from (56), the system of equations is overdetermined with respect to the sought function “y”. In connection with this, a reasonable question arises: under what conditions are these equations equivalent?

Using (56), it is easy to find expressions for two different derivatives of the off-diagonal component of the metric tensor:

$$y_1={\displaystyle \frac{\partial y}{\partial u_1}}={\displaystyle \frac{{ϵ}^{\left(i\right)}{k}_1\chi +k_2(1+y\chi )}{a{\chi }^2+{(1+y\chi )}^2}},y_2={\displaystyle \frac{\partial y}{\partial u_2}}={\displaystyle \frac{{ϵ}^{\left(r\right)}{k}_2\chi -k_1(1+y\chi )}{a{\chi }^2+{(1+y\chi )}^2}}.$$

(57)

Now we can calculate expressions for the mixed second derivatives of an off-diagonal element of the metric tensor using the expressions (57):

$$\begin{array}{c}\hfill y_{12}={\displaystyle \frac{{\partial }^{2}y}{\partial u_2\partial u_1}}={\displaystyle \frac{{ϵ}^{\left(i\right)}(k_{1;2}\chi +k_1{\chi }_{;2})+k_{2;2}(1+y\chi )+k_2(y_2\chi +y{\chi }_{;2})}{a{\chi }^2+{(1+y\chi )}^2}}\\ \hfill -2{\displaystyle \frac{{ϵ}^{\left(i\right)}{k}_1\chi +k_2(1+y\chi )}{{[a{\chi }^2+{(1+y\chi )}^2]}^2}}\left[a\chi {\chi }_{;2}+(1+y\chi )(y_2\chi +y{\chi }_{;2})\right],\\ \hfill y_{21}={\displaystyle \frac{{\partial }^{2}y}{\partial u_1\partial u_2}}={\displaystyle \frac{{ϵ}^{\left(r\right)}(k_{2;1}\chi +k_2{\chi }_{;1})-k_{1;1}(1+y\chi )-k_1(y_1\chi +y{\chi }_{;1})}{a{\chi }^2+{(1+y\chi )}^2}}\\ \hfill -2{\displaystyle \frac{{ϵ}^{\left(r\right)}{k}_2\chi -k_1(1+y\chi )}{{[a{\chi }^2+{(1+y\chi )}^2]}^2}}\left[a\chi {\chi }_{;1}+(1+y\chi )(y_1\chi +y{\chi }_{;1})\right],\end{array}$$

(58)

where ${\chi }_{;\nu }=\partial \chi /\partial u_{\nu }$ and $k_{\nu ;\mu }=\partial k_{\nu }/\partial u_{\mu },(\nu ,\mu =1,2).$

Based on the basic requirement of mathematical analysis, for an analytic function there must be symmetry of mixed second derivatives (Schwarz’s theorem, see [33]):

$$y_{12}={\displaystyle \frac{{\partial }^{2}y}{\partial u_1\partial u_2}}=y_{21}={\displaystyle \frac{{\partial }^{2}y}{\partial u_2\partial u_1}},$$

(59)

which is a necessary condition for a twice continuously differentiable function. In the context of partial differential equations, it is called the Schwarz integrability condition. Now if we write this equality (59) explicitly, taking into account the expressions (58), it will look like this:

$$2{\displaystyle \frac{{ϵ}^{\left(r\right)}{k}_2\chi -k_1(1+y\chi )}{a{\chi }^2+{(1+y\chi )}^2}}\left[a\chi {\chi }_{;1}+(1+y\chi )(y_1\chi +y{\chi }_{;1})\right]+(k_1+k_{2}y){\chi }_{;2}=$$

$$2{\displaystyle \frac{{ϵ}^{\left(i\right)}{k}_1\chi +k_2(1+y\chi )}{a{\chi }^2+{(1+y\chi )}^2}}\left[a\chi {\chi }_{;2}+(1+y\chi )(y_2\chi +y{\chi }_{;2})\right]+(k_2-k_{1}y){\chi }_{;1}$$

$$-(k_{1;1}+k_{2;2})(1+y\chi )-(k_1{y}_1+k_2{y}_2+({ϵ}^{\left(i\right)}{k}_{1;2}-{ϵ}^{\left(r\right)}{k}_{2;1})\chi =0.$$

(60)

Finally, taking into account the expressions (57) from equation (60), we obtain the following algebraic equation of the 4th degree with coefficients depending on the coordinates of the additional subspace:

$$\sum _{n=0}^4{A}_n(u_1,u_2,\tau )y^n=0,$$

(61)

where the coefficients of the algebraic equation $A_n(u_1,u_2,\tau )$ are defined by the expressions:

$$\begin{array}{c}\hfill A_0=a\left[4a{u}_1-4{ϵ}^{\left(r\right)}{u}_1^2{u}_2^2-\mathcal{X}(u_1,u_2,\tau )\right],A_1=-2a{u}_2\left[{ϵ}^{\left(r\right)}+{ϵ}^{\left(i\right)}\right],\\ \hfill A_2=\left[24a{u}_1+8{ϵ}^{\left(r\right)}{u}_1^2{u}_2^2+2\mathcal{X}(u_1,u_2,\tau )\right],A_3=4A_1,A_4=32u_1,\end{array}$$

(62)

where $\mathcal{X}(u_1,u_2,\tau )={ϵ}^{\left(i\right)}{\left[u_1^2-u_2^2+{\omega }_1^2\left(\tau \right)\right]}^2$.

It should be noted that a similar algebraic equation for determining the off-diagonal term is also obtained for the case of the manifold ${\mathsf{\Xi }}_{\{{\mathbf{u}}^{\left(2\right)},\tau \}}^{\left(2\right)}$ with the only difference that in the coefficients of the equation (62) instead of the frequency ${\omega }_1\left(\tau \right)$, the frequency ${\omega }_2\left(\tau \right)$ appears.

Now, if we specify the necessary parameters ${ϵ}^{\left(i\right)},{ϵ}^{\left(r\right)}$ and ${\omega }_l\left(\tau \right)$ included in the algebraic equation (61), this will allow us to explicitly solve the equation and visualize it as a 2D-surface (submanifold ${\mathsf{\Xi }}_{\left\{{\mathbf{u}}^{\left(l\right)}\right\}}^{\left(2\right)})$. As the calculations performed show (see Figure 1), depending on the set of problem parameters, these submanifolds can possess nontrivial geometric and topological properties. In particular, they can have different Betti numbers, which characterize topological invariants, representing the rank of their, in this case, 2D homology group, and counting the number of 2D “holes".

5. The Mathematical Expectation of the Wave Function of a Single Photon

5.1. The Measure of the Functional Space

To continue the analytical constructions of the theory, it is necessary to determine the distance between functions in the space ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$, or more precisely, to determine the measure of the functional space (see [30,32]).

Let the conditional probability $P(\mathbf{u},\tau |{\mathbf{u}}^{\prime },{\tau }^{\prime })$ satisfy the following limiting condition:

$$\begin{array}{c}\hfill \underset{\tau \to {\tau }^{\prime }}{lim}P(\mathbf{u},\tau |{\mathbf{u}}^{\prime },{\tau }^{\prime })=\delta (\mathbf{u}-{\mathbf{u}}^{\prime }),\tau ={\tau }^{\prime }+\mathsf{\Delta }\tau .\end{array}$$

(63)

The latter means that for small intervals, i.e. for $\mathsf{\Delta }\tau =\tau -{\tau }^{\prime }\ll 1$, we can represent the solution of equations (46)-(47) in the form:

$$P(\mathbf{u},\tau |{\mathbf{u}}^{\prime },{\tau }^{\prime })={\displaystyle \frac{1}{2\pi \sqrt{|det\mathsf{{\rm Y}}|}\mathsf{\Delta }\tau }}exp\left\{-{\displaystyle \frac{{\left[\mathbf{u}-{\mathbf{u}}^{\prime }-\mathbf{K}(\mathbf{u},\tau )\mathsf{\Delta }\tau \right]}^T{\mathsf{{\rm Y}}}^{-1}\left[\mathbf{u}-{\mathbf{u}}^{\prime }-\mathbf{K}(\mathbf{u},\tau )\mathsf{\Delta }\tau \right]}{2\mathsf{\Delta }\tau }}\right\},$$

(64)

where $\mathsf{{\rm Y}}$ is the second-rank matrix with elements; ${ϵ}_{11}={ϵ}^{\left(r\right)},{ϵ}_{22}={ϵ}^{\left(i\right)}$ and ${ϵ}_{12}={ϵ}_{21}=0$, while ${[···]}^T$ denotes a vector transposition.

Using (64), we can write an explicit expression for the probability distribution of the transition between two very close points:

$$\begin{array}{c}\hfill P(u_1,u_2,\tau |u_1^{\prime },u_2^{\prime },{\tau }^{\prime })={\displaystyle \frac{1}{2\pi \sqrt{{ϵ}^{\left(r\right)}{ϵ}^{\left(i\right)}}\mathsf{\Delta }\tau }}\times \\ \hfill exp\left\{-{\displaystyle \frac{{\left[u_1-u_1^{\prime }-k_1(u_1,u_2,\tau )\mathsf{\Delta }\tau \right]}^2}{2{ϵ}^{\left(r\right)}\mathsf{\Delta }\tau }}-{\displaystyle \frac{{\left[u_2-u_2^{\prime }-k_2(u_1,u_2)\mathsf{\Delta }\tau \right]}^2}{2{ϵ}^{\left(i\right)}\mathsf{\Delta }\tau }}\right\},\end{array}$$

(65)

where the coefficients $k_1$ and $k_2$ are defined from expression (46).

In the case where dissipation in the medium is absent or extremely small, for example when light passes through a vacuum, it is necessary to set ${ϵ}^{\left(i\right)}=0$, and accordingly the distribution (65) takes the following form:

$$\begin{array}{c}\hfill P(u_1,u_2,\tau |u_1^{\prime },u_2^{\prime },{\tau }^{\prime })={\displaystyle \frac{1}{2\pi \sqrt{{ϵ}^{\left(r\right)}\mathsf{\Delta }\tau }}}\times \\ \hfill exp\left\{-{\displaystyle \frac{{\left[u_1-u_1^{\prime }-k_1(u_1,u_2,\tau )\mathsf{\Delta }\tau \right]}^2}{2{ϵ}^{\left(r\right)}\mathsf{\Delta }\tau }}\right\}\delta \left[u_2-u_2^{\prime }+k_2(u_1,u_2)\right].\end{array}$$

(66)

Finally, taking into account the above, we can write down the probability of transition between any two points along any predetermined trajectory in an uncountable-dimensional functional space ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$:

$$\begin{array}{c}\hfill D\mu \left(\mathbf{u}\right)=d\mu \left({\mathbf{u}}_0\right)\underset{N\to \infty }{lim}\{{\left({\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{N/\tau }{\sqrt{{ϵ}^{\left(r\right)}{ϵ}^{\left(i\right)}}}}\right)}^N\prod _{l=0}^{N}d{u}_{1(l+1)}d{u}_{2(l+1)}exp[-{\displaystyle \frac{N/\tau }{2{ϵ}^{\left(r\right)}}}\times \\ \hfill {\left(u_{1(l+1)}-u_{1\left(l\right)}-k_{1(l+1)}{\displaystyle \frac{{\tau }_{l+1}}{N}}\right)}^2-{\displaystyle \frac{N/\tau }{2{ϵ}^{\left(i\right)}}}{\left(u_{2(l+1)}-u_{2\left(l\right)}-k_{2(l+1)}{\displaystyle \frac{{\tau }_{l+1}}{N}}\right)}^2\left]\right\},\end{array}$$

(67)

where $d\mu \left({\mathbf{u}}_0\right)=\delta (u_1-u_{1\left(0\right)})\delta (u_2-u_{2\left(0\right)})d{u}_{1}d{u}_2$ denotes the measure of the initial distribution; in addition, the following notations are made:

$$u_{1\left(l\right)}=u_1\left({\tau }_l\right),u_{2\left(l\right)}=u_2\left({\tau }_l\right),k_{1\left(l\right)}=k_1(u_{1\left(l\right)},u_{2\left(l\right)},{\tau }_l),k_{2\left(l\right)}=k_2(u_{1\left(l\right)},u_{2\left(l\right)},{\tau }_l).$$

Note that the expression (67) is essentially the Fokker-Planck measure of the functional space ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$, both in the case when the fluctuation powers are small and when they are large and the corresponding geometry in the limit of statistical equilibrium is non-commutative. Recall that in the limit of statistical equilibrium, the functional space ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$ is compactified into 4D additional subspace ${\mathbb{R}}_{\left\{\mathbf{u}\right\}}^4$ or ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right\}}^{\left(4\right)}\left(\tau \right)$, respectively, depending on the values of the fluctuation powers ${ϵ}^{\left(i\right)}$ and ${ϵ}^{\left(r\right)}$.

5.2. The Wave Function of a Single Photon in the Limit of Statistical Equilibrium

Now that we have constructed the Fokker-Planck measure of the functional space (67), in the limit of statistical equilibrium, it is easy to calculate the mathematical expectation of the photon wave function in the plane perpendicular to the direction of photon propagation.

Definition 2. The mathematical expectation of a random variable will be by the following integral over the functional space ${\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$:

$$\mathbb{E}\left[\mathcal{G}(q_1,q_2,\tau ;\left\{f\right\})\right]={\int }_{{\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}}\mathcal{G}(q_1,q_2,\tau ;\left\{f\right\})D\mu \left(\mathbf{u}\right).$$

(68)

where $\mathcal{G}(q_1,q_2,\tau ;\left\{f\right\})$ is a random variable, and $\mathbb{E}\left[···\right]$ denotes the mathematical expectation.

The solution of the Langevin-Schrödinger equation (35) for a photon moving in a nanowaveguide with random influences can be represented as:

$$\mathcal{F}\left(\mathbf{n}\right|\mathbf{q},\tau ;\left\{f\right\})=\prod _{l=1}^2{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\}).$$

(69)

Taking into account (68) and the fact that the additional subspace in general is reduced to the direct product of two 2D subspaces (see subsections 4.2 and 4.3), it is easy to prove that the mathematical expectation of the photon wave function can be written in factorized form:

$$\mathbb{E}\left[\mathcal{F}\left(\mathbf{n}\right|\mathbf{q},\tau ;\left\{f\right\})\right]=\prod _{l=1}^2\mathbb{E}\left[{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})\right],$$

(70)

where

$$\mathbb{E}\left[{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})\right]={\int }_{{\mathsf{\Xi }}_{\left\{{\mathbf{u}}^{\left(l\right)}\left(\tau \right)\right\}}}{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})D\mu \left({\mathbf{u}}^{\left(l\right)}\right).$$

(71)

Before we start calculating the functional integral (71), it is important to write down the explicit form of the photon random wave function for an arbitrary quantum state. In particular, taking into account (42) and (44), we can write:

$${\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})=\sqrt{{\displaystyle \frac{1}{2^{n_l}{n}_l!\sqrt{\pi }}}}exp\left\{{\displaystyle \frac{1}{2}}\left[i{\mathsf{\Pi }}^{\prime }\left(\tau \right)q_l^2-\mathsf{\Pi }\left(\tau \right)\right]\right\}{H}_{n_l}\left(q_l{e}^{-n_l\mathsf{\Pi }\left(\tau \right)}\right).$$

(72)

where ${\mathsf{\Pi }}_l\left(\tau \right)={\int }_0^{\tau }\left[u_1^{\left(l\right)}\left({\tau }^{\prime }\right)+i{u}_2^{\left(l\right)}\left({\tau }^{\prime }\right)\right]d{\tau }^{\prime }$ (see (42), definition of the function $\mathsf{\Pi }\left(\tau \right)$), in addition, ${\mathsf{\Pi }}^{\prime }\left(\tau \right)=d\mathsf{\Pi }\left(\tau \right)/d\tau$. Below, to simplify notation, the index “l” above the variables $u_1$ and $u_2$ will be omitted.

From expression (72) it is easy to find the wave functions of the first two quantum states:

$$\begin{array}{c}\hfill {\mathcal{F}}_l\left(0\right|q_l,\tau ;\left\{f_l\right\})=\sqrt{{\displaystyle \frac{1}{\sqrt{\pi }}}}exp\left\{-{\displaystyle \frac{1}{2}}\underset{0}{\overset{\tau }{\int }}\right(u_1\left({\tau }^{\prime }\right)+i{u}_2\left({\tau }^{\prime }\right))d{\tau }^{\prime }+{\displaystyle \frac{1}{2}}(i{u}_1-u_2)q_l^2\},\\ \hfill {\mathcal{F}}_l\left(1\right|q_l,\tau ;\left\{f_l\right\})=q_l\sqrt{{\displaystyle \frac{1}{2\sqrt{\pi }}}}exp\left\{-{\displaystyle \frac{3}{2}}\underset{0}{\overset{\tau }{\int }}\left(u_1\left({\tau }^{\prime }\right)+i{u}_2\left({\tau }^{\prime }\right)\right)d{\tau }^{\prime }+{\displaystyle \frac{1}{2}}(i{u}_1-u_2)q_l^2\right\}.\end{array}$$

(73)

Substituting (73) into (71) and using the generalized Feynman–Kac theorem [29], we can calculate the corresponding functional integrals and find the ME of the photon wave function in various quantum states, representing them as double integrals:

$$\begin{array}{c}\hfill {\overline{\mathcal{F}}}_l\left(0\right|q_l,\tau ;\left\{f_l\right\})=\sqrt{{\displaystyle \frac{1}{\sqrt{\pi }}}}\underset{-\infty }{\overset{+\infty }{\int }}d{u}_1\underset{0}{\overset{+\infty }{\int }}d{u}_2{Q}_l(0;u_1,u_2,\tau )e^{{\mathcal{A}}_l(u_1,u_2;q_l)},\\ \hfill {\overline{\mathcal{F}}}_l\left(1\right|q_l,\tau ;\left\{f_l\right\})=q_l\sqrt{{\displaystyle \frac{1}{2\sqrt{\pi }}}}\underset{-\infty }{\overset{+\infty }{\int }}d{u}_1\underset{0}{\overset{+\infty }{\int }}d{u}_2{Q}_l(1;u_1,u_2,\tau )e^{{\mathcal{A}}_l(u_1,u_2;q_l)},\end{array}$$

(74)

where ${\mathcal{A}}_l(u_1,u_2;q_l)=(i{u}_1-u_2)q_l^2/2$ and ${\overline{\mathcal{F}}}_l\left(n\right|q_l,\tau ;\left\{f_l\right\})=\mathbb{E}\left[{\mathcal{F}}_l\left(n\right|q_l,\tau ;\left\{f_l\right\})\right]$ denotes the ME of the wave function.

As for the integrand $Q_l(n;u_1,u_2,\tau )$, it is a solution to the following second-order complex PDE:

$${\partial }_{\tau }{Q}_l(n;u_1,u_2,\tau )=\left\{{\widehat{L}}_l-(n+1/2)(u_1+i{u}_2)\right\}{Q}_l(n;u_1,u_2,\tau ),n=0,1,2,...,$$

(75)

where the operator ${\widehat{L}}_l$ is defined by the formula (47).

By representing the solution of equation (75) as a sum:

$$Q_l(n;u_1,u_2,\tau )=Q_l^{\left(r\right)}(n;u_1,u_2,\tau )+i{Q}_l^{\left(i\right)}(n;u_1,u_2,\tau ),$$

(76)

we can obtain the following system of PDEs:

$$\{\begin{array}{c}{\partial }_{\tau }{Q}_l^{\left(r\right)}=\left\{{\widehat{L}}_l-(n+1/2)u_1\right\}{Q}_l^{\left(r\right)}+(n+1/2)u_2{Q}_l^{\left(i\right)},\hfill \\ {\partial }_{\tau }{Q}_l^{\left(i\right)}=\left\{{\widehat{L}}_l-(n+1/2)u_1\right\}{Q}_l^{\left(i\right)}-(n+1/2)u_2{Q}_l^{\left(r\right)}.\hfill \end{array}$$

(77)

The functions $Q_l^{\left(r\right)}$ and $Q_l^{\left(i\right)}$ can be given a probabilistic meaning and normalized to unity. In particular, in what follows we will use the normalized distributions ${\overline{Q}}_l^{\left(r\right)}$ and ${\overline{Q}}_l^{\left(i\right)}$, which are defined by the expressions:

$$\begin{array}{c}\hfill {\overline{Q}}_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau )=Q_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau )/J_l(n;\tau ),\upsilon =i,r,\end{array}$$

(78)

where $J_l(n;\tau )$ is the probability normalization coefficient, determined by the expression:

$$J_l(n;\tau )={\int }_0^{+\infty }{\int }_{-\infty }^{+\infty }\left[Q_l^{\left(r\right)}(n;u_1,u_2,\tau )+Q_l^{\left(i\right)}(n;u_1,u_2,\tau )\right]d{u}_{1}d{u}_2.$$

It is obvious that after normalization of the distributions ${\overline{Q}}_l^{\left(r\right)}(u_1,u_2,\tau )$ and ${\overline{Q}}_l^{\left(i\right)}(u_1,u_2,\tau )$ the probability in JS is preserved and, accordingly, the equality is satisfied:

$${\int }_0^{+\infty }{\int }_{-\infty }^{+\infty }\left[{\overline{Q}}_l^{\left(r\right)}(n;u_1,u_2,\tau )+{\overline{Q}}_l^{\left(i\right)}(n;u_1,u_2,\tau )\right]d{u}_{1}d{u}_2=1.$$

(79)

Now, returning to the ME of wave function (74), it is necessary to note that it, in essence, describes a “closed", self-consistent, JS “photon + environment” and, therefore, must be normalized to unity. In other words, the normalized ME of the photon’s wave function will be represented as:

$${\underline{\overline{\mathcal{F}}}}_l\left(n\right|q_l,\tau ;\left\{f_l\right\})={\overline{\mathcal{F}}}_l\left(n\right|q_l,\tau ;\left\{f_l\right\})/{\mathcal{J}}_l(n,\tau ),$$

(80)

where ${\mathcal{J}}_l(n,\tau )={\int }_{-\infty }^{+\infty }{\overline{\mathcal{F}}}_l\left(n\right|q_l,\tau ;\left\{f_l\right\}){\overline{\mathcal{F}}}_l^{*}\left(n\right|q_l,\tau ;\left\{f_l\right\})d{q}_l$ is the normalization factor.

To solve the SDEs system, it is convenient to formulate the Neumann initial-boundary value problem. In particular, we will represent the initial distributions as the product of two Dirac delta-functions over the perpendicular coordinates $u_1$ and $u_2$, accordingly:

$${\overline{Q}}_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau ){|}_{\tau =0}=\prod _{j=1}^2\delta (u_j-u_{0j}),\upsilon =i,r.$$

(81)

Recall that the product of two Dirac delta-functions should be understood in such a way that when $u_1=u_{01}$ and $u_2=u_{02}$ then it is infinite, in all other cases it is zero.

Based on the conditions of the problem under consideration, we can write the normalization condition for the total probability as follows:

$${\int }_0^{+\infty }{\int }_{-\infty }^{+\infty }\left[p{\overline{Q}}_l^{\left(r\right)}(n;u_1,u_2,{\tau }_0)+(1-p){\overline{Q}}_l^{\left(i\right)}(n;u_1,u_2,{\tau }_0)\right]d{u}_{1}d{u}_2=1/2,$$

(82)

where $p\in [0,1]$ denotes the weight of the distribution associated with elastic quantum processes accompanying a photon in the environment, and $(1-p)$ is the weight of the distribution associated with inelastic quantum processes. Thus, the probability distributions ${\overline{Q}}_l^{\left(r\right)}(n;u_1,u_2,\tau )$ and ${\overline{Q}}_l^{\left(i\right)}(n;u_1,u_2,\tau )$ are combined into a common probability, starting from the initial moment of evolution time $ta{u}_0$ for a single photon. If we check the probability normalization condition (79) for the initial time $\tau =0$, then it is easy to verify that the initial distributions (81) satisfy the normalization to unity, since the integration over the coordinate $u_2$ is carried out over the segment of the numerical axis $[0,\infty )$, i.e. ${\int }_0^{\infty }\delta \left(u_2\right)d{u}_2=1/2.$

Below there are visualizations (see graphs in Figure 3 and Figure 4) of calculations of the evolution of normalized distributions of the environmental fields ${\overline{Q}}_1^{\left(r\right)}$ and ${\overline{Q}}_1^{\left(i\right)}$ in the ground and first excited states. As can be seen from the graphs, in the process of evolution of QS under certain conditions of the problem, a fairly rapid transfer of probability occurs from the distribution ${\overline{Q}}_1^{\left(i\right)}$ to ${\overline{Q}}_1^{\left(r\right)}$.

As for the boundary conditions for solving the system of partial differential equations (77), they should be specified on the axis $u_1$, the natural boundary of the problem, as follows:

$${\displaystyle \frac{\partial {\overline{Q}}_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau )}{\partial u_2}}{|}_{u_2=0}={\overline{Q}}_l^{\left(\upsilon \right)}(n;u_1,0,\tau )={\mathbb{Q}}_l^{\left(\upsilon \right)}(n;u_1,\tau ).$$

(83)

Finally, using expressions (70) and (74), we can construct the mathematical expectations of the first four quantum states of a photon in a nanowaveguide:

$$\begin{array}{c}\hfill \mathbb{E}\left[\mathcal{F}(0,0|\mathbf{q},\tau ;\left\{f\right\})\right]={\overline{\mathcal{F}}}_1\left(0\right|q_1,\tau ;\left\{f_1\right\})·{\overline{\mathcal{F}}}_2\left(0\right|q_2,\tau ;\left\{f_2\right\}),\\ \hfill \mathbb{E}\left[\mathcal{F}(0,1|\mathbf{q},\tau ;\left\{f\right\})\right]={\overline{\mathcal{F}}}_1\left(0\right|q_l,\tau ;\left\{f_1\right\})·{\overline{\mathcal{F}}}_2\left(1\right|q_2,\tau ;\left\{f_2\right\}),\\ \hfill \mathbb{E}\left[\mathcal{F}(1,0|\mathbf{q},\tau ;\left\{f\right\})\right]={\overline{\mathcal{F}}}_1\left(1\right|q_1,\tau ;\left\{f_1\right\})·{\overline{\mathcal{F}}}_2\left(0\right|q_2,\tau ;\left\{f_2\right\}),\\ \hfill \mathbb{E}\left[\mathcal{F}(1,1|\mathbf{q},\tau ;\left\{f\right\})\right]={\overline{\mathcal{F}}}_1\left(1\right|q_1,\tau ;\left\{f_1\right\})·{\overline{\mathcal{F}}}_2\left(1\right|q_2,\tau ;\left\{f_2\right\}).\end{array}$$

(84)

Recall that the quantum states of the photon $\mathbb{E}\left[\mathcal{F}(0,1|\mathbf{q},\tau ;\left\{f\right\})\right]$ and $\mathbb{E}\left[\mathcal{F}(1,0|\mathbf{q},\tau ;\left\{f\right\})\right]$ are, in the general case, obviously different and do not coincide.

Based on the obtained formulas for the mathematical expectation of the wave function of a photon in various quantum states (84), one can write the corresponding distributions of the photon energy density in 3D space as a function of time:

$$\begin{array}{c}\hfill {\mathcal{W}}_{nm}={\left|\mathbb{E}\left[\mathcal{F}(n,m|\mathbf{q},\tau ;\left\{f\right\})\right]\right|}^2={\left|{\overline{\mathcal{F}}}_1\left(n\right|q_1,\tau ;\left\{f_1\right\})\right|}^2·{\left|{\overline{\mathcal{F}}}_2\left(m\right|q_2,\tau ;\left\{f_2\right\})\right|}^2,\end{array}$$

(85)

where n and m are vibrational quantum numbers. In order for ${\mathcal{W}}_{nm}$ to acquire the meaning of the probability distribution of the energy of electromagnetic fields in a single photon, it is necessary to replace the wave functions in the expressions (85) on normalized functions (80).

Below in subsection 9.2, for one particular set of initial data, see Table 1, the spatial distribution of the electromagnetic fields energy of one photon for the first four quantum states is calculated and visualized, see Figure 5.

5.3. Reduced Density Matrix: Quantized Small Environment

Now it is important to construct the elements of the reduced density matrix, which play an important role in calculating various QS parameters. Taking into account (72), we can write the following expression for the element of the stochastic density matrix:

$${\rho }_{nm}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})=\mathcal{F}(n,m|\mathbf{q},\tau ;\left\{f\right\})·{\mathcal{F}}^{*}(n,m|{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\}).$$

Definition 3. An element of the reduced density matrix (RDM) ${\rho }_{nm}(\mathbf{q},\tau )$ will be called the mathematical expectation of an element of a random density matrix ${\rho }_{nm}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };{\left\{f\right\}}^{\prime })$, obtained after functional integration:

$${\rho }_{nm}(\mathbf{q},\tau )=T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}\left[{\rho }_{nm}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})\right]{|}_{\mathbf{q}={\mathbf{q}}^{\prime };\tau ={\tau }^{\prime }},$$

(86)

where the operation, $T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}$ denotes the functional integration over the environment fields:

$$\begin{array}{c}\hfill T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}\left[\rho (\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})\right]={\int }_{{\mathsf{\Xi }}_{\left\{\mathbf{u}\right(\tau \left)\right\}}}\rho (\mathbf{q},\tau ;\left\{f\right\}|{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})D\left\{\mathbf{u}\left(\tau \right)\right\}.\end{array}$$

(87)

By substituting the corresponding random wave function from (72) into the expression (86)-(87) and performing functional integration, one can find the corresponding element of the RDM. In particular, after simple calculations for the first four elements of the RMD, the following expressions can easily be found:

$$\begin{array}{c}\hfill {\varrho }_{00}(\mathbf{q},\tau )=\prod _{l=1}^2{\varrho }_0^{\left(l\right)}={\displaystyle \frac{1}{\pi }}\prod _{l=1}^2{\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_l(0;u_1,u_2,\tau )e^{-u_2{q}_l^2},\\ \hfill {\varrho }_{11}(\mathbf{q},\tau )=\prod _{l=1}^2{\varrho }_1^{\left(l\right)}={\displaystyle \frac{q_1^2{q}_2^2}{4\pi }}\prod _{l=1}^2{\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_l(1;u_1,u_2,\tau )e^{-u_2{q}_l^2},\\ \hfill {\varrho }_{01}(\mathbf{q},\tau )={\rho }_0^{\left(1\right)}{\rho }_1^{\left(2\right)}={\displaystyle \frac{q_2^2}{2\pi }}{\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_1(0;u_1,u_2,\tau )e^{-u_2{q}_1^2}\times \\ \hfill {\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_2(1;u_1,u_2,\tau )e^{-u_2{q}_2^2},\\ \hfill {\varrho }_{10}(\mathbf{q},\tau )={\rho }_1^{\left(1\right)}{\rho }_0^{\left(2\right)}={\displaystyle \frac{q_1^2}{2\pi }}{\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_1(1;u_1,u_2,\tau )e^{-u_2{q}_1^2}\times \\ \hfill {\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_2(0;u_1,u_2,\tau )e^{-u_2{q}_2^2}.\end{array}$$

(88)

It is not difficult to verify from (88) that the matrix elements ${\varrho }_{01}(\mathbf{q},\tau )$ and ${\varrho }_{10}(\mathbf{q},\tau )$ are in general asymmetric.

As for the integrand ${\mathcal{P}}_l(n;u_1,u_2,\tau )$ in equations (88), it represents the probability density of the environmental fields after self-organization of JS. It is easy to verify that this function satisfies the following PDE:

$${\partial }_{\tau }{\mathcal{P}}_l\left(n\right|u_1,u_2,\tau )=\left\{{\widehat{L}}_l-(2n+1)u_1\right\}{\mathcal{P}}_l\left(n\right|u_1,u_2,\tau ),n=0,1,2...,$$

(89)

where the operator ${\widehat{L}}_l$ has the form (47).

To solve the equation (89), one can formulate the Neumann initial-boundary value problem similarly to Section 5.2, only in this case for one equation:

$${\mathcal{P}}_l(n;u_1,u_2,\tau ){|}_{\tau =0}=\prod _{j=1}^2\delta (u_j-u_{0j}),{\displaystyle \frac{\partial {\mathcal{P}}_l(n;u_1,u_2,\tau )}{\partial u_2}}{|}_{u_2=0}={\mathbb{P}}_l(n;u_1,\tau ).$$

(90)

Note that the distribution ${\mathcal{P}}_l(n;u_1,u_2,\tau )$ finally takes on the meaning of a probability density after normalization:

$${\overline{\mathcal{P}}}_l(n;u_1,u_2,\tau )={\mathcal{P}}_l(n;u_1,u_2,\tau )/{\mathcal{J}}_l^0(n,\tau ),$$

where ${\mathcal{J}}_l^0(n,\tau )={\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_l(n;u_1,u_2,\tau )$ denotes the normalization coefficient.

Finally, we note a very important fact: as follows from equation (89), the fields of the environment form the so-called small environment (SE), which, as a result of interaction with the QS, itself becomes quantized.

Below in subsection 9.3, calculations of the evolution of the elements of the reduced density matrix and their visualization in the form of 3D graphs (see Figure 6) are given for a specific set of parameters of the problem under consideration.

6. Entropy of a Single Photon Propagating in a Medium

The propagation of a single photon in a nanowaveguide with random influences is a typical irreversible quantum process. Let us recall that such processes are usually studied within the framework of the representation of a non-stationary density matrix, using the Liouville – von Neumann equation [34] or its various generalizations [35,36]. However, there is an important limitation that makes the application of these representations in the case under consideration unsuitable. It should be noted that standard representations of the density matrix take into account the influence of the environment on a quantum system, while the reverse influence, namely the influence on the environment, is not taken into account or is taken into account insufficiently consistently, which in some cases leads to a noticeable loss of information of JS. In other words, these difficulties are insurmountable within the framework of the concept of open quantum systems, since they in any case allow approximate approaches in deriving the master equation.

In the problem under consideration, the photon undergoes multiple scattering, during which it can also be absorbed and emitted by the medium, but in the form of one or two photons with different parameters. In this case, it is obviously necessary to use a representation that takes into account the influence of the quantum system on the environment and does not allow the loss of information about the “QS + environment" joint system (JS).

Following the articles [29,30,32], where the particle and its random environment are considered as a single closed system within the framework of the random matrix method, we can write the form of the stochastic density matrix (SDM):

$$\rho (\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})=\sum _{n,m}{W}_{nm}{\rho }_{nm}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\}).$$

(91)

In the representation (91) the term $W_{nm}$ denotes the population of the levels of two non-interacting quantum harmonic oscillators with energies; $E_n=(n+1/2){\omega }_1^{-}$ and $E_m=(m+1/2){\omega }_2^{-}$. Integrating expression (91) over the coordinates $q_1,q_2\in {\mathbb{R}}^2$, where ${\mathbb{R}}^2$ is a 2D Euclidean subspace, we find the following normalization condition:

$$\sum _{n,m}{W}_{n,m}=1,W_{n,m}\ge 0.$$

Using the definition of a random density matrix (91), we can construct a regular reduced density matrix by performing functional integration.

For simplicity, we assume that the population of the quantum levels is such that $W_{00}=1$ and $W_{nm}=0$, for all integers $n,m\ge 1$.

Using the first expression in (88) for the RDM, we can construct the von Neumann entropy, which characterizes the measure of disorder of a QS, i.e., two linearly coupled 1D quantum oscillators in the ground state. Recall that the standard definition of von Neumann entropy is as follows [34]:

$${\mathsf{\Lambda }}_N\left(\tau \right)=-T{r}_{\mathbf{q}}\{{\rho }_{00}(\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime })ln{\rho }_{00}(\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime })\}{|}_{{\mathbf{q}}^{\prime }=\mathbf{q};{\tau }^{\prime }=\tau }.$$

(92)

In particular, substituting from the (88) first expression for RDM into the equation (92), we can obtain the following expression:

$${\mathsf{\Lambda }}_N\left(\tau \right)=-{\mathcal{N}}^{\left(1\right)}\left(\tau \right){\mathsf{\Lambda }}_N^{\left(2\right)}\left(\tau \right)-{\mathcal{N}}^{\left(2\right)}\left(\tau \right){\mathsf{\Lambda }}_N^{\left(1\right)}\left(\tau \right),$$

(93)

where ${\mathcal{N}}^{\left(l\right)}\left(\tau \right)$ and ${\mathsf{\Lambda }}_N^{\left(l\right)}\left(\tau \right),(l=1,2)$ denote the number of states and the entropy of the l-th single oscillator, respectively. Note that the number of states is determined by the expression:

$$\begin{array}{c}\hfill {\mathcal{N}}^{\left(l\right)}\left(\tau \right)={\int }_{-\infty }^{+\infty }{\rho }_0^{\left(l\right)}(q_l,\tau )d{q}_j={\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }{\displaystyle \frac{d{u}_2}{\sqrt{u_2}}}{\mathcal{P}}_l(0;u_1,u_2,\tau ),\end{array}$$

(94)

where ${\rho }_0^{\left(l\right)}(q_l,\tau )$ is the RDM of the l-th 1D-th oscillator, which has the form:

$${\rho }_0^{\left(l\right)}(q_l,\tau )={\displaystyle \frac{1}{\sqrt{\pi }}}{\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\mathcal{P}}_l(0;u_1,u_2,\tau )exp\left\{-u_2{q}_l^2\right\}.$$

As for the entropy of the l-th 1D-th oscillator, it is determined by the expression:

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_N^{\left(l\right)}\left(\tau \right)={\int }_{-\infty }^{+\infty }{\rho }_0^{\left(l\right)}(q_l,\tau )ln{\rho }_0^{\left(l\right)}(q_l,\tau )d{q}_l,\end{array}$$

(95)

As can be seen from the expression (95), the von Neumann entropy has a rather complex form, which does not allow for analytical integration of the expression over the coordinate $q_l$. However, it should be noted that the definition (95), unlike the standard definition of von Neumann entropy [34], does not fully account for the influence of a QS, in this case a photon, on its random environment. This is due to the fact that the averaging procedure of a quantum system is carried out only by the distribution ${\mathcal{P}}_l(n;u_1,u_2,\tau )$, which does not preserve the full probability in the JS.

Another, more elegant definition of the entropy of a quantum system, where integration over spatial coordinates $q_l$ can be performed, can be implemented within the framework of the random density matrix method (SDM).

Definition 4. The time-dependent entropy of a quantum subsystem immersed in a thermostat can be determined using the mathematical expectation of a random density matrix:

$${\mathsf{\Lambda }}_G\left(\tau \right)=-T{r}_{\mathbf{q}}\{T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}\left[\rho (\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})ln\rho (\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})]\right\}.$$

(96)

For definiteness, we will consider the entropy of a photon in the ground state. Given the expression for the random wave function of a single 1D oscillator (73), SDM in the ground state can be written in factorized form:

$$\begin{array}{c}\hfill {\rho }_{00}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})=\prod _{l=1}^2{\rho }_0^{\left(l\right)}(q_l,\tau ;\left\{f\right\}:q_l^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\}),\end{array}$$

(97)

where

$${\rho }_0^{\left(l\right)}(q_l,\tau :q_l^{\prime },{\tau }^{\prime })={\displaystyle \frac{1}{\sqrt{\pi }}}exp\left\{-{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }-{\displaystyle \frac{1}{2}}\left[u_2\left(\tau \right)q_l^2+u_2\left({\tau }^{\prime }\right){q_l^{\prime }}^2\right]\right\}.$$

(98)

Substituting (97) into (96) for entropy, we obtain the following expression:

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_G\left(\tau \right)=-T{r}_{\mathbf{q}}\{T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}\left[{\rho }_{00}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})ln{\rho }_{00}(\mathbf{q},\tau ;\left\{f\right\}:{\mathbf{q}}^{\prime },{\tau }^{\prime };\left\{f^{\prime }\right\})\right]\\ \hfill =-T{r}_{\mathbf{q}}\left\{T{r}_{\left\{\mathbf{u}\right(\tau \left)\right\}}{\rho }_0^{\left(1\right)}{\rho }_0^{\left(2\right)}ln\left[{\rho }_0^{\left(1\right)}{\rho }_0^{\left(2\right)}\right]\right\}=-{\mathcal{N}}^{\left(1\right)}\left(\tau \right){\mathsf{\Lambda }}_G^{\left(2\right)}\left(\tau \right)-{\mathcal{N}}^{\left(2\right)}\left(\tau \right){\mathsf{\Lambda }}_G^{\left(1\right)}\left(\tau \right),\end{array}$$

(99)

where ${\mathcal{N}}^{\left(l\right)}=T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left[{\rho }_0^{\left(l\right)}\right]$ and ${\mathsf{\Lambda }}_G^{\left(l\right)}=T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left[{\rho }_0^{\left(l\right)}ln{\rho }_0^{\left(l\right)}\right]$ is the number of states and the entropy of the l-th quantum harmonic oscillator in the ground state, respectively.

In the expression (99) a term of the type $T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left[{\rho }_0^{\left(l\right)}\right]$ can be easily calculated. It describes the number of states of 1D QHO in the ground state, which is immersed in a random environment:

$${\mathcal{N}}^{\left(l\right)}\left(\tau \right)=T{r}_{q_l}\left\{T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}{\rho }_0^{\left(l\right)}\right\}={\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{+\infty }d{u}_2{\displaystyle \frac{1}{\sqrt{u_2}}}{\mathcal{P}}_l(0;u_1,u_2,\tau ).$$

(100)

The term of the type ${\mathsf{\Lambda }}^{\left(l\right)}\left(\tau \right)=T{r}_{q_l}\left\{T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}{\rho }_0^{\left(l\right)}ln{\rho }_0^{\left(l\right)}\right\}$ in expression (99) represents the entropy of one 1D quantum oscillator, which can be written as follows:

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}^{\left(l\right)}\left(\tau \right)=-{\displaystyle \frac{1}{\sqrt{\pi }}}T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left\{e^{-{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }-u_2{q}_l^2}\left[u_2{q}_l^2+{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }\right]\right\}\\ \hfill =-{\displaystyle \frac{1}{2}}{\mathcal{N}}^{\left(l\right)}\left(\tau \right)-{\displaystyle \frac{1}{\sqrt{\pi }}}T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left\{e^{-{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }}\left[{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }\right]\right\}.\end{array}$$

(101)

Now let us move on to the question of calculating the second term in the expression (101), namely the function:

$${\mathcal{I}}^{\left(l\right)}\left(\tau \right)={\displaystyle \frac{1}{\sqrt{\pi }}}T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left\{e^{-{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }}\left[{\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }\right]\right\}.$$

Note that ${\mathcal{I}}^{\left(l\right)}\left(\tau \right)={\partial }_{\lambda }{\mathcal{I}}^{\left(l\right)}(\lambda ;\tau ){|}_{\lambda =0}$, where ${\mathcal{I}}^{\left(l\right)}(\lambda ;\tau )$ has the form:

$$\begin{array}{c}\hfill {\mathcal{I}}^{\left(l\right)}(\lambda ;\tau )={\displaystyle \frac{1}{\sqrt{\pi }}}T{r}_{q_l}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left\{e^{-(1+\lambda ){\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }-u_2{q}_l^2}\right\}.\end{array}$$

(102)

Integrating over the $q_l$ coordinate into expression (102), we find:

$${\mathcal{I}}^{\left(l\right)}(\lambda ;\tau )={\displaystyle \frac{1}{\sqrt{u_2}}}T{r}_{\left\{{\mathbf{u}}^l\left(\tau \right)\right\}}\left\{e^{-(1+\lambda ){\int }_0^{\tau }{u}_1\left({\tau }^{\prime }\right)d{\tau }^{\prime }}\right\}.$$

(103)

Finally, performing functional integration in (103) we obtain:

$${\mathcal{I}}^{\left(l\right)}(\lambda ;\tau )={\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{\infty }d{u}_2{\displaystyle \frac{1}{\sqrt{u_2}}}{\mathcal{P}}_l(\lambda ;u_1,u_2,\tau ),$$

(104)

where the integrand function ${\mathcal{P}}_l(\lambda ;u_1,u_2,\tau )$ satisfies the following PDE:

$$\left\{{\partial }_{\tau }-{\widehat{L}}_l+(1+\lambda )u_1\right\}{\mathcal{P}}_l(\lambda ;u_1,u_2,\tau )=0.$$

(105)

Taking into account the definition of the function ${\mathcal{I}}^{\left(l\right)}\left(\tau \right)$, we can write:

$${\mathcal{I}}^{\left(l\right)}\left(\tau \right)={\int }_{-\infty }^{+\infty }d{u}_1{\int }_0^{\infty }d{u}_2{\displaystyle \frac{1}{\sqrt{u_2}}}{\mathcal{D}}_l(0;u_1,u_2,\tau ),$$

(106)

where ${\mathcal{D}}_l(0;u_1,u_2,\tau )={\partial }_{\lambda }{\mathcal{P}}_l(\lambda ;u_1,u_2,\tau ){|}_{\lambda =0}.$

Now, taking into account (105), it is easy to prove that the function $\mathcal{D}(0;u_1,u_2,\tau )$ is a solution to the following non-homogeneous PDE:

$$\{{\partial }_{\tau }-{\widehat{L}}_l+u_1\}{\mathcal{D}}_l(0;u_1,u_2,\tau )+u_1{\mathcal{P}}_l(0;u_1,u_2,\tau )=0.$$

(107)

As can be easily verified, equation (107) involves a source ${\mathcal{P}}_l(0;u_1,u_2,\tau )$ that obeys the wave equation (89), and so to find a solution ${\mathcal{D}}_l(0;u_1,u_2,\tau )$ we need to solve (107) and (89) together, self-consistently. The entropy of the l-th 1D quantum harmonic oscillator will have the following form:

$${\mathsf{\Lambda }}^{\left(l\right)}\left(\tau \right)=-{\displaystyle \frac{1}{2}}{\mathcal{N}}^{\left(l\right)}-{\mathcal{I}}^{\left(l\right)}\left(\tau \right).$$

(108)

Thus, we have defined two different types of time-dependent entropies for a QS, ${\mathsf{\Lambda }}_N\left(\tau \right)$ and ${\mathsf{\Lambda }}_G\left(\tau \right)$, which do not take into account the influence of the quantum system on the environment. In other words, these entropies describe open systems and are inaccurate in the case of strong interactions between the QS and its environment. This problem can be solved — that is, the question of conserving total probability in a self-organizing closed system (see definition of JS) —by radically changing the definition of entropy, using the mathematical expectation of the photon’s wave function.

Definition 5. The time-dependent entropy of a quantum subsystem self-organizing in a random environment will be denoted by the following expression:

$${\mathsf{\Lambda }}_S\left(\tau \right)=-T{r}_{\mathbf{q}}\left\{\varrho \left(\mathbf{n}\right|\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime })ln\varrho \left(\mathbf{n}\right|\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime })\right\},$$

(109)

where $\varrho \left(\mathbf{n}\right|\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime })$ denotes the mathematical expectation of RDM, which is defined by the following bilinear form (see expressions (74)-(77)), in addition, $\mathbf{n}=n_1,n_2$ are quantum numbers:

$$\varrho \left(\mathbf{n}\right|\mathbf{q},\tau :{\mathbf{q}}^{\prime },{\tau }^{\prime }){|}_{\mathbf{q}={\mathbf{q}}^{\prime };\tau ={\tau }^{\prime }}=\prod _{l=1}^2{\overline{\mathcal{F}}}_l\left(n_l\right|q_l,\tau )·{\overline{\mathcal{F}}}_l^{*}\left(n_l\right|q_l,\tau ).$$

(110)

Below, for the definiteness, we will consider the entropy of a single photon in the ground state during the relaxation process in the random environment.

Substituting the representation for RDM (110) into equation (109), we obtain the following expression for the entropy QS:

$${\mathsf{\Lambda }}_S\left(\tau \right)=-{\mathcal{N}}_S^{\left(1\right)}\left(\tau \right){\mathsf{\Lambda }}_S^{\left(2\right)}\left(\tau \right)-{\mathcal{N}}_S^{\left(2\right)}\left(\tau \right){\mathsf{\Lambda }}_S^{\left(1\right)}\left(\tau \right),$$

(111)

where ${\mathcal{N}}_S^{\left(l\right)}\left(\tau \right)$ and ${\mathsf{\Lambda }}_S^{\left(l\right)}\left(\tau \right)$ denote, respectively, the density of states and the entropy of the l-th 1D QHO, which are defined by following expressions:

$${\mathcal{N}}_S^{\left(l\right)}\left(\tau \right)=\underset{-\infty }{\overset{+\infty }{\int }}d{q}_l{\varrho }_l\left(0\right|q_l,\tau )=\underset{-\infty }{\overset{+\infty }{\int }}{\overline{\mathcal{F}}}_l\left(0\right|q_l,\tau )·{\overline{\mathcal{F}}}_l^{*}\left(0\right|q_l,\tau )d{q}_l,$$

(112)

and accordingly:

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_S^{\left(l\right)}\left(\tau \right)=\underset{-\infty }{\overset{+\infty }{\int }}d{q}_l{\varrho }_l\left(0\right|q_l,\tau )ln{\varrho }_l\left(0\right|q_l,\tau )=\\ \hfill \underset{-\infty }{\overset{+\infty }{\int }}{\overline{\mathcal{F}}}_l\left(0\right|q_l,\tau )·{\overline{\mathcal{F}}}_l^{*}\left(0\right|q_l,\tau )ln\left\{{\overline{\mathcal{F}}}_l\left(0\right|q_l,\tau )·{\overline{\mathcal{F}}}_l^{*}\left(0\right|q_l,\tau )\right\}d{q}_l.\end{array}$$

(113)

Using the expression for the mathematical expectation of the photon wave function (74), we can calculate the coefficient ${\mathcal{N}}_S^{\left(l\right)}\left(\tau \right)$ and represent it as a four-dimensional integral representation:

$$\begin{array}{c}\hfill {\mathcal{N}}_S^{\left(l\right)}\left(\tau \right)=\sqrt{2}\underset{-\infty }{\overset{+\infty }{\int }}d{u}_1\underset{0}{\overset{+\infty }{\int }}d{u}_2\underset{-\infty }{\overset{+\infty }{\int }}d{\overline{u}}_1\underset{0}{\overset{+\infty }{\int }}d{\overline{u}}_2{\displaystyle \frac{Q_l(0;u_1,u_2,\tau )Q_l^{*}(0;{\overline{u}}_1,{\overline{u}}_2,\tau )}{\sqrt{(u_2+{\overline{u}}_2)+i(u_1-{\overline{u}}_1)}}},\end{array}$$

(114)

where

$$Q_l(0;u_1,u_2,\tau )Q_l^{\left(i\right)}{Q}_l^{*}(0;{\overline{u}}_1,{\overline{u}}_2,\tau )=$$

$$Q_l^{\left(r\right)}(0;u_1,u_2,\tau )Q_l^{\left(r\right)}(0;{\overline{u}}_1,{\overline{u}}_2,\tau )+Q_l^{\left(i\right)}(0;u_1,u_2,\tau )Q_l^{\left(i\right)}(0;{\overline{u}}_1,{\overline{u}}_2,\tau )-$$

$$i{Q}_l^{\left(r\right)}(0;u_1,u_2,\tau )Q_l^{\left(i\right)}(0;{\overline{u}}_1,{\overline{u}}_2,\tau )+i{Q}_l^{\left(i\right)}(0;u_1,u_2,\tau )Q_l^{\left(r\right)}(0;{\overline{u}}_1,{\overline{u}}_2,\tau ).$$

Note that the expression (114) defining the density of states is, by definition, a positive real quantity.

Unfortunately, it is not possible to carry out analytical integration of the expression for the entropy of an individual photon (112) over the $q_l$ coordinate, which greatly complicates its numerical calculation.

Thus, we have constructed the entropy of a system of photons using the standard von Neumann definition (95), the SDM method, and in the framework of the mathematical expectation of the photon wave function (109)-(110). Thus, we have constructed the entropy of a system of photons using the standard von Neumann definition (95), the SDM method, and within the ME method of the photon wave function (109)-(110). Under certain conditions, a single-photon system can decay into two photons, and the state of each individual photon will be described by the wave function of a 1D QHO immersed in a random environment. The difference between these representations is that in the framework of the ME method for the wave function, the influence of the quantum system on the environment is correctly taken into account, which is fundamentally important in the case of strong interaction between the quantum system and its environment.

An important feature of expression (95) and (112) is that, as can be seen, if they characterize the evolution of a system of two single photons, then these photons are obviously entangled through the entropies of the individual photons.

7. Transition Probabilities to Bell States

Let us consider a problem where a photon, moving in a nanowaveguide, is scattered by quantum dots randomly located along its trajectory. Clearly, in this case, as a result of inelastic scattering, there is a possibility of the photon being absorbed by the quantum dot, which will eventually emit two entangled photons. For the organization of quantum communications using photons, the degree of their entanglement plays a very important role. In particular, maximum entanglement between photons is achieved using four so-called Bell states (see for example [37,38]).

Mathematically, single-photon states, taking into account their helicity, in a plane perpendicular to the direction of photon propagation, can be represented as two-component vector wave functions of the form:

$$|1\rangle =\left(\begin{array}{c}{|1\rangle }_{\uparrow }\\ {|1\rangle }_{\downarrow }\end{array}\right),|2\rangle =\left(\begin{array}{c}{|2\rangle }_{\uparrow }\\ {|2\rangle }_{\downarrow }\end{array}\right),$$

(115)

where $|1\rangle =|{\mathbf{r}}_1\rangle$ and $|2\rangle =|{\mathbf{r}}_2\rangle$ denote the singlet wave functions of the first and second photons, respectively, ${\mathbf{r}}_1,{\mathbf{r}}_2\in {\mathbb{R}}^2$ denote the radius vectors of localization of the corresponding photons, ${|l\rangle }_{\uparrow }$ and ${|l\rangle }_{\downarrow }(l=1,2)$ denote photons with right and left polarization, respectively. Taking into account the above, we can write down all four time-like Bell states that describe the maximally entangled state of two qubits and form an orthogonal basis in the Hilbert space:

$$\begin{array}{c}\hfill |{\psi }_{\pm }\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left\{{|1\rangle }_{\uparrow }\otimes {|2\rangle }_{\uparrow }\pm {|1\rangle }_{\downarrow }\otimes {|2\rangle }_{\downarrow }\right\},\\ \hfill |{\varphi }_{\pm }\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left\{{|1\rangle }_{\uparrow }\otimes {|2\rangle }_{\downarrow }\pm {|1\rangle }_{\downarrow }\otimes {|2\rangle }_{\uparrow }\right\}.\end{array}$$

(116)

Using the Riemann-Silberstein vector representation [26,27] for two right-handed (left-handed) helicities of photons simultaneously localized in two different spatial positions ${\mathbf{r}}_1={\mathbf{r}}_1(x_1,x_2)$ and ${\mathbf{r}}_2={\mathbf{r}}_2(x_1,x_2)$ at time “$\tau$" (see also expression (10)), the wave functions can be written as follows:

$$\begin{array}{c}\hfill {|1\rangle }_{\uparrow \downarrow }=\left\{{\displaystyle \frac{{\mathbf{D}}_1\left({\mathbf{r}}_1\right)}{\sqrt{\epsilon \left({\mathbf{r}}_1\right)}}}\pm i{\displaystyle \frac{{\mathbf{B}}_1\left({\mathbf{r}}_1\right)}{\sqrt{\mu \left({\mathbf{r}}_1\right)}}}\right\},{|2\rangle }_{\uparrow \downarrow }=\left\{{\displaystyle \frac{{\mathbf{D}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_2\right)}}}\pm i{\displaystyle \frac{{\mathbf{B}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\mu \left({\rho }_2\right)}}}\right\}.\end{array}$$

(117)

Using (117), we can easily construct the wave functions of all two-photon states, namely: $\left\{{|1\rangle }_{\uparrow }\otimes {|2\rangle }_{\uparrow }{;|1\rangle }_{\downarrow }\otimes {|2\rangle }_{\downarrow }{;|1\rangle }_{\uparrow }\otimes {|2\rangle }_{\downarrow }{;|1\rangle }_{\downarrow }\otimes {|2\rangle }_{\uparrow }\right\}$ and, with their help, the corresponding Bell states (116).

Taking into account above, the Bell states (116) can be rewritten in the following form:

$$\begin{array}{c}\hfill |{\psi }_{+}\rangle =\sqrt{2}\left\{{\displaystyle \frac{{\mathbf{D}}_1\left({\rho }_1\right)\otimes {\mathbf{D}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_1\right)\epsilon \left({\mathbf{r}}_2\right)}}}-{\displaystyle \frac{{\mathbf{B}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{B}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\mu \left({\mathbf{r}}_1\right)\mu \left({\mathbf{r}}_2\right)}}}\right\},\\ \hfill |{\psi }_{-}\rangle =i\sqrt{2}\left\{{\displaystyle \frac{{\mathbf{D}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{B}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_1\right)\mu \left({\mathbf{r}}_2\right)}}}+{\displaystyle \frac{{\mathbf{B}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{D}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_2\right)\mu \left({\mathbf{r}}_1\right)}}}\right\},\\ \hfill |{\varphi }_{+}\rangle =\sqrt{2}\left\{{\displaystyle \frac{{\mathbf{D}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{D}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_1\right)\epsilon \left({\mathbf{r}}_2\right)}}}+{\displaystyle \frac{{\mathbf{B}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{B}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\mu \left({\mathbf{r}}_1\right)\mu \left({\mathbf{r}}_2\right)}}}\right\},\\ \hfill |{\psi }_{-}\rangle =i\sqrt{2}\left\{{\displaystyle \frac{{\mathbf{D}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{B}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_1\right)\mu \left({\mathbf{r}}_2\right)}}}-{\displaystyle \frac{{\mathbf{B}}_1\left({\mathbf{r}}_1\right)\otimes {\mathbf{D}}_2\left({\mathbf{r}}_2\right)}{\sqrt{\epsilon \left({\mathbf{r}}_2\right)\mu \left({\mathbf{r}}_1\right)}}}\right\}.\end{array}$$

(118)

The electric and magnetic field strengths of photons can be represented in vector form. In particular, for the first photon, this will look like:

$${\mathbf{D}}_1\left({\mathbf{r}}_1\right)=\left(\begin{array}{c}{D}_1\left(x_1\right)\\ 0\end{array}\right),{\mathbf{B}}_1\left({\mathbf{r}}_1\right)=\left(\begin{array}{c}{B}_1\left(x_1\right)\\ 0\end{array}\right),$$

(119)

and for the second photon, respectively:

$${\mathbf{D}}_2\left({\mathbf{r}}_2\right)=\left(\begin{array}{c}{D}_2\left(x_2\right)\\ 0\end{array}\right),{\mathbf{B}}_2\left({\mathbf{r}}_2\right)=\left(\begin{array}{c}{B}_2\left(x_2\right)\\ 0\end{array}\right).$$

(120)

Now, using the expressions (119)-(120), the equations (118) can be written explicitly:

$$\begin{array}{c}\hfill |{\psi }_{+}\rangle =\sqrt{2}\left\{{\displaystyle \frac{D_1\left(x_1\right)D_2\left(x_2\right)}{\sqrt{\epsilon \left(x_1\right)\epsilon \left(x_2\right)}}}-{\displaystyle \frac{B_1\left(x_1\right)B_2\left(x_2\right)}{\sqrt{\mu \left(x_1\right)\mu \left(x_2\right)}}}\right\},\\ \hfill |{\psi }_{-}\rangle =i\sqrt{2}\left\{{\displaystyle \frac{D_1\left(x_1\right)B_2\left(x_2\right)}{\sqrt{\epsilon \left(x_1\right)\mu \left(x_2\right)}}}+{\displaystyle \frac{B_1\left(x_1\right)D_2\left(x_2\right)}{\sqrt{\mu \left(x_1\right)\epsilon \left(x_2\right)}}}\right\},\\ \hfill |{\varphi }_{+}\rangle =\sqrt{2}\left\{{\displaystyle \frac{D_1\left(x_1\right)D_2\left(x_2\right)}{\sqrt{\epsilon \left(x_1\right)ϵ\left(x_2\right)}}}+{\displaystyle \frac{B_1\left(x_1\right)B_2\left(x_2\right)}{\sqrt{\mu \left(x_1\right)\mu \left(x_2\right)}}}\right\},\\ \hfill |{\varphi }_{-}\rangle =i\sqrt{2}\left\{{\displaystyle \frac{D_1\left(x_1\right)B_2\left(x_2\right)}{\sqrt{\epsilon \left(x_1\right)\mu \left(x_2\right)}}}-{\displaystyle \frac{B_1\left(x_1\right)D_2\left(x_2\right)}{\sqrt{\mu \left(x_1\right)\epsilon \left(x_2\right)}}}\right\},\end{array}$$

(121)

where $D_{1,2}\left(x_1\right)$ and $B_{1,2}\left(x_2\right)$ are the scalar components of the electric and magnetic fields of the corresponding photons.

Now, we can move on to the question of constructing the mathematical expectation of Bell states. For definiteness, let us consider transitions from the ground state (70) and (74) to the maximally entangled Bell states (116). Assuming that the Bell states (118) describe $\left(out\right)$ asymptotic states, we can define the transition $\mathrm{S}$-matrix elements by projecting the evolving full wave function of the single-photon state (70) onto the various Bell states (121):

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{\pm }\rangle }=\underset{\tau \to +\infty }{lim}\underset{-\infty }{\overset{+\infty }{\int }}d{x}_1\underset{-\infty }{\overset{+\infty }{\int }}d{x}_2\left\{\prod _{l=1}^2\mathbb{E}\left[{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})\right]\right\}·{|{\psi }_{\pm }\rangle }^{*},\\ \hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\varphi }_{\pm }\rangle }=\underset{\tau \to +\infty }{lim}\underset{-\infty }{\overset{+\infty }{\int }}d{x}_1\underset{-\infty }{\overset{+\infty }{\int }}d{x}_2\left\{\prod _{l=1}^2\mathbb{E}\left[{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})\right]\right\}·{|{\varphi }_{\pm }\rangle }^{*},\end{array}$$

(122)

where “*” denotes complex conjugation.

For further calculations, it is convenient to represent the dependences of the electric and magnetic fields in the form of Gaussian distributions:

$$\begin{array}{c}\hfill {\mathbb{G}}_l\left(x_l\right)={\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{\left(h\right)l}}}}exp\left\{-{\displaystyle \frac{{(x_l-{\mu }_{\left(h\right)l})}^2}{2{\sigma }_{\left(h\right)l}^2}}\right\},\mathbb{G}=B,D,h=b,d,\end{array}$$

(123)

where ${\mu }_{\left(h\right)l}$ denotes the center of the corresponding distribution, whereas ${\sigma }_{\left(h\right)l}$ is the standard deviation of the distribution.

Now we can write down explicit forms of all elements of the transition $\mathrm{S}$-matrix, representing them in the form of four-fold integral representations:

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{-}\rangle ;|{\varphi }_{-}\rangle \}}={\displaystyle \frac{1}{\sqrt{\pi }}}\underset{\tau \to +\infty }{lim}\prod _{l=1}^2{\int }_{-\infty }^{+\infty }d{u}_1^{\left(l\right)}{\int }_0^{+\infty }d{u}_2^{\left(l\right)}{Q}_l(0;u_1^{\left(l\right)},u_2^{\left(l\right)},\tau )\times \\ \hfill \{{\displaystyle \frac{1}{\epsilon }}{\left\{{\sigma }_{\left(d\right)1}{\sigma }_{\left(d\right)2}\left[{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(d\right)l}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right]\right\}}^{-1/2}\\ \hfill -{\displaystyle \frac{1}{\mu }}{\left\{{\sigma }_{\left(b\right)1}{\sigma }_{\left(b\right)2}\left[{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(b\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right]\right\}}^{-1/2}\},\\ \hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{+}\rangle ;|{\varphi }_{+}\rangle \}}={\displaystyle \frac{-i}{\sqrt{\pi }}}\underset{\tau \to +\infty }{lim}\prod _{l=1}^2{\int }_{-\infty }^{+\infty }d{u}_1^{\left(l\right)}{\int }_0^{+\infty }d{u}_2^{\left(l\right)}{Q}_l(0;u_1^{\left(l\right)},u_2^{\left(l\right)},\tau )\times \\ \hfill {\displaystyle \frac{1}{\sqrt{\epsilon \mu }}}\{{\left[{\sigma }_{\left(d\right)1}{\sigma }_{\left(b\right)2}\left\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+(2{\sigma }_{\left(d\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)(2{\sigma }_{\left(b\right)2}^{-2}+{\mathcal{A}}_1{\mathcal{A}}_2)\right\}\right]}^{-1/2}+\\ \hfill +\left[{\sigma }_{\left(b\right)1}{\sigma }_{\left(d\right)2}\right\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+(2{\sigma }_{\left(b\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\end{array}$$

(124)

(2-2(d)2+A1+A2)}]-1/2}. where coefficients ${\mathcal{A}}_1$ and ${\mathcal{A}}_2$ are defined by the formula, see (74).

Finally, it should be noted that the probability of transition to a specific Bell state will be determined by the square of the modulus of the transition $\mathrm{S}$-matrix element, i.e., ${\mathfrak{W}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{+}\rangle ;|{\varphi }_{+}\rangle \}}={\left|{\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{+}\rangle ;|{\varphi }_{+}\rangle \}}\right|}^2$ and ${\mathfrak{W}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{-}\rangle ;|{\varphi }_{-}\rangle \}}={\left|{\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to \left\{\right|{\psi }_{-}\rangle ;|{\varphi }_{-}\rangle \}}\right|}^2$. It is important to note that the impossibility of factorizing a fourfold integration over the extended 4D subspace into two twofold integrations reflects the fact that the two photon states are entangled.

8. Initial-Boundary Conditions for a System of Two Second-Order PDEs

To simplify the numerical modeling of the system of PDEs (77), we will define and explicitly formulate the initial-boundary value problem (81)-(83) in 2D Euclidean space:

$${\mathsf{\Xi }}_{\left\{\mathbf{u}\right\}}^{\left(2\right)}\cong {\mathbb{R}}^2\equiv (-\infty ,+\infty )\times [0,+\infty ),\mathbf{u}=\{u_1,u_2\}.$$

Recall that although the initial conditions of the problem (81) (see also the extended condition of its normalization (79)) is clearly defined, nevertheless, in this form it is not suitable for numerical calculation. Based on physical considerations, it is more convenient to approximate the 2D Dirac delta-function with a 2D Gaussian function:

$$Q_l^{\left(\upsilon \right)}(u_1,u_2,\tau ){|}_{\tau =0}\approx {\mathfrak{B}}_{l}exp\left\{-{\mathfrak{D}}_{\upsilon }\left[{\left({\displaystyle \frac{u_1}{a_{\upsilon }}}\right)}^2+{\left({\displaystyle \frac{u_2}{b_{\upsilon }}}\right)}^2\right]\right\},\upsilon =i,r.$$

(125)

In particular, for the sake of definiteness we will use the following set of values of the parameters included in the Gaussian function; ${\mathfrak{B}}_l=500,{\mathfrak{D}}_{\upsilon }=\pi {\mathfrak{B}}_l{a}_{\upsilon }^2$ and $a_{\upsilon }=b_{\upsilon }=1/2$. It should be noted that this set of parameter values ensures that each individual 1D Gaussian function is normalized to unity (see Figure 2). Based on the conditions of the problem under consideration, other combinations of parameters can of course be used.

8.1. Boundary Conditions on One Axis

Turning to the question of boundary conditions, let us first consider the behavior of the system of equations (77) near the axis $u_1\in (-\infty ,+\infty )$. In this case, the equations are separated and take the same form, which allows solutions to be sought in the form:

$$Q_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau ){|}_{u_2\sim 0}\approx \left(1+{\displaystyle \frac{1}{2}}{\alpha }_0{u}_2^2\right)exp\left\{-{\displaystyle \frac{1}{2}}{\alpha }_0{u}_2^2\right\}{\mathbb{Q}}_l^{\left(\upsilon \right)}(n;u_1,\tau ),\upsilon =i,r,$$

(126)

where ${\alpha }_0>0$ is a constant that will be determined based on physical considerations, and n is the vibrational quantum number.

From expression (126) it obviously follows that in the limit of probabilistic equilibrium for the solution $Q_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau )$ the limit transition is valid: ${lim}_{u_2\to 0}{Q}_l^{\left(\upsilon \right)}(u_1,u_2,\tau )={\mathbb{Q}}_l^{\left(\upsilon \right)}(n;u_1,\tau )$ (see also expression (83), which defines the Neumann boundary condition).

Now, substituting (126) into (77), that is, into the PDE of two variables, we obtain the following second-order PDE, but with one variable:

$${\displaystyle \frac{\partial }{\partial \tau }}{\mathbb{Q}}_l^{\left(\upsilon \right)}={\displaystyle \frac{\partial }{\partial u_1}}\{{ϵ}^{\left(r\right)}{\displaystyle \frac{\partial }{\partial u_1}}{\mathbb{Q}}_l^{\left(\upsilon \right)}+\left(u_1^2+{\omega }_l^2\left(\tau \right)\right){\mathbb{Q}}_l^{\left(\upsilon \right)}\}-(n-3/2)u_1{\mathbb{Q}}_l^{\left(\upsilon \right)}.$$

(127)

Recall that the solution ${\mathbb{Q}}_l^{\left(\upsilon \right)}(u_1,\tau )$ is a boundary condition for the solution of the system of PDEs (83), which we define on the $u_1$ axis. In addition, it is easy to see that the boundary conditions reflect the fact that the environment is quantized.

Note that in a similar way one can define the Neumann initial-boundary value problem (90) for solving a single PDE (see (89)). In particular, a Gaussian distribution type (125) can be chosen as the initial condition. As for the boundary condition for solving the PDE (89) for the probability distribution ${\mathcal{P}}_l(n;u_1,u_2,\tau )$, then, substituting the representation:

$${\mathcal{P}}_l(n;u_1,u_2,\tau ){|}_{u_2\sim 0}\approx \left(1+{\displaystyle \frac{1}{2}}{\alpha }_0{u}_2^2\right)exp\left\{-{\displaystyle \frac{1}{2}}{\alpha }_0{u}_2^2\right\}{\mathbb{P}}_l(n;u_1,\tau ),$$

(128)

into equation (89), in the limit $u_2\to 0$ one can find the following PDE of one variable (see also the second condition in (90)):

$${\displaystyle \frac{\partial }{\partial \tau }}{\mathbb{P}}_l={\displaystyle \frac{\partial }{\partial u_1}}\{{ϵ}^{\left(r\right)}{\displaystyle \frac{\partial }{\partial u_1}}{\mathbb{P}}_l+\left(u_1^2+{\omega }_l^2\left(\tau \right)\right){\mathbb{P}}_l\}-(2n-1)u_1{\mathbb{P}}_l.$$

(129)

After normalizing the solution ${\mathbb{P}}_l^{\left(\upsilon \right)(n;u_1,\tau )}$, it acquires the meaning of a probability density.

It should be noted that this equation is a special case of the well-known continuity equation (see [40,41]), which has an exact solution only for $n=3/2$ and at large times. However, in the case under consideration, n is an vibrational quantum number and, accordingly, can only take integer values:

$${\mathbb{Q}}_l^{\left(\upsilon \right)}(n=3/2;u_1,\infty )={\displaystyle \frac{\mathcal{Z}}{\mathcal{C}\left({\mathcal{E}}_l\right)}}exp\left\{-{\displaystyle \frac{u_1^3}{3\mathcal{Z}}}-{\displaystyle \frac{{\mathcal{E}}_l{u}_1}{\mathcal{Z}}}\right\}{\int }_{-\infty }^{u_1}exp\left\{{\displaystyle \frac{t^3}{3\mathcal{Z}}}+{\displaystyle \frac{{\mathcal{E}}_{l}t}{\mathcal{Z}}}\right\}dt,$$

(130)

where ${\mathbb{Q}}_l^{\left(\upsilon \right)}(u_1,\infty )={lim}_{\tau \to +\infty }{\mathbb{Q}}_l^{\left(\upsilon \right)}(u_1,\tau )$, in addition, ${\mathcal{E}}_l={\left({\omega }_l^{+}\right)}^2$ and $\mathcal{Z}=1/{ϵ}^{\left(r\right)}$.

As for the normalization constant $\mathcal{C}\left({\mathcal{E}}_l\right)$, it is determined from the condition of normalization of the distribution ${\mathbb{Q}}_l^{\left(\upsilon \right)}(u_1,\infty )$ by one and has the form:

$${\mathcal{C}}^{-1}\left({\mathcal{E}}_l\right)={\displaystyle \frac{\sqrt{\pi }}{{\mathcal{Z}}^{1/3}}}{\int }_0^{\infty }{z}^{-1/2}exp\left\{-{\displaystyle \frac{z^3}{12}}-{\displaystyle \frac{{\mathcal{E}}_{l}z}{{\mathcal{Z}}^{2/3}}}\right\}dz.$$

Thus, we have defined natural boundary conditions on the $u_1$ axis satisfying equation (127) (see also condition (83)), which together with the initial conditions (81) defines the Neumann initial-boundary value problem for the system of PDEs (77). As shown above, for the solution of PDE (89) we obtain an initial-boundary conditions similar to (127), only now for the distribution function ${\mathbb{P}}_l(n;u_1,\tau )$. It is worth noting that the initial and boundary conditions of the problem can be modified depending on the specifics of the problem under consideration. This will allow us to more accurately describe various quantum processes in the environment with involving a photon, including phenomena such as the absorption and emission of photons in the environment.

8.2. Boundary Conditions on Two Perpendicular Axes

Now let us consider the behavior of the system of PDEs (77) near two perpendicular axes $u_1\in (-\infty ,+\infty )$ and $u_2\in [0,+\infty )$. As we have seen above, near the $u_1$ axis, both equations of the system (77) are separated and take the same form, in connection with which the same boundary conditions are considered for them. Below we will use the solution ${\mathbb{Q}}_l^{\left(r\right)}(n;u_1,\tau )$ (see expression (127)) as a boundary condition on the $u_1$ axis, but we still need to find a second boundary condition on the $u_2$ axis, which we denote by the function ${\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )$.

Let us consider the solutions $Q_l^{\left(r\right)}$ and $Q_l^{\left(i\right)}$, representing them near the $u_2$ axis in the form:

$$Q_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau ){|}_{u_1\sim 0}\approx \left(1+{\displaystyle \frac{1}{2}}{\beta }_0{u}_1^2\right)exp\left\{-{\displaystyle \frac{1}{2}}{\beta }_0{u}_1^2\right\}{\mathbb{Q}}_l^{\left(\upsilon \right)}(n;u_1,u_2,\tau ),\upsilon =i,r,$$

(131)

where ${\beta }_0>0$ is a some constant that will be determined based on physical considerations. Substituting the representation (131) into the system of PDEs (77), we obtain:

$$\{\begin{array}{c}{\partial }_{\tau }{\mathbb{Q}}_l^{\left(r\right)}={ϵ}^{\left(i\right)}{\partial }_{u_2{u}_2}^2{\mathbb{Q}}_l^{\left(r\right)}+(n+1/2)u_2{\mathbb{Q}}_l^{\left(i\right)},\hfill \\ {\partial }_{\tau }{\mathbb{Q}}_l^{\left(i\right)}={ϵ}^{\left(i\right)}{\partial }_{u_2{u}_2}^2{\mathbb{Q}}_l^{\left(i\right)}-(n+1/2)u_2{\mathbb{Q}}_l^{\left(r\right)},\hfill \end{array}$$

(132)

where ${\partial }_{u_2{u}_2}^2={(\partial /\partial u_2)}^2.$

As can be seen from the system (132), the equations on the $u_2$ axis are not separated and remain connected, despite the obvious property of symmetry. In particular, when replacing the coordinates $u_2\to -u_2$ in the first equation, we obtain:

$${\partial }_{\tau }{\mathbb{Q}}_l^{\left(r\right)}(n;-u_2,\tau )={ϵ}^{\left(i\right)}{\displaystyle \frac{{\partial }^2}{\partial u_2^2}}{\mathbb{Q}}_l^{\left(r\right)}(n;-u_2,\tau )-(n+1/2)u_2{\mathbb{Q}}_l^{\left(i\right)}(n;-u_2,\tau ).$$

Now, if in the equation we make the following substitutions:

$${\mathbb{Q}}_l^{\left(r\right)}(n;-u_2,\tau )={\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau ),{\mathbb{Q}}_l^{\left(i\right)}(n;-u_2,\tau )={\mathbb{Q}}_l^{\left(r\right)}(n;u_2,\tau ),$$

then, obviously, from the first equation of the system (132) we obtain the second equation. In other words, the system of equations (132) can be split into two uncoupled equations of the form:

$$\{\begin{array}{c}{\partial }_{\tau }{\mathbb{Q}}_l^{\left(r\right)}(n;u_2,\tau )={ϵ}^{\left(i\right)}{\partial }_{u_2{u}_2}^2{\mathbb{Q}}_l^{\left(r\right)}(n;u_2,\tau )+(n+1/2)u_2{\mathbb{Q}}_l^{\left(r\right)}(n;-u_2,\tau ),\hfill \\ {\partial }_{\tau }{\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )={ϵ}^{\left(i\right)}{\partial }_{u_2{u}_2}^2{\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )-(n+1/2)u_2{\mathbb{Q}}_l^{\left(i\right)}(n;-u_2,\tau ).\hfill \end{array}$$

(133)

It is obvious that these equations cannot be solved separately due to the shift of the argument of the desired function, and therefore such a separation can be considered conditional. Nevertheless, this representation of the equation allows us to perform additional reasoning to find new, nontrivial solutions.

Now for definiteness, let us assume that the function ${\mathbb{Q}}_l^{\left(i\right)}(n;-u_2,\tau )$ is even. Then the second equation in the (133) can be written as follows:

$${\partial }_{\tau }{\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )={ϵ}^{\left(i\right)}{\partial }_{u_2{u}_2}^2{\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )-(n+1/2)u_2{\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau ).$$

(134)

Representing the solution ${\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )$ in the form:

$${\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )=e^{-E\tau }{\tilde{\mathbb{Q}}}_l^{\left(i\right)}(n;u_2),$$

(135)

from (134) it is easy to find the following stationary ordinary differential equation (ODE):

$$\left\{{ϵ}^{\left(i\right)}{d}_{u_2{u}_2}^2+\left[E-(n+1/2)u_2\right]\right\}{\tilde{\mathbb{Q}}}_l^{\left(i\right)}(n;u_2)=0,$$

(136)

where the coordinate $u_2>0$ and E is a constant like energy.

By requiring the boundary conditions:

$${\tilde{Q}}_l^{\left(i\right)}(n;u_2){{|}_{u_2=0}=0,{\tilde{Q}}_l^{\left(i\right)}(n;u_2)|}_{u_2\to \infty }\to 0$$

(137)

to be satisfied for solving ODE (136), we can obtain the following analytical expression (see for example [44]):

$${\tilde{\mathbb{Q}}}_l^{\left(i\right)}(n;u_2)=\mathrm{Ai}\left(\zeta \right)={\displaystyle \frac{1}{\pi }}\sqrt{{\displaystyle \frac{\zeta }{3}}}{K}_{1/3}\left({\displaystyle \frac{2}{3}}\zeta \right),\zeta =\left(u_2-{\displaystyle \frac{2E}{2n+1}}\right){\left({\displaystyle \frac{2n+1}{ϵ}}\right)}^{1/3},$$

(138)

where $\zeta >0$, in addition, $\mathrm{Ai}\left(x\right)$ and $K_{1/3}\left(x\right)$ are the Airy function and the modified Hankel function, respectively.

For negative values of the variable $\zeta <0$, the Airy function is represented as a sum (see [45]):

$$\mathrm{Ai}(-\zeta )={\displaystyle \frac{\sqrt{\zeta }}{3}}\left[J_{1/3}\left({\displaystyle \frac{2}{3}}{\zeta }^{3/2}\right)+J_{-1/3}\left({\displaystyle \frac{2}{3}}{\zeta }^{3/2}\right)\right],\zeta <0,$$

where $J_{1/3}\left(x\right)$ denotes the Bessel function of the first kind of order $1/3.$

Due to the second boundary condition (137), the solution $\mathrm{Ai}(-\zeta )$ must vanish at values $\zeta =-2E{\left[ϵ{(2n+1)}^2\right]}^{-1/3}$ that have an energetic meaning and represent a discrete set.

Finally, we can set the boundary conditions for the solution of problem (77) on two perpendicular axes, namely, ${\mathbb{Q}}_l^{\left(r\right)}(n;u_1,\tau )$ satisfying PDE (127) on the $u_1$ axis and ${\mathbb{Q}}_l^{\left(i\right)}(n;u_2,\tau )$ being the solution of PDE system (132) on the $u_2$ axis. It is easy to see that the new boundary conditions, like the previous ones in subsection 8.1, are quantized. As for the initial conditions, they remain unchanged as defined in subsection 8.1, see expression (125).

9. Algorithms and Numerical Methods. Visualization of Numerical Calculations

A complete quantum description of photon motion in a nanowaveguide with randomly located two-level quantum dots inside, at which the photon can be absorbed and then emitted, is mathematically correctly described using integral representations (74), (88), and (Section 7). Recall that in these representations the integrand is a solution of a second-order PDE (89) or a system of two such PDEs (77), which are quantized, i.e. depend on the vibrational quantum numbers $n_1$ and $n_2$. Further study of the problem reduces to a rather complex computational task, where the above-mentioned parabolic PDEs play a key role.

9.1. Algorithms for Solving PDE and System Consisting of Such Equations

Based on preliminary analysis and test calculations, we found it appropriate to choose an explicit finite-difference scheme for the numerical solution of the problem that is the convection–diffusion type, see for example [43]. Note that its accuracy in coordinates is second-order and its accuracy in time is first-order.

Recall that in ([32]) we have already considered equations of the same type. In particular, various numerical solution methods using finite-difference schemes were tested. In this paper, we will repeat some of the indicated schemes for the solved PDEs (77) and (89). Note that, despite the apparent simplicity of this method, based on our experience in previous studies, we consider its use to be fully justified in terms of efficiency, accuracy and stability.

Listing 1. Numerical algorithm for solving PDEs (89).

1. The continuous domain ${\mathbb{R}}^2$ for equation (89) is replaced by the computational domain $[{\left(u_1\right)}_{min},{\left(u_1\right)}_{max}]\times [{\left(u_2\right)}_{min},{\left(u_2\right)}_{max}]\times [0,T]$. In the computational domain, a uniform difference grid is specified over time t and spatial coordinates $u_1$ and $u_2$:

$$\begin{array}{c}\hfill {\left(u_1\right)}_j=j\mathsf{\Delta }{u}_1,j\in [1,M],{\left(u_1\right)}_{min}={\left(u_1\right)}_{j=1},{\left(u_1\right)}_{max}={\left(u_1\right)}_{j=M},\\ \hfill {\left(u_2\right)}_k=k\mathsf{\Delta }{u}_2,k\in [1,L],{\left(u_2\right)}_{min}={\left(u_2\right)}_{k=1},{\left(u_2\right)}_{max}={\left(u_2\right)}_{k=L},\\ \hfill {\tau }^m=\mathsf{\Delta }\tau ,m=0,1,2,...T/\mathsf{\Delta }\tau -1,\end{array}$$

(139)

where $\mathsf{\Delta }{u}_1$ and $\mathsf{\Delta }{u}_2$ denote steps on the spatial coordinates, while $\mathsf{\Delta }\tau$ is the step of time. The function f on the grid is written as follows: $f_{j,k}^m=f[{\left(u_1\right)}_j,{\left(u_2\right)}_k,{\tau }^m]$, where $f={\mathcal{P}}_l^{\left(n\right)}={\mathcal{P}}_l(n;u_1,u_2,\tau )$.

In these notations, equation (89) can be written in the following difference form:

$$\begin{array}{c}\hfill {\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k}^{m+1}={\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k}^m+g_1\left\{{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j+1,k}^m-2{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k}^m+{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j-1,k}^{m+1}\right\}+g_2\{{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k+1}^m\\ \hfill -2{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k}^m+{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k-1}^m\}+{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_1}}\left[{\left(u_1^2\right)}_j-{\left(u_2^2\right)}_k+{\omega }_l^m\right]\left\{{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j+1,k}^m-{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j-1,k}^m\right\}\\ \hfill +{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_2}}{\left(u_1\right)}_j{\left(u_2\right)}_k\left\{{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k+1}^m-{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k-1}^m\right\}+\mathsf{\Delta }\tau [4-(2n+1)]{\left(u_1\right)}_j{\left[{\mathcal{P}}_l^{\left(n\right)}\right]}_{j,k}^m,\end{array}$$

(140)

where the following notations are made:

$${\omega }_l^m={\omega }_l\left({\tau }^m\right),g_1={ϵ}_l^{\left(r\right)}{\displaystyle \frac{\mathsf{\Delta }\tau }{{\left(\mathsf{\Delta }{u}_1\right)}^2}},g_2={ϵ}_l^{\left(i\right)}{\displaystyle \frac{\mathsf{\Delta }\tau }{{\left(\mathsf{\Delta }{u}_2\right)}^2}}.$$

2. Instead of the Dirac delta-function, the Gaussian distribution is used as the initial condition in this work, see Section 8, expression (125).

3. To calculate PDE (89), it is necessary to set a boundary condition in the form of difference equation on the coordinate axis $u_1$, which can be easily found by approximating Equation (89), see Section 8, expression (129). Note that these difference equations must be solved taking into account the Dirichlet boundary condition; ${\mathbb{P}}_l{|}_{u_1=\pm \left|{\lambda }_0\right|}=0,\left(\right|{\lambda }_0|>>1,l=1,2).$

Listing 2. Numerical algorithm for solving the system of PDEs (77).

1. The continuous region ${\mathbb{R}}^2$ for the system of PDEs (77) is replaced by a discrete grid, as described in Listing 1. Similarly, using (77), we can write the difference equations for this system of PDEs:

$$\begin{array}{c}\hfill {\left[Q_l^{\left(r\right)}\right]}_{j,k}^{m+1}={\left[Q_l^{\left(r\right)}\right]}_{j,k}^m+g_1\left\{{\left[Q_l^{\left(r\right)}\right]}_{j+1,k}^m-2{\left[Q_l^{\left(r\right)}\right]}_{j,k}^m+{\left[Q_l^{\left(r\right)}\right]}_{j-1,k}^{m+1}\right\}+g_2\{{\left[Q_l^{\left(r\right)}\right]}_{j,k+1}^m\\ \hfill -2{\left[Q_l^{\left(r\right)}\right]}_{j,k}^m+{\left[Q_l^{\left(r\right)}\right]}_{j,k-1}^m\}+{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_1}}\left[{\left(u_1^2\right)}_j-{\left(u_2^2\right)}_k+{\omega }_l^m\right]\left\{{\left[Q_l^{\left(r\right)}\right]}_{j+1,k}^m-{\left[Q_l^{\left(r\right)}\right]}_{j-1,k}^m\right\}\\ \hfill +{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_2}}{\left(u_1\right)}_j{\left(u_2\right)}_k\left\{{\left[Q_l^{\left(r\right)}\right]}_{j,k+1}^m-{\left[Q_l^{\left(r\right)}\right]}_{j,k-1}^m\right\}-\mathsf{\Delta }\tau (n-7/2){\left(u_1\right)}_j{\left[Q_l^{\left(r\right)}\right]}_{j,k}^m+\\ \hfill \mathsf{\Delta }\tau (n+1/2){\left(u_2\right)}_k{\left[Q_l^{\left(i\right)}\right]}_{j,k}^m,\end{array}$$

(141)

and, accordingly:

$$\begin{array}{c}\hfill {\left[Q_l^{\left(i\right)}\right]}_{j,k}^{m+1}={\left[Q_l^{\left(i\right)}\right]}_{j,k}^m+g_1\left\{{\left[Q_l^{\left(i\right)}\right]}_{j+1,k}^m-2{\left[Q_l^{\left(i\right)}\right]}_{j,k}^m+{\left[Q_l^{\left(i\right)}\right]}_{j-1,k}^{m+1}\right\}+g_2\{{\left[Q_l^{\left(i\right)}\right]}_{j,k+1}^m\\ \hfill -2{\left[Q_l^{\left(i\right)}\right]}_{j,k}^m+{\left[Q_l^{\left(i\right)}\right]}_{j,k-1}^m\}+{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_1}}\left[{\left(u_1^2\right)}_j-{\left(u_2^2\right)}_k+{\omega }_l^m\right]\left\{{\left[Q_l^{\left(i\right)}\right]}_{j+1,k}^m-{\left[Q_l^{\left(i\right)}\right]}_{j-1,k}^m\right\}\\ \hfill +{\displaystyle \frac{\mathsf{\Delta }\tau }{2\mathsf{\Delta }{u}_2}}{\left(u_1\right)}_j{\left(u_2\right)}_k\left\{{\left[Q_l^{\left(i\right)}\right]}_{j,k+1}^m-{\left[Q_l^{\left(i\right)}\right]}_{j,k-1}^m\right\}-\mathsf{\Delta }\tau (n-7/2){\left(u_1\right)}_j{\left[Q_l^{\left(i\right)}\right]}_{j,k}^m-\\ \hfill \mathsf{\Delta }\tau (n+1/2){\left(u_2\right)}_k{\left[Q_l^{\left(r\right)}\right]}_{j,k}^m.\end{array}$$

(142)

Note that similarly for the functions ${\mathbb{P}}_l^{\left(n\right)},{\mathbb{Q}}_l^{\left(r\right)}$ and ${\mathbb{Q}}_l^{\left(i\right)}$ describing the boundary conditions, the corresponding PDEs on the coordinate axis $u_1$ (see Section 8 are derived.

3. An important property of the probability densities $Q_l^{\left(r\right)}(u_1,u_2,\tau )$ and $Q_l^{\left(i\right)}(u_1,u_2,\tau )$, which is very important in numerical calculations, is that at the center of the coordinate axes they must satisfy the following conditions, namely:

$$Q_l^{\left(r\right)}(u_1,u_2=0,\tau ){|}_{u_1=0}=Q_l^{\left(r\right)}(u_1=0,u_2,\tau ){|}_{u_2=0},$$

$$Q_l^{\left(i\right)}(u_1,u_2=0,\tau ){|}_{u_1=0}=Q_l^{\left(i\right)}(u_1=0,u_2,\tau ){|}_{u_2=0}.$$

4. In addition, we will proceed from the assumption that at the boundary of the computational domain, solutions are subject to the following conditions:

$$Q_l^{\left(r\right)}(u_1,u_2,\tau ){|}_{\partial \mathsf{\Omega }}=0,Q_l^{\left(i\right)}(u_1,u_2,\tau ){|}_{\partial \mathsf{\Omega }}=0,$$

where $\partial \mathsf{\Omega }$ denotes the boundaries of the computational domain.

5. As the initial condition for solving the system of PDEs (77), a Gaussian distribution with parameters as in Listing 1 is chosen for both functions $Q_l^{\left(r\right)}(u_1,u_2,\tau )$ and $Q_l^{\left(i\right)}(u_1,u_2,\tau )$.

Finally, we note that the calculation was carried out for the integration region ${\mathbb{R}}^2\cong (-\infty ,+\infty )\times [0,+\infty )$ on a grid of $1100\times 900$ nodes for $u_1\times u_2$, respectively, steps in space: $\mathsf{\Delta }{u}_1=\mathsf{\Delta }{u}_2=0.025$, time step $\mathsf{\Delta }\tau ={10}^{-5}$.

9.2. Numerical Modeling of Photon’s Parameters and Their Visualization

As we have seen above, for constructing mathematical expectations of different parameters of a photon, the key role is played by the distributions of the fields of the photon’s environment in the limit of statistical equilibrium described by the system of PDEs (77).

To numerically solve the system of PDEs (77), the following values of frequencies and diffusion powers parameters were used. Specifically, to determine the frequencies (see equation (39)) we selected the parameters listed in Table 1:

As for the powers of fluctuations, for definitency we carried out calculations for specific values ${ϵ}^{\left(r\right)}={ϵ}^{\left(i\right)}=1$.

9.3. Distribution of Electromagnetic Fields Energy in a Plane Perpendicular to the Motion of a Photon

As noted above, the structure of a photon is a very important characteristic of a particle, which ultimately determines the cross-sections of the processes in which it participates. Based on the expressions obtained for a photon (85) that has undergone multiple elastic and inelastic collisions in a medium with randomly located two-level quantum dots, we calculated the energy distribution of electromagnetic fields in the photon using the data in Table 1. As can be seen from the Figure 4, photons of the same frequency can have completely different structures depending on the vibrational quantum numbers $n_1$ and $n_2$. Calculations also show that their shape, depending on the evolution time, practically stops changing from the moment $\tau =3$ and all distributions are saturated and become stable.

9.4. Reduced Density Matrix Elements

As is known, the reduced density matrix plays an important role in studying the statistical properties of a quantum system, therefore, to deeply understand the evolution of the system, we carried out test calculations of the first four elements of the density matrix. It is important to note that the developed mathematical approach, conventionally called quantum mechanics with its environment, continuously describes the state of a quantum system in various conditions, including its transition from a free to a bound state. In other words, even if a photon is absorbed by the environment, the approach allows us to track its evolution. To calculate the first four matrix elements, we chose the initial boundary conditions so that at the beginning of the evolution they described a photon that was absorbed by the environment with a high probability (about 90 percent). In other words, we consider the entangled quantum state of a photon with the environment, when 90 percent of it is absorbed by the environment and 10 percent is not absorbed, and track its evolution (see Figure 5).

The first row shows the time-dependent element of the reduced density matrix in the ground state, i.e., in ${\rho }_{00}(x_1,x_2,\tau )$. As can be seen from the graphs, the distribution is quickly established and changes very slowly and little with time. The second row shows the time-dependent element of the reduced density matrix ${\rho }_{01}(x_1,x_2,\tau )$, describing the excited photon absorbed by the environment at the first instant of time. As can be seen from the graphs, the distribution changes rapidly and, over time, is reconstructed into the state of a free photon, which can be interpreted as the emission of a photon by the environment. It is interesting to note that when a photon was absorbed by the environment and is in the ground state, its wave state does not change, i.e., over time, the environment does not emit it again, as it happens in the case of an excited photon in a captured state. In the third and fourth lines, where the plots of the matrix elements ${\rho }_{10}(x_1,x_2,\tau )$ and ${\rho }_{11}(x_1,x_2,\tau )$ are given, we again observe the rapid evolution of the excited photon state captured by the environment and its transition to the state of a free photon. It is also clear that, depending on the specific excited states being considered, the distributions and their evolution are very different. In any case, the plots in the right column, for times $\tau \ge 10$, describe the states of a free photon described in the paper [30].

10. Conclusions

It is a priori obvious that the probability of various atomic-molecular and nuclear processes initiated by photons is closely related to the spatial structure of the photons themselves. Moreover, it has long been known that electromagnetic fields consist of elementary massless Bose particles – photons, but only now has the sensitivity of the instruments used reached such a level that it allows us to conduct a variety of experiments with single photons. In this regard, an accurate quantum description of photons is of fundamental importance not only for the completeness of the foundations of quantum mechanics, but also for the further development of such a modern scientific and technical field as quantum photonics. The difficulties of constructing a quantum theory of the photon and its consistent correct description are analyzed in sufficient detail in a series of works by Bialynicki-Birula [1,21,22]. Here we will highlight two of them, which, in our opinion, are the most difficult to consider. The first is the lack of connection between the argument of the photon wave function and the photon coordinate operator, which simply does not exist. The second is the non-conservation of the number of photons, which is especially important when considering problems in a medium where processes with the absorption and emission of photons are possible.

To avoid formal complications and non-fundamental analogies, especially with massive particles, which could further complicate the quantum description of the photon, we considered the problem in the framework of the Yang-Mills equations for Abelian fields using the $SU\left(2\right)xU\left(1\right)$ gauge symmetry group. We investigated the most general formulation of the problem, in which a photon moves in an arbitrary, including in a random, environment. The 4D Minkowski space was represented as a domain with fabrics, characterized by the fact that the speed of light in 3D space is not constant and depends on the coordinates. Using the first-order Y-M vector equation (1) as an identity, we obtain a system of three second-order PDEs (3) describing the region of localization of the photon’s electromagnetic field energy in 3D space, which we call the photon’s wave function. It is important to note that the wave equations (3) taking into account the Riemann–Silberstein vector representation (10) also allow us to describe the photon spin, which is very important for studying subtle quantum phenomena such as photons entanglement, etc. Recall that the full wave function of a photon (9) is defined as the sum of the components of the vector wave function (2).

The paper examines in detail the case where the medium is homogeneous along the direction of photon propagation. For this case, an equation for the total photon wave function (15) was derived, taking into account the dispersion of the speed of light in a plane perpendicular to the photon’s direction of motion. It is shown that in the absence of dispersion of the speed of light, the wave motion of a photon is described by the Klein-Gordon-Fock equation (18).

The problem of single photon motion in a nanowaveguide is considered within the framework of the wave equation (18). The mathematical problem is reduced with high accuracy to the 2D QHO problem with variable frequencies, which in the case of equal frequencies is solved exactly (30) (see also [30]), and the set of wave functions forms an orthonormal basis in the Hilbert space (see expression (34)). In other words, a photon in a nanowaveguide is now described not only by its own frequency ${\omega }_0$ (see expression (19)) and the spin projection $\pm 1$, but also by two vibrational quantum numbers $n_1$ and $n_2$, which actually determine its spatial structure.

In this paper, the problem under consideration is significantly expanded and generalized to the case of a dissipative medium. In particular, it is assumed that a photon passing through a nanowaveguide experiences random elastic and inelastic collisions with two-level quantum dots, which leads to the probability of absorption of one photon and emission of two entangled photons by the medium. To describe this complicated process, where the law of conservation of the number of photons is clearly violated, a mathematical representation of complex probabilistic processes is used. In the framework of the Langevin-Schrödinger type SDE (35), the mathematical expectation of the photon wave function is studied in detail and a closed integral representation was constructed for it, where the integrand is the solution to a system of two second-order PDEs. This allows us to calculate and visualize the electromagnetic field distribution of a photon of a given frequency in various quantum states and verify that photons of the same frequency can have completely different structures (see Figure 4). Obviously, the efficiency of interaction of a photon with any other quantum system will be different depending on its spatial structure.

As is known, an important indicator of a dynamical system, which allows us to determine the features of its motion, is entropy. In this regard, we constructed expressions for the entropy of a single photon in the ground state using two different methods (see equations (99)-(95) and (98)), which allows us to study in detail the evolution process of a single photon in a nanowaveguide. As follows from both definitions, two different quantum states of the 1D QHO are entangled, which characterize the entanglement of two single photons through their individual entropies. Recall that the authors recently showed [32], that the entropy of a classical 1D oscillator immersed in a thermostat takes on a negative value in the form of a peaks at different stages of its evolution, which may serve as evidence for the Schrödinger’s negentropy hypothesis [39]. Based on the significant similarities between classical and quantum self-organizing systems, we hypothesize that calculating the quantum entropy of entangled photons (see Definition 5) will also allow us to identify situations where a quantum system generates negentropy. Note that this would be a fundamental proof of Schrödinger’s negentropy hypothesis for quantum self-organizing systems.

For quantum communications, Bell states play a key role, so calculating the probabilities of transitions to these states in the problem under consideration is a very important task. In this regard, we have studied the problem in detail and constructed mathematical expectations of the elements of the S-matrix of transitions for various Bell states in the form of integral representations (see (Section 7)), which also allows us to construct mathematical expressions for the probabilities of the corresponding transitions.

Note that studying the structure of the photon and its evolution could also be of great interest for cosmology in light of the Lorentz symmetry violation associated with extensions of the Standard Model, due to the non-conservation of the energy-momentum tensor of a light wave when crossing an electromagnetic background field, even if the field is constant [42]. Such studies are also relevant for obtaining additional important information about the parameters of massive astrophysical objects and the space surrounding them. In this sense, the problem of photon motion in cosmic space is ideally suited for description within the framework of the mathematical problem of quantum motion of photons in Minkowski space with fabrics.

In conclusion, we would like to note that the developed representation for the wave function of a single photon sheds light on previously unknown properties of one of the fundamental particles of nature, widely present in various physical, chemical and biological processes, which may generate new ideas for the development of fine quantum technologies and quantum photonics in general.

11. Appendix

Let us write the mathematical expectation of the wave function of a photon in the ground state (74) in coordinates $x_1$ and $x_2$ (see coordinate transformation (25)):

$$\begin{array}{c}\hfill \prod _{l=1}^2\mathbb{E}\left[{\mathcal{F}}_l\left(n_l\right|q_l,\tau ;\left\{f_l\right\})\right]={\displaystyle \frac{1}{\sqrt{\pi }}}\prod _{l=1}^2{\widehat{\mathcal{K}}}_l\left(\tau \right)e^{{\mathcal{A}}_l(u_1^{\left(l\right)},u_2^{\left(l\right)};q_l)}={\displaystyle \frac{1}{\sqrt{\pi }}}\left\{{\widehat{\mathcal{K}}}_1\left(\tau \right){\widehat{\mathcal{K}}}_2\left(\tau \right)\right\}\times \\ \hfill exp\left\{-{\displaystyle \frac{1}{4}}\left[\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_1^2-2\left({\mathcal{A}}_1-{\mathcal{A}}_2\right)x_1{x}_2-\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_2^2\right]\right\},\end{array}$$

(143)

where ${\mathcal{A}}_l\left(u_1^{\left(l\right)},u_2^{\left(l\right)};q_l\right)=\left(i{u}_1^{\left(l\right)}-u_2^{\left(l\right)}\right)q_l^2/2$, in addition, ${\widehat{\mathcal{K}}}_l\left(\tau \right)$ is an integral operator of the following form:

$${\widehat{\mathcal{K}}}_l\left(\tau \right)={\int }_{-\infty }^{+\infty }d{u}_1^{\left(l\right)}{\int }_0^{+\infty }d{u}_2^{\left(l\right)}{Q}_l(0;u_1^{\left(l\right)},u_2^{\left(l\right)},\tau ).$$

Substituting the representation (143) into the expression (122) and taking into account (123), we can perform simple calculations by the coordinates $x_1$ and $x_2$ and find the transition $\mathrm{S}$-matrix elements (see (122)) in the form of fourfold integral representation. In particular, detailed calculations for the first element lead to the following result:

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{+}\rangle }={\displaystyle \frac{1}{2\pi \sqrt{\pi }}}\underset{\tau \to +\infty }{lim}{\widehat{\mathcal{K}}}_1\left(\tau \right){\widehat{\mathcal{K}}}_2\left(\tau \right)\times \\ \hfill \underset{-\infty }{\overset{+\infty }{\int }}d{x}_1\underset{-\infty }{\overset{+\infty }{\int }}d{x}_{2}exp\left\{-{\displaystyle \frac{1}{4}}\left[\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_1^2-2\left({\mathcal{A}}_1-{\mathcal{A}}_2\right)x_1{x}_2+\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_2^2\right]\right\}\times \\ \hfill \left[{\displaystyle \frac{1}{\epsilon }}{\left({\displaystyle \frac{2}{{\sigma }_{\left(d\right)1}{\sigma }_{\left(d\right)2}}}\right)}^{1/2}exp\left\{-\sum _{l=1}^2{\displaystyle \frac{x_l^2}{2{\sigma }_{\left(d\right)l}^2}}\right\}-{\displaystyle \frac{1}{\mu }}{\left({\displaystyle \frac{2}{{\sigma }_{\left(b\right)1}{\sigma }_{\left(b\right)2}}}\right)}^{1/2}exp\left\{-\sum _{l=1}^2{\displaystyle \frac{x_l^2}{2{\sigma }_{\left(b\right)l}^2}}\right\}\right].\end{array}$$

(144)

It should be noted that for simplicity we assume that the shift of the center of the Gaussian distribution is zero, i.e. ${\mu }_{\left(g\right)l}=0$ (see expression (123)), also it is assumed that the permittivity and magnetic permeability are constant and do not depend on the coordinates.

Now we can move on to calculating the double integrals in the expression (144). The first double integral can be rewritten in the following form:

$$\begin{array}{c}\hfill {\mathcal{X}}_{\epsilon ,d}^{\left(1\right)}={\mathcal{V}}_{\epsilon ,d}\underset{-\infty }{\overset{+\infty }{\int }}d{x}_{2}exp\left\{-{\displaystyle \frac{1}{4}}\left[{\displaystyle \frac{2}{{\sigma }_{\left(d\right)2}^2}}+{\mathcal{A}}_1+{\mathcal{A}}_2\right]x_2^2\right\}\underset{-\infty }{\overset{+\infty }{\int }}exp\{-{\displaystyle \frac{1}{4}}[{\displaystyle \frac{2}{{\sigma }_{\left(d\right)1}^2}}+\\ \hfill {\mathcal{A}}_1+{\mathcal{A}}_2]x_1^2+{\displaystyle \frac{1}{2}}\left({\mathcal{A}}_1-{\mathcal{A}}_2\right)x_1{x}_2\}d{x}_1,{\mathcal{V}}_{\epsilon ,d}={\displaystyle \frac{1}{2\epsilon \pi \sqrt{\pi }}}\sqrt{{\displaystyle \frac{2}{{\sigma }_{\left(d\right)1}{\sigma }_{\left(d\right)2}}}}.\end{array}$$

(145)

By integrating over $x_1$ in expression (145), we can easily obtain the following expression:

$${\mathcal{X}}_{\epsilon ,d}^{\left(1\right)}={\overline{\mathcal{V}}}_{\epsilon ,d}\underset{-\infty }{\overset{+\infty }{\int }}d{x}_{2}exp\left\{-{\displaystyle \frac{1}{4}}\left[{\displaystyle \frac{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2}{2{\sigma }_{\left(d\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2}}+\left(2{\sigma }_{\left(d\right)2}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2\right)\right]x_2^2\right\},$$

(146)

where ${\overline{\mathcal{V}}}_{\epsilon ,d}=\sqrt{4\pi }{(2{\sigma }_{\left(d\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)}^{-1/2}{\mathcal{V}}_{\epsilon ,d}.$

Having calculated the Gaussian integral (146), we finally obtain the following expression:

$${\mathcal{X}}_{\epsilon ,d}^{\left(1\right)}=2\pi {\mathcal{V}}_{\epsilon ,d}\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(d\right)l}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2){\}}^{-1/2}.$$

(147)

The second double integral in (144) can be calculated in a similar way:

$${\mathcal{X}}_{\mu ,b}^{\left(2\right)}=2\pi {\mathcal{V}}_{\mu ,b}{\left\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(b\right)l}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right\}}^{-1/2}.$$

(148)

where ${\mathcal{V}}_{\mu ,b}={\displaystyle \frac{1}{\mu \pi \sqrt{2\pi }}}{\left\{{\sigma }_{\left(b\right)1}{\sigma }_{\left(b\right)2}\right\}}^{-1/2}$.

Finally, using (122) and (147)-(148), we can obtain an expression for the first element of the transition $\mathrm{S}$- matrix to Bell states:

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{+}\rangle }={\displaystyle \frac{1}{\sqrt{\pi }}}\underset{\tau \to +\infty }{lim}{\widehat{\mathcal{K}}}_1\left(\tau \right){\widehat{\mathcal{K}}}_2\left(\tau \right)\times \\ \hfill \{{\displaystyle \frac{1}{\epsilon }}{\left\{{\sigma }_{\left(d\right)1}{\sigma }_{\left(d\right)2}\left[{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(d\right)l}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right]\right\}}^{-1/2}\\ \hfill -{\displaystyle \frac{1}{\mu }}{\left\{{\sigma }_{\left(b\right)1}{\sigma }_{\left(b\right)2}\left[{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\prod _{l=1}^2(2{\sigma }_{\left(b\right)l}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right]\right\}}^{-1/2}\}.\end{array}$$

(149)

Based on the symmetry between the elements of the transition $\mathrm{S}$ -matrix (see (122)), using (149), we can write the expression for the matrix element ${\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\varphi }_{+}\rangle }$ by simply changing the sign $“-”$ to $“+”$ before the second term in curly brackets $“\left\{\dots \right\}”$.

Let us consider the matrix element ${\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{-}\rangle }$, which will be conveniently written in the following form:

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{-}\rangle }={\displaystyle \frac{-i}{2\pi \sqrt{\pi \epsilon \mu }}}\underset{\tau \to +\infty }{lim}{\widehat{\mathcal{K}}}_1\left(\tau \right){\widehat{\mathcal{K}}}_2\left(\tau \right)\times \\ \hfill \underset{-\infty }{\overset{+\infty }{\int }}d{x}_1\underset{-\infty }{\overset{+\infty }{\int }}d{x}_{2}exp\left\{-{\displaystyle \frac{1}{4}}\left[\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_1^2-2\left({\mathcal{A}}_1-{\mathcal{A}}_2\right)x_1{x}_2+\left({\mathcal{A}}_1+{\mathcal{A}}_2\right)x_2^2\right]\right\}\times \\ \hfill \left[\sqrt{{\displaystyle \frac{2}{{\sigma }_{\left(d\right)1}{\sigma }_{\left(b\right)2}}}}exp\left\{-{\displaystyle \frac{x_1^2}{2{\sigma }_{\left(d\right)1}^2}}-{\displaystyle \frac{x_2^2}{2{\sigma }_{\left(b\right)2}^2}}\right\}+\sqrt{{\displaystyle \frac{2}{{\sigma }_{\left(b\right)1}{\sigma }_{\left(d\right)2}}}}exp\left\{-{\displaystyle \frac{x_1^2}{2{\sigma }_{\left(b\right)1}^2}}-{\displaystyle \frac{x_2^2}{2{\sigma }_{\left(d\right)2}^2}}\right\}\right].\end{array}$$

(150)

By performing simple calculations in (150), we eventually obtain the following expression:

$$\begin{array}{c}\hfill {\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\psi }_{-}\rangle }={\displaystyle \frac{-i}{\sqrt{\pi \epsilon \mu }}}\underset{\tau \to +\infty }{lim}{\widehat{\mathcal{K}}}_1\left(\tau \right){\widehat{\mathcal{K}}}_2\left(\tau \right)\times \\ \hfill \{{\left[{\sigma }_{\left(d\right)1}{\sigma }_{\left(b\right)2}\left\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+(2{\sigma }_{\left(d\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)(2{\sigma }_{\left(b\right)2}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right\}\right]}^{-1/2}+\\ \hfill {\left[{\sigma }_{\left(b\right)1}{\sigma }_{\left(d\right)2}\left\{{({\mathcal{A}}_1-{\mathcal{A}}_2)}^2+\left(2{\sigma }_{\left(b\right)1}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right(2{\sigma }_{\left(d\right)2}^{-2}+{\mathcal{A}}_1+{\mathcal{A}}_2)\right\}\right]}^{-1/2}\}.\end{array}$$

(151)

Note that in a similar way we can calculate the transition matrix element ${\mathrm{S}}_{\mathbb{E}\left[\mathcal{F}\right]\to |{\varphi }_{-}\rangle }$.

Thus, we have calculated all the necessary transition $\mathrm{S}$-matrix elements, which allows us to construct the transition probabilities to Bell states (see expressions (124)). Further calculation of transitions to Bell states can only be carried out by numerical simulation.