Universal non-Markovianity detection in hybrid open quantum systems

A universal characterization of non-Markovianity for any open hybrid quantum systems is presented. This formulation is based on the negativity volume of the generalized Wigner function, which serves as an indicator of the quantum correlations in any composite quantum systems. It is shown, that the proposed measure can be utilized for any single or multi-partite quantum system, containing any discrete or continuous variables. To demonstrate its power in revealing non-Markovianity in such quantum systems, we additionally consider a few illustrative examples.

In this contribution, we employ the definition of GWF for arbitrary quantum system to introduce a new measure of non-Markovianity. We demonstrate the power of the proposed measure to reveal non-Markovian dynamics in composite systems on two canonical examples, namely, on qubit-qubit and qubit-squeezed coherent states. The introduced measure is based on the changes of the negativity volume (NV) of GWF when the system is interacting with its environment 34,35 . As we show below, the Markovian dynamics is reflected by a pure monotonic behaviour (decreasing) of the NV, as the system and the environment became more and more entangled. On the other side, the non-Markovian dynamics shows reduced (and also increasing) changes of NV. This feature offers us to quantify a degree of non-Markovianity based on the negativity volume of the GWF. Most importantly, such behaviour is observed for any CV, DV, or hybrid quantum systems, all being either single or multipartite. We also note that the concept of negativity volume of the standard Wigner function in a CV domain has been already utilized in detection of non-Markovianity in bosonic systems 36 .

Hybrid state evolution
A generic open quantum system dynamics can be described by a master equation: where L is, in general, a time-dependent Liouvillian superoperator, and for a weakly coupled bosonic Markovian environment the Liouvillian L attains the Gorini-Kossakowski-Sudarshan-Lindblad form 37,38 . For a hybrid system, consisting of the discrete and continuous parts, the Liouvillian L can be decomposed as follows where The Liouvillian superoperators L 0 , L d , and L c describe the coherent, incoherent discrete and incoherent continuous system evolution, respectively. Such superoperators can induce either dissipation or amplification in reduced quantum systems 39 . L is a Lindblad operator with L † being its Hermitian conjugate; the square [ ] and curl { } brackets denote commutator and anticommutator, respectively; the symbol · is a placeholder to indicate where an operator should be applied; and the damping coefficients γ (t) are, in general, considered to be timedependent. However, a weakly coupled bosonic non-Markovian environments can be as well introduced in Eq.
(3) as the damping coefficients γ (t) are allowed to flip their signs 25 .
The density matrix ρ of the hybrid state in Eq. (1) can be, thus, written as following 1,40 following 1,40 , and ρ ′ i (ρ i ) are left(right) Liouvillian eigenmatrices. We note that because the Liouvillian is a non-Hermitian superoperator, it can attain both left and right eigenmatrices. Eventually, for the time-dependend coefficients in Eq. (3) is beneficial to express Eq. (4) in terms of more suitable states, whose evolution can be easily calculated. The total Hilbert space of a system is Hilbert space of discrete (continuous) part of a system. Moreover, due to the CV part, the total Hilbert space is infinite, i.e., N → ∞ , since N c → ∞ . Nevertheless, it is always possible to make a truncation of Hilbert space at some finite N without affecting much the very description of a system, though N can still be large. Indeed, in real experiments one always deals with a finite number of photons in a system, thus, the truncation assumption, in practice, can be justified [41][42][43] .
The generalized Wigner function W for any hybrid system, is defined as follows 30 : where ˆ is a kernel operator, which is represented by a tensor product of kernels corresponding to DV and CV subsystems. In particular, for a bipartite hybrid system composed of a qubit and photonic part, the kernel operator in Eq. (5) can be written as ˆ =ˆ q ⊗ˆ p , where ˆ q and ˆ p are kernel operators corresponding to the qubit and photonic subsystems, respectively, and which read as 30 : The operator Û is a rotational operator in SU(2) algebra, namely Û = e iσ 3 φ e iσ 2 θ e iσ 3 � with Pauli operators σ i , i = 1, 2, 3 , and angles φ, � ∈ [0, 2π] , θ ∈ [0, π] . � q =Î 2 − √ 3σ 3 is a parity operator of the qubit. The operator D is the displacement operator of the coherent state, i.e., D = eâ www.nature.com/scientificreports/ boson operator. The corresponding bosonic parity operator reads as � p = e iπâ †â . The form of a kernel operator =ˆ q ⊗ˆ p ⊗ · · · for a more complex hybrid states, e.g., consisting of N-level systems, can be found in Refs. 30,31 . The normalization condition Wd = 1 is obtained by means of an appropriate integral measure d , which is a product of normalized differential volume of SU(2) space of a qubit with a Haar measure dν 31,44,45 . Now, by combining Eqs. (4) and (5), one arrives at the following expression for the Wigner function of the initial state ρ: where W i are Wigner functions corresponding to ρ i obtained using Eq. (5).
The GWF negativity volume is found as 34,35 The usefulness of using NV of the GWF lies in its ability to detect quantum correlations in hybrid systems 34 , and thus its potential in identification of non-Markovian backflow of nonclassicality from environment to such systems. It is worth mentioning that the GWF can already exhibit negativity for single qubit states, which do not possess any quantum correlations 30 , that is in a striking contrast to ordinary Wigner function, defined for continuous states, where the negative values signal the presence of nonclassicality 39 .
The divisibility condition for quantum maps 12 mentioned earlier implies that the NV of the GWF remains a monotonic function of time and such the vector v , in Eq. (7), can only contract. Indeed, for the Markovian dynamics the components of the vector v exponentially decay. This means that any deviation of system dynamics from, in general, the time-inhomogeneous Markovianity, will be reflected in the non-monotonical behaviour of the NV of the GWF. Therefore, any non-Markovian processes can induce the increase of the vector of states v , namely w , and, thus, a sudden increase of the Wigner function negativity volume, which allows to detect a backflow of nonclassicality from an environment to a system. This is related to the flip of the sign in the decaying rates γ < 0 , will be detected by the change in the sign of the speed of the NV N , which for Markovian dynamics is always negative, i.e, dN M /dt < 0 . Thus, in order to quantify the degree of non-Markovianity, we can define the following measure: which can vary in the range D N ∈ [0, 1] , and, therefore, the Markovianity is determined whenever D N = 0.
We will demonstrate the ability of the measure D N , given in Eq.

Environment models
In what follows, in our demonstration of the ability of the NV of the GWF to detect the non-Markovian dynamics of hybrid systems, without loss of generality, we will focus on hybrid systems whose discrete part is formed by a qubit and the continuous part by a photonic field.
As such, we consider the following incoherent Lindblad operators, in Eq. (3), for the qubit and photonic parts, respectively: Namely, the qubit decoherence is encoded by the Pauli matrix σ z . Whereas the decoherence of the photonic part is represented by the amplitude damping L p,1 ∝â and photon dephasing L p,2 ∝n processes, where â and n are the annihilation and photon number operators, respectively. The damping functions of the qubit and photonic parts, in Eq. (10), are denoted as γ q (t) and κ i (t) , i = 1, 2 , respectively. Also, we assume that the Hamiltonian Ĥ in Eq. (3) describes the free evolution of the qubit and photonic field. It attains the form Ĥ = ω q0 2σ z + ω p0n , where ω q0 and ω p0 is the qubit transition frequency and the photonic field frequency, respectively. We will work in the interaction picture, and, for simplicity, we set ω q0 = ω p0 = ω 0 , as we consider both decoherences acting separately in examined cases.
The non-Markovian evolution is introduced via the time convolution-less (TCL) projection operator techniques employing up to the fourth-order correction terms representing the damping functions γ q (t) , as well as κ 1 (t) and κ 2 (t) , using the following expression originating in a perturbation expansion 1 : and c 2 , c 4 are non-negative strength coefficients 1 , which are employed to simply control the interaction, as they can be absorbed to coefficients for a reservoir. Again, here, the form of the damping function γ (t) formally embodies all the damping functions given in Eq. (10). That is, the functions γ q (t) , κ 1 (t) and κ 2 (t) are all decomposed in the fashion presented in Eq. (11). As we are only interested in the evolution of NV, the frequency Lamb shifts induced by coupling to environments, which are time dependent, are ignored. One can also consider exactly solvable models of non-Markovian environments [46][47][48] . For the sake of simplicity, we consider all damping functions in the form of Eq. (11) and γ q (t) = κ 1 (t) = κ 2 (t) as we always select only one type of open quantum dynamics, given by Eq. (3), being enabled during a time evolution (for more details see 39 ). Nevertheless, all conclusions drawn in the following parts are applicable as well to any model of system-environment interactions, even with all open quantum dynamics combined together.
Typical representative environments are described by the Lorentzian and Ohmic-type spectral functions. The first type of environment is related to an atom placed inside of a detuned optical cavity, which is represented by the Jaynes-Cummings model 39 . Similar spectral function also describes a qubit, represented by a spin-1/2 system, coupled to a spin environment 49 . Such reservoir is characterized by the Lorentzian function 50 : The parameter is related to the relaxation rate and δ 0 to the coupling. The detuning of the cavity mode from the system transition frequency ω 0 is denoted by . The corresponding correlation function for Eq. (14) is obtained as: Illustrative profiles of damping functions are shown in Fig.1a. The second type is represented by the Ohmictype spectral function, which can imitate many different thermal baths and is applicable, for instance, to superconducting Josephson junctions 51 and nano-mechanical resonators 52 . The corresponding spectral function is given by 5 : where ω c is the critical wavelength of the environment spectra and η is the coupling strength. The parameter s denotes the Ohmic for s = 1 , sub-Ohmic for 0 < s < 1 , and super-Ohmic for s > 1 cases. Three illustrative realizations are shown in Fig.1b. The corresponding correlation function is obtained as: The symbol Ŵ marks the gamma function. www.nature.com/scientificreports/

Numerical simulations
In this section, we numerically compute the NV of the GWF of the qubit-qubit and qubit-squeezed coherent states subjected to different non-Markovian environments. A photonic part in the second hybrid system is represented by a squeezed coherent state. For that we numerically integrate the master equation in (1) by applying the fourth-order Runge-Kutta method for different time-dependent damping functions, considered in the previous section. The choice of the studied states is dictated by the frequent use of the qubit and squeezed coherent states in entanglement-based quantum information protocols utilizing the DV and CV states, respectively 53,54 . By knowing the dynamics of a density matrix ρ(t) , we then calculate the non-Markovianity measure D N , in Eq. (9), for studied states to demonstrate its ability to quantitatively capture the deviation of the time evolution of hybrid states from the ordinary Markovian dynamics, which is also reflected in the qualitative non-monotonic behaviour of the NV of the GWF.
Qubit-qubit state. Let us start, first, from the qubit-qubit state, written in the following form: The state in Eq. (18) is one of the maximally entangled Bell states 55 . Now, lets assume that the first qubit is coupled to a non-Markovian dephasing channel, modelled by the Lindblad operator L q in Eq. (10). In Fig. 2a we plot the time evolution of the NV N of the GWF along with corresponding values of D N for the qubit-qubit state for three different cases of Lorentzian environments, which have been depicted in Fig. 1a. The blue curve corresponds to the Markovian case, as the NV only gradually decreases in time. The other two curves (orange and green) exhibit a partial revival in the negativity volumes indicating the non-Markovian evolution. The similar trends in the behaviour of the NV for the qubit-qubit state can be observed for the Ohmic-type environments shown in Fig. 2b. The different kinds of the Ohmic bath have been displayed in Fig. 1b.
Note also that in both Fig. 2a,b, the NV tends to the value N = 1/ √ 3 − 1/2 ≈ 0.08 with t → ∞ , as the initial qubit-qubit entanglement eventually completely degrades. This is so-called critical value of the NV, namely, it is the maximal value of the NV which is generated by non-entangled, i.e., 'classical' qubits. Above that critical value, the NV indicates the presence of the quantum correlations in the system 34 .
Qubit-squeezed coherent state. The second example, we would like to consider, is an entangled hybrid state 56 , where a DV part is presented by a qubit and the CV part by a squeezed coherent state (SCS). This entangled hybrid state reads as The SCS is defined as |ξ i , α i � =D(α i )Ŝ(ξ i )|0� c , where D (α) = exp αâ † − α * â is the displacement operator, the single-mode squeezing operator is Ŝ (ξ ) = exp ξâ †2 2 − ξ * â 2 2 , and |0� c is vacuum state. Moreover, we assume here that only the photonic part, i.e., SCS, is exposed to a non-Markovian environment. For our numerical calculation purposes, the Hilbert space for the photonic part can be effectively truncated as we consider small values of |α| , which result in vanishing amplitudes for large photon numbers 43 .
To visualize the effect of the photon decoherence, we plot the GWF of the reduced photonic part for three different scenarios, namely, before and after the phase and amplitude damping in Fig. 3a-c, respectively. The GWF W p (β) of the reduced photonic state is found by integrating out the degrees of freedom θ of the qubit state, that is W p (β) = dθW(θ, β) (see Ref. 34 for details). In case of the phase damping [presented by the Lindblad operator L p,2 in Eq. (10)], the reduced GWF of the initial SCS transforms as depicted in Fig. 3b. As one can see, the phase of the initial reduced SCS is smeared in time. Although the real amplitude remains unchanged, which is interpreted  Fig. 3c, that is, it just degrades to the vacuum state |0� c . The Lorentzian environment is presented in the Fig. 4 for the both phase and amplitude damping channels and for varying parameters of environments. The blue curve in both plots marks the Markovian case, and the yellow and green curves are related to the non-Markovian case. One can see the related degrees of non-Markovianity, given by Eq. (9), in corners of each plot. Additionally, the main distinctiveness between curves in Fig. 4 is the minimum N achieved in each case. In the case of amplitude damping Fig. 4b, the Wigner function reaches the minimal N in the Markovian dynamics, in a given time interval, whereas other curves are approaching this value only slowly. This means that the squeezed coherent state turns into the vacuum state |0� c represented in Fig. 3c and a source of negativity is hidden in the qubit part. For the phase dumping, the resulting state, which is a mixed state, contains still more correlations than the previous case, which is reflected by higher values of N , see Fig. 3b. A quite similar behavior can be observed as well in Fig. 5 for the Ohmic-type environments.

Conclusion
We have introduced a new measure of non-Markovianity based on the negativity volume of the generalized Wigner function, which can be applied to any quantum composite system, i.e., hybrid system. We have illustrated its power to reveal non-Markovian dynamics by considering a few canonical examples. Namely, we have applied it to the qubit-qubit and qubit-squeezed-coherent hybrid states. Such composite systems have been exposed to different bosonic environments. Despite the fact that the photonic part was only considered to be a single-mode field, all conclusions remain the same for the multi-mode photonic fields as only computation demands increase.
Most importantly, our proposed measure is state-independent and, therefore, is universal, as it is solely based on the GWF defined for any hybrid state. Moreover, the experimental advances 57 allow us to expect that the experimental measurements of the GWF will become a standard tool in characterizing complex composite quantum systems and their dynamics in the near future. Additionally, the knowledge of the GWF offers the opportunity to perform not only quantitative but also qualitative analysis of the time evolution of any quantum hybrid systems 32,33 .
We note, however, that the knowledge of the NV of the GWF, in general, requires an implementation of a state tomography, which can be source demanding. Nonetheless, the task can be substantially simplified if the initial state is known 58 . In this case, one can focus on detecting a specific region of the GWF, which demonstrates the largest concentration of the NV in a system, and, especially, when such a region exhibits an interference-like   . Negativity volume N for the qubit-photonic system, where the photonic part is coupled to the Ohmic-type environments while undergoes of (a) phase and (b) amplitude damping. The blue curves correspond to s = 0.5 , the orange curves to s = 1 and the green curves to s = 3 . In (a) for s = 0.5 η = 0.005 . The remaining parameters are the same as for the Fig. 4.