11institutetext: Department of Computer Science, University of Copenhagen,
Denmark 11email<EMAIL_ADDRESS>
22institutetext: Department of Biomedical Sciences, University of Copenhagen,
Denmark

# Extremely weakly-supervised blood vessel segmentation with physiologically
based synthesis and domain adaptation

Peidi Xu 11 Olga Sosnovtseva 22 Charlotte Mehlin Sørensen 22 Kenny Erleben 11
Sune Darkner 11

###### Abstract

Accurate analysis and modeling of renal functions require a precise
segmentation of the renal blood vessels. Micro-CT scans provide image data at
higher resolutions, making deeper vessels near the renal cortex visible.
Although deep-learning-based methods have shown state-of-the-art performance
in automatic blood vessel segmentations, they require a large amount of
labeled training data. However, voxel-wise labeling in micro-CT scans is
extremely time-consuming given the huge volume sizes. To mitigate the problem,
we simulate synthetic renal vascular trees physiologically while generating
corresponding scans of the simulated trees by training a generative model on
unlabeled scans. This enables the generative model to learn the mapping
implicitly without the need for explicit functions to emulate the image
acquisition process. We further propose an additional segmentation branch over
the generative model trained on the generated scans. We demonstrate that the
model can directly segment blood vessels on real scans and validate our method
on both 3D micro-CT scans of rat kidneys and a proof-of-concept experiment on
2D retinal images. Code and 3D results are available at
111https://github.com/miccai2023anony/RenalVesselSeg

###### Keywords:

Blood vessel Renal vasculature Semantic segmentation Physiological simluation
Generative model Domain adaptation.

##  1 Introduction

The vasculature in each organ has a characteristic structure tailored to
fulfil the particular requirements of the organ. The renal vasculature, which
serves as a resource distribution network, plays a significant part in the
kidney’s physiology and pathophysiology. Not only does it distribute blood and
nutrients to individual nephrons and regulates the filtration of blood in the
kidney, but it also acts as a communication network, allowing neighboring
nephrons to interact through electrical signals transmitted along the vessels
[16]. Automatic segmentation of renal blood vessels from medical scans is
usually the essential first step for developing realistic computer simulations
of renal functions.

General deep-learning based segmentation of blood vessels: Deep learning
models have been widely used for automatic blood vessel segmentations and have
shown state-of-the-art performances applied on lungs, liver, and eyes [11, 5,
12, 27]. However, only a few efforts were made for renal blood segmentation.
Recently, He et al. proposed Dense biased networks [9] to segment renal
arteries from abdominal CT angiography (CTA) images by fusing multi-receptive
fields and multi-resolution features for the adaptation of scale changes.
However, the limited resolution of CTA images only allows the models to reach
interlobar arteries and enables the estimation of blood-feeding regions, which
is useful for laparoscopic partial nephrectomy but not for analyzing realistic
renal functions. Therefore, there is a need for imaging with higher
resolution, e.g., micro-computed tomography (micro-CT) scans.

Micro-computed tomography and related deep-learning works: Micro-CT shares the
same imaging principle as conventional clinical CT, but with a much smaller
field of view, such that microscale images of high spatial resolution can be
obtained [20]. Micro-CT scans are commonly used to study various
microstructures including blood vessels [1]. Few existing research on the
auto-segmentation of micro-CT scans focuses on segmenting internal organs of
the heart, spinal cord, right and left lung [15], and blood vessels on
colorectal tissue [19] using either nn-UNet [10] or variant 3DUNet [21, 4].
There is, however, no prior work in segmenting vasculatures in organs like
kidneys from micro-CT scans. Crucially, most of the above deep learning
methods require a large number of label maps to train the segmentation
network. Manual labeling of micro-CT scans is extremely time-consuming given
the huge volume size. Therefore, in our case, we do not have any clean label
maps to train a segmentation model.

Synthetic training data for blood vessel segmentation: Transfer learning from
artificially generated data is one possible technique to train deep learning
models in a data scarcity setting. The process involves pre-training models on
synthetic data, which are then fine-tuned on a small set of labeled real data.
In medical image segmentation, this strategy has been widely applied to tumor
segmentations [13]. Since blood vessels do follow certain physiological and
anatomical properties, e.g., Murray’s law [18], this approach has also been
applied to train segmentation models for mouse brain vasculature with
physiologically synthesized vessels [25]. However, these works only pre-train
the models on synthetic data and still require real labeled data for fine-
tuning. Recently, Menten et al. [17] synthesize retinal vasculature and then
emulate the corresponding optical coherence tomography angiography (OCTA)
images. They show that a UNet trained on these emulated image-label pairs can
directly segment real OCTA images. However, the way they generate scans from
synthesized labels is completely explicit, which includes a series of physics-
based image transformation functions that emulate the image acquisition
process (e.g., OCTA). These functions clearly require expert knowledge and do
not translate to micro-CT settings.

Generative models for domain adaptation: In practice, a relatively large
number of unlabeled scans are usually available. Thus, a more general way to
generate scans without emulating image acquisition explicitly is to utilize
these unlabeled scans via a generative model, i.e., domain adaptation. Medical
image segmentation with domain adaptation is an active research area, and a
popular method is Generative Adversarial Networks (GANs). In particular,
CycleGAN [29] has been used to perform domain adaptation for medical image
segmentation such as liver and tumor segmentation in CT images, where the goal
is to train a segmentation model on one domain and then apply it to another
domain with missing or scarce annotated data [2, 8] . Recently, Chen et al.
[3] segment cerebral vessels in 2D Laser speckle contrast imaging (LSCI)
images using public fundus images with segmentation labels as the source
domain and a CycleGAN for domain adaptation. The ability of CycleGAN to
perform translation without paired training data makes it a powerful tool for
domain adaptation.

Our contribution: We propose a framework with two main components: 1) a
physiology-based simulation that synthesizes renal vascular trees and 2) a
generative model with an additional segmentation branch that adapts the
synthesized vascular trees to real scans while performing segmentation
simultaneously. For 1), we extend the work [26] from physiologically
synthesizing renal arterial tree to venous tree. Since a small prebuilt tree
needs to be manually provided in the initialization step of the process, we
call our method extremely weakly-supervised. For 2), we aim to “emulate”
corresponding scans using CycleGANs. Specifically, we follow the idea in [3]
to train a vessel segmentation network over the output from CycleGAN while
extending to 3D on kidney micro-CT images. Notably, although 3D CycleGAN has
been adopted for segmenting brain tissues and heart chambers [2, 8, 28], no
similar work exists on subtle structures like blood vessels in 3D. Moreovoer,
these works still require scans or labels from other sources, modalities, or
time points as the source domain. Instead, our source domain is purely
physiologically synthesized vascular trees. We show that our combined model
can directly segment blood vessels in real scans from the target domain and
demonstrate the validity of our approach on segmenting vasculatures both on 3D
micro-CT scans of rat kidneys and 2D retinal images.

##  2 Method

###  2.1 Physiologically-based simulation of renal blood vessels

Figure 1: Physiologically based vessel synthesizing pipeline. Details are
given in [26]. The last two subfigures are shown with maximum intensity
projection (MIP).

Constraint Constructive Optimization (CCO) [22] and its variant Global
Constructive Optimization (GCO) [7] are widely used angiogenesis-based methods
that simulate the growth of vascular trees. These methods turn tree growth
into an optimization problem based on the biological, physiological, and
hemodynamic factors involved in the process. Here, the vascular tree is
modeled by a directed acyclic graph
$\mathcal{G}\equiv(\mathcal{V},\mathcal{E})$ where $\mathcal{V}$ is a set of
nodes in the two endpoints of each vessel centerline with its coordinates in
Euclidean space as node features, and $\mathcal{E}$ is a set of directed edges
representing each vessel segment as a cylindrical tube with its specific
radius as edge features. Boundary conditions such as terminal radius and flow
distributions are imposed to represent physiologic conditions. The algorithms
then find a tree that minimizes the system’s overall cost function while
fulfilling several constraints.

Here we follow [26] which adopts GCO as the backbone model for generating the
renal vascular trees with optimal branching structures by performing multi-
scale optimizations through iterating several operations such as splitting and
relaxation [7]. However, instead of the arterial tree presented in [26], we
only focus on the venous tree because it constitutes most of the vessel
foreground. In summary, veins follow a similar pattern to arteries but are
thicker, which is accomplished by sampling more terminal nodes with large
radii (Fig. 1a). Detailed modifications over boundary conditions from arterial
to venous trees are given in supplementary material, including tuning the
weighting factors of the loss defined in [26]. Together with inherent
randomnesses in the GCO process itself, the generated tree (Fig. 1d) will look
different each run. This enables a variety of synthesized vascular trees to
train the later deep-learning model.

Note that though a prebuilt tree
$\mathcal{G}_{0}\equiv(\mathcal{V}_{0},\mathcal{E}_{0})$ is required to guide
the GCO process as noted in [26], $\mathcal{G}_{0}$ involves less than 20
nodes, as shown in Fig. 1b, which can be manually selected. This node
selection process should take much less time than a voxel-wise manual
annotation of the whole blood vessels. Therefore, we call our pipeline
extremely weakly-supervised with a partially annotated tree structure but
without any real segmentation label maps.

###  2.2 CycSeg: CycleGAN with additional segmentation branch

To create a synthetic image dataset, the reconstructed vascular tree
structures $\mathcal{G}\equiv(\mathcal{V},\mathcal{E})$ (Fig. 1d) are then
remapped to a 3D binary label map (Fig. 1e) by voxelization, the detail of
which is given in [26]. We then generate the corresponding gray-scale
synthetic images (Fig. 1f) by simply assigning vessel foreground and
background with random integers in $[128,255]$ and $[0,127]$ respectively. Of
course, a segmentation model trained on these image-label pairs (Fig. 1f & e)
will not work on real scans because of this oversimplified scan construction.

In order to adapt corresponding scans that emulate the micro-CT acquisition
process out of the label maps from the previous step, unlabeled real micro-CT
scans are utilized to train a generative model. Our backbone model is CycleGAN
[29], which we extend to 3D while integrating an additional segmentation
branch with segmenter $S$ using a standard 3D UNet[4], as shown in Fig.
2(Left). We refer to our model as CycSeg later. The only modification we make
in 3D is that 32 filters are used in the first layer instead of 64 to ease the
computational load. Please refer to the supplementary material for detailed
model architectures.

We strictly follow the original design [29] while extending to 3D for the loss
functions definition and training procedure. Briefly, the two-way GAN loss
$\mathcal{L}_{\mathrm{GAN}}\left(G_{A\rightarrow B},D_{B}\right)$ and
$\mathcal{L}_{\mathrm{GAN}}\left(G_{B\rightarrow A},D_{A}\right)$ is to
encourage the generator network to produce samples that are indistinguishable
from real samples. Besides, cycle-consistency loss
$\mathcal{L}_{\mathrm{cyc}}(G_{A\rightarrow B},G_{B\rightarrow A})$ is to
assure the property ${a\sim p_{\text{data }}(a)}:G_{B\rightarrow
A}(G_{A\rightarrow B}(a))=a$ and ${b\sim p_{\text{data }}(b)}:G_{A\rightarrow
B}(G_{B\rightarrow A}(b))=b$. In our case, $p_{\text{data }}(a)$ and
$p_{\text{data }}(b)$ denotes the real micro-CT scans distribution and the
synthetic scans distribution from the physiologically generated trees,
respectively. A final identity loss $\mathcal{L}_{\text{id }}(G_{A\rightarrow
B},G_{B\rightarrow A})$ is to stabilize the two generators. Please refer to
[29] for the exact definition and computation of the above losses.

In our CycSeg, the additional segmenter ($S$) is trained on fake A (output
from $G_{B\rightarrow A}$) with the goal to work on real A, as shown in Fig.
2. This introduces the segmentation loss
$\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow A})$. Thus, the final loss
$\mathcal{L}_{\text{tot}}$ is defined as

$\displaystyle\mathcal{L}_{\text{tot}}$
$\displaystyle=\mathcal{L}\left(G_{A\rightarrow B},G_{B\rightarrow
A},D_{A},D_{B},S\right)$ (1)
$\displaystyle=\mathcal{L}_{GAN}\left(G_{A\rightarrow
B},D_{B}\right)+\mathcal{L}_{GAN}\left(G_{B\rightarrow A},D_{A}\right)$
$\displaystyle+\lambda_{1}\mathcal{L}_{cyc}(G_{A\rightarrow B},G_{B\rightarrow
A})+\lambda_{2}\mathcal{L}_{\text{id}}(G_{A\rightarrow B},G_{B\rightarrow
A})+\lambda_{3}\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow A})$

Here $\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow A})$ is an unweighted
combination of a dice loss and standard cross-entropy loss. Specifically, the
segmenter $S$ takes the output from the generator $G_{B\rightarrow A}$ as
input. Thus, given the physiologically generated label map $g$ and its
corresponding synthetic gray-scale image $x\sim p_{\text{data }}(b)$, the
segmenter outputs $p=S(G_{B\rightarrow A}(x))$, and
$\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow A})$ is defined as

$\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow
A})=\frac{1}{N}\left(-\sum\nolimits_{i}^{N}g_{i}\log
p_{i}+1-\frac{2\sum_{i}^{N}p_{i}g_{i}}{\sum_{i}^{N}p_{i}+\sum_{i}^{N}g_{i}}\right)$
(2)

where $N$ is the total number of voxels in each sampled 3D patch.

We follow [29] by setting the weights $\lambda_{1}=10$ and $\lambda_{2}=5$
while setting $\lambda_{3}=3$ experimentally. Note that although
$G_{B\rightarrow A}$ is one of the input models to compute
$\mathcal{L}_{\text{seg}}(S,G_{B\rightarrow A})$, all the CycleGAN components
including $G_{B\rightarrow A}$ are frozen when training the segmenter $S$ by
blocking backpropagation. Moreover, since the generator, discriminator, and
segmenter are trained alternately, $\lambda_{3}$ does not strongly impact the
training and only affects early stopping.

During inference, all CycleGAN components are discarded, while real scans
(domain A) are directly passed to segmenter $S$ to output segmentation maps.

Figure 2: Left: An illustration of CycleGAN with an additional segmentation
branch working on adapted (fake) domain A images by $G_{B\rightarrow A}$. All
3D patches are shown with MIP. Computations for $\mathcal{L}_{\text{GAN}}$,
$\mathcal{L}_{\text{id}}$, $\mathcal{L}_{\text{cyc}}$ are not shown for
simplicity. Middle: An example result in 3D. Right: A sample slice overlaid
with segmentation.

##  3 Experiments and Results

Dataset: The kidney cast was prepared as described in [1] in agreement with
approved protocols (approval granted from the Danish Animal Experiments
Inspectorate under the Ministry of Environment and Food, Denmark). The rat
kidneys were ex vivo scanned in a ZEISS XRadia 410 Versa $\mu$CT scanner with
an isotropic voxel size of 22.6 $\mu$m [1], with a fixed dimension of
$1000\times 1024\times 1014$. To ease the computational load, scans are auto-
cropped to $(692\pm 33)\times(542\pm 28)\times(917\pm 28)$ by intersected
bounding cubes from Otsu’s thresholding in each dimension [26]. Here, we use 7
unlabeled scans (domain A) for training and 1 labeled scan for testing, while
other 3 unlabeled scans are only validated visually. The synthesized dataset
(domain B) has 15 image-label pairs by tuning parameters used in GCO and
running multiple times as discussed in Section 2.1.

Pre-processing: Each generated patch is only preprocessed by simple min-max
normalization $X_{\text{scale }}=(x_{i}-x_{\text{min}})/(x_{max}-x_{min})$.

Experimental Setup and Training process: The network is implemented in PyTorch
and trained on NVIDIA A100 with a batch size of 1 and patch size of 208 for
200 epochs. All three components are optimized using the Adam optimizer with
the same learning rate of $2\times 10^{-4}$ and reduced by 1% for every epoch.
We apply early stopping if $\mathcal{L}_{\text{tot}}$ (cf. Eq.(1)) of ten
consecutive epochs does not decrease. Training takes approximately three days
to reach convergence, while segmentation during inference takes only around
two minutes per scan.

Table 1: Segmentation result on private renal and public retina dataset (in $\%$).  Training data | Model | Acc | DICE | clDICE
---|---|---|---|---
Renal | Synthetic label + Raw scan | CycSeg | 99.2$\pm$0.0 | 76.8$\pm$0.3 | 65.0$\pm$0.8
CycSegP | 99.2$\pm$0.0 | 77.6$\pm$0.2 | 72.4$\pm$1.0
Retinal | CHASE label + DRIVE image | CycSeg | 96.1$\pm$0.1 | 74.8$\pm$0.4 | 76.3$\pm$0.1
DRIVE label + DRIVE image | UNet2D | 96.5$\pm$0.1 | 79.6$\pm$0.5 | 78.9$\pm$0.8

###  3.1 Results

As shown in Fig. 2, the CycleGAN successfully adapts realistic noise during
micro-CT acquisition to the synthesized images (from real B to fake A), while
a UNet trained over the adapted image with the corresponding synthesized label
maps can segment real micro-CT scans. Although the ex-vivo scan separates the
organ of interest from other parts, the segmentation of venous trees is still
challenging due to various noises during the micro-CT acquisition and efficacy
of the contrast media. Despite of few false positives, the 3D results in Fig.
2 show a smooth and clear venous tree structure, which indicates that the
segmenter is trained to recognize veins from noises in real scans by the
adapted images from $G_{B\rightarrow A}$. Since the result should anatomically
be a connected tree, we also apply a simple connected component post-
processing to the model’s output (CycSegP), which successfully removes most of
the floating points and gives better results.

Table 1 shows the numerical evaluation of the segmentation performance on one
test image including accuracy, DICE score, and the topology-aware centerline
DICE (clDICE) [23]. Although we cannot compare our results with standard
segmentation models due to the lack of any labeled training data, we believe
that the visual inspection in Fig. 2 and relatively high quantitative results
in Table 1 over 3D vessel segmentation task demonstrate the potential of our
method in segmenting and building 3D renal models from micro-CT scans.

###  3.2 Proof-of-concept on retinal blood vessel segmentation

We acknowledge that the quantitative analysis above is not thorough with only
one test data. Thus, we conduct a proof-of-concept on a 2D retinal blood
vessel segmentation task. We follow the previous experimental setups but with
a patch size of 256 and batch size of 32 due to the ease of computational load
in 2D.

Dataset and domain construction: We adopt the popular DRIVE dataset [24] as
the target domain, which includes 40 digital fundus images captured by a Canon
CR5 3CCD camera with a resolution of $584\times 565$ pixels. However, the 3D
vascular structure will not be physiologically correct when projecting to 2D
images, and some works argue against compliance with Murray’s law in retinal
blood vessels [14]. Thus, we do not synthesize retinal blood vessels using a
similar physiologically based strategy in kidneys and only focus on the
validity of the segmentation power with domain adaptation. Thus, we directly
adopt label maps from another dataset (CHASE [6]) as the fake source domain.
Specifically, we use 20 from 40 images of DRIVE dataset without label maps,
together with 28 label maps from CHASE dataset without corresponding images to
jointly train the CycSeg. The remaining 20 labeled images from DRIVE are used
for testing. Note that we only work on gray-scale images, and label maps from
CHASE are resized to $584\times 565$ using nearest neighbor interpolation.

Figure 3: Segmentation results on a DRIVE test image of different strategies.

As shown in Fig. 3, CycSeg trained on unpaired images is able to segment blood
vessels on real images from DRIVE directly. An illustration of retinal domain
adaptation in analogue to Fig. 2(Left) is in the supplementary material. As an
ablation study, we adopt the same train-test split to directly train a UNet
[21] using real image-label pairs from the DRIVE training set. From Fig. 3,
both UNet and CycSeg can segment large vessels well while (especially CycSeg)
having difficulty detailing small vessels. Numerical results in Table 1 show
that the performance of CycSeg trained purely on adapted images is still
lower, but the difference is acceptable, as our goal is not to outperform the
state-of-the-art but to propose a pipeline that can do segmentation without
any labeled images from the target domain. Future work would be to train a
UNet using the same generative strategy together with real image-label pairs
to see how much information gains it can offer compared with UNet trained only
on real image-label pairs.

##  4 Conclusion

We have presented a pipeline that segments blood vessels from real scans
without any manually segmented training data. The pipeline first synthesizes
label maps of renal vasculatures using a physiologically based model and then
generates corresponding scans via a 3D CycleGAN. The CycleGAN emulates image
features and artifacts of micro-CT acquisition processes implicitly from
unlabeled scans. Simultaneously, an additional segmentation branch on top of
CycleGAN trained over the adapted scans enables the segmentation of blood
vessels on real scans. This removes any need for expert knowledge of scanning
settings like proposed in [17]. We believe our pipeline can crucially reduce
annotations needed to segment blood vessels and easily adapt to other organs
or modalities and enable computerized diagnosis of vessel-related diseases in
clinical practice.

Since segmentation is the final objective, the intermediate image adaptation
task is only visually inspected. Future work could include numerical tests
like image structure clustering and Amazon Mechanical Turk perceptual studies
[29]. A modification over [26] is necessary to model the pair-wise coupling of
veins and arteries. This will enable physiologically correct label maps for an
Artery/Vein multi-class segmentation task in the future, which is far beyond
the scope of this work. Nonetheless, more clean ground truth labels should
enable a more thorough validation and benchmark test with standard
segmentation models.

## References

  * [1] Andersen, S.B., Taghavi, I., Kjer, H.M., Søgaard, S.B., Gundlach, C., Dahl, V.A., Nielsen, M.B., Dahl, A.B., Jensen, J.A., Sørensen, C.M.: Evaluation of 2d super-resolution ultrasound imaging of the rat renal vasculature using ex vivo micro-computed tomography. Scientific reports 11(1), 1–13 (2021)
  * [2] Bui, T.D., Wang, L., Lin, W., Li, G., Shen, D.: 6-month infant brain mri segmentation guided by 24-month data using cycle-consistent adversarial networks. In: 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI). pp. 359–362. IEEE (2020)
  * [3] Chen, H., Shi, Y., Bo, B., Zhao, D., Miao, P., Tong, S., Wang, C.: Real-time cerebral vessel segmentation in laser speckle contrast image based on unsupervised domain adaptation. Frontiers in Neuroscience p. 1523 (2021)
  * [4] Çiçek, Ö., Abdulkadir, A., Lienkamp, S.S., Brox, T., Ronneberger, O.: 3d u-net: learning dense volumetric segmentation from sparse annotation. In: International conference on medical image computing and computer-assisted intervention. pp. 424–432. Springer (2016)
  * [5] Cui, H., Liu, X., Huang, N.: Pulmonary vessel segmentation based on orthogonal fused u-net++ of chest ct images. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part VI 22. pp. 293–300. Springer (2019)
  * [6] Fraz, M.M., Remagnino, P., Hoppe, A., Uyyanonvara, B., Rudnicka, A.R., Owen, C.G., Barman, S.A.: An ensemble classification-based approach applied to retinal blood vessel segmentation. IEEE Transactions on Biomedical Engineering 59(9), 2538–2548 (2012)
  * [7] Georg, M., Preusser, T., Hahn, H.K.: Global constructive optimization of vascular systems (2010), https://openscholarship.wustl.edu/cse˙research/36
  * [8] Gilbert, A., Marciniak, M., Rodero, C., Lamata, P., Samset, E., Mcleod, K.: Generating synthetic labeled data from existing anatomical models: an example with echocardiography segmentation. IEEE Transactions on Medical Imaging 40(10), 2783–2794 (2021)
  * [9] He, Y., Yang, G., Yang, J., Chen, Y., Kong, Y., Wu, J., Tang, L., Zhu, X., Dillenseger, J.L., Shao, P., et al.: Dense biased networks with deep priori anatomy and hard region adaptation: Semi-supervised learning for fine renal artery segmentation. Medical image analysis 63, 101722 (2020)
  * [10] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18(2), 203–211 (2021)
  * [11] Jia, D., Zhuang, X.: Learning-based algorithms for vessel tracking: A review. Computerized Medical Imaging and Graphics 89, 101840 (2021)
  * [12] Keshwani, D., Kitamura, Y., Ihara, S., Iizuka, S., Simo-Serra, E.: Topnet: Topology preserving metric learning for vessel tree reconstruction and labelling. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part VI 23. pp. 14–23. Springer (2020)
  * [13] Lindner, L., Narnhofer, D., Weber, M., Gsaxner, C., Kolodziej, M., Egger, J.: Using synthetic training data for deep learning-based gbm segmentation. In: 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). pp. 6724–6729. IEEE (2019)
  * [14] Luo, T., Gast, T.J., Vermeer, T.J., Burns, S.A.: Retinal vascular branching in healthy and diabetic subjects. Investigative ophthalmology & visual science 58(5), 2685–2694 (2017)
  * [15] Malimban, J., Lathouwers, D., Qian, H., Verhaegen, F., Wiedemann, J., Brandenburg, S., Staring, M.: Deep learning-based segmentation of the thorax in mouse micro-ct scans. Scientific Reports 12(1), 1822 (2022)
  * [16] Marsh, D.J., Postnov, D.D., Sosnovtseva, O.V., Holstein-Rathlou, N.H.: The nephron-arterial network and its interactions. American Journal of Physiology-Renal Physiology 316(5), F769–F784 (2019)
  * [17] Menten, M.J., Paetzold, J.C., Dima, A., Menze, B.H., Knier, B., Rueckert, D.: Physiology-based simulation of the retinal vasculature enables annotation-free segmentation of oct angiographs. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part VIII. pp. 330–340. Springer (2022)
  * [18] Murray, C.D.: The physiological principle of minimum work: I. the vascular system and the cost of blood volume. Proceedings of the National Academy of Sciences 12(3), 207–214 (1926)
  * [19] Ohnishi, T., Teplov, A., Kawata, N., Ibrahim, K., Ntiamoah, P., Firat, C., Haneishi, H., Hameed, M., Shia, J., Yagi, Y.: Three-dimensional vessel segmentation in whole-tissue and whole-block imaging using a deep neural network: Proof-of-concept study. The American Journal of Pathology 191(3), 463–474 (2021)
  * [20] Ritman, E.L.: Current status of developments and applications of micro-ct. Annual review of biomedical engineering 13, 531–552 (2011)
  * [21] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [22] Schreiner, W., Karch, R., Neumann, F., Neumann, M.: Constrained constructive optimization of arterial tree models. Scaling in biology 145, 65 (2000)
  * [23] Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P., Bauer, U., Menze, B.H.: cldice-a novel topology-preserving loss function for tubular structure segmentation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 16560–16569 (2021)
  * [24] Staal, J., Abràmoff, M.D., Niemeijer, M., Viergever, M.A., Van Ginneken, B.: Ridge-based vessel segmentation in color images of the retina. IEEE transactions on medical imaging 23(4), 501–509 (2004)
  * [25] Todorov, M.I., Paetzold, J.C., Schoppe, O., Tetteh, G., Shit, S., Efremov, V., Todorov-Völgyi, K., Düring, M., Dichgans, M., Piraud, M., et al.: Machine learning analysis of whole mouse brain vasculature. Nature methods 17(4), 442–449 (2020)
  * [26] Xu, P., Holstein-Rathlou, N.H., Søgaard, S.B., Gundlach, C., Sørensen, C.M., Erleben, K., Sosnovtseva, O., Darkner, S.: A hybrid approach to full-scale reconstruction of renal arterial network. Scientific Reports 13(1), 7569 (2023)
  * [27] Zhang, J., Zhang, Y., Xu, X.: Pyramid u-net for retinal vessel segmentation. In: ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). pp. 1125–1129. IEEE (2021)
  * [28] Zhang, Z., Yang, L., Zheng, Y.: Translating and segmenting multimodal medical volumes with cycle-and shape-consistency generative adversarial network. In: Proceedings of the IEEE conference on computer vision and pattern Recognition. pp. 9242–9251 (2018)
  * [29] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE international conference on computer vision. pp. 2223–2232 (2017)

Supplementary Material

Figure 4: Process over 2D Retinal images. Real Domain B images are simply generated from Domain B labels by assigning vessel foreground and background with random values in $[0,127]$ and $[128,255]$ respectively. Table 2: 3D Generator structure. pad and out_pad represent the amount of padding applied in all three dimensions to the input and the output, respectively. Please refer to Github Repo for details such as the usage of InstanceNorm3d and activation functions. Layer | Output feature map size
---|---
Input | $1\times 208\times 208\times 208$
$32\times 7\times 7\times 7$ conv, stride 1, pad 3 | $32\times 208\times 208\times 208$
$64\times 3\times 3\times 3$ conv, stride 2, pad 1 | $64\times 104\times 104\times 104$
$128\times 3\times 3\times 3$ conv, stride 2, pad 1 | $128\times 52\times 52\times 52$
3 consecutive Residual Blocks $128\times 3\times 3\times 3$ conv | $128\times 52\times 52\times 52$
$64\times 3\times 3\times 3$ conv, stride 2, pad 1 out_pad 1 | $64\times 104\times 104\times 104$
$32\times 3\times 3\times 3$ conv, stride 2, pad 1 out_pad 1 | $32\times 208\times 208\times 208$
$1\times 7\times 7\times 7$ conv, stride 1, pad 3 | $1\times 208\times 208\times 208$
Table 3: 3D Discriminator structure. Layer | Output feature map size
---|---
Input | $1\times 208\times 208\times 208$
$32\times 4\times 4\times 4$ conv, stride 2, pad 1 | $32\times 104\times 104\times 104$
$64\times 4\times 4\times 4$ conv, stride 2, pad 1 | $64\times 52\times 52\times 52$
$128\times 4\times 4\times 4$ conv, stride 2, pad 1 | $128\times 26\times 26\times 26$
$256\times 4\times 4\times 4$ conv, stride 1, pad 1 | $256\times 25\times 25\times 25$
$1\times 4\times 4\times 4$ conv, stride 1, pad 1 | $1\times 24\times 24\times 24$
Table 4: Boundary conditions and physics parameters involved in the GCO process for constructing renal venous trees. Note that $Q_{t}$ is derived by equal terminal flow assumption: $Q_{t}=\frac{Q}{N}\approx\frac{7ml/min}{68564}\approx 1.7\times 10^{6}\mu m^{3}\,s^{-1}$. Other parameters are identical from [2]. Parameter | Meaning | Value/Sampling distribution
---|---|---
$r_{0}$ | radius of terminal veins | $\sim\mathcal{N}(10.79,2.41)$
$N$ | number of terminal veins | $\sim\mathcal{N}(68564,16647)$
$Q$ | inlet flow to renal venous tree | $7ml/min=1.167\times 10^{11}\mu m^{3}\,s^{-1}$
$Q_{t}$ | terminal flow out of renal venous tree | $1.7\times 10^{6}\mu m^{3}\,s^{-1}$
$\mu$ | viscosity | $3.6\times 10^{-3}\mathrm{~{}Pa}\,\mathrm{~{}s}$
$w_{c}$ | weight factor for material cost | $\\{5,6,7\\}\times 10^{-8}N\,{\mu m^{-2}}s^{-1}$


# Modeling the effects of perturbations and steepest entropy ascent on the
time evolution of entanglement

Cesar Damian  Department of Mechanical Engineering, Universidad de Guanajuato,
Salamanca, GTO 36885, Mexico<EMAIL_ADDRESS>(sabbatical at Virginia Tech)
Robert Holladay  Mechanical Engineering Department, Virginia Tech, Blacksburg,
VA 24061  Adriana Saldana  Agricultural Engineering Department, Universidad de
Guanajuato, Salamanca, GTO 36885, Mexico  Michael von Spakovsky
<EMAIL_ADDRESS>Mechanical Engineering Department, Virginia Tech, Blacksburg,
VA 24061

###### Abstract

This work presents an analysis of the evolution of perturbed Bell diagonal
states using the equation of motion of steepest-entropy-ascent quantum
thermodynamics (SEAQT), the Lindblad equation, and various measures of loss of
entanglement. First, a brief derivation is presented showing that Bell
diagonal states are stationary states that are not stable equilibrium states
relative to the SEAQT equation of motion, highlighting the need for the
development of perturbation methods to study the evolutions of nearby states.
This contrasts with the Lindblad equation of motion for which only some of the
Bell diagonal states are stationary. Next, two perturbation methods are
presented. The first is a general method for perturbing bipartite systems and
the second is a method based on a set of unitary operations that are
constrained to hold the system energy and system entropy constant. Sets of
density operators are randomly generated with each method and the resulting
time-varying characteristics of the system’s entanglement are analyzed. The
findings reveal that the evolutions associated with the constrained
perturbations accurately predict the loss of non-locality and align well with
the measured concurrence.

††preprint: APS/123-QED

##  I Introduction

Quantum information and computing systems rely on the quantum phenomena of
superposition and entanglement to implement useful tasks. In general, two
qubits may become entangled if they interact, and certain applications may
require that this entanglement persist for some specified amount of time even
if the qubits are no longer interacting (such as in the case of quantum
teleporting [1]). Thus, as dissipation occurs throughout a process, it is
necessary to understand the effects it has on the entangled state of the
system. While various dissipative models predict different effects on the
evolution of system state [2], this work focuses on two dynamics of system
state evolution. The first is based on the principle of steepest entropy
ascent (SEA) as implemented in the SEA quantum thermodynamic (SEAQT) equation
of motion. The second is the linear dynamics modeled by the Lindblad equation
of motion for open quantum systems.

Liu $et\;al.$ [3] theoretically predict the loss of entanglement between two
non-interacting qubits in a Bell diagonal state using the Kraus-Operator-Sum
approach. This approach does not assume a specific physical realization for
the qubits but does assume that both of the qubits interact with an external
global dephasing (decoherence) field, which is the cause of the dissipation
and loss of entanglement. They show that the concurrence (a measure of
entanglement) of the system remains unaffected by the dephasing field-induced
dissipation, resulting in what is referred to as time-invariant entanglement.
However, they also show that the maximum expectation value of the Clauser-
Horne-Shimony-Holt (CHSH) operator,
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$, which is another measure of
entanglement, exhibits decay due to the assumed interaction with the dephasing
field. Bell [4] and Clauser $et\;al.$ [5] demonstrate that Local-Hidden-
Variable theories can only produce correlations of the
$\langle{\hat{B}_{CHSH}\rangle}$ having a maximum value of 2. Cirel’son [6]
later shows that when accounting for quantum phenomena, a value as high as
$2\sqrt{2}$ can be obtained.

Liu $et\;al.$ [3] show that some entangled quantum states having values of
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}>2$ undergo an evolution in this
global dephasing (or decoherence) environment that results in
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}<2$ whereas others exhibit a decay
but maintain $\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}>2$ and, thus, always
produce non-local correlations. A drop in value of
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ below 2 implies that the
correlations resulting from measurements of a quantum system can no longer be
attributed with certainty to the entangled nature of the quantum system, as it
is possible for unentangled systems to exhibit values of
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ as high as 2. This is what is
referred to as the sudden death of non-locality.

The loss of non-locality is validated by Liu $et\;al.$ [3] experimentally in a
quantum optics experiment using pairs of entangled photons that demonstrate
that certain Bell diagonal states exhibit sudden death of non-locality (i.e.,
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ drops below 2), whereas other
Bell diagonal states maintain non-local correlations. In all cases, however,
the entanglement as measured by the concurrence remains constant.

The present work quantifies the effects of the SEA principle as well as the
model of open quantum systems represented by the Lindblad equation on the
dynamics of an entangled system and, consequently, on how the entanglement as
measured by concurrence and $\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$
changes with system state. In particular, it is demonstrated that under the
SEA principle, certain states may also exhibit decay in
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ but never drop below 2 whereas
others undergo sudden-death of non-locality depending on which of the two
types of perturbations are used to generate the initial state of the
evolution. As will also be seen, the type of perturbation used to generate the
initial state determines whether or not the concurrence under the Lindblad
formulation remains constant.

To accomplish this, Section II first presents the measures used to quantify
system entanglement. This is followed by a description and discussion of the
SEAQT and Lindblad equations and the Hamiltonian used to model the experiment.
The behavior of the perturbed Bell diagonal states relative to the two
equations of motion is also discussed, and because Bell diagonal states are
non-dissipative in the SEAQT framework, the perturbation methods used here to
perturb the non-dissipative baseline Bell diagonal state to a dissipative one
are presented. Specifically, two types of perturbation methods are used. The
resulting dynamics relative to the SEAQT and Lindblad equations of motion are
presented in Section III. Finally, Section IV provides the primary conclusions
of this work.

##  II Model Development

To model the evolution of entanglement due to the principle of steepest
entropy ascent, three entanglement measures and the SEAQT equation of motion
are first reviewed. The entanglement measures are reviewed in Section II.1 and
are used to examine how the entanglement changes throughout the system’s
evolution. Next, Section II.2 presents the SEAQT equation of motion for a
general quantum system and how it is used to predict the dissipation and
associated change in entanglement characteristics. The class of states of
interest in this work are known as Bell-diagonal states and turn out to be
non-dissipative relative to the SEAQT equation of motion (a brief proof of
which is also presented in Section III). Thus, to examine the dynamics of
dissipation associated with SEAQT, it is necessary to generate dissipative,
perturbed states from a baseline Bell diagonal state. Section II.3 then
presents the Lindblad formulation where Lindblad operators are proposed to
represent a pure dephasing environment independent of the state of the system
and interpreted as a bath of photons interacting with the system. The two
methods for creating the required perturbed states are presented in Section
II.5 and common features relating the perturbed state to that of the baseline
Bell diagonal state are highlighted. The evolutions of these perturbed states
using both the SEAQT and Lindblad equations of motion and the associated
changes in entanglement are then presented in Section III.

###  II.1 Entanglement Measures

To analyze the evolution of the entanglement within the SEAQT framework, two
different metrics are used: i) the concurrence, $E(\hat{\rho})$, and ii) the
maximum expectation value of the CHSH operator,
$\langle{\hat{B}_{CHSH}\rangle}$, were $\hat{\rho}$ is the density operator.
Many others exist (for a rather comprehensive review see [7]). The relative
entropy $D(\hat{\rho}||\hat{\rho}_{0})$, which is related to the difference
between two states, is also explored to measure the distance between the
baseline Bell diagonal state and the perturbed state.

The concurrence, $E(\hat{\rho})$, proposed by Hill and Wooters [8] and Wooters
[9] is intended to be a faithful representation of the entanglement and is
computed as

$E(\hat{\rho})=\max(0,\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4})$ (1)

where the $\lambda_{i}$ are the eigenvalues of the operator $\hat{R}$ in
decreasing order as $i$ increases and $\hat{R}$ is determined using the
relation

$\hat{R}=\sqrt{\sqrt{\hat{\rho}}\hat{\tilde{\rho}}\sqrt{\hat{\rho}}}$ (2)

Here the matrix $\hat{\tilde{\rho}}$ is computed as

$\hat{\tilde{\rho}}=\left(\hat{\sigma}_{Y}\otimes\hat{\sigma}_{Y}\right)\hat{\rho}^{*}\left(\hat{\sigma}_{Y}\otimes\hat{\sigma}_{Y}\right)$
(3)

where $\hat{\rho}^{*}$ is the complex conjugate of $\hat{\rho}$ and
$\hat{\sigma_{Y}}$ is the Pauli-y matrix.

The second common entanglement measure used here is the maximum expectation
value of the CHSH operator,
$\langle{\hat{B}_{CHSH}(\hat{\rho})\rangle}_{\text{max}}$, which is expressed
as

$\langle{\hat{B}_{CHSH}(\hat{\rho})\rangle}_{\text{max}}=2\sqrt{h_{i}+h_{j}}$
(4)

where $h_{i}$ and $h_{j}$ are the two largest eigenvalues of the matrix
$T_{\hat{\rho}}^{T}T_{\hat{\rho}}$ with $T_{\hat{\rho}}^{T}$ the transpose of
$T_{\hat{\rho}}$ and the elements of $T_{\hat{\rho}}$ given by
$t_{ij}=\text{Tr}(\hat{\rho}\left(\hat{\sigma}_{i}\otimes\hat{\sigma}_{j}\right))$.
Here again the $\hat{\sigma}_{i}$ and $\hat{\sigma}_{j}$ are the Pauli
matrices used to construct $T_{\hat{\rho}}$.

Now, for a composite system in a dissipative state consisting of subsystems
$A$ and $B$, the state of each subsystem in the SEAQT framework evolves to a
unique local equilibrium (canonical) state given by
$\hat{\rho}_{\text{eq},A(B)}=e^{-\beta_{A(B)}\hat{H}_{A(B)}}/Z_{A(B)}$ with
each canonical state in mutual stable equilibrium with the other. In this
expression, $\hat{H}$ is the Hamiltonian, $Z$ is the partition function and
$\beta$ is a Lagrange multiplier inversely proportional to the temperature of
the state. In the Lindblad formulation, the subsystem states evolve towards
stationary states that are characterized by a time-invariant density matrix
under the dynamics governed by the Lindblad equation. However, these
stationary states (i.e., steady states) are not equilibrium (or canonical)
states since under the SEAQT dynamics they are dissipative and continue to
evolve toward stable equilibrium.

Finally, to analyze the effects of a perturbation on the state of a system,
the quantum relative entropy, which measures the distinguishability between
two states, is used. It is given by

$D(\hat{\rho}||\hat{\rho}_{0})=\text{Tr}\hat{\rho}\ln\hat{\rho}-\text{Tr}\hat{\rho}\ln\hat{\rho}_{0}\,.$
(5)

When there is no perturbation, the relative entropy is zero. Otherwise, it is
positive, sometimes decreasing and then increasing or always increasing as the
state moves away from the Bell diagonal state employed as our baseline state.

###  II.2 SEAQT Equation of Motion for a General Quantum System

Because entanglement is inherently a feature of composite system states,
analyzing the change in entanglement requires analyzing the dynamics of the
entire composite system. The associated Hilbert space for these composite
systems is $\mathcal{H}=\bigotimes_{J=1}^{M}\mathcal{H}_{J}$ where $M$ is the
number of subsystems in the composite system. Thus, to understand the effects
that the dissipation associated with the SEA principle has on the time
evolution of the entanglement, the SEAQT equation of motion for a general
quantum system is required and is expressed as

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}\left[\hat{H},\hat{\rho}\right]-\frac{D\left(\hat{\rho}\right)}{Dt}$
(6)

where the dissipation operator $D\left(\hat{\rho}\right)/Dt$ for a composite
of $M$ subsystems is

$\frac{D\left(\hat{\rho}\right)}{Dt}=\sum_{J=1}^{M}\frac{1}{\tau_{J}}\hat{D}_{J}\otimes\hat{\rho}_{\bar{J}}$
(7)

Here the partial trace of the composite density operator onto the
$J^{\text{th}}$ subsystem is given by
$\hat{\rho}_{J}=\text{Tr}_{\bar{J}}(\hat{\rho})$ where $\bar{J}$ represents
the composite of all subsystems $\neq J$. Likewise, the density operator for
the composite of all subsystems except the $J^{\text{th}}$ can be computed as
$\hat{\rho}_{\bar{J}}=\text{Tr}_{J}(\hat{\rho})$. The dissipation operator for
the $J^{\text{th}}$ subsystem, $\hat{D}_{J}$, is then given by

$\hat{D}_{J}=\frac{1}{2}\left(\hat{\rho}_{J}\tilde{D}_{J}+\left(\hat{\rho}_{J}\tilde{D}_{J}\right)^{\dagger}\right)$
(8)

where $\tilde{D}_{J}$ represents a ratio of determinants for the
$J^{\text{th}}$ subsystem that ensures that the expectation values associated
with the generators of the motion, $\hat{R}_{iJ}$, of the $J^{\text{th}}$
subsystem remain constant. $\tilde{D}_{J}$ is then computed as

$\tilde{D}_{J}=\sqrt{\hat{\rho}_{J}}\frac{\left|\begin{matrix}\left(\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{N(J)J}\right)^{J}\\\
\left(\hat{R}_{1J},\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{R}_{1J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{1J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{1J},\hat{R}_{N(J)J}\right)^{J}\\\
\left(\hat{R}_{2J},\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{R}_{2J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{2J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{2J},\hat{R}_{N(J)J}\right)^{J}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\left(\hat{R}_{N(J)J},\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{R}_{N(J)J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{N(J)J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{N(J)J},\hat{R}_{N(J)J}\right)^{J}\end{matrix}\right|}{\Gamma^{J}}$
(9)

where the subscript $N(J)$ is the number of generators of the motion for
subsystem $J$ and the Gram determinant $\Gamma^{J}$ is computed as

$\Gamma^{J}=\left|\begin{matrix}\left(\hat{R}_{1J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{1J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{1J},\hat{R}_{N(J)J}\right)^{J}\\\
\left(\hat{R}_{2J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{2J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{2J},\hat{R}_{N(J)J}\right)^{J}\\\
\vdots&\vdots&\ddots&\vdots\\\
\left(\hat{R}_{N(J)J},\hat{R}_{1J}\right)^{J}&\left(\hat{R}_{N(J)J},\hat{R}_{2J}\right)^{J}&\cdots&\left(\hat{R}_{N(J)J},\hat{R}_{N(J)J}\right)^{J}\end{matrix}\right|$
(10)

For the remainder of this work, all subsystems $J$ will have the same two (and
only two) generators of the motion. The first, $\hat{R}_{1J}=\hat{I}$, ensures
that $\text{Tr}(\hat{\rho})=1$, and the second, $\hat{R}_{2J}=\hat{H}$,
ensures that the system evolution conserves the energy. Eq. (9), thus, reduces
to

$\tilde{D}_{J}=\sqrt{\hat{\rho}_{J}}\frac{\left|\begin{matrix}\left(\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{I}\right)^{J}&\left(\hat{H}\right)^{J}\\\
\left(\hat{I},\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{I},\hat{I}\right)^{J}&\left(\hat{I},\hat{H}\right)^{J}\\\
\left(\hat{H},\hat{B}\ln(\hat{\rho})\right)^{J}&\left(\hat{H},\hat{I}\right)^{J}&\left(\hat{H},\hat{H}\right)^{J}\end{matrix}\right|}{\left|\begin{matrix}\left(\hat{I},\hat{I}\right)^{J}&\left(\hat{I},\hat{H}\right)^{J}\\\
\left(\hat{H},\hat{I}\right)^{J}&\left(\hat{H},\hat{H}\right)^{J}\end{matrix}\right|}$
(11)

where each term $\left(\hat{F},\hat{G}\right)^{J}$ is defined as

$\left(\hat{F},\hat{G}\right)^{J}=\frac{1}{2}\text{Tr}_{J}\left(|\hat{\rho}_{J}|\\{(\hat{F})^{J},(\hat{G})^{J}\\}\right)$
(12)

and the so-called locally perceived operator $\left(\hat{F}\right)^{J}$ is
computed as

$\left(\hat{F}\right)^{J}=\text{Tr}_{\bar{J}}\left(\left(\hat{I}_{J}\otimes\hat{\rho}_{\bar{J}}\right)\hat{F}\right)$
(13)

and in a similar fashion $\left(\hat{G}\right)^{J}$. The above formulation has
been shown to effectively simulate pure dephasing, which in the SEAQT
framework is an effect only of the irreversibilities of the system, i.e., of
internal dissipation [10].

Now, for a composite system whose subsystems $A$ and $B$ are said to be weakly
interacting, the composite system Hamiltonian can be written as the sum of
Hamiltonians for each subsystem $A$ and $B$, i.e., $\hat{H}_{A}$ and
$\hat{H}_{B}$, and an interaction term $\hat{V}_{AB}$ such that

$\hat{H}=\hat{H}_{A}\otimes\hat{I}_{B}+\hat{I}_{A}\otimes\hat{H}_{B}+\hat{V}_{AB}$
(14)

where the $\hat{I}_{A(B)}$ are the identity operators in $\mathcal{H}_{A(B)}$
and the composite system identity operator can be computed as
$\hat{I}=\hat{I}_{A}\otimes\hat{I}_{B}$. While it is well known that two
systems must interact to become entangled, this work examines the evolution of
two systems that have already interacted to become entangled and are no longer
interacting, implying that $\hat{V}_{AB}=0$. With this model, it is possible
to attribute any change in the entanglement characteristics of the composite
system solely to internal dissipation, since the subsystems no longer interact
to change the entanglement characteristics through the reversible dynamics.
The dynamics of the SEAQT equation of motion are examined below for a specific
class of entangled states, known as Bell diagonal states, to better understand
how the entanglement changes. However, a discussion of the Lindblad equation
of motion is provided first in the next section.

###  II.3 Lindblad Equation of Motion

The theory of open quantum systems offers a way to determine the evolution of
a system subject to external perturbations by an environment. The equation of
motion for this framework is the Lindblad equation. To compare the loss of
non-locality as well as the loss of concurrence obtained by SEAQT formulation
with that predicted by the Lindblad equation, the latter is formulated based
on a Markovian dynamics that can incorporate the loss of locality due to pure
phase damping via the jump operators

$\hat{L}_{1}^{A(B)}=\left(\begin{matrix}0&0\\\
0&1\end{matrix}\right)\,,\quad\hat{L}_{2}^{A(B)}=\left(\begin{matrix}1&0\\\
0&0\end{matrix}\right)$ (15)

The use of these operators results in dephasing irrespective of whether or not
the system is in state $|0\rangle$ or state $|1\rangle$ of the Hamiltonian
basis or in some superposition of these states. In other words, each operator
contributes to the decay of the off-diagonal elements of the density matrix
without energy exchange. This choice of operators leads to the equation of
motion for the reduced-density matrix given by

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}\left[\hat{H},\hat{\rho}\right]+\gamma\mathcal{L}(\hat{\rho})$
(16)

where the coefficient $\gamma$ represents the strength of the interaction and
the operator $\mathcal{L}(\hat{\rho})$ depends on the jump operators as
follows:

$\mathcal{L}(\hat{\rho})=\sum_{m=1}^{2}\hat{L}_{m}\hat{\rho}\hat{L}_{m}^{\dagger}-\frac{1}{2}\\{\hat{L}_{m}^{\dagger}\hat{L}_{m},\hat{\rho}\\}\,.$
(17)

To incorporate the dynamics of the two subsystems A and B, the general density
matrix is given by $\hat{\rho}=\hat{\rho}_{A}\otimes\hat{\rho}_{B}$ and
initially assumed to be a Bell diagonal state that once perturbed evolves in
time without generating new correlations. The density operators for each
subsystem are then found from partial traces of $\hat{\rho}$ and the jump
operators for each subsystem constructed. The global jump operator is given by
the Kronecker product $\hat{L}_{i}^{A}\otimes\hat{L}_{i}^{B}$.

###  II.4 Equilibrium and Stationary States

In the SEAQT framework, any state whose entropy does not evolve in time (i.e.,
$dS/dt=0$) is a non-dissipative state. In Eq. (6), this corresponds to
$D\left(\hat{\rho}\right)/Dt=0$ provided there are no reversible dynamics,
i.e., $\hat{\rho}$ commutes with the Hamiltonian $\hat{H}$. This is the case
for every stable equilibrium state. When the only generators of the motion are
the identity and Hamiltonian operators, the density operator for such states
is defined by the canonical distribution
$\hat{\rho}_{eq}=\textrm{exp}(-\beta\hat{H})/Z$.

Other states in which there is no dissipation are so-called partially
canonical states of a simple quantum system [11], which are not stable
equilibrium states, i.e., they are not globally stable but instead only
locally so. This is the case for unstable or metastable equilibrium states
[12] or for maximally entangled states such as Bell diagonal states. Composite
quantum systems may experience partially canonical states where the reversible
dynamics of the system are still present but the dissipation is zero. As an
example, consider the spectral representation of the composite system density
operator given by

$\hat{\rho}=\sum_{i}\lambda_{i}\hat{P}_{i}$ (18)

where the $\lambda_{i}$ are the density operator eigenvalues and the
$\hat{P}_{i}$ are the associated projectors. Now, if the system is assumed to
be bipartite where both subsystems $A$ and $B$ have two energy eigenlevels and
the composite system state is a Bell diagonal state, the four projectors can
be written as

$\hat{P}_{\Phi^{\pm}}=|\Phi^{\pm}\rangle\langle{\Phi^{\pm}}|=\frac{1}{2}\begin{bmatrix}1&0&0&\pm
1\\\ 0&0&0&0\\\ 0&0&0&0\\\ \pm 1&0&0&1\\\ \end{bmatrix}$ (19)
$\hat{P}_{\Psi^{\pm}}=|\Psi^{\pm}\rangle\langle\Psi^{\pm}|=\frac{1}{2}\begin{bmatrix}0&0&0&0\\\
0&1&\pm 1&0\\\ 0&\pm 1&1&0\\\ 0&0&0&0\\\ \end{bmatrix}$ (20)

To obtain a representation of the state for each subsystem, the partial trace
of the composite can be found such that

$\hat{\rho}_{J}=\text{Tr}_{\bar{J}}\left(\hat{\rho}\right)=\sum_{i}\lambda_{i}\text{Tr}_{\bar{J}}\left(\hat{P}_{i}\right)=\sum_{i}\lambda_{i}\hat{I}_{J}/2=\hat{I}_{J}/2$
(21)

where $J=A,B$ and $\bar{J}=B,A$. Thus, for any Bell diagonal state, the
resulting density operator for subsystem $A$ or $B$ is maximally mixed. Noting
that any function of a density operator can be written as a function of its
eigenvalues, i.e.,

$f\left(\hat{\rho}\right)=\sum_{i}f\left(\lambda_{i}\right)\hat{P}_{i}$ (22)

it follows that for any Bell diagonal state, the partial trace of any function
of the density operator (e.g., $\hat{B}\ln(\hat{\rho})$ where $\hat{B}$ is the
idempotent operator defined in [13]) will be proportional to the identity
operator. Next, theorem 7 of [14] states that any density operator satisfying
the following relationship for all $J$ is a non-dissipative state of the
composite system for which the generators of the motion are $\hat{I}$ and
$\hat{H}$:

$\hat{\rho}_{J}\left(\hat{B}\ln(\hat{\rho})\right)^{J}=\hat{\rho}_{J}\left[\lambda_{IJ}\hat{I}_{J}+\lambda_{HJ}\hat{H}\right]$
(23)

where $\lambda_{IJ}$ and $\lambda_{HJ}$ must be real scalars. Evaluating
$\left(\hat{B}\ln(\hat{\rho})\right)^{J}$ for a Bell diagonal state results in

$\displaystyle\left(\hat{B}\ln(\hat{\rho})\right)^{J}=\text{Tr}_{\bar{J}}\left(\left(\hat{I}_{J}\otimes\hat{\rho}_{\bar{J}}\right)\hat{B}\ln(\hat{\rho})\right)$
$\displaystyle=\text{Tr}_{\bar{J}}\left(\hat{B}\ln(\hat{\rho})\right)=\alpha\hat{I}_{J}$
(24)

where $\alpha$ is a real scalar dependent upon the eigenvalues of the Bell
diagonal state. Thus, all Bell states satisfying Eq. (23) are non-dissipative,
yet they are not stable equilibrium states. Therefore, to analyze the change
of entanglement of states related to Bell diagonal states in particular using
the SEAQT equation of motion, it is necessary to use perturbation methods to
find dissipative states that have features similar to Bell diagonal states.

In the case of the Lindblad framework, the stationary or steady state to which
the equation of motion evolves, whether thermal or not, may be independent of
the initial state, which because of our choice of jump operators here is.
Thus, the stationary state is unique provided

$\mathcal{L}(\hat{\rho}_{SS})=0$ (25)

which implies that the semi-group relaxes as long as the zero eigenvalues of
$\mathcal{L}$ are non-degenerate and the remaining ones have a negative real
part. If not, the final stationary state may depend on the initial state [15].
In either case, the stationary state is dissipative and not canonical.

###  II.5 Perturbation Approaches

The baseline density operator $\hat{\rho}_{0}$ that is perturbed here is the
one for the Bell diagonal state in the c-configuration space given by

$\hat{\rho}_{0}=\frac{1}{4}\left(\hat{I}_{A}\otimes\hat{I}_{B}+\sum_{i=1}^{3}c_{i}\hat{\sigma}_{iA}\otimes\hat{\sigma}_{iB}\right)$
(26)

where the $\hat{\sigma}_{iA(B)}$ are the Pauli x, y, and z matrices,
respectively, $|c_{i}|\leq 1$ are real scalar coefficients. In the following,
the scalar coefficients defining the state are $c_{1}=0.996$, $c_{2}=0.4$, and
$c_{3}=-0.4$. This section presents two distinct approaches for generating
perturbed states from $\hat{\rho}_{0}$ where each perturbation approach
produces perturbed states that share some particular feature with
$\hat{\rho}_{0}$.

####  II.5.1 Pure State Weighted-Average Approach

The first proposed approach is a simple weighted average of $\hat{\rho}_{0}$
with the density operator, $\hat{\rho}_{\text{pure}}$, of a pure state of zero
energy whose subsystems both also have zero energy. In a sense, this
perturbation approach is compatible with a constrained perturbation with only
energy as the constraint. The following expression is used to generate a
perturbed state $\hat{\rho}_{0}^{\prime}$ with this method:

$\hat{\rho}_{0}^{\prime}=\zeta\hat{\rho}_{0}+\left(1-\zeta\right)\hat{\rho}_{\text{pure}}$
(27)

where

$\hat{\rho}_{\text{pure}}\equiv\hat{\rho}_{A,\text{pure},x}\otimes\hat{\rho}_{B,\text{pure},x}=\frac{1}{4}\begin{bmatrix}1&1\\\
1&1\end{bmatrix}\otimes\begin{bmatrix}1&1\\\ 1&1\end{bmatrix}$ (28)

Here the subscript $x$ denotes that the polarization vector of
$\hat{\rho}_{A(B),\text{pure},x}$ points along the x-axis of the Bloch sphere.
The weighting parameter $\zeta$ ranges from 0 to 1. This method is used
because the resulting perturbed state density operator,
$\hat{\rho}_{0}^{\prime}$, has the same energy as $\hat{\rho}_{0}$ (though
only one particular value of $\zeta$ ensures that $\hat{\rho}_{0}^{\prime}$
has the same entropy as $\hat{\rho}_{0}$) and also because no new entanglement
is introduced into the perturbed state. Thus, while $\hat{\rho}_{0}^{\prime}$
is not a Bell diagonal state, this guarantees that the stable equilibrium
state to which $\hat{\rho}_{0}^{\prime}$ evolves is a Bell diagonal state
although one different than the original Bell diagonal state.

####  II.5.2 General Bipartite System Approach

The second perturbation approach is the general bipartite perturbation scheme
referenced in [16, 17]. This perturbation method enables the perturbation of
an arbitrary density operator, while preserving its unit trace, energy,
entropy, and other physical constraints. The algorithm begins by calculating
the non-negative square root of the baseline density operator as

$\hat{\gamma}_{0}=\sqrt{\hat{\rho}_{0}},$ (29)

after which the perturbed square root is constructed by perturbing
$\hat{\gamma}_{0}$ according to

$\hat{\gamma}_{\epsilon}=\hat{\gamma}_{0}+\frac{1}{2}\sum_{i,j=0}^{3}\eta_{ij}\hat{\sigma}_{i}^{A}\otimes\hat{\sigma}_{j}^{B},$
(30)

where $\hat{\sigma}_{i}$ represents the set of Pauli matrices $\hat{I}$,
$\hat{\sigma}_{x}$, $\hat{\sigma}_{y}$ and $\hat{\sigma}_{z}$, the
superscripts $A$ and $B$ denote the subsystems to which the operators pertain,
and $\eta\in$ GUE with mean zero and standard deviation $\sigma$. The
perturbed density matrix is then constrained to maintain the required
constraints as

$\text{Tr}\left(\hat{\gamma}_{r}^{2}\hat{C}_{i}(\hat{\gamma}_{r})\right)=\text{Tr}\left(\hat{\rho}_{0}\hat{C}_{i}(\hat{\gamma}_{0})\right),$
(31)

where the $C_{i}$ are the
${\hat{I},\hat{H},-\hat{B}(\hat{\gamma}^{2}-\hat{I})}$ operators. Eq. (31)
thus, maintains constant unit trace, energy, and linear entropy respectively.
Furthermore,

$\hat{\gamma}_{r}=\hat{\gamma}_{\epsilon}-\sum_{i=1}^{3}\\{\hat{G}_{i}(\hat{\gamma}_{\epsilon}),\hat{\gamma}_{\epsilon}\\}\lambda_{i},$
(32)

where the $\lambda_{i}$ coefficients are derived from the constraint
equations, Eq. (31), and the operators $\hat{G}_{i}$ are part of the
symmetrized gradients $\\{\hat{G}_{i}(\hat{\gamma}^{2}),\hat{\gamma}\\}$ with
$\hat{G}_{1}=\hat{I}$, $\hat{G}_{2}=\hat{H}$ and
$\hat{G}_{3}(\hat{\gamma})=-(2\hat{\gamma}^{2}-I)$. The resulting set of
equations contains 12 roots according to the classical Bezout theorem, and
only the solution with $E\left(\hat{\rho}\right)$ close to the Bell diagonal
state is kept, i.e., that perturbed state with entanglement and separability
close to the original Bell diagonal state.

##  III Results

This section presents results for the evolution of system entanglement using
the SEAQT and Lindblad equations of motion starting from perturbed states. In
addition, various correlations between the properties of the perturbed and
stable equilibrium states are presented for the SEAQT framework. First,
Section III.1 shows entanglement evolutions using the weighted average
perturbation approach to generate the initial states used by the SEAQT
equation of motion as well as the Lindblad equation of motion. Next, Section
III.2 does likewise except that the initial states used are generated by the
general bipartite perturbation approach with unit trace and constant energy
and entropy constraints and a perturbation matrix with a standard deviation of
$\sigma=10^{-1}$. This arbitrary value is chosen to align with experimental
values, but alternative selections could be employed.

a)
b)

Figure 1: Theoretical and experimental a) concurrence and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ evolutions for the Bell diagonal
state characterized by the parameters $c_{1}=1$, $c_{2}=0.4$, and
$c_{3}=-0.4$; the experimental (black dots) and Markovian theoretical (red
curve) results of [3] were generated using a quantum optics experiment and the
Kraus-Operator-Sum approach; the initial state for the SEAQT and Lindblad
evolutions is generated using the first perturbation method with $\zeta=0.68$.

To gain a qualitative understanding of the effects of dissipation due to an
external environment on the entanglement properties of the system, theoretical
and experimental results developed by Liu $et\;al.$ [3] for the evolution of a
Bell diagonal state in a stochastic global dephasing (or decoherence) field
are shown in Fig. 1. Their theoretical results generated using the Kraus-
Operator-Sum approach predict that the concurrence remains constant, while
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ decays. Thus, after a finite
time (dimensionless time of $t\approx 0.58$) the theoretical results suggest
that the system undergoes so-called sudden death of non-locality, i.e.,
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ goes below 2. The experimental
results produced by Liu $et\;al.$ qualitatively validate the theoretical
results since the measured concurrence remains roughly constant, while the
measured $\langle\hat{B}_{CHSH}\rangle_{\text{max}}$ decays through time
although sudden death of non-locality occurs at a dimensionless time of
$t\approx 0.2$.

Fig. 1 also presents results for the evolution of $E(\hat{\rho})$ and $\langle
B_{CHSH}\rangle_{\text{max}}$ using the SEAQT and Lindblad equations of motion
and a weighting parameter of $\zeta=0.68$ to generate the initial state using
the first perturbation method. As can be seen, the sudden death of non-
locality predicted with the SEAQT as well as the Lindblad equation matches the
experimental result better than the model proposed by Liu $et\;al.$,
predicting this death of non-locality at times near 0.2. However, using the
weighted-average perturbation method, the $E(\hat{\rho})$ is significantly
below the 0.4 reported by the experimental results. The discrepancy is due to
the weighting parameter which moves the perturbed state away from the Bell
diagonal state by nearly 30% in the direction of a pure state. In Section
III.1 the evolution of this entanglement measure and the other and their
relation to the entropy as well as the entropy generation are explored more
systematically, indicating a correlation between thermodynamics and both
entanglement measures.

###  III.1 Parametric Study of Entanglement Evolution Using Initial States
Generated with the Weighted-Average Perturbation Method

####  III.1.1 SEAQT framework

This section presents evolutions for the i) system concurrence,
$E(\hat{\rho})$, and the maximum expectation value of the CHSH operator,
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$; ii) correlations between the
relative entropy of the perturbed initial state and final stable equilibrium
state and the weighting parameter $\zeta$ and the coefficient $c_{1}$; iii)
correlations between the perturbed state and the final state entropy as a
function of the coefficient $c_{1}$; iv) correlations between the equilibrium
state entropy and the weighting parameter $\zeta$ and the coefficient $c_{1}$;
and v) correlations between the total entropy generation and the weighting
parameter $\zeta$ and the coefficient $c_{1}$.

Fig. 2 shows the evolutions of $E(\hat{\rho})$ and
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ for nine different values of the
weighting parameter $\zeta$. For both properties, initial decay in their value
occurs as entropy is produced before becoming constant as the system reaches a
stable equilibrium state. In addition to the initial values of these
properties dropping as $\zeta$ decreases, the overall decay in the value of
each property increases as $\zeta$ decreases. Fig. 2 also shows that for
values of $\zeta\approx 0.76$ and below, the composite system undergoes sudden
death of non-locality in the way that the Bell diagonal state predicted by Liu
$et\;al.$ does. For values above $\zeta\approx 0.76$, the correlations due to
quantum entanglement remain as the system relaxes to a stable equilibrium
state.

a)
b)

Figure 2: a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ evolutions in the SEAQT
framework for values of the weighting parameter $\zeta=1$, $0.96$, $0.92$,
$0.88$, $0.84$, $0.8$, $0.76$, $0.72$, and $0.68$ where $\zeta=1$ corresponds
to the uppermost dashed line and each lower solid line corresponds to the next
lower value of $\zeta$; the experimental results (black dots) of [3] were
generated using a quantum optics experiment.

a)
b)

Figure 3: Plots of a) relative entropy between the perturbed and the Bell
diagonal and b) the perturbed and equilibrium state entropy for
$c_{1}=0.4,0.8,0.9,0.996$ in the SEAQT framework. The red line indicates the
states where the final state equals the perturbed initial entropy whereas the
$\blacksquare$ marker indicates the pure state entropy and the $\times$ marker
indicates the Bell diagonal state with $c_{1}=1$.

a)
b)

Figure 4: Correlations in the SEAQT framework between the weighting parameter
$\zeta$ and a) the stable equilibrium state entropy and b) the total entropy
generation for $c_{1}=0.1,0.4,0.8,0.9$ and $0.996$.

Fig. 3 a) shows the change of the relative entropy between the initial
perturbed state and the final stable equilibrium state predicted by the SEAQT
framework versus the weighting parameter $\zeta$, which varies between zero
and one in increments of $0.01$ and for different values of the coefficient
$c_{1}$. The Bell diagonal state corresponds to a $\zeta=1$ whereas a
$\zeta=0$ corresponds to the pure state. For the perturbed states, the
relative entropy decreases monotonically initially with increasing $\zeta$
after which it plateaus or increases. In addition, different baseline Bell
diagonal states defined here by changes in the scalar coefficient $c_{1}$
result in smaller decreases in the relative entropy as the value of $c_{1}$
decreases.

Fig. 3 b) provides the final state entropy of each evolution versus the
initial perturbed state entropy for different choices of the $c_{1}$
parameter. The right part of the curves shown in Fig. 3 b) represents states
with $\zeta\approx 1$ whereas the left part corresponds to states with
$\zeta\approx 0$. In addition, the red line shown represents the states whose
initial perturbed entropy is equal to the final state entropy. Thus, states
close to the Bell diagonal states remain close throughout the evolution of the
state in the SEAQT framework. This is an indication of the Lyapunov stability
sought by the SEAQT equation of motion [18] where the equilibrium solution of
the final state found by this equation is indeed the thermodynamically stable
state which is non-dissipative. Note that SEAQT also admits metastable
solutions that are related to the elements of the kernel of the dissipative
operator $D(\hat{\rho})/Dt$. In general, these solutions take the form

$\hat{\rho}=\frac{\hat{B}\exp\left(-\beta
H\right)\hat{B}}{\text{Tr}\left(\hat{B}\exp\left(-\beta H\right)\right)}$ (33)

where $\hat{B}=I-P_{\text{ker}\hat{\rho}}$ is the idempotent operator or
projector onto the range of $\hat{\rho}$, while $I-\hat{B}$ is the projector
onto the kernel of $\hat{\rho}$. Because of this, the von Neumann entropy can
be expressed as

$\displaystyle S$ $\displaystyle=$
$\displaystyle-\text{Tr}\left(\hat{\rho}\hat{B}\ln\hat{\rho}\right)$ (34)
$\displaystyle=$
$\displaystyle-\text{Tr}\left(\hat{\rho}\ln(\hat{\rho}+P_{\text{ker}\rho})\right)\,.$
(35)

and it can readily be shown that $S$ is well defined for any $\hat{\rho}$,
even one that is singular [19]. This, in particular, is the case for the pure
state entropy reported in Fig. 3 b). Thus, for pure states, the null vector,
$|\psi\rangle$, can easily be calculated, i.e.,
$\hat{\rho}_{0}|\psi\rangle=0$, which leads to the projector
$P_{\text{ker}\hat{\rho}}=|\psi\rangle\langle\psi|$ and a value for the von
Neumann entropy of the pure state of zero, i.e.,

$\displaystyle S_{\text{Pure}}$
$\displaystyle=-\text{Tr}(\hat{B}\hat{\rho}_{\text{pure}}\ln(\hat{\rho}_{\text{pure}}))=0$
(36)

In other words, since the pure state corresponds to a completely known state,
the projector projects out the zero eigenvalues leaving only the eigenvalue
$\lambda=1$.

The other limiting case reported in Fig. 3 b) is that of the Bell diagonal
state for which the null space depends on the selection of the coefficients
$c_{i}$, and for which all the null vectors are sums of the projectors

$\hat{P}_{\Phi^{\pm}}=|\Phi^{\pm}\rangle\langle\Phi^{\pm}|\,,\hat{P}_{\Psi^{\pm}}=|\Psi^{\pm}\rangle\langle\Psi^{\pm}|$
(37)

Furthermore, the null space can be classified according to its cardinality,
which accounts for the number of null vectors for a given selection of the
$c_{i}$. In general, the cardinality of the kernel of $\hat{\rho}_{0}$ can be
written as the Binomial coefficient, i.e.,

$\text{card}\,\left(\text{ker}\hat{\rho}\right)=\left(\begin{matrix}\text{dim}\mathcal{H}\\\
\text{dim}\left(\text{ker}(\hat{\rho})\right)\end{matrix}\right)$ (38)

Thus, the limiting case of $c_{1}=1$ and $c_{3}=-c_{2}$ seen in Fig. 3 b) is
an element of a set with cardinality 6 whose projector is
$\hat{B}=I-|\Phi^{-}\rangle\langle\Phi^{-}|-|\Psi^{-}\rangle\langle\Psi^{-}|$
and for which the von Neumann entropy is found to be

$\displaystyle S_{\text{eq},0}$
$\displaystyle=-\text{Tr}(\hat{\rho}_{0}\hat{B}\ln\rho_{0})$ (39)
$\displaystyle=\ln
2-\frac{1}{2}c_{2}\text{arctanh}\,c_{2}-\ln(1-c_{2}^{2})^{1/2}\,,$ (40)
$\displaystyle=0.6108\,.$ (41)

Now, Fig. 4 a) shows the equilibrium state entropy versus the weighting
parameter for different values of $c_{1}$. As can be seen, the larger stable
equilibrium entropy is associated with Bell diagonal states with $c_{1}=0.1$,
and as $c_{1}$ approaches 1, the stable equilibrium state entropy decreases
towards the limiting stationary (steady-state) entropy of
$S_{\text{eq},0}=0.6108$.

For different values of $c_{1}$, Fig. 4 b) shows correlations between the
weighting parameter $\zeta$ and the total entropy generation, which in this
case is equivalent to the change of entropy between the initial perturbed
state and the final stable equilibrium state entropy, i.e.,

$S_{\text{Gen}}=S(t=t_{\text{eq}})-S(t=0)\,,$ (42)

As can be seen, the entropy generation diminishes as the magnitude of the
weighting parameter, $\zeta$, increases, which means that the closer the
initial perturbed state is to the Bell diagonal state, the closer the final
stable equilibrium state is to the initial state. Again, this behavior is
related to the Lyupanov stability of the equation of motion, since the Bell
diagonal states are in the kernel of the dissipative operator, and, thus, the
states closer to the Bell diagonal state remain closer throughout the dynamic
evolution of the equation of motion. In contrast, initial perturbed states
closer to the pure state result in a higher amount of entropy generation.

####  III.1.2 Lindblad framework

In this section, the dynamics of the Linblad framework are considered. How
these dynamics affect the information measures $C(\hat{\rho})$ and $\langle
B_{CHSH}\rangle_{\text{max}}$ is examined first. Fig. 5 shows nine
perturbations employing the weighting parameter. As is seen, in contrast with
the SEAQT formulation, the Lindblad dynamics predicts a constant value for the
concurrence but like the SEAQT framework a decay of the non-locality measure
where the sudden death of non-locality for the $\zeta=0.68$ curve occurs at
$t\approx 0.23$ (as compared to $t\approx 0.19$ for the SEAQT evolution for
this value of the weighting parameter).

a)
b)

Figure 5: Theoretical and experimental a) concurrence and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ evolutions in the Lindblad
framework for values of the weighting parameter $\zeta=1$, $0.96$, $0.92$,
$0.88$, $0.84$, $0.8$, $0.76$, $0.72$, and $0.68$ where $\zeta=1$ corresponds
to the uppermost dashed line and each lower solid line corresponds to the next
lower value of $\zeta$; the experimental results (black dots) of [3] were
generated using a quantum optics experiment.

In general, the entropy changes in the Lindblad framework depend on the
Lindblad operator, $\mathcal{L}$, such that

$\frac{dS}{dt}=\gamma\text{Tr}\,\left(\mathcal{L}\left(\hat{\rho}\right)\ln\hat{\rho}\right)$
(43)

where the Lindblad operator for the case of Bell diagonal states can be
written as

$\mathcal{L}\left(\hat{\rho}_{0}\right)=-\frac{1}{4}\gamma^{2}\left(\begin{array}[]{cccc}0&0&0&c_{1}-c_{2}\\\
0&0&0&0\\\ 0&0&0&0\\\ c_{1}-c_{2}&0&0&0\\\ \end{array}\right)$ (44)

A state is not a stationary state of the Lindblad equation except when
$c_{1}=c_{2}$. In such cases, the entropy remains constant, and the density
matrix’s evolution is influenced only by the symplectic term.

Fig. 6 a) illustrates the relative entropy as a function of the weighting
parameter $\zeta$ and various values of $c_{1}$. A decrease in $c_{1}$ leads
to a greater deviation of the stationary (steady-state) entropy from the
initial perturbed state entropy. In contrast, when $c_{1}$ is approximately 1,
the difference between the perturbed Bell diagonal and stationary states
reduces significantly. In Fig. 6 b), the change in final entropy relative to
the initial perturbed state entropy is shown. The left most point of the
curves corresponds to $\zeta=0$, and the right most point to $\zeta=1$. As is
seen, for $c_{1}=0.4$ and $\zeta=1$, the final entropy equals the initial
perturbed entropy since $c_{1}=c_{2}$ in Eq. (44). However, for larger values
of $c_{1}$, this is not the case, indicating that, in general, Bell diagonal
states are not stationary states within the Lindblad framework.

a)
b)

Figure 6: Plots of a) Relative entropy between the perturbed Bell diagonal and
stationary states and b) the perturbed Bell diagonal and stationary state
entropies for $c_{1}=0.4,0.8,0.9,0.996$ in the Lindblad framework. The red
line indicates the states where the final state entropy equals the perturbed
initial entropy whereas the $\blacksquare$ marker indicates the pure state
entropy and the $\times$ marker indicates the Bell diagonal state entropy with
$c_{1}=1$.

In Fig. 7, both the stationary state entropy and the entropy generation as
functions of the weighting parameter $\zeta$ are depicted. Notably, the
highest stationary state entropies correspond to the perturbed states with
$c_{1}=0.4$, while the entropy generation is the least for this caae.
Additionally, for $\zeta=1$ and $c_{1}=0.4$, the entropy generation is zero,
which, according to Eq. (44), implies the absence of entropy changes. This
observation suggests a unique stationary (steady state) under these specific
parameter conditions.

a)
b)

Figure 7: Correlationsin the Lindblad framework between the weighting
parameter $\zeta$ and a) the stationary state entropy and b) the total entropy
generation for $c_{1}=0.1,0.4,0.8,0.9$ and $0.996$.

###  III.2 Entanglement Evolution using Initial States generated with the
General Bipartite Constant Energy and Entropy Constrained Perturbations

####  III.2.1 SEAQT framework

This section presents results using the general bipartite system perturbation
approach with constant energy and constant linear entropy constraints and a
variance of $\sigma=10^{-1}$. Specifically, results are presented for the
evolutions of $E(\hat{\rho})$ and
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ as are results for the
correlations between the perturbed and equilibrium state $E(\hat{\rho})$ and
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ and the total entropy generation
and the total entropy change.

Fig. 8 shows the evolution within the SEAQT framework of $E(\hat{\rho})$ and
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ for five random perturbations.
For these perturbations, both entanglement measures exhibit little decay
before the system reaches an equilibrium state. Thus, the perturbation
approach with constant entropy and energy shows that the SEAQT evolution of
the density matrix for perturbations 1 (GP 1), 3 (GP 3), and 5 (GP5) match
both the experimental values for $E(\rho)$ and those for the $\langle
B_{CHSH}\rangle$ measure. This was not the case for any of the evolutions
associated with the weighted-average perturbations, i.e., either the
experimental concurrence was matched or the $\langle B_{CHSH}\rangle$ measure
but not both. In general, these evolutions of the concurrence as well as of
the $\langle B_{CHSH}\rangle$ measure will decay further since stable
equilibrium has not been reached for the final values shown in the figure.

a)
b)

Figure 8: The a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ SEAQT evolutions for the five
initial states generated with the general perturbation method.

Fig. 8 also shows that the evolutions for GP1 to GP3 and for GP5 maintain
values of $\langle B_{CHSH}\rangle_{\text{max}}$ below the non-locality
condition at all times and all but the GP2 evolution match the measured
concurrence well. On the other hand, the beginning values of the GP4 evolution
show the presence of quantum non-locality with the loss of non-locality only
occurring at about $t=0.24$.

Fig. 9 a) shows histograms of the concurrence for the initial and final states
obtained from 1000 random perturbations with $\eta\in\text{GUE}$ and mean and
standard deviations of 0 and $10^{-1}$, respectively. It is observed that the
initial concurrence has a mean value of 0.38, closely matching the concurrence
of the pure Bell diagonal state. This is expected, as the constrained
perturbation assures that the initial perturbed state has a concurrence
similar to that of the original Bell diagonal state. Moreover, at the end of
the dynamics, the average concurrence has decreased somewhat, indicating a
reduction in entanglement. In contrast, Fig. 9 b) shows that the $\langle
B_{\text{CHSH}}\rangle_{\text{max}}$ exhibits higher dispersion for both the
initial and final states. Furthermore, considering the limit of non-locality
in states where $\langle B_{\text{CHSH}}\rangle_{\text{max}}\leq 2$, it is
found that the random perturbations result in nearly 32 % of states exhibiting
non-locality properties at the end of the dynamics.

a) b)

Figure 9: Initial perturbed state and final stable equilibrium state
histograms for a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ in the SEAQT framework. The gray
bars correspond to the initial state values whereas the orange bars to the
final state values.

Fig. 10 shows the correlation between the total entropy generation (or
production) and the total change of $E(\hat{\rho})$ and
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ from the initial perturbed state
to the final stable equilibrium state. Both properties show clear trends of
decrease with increasing entropy generation. One interesting aspect to note is
that while the average $E(\hat{\rho})$ is seen to only decrease relative to
the SEAQT dynamics, for a few of the initial perturbed states the evolution
resulting from this dynamics does cause an increase in
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$. Since the SEAQT equation of
motion and the resulting dissipation do not increase the entanglement for two
non-interacting subsystems, this increase indicates that these few states are
sufficiently close to the Bell diagonal state so that there is no decay in
entanglement or non-locality. Thus, for non-monotonic measures, like
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$, a less entangled state can have
a higher value of $\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ than a more
entangled state (e.g., a completely unentangled state can have a value as high
as 2 whereas entangled states can have values less than or greater than 2).

Additionally, Fig. 10 illustrates that an increase in the information measure
corresponds to an increase in entropy generation. This correlation can be
measured by considering the entropy generation as well as the decay of
$E(\hat{\rho})$ or $\langle B_{\text{CHSH}}\rangle_{\text{max}}$ as random
variables with Pearson coefficient

$r_{E(B)}=\frac{\text{Cov}(\delta
E(\hat{\rho}),S_{\text{Gen}})}{\sqrt{\text{Var}(\delta
E(\hat{\rho}))\cdot\text{Var}(S_{\text{Gen}})}}$ (45)

where $\delta E(\hat{\rho})=E(\hat{\rho})_{\text{eq}}-E(\hat{\rho})_{0}$ and
similarly for $\langle B_{\text{CHSH}}\rangle_{\text{max}}$. In this case, the
Pearson coefficients are $r_{E}\approx-0.94$ and $r_{B}\approx-0.91$
indicating a strong negative correlation between the information measures and
the entropy generation, i.e., that as the entropy generation increases the
change in $E(\hat{\rho})$ as well as $\langle
B_{\text{CHSH}}\rangle_{\text{max}}$ reduces.

a)
b)

Figure 10: Correlations of changes in a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ versus the total entropy
generation in the SEAQT framework.

####  III.2.2 Lindblad framework

Results for the evolutions of $E(\hat{\rho})$ and
$\langle\hat{B}_{CHSH}\rangle_{\text{max}}$ based on the Lindblad framework
are given in this section as are correlations between $E(\hat{\rho})$ and
$\langle\hat{B}_{CHSH}\rangle_{\text{max}}$ and the initial perturbed states
generated with the general perturbation method and the stationary states found
using the Lindblad framework. As observed in Fig. 11 the evolution of
concurrence for the five random perturbations is almost constant, whereas
$\langle\hat{B}_{CHSH}\rangle_{\text{max}}$ shows a decay in time which is
within the error bars of the experimental results. In general, all the
perturbations, predict a sudden death of non-locality for values of the
dimensionless time between about 0.05 and about 0.3 whereas $C(\hat{\rho})$
remains close to the 0.4 value reported in the experiments.

a)
b)

Figure 11: a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ evolutions based on the Lindblad
framework for the five random initial states generated using the general
perturbation method.

Fig. 12 presents histograms based on 1000 constrained state perturbations
generated using the general perturbation method and the corresponding
stationary states found using the Lindblad dynamics. In Fig. 12 a), the mean
concurrence value for the initial perturbed states (gray bars) is
approximately 0.4. This is a consequence of selecting specific roots from the
set of constraints for initializing the density matrix close to the
concurrence of the reference Bell diagonal state. Under the Lindblad dynamics,
the average concurrence tends to decrease, exhibiting greater variance than
that predicted by the SEAQT dynamics. Conversely, as shown in Fig. 12 b), the
initial and final state values of the $\langle B_{CHSH}\rangle_{\text{max}}$
measure display a less broad dispersion compared to the concurrence
histograms. Remarkably, at the end of the dynamics, $\langle
B_{CHSH}\rangle_{\text{max}}$ is less than 2.0 in 75 % of the cases,
indicating that most of the perturbed cases have a loss of non-locality.

a) b)

Figure 12: Initial perturbed state and final stationary state histograms for
a) $E(\hat{\rho})$ and b) $\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ in the
Lindblad framework. The gray bars correspond to the initial state values and
the orange bars to the final state values.

Finally, the scatter plots in Figure 13 depict the correlations of the changes
in the concurrence and the $\langle B_{CHSH}\rangle_{\text{max}}$ measure
relative to the entropy generation with Pearson coefficients of
$r_{E}\approx-0.83$ and $r_{B}\approx 0.91$. The former indicates a weaker
correlation between the change in the concurrence and the entropy generation
than that exhibited with the SEAQT dynamics, while the latter indicates that
the correlation between the change of the $\langle
B_{CHSH}\rangle_{\text{max}}$ measure and the entropy generation is almost the
same as that in the SEAQT dynamics. Notably, as the entropy generation
increases, the spread in the concurrence and the $\langle
B_{CHSH}\rangle_{\text{max}}$ measure values increases. This is true for the
SEAQT dynamics as well although the increase in the spread is not as
pronounced.

a)
b)

Figure 13: Correlations of change in a) $E(\hat{\rho})$ and b)
$\langle{\hat{B}_{CHSH}\rangle}_{\text{max}}$ versus the total entropy
generation in the Lindblad framework.

###  III.3 Comparison between the SEAQT and Lindblad frameworks

Fig. 14 illustrates the dynamics of the entropy, time rate of entropy change,
$E(\hat{\rho})$, and $\langle B_{CHSH}\rangle_{\text{max}}$ under both the
SEAQT and Lindblad frameworks. As seen in Fig. 14 a), the Lindblad approach
consistently results in a higher entropy than that predicted by the SEAQT
equation of motion at each instant of time. In contrast, the time rate of
change of the entropy is lower in the Lindblad framework, indicating a less
rapid approach to the final stationary state in this framework than to the
final stable equilibrium state in the SEAQT framework. As to the concurrence,
Fig. 14 b) shows that both models exhibit a relatively stable value over time
with only slight fluctuations. Intriguingly, in the context of the $\langle
B_{CHSH}\rangle_{\text{max}}$ measure, Fig. 14 c) shows that both frameworks
concur on an abrupt disappearance of non-locality, occurring at a
dimensionless time less than about 0.3, a finding borne out by the
experimental data. This simultaneous loss of non-locality in both frameworks
highlights the robustness of this phenomenon across different theoretical
approaches and reinforces its empirical validity.

a)
b)
c)

Figure 14: Evolution of a) the entropy and the time rate of change of the
entropy, b) the $E(\rho)$ and c) the $\langle B_{CHSH}\rangle_{\text{max}}$
measure in the SEAQT and Lindblad frameworks. The experimental data is include
in figures b) and c).

##  IV Conclusions

This paper has presented a study of the evolution of a bipartite system in an
entangled state using the dynamics of the SEAQT and Lindblad frameworks and
provided insight into the entanglement dynamics within both frameworks for
non-interacting subsystems. Results for these evolutions are indirectly
compared to evolutions predicted by the Kraus-Operator-Sum approach as well as
to data obtained from a quantum optics experiment. While the dynamics
predicted by the SEAQT framework differ from those predicted by the Kraus-
Operator-Sum approach, both predict the possibility of almost constant
concurrence and the presence of quantum entanglement or non-locality (via the
CHSH operator maximum expectation value being above 2). In addition, SEAQT and
Lindblad frameworks predict the sudden death of non-locality via a drop below
2 of the maximum CHSH operator expectation value at dimensionless times
consistent with the experimental data but much smaller than the times
predicted by the Kraus-Operator-Sum approach.

Using the weighted-perturbation approach to generate the initial states, a
good match between the experimental bounds and the predictions made with the
SEAQT and Lindblad frameworks for the behavior of the $\langle
B_{CHSH}\rangle_{\text{max}}$ measure is found. At the same time, a
significant discrepancy with respect to $E(\hat{\rho})$ occurs. In contrast,
using the general perturbation approach to generate the initiate states,
predictions by the dynamics of both the SEAQT and Lindblad frameworks of the
behavior of $\langle B_{CHSH}\rangle_{\text{max}}$ and $E(\hat{\rho})$ match
quite well with the experimental results. These results furthermore suggest
that the experimental sets are not able to create a diagonal state, but a
state close to it with the same energy and entropy.

Finally, both information measures show a strong negative correlation,
represented by Pearson coefficients greater than 0.8, with the entropy
generation, suggesting a deep connection between thermodynamics and quantum
information. Remarkably, the SEAQT shows a higher correlation coefficient for
$E(\hat{\rho})$ than the Lindblad formulation. This could be related to the
fact that SEAQT framework admits of a fluctuation-dissipation formulation,
which relates the entropy generation to the covariance (fluctuation) of the
generalized Massieu function. The intriguing connection between the decay of
these information measures in the SEAQT framework is left for future work.

##  V Acknowledgment

C. D. and A. S.-R., report financial support was provided by National Council
on Science and Technology. C. D. was supported by Universidad de Guanajuato
through the project CIIC 039/2023. C. D. and A. S.-R. thanks to Virginia Tech
for its kind hospitality in the development of this project. Finally, R. T.
Holladay would like to thank the U.S. Department of Defense for its Science,
Mathematics, and Research for Transformation (SMART) scholarship.

## References

  * Bennett  _et al._ [1993] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, Teleporting an unknown quantum state via dual classical and einstein-podolsky-rosen channels, Phys. Rev. Lett. 70, 1895 (1993).
  * Yu and Eberly [2003] T. Yu and J. H. Eberly, Qubit disentanglement and decoherence via dephasing, Phys. Rev. B 68, 165322 (2003).
  * Liu  _et al._ [2016] B.-H. Liu, X.-M. Hu, J.-S. Chen, C. Zhang, Y.-F. Huang, C.-F. Li, G.-C. Guo, G. b. u. Karpat, F. F. Fanchini, J. Piilo, and S. Maniscalco, Time-invariant entanglement and sudden death of nonlocality, Phys. Rev. A 94, 062107 (2016).
  * Bell [1964] J. S. Bell, On the einstein podolsky rosen paradox, Physics Physique Fizika 1, 195 (1964).
  * Clauser  _et al._ [1969] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Proposed experiment to test local hidden-variable theories, Phys. Rev. Lett. 23, 880 (1969).
  * Cirel’son [1980] B. S. Cirel’son, Quantum generalizations of bell’s inequality, Letters in Mathematical Physics 4, 93 (1980).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Hill and Wootters [1997] S. Hill and W. K. Wootters, Entanglement of a pair of quantum bits, Phys. Rev. Lett. 78, 5022 (1997).
  * Wootters [1998] W. K. Wootters, Entanglement of formation of an arbitrary state of two qubits, Phys. Rev. Lett. 80, 2245 (1998).
  * Montañez Barrera  _et al._ [2020] J. A. Montañez Barrera, C. E. Damian-Ascencio, M. R. von Spakovsky, and S. Cano-Andrade, Loss-of-entanglement prediction of a controlled-phase gate in the framework of steepest-entropy-ascent quantum thermodynamics, Phys. Rev. A 101, 052336 (2020).
  * Beretta [2006] G. P. Beretta, Nonlinear model dynamics for closed-system, constrained, maximal-entropy-generation relaxation by energy redistribution, Phys. Rev. E 73, 026113 (2006).
  * Gyftopoulos and Beretta [2005] E. P. Gyftopoulos and G. P. Beretta,  _Thermodynamics: Foundations and Applications: 2nd Edition_ (Dover Publication, 2005).
  * Beretta [2009a] G. P. Beretta, Nonlinear quantum evolution equations to model irreversible adiabatic relaxation with maximal entropy production and other nonunitary processes, Reports on Mathematical Physics 64, 139 (2009a).
  * Beretta  _et al._ [1985] G. P. Beretta, E. P. Gyftopoulos, and J. L. Park, Quantum thermodynamics. a new equation of motion for a general quantum system, Il Nuovo Cimento B (1971-1996) 87, 77 (1985).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Dynamics of non-markovian open quantum systems, Rev. Mod. Phys. 89, 015001 (2017).
  * Montanez-Barrera  _et al._ [2022] J. Montanez-Barrera, R. Holladay, G. Beretta, and M. R. von Spakovsky, Method for generating randomly perturbed density operators subject to different sets of constraints, Quantum Information Processing 21, 314 (2022).
  * Holladay [2019] R. T. Holladay,  _Steepest-Entropy-Ascent Quantum Thermodynamic Modeling of Quantum Information and Quantum Computing Systems_ , Ph.D. thesis, M.E. Dept., Virginia Tech, Blacksburg, VA (2019).
  * Beretta [2009b] G. P. Beretta, Nonlinear quantum evolution equations to model irreversible adiabatic relaxation with maximal entropy production and other nonunitary processes, Reports on Mathematical Physics 64, 139 (2009b).
  * Beretta [2010] G. P. Beretta, Maximum entropy production rate in quantum thermodynamics, Journal of Physics: Conference Series 237, 012004 (2010).