Coupled Harmonic Oscillator in a System of Free Particles

Abstract

The coupled quantum harmonic oscillator is one of the most researched and important model systems in quantum optics and quantum informatics. This system is often investigated for quantum entanglement in the environment. As a result, such studies are complex and can only be carried out using numerical methods that do not reveal the general pattern of such systems. In this work, the external environment is considered to be two independent particles interacting with coupled harmonic oscillators. It is shown that such a system has an exact analytical solution to the dynamic Schrödinger equation. The analysis of this solution is carried out, and the main parameters of this system are revealed. The solutions obtained can be used to study more complex systems and their quantum entanglement.

1. Introduction

A coupled harmonic oscillator is the most important model system in quantum optics and computer science. For example, a model of linear beam splitter in quantum optics can be represented by two coupled harmonic oscillators [1]. This model is also used to explain the problem of photosynthesis based on quantum entangled states [2,3,4]. Study of properties of coupled harmonic oscillators, mainly quantum entanglement, is a separate direction in quantum physics. First of all, this is due to the fact that bound harmonic oscillators are a good model of real physical objects. Such objects may include thermal vibrations of coupled atoms, photons in cavities, ions in traps, and many more. The study of coupled harmonic oscillators is one of the main theoretical methods to investigate quantum decoherence [5,6,7]. It is known that quantum decoherence leads to violation of quantum coherence of a system. When decoherence occurs, the system itself has classical features that correspond to the information available in the environment. This means that the system is mixed or entangled with the environment. Reducing the quantum decoherence of a system when it interacts with a classical system is one of the important problems in quantum informatics on the way to creating quantum computers. Various methods are used to reduce decoherence, such as isolating the quantum system (using extremely low temperatures and high vacuum) [8,9], introducing error-resistant codes into quantum computing [9,10], the use of intense sources of quantum entangled particles [7,11,12,13], or the application of various complex mechanisms (e.g., Reference [14]). Theoretical study of quantum decoherence is mainly related to the study of quantum entanglement of a system in its interaction with a system which is a model of a classical medium [7,15,16,17]. Usually, oscillators are used as quantum entangled particles, which interact with the system of oscillators (thermal bath) modeling the classical medium (e.g., References [16,17,18] ). This choice is justified because oscillators are the simplest model for studying complex systems. In order to reduce quantum decoherence, one usually selects various model parameters, at which used measures of quantum entanglement become minimal. It should be specified that such parameters are chosen by numerically solving the Schrödinger equation for the chosen model. Moreover, differential equations in such models are simplified to a Quantum Master Equation (QME) (including Lindblad equation, Redfield equation, etc.) to enable numerical solution of the system under study. Of course, it is very difficult to find the best and correct parameters of the chosen model, when decoherence will be minimal, in the numerical solution. Besides, it is rather difficult to study the model itself, i.e., what can be changed in the model itself so that decoherence would be minimal. Such analysis requires analytical solutions of the chosen model, in which it is possible to allocate the basic parameters, at which the decoherence will be minimal.

In this paper, a solution of a dynamic Schrödinger equation for a model of bound harmonic oscillators interacting with 2 free particles is presented. The obtained solution has an analytical solution, which shows the terms responsible for decoherence. The obtained solution can be used to analyze quantum entanglement and decoherence of bound harmonic oscillators.

2. Model and Solution

Consider the model of coupled harmonic oscillators interacting with two independent particles (see Figure 1). The Hamiltonian of such a system is $\widehat{H}={\widehat{H}}_1+{\widehat{H}}_2+{\widehat{H}}_{int}$, where

$${\widehat{H}}_1=\frac{1}{2}\left(\frac{1}{m_1}{\widehat{p}}_1^2+\frac{1}{m_2}{\widehat{p}}_2^2+a_1{x}_1^2+a_2{x}_2^2+a_3{x}_1{x}_2\right),$$

(1)

$${\widehat{H}}_2=\frac{1}{2}\left(\frac{1}{m_3}{\widehat{p}}_3^2+\frac{1}{m_4}{\widehat{p}}_4^2\right),$$

(2)

$${\widehat{H}}_{int}=\left(a_4{x}_1+a_5{x}_2\right){\widehat{p}}_3+\left(a_6{x}_1+a_7{x}_2\right){\widehat{p}}_4,$$

(3)

where ${\widehat{p}}_i=-i\hslash \partial /\partial x_i$ is the momentum operator, $m_i$ is the mass for i particle, respectively, $a_1$ and $a_2$ are stiffness coefficients of 1 and 2 oscillators, respectively, $a_3$ is the coupling factor between 1 and 2 oscillators, and $a_i$ (i = 4, 5, 6, 7) are the factors determining the value of interaction of free particles with each of the oscillators.

In Equation (1), ${\widehat{H}}_1$ is the Hamiltonian of two bound oscillators of mass, $m_1,m_2$, respectively. In Equation (2), ${\widehat{H}}_2$ is the Hamiltonian of two free particles with masses, $m_3,m_4$, respectively. In Equation (3), ${\widehat{H}}_{int}$ is the Hamiltonian responsible for the interaction of bound oscillators with two particles. It should be added that this type of interaction ${\widehat{H}}_{int}$ is primarily defined in possible application cases [1,19,20], as well as the possibility to further obtain an analytical solution. This type of interaction can be implemented, for example, in a two-port lossy waveguide beam splitter. Indeed, the two modes of the electromagnetic field can be represented by the Hamiltonian $H_1$, e.g., as in References [1,21]. Each mode of the electromagnetic field interacts with the environment (Hamiltonian $H_2$) through the electromagnetic interaction $H_3$, which leads to an interaction operator with a charged particle in the well-known form ${(\mathbf{p}-e/c\mathbf{A})}^2/\left(2m\right)$, where $\mathbf{p}$ is the momentum operator of a particle of mass m and charge e, and $\mathbf{A}$ is the electromagnetic field potential vector. In this case, the interaction of two modes with particles of mass $m_3$ and $m_4$ will be described by the Hamiltonian equation $\widehat{H}$, similar to References [1,19,20] (of course, considering $m_1$ and $m_2$ as just some constants, not masses).

A system consisting of two unbound particles will be classical only when the masses $m_3,m_4$ have large values (more precisely, the De Broglie wavelength for particles with masses $m_3,m_4$ should ${\lambda }_B\ll a$, where a is the characteristic size of the system in question). Since, as will be shown later, the model under consideration has an exact analytic solution, the masses $m_3,m_4$ can be chosen arbitrarily. That is why we will search for the solution in general form, not limiting ourselves to large values of $m_3,m_4$. We should also add that this choice of the classical system is one of a large number of variants of models that are now being actively studied. For example, References [16,17] investigated quantum entanglement of coupled oscillators in thermal bath. Harmonic oscillators interacting with coupled harmonic oscillators acted as thermal baths. A more detailed description of such systems can be found in References [22,23,24]. However, only in this model was it possible to obtain an analytical solution.

For further consideration of the presented model, it is convenient to make substitutions of variables in the system under consideration:

$$\begin{array}{c}\hfill {\tilde{x}}_1=x_1{\left(\frac{m_1}{m_2}\right)}^{1/4},{\tilde{x}}_2=x_2{\left(\frac{m_2}{m_1}\right)}^{1/4},{\tilde{x}}_3=x_3{\left(\frac{m_3}{m_4}\right)}^{1/4},{\tilde{x}}_4=x_4{\left(\frac{m_4}{m_3}\right)}^{1/4},\\ \hfill {\tilde{a}}_1=a_1{\left(\frac{m_2}{m_1}\right)}^{1/2},{\tilde{a}}_2=a_2{\left(\frac{m_1}{m_2}\right)}^{1/2},{\tilde{a}}_4=a_4{\left(\frac{m_2{m}_3}{m_1{m}_4}\right)}^{1/4},{\tilde{a}}_5=a_5{\left(\frac{m_1{m}_3}{m_2{m}_4}\right)}^{1/4},\\ \hfill {\tilde{a}}_6=a_6{\left(\frac{m_2{m}_4}{m_1{m}_3}\right)}^{1/4},{\tilde{a}}_7=a_7{\left(\frac{m_1{m}_4}{m_2{m}_3}\right)}^{1/4},M_1={\left(m_1{m}_2\right)}^{1/2},M_2={\left(m_3{m}_4\right)}^{1/2}.\end{array}$$

(4)

Next, it is easier to use expressions without the tilde $\tilde{x}$, so we remove this sign from the top. We end up with (3) in the same form, and (1) and (2) in simpler forms:

$${\widehat{H}}_1=\frac{1}{2M_1}\left({\widehat{p}}_1^2+{\widehat{p}}_2^2\right)+\frac{1}{2}\left(a_1{x}_1^2+a_2{x}_2^2+a_3{x}_1{x}_2\right),$$

(5)

$${\widehat{H}}_2=\frac{1}{2M_2}\left({\widehat{p}}_3^2+{\widehat{p}}_4^2\right).$$

(6)

Next, consider the nonstationary Schrödinger equation $\widehat{H}\Psi =i\hslash \frac{\partial \Psi }{\partial t}$. Consider the case of initial conditions, which is often studied in the literature, e.g., in References [1,25]. At the initial moment of time, $t=0$, the bound oscillators have not interacted with the particles, i.e., $\widehat{H}\left(t=0\right)={\widehat{H}}_1+{\widehat{H}}_2$. In this case, we know that, if we make variable substitutions in the form $y_1={(a_1^{{}^{\prime }}{M}_1/{\hslash }^2)}^{1/4}(x_{1}cos\theta -x_{2}sin\theta ),y_2={(a_2^{{}^{\prime }}{M}_1/{\hslash }^2)}^{1/4}(x_{1}sin\theta +x_{2}cos\theta )$, where $tan\theta =\frac{ϵ}{\left|ϵ\right|}\sqrt{{ϵ}^2+1}-ϵ$, where $ϵ=\frac{a_2-a_1}{a_3}$, $a_1^{{}^{\prime }}=a_1-\frac{a_3}{2}tan\theta ,a_2^{{}^{\prime }}=a_2+\frac{a_3}{2}tan\theta$ (see e.g., Reference [1]), the Hamiltonian (5) can be reduced to a diagonal form. As a result, at $t=0$, the wave function of the system in $|{\Psi }_{n_0,m_0,{\mathbf{k}}_0}(y_1,y_2,\mathbf{r},t=0)\rangle =c_{n_0}{c}_{m_0}C{e}^{-\frac{1}{2}{y}_1^2}{e}^{-\frac{1}{2}{y}_2^2}{H}_{n_0}\left(y_1\right)H_{m_0}\left(y_2\right)e^{i{\mathbf{k}}_0\mathbf{r}}$, where $c_{n_0}=1/\sqrt{\sqrt{\pi }{n}_0!2^{n_0}}$, $c_{m_0}=$$1/\sqrt{\sqrt{\pi }{m}_0!2^{m_0}}$,$n_0,m_0$ are quantum numbers, $H_n\left(z\right)$ are Hermite polynomials, $C={\left(L_3{L}_4\right)}^{-1/2}$ is the normalization factor for particles with masses $m_3,m_4$ (constrained in $L_3$ and $L_4$ space, respectively), ${\mathbf{k}}_0=k_3\mathbf{i}+k_4\mathbf{j}$ is wave vector of particles ($\mathbf{i}$ and $\mathbf{j}$ are unit perpendicular vectors introduced to simplify writing), and $\mathbf{r}=x_3\mathbf{i}+x_4\mathbf{j}$. It should be added that quantum entanglement of bound oscillators, in the case of the initial state (without wave function $e^{i{\mathbf{k}}_0\mathbf{r}}$ ) has been studied in detail in many papers, for example, in ref. [25], and does not need a detailed analysis.

As a result, the Hamiltonian of the system, using the new variables $y_1,y_2$, will be

$$\begin{array}{c}\hfill \widehat{H}=\frac{\hslash {\omega }_1^{{}^{\prime }}}{2}\left(-\frac{{\partial }^2}{\partial y_1^2}+y_1^2\right)+\frac{\hslash {\omega }_2^{{}^{\prime }}}{2}\left(-\frac{{\partial }^2}{\partial y_2^2}+y_2^2\right)+\left(y_1\mathbf{A}+y_2\mathbf{B}\right)\widehat{\mathbf{p}}+\frac{{\widehat{p}}^2}{2M_2},\end{array}$$

(7)

where ${\omega }_1^{{}^{\prime }}=\sqrt{a_1^{{}^{\prime }}/M_1}$, ${\omega }_2^{{}^{\prime }}=\sqrt{a_2^{{}^{\prime }}/M_1}$,

$$\mathbf{A}=\sqrt{\hslash }{\left(a_1^{{}^{\prime }}{M}_1\right)}^{-\frac{1}{4}}\left((a_{4}cos\theta -a_{5}sin\theta )\mathbf{i}+(a_{6}cos\theta -a_{7}sin\theta )\mathbf{j}\right),$$

$$\mathbf{B}=\sqrt{\hslash }{\left(a_2^{{}^{\prime }}{M}_1\right)}^{-\frac{1}{4}}\left((a_{4}sin\theta +a_{5}cos\theta )\mathbf{i}+(a_{6}sin\theta +a_{7}cos\theta )\mathbf{j}\right),$$

and $\widehat{\mathbf{p}}={\widehat{p}}_3\mathbf{i}+{\widehat{p}}_4\mathbf{j}$. In order to find a solution of the nonstationary Schrödinger equation with Hamiltonian (7), it is necessary to diagonalize the Hamiltonian (7) for the variables $y_1,y_2$. In this case, the nonstationary Schrödinger equation will be in the form ${\widehat{H}}^{{}^{\prime }}{\Psi }^{{}^{\prime }}=i\hslash \frac{\partial {\Psi }^{{}^{\prime }}}{\partial t}$, where ${\widehat{H}}^{{}^{\prime }}=\widehat{S}\widehat{H}{\widehat{S}}^{-1}$ and ${\Psi }^{{}^{\prime }}=\widehat{S}\Psi$. The operator $\widehat{S}$ is some operator transforming the Hamiltonian (7) to a diagonal form for the variables $y_1,y_2$. To find it, we represent $\widehat{S}$ in the form $\widehat{S}=e^{\widehat{\alpha }\frac{\partial }{\partial y_1}}{e}^{\widehat{\beta }\frac{\partial }{\partial y_2}}$, where $\widehat{\alpha }$ and $\widehat{\beta }$ are some operators independent of the variables $y_1,y_2$, but which can depend on $\frac{\partial }{\partial }\partial x_1,\frac{\partial }{\partial x_2}$. Using the well-known Magnus expansion

$$e^{\widehat{X}}\widehat{Y}{e}^{-\widehat{X}}=\widehat{Y}+\left[\widehat{X},\widehat{Y}\right]+\frac{1}{2!}\left[\widehat{X},\left[\widehat{X},\widehat{Y}\right]\right]+\frac{1}{3!}\left[\widehat{X},\left[\widehat{X},\left[\widehat{X},\widehat{Y}\right]\right]\right]+\dots ,$$

it is not difficult to show that, to diagonalize the Hamiltonian (7):

$$\begin{array}{c}\hfill \widehat{\alpha }=-\frac{\mathbf{A}\widehat{\mathbf{p}}}{\hslash {\omega }_1^{{}^{\prime }}},\widehat{\beta }=-\frac{\mathbf{B}\widehat{\mathbf{p}}}{\hslash {\omega }_2^{{}^{\prime }}}.\end{array}$$

(8)

We end up with a Hamiltonian in the form

$$\begin{array}{c}\hfill {\widehat{H}}^{{}^{\prime }}=\frac{\hslash {\omega }_1^{{}^{\prime }}}{2}\left(-\frac{{\partial }^2}{\partial y_1^2}+y_1^2\right)+\frac{\hslash {\omega }_2^{{}^{\prime }}}{2}\left(-\frac{{\partial }^2}{\partial y_2^2}+y_2^2\right)-\frac{{\left(\mathbf{A}\widehat{\mathbf{p}}\right)}^2}{2\hslash {\omega }_1^{{}^{\prime }}}-\frac{{\left(\mathbf{B}\widehat{\mathbf{p}}\right)}^2}{2\hslash {\omega }_2^{{}^{\prime }}}+\frac{{\widehat{p}}^2}{2M_2}.\end{array}$$

(9)

As a result, the desired solution can be represented as $\Psi ={\widehat{S}}^{-1}{\Psi }^{{}^{\prime }}$. To correctly account for the initial conditions, the solution will be sought in the form $\Psi =e^{-\frac{i}{\hslash }\widehat{H}t}{\Psi }_{n_0,m_0,{\mathbf{p}}_0}\left(t=0\right)$, where ${\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)$ is the wave function of the initial state presented above. We then obtain ${\Psi }^{{}^{\prime }}=\widehat{S}{e}^{-\frac{i}{\hslash }\widehat{H}t}{\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)$. Using the easily verifiable expression $\widehat{S}{e}^{-\frac{i}{\hslash }\widehat{H}t}{\widehat{S}}^{-1}=e^{-\frac{i}{\hslash }{\widehat{H}}^{{}^{\prime }}t}$, we get ${\Psi }^{{}^{\prime }}=e^{-\frac{i}{\hslash }{\widehat{H}}^{{}^{\prime }}t}\widehat{S}{\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)$. Then, we decompose the expression $\widehat{S}{\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)$ in the basis of the eigenfunctions ${\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}$ of the Hamiltonian (9) for the variables $y_1,y_2$, i.e., $\widehat{S}{\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)={\sum }_{n,m}{A}_{n,m}^{n_0,m_0}{\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}$, where

$$\begin{array}{c}\hfill {\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}=c_n{c}_{m}C{e}^{-\frac{1}{2}{y}_1^2}{e}^{-\frac{1}{2}{y}_2^2}{H}_n\left(y_1\right)H_m\left(y_2\right)e^{i\frac{{\mathbf{p}}_0\mathbf{r}}{\hslash }},\\ \hfill A_{n,m}^{n_0,m_0}=\langle {\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}|\widehat{S}|{\Psi }_{n_0,m_0,{\mathbf{p}}_0}(t=0)\rangle .\end{array}$$

(10)

Then, we get ${\Psi }^{{}^{\prime }}={\sum }_{n,m}{e}^{-\frac{i}{\hslash }{E}_{n,m,{\mathbf{p}}_0}t}{A}_{n,m}^{n_0,m_0}{\Psi }_{n,m}^{{}^{\prime }}$, where $E_{n,m,{\mathbf{p}}_0}$ is the Hamiltonian energy eigenvalue (9), which will be

$$\begin{array}{c}\hfill E_{n,m,{\mathbf{p}}_0}=\hslash {\omega }_1^{{}^{\prime }}\left(n+\frac{1}{2}\right)+\hslash {\omega }_2^{{}^{\prime }}\left(m+\frac{1}{2}\right)+\frac{p_0^2}{2M_2}-\frac{{\left(\mathbf{A}{\mathbf{p}}_0\right)}^2}{2\hslash {\omega }_1^{{}^{\prime }}}-\frac{{\left(\mathbf{B}{\mathbf{p}}_0\right)}^2}{2\hslash {\omega }_2^{{}^{\prime }}}.\end{array}$$

(11)

Then, the wave function is $\Psi ={\sum }_{n,m}{e}^{-\frac{i}{\hslash }{E}_{n,m,{\mathbf{p}}_0}t}{A}_{n,m}^{n_0,m_0}{\widehat{S}}^{-1}{\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}$. Here, we can see that ${\widehat{S}}^{-1}{\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}=e^{-\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}}\frac{\partial }{\partial y_1}}{e}^{-\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}\frac{\partial }{\partial y_2}}{\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}$. You can also see that the operators in the exponent are offset operators on coordinates $y_1$ and $y_2$, respectively. As a result, ${\widehat{S}}^{-1}{\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}={\Psi }_{n,m,{\mathbf{p}}_0}^{{}^{\prime }}\left(y_1+\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}},y_2+\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}\right)$. Then, we get

$$\begin{array}{c}\hfill \Psi =\sum _{n,m}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{\Psi }_{n,m}^{{}^{\prime }}\left(y_1+\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}},y_2+\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}\right),\\ \hfill E_{n,m}=\hslash {\omega }_1^{{}^{\prime }}\left(n+\frac{1}{2}\right)+\hslash {\omega }_2^{{}^{\prime }}\left(m+\frac{1}{2}\right).\end{array}$$

(12)

In Equation (12), the term $C{e}^{i\frac{{\mathbf{p}}_0\mathbf{r}}{\hslash }}$ is omitted since it is not of interest for calculating probabilities, respectively, ${\Psi }_{n,m}^{{}^{\prime }}$ is the wave function (10), but without $C{e}^{i\frac{{\mathbf{p}}_0\mathbf{r}}{\hslash }}$. In addition, $E_{n,m,{\mathbf{p}}_0}$ is replaced by $E_{n,m}$, since this replacement is not essential for probability calculations, since the wave function only changes by an immaterial phase. In other words, the wave function (12) takes into account only states of bound oscillators, and we are not interested in states (energy and momentum) of particles interacting with these oscillators (summed over all finite states). In principle, it is possible to limit ourselves to Equation (12) for the solution of the dynamic Schrödinger equation, with the condition of finding the coefficient $A_{n,m}^{n_0,m_0}$ calculated by the Formula (10). In reality, we are interested in the probabilities of certain processes, which have a very definite physical meaning. To find such a probability, it is necessary to decompose the wave function (12) by eigenfunctions of the unperturbed system, i.e. the state under study. Then, the square of the modulus of the expansion coefficient will have the meaning of detecting the system at time t in the given state. Based on the initial condition of our problem, the state under study will be in the form $|{\Psi }_{k,p}(y_1,y_2)\rangle =|{\Psi }_k\left(y_1\right)\rangle |{\Psi }_p\left(y_2\right)\rangle =c_k{c}_p{e}^{-\frac{1}{2}{y}_1^2}{e}^{-\frac{1}{2}{y}_2^2}{H}_k\left(y_1\right)H_p\left(y_2\right)$, where $c_k=1/\sqrt{\sqrt{\pi }k!2^k}$, $c_p=1/\sqrt{\sqrt{\pi }p!2^p}$, $k,p$ are quantum numbers. So, we get decompositions in the form ${\Psi }_{n,m}^{{}^{\prime }}\left(y_1+\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}},y_2+\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}\right)={\sum }_{k,p}{B}_{k,p}^{n,m}|{\Psi }_{k,p}(y_1,y_2)\rangle$, where

$$\begin{array}{c}\hfill B_{k,p}^{n,m}=\langle {\Psi }_{k,p}(y_1,y_2){\Psi }_{n,m}^{{}^{\prime }}\left(y_1+\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}},y_2+\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}\right)\rangle .\end{array}$$

(13)

You can see that the coefficient $B_{k,p}^{n,m}$ coincides with the coefficient $A_{k,p}^{n,m}$ at ${\mathbf{p}}_0\to -{\mathbf{p}}_0$, i.e., $B_{k,p}^{n,m}=A_{k,p}^{n,m}({\mathbf{p}}_0\to -{\mathbf{p}}_0)$. It is not difficult to obtain using the table values of integrals [26]

$$\begin{array}{c}\hfill B_{k,p}^{n,m}={(-1)}^{(k-n)\theta (k-n)}{(-1)}^{(p-m)\theta (p-m)}{\sqrt{\frac{k!}{n!}}}^{sgn(n-k)}{\sqrt{\frac{p!}{m!}}}^{sgn(m-p)}{e}^{-\frac{1}{4}\left({\alpha }^2+{\beta }^2\right)}\\ \hfill \times {\left(\frac{\alpha }{\sqrt{2}}\right)}^{|n-k|}{\left(\frac{\beta }{\sqrt{2}}\right)}^{|m-p|}{L}_{(n+k-|n-k\left|\right)/2}^{|n-k|}\left({\left(\frac{\alpha }{\sqrt{2}}\right)}^2\right)L_{(m+p-|m-p\left|\right)/2}^{|m-p|}\left({\left(\frac{\beta }{\sqrt{2}}\right)}^2\right),\\ \hfill A_{n,m}^{n_0,m_0}={(-1)}^{(n_0-n)\theta (n_0-n)}{(-1)}^{(m_0-m)\theta (m_0-m)}{\sqrt{\frac{n!}{n_0!}}}^{sgn(n_0-n)}{\sqrt{\frac{m!}{m_0!}}}^{sgn(m_0-m)}{e}^{-\frac{1}{4}\left({\alpha }^2+{\beta }^2\right)}\\ \hfill \times {\left(\frac{\alpha }{\sqrt{2}}\right)}^{|n_0-n|}{\left(\frac{\beta }{\sqrt{2}}\right)}^{|m_0-m|}{L}_{(n_0+n-|n_0-n\left|\right)/2}^{|n_0-n|}\left({\left(\frac{\alpha }{\sqrt{2}}\right)}^2\right)L_{(m_0+m-|m_0-m\left|\right)/2}^{|m_0-m|}\left({\left(\frac{\beta }{\sqrt{2}}\right)}^2\right),\end{array}$$

(14)

where $L_a^b\left(z\right)$ are Laguerre polynomials, $\theta \left(x\right)$ is the Haveside tetta function, $sgn\left(x\right)$ is the signum function, and the coefficients $\alpha =\frac{\mathbf{A}{\mathbf{p}}_0}{\hslash {\omega }_1^{{}^{\prime }}}$, $\beta =\frac{\mathbf{B}{\mathbf{p}}_0}{\hslash {\omega }_2^{{}^{\prime }}}$. We end up with a solution to the dynamic Schrödinger equation in the form

$$\begin{array}{c}\hfill \Psi =\sum _{n,m,k,p}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{B}_{k,p}^{n,m}|{\Psi }_{k,p}(y_1,y_2)\rangle .\end{array}$$

(15)

It is also convenient to represent the Equation (15) as

$$\begin{array}{c}\hfill \Psi =\sum _{n,m,k,p}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{B}_{k,p}^{n,m}|{\Psi }_k\left(y_1\right)\rangle |{\Psi }_p\left(y_2\right)\rangle =\sum _{n,m,k,p}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{B}_{k,p}^{n,m}|k\rangle |p\rangle .\end{array}$$

(16)

The probability of finding the system in the state $|k\rangle |p\rangle$ is then

$$\begin{array}{c}\hfill P_{k,p}={\left|\sum _{n,m}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{B}_{k,p}^{n,m}\right|}^2.\end{array}$$

(17)

From Equation (16), as well as (17), we can see that, at $\alpha \ll 1$ and $\beta \ll 1$, the wave function $\Psi =e^{-\frac{i}{\hslash }{E}_{n_0,m_0}t}|{\Psi }_{n_0,m_0}(y_1,y_2)\rangle$, that is, the initial state wave function. Thus, we have two main coefficients $\alpha$ and $\beta$, indicating how quantum decoherence of coupled harmonic oscillators occurs. The coefficient $\alpha$ indicates quantum decoherence of the state $|{\Psi }_{n_0}\left(y_1\right)\rangle$, and $\beta$ indicates decoherence of the state $|{\Psi }_{m_0}\left(y_2\right)\rangle$. How quantum entanglement of bound harmonic oscillators in the state $|{\Psi }_{n_0}\left(y_1\right)\rangle |{\Psi }_{m_0}\left(y_2\right)\rangle$ depends on the main parameters of the system has been well studied previously; see, e.g., Reference [25]. Next, we present the results of the calculations, in Figure 2, probabilities $P_{k,p}$ as a functional dependence on time (more precisely ${\omega }_1^{{}^{\prime }}t$ and ${\omega }_2^{{}^{\prime }}t$ ) for given parameters $\alpha =1$ and $\beta =1$ and initial states $n_0=m_0=2$.

Figure 2 shows that the influence of the parameters $\alpha ,\beta$ is significant. First, this is due to the fact that, if $\alpha ,\beta$ is small, i.e., if we assume that $\alpha =\beta \to 0$, then the nonzero probability will be only at $n_0=k,m_0=p$. In other words, at $\alpha =\beta \to 0$, the number of states is conserved as $n_0+m_0=const$ during the evolution of the quantum system (this is shown in References [1,25]). For example, for a waveguide beam splitter, which can be described as a coupled harmonic oscillator, this means that the number of photons in the lossless system is always the same, i.e., as in Reference [1], is conserved. In our case, at $\alpha =\beta =1$, we do not observe this conservation, i.e., $k+p\ne const$. This means that part of the oscillator’s energy is absorbed by the external environment, or vice versa. It should be added that, in the general case, there is not only absorption of energy by the external environment but also receiving of energy by the oscillators from the external environment, i.e., $k+p>n_0+m_0$.

Of course, the general analysis of quantum entanglement in the presented model is a separate topic for study. Nevertheless, it is possible to consider the special case when $\alpha$ and $\beta$ are small. Before that, we will obtain expressions which can be used in the future for calculations of quantum entanglement. Let us represent the wave function (16) as a expansion on pure separable states $|i\rangle ={\Psi }_i\left(x_1\right)$ and $|l\rangle ={\Psi }_l\left(x_2\right)$ of harmonic oscillators. This requires the wave function of the pure non-separable state $|k\rangle |p\rangle$ (represented in Equation (16)) to decompose the pure separable states $|i\rangle$$|l\rangle$. Such a decomposition was carried out in Reference [25]; we obtain

$$\begin{array}{c}\hfill |k\rangle |p\rangle =\sum _{i=0}^{k+p}{C}_{k,p}^{i,k+p-i}|i\rangle |k+p-i\rangle ,\\ \hfill C_{k,p}^{i,l}=\frac{{\mu }^{i+k}\sqrt{p!k!}}{{(1+{\mu }^2)}^{\frac{k+p}{2}}\sqrt{i!l!}}{P}_k^{(-(1+p+k),p-i)}\left(-\frac{2+{\mu }^2}{{\mu }^2}\right),\end{array}$$

(18)

where $P_a^{(b,c)}\left(x\right)$ is the Jacobi polynomial, and $\mu =tan\theta$. In Equation (18), as was shown in Reference [11], in calculating the integral, the condition $k+p=i+l$ is fulfilled. Taking into account the above, the wave function (16) can be represented as a decomposition over pure separable states in the form

$$\begin{array}{c}\hfill \Psi =\sum _{i,l}{D}_{i,l}|i\rangle |l\rangle ,D_{i,l}=\sum _{k,n,m}{C}_{k,i+l-k}^{i,l}{e}^{-\frac{i}{\hslash }{E}_{n,m}t}{A}_{n,m}^{n_0,m_0}{B}_{k,i+l-k}^{n,m}.\end{array}$$

(19)

From Equation (19) follows an interesting property: when $\alpha =\beta \to 0$, it turns out that $D_{i,l}=C_{n_0,m_0}^{i,n_0+m_0-i}$ (it turns out that $l=m_0+n_0-i$). As a result, with $\alpha =\beta \to 0$, it turns out that $\Psi ={\sum }_{i=0}^{n_0+m_0}{C}_{n_0,m_0}^{i,n_0+m_0-i}|i\rangle |n_0+m_0-i\rangle$. Then, $C_{n_0,m_0}^{i,n_0+m_0-i}$ is the known Schmidt mode in the Schmidt decomposition. In this case, as we know, quantum entanglement can be calculated using von Neumann entropy $S_N=-{\sum }_{i=0}^{n_0+m_0}{\lambda }_{i}ln\left({\lambda }_i\right)$, where ${\lambda }_i={\left|C_{n_0,m_0}^{i,n_0+m_0-i}\right|}^2$ [25,27]. Quantum entanglement in this particular case is studied in the work [25]. For quantum entanglement calculations in the general case, where $\alpha \ne 0$ and $\beta \ne 0$, von Neumann entropy is not a suitable measure of quantum entanglement. In this case, we can use such measures of quantum entanglement as Concurrence or Negativity [28].

3. Discussion and Conclusions

Thus, we obtained an analytical expression for the wave function of a bound harmonic oscillator under the action of two free particles (16), as well as the probability of detecting bound oscillators in unperturbed states $|k\rangle |p\rangle$; see Equation (17). Two main parameters of the system responsible for quantum decoherence have been found: $\alpha$ and $\beta$. The larger parameters $\alpha$ and $\beta$ are, the stronger the bound harmonic oscillators interact with particles. At such interaction, the initial state $|n_0\rangle |m_0\rangle$ is destroyed, and the system itself has features which correspond to information available in the environment. Indeed, at $\alpha \gtrsim 1$ or $\beta \gtrsim 1$ or $(\alpha ,\beta )\gtrsim 1$, the states $|k\rangle |p\rangle$ may differ significantly from the initial $|n_0\rangle |m_0\rangle$, that is, the oscillator may be in many states, depending on the probability $P_{k,p}$, which, in turn, depends on the degree of interaction with the environment. The study of obtained expressions on quantum entanglement is a separate topic, which may be investigated in the future. It should be added that the coefficients $\alpha$ and $\beta$, responsible for the degree of interaction of coupled oscillators with the medium, depend on many parameters. First of all, it depends on energy of free particles or their momentum ${\mathbf{p}}_0$; secondly, it depends on the degree of interaction of particles with the oscillator these are coefficients $a_4,a_5,a_6,a_7$. It also depends on the degree of coupling of the oscillators, i.e., the coefficient $\theta$, and on their frequencies, ${\omega }_1,{\omega }_2$. In other words, the coefficients $\alpha$ and $\beta$ have complex dependence on many parameters of the system under study, but, as a result, the probability $P_{k,p}$ depends on 4 dimensionless quantities ${\omega }_1^{{}^{\prime }}t,{\omega }_2^{{}^{\prime }}t,\alpha ,\beta$ (without taking into account the initial states $n_0,m_0$).

The obtained expressions have a simple analytical form, which contributes to a simple analysis. This model has good prospects for further study and identification of the basic regularities of quantum entanglement of coupled oscillators. In addition, this model can be extended to the case of interaction of bound oscillators with many free particles.