3.1.1. Coupling to a Thermal Reservoir

Let consider the spin-S systems interacting with a thermal reservoir modeled by an infinite chain of quantum harmonic oscillators with *ωk*, *bk*, and *b*† *<sup>k</sup>* being, respectively, the frequency, annihilation, and creation operators for the *k*-th oscillator. The total Hamiltonian of the system is given by

$$H = \omega\_0 \mathbf{S}\_z + \sum\_k \omega\_k b\_k^\dagger b\_k + \sum \mathbf{S}\_z (\mathcal{g}\_k b\_k^\dagger + \mathcal{g}\_k^\* b\_k)\_\prime \tag{13}$$

in which *ω*<sup>0</sup> denote the transition frequency between any neighboring energy states of the spin, and *Sz*, the *z* component of spin operator, can be represented by a diagonal matrix *Sz* = diag[*s*,*s* − 1, ... , −*s*] in the eigen-basis {|*i*, *i* = *s*, ... , −*s*}. In the interaction picture Equation (13) into is expressed as

$$H\_I = \sum S\_z \left( g\_k b\_k^\dagger e^{i\omega\_k t} + g\_k^\* b\_k e^{-i\omega\_k t} \right), \tag{14}$$

where *gk* denotes the coupling strength between the spin and the environment through the dephasing interaction. Up to an overall phase factor, the corresponding unitary propagator is obtained as

$$V(t) = \exp\left[\frac{1}{2} S\_z \sum\_k \left(a\_k b\_k^\dagger - a^\* b\_k\right)\right],\tag{15}$$

where *α<sup>k</sup>* = 2*gk* <sup>1</sup> <sup>−</sup> *<sup>e</sup>iω<sup>k</sup> <sup>t</sup>* /*ωk*.

It is assumed that the initial state of the spin-bath system is in a product state *ρT*(0) = *ρ*(0) ⊗ *ρ<sup>B</sup>* in which *ρ*(0) denotes the initial state of spin, and

$$
\rho\_B = \frac{1}{Z\_B} e^{-\beta \sum\_k \omega\_k b\_k^\dagger b\_k} \tag{16}
$$

represents the thermal equilibrium state of the bath with partition function *ZB* and inverse temperature *β* = <sup>1</sup> *kBT* . The evolved state of the system can be calculated by [83]

$$\rho\_{nm}(t) = \rho\_{nm}(0) \exp\left[-(n-m)^2 \Gamma(t)\right],\tag{17}$$

where *n*, *m* = −*s*, −*s* + 1, ... , 0, ... ,*s* − 1,*s* and, in the continuum-mode limit, the decoherence function is given by

$$
\Gamma(t) = \int\_0^\infty f(\omega) \coth\left(\frac{\omega}{2k\_b T}\right) \frac{1 - \cos(\omega t)}{\omega^2} d\omega,\tag{18}
$$

with spectral density *J*(*ω*) = ∑*<sup>k</sup>* |*gk*| <sup>2</sup>*δ*(*<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>k*).

The Γ(*t*) behavior closely depends on the characteristics of the environment. Here we consider the Ohmic-like reservoirs with spectral density

$$J(\omega) = a \frac{\omega^{\circ}}{\omega\_{\circ}^{\circ - 1}} \exp\left(\frac{-\omega}{\omega\_{\circ}}\right),\tag{19}$$

where *α* represents a dimensionless coupling strength, and *ω<sup>c</sup>* denotes the cutoff frequency of the bath. Changing the Ohmic parameter *s*, one can obtain sub-Ohmic (0 < *s* < 1), Ohmic (*s* = 1) and super-Ohmic (*s* > 1) reservoirs.

#### 3.1.2. Coupling to a Squeezed Vacuum Reservoir

In the case that the spin system is coupled to a squeezed vacuum reservoir, the reduced density-matrix elements are similar to the ones presented in Equation (17) when the decoherence function Γ(*t*) is replaced by

$$\gamma(t) = \int\_0^\infty f(\omega) \frac{(1 - \cos(\omega t))}{\omega^2} [\cosh(2r) - \sinh(2r)\cos(\omega t - \theta)] d\omega,\tag{20}$$

where *r* is the squeezed amplitude parameter, and *θ* denotes the squeezed angle.

Because the structures of the density matrices are the same in both scenarios (coupling to thermal and squeezed vacuum reservoirs), we only focus on the interaction of the system with the squeezed vacuum reservoir, noting that the general results also holds for the thermal reservoir.

We take the qudit in the pure initial state

$$|\psi\rangle = \frac{1}{\sqrt{2s+1}} (e^{i\phi}|s\rangle + |s-1\rangle + |s-2\rangle + \cdots + |-s\rangle),\tag{21}$$

which leads to the evolved state *ρ*(*t*) given by

$$\rho(t) = \frac{1}{2s+1} \begin{pmatrix} 1 & e^{-\gamma(t)}e^{i\phi} & \cdots & e^{-(2s)^2\gamma(t)}e^{i\phi} \\ e^{-\gamma(t)}e^{-i\phi} & 1 & \cdots & e^{-(2s-1)^2\gamma(t)} \\ e^{-4\gamma(t)}e^{-i\phi} & e^{-\gamma(t)} & \cdots & e^{-(2s-2)^2\gamma(t)} \\ \vdots & 1 & \ddots & \\ e^{-(2s)^2\gamma(t)}e^{-i\phi} & e^{-(2s-1)^2\gamma(t)} & \cdots & 1 \end{pmatrix}. \tag{22}$$

Therefore, the time derivative of the HSS-based witness is obtained as

$$\chi(t) = -\frac{1}{2s+1} \frac{\partial \gamma(t)}{\partial t} \frac{\sum\_{k=1}^{2s} k^2 e^{-2k^2 \gamma(t)}}{\sum\_{k=1}^{2s} e^{-2k^2 \gamma(t)}}. \tag{23}$$

The HSS-based witness *χ*(*t*) > 0 tells us that the process is non-Markovian whenever *∂γ*(*t*) *<sup>∂</sup><sup>t</sup>* < 0, which corresponds to time intervals in which the decoherence function decreases, leading to the re-coherence phenomenon. As known, in this system the non-Markovian effects, originating from the non-divisible maps, appear when the decoherence function temporarily decays with time [84]. Therefore, our witness correctly predicts the intervals at which the memory effects arise in this single-qudit system. Moreover, when *γ*(*t*) is a monotonous increasing function of time, the dynamics is Markovian because the coherence decays monotonously with time.

#### *3.2. Hybrid Qubit–Qutrit System Interacting with Various Quantum and Classical Environments*

The composite hybrid qubit(*A*)–qutrit(*B*) system consists of a spin–<sup>1</sup> <sup>2</sup> subsystem (qubit A) and a spin-1 subsystem (qutrit B). In the following, we study the interaction of this composite system with local non-Markovian environments *A* and *B*, or with a common environment *C* modeling quantum or classical noises. The theoretical schematic of this system is depicted in Figure 1.

#### 3.2.1. Coupling to Independent Squeezed Vacuum Reservoirs

Now we investigate the scenario in which each of the subsystems, i.e., the qubit *A* (*sA* = <sup>1</sup> <sup>2</sup> ) and qutrit *B* (*sB* = 1), interacts independently with its local squeezed vacuum reservoir. For simplicity we assume that the characteristics of the reservoirs are similar. Equation (17), with the decoherence factor introduced in Equation (20), gives the reduced density matrices of the subsystems. Computing them and applying the method presented in [81], one can obtain the elements of the evolved density matrix of the composite system as [85]

$$\rho\_{AB\_{nm}}(t) = \rho\_{AB\_{nm}}(0) \exp\left[-(n\_A - m\_A)^2 - (n\_B - m\_B)^2\right] \gamma(t),\tag{24}$$

**Figure 1.** Illustration of the composite qubit(*A*)-qutrit(*B*) system; Blue dashed lines represent entanglement between the subsystems. The bipartite system can interact either with independent local environments *EA*, *EB* or with a common environment *EC*.

where *nA*, *mA* = −*sA*,...,*sA* and *nB*, *mB* = −*sB*,...,*sB*.

**Pure initial state.** We take the hybrid qubit–qutrit system initially in a pure state given by [61]

$$|\psi\rangle = \frac{1}{\sqrt{6}} \left( \epsilon^{i\phi} |00\rangle + |01\rangle + |02\rangle + |10\rangle + |11\rangle + |12\rangle \right),\tag{25}$$

which leads to a dynamics of the system described by the evolved reduced density matrix *ρ*(*t*) whose elements are presented in Appendix A.1. Then, the HSS is obtained as

$$H\text{SS} = \frac{1}{6}\sqrt{2e^{-2\gamma(t)} + e^{-4\gamma(t)} + e^{-8\gamma(t)} + e^{-10\gamma(t)}}.\tag{26}$$

The dynamics of negativity, MID and HSS computed by the evolved state of the system are plotted in Figure 2. We find that each of the measures initially decreases with time, then starts to increase, and finally remains approximately constant over time, a behavior known as the freezing phenomenon [86–92]. As discussed, the revival of the quantum correlation measures can be attributed to the non-Markovian evolution of the system [47]. We see that the behaviors of the HSS, negativity and quantum correlation exhibit an excellent qualitative agreement. Consequently, the HSS-based witness can precisely capture the non-Markovian dynamics of the composite system.

**Figure 2.** Evolution of the negativity, MID and HSS as a function of dimensionless time *τ* = *ω*0*t* when each subsystem of the hybrid qubit–qutrit system, starting from the initial pure state, is independently subject to a squeezed vacuum reservoir. The values of the other parameters are *α* = 0.1, *ω<sup>c</sup>* = 20*ω*0, *r* = 0.3, *φ* = *π* and *s* = 3.

**Mixed initial state.** The non-Markovianity of the system, as faithfully individuated by quantum correlation measures, may in general depend on the initial state. It is thus important to investigate whether the HSS witness, obtained from the initial pure state of Equation (25) by definition, is capable to identify the non-Markovian character of the system dynamics also when the system starts from a mixed state. We shall study this aspect here and in all the other environmental conditions considered hereafter (see sections below devoted to a mixed initial state).

We consider the one-parameter mixed entangled state as the initial state of the hybrid qubit–qutrit system [93]

$$\rho\_0(p) = \frac{p}{2}(|01\rangle\langle 01| + |11\rangle\langle 11|) + p|\psi^+\rangle\langle \psi^+| + (1 - 2p)|\psi^-\rangle\langle \psi^-|,\tag{27}$$

where

$$\begin{aligned} \vert \psi^+ \rangle &= \frac{1}{\sqrt{2}} (\vert 00 \rangle + \vert 12 \rangle), \\ \vert \psi^- \rangle &= \frac{1}{\sqrt{2}} (\vert 02 \rangle + \vert 10 \rangle), \end{aligned} \tag{28}$$

in which the entanglement parameter *p* varies from 0 to 1 such that *ρ*(*p*) is entangled except for *p* = <sup>1</sup> <sup>3</sup> . We point out that such a state is taken as the initial state of the system for the dynamics of the quantum correlation quantifiers, namely negativity and MID. We find that Equation (27) leads to the evolved state of the system

$$
\rho(t) = \begin{pmatrix}
\frac{p}{2} & 0 & 0 & 0 & 0 & \frac{p}{2}\mathcal{F} \\
0 & \frac{p}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2}\mathcal{F} & 0 & 0 \\
0 & 0 & \frac{1-2p}{2}\mathcal{F} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{p}{2} & 0 \\
\frac{p}{2}\mathcal{F} & 0 & 0 & 0 & 0 & \frac{p}{2}
\end{pmatrix},\tag{29}
$$

where <sup>F</sup> <sup>=</sup> *<sup>e</sup>*−5*γ*(*t*). Then, the negativity is given by [71]

$$\begin{split} \mathcal{N} &= \frac{(p-1)}{2} + \frac{1}{4}|p + (1-p)\mathcal{F}| + \frac{1}{4}|p - (1-p)\mathcal{F}| + \frac{1}{4}|p - (1-2p)\mathcal{F}| + \\ &\frac{1}{4}|p + (1-2p)\mathcal{F}|. \end{split} \tag{30}$$

Moreover, using Equation (10) we can compute the MID as

$$\mathcal{M} = \frac{(1-p)}{2} [(1+\mathcal{F})\log\left(1+\mathcal{F}\right) + (1-\mathcal{F})\log\left(1-\mathcal{F}\right)].\tag{31}$$

In Figure 3, we compare the evolution of HSS, obtained from the initial pure state of Equation (25), with the dynamics of negativity and MID, computed for the mixed initial state of Equation (15), for different values of *p*. The dynamics of the HSS is again in perfect agreement with that observed for the entanglement and quantum correlations as quantified by the negativity and MID, respectively. Therefore, the HSS-based witness, computed versus the phase parameter encoded into an initial pure state of the system, can efficiently detect the non-Markovian dynamics even in the case when the initial state of our highdimensional system is not pure. It should be noted that in the presence of sudden death of entanglement, which occurs for some values of the entanglement parameter (for example, for *p* = 0.4), only the HSS and MID show the same dynamics. Hence, the negativity cannot be used as a faithful witness of non-Markovianity when it exhibits the sudden death phenomenon.

In the case of initially entangled noninteracting qubits in independent non-Markovian quantum environments, entanglement or quantum correlation revivals can be explained in terms of transfer of correlations back and forth from the composite system to the various parts of the total system. This is due to the back-action via the environment on the system, which creates correlations between qubits and environments and between the environments themselves. Accordingly, in this case the non-Markovianity is defined as backflow of information from the environment(s) to the system(s).

**Figure 3.** Comparing the evolution of negativity and MID computed for the initial mixed state of the hybrid qubit–qutrit system, when each subsystem is independently coupled to a squeezed vacuum reservoir, with HSS (obtained from the initial pure state) for different values of the entanglement parameter *p*. In all plots the remaining parameters are *α* = 0.1, *s* = 3, *ω<sup>c</sup>* = 20 *ω*0, *r* = 0.3.

#### 3.2.2. Coupling to Classical Environments

Here we assume that the hybrid qubit-qutrit system is affected by a classical environment implemented by random telegraph noise (RTN) with a Lorentzian spectrum. It is a famous class of non-Gaussian noises used to generate the low-frequency <sup>1</sup> *<sup>f</sup> <sup>α</sup>* noise both theoretically and experimentally. It is also responsible for coherent dynamics in quantum solid-state nanodevices [94–96]. Physically, the RTN may result from one of the following scenarios: (i) charges flipping between two locations in space (charge noise); (ii) electrons trapping in shallow subgap formed at the boundary between a superconductor and an insulator (noise of critical current); and (iii) spin diffusion on a superconductor surface generated by the exchange mediated by the conduction electrons (flux noise) [97,98]. The Hamiltonian of the qubit–qutrit system under the RTN is given by

$$\begin{aligned} \mathcal{H}(t) &= \mathcal{H}\_0 + \mathcal{H}\_I \\ \mathcal{H}\_0 &= \sum\_{k=A,B} \epsilon\_k S\_k^Z, \mathcal{H}\_I = \sum\_{k=A,B} \left[ f\_k L\_k(t) + f\_c \mathcal{C}(t) \right] S\_{z'}^k \end{aligned} \tag{32}$$

where *<sup>k</sup>* denote the energy of an isolated qubit (qutrit), *S<sup>A</sup> <sup>z</sup>* = *σ<sup>z</sup>* and *S<sup>B</sup> <sup>z</sup>* represent the spin operators of, respectively, the qubit and the qutrit in the *z*-direction. Moreover, *Jk* and *Jc* represent the coupling strengths of each marginal system to the local and non-local RTN, such that we consider two types of system-environment interactions, namely


Furthermore, *Lk*(*t*) and *C*(*t*) denote the random variables used to introduce the stochastic processes. They are used to describe the different conditions under which the subsystems undergo decoherence due to the environment. Here, they represent classical random fluctuating fields such as bistable fluctuators flipping between two fixed values ±*m* at rates *γ<sup>k</sup>* and *γ*, respectively. For simplicity, we assume that *γ<sup>k</sup>* = *γ*. For the *autocorrelation function* of the random variable *η*(*t*) = {*Lk*(*t*); *C*(*t*)} we have *δη*(*t*)*δη*(*t* ) = exp[−2*γ*|*t* − *t* <sup>|</sup>] with a Lorentzian power spectrum *<sup>S</sup>*(*ω*) = <sup>4</sup>*<sup>γ</sup> <sup>ω</sup>*2+*γ*<sup>2</sup> . Defining the parameter *q* = *<sup>γ</sup> <sup>ν</sup>* , we can identify two regimes for the dynamics of quantum correlations: the Markovian regime (*q* 1: fast RTN), and the non-Markovian regime (*q* 1: slow RTN). The time-evolving state of the system under the influence of the RTN is given by

$$
\rho(\{\eta\},t) = \mathcal{U}(\{\eta\},t)\rho(0)\mathcal{U}^\dagger(\{\eta\},t).\tag{33}
$$

in which the time-evolution operator *U*({*η*}, *t*) called the stochastic unitary operator in the interaction picture is given by

$$\mathcal{U}(\{\eta\},t) = \exp\left[-i\int\_0^t \mathcal{H}\_I(t')dt'\right].\tag{34}$$

where *η*(*t*) = {*Lk*(*t*); *C*(*t*)} stands for the different realizations of the stochastic process. Because *U*({*η*}, *t*) depends on the noise, we should perform the ensemble average over the noise fields to obtain the reduced density matrix of the open system, i.e.,

$$
\rho\_{\dot{\rm ic}(\varepsilon t)} = \langle \rho(\{\eta\}, t) \rangle\_{\eta(t)}.\tag{35}
$$

The evolved state of the system in the presence of independent environments (ie) and collective environments (ce) is obtained as

$$\begin{aligned} \rho\_{i\varepsilon}(t) &= \langle \langle \rho(\theta\_A(t), \theta\_B(t), t) \rangle\_{\theta\_A} \rangle\_{\theta\_B} \\ \rho\_{\text{cr}}(t) &= \langle \rho(\theta(t), t) \rangle\_{\theta} \end{aligned} \tag{36}$$

where *θk*(*t*) = *ν* \$ *t* <sup>0</sup> *Lk*(*t* )*dt* (*k* = *A*, *B*) and *θ*(*t*) = *ν* \$ *t* <sup>0</sup> *C*(*t* )*dt*. Calculation of the above terms requires the computation of averaged terms of the type *e*±*in<sup>θ</sup>* (*<sup>n</sup>* <sup>∈</sup> *<sup>N</sup>*) given by [99]

$$\begin{aligned} \langle e^{in\theta} \rangle = D\_{\text{ll}}(\tau) = \langle \cos \left( n\theta \right) \rangle \pm i \langle \sin \left( n\theta \right) \rangle, \\ \langle \sin \left( n\theta \right) \rangle = 0, \end{aligned} \tag{37}$$

$$\langle \cos \left( n \theta \right) \rangle = \begin{cases} e^{-q \tau} \left[ \cosh \left( \xi\_{qn} \tau \right) + \frac{q}{\xi\_{qn}} \sinh \left( \xi\_{qn} \tau \right) \right], & q > n \\ e^{-q \tau} \left[ \cos \left( \xi\_{nq} \tau \right) + \frac{q}{\xi\_{nq}} \sin \left( \xi\_{nq} \tau \right) \right], & q < n \end{cases}$$

where *<sup>ξ</sup>ab* <sup>=</sup> <sup>√</sup>*a*<sup>2</sup> <sup>−</sup> *<sup>b</sup>*<sup>2</sup> ((*a*, *<sup>b</sup>*) = *<sup>n</sup>*, *<sup>q</sup>*), and *<sup>τ</sup>* <sup>=</sup> *<sup>ν</sup><sup>t</sup>* denotes the scaled (dimensionless) time [71].

**Pure initial state in the presence of independent classical environments.** Here, we assume that each of the qubits and qutrits interact locally with local RTN, while the composite system starts with the pure initial state in Equation (25). For this case, the elements of evolved density matrix are given in Appendix A.2. Then the HSS is obtained as

$$H\text{SS} = \frac{1}{6}\sqrt{D\_1^2(\tau) + 2D\_2^2(\tau) + D\_2^2(\tau)D\_1^2(\tau) + D\_2^4(\tau)}.\tag{38}$$

In Figure 4, we illustrate the time behaviors of the negativity, MID and HSS in the non-Markovian regime as a function of the dimensionless time. It is clear that when the entanglement sudden death occurs, the HSS and MID synchronously oscillate with time as they are suppressed to the minimum value and then rise. Moreover, at the first revival of the measures, the minimum point of the HSS exactly coincides with that of the negativity. After that moment we see that maximum (minimum) points of the HSS are in complete coincidence with maximum (minimum) points of the negativity as well as the MID. This perfect qualitative agreement between HSS and entanglement or quantum correlations is evidence that the HSS-based witness can precisely detect non-Markovianity in the presence of classical noises.

**Figure 4.** Evolution of negativity, MID and HSS as a function of dimensionless time *τ* = *νt* when each subsystem of the hybrid qubit–qutrit system, starting from the initial pure state, is independently subject to a random telegraph noise in non-Markovian regime *q* = 0.1.

**Mixed initial state in the presence of independent classical environments.** Now we compare the dynamics of the HSS, obtained from the initial pure state of Equation (25), with the evolution of the negativity and quantum correlation computed for the initial mixed state of Equation (27). The evolved density matrix, the corresponding negativity and quantum correlation are obtained from, respectively, Equations (29)–(31) replacing F with *D*2(*τ*) 2 .

Figure 5 exhibits this comparison for different values of the entanglement parameter *p*. Not considering the periods when the sudden death of the entanglement occurs, we observe that the maximum and minimum points of the measures are very close to each other and small deviations originate from the fact that the initial state, used for computation of the HSS-based measure, should be optimized over all possible parametrizations. Therefore, the HSS-based measure remains as a valid non-Markovianity identifier in the presence of the classical noises.

**Mixed initial state in the presence of a common classical environment.** Let us now compare the dynamics of the HSS, obtained as usual from the initial pure state of Equation (25) by definition, with the evolution of the negativity and quantum correlation computed for the initial mixed state of Equation (27), when both the qubit and the qutrit are embedded into a common RTN source in the non-Markovian regime. The elements of the evolved dynamical density matrix are given in Appendix A.3. Then, one can easily determine the HSS as

$$HSS = \frac{1}{6}\sqrt{D\_1(\tau)^2 + 2D\_2(\tau)^2 + D\_3(\tau)^2 + D\_4(\tau)^2}.\tag{39}$$

**Figure 5.** Comparing the evolution of negativity and MID computed for the initial mixed state of the hybrid qubit–qutrit system, when each subsystem is independently coupled to a random telegraph noise, with HSS (obtained from the initial pure state) for different values of the entanglement parameter *p* in the non-Markovian regime: *q* = 0.1.

Moreover, the evolved density matrix of the hybrid qubit–qutrit system for the initial mixed state of Equation (27) is obtained as

$$
\rho(t) = \begin{pmatrix}
\frac{p}{2} & 0 & 0 & 0 & 0 & \frac{p}{2}\mathcal{F}e^{i\phi} \\
0 & \frac{p}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{p}{2} & 0 \\
\frac{p}{2}\mathcal{F}e^{-i\phi} & 0 & 0 & 0 & 0 & \frac{p}{2}
\end{pmatrix},\tag{40}
$$

where F = *D*4(*τ*).

As a consequence, we find that the negativity and MID are, respectively,

$$\mathcal{N} = \frac{1}{4} [(p - 1) + |3p - 1| + |(1 - 2p) - p\mathcal{F}| + |(1 - 2p) + p\mathcal{F}|].\tag{41}$$

$$\mathcal{M} = (1 - 2p) + \frac{p}{2}(1 + \mathcal{F})\log\left(1 + \mathcal{F}\right) + \frac{p}{2}(1 - \mathcal{F})\log\left(1 - \mathcal{F}\right). \tag{42}$$

For common environments, we know that mutual interaction between subsystems, induced by the common environment, may lead to the preservation of correlations or even result in creation of quantum correlations between the subsystems [82,100–102]. Therefore, revivals of the quantum correlations cannot be necessarily linked to pure non-Markovianity effects and hence we do not expect complete consistency between the HSS and quantum correlations behaviors (see Figure 6 demonstrating this feature of common environments causing the MID to fail in detecting non-Markovianity). Except for these situations, we see that the maximum (minimum) points of the HSS computed for the initial pure state are very close to those of the MID calculated for the initial mixed state.

It should be noted that the classical environments cannot store any quantum correlations on their own, and hence they do not become entangled with their respective quantum systems. Accordingly, common interpretation of non-Markovianity in accordance with inflow (outflow) of information to (from) the system may be problematic in the presence of

the RTN and other similar classical noises [47,103]. In other words, it is somewhat misleading to talk about information flow from the system(s) to the environment(s) or information backflow from the environment(s) to the system(s). The better interpretation is to say that the quantum system has a recording memory of the events affecting its dynamics. When the quantum memory starts remembering, the information about the past events becomes accessible, leading to revival of the quantum correlations and hence to the appearance of quantum non-Markovianity [104].

**Figure 6.** Comparing the evolution of negativity and MID computed for the initial mixed state of the hybrid qubit–qutrit system, when its subsystems are subject to a common RTN source, with HSS (obtained from the initial pure state) for different values of the entanglement parameter *p* in the non-Markovian regime: *q* = 0.1.

#### 3.2.3. Composite Classical-Quantum Environments

Here we investigate a hybrid system formed by a qubit subjected to a random telegraph noise and a qutrit independently subjected to a squeezed vacuum reservoir. The Hamiltonian of such a system can be written as

$$\mathcal{H} = \mathcal{H}\_{\text{qb}}(t) \otimes \mathcal{Z}\_{\text{qt}} + \mathcal{Z}\_{\text{qb}} \otimes \mathcal{H}\_{\text{qt}}(t). \tag{43}$$

where Iqb(qt) denotes the identity operator acting on the subspace of the qubit (qutrit). Moreover, the Hamiltonians of the local interaction of the qubit and qutrit, H*qb*(*t*) and H*qt*(*t*), as well as their corresponding evolution operators, U*qb*(*θ*, *t*) and U*qt*(*θ*, *t*) can be extracted from Sections 3.2.2 and 3.1. In addition, one can consider the unitary evolution operator of the system as U = U*qb*(*θ*, *t*) ⊗ U*qt*(*t*). Then, the evolved density matrix of the this system can then be obtained by averaging the unitary evolved density matrix over the stochastic process induced by the RTN.

**Pure initial state.** The elements of the evolved density matrix when starting from the pure state of Equation (25) are given in Appendix A.4, leading to the following expression for the HSS:

$$HSS = \frac{1}{6}\sqrt{\left(e^{-2\gamma(t)} + e^{-8\gamma(t)}\right)\left(1 + D\_2(\tau)^2\right) + D\_2(\tau)^2}.\tag{44}$$

The time behaviors of negativity, MID and HSS are shown in Figure 7 illustrating that all measures exhibit simultaneous oscillations with time such that their maximum and minimum points exactly coincide. This excellent agreement confirms the faithfulness of the HSS-based measure to detect memory effects.

**Figure 7.** Evolution of negativity, MID and HSS as a function of dimensionless time *τ* when the subsystems of the hybrid qubit–qutrit system, starting from the initial pure state, are independently subject to composite classical-quantum environments. The values of the other parameters are given by *α* = 0.1, *ω<sup>c</sup>* = 20 *ω*0, *r* = 0.3, and *ν* = 100.

**Mixed initial state.** Using Equation (27) as the initial state and computing the evolved state of the system (see Appendix B.4), we find that the the negativity and MID, respectively, are in the form of Equations (30) and (31) with <sup>F</sup> <sup>=</sup> *<sup>D</sup>*2(*τ*)*e*−4*γ*(*t*). In Figure 8, the dynamics of negativity and MID, obtained for the initial mixed state, has been compared with that of the HSS (computed for the initial pure state) in the non-Markovian regime.

**Figure 8.** Comparing the evolution of the negativity and MID, computed for the initial mixed state of the hybrid qubit–qutrit system, when the subsystems are independently subject to composite classical-quantum environments, with the HSS obtained from the initial pure state for different values of the entanglement parameter *p* in the non-Markovian regime: *q* = 0.1. The values of the other parameters are given by *α* = 0.1, *s* = 3, *ω<sup>c</sup>* = 20 *ω*0, *p* = 0 and *v* = 100.

The related analyses are similar to those in the above discussed scenarios, showing that the HSS-based witness may be a proper non-Markovianity identifier even if the initial state of high-dimensional systems is not pure.

#### **4. Conclusions**

Recently, the HSS-based witness, a quantifier of quantum statistical speed which has the advantage of avoiding the diagonalization of the evolved density matrix, has been introduced as a trustful witness of non-Markovianity in low-dimensional systems [61]. In this work, we have generalized this result showing that the proposed witness is a bonafide identifier of non-Markovianity for high-dimensional and multipartite open quantum systems with finite Hilbert spaces. This result stems from the observation that the HSSbased witness is in perfect agreement with established non-Markovianity identifiers based on the dynamical breakdown of monotonicity for quantum information resources, such as negativity and measurement-induced disturbance. We have found that, despite the common interpretation of non-Markovianity in terms of backflow of information from the environment to the system may be problematic [6], the HSS-based witness is capable to detect memory effects of the evolved quantum system.

In order to construct a non-Markovianity measure on the basis of a geometric distance between two quantum states, one of desirable properties is that the distance is contractive, i.e., nonincreasing under any completely positive trace preserving (CPTP) map. It has been shown that the HSS is contractive under CPTP maps in low-dimensional Hermitian systems [61]. Checking all of the dynamical cases presented here, we have found that the contractivity of the HSS holds not only in low dimensional systems but also in finite highdimensional ones. Recently, an HSS-like measure has been used to analyze the quantum speed limit for continuous-variable systems following Gaussian preserving dynamics [105]. Therefore, our results also motivate further studies about HSS applications in detecting non-Markovianity in continuous variable systems.

By definition, the HSS-based witness of memory effects is obtained by maximizing the speed of a classical distance measure between the probability distributions, over all quantum measurements. This, as a prospect, may induce the idea of the the possibility to use classical-like description of density matrix properties in probability representation of quantum mechanics.

Recently, K. Goswami et al. [106] have reported a quantum-optics experimental setup to implement a non-Markovian process—specifically, a process with initial classical correlations between system and environment. It should be noted that in all systems investigated in this paper we have adopted the usual assumption that the system and its environment are initially uncorrelated. It would be interesting to generalize the application of the HSSbased non-Markovianity witness to scenarios in which initial correlations between the system and environment rise. This will be studied in detail in our future work.

**Author Contributions:** Conceptualization, H.R.J.; methodology, K.M., M.K.S., H.R.J. and R.L.F.; formal analysis, K.M. and M.K.S.; investigation, K.M., M.K.S., H.R.J. and R.L.F.; writing—original draft preparation, K.M.; writing—review and editing, H.R.J., R.M. and R.L.F.; supervision, R.M. and R.L.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** R.M. acknowledges support from NSERC, MEI and the CRC program in Canada. R.L.F. acknowledges support from "Sistema di Incentivazione, Sostegno e Premialità della Ricerca Dipartimentale" of the Department of Engineering, University of Palermo.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **Appendix A. Pure Hybrid Qubit–Qutrit Evolved Density Matrix**

This appendix presents the elements of the evolved density matrix of hybrid qubit– qutrit system, starting from the initial pure state of Equation (25), in the presence of quantum and classical noises. This evolved state is required for the assessment of non-Markovianity via the HSS-based witness.

#### *Appendix A.1. Squeezed Vacuum Reservoirs*

The elements of the evolved density matrix, when each subsystem of the hybrid qubit– qutrit system is independently subject to a squeezed vacuum reservoir, in the computational basis |00, |01, |02, |10, |11, |12 are given by

$$\rho\_{11}(t) = \rho\_{22}(t) = \rho\_{33}(t) = \rho\_{44}(t) = \rho\_{55}(t) = \rho\_{66}(t) = \frac{1}{6},$$

$$\rho\_{12}(t) = \rho\_{14}(t) = \rho\_{21}^\*(t) = \rho\_{41}^\*(t) = \frac{1}{6}e^{i\theta}e^{-\gamma(t)},$$

$$\rho\_{13}(t) = \rho\_{31}^\*(t) = \frac{1}{6}e^{i\theta}e^{-4\gamma(t)}, \quad \rho\_{15}(t) = \rho\_{51}^\*(t) = \frac{1}{6}e^{-2\gamma(t)},$$

$$\rho\_{16}(t) = \rho\_{61}^\*(t) = e^{i\theta}e^{-5\gamma(t)},$$

$$\rho\_{23}(t) = \rho\_{25}(t) = \rho\_{36}(t) = \rho\_{45}(t) = \rho\_{52}(t) = \rho\_{53}(t) = \rho\_{54}(t) = \rho\_{56}(t) \tag{A1}$$

$$= \rho\_{63}(t) = \rho\_{65}(t) = \frac{1}{6}e^{-\gamma(t)},$$

$$\rho\_{46}(t) = \rho\_{26}(t) = \rho\_{35}(t) = \rho\_{42}(t) = \rho\_{53}(t) = \rho\_{62}(t) = \frac{1}{6}e^{-2\gamma(t)},$$

$$\rho\_{34}(t) = \rho\_{43}(t) = \frac{1}{6}e^{-5\gamma(t)}.$$

#### *Appendix A.2. Independent Random Telegraph Noise*

The elements of the evolved density matrix, when each subsystem of the hybrid qubit–qutrit system is independently subject to the classical random telegraph noise, can be obtained as

$$\rho\_{11}(t) = \rho\_{22}(t) = \rho\_{33}(t) = \rho\_{44}(t) = \rho\_{55}(t) = \rho\_{66}(t) = \frac{1}{6}$$

$$\rho\_{12}(t) = \rho\_{21}^\*(t) = \frac{1}{6}e^{i\phi}D\_1(\tau)$$

$$\rho\_{23}(t) = \rho\_{32}(t) = \rho\_{45}(t) = \rho\_{54}(t) = \rho\_{56}(t) = \rho\_{65}(t) = \frac{1}{6}D\_1(\tau)$$

$$\rho\_{13}(t) = \rho\_{14}(t) = \rho\_{31}^\*(t) = \rho\_{41}^\*(t) = \frac{1}{6}e^{i\phi}D\_2(\tau)$$

$$\rho\_{25}(t) = \rho\_{36}(t) = \rho\_{46}(t) = \rho\_{52}(t) = \rho\_{63}(t) = \rho\_{64}(t) = \frac{1}{6}D\_2(\tau) \tag{A2}$$

$$\rho\_{15}(t) = \rho\_{51}^\*(t) = \frac{1}{6}e^{i\phi}D\_2(\tau)D\_1(\tau)$$

$$\rho\_{24}(t) = \rho\_{26}(t) = \rho\_{35}(t) = \rho\_{42}(t) = \rho\_{53}(t) = \rho\_{62}(t) = \frac{1}{6}D\_2(\tau)D\_1(\tau)$$

$$\rho\_{16}(t) = \rho\_{61}^\*(t) = \frac{1}{6}e^{i\Phi}D\_2^2(\tau)$$

$$\rho\_{34}(t) = \rho\_{43}(t) = \frac{1}{6}D\_2^2(\tau).$$

#### *Appendix A.3. Common Random Telegraph Noise*

The elements of the evolved density matrix, when the qubit and qutrit are subject to a common RTN source, are given by

$$\rho\_{11}(t) = \rho\_{22}(t) = \rho\_{33}(t) = \rho\_{44}(t) = \rho\_{55}(t) = \rho\_{66}(t) = \frac{1}{6}$$

$$\rho\_{12}(t) = \rho\_{21}^\*(t) = \frac{1}{6}e^{j\phi}D\_1(\tau)$$

$$\rho\_{23}(t) = \rho\_{32}(t) = \rho\_{24}(t) = \rho\_{42}(t) = \rho\_{55}(t) = \rho\_{53}(t) = $$

$$\rho\_{45}(t) = \rho\_{54}(t) = \rho\_{56}(t) = \rho\_{65}(t) = \frac{1}{6}D\_1(\tau)$$

$$\rho\_{13}(t) = \rho\_{14}(t) = \rho\_{31}^\*(t) = \rho\_{41}^\*(t) = \frac{1}{6}e^{j\phi}D\_2(\tau)$$

$$\rho\_{25}(t) = \rho\_{36}(t) = \rho\_{65}(t) = \rho\_{35}(t) = \rho\_{64}(t) = \frac{1}{6}D\_2(\tau) \tag{A3}$$

$$\rho\_{15}(t) = \rho\_{56}^\*(t) = \rho\_{55}^\*(t) = \rho\_{52}(t) = \rho\_{62}(t) = \rho\_{62}(t) = \frac{1}{6}D\_2(\tau)D\_1(\tau)$$

$$\rho\_{16}(t) = \rho\_{61}^\*(t) = \frac{1}{6}e^{j\phi}D\_4(\tau)$$

$$\rho\_{34}(t) = \rho\_{43}(t) = \frac{1}{6}.$$

*Appendix A.4. Composite Classical-Quantum Environments*

The elements of the evolved density matrix, when the qubit and qutrit are independently subject to, respectively, random telegraph noise channel and squeezed vacuum reservoirs, can be obtained as

$$\rho\_{11}(t) = \rho\_{22}(t) = \rho\_{33}(t) = \rho\_{44}(t) = \rho\_{55}(t) = \rho\_{66}(t) = \frac{1}{6}$$

$$\rho\_{12} = \rho\_{21}^\* = \frac{1}{6}e^{i\rho}e^{-\gamma(t)}$$

$$\rho\_{23}(t) = \rho\_{32}(t) = \rho\_{45}(t) = \rho\_{54}(t) = \rho\_{56}(t) = \rho\_{65}(t) = \frac{1}{6}e^{-\gamma(t)}$$

$$\rho\_{13}(t) = \rho\_{31}^\*(t) = \frac{1}{6}e^{i\phi}e^{-4\gamma(t)}$$

$$\rho\_{14}(t) = \rho\_{41}^\*(t) = \frac{1}{6}e^{i\phi}D\_{2}(\tau)$$

$$\rho\_{15}(t) = \rho\_{51}^\*(t) = \frac{1}{6}e^{i\phi}D\_{2}(\tau)e^{-\gamma(t)} \tag{A4}$$

$$\rho\_{16}(t) = \rho\_{56}^\*(t) = \rho\_{52}^\*(t) = \rho\_{63}(t) = \frac{1}{6}D\_{2}(\tau)$$

$$\rho\_{24}(t) = \rho\_{26}(t) = \rho\_{35}(t) = \rho\_{24}(t) = \rho\_{62}(t) = \frac{1}{6}D\_{2}(\tau)e^{-\gamma(t)}$$

$$\rho\_{34}(t) = \rho\_{43}(t) = \frac{1}{6}D\_{2}(\tau)e^{-4\gamma(t)}$$

$$\rho\_{45}(t) = \rho\_{64}(t) = \frac{1}{6}e^{-4\gamma(t)}.$$

#### **Appendix B. Mixed Hybrid Qubit–Qutrit Evolved Density Matrix**

This appendix presents the elements of the evolved density matrix of hybrid qubit– qutrit system, starting from the initial mixed state of Equation (27), in the presence of quantum and classical noises.

#### *Appendix B.1. Squeezed Vacuum Reservoirs*

The elements of the evolved density matrix, when each subsystem of the hybrid qubit–qutrit system is independently subject to a squeezed vacuum reservoir, are given by

$$
\rho(t) = \begin{pmatrix}
\frac{p}{2} & 0 & 0 & 0 & 0 & \frac{p}{2}\mathcal{F} \\
0 & \frac{p}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2}\mathcal{F} & 0 & 0 \\
0 & 0 & \frac{1-2p}{2}\mathcal{F} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{p}{2} & 0 \\
\frac{p}{2}\mathcal{F} & 0 & 0 & 0 & 0 & \frac{p}{2}
\end{pmatrix},\tag{A5}
$$

and the partial transpose with respect to the subsystem *A* is

$$\left(\rho(t)^{AB}\right)^{T\_A} = \begin{pmatrix} \frac{p}{2} & 0 & 0 & 0 & 0 & \frac{1-2p}{2}\mathcal{F} \\ 0 & \frac{p}{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1-2p}{2} & \frac{p}{2}\mathcal{F} & 0 & 0 \\ 0 & 0 & \frac{p}{2}\mathcal{F} & \frac{1-2p}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{p}{2} & 0 \\ \frac{1-2p}{2}\mathcal{F} & 0 & 0 & 0 & 0 & \frac{p}{2} \end{pmatrix},\tag{A6}$$

where the <sup>F</sup> <sup>=</sup> *<sup>e</sup>*−5*γ*((*t*)).

#### *Appendix B.2. Independent Random Telegraph Noise*

The elements of the evolved density matrix, when each subsystem of the hybrid qubit– qutrit system is independently subject to the classical random telegraph noise, are given by Equation (A5) with F = *D*2(*τ*) 2 .

#### *Appendix B.3. Common Random Telegraph Noise*

The evolved density matrix, when the qubit and qutrit are subject to a common RTN source, is given by

$$
\rho(t) = \begin{pmatrix}
\frac{p}{2} & 0 & 0 & 0 & 0 & \frac{p}{2}\mathcal{F} \\
0 & \frac{p}{2} & 0 & 0 & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & \frac{1-2p}{2} & \frac{1-2p}{2} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{p}{2} & 0 \\
\frac{p}{2}\mathcal{F} & 0 & 0 & 0 & 0 & \frac{p}{2}
\end{pmatrix},\tag{A7}
$$

where F = *D*4(*τ*).

#### *Appendix B.4. Composite Classical-Quantum Environments*

The elements of the evolved density matrix, when the qubit and qutrit are independently subject to, respectively, random telegraph noise channel and squeezed vacuum reservoirs, are given by Equation (A5) with <sup>F</sup> <sup>=</sup> *<sup>D</sup>*2(*τ*)*e*−4*γ*(*t*).

#### **References**

