Problem with TEBD on Fermi-Hubbard model with Electron indices

I am trying to prepare the ground state of a simple 1d Fermi-Hubbard model and then I want to TEBD it. Ideally I would quench it with the evolution given by the Hamiltonian at different repulsion $U$, but at the moment I am trying to do it for the simple case of same $U$ where the ground state should remain unchanged by Hamiltonian dynamics (up to small errors in Trotter).
However, if I use the dagger operators with the Jordan-Wigner strings attached (like Cdagup, Cdn etc), I don't get this conservation of energy and fidelity wrt. to the initial state.
Changing this to the JW-less operators A... fixes the issue, but this leaves me very suspicious. What am I doing wrong in my implementation?
I am attaching a MWE of my issue, where changing ladder_op_type = "C" to ladder_op_type = "A" reproduces what I described.

using ITensors, ITensorMPS

function hopping_gate(sites, i, j, spin, theta)
    ladder_op_type = "C"
    hj = op("$(ladder_op_type)dag$spin", sites[i])* op("$(ladder_op_type)$spin", sites[j]) + op("$(ladder_op_type)dag$spin", sites[j])* op("$(ladder_op_type)$spin", sites[i])
    return exp(-im * theta * hj)
end

let 

#number of sites, interactions etc 
N_spinful_sites = 8
U = 1.
t = 1.
dτ = 0.04

#choose between conserving quanties or not 
conserve = true
if conserve
    sites = siteinds("Electron", N_spinful_sites; conserve_qns=true)
else
    sites = siteinds("Electron", N_spinful_sites; conserve_qns=false)
end


#Build Hamiltonian and Trotter layer 

ladder_op_type = "C"

os = OpSum()
trotter_layer::Vector{ITensor} = []

#up hopping 
for j in 1:(N_spinful_sites - 1)
    os += -t, "$(ladder_op_type)dagup", j, "$(ladder_op_type)up", j + 1
    os += -t, "$(ladder_op_type)dagup", j + 1, "$(ladder_op_type)up", j
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "up", -dτ*t/2))
end

#down hopping
for j in 1:(N_spinful_sites - 1)
    os += -t, "$(ladder_op_type)dagdn", j, "$(ladder_op_type)dn", j + 1
    os += -t, "$(ladder_op_type)dagdn", j + 1, "$(ladder_op_type)dn", j
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "dn", -dτ*t/2))
end

#onsite repulsion 
for j in 1:(N_spinful_sites)
    os += U, "Nupdn", j 
    push!(trotter_layer, exp(-im * dτ * U * op("Nupdn", sites[j])))
end

#down hopping
for j in 1:(N_spinful_sites - 1)
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "dn", -dτ*t/2))
end 
#up hopping
for j in 1:(N_spinful_sites - 1)
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "up", -dτ*t/2))
end

#build Hamiltonian 
H = MPO(os, sites)

#build initial state for DMRG
if conserve 
    Nup = 3
    Ndown = 3
    state = fill("Emp", N_spinful_sites)
    @assert Nup + Ndown <= N_spinful_sites
    for j=1:Nup
        state[j] = "Up"
    end
    for j=(Nup+1):(Nup+Ndown)
        state[j] = "Dn"
    end

    psi0 = random_mps(sites, state)
else
    psi0 = random_mps(sites)
end
@show sum(expect(psi0, "Nup")), sum(expect(psi0, "Ndn"))

#DMRG
maxdim = [20, 40, 80, 160, 500]
cutoff=1.e-12
nsweeps = 2*length(maxdim)
energy, psi_gs = dmrg(H, psi0; nsweeps, maxdim, cutoff)
@show sum(expect(psi_gs, "Nup")), sum(expect(psi_gs, "Ndn"))


#dynamis 
χ = 2000
nlayers = 10
psi = deepcopy(psi_gs)
@show typeof(psi)
for k=1:nlayers
    psi = ITensorMPS.apply(trotter_layer, psi; maxdim=χ, cutoff=1.e-18)
    energy_val = real(inner(psi, Apply(H, psi)))
    overlap_with_gs = abs(inner(psi_gs, psi))^2
    println("<psi|H|psi> = $energy_val, <psi_gs|psi> = $overlap_with_gs, maxlinkdim(psi) = $(maxlinkdim(psi))")
end


end

(The only similar issue I could find was this older comment that never got solved)

Run the attached MWE with ladder_op_type set to "C" and "A"; start at hopping_gate and the construction of H and trotter_layer, then compare the TEBD energy and overlap with psi_gs. Done means the C-operator evolution matches the stated same-U expectation without relying on the A operators.