W-state as MPS

[Robin123462]

Hello,

I'm trying to implement the W-state (https://en.wikipedia.org/wiki/W_state) as an MPS. My first approach was the MPS.from_Bflat() function, but I'm not quite sure which shape the Matrix should have. I've tried for example the State 1/sqrt(3) (|100> + |010> + |001>) by the three matrices
A_1 = [[[0,1],[0]],[[0,1],[1,0]]] A_2 = ... (as in https://en.wikipedia.org/wiki/Matrix_product_state) but that did not work.

My second approach is to run DMRG on the XX Model (I used the XXZChain with Jz = 0, h = 0) since the groundstate of the XX Model should be the W state, but then I'm not quite sure how to check if this is the state I'm looking for. Is there any method to get the resulting state in the form of basis states like 1/sqrt(3) (|100> + |010> + |001>)? I tried psi.get_theta(0,3) but wasn't really sure how to interpret the output.

Thanks a lot.

Robin

[Johannes]

You can write this analytically as a bond-dimension 2 MPS.

You can see that you need bond dimension 2, since you always have either exactly one or no 1 on the left/right of any of the two cuts (between the first two sites, or the second two sites).
The Schmidt decomposition between the last two sites is given as
\[\ket{\psi} = \frac{\sqrt{2}}{\sqrt{3}} \frac{\ket{1 0} + \ket{01}}{\sqrt{2}} \otimes \ket{0} + \frac{1}{\sqrt{3}} \ket{0 0} \otimes \ket{1} \]
Here, I've already normalized the left/right Schmidt states, hence the \(\sqrt{2}\) which cancels, and ordered Schmidt values (left-most constants in the sum) large to small. This should correspond to the representation of the MPS in A_0 A_1 S_2 B_2 canonical form, so we can read off the Schmidt values
S_2 = [sqrt(2)/sqrt(3), 1/sqrt(3)] and the last B_2 as identity from the left virtual to the physical leg with trivial right leg,
B_2[p=0] = [[1], [0]]; B_2[p=1] = [[0], [1]].

By symmetry, the Schmidt decompostion cutting between the first two sites reads
\[\ket{\psi} = \frac{\sqrt{2}}{\sqrt{3}} \ket{0} \otimes \frac{\ket{1 0} + \ket{0 1}}{\sqrt{2}} + \frac{1}{\sqrt{3}} \ket{1}\otimes \ket{0 0} \]
so again, the Schmidt values are
S_1 = [sqrt(2)/sqrt(3), 1/sqrt(3)], and we can read off the left-most A_0 as
A_0[p=0] = [[1, 0]]; A_0[p=1] = [[0, 1]].

The remaining, somewhat tricky part is then to read of the B_1, which is a 2x2 matrix for each fixed physical p on site 0.
It should be an isometry writing the two Schmidt states \(\frac{\ket{1 0} + \ket{0 1}}{\sqrt{2}}\) and \(\ket{00}\) from the left decomposition in terms of the local basis and the right Schmidt basis \(\ket{0}\) and \(\ket{1}\).
Hence, it reads:
B_1[p=0] = [[0, 1/sqrt(2)], [1, 0]]; B_1[p=1] = [[1/sqrt(2), 0], [0, 0]]

That gives us the MPS in form A_0 S_1 B_1 B_2. If you want to write it with MPS.from_B, you can also transform the first A tensor into an B tensor by using
B_0 = S_0^{-1} A_0 S_1. Since we have trivial bond dimension 1 on the left of A_0, the Schmidt value S_0 = [1] is trivially 1 (the norm of the state), and we get B_0[p=0] = [[sqrt(2)/sqrt(3), 0]]; B_0[p=1] = [[0, 1/sqrt(3)]].

Putting things together in TeNPy, we get

Code: Select all

from numpy import sqrt
from tenpy.networks.mps import MPS
from tenpy.networks.site import SpinHalfSite

B_0 = np.zeros([2, 1, 2])  # p vL vR
B_1 = np.zeros([2, 2, 2])
B_2 = np.zeros([2, 2, 1])
B_0[0] = [[sqrt(2)/sqrt(3), 0]]
B_0[1] = [[0, 1/sqrt(3)]]
B_1[0] = [[0, 1/sqrt(2)], [1, 0]]
B_1[1] = [[1/sqrt(2), 0], [0, 0]] 
B_2[0] = [[1], [0]]
B_2[1] = [[0], [1]]

S = [[1.], [sqrt(2)/sqrt(3), 1/sqrt(3)], [sqrt(2)/sqrt(3), 1/sqrt(3)], [1.]]
psi = MPS.from_Bflat([SpinHalfSite(conserve=None)]*3, [B_0, B_1, B_2], S, form='B')

Indeed, you can check the output of get_theta:

Code: Select all

psi.get_theta(0, 3).to_ndarray().reshape(2**3)
# gives
# array([0.        , 0.57735027, 0.57735027, 0.        , 0.57735027,
#        0.        , 0.        , 0.        ])

The 8 states here are orderd due to numpy reshaping rules as \(\ket{000}, \ket{0 0 1} \ket{010} \ket{011} \ket{100} \ket{101} ...\), so we see that we constructed indeed the correct W state.

PS: Note that if you add charge conservation as SpinHalfSite(conserve='Sz', sort_charge=True), the option permute=True of from_Bflat will cause it to permute the local basis, and the theta will look the wrong way arround (exchanging 0 <->1 states on each sites). Note for the future: This issue will be fixed in TenPy version 2.0 with the rewrite of the np_conserved part.

[Robin123462]

Hello Johannes,

thanks a lot for the quick and detailed reply, it helps me a lot.
Unfortunately, another error occurred and I'm unsure whether it's a mistake on my part, maybe you can help me there too.
I've run the following model:

Code: Select all

    class XYModel(CouplingModel, MPOModel): 
        def __init__(self, model_params):
            L = model_params.get('L')
            site = SpinHalfSite(conserve='parity',sort_charge=False)
            lat = Chain(L, site,bc='periodic', bc_MPS='finite')
            CouplingModel.__init__(self, lat)
            for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
                self.add_coupling(-1., u1, 'Sigmax', u2, 'Sigmax', dx)
                self.add_coupling(-1., u1, 'Sigmay', u2, 'Sigmay', dx)
            MPOModel.__init__(self, lat, self.calc_H_MPO())

and for L=3 I get the right result

Code: Select all

psi[1].get_theta(0,3).to_ndarray().reshape(2**3))
#[ 0.00000000e+00 -5.77350269e-01 -5.77350269e-01  0.00000000e+00
#-5.77350269e-01  0.00000000e+00  0.00000000e+00 -1.28197512e-16]

but for L=4 my result is:

Code: Select all

psi[1].get_theta(0,4).to_ndarray().reshape(2**4))
#[ 0.          0.          0.         -0.35355339  0.         -0.5
#-0.35355339  0.          0.         -0.35355339 -0.5         0.
# -0.35355339  0.          0.          0.        ]

which seems different to the W-state for 4 qubits (1/sqrt(4) (|0001> + |0010> ...)), do you have any Idea why is that?
Again thanks for your help.

[Johannes]

The XX model has charge conservation of total Sz. In the zero-particle sector, you can map it exactly to a single-particle fermion hopping, and the ground state is the W state. However, this is not true in the sector of higher filling when you have peridoic boundary conditions: for the coupling going from site 0 to site L-1, you need the full Jordan-wigner string inbetween.

The effect of this is that the true ground state of the XX chain is in the half-filling sector (same number of up spins as down spins), it has lower energy than the W state. You can check that the W state still is an eigenstate, but at higher energies.

To get the W state for larger L, you can either continue to follow your strategy and use the ground state of the XX model in the sector with just one up spin, for which you just need to turn on Sz conservation and use an initial MPS with 1 up spin only.
Alternatively, you can generalize the analytic construction of the MPS.

In general, for a single-particle wave function you can construct the MPS analytically.
A simple way for the general-L construction would be as follows (this is a generic construction for single-site wave functions for finite systems):
the state we want to write it is \((\frac{1}{\sqrt(L)} \sum_i S^+_i ) \ket{\downarrow \cdots downarrow}\),
look up how to construct the operator \((\frac{1}{\sqrt(L)} \sum_i S^+_i )\) as an MPO with bond dimension 2, namely with
W matrices W = [[Id, Sp], [0, Id]] . The initial down product state trivially removes one physical leg of the MPO, so we see that we can directly take the W matrices as MPS tensors, although not in canonical form. But you can still initialize that as an MPS in tenpy, passing form=None to MPS.from_Bflat, and then simply call psi.canonical_form().

Actually, it would be nice to have this as a function in the MPS, something like MPS.from_single_particle_wave_function() - I'd be happy to merge a pull request for this

Edit: See Issue #312