Quantum and Topological Dynamics of the GKSL Equation in the Camel-Like Framework

1. Introduction

In order to analyze the interactions between two or more elements in a quantum information system, one can represent these elements or events as states of information, termed infons. These infons might describe, for example, the interaction between two bosons [1], the exchange of cognitive signals between neurons in the brain [2], or the quantum tunneling of electron pairs across a potential barrier [3]. The dynamics of these interacting elements can be effectively modeled within the framework of open quantum systems [1,4]. This formalism provides a fully quantum-mechanical description of a system $\mathcal{S}$ as it evolves under the influence of a local environment $ϵ$ [5,6,7]. The environment $ϵ$ accounts for effects such as decoherence and dephasing, which impact the state of the system $\mathcal{S}$.

The system-environment relationship can be represented hierarchically as:

$$\mathcal{S}\subset ϵ\subset \mathcal{U},$$

(1)

where $\mathcal{U}$ represents the highest containing state (e.g., an isolated entangled system subjected to an electromagnetic field).

The evolution of the system $\mathcal{S}$ under the influence of $ϵ$ can be described using the density matrix formalism. Let ${\rho }_{\mathrm{initial}}$ denote the initial state of $\mathcal{S}$. Under the action of a Hamiltonian $\widehat{H}$ (representing a local perturbation within $ϵ$), the system evolves to a final state ${\rho }_{\mathrm{final}}$:

$${\rho }_{\mathrm{initial}}\stackrel{\widehat{H}}{\to }{\rho }_{\mathrm{final}}.$$

(2)

The environment $ϵ$ can be an electromagnetic field for particles [8,9,10,11] or even a cognitive environment for neuronal cells, including electric currents, neurotransmitters, and other stimuli [7,12,13].

1.1. Quantum-like Formalism for the GKSL Equation

The foundational framework of open quantum systems can be extended beyond traditional physical scenarios to model quantum-like processes in biological and cognitive systems. Such models have been increasingly applied to represent quantum phenomena across natural sciences [7,14,15].

For instance, neuronal signals (treated as infons) can be described by quantum observables and probabilities, where information processing is governed by specific dynamics between decoherence and non-unitary evolution of quantum states. The transition of a neuronal state under environmental perturbations can be represented as a superposition of wavefunctions [7], where the wavefunction of the neuron exposed to some perturbation undergoes the following transition:

$${\psi }_{ss}\stackrel{\widehat{H}}{\to }{C}_0{\psi }_{ex}+C_1{\psi }_{qs},$$

(3)

where ${\psi }_{ss}$ denotes the wavefunction at steady state with the environment, ${\psi }_{ex}$ the excited state, and ${\psi }_{qs}$ the quiescent state. Here, $C_0$ and $C_1$ are complex coefficients satisfying the normalization condition $|C_0{|}^2+{\left|C_1\right|}^2=1$. $\widehat{H}$ is a locally acting perturbation within $ϵ$, and $ϵ$ represents the totality of electrochemical signals in the environment of the neuron, including electric currents, EMF, and neurotransmitters generated by other neurons [12].

In line with quantum theory, the transition induced by $\widehat{H}$ in (3) can be represented within the density matrix formalism:

$${\rho }_{\mathrm{initial}}\stackrel{\widehat{H}}{\to }{\rho }_{\mathrm{final}},$$

(4)

where ${\rho }_{\mathrm{initial}}$ transitions to ${\rho }_{\mathrm{final}}$, describing the probabilistic nature of superposition. Now, in order to construct the density matrix from a state vector we define $\mathcal{H}$ as a complex Hilbert space, where a pure quantum state is represented by a unit vector $|\psi \rangle \in \mathcal{H}$, and its dual vector $\langle \psi |\in {\mathcal{H}}^{*}$ by:

$$\langle \psi |:|\varphi \rangle \mapsto \langle \psi |\varphi \rangle \in \mathbb{C}.$$

(5)

Then, the fhe final density matrix ${\rho }_{\mathrm{final}}$ corresponds to the pure state described by the outer product:

$${\rho }_{\mathrm{final}}=|{\psi }_{\mathrm{final}}\rangle \langle {\psi }_{\mathrm{final}}|,$$

(6)

with the pure state ${\psi }_{\mathrm{final}}$ given by:

$$|{\psi }_{\mathrm{final}}\rangle =C_0|0\rangle +C_1|1\rangle ,$$

(7)

and

$$\langle {\psi }_{\mathrm{final}}|=C_0^{*}\langle 0|+C_1^{*}\langle 1|,$$

(8)

where $C_0$ and $C_1$ are complex coefficients, and $C_0^{*}$ and $C_1^{*}$ are their complex conjugates, respectively.

By examining equation (3), we can see that ${\psi }_{ex}$ and ${\psi }_{qs}$ are represented by ${\rho }_{final}$ in (6), with arbitrary coefficients in both real and complex form. Thus, we have demonstrated how the biological state transition in (3) can be represented by a superposition of quantum states in an open quantum system via (6). This formulation automatically satisfies the Born rule, which gives the probability of measuring the system:

$$P\left(\alpha \right)=|\langle \alpha |{\psi }_{\mathrm{final}}\rangle {|}^2.$$

(9)

The time-dependent interaction of such a measurable open quantum system with a specific environment can be studied using the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation [16,17,18]. This equation describes how the system stabilizes to a steady state as time approaches infinity [2]:

$$\underset{t\to \infty }{lim}\widehat{\rho }\left(t\right)={\widehat{\rho }}_{steady}$$

(10)

1.2. GKSL Operator Structure and Dynamics

Let us first introduce the space of states, which we devise using the finite-dimensional Hilbert space $\mathcal{H}$ endowed with the scalar product $\langle ·,·\rangle$. Let the space of density operators (which represent the states of the system) be denoted by $\mathcal{D}\left(\mathcal{H}\right)$, and the space of all linear operators be given by $\mathcal{L}\left(\mathcal{H}\right)$. Furthermore, $\mathcal{L}\left(\mathcal{H}\right)$ is a complex Hilbert space with the scalar product $\langle A|B\rangle =\mathrm{Tr}\left(A^{†}B\right)$, with A and B being operators. We have furthermore that the environment (i.e. dephasing/decoherence) is represented by superoperators in $\mathcal{L}\left(\mathcal{H}\right)$ that act on the density matrix $\rho$. The GKSL equation is the time-dependent master equation for open quantum systems, developed in [17,19] from the Schrödinger equation. It is described by

$${\displaystyle \frac{\partial \rho \left(t\right)}{\partial t}}=-i[H,\rho \left(t\right)]+\mathcal{L}\left[\rho \left(t\right)\right].$$

(11)

Here, $\rho \left(t\right)$ is the unknown probability density matrix, which evolves with time and is used to resolve properties of an open quantum system, such as entropy, entanglement, and decoherence [17,19]. The first term on the RHS in (11), the Hamiltonian, is the unitary and deterministic part of the evolution, while the second term on the RHS in (11), the dissipation operator, represents the non-unitary effects on the density matrix that ultimately represent non-deterministic evolution of the system. The Hamiltonian represents the internal dynamics of the system at hand, and are often given by the Pauli matrices, however in some cases it can also include component related to the interaction with the environment [20]. The dissipation operator, $\mathcal{L}$, is the GKSL superoperator given by:

$$\mathcal{L}\left[F\left(\rho \right)\right]=\sum _k\left(L_k\rho L_k^{†}-{\displaystyle \frac{1}{2}}\{L_k^{†}{L}_k,\rho \}\right),$$

(12)

where $L^{†}$ is the adjoint of the jump operator L, defined by the relation:

$$\langle \psi ,L\varphi \rangle =\langle L^{†}\psi ,\varphi \rangle .$$

(13)

This relation holds for any linear operator L, Hermitian or not, and defines $L^{†}$ with respect to the inner product. In general, the jump operators L in the Lindblad formalism are non-Hermitian, representing dissipative processes such as decoherence or decay. We have hence the GKSL equation given by

$${\displaystyle \frac{d\rho }{dt}}=-i[H,\rho ]+\gamma \sum _k\left(L_k\rho L_k^{†}-{\displaystyle \frac{1}{2}}\{L_k^{†}{L}_k,\rho \}\right),$$

(14)

where $\gamma$ is the relaxation factor for decoherence processes to the dissipation operator. Each jump operator represents the special type of interaction of the system $\mathcal{S}$ by the density matrix $\rho$ in (10), with the environment $ϵ$ (described by $\mathcal{L})$, and vary according to the decoherence and dissipation resulting by the formalism of H and $\mathcal{L}$. This forms the main difference between the original Schrödingers equation and the GKSL equation: the GKSL equation produces stabilization of its solution to the steady-state (10).

1.3. Camel-Like Framework in Entropy Evolution of Open Quantum Systems

In the study of complex systems, entropy serves as a fundamental measure of disorder and information content. We introduce the concepts of N-hump camel-like entropy to describe distinct temporal profiles of entropy evolution in light of the studies by Khrennikov [7,12,13,15,21]. A 1-hump camel-like entropy is characterized by a single local maximum and minimum, followed by damped oscillations, while a 2-hump camel-like entropy exhibits two such maxima and minima before stabilizing. These definitions are inspired by the analogy of a camel’s humps, providing an intuitive visualization of entropy dynamics. In quantum information research, these profiles could prove invaluable for analyzing the behavior of quantum systems undergoing decoherence, entanglement dynamics, or information scrambling. For instance, a 2-hump entropy profile might indicate the presence of intermediate metastable states in a quantum process, while a 1-hump profile could describe simpler, single-phase transitions. Furthermore, we can also obtain N-humps and obtain Rabi-like oscillations [22] in the entropy evolution which have similarities to the study of spin-state fractions in single atom system [23]. By categorizing entropy evolution in this way, we may gain a deeper understanding of the underlying mechanisms governing quantum systems, paving the way for improved control and optimization in quantum computing and communication protocols using the hereafter defined camel-like framework presented by Khrennikov [2,7,12,13,15,21].

Definition 1.

(Camel-like behaviour of entropy)

One-Hump Camel-Like behaviour of the Entropy

Consider the following time intervals:

$$I_1=[0,t_{max}),I_2=(t_{max},t_{min}),I_3=(t_{min},+\infty )$$

We characterize these intervals using the derivative of entropy $E^{\prime }\left(t\right)$ as follows:

• $E^{\prime }\left(t\right)>0$ on $I_1$, and $E^{\prime }\left(t_{max}\right)=0$
• $E^{\prime }\left(t\right)<0$ on $I_2$, and $E^{\prime }\left(t_{min}\right)=0$
• On $I_3$, $E\left(t\right)$ exhibits damped oscillations in amplitude.

This behavior corresponds to a one-hump camel analogy, where the entropy profile has a single local maximum at $t_{max}$ and a single local minimum at $t_{min}$.

Two-Hump Camel-Like behaviour of the Entropy

Consider the following time intervals:

$$I_1=[0,t_{{max}_1}),I_2=(t_{{max}_1},t_{{min}_1}),I_3=(t_{{min}_1},t_{{max}_2}),I_4=(t_{{max}_2},t_{{min}_2}),I_5=(t_{{min}_2},+\infty )$$

We characterize these intervals using the derivative of entropy $E^{\prime }\left(t\right)$ as follows:

• $E^{\prime }\left(t\right)>0$ on $I_1$, and $E^{\prime }\left(t_{{max}_1}\right)=0$
• $E^{\prime }\left(t\right)<0$ on $I_2$, and $E^{\prime }\left(t_{{min}_1}\right)=0$
• $E^{\prime }\left(t\right)>0$ on $I_3$, and $E^{\prime }\left(t_{{max}_2}\right)=0$
• $E^{\prime }\left(t\right)<0$ on $I_4$, and $E^{\prime }\left(t_{{min}_2}\right)=0$
• On $I_5$, $E\left(t\right)$ exhibits damped oscillations in amplitude.

This behavior corresponds to a two-hump camel analogy, where the entropy profile has two local maxima at $t_{{max}_1}$ and $t_{{max}_2}$, and two local minima at $t_{{min}_1}$ and $t_{{min}_2}$.

N-hump camel-like behaviour of the Entropy

Let $E\left(t\right)$ be the entropy of a quantum system evolving under the GKSL equation. We define $E\left(t\right)$ as having Nhumps (local maxima) if there exists a sequence of critical points $\{t_{{\max }_1},t_{{\min }_1},\dots ,t_{{\max }_N},t_{{\min }_N}\}$ such that the time domain ${\mathbb{R}}^{+}$ is partitioned into $2N+1$ intervals:

$$\begin{array}{cc}\hfill I_0& =[0,t_{{\max }_1}),\hfill \\ \hfill I_{2k-1}& =(t_{{\min }_{k-1}},t_{{\max }_k})\mathrm{for}k=1,\dots ,N,\hfill \\ \hfill I_{2k}& =(t_{{\max }_k},t_{{\min }_k})\mathrm{for}k=1,\dots ,N,\hfill \\ \hfill I_{2N+1}& =(t_{{\min }_N},+\infty ),\hfill \end{array}$$

where $t_{{\min }_0}\equiv 0$, and the entropy derivative $E^{\prime }\left(t\right)$ satisfies:

• Growth Phases: $E^{\prime }\left(t\right)>0$ on $I_{2k-1}$ for all k, with $E^{\prime }\left(t_{{\max }_k}\right)=0$ (local maxima).
• Decay Phases: $E^{\prime }\left(t\right)<0$ on $I_{2k}$ for all k, with $E^{\prime }\left(t_{{\min }_k}\right)=0$ (local minima).
• Asymptotic Damping: On $I_{2N+1}$, $E\left(t\right)$ exhibits underdamped oscillations converging to a steady-state value $E_{\mathrm{ss}}$.

From this definition, we see that the defined camel-like entropy profiles, characterized by transient maxima and minima in disorder, evoke parallels to barrier-crossing phenomena across diverse systems, based on their number of humps. For instance, in enzyme catalysis [24], the one-hump entropy profile mirrors the transient disordered transition state (hump) that precedes product stabilization—akin to overcoming a single entropic barrier during catalysis. In contrast, two-hump entropy aligns with multi-phase processes like protein folding [25], where disordered intermediates and misfolded states generate sequential entropy peaks, or quantum tunneling [26] with successive decoherence events. Moreover, spin-fraction measurement for single atom systems [23] show Rabi-oscillations, which are in a damped form here described as N-hump camel like behaviour.

1.4. Operators

Now that we have formalized the definition of the camel-like entropy, we can define the Hamiltonian H and the dissipation operator $\mathcal{L}$ that governs the dynamics of the system following the camel-like formalism proposed by Khrennikov [2,7,12,13,15,21]. By he GKSL equation:

$${\rho }^{\prime }\left(t\right)=-i[H,\rho \left(t\right)]+\gamma \left(C\rho \left(t\right)C^{†}-{\displaystyle \frac{1}{2}}\{C^{†}C,\rho \left(t\right)\}\right),$$

(15)

we have, $\gamma$ which is the relaxation rate representing the strength of interaction between a quantum system and its environment. We present the model example of the order-stable dynamics that matches the biological order-preserving behavior and select the Hamiltonian H and the interaction operator C as follows [2]:

$$H={\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),C=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right).$$

(16)

We calculate these operators into their four-dimensional form:

$$H={\sigma }_x\otimes {\sigma }_x=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\end{array}\right)$$

(17)

The dissipation operator redimensioned to $N=4$ gives the collapse operator C, which we extend to act on the two-qubit system. By choosing an operator that acts locally on one qubit, specifically $C={\sigma }_{+}\otimes I$, where ${\sigma }_{+}$ is the raising operator for the first qubit and I is the identity matrix for the second qubit, we obtain:

$$C={\sigma }_{+}\otimes I=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right).$$

(18)

Based on analysis done in [27] we set the relaxation factor $\gamma ={\displaystyle \frac{1}{2}}$.

Having set the premises for the study of the camel-like framework of the GKSL equation applied on a set of quantum states, we now define our aims. The first aim is to study the evolution of the the von Neumann entropy for the set of selected quantum states, and the quantum correlations of quantum discord, entanglement entropy, mutual information and the Bell-CHSH parameter, so that we show the performance of the camel-like framework on pure states, product states and mixed states. Second, we wish to analyse the properties of the camel-like framework by using a gradient field and topological basin analysis, so that we can study the evolution of solutions on the Bloch sphere and their stability. At last, we provide analytical solutions which satisfy the definition of the 1. Using this combination of analytical tools, we thus aim to present how the camel-like framework works and to which physical experiment it has strong similarities to.

2. Methods of Analysis

2.1. State Representation

We can form the given (normalized) states for the three classes of quantum states, pure states, product states and mixed states: The states of interest are represented thus by the three categories:

Bell-states.

We include two Bell states in our calculations, namely ${\varphi }^{+}$ and ${\varphi }^{-}$.

Product states.

$$|{\psi }_{\mathrm{prod}}\rangle =\sqrt{{\displaystyle \frac{3}{5}}}|00\rangle +\sqrt{{\displaystyle \frac{3}{20}}}|01\rangle +\sqrt{{\displaystyle \frac{1}{5}}}|10\rangle +\sqrt{{\displaystyle \frac{1}{20}}}|11\rangle .$$

(19)

In tensor product form:

$$|{\psi }_{\mathrm{prod}}\rangle =\left(\begin{array}{cc}\sqrt{{\displaystyle \frac{3}{5}}}& {\displaystyle \frac{\sqrt{\frac{3}{5}}}{2}}\\ {\displaystyle \frac{1}{\sqrt{5}}}& {\displaystyle \frac{1}{2\sqrt{5}}}\end{array}\right)$$

(20)

$$|{\psi }_{\mathrm{prod}}^{\prime }\rangle ={\displaystyle \frac{\sqrt{6}}{6}}|00\rangle +{\displaystyle \frac{\sqrt{6}}{6}}|01\rangle +{\displaystyle \frac{\sqrt{3}}{3}}|10\rangle +{\displaystyle \frac{\sqrt{3}}{3}}|11\rangle ,$$

(21)

in tensor product form:

$$|{\psi }_{\mathrm{prod}}^{\prime }\rangle =\left(\begin{array}{cc}{\displaystyle \frac{1}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{6}}}\\ {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{3}}}\end{array}\right)$$

(22)

Mixed states.

The mixed states are defined with probabilities $p=q=0.5$ using the pure states:

$$\begin{array}{cc}\hfill {\rho }_1& =p{\rho }_{{\psi }_3}+(1-p){\rho }_{{\psi }_4}\hfill \\ \hfill & =0.5|{\psi }_3\rangle \langle {\psi }_3|+0.5|{\psi }_4\rangle \langle {\psi }_4|\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\rho }_2& =q{\rho }_{{\psi }_5}+(1-q){\rho }_{{\psi }_6}\hfill \\ \hfill & =0.5|{\psi }_5\rangle \langle {\psi }_5|+0.5|{\psi }_6\rangle \langle {\psi }_6|\hfill \end{array}$$

(24)

where:

$$\begin{array}{cc}\hfill |{\psi }_3\rangle & ={\displaystyle \frac{1}{\sqrt{2}}}|00\rangle +{\displaystyle \frac{1}{\sqrt{2}}}|01\rangle ,\hfill \\ \hfill {\rho }_{{\psi }_3}& =|{\psi }_3\rangle \langle {\psi }_3|\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill |{\psi }_4\rangle & =|10\rangle ,\hfill \\ \hfill {\rho }_{{\psi }_4}& =|{\psi }_4\rangle \langle {\psi }_4|\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill |{\psi }_5\rangle & =|11\rangle ,\hfill \\ \hfill {\rho }_{{\psi }_5}& =|{\psi }_5\rangle \langle {\psi }_5|\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill |{\psi }_6\rangle & ={\displaystyle \frac{\sqrt{3}}{2}}|00\rangle +{\displaystyle \frac{1}{2}}|10\rangle ,\hfill \\ \hfill {\rho }_{{\psi }_6}& =|{\psi }_6\rangle \langle {\psi }_6|\hfill \end{array}$$

(28)

2.2. Numerical Calculations

We employ numerical methods to solve the Lindblad equation, specifically utilizing matrix exponential techniques where we use the forward-Euler method [28] for solving the Lindblad equation to calculate purity, entropy, quantum discord and mutual information. The calculation of von Neumann entropy uses logarithm base 2 to align with information theory, where entropy measures uncertainty in bits. When we calculate the entanglement entropy, we use the mesolve method as part of the qutip package [29]. Finally, for the calculations of the CHSH parameter,, we use the Runge-Kurra method of 4th order. These approaches allows us to study the evolution of the density matrix through the matrix exponential. We analyze the numerical solutions for the density matrices over time to observe the evolution of the mentioned quantum relations as a result of the Lindblad dynamics. All simulations are carried out with Python.

2.3. Quantum Correlations Computational Code

We form Python codes for the study of the time evolution of the von Neumann entropy for the considered pair of states (Bell states, product states and mixed states). All states are represented as density matrices, computed as the outer products of the states with their conjugates. The Hamiltonian governing the system is defined as $H={\sigma }_x\otimes {\sigma }_x$, where ${\sigma }_x$ is the Pauli-X matrix. The collapse operator is defined as $\mathcal{C}={\sigma }_{+}\otimes I$, with ${\sigma }_{+}$ being the raising operator for a single qubit and I the identity matrix. The evolution of the system is modeled using the Lindblad equation, which accounts for both the Hamiltonian and dissipative processes. The Lindblad evolution step is implemented by iterating over small time steps $\Delta t$, updating the density matrix at each step, and normalizing the result to ensure it remains a valid density matrix. The von Neumann entropy is computed at each time step using the eigenvalues of the density matrix. To avoid numerical issues, eigenvalues are clipped to a small value before taking their logarithm. The time evolution is simulated for 1000 steps, with a time step of $\Delta t=0.01$.

2.3.1. von Neumann and Entanglement Entropy

The entanglement entropy is computed using the von Neumann entropy, which quantifies the amount of quantum entanglement in a system. In the code, we compute this as:

$$S\left(\rho \right)=-\mathrm{Tr}\left(\rho {log}_2\rho \right)$$

where the logarithm is taken in base 2, and the trace operation sums over the eigenvalues of the matrix. To calculate the entanglement entropy for the two-qubit systems, the following procedure is applied in the Python code:

1. Partial Trace: The first step is to compute the partial trace over one of the qubits (in this case, the second qubit) to obtain the reduced density matrix for the first qubit. The partial trace operation effectively “traces out” the second qubit, leaving a reduced matrix that describes the subsystem of interest. For a two-qubit system described by the density matrix $\rho$, the partial trace over the second qubit is computed as:

$${\rho }_{\mathrm{red}}={\mathrm{Tr}}_2\left(\rho \right)$$

For the two-qubit Bell states, this results in a 2x2 matrix representing the reduced density matrix of the first qubit.

2. Eigenvalue Calculation: After obtaining the reduced density matrix, its eigenvalues are computed. These eigenvalues represent the probabilities associated with the possible states of the reduced system.

3. Von Neumann Entropy: Once the eigenvalues of the reduced density matrix are obtained, the von Neumann entropy is computed by applying the formula:

$$S\left({\rho }_{\mathrm{red}}\right)=-\sum _i{\lambda }_i{log}_2\left({\lambda }_i\right)$$

where ${\lambda }_i$ are the eigenvalues of the reduced density matrix ${\rho }_{\mathrm{red}}$. To prevent numerical issues, the eigenvalues are clipped to a small value before taking the logarithm (e.g., ${\lambda }_i\ge {10}^{-14}$) to avoid computing ${log}_2\left(0\right)$. This process is repeated for each time step during the evolution, allowing us to track how the entanglement entropy changes over time due to the effects of the dissipative dynamics introduced by the Lindblad evolution.

2.3.2. Mutual Information

The mutual information is computed in the Python script as:

$$I(A:B)=S\left(A\right)+S\left(B\right)-S\left(AB\right)$$

where $S\left(A\right)$, $S\left(B\right)$, and $S\left(AB\right)$ are the von Neumann entropies of the reduced density matrices ${\rho }_A$, ${\rho }_B$, and the total density matrix $\rho$, respectively. The mutual information is calculated at each time step for both states by tracing out subsystems and computing the von Neumann entropies. The results are stored and plotted as a function of time. This process tracks how the mutual information between the subsystems evolves over time, showing the dynamics of quantum correlations in the Bell states.

2.3.3. Derivative of the Entanglement Entropy

The first derivative of the entanglement entropy with respect to time ${\displaystyle \frac{d{S}_A}{dt}}$ is computed numerically using the ‘np.gradient’ function on the entanglement entropy values, which are computed as in the entropy-code.

2.3.4. Quantum Discord

Quantum discord is computed as the difference between mutual information and classical correlation. The quantum discord is computed at each time step for both states. The function $‘quantum-discord’$ function calculates this using the defined mutual information and classical correlation functions. The quantum discord at $t=0$ is computed and printed for both states, providing insight into the initial correlations. The eigenvalues of the initial density matrices are computed to analyze the purity and state properties of the system at the start.

2.3.5. CHSH Parameter

In the Python code used to calculate the time evolution of the Clauser-Horne-Shimony-Holt (CHSH) parameter for both states, we calculate it by evaluating the expectation values of four measurement settings, $E\left(AB\right)$, $E\left(A^{\prime }B\right)$, $E\left(A{B}^{\prime }\right)$, and $E\left(A^{\prime }{B}^{\prime }\right)$, for the two qubit system at each timestep. The code then computes the CHSH parameter over time and plots its evolution, comparing it to the classical bound of 2 and the quantum maximum of $2\sqrt{2}$. This allows for the analysis of quantum correlations and the violation of local realism as a function of time for both states. All calculations are carried out using the following Python packages:

• numpy: A fundamental package for numerical computing in Python, used for array manipulation and linear algebra operations such as matrix products and Kronecker products.
• scipy.linalg: Specifically, the expm function from this module is used for matrix exponentiation.
• matplotlib.pyplot: A plotting library used for generating visualizations, such as graphs for quantum discord and CHSH parameters.
• qutip: A quantum toolbox used for simulating quantum systems. Key functions from qutip include:

  −

  sigmax(), sigmay(), sigmaz(): Pauli matrices for qubits.

  −

  qeye(): Identity matrix for qubits.

  −

  basis(): Creates the computational basis states.

  −

  tensor(): Computes tensor products of quantum states and operators.

  −

  mesolve(): Solves the Schrödinger equation with dissipation (Lindblad evolution).

  −

  entropy_vn(): Computes the von Neumann entropy.

2.4. Analytical Solutions

The analytical solutions to the GKSL equation under the camel-like framework were found using a Python script that uses SymPy to model and solve the time evolution of a two-level quantum system described by a density matrix. The time variable t is defined as a real-valued symbol, and the matrix components ${\rho }_{11}\left(t\right)$, ${\rho }_{12}\left(t\right)$, and ${\rho }_{21}\left(t\right)$ are declared as time-dependent functions. These represent, respectively, the population of the first state and the coherences of the system. The system of differential equations encodes the dynamics of the density matrix, possibly derived from a Lindblad-type master equation or an effective non-Hermitian Hamiltonian. Specifically, ${\dot{\rho }}_{11}$ evolves according to a population exchange term and a term involving the imaginary part of the coherences, while ${\dot{\rho }}_{12}$ and ${\dot{\rho }}_{21}$ evolve under damping and driving terms that couple to the population difference. Initial conditions are provided in the form of a dictionary, specifying the values of the components of $\rho \left(t\right)$ at time $t=0$; for example, setting ${\rho }_{11}\left(0\right)=0$, ${\rho }_{12}\left(0\right)=0$, and ${\rho }_{21}\left(0\right)=0$ corresponds to the system being initially in the excited state, that is, $\rho \left(0\right)=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)$. The script then attempts to solve the system symbolically using dsolve, and generates the solution in the command prompt. The program also checks the hermiticity of the resulting matrix by verifying whether ${\rho }_{12}={\rho }_{21}^{*}$, whether ${\rho }_{11}$ remains real-valued, and whether the trace is preserved. Finally, the steady-state limit ${\rho }_{\infty }={lim}_{t\to \infty }\rho \left(t\right)$ is computed and formatted in the command prompt.

2.5. Gradient Map Analysis on the ${\mathbb{S}}^2$ Manifold Represented by the Bloch Sphere

The gradient flow analysis was implemented in MATLAB through a custom script that visualizes critical dynamics on the Bloch sphere representation of quantum states. The script defines the mathematical components $g_1\left(\theta \right)=(cos\theta +\sqrt{2}/2)(cos\theta +1)$ and $g_2(\theta ,\varphi )={sin}^2\varphi {sin}^2(\varphi -\pi /2)(cos\theta +1)$ from (70), along with the squared magnitude function $h(\theta ,\varphi )=g_1^2+g_2^2$ whose gradient drives the flow analysis. Using numerical differentiation with a delta of ${10}^{-6}$, the gradient is computed and incorporated into a differential equation system modeling the negative gradient flow, solved via MATLAB’s ode45 Runge-Kutta solver with high precision tolerances (${10}^{-6}$). The Bloch sphere is rendered as a surface mesh colored by the z-coordinate, with the southern quarter-circle at $\theta =3\pi /4$ explicitly highlighted in red as the critical manifold. Six unstable initial states are marked with distinct colors and enlarged markers, representing the singularities in (53) where the map $f:{\mathbb{S}}^2\to {\mathbb{R}}^2$ is undefined (two of the six states overlap on the South Pole - hence yielding 5 visible states as points on the Bloch sphere). White gradient flow vectors are superimposed on the sphere using normalized arrows, and cyan trajectories from the unstable points visualize convergence behavior toward stable attractors. The analysis includes seven viewing angles, including a dedicated south pole view, along with temporal evolution plots of the h-function to demonstrate the system’s dynamical stability and convergence properties around the critical southern ring where states share a common instability towards decoherence.

2.6. Topological Basin Analysis on the ${\mathbb{S}}^2$ Manifold Represented by the Bloch Sphere

The MATLAB script for analyzing the topological dynamics of the GKSL equation simulates gradient flows on a spherical coordinate grid $(\theta ,\phi )$ across the Bloch sphere, with the singularities in the set (53) corresponding to the zeros of the functions $g_1$ and $g_2$. The system evolves according to the negative gradient of the scalar function $h(\theta ,\phi )=g_1{\left(\theta \right)}^2+g_2{(\theta ,\phi )}^2$, where $g_1\left(\theta \right)=(cos\theta +\sqrt{2}/2)(cos\theta +1)$ and $g_2(\theta ,\phi )={sin}^2\phi {sin}^2(\phi -\pi /2)(cos\theta +1)$. The singularities occur where both $g_1=0$ and $g_2=0$: at the north pole ($\theta =0$), the south pole ($\theta =\pi$), and along the critical latitude $\theta =3\pi /4$ (where $cos\theta =-\sqrt{2}/2$). The gradient flow equations incorporate regularization to handle the coordinate singularity at the poles. Trajectories are integrated via MATLAB’s ode45 solver from initial conditions sampled on a $30\times 60$ grid, and classified into three basins: (1) the northern hemisphere ($\theta <3\pi /4$), (2) the south pole region ($3\pi /4\le \theta \le \pi$), and (0) elsewhere. The critical boundary at $\theta =3\pi /4$ acts as a separatrix between basins. Visualizations in both 2D spherical coordinates and 3D Cartesian projections explicitly mark these singularities and basin boundaries, with multiple viewing angles emphasizing the southern hemisphere structure to reveal the system’s global topological landscape.

3. Results and Discussion

3.1. Entropy for Bell-States, Product States and Mixed States

The numerical calculations allow us to calculate the most fundamental property of the system, defined by the von Neumann entropy $S\left(\rho \right)$. The von Neumann entropy (also defined as quantum entropy) quantifies the uncertainty of a quantum state represented by the density operator $\rho$, defined as:

$$S\left(\rho \right)=-\mathrm{Tr}(\rho log\rho ).$$

(29)

Following Figure 1, both Bell states $|{\varphi }^{+}\rangle$ and $|{\varphi }^{-}\rangle$ display initially an entropy of approximately zero, indicating that the system starts in a pure state. Over time, the entropy increases to about 1.85 bits, suggesting a process of mixing, likely due to decoherence. Notably, the identical evolution of these states implies that the dynamics are phase-insensitive, meaning that the system’s evolution does not depend on the specific phase relationship between the states. For the product states $|{\psi }_{prod}\rangle$ and $|{\psi }_{prod}^{\prime }\rangle$, the initial entropy of the full system is also zero, but the reduced entropies are approximately 0.72 bits and 0.97 bits, respectively. As the system evolves, the entropy increases, reaching values of 1.3 bit and 0.8 bits. This growth in entropy suggests that entanglement is generated over time, transforming the initially separable states into partially entangled states. The differences in the entropy values indicate that the two product states undergo different degrees of entanglement formation, possibly due to variations in their initial superposition coefficients.

For the mixed states ${\rho }_1$ and ${\rho }_2$, the initial reduced entropies are 1 bit and approximately 0.65 bits. The entropy of ${\rho }_1$ fluctuates between 0.75 and 1.75 bits before stabilizing at 1.6 bits, while ${\rho }_2$ reaches a peak of 1.9 bits and settling there. This suggests more complex dynamical behavior compared to the Bell and product states, possibly due to the probabilistic mixture of pure states in their initial configurations. The long-term stabilization of entropy values indicates that the system reaches a steady state, reflecting a balance between entanglement generation and dissipation effects. In overall, all simulations generate steady states with a maximal value of $S\approx 1.85$, it is however noteworthy that the interval $t\in [0,2]$ represents a transition phase with decoherence and coherence which generate one-hump camel-like entropy behaviour, particularly of the product and mixed states. The Bell states display a lesser undulation in this interval, which however, is more pronounced in the entanglement entropy seen in next section.

Regarding the entropy plots (see Figure 1), we observe a striking similarity to the spin-up electron fraction measured for a single phosphorus donor atom in silicon subjected to microwave radiation, as reported by Pla et al. [23], however, with a considerable damping effect. Specifically, their results show Rabi oscillations of the spin-up fraction, which we find mirrored in the non-monotonic time evolution of the von Neumann entropy in our system for both product states ($\psi$, ${\psi }^{\prime }$) and mixed states (${\rho }_1$, ${\rho }_2$). For the product states in Figure 1, particularly ${\psi }^{\prime }$ and ${\rho }_1$, we observe the entropy behavior as “camel-like” due to the characteristic rise, dip, and subsequent increase, resembling a double-humped structure in the entropy curves. Indeed, for the mixed states in Figure 1, both ${\rho }_1$ and ${\rho }_2$ exhibit a rise-dip pattern, with ${\rho }_1$ (blue curve) showing a more pronounced ‘camel-like” shape, 1.75 bits, while ${\rho }_2$ (red curve) reaches a higher entropy of 1.85 (as the Bell states) bits with a subtler dip. This non-monotonic pattern in both cases can be attributed to coherent dynamics driven by noncommuting terms in the Hamiltonian which induce damped Rabi oscillations [10,22]. These oscillations cause the qubit’s state to evolve coherently, leading to a transient oscillatory profile in the entropy when traced over a subsystem; however, weak dissipation, as introduced by the dissipation term $\gamma \mathcal{L}\left(\rho \right)$, damps these oscillations after a single cycle, resulting in the observed “camel-like” shape before the entropy stabilizes at a steady value.

We add here an additional analysis of the evolution of the entropy for the Bell states using a parameterized coefficient of $\gamma$. We parameterize $\gamma ={cos}^2\varphi$ on the domain $\varphi \in [0,\pi /2]$. By this we obtain a map of the effect of the relaxation factor $\gamma$ on the unit interval, providing a smooth variation for the entropy evolution analysis, which is shown in Figure 2.

We find an interesting result in the plot of the parameterized under the framework given in [2] with relaxation factor $\gamma$ (Figure 2), where the maximal number of undulations in the entropy (between 1 and 2 bits) occurs for $\varphi$ values approximately between $\varphi \approx 1.0\approx {\displaystyle \frac{1}{3}}\pi$ radians and $\varphi \approx 1.2\approx {\displaystyle \frac{2}{5}}\pi$ radians. The special behavior of the entropy observed in Figure 2 is aligned with the damped Rabi oscillations detected for spin-fractions of single silicon atoms studied by Pla et al. [23]. In this study, several different spin oscillations were observed after various input powers were applied to a single silicon atom device. The plots of these different spin-oscillations result as highly similar to the various plots represented in Figure 2, indicating that the Khrennikov set-up of the GKSL equation may be suitable for the calculation of properties of such spin-fractions from single atom systems and entropy properties therein.

The Rabi-oscillations we observe in Figure 2 signify intermittent entanglement degradation and revival, arising from the competition between the Hamiltonian’s coherent dynamics and the dissipative action of the ${\sigma }_{+}\otimes I$ operator. Furthemore, the undulations in the entropy evolution for $|{\mathrm{\Phi }}^{-}\rangle$ under $\gamma ={cos}^2\varphi$ are signatures of non-Markovian dynamics [30], violating the Markovian monotonicity condition ${\partial }_{t}S\left(t\right)\ge 0$ where these oscillations lead to an increasing number of metastable entanglement revival via Liouvillian exceptional points [31]. These points are degeneracies in quantum open systems (and classical systems as well) and have significant relevance to physics and optics [31]. Such exceptional points define an entropy evolution where two or more eigenvalues coalesce [32] and thus revoke that $\mathrm{\Phi }$-type Bell states could be optimally protected against external influence (i.e. decoherence channels) when $\gamma =0.1-0.3$, making them promising for quantum devices to reduce noise and error.

3.2. Entanglement Entropy

The entanglement entropy is defined by the von Neumann entropy of the reduced density matrix and is employed by being computed for one of the states, working as a key metric for understanding the coherence and mixedness of quantum states over time. Following the top plot in Figure 3, we focus on the Bell states $|{\mathrm{\Phi }}^{+}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left(\right|00\rangle +|11\rangle )$ and $|{\mathrm{\Phi }}^{-}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left(\right|00\rangle -|11\rangle )$, both initially maximally entangled with $S\left(A\right)=1$ (due to normalization, scaling to a maximum value of 1 for clarity). As time progresses, the entropy (represented by a single red curve, as both states evolve identically) dips slightly around $t=1$, reflecting the dissipative influence of the raising operator, which nudges the system toward the ground state $|00\rangle$. However, the entropy stabilizes near 1, significantly above zero, indicating robust residual entanglement. This persistence suggests that the Hamiltonian H, by facilitating transitions between $|00\rangle$ and $|11\rangle$, effectively counteracts the disentangling effects of dissipation, preserving quantum correlations in the long-time limit.

For the product states $|{\psi }_{prod}\rangle$ and $|{\psi }_{prod}^{\prime }\rangle$, the entanglement entropy dynamics reveal three distinct phases. Initially, both states show zero entanglement entropy (as expected for separable states), followed by a sharp increase to approximately 0.8-0.9 bits, indicating significant entanglement generation under the system’s dynamics. The subsequent plateau around 1 bits demonstrates stable entanglement preservation, with the states remaining highly (though not maximally) entangled. The close similarity between both curves suggests the dynamics affect different product states in qualitatively similar ways, while the slight difference in their steady-state values may reflect varying degrees of correlation in their initial configurations. The mixed states ${\rho }_1$ and ${\rho }_2$ (Figure 3) exhibit more complex behavior, following a damped Rabi-oscillation mode, beginning with zero entanglement despite their mixed nature, which highlights their preparation using mixtures of product states. Both states rapidly develop entanglement, with both reaching maximal entanglement (1 bit) and also dipping like the Bell states in the interval $I=[1,2]$. This non-monotonic evolution suggests competing dynamics with initial entanglement generation through unitary evolution followed by partial disentanglement, due to decoherence or specific Hamiltonian terms. We pay particular attention to the interval $I=[0-2]$, which aligns with the definition of one-hump in 1 and study it further by the mutual information with regards to all states.

3.3. Mutual Information

We study the physical implications that the camel like framework [2] has for the Bell-states, the product states and the mixed states with emphasis to definition 1. For this, we calculate the mutual information, which is a measure of the amount of information exchanged between the states during the Lindbladian evolution of the density matrices. In Figure 4 we see in the top plot the mutual information evolution of the Bell states $|{\mathrm{\Phi }}^{+}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left(\right|00\rangle +|11\rangle )$ (red) and $|{\mathrm{\Phi }}^{-}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}\left(\right|00\rangle -|11\rangle )$ (blue). Initially, both state. Over time, $I(A:B)$ decreases monotonically to 0.15 bits by $t=10$, driven by decoherence from the Lindblad term with collapse operator $C={\sigma }_{+}\otimes I$. Despite this, the entanglement entropy $S\left({\rho }_A\right)\approx 1$ bit (see Figure 4) indicates persistent quantum entanglement, maintained by the Hamiltonian $H={\sigma }_x\otimes {\sigma }_x$. The dissipation drives the first qubit toward $|0\rangle$, increasing $S\left({\rho }_{AB}\right)$ to 1.75 bits by $t\approx 2$ (data not shown due to graphical limitations), while the second qubit’s entropy decreases to $S\left({\rho }_B\right)\approx 0.75$ bits due to asymmetric dissipation. This yields $I(A:B)=S\left({\rho }_A\right)+S\left({\rho }_B\right)-S\left({\rho }_{AB}\right)\approx 1+0.75-1.75=0$, though the plot shows a residual 0.1 bits, possibly due to incomplete decoherence. The decline in mutual information suggests that the GKSL framework induces classical-like information exchange, despite persistent quantum entanglement.

The middle plot in Figure 4 compares product states where we have an initially entangled state $|{\psi }_{prod}\rangle$ (blue) and a separable state $|{\psi }_{prod}^{\prime }\rangle =|{\varphi }_A\rangle \otimes |{\varphi }_B\rangle$ (red dashed). For $|{\psi }_{prod}\rangle$, the mutual information starts at 0 bits, peaks at 0.35 bits around $t=1$, and stabilizes at 0.08 bits, reflecting a decay of $\sim 0.3$ bits. This behavior arises from the competition between entanglement generation by $H={\sigma }_x\otimes {\sigma }_x$ and correlation loss to the environment at a rate $\gamma =1$. In contrast, $|{\psi }_{prod}^{\prime }\rangle$ remains at 0 bits, confirming the absence of quantum correlations throughout the evolution.

The bottom plot in Figure 4 examines mixed states ${\rho }_1$ (red) and ${\rho }_2$ (blue). For ${\rho }_1$, $I(A:B)$ starts at 0 bits, peaks at 0.75 bits around $t=1$, and stabilizes at 0.1 bits. For ${\rho }_2$, it begins at 0.8 bits, dips to 0.2 bits, before increasing to 0.35 bits, and also stabilizes at 0.1 bits. Both mixed states exhibit a transition around $t\in [1,1.5]$, with ${\rho }_2$ showing a higher initial mutual information, suggesting stronger initial dependence between subsystems.

3.4. Rate of Information Exchange

The low steady-state mutual information across all states (0.1–0.25 bits) which we observed in Figure 4 indicates that the Lindblad framework described in Section 1.3 generally drives systems toward classical-like information exchange. However, comparing the transition period $t\in [1,2]$ in Figure 4 with the entanglement entropy plots in Figure 3, we identify a key transition state for $|{\psi }_{prod}\rangle$, $|{\psi }_{prod}^{\prime }\rangle$, ${\rho }_1$, and ${\rho }_2$, where the rate of change of entanglement entropy is maximized. This rate is given by

$${\displaystyle \frac{d{S}_A}{dt}}=-\mathrm{Tr}\left({\displaystyle \frac{d{\rho }_A}{dt}}(log{\rho }_A+1)\right),$$

(30)

where A represents the reduced state of $|{\psi }_{prod}\rangle$, $|{\psi }_{prod}^{\prime }\rangle$, ${\rho }_1$, or ${\rho }_2$. The maximal rate of information exchange is thus defined as

$$\mathrm{maximal}\mathrm{rate}\mathrm{of}\mathrm{information}\mathrm{exchange}=\underset{\lambda }{max}\left({\displaystyle \frac{d{S}_A}{dt}}\right).$$

(31)

This period of maximal rate corresponds to the largest increase in entanglement entropy, contrasting with the overall decrease driven by dissipation. The plot of ${\displaystyle \frac{d{S}_A}{dt}}$ is shown in Figure 5 and highlights this transition, aligning with the mutual information peaks around $t\in [1,2]$, which is a period contained in the one-hump interval defined in 1, and defined in definition 1. Thus, we observe that the transition state ($f\left(x_b\right)$) from 1 corresponds to the maximum of the first derivative of the entanglement entropy for all states. This is particularly evident in the rate of change of the information exchange for the Bell states (Figure 5), which exhibit only a minor “hump” in the entanglement entropy within the time interval $t\in [1,2]$ (see Figure 3). This feature is reflected in Figure 5 as an immediate maximum followed by a minimum in the rate of change of the entanglement entropy.

The maximum in the derivative highlights the period $t\in [1,2]$ as the interval of the most rapid rate of information exchange across all simulated systems. This observation suggests that the maximum of the rate of change of the entanglement entropy may serve as a key indicator for identifying whether the information exchange is transitioning toward classical or quantum-like behavior. Specifically, if the information exchange is developing towards classical, we expect the rate of change to decrease, while if it develops towards quantum-like, we anticipate a maximum in the rate of change of the entanglement entropy. This behavior is further contextualized by the mutual information, which decreases more slowly for the Bell states during this period, as shown in Figure 5. A more pronounced effect is observed for the product states and mixed states in Figure 5 (compare the y-axis values between the three plots), where the maxima in the rate of information exchange align with those in Figure 4. Thus, calculating the rate of information exchange emerges as a critical method for predicting whether a system is transitioning from classical to quantum information exchange, or vice versa.

Notably, quantum correlations, such as those measured by quantum discord or CHSH parameters (discussed in the next sections), are often insufficient to determine whether a system is evolving toward classical or quantum information exchange. Instead, these correlations primarily indicate whether the system currently exhibits quantum-like or classical information exchange. In contrast, when the entanglement entropy can be resolved, its rate of change provides a robust prediction of the system’s evolution toward either quantum-like or classical information exchange.

In Figure 5, we observe that the transition state ($f\left(x_b\right)$) from definition 1 corresponds to the maximum of the first derivative of the entanglement entropy for all states. This is particularly evident in the Bell states, which exhibit a minor “hump” in the entanglement entropy within the time interval $t\in [1,2]$ (see Figure 3). This feature is reflected in Figure 5 as an immediate maximum followed by a minimum in the rate of change of the entanglement entropy. The maximum in the derivative highlights the period $t\in [1,2]$ as the interval of the most rapid rate of information exchange across all simulated systems. This observation suggests that the maximum of the rate of change of the entanglement entropy may serve as a key indicator for identifying whether the information exchange is transitioning toward classical or quantum-like behavior. Specifically, if the information exchange is classical, we expect the rate of change to decrease, whereas if it is quantum-like, we anticipate a maximum in the rate of change of the entanglement entropy. This behavior is further contextualized by mutual information, which decreases more slowly for the Bell states during this period, as shown in Figure 4. A more pronounced effect is observed for the product states and mixed states in Figure 4, where the maxima in the rate of information exchange align with those in Figure 5. Thus, calculating the rate of information exchange emerges as a critical method to predict whether a system is transitioning from classical to quantum information exchange. Interestingly, these results are supported by recent studies [33], which show that quantum-to-classical transitions are not derived from operator non-commutativity but from the dynamical properties of the systems in Hilbert space, as we consider here in our calculations.

Notably, quantum correlations, such as those measured by quantum discord or CHSH parameters (discussed in the next sections), are often insufficient to determine whether a system is evolving toward classical or quantum information exchange. Instead, these correlations primarily indicate whether the system currently exhibits quantum-like or classical information exchange. In contrast, when the entanglement entropy can be resolved, its rate of change provides a robust prediction of the system’s evolution toward either quantum-like or classical information exchange. The former is associated with phenomena such as superdense coding or quantum teleportation, while the latter corresponds to classical phenomena, such as electromagnetic effects on electrons or the Compton effect [34].

3.5. Quantum Discord

Quantum discord gives a measure on the level of “quantumness” of a system, hence to which degree the behaviour of the qubits is classical or quantistic. This is important to be determined, since the relationship between “quantumness” and quantum relations is not always as expected. For instance, a system can behave classically, but still be entangled, which indicates that quantum discord is an important property of systems, and can give better classifications of their behaviour in an open quantum system simulation. We calculate and plot the quantum discord of our model states, the Bell states, the product states in (19), (21) and the mixed states in (23), (24) by the Khrennikov-picture, using the Hamiltonians and dissipation operators given by (16) (Figure 6). In this figure, we see that we obtain a stabilization of the quantum discord of the Bell-states to a stable steady state after experiencing a small “hump” at the beginning of the simulation. This is in agreement with the evolution of the mutual information (Figure 4). The fact that the Bell-states start from quantum discord of D=1 is correct as quantum discord measure non-classical correlations in mixed states, which Bell states initially are not. Clearly, the Bell-states lose coherence and their correlation turns into a classical correlation, to a higher extent, by attaining a steady state with a discord of D=0.1. It is noteworthy that the decrease in discord is showing a “hump” during the transition state period of $t\in [1,2]$, as for the previous correlations and as defined by definition 1.

The quantum discord evolution in Figure 6 reveals distinct dynamics for the entangled ($|{\psi }_{\mathrm{prod}}\rangle$) and separable ($|{\psi }_{\mathrm{prod}}^{\prime }\rangle$) states. The entangled state displays initial discord $D=0$ bits, confirming non-classical correlations, which suddenly increases to $D\approx 0.2$ bits in the interval $t\in [1,2]$ and by $t=10$ it stabilizes to 0.05 bits. In contrast, the separable state maintains $D=0$ throughout, as expected for a classically correlated system. Notably, the non-zero steady-state discord for $|{\psi }_{\mathrm{prod}}\rangle$ implies persistent quantum correlations despite environmental interaction.

The mixed states, ${\rho }_1$ and ${\rho }_2$ behave similarly to the product states, however with a stronger synchronicity during the evolution. Their synchronous behaviour rests in their higher purity than the product states, however, their complete loss of quantum discord into the classical realm indicate that however pure, their combination generates classical relations for their substates (Figure 6).

From Figure 6 we see a further confirmation that the period $t=[1,2]$ is the critical transition which delineates that the system is moving towards quantum-like information exchange by increasing the quantum discord for the product states and inducing a slower rate of negative change of the discord for the Bell state and the the mixed state ${\rho }_1$. Indeed, this is equally well-confirmed by the rate of change of the entanglement entropy of all the states, as an indicator of whether a system moves towards quantum-like information exchange or classical information exchange.

From Figure 6 we see a further confirmation that the period $t=[1,2]$ is the critical transition which delineates that the system is moving towards quantum-like information exchange (by increasing discord for the product states and experiencing a slower rate of negative change of the discord for the Bell state and the the mixed state ${\rho }_1$ (Figure 6). Indeed, this is equally well-confirmed by the rate of change of the entanglement entropy of all the states, as an indicator of whether a system moves towards quantum-like information exchange or classical information exchange.

3.6. Bell-CHSH Parameter

For a proper quantum steering assessment, we evaluate the evolution of the CHSH parameter. Here we study the evolution of the parameter that allows for identifying the level of steering that occurs for the Bell states ${\mathrm{\Phi }}^{+}$, ${\mathrm{\Phi }}^{-}$, the product states in (19) and the mixed states in (23). We solve the Lindblad equation by the Runge-Kutta method of 4th order, under Lindblad dynamics governed by a Hamiltonian $H={\sigma }_x$ and the collapse operator $C={\sigma }_{+}\otimes \mathbb{I}$ and plot the evolution in Figure 7.

The calculations show that the entangled Bell states $|{\mathrm{\Phi }}^{+}\rangle$ and $|{\mathrm{\Phi }}^{-}\rangle$ exhibit distinct behaviors in their CHSH parameter (Figure 7). The CHSH parameter for both pure states starts at the quantum upper bound $2\sqrt{2}$ which correctly reflects property of pure states by the CHSH parameter, which however decays rapidly within a few time steps, reflecting the loss of quantum coherence due to dissipation by the model dissipation operator and Hamiltonian. This describes that the presence of locally hidden variables, rapidly increases during the evolution of the GKSL equation with the given framework by [2]. The critical transition state, referred to as the “camel-hump” in definition 1, demonstrates that surpassing the maximum leads to decoherence for the pure states ${\mathrm{\Phi }}^{+}$ and ${\mathrm{\Phi }}^{-}$ specifically within the interval $t\in [0,1]$. This interval corresponds to the period $t\in [0,1]$ in Figure 3, where the entanglement entropy rate, ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}$, is significantly negative. Clearly, by the case of ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}<0$ in this interval leads to classical evolution and hence, as we see here for the CHSH parameter, a loss of coherence for pure states.

The evolution of the CHSH parameter under Lindblad dynamics for the product states $|{\psi }_{prod}\rangle$ and $|{\psi }_{prod}^{\prime }\rangle$ highlights the distinct behavior of these states under dissipation. For $|{\psi }_{prod}\rangle$, given by the initial CHSH parameter is $\sim 1.4$, indicating classical correlations with ${\psi }_{prod}^{\prime }$. However, dissipation causes a rapid decrease in the CHSH parameter to approximately 0, reflecting a significant loss of coherence. This final value indicates that both product state retain weak correlations. When we consider the mixed states we find interesting results. For ${\rho }_2$, the mixed state defined as

$${\rho }_2=q|{\psi }_5\rangle \langle {\psi }_5|+(1-q)|{\psi }_6\rangle \langle {\psi }_6|,$$

with $q=0.5$ which starts with a CHSH value of $\sim 1$, reflecting classical correlations between the subsystems under interaction with ${\rho }_1$. The CHSH parameter decreases to a minimum of 0.22 by time step 2, suggesting dissipation suppresses these correlations. It then increases and stabilizes at 0.4 in the steady state, indicating the state retains weak classical correlations, but the value remains below the quantum threshold of 2. For ${\rho }_1$, defined as

$${\rho }_1=p|{\psi }_3\rangle \langle {\psi }_3|+(1-p)|{\psi }_4\rangle \langle {\psi }_4|,$$

with $p=0.5$, the CHSH parameter starts at -0.5, which is unusual. This negative value suggests that classical correlations or interference between the components $|{\psi }_3\rangle$ and $|{\psi }_4\rangle$ dominate initially. Over time, the CHSH value increases to -0.34, then asymptotically settles at -0.6, indicating that ${\rho }_1$ remains in a classical regime with no significant quantum steering potential in the steady state. The negative CHSH value for ${\rho }_1$ suggests interference or classical correlations, rather than quantum coherence. As ${\rho }_1$ is a mixture of $|{\psi }_3\rangle$ and $|{\psi }_4\rangle$, the state lacks the necessary alignment for generating positive quantum correlations, leading to a negative value. Indeed, the CHSH components of ${\rho }_1$ were calculated by the Python script as:

$$E_{Z\otimes B}=0.1414(\mathrm{weakly}\mathrm{positive}),$$

$$E_{Z\otimes B^{\prime }}=-0.7071(\mathrm{strongly}\mathrm{negative}),$$

$$E_{X\otimes B}=0.0,$$

$$E_{X\otimes B^{\prime }}=0.0,$$

where X and Z are the respective Pauli matrices.

These values confirm that the simulation proceeded correctly. We can therefore affirm the negative values of ${\rho }_1$, and add that the negative CHSH parameter indicates that the correlations in the state are anti-aligned with the measurement settings used in the CHSH inequality. This does not mean the state is unphysical or lacks entanglement; it simply means that the specific measurement settings used are not optimal for detecting the nonlocality of the state.

In overall, we can see that the Lindbladian evolution by the camel-like behaviour in [2] gives no quantum steering between any steady states of the considered states, as their respective CHSH parameter fall below the classical threshold. A final word for the CHSH parameter is worth mentioning. The CHSH parameter cannot be used to predict whether the system evolves towards quantum-like or classical information exchange, as its evolution must violate the upper bound of 2 before the system has a direct quantum steering relationship between two states. The increase of a systems’ CHSH parameter over time does not warrant violation of the the upper classical bound, therefore, we cannot conclude whether the system is evolving towards a classical or quantum like information exchange by relying on the CHSH parameter. Such an evolution can only be predicted using the rate of the change of the entanglement entropy and assessed additionally by calculating the quantum discord. This is so, since the former highlights an evolution towards classicality (if ${\displaystyle \frac{d{S}_e\left(t\right)}{dt}}<0$), and an evolution towards quantumness (if ${\displaystyle \frac{d{S}_e\left(t\right)}{dt}}>0$), over a particular period of time $I\subset J$, where J is the entire period of the evolution of the quantum system.

3.7. Entropy and Entanglement Entropy by Parametrized Coefficients of Mixed States

The goal of this final part of the analysis is to study the camel-hump-landscape defined by definition 1 for the entropy and entanglement entropy for mixed quantum states parameterized by $\varphi$ under Lindblad evolution. As shown in the method section, mixed states are constructed as convex combinations of pure states, where the coefficients $\alpha ,\beta ,\gamma ,$ and $\delta$ are defined by trigonometric functions of $\varphi$ (see Method section). The mixed state that is considered for this final calculation is defined as:

$${\rho }_0=0.5|{\psi }_1\rangle \langle {\psi }_1|+0.5|{\psi }^{\prime }\rangle \langle {\psi }^{\prime }|,$$

(32)

where the pure states $|{\psi }_1\rangle$ and $|{\psi }^{\prime }\rangle$ are given by:

$$|{\psi }_1\rangle =\alpha |00\rangle +\beta |01\rangle ,|{\psi }^{\prime }\rangle =\gamma |10\rangle +\delta |11\rangle .$$

(33)

So the convex combination of these two pure states forms the mixed state as:

$${\rho }_0=0.5\left(\alpha \right|00\rangle +\beta |01\rangle )({\alpha }^{*}\langle 00|+{\beta }^{*}\langle 01\left|\right)+0.5\left(\gamma \right|10\rangle +\delta |11\rangle )({\gamma }^{*}\langle 10|+{\delta }^{*}\langle 11\left|\right).$$

(34)

These states evolve according to the Lindblad master equation, governed by the camel-like evolution Hamiltonian as given by [2]. The von Neumann entropy, $S\left(\rho \right)=-\mathrm{Tr}(\rho log\rho )$, is calculated at each time step to quantify the uncertainty or mixedness of the state. By varying $\varphi \in [-\pi ,\pi ]$, the dependence of steady-state entropy on the parameterization is analyzed, providing insight into how the structure of the coefficients influences the quantum system’s dynamics.

By Figure 8 (left), we see that the entropy entropy dynamics observed in the simulations reveal a periodic behavior with a period of $2\pi$ determined by the trigonometric parametrization of the coefficients $\alpha ,\beta ,\gamma ,$ and $\delta$ using $cos\left(\varphi \right)$ and $sin\left(\varphi \right)$. This periodicity, inherent to the symmetry of the parameterization, is reflected in the steady-state entropy values, which exhibit maxima at $\varphi =-\pi /4$ and $\varphi =3\pi /2$, and minima at $\varphi =\pi /4$ and $\varphi =-3\pi /2$, accordingly with the nature of ${cos}^2\varphi$ and ${sin}^2\varphi$ by the normalization condition ${cos}^2\varphi +{sin}^2\varphi =1$. These points correspond to configurations where the quantum system achieves the highest and lowest levels of mixedness, respectively (Figure 8).

Figure 8 (left) shows that the periodic high-entropy states yield an increased uncertainty in the system, while the low-entropy states imply closer alignment to pure states (compare with Figure 2). However, a key-insight is that the camel-like landscape gradually vanishes as the entropy grows by the periodic change of the value of $\varphi$ (Figure 8, left). By considering definition 1, we see in Figure 8 (left) that as the value of the coefficients evolves towards $\varphi =-\pi /4$ and $\varphi =3\pi /2$ as the transition state. This implies that the transition state is no longer a transition state followed by a relaxed steady-state, but a step towards a steady state of higher entropy than the transition state itself. By this we can conclude that the camel-like behaviour of the entropy given by [2] given in Figure 2 is sensitive to the parametrization of the coefficients of the quantum states.

Nevertheless, we find an equally interesting result in the evolution of the entanglement entropy (Figure 8, left). The parametrization of the quantum state coefficients by $\varphi$ has no effect on the entropy of the reduced density matrices. This is expected, as entanglement entropy depends only on the structure of the total state and not on any direct interaction between ${\psi }_1$ and ${\psi }^{\prime }$. As a result, the entanglement entropy for the mixed states follows the same transition-state-based landscape as when the coefficients are fixed (Figure 3).

We can now construct the following theorem, relating the behaviour of the entropy in the entropy, entanglement entropy, mutual information, rate of entanglement entropy, the mutual information and the quantum discord in Figure 2, Figure 3, Figure 4, Figure 6, Figure 8.

Theorem 1.

(Classical versus non-classical evolution)

Let $\rho \left(t\right)$ be the density matrix of a quantum system at time t, and let $S_E\left(t\right)$ denote the entanglement entropy of a bipartite subsystem $A\cup B$ of the system, where $S_E\left(t\right)=-Tr({\rho }_A\left(t\right)log{\rho }_A\left(t\right))$ and ${\rho }_A\left(t\right)={Tr}_B\left(\rho \left(t\right)\right)$ is the reduced density matrix of subsystem A. Let furthermore J be the entire period of the evolution of the states under a completely positive trace-preserving (CPTP) map, and let $I\subset J$.

• Assume:

• The system evolves under a completely positive trace-preserving (CPTP) map (i.e. the GKSL equation), which ensures the physicality of the evolution.
• The entanglement entropy $S_E\left(t\right)$ is differentiable with respect to time t.

• Then:

• If ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}<0$ for all t in some interval I, the system’s information exchange evolves towards classical-like behavior in I. This implies a reduction in quantum correlations, such as entanglement, and a shift towards classical correlations, where the system behaves more classically.
• If ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}>0$ for all t in some interval I, the system’s information exchange evolves towards quantum-like behavior in I. This implies an increase in quantum correlations, such as entanglement, which enables non-classical phenomena like quantum teleportation.
• The measured substate shifts the quantum discord toward 0 by adding classical contributions, while the unmeasured substate shifts it toward 1 if ${\displaystyle \frac{d{S}_E}{dt}}>0$, and conversely if ${\displaystyle \frac{d{S}_E}{dt}}<0$.

Proof. 

Since the system evolves under a CPTP map (completely positive trace-preserving map), the density matrix $\rho \left(t\right)$ evolves as $\rho \left(t\right)={\mathcal{E}}_t\left(\rho \left(0\right)\right)$, where ${\mathcal{E}}_t$ is a CPTP map. This guarantees that ${\rho }_A\left(t\right)$ remains a valid density matrix. Furthermore, we base our proof on the universal relation:

$$\text{Bell-nonlocality}\supset \mathrm{Quantum}\mathrm{steering}\supset \mathrm{Quantum}\mathrm{entanglement}\supset \mathrm{Quantum}\mathrm{discord}$$

(35)

The entanglement entropy $S_E\left(t\right)$ is defined as the von Neumann entropy of the reduced density matrix ${\rho }_A\left(t\right)$:

$$S_E\left(t\right)=-\mathrm{Tr}({\rho }_B\left(t\right)log{\rho }_B\left(t\right)).$$

By assumption, $S_E\left(t\right)$ is differentiable, and we analyze the implications of ${\displaystyle \frac{d{S}_E}{dt}}$ and prove point 1. and 2.

Case 1: If ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}<0$ for all t in some interval I, this implies that the entanglement entropy is decreasing over time. Since entanglement entropy measures quantum correlations, its decrease suggests that quantum information is being lost, likely due to decoherence. This results in a reduction of quantum correlations, leaving only classical correlations dominant, signifying classical-like behavior.

Case 2: If ${\displaystyle \frac{d{S}_E\left(t\right)}{dt}}>0$ for all t in some interval I, then entanglement entropy is increasing. Since increasing entropy indicates an increase in quantum correlations, the system is evolving towards more non-classical correlations, including entanglement and quantum discord. This suggests that the system’s information exchange is dominated by quantum effects, signifying quantum-like behavior.

In order to prove point 3., we analyze explicitly the implications that the theorem has when we apply a measurement to either of the substates of the system, A or B, under the evolving CPTP-map.

We start with the definition of the quantum discord:

$$\mathcal{D}(A:B)=I(A:B)-C(A:B),$$

(36)

where

$$I(A:B)=S_A+S_B-S_{AB},$$

(37)

and (using an equivalent formulation)

$$C(A:B)=S_A-\underset{\left\{M_a\right\}}{max}\sum _a{p}_{a}S\left({\rho }_{B|a}\right).$$

(38)

Thus,

$$\mathcal{D}(A:B)=S_B-S_{AB}+\underset{\left\{M_a\right\}}{max}\sum _a{p}_{a}S\left({\rho }_{B|a}\right).$$

(39)

Taking the time derivative (and assuming that the optimal measurement varies smoothly so that the derivative passes through the maximization) yields

$${\displaystyle \frac{d\mathcal{D}(A:B)}{dt}}={\displaystyle \frac{d{S}_B}{dt}}-{\displaystyle \frac{d{S}_{AB}}{dt}}+{\displaystyle \frac{d}{dt}}\left\{\underset{\left\{M_a\right\}}{max}\sum _a{p}_{a}S\left({\rho }_{B|a}\right)\right\}.$$

(40)

We now consider two cases,

Variation in $S_A$: In the definitions above, $S_A$ appears in both $I(A:B)$ and $C(A:B)$ and cancels exactly. Thus, variations in $S_A$ do not contribute to ${\displaystyle \frac{d\mathcal{D}(A:B)}{dt}}$.

Variation in $S_B$: In Eq. (40), if we assume that the contributions from ${\displaystyle \frac{d{S}_{AB}}{dt}}$ and

$${\displaystyle \frac{d}{dt}}\left\{\underset{\left\{M_a\right\}}{max}\sum _a{p}_{a}S\left({\rho }_{B|a}\right)\right\}$$

are constant or negligible, then

$${\displaystyle \frac{d\mathcal{D}(A:B)}{dt}}\approx {\displaystyle \frac{d{S}_B}{dt}}.$$

Consequently,

• If ${\displaystyle \frac{d{S}_B}{dt}>0}$, then ${\displaystyle \frac{d\mathcal{D}(A:B)}{dt}>0}$.
• If ${\displaystyle \frac{d{S}_B}{dt}<0}$, then ${\displaystyle \frac{d\mathcal{D}(A:B)}{dt}<0}$.

Therefore, when performing a measurement on a bipartite system, the sign of ${\displaystyle \frac{d{S}_B}{dt}}$ directly contributes to the sign of the change of the quantum discord. We can interchange the roles of the substates A and B, such that A does not vanish but B does instead in the initial equation of the quantum discord. This would change the substate that is subjected to a measurement. Hence, the substate that is measured contributes to the change of the sign of the quantum discord towards 0 (by addition classical contributions to the system), while the substate that is not measured contributes to the change of the sign of the quantum discord towards 1, if ${\displaystyle \frac{d{S}_E}{dt}}>0$, and conversely if ${\displaystyle \frac{d{S}_E}{dt}}<0$. Hence, we have completed the proof which holds for a bipartite system in the absence and presence of a measurement.

   □

Remarks: 

• The theorem assumes that entanglement entropy is a valid measure of quantum entanglement, which holds for pure and product states, by considering the results from this study. The theorem holds also for mixed states since the parametrization by $\varphi$ of the coefficients of the subsystems does not affect the entanglement entropy landscape.
• The theorem applies to bipartite systems. For multipartite systems, the behavior of entanglement entropy can be more complex, and additional considerations may be necessary.
• The theorem assumes differentiable entanglement entropy, which may not hold in all cases (e.g., during quantum phase transitions or in open quantum systems with non-Markovian dynamics or under discrete Markovian quantum dynamics).

3.8. Analytical Approach to the GKSL Equation

As we have now introduced a set of methods for solving the GKSL equation analytically, and provided an example, we shall enter into the calculations for solving the GKSL equation under a specific framework, namely the Khrennikov framework. Analytical solutions of this open quantum system are important as these functions can provide crucial insights into how entropy scales in different dissipation regimes and as we obtain closed-form expressions for the entropy as a function of time, they enable us to study the entropy growth or decay [17]. The solutions can thus reveal how the entropy behaves near or at equilibrium, particularly identifying the steady-state entropy and revealing its dependence on the system’s constraints [35]. We can also derive the exact description of how entropy scales with time, by using the derived analytical solutions which can evolve according to exponential decay or power-law growth, depending on the nature of the dissipation [30]. Furthermore, the analytical solutions can aid to facilitate the study of quantum-to-classical information exchange transitions, and one can identify how this can be caused by change in decoherence or quantum coherence affecting the system’s entropy [36]. Indeed, the dissipation operators govern this role on the transition from classical-to-quantum or vice versa, and hence, the analytical solutions will prove valuable in concert with the structure of the dissipation operators to reveal this mechanism. This is essential in understanding the role of quantum effects, where entropy might increase at a different rate compared to classical systems, and has also with strong relevance to theorem 1.

3.8.1. Formulation of the Problem

We consider the evolution of a two-level quantum system (qubit) under the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation:

$${\rho }^{\prime }\left(t\right)=-i[H,\rho \left(t\right)]+\gamma \left(C\rho \left(t\right)C^{†}-{\displaystyle \frac{1}{2}}\{C^{†}C,\rho \left(t\right)\}\right),$$

(41)

where: $\rho \left(t\right)$ is the density matrix of the system, H is the Hamiltonian governing unitary evolution, C is the Lindblad operator representing system-environment interaction, and $\gamma$ is a coupling constant controlling the strength of dissipation, which we calculated to $\gamma ={\displaystyle \frac{1}{2}}$ in the Spectral analysis section.

Following [2], we apply:

$$H={\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),C=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),\gamma ={\displaystyle \frac{1}{2}}.$$

(42)

in the spirit of the transition-entropy regime by Khrennikov [2,7,12,13,15,21], which we studied numerically in the first part of this thesis.

Structure of the Density Matrix

We assume a general initial density matrix of the form:

$$\rho \left(t\right)=\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{12}\\ {\rho }_{21}& 1-{\rho }_{11}\end{array}\right),$$

(43)

where all diagonal entries are real functions, while off-diagonal entries are complex conjugates of one another.

3.8.2. Derivation of the Evolution Equations

Substituting (42) in (41) we solve for each component:

$$-i[H,\rho ]=\left(\begin{array}{cc}i({\rho }_{12}-{\rho }_{21})& i(2{\rho }_{11}-1)\\ i-2i{\rho }_{11}& -i({\rho }_{12}-{\rho }_{21}),\end{array}\right)$$

(44)

$$C\rho C^{†}=\left(\begin{array}{cc}1-{\rho }_{11}& 0\\ 0& 0\end{array}\right),$$

(45)

$$C^{†}C\rho =\left(\begin{array}{cc}0& 0\\ {\rho }_{21}& 1-{\rho }_{11},\end{array}\right)$$

(46)

$$\rho C^{†}C=\left(\begin{array}{cc}0& {\rho }_{12}\\ 0& 1-{\rho }_{11}.\end{array}\right)$$

(47)

Then we get the anti-commutator term:

$${\displaystyle \frac{1}{2}}\left\{C^{†}C,\rho \right\}=\left(\begin{array}{cc}0& {\displaystyle \frac{{\rho }_{12}}{2}}\\ {\displaystyle \frac{{\rho }_{21}}{2}}& 1-{\rho }_{11}\end{array}\right)$$

(48)

combining all terms and we get:

$${\displaystyle \frac{d\rho }{dt}}=\left(\begin{array}{cc}{\displaystyle \frac{1}{2}}\left(-{\rho }_{11}+2i{\rho }_{12}-2i{\rho }_{21}+1\right)& 2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{12}-i\\ -2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{21}+i& {\displaystyle \frac{1}{2}}\left({\rho }_{11}-2i{\rho }_{12}+2i{\rho }_{21}-1\right)\end{array}\right)$$

(49)

3.8.3. Basis States and Initial Conditions

We now investigate the determination of optimal initial conditions for the analytical solution of the GKSL equation for $\rho \left(t\right)$. Guided by Definition 1, we seek initial conditions yielding 3-, 4-, or N-hump solutions in the entropy plots. Considering Definition 1, we first perform numerical calculations for a parametrized superposition of pure states, using one and then two parameters (corresponding to a single point and a curve on the Bloch sphere, respectively). These numerical results are then used to identify which parametrized states best reflect Definition 1 and to determine the relevant intervals for constructing appropriate models of N-hump camel-like entropy. Essentially, camel-like entropy represents a damped Rabi oscillation due to decoherence, eventually reaching a steady state. Based on this, we begin by identifying the initial conditions for the single-parameter states.

3.8.4. Initial Conditions from Parametrization on the Bloch Sphere

We now parametrize the pure states using two parameters, representing a curve on the Bloch sphere. To achieve this, we parametrize into the representation of the pure state superposition on the Bloch sphere.

$$|\psi (\theta ,\varphi )\rangle =cos\theta |0\rangle +sin\theta e^{i\varphi }|1\rangle$$

(50)

which has the corresponding density matrix:

$${\rho }_0=|\psi (\theta ,\varphi )\rangle \langle \psi (\theta ,\varphi )\rangle$$

(51)

explicitly this gives the initial density matrix for the numerical calculation by:

$$\rho =|\psi \rangle \langle \psi |=\left(\begin{array}{cc}{cos}^2\left({\displaystyle \frac{\theta }{2}}\right)& cos\left({\displaystyle \frac{\theta }{2}}\right)sin\left({\displaystyle \frac{\theta }{2}}\right)e^{-i\varphi }\\ cos\left({\displaystyle \frac{\theta }{2}}\right)sin\left({\displaystyle \frac{\theta }{2}}\right)e^{i\varphi }& {sin}^2\left({\displaystyle \frac{\theta }{2}}\right)\end{array}\right).$$

(52)

We proceed to calculate numerical solutions of the GKSL equation (41) with $\gamma ={\displaystyle \frac{1}{2}}$ and the operators defined in (16). Using the Bloch sphere density matrix from (52) we vary $\varphi$ and $\theta$ and identify the suitable initial conditions by the resulting spectrum of numerical entropy solutions is shown in Figure 9 as overlaid plots.

From Figure 9, we can see that six unstable states immediately decohere forming an important class of quantum states with importance for quantum coding purposes. We identify these on the left-most part of the plot and find that these states on the Bloch sphere form the set:

$$\mathcal{S}=\left\{\left({\displaystyle \frac{3\pi }{4}},0\right),\left({\displaystyle \frac{3\pi }{4}},{\displaystyle \frac{3\pi }{2}}\right),\left({\displaystyle \frac{3\pi }{4}},\pi \right),\left({\displaystyle \frac{3\pi }{4}},{\displaystyle \frac{\pi }{2}}\right),\left(\pi ,0\right),\left(\pi ,{\displaystyle \frac{3\pi }{2}}\right)\right\}$$

(53)

which gives the specific states:

$$\begin{array}{cc}\hfill |\psi ({\displaystyle \frac{3\pi }{4}},0)\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle +sin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \\ \hfill |\psi ({\displaystyle \frac{3\pi }{4}},{\displaystyle \frac{3\pi }{2}})\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle -isin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \\ \hfill |\psi ({\displaystyle \frac{3\pi }{4}},\pi )\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle -sin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \\ \hfill |\psi ({\displaystyle \frac{3\pi }{4}},{\displaystyle \frac{\pi }{2}})\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle +isin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \\ \hfill |\psi (\pi ,0)\rangle & =|1\rangle \hfill \\ \hfill |\psi (\pi ,{\displaystyle \frac{3\pi }{2}})\rangle & =-i|1\rangle \hfill \end{array}$$

(54)

which are states that are highly unstable towards decoherence and dephasing. We have a total of six of such states, where the last two in (54) overlap on the South Pole. In order to find how these states appear on the Bloch sphere, we plot the Bloch sphere with these six states arranged as spherical coordinates (Figure 10).

From the plot in Figure 10, we find that the unstable states arise around the South pole of the Bloch sphere, in the well-defined third quarter-circle for $\theta$, at each quadrant intersection for the value of $\varphi$. Notably, the four states on this quarter circle correspond to equal superpositions with distinct relative phases, aligned with the principal axes of the Bloch sphere. These geometrically significant points delineate a partition of the Bloch sphere where entropic instability becomes prominent, particularly through the emergence of rapid decoherence with no oscillatory components, which then decays oscillatory towards the steady state (see Figure 9).

3.8.5. Behaviour of Decoherence-Sensitive States on the Bloch Sphere

The specific distribution of points in $\mathcal{S}$ in (53) consists of four points at the three-quarter circle ($\theta ={\displaystyle \frac{3\pi }{4}}$) and two that overlap completely at the South pole ($\theta =\pi$), and this set forms a $C_4$-symmetry around the z-axis when viewed purely geometrically (Figure 10). These points represent distinct quantum states: $|1\rangle$ (the second point at the South pole in (53) is omitted for simplicity) at the South pole, and four states on the third quarter circle with amplitudes shown in (54). The four states around the third quarter circle have the same superposition amplitudes but differ in their relative phases, aligned with the principal axes in the x-y plane of the Bloch sphere. The behaviour for each point is evaluated for the GKSL operators under the camel-like framework which we analyse one by one:

$$\begin{array}{c}\hfill H={\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),C=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),\gamma ={\displaystyle \frac{1}{2}}\end{array}$$

(55)

For the states on the three-quarter circle at $\theta ={\displaystyle \frac{3\pi }{4}}$, we have superposition states with amplitudes $cos\left({\displaystyle \frac{3\pi }{8}}\right)\approx 0.383$ and $sin\left({\displaystyle \frac{3\pi }{8}}\right)\approx 0.924$:

$$\begin{array}{cc}\hfill \mathrm{At}\varphi =0\text{:}|{\psi }_1\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle +sin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill H|{\psi }_1\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle +sin\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle \hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill C|{\psi }_1\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill C^{†}C|{\psi }_1\rangle & =sin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(59)

Notably, these states are characterized by the special property that $C|\psi \rangle$ produces a state proportional to $|1\rangle$ with amplitude $cos\left({\displaystyle \frac{3\pi }{8}}\right)$, while $C^{†}C|\psi \rangle$ yields the $|1\rangle$ component of the original state. This creates a specific balance between the unitary and dissipative dynamics that leads to maximal decoherence rates.

The states at $\varphi ={\displaystyle \frac{\pi }{2}},\pi ,{\displaystyle \frac{3\pi }{2}}$ follow the same pattern but with different relative phases:

$$\begin{array}{cc}\hfill \mathrm{At}\varphi =\pi /2\text{:}|{\psi }_2\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle +isin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill \mathrm{At}\varphi =\pi \text{:}|{\psi }_3\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle -sin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill \mathrm{At}\varphi =3\pi /2\text{:}|{\psi }_4\rangle & =cos\left({\displaystyle \frac{3\pi }{8}}\right)|0\rangle -isin\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle \hfill \end{array}$$

(62)

In each case, the action of the collapse operator C produces the same output $cos\left({\displaystyle \frac{3\pi }{8}}\right)|1\rangle$, while $C^{†}C$ preserves the relative phase. The specific angle ${\displaystyle \frac{3\pi }{8}}$ appears to maximize the interplay between coherent evolution under H and decoherence under C, creating instability points where the quantum state cannot maintain coherence.

For the South pole state:

$$\begin{array}{c}\hfill \mathrm{South}\mathrm{Pole}(\pi ,0):|1\rangle \end{array}$$

(63)

$$\begin{array}{cc}\hfill H|1\rangle & =|0\rangle \hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill C|1\rangle & =0\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill C^{†}C|1\rangle & =|1\rangle \hfill \end{array}$$

(66)

This state experiences maximal dissipation with $C^{†}C|1\rangle =|1\rangle$, while the collapse operator C annihilates it. The dissipative term dominates completely, making this a strong decoherence point where the state is “trapped” by the dissipation, unable to evolve coherently under the combined dynamics.

A common feature of all five unstable states (or six, if including the final state in (54), which we omit here for simplicity) is that they represent critical configurations where the delicate balance between unitary evolution and dissipation gives rise to either maximal decoherence rates or degenerate dynamics. These properties render them unsuitable as stable quantum states under this particular dissipative channel. A similar scenario arises in the spectral analysis of the governing operator, where the parameter $\gamma$ modulates the system’s oscillatory and exponential components. However, in the current context, the imbalance between unitary evolution and decoherence is especially pronounced for the rapidly decohering states in the set $\mathcal{S}$.

This issue motivates a deeper investigation, specifically, an analysis of the evolution of quantum states on the Bloch sphere. We aim to study the transition from the unstable states in $\mathcal{S}$ to the remaining states on the sphere. This evolution can be rigorously analyzed by constructing a Lyapunov function, derived from an analytic function whose zero set corresponds to the unstable quantum states on the Bloch sphere, with all other points lying in its complement. In the following section, we derive this analytic function, which characterizes the relevant sets of quantum states on the Bloch sphere.

3.9. Topological Analysis on the Bloch Sphere

The set $\mathcal{S}$ in (53) can be classified topologically on the Bloch sphere ${\mathbb{S}}^2$. We establish coordinates on this manifold using the standard spherical parameterization:

$$(x,y,z)=(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta )$$

(67)

This allows us to combine topologically all stable solutions (which decohere slowly under the camel-like framework) as the complement ${\mathbb{S}}^2\setminus X$, forming a non-compact topological space homeomorphic to a 2-sphere with 5 punctures, where each puncture corresponds to one of the highly unstable states in $\mathcal{S}$ in (53) (recall two of the six states are the same on the Bloch spheres’ South pole). The set $\mathcal{S}$ of unstable states in (53) can be interpreted as the zero-set of a map

$$f:{\mathbb{S}}^2\to {\mathbb{R}}^2$$

(68)

where f vanishes precisely at these six specific points. We can express this map as:

$$f(\theta ,\varphi )=(g_1(\theta ,\varphi ),g_2(\theta ,\varphi ))$$

(69)

where $g_1$ captures the pole structure (with zeros at $\theta \in \{0,\pi /4,\pi \}$), and $g_2$ represents the angular dependence (with zeros at specific $\varphi$ values depending on $\theta$). Here we can form the function:

$$\begin{array}{cc}\hfill g_1\left(\theta \right)& =\left(cos\theta +{\displaystyle \frac{\sqrt{2}}{2}}\right)·(cos\theta +1)\hfill \\ \hfill g_2(\theta ,\varphi )& ={sin}^2\left(\varphi \right)·{sin}^2\left(\varphi -{\displaystyle \frac{\pi }{2}}\right)·(cos\theta +1)\hfill \end{array}$$

(70)

This function vanishes precisely at the eight points of set $\mathcal{S}$ in (53). The component $g_1$ captures the latitude dependence, vanishing at $\theta \in \{{\displaystyle \frac{3\pi }{4}},\pi \}$, while $g_2$ encodes the longitude constraints, vanishing at the south pole and at four specific longitudes $\varphi \in \{0,{\displaystyle \frac{\pi }{2}},\pi ,{\displaystyle \frac{3\pi }{2}}\}$ when $\theta ={\displaystyle \frac{3\pi }{4}}$. The factor $(cos\theta +1)$ in $g_2$ ensures it vanishes at the south pole regardless of $\varphi$. The reason we construct this function which vanishes for the unstable points selected from Figure 9 is based on forming a topological analysis of where the stable states arise, which are states that satisfy a gradual decoherence and dephasing under the camel-like framework. These decoherence- instable states in the set (53) are indeed still valid as solutions to the GKSL equation, however, as they decohere immediately, they form a special case of states that pose specific consequences for quantum computation purposes,a s we briefly mentioned above. Therefore, isolating these states on the sphere can give us a topologically intuitive view of how states develop stability and instability to decoherence effects by the dissipation operator in the GKSL equation.

Thus, the function (70) characterizes the set $\mathcal{S}$ of the unstable states as the zero set where $f(\theta ,\varphi )=(g_1\left(\theta \right),g_2(\theta ,\varphi ))=(0,0)$, while (more) stable states correspond to points where $f(\theta ,\varphi )\ne (0,0)$. We can use this framework on the Bloch sphere to start by performing a gradient flow analysis to study the evolution of states as they approach the highly unstable states in $\mathcal{S}$ on the Bloch sphere.

3.9.1. Gradient Flow Analysis on the Bloch Sphere

The function

$$h(\theta ,\varphi )=g_1{\left(\theta \right)}^2+g_2{(\theta ,\varphi )}^2$$

where $g_1{g}_2$ are defined in (70) is constructed as a scalar potential function (or a Lyapunov-like function) on the Bloch sphere, ${\mathbb{S}}^2$ , whose global minima are precisely the five unstable points in the set $\mathcal{S}$ in (53). These minima occur at the four points on the three-quarter circle ($\theta ={\displaystyle \frac{3\pi }{4}}$) with $\varphi \in \{0,{\displaystyle \frac{\pi }{2}},\pi ,{\displaystyle \frac{3\pi }{2}}\}$, and at the south pole ($\theta =\pi$). By computing the negative gradient, $-\nabla h$ , we define a gradient dynamical system ${\displaystyle \frac{d}{dt}}(\theta ,\varphi )=-\nabla h(\theta ,\varphi ).$ The integral curves (trajectories or flow lines) of this vector field represent paths of steepest descent on the potential surface h . Plotting these trajectories visualizes the basins of attraction for each of the five special points under this gradient flow. This analysis provides topological insights into the phase portrait of the original system by partitioning the Bloch sphere according to the influence of these unstable points and their effect on solving the system using the ODE45 method. By visualizing the trajectories that follow ${\displaystyle \frac{d}{dt}}(\theta ,\varphi )=-\nabla h(\theta ,\varphi )$ on the Bloch sphere, we can identify the resulting in a gradient field structure where the set $\mathcal{S}$ of unstable points form a $C_4$ symmetry in the $\varphi$ direction is preserved by the gradient flow due to the periodic structure of $g_2$. This partitioning reveals how the Bloch sphere is divided into the two distinct regions, Northern and Southern hemisphere, where quantum states in each region are attracted to their corresponding decoherence hotspot, either the South pole, or the four concentric points on the Braiding ring $(\theta ={\displaystyle \frac{3\pi }{4}}$ (see Figure 11).

The flow lines in Figure 11 indicate that unstable points near $\theta ={\displaystyle \frac{3\pi }{4}}$ act as attractors for certain states descending from the North Pole. As shown in the figure, the gradient structure around the South Pole reveals a complex flow pattern towards the attractor at the South Pole. The deep red region near the North Pole, corresponding to the highest gradient values, shows that decoherence and dephasing drive states as a source toward the South Pole, while the gradient flow (white arrows) reveals this as the dominant evolution pathway. The gradient lines evolve globally away from the third quarter circle ($\theta ={\displaystyle \frac{3\pi }{4}}$), exhibiting bidirectional flow, toward the South Pole (dominant pathway), and back toward the North-West and North-East near four specific points on the third quarter circle. This behavior, clearly visible in Figure 11 (South pole image), reveals that the third quarter circle acts as a topological transition region, by the specific Lyapunov function. The gradients around both the South Pole and this quarter circle show bifurcating flows, creating what may be described as an open barrier that introduces disorder in the evolution of states.

These observations demonstrate a topological tendency of the solutions of the GKSL equation where the third quarter-circle at $\theta ={\displaystyle \frac{3\pi }{4}}$ acts as a transition region and the states can evolve both slightly Northward and due Southward. The system thus exhibits a topological transition containing unstable states and trajectories, with gradient-driven evolution creating complex flow patterns around these selected unstable points in (53).

It is worth noting that these critical points simply emerge from the global shape of h and do not affect the form of the gradient field itself. The ‘camel-like” geometry of the entropy thus arises solely from the interplay of the Hamiltonian and dissipative terms encoded in h, rather than from any special, localized influence of the unstable points. This implies that the gradient flow trajectories around the unstable point is real also in a physical experiment, as long as we consider the points in $\mathcal{S}$ in (53) as “special” or “different”, and hence are evaluated in h by being outside the desired stability against decoherence of states under the camel-like framework.

Finally, we note that the trajectories on the Bloch sphere represent a multitude of different transitions between pure states under an entirely unitary evolution, as they remain on the surface of the Bloch sphere. In a physical case, this would be represented by a series of diffractometers for beamed photons which preserves the purity while changing their polarization. Conclusively, in this section, we have thus constructed a Lyapunov function based on five highly unstable quantum states under the camel-like framework by their entropy plot in Figure 9, in order to design an trajectory path for the evolution of states under a completely unitary and real experiment in quantum information, so that we can predict the evolution of quantum states under the camel-like framework. It is also worth noting, that the North Pole represents the most stable point in terms of decoherence-resistance, and the South, represents the most unstable. By the gradient field we have designed here, the instability of states is thus not related to instability/stability towards decoherence, but stability and instability based on the gradient flow values that $-\nabla h$ attains on the Bloch sphere, which are useful for the mapping we have presented in Figure 11.

3.9.2. Topological Constraints on the Bloch Sphere

With the notion that the North pole acts as a source and the South pole’s as an attractor by considering the gradient field dynamics, we can form a basic basin map and analyze the evolution of the stability of the quantum states subdivided into topological basins. We develop this analysis from the quantum system’s fundamental functions given in (70). The gradient flow dynamics are derived from the potential function $h(\theta ,\varphi )=g_1{\left(\theta \right)}^2+g_2{(\theta ,\varphi )}^2$, with the evolution equations

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\theta }{dt}}& =-{\displaystyle \frac{\partial h}{\partial \theta }}\hfill \\ \hfill {\displaystyle \frac{d\varphi }{dt}}& =-{\displaystyle \frac{1}{sin\theta }}{\displaystyle \frac{\partial h}{\partial \varphi }}\hfill \end{array}$$

(71)

which we derive from the Riemann flow on the ${\mathcal{S}}^2$ sphere [37]. Following these dynamics, we classify the quantum state space into distinct basins:

$$\mathrm{Basin}(\theta ,\varphi )=\left\{\begin{array}{cc}\mathrm{BASIN}1,\hfill & \mathrm{if}\theta <{\displaystyle \frac{3\pi }{4}}\hfill \\ \mathrm{BASIN}2,\hfill & \mathrm{if}{\displaystyle \frac{3\pi }{4}}<\theta \le \pi \hfill \end{array}\right.$$

(72)

The boundaries of each basin are defined as a manifold with boundary and can thus be viewed as compact topological sets on the Bloch sphere, which can be assigned well-defined closures and where any calculation will converge within the established boundaries. We can thus define these basins topologically as:

$$\partial B=\{(\theta ,\varphi )\in {\mathbb{S}}^2\mid \forall ϵ>0,\exists ({\theta }_1,{\varphi }_1),({\theta }_2,{\varphi }_2)\in B_{ϵ}(\theta ,\varphi )$$

(73)

where $({\theta }_1,{\varphi }_1)\in B_1$ and $({\theta }_2,{\varphi }_2)\in B_2$. Here, $B_{ϵ}(\theta ,\varphi )$ represents the $ϵ$-neighborhood around the point $(\theta ,\varphi )$, and $B_1,B_2$ are distinct basins. The unstable points (in terms of decoherence) of the system given in (53) are characterized by:

$$\mathrm{Unstable}\mathrm{points}=\{(\theta ,\varphi )\in {\mathbb{S}}^2\mid g_1\left(\theta \right)=0\mathrm{and}{g}_2(\theta ,\varphi )=0\}$$

(74)

We plot the subdivision of basins on the Bloch sphere in Figure 12 based on this Riemannian flow geometry, where the topological ridge sinks toward the third quarter circle on periodic longitudes where we have the Braiding ring, composed of unstable states. This gives two basins, displayed in red and blue, where field lines move towards the South Pole from the North pole, north of $\theta ={\displaystyle \frac{3\pi }{4}}$.

In Section 3.9.4 we shall discuss the physical implications of this topological manifold. Meanwhile, we form a theorem summarizing the results from topological analysis of the solutions of the GKSL equation under the camel-like framework.

Theorem 2

(Global Convergence on the Bloch Sphere).

Consider the gradient flow $\dot{\mathbf{x}}=-\nabla h$ of the function $h(\theta ,\varphi )=g_1{\left(\theta \right)}^2+g_2{(\theta ,\varphi )}^2$ on ${\mathbb{S}}^2$, where

$$\begin{array}{cc}\hfill g_1\left(\theta \right)& =(cos\theta +\frac{\sqrt{2}}{2})(cos\theta +1),\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill g_2(\theta ,\varphi )& ={sin}^2\varphi ·{sin}^2(\varphi -\frac{\pi }{2})·(cos\theta +1).\hfill \end{array}$$

(76)

Then:

• The unstable points in terms of rapid decoherence consist of the critical points by the south pole $(\pi ,\varphi )$, and four saddle points at $\{(\frac{3\pi }{4},k\frac{\pi }{2}):k=0,1,2,3\}$. Also, the North pole is a critical point, however stable in terms of decoherence.
• The south pole is a global attractor: for all initial conditions $({\theta }_0,{\varphi }_0)\in {\mathbb{S}}^2\setminus \left\{(0,\varphi )\right\}$, the flow converges to $(\pi ,\varphi )$.

Proof. 

(1) Critical points satisfy $\nabla h=0$. Since $h=g_1^2+g_2^2$, we have $\nabla h=2g_1\nabla g_1+2g_2\nabla g_2$. At $\theta =0$: $g_1\left(0\right)\ne 0$ but $\partial g_1{/\partial \theta |}_{\theta =0}=0$, and $g_2(0,\varphi )=0$, giving isolated critical points. At $\theta =\pi$: both $g_1\left(\pi \right)=g_2(\pi ,\varphi )=0$ due to the factor $(cos\theta +1)$. At $\theta =3\pi /4$: $\partial g_1/\partial \theta =0$ and $g_2=0$ when $sin\varphi =0$ or $cos\varphi =0$, yielding four saddle points.

(2) Since $h(\pi ,\varphi )=0$ and $h(\theta ,\varphi )>0$ for all $\theta \ne \pi$, the south pole is the unique global minimum. Furthermore, along any trajectory of the gradient flow $\dot{\mathbf{x}}=-\nabla h$, the function h decreases:

$${\displaystyle \frac{d}{dt}}h\left(\mathbf{x}\left(t\right)\right)=\nabla h·\dot{\mathbf{x}}=-{|\nabla h|}^2\le 0.$$

This means that h always decreases (or stays constant) along trajectories, which can be seen in the gradient plot. Since the sphere is bounded and $h\ge 0$ everywhere, h must approach some limiting value. The only points where trajectories can stop are where $\nabla h=0$ (the critical points). At the north pole and saddle points on the third quarter circle, small perturbations cause trajectories to move away (since h decreases from these points). However, at the south pole where $h=0$ (the minimum), trajectories cannot decrease h further, so they must stop. Therefore, all trajectories eventually reach the south pole, making it the global attractor.    □

Remark 1.

It is important to note that the notion of instability described by the gradient flow differs from the instability observed in the von Neumann entropy plots in Figure 9. The gradient flow instability arises from the rate at which the Lyapunov-like function changes over the Bloch sphere, whereas the decoherence-related instability of the states in the set $\mathcal{S}$ is due to their susceptibility to the dissipation operator in the GKSL equation under the camel-like entropy framework, which we precisely observe in Figure 9. The gradient field of the Lyapunov function thus illustrates a descent from decoherence-stable states (North Pole) toward decoherence-sensitive ones (South Pole).

3.9.3. Evolution by Purity in the Bloch Ball

In this section, we investigate the evolution of pure states on the Bloch ball (where the entire space inside ${\mathbb{S}}^2$ is considered as a vector space of solutions), following the properties of the GKSL equation under the camel-like framework. The evolution starts from pure states (on the Bloch sphere - ${\mathbb{S}}^2$) to mixed states in the center of the Bloch ball (a dense globular manifold). The evolution on the Bloch ball is tracked to simulate a pair of entangled particles that are passed through a series of diffractometers and electromagnetic fields, rearranging their polarity and reducing their state of purity into mixed state. Naturally, this process is entirely irreversible due to the uncertainty principle, and evolves from the Bloch surface to the center, following the paths allowed by the GKSL equation under a camel-like framework. These paths should not be confused with the gradient paths in the ${\mathbb{S}}^2$ manifold from Figure 11, where the gradient field lines determine the evolution paths from pure metastable states to stable pure states (stability in the sense of being an admissible solution to the GKSL equation).

The results are shown in Figure 13 below.

The universal steady state, the great attractor of all states under the GKSL equation in the camel-like framework, is computed in Python and presented as the density matrix shown in Figure 13.

$${\rho }_{Universal}=\left(\begin{array}{cc}0.515041& 0.00134759+0.121413i\\ 0.00134759-0.121413i& 0.484959.\end{array}\right)$$

(77)

3.9.4. Relevance for Physics

The Lyapunov function with four singularities on an unstable equilibrium circle represents a geometric quantum control framework that actively avoids the most decoherence-prone states on the Bloch sphere. These four points at $\theta =3\pi /4$ are fully admissible pure quantum states, but they exhibit the fastest decoherence rates in the system, making them the most fragile states to maintain coherence by the notion of stability of states in quantum information experiments [38,39]. By strategically placing singularities of $g_2(\theta ,\varphi )$ at these four maximally-decohering points, the gradient flow creates a control landscape that steers quantum trajectories away from these vulnerable regions, effectively implementing a decoherence-avoiding quantum control protocol [40]. This approach differs from traditional geometric phase schemes by explicitly incorporating decoherence information into the Lyapunov function design, where the unstable manifold at $\theta =3\pi /4$ represents the locus of states with maximal environmental coupling [41]. The resulting flow pattern, with its characteristic basin structure, ensures that quantum states rapidly escape the high-decoherence region and converge to the more stable states, thereby minimizing exposure to decoherence throughout the evolution [38,40]. This framework can be experimentally realized using polarization qubits passing through engineered diffractometers, where polarization-dependent diffraction losses implement the Lyapunov function $h(\theta ,\varphi )$. Specifically, using a first diffractometer we create the $g_1\left(\theta \right)$ dependence through selective polarization, while a second holographic diffractometer implements $g_2(\theta ,\varphi )$ by creating null diffraction points (the unstable states we call “singularities”) at the four decoherence-prone states [42,43]. The measurement backaction from photons lost to higher diffraction orders drives the remaining zero-order photons along the gradient flow trajectories, effectively steering them away from the fragile states at $\theta =3\pi /4$ toward the decoherence-protected south pole [44].

4. Conclusions

This study provided novel results of simulating the GKSL-equation in modeling open quantum systems, emphasizing an analysis of entropy, entanglement entropy and its rate of change, quantum discord, mutual information, and quantum steering for Bell, product, and mixed states. The results revealed the evolution under environmental decoherence, with entropy exhibiting a characteristic “camel-like” behavior, a barrier-landscape framework for the GKSL equation developed by Khrennikov and Basileva [2]. Most importantly, these results were based on a thorough analysis of the spectrum of the GKSL operator, where the optimal value of the relaxation factor $\gamma$ was found, and implemented in the calculations for the open quantum system. We have used topology and differential geometry to study the properties of the GKSL equation under the camel-like framework, on the Bloch sphere. Beyond its theoretical significance, the findings of this study may have practical implications in fields such as quantum computing and quantum thermodynamics. The ability to identify transitions towards a classical or a quantum-like information exchange during a systems’ evolution could be useful for error correction in quantum computers or understanding dissipation in quantum heat engines. Future work could explore how the camel-like entropy behavior influences quantum circuits’ stability or thermodynamic efficiency in quantum systems.

Appendix A.1. Analytical Expressions with Initial Condition from Single Parameterization Procedure at ϕ=π

We turn now our attention to extracting valid initial conditions following the rationale in [27]. We therefore insert for $\varphi =\pi$ and obtain the initial conditions for the analytical solutions of the GKSL equation as:

$${\rho }_0=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)$$

(A1)

which is the pure state $\mid 0\rangle$.

We form the derivative of the density matrix:

$${\displaystyle \frac{d\rho }{dt}}=\left(\begin{array}{cc}{\dot{\rho }}_{11}& {\dot{\rho }}_{12}\\ {\dot{\rho }}_{21}& {\dot{\rho }}_{22}\end{array}\right)=\left(\begin{array}{cc}{\dot{\rho }}_{11}& {\dot{\rho }}_{12}\\ {\dot{\rho }}_{21}& -{\dot{\rho }}_{11}\end{array}\right)$$

(A2)

and form the GKSL equation using the GKSL operator in (49):

$$\left(\begin{array}{cc}{\dot{\rho }}_{11}& {\dot{\rho }}_{12}\\ {\dot{\rho }}_{21}& -{\dot{\rho }}_{11}\end{array}\right)=\left(\begin{array}{cc}{\displaystyle \frac{1}{2}}\left(-{\rho }_{11}+2i{\rho }_{12}-2i{\rho }_{21}+1\right)& 2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{12}-i\\ -2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{21}+i& {\displaystyle \frac{1}{2}}\left({\rho }_{11}-2i{\rho }_{12}+2i{\rho }_{21}-1\right)\end{array}\right)$$

(A3)

which encompasses the system of ODEs :

$$\begin{array}{cc}\hfill {\dot{\rho }}_{11}& ={\displaystyle \frac{1}{2}}\left(-{\rho }_{11}+2i{\rho }_{12}-2i{\rho }_{21}+1\right)\hfill \\ \hfill {\dot{\rho }}_{12}& =2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{12}-i\hfill \\ \hfill {\dot{\rho }}_{21}& =-2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{21}+i,\hfill \end{array}$$

(A4)

which we solve computationally with the SymbPy method.

Analytical Solutions

The system was solved symbolically yielding:

$$\begin{array}{cc}\hfill {\rho }_{11}\left(t\right)& ={\displaystyle \frac{1}{2805}}\left(1445e^{{\displaystyle \frac{3t}{4}}}+8\left(85+\sqrt{255}i\right)e^{{\displaystyle \frac{t\left(3-\sqrt{255}i\right)}{8}}}+8\left(85-\sqrt{255}i\right)e^{{\displaystyle \frac{t\left(3+\sqrt{255}i\right)}{8}}}\right)e^{-{\displaystyle \frac{3t}{4}}}\hfill \\ \hfill {\rho }_{12}\left(t\right)& ={\displaystyle \frac{2}{2805}}\left(170i{e}^{{\displaystyle \frac{3t}{4}}}-\left(21\sqrt{255}+85i\right)e^{{\displaystyle \frac{t\left(3-\sqrt{255}i\right)}{8}}}+\left(21\sqrt{255}-85i\right)e^{{\displaystyle \frac{t\left(3+\sqrt{255}i\right)}{8}}}\right)e^{-{\displaystyle \frac{3t}{4}}}\hfill \\ \hfill {\rho }_{21}\left(t\right)& ={\displaystyle \frac{2}{2805}}\left(-170i{e}^{{\displaystyle \frac{3t}{4}}}+\left(21\sqrt{255}+85i\right)e^{{\displaystyle \frac{t\left(3-\sqrt{255}i\right)}{8}}}+\left(-21\sqrt{255}+85i\right)e^{{\displaystyle \frac{t\left(3+\sqrt{255}i\right)}{8}}}\right)e^{-{\displaystyle \frac{3t}{4}}}\hfill \end{array}$$

(A5)

Appendix A.1.1. Properties of the Single-Parameter Analytical Solution

The density matrix components ${\rho }_{11}\left(t\right)$, ${\rho }_{12}\left(t\right)$, and ${\rho }_{21}\left(t\right)$ in (A5) constitutes a physically valid description of the quantum system’s evolution. The solution exhibits three key features that ensure its physical consistency. The Hermitian nature of the density matrix is preserved exactly through the conjugate relationship ${\rho }_{12}\left(t\right)={\rho }_{21}^{*}\left(t\right)$, which was verified computationally with zero residual difference. Furthermore, the reality condition for diagonal elements holds strictly, as $\Im \left({\rho }_{11}\left(t\right)\right)\equiv 0$ at all times. This is directly verified from the analytical expression in (A5) where the diagonal element ${\rho }_{11}\left(t\right)$ is guaranteed to be real-valued for all $t\ge 0$. This reality condition can be verified by analyzing the complex conjugate symmetry of the exponentials. Using Euler’s identity, the imaginary terms cancel exactly:

$$\begin{array}{cc}\hfill e^{-{\displaystyle \frac{t(3\pm \sqrt{255}i)}{8}}}& =e^{-3t/8}\left[cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\mp isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right]\hfill \end{array}$$

(A6)

Substituting into ${\rho }_{11}\left(t\right)$ yields:

$$\begin{array}{cc}\hfill {\rho }_{11}\left(t\right)& ={\displaystyle \frac{17}{33}}+{\displaystyle \frac{8e^{-3t/8}}{2805}}[85\left(e^{i\sqrt{255}t/8}+e^{-i\sqrt{255}t/8}\right)\hfill \\ \hfill & +\sqrt{255}i\left(e^{i\sqrt{255}t/8}-e^{-i\sqrt{255}t/8}\right)]\hfill \\ \hfill & ={\displaystyle \frac{17}{33}}+{\displaystyle \frac{16e^{-3t/8}}{2805}}\left[85cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)-\sqrt{255}sin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right]\hfill \end{array}$$

(A7)

The final form contains only real-valued trigonometric functions, confirming $\Im \left({\rho }_{11}\left(t\right)\right)\equiv 0$. Moreover, the time dependence of (A5) exhibits underdamped oscillations of the form $e^{-t(3\pm i\sqrt{255})/8}$, with decay rate $\gamma =3/8$ and oscillation frequency $\omega =\sqrt{255}/8$. In the long-time limit, these terms vanish, leaving a calculated steady-state solution:

$${\rho }_{\infty }=\left(\begin{array}{cc}{\displaystyle \frac{17}{33}}& {\displaystyle \frac{4i}{33}}\\ -{\displaystyle \frac{4i}{33}}& {\displaystyle \frac{16}{33}}\end{array}\right)$$

(A8)

This steady state is the exact form of the formula given in (10), respective for the analytic solution considered by (A5). This steady state density matrix is Hermitian, has unit trace, and is positive semi-definite, preserving all quantum mechanical consistency conditions.

We now plot the von Neumann entropy of the solution in (A5).

Figure A1. von Neumann entropy evolution of the analytical solution in (A5) and numerical solutions to the GKSL equation in (41) under the camel-like framework, with initial conditions in (A1).

Figure A1. von Neumann entropy evolution of the analytical solution in (A5) and numerical solutions to the GKSL equation in (41) under the camel-like framework, with initial conditions in (A1).

The damped oscillations in the von Neumann entropy (Figure A1) show a qualitative resemblance to both the three-hump analytical solutions of the GKSL equation derived from single-parameter approach (see equation (A5)). The plot in Figure A1 shows familiarities with the Rabi oscillations reported by Pla et al. [23]. However, the dissipative nature of our GKSL model implies that a direct comparison to the purely coherent Rabi oscillations in [23] requires careful consideration and is better termed as “damped Rabi oscillations”. This model may thus provide insights into spin-fraction evolution under damped electromagnetic perturbation in single atoms, but further analysis is warranted experimentally.

Appendix A.2. Analytical Expressions of a Pure and Stable Point on the Bloch Sphere: $(\theta ,\varphi )=({\displaystyle \frac{\pi }{4}},{\displaystyle \frac{3\pi }{2}})$ as Initial Condition

We now seek to obtain an analytical solution using an admissible point, $(\theta ,\varphi )=({\displaystyle \frac{\pi }{4}},{\displaystyle \frac{3\pi }{2}})$, which lies on the quarter circle on the Northern hemisphere of the Bloch sphere, between singularities (see Section 3.9). This gives the following state:

$$\begin{array}{cc}\hfill |\psi \left({\displaystyle \frac{\pi }{4}},{\displaystyle \frac{3\pi }{2}}\right)\rangle & =cos\left({\displaystyle \frac{\pi /4}{2}}\right)|0\rangle +e^{i(3\pi /2)}sin\left({\displaystyle \frac{\pi /4}{2}}\right)|1\rangle \hfill \\ \hfill & =cos\left({\displaystyle \frac{\pi }{8}}\right)|0\rangle +e^{i(3\pi /2)}sin\left({\displaystyle \frac{\pi }{8}}\right)|1\rangle \hfill \\ \hfill & =cos\left({\displaystyle \frac{\pi }{8}}\right)|0\rangle +(-i)sin\left({\displaystyle \frac{\pi }{8}}\right)|1\rangle \hfill \\ \hfill & ={\displaystyle \frac{\sqrt{2+\sqrt{2}}}{2}}|0\rangle -i{\displaystyle \frac{\sqrt{2-\sqrt{2}}}{2}}|1\rangle ,\hfill \end{array}$$

(A9)

which is in density matrix form:

$${\rho }_0=\left(\begin{array}{cc}{\displaystyle \frac{2+\sqrt{2}}{4}}& i{\displaystyle \frac{\sqrt{2}}{4}}\\ -i{\displaystyle \frac{\sqrt{2}}{4}}& {\displaystyle \frac{2-\sqrt{2}}{4}}\end{array}\right)$$

(A10)

We have system of ODEs from above which we solve with MATLAB, using ${\rho }_0$ as initial condition:

$$\begin{array}{cc}\hfill {\dot{\rho }}_{11}& ={\displaystyle \frac{1}{2}}\left(-{\rho }_{11}+2i{\rho }_{12}-2i{\rho }_{21}+1\right)\hfill \\ \hfill {\dot{\rho }}_{12}& =2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{12}-i\hfill \\ \hfill {\dot{\rho }}_{21}& =-2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{21}+i.\hfill \end{array}$$

(A11)

The final analytical expression results as:

$$\begin{array}{cc}\hfill {\rho }_{11}\left(t\right)& =e^{-{\displaystyle \frac{3t}{8}}}\left({\displaystyle \frac{\sqrt{255}}{16}}-{\displaystyle \frac{i}{16}}\right)\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)-isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \\ \hfill & ·\left({\displaystyle \frac{\sqrt{2}\sqrt{255}}{2040}}+e^{{\displaystyle \frac{3t}{8}}+{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}-{\displaystyle \frac{2i}{33}}\right)-{\displaystyle \frac{i\sqrt{2}}{8}}+{\displaystyle \frac{16388992192603249\sqrt{255}}{1579074619346780160}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & +e^{-{\displaystyle \frac{3t}{8}}}\left({\displaystyle \frac{\sqrt{255}}{16}}+{\displaystyle \frac{i}{16}}\right)(e^{{\displaystyle \frac{3t}{8}}-{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\sqrt{255}}{4737223858040340480}}\left(i\sqrt{2}\sqrt{255}·2322168557862912+\sqrt{2}·2322168557862912\right.\hfill \\ \hfill & \left.-i·\sqrt{255}·1125899906842624+49166976577809747\right))\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)+isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \\ \hfill {\rho }_{12}\left(t\right)& =-e^{-{\displaystyle \frac{3t}{8}}}\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)-isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \\ \hfill & ·\left({\displaystyle \frac{\sqrt{2}\sqrt{255}}{2040}}+e^{{\displaystyle \frac{3t}{8}}+{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}-{\displaystyle \frac{2i}{33}}\right)-{\displaystyle \frac{i\sqrt{2}}{8}}+{\displaystyle \frac{16388992192603249\sqrt{255}}{1579074619346780160}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & +e^{-{\displaystyle \frac{3t}{8}}}(e^{{\displaystyle \frac{3t}{8}}-{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\sqrt{255}}{4737223858040340480}}\left(i\sqrt{2}\sqrt{255}·2322168557862912+\sqrt{2}·2322168557862912\right.\hfill \\ \hfill & \left.-i·\sqrt{255}·1125899906842624+49166976577809747\right))\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)+isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \\ \hfill {\rho }_{21}\left(t\right)& =e^{-{\displaystyle \frac{3t}{8}}}\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)-isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \\ \hfill & ·\left({\displaystyle \frac{\sqrt{2}\sqrt{255}}{2040}}+e^{{\displaystyle \frac{3t}{8}}+{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}-{\displaystyle \frac{2i}{33}}\right)-{\displaystyle \frac{i\sqrt{2}}{8}}+{\displaystyle \frac{16388992192603249\sqrt{255}}{1579074619346780160}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & -e^{-{\displaystyle \frac{3t}{8}}}(e^{{\displaystyle \frac{3t}{8}}-{\displaystyle \frac{i\sqrt{255}t}{8}}}\left({\displaystyle \frac{46\sqrt{255}}{2805}}+{\displaystyle \frac{2i}{33}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\sqrt{255}}{4737223858040340480}}\left(i\sqrt{2}\sqrt{255}·2322168557862912+\sqrt{2}·2322168557862912\right.\hfill \\ \hfill & \left.-i·\sqrt{255}·1125899906842624+49166976577809747\right))\left(cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)+isin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\right)\hfill \end{array}$$

(A12)

Appendix A.2.1. Properties of the Analytical Solution on the Bloch Sphere

The analytical solution given above maintains exact Hermiticity through the conjugate relationship ${\rho }_{21}\left(t\right)={\rho }_{12}^{*}\left(t\right)$, with all imaginary components in the transient terms $e^{-t(3\pm i\sqrt{255})/8}$ canceling exactly (data not shown for limitation purposes). Trace preservation emerges naturally via ${\rho }_{22}\left(t\right)=1-{\rho }_{11}\left(t\right)$, where the diagonal elements’ time dependence guarantees unit trace for all $t\ge 0$. Moreover, the solution approaches the non-equilibrium steady state

$$\underset{t\to \infty }{lim}\rho \left(t\right)={\rho }_{steady}=\left(\begin{array}{cc}{\displaystyle \frac{17}{33}}& {\displaystyle \frac{4i}{33}}\\ -{\displaystyle \frac{4i}{33}}& {\displaystyle \frac{16}{33}}\end{array}\right),$$

(A13)

which yields a preservation of quantum coherence through non-vanishing off-diagonal elements. This reflects the balance between the unitary evolution frequency $\sqrt{255}/8$ and the dissipative decay rate $\gamma =3/8$ (the same as for the solution in (A5)). The transient dynamics is governed by oscillatory terms $cos(\sqrt{255}t/8)$ and $sin(\sqrt{255}t/8)$, exponentially damped by $e^{-3t/8}$, while the initial coherence $i\sqrt{2}/4$ determines both the transient behavior and the properties of the steady state. We note that the steady state density matrix given in (A13) is Hermitian and has unit trace and it is positive semi-definite, hence satisfying the necessary requirements for a solution in open quantum systems.

The steady state in (A13) is now visualized on the Bloch sphere (Figure A2). The plot shows that this is a mixed state of very low purity, being close to the center of the Bloch sphere (pure states are on the surface of the sphere).

Figure A2. Steady state in (A13) of the analytical solution on the Bloch sphere (A12), derived by solving the system of equations (49) with the initial conditions from the Bloch parametrization in (A10). Final steady-state bloch state-vector: ${\overrightarrow{\rho }}_{\mathrm{steady}}=[0,0.24,0.03]$.

Figure A2. Steady state in (A13) of the analytical solution on the Bloch sphere (A12), derived by solving the system of equations (49) with the initial conditions from the Bloch parametrization in (A10). Final steady-state bloch state-vector: ${\overrightarrow{\rho }}_{\mathrm{steady}}=[0,0.24,0.03]$.

The Bloch vector ${\overrightarrow{\rho }}_{steady}=[0,0.24,0.03]$ (Figure A2) represents a quantum state with purity $\mathrm{Tr}\left({\rho }^2\right)={\displaystyle \frac{577}{1089}}\approx 0.53$, exhibiting partial coherence, displaying slightly unequal populations in computational basis states (${\rho }_{00}={\displaystyle \frac{17}{33}}\approx 0.515$ and ${\rho }_{11}={\displaystyle \frac{16}{33}}\approx 0.485$). The non-zero y-component with magnitude ${\displaystyle \frac{8}{33}}\approx 0.242$ indicates moderate quantum superposition, while its location well inside the Bloch sphere (Bloch vector length $\approx 0.244$) confirms it cannot be represented by a pure-state vector. This state exhibits significant mixedness with modest coherence, making it relevant for quantum information theory as it represents a typical steady state that might arise in open quantum systems subject to decoherence and dissipation processes. The density matrix for this steady state displays a similar von Neumann entropy of the analytical solution in (A12) as the previous state in (A5).

Figure A3. von Neumann entropy evolution of the analytical solution in (A12) to the GKSL equation in (41) under the camel-like framework, with the initial conditions in (A10).

Figure A3. von Neumann entropy evolution of the analytical solution in (A12) to the GKSL equation in (41) under the camel-like framework, with the initial conditions in (A10).

The entropy evolution in Figure A3 shows oscillations similar to those found in the single-parametrization approach (Section A.1). While these oscillations resemble Rabi oscillations, they differ crucially in their decay to a steady state, a hallmark of open quantum systems, unlike the undamped Rabi oscillations reported by Pla et al. [23]. This damped oscillatory behavior finds parallels in other physical systems, such as neutrino oscillations (exhibiting decaying sinusoidal patterns due to interactions and detection limitations) and neutral kaon ($K^{\theta }$) oscillations (where damping arises from kaon decay) [45].

Appendix A.3. Analytical Expressions of a Metastable Point by Bloch Parameter (θ,ϕ)=(0.05,0) as Initial Condition

In this final section, we shall derive the analytic expression of a metastable point belonging to the $B_N$ basin in the Bloch sphere, reported in the topological analysis section. The quantum state corresponding to $(\theta ,\varphi )=(0.05,0)$ is

$$|\psi \rangle =cos\left(0.05\right)|0\rangle +sin\left(0.05\right)e^{i·0}|1\rangle \approx 0.999|0\rangle +0.05|1\rangle$$

(A14)

which we calculate in density matrix form to:

$$\rho =|\psi \rangle \langle \psi |=\left(\begin{array}{cc}{cos}^2\left(0.05\right)& cos\left(0.05\right)sin\left(0.05\right)\\ sin\left(0.05\right)cos\left(0.05\right)& {sin}^2\left(0.05\right)\end{array}\right)\approx \left(\begin{array}{cc}0.998& 0.05\\ 0.05& 0.0025\end{array}\right)$$

(A15)

We then use the same system of ODEs as for the other initial states,

$$\begin{array}{cc}\hfill {\dot{\rho }}_{11}& ={\displaystyle \frac{1}{2}}\left(-{\rho }_{11}+2i{\rho }_{12}-2i{\rho }_{21}+1\right)\hfill \\ \hfill {\dot{\rho }}_{12}& =2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{12}-i\hfill \\ \hfill {\dot{\rho }}_{21}& =-2i{\rho }_{11}-{\displaystyle \frac{1}{4}}{\rho }_{21}+i,\hfill \end{array}$$

(A16)

which we solve using MATLAB code with (A15) as initial conditions.

We obtain:

$$\begin{array}{cc}\hfill {\rho }_{11}\left(t\right)& ={\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left({\displaystyle \frac{7967}{8250}}+{\displaystyle \frac{1600}{33}}\mathrm{i}\right)\hfill \\ \hfill & +\sqrt{255}{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}sin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left({\displaystyle \frac{8011}{701250}}+{\displaystyle \frac{320}{561}}\mathrm{i}\right)\hfill \\ \hfill & +{\displaystyle \frac{1}{33}}-{\displaystyle \frac{1600}{33}}\mathrm{i}\hfill \\ \hfill {\rho }_{12}\left(t\right)& ={\displaystyle \frac{{\mathrm{e}}^{-\frac{t}{4}}}{20}}+{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left({\displaystyle \frac{400}{33}}-{\displaystyle \frac{8}{33}}\mathrm{i}\right)\hfill \\ \hfill & +\sqrt{255}{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}sin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left(-{\displaystyle \frac{560}{187}}+{\displaystyle \frac{20912}{350625}}\mathrm{i}\right)\hfill \\ \hfill & -{\displaystyle \frac{400}{33}}+{\displaystyle \frac{8}{33}}\mathrm{i}\hfill \\ \hfill {\rho }_{21}\left(t\right)& ={\displaystyle \frac{{\mathrm{e}}^{-\frac{t}{4}}}{20}}+{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left(-{\displaystyle \frac{400}{33}}+{\displaystyle \frac{8}{33}}\mathrm{i}\right)\hfill \\ \hfill & +\sqrt{255}{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}sin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left({\displaystyle \frac{560}{187}}-{\displaystyle \frac{20912}{350625}}\mathrm{i}\right)\hfill \\ \hfill & +{\displaystyle \frac{400}{33}}-{\displaystyle \frac{8}{33}}\mathrm{i}\hfill \\ \hfill {\rho }_{22}\left(t\right)& ={\displaystyle \frac{32}{33}}+\sqrt{255}{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}sin\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left(-{\displaystyle \frac{8011}{701250}}-{\displaystyle \frac{320}{561}}\mathrm{i}\right)\hfill \\ \hfill & +{\mathrm{e}}^{-{\displaystyle \frac{3t}{8}}}cos\left({\displaystyle \frac{\sqrt{255}t}{8}}\right)\left(-{\displaystyle \frac{7967}{8250}}-{\displaystyle \frac{1600}{33}}\mathrm{i}\right)\hfill \\ \hfill & +{\displaystyle \frac{1600}{33}}\mathrm{i}\hfill \end{array}$$

(A17)

We analyse the entropy evolution of this analytical solution in Figure A4.

The von Neumann entropy plot reveals the evolution of our quantum state under the specified dynamics, starting from a nearly pure state ($S\approx 0.016$) in basin $B_N$, reaching maximum mixedness ($S\approx 0.877$) at $t\approx 2.98$, and eventually returning to purity. This behavior reflects the specific trajectory determined by our differential equations, showing substantial mixing before eventual purification. States initialized in $B_N$ exhibit prolonged dwelling times before eventually finding escape channels toward their final attractors, making them physically significant for quantum information processing. This topological obstruction creates effective barriers in the quantum state space, allowing for the existence of these metastable states that remain coherent for extended periods despite being energetically unfavorable, a property that emerges naturally from the gradient flow structure we analyzed.

The model of the metastable state has a entropy behavior that is experiencing something akin to quantum recurrences or resonances, but with very precise returns to pure states in between periods of mixing. Such behavior would be extremely difficult to realize in natural physical systems, as it requires highly controlled and isolated quantum dynamics that resist the typical irreversible decoherence processes, therefore we can conclude that metastable state in the respective topological basin do not reflect sound physical phenomena.

Figure A4. Evolution of the von Neumann entropy for a quantum state initialized in the stable basin $B_N$. The system begins in a nearly pure state (entropy $\approx 0.016$), reaches maximum mixedness (entropy $\approx 0.877$) at $t\approx 2.98$, and eventually returns to purity. This behavior demonstrates the characteristic dynamics of states in the $B_N$ basin, showing temporary confinement followed by escape to the final attractor. The non-monotonic entropy evolution reflects the topological structure of the gradient flow defined by the GKSL master equation.

Figure A4. Evolution of the von Neumann entropy for a quantum state initialized in the stable basin $B_N$. The system begins in a nearly pure state (entropy $\approx 0.016$), reaches maximum mixedness (entropy $\approx 0.877$) at $t\approx 2.98$, and eventually returns to purity. This behavior demonstrates the characteristic dynamics of states in the $B_N$ basin, showing temporary confinement followed by escape to the final attractor. The non-monotonic entropy evolution reflects the topological structure of the gradient flow defined by the GKSL master equation.