I am trying to compute the first excited state of the Ising model in Julia and the DMRG converges without a problem for the ground state for any system size N but fails to do so for N > 10 for the first excited state. It does converge for the first excited state for up to N = 10. I have not tested N between 10 and 20 but then for N = 20 for example the ground state converges as expected and then the first excited state energy blows up immediately and is much more negative than the ground state.

I am having the same problem for the Schwinger model and there the situation is even worse, maybe because the Hamiltonian is more complicated, but there the first excited state algorithm fails to converge for N > 2 and in fact there is such a big numerical instability that I get errors of the form: ERROR: LoadError: operator does not appear to be hermitian: 17.330705523490913 vs 2.772575076898464e9.

I am checking my ansatz for the first excited state and has very small overlap with the ground state and the same is true for the result of the first excited state with the ground state even if the energy of the first excited state blew up beyond the energy of the ground state.

I give below a full working code for the Ising model which ends up giving a ground state energy of -4.471756355671585 and a first excited energy of -81856.75774379515. I have tried increasing the bond dimension, or adding some noise or increasing the maximum number of sweeps but nothing worked.

Is there anything you can suggest?

```
using ITensors
function get_Ising_OpSum(N, J, g_z, g_x)
ampo = OpSum()
for n=1:N-1
ampo += 4*J,"Sz",n,"Sz",n+1
end
for n=1:N
ampo += 2*g_x,"Sx",n
ampo += 2*g_z,"Sz",n
end
return ampo
end
function get_MPO_from_OpSum(OpSum, sites)
return MPO(OpSum, sites)
end
N = 20
J = 0.0001
g_z = 0.1
g_x = 0.2
ns = 500
D = 20
sites = siteinds("S=1/2", N)
opsum_ising = get_Ising_OpSum(N, J, g_z, g_x)
H = get_MPO_from_OpSum(opsum_ising, sites)
initial_ansatz_0 = randomMPS(sites, D)
sweeps = Sweeps(ns, maxdim = D)
energy_0, psi_0 = dmrg(H, initial_ansatz_0, sweeps, ishermitian = true, maxdim = D)
# Ms = [psi_0]
# w = energy_0
# initial_ansatz_1 = randomMPS(sites, D)
# energy_1, psi_1 = dmrg(H, Ms, initial_ansatz_1, sweeps, weight = w, ishermitian = true, maxdim = D)
Heff = H + energy_0.*outer(psi_0', psi_0)
initial_ansatz_1 = randomMPS(sites, D)
println("Overlap of ansatz with gs: ", inner(psi_0, initial_ansatz_1))
initial_noise = 1e-2
noise_vector = LinRange(initial_noise, 0.0, ns)
energy_1, psi_1 = dmrg(Heff, initial_ansatz_1, sweeps, ishermitian = true, maxdim = D, noise = noise_vector)
println("Overlap of 1st excited state with gs: ", inner(psi_0, psi_1))
println(energy_0)
println(energy_1)
```