*Article* **The Entanglement Generation in** *PT* **-Symmetric Optical Quadrimer System**

#### **Joanna K. Kalaga 1,2**


Received: 12 July 2019; Accepted: 29 August 2019; Published: 3 September 2019

**Abstract:** We discuss a model consisting of four single-mode cavities with gain and loss energy in the first and last modes. The cavities are coupled to each other by linear interaction and form a chain. Such a system is described by a non-Hermitian Hamiltonian which, under some conditions, becomes PT -symmetric. We identify the phase-transition point and study the possibility of generation bipartite entanglement (entanglement between all pairs of cavities) in the system.

**Keywords:** PT -symmetry; entanglement; negativity

#### **1. Introduction**

One of the main principles of quantum mechanics is the assumption that the Hamiltonian describing a quantum system must be Hermitian. In consequence, all of Hamiltonian's eigenvalues are real. In 1998, Bender and Boettcher [1] showed that the Hermiticity of a Hamiltonian (*H*ˆ = *H*ˆ †) is not the only condition to obtain its real eigenvalues. When non-Hermitian Hamiltonian has PT symmetry, its eigenvalues are real. PT symmetry means that the Hamiltonian satisfies commutation relations *H*ˆ , *P*ˆ*T*ˆ = 0, where *P*ˆ is a linear parity operator which changes the sign of the momentum operator and the position operator; whereas *T*ˆ is the antilinear time reversal operator.

The first of the studied PT -symmetric Hamiltonians with a real spectrum was that described by the equation *<sup>H</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>p</sup>*ˆ<sup>2</sup> <sup>+</sup> *ix*ˆ<sup>3</sup> [1–3]. The most commonly discussed quantum mechanical systems with PT symmetry are those described by Hamiltonians including a complex potential. When the imaginary part has a plus sign, the system obtains energy from the environment; whereas the minus sign means that the system gives energy into the environment. When the system satisfies PT symmetry, the loss and gain must be balanced. For optical PT -symmetric systems, the refractive index *n* can play a role of the potential energy *V* [4]. Quantum optical systems should also exhibit such balance between the gain and loss of energy to satisfy PT symmetry condition. For such a case, the refractive index obeys the symmetry relation *n*(*x*) = *n*∗(−*x*).

In recent years, the different kinds of PT -symmetric systems have been considered. For instance, there has been the pair of coupled *LC* circuits [5], optical lattice [6,7], optomechanical system [8], optical waveguides [9], quantum-dot [10], and others. The PT -symmetric systems can exhibit numerous interesting phenomena such as invisibility [11], chaos induced by the PT symmetry breaking [12], and many others.

For the PT -symmetric system, we can observe a phase-transition point which is the point where the system loses its PT -symmetric properties. If the system is in the PT -symmetric phase, all eigenvalues of Hamiltonian are real. When the system is in PT symmetry broken phase, it has complex eigenvalue spectra. Such a transition point from the unbroken PT -symmetric phase to the PT symmetry broken phase is called exceptional point [13–15]. In general, this singular point

occurs when eigenvalues and corresponding eigenvectors of the system depend on some parameters. When those parameters reach a critical value, then eigenvalues of the system coalesce and the spectrum becomes complex. If two eigenvalues coalescence, we have second-order exceptional points. During recent years, various type of singularities and their features have been studied [16–22]. The PT -symmetric systems exhibit many interesting phenomena related to the presence of an exceptional point. For instance, enhancing spontaneous emission [23], unidirectional invisibility in fiber networks [24], and loss-induced lasing [25].

In this paper, we concentrate on the entangled states generation in the PT -symmetric quadrimer system. In particular, we are interested in producing bipartite entangled states, and the influence of gain and loss rate on producing such states. The generation of entanglement in various quantum systems is one of the fundamental areas of interest in quantum information theory. Entanglement can be observed in various physical systems such as Bose–Einstein condensates [26,27], cavity QED [28], quantum dots [29,30], trapped ions [31], and many others [32–37]. The research related to the production of entanglement in open systems is of particular importance. For such a system, a crucial role is played by the decoherence processes, which are the consequence of the interaction of the system with the environment. One should realize that the entanglement is very sensitive to the noise processes. The interaction with the environment leads to losses of coherence, and, in consequence, destroys entanglement. Very interesting is sudden death and the rebirth of entanglement observed in various systems interacting with the environment [38–41].

The paper is organized as follows. In Section 2, we describe the PT -symmetric system. In particular, we derive the formulas determining the location of the exceptional point, the point where the system loses its PT symmetry, then, the eigenvalues of Hamiltonian become complex. For the different values of gain/loss rate, we analyze the possibility of generation entanglement state in Section 3. We show how the gain/loss rate influences the entanglement production process.

#### **2. The Model**

The considered system is composed of four identical single-mode cavities of the resonance frequency *ω*. The first cavity is passive (it loses the energy), the second and third ones are neutral (no losses and no gain of the energy), and the last one is active (it gains the energy). All cavities are coupled mutually by linear interaction and form a chain (see Figure 1). Such a system can be described by the following Hamiltonian:

$$\begin{aligned} \hat{H} &= \left(\omega - i\gamma\right) \mathfrak{d}\_1^\dagger \mathfrak{d}\_1 + \omega \mathfrak{d}\_2^\dagger \mathfrak{d}\_2 + \omega \mathfrak{d}\_3^\dagger \mathfrak{d}\_3 + \left(\omega + i\gamma\right) \mathfrak{d}\_4^\dagger \mathfrak{d}\_4 \\ &+ \beta \left(\mathfrak{d}\_1^\dagger \mathfrak{d}\_2 + \mathfrak{d}\_2^\dagger \mathfrak{d}\_1 + \mathfrak{d}\_2^\dagger \mathfrak{d}\_3 + \mathfrak{d}\_3^\dagger \mathfrak{d}\_2 + \mathfrak{d}\_3^\dagger \mathfrak{d}\_4 + \mathfrak{d}\_4^\dagger \mathfrak{d}\_3\right), \end{aligned} \tag{1}$$

where *a*ˆ*<sup>i</sup>* and *a*ˆ† *<sup>i</sup>* are bosonic annihilation and creation operators, respectively, whereas the indices 1, 2, 3, 4 label four subsystems. The parameter *β* describes the strength of linear interaction between two nearest cavities and *γ* denotes the strength of decay or the gain of cavities. We assume here that *h*¯ = 1 and the parameters *ω*, *γ*, and *β* are real.

**Figure 1.** Scheme of the model.

The Hamiltonian defined in Equation (1) is non-Hermitian, but for some situations it becomes PT -symmetric and its eigenvalues are real. To find the PT phase transition point, we write the Hamiltonian *H*ˆ in this form:

$$
\hat{H} = \left[ \begin{array}{cccc} \hat{a}\_1^\dagger & \hat{a}\_2^\dagger & \hat{a}\_3^\dagger & \hat{a}\_4^\dagger \end{array} \right] \mathcal{H} \begin{bmatrix} \hat{a}\_1 \\ \hat{a}\_2 \\ \hat{a}\_3 \\ \hat{a}\_4 \end{bmatrix} \,\tag{2}
$$

where

$$
\mathcal{H} = \begin{bmatrix}
\omega - i\gamma & \beta & 0 & 0 \\
\beta & \omega & \beta & 0 \\
0 & \beta & \omega & \beta \\
0 & 0 & \beta & \omega + i\gamma
\end{bmatrix}.
\tag{3}
$$

Next, we can calculate the eigenvalues of the Hamiltonian H, and they are

$$\begin{array}{rcl} E\_1 &=& \frac{1}{2} \left( 2\omega - \sqrt{2} \sqrt{3\beta^2 - \sqrt{5\beta^4 - 2\beta^2 \gamma^2 + \gamma^4} - \gamma^2} \right), \\ E\_2 &=& \frac{1}{2} \left( 2\omega + \sqrt{2} \sqrt{3\beta^2 - \sqrt{5\beta^4 - 2\beta^2 \gamma^2 + \gamma^4} - \gamma^2} \right), \\ E\_3 &=& \frac{1}{2} \left( 2\omega - \sqrt{2} \sqrt{3\beta^2 + \sqrt{5\beta^4 - 2\beta^2 \gamma^2 + \gamma^4} - \gamma^2} \right), \\ E\_4 &=& \frac{1}{2} \left( 2\omega + \sqrt{2} \sqrt{3\beta^2 + \sqrt{5\beta^4 - 2\beta^2 \gamma^2 + \gamma^4} - \gamma^2} \right). \end{array} \tag{4}$$

We see that all eigenvalues depend on the frequency *ω*, the strength of coupling *β*, and the gain and loss rate *γ*. The eigenvalues are real when two relations are satisfied simultaneously: <sup>3</sup>*β*<sup>2</sup> <sup>−</sup> 5*β*<sup>4</sup> <sup>−</sup> <sup>2</sup>*β*2*γ*<sup>2</sup> <sup>+</sup> *<sup>γ</sup>*<sup>4</sup> <sup>−</sup> *<sup>γ</sup>*<sup>2</sup> <sup>≥</sup> 0 and <sup>3</sup>*β*<sup>2</sup> <sup>+</sup> 5*β*<sup>4</sup> <sup>−</sup> <sup>2</sup>*β*2*γ*<sup>2</sup> <sup>+</sup> *<sup>γ</sup>*<sup>4</sup> <sup>−</sup> *<sup>γ</sup>*<sup>2</sup> ≥ 0.

In Figure 2a,b the real and imaginary parts of *Ei* are presented, respectively. We see there that the phase-transition point is localized when *γ* = *β*. If the gain/loss rate *γ* is smaller than the coupling parameter *β*, the spectrum is real and the system is in the PT -symmetric phase. As the parameter *γ* exceeds *β*, the system passes into the broken PT -symmetric phase, and the eigenvalues become complex. Observed here, the phase-transition point is the second-order exceptional point at which the two eigenvalues coalesce.

**Figure 2.** The values of (**a**) the real part of *Ei* and (**b**) the imaginary part of *Ei*, as a function of ratio *γ*/*β*. All parameters are scaled in *ω* units and *β* = *ω*/2.

To study the system's dynamics, we analyze the time evolution of the density matrix *ρ*ˆ. We can note that PT -symmetric system is an open system which exchanges energy with the environment. The energy is gained and lost by the system. The evolution of such a system is determined by a master equation for the operator *ρ*ˆ

$$\frac{d}{dt}\mathfrak{\dot{\rho}} = -i\left[\hat{H}\_{0\prime}\mathfrak{\dot{\rho}}\right] + \mathcal{L}\mathfrak{\dot{\rho}},\tag{5}$$

where

$$
\hat{H}\_0 = \omega \left( \hat{a}\_1^\dagger \hat{a}\_1 + \hat{a}\_2^\dagger \hat{a}\_2 + \hat{a}\_3^\dagger \hat{a}\_3 + \hat{a}\_4^\dagger \hat{a}\_4 \right) + \beta \left( \hat{a}\_1^\dagger \hat{a}\_2 + \hat{a}\_2^\dagger \hat{a}\_1 + \hat{a}\_2^\dagger \hat{a}\_3 + \hat{a}\_3^\dagger \hat{a}\_2 + \hat{a}\_3^\dagger \hat{a}\_4 + \hat{a}\_4^\dagger \hat{a}\_3 \right), \tag{6}
$$

and L is the Liouvillian superoperator, which in our case it takes the form:

$$\mathcal{L}\dot{\rho} = \gamma \left( 2\dot{a}\_1 \rho \dot{a}\_1^\dagger - \dot{a}\_1^\dagger \dot{a}\_1 \dot{\rho} - \rho \dot{a}\_1^\dagger \dot{a}\_1 \right) + \gamma \left( 2\dot{a}\_4^\dagger \rho \dot{a}\_4 - \dot{a}\_4 \dot{a}\_4^\dagger \dot{\rho} - \rho \dot{a}\_4 \dot{a}\_4^\dagger \right). \tag{7}$$

In Equation (7), the first term describes the loss and the second term is related to the gain of energy.

In further considerations, we assume that initially only one of the four cavities is in the one-photon state |1, and the other three are in vacuum state |0. Therefore, we will discuss here four cases:


where |*ijkl* = |*i*<sup>1</sup> ⊗ |*j*<sup>2</sup> ⊗ |*k*<sup>3</sup> ⊗ |*l*<sup>4</sup> are the four-mode states. For the first of them, only the passive cavity is in the one-photon state at the initial time, then, the density matrix describing the system can be written as |10001000|. The next two cases analyzed correspond to the situation when one of the neutral cavities is in the one-photon state. Finally, for the last situation discussed here, the system starts evolution from the state |00010001|. For such a case, only the active cavity is in the one-photon state at the initial time.

#### **3. The Entanglement Generation**

In further considerations, we discuss the generation of two-mode entangled states. For the analysis of the degree of bipartite entanglement between two cavities, we apply the entanglement measure which is based on the positive partial transposition criterion, the *bipartite negativity*. It was defined in [42,43] as

$$N\_{i\bar{j}}(\rho\_{i\bar{j}}) = \frac{1}{2} \sum\_{i} |\lambda\_{\mathbb{R}}| - \lambda\_{\mathbb{H}\_{\prime}} \tag{8}$$

where *ρij* = *Trk*,*<sup>l</sup> <sup>ρ</sup>ijkl* is the two-mode density matrix, *λ<sup>n</sup>* is *n*-th eigenvalue of the matrix *ρTi ij* , and *<sup>ρ</sup>Ti ij* describes the partial transposition (with respect to the *i*-th mode) of the two-mode density matrix *ρij*.

In our consideration, we analyze only maximal values of the bipartite negativity *N*<sup>0110</sup> defined in the space of four two-mode states: |0|0; |0|1; |1|0; and |1|1. Therefore, to quantify the bipartite entanglement, we chose the negativity because this quantity can clearly differentiate between entangled and unentangled states when it is applied to two-qubit or qubit—qutrit systems. The negativity takes values between 0 for separable states and 1 for maximally entangled ones. To find values of *N*0110, we generate the time evolution of the density matrix for the whole system, next, we find the reduced density matrix *ρij* and calculate negativities for the subsystems spanned in the four states. We note that for the system described by PT -symmetric Hamiltonian, its evolution is nonunitary. Therefore, we should normalize the density matrix during the all evolution-time, and then, calculate the negativity with the application of such normalized density matrix *ρijkl*.

As it was mentioned before, we study the possibility of generation of the two-mode entangled state for four cases related to the four initial states (see Figure 3).

First, we consider the situation when the first cavity (passive cavity) is excited, and the initial state is *ρ*ˆ(*t* = 0) = |10001000|. In Figure 3a, we show the dependence of the maximal values of negativities *N*<sup>0110</sup> on *γ*/*β*. One can see that only entanglement between modes 1-2 can be produced for all situations considered here. What is relevant is that the degree of entanglement appearing in the system strongly depends on the value of the gain/loss parameter. With the increased value of parameter *γ*, we can observe decreasing entanglement in all pairs modes. For example, when *γ* = 0.99*β*, the maximal value of *N*<sup>12</sup> becomes closed to 0.14. For small values of *γ*/*β*, the entanglement between all pairs of cavities is generated. Whereas for large values of *γ*/*β*, the entanglement corresponding to pairs of cavities 1-3, 1-4, 2-3, 2-4, 3-4 is not produced.

Next, we check the system's ability to produce entanglement when the evolution of our system starts from the state *ρ*ˆ(*t* = 0) = |01000100|. For such a case, the excitation initially appears in the second cavity (the cavity without loss and gain). In Figure 3b, we see that analogously as in the previous case, the entanglement between all cavities is generated only for small values of *γ*/*β* and the strength of maximal possible bipartite entanglement depends on the gain/loss rate again. For weak losses/gains (*γ* < 0.03*β*), the entanglement between the cavities 1-2, 2-3, 1-3 is significant, but the strongest entanglement can be observed between subsystems 2 and 3. For strong losses/gain, the entanglement between the cavities 1-2, 2-3 plays the main role, whereas the entanglement between modes 3-4, 2-4, 1-4 is not produced.

In the next analyzed case, for *t* = 0, the third cavity is in the one-photon state and the evolution of the system starts from the state *ρ*ˆ(*t* = 0) = |00100010|. Figure 3c shows that with increasing parameter *γ*, all bipartite entanglements become weaker and entanglement between subsystems 1-4 and 2-4 disappears. For small values of *γ*/*β*, when the strength of the loss/gain rate increases, the maximal values of all negativities significantly decrease. Whereas, for greater values of *γ*, the values of negativities describing entanglement between modes 1-2, 2-3, 1-3 practically remain constant. For *γ* > 0.2*β*, the strongest entanglement we observe is between two neutral cavities.

In the last case, the system's evolution starts from state *ρ*ˆ(*t* = 0) = |00010001|, and the excitation initially appears in the active cavity. In Figure 3d, we see that analogously as in the previous case, with increasing *γ* the maximal values of all negativities decrease. What is interesting is that in such a case, all negativities reach nonzero values. The entanglement between cavities 3 and 4 plays the main role. For subsystems 1-2, 2-3, 1-3, 2-4, the entanglement is weaker, and the weakest correlation we observe is between passive and active cavities 1-4. For small values of *γ*, the negativity *N*<sup>14</sup> significantly decreases; and for *γ* > 0.2*β* reaches values close to zero (*N*<sup>14</sup> < 0.01).

**Figure 3.** *Cont.*

**Figure 3.** Maximal values of the negativity versus *γ*/*β* for *ω* = 2*β* and for various initial states: (**a**) |10001000|, (**b**) |01000100|, (**c**) |00100010|, (**d**) |00010001|.

#### **4. Conclusions**

Obtaining an entanglement quantum state plays a crucial role in the quantum information theory. The entangled states have many possible applications, such as teleportation and dense coding. Quantum correlations, including entanglement, play a significant role in the development of quantum computation and quantum information processing. Therefore, it is important to understand the entanglement nature of quantum systems. In recent years, a lot of papers have been presented on characterizing entanglement in bipartite systems. The entanglement in such systems can be analyzed relatively easily. The quantum correlations in systems consisting of more than two subsystems have not been understood well, so they are extensively studied. Therefore, we were interested in studying the ability of a PT -symmetric system to generate the entangled states, which are especially interesting

from the point of view of quantum information theory. The purpose of this article was to contribute to studies by addressing the entanglement generation.

In this paper, the optical quadrimer system containing four cavities characterized by frequency *ω* was discussed. The cavities were coupled with each other in such a way that the system formed a chain. Additionally, the first and the last cavity were a subject of losses and gains of energy, respectively, and the rate of losses and gain was balanced.

For such a system, we have found a phase-transition point, and we have shown that when the gain/loss rate is equal or smaller than coupling parameter, all eigenvalues of the system's Hamiltonian are real. In other words, the system is PT -symmetric, contrary to the situation when the gain/loss rate reaches values greater than that of coupling strength and the optical quadrimer is in the broken PT -symmetric phase. For such a case, the system has complex eigenvalue spectra.

We have analyzed four situations corresponding to the four various initial states of the system. We were interested in the possibility of generation of entangled states and the influence of the gain/loss rate strength on the effectiveness of production of such states. We have shown that the degree of entanglement strongly depends on the value of the parameter *γ*. We have estimated the range of such values for which the strongest bipartite entanglement can be created. For all analyzed cases, when the values of loss/gain rate are small, with increasing *γ*, the entanglements significantly decrease. For the greater values of loss/gain rate, if the initially active or passive cavity is excited, the strongest entanglement appears between the cavity which, for *t* = 0, is in one-photon state and its nearest neighbor. For the initial time, when one of the neutral cavities is excited, we can observe the strongest entanglement between subsystems 2 and 3.

**Funding:** The results described in the present work were achieved with the financial support of the ERDF/ESF project "Nanotechnologies for Future" (CZ.02.1.01/0.0/0.0/16\_019/0000754) and the program of the Polish Minister of Science and Higher Education under the name "Regional Initiative of Excellence" in 2019–2022, project no. 003/RID/2018/19, funding amount 11 936 596.10 PLN.

**Conflicts of Interest:** The author declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


c 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Single-Qubit Driving Fields and Mathieu Functions**

#### **Marco Enríquez 1,\*, Alfonso Jaimes-Nájera 2,3 and Francisco Delgado <sup>1</sup>**


Received: 13 July 2019; Accepted: 7 September 2019; Published: 16 September 2019

**Abstract:** We report a new family of time-dependent single-qubit radiation fields for which the correspondent evolution operator can be disentangled in an exact way via the Wei–Norman formalism. Such fields are characterized in terms of the Mathieu functions. We show that the regions of stability of the Mathieu functions determine the nature of the driving fields: For parameters in the stable region, the fields are oscillating, being able to be periodic under certain conditions. Whereas, for parameters in the instability region, the fields are pulse-like. In addition, in the stability region, this family admits solutions for evolution loops in quantum control. We obtain some prescriptions to reach such a control effect. Geometric phases in the evolution are also analyzed and discussed.

**Keywords:** quantum control; Mathieu functions; time-dependent driving fields

#### **1. Introduction**

Exactly solvable models are valuable in quantum mechanics even though the number of cases is very limited. In the case of two-level systems (or *qubits*), some solutions are well-known and longstanding [1–3]. However, the interest in the control of these systems has recently increased due to its potential applications in quantum computation and quantum information [4]. Since the goal is to achieve precise control operations, one should deal with time-dependent driving fields in general. In this context, it is worth mentioning that the design of high-fidelity control operations is required [5,6] and some techniques, such as the adiabatic rapid passage, have been proposed [7]. Nevertheless, the analytical solution for an arbitrary driving field is still an open problem.

The dynamics of a system are determined by the correspondent dynamical law, the interactions to which it is subjected, as well as the initial conditions. However, an alternative approach to obtain exact solutions to such a dynamical law is constituted by the so-called inverse techniques in which some aspects of the dynamics are prescribed and then the interactions are found. This scheme has been widely used to deal with the control of two-level systems through time-dependent radiation fields [8–11]. Recently, a method to generate new solutions to the one qubit dynamical problem has been proposed within this framework. By requiring that the time-evolution operator be exactly factorized as a product of independent exponential factors involving only one Lie algebra generator, new families of time-dependent driving fields are obtained [12]. The evolution-operator disentangling problem lies in the Wei–Norman theorem context [13]. In the case of the *su*(2) algebra, it is well-known that the direct factorization problem is equivalent to solving a parametric oscillator-like equation whose solutions are known in a limited number of cases [14–17]. Thus, in the inverse solution proposed in [12], one departs from certain known solutions of such an equation and then the driving fields are determined. Accordingly, the dynamics is sensitive to the nature of such solutions.

The main purpose of the present paper is twofold. First, we obtain new families of analytically solvable driving fields using the formalism presented in [12]. Within this framework we show that the dynamics of one qubit interacting with such a family is closely related to the Mathieu functions properties. Second, we present an analysis on the correspondent dynamics. The interest in Mathieu functions lies in the fact that they are widely used in several branches of physics. They have emerged throughout physics mainly in systems with elliptic symmetry or those involving periodic or oscillating behavior. Initially, they were found by Mathieu in 1868 when he studied the oscillating modes of an elliptical membrane [18]. Since then, extensive theoretical work has been done studying their mathematical properties [19]. This has permitted the application of Mathieu functions to study diverse kinds of systems such as elliptical optical waveguides [20], propagation of invariant optical beams [21], the dynamics of electrons in periodic lattices [22], spectral singularities in PT symmetric potentials in quantum mechanics [23], and Aharonov–Bohm oscillations [24], among others.

The present work is organized as follows. In Section 2, the method to generate exactly solvable driving fields is revisited. We discuss in Section 3 the dynamics of a two-level system interacting with a family of precessing fields with oscillating amplitude generated using this formalism. When the dynamics is solved, we can request the condition to get a cyclic evolution for any initial state in the form of an evolution loop [25,26]. In Section 4, we get the prescriptions of such effect in terms of the dynamical parameters warranting the cyclic behavior for any initial state under the dynamics. An additional analysis on the geometric phase behavior is then conducted. Finally, conclusions and some perspectives are presented in Section 5.

#### **2. Analytically Solvable Driving Fields**

This Section is devoted to presenting the formalism to construct driving fields for which the time-evolution operator is exactly factorized. We first consider the two-level system ruled by the time-dependent Hamiltonian

$$H\_2(t) = \Delta \sigma\_0 + V(t)\sigma\_+ + \overline{V}(t)\sigma\_-,\tag{1}$$

where <sup>Δ</sup> <sup>∈</sup> <sup>R</sup> and {*σ*0, *<sup>σ</sup>*±} are the three generators of the *su*(2) algebra defined in terms of the three traditional Pauli operators as follows: *σ*<sup>0</sup> = <sup>1</sup> <sup>2</sup>*σz*, *<sup>σ</sup>*<sup>±</sup> <sup>=</sup> <sup>1</sup> <sup>2</sup> (*σ<sup>x</sup>* ± *iσy*). Furthermore, the interaction of the qubit with a classical field is described by the complex-valued function *V*. The correspondent time-evolution operator *U* is a solution of the equation

$$i\frac{d\mathcal{U}(t)}{dt} = H\_2(t) \cdot \mathcal{U}(t), \quad \mathcal{U}(0) = \mathbb{I}.\tag{2}$$

As the Hamiltonian (1) is a linear combination of the *su*(2) algebra generators, the Wei–Norman theorem [13] establishes that the time-evolution operator can be written in the form

$$
\Delta I(t) = e^{a(t)\sigma\_+} e^{\Delta f(t)\sigma\_0} e^{\oint (t)\sigma\_-},\tag{3}
$$

where the factorization functions *α*, *f* , and *β* satisfy the following coupled system of non-linear equations:

$$\begin{aligned} \alpha' - \alpha \Delta f' - \beta' e^{-\Delta f} &=& -iV, \\ \Delta f + 2\alpha \beta' e^{-\Delta f} &=& -i\Delta, \\ \beta' e^{\Delta f} &=& -i\overline{V}, \end{aligned} \tag{4}$$

with the initial conditions *α*(0) = *f*(0) = *β*(0) = 0. By developing the *su*(2) algebra in the exponential operators, it is an easy task to demonstrate that

$$\begin{split} \mathcal{U}(t) &= \, ^t\mathbb{I}(\cosh\frac{\Delta f(t)}{2} + \frac{a(t)\beta(t)}{2}e^{-\frac{\Delta f(t)}{2}}) + 2\sigma(\sinh\frac{\Delta f(t)}{2} + \frac{a(t)\beta(t)}{2}e^{-\frac{\Delta f(t)}{2}}) \\ &+ \sigma\_+ a(t)e^{-\frac{\Delta f(t)}{2}} + \sigma\_- \beta(t)e^{-\frac{\Delta f(t)}{2}}, \end{split} \tag{5}$$

where the operator *U*(*t*) in (5) has linked conditions fulfilled by *α*(*t*), *β*(*t*), and Δ*f*(*t*). In fact, the unitary condition on *U*(*t*): *U*(*t*)*U*†(*t*) = I implies that conditions *α*(*t*) = 0, *β*(*t*) = 0, and Δ*f*(*t*) are pure imaginary and fulfill together if at least two of them are requested, which can be proved directly.

Analytical solutions of System (4) are found in a limited number of cases [14]. However, in [12], a formalism to generate exact solutions for the disentangling problem has been developed. In such an approach, the radiation field is given by *V*(*t*) = *e*−*i*Δ*<sup>t</sup> R*(*t*), where the function *R* reads

$$R(t) = \frac{R\_0}{\mu^2(t)} \exp\left[i\lambda \int\_0^t \frac{ds}{\mu^2(s)}\right];\tag{6}$$

here, the real parameter *λ* is going to be fixed by the initial conditions and the real-valued function *μ* satisfies the Ermakov equation

$$
\mu''(t) + \Omega^2(t)\mu(t) = \frac{\Omega\_0^2}{\mu^3(t)}, \quad \Omega\_0 = \left[|R\_0|^2 + \frac{\lambda^2}{4}\right]^{1/2}.\tag{7}
$$

For reasons to be clarified later, the real number Ω<sup>0</sup> is called the generalized Rabi frequency. Furthermore, the initial conditions are determined by the logarithm derivative of (6) at *t* = 0

$$\frac{R\_0'}{R\_0} = i\frac{\lambda}{\mu\_0^2} - 2\frac{\mu\_0'}{\mu\_0} \,\tag{8}$$

where *R* <sup>0</sup> := *R* (0), *μ* <sup>0</sup> = *μ* (0) and without loss of generality we take *μ*<sup>0</sup> := *μ*(0) = 1. Thus, once [ln *R*(0)] and *R*<sup>0</sup> are provided as free parameters, the initial conditions are established. Indeed, given *<sup>δ</sup>*1, *<sup>δ</sup>*<sup>2</sup> <sup>∈</sup> <sup>R</sup> such that (ln *<sup>R</sup>*0) <sup>=</sup> *<sup>δ</sup>*<sup>1</sup> <sup>+</sup> *<sup>i</sup>δ*2, it is found that *<sup>λ</sup>* <sup>=</sup> *<sup>δ</sup>*<sup>2</sup> and *<sup>μ</sup>* <sup>0</sup> = −*δ*1/2. It is a well-known fact that a particular solution to Equation (7) is related to the parametric-like equation [27,28]:

$$
\varphi''(t) + \Omega^2(t)\varphi(t) = 0,\tag{9}
$$

where Ω<sup>2</sup> is a real-valued function and the initial conditions read

$$\lim\_{t \to 0} \frac{1}{R(t)} \left[ \frac{\varphi'(t)}{\varphi(t)} + \frac{1}{2} \frac{R'(t)}{R(t)} \right] = 0, \quad \varphi(0) = 1. \tag{10}$$

The former expression is equivalent to the initial condition *α*(0) = 0 (see [14] for details). According to Pinney (alternatively Ermakov), if *u* and *v* are two linear independent solutions of (9), then

$$\mu(t) = \left[\mu^2(t) - \Omega\_0^2 \mathcal{W}^{-2} \upsilon^2(t)\right]^{1/2} \tag{11}$$

is a solution of (7). Furthermore, *u*(*t*0) = *μ*(*t*0), *u* (*t*0) = *μ* (*t*0), *v*(*t*0) = 0, *v* (*t*0) = 0, and *W* stands for the Wronskian of *u* and *v*. In addition, note that *R*(0) = *R*<sup>0</sup> is also a free parameter to be specified and the function *ϕ* solves the factorization problem [12,14]. Indeed,

$$a(t) = \frac{i}{R\_0} \mu^2(t) \exp\left[-i\Delta t - i\lambda \int\_0^t \frac{ds}{\mu^2(s)}\right] \left[\frac{q^\prime(t)}{q(t)} - \frac{\mu^\prime(t)}{\mu(t)} + \frac{i\lambda}{2\mu^2(t)}\right],\tag{12}$$

$$\beta(t) = -i\mathbb{R}\_0 \int\_0^t \frac{ds}{q^2(s)},\tag{13}$$

$$
\Delta f(t) = \ln \left[ \frac{\mu^2(t)}{\rho^2(t)} \right] - i\lambda \int\_0^t \frac{ds}{\mu^2(s)} - i\Delta t. \tag{14}
$$

The time-evolution of a state can be straightforwardly computed. Because an instance if the qubit's initial state is |0, at any time the state of the system reads

$$|\psi(t)\rangle = e^{-\Lambda f(t)/2} [e^{\Lambda f(t)} + a(t)\beta(t)]|0\rangle + e^{-\Lambda f(t)/2} \beta(t)|1\rangle. \tag{15}$$

An important physical observable is the atomic population inversion defined as *P*(*t*) := *σ*0. For a general state |*ψ*(*t*) = *c*0|0 + *c*1|1, the population inversion reads *P*(*t*) = |*c*0| <sup>2</sup> − |*c*1<sup>|</sup> 2. This quantity is defined on the interval [−1, 1]. In fact, for the largest (lowest) possible value of *P*(*t*), the state is in the excited (ground) state with certainty. In general, for positive (negative) values of the population inversion the probability of finding the state in the upper (lower) energy state is larger. In the Rabi model, the population inversion is a periodic function of time where field amplitude and detuning (the gap between the field frequency and the spacing level energy of the atom) determine the oscillation period [12,14,29].

We finish this section by summarizing the method: One should start with the second-order differential Equation (9) choosing Ω such that the solutions are known. The function *μ* is then obtained via the Ermakov–Pinney solution and so it is possible to generate families of analytically solvable driving fields (6) for which the initial conditions constitute free parameters to control the dynamics. Some examples have been reported in [12]. For instance, if Ω<sup>2</sup> is a negative real constant, the correspondent driving field is a decaying one. However, a precessing field with oscillating amplitude is achieved if Ω<sup>2</sup> is a positive real constant Ω1. Such a family of control fields constitutes a generalization of the circularly polarized field, which is obtained as a particular case of this model.

#### **3. Dynamics in a Precessing Field with Oscillating Amplitude**

In this section, we report a new family of driving fields for which the time-evolution operator is exactly disentangled. We consider the function <sup>Ω</sup>2(*t*) = *<sup>ω</sup>*<sup>0</sup> <sup>−</sup> *<sup>ω</sup>*<sup>1</sup> cos2 (*t*) in the parametric oscillator-like Equation (9), which can be written as a Mathieu equation [19,30] as follows:

$$[\varphi''(t) + [a - 2q \cos(2t)]\varphi(t) = 0,\tag{16}$$

where the identity 2 cos2 *t* = 1 + cos(2*t*) has been used and

$$a = \omega \upsilon - \frac{1}{2}\omega\_1, \qquad \qquad q = \frac{1}{4}\omega\_1. \tag{17}$$

Furthermore, the following initial conditions read [ln *<sup>R</sup>*(0)] <sup>=</sup> *<sup>i</sup>δ*, where *<sup>δ</sup>* <sup>∈</sup> <sup>R</sup> and *<sup>ϕ</sup>*<sup>0</sup> <sup>=</sup> 1. According to (11), the function *μ* can be written as

$$
\mu^2(t) = \mu^2(t) + \left(\frac{\Omega\_0 v(t)}{W}\right)^2,\tag{18}
$$

where Ω<sup>2</sup> <sup>0</sup> = |*g*| <sup>2</sup> <sup>+</sup> *<sup>δ</sup>*2/4 has been determined considering *<sup>R</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup>*ig*, *<sup>g</sup>* <sup>∈</sup> <sup>C</sup>. Note that these particular expressions for the initial conditions retrieve the well-known Rabi frequency in the problem of a qubit interacting with a circularly polarized field. Furthermore, it is easy to prove that the Wronskian *W* does not depend on *t*.

In most cases, the two linearly independent solutions of the Mathieu Equation (16) can be written in terms of the even and odd Mathieu functions *MC*(*a*, *q*, *t*) and *MS*(*a*, *q*, *t*) [19]. As is clarified in the next section, we restrict ourselves to this case.

#### *3.1. The Dynamics of Driving Fields and the Theory of Mathieu Functions*

In this section, we study the behavior of the solution *μ* (18) of the Ermakov Equation (7) as a function of *t*, but also as a function of its parameters. The behavior of the driving fields can be very different depending on the values of the frequencies *ω*<sup>0</sup> and *ω*1. Moreover, *R*(*t*) can be a periodic, non-periodic, oscillating, or pulsed-like function of time. In this respect, the theory of Mathieu functions will be useful.

Given a value of the parameter *q* of the Mathieu Equation (16), there exist values of *a* for which the Mathieu functions can be periodic, bounded, or unbounded [19]. When the Mathieu function *MC* (*MS*) is periodic, the values of *a* are called *characteristic values* and they are denoted by *ar* (*br*), becoming functions of the parameter *q*

$$a\_{\varGamma} = r^2 + \sum\_{i=1}^{\infty} a\_i(r) \, q^i, \qquad \qquad b\_{\varGamma} = r^2 + \sum\_{i=1}^{\infty} \beta\_i(r) \, q^i. \tag{19}$$

where *<sup>α</sup>i*(*r*), *<sup>β</sup>i*(*r*) <sup>∈</sup> <sup>R</sup> are given in [19]. The parameter *<sup>r</sup>* <sup>≥</sup> 0 is called the *characteristic exponent* and it determines the periodic properties of the Mathieu functions. The even Mathieu function *MC*(*ar*, *q*, *t*) is periodic with period *π* or 2*π* if the characteristic exponent *r* can be written as 2*n*, or 2*n* + 1, respectively, where *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> ∪ {0}. The odd Mathieu function *MS*(*br*, *<sup>q</sup>*, *<sup>t</sup>*) is periodic with period *<sup>π</sup>* or 2*<sup>π</sup>* if *<sup>r</sup>* can be written as 2*n* + 2, or 2*n* + 1, respectively [19]. When *r* is an integer *ar* = *br*, then the Mathieu functions *MC*(*ar*, *q*, *t*) and *MS*(*br*, *q*, *t*) are not solutions to the same Mathieu equation. On the other hand, when *r* is not an integer, *ar* = *br* holds. Then, *MC*(*ar*, *q*, *t*) and *MS*(*ar*, *q*, *t*) are solutions to the same Mathieu equation, and moreover, they are linearly independent [19]. Considering *r* as a positive rational number, then, it can be written as *r* = *m* + *p*/*s*, where *m* ≥ 0 is an integer and *p* = 0, *s* are relative prime integers. Hence, Mathieu functions *MC*(*ar*, *q*, *t*) and *MS*(*ar*, *q*, *t*) are periodic with period 2*πs*, for *s* > 2. On the other hand, if *r* is an irrational number, Mathieu functions are not periodic but bounded functions which do not decay to zero as *t* → ∞ [19].

We define the *q* − *a* space as the two dimensional space in which *q* and *a* are the abscissa and ordinate, respectively. The behavior of the Mathieu functions is dictated by the region of the *q* − *a* space in which the point (*a*, *q*) lies, and so the discussion of the previous paragraph can be summarized into a picture in the *q* − *a* space as follows. As mentioned above, the characteristic values are functions of *q*, therefore in the *q* − *a* space, the curves of *ar* and *br* lie there as functions of *q*. Such curves are from now on called the *characteristic curves*. In Figure 1a, the behavior of the characteristic curves can be observed in the space *q* − *a* for different integer values of *r*, while in Figure 1b, the case is shown with several non-integer values of *r*. The characteristic curves for all non-integer positive values of *r*, a case in which *ar* = *br*, saturate a certain region of the *q* − *a* space [19] which is shown shaded in Figure 1b. This region, along with the characteristic curves of integer values of *r*, is known as the region of *stability* of the Mathieu functions. A Mathieu function is called *stable* if it is bounded or tends to zero as *t* → +∞, as well as it is called *unstable* if it diverges as *t* → +∞. Therefore, if the point (*a*, *q*) lies inside (outside) the region of stability, the corresponding Mathieu function is a bounded (unbounded) function of *t*. Regarding Equation (17), it can be observed that *ω*<sup>0</sup> can be written as *ω*<sup>0</sup> = *a* + *ω*1/2. Taking *a* as a characteristic value *ar*(*q*), *ω*<sup>0</sup> can be considered as a characteristic frequency (value), in the sense that

$$
\omega\_0(\omega\_1, r) = a\_r(\omega\_1) + \omega\_1/2,\tag{20}
$$

in the analog *ω*<sup>1</sup> − *ω*<sup>0</sup> space, since *q* = *ω*1/4. Figure 2 shows the behavior of the characteristic curves *ω*0(*ω*1,*r*) in the *ω*<sup>1</sup> − *ω*<sup>0</sup> space.

**Figure 1.** (**a**) The characteristic values *ar*(*q*) (blue-dashed curves) for *r* = 0, 1, 2, 3, 4, 5, and *br* (red curves) for *<sup>r</sup>* <sup>=</sup> 1, 2, 3, 4, 5, as functions of *<sup>q</sup>*, for several values of *<sup>r</sup>* <sup>∈</sup> <sup>N</sup>. Note that *ar*(0) = *br*(0) = *<sup>r</sup>*2. (**b**) The shaded region is produced by all the characteristic curves with non-integer characteristic values *r*, case in which *ar* = *br*.

**Figure 2.** The characteristic frequency *ω*0(*ω*1,*r*) = *ar*(*ω*1) + *ω*1/2 as functions of *ω*1, for several rational values of *r*.

Nevertheless, we are interested in the two linearly independent solutions of the Mathieu equation conforming to the Ermakov equation solution *μ*. Then, we must separate the *ω*<sup>1</sup> − *ω*<sup>0</sup> space into two regions: the one in which *both* linearly independent solutions are stable, which will be denoted as A; and the one in which *at least* one of the two solutions is unstable, represented by <sup>A</sup>*C*. In region <sup>A</sup>, the function *<sup>μ</sup>* is a bounded function of *<sup>t</sup>*, while in the region <sup>A</sup>*C*, *<sup>μ</sup>* is unbounded, as can be observed from Equation (18).

As mentioned above, within the stability region, there exist two cases: (*i*) the characteristic exponent *r* is not an integer; and (*ii*) it is an integer. In case (*i*), both Mathieu functions *MC* and *MS* are linearly independent stable solutions of the Mathieu Equation (16). In case (*ii*), the second linearly independent solution *MC* as well as the second linearly independent *MS* are unstable [19]. On the other hand, in the instability region, both linearly independent solutions of the Mathieu equation are unbounded. Therefore, A can be described as the region of the *ω*<sup>1</sup> − *ω*<sup>0</sup> space in which the characteristic values *r* are positive non-integers (see Figure 2), namely

$$\mathcal{A} = \{ (\omega\_1, \omega\_0) \in \mathbb{R}^2 \mid \omega\_0(\omega\_1, r) = a\_\mathcal{I}(\omega\_1) + \omega\_1/2, r > 0, r \notin \mathbb{Z} \}, \tag{21}$$

where *ar*(*ω*1) is the characteristic value as function of *ω*1(*q*), see Equation (17). Let us denote by B the set of points of the characteristic curves for integer values of *r* > 0 (denoted in Figure 1a), so the stability region in the *ω*<sup>1</sup> − *ω*<sup>0</sup> space can be written as A∪B (note that the set B is part of the boundary of A (see Figure 1)). Since this is case (*ii*), one of the solutions to the Mathieu equation is unstable, then B⊂A*C*. Hence, within <sup>A</sup>*C*, there are two cases: the point (*ω*1, *<sup>ω</sup>*0) belongs to the unstable region or to B.

For simplicity, we focus on the case in which *μ* can be written in terms of the Mathieu functions *MC* and *MS*. This case includes the whole region <sup>A</sup> as well as <sup>A</sup>*<sup>C</sup>* − B (unstable region case), that is to say, we are excluding the points in B. Nevertheless, as the second solutions in B are unstable, the asymptotic behavior of *<sup>μ</sup>* for *<sup>t</sup>* <sup>→</sup> <sup>∞</sup> for (*ω*1, *<sup>ω</sup>*0) ∈ A*<sup>C</sup>* is similar whether (*ω*1, *<sup>ω</sup>*0) belongs to <sup>B</sup> or not.

Since Equation (16) is invariant under a complex conjugate operation the real or imaginary part of any of its solutions is also a solution and being *μ* a real-valued function, it can be written as

$$
\mu^2(t) = c^2 \text{Re}[M\_\mathbb{C}(t)]^2 + \frac{\Omega\_0^2}{\mathcal{W}^2} \text{Re}[M\_\mathbb{S}(t)]^2,\tag{22}
$$

where we have defined *u*(*t*) = *c*Re[*MC*(*t*)], *c* = 1/Re[*MC*(0)], and *v*(*t*) = Re[*MS*(*t*)] (the dependence on the frequencies *ω*0, *ω*<sup>1</sup> has been omitted for shortness) in order to satisfy the initial conditions stated in Equation (10) and the ones mentioned after Equation (18). Furthermore, *W* denotes the Wronskian of two functions. Then, the driving field *R*(*t*) can be written as

$$R(t) = \frac{-i\vec{g}}{c^2 \text{Re}[M\_\mathbb{C}(t)]^2 + \frac{\Omega\_0^2}{W^2} \text{Re}[M\_\mathbb{S}(t)]^2} \exp\left(i\delta \int\_0^t \frac{ds}{c^2 \text{Re}[M\_\mathbb{C}(s)]^2 + \frac{\Omega\_0^2}{W^2} \text{Re}[M\_\mathbb{S}(s)]^2}\right). \tag{23}$$

In order to write the factorizing functions, we have to obtain the solution to the parametric oscillator-like Equation (16), satisfying the initial conditions (8), which can be written as

$$\varphi(t) = \frac{1}{M\_{\mathcal{C}}(0)} M\_{\mathcal{C}}(t) - \frac{i\delta}{2M\_{\mathcal{S}}'(0)} M\_{\mathcal{S}}(t),\tag{24}$$

where *M <sup>S</sup>*(*t*) = *dMS*(*t*)/*dt*. The factorizing functions can be written as

$$a(t) = \frac{\mu^2(t)}{-i\bar{\mathfrak{g}}} \left[ \frac{\varrho'(t)}{\varrho(t)} - \frac{\mu'(t)}{\mu(t)} + \frac{i\delta}{2\mu^2(t)} \right] e^{-i(\Delta t + \mu\_1(t))},\tag{25}$$

$$\beta(t) = R\_0 \int\_0^t \frac{ds}{\varphi^2(s)},\tag{26}$$

$$
\Delta f(t) = \ln \left[ \frac{\mu^2(t)}{\varphi^2(t)} \right] - i\mu\_1(t) - i\Delta t\_\prime \tag{27}
$$

where

$$
\mu\_1(t) = \delta \int\_0^t \frac{ds}{\mu^2(s)}. \tag{28}
$$

Now, we are able to calculate the time-evolution of the initial state *ψ*(0) = |1 given in (15). As commented previously, population inversion, *P*(*t*) ∈ [−1, 1] is a physical quantity very useful and descriptive for two-level systems with respect to the system dynamics. It depicts, for the quantum state, the dynamical variation between the excited state |0 (*P*(*t*) = 1) and the base state |1 (*P*(*t*) = −1). Thus, in our case, the population inversion can be written after some algebra as

$$P(t) = \mu^2(t) |\varphi(t)\beta(t)|^2 \left\{ \left| \frac{1}{\beta(t)\varphi^2(t)} + \frac{a(t)}{\mu^2(t)} e^{i[\Delta t + \mu\_1(t)]} \right|^2 - \frac{1}{\mu^4(t)} \right\}.\tag{29}$$

#### 3.1.1. Driving Fields in the Region A

In this case, the solution *μ* to Ermakov equation is a time oscillating function that can be periodic if the characteristic exponent *r* is a rational number. Hence, the driving field *R*(*t*), as well as *V*, is an oscillating function which becomes periodic when the parameters *g*, *δ*, Δ, and the characteristic exponent *r* satisfy the condition [12]

$$p\delta \int\_0^\pi \frac{ds}{\mu^2(s)} + \Delta \tau \equiv 0 \pmod{2\pi},\tag{30}$$

where *p* is a natural number and *τ* is the period of *μ*. Solutions exhibit lots of possibilities due to the set of physical parameters *a*, *q*, *g*, *δ*, Δ, and *τ* (as well as *ω*1, *ω*<sup>2</sup> instead of *a*, *q*). They could be selected or highlighted to get several possible control effects: periodicity in the field *R*(*t*), the reaching of evolution loops, or simply the control of population inversion. We analyze several of these effects separately in the following sections of the article.

Figure 3 shows three examples of periodic driving fields and their associated population inversion. In these cases, the characteristic exponent *r* is a rational number, so the associated *μ* is periodic and the parameters *g* and *δ* are such that *R*(*t*) is a periodic function of time. Furthermore, in order to show the time flow correspondence between the left plots with the right ones, a color scale from blue to red has been used. As can be observed from Figure 3b,e, each minimum of the population inversion barely corresponds to a local maximum of the transverse driving field amplitude |*R*(*t*)|. Furthermore, the non-negative parameter *p* is related to the number of "petals" in each flower-like structure of the driving field. In [12], it is shown that our model can describe a two-level system interacting with a magnetic transverse field *B*⊥. It is easy to show that real and imaginary parts of *R*(*t*) are nothing other than the field components performing a certain rotation. This fact is due to the transverse field *B*<sup>2</sup> <sup>⊥</sup> <sup>=</sup> *<sup>B</sup>*<sup>2</sup> *<sup>x</sup>* + *B*<sup>2</sup> *<sup>y</sup>* = Im[*R*(*t*)]<sup>2</sup> + Re[*R*(*t*)]<sup>2</sup> causing the state to instantaneously precess around the direction of the entire field, the direction of such precession becomes more horizontal while <sup>|</sup>*B*⊥||*Bz*<sup>|</sup> <sup>=</sup> <sup>|</sup>Δ<sup>|</sup> <sup>2</sup> . Thus, when the relative orientation between the field direction and the state direction on the Bloch sphere becomes more perpendicular, the states rotates more extremely on the sphere, being able to reach a more extreme position there, just as it happens in the case of a circularly polarized field [31]. These rotations becomes more effective the bigger the field is because the instantaneous angular frequency depends on it. When the value of the parameter *δ* is changed, with the others held constant, an oscillatory non-periodic driving field is obtained (see the continuous blue curves in Figure 4). In this case, the changes in the amplitude of the driving field are not fast enough for the population inversion to reach its minimum value −1. Now, by only increasing the parameter *g* and holding the others fixed (the magnitude of *B*<sup>⊥</sup> is proportional to |*g*|), we observe in Figure 4 (red dashed curves) that the changes in the amplitude of the driving field are increased with respect to the previous case, pushing the increment of the range of the peaks of the population inversion towards (−1, 1). The positions of the maxima of the population inversion remain unchanged. We can conclude, in this case, that as the parameter *g* increases, the peaks of the population inversion are sharpened, presenting a resonant-like behavior. Furthermore, in any of the previous cases, as shown at the end of Section 3.1, the population inversion is periodic.

**Figure 3.** The driving field (23) in the complex plane for *g* = 1, Δ = 1, (**a**) *r* = 1.5, *ω*<sup>1</sup> = 4, *ω*<sup>0</sup> = 4.53718, *δ* = 0.7071, (**b**) *r* = 3.5, *ω*<sup>1</sup> = 4, *ω*<sup>0</sup> = 14.2946, *δ* = 0.28867, (**c**) *r* = 3.5, *ω*<sup>1</sup> = 40, *ω*<sup>0</sup> = 36.2607, *δ* = 0.28867. In (**d**–**f**), the corresponding population inversion is shown. All of the previous values of *ω*<sup>0</sup> correspond to characteristic values of the Mathieu functions (16) with rational values of *r*, that is, (*ω*1, *ω*0) ∈ A. Time flows from blue to red matching the scale between the types of plots.

**Figure 4.** The driving fields (23) in the complex plane (left column) and the associated population inversion (right column) for the cases of Figure 3, except for the values of *g* and *δ*, which has been increased as follows: (**a**,**d**) *g* = 1, *δ* = 2 (blue), and *g* = 5, *δ* = 2 (red dashed), (**b**,**e**) *g* = 1, *δ* = 5 (blue), and *g* = 7, *δ* = 5 (red dashed), (**c**,**f**) *g* = 1, *δ* = 4 (blue), and *g* = 7, *δ* = 4 (red dashed). The driving field precesses counterclockwise.

### 3.1.2. Driving Fields in the Region <sup>A</sup>*<sup>C</sup>*

As mentioned above, at least one of the Mathieu functions, which conform the solution *μ* to the Ermakov equation, is an unbounded function diverging as *t* → ∞. In this case, we have observed that the zeros of the Mathieu functions *MC*(*t*) and *MS*(*t*) tend to approach each other more and more as *t* increases, a phenomenon known as condensation of zeros [19]. This gives rise to pulsed-like peaks in the corresponding driving field *<sup>R</sup>*(*t*) (|*R*(*t*)<sup>|</sup> <sup>∝</sup> 1/*μ*2), as can be observed in Figure 5a. Therefore, as time passes, the pulses increase their size and become narrower causing a pulsed-like behavior of the population inversion, as can be observed in Figures 5b and 5c.

**Figure 5.** (**a**) The driving field (23) in the complex plane. (**b**) Its square modulus as a function of time. The inset is a zooming of the plot for *t* ∈ [0, 10] to show the correspondent local maxima, and (**c**) the associated population inversion with *ω*<sup>1</sup> = 1, *ω*<sup>0</sup> = 1.5, *g* = 1, *δ* = 0.2887, where (*ω*0, *ω*1) ∈ A/ . The driving field precesses counterclockwise.

#### **4. Evolution Loops, Cyclic Evolution and Phases**

Evolution loops are a control phenomena present in certain quantum systems admitting dynamic solutions to the restriction *U*(*τ*) = *eiφ*I for some *τ* > 0, which clearly revert the evolution after time *τ*. Furthermore, the parameter *φ* is the global phase acquired during the process. It is relevant because that reduction becomes independent from the initial states, thus, any initial state evolves after such time. These phenomena boost other intermediate effects which have value in quantum control. In this section, we show that the field being considered admits this kind of effects. In the current case, because if *α*(*τ*) = 0 and *β*(*τ*) = 0, then Δ*f*(*τ*) is automatically pure imaginary, then, if in addition, *<sup>φ</sup>* <sup>=</sup> <sup>Δ</sup>*f*(*τ*) = <sup>4</sup>*nπi*, *<sup>n</sup>* <sup>∈</sup> <sup>Z</sup>, *<sup>U</sup>*(*τ*) reduces exactly to <sup>I</sup>. Such is the condition to reach evolution loops for some *t* = *τ* (otherwise, if only *α*(*τ*) = 0 = *β*(*τ*) are fulfilled, the evolution reaches the initial state with a non-zero phase). Then, the procedure to find such prescriptions is based in the finding of *τ* satisfying *<sup>α</sup>*(*τ*) = <sup>0</sup> <sup>=</sup> *<sup>β</sup>*(*τ*), then by using (27), we can adjust the value of <sup>Δ</sup> to reduce <sup>Δ</sup>*f*(*t*) to 4*nπi*, *<sup>n</sup>* <sup>∈</sup> <sup>Z</sup>. This is always possible due to the first two terms on the right side being independent of Δ. Despite this, there are a lot of solutions because there are five free parameters: *a*, *q* (otherwise *ω*0, *ω*1), *δ*, *g*, *t*.

Figure 6 shows three examples of such an effect from a large variation in the initial state |1 until a smaller one in terms of parameters *a*, *q*, *g*, *δ*, Δ, and *τ*. Those parameters work in a combined way to give dynamical effects, thus, they are not related with the periodicity of *R*(*t*), *P*(*t*) of *U*(*t*) in an exclusive way. In this sense, examples being considered are only illustrative among a wide variety of possibilities (obtained numerically upon the conditions *α*(*τ*) = *β*(*τ*) = 0), but exhibiting different ranges of variation for *P*(*t*). In any case, note that *g* and *δ* play an important role in such variation, as is expected. Plots on the left exhibit the value of *P*(*t*) while plots on the right show the dynamics represented on the Bloch sphere. Nevertheless, evolution loops are independent from the initial states, examples consider the initial state as the ground state, |1. The color changes from blue to red in agreement with the *P*(*t*) value. Note that the intermediate reaching of unitary multiples of |1 with non-zero phases until the zero phase is reached in *τ*. The last example in Figure 6c meets the evolution loop under periodic driving field, generating a more regular dynamics around |1.

**Figure 6.** Population inversion together with the dynamics represented on Bloch spheres for several prescriptions generating evolution loops. In any case, the color depicts the value of *P*(*t*) in agreement with the left plot. (**a**) *a* = 3.37813, *q* = 2.6118, *g* = 0.53635, *δ* = 3.30558, Δ = 0.190982, *t* = 11.3718. (**b**) *a* = 3.13268, *q* = 2.70972, *g* = 1.93404, *δ* = 1.41305, Δ = 0.74138, *t* = 18.7095. (**c**) *a* = 3.01588, *q* = 0.022235, *g* = 0.0123357, *δ* = 1.0189, Δ = 0.00100538, *t* = 18.3134 depicting a tiny regular evolution loop with a periodic field around of |1.

#### *Dynamical and Geometric Phases*

It is well-known that a geometric phase exists related to a quantum system cyclic evolution governed by a slow change of parameters [32,33]. The relevance of such a phase factor lies in the foundations of quantum theory; however, it has also recently found important applications in quantum information and computation, as a part of the system evolution [34,35].

Consider the dynamic phase *φd*, as well as the Aharonov–Anandan geometric phase *φ<sup>g</sup>* in evolution loop dynamics [36,37], which depicts the phase when the system returns to its initial physical state at the time *τ* (in the current case to the state |0):

$$\phi\_d = -\int\_0^\tau \langle \Psi(t) | H\_2(t) | \Psi(t) \rangle dt \quad \rightarrow \quad \phi\_\S = \arg(\langle \Psi(0) | \Psi(\tau) \rangle) - \phi\_d \equiv \phi - \phi\_d. \tag{31}$$

which even in the simplest case of a time-independent magnetic field, a non-trivial value is obtained: *φ<sup>d</sup>* = −*π* cos *θ*<sup>0</sup> and *φ<sup>g</sup>* = *π*(cos *θ*<sup>0</sup> − 1), where *θ*<sup>0</sup> is the relative angle between the magnetic field and the initial Bloch vector. These values are indeed independent of the field amplitude [38]. Conversely, in time-dependent cases, the dynamical trajectory on the Bloch sphere can be manipulated with respect to the field parameters, yielding different values of the geometric phase. Such is our case, as shown in the examples in the Figure 6: (a) *φ<sup>d</sup>* = −38.057, *φ<sup>g</sup>* = 126.022, *φ* = 28*π*, (b) *φ<sup>d</sup>* = −2.19494, *φ<sup>g</sup>* = 52.460, *φ* = 16*π*, and (c) *φ<sup>d</sup>* = −13.4832, *φ<sup>g</sup>* = 88.881, *φ* = 24*π*. In these examples, the global phase *φ* has been selected as the minimum value assuring the positivity for Δ. The present analysis represents a first step in the construction of logical gates for holonomic quantum computation [34,35]. These results will be reported elsewhere.

#### **5. Conclusions**

We have presented the study of driving fields generated by the inverse approach given in [12], departing from a sinusoidal parametric oscillator-like equation. Regarding the solution of the associated Ermakov equation, the driving fields and the population inversion can be written in terms of the Mathieu functions. We have shown that the theory of Mathieu functions is determinant in the dynamical analysis of such driving fields, as a clear split in two regions for the space of the frequencies *ω*<sup>0</sup> and *ω*<sup>1</sup> is acquired. Thus, different dynamical behaviors of the driving fields are shown. In the region A, we have shown that the driving fields have an oscillatory nature, but they can still be periodic or non-periodic depending on the value of the parameters *g* and *δ* given by a precise prescription obtained analytically. In this same region, the population inversion is periodic provided that the associated Mathieu functions are also periodic. In the complement of the previous region, for example <sup>A</sup>*C*, the behavior of the driving fields is pulse-like; with the amplitude and sharpness of their peaks increasing in time, yielding a resonant-like behavior for the associated population inversion. We consider that this system presents the advantage of having several kinds of dynamics that can be prescribed on demand by fine tuning the frequencies *ω*<sup>0</sup> and *ω*<sup>1</sup> (alternatively *a* or *q*) of the parametric oscillator-like potential. Together, we showed that evolution loop solutions are numerically achievable and the correspondent geometric phases can be calculated. Our study represents the first stage in the development of technological applications involving single-qubit devices. For instance, in some many-body models, precise selective control operations on single qubits are required, which could be achievable through the driving fields obtained in this work. Such analysis will be reported elsewhere. We hope our results shed some light on the matter.

**Author Contributions:** Investigation, M.E., A.J.-N. and F.D.; Writing—original draft, M.E., A.J.-N. and F.D.; Writing—review & editing, M.E., A.J.-N. and F.D.

**Funding:** The support of CONACyT through project A1-S-24569 is acknowledged.

**Acknowledgments:** M.E. and F.D. would like to acknowledge to School of Engineering and Science of Tecnologico de Monterrey as well as the financial support of Novus Grant with PEP no. PHHT023\_17CX00001, TecLabs, Tecnologico de Monterrey, Mexico, in the production of this work. A.J.-N. thanks the support of INAOE and CICESE.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Asymmetry of Quantum Correlations Decay in Nonlinear Bosonic System**

#### **Anna Kowalewska-Kudłaszyk \*,† and Grzegorz Chimczak †**

Nonlinear Optics Division, Faculty of Physics, Adam Mickiewicz University, Uniwersytetu Pozna ´nskiego 2, 61-614 Pozna ´n, Poland

**\*** Correspondence: annakow@amu.edu.pl

† These authors contributed equally to this work.

Received: 12 July 2019; Accepted: 6 August 2019; Published: 8 August 2019

**Abstract:** We study the problem of the influence of one-sided different noisy channels to the quantum correlations decay in a symmetric bosonic system. We concentrate on one type of these correlations—the entanglement. The system under consideration is composed of two nonlinear oscillators coupled by two-boson interactions and externally driven by a continuous coherent field. Our low-dimensional system can be treated as 2-qutrit one. Two different noisy channels (the amplitude and the phase-damping reservoirs) are applied to both of the system's modes. We show that there is a noticeable difference in the quantum entanglement in 2-qubit subspaces of the whole system decrease after swapping the reservoirs between the modes of the considered symmetric system. It appears also that the degree of obtained entanglement depends crucially on the position of the appropriate type of reservoir. The origin of the observed asymmetry is also explained.

**Keywords:** nonlinear oscillator; quantum entanglement; open system

#### **1. Introduction**

Obtaining low-dimensional optical systems is still one of the challenges in contemporary quantum optics. Systems evolving among a strictly limited number of quantum states play an important role, e.g., in quantum cryptography [1,2], state teleportation [3,4], and quantum information processing [5].

The methods leading to the state truncation rely mainly on the applications of nonlinear media of various types. In the literature we may find both theoretical and experimental examples of low-dimensional bosonic systems obtained with the help of effective nonlinearities of various orders and strengths [6–9]. One of such methods makes use of quantum effects leading to the phenomena known as photon blockade (or more generally, boson blockade). This is the phenomenon in which the Hilbert space of the considered system is significantly truncated in such a way that only a very limited number of bosonic states is available for the whole system. Such a truncation can be obtained in two different ways, and two different physical mechanisms lead to photon blockades.

In the first mechanism, which is known as *conventional photon blockade*, the origin of truncation lies in obtaining an effective nonlinearity in the Hamiltonian describing the system [10–13]. This nonlinearity is responsible for the non-equidistant eigenstates spectrum of the system. An external driving resonant with one transition is not resonant with others, and thus, the number of allowed states is limited (truncated). For example, the effective nonlinearity of a cavity field–atom interaction leads to the anharmonicity of the Jaynes–Cummings ladder of eigenstates. The presence of single photon in this optical cavity suppresses the appearance of the second one in the system, because a second photon at the same frequency is not

resonant with the next transitions in the ladder [10]. This mechanism works properly only if the second photon is detuned from each of those next transitions by amounts that are much larger than the respective linewidths. Therefore, strong nonlinearities, and thus strong light–matter regimes, are necessary in this first mechanism. Large nonlinearities may be problematic from an experimental point of view, but many examples have already been presented in which strong effective nonlinearities are generated, e.g., due to: cavity field–atom interactions in the dispersive limit [10], optomechanical interactions [14], interactions between the atoms in optical lattices [15], nanoresonators [16], and interactions in superconducting circuit-QED systems [17,18].

In the weak light–matter regime, where the effective nonlinearity is weak, the second mechanism (known as *unconventional photon blockade*) is necessary [19–21]. It is a strongly resonant effect and requires a minimum input intensity to operate [22]. The physical origin of this blockade effect is the result of multipath destructive interference, which is responsible for uncoupling of some of the states.

In the considerations presented in this paper, we apply the systems characterized by strong nonlinearities to make use of the physical mechanism leading to the conventional photon blockade. Such systems are also known in the literature as *nonlinear quantum scissors*, due to their properties resulting in truncation of the Hilbert space [23–25]. Due to the presence of strong nonlinearity, which changes the energy spectrum of the considered system, there is only a limited number of resonant states, and finally the effective state truncation is possible.

In many papers the limited number of two-mode states is the starting point of the discussion of entanglement formation. Such possibility was reported for example in papers dealing with nonlinear quantum scissors, which exploited strongly nonlinear media (with *χ*(3) much larger than other parameters describing the interactions) to obtain 2-qubit, qutrit-qubit or 2-qutrit systems [26–28]. Various external driving mechanisms were also analyzed: continuous excitation, performed by a laser field [24,26]; intense optical field generating two modes by parametric down-conversion process [27]. Pulsed laser driving with constant or different pulse widths was also analyzed [29]. In all these cases, generation of maximally or almost maximally entangled states were reported. It was also shown that changes in the duration of laser pulse or even random perturbations in pulse duration can improve the ability of the nonlinear system to produce entangled states [29]. It may be of special worth especially for real experimental situations when some perturbations in maintaining the precise and stable timing of pulse laser sources can occur.

In the present paper, a coupled oscillatory system evolving among limited number of quantum states externally driven in both modes by a coherent field is presented. We concentrate on an open system evolution of coupled nonlinear media and apply different (amplitude and phase) noisy channels to both modes of the considered system. We show that we may influence the rate of disentanglement and the values of entanglement created in the qubit subspaces of the whole system by swapping the noisy channels between the modes. Such changes in the environment type, which influence the entangled qubits, may occur when both parts of the entangled pair are sent through different paths and therefore may encounter various local disturbances. The observed asymmetry in the response of the analyzed 2-mode system may give some information on which type of environment the system's parts were exposed to and which of the qubits was affected by a specific type of noise. The influence of the reservoir and the reservoir combinations on the symmetry of quantum correlation were analyzed for example in [30] for some mixed states under the influence of depolarizing channel. For the specified X state, the problem of symmetry in correlations decay was reported for example in [31].

#### **2. The Model**

The systems we are dealing with consist of coupled quantum nonlinear media with external excitations, which in general can be of various types. Such systems are well known as the potential

sources of maximally or almost maximally entangled states. Various types of coupling and external driving have already been discussed [26–28,32]. In all of those cases it appeared that it was possible to restrict the evolution of the whole coupled system to several bosonic states only—therefore, under the assumptions of weak interactions (as compare with the nonlinearity), truncation of the appropriate wave function was possible. In consequence it has also been shown that such nonlinear systems can be treated as pairs of qubits, or qutrits or qutrit-qubit pairs. It has also been shown that under some assumptions, creation of for example Bell-type [24,26] or W-states [32] was possible.

Such models can be realized either in the quantum optical domain with media characterized by ultra-strong Kerr nonlinearities, or more often in the systems described by effective Hamiltonians including nonlinearities of the Kerr type.

Each of the parts forming the coupled system is characterized by the following nonlinear Hamiltonian in the interaction picture:

$$\hat{H}\_{NL} = \frac{\chi\_{a}^{(3)}}{2} \hbar^{\dagger 2} \hbar^2 + \frac{\chi\_{b}^{(3)}}{2} \hbar^{\dagger 2} \hbar^2 = \frac{\chi\_{a}^{(3)}}{2} \hbar\_{a} (\hbar\_{a} - 1) + \frac{\chi\_{b}^{(3)}}{2} \hbar\_{b} (\hbar\_{b} - 1) \,. \tag{1}$$

where *a* and *b* label two modes of the system, *n*ˆ *<sup>a</sup>* and *n*ˆ *<sup>b</sup>* are the photon number operators in those modes. In general, both parts can have different nonlinear properties: *χ<sup>a</sup>* = *χb*.

In our considerations of the decay asymmetry we will concentrate on the nonlinear exchange between the bosonic modes. In the optical counterparts of such a system that type of interaction is realized by two-photon processes:

$$
\hat{H}\_{\text{int}\_{nml}} = \, \_\epsilon \epsilon (\mathfrak{d}^\dagger)^2 \hat{b}^2 + \epsilon \, ^\*(\hat{b}^\dagger)^2 \hat{a}^2 \,, \tag{2}
$$

where  describes the strength of the mutual interaction. Apart from the interactions between nonlinear parts, we assume that the system is driven by coherent laser fields with intensities *α* and *β*. Therefore, the appropriate Hamiltonian has the following form:

$$
\hat{H}\_{\rm ext} = a\mathfrak{d}^{\dagger} + a^\*\mathfrak{d} + \beta\mathfrak{b}^{\dagger} + \beta^\*\mathfrak{b} \,. \tag{3}
$$

Under assumptions of weak internal and external interactions (if compared to Kerr-type nonlinearity), the evolution of the discussed system is closed within a few two-mode states only, despite the continuous excitations. To make use of the nonlinear bosonic exchange between the modes of a coupled oscillatory system we prepare it in the initial excited state |Ψ(*t* = 0) = |0*a*|2*b*, and when symmetrically exciting both modes the adequate wave function takes the following form:

$$|\Psi(t)\rangle\_{\rm cut\_{nml}} = c\_{02}(t)|0\rangle\_{a}|2\rangle\_{b} + c\_{12}(t)|1\rangle\_{a}|2\rangle\_{b} + c\_{20}(t)|2\rangle\_{a}|0\rangle\_{b} + c\_{21}(t)|2\rangle\_{a}|1\rangle\_{b} \,. \tag{4}$$

and forms a 2-qutrit system, as in both modes only states |0, |1 and |2 are populated. The appropriate fidelity for obtaining the wave function (4) is plotted in Figure 1. We calculate fidelities between the two wave functions of the system: the truncated wave function for the specified subspace and the whole numerically obtained wave function in an extended basis. We can easily see that our nonlinear coupled system evolves only between the states belonging to two subspaces: {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; } (solid line) and {|1*a*|1*b*; |1*a*|2*b*; |2*a*|1*b*; |2*a*|2*b*; } (dashed line). Therefore, the influence of the states other than (4) is negligible and we can treat the state truncation as fully justified.

The analysis of formation of maximally entangled states within the system, given in [25,26], enables us to claim that the coupled oscillators form the following entangled states with probabilities equal to unity (|*B*1(2)) or slightly less than unity (|*B*<sup>3</sup> for which fidelity ≈ 0.98). As also the mode *b* is externally driven by a coherent field, additionally the state |2*a*|1*<sup>b</sup>* is included in the evolution. Consequently, also states |*B*3(4) can be obtained, but with lower probabilities. These maximally entangled states are defined by:

$$|\langle B \rangle\_{1(2)}| = \frac{1}{\sqrt{2}} \left( |2\rangle\_{a} |0\rangle\_{b} \pm i |0\rangle\_{a} |2\rangle\_{b} \right),\tag{5}$$

$$|B\rangle\_{\mathfrak{A}(4)} = \frac{1}{\sqrt{2}} \left( |2\rangle\_{\mathfrak{a}} |1\rangle\_{\mathfrak{b}} \pm i |1\rangle\_{\mathfrak{a}} |2\rangle\_{\mathfrak{b}} \right). \tag{6}$$

**Figure 1.** Fidelity between the numerically obtained wave function describing the whole two-mode system in an extended basis and the truncated wave function for the subspace {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; } — solid line, and for the subspace {|1*a*|1*b*; |1*a*|2*b*; |2*a*|1*b*; |2*a*|2*b*; } — dashed line. The initial state of the system is a two-photon state |0*a*|2*b*, *γ* = 0; *-*= *α* = *χ*/100.

#### **3. Entanglement Decay under Dissipation Channel Combinations**

The analysis of the influence of various dissipation channels on the system's dynamics can be described by standard techniques when considering master equation for the density matrix of the system. For our purpose we will assume two-sided different noisy channels. We will focus on the amplitude and phase channels applied to both of the system's modes. In particular, we will concentrate on the problem of the decay of the created entangled states for the symmetrically driven system. The amplitude damping in Born and Markov approximations for a specified mode can be obtained by adding an appropriate Liouvillian to the master equation:

$$\frac{d}{dt}\not{p}\_{\!\!\!} = \left. -i(\not{H}\not{p} - \not{p}\not{H}) + \sum\_{k=a,b} \gamma\_{k} \left[ \not{k}\not{p}\not{k}^{\sharp} - \frac{1}{2} \left( \not{p}\not{k}^{\sharp}\not{k} + \not{k}^{\sharp}\not{k} \not{p} \right) \right] \right. \tag{7}$$

For the phase-damping channel, the adequate master equation has the following form:

$$\frac{d}{dt}\boldsymbol{\hat{\rho}}\quad =\quad -i(\boldsymbol{\hat{H}}\boldsymbol{\hat{\rho}}-\boldsymbol{\rho}\boldsymbol{\hat{H}})+\frac{1}{2}\sum\_{k=a,b}\gamma\_{k}\left[2\boldsymbol{\hat{k}}^{\dagger}\boldsymbol{\hat{k}}\boldsymbol{\hat{\rho}}\boldsymbol{\hat{k}}^{\dagger}\boldsymbol{\hat{k}}-\left(\boldsymbol{\hat{k}}^{\dagger}\boldsymbol{\hat{k}}\right)^{2}\boldsymbol{\hat{\rho}}-\boldsymbol{\hat{\rho}}\left(\boldsymbol{\hat{k}}^{\dagger}\boldsymbol{\hat{k}}\right)^{2}\right],\tag{8}$$

where ˆ *k* denotes the specified mode of the system, and we additionally assume zero temperature reservoirs. The Hamiltonian *H*ˆ consists of the nonlinear part (1), the Hamiltonian (2) describing the interactions within the oscillatory system, and the Hamiltonian (3) with excitations performed in both modes. The whole coupled system is fully symmetric, as depicted in Figure 2, and for our analysis we will assume switching the location of the reservoir with specified type of noise between the qubits which we can create the entangled states and look for the symmetry in the quantum correlations dynamics.

We will focus on one type of quantum correlations, namely the entanglement. For the analysis of entanglement, we can apply negativity [33,34], which is a measure of entanglement degree in 2-qubit and qubit-qutrit systems

$$N(\rho) = \frac{1}{2} \sum\_{i} |\lambda\_i| - \lambda\_{i\prime} \tag{9}$$

where *λ<sup>i</sup>* is the *i*-th eigenvalue of a matrix *ρTk* . This matrix is obtained by performing partial transpose of the matrix *ρ* of the whole quantum system with respect to one of the subsystems, *a* or *b*. Using that measure we can identify the maximally entangled states (in the systems of 2 qubits and qubit-qutrit) as those for which the negativity equal to unity is obtained. *N*(*ρ*) = 0 holds for separable states. In further considerations we will analyze the process of entanglement creation within the 2-qubit subsystem of the whole system. Therefore we restrict ourselves to the subspace spanned by the following two-mode states: {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; }; *N*<sup>0220</sup> is the abbreviation for the negativity defined according to (9), with the matrix *ρ* describing this 2-qubit subsystem. The second 2-qubit subsystem is related to the states: {|1*a*|1*b*; |1*a*|2*b*; |2*a*|1*b*; |2*a*|2*b*; } and the appropriate negativity *N*1221.

**Figure 2.** Visualization of different location of noise in the model of the symmetric coupled oscillatory system.

It is well known that the amplitude reservoir describes the process of energy dissipation induced by the environment. Therefore, the decrease in the states' populations leads consequently to the decay of quantum correlations such as for example quantum entanglement. Phase-damping reservoir, on the other hand, describes the process of a random change in time of relative phases of the superposed states. In that process the number of photons remains unchanged and the energy is preserved. Nevertheless, the dephasing process is also responsible for quantum correlations destruction, but without changes of the total energy of the system. When looking at the density matrix—the phase-damping channel affects the off-diagonal elements only, while the amplitude damping channel influences the diagonal as well as the off-diagonal elements.

In [35] the influence of symmetrically applied amplitude and phase-damping channels was analyzed and it was shown that the dephasing channel, although it leads to the decrease of entanglement (measured by the negativity) but the rate of such a process is much slower than for the system exposed to the amplitude damping environment. Also, the entanglement sudden death, observed when energy dissipation is allowed, was not possible due to the dephasing process only. Another interesting feature of the analyzed qutrit-qubit system was also noticed. Specifically, the state |1*a*|2*<sup>b</sup>* appeared to be helpful in restoring the degree of entanglement when the system is exposed to the amplitude damping process. It is possible because of constant external driving and formation of the entangled state <sup>|</sup>*B* <sup>=</sup> <sup>√</sup><sup>1</sup> 2 |2*a*|0*<sup>b</sup>* + *i*|1*a*|2*<sup>b</sup>* , which involves |2*a*|0*b*. Therefore, the negativity can be periodically increased due to the presence of correlations, which depend on the external coherent field of the amplitude *α*. Correlations can therefore flow out of the considered 2-qubit subspace and return, due to the interactions with the additional two-mode state. Actually, our coupled nonlinear system is a higher-dimensional one. When only unitary evolution is assumed, the wave function can be expressed by (4), therefore we must consider 2-qutrits instead of just a qutrit-qubit system. Apart from the state |1*a*|2*<sup>b</sup>* there is also nonzero probability of obtaining the state |2*a*|1*b*, and other entangled states |*B*3(4) are also possible to obtain.

The state |2*a*|2*<sup>b</sup>* is a common state for both of the 2-qubit subspaces, therefore it is evident that they are not independent and correlations obtained in one of them may influence the correlations observed in the other one.

We assume that the whole system is fully symmetrical — nonlinear media are of the same type, *χ<sup>a</sup>* = *χb*, they are symmetrically driven with the same strengths *α* = *β* (see Figure 2). We apply different damping channels to both modes *a* and *b* and look for any differences under swapping the channels, which affect both parts of the whole considered system.

Still we concentrate mainly on the entanglement created only between the states {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; } (*N*0220) and between the states {|1*a*|1*b*; |1*a*|2*b*; |2*a*|1*b*; |2*a*|2*b*; } (*N*1221). We are going to ask ourselves the question, whether it is possible in a real physical system composed of these states to see any asymmetrical behavior of negativity decay by switching the damping channels between the modes *a* and *b*. The appropriate figures (Figure 3) presenting the time dependence of negativities calculated for 2-qubit subspace for various relations between the mutual (*ε*) and the external (*α*) interactions, are presented. First of all, it can be easily seen that there is an evident difference between the values of negativities obtained and the rate of decay under the symmetric channel swapping. Therefore, a fully symmetrical quantum system can exhibit a clear asymmetry in evolution of quantum correlations when applying and swapping different types of reservoirs.

As seen in Figure 3, for all the considered cases larger values of entanglement characterized by *N*<sup>0220</sup> and *N*1221, are possible, when the mode *a* is exposed to the amplitude damping channel and the mode *b* decays to the phase reservoir (Figure 3a,c), than when the reservoirs are swapped between the modes (Figure 3b,d). There is also a noticeable difference in the time span for which the system can be considered to be an entangled one. Again, longer times for disentangling of the qubits are necessary if the mode *a* decays to the reservoir with amplitude damping and the mode *b* to the phase than for the case with the mode *a* decaying to the dephasing channel and the mode *b* to the amplitude channel. Therefore, there is a certain asymmetry in the response of the coupled system to the exchange of the damping channels between the parts of the whole system. We can also see that the presence of correlations is connected mainly with the subspace {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; }. In all the considered cases (in various damping channels configurations) the values of negativities *N*<sup>0220</sup> are noticeably higher than appropriate negativities *N*<sup>1221</sup> for the second considered subspace. Additionally, the decrease in the value of *N*<sup>0220</sup> is obviously connected with the increase of *N*1221. Thus, there is some kind of flow of entanglement between the considered subspaces of the whole system.

Next we will address the problem of the specific entangled states (5)–(6) formation to decide, which of them are responsible for nonzero negativities and which influence the observed asymmetry in *N*0220(*t*) and *N*1221(*t*) dependencies after exchanging the reservoirs between the modes *a* and *b*.

It is evident that all the states (5)–(6) should be possible to observe during the time evolution of the system. States |*B*1(2) are created due to the 2-photon interactions, states |*B*3(4) are obtained as a consequence of additional external excitation, which is responsible for populating the states |1*a*|2*<sup>b</sup>* and |2*a*|1*b*.

In Figure 4 the evolutions of appropriate entangled states are plotted. When the mode *a* is decaying into the amplitude damping channel, the states |*B*1(2) oscillate with almost the same frequencies, simultaneously obtaining maxima and minima. When the mode *a* is exposed to the dephasing

process, the evolution of |*B*<sup>2</sup> state is slightly perturbed by another frequency. From the behavior of entangled states presented, it comes out that a complete decay of negativity is related to the higher probability of obtaining the states |*B*3(4) —in such a case the correlations initially created in the subspace {|0*a*|0*b*; |0*a*|2*b*; |2*a*|0*b*; |2*a*|2*b*; } are transferred out of the considered subspace. Due to the presence of one common state belonging to both subspaces of the whole system, correlations can be significantly transferred to the subspace {|1*a*|1*b*; |1*a*|2*b*; |2*a*|1*b*; |2*a*|2*b*; }. It happens in shorter time for the situation when the mode *a* is exposed to the dephasing process, as it is easier to populate the state |1*a*|2*<sup>b</sup>* by an external field, while the mode *a* does not lose its population. The presence of the state |1*a*|2*<sup>b</sup>* is therefore crucial in obtaining the asymmetry in correlations, measured by the negativity decay introduced by swapping the reservoirs between the qubits forming the entangled states. As seen from Figure 4, the evolution of the created entangled states strongly depends on the type of reservoir applied to each of the modes. Such a difference is obviously seen in the values of negativities for the qubit subspaces and results in the asymmetry in correlations decay.

**Figure 3.** Negativities for 2-qubit subspaces—*N*<sup>0220</sup> (solid line) and *N*<sup>1221</sup> (dashed line) versus time scaled in 1/*χ* units for combinations of damping channels. In (**a**) and (**c**) the mode *a* is exposed to the amplitude damping channel and mode *b* into the phase-damping channel. In (**b**) and (**d**) the channels are swapped. The initial state of the system is a two-photon one |0*a*|2*b*, *γ* = *-* = *α* = *χ*/100 for (**a**) and (**b**); *γ* = *-*= *χ*/100 and *α* = *χ*/25 for (**c**) and (**d**).

It appears that the state |1*a*|2*<sup>b</sup>* plays a crucial role in formation and behavior of the entanglement in the coupled system. Its presence influences significantly the time of disentanglement and indirectly the values of *N*<sup>0220</sup> — higher probabilities of |*B*3(4) decrease the probabilities of |*B*1(2). When the mode *a* of the state |1*a*|2*<sup>b</sup>* decays by the amplitude channel to the state |0*a*|2*b*, it increases the probability of generating the states |*B*1(2) and results in higher values of *N*0220. On the other hand, the mode *a* is externally pumped by the linear interaction and again the state |1*a*|2*<sup>b</sup>* can be formed.

**Figure 4.** Fidelities for obtaining the entangled states for *γ* = *-* = *χ*/100, *α* = *χ*/25 and the initial state |0*a*|2*b*. In (**a**) the mode *a* is exposed to the amplitude damping channel and the mode *b* to the phase-damping channel. In (**b**) the channels are swapped. |*B*1—solid line; |*B*2—dashed line; |*B*3—dashed-dotted line and |*B*4—dotted line.

#### **4. Conclusions**

We discussed a coupled nonlinear system whose parts are driven by coherent laser fields. We additionally assumed that interactions between the nonlinear media are performed via two-boson exchange. Such systems, when appropriate conditions describing relations between the strengths of the interactions and the nonlinearity are fulfilled, evolve only between a limited number of states despite the constant energy supply. We have shown that the model we are dealing with can be treated as a 2-qutrit system and the unitary evolution can be described by a truncated wave function. As the nonlinearity describing the interacting media must be large (as compared to other interactions), in that sense the possibility of the Hilbert space truncation may be regarded as the conventional photon blockade. We have discussed the problem of the entanglement decay under the influence of different noisy channels. Both parts of the considered coupler (the mode *a* and the mode *b*) were exposed to the local amplitude or phase-damping environments. Our goal was to show that despite the symmetry of our system, it is possible to differentiate the decay rates and the values of the entanglement obtained within 2-qubits subspaces by appropriately choosing the type of system-environment interactions. For that purpose, we have analyzed the time evolution of negativities after first applying and then swapping different noisy channels to the modes *a* and *b*. We have identified the asymmetry in the response of the system to the interactions with the reservoirs. There is a noticeable difference in the rate at which entanglement decreases, due to the location of the noise between the qubits.

**Author Contributions:** Conceptualization and Formal analysis, A.K.-K.; Software, Validation and Investigation, A.K.-K. and G.C.; Visualization, G.C.; Writing—original draft, A.K.-K.; Writing—review & editing, G.C.

**Funding:** This research received no external funding

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**




c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
