Hi,
I’m writing a density matrix as a supervector, e.g. density matrix of a N-site system is described by a 2N-site array where half of the array is the bra and the other is ket. I want to compute the trace of this density matrix by contracting the physical indices of one-half of the supervector with the other half. I would appreciate if someone can show me how to implement such action in C++, I’m new to Itensor. Thank you.
1 Like
Hi, thanks for the interesting question. Here are some thoughts and helpful suggestions:
If it’s possible for you to use the Julia version of ITensor instead, it has a lot more tools at this point to do the kind of thing you are asking about. In particular there are nice tools to work with MPOs, which is equivalent to the kind of representation you are talking about.
To clarify your question, have you already made this supervector in code and you are just asking about how to perform the trace? If so, how are you representing this supervector? As an MPS?
Or is your question about doing all of the above?
Thanks,
Miles
My strategy is to represent the supervector as an MPS, the Liouvillian as a non-Hermitian MPO and apply TDVP. I already managed to complete all those steps but am still thinking how to contract the MPS with itself, which is equivalent to taking the trace.
I also appreciate if there are complete tools for solving the Lindblad equations for dissipative interacting 1D systems.
Best,
Duy
Hi Duy,
Regarding tools to solve the Lindblad equations, would you mind posting the equation you are trying to solve? I haven’t worked with open systems much myself so I often forget the formulation. If I recall correctly, the steady state can be found by evolving by the L superoperator so maybe that’s something that can be done easily enough with existing ITensor tools.
If you are representing the density matrix as an MPS and you want to trace it, you can just trace the pair of physical indices on neighboring MPS tensors (I assume that bra and ket indices describing the same “site” or local Hilbert space are on neighboring MPS tensors?).
To do that, first contract the neighboring MPS tensors together. Then contract the result with a “delta” ITensor, like in this code example about tracing ITensors:
https://itensor.github.io/ITensors.jl/stable/examples/ITensor.html#Tracing-an-ITensor
Finally, you can contract all of the resulting ITensors together to get a scalar ITensor and call the scalar function on this result to extract a scalar value.
Please let me know if that solves your issue & if it works well in some test cases.
Hi Miles,
Thanks for your response. I am working on this equation of an open system under the dephasing noise
d\rho/dt=i[H,\rho] + \Gamma\sum_i ( S_z^i\rho S_z^i - \rho/4)
Here, i is the physical site index. So, by making \rho into a supermatrix, we have a superoperator
H_{supper}=iH\otimes1-i1\otimes H^T + \Gamma\sum_i (S_z^i\otimes S_z^i - 1/4)
The \otimes represents the cross product between ket and bra. In this problem, I am more interested in the evolution at intermediate time rather than the steady state at long time, so I use the non-Hermitian TDVP engine of Itensor. One thing I notice is that if I place all the bra tensor in one half and all the ket tensor in the other half of the chain, the performance is much better than placing the bra and associated ket tensor next to each other. Probably because my noise is only perturbative, while the main drive comes from the Hermitian part.
I already managed to run the code but I welcome if anyone has a better idea to improve the performance.
1 Like
I see, this is interesting. Glad you are already using TDVP and finding that to be working for your case.
I agree it’s non-trivial about the location or ordering of the bra and ket sites. It makes sense that for small \Gamma that separating them is best, whereas for large \Gamma having them next to each other is likely best.
Since they are separated spatially in your MPS, it could be more challenging to get the trace of your density matrix. Naively one would do it by either “swapping” the sites until bra/ket pairs are next to each other or by acting with an identity operator connecting the two far away sites. But those approaches could be expensive as a function of system size.
So here’s an idea for your setup: if you actually reverse the order of the bra sites relative to the ket sites, you may be able to efficiently trace out the pairs one at a time. So like if the order is 1,2,3,…, N-1, N, N, N-1 , …, 3,2,1 then you could first trace the bra index for Hilbert space number N with the ket index for numer N since they are already next to each other. This just leaves a tensor in the middle of the MPS on a bond that you can absorb into the MPS tensor next to it. Now the bra and ket indices associated with (N-1) are next to each other and you can repeat the process. Does that make sense?
Miles
Hi Duy, i’m doing the similar problem.could you explain more about the supervectoriztion of the MPS \rho in itensor language? thank you.
Hi Trudy,
So the idea is squeezing a 2^N x 2^N density matrix into a 2^(2N) column (notice the Hilbert space dimension is the same). Because of this, you can describe the density matrix by an MPS of 2N sites, half of them is for ket and the other half is for bra. You also need rewrite the Linbladian as well. I found this paper http://arxiv.org/abs/1804.09796 provide a very pedagogical introduction, hope it can help.
2 Likes
Hi Duy. Thank you so much for your help! I benefit a lot from this paper.
i followed the construction process of MPS |\rho\rangle\rangle, which is still a MPS of length N. The physical dimension and bond dimenion of each site is changed to its square power. Now i need to map the Hamilton H = -I \sum_{j=1}^{L-1}\sigma_j^z \sigma_{j+1}^z -g \sum_{j=1}^{L}\sigma_j^x to Liouville space too, which would need a construction of customised ITensor SiteType , where the operators are in form of I \otimes \sigma_x, I\otimes \sigma_z.
I noticed this process is a bit cumbersome. And it is different from your strategy
that describe |\rho\rangle\rangle by an MPS of 2N sites. So i was wondering if my understanding is wrong? Any of your suggestion would be appreciated. Thank you.
Hi Trudy,
Changing the local site dimension from 2 to 4 is equivalent to using two sites of dimension 2 (2 x 2 = 4). This is what I mean by doubling the number of site. This way you don’t need to define customized SiteType.
Now with 2N sites, you will label half of them to be “bra” spins and another half to be “ket” spins. Operators than only act on the density matrix from the left (H\rho) only acts on bra spin, operators act from the right side of \rho (\rho H) acts on ket spin, and quantum channels (L^\dagger \rho L) acts on both ket and bra spins. Hope this help.
Duy
1 Like
Hi Duy,
Thank you so much for these valuable advice!