Normalization Issue in the semiclassical solvers

[Zhangpengfeiscu]

I am attempting to calculate the ground state of a semiclassical system using imaginary-time evolution within the semiclassical.schroedinger_dynamic function. I am facing difficulties ensuring the normalization of the wavefunction.

I have referenced the examples provided in the QOJulia documentation for vortices (https://docs.qojulia.org/examples/vortex/) and spin-orbit-coupled BEC (https://docs.qojulia.org/examples/spin-orbit-coupled-BEC1D/). Based on these examples, I have attempted two methods to incorporate normalization operations into my function:

1. Direct Normalization in Function Definitions:

function fquantum(t, psi, u) 
    normalize!(psi)
    # build Hamiltonian H
    return -1im*H
end

function fclassical!(du, u, psi, t) 
    normalize!(psi)
    # build classical equation
end

tout, ρt = semiclassical.schroedinger_dynamic(tlist, ψsc0, fquantum, fclassical!);

2. Using a Renormalization Callback:

# renormalization callback
norm_func(ψsc0, t, integrator) = normalize!(ψsc0)
ncb = FunctionCallingCallback(norm_func; func_everystep = true, func_start = false)

abstol_int = 1e-4
reltol_int = 1e-4
maxiters_int = 1e8

tout, ρt = semiclassical.schroedinger_dynamic(tlist, ψsc0, fquantum, fclassical!, callback=ncb, abstol=abstol_int, reltol=reltol_int, maxiters=maxiters_int);

Despite these efforts, the wavefunction obtained from both methods does not appear to be normalized.
Could someone please advise on how to properly implement wavefunction normalization within the semiclassical.schroedinger_dynamic function?
Thank you for your assistance.

[david-pl]

Internally, the solver works with just a vector of complex numbers. That's also true for a Schrödinger equation without using semiclassical dynamics. However, there is a slight difference in how data is copied over to that vector of numbers. Essentially, when you're doing normalize!(psi) in your first approach, you're normalizing a copy of the state rather than the actual state. This is very unfortunate since there's no way for you to know that, so I'm sorry that this tripped you up, but I don't think anyone had this use case before.

I'm not sure what the best way to fix it yet is, but in the meantime your second approach should actually work. But you have to be aware of how the data is structured since you're working with raw data, not quantum optics objects. Essentially a semiclassical state just consists of a Ket and a vector of numbers representing the classical part. When copied over, the state in the integrator is then just a vector of numbers, where the data is joined. So you have a vector of the form such as [ψsc.quantum.data; ψsc.classical].

You can adapt your normalization callback to only normalize the quantum part. Something like this should work:

const n = length(ψsc0.quantum.data)  # length of the quantum part of the data
function norm_func(integrator_state, t, integrator)
    quantum_part = @view integrator_state[1:n]
    normalize!(quantum_part)
end

Can you try this and let me know if it works?

[Zhangpengfeiscu]

@david-pl Thank you for your helpful advice regarding the normalization issue in our solver. Your explanation was clear and the proposed solution with the normalization callback was effective.
I've implemented the changes, and everything is now working as expected. Thank you once again for your valuable assistance