Numerical analysis of bipartite entanglement evolution in simple cubic 1/2-spin system with additional spin 1 dopant

The system in consideration is a singular simple cubic cell of spins s=12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s = \frac{1}{2}$$\end{document} doped in centre with additional spin S=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S = 1$$\end{document}. The structure is located in external magnetic field and is undergoing dissipative, Markovian evolution. The analysis of system time evolution is focused on the entanglement evaluation and determination of its behaviour over time. We found out the distinct effect of doping the structure with additional spin S=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S = 1$$\end{document} in the time evolution of the bipartite entanglement in |W⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vert W \rangle $$\end{document} state. We also distinguished the raising and lowering spin interaction with environment and determined its impact on decoherence. Furthermore, we have shown that the entanglement between spins is in the form of a damped superposition of Gaussian functions. Its pulsed nature results from the unitary evolution in spin structure, while the damping is caused by the influence of the environment.


Introduction
The entanglement, specific correlation between the elements of composite quantum system, turned out to be one of the most surprising features of a quantum theory, alongside with its non-determinism and discretisation of physical parameters. Since the famous papers by Schrödinger [1] and Einstein, Podolsky and Rosen [2] it remains one of the most frequently discussed topics of modern physics. Development of quantum information theory drew attention to possibility of utilising the entanglement as a resource: in quantum cryptography [3][4][5], quantum computation [6,7] and quantum information processing [8][9][10]. Due to this reason, the gathering knowledge on the entanglement formation and behaviour over time is highly anticipated, and determining the set of parameters for which the entanglement in particular quantum system will be high and stable over time may be used as a guideline for manufacturing the real layouts to implement the quantum informatics protocols on.
There exist quantum states that are highly entangled and states with low entanglement or even without any entanglement at all. There are multiple ways to quantitatively express the amount of entanglement in system: one of them is to determine the conditions on the elements of density matrix to represent separable states in accordance to the Peres-Horodecki criterion [31,32] and then calculating the distance between the analysed entangled state and the set of all separable states, utilising the Schmidt-Hilbert norm [33,34]-the further the state is from the set, the more entangled it is. Other measures include among others: negativity [32], relative entropy of entanglement [33], entanglement of formation [35] or, probably the most common, concurrence [36]. This last measure proven to be a reliable mean to evaluate the amount of entanglement for coupled systems of dimensionality 2 ⊗ 2 and 2 ⊗ 3 [37] and therefore was utilised in presented analysis.
The aim of this work is to determine the time behaviour of the bipartite entanglement in a three-dimensional finite spin structure coupled with Markovian environment. To this end, the Lindblad master equation was applied, and it was solved by numerical finite step method, using the Runge-Kutta algorithm. Entanglement concurrence time evolution was computed for 1 2 -1 2 spins and for 1 2 and 1 spins using initial state, which is equivalent of |W state [38,39] for the system without the S = 1 dopant.
The paper is organised as follows: In Sect. 2, we present the analysed structure and describe the mathematical formalism behind the system Hamiltonian and time evolution of the structure state; we also present the tools to quantify the entanglement in the system. Section 3 covers the presentation and discussion of obtained results. Finally, Sect. 4 serves as a summary of the research.

Model
The analysed model is a finite, three-dimensional structure similar to simple cubic cell with spins s = 1 2 , placed in the corners. The structure consists also an additional dopant with spin S = 1 in the centre. The entire system is located in an external  Fig.1.
If we assume that spins coupling range does not exceed the lattice constant, and they interact according to Heisenberg XXX model, the Hamiltonian of the analysed structure can be written as: where j is exchange integral for spins s = 1 2 and J is the integral for exchange interaction between spin S = 1 and s = 1 2 . The superscript z indicates the component of spin in the direction of the external magnetic field. We are primarily interested in time evolution of system coupled to dissipative environment satisfying the Born-Markov condition, which can be realized for example in a form of a large bosonic reservoir. In aim to determine the change over time of density matrix ρ(t) describing the overall state of the system, we utilize the Lindblad equation: In (2), the first term (commutator) describes coherent time evolution governed by the Hamiltonian, while the second D(ρ) indicates the dissipator, assumed here as: 168 Page 4 of 17 M. Kaczor, P. Jakubczyk with L being Lindblad operator: Parameters γ P and γ N are decoherence rates connected with spin raising and lowering effect of the system interaction with the reservoir, and the operatorsŝ + q,i (ŝ − q,i ) andŜ + (Ŝ − ) denote the raising (lowering) operators for spins 1/2 and 1, respectively. Individual forms of these operators result from formulas: and in calculational basis they can be represented as: To quantitatively describe the strength of bipartite entanglement in the system, the concurrence measure is applied [35,36]. For the entanglement between the two adjacent spins s = 1 2 , this measure must be calculated in respect to subsystem consisting only the two examined spins (be it subsystem A, in no way related to spins s a,i ), so the overall density matrix has to be reduced to a smaller form, in accordance with the formula:ρ where Tr B ρ denotes partial trace of the density matix ρ over the subsystem B consisting of six remaining spins s = 1 2 and spin S = 1; states {|φ B j } j span complete basis in subsystem B. Concurrence of the entanglement is calculated for the reduced density matrix as: In equation above λ 4 ≥ λ 3 ≥ λ 2 ≥ λ 1 are the eigenvalues of a matrix resulted from following transformation on the reduced density matrix: where σ Y 2 is Y-Pauli matrix-lower index to describe the dimensionality of it and does not correspond to the position of spin in the spin structure.
To calculate the concurrence value between the spins S = 1 and s = 1 2 , the density matrix also needs to be reduced, in this case, however, the reduction takes place in (a) Time dependence of entanglement concurrence without interaction with dopant spin, external magnetic field and with no influence from environment (γ P = γ N = 0). Exchange integral is set to be j = 1.
(b) Time dependence of entanglement concurrence without external magnetic field and with no influence from environment (γ P = γ N = 0). Exchange integrals are set as j = 1 and J = 1. respect to subsystem of dimensionality 2 ⊗ 3 (therefore B in equation below is the subsystem consisting the seven spins s = 1 2 that are not taken into consideration): Concurrence for the systems of mentioned dimensionality is calculated similarly to (6): where λ 6 ≥ λ 5 ≥ λ 4 ≥ λ 3 ≥ λ 2 ≥ λ 1 are the eigenvalues of the matrix [40,41]:    [42,43]:

Time evolution of entanglement concurrence
The time behaviour of the entanglement in the system was calculated by numerical methods. Analytical solution for the system without the environment taken into account and cursory solution of this problem with environment included is provided in appendix; however, due to the scale of the problem the results presented below were obtained by computer-assisted calculations. Equation (2) is solved approximately by multi-step method using fourth-order Runge-Kutta algorithm, and for each step the concurrences were computed in accordance to the formulas above. As an initial state, we choose partially entangled state |α , which can be presented in the form: where the first term of the tensor product corresponds to one excitation in a ferromagnetically ordered frame of spins s = 1 2 , while the second represents the state of the S = 1 spin. In other words, the first term is equivalent of |W state for the frame, in notation below the arrows refers to values of spins s = 1 2 z-projection-↑: m S = 1 2 , ↓: m S = − 1 2 , the second term can in general take three values; −1, 0 and 1 and refer to value of spin S = 1 z-projection.
Since our system is highly symmetric, to determine the structure of bipartite entanglement it is sufficient to calculate the entanglement between two spins 1/2 and between any spin 1/2 and spin 1 of the considered spin system. In Fig. 2a, we present evolution of the bipartite entanglement for spins a1 and b1 for the case of no interaction with spin S = 1 (J = 0). If we then turn on the interaction of spins 1/2 with spin 1 (J = 0), then the evolution of bipartite entanglement proceeds as shown in Fig. 2b.
Without the S = 1 dopant, the |α state is an eigenstate of the frame of spins s = 1 2 , therefore in isolated system the time evolution does not change the initial bipartite entanglement; in other words, during unitary evolution all 1/2 spins are equally entangled with concurrence equal to 0.25 (see Fig. 2a). Introduction of the dopant into the structure changes the behaviour of the bipartite entanglement which becomes a time dependent and periodic, both for C 1 2 , 1 2 and for C 1 2 ,1 (see Fig. 2b). Within one period, peaks of highest concurrence can be observed: C 1 2 , 1 2 reaches 0.25 at the begging and the end of the period, while C 1 2 ,1 exceeds the value 0.3 twice during that period. The In contrast to the interaction integral between the spins located at the vertices of the cubic structure under consideration, the dopant interaction integral is relevant. By varying the value of the J integral, we can adjust the period of the entanglement variation by changing its frequency, as shown in Fig. 6. The change in frequency of entanglement evolution holds for both kinds of entanglement: for C 1 2 , 1 2 and for C 1 2 ,1the periods for the former are longer. It is, however, noteworthy that without dopant time of total concurrence loss for C 1 2 , 1 2 is extended in comparison with situations with dopant, but the price is impossibility of utilising the entanglement between spins s = 1 2 and S = 1. The next step of our simulations was to analyse the different characteristics of dissipative environment. We prepared three sets of various dissipative parameters and conducted a comparative analysis on them. The sets correspond to situations, where there is an equal probability that environment raise or lower any spin in structure (γ N = γ P ); the environment prefers to lower spins than to raise them (γ N > γ P ); the environment prefers to raise spins than lower them (γ N < γ P ). The results are presented in Figs. 7, 8 and 9.
Looking at Fig. 7, it is easy to see that an increase of γ P and γ N parameters has a strongly dampening influence on the entanglement of both C 1 2 , 1 2 and C 1 2 ,1 . However, if we compare the red and blue plots in Figs. 8 and 9, it can be seen that the decisive influence on the decoherence rate is exerted by the spin-lowering factors from the environment. Difference between the plots in Fig. 8 is minimal (due to identical values of γ N ) and is difficult to notice even in close-up (see Fig. 10a), while similar difference for Fig. 9 can be spotted on the original figure-for sake of clarity; we nevertheless pictured the close-up in Fig. 10b-here, the γ P have identical values, so the results confirm the previously stated observation.

Conclusion
We have conducted a numerical analysis of bipartite entanglement evolution in simple cubic cell of spins s = 1 2 with addition of spin S = 1 in the centre, placed in external magnetic field and in dissipative, Markovian environment. We successfully applied the open system formalism, namely Lindblad master equation, in research focused on determining the impact of additional spin on the behaviour of bipartite entanglement between spins s = 1 2 and between spin s = 1 2 and spin S = 1. As an initial state, we assumed the equivalent of W state for structure devoid of an additional node. We compared the behaviour of concurrence in isolated system with and without additional spin, and for different exchange integrals and applied magnetic fields. Moreover, we analysed various environment parameters to determine its impact on the decoherence. We found out that inclusion of the additional spin introduces the state into a behaviour where entanglement pulses can be distinguished (the periods of increased entanglement), which is a direct result of getting the structure out of its eigenstate by expanding the Hilbert space of the system (it gains additional degrees of freedom). We would especially like to draw the readers attention to two periods of entanglement between s = 1 2 and S = 1, where, in isolated system, the concurrence C 1 2 ,1 , a resource available only in doped system, reaches its maxima, which, in turn, are greater than the maximal value of the concurrence between two spins s = 1 2 . We found out that for initial state |α = |W ⊗ |0 neither the external magnetic field, nor the exchange integral between spins s = 1 2 has any impact on the concurrence behaviour in time. However, the exchange integral between spins of the frame and spin S = 1 has-the increase of exchange interaction strength results in higher frequency of the entanglement pulses. In open systems, the analysis showed that both concurrences decrease exponentially; decrease rate is dictated by the environment parameters. However, for our definition of W state it turned out that the rate of lowering spins γ N has a bigger impact on decoherence time than its raising spins counterpart γ P , yet this conclusion appears to result solely from the definition on initial state.
Author Contributions MK and PJ did conceptualisation, formal analysis, visualisation, supervision and writing (review and editing); MK done data curation and provided methodology, software and writing (original draft); PJ validated the study. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement
The data presented in this study are available on request from the corresponding author.

Conflict of interest
The authors declare no conflict of interests or competing interests that can be relevant in the context of this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. m-th row and n-th column of density matrix in this calculational basis, denoted by ρ m,n , is expressed by: where |m and |n are the base vectors of calculational basis. Further calculations (reduction of he Hilbert space, determination of concurrence using (6) or (9)) do not require any complex calculation methods; nevertheless, they are extremely labourintensive and their determination is best entrusted to the computational software.
In order to attempt an analytical solution for an open system, we start the analysis with determining the effect of the dissipator on the p k,l element of density matrix, initially for only one dissipative operatorÂ in Lindblad equation: All the operators we used in dissipator act on the space spanned by the computational basis vectors; therefore, they can be decomposed to the form: which allows to write the k,l-th element of dissipator as: Proceeding as in the case of a closed system, one can derive the result: (A7) As one can notice, a complete solution to determine the explicit form of the density matrix requires solving a system of interconnected differential equations. However, the case discussed in paper allows for two more simplification: firstly, the operators in dissipator are of two kinds-raising and lowering spin operators. In their calculation basis matrix form, they have only two types of elements-0 and 1. The only terms that will not be zeroed out in the above equation are those that relate to the two elements of the matrixÂ that both are equal to one. What is more, the initial state consists of several calculational basis vectors, all of them with the same probability amplitude (for |α this amplitude is ρ 0 = 1 8 ). Assuming operationally that Ĥ , ρ(t) = 0, the Lindblad equation for the situation considered in paper simplifies to: where η can be equal to either +2 or −1, depending on what part of Lindblad operator this term originates from, x can be P (for spin raising operator) or N (for spin lowering operator), and ρ α,β are the terms from (A7) that were associated with two nonzero elements of operator matrix. The selection of these elements can be laborious, yet it is perfectly doable knowing the matrix forms of the operators in question. The implicit form of particular density matrix element reads: ρ 0 · e ηγ X t .
While solution (A9) is far from rigorous, it allows, along with (A5), to draw some conclusions: the unitary density matrix evolution for discussed system without the bosonic environment consists of periodic changes in the population of states that are not eigenstates, while the inclusion of the dissipative environment into account will result in exponential change in time of the density matrix elements. Both of these findings have their confirmation in numerical results presented in the main paper.