Mixing of LockedAmplitude and ShapedAmplitude

[mabuni1998]

I have an optimization problem where I have three time-dependent fields. The first field is my control, the second is a LockedAmplitude (not optimized), and the third should be the product of the two. However, it is not clear to me how I should define that?

As an example, it could be:

Tend = 10
times = range(0, Tend, step=0.01)
shape(t) =  QuantumControl.Shapes.flattop(t, T=Tend, t_rise=0.3, func=:blackman);
optimized_field = QuantumControl.Controls.ShapedAmplitude(t -> 1, times; shape) #Should be optimized

constant_ampl = QuantumControl.Controls.LockedAmplitude(t->t^2, times) #Amplitude that should not be optimized. In reality the time dependence is more complicated

#Cross term that is the product of the two fields
cross_term = optimized_field*constant_ampl #This obviously does not work, but I am unsure how to define an equivalent function.

#The Hamiltonian should look like
H = hamiltonian((h1,constant_ampl),(h2,optimized_field),(h3,cross_term)) 

I considered whether I should refer to optimized_field.control, since this vector is mutated in the optimization, but it is not entirely clear to me whether this will work.

[goerz]

If we call your constant_ampl $\Gamma(t)$, the constant shape in the optimized_field $S(t)$ and the actual control t -> 1 $\epsilon(t)$, it seems to me your Hamiltonian is

$$ \hat{H} = \Gamma(t) \hat{H}_1 + S(t) \epsilon(t) \hat{H}_2 + \Gamma(t) S(t) \epsilon(t) \hat{H}_3 $$

The product $\Gamma(t) S(t)$ is just another shape function for a ShapedAmplitude. So this should be totally fine, as something like:

ϵ(t) = 1
S(t) =  QuantumControl.Shapes.flattop(t, T=T, t_rise=0.3, func=:blackman)
Γ(t) = t^2
optimized_field = QuantumControl.Controls.ShapedAmplitude(ϵ; shape=S)
constant_ampl = QuantumControl.Controls.LockedAmplitude(Γ)
cross_term = QuantumControl.Controls.ShapedAmplitude(ϵ; shape=(t -> S(t) * Γ(t)))
H = hamiltonian((h1,constant_ampl),(h2,optimized_field),(h3,cross_term))

Note that I didn't include times when instantiating the amplitudes, so that the time discretization is delayed. Otherwise, you run into a problem that ϵ is discretized twice and you end up with two controls (always double-check your full Hamiltonian with get_controls to make sure you have the expected numbers of controls).

[will-atom]

Answering my own question: yes

shape(t) = flattop(t, T=tmax, t_rise=t_rise)
ϵ = ShapedAmplitude(t -> 1.0; shape=shape)
ϵ_sq = LockedAmplitude(t->ϵre(t)^2)

[will-atom]

A question more related to the one from @mabuni1998: What if Gamma(t) was also a control, so the cross-term is the product of 2 controls? Can that be done?

[goerz]

Yes! Define an amplitude that is the product of two controls and a static shape.

Everything is based on the distinction between amplitudes and controls, so make sure you understand the definitions in the glossay. There are absolutely no restrictions on how any number of controls or other time-dependencies can appear in the Hamiltonian.

[mabuni1998]

I think it would be very helpful with some examples in the documentation where it is shown how some of these more exotic examples of control amplitudes can be constructed. Even knowing the glossary, and reading that the control amplitude can have any dependence on one or more controls, I would not find it trivial to define such a control amplitude. The above example from @will-atom helps, but from the documentation of LockedAmplitude, it seems that it should not be used together with controls: "A time-dependent amplitude that is not a control." The only other control amplitude in the documentation is the ShapedAmplitude, but from here a control amplitude where two controls are mixed is for example not shown (probably it is trivial to you, but it might not be to a new user).

Anyway, I think the functionality of the package is fantastic, but I would love to be better at using it, which currently is a steep learning curve and that is why I am suggesting having a section in the documentation where some of these examples of more complicated control amplitudes is shown.

Finally, for the case where Gamma(t) is also a control,I would also be very interested in how this could be defined.

[goerz]

I think it would be very helpful with some examples in the documentation where it is shown how some of these more exotic examples of control amplitudes can be constructed.

I completely agree! I would love to have some tutorials of some more exotic examples. But this will, at least in part, depend on contributions, since I don't have any such examples myself.

it seems that it should not be used together with controls: "A time-dependent amplitude that is not a control."

Well, that's certainly not the intention. A time-dependent-amplitude not being a control (or probably that should be "containing") has absolutely no implications for controls that are present in other amplitudes, in other terms of the Hamiltonian. I'm always looking for ways to clarify the documentation, so any suggestions would be helpful.

probably it is trivial to you, but it might not be to a new user

No, it's not trivial. In fact, quite the opposite: I don't have any examples myself for such exotic controls. That's why there are no tutorials or examples.

It's also not trivial from the technical perspective: the distinction between amplitudes and controls is what I would consider one of the major "innovations" of this framework, compared to existing quantum control packages. But the motivation for this is "flexibility", in that it should be possible to describe any possible control system. That does not mean that I have any examples for specific "exotic" schemes myself. The guiding principle is "bring your own Hamiltonian" (also, "bring your own functionals", and maybe even bring your own propagators. There is definitely the expectation that people implement their own problem-specific Hamiltonians or amplitudes. That does imply a learning curve, which I'm okay with, up to a point. There should still be plenty of documentation, and ideally, examples. But I'm pretty limited in the examples I can provide, so these must be contributed. It would be extremely helpful if you had a write-up with a description of the physical system (the Hamiltonian) and the control problem you are considering.

Finally, for the case where Gamma(t) is also a control,I would also be very interested in how this could be defined.

This involves some significant work. You'd start by defining a custom struct like

struct ShapedProductAmplitude
    controls # vector of controls
    shape
end

to be instantiated as ShapedProductAmplitude((Γ, ϵ), S), and then all the methods specified in QuantumPropagators.Interfaces.check_amplitude and QuantumControl.Interfaces.check_amplitude need to be implemented (and verified via by calling check_amplitude).

[mabuni1998]

I see that makes a lot of sense! The Hamiltonian I am currently interested in optimizing has the following form:

$H_{\mathrm{red}} / \hbar=\sum_m \omega_m^c(t) a_m^{\dagger} a_m+\sum_{n, m} C_{n m} \sqrt{\kappa_n^*(t) \kappa_m(t)} a_n^{\dagger} a_m+\chi_3 \sum_m a_m^{\dagger} a_m^{\dagger} a_m a_m,$

We are thus considering a certain amount of photonic modes (with the associated creation and annihilation operator $a_m$). We vary their detuning or frequency $\omega_m(t)$ as a control signal as well as their coupling rates $\kappa_n(t)$. Notice the form of the coupling term, where we have a constant $C_{mn}$, which ideally could also be optimized (see more on this point below). But also the two control signals $\kappa_m(t)$ and $\kappa_n(t)$ combine under a square-root. Finally each cavity contains a chi-three nonlinearity which is not optimized. An example goal could then be to consider four cavities with a dual rail encoding so that $|00\rangle_L = |1,0,1,0\rangle$, $|10\rangle_L = |0,1,1,0\rangle$, $|01\rangle_L = |1,0,0,1\rangle$, $|11\rangle_L = |0,1,0,1\rangle$, where $|0,0\rangle_L$ denotes the logical state and $|1,0,1,0\rangle$ the occupation of the four physical cavities. One could then try to implement a C-phase gate, so that $|11\rangle_L \rightarrow -|11\rangle_L$ and all other states are unchanged.

I have sketched out an implementation of this squareroot product of two functions below, which passes both checks. I haven't tested whether it works in an optimization yet. Do you have any comments on how it was implemented?

using QuantumControl
using QuantumPropagators

abstract type SqrtProductAmplitude <: QuantumPropagators.Amplitudes.ControlAmplitude end


function SqrtProductAmplitude(control1,control2)
    if (control1 isa Vector{Float64}) && (control2 isa Vector{Float64})
        return SqrtProductAmplitudePulse((control1, control2))
    else
        try
            ϵ_t = control1(0.0)
        catch
            error(
                "SqrtProductAmplitude control 1 must either be a vector of values or a callable"
            )
        end
        try
            ϵ_t = control2(0.0)
        catch
            error("SqrtProductAmplitude control 2 must either be a vector of values or a callable")
        end
        return SqrtProductAmplitudeContinuous((control1, control2))
    end
end

struct SqrtProductAmplitudePulse <: SqrtProductAmplitude
    controls # vector of controls
end

struct SqrtProductAmplitudeContinuous <: SqrtProductAmplitude
    controls # vector of controls
end



Base.Array(ampl::SqrtProductAmplitudePulse  ) = sqrt.(ampl.controls[1] .*ampl.controls[2] )
(ampl::SqrtProductAmplitudeContinuous)(t::Float64) = sqrt(ampl.controls[1](t) * ampl.controls[2](t))

function QuantumPropagators.Controls.evaluate(ampl::SqrtProductAmplitude, args...; vals_dict=IdDict())
    S_t = QuantumPropagators.Controls.evaluate(ampl.controls[1], args...; vals_dict)
    ϵ_t = QuantumPropagators.Controls.evaluate(ampl.controls[2], args...; vals_dict)
    return sqrt.(S_t .* ϵ_t)
end

function Base.show(io::IO, ampl::SqrtProductAmplitude)
    print(io, "SqrtProductAmplitude(::$(typeof(ampl.controls[1])),::$(typeof(ampl.controls[2]))")
end


function QuantumControl.Controls.get_control_deriv(ampl::SqrtProductAmplitudePulse,control)
    if control ≡ ampl.controls[1]
        return QuantumControl.Amplitudes.LockedAmplitude(.- 1 ./ (2 .* sqrt(ampl.controls[1] .* ampl.controls[2])) .* ampl.controls[2])
    elseif control ≡ ampl.controls[2]
         return QuantumControl.Amplitudes.LockedAmplitude(.- 1 ./ (2 .* sqrt(ampl.controls[1] .* ampl.controls[2])) .* ampl.controls[1])
    else
        0
    end
end

function QuantumControl.Controls.get_control_deriv(ampl::SqrtProductAmplitudeContinuous,control)
    if control ≡ ampl.controls[1]
        return QuantumControl.Amplitudes.LockedAmplitude(t->-1/(2*sqrt(ampl.controls[1](t)*ampl.controls[2](t)))*ampl.controls[2](t))
    elseif control ≡ ampl.controls[2]
         return QuantumControl.Amplitudes.LockedAmplitude(t->-1/(2*sqrt(ampl.controls[1](t)*ampl.controls[2](t)))*ampl.controls[1](t))
    else
        0
    end
end



QuantumControl.Controls.get_controls(a::SqrtProductAmplitude) = a.controls

function QuantumPropagators.Controls.substitute(ampl::CT, replacements) where {CT<:SqrtProductAmplitude}
    if ampl in keys(replacements)
        return replacements[ampl]
    else
        control1 = QuantumPropagators.Controls.substitute(ampl.controls[1], replacements)
        control2 = QuantumPropagators.Controls.substitute(ampl.controls[2], replacements)
        CT((control1, control2))
    end
end


ϵ1(t) = 1.0
ϵ2(t) = 1.0

field = SqrtProductAmplitude(ϵ1, ϵ2)

dt = 0.1
T = 10
tlist = collect(range(0, T, step=dt))

QuantumPropagators.Interfaces.check_amplitude(field; tlist, quiet=false)
QuantumControl.Interfaces.check_amplitude(field; tlist, quiet=false)

As for the constant factor in a Hamiltonian, lets say we have $H= k * O$, where k is a constant we want to optimize and $O$ is an operator. If we then want to compute the derivate of some fidelity $F = \langle \psi_T | \psi_F \rangle$, where $\psi_T$ is the target state and $\psi_F = \prod_n U_n |\psi_i \rangle$ is the final state after some evolution of an initial state via the evolution operators $U_n$. If we look at the derivative with regards to k, we then have (I might have the numbering of the U's backward below):

$$ \begin{aligned} \frac{\partial F}{\partial k}= & \left\langle\psi_T \left\lvert\ \frac{\partial U_1}{\partial k} U_2 \ldots\right \lvert \psi_i \right \rangle+\left\langle\psi_T\right| U_1 \frac{\partial U_2}{\partial k} U_3 \ldots .|\psi_i\rangle \\ & +\left\langle\psi_T\right| U_1 U_2 \frac{\partial U_3}{\partial k} U_4 \ldots\left|\psi_i\right\rangle + \cdots \end{aligned} $$

This is exactly the sum over n of $\nabla\tau_{ln}^{(k)}$ of eq. 3 in your Quantum Journal paper of 2022 about the framework. So if we have a "constant control" we should simply sum the contribution of the derivative of each of the discretizations. I don't have a complete overview of the framework, but I would imagine that this is not a fundamental challenge but rather a software implementation challenge. I would be very interested if you have any inputs or ideas on how the above could be included in the framework. Thanks!