[ITensorMPS] [BUG] Default algorithm for adding MPSs leading to large bond dimensions

Thank you for your efforts in the development.

Description of bug

Adding an MPS with an MPS representing 0 with the default algorithm leads to a larger bond dimension than that obtained with the directsum algorithm

Minimal code demonstrating the bug or unexpected behavior

using ITensors
using Random

function zeromps(::Type{T}, sites) where {T<:Number}
    M = MPS(T, sites; linkdims=1)
    l = linkinds(M)
    for n in eachindex(M)
        if n == 1
            M[n] = ITensor(T, sites[n], l[n])
        elseif n == length(M)
            M[n] = ITensor(T, l[n - 1], sites[n])
        else
            M[n] = ITensor(T, l[n - 1], sites[n], l[n])
        end
        M[n] .= zero(T)
    end
    return M
end

L = 10
sites = [Index(2, "s=$n") for n in 1:L]

Random.seed!(1234)
A = zeromps(ComplexF64, sites)
B = randomMPS(ComplexF64, sites; linkdims=10)
AB = +(A, B; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB)
truncate!(AB)
@show maxlinkdim(AB)

AB_ref = +(A, B; alg="directsum")
truncate!(AB_ref; cutoff=1e-25, maxdim=typemax(Int))
@show maxlinkdim(AB_ref)

Expected output or behavior

Since A=0, we should be able to restore the original bond dimension (=10).

Actual output or behavior

The bond dimension of the result is much larger than the expected one, even after an explicit truncation.

maxlinkdim(AB) = 51
maxlinkdim(AB) = 25
maxlinkdim(AB_ref) = 10

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M2
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_NUM_THREADS = 1
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_PROJECT = @.

• Output from using Pkg; Pkg.status("ITensors"):

julia> using Pkg; Pkg.status("ITensors")
  [9136182c] ITensors v0.6.21

Start by running the minimal Julia example and comparing the default +(A, B) result with +(A, B; alg="directsum"). Trace the MPS addition path used by the default algorithm; done means adding the zero MPS preserves B's bond dimension of 10 and explicit truncation agrees with the directsum result.