What's even more worrying is that the expectation values of say "Sz" are also changing
However, it's not a bug of TenPy, but there is a simple explanation for this behaviour
The problem is that the ground state is degenerate (up to a finite size gap below the precision).
You chose the limit Jz=4, which corresponds to the strong Ising anisotropy, so the exact ground state is the Neel state with some very minor fluctuations.
But there are two Neel states |N1> = |ududud> and |N2> = |dududu>, relatet by the symmetry of flipping all spins.
For a finite system, only a superposition of the both states (+ fluctuations) is the exact ground states (consider the extreme case L=2, where you get the singlet |psi_0> = 1/\sqrt{2}(|N1> - |N2>) as exact ground state, and the triplet |psi_1> = 1/\sqrt{2}(|N1> + |N2>) has a higher energy.
However, the gap between the different superpositions of N1 and N2 goes down exponentially with system size, such that at some point the energy gap is below machine precision. Hence, DMRG can give you *any* superposition of the two Neel states as a "correct" answer.
The two superpositions have clearly different expectation values for "Sz", and also different entropies - you can add an arbitrary value between 0 and log(2) to the entropy by forming an appropriate superposition.
(If you have larger degeneracy, you can even add log(omega), where omega is the degeneracy of the ground state.)
One way to avoid these superpositions is to limit chi stronger.
In your case, you allowed a chi=512, which is much larger than necessary; the limititation of chi didn't contribute to the truncation at all.
If you limit e.g. chi=50, you get a single one of the two superposition states only. Let's say you need a bond dimension chi=M to represent (|N1>+fluctuations) being a ground state, and the same chi=M to represent (|N2>+fluctuations), then the superposition 1/sqrt(2)(|N1>+fluctuations + |N2>+fluctuations) requires a bond dimension chi=M*M = M^2.
In that way, the superpositions are suppressed by the truncation, and you end up with a single state.
The DMRG parameter
chi_list
allows to increase chi gradually during the DMRG simulation, to initially truncate strongly (suppressing superpositions), and later increase chi to reduce the truncation error (with an already symmetry broken state).
In particular for infinite DMRG calculations, this is a very good technique because the "environments" give a symmetry-broken boundary and hence a superposition is no longer possible. In the finite DMRG, this might not work so well, because the complete "environment" can be changed from one sweep to the other.
Alternatively, you can add a small term to the Hamiltonian to explicitly break the symmetry; e.g. in this case a staggered magnetic field. Playing around, I found that in this case a field
hz = np.array([1.e-6, -1.e-6]* (L//2))
was enough to fix the direction of the Neel state.
Of course, this requires to know the symmetry leading to the degeneracy. It might also change the physics, so be careful whether your "perturbation" still describes the same physical system