Suppress population in a state throughout the trajectory

[will-atom]

I have a problem with a lossy intermediate state that I usually represent with a non-hermitian hamiltonian. I don't know if this problem warrants the full complexity of optimization of a dissipative gate as described in this example, but I would like to minimize the population in the lossy state throughout the optimization. I have tried creating a non-Hermitian Hamiltonian, but this resulted in errors from the optimizer, so I assume it's not compatible. Is it best to just use a Liouvillian approach?

[goerz]

I have tried creating a non-Hermitian Hamiltonian, but this resulted in errors from the optimizer, so I assume it's not compatible

No, that’s the best way to do it, and it absolutely should be supported

[will-atom]

Ok, thanks, I'll try that again and post my error message.

[will-atom]

Here is a typical result. If I omit the imaginary elements from the diagonal of the Hamiltonian (i.e. if it is hermitian) then the problem converges.

iter.        J_T     |∇J_T|       ΔJ_T   FG(F)    secs
     0   9.49e-01   7.80e-02        n[/a](https://file+.vscode-resource.vscode-cdn.net/a)    1(0)     0.2
     1   9.21e-01   2.22e-01  -2.77e-02    3(0)     0.7
     2   9.20e-01   3.40e-01  -1.26e-03    3(0)     0.7
     3   8.67e-01   4.68e-01  -5.26e-02    6(0)     1.6
GRAPE Optimization Result

- Started at 2023-05-21T14:52:24.362
- Number of objectives: 2
- Number of iterations: 3
- Number of pure func evals: 0
- Number of func[/grad](https://file+.vscode-resource.vscode-cdn.net/grad) evals: 15
- Value of functional: 8.67446e-01
- Reason for termination: Exception: TaskFailedException

    nested task error: AssertionError: Divided differences too small
    Stacktrace:
     [1] extend_newton_coeffs!(a::OffsetArrays.OffsetVector{ComplexF64, Vector{ComplexF64}}, n_a::Int64, leja::OffsetArrays.OffsetVector{ComplexF64, Vector{ComplexF64}}, func::QuantumPropagators.var"#33#35", n_leja::Int64, radius::Float64)
       @ QuantumPropagators.Newton [~/.julia/packages/QuantumPropagators/eKOEy/src/newton.jl:204](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/~/.julia/packages/QuantumPropagators/eKOEy/src/newton.jl:204)
     [2] newton!(Ψ::QuantumGradientGenerators.GradVector{5, Vector{ComplexF64}}, H::QuantumGradientGenerators.GradgenOperator{5, QuantumPropagators.Generators.Operator{SparseMatrixCSC{ComplexF64, Int64}, Float64}, QuantumPropagators.Generators.Operator{SparseMatrixCSC{ComplexF64, Int64}, Float64}}, dt::Float64, wrk::QuantumPropagators.Newton.NewtonWrk{QuantumGradientGenerators.GradVector{5, Vector{ComplexF64}}}; kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:func, :norm_min, :relerr, :max_restarts), Tuple{QuantumPropagators.var"#33#35", Float64, Float64, Int64}}})
       @ QuantumPropagators.Newton [~/.julia/packages/QuantumPropagators/eKOEy/src/newton.jl:300](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/~/.julia/packages/QuantumPropagators/eKOEy/src/newton.jl:300)
     [3] prop_step!(propagator::QuantumPropagators.NewtonPropagator{QuantumGradientGenerators.GradGenerator{QuantumPropagators.Generators.Generator{SparseMatrixCSC{ComplexF64, Int64}, ParametrizedAmplitude}, QuantumPropagators.Generators.Generator{SparseMatrixCSC{ComplexF64, Int64}}, Vector{Float64}}, QuantumGradientGenerators.GradgenOperator{5, QuantumPropagators.Generators.Operator{SparseMatrixCSC{ComplexF64, Int64}, Float64}, QuantumPropagators.Generators.Operator{SparseMatrixCSC{ComplexF64, Int64}, Float64}}, QuantumGradientGenerators.GradVector{5, Vector{ComplexF64}}})
       @ QuantumPropagators [~/.julia/packages/QuantumPropagators/eKOEy/src/newton_propagator.jl:121](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/~/.julia/packages/QuantumPropagators/eKOEy/src/newton_propagator.jl:121)
     [4] macro expansion
       @ [~/.julia/packages/GRAPE/uuiGq/src/optimize.jl:244](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/~/.julia/packages/GRAPE/uuiGq/src/optimize.jl:244) [inlined]
     [5] (::GRAPE.var"#151#threadsfor_fun#68"{GRAPE.var"#151#threadsfor_fun#53#69"{UnitRange{Int64}, Vector{Matrix{ComplexF64}}, Int64, Int64, Vector{Matrix{ComplexF64}}, Vector{Vector{ComplexF64}}, GrapeWrk{}}})(tid::Int64; onethread::Bool)
       @ GRAPE [./threadingconstructs.jl:84](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/threadingconstructs.jl:84)
     [6] #151#threadsfor_fun
       @ [./threadingconstructs.jl:51](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/threadingconstructs.jl:51) [inlined]
     [7] (::Base.Threads.var"#1#2"{GRAPE.var"#151#threadsfor_fun#68"{GRAPE.var"#151#threadsfor_fun#53#69"{UnitRange{Int64}, Vector{Matrix{ComplexF64}}, Int64, Int64, Vector{Matrix{ComplexF64}}, Vector{Vector{ComplexF64}}, GrapeWrk{}}}, Int64})()
       @ Base.Threads [./threadingconstructs.jl:30](https://file+.vscode-resource.vscode-cdn.net/home/will/us2qc_calculations/threadingconstructs.jl:30)
- Ended at 2023-05-21T14:52:27.892 (3 seconds, 530 milliseconds)

[goerz]

Ok, so your time propagation is unstable, which means there's something wrong with your Hamiltonian.

Are you using the :newton propagation method for the Hermitian case as well? If you're not specifying it explicitly, then you're probably using the default :cheby, not :newton. Make sure you can successfully run a propagation with method=:newton for the Hermitian case.

In general, before you start doing any optimal control, it's important to have a good understanding of the time propagation first, both in terms of physical behavior and numerical stability. You should analyze the dynamics of your Hamiltonian with a range of reasonable control fields.

How big are the imaginary elements? They should be on a similar order of magnitude as the non-imaginary elements. Try to vary them a little.

What is the spectral radius of the non-Hermitian Hamiltonian? What is the size of the time step? Ideally, both individually, and their product, should be of order 1e-3 to 1e1.

Are the control fields after iteration 3 still of reasonable amplitude?

Maybe try Krotov's method, too, which gives you more control over the step size in each iteration.

What is your system dimension?

A potentially more stable alternative to a non-Hermitian Hamiltonian is to write a custom propagator that does a normal Hermitian propagation step, but then simply removes the population in the levels you want to suppress. At some point, I should add support for the propagation callback to the optimal control methods, which would make this easier.

[wcairnx]

Specifying method=:cheby explicitly for both the Hermitian and non-Hermitian cases resolved my error, though I have not explored in detail why this is. Unfortunately I don't have bandwidth to dig into this now, but may be able to follow up at a later date. Thanks for your help!

[goerz]

Cheby requires a Hermitian Hamiltonian. I’m pretty sure it should blow up (most likely give you NaNs) if you try to use it with a non-Hermitian Hamiltonian. If you’re not seeing it blow up, maybe you’re not doing what you think you’re doing