document: entanglement in quantum multipartite systems is a unique property in quantum world . 
   it plays an important role in quantum information processing @xcite . 
   therefore , the study of its essential features and dynamical behavior under the ubiquitous decoherence of relevant quantum system has attracted much attention in recent years @xcite . 
   for example , it was found that the entanglement of qubits under the markovian decoherence can be terminated in a finite time despite the coherence of single qubit losing in an asymptotical manner @xcite . the phenomenon called as entanglement sudden death ( esd 
   ) @xcite has been observed experimentally @xcite . 
   this is detrimental to the practical realization of quantum information processing using entanglement . 
   surprisingly , some further studies indicated that esd is not always the eventual fate of the qubit entanglement . 
   it was found that the entanglement can revive again after some time of esd @xcite , which has been observed in optical system @xcite . 
   it has been proven that this revived entanglement plays a constructive role in quantum information protocols @xcite . even in some occasions 
   , esd does not happen at all , instead finite residual entanglement can be preserved in the long time limit @xcite . 
   this can be due to the structured environment and physically it results from the formation of a bound state between the qubit and its amplitude damping reservoir @xcite . 
   these results show rich dynamical behaviors of the entanglement and its characters actually have not been clearly identified . 
   recently , lpez _ et al . 
   _ asked a question about where the lost entanglement of the qubits goes @xcite . 
   interestingly , they found that the lost entanglement of the qubits is exclusively transferred to the reservoirs under the markovian amplitude - damping decoherence dynamics and esd of the qubits is always accompanied with the entanglement sudden birth ( esb ) of the reservoirs . 
   a similar situation happens for the spin entanglement when the spin degree of freedom for one of the two particles interacts with its momentum degree of freedom @xcite . 
   all these results mean that the entanglement does not go away , it is still there but just changes the location . 
   this is reminiscent of the work of yonac _ et al . _ 
   @xcite , in which the entanglement dynamics has been studied in a double jaynes - cummings ( j - c ) model . 
   they found that the entanglement is transferred periodically among all the bipartite partitions of the whole system but an identity ( see below ) has been satisfied at any time . 
   this may be not surprising since the double j - c model has no decoherence and any initial information can be preserved in the time evolution . 
   however , it would be surprising if the identity is still valid in the presence of the decoherence , in which a non - equilibrium relaxation process is involved . in this paper 
   , we show that it is indeed true for such a system consisted of two qubits locally interacting with two amplitude - damping reservoirs . 
   it is noted that although the infinite degrees of freedom of the reserviors introduce the irreversibility to the subsystems , this result is still reasonable based on the fact that the global system evolves in a unitary way . 
   furthermore , we find that the distribution of the entanglement among the bipartite subsystems is dependent of the explicit property of the reservoir and its coupling to the qubit . 
   the rich dynamical behaviors obtained previously in the literature can be regarded as the special cases of our present result or markovian approximation . 
   particularly , we find that , instead of entirely transferred to the reservoirs , the entanglement can be stably distributed among all the bipartite subsystems if the qubit and its reservoir can form a bound state and the non - markovian effect is important , and the esd of the qubits is not always accompanied with the occurrence of esb of reservoirs . 
   irrespective of how the entanglement distributes , it is found that the identity about the entanglement in the whole system can be satisfied at any time , which reveals the profound physics of the entanglement dynamics under decoherence . 
   this paper is organized as follows . in sec . 
   [ model ] , the model of two independent qubits in two local reservoirs is given . and 
   the dynamical entanglement invariance is obtained based on the exact solution of the non - markovian decoherence dynamics of the qubit system . in sec . 
   [ edd ] , the entanglement distribution over the subsystems when the reservoirs are pbg mediums is studied explicitly . 
   a stable entanglement - distribution configuration is found in the non - markovian dynamics . 
   finally , a brief discussion and summary are given in sec . 
   we consider two qubits interacting with two uncorrelated vacuum reservoirs . due to the dynamical independence between the two local subsystems 
   , we can firstly solve the single subsystem , then apply the result obtained to the double - qubit case . 
   the hamiltonian of each local subsystem is @xcite @xmath0 where @xmath1 and @xmath2 are the inversion operators and transition frequency of the qubit , @xmath3 and @xmath4 are the creation and annihilation operators of the @xmath5-th mode with frequency @xmath6 of the radiation field . 
   the coupling strength between the qubit and the reservoir is denoted by @xmath7 , where @xmath8 and @xmath9 are the unit polarization vector and the normalization volume of the radiation field , @xmath10 is the dipole moment of the qubit , and @xmath11 is the free space permittivity . 
   for such a system , if the qubit is in its ground state @xmath12 and the reservoir is in vacuum state at the initial time , then the system does not evolve to other states . 
   when the qubit is in its excited state @xmath13 , the system evolves as @xmath14 here @xmath15 denotes that the qubit jumps to its ground state and one photon is excited in the @xmath5-th mode of the reservoir . 
   @xmath16 satisfies an integro - differential equation @xmath17 where the kernel function @xmath18 is dependent of the spectral density @xmath19 . introducing the normalized collective state of the reservoir with one excitation as @xmath20 and with zero excitation as @xmath21 @xcite , eq . 
   ( [ t9 ] ) can be written as @xmath22 , where @xmath23 . 
   it should be emphasized that the introducing of normalized collective state is not a reduction of present model to the j - c model @xcite , as noted in @xcite . 
   the dynamics is given by eq . 
   ( [ t7 ] ) , which is difficult to obtain analytically since its non - markovian nature . in general the numerical integration should be used . 
   it is emphasized that our treatment to the dynamics of the system is exact without resorting to the widely used born - markovian approximation . to compare with the conventional approximate result 
   , we may derive straightforwardly the master equation from eq . 
   ( [ t9 ] ) after tracing over the degree of freedom of the reservoir @xcite , @xmath24+\gamma ( t)[2\sigma _ { -}\rho ( t)\sigma _ { + }   \nonumber \\ & & -\sigma _ { + } \sigma _ { -}\rho ( t)-\rho ( t)\sigma _ { + } \sigma _ { - } ] , \label{mstt}\end{aligned}\]]where the time - dependent parameters are given by @xmath25,~\gamma ( t)=-\text{re}[\frac{% \dot{b}(t)}{b(t)}]$ ] . 
   the time - dependent parameters @xmath26 and @xmath27 play the roles of lamb shifted frequency and decay rate of the qubit , respectively . 
   the integro - differential equation ( [ t7 ] ) contains the memory effect of the reservoir registered in the time - nonlocal kernel function and thus the dynamics of qubit displays non - markovian effect . 
   if the time - nonlocal kernel function is replaced by a time - local one , then eq . ( [ mstt ] ) recovers the conventional master equation under born - markovian approximation @xcite .    according to the above results , 
   the time evolution of a system consisted of two such subsystems with the initial state @xmath28 is given by @xmath29 where @xmath30 and @xmath31 are the coefficients to determine the initial entanglement in the system . 
   from @xmath32 , one can obtain the time - dependent reduced density matrix of the bipartite subsystem qubit1-qubit2 ( @xmath33 ) by tracing over the reservoir variables . 
   it reads @xmath34 where @xmath35 and @xmath36 . 
   similarly , one can obtain the corresponding reduced density matrices for other subsystems like reservoir1-reservoir2 ( @xmath37 ) and qubit - reservoir ( @xmath38 , @xmath39 , @xmath40 , @xmath41 ) .    using the concurrence @xcite to quantify entanglement 
   , we can calculate the entanglement of each subsystem as @xmath42 with @xmath43 for different bipartite partitions labeled by @xmath44 as @xmath45 one can verify that @xmath46 in eqs . 
   ( [ qq])-([qr2 ] ) satisfy an identity @xmath47 where @xmath48 is just the initial entanglement present in @xmath33 . 
   ( [ t15 ] ) recovers the explicit form derived in a double j - c model @xcite when each of the reservoirs contains only one mode , i.e. @xmath49 , where the decoherence is absent and the dynamics is reversible . 
   it is interesting that this identity is still valid in the present model because the reservoirs containing infinite degrees of freedom here lead to a completely out - of - phase interaction with qubit and an irreversibility . 
   furthermore , one notes that the identity is not dependent of any detail about @xmath16 , which only determines the detailed dynamical behavior of each components in eq . 
   ( [ t15 ] ) . 
   this result manifests certain kind of invariant nature of the entanglement . 
   ( [ t15 ] ) can be intuitively understood by the global multipartite entanglement of the whole system . 
   the global entanglement carried by the subsystem @xmath50 can be straightforwardly calculated from eq . 
   ( [ phit ] ) by generalized concurrence @xcite as @xmath48 , which , coinciding with the bipartite entanglement initially present in @xmath33 , just is the right hand side of eq . 
   ( [ t15 ] ) . 
   since there is no direct interaction between @xmath51 and @xmath52 , this global entanglement is conserved during the time evolution . from this point , 
   our result is consistent with the one in refs . @xcite and @xcite . 
   another observation of eq . 
   ( [ t15 ] ) is that the different coefficients in the left hand side are essentially determined by the energy / information transfer among the local subsystems . 
   explicitly , in our model the total excitation number is conserved , so the energy degradation in @xmath53 with factor @xmath16 is compensated by the energy enhancement in @xmath54 with factor @xmath55 . 
   this causes that @xmath56 , @xmath57 , and @xmath58 , in all of which the double excitation is involved , have similar form except for the different combinations of @xmath16 and @xmath55 in eqs . 
   ( [ qq ] ) , ( [ rr ] ) , and ( [ qr2 ] ) . 
   the dynamical consequence of the competition of the two terms in these equations causes the sudden death / birth of entanglement characterized by the presence of negative @xmath59 . 
   a different case happens for @xmath60 , where only single excitation is involved and no sudden death is present . with these observation , one can roughly understand why such combination in left hand side of eq . 
   ( [ t15 ] ) gives the global entanglement .    the significance of eq . 
   ( [ t15 ] ) is that it gives us a guideline to judge how the entanglement spreads out over all the bipartite partitions . 
   it implies that entanglement is not destroyed but re - distributed among all the bipartite subsystems and this re - distribution behavior is not irregular but in certain kind of invariant manner . 
   the similar invariant property of entanglement evolution has also been studied in ref . 
   @xcite .    in the following 
   we explicitly discuss the entanglement distribution , especially in the steady state , by taking the reservoir as a photonic band gap ( pbg ) medium @xcite and compare it with the previous results . 
   we will pay our attention mainly on the consequence of the non - markovian effect on the entanglement distribution and its differences to the results in refs . @xcite and @xcite . 
   for the pbg medium , the dispersion relation near the upper band - edge is given by @xcite @xmath61 where @xmath62 , @xmath63 is the upper band - edge frequency and @xmath64 is the corresponding characteristic wave vector . in this case 
   , the kernel function has the form @xmath65 where @xmath66 is a dimensionless constant . in solving eq . 
   ( [ t7 ] ) for @xmath16 , eq . 
   ( [ kn ] ) is evaluated numerically . 
   here we do not assume that @xmath5 is replaced by @xmath64 outside of the exponential @xcite . 
   so our result is numerically exact . in the following 
   we take @xmath67 as the unit of frequency .    [ 
   cols="^,^ " , ]     in figs . 
   [ bst ] and [ nbst ] , we show the entanglement evolutions of each subsystem for two typical cases of @xmath68 and @xmath69 , which correspond to the atomic frequency being located at the band gap and at the upper band of the pbg medium , respectively . in the both cases the initial entanglement in @xmath33 begins to transfer to other bipartite partitions with time but their explicit evolutions , in particular the long time behaviors , 
   are quite different . in the former case 
   , the entanglement could be distributed stably among all possible bipartite partitions . 
   [ bst](a ) shows that after some oscillations , a sizeable entanglement of @xmath33 is preserved for the parameter regime of @xmath70 . 
   remarkably , the entanglement in @xmath71 forms quickly in the full range of @xmath30 [ fig . [ bst](c ) ] and dominates the distribution . on the contrary , only slight entanglement of @xmath37 is formed in a very narrow parameter regime @xmath72 , as shown in fig . 
   [ bst](b ) . 
   however , when @xmath73 is located at the upper band of the pbg medium , the initial entanglement in @xmath33 is transferred completely to the @xmath37 in the long - time limit , as shown in fig . 
   [ nbst ] . at the initial stage , @xmath74 and @xmath75 
   are entangled transiently , but there is no stable entanglement distribution . 
   this result is consistent with that in refs . 
   it is noted that the entanglement in @xmath76 comes from two parts : one is transferred from @xmath33 , the other is created by the direct interaction between @xmath53 and @xmath54 . this can be seen clearly from fig . 
   [ bst](c ) when @xmath30 is very small , the initial entanglement of @xmath33 is very small , while that of @xmath76 is rather large , which just results from the interaction between @xmath53 and @xmath54 . 
   another interesting point is a stable entanglement can even be formed for the non - interacting bipartite system @xmath39 [ see [ bst](d ) when @xmath77 . 
   this entanglement transfer also results from the local interaction between @xmath53 and @xmath54 @xcite . 
   it is not difficult to understand these rich behaviors of entanglement distribution according to eqs . 
   ( [ qq])-([qr2 ] ) and its invariance ( [ t15 ] ) . from these equations 
   , one can clearly see that the entanglement dynamics and its distributions in the bipartite partitions are completely determined by the time - dependent factor @xmath78 of single - qubit excited - state population . 
   [ tong0 ] shows its time evolutions for the corresponding parameter regimes presented above . 
   we notice that @xmath79 when @xmath73 is located at the band gap , which means that there is some excited - state population in the long - time limit . 
   this phenomenon known as population trapping @xcite is responsible for the suppression of the spontaneous emission of two - level system in pbg reservoir and has been experimentally observed @xcite . 
   such population trapping just manifests the formation of bound states between @xmath53 and @xmath54 @xcite , which has been experimentally verified in @xcite . 
   consequently , @xmath53 and @xmath54 are so correlated in the bound states that the initial entanglement in @xmath33 can not be fully transferred to @xmath37 . 
   the oscillation during the evolution is just the manifestation of the strong non - markovian effect induced by the reservoirs . 
   on the contrary , if @xmath73 is located in the upper band , then @xmath80 and the qubits decay completely to their ground states . in this case 
   the bound states between @xmath53 and @xmath54 are absent and , according to eq . 
   ( [ t15 ] ) , the initial entanglement in @xmath33 is completely transferred to @xmath37 , as clearly shown in eq . 
   ( [ rr ] ) . 
   in addition , in refs . @xcite it was emphasized that esd of @xmath33 is always accompanied with esb of @xmath37 . 
   however , this is not always true . 
   to clarify this , we examine the condition to obtain esd of the qubits and the companying esb of the reservoirs . from eqs . 
   ( [ qq ] ) and ( [ rr ] ) it is obvious that the condition is @xmath81 and @xmath82 at any @xmath83 and @xmath84 , which means @xmath85 when the bound states is absent , @xmath86 , the condition ( [ dd ] ) can be satisfied when @xmath87 . 
   so one can always expect esd of the qubits and the companying esb of the reservoirs in the region @xmath88 , as shown in fig . 
   [ nbst ] and refs . 
   however , when the bound states are available , the situation changes . 
   in particular , when @xmath89 in the full range of time evolution , no region of @xmath90 can make the condition ( [ dd ] ) to be satisfied anymore . 
   for clarification , we present three typical behaviors of the entanglement distribution in fig . 
   [ td ] . in all these cases 
   the bound states are available . 
   [ td](a ) shows the situation where the entanglement is stably distributed among all of the bipartite subsystems . 
   [ td](b ) indicates that the entanglement of @xmath91 shows esb and revival , while the entanglement of @xmath33 does not exhibit esd . 
   [ td](c ) shows another example that while the entanglement of @xmath33 has esd and revival @xcite , the entanglement of @xmath37 does not show esb but remains to be zero . 
   [ td](b , c ) reveal that esd in @xmath33 has no direct relationship with esb in @xmath37 . 
   the above discussion is not dependent of the explicit spectral density of the individual reservoir . to confirm this 
   , we consider the reservoir in free space . 
   the spectral density has the ohmic form @xmath92 , which can be obtained from the free - space dispersion relation @xmath93 . 
   one can verify that the condition for the formation of bound states is : @xmath94 @xcite . in fig . 
   [ osd ] , we plot the results in this situation . 
   the previous results can be recovered when the bound states are absent @xcite . on the contrary , when the bound states are available , a stable entanglement is established among all the bipartite partitions . 
   therefore , we argue that the stable entanglement distribution resulted from the bound states is a general phenomenon in open quantum system when the non - markovian effect is taken into account .    in summary , we have studied the entanglement distribution among all the bipartite subsystems of two qubits embedded into two independent amplitude damping reservoirs . 
   it is found that the entanglement can be stably distributed in all the bipartite subsystems , which is much different no matter to the markovian approximate result @xcite or to the decoherenceless double j - c model result @xcite , and an identity about the entanglement in all subsystems is always satisfied . 
   this identity is shown to be independent of any detail of the reservoirs and their coupling to the qubit , which affect only the explicit time evolution behavior and the final distribution . 
   the result is significant to the study of the physical nature of entanglement under decoherence . 
   it implies an active way to protect entanglement from decoherence by modifying the properties of the reservoir via the potential usage of the newly emerged technique , i.e. quantum reservoir engineering @xcite . 
   this work is supported by the fundamental research funds for the central universities under grant no . 
   lzujbky-2010 - 72 , gansu provincial nsf under grant no . 
   0803rjza095 , the national nsf of china , the program for ncet , and the cqt wbs grant no . r-710 - 000 - 008 - 271 .
summary: we study the entanglement dynamics of two qubits , each of which is embedded into its local amplitude - damping reservoir , and the entanglement distribution among all the bipartite subsystems including qubit - qubit , qubit - reservoir , and reservoir - reservoir . 
   it is found that the entanglement can be stably distributed among all components , which is much different to the result obtained under the born - markovian approximation by c. e. lpez _ 
   et al . _ 
   [ phys . 
   rev . 
   lett . 
   * 101 * , 080503 ( 2008 ) ] , and particularly it also satisfies an identity . 
   our unified treatment includes the previous results as special cases . 
   the result may give help to understand the physical nature of entanglement under decoherence .