Научная статья на тему 'Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time'

Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time Текст научной статьи по специальности «Физика»

CC BY
62
13
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Компьютерная оптика
Scopus
ВАК
RSCI
ESCI
Область наук
Ключевые слова
non-Markovian evolution / bosonic reservoir / stochastic field / two-level system.

Аннотация научной статьи по физике, автор научной работы — V.A. Mikhailov, N.V. Troshkin

In this paper we investigate non-Markovian evolution of a two-level system (qubit) in a bosonic bath under influence of an external classical fluctuating environment. The interaction with the bath has the Lorentzian spectral density, and the fluctuating environment (stochastic field) is represented by a set of Ornstein-Uhlenbeck processes. Each of the subenvironments of the composite environment is able to induce non-Markovian dynamics of the two-level system. By means of the numerically exact method of hierarchical equations of motion, we study steady states of the two-level system, evolution of the reduced density matrix and the equilibrium emission spectra in dependence on the frequency cutoffs and the coupling strengths of the subenvironments. Additionally, we investigate the impact of the rotating wave approximation (RWA) for the interaction with the bath on accuracy of the results.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time»

Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time

V.A. Mikhailov ', N.V. Troshkin 1 'Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34

Abstract

In this paper we investigate non-Markovian evolution of a two-level system (qubit) in a bosonic bath under influence of an external classical fluctuating environment. The interaction with the bath has the Lorentzian spectral density, and the fluctuating environment (stochastic field) is represented by a set of Ornstein-Uhlenbeck processes. Each of the subenvironments of the composite environment is able to induce non-Markovian dynamics of the two-level system. By means of the numerically exact method of hierarchical equations of motion, we study steady states of the two-level system, evolution of the reduced density matrix and the equilibrium emission spectra in dependence on the frequency cutoffs and the coupling strengths of the subenvironments. Additionally, we investigate the impact of the rotating wave approximation (RWA) for the interaction with the bath on accuracy of the results.

Keywords: non-Markovian evolution, bosonic reservoir, stochastic field, two-level system.

Citation: Mikhailov VA, Troshkin NV. Non-Markovian decoherence of a two-level system in a Lorentzian bosonic reservoir and a stochastic environment with finite correlation time. Computer Optics 2021; 45(3): 372-381. DOI: 10.18287/2412-6179-CO-776.

Introduction

Almost every quantum system interacts with its surroundings in a way that makes the system nearly impossible to isolate completely. The interaction gives rise to processes of decoherence and appears to be harmful in some circumstances, e.g. in quantum information processing [1 - 4], and to be a valuable resource in others [5 - 7]. In some cases, it is possible to capture evolution of the quantum system by means of a master equation of Lindblad form [8 - 10]. Systems of that type are called Markovian, and their evolution has the form of a quantum Markov process. In the Markovian approximation the environment is assumed memoryless owing to sufficiently fast relaxation processes that restore its equilibrium state almost instantly. Markovian systems are common in quantum optics, where a quantum system is often weakly coupled to an environment characterized by negligibly small correlation time.

Systems characterized by memoryful environments belong to the much wider category of non-Markovian systems [9, 11]. Due to increased experimental and computational capabilities, non-Markovian systems are of great interest today. Among them there are such well-known systems as quantum dots [12 - 15], micromechan-ical resonators [16, 17], superconducting qubits [18 - 20], and many others. Non-Markovian effects are ubiquitous in physics, chemistry, and biology and for systems interacting with either bosonic or fermionic reservoirs, like photosynthetic systems [21 - 25], molecular aggregates [26], molecular magnets [27], and solar cells [28]. Recently non-Markovian environments started to gain attention in quantum information processing [29 - 31], where attempts are made to utilize the backflow of information from the environment.

Accurate description of non-Markovian systems is a more complicated problem in comparison with description of Markovian systems. There are a plenty of methods known, but none of them are generally applicable, i.e. each one has its strong and weak sides, and the number of systems it can describe efficiently is often limited. Analytical methods are represented mostly by perturbative expansions for some special parameter regimes, e.g. effective weak coupling theories or the projection operator techniques [32, 33]. Numerically exact methods include the ones based on enlarging of the system state space, e.g. by extending the system space by the most relevant environmental modes [34 - 36], the ones utilizing tensor network approximations in propagation of influence func-tionals [37 - 39] and in mappings on effective 1D fermi-onic and bosonic chains [40 - 42], and etc. One of the most well-established numerically exact methods is the method of hierarchical equations of motion (HEOM) [43 - 45]. HEOMs utilize infinite systems of recurrent differential equations to encode the memory kernel of system-environment interaction and are able to handle a great variety of environmental spectral densities.

Switching from the Markovian approximation to a full non-Markovian description reveals many interesting phenomena. The most common one is the emergence of oscillations [46 - 50]. Non-Markovianity is known to be able to affect steady states (equilibrium states) of a system, e.g. it causes non-canonical steady states to appear [51], and also it affects the correlation functions. Actually, if non-Markovianity of a quantum process is sufficiently high, the quantum regression theorem (QRT) stops giving reliable correlation functions [52 - 54], which leads, for example, to significant differences between the predicted and the actual emission spectra. It is a

common practice to use the rotating wave approximation (RWA) for interaction with Markovian environments. Often the RWA is still used where the evolution is already non-Markovian, which can cause problems if non-Markovianity of the evolution is sufficiently high [11, 55, 56]. For example, wrongly used RWA may lead to incorrect shifts of the system frequencies [57] or may cutoff all non-Markovianity [56].

When an environment is composite, i.e. consists of several subenvironments, it is possible to utilize one of the subenvironments to control decoherence of a quantum system. The case of a stochastic subenvironment as a control tool is rather popular in literature [58 - 63] and has similarities with dynamical decoupling schemes that alter the environment spectral density via filtering functions realized in sequences of laser impulses [64].

In the paper we investigate non-Markovian evolution of a two-level system (TSL) in a bosonic bath under influence of an additional external fluctuating environment. The bath spectral density function is chosen to be Lorentzian [8, 65, 66]. The Lorentzian spectral density is suitable, for example, for interaction between a Jaynes-Cummings cell and a zero-temperature bosonic bath [67]. The stochastic environment is represented by a set of Ornstein-Uhlenbeck random processes. Following [68], we derive a HEOM capable of handling both RWA and non-RWA couplings with the bath equally accurate. Assuming the bath and the fluctuating environment to be independent, we analyze steady states of the TLS, evolution of the reduced density matrix, and equilibrium emission spectra. We investigate the dependence on frequency cutoffs and coupling strengths of spectral densities of the subenvironments in both RWA and non-RWA cases and provide a comparison with the ones obtained in the Mar-kovian approximation [69].

The paper is organized as follows. In Sec. 1 we introduce the model, in Sec. 2 we present the hierarchical equations of motion. Next, we study the TLS evolution numerically. In Sec. 3 we study steady states of the TLS, in Sec. 4 we investigate evolution of the reduced density matrix, and in Sec. 5 we investigate the emission spectra. Finally, we draw conclusions.

1. Model

The full Hamiltonian for the system can be written as

H = HA + £+ HIB + Hif , (1)

k=1

where Ha = h wo C+ or_ is the Hamiltonian for the free TLS, o>o is the TLS frequency, c+ and cc- are the rising and the lowering operators of the TLS, respectively; bk+ and bk form a set of creation and annihilation operators describing the modes of the bosonic bath; Hif is the Hamiltonian for the interaction of the TLS with the fluctuating environment (stochastic field) and Hib describes the interaction between the TLS and the bath.

The Hamiltonian for interaction with the stochastic field is defined by the next expression

HIF = hQ.(t )ct+CT _ + )CT + )CT _

(2)

where Q (t), £ (t) are random functions, the overbar means complex conjugation. The interaction gives rise to two decoherence channels, a dephasing channel and a relaxation channel, originating from the first and the second terms in (2), respectively.

The random function Q (t) is a real random process and £ (t) is a complex random process. We assume that Q (t) and £ (t) are Markov processes of Ornstein-Uhlenbeck (OU) type [70], and consider the real and imaginary parts of £ (t) as two independent real OU processes £1 (t) and £2 (t), respectively. Correlation functions of the random processes have the same form

(v(t )v(t ')) = ^ e

Yv

_Yv\t _/|

(3)

where v e {Q, £1, £2}, Av is the standard deviation of the OU process and defines the coupling strength with the stochastic field, 1 / Yv is the correlation time of the OU process and Yv has the physical meaning of the cutoff frequency of the environment.

The Hamiltonian Hib is used in two forms, the form corresponding to the full electric-dipole interaction (non-RWA) and the form representing the interaction in the rotating wave approximation (RWA). Introducing a new auxiliary TLS operator a, we can combine both forms of Hib in one expression

Hib = X (gkâbk + gkâ+bk ) ;

(4)

where gk are the TLS-bath coupling constants, â+ is the adjoint of â. If â = â+ + à-, (4) describes the full interaction, for â = à - it corresponds to the RWA interaction.

For path integral methods it is more naturally to define interaction with environment in the continuous form via the spectral density function, instead of utilizing the coupling constants gk directly. In the paper we consider the Lorentzian spectral density [8, 65 - 67]

r ( ) 1 ¿2ABYB Jl (ro) = —

n (ro_roo)2 + YB

(5)

where Ab defines the coupling strength, yb is the bath spectral width, also called the environment cutoff frequency. For the noise induced by the bath, 1 / yb represents the correlation time. The parameters have close relation to the corespondent parameters of the stochastic field Av and Yv (3) and generally have the same physical meaning in terms of impact on the TLS dynamics.

2. Hierarchical equations of motion

Let us suppose that before the initial moment of time the TLS does not interact with the bath and the stochastic field

k=1

KoMntrorepHaa omma, 2021, tom 45, №3 DOI: 10.18287/2412-6179-CO-776

373

and both of the subenvironments are at equilibrium. Then the total density matrix at t=to has the factorized form

Pot (Q f, f, to) = Peq (Q, Sb ^2)p('4)(to) »^(to) , (6)

where Peq (Q, §1, §2) is the factorizable Gaussian equilibrium distribution function of the stochastic field, p(A) (to) denotes the initial density matrix of the TLS, and p(B)eq(to) is the equilibrium bath density matrix at zero temperature.

For the factorized initial conditions (6) and the Lorentzian spectral density (5), we can obtain the HEOM by the steps presented in [68], where we replace the high-temperature Drude spectral density with the Lorentzian spectral density and take the limit в ^ <x> for the bath subenvironment. The HEOM expression has the same form and can be written as

д i --fpm4 ч t) = --и ap m4 ч t) +

д t n

+ X [ mk a [F> p m4>(t) + &Z P m4L( t) +

k e {field}

+ mkф£кPLAL (t)] + (7)

+ X [mka\B>pLA>(t) + Ф%PLAi (t) +

кe {bath }

+ m*<ФВ')кPLal (t)],

where m combines the field and the bath indexes into one composite index. We assume that the first three components of m index the recursion relations for the stochastic field in the order {Q, §1, §2}. The last two components index the bath recursion relations. The special notation m|k+i is used for the index m with the k-th component increased by 1

m = (m1, m2,...),

I ( 1 ) (8)

m I*+1 = (mi,m2,...,m* +1,...),

p(A)m(t) denotes the m-th auxiliary density matrix, Hxa is the commutator superoperator for the free TLS, HAp = HaP - PHa. The actual TLS density matrix starts the recursion and has index m = 0, p(A)(t) = p(A)o(t). The constants a(F)k and the operators Ф(\к and Ф(1)г,* belong to the stochastic field part of the recursion relations and can be expressed via parameters of the stochastic field and the field coupling operators

akF > =-Tvk,

ф% = -Avk (i / nVFvk, (9)

ф% =-Avk (i / A)VF,vk,

where v* denotes the k-th element in {Q, §1, §2}, e.g. V1 = Q, i^F.v* are commutator superoperators, i.e. Vхf.v* p = Vf.v* p - p Vf ,vk, the corresponding operators on the TLS subspace

Vr Q = Йст+< -,

Vf ,f1 = Й(а++ст -), (Ю)

Ф% = /й(а + -< -).

The remaining coefficients a(B)k, Ф(())г, k and Ф(1)г. k define the bath part of the recurrence relations and depend on the bath spectral density. If the spectral density is Lorentzian (5), we have the next expressions for the coefficients

akf =-( У в - i®o),

akf' =-(У в + ifflo),

<bB0i = bct ,

(11)

<B°l bc2, <Bk =-(AВ /2)cR, <B\ = (AВ / 2)c ,

where 61 = a and 62 = a+, the superscript "X" denotes the commutator superoperator defined earlier, and the superscript "R" denotes the superoperator for action from the right, i.e. 6R2 p = P62.

From (11) it follows that the bath part of the HEOM cannot be transformed into the one-indexed form, in contrast with the stochastic field part [44] and the case of non-RWA interaction with the high-temperature Drude bath [68], because there are two unequal constants a(B)k1 and a(B)k2. While the infinite-temperature Drude bath appears to be quite similar to the stochastic field in terms of the TLS dynamics [43], we expect that the zero-temperature Lorentzian bath and the stochastic field act on the TLS in qualitatively different ways.

The HEOM (7) is an infinite system of ordinary differential equations for elements of the main (m = 0) and the auxiliary (m Ф 0) density matrices over the TLS state space. The auxiliary density matrices improve accuracy and should be cut at some m, where the accuracy goal is satisfied. It is possible to truncate the recurrence relations in such a way that the zero-order approximation will coincide with the quantum master equation for relevant parameter regimes [43]. In practice a simpler way is usually sufficient, where all further auxiliary density matrices are considered zeros. The truncated system of ODEs can be solved by any numerical method for stiff systems of ODEs. In our case the best performance was achieved with explicit Runge-Kutta methods.

3. Steady states

The steady states of the TLS reachable from the selected initial state (6) can be obtained by propagating the state forward in time until the reduced density matrix stops changing. Typical response of the steady states on changes of the environment parameters is presented in fig. 1.

For a wide range of frequency cutoffs and coupling strengths tested, the stochastic field acting alone brings the TLS to the steady state where both the excited and the ground states are equally possible. Qualitatively different picture is observed for the TLS in the Lorentzian bath. Here, if one of the approximations is used, either the RWA or the Markovian, the TLS equilibrates in its ground state. Otherwise, in case of the non-RWA interac-

tion with the bath, the probability of the excited state is above zero and gradually rises with the frequency cutoff, starting near zero and then approaching 1 / 2 from below (fig. 1a). Also in the non-RWA case the stationary state exhibits a weak dependence on the coupling strength, the

Pw(°°) 0.5

Pw(°°) 0.50.40.30.20.1

excited state probability rises from zero to 1 / 2 from below, but much slower, then it is in the frequency cut-off case. Thus, the RWA, same as the Markovian approximation that intrinsically utilizes the RWA, alters the stationary states behavior drastically.

P«(°°) 0.5 i

Ys/coo

0.4 0.3 0.2 0.1

4

P/i(œ)

0.5

0.4

0.3- It

/

0.2- 1

0.1■

(b)

(c)

Yfi/W

1/2

0.4

0.8

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

1.2

Pw(°°) 0.50.40.30.20.1-

Ав/(о о

Ap/m'o'

(d) 1 2 3 4 5 (e) 0.4 0.8 1.2

Fig. 1. (a) The excited state population of the TLS in a steady state for interaction with the Lorentzian bath only (the field is off) as a function of the bath frequency cutoff yb , Ab /mo1/2 = 0.4. Solid curve denotes the non-RWA coupling with the bath. The RWA and the Markovian approximation curves coincide with the horizontal axis. (b, c) The excited state population of the TLS in a steady state for simultaneous interaction with the stochastic field and the Lorentzian bath as a function of (b) the bath frequency cutoff yb , Ab /mo1'2 = 0.4, and (c) the bath coupling strength Ab , yb /mo = 0.8. The dashed and the solid black curves denote the RWA and the non-RWA couplings, respectively, the gray curve stands for the Markovian approximation. The stochastic field is characterized by Yv /mo = 0.2 and Av /mo1'2 = 0.4. (d, e) The excited state population of the TLS in a steady state for simultaneous interaction with the stochastic field and the Lorentzian bath as a function of'(d) the field frequency cutoff Yv = yf , Av /mo1/2 = 0.4, and (e) the field coupling strength Av = Af, Yv /mo = 0.2. The dashed and the solid black curves denote the RWA and the non-RWA couplings, respectively, the gray curve stands for the Markovian approximation. The bath is characterized by yb /mo = 0.8, Ab /mo"2 = 0.4

When the TLS interacts with both subenvironments (fig. 1b - e), the difference between the stationary states for each of the subenvironments becomes visible. Because the difference may be big, e.g. for low bath frequency cutoffs in the non-RWA case or for any parameters in the RWA case, the impact of the stochastic field can be significant.

If we fix parameters of the stochastic field and start increasing the bath frequency cutoff yb (fig. 1b) or the coupling strength with the bath Ab (fig. 1c) starting from zero, the bath contribution in a steady state grows and the excited state population in the steady state decreases, because the steady states for decoherence in the bath always lie lower. Because in the RWA case the steady states are all completely unexcited (fig. 1a), the RWA curve is always below the non-RWA one. The distance between them constantly increases as the non-RWA steady state for decoherence in the bath shifts up with both the bath frequency cutoff yb and the bath coupling strength Ab. The RWA curves exhibit saturation in both figures, but the non-RWA curves reach minimums, go up and approach the Markovian curves from below (more pronounced in fig. 1b). In the Markovian approximation the steady states seem insensitive to any changes of the bath spectral density parameters.

Now let us fix the bath parameters and change the stochastic field instead. If we begin to gradually increase the stochastic field frequency cutoff from zero changing all the random processes simultaneously Yv = Yf (fig. 1e), the field contribution in the resulting steady state starts growing and the steady excited state population of the TLS starts growing either, because the steady states for decoherence in the stochastic field always lie higher. At some value of yf we reach the maximum, and after it the excited state population starts decreasing. The RWA applied to the interaction with the bath shifts the steady state down with respect to the one of the non-RWA curve, while the Markovian approximation gives much higher steady excited states populations for small frequency cutoffs and significantly overestimates the rate at which they decrease. The situation is similar if we change the field coupling strength in the same way Av = Af (fig. 1d), but there is no maximum, and the curves continues to rise, approaching the Markovian curve from below.

The observed behavior can be explained via the magnitude of the environment spectral density (for the stochastic field it is the spectrum of corresponding random processes) in the vicinity of the TLS resonant frequency. The OU random processes and the Lorentzian bath have

Компьютерная оптика, 2021, том 45, №3 DOI: 10.18287/2412-6179-CO-776

375

spectral densities with one peak. The OU process peak is located at œ = 0 and becomes wider and lower when the correspondent frequency cutoff increases. In fig. Id we see how the stochastic field contribution in the steady states grows at first, causing the rise of the steady states curves, because the spectral density peak widens and the magnitude of the spectral density near the resonant frequency increases, then the contribution falls, because the process of the spectral density peak declining becomes dominating, and the curves go down too. The switch from the growth to the decline of the field spectral density near the resonance frequency determines the maximum of the RWA and the non-RWA curves.

In contrast, the peak of the bath spectral density is always located at resonance and widens when the bath frequency cutoff increases. In fig. lb we see how the widening increases the magnitude of the bath spectral density in the vicinity of œ0 and the contribution of the bath in the steady states start growing, causing the curves to go down, because the TLS steady states for interaction only with the bath are located lower. Then the magnitude reaches its maximum and the contribution saturates. The subsequent rise of the non-RWA curves can be explained by the rise of the non-RWA curve for decoherence in the bath only (fig. 1a).

Coupling strengths of both subenvironments affect only heights of the corresponding spectral density peaks, when a coupling strength grows, the peak grows either. It results in gradual increase of the magnitude of a spectral density near the TLS resonant frequency and can be seen in fig. 1c, e. In fig. 1c the bath contribution increases, shifting the curves down, while in fig. 1e the field contribution increases, shifting the curves up. The situation becomes more complex if we stop using the restrictions Yv = Yf with Av = Af and allow arbitrary changes for each of the underlying random processes of the stochastic field. The detailed study of impact of each of the processes on the steady states will be presented elsewhere.

4. Density matrix evolution

An initially excited TLS placed in an equilibrium non-Markovian environment loses coherence due to the interaction with the environment, but the process is not monotone. At some point during the evolution the information backflow from the environment begins restoring the coherence, then the backflow weakens and the TLS starts losing coherence again. As a result, the oscillatory behavior emerges.

The evolution of the reduced density matrix from the factorized initial state (6) for different parameter regimes is presented in figs. 2, 3. The curves corresponding to the non-Markovian evolution, both for the RWA and the non-RWA types of coupling with the bath, exhibit rapidly vanishing oscillations, which are more evident for the RWA curves and for all the non-Markovian curves corresponding to decoherence in the stochastic field (fig. 2).

P«(0 1.0

0.90.8 0.7 0.60.50.4-

(a)

CO ot

Pi lit)

(O0t

Fig. 2. Evolution of the TLS excited state population in the stochastic field in dependence on (a) the field frequency cutoff YF and (b) the field coupling strength Af. Black denotes the non-Markovian curves, gray stands for the Markovian ones. In (a) Yv /mo = yf /mo = {0.2, 0.4, 0.8} and Av /mo1'2 = Af /mo1'2 = 1.6, for {dotted, dashed, solid} curves respectively, and in (b)

Yv /mo = yf /mo = 0.4, Av /mo1/2 = Af /mo1/2 = {0.4, 0.8,1.6}

The amplitude of the oscillations has a clear relation to the shapes of spectral densities of the subenvironments in the vicinity of the TLS resonance frequency rao. If the environment spectral density is flat in the vicinity of rao, which is the case of large frequency cutoffs, there are no oscillations. For example, the Markovian approximation curves in figs. 2a, 3a exhibit no oscillations, because the Markovian approximation assumes that environment correlation times are small, which corresponds to large frequency cutoffs. If the environment spectral density is not flat in the vicinity of the TLS resonance, the oscillations appear. The dependence of the oscillations amplitude on the frequency cutoff value in figs. 2a, 3a is more evident in case of decoherence in the bath, because the peak of its spectral density is located at the TLS resonance frequency. The coupling strength of the environment impacts the amplitude of the oscillations in the opposite way (figs. 2b, 3b).

The evolution becomes faster if the coupling strength increases, i.e. the minimums are located closer and the steady states are reached earlier (fig. 2b, 3). The frequency cutoff impacts the speed of the evolution in a more complex way, there is a cutoff frequency for which the evolution speed is maximum (not shown).

The main difference between the RWA and the non-RWA curves in fig. 3 resides in values of the minimums. The RWA curves have its minimums placed at the hori-

zontal axis where also the stationary values are located, while the minimums of all the non-RWA curves are placed strictly above the corresponding stationary values. The stochastic field curves for sufficiently large couplings have their first minimum located below the stationary value (fig. 2), which resembles the behavior of the non-RWA curves for decoherence in the bath. In overall, the non-RWA evolution of the reduced density matrix reminds the smoothed version of the RWA evolution. When the TLS interacts with both subenvironments simultaneously (fig. 4), the evolution becomes faster due to the presence of an additional decoherence channel, the first minimum is moved to the left and the stationary value is reached earlier. Because the steady state is shifted up by the stochastic field, the curves lie above the corresponding curves for decoherence in the bath only. Also the stochastic field significantly damps the oscillations, which is evident even for rather weak coupling strengths with the field in comparison with the coupling strength with the bath. The Markovian approximation shows the fastest decoherence among all the curves.

Pn{t)

the bath frequency cutoff yb and (b) the bath coupling strength Ab. The thin and the thick black curves (any stroke style) denote the RWA and the non-RWA couplings, respectively, the gray curves stand for the Markovian approximation.

In (a) yb /mo = {0.2, 0.4, 0.8}, Ab /mo1'2 = 1.6, for {dotted, dashed, solid} curves respectively, and in (b) yb /mo = 0.4mo, Ab /mo1'2 = {0.4, 0.8,1.6}

5. Emission spectrum

We obtain the equilibrium emission spectra of the TLS by applying the Fourier transform to the two-time correlation function <&+ (t2) a_ (ti)>, where t2 > ti. The time

ti is selected sufficiently big for the reduced density matrix evolution to reach its stationary phase, thereby the correlation function may be considered stationary. We calculate the stationary correlation function in the following way. First we propagate the initial state to the steady state, then apply the operator a_ to all density matrices p(A)m (ti), lying in the TLS subspace, next the result is propagated to t2, where &+ is applied.

P n(t)

simultaneous interaction with the stochastic field and the Lorentzian bath (solid curves) in comparison with the case of interaction with the bath only (dashed curves). The thin and the

thick black curves (any stroke style) denote the RWA and the non-RWA couplings, respectively, the gray curves stand for the Markovian approximation, Yv/mo = yf/mo = 0.4, Av /mo112 = Af /mo"2 = 0.4, and yb /mo = 0.4

The equilibrium emission spectra for interaction with the Lorentzian bath can be obtained only in the case of the non-RWA coupling, because the steady states in the Markovian and the RWA approximations are completely unexcited and do not emit (fig. 5). If the frequency cutoff is large, the non-RWA emission spectrum has one peak, which shifts to the right when the cutoff becomes smaller (not shown). At some cutoff value the top of the peak becomes a plateau and splits in two practically non-distinguishable peaks of non-equal intensity which are placed symmetrically with respect to ra = rao (not shown). Then the peaks move in the opposite directions, slightly declining, but stop at some value of the cutoff and start moving backwards while becoming more and more distinct (fig. 5a). At the same time the overall intensity of the spectrum gradually decreases. As a result, the two-peaked spectra are fairly weak in comparison with the one-peaked spectra. In fig. 5 we use the normalization by the maximum value, so the actual spectrum intensity is not shown. Potentially the appearance of the two peaks is explained by the presence of the two complex-conjugated coefficients a(B)ki and a(B)k2 in (ii) instead of one, e.g. for decoherence in the stochastic field. Also there is a zero-intensity point for ra = - rao and two peaks on both sides of it. The peaks move towards each other when the frequency cutoff decreases.

The dependence on the coupling strength with the bath is slightly more simple. If we increase it, the main peak widens and moves to the right, then splits in two peaks of unequal intensity (fig. 5b). If we increase the

KoMntroTepHaa onTma, 202i, tom 45, №3 DOI: i0.i8287/24i2-6i79-CO-776

377

coupling strength further, the right peak, which is the main peak, decreases, and the left (the side peak) grows, while moving in the opposite directions, the peaks resolution becomes better. For large coupling strengths the side peak becomes the dominant and approaches ra = 0.

1.0

(a) -3

(b) -3-2-1 0 1 2 (B/ra o

Fig. 5. Emission spectra of the TLS in the Lorentzian bath, normalized by maximum values, in dependence on (a) the bath frequency cutoff yb and (b) the bath coupling strength Ab.

In (a) yb /mo = {0.1, 0.2, 0.4}, Ab /mo1/2 = 0.6, for {dotted, dashed, solid} curves respectively, and in (b) yb /mo = 0.2, Ab /mo1/2 = {0.3, 0.6,1.2}

For comparison, we show the equilibrium emission spectra for decoherence in the stochastic field in fig. 6. The spectrum differs qualitatively from the case of deco-herence in the Lorentzian bath. The stochastic environment spectral density has a peak located at ra = 0 for any value of the frequency cutoff and the coupling strength. The peak becomes more distinct if the frequency cutoff lowers or the coupling strength rises. The resonance is clearly visible in fig. 6, where a side peak located at ra=0 appears if the cutoff frequency is sufficiently small or the coupling strength is sufficiently large. At the same time, the main peak shifts to the right and widens. The spectrum energy redistributes from the main peak to the side peak, and the side peak grows while the main peak decreases.

The Markovian approximation effectively considers the environmental spectral density flat (large frequency cutoffs), or, equivalently, it considers only a small region of the spectral density in the vicinity of the TLS resonance. If the spectral density resonance is located sufficiently far from the TLS resonance, the Markovian approximation loses the essential information about the peak existence. In fig. 6 it results in wide one-peaked

spectra, which peaks are located in the vicinity of the TLS resonance frequency.

/'.I \ \ i \ \

jV.Zl

V

(b) -3-2-10 1 2 3 CO/GOO

Fig. 6. Emission spectra of the TLS in the stochastic field, normalized by maximum values, in dependence on (a) the field frequency cutoff yf and (b) the field coupling strength Af. Black denotes the non-Markovian curves, gray stands for the Markovian ones. In (a) yf /mo = {0.1, 0.2, 0.4}, Af/mo1/2 = 0.6, for {dotted, dashed, solid} curves respectively, and in (b) yf /mo = 0.2, Af /mo1/2 = {0.4, 0.6, 0.8}

In fig. 7 we show the impact of the stochastic field on the equilibrium emission spectra for decoherence in the bath. The stochastic field causes the emergence of the zero-frequency peak, like it does in fig. 6, so the spectrum obtains the three-peaked form. Also it lowers the left peak of the doublet (the side peak in fig. 5), widens it and shifts the main peak (the rightmost peak) to the right. The field smoothes the negative frequencies spectrum, making the zero point disappear, and increases intensity of the RWA curve so that it can be observed. The Markovian approximation is clearly inaccurate in the parameter regions selected.

Conclusion

We have studied non-Markovian evolution of a TLS in a composite environment consisting of two subenvironments, a zero-temperature bosonic bath characterized by the Lo-rentzian spectral density and a stochastic field of the Ornstein-Uhlenbeck type, and analyzed the impact of the ro-tating-wave approximation for the interaction with the bath.

It was shown that the steady states for decoherence in the bath depend on the TLS-bath coupling type: the full interaction leads to different steady states in comparison with the cases when either the RWA or the Markovian

approximation is used. We investigated the joint influence of the subenvironments on the steady states and found connections with the shape of the environment spectral density in the vicinity of the TLS resonant frequency.

Fig. 7. Emission spectra of a TLS for simultaneous interaction with the stochastic field and the Lorentzian bath (solid curves) in comparison with the case of interaction with the bath only (dashed curves), normalized by maximum values. The thin and the thick black curves (any stroke style) denote the RWA and the non-RWA couplings, respectively, the gray curves stand for the Markovian approximation. Yv/mo = YF/mo = 0.2,

Av/mo1/2 = Af/mo1/2 = 0.2, and yb /mo = 0.2, Ab /mo1/2 = 0.6

We demonstrated how the reduced density matrix evolution depends on the frequency cutoff and the coupling strength of the environment. In all cases, except the Markovian approximation ones, the reduced density matrix exhibits the oscillatory behavior, the amplitude of which can be explained via shapes of spectral densities of the subenvironments in the vicinity of the TLS resonance frequency. We showed that an increase of the frequency cutoff smooths the oscillations and shifts the first minimum to the right while an increase of the coupling strength acts in the opposite way. Comparing the cases of the full and the RWA couplings to the bath, we found that minimums of the oscillations in the RWA can be located below the stationary value in contrast to the full interaction case where this is not observed.

We investigated the dependece of the TLS equilibrium emission spectra on the frequency cutoffs and the coupling strengths of the subenvironments and found that the spectra can have the doublet form in the positive frequencies domain and the doublet form in the negative frequencies domain at the same time, if the TLS interacts with the bath only, but the intensity of the spectrum is relatively low. For interaction with both subenvironments, the spectrum can have three distinct peaks, one of which is located at the zero frequency, and the other two are located at the opposite sides of the TLS resonance.

References

[1] Koch CP. Controlling open quantum systems: tools, achievements, and limitations. J Phys Condens Matter 20i6; 28(2 i): 2i300i. DOI: i0.i088/0953-8984/28/2i/2i300i.

[2] Khurana D, Agarwalla BK, Mahesh TS. Experimental emulation of quantum non-Markovian dynamics and co-

herence protection in the presence of information backflow. Phys Rev A 2019; 99: 022107. DOI: 10.1103/PhysRevA.99.022107.

[3] D'Arrigo A, Falci G, Paladino E. Quantum zeno and antizeno effect on a two-qubit gate by dynamical decoupling. Eur Phys J Spec Top 2019; 227(15): 2189-2194. DOI: 10.1140/epjst/ e2018-800081-0.

[4] Jing J, Wu L-A. Decoherence and control of a qubit in spin baths: an exact master equation study. Sci Rep 2018; 8(1): 1471. DOI: 10.1038/s41598-018-19977-9.

[5] Ban M. Decoherence of a two-qubit system interacting with initially correlated random telegraph noises. Quantum Inf Process 2020; 19(2): 46. DOI: 10.1007/s11128-019-2539-4.

[6] Moreira S, Marques B, Paiva R, Cruz L, Soares-Pinto D, Semiao F. Enhancing quantum transport efficiency by tuning non-Markovian dephasing. Phys Rev A 2020; 101(1): 012123. DOI: 10.1103/PhysRevA.101.012123.

[7] Maier C, Brydges T, Jurcevic P, Trautmann N, Hempel C, Lanyon B, Hauke P, Blatt R, Roos C. Environmentassisted quantum transport in a 10-qubit network. Phys Rev Lett 2019; 122(5): 050501. DOI: 10.1103/PhysRevLett.122.050501.

[8] Breuer H-P, Petruccione F, et al. The theory of open quantum systems. Oxford: Oxford University Press; 2002.

[9] Rivas A, Huelga SF, Plenio MB. Quantum non-markovianity: characterization, quantification and detection. Rep Progr Phys 2014; 77(9): 094001. DOI: 10.1088/0034-4885/77/9/094001.

[10] Lindblad G. On the generators of quantum dynamical semigroups. Commun Math Phys 1976; 48(2): 119-130. DOI: 10.1007/ BF01608499.

[11] de Vega I, Alonso D. Dynamics of non-Markovian open quantum systems. Rev Mod Phys 2017; 89: 015001. DOI: 10.1103/RevModPhys.89.015001.

[12] Wu J, Chen S, Seeds A, Liu H. Quantum dot optoelectronic devices: lasers, photodetectors and solar cells. J Phys D Appl Phys 2015; 48(36): 363001. DOI: 10.1088/00223727/48/36/363001.

[13] Meden V. The Anderson-Josephson quantum dot—a theory perspective. J Phys Cond Matter 2019; 31(16): 163001. DOI: 10. 1088/1361-648x/aafd6a.

[14] Tahara H, Ogawa Y, Minami F, Akahane K, Sasaki M. Longtime correlation in non-Markovian dephasing of an exciton-phonon system in inas quantum dots. Phys Rev Lett 2014; 112: 147404. DOI: 10.1103/PhysRevLett.112.147404.

[15] Bera D, Qian L, Tseng T-K, Holloway P. Quantum dots and their multimodal applications: A review. Materials 2010; 3(4): 2260-2345. DOI: 10.3390/ma3042260.

[16] Aspelmeyer M, Kippenberg T, Marquardt F. Cavity opto-mechanics. Rev Mod Phys 2014; 86(4): 1391-1452. DOI: 10.1103/ RevModPhys.86.1391.

[17] Gröblacher S, Trubarov A, Prigge N, Cole GD, Aspelmeyer M, Eisert J. Observation of non-Markovian microme-chanical brownian motion. Nat Commun 2015; 6(1): 7606. DOI: 10.1038/ncomms8606.

[18] Andersson G, Suri B, Guo L, Aref T, Delsing P. Non-exponential decay of a giant artificial atom. Nature Phys 2019; 15(11): 1123-1127. DOI: 10.1038/s41567-019-0605-6.

[19] Potocnik A, Bargerbos A, Schröder FAYN, Khan SA, Collo-do MC, Gasparinetti S, Salathe Y, Creatore C, Eichler C, Tü-reci HE, Chin AW, Wallraff A. Studying light-harvesting models with superconducting circuits. Nature Commun 2018; 9(1): 904. DOI: 10.1038/s41467-018-03312-x.

[20] Yu D, Dumke R. Open ising model perturbed by classical colored noise. Phys Rev A 2019; 100(2): 022124. DOI: 10.1103/PhysRevA.100.022124.

Компьютерная оптика, 2021, том 45, №3 DOI: 10.18287/2412-6179-CO-776

379

[21] Pfalzgraff W, Montoya-Castillo A, Kelly A, Markland T. Efficient construction of generalized master equation memory kernels for multi-state systems from nonadiabatic quantum-classical dynamics. J Chem Phys 20i9; i50(24): 244i09. DOI: i0.i063/i.50957i5.

[22] Hwang-Fu Y-H, Chen W, Cheng Y-C. A coherent modified redfield theory for excitation energy transfer in molecular aggregates. Chem Phys 20i5; 447: 46-53. DOI: i0.i0i6/j.chemphys.20i4.ii.026.

[23] Chin AW, Prior J, Rosenbach R, Caycedo-Soler F, Huelga SF, Plenio MB. The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes. Nat Phys 20i3; 9(2): ii3-ii8. DOI: i0.i038/nphys25i5.

[24] Lee MK, Huo P, Coker DF. Semiclassical path integral dynamics: Photosynthetic energy transfer with realistic environment interactions. Annu Rev Phys Chem 20i6; 67(i): 639-668. DOI: i0.ii46/annurev-physchem-0402i5-ii2252.

[25] Segal D, Agarwalla BK. Vibrational heat transport in molecular junctions. Annu Rev Phys Chem 20i6; 67(i): i85-209. DOI: i0.ii46/annurev-physchem-0402i5-ii2i03.

[26] Plenio MB, Almeida J, Huelga SF. Origin of long-lived oscillations in 2d-spectra of a quantum vibronic model: Electronic versus vibrational coherence. J Chem Phys 20i3; i39(23): 235i02. DOI: i0.i063/i.4846275.

[27] Coish W, Baugh J. Nuclear spins in nanostructures. Phys Status Solidi B Basic Res 2009; 246(i0): 2203-22i5. DOI: i0.i002/pssb. 200945229.

[28] Barford W. Electronic and optical properties of conjugated polymers. Oxford: Oxford University Press; 20i3.

[29] Latune CL, Sinayskiy I, Petruccione F. Quantum force estimation in arbitrary non-Markovian gaussian baths. Phys Rev A 20i6; 94: 052ii5. DOI: i0.ii03/PhysRevA.94.052ii5.

[30] Bylicka B, Chruscinski D, Maniscalco S. Non-Markovianity and reservoir memory of quantum channels: a quantum information theory perspective. Sci Rep 20i4; 4(i): 5720. DOI: i0.i038/srep05720.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

[31] Xiang G-Y, Hou Z-B, Li C-F, Guo G-C, Breuer H-P, Laine E-M, Piilo J. Entanglement distribution in optical fibers assisted by nonlocal memory effects, EPL 20i4; i07(5): 54006. DOI: i0.i209/0295-5075/i07/54006.

[32] McCutcheon DPS, Dattani NS, Gauger EM, Lovett BW, Nazir A. A general approach to quantum dynamics using a variational master equation: Application to phonon-damped rabi rotations in quantum dots. Phys Rev B 20ii; 84: 08i305. DOI: i0.ii03/PhysRevB.84.08i305.

[33] Jang S. Theory of coherent resonance energy transfer for coherent initial condition. J Chem Phys 2009; i3i(i6): i64i0i. DOI: i0.i063/i.3247899.

[34] Garraway BM. Nonperturbative decay of an atomic system in a cavity. Phys Rev A i997; 55: 2290-2303. DOI: i0.ii03/PhysRevA.55.2290.

[35] Mascherpa F, Smirne A, Somoza AD, Fernández-Acebal P, Donadi S, Tamascelli D, Huelga SF, Plenio MB. Optimized auxiliary oscillators for the simulation of general open quantum systems. Phys Rev A 2020; i0i: 052i08. DOI: i0.ii03/PhysRevA.i0i.052i08.

[36] Tamascelli D, Smirne A, Huelga SF, Plenio MB. Nonper-turbative treatment of non-Markovian dynamics of open quantum systems. Phys Rev Lett 20i8; i20: 030402. DOI: i0.ii03/PhysRevLett.i20.030402.

[37] Makri N, Makarov DE. Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory. J Chem Phys i995; i02(ii): 4600-46i0. DOI: i0.i063/i.469508.

[38] Makri N, Makarov DE. Tensor propagator for iterative quantum time evolution of reduced density matrices. II. Numerical methodology. J Chem Phys 1995; 102(11): 4611-4618. DOI: 10.1063/1. 469509.

[39] Strathearn A, Kirton P, Kilda D, Keeling J, Lovett BW. Efficient non-Markovian quantum dynamics using time-evolving matrix product operators. Nat Commun 2018; 9(1): 3322. DOI: 10.1038/ s41467-018-05617-3.

[40] Prior J, Chin AW, Huelga SF, Plenio MB. Efficient simulation of strong system-environment interactions. Phys Rev Lett 2010; 105: 050404. DOI: 10.1103/PhysRevLett.105.050404.

[41] Tamascelli D, Smirne A, Lim J, Huelga SF, Plenio MB. Efficient simulation of finite-temperature open quantum systems. Phys Rev Lett 2019; 123: 090402. DOI: 10.1103/PhysRevLett.123.090402.

[42] Nüßeler A, Dhand I, Huelga SF, Plenio MB. Efficient simulation of open quantum systems coupled to a fermionic bath. Phys Rev B 2020; 101: 155134. DOI: 10.1103/PhysRevB.101.155134.

[43] Tanimura Y. Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissi-pative systems. J Phys Soc Japan 2006; 75(8): 082001. DOI: 10.1143/JPSJ.75.082001.

[44] Tanimura Y, Kubo R. Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath. J Phys Soc Japan 1989; 58(1): 101-114. DOI: 10.1143/JPSJ.58.101.

[45] Tanimura Y. Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and ther-modynamic quantities. J Chem Phys 2014; 141(4): 044114. DOI: 10.1063/1. 4890441.

[46] Semin V, Sinayskiy I, Petruccione F. Arbitrary spin in a spin bath: Exact dynamics and approximation techniques, Phys Rev A 2014; 89: 012107. DOI: 10.1103/PhysRevA.89.012107.

[47] Rossi MAC, Paris MGA. Non-Markovian dynamics of single- and two-qubit systems interacting with Gaussian and non-Gaussian fluctuating transverse environments. J Chem Phys 2016; 144(2): 024113. DOI: 10.1063/1.4939733.

[48] Mwalaba M, Sinayskiy I, Petruccione F. Dynamics and thermalization in a simple mesoscopic fermionic bath, Phys Rev A 2019; 99: 052102. DOI: 10.1103/PhysRevA.99.052102.

[49] Pavelev A, Semin V. Investigation of non-Markovian dynamics of two dipole-dipole interacting Qubits based on numerical solution of the non-linear stochastic Schrödinger equation. Computer Optics 2019; 43(2): 168-173. DOI: 10.18287/2412-6179-2019-43-2-168-173.

[50] Vasilev D, Semin V. Qubit dynamics in extern laser field. Computer Optics 2019; 43(4): 562-566. DOI: 10.18287/ 2412-6179-2019-43-4-562-566.

[51] Iles-Smith J, Lambert N, Nazir A. Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems. Phys Rev A 2014; 90: 032114. DOI: 10.1103/PhysRevA.90.032114.

[52] De Santis D, Johansson M, Bylicka B, Bernardes NK, Acin A. Correlation measure detecting almost all non-markovian evolutions. Phys Rev A 2019; 99: 012303. DOI: 10.1103/PhysRevA.99.012303.

[53] Ali MM, Lo P-Y, Tu MW-Y, Zhang W-M. Non-Markovianity measure using two-time correlation functions. Phys Rev A 2015; 92: 062306. DOI: 10.1103/PhysRevA.92.062306.

[54] McCutcheon DPS. Optical signatures of non-Markovian behavior in open quantum systems. Phys Rev A 2016; 93: 022119. DOI: 10.1103/ PhysRevA.93.022119.

[55] Fleming C, Cummings NI, Anastopoulos C, Hu BL. The rotating-wave approximation: consistency and applicability from an open quantum system analysis. J Phys A Math Theor 2010; 43(40): 405304. DOI: 10.1088/17518113/43/40/405304.

[56] Makela H, Mottonen M. Effects of the rotating-wave and secular approximations on non-Markovianity. Phys Rev A 2013; 88: 052111. DOI: 10.1103/PhysRevA.88.052111.

[57] Eastham PR, Kirton P, Cammack HM, Lovett BW, Keeling J. Bath-induced coherence and the secular approximation. Phys Rev A 2016; 94: 012110. DOI: 10.1103/PhysRevA.94.012110.

[58] Jing J, Yu T, Lam C-H, You JQ, Wu L-A. Control relaxation via dephasing: A quantum-state-diffusion study. Phys Rev A 2018; 97: 012104. DOI: 10.1103/PhysRevA.97.012104.

[59] Jing J, Li R, You JQ, Yu T. Nonperturbative stochastic dynamics driven by strongly correlated colored noise. Phys Rev A 2015; 91: 022109. DOI: 10.1103/PhysRevA.91.022109.

[60] Jing J, Wu L-A. Control of decoherence with no control, Sci Rep 2013; 3(1): 2746. DOI: 10.1038/srep02746.

[61] Brian Walton D, Visscher K. Noise suppression and spectral decomposition for state-dependent noise in the presence of a stationary fluctuating input. Phys Rev E 2004; 69: 051110. DOI: 10.1103/PhysRevE.69.051110.

[62] Wang ZH, Ji YJ, Li Y, Zhou DL. Dissipation and decoherence induced by collective dephasing in a coupled-qubit

system with a common bath. Phys Rev A 2015; 91: 013838. DOI: 10.1103/PhysRevA.91.013838.

[63] Semin V. Non-Markovian relaxation of a three-level atom in two laser fields with noise. Laser Phys 2020; 30(2): 025204. DOI: 10.1088/1555-6611/ab65c3.

[64] Biercuk MJ, Doherty AC, Uys H. Dynamical decoupling sequence construction as a filter-design problem. J Phys BAt Mol Opt 2011; 44(15): 154002. DOI: 10.1088/09534075/44/15/154002.

[65] Wu W, Lin H-Q. Quantum zeno and anti-zeno effects in quantum dissipative systems. Phys Rev A 2017; 95: 042132. DOI: 10.1103/PhysRevA.95.042132.

[66] Wu W. Realization of hierarchical equations of motion from stochastic perspectives. Phys Rev A 2018; 98: 012110. DOI: 10.1103/PhysRevA.98.012110.

[67] Li J-G, Zou J, Shao B. Non-Markovianity of the damped jaynescummings model with detuning. Phys Rev A 2010; 81: 062124. DOI: 10.1103/PhysRevA.81.062124.

[68] Mikhailov VA, Troshkin NV. Non-Markovian dynamics of a two-level system in a bosonic bath and a gaussian fluctuating environment with finite correlation time. Phys Rev A 2021; 103: 012208. DOI: 10.1103/PhysRevA.103.012208.

[69] Mikhailov VA, Troshkin NV. Master equation averaged over stochastic process realizations for the description of a three-level atom relaxation. Computer Optics 2016; 40(5): 649-653. DOI: 10.18287/2412-6179-2016-40-5-649-653.

[70] Risken H, Frank T. The Fokker-Planck equation. Berlin, Heidelberg: Springer-Verlag; 1996. DOI: 10.1007/978-3642-61544-3.

Authors' information

Victor Alexandrovich Mikhailov (b. 1952) graduated from Kuibyshev Aviation Institute in 1975 (KuAI; presently Samara National Research University named after academician S.P. Korolyov), majoring in Aircraft Engines, graduated from Kuibyshev State University, the faculty of physics, in 1981, majoring in Theoretical Physics, Candidate of Physical and Mathematical Sciences, associate professor of the Physics department at Samara National Research University. Research interests are in generalized coherent states, dynamics and relaxation of quantum systems, Fokker-Planck equations. E-mail: va_mikhailov@mail.ru .

Nikolay Vyacheslavovich Troshkin (b. 1989) graduated from Samara State Aerospace University in 2012 (SSAU; presently Samara National Research University named after Academician S.P. Korolev), majoring in Applied Physics, postgraduate student at Samara National Research University. Research interests are in quantum many-body physics, low-rank approximations, algorithms. E-mail: nick.troshkin@smail.com .

Received July 3, 2016. The final version - December 7, 2020.

KoMntrorepHaa omma, 2021, tom 45, №3 DOI: 10.18287/2412-6179-CO-776

381

i Надоели баннеры? Вы всегда можете отключить рекламу.