Refactor indexing

[david-pl]

This lays the foundation for refactoring the entire indexing implementation. The basic idea is one we already discussed in the very beginning. Essentially, we want to use the following commutation relation for indexed bosonic operators

$$ a_i a_j^\dagger = a_j^\dagger a_i + \delta_{ij} $$

Similarly, for transition operators, we have

$$ \sigma_r^{ij} \sigma_s^{k\ell} = \delta_{jk}\delta_{rs}\sigma_r^{i\ell} + (1 - \delta_{rs})\sigma_r^{ij}\sigma_s^{k\ell} $$

The problem here is the fact that the original product of transition operators on the left-hand side appears again on the right-hand side. This can potentially lead to infinite recursion. Currently, the solution to this problem is to store the fact that in the latter product $r \neq s$ on any sum expression that may occur. This has other problems, mostly the complex implementation and handling all kinds of special cases.

Instead, here I propose to implement the commutation relations using the above relations directly. To prevent infinite recursion on the transition operators, we store an additional field on indexed operators called merge_events, which is just a Vector{UUID}. Whenever two indexed transitions get rewritten according to the above, we generate a random UUID, currently using UUIDs.uuid4(), and store that on two copies of the original transition operators. Therefore, the operators behave just as they should when multiplied together with other operators, except when we find a shared UUID in the merge_events. That can only occur when a product of indexed transitions was returned from the * implementation and thus they have already been merged.
To ensure that in products involving three or more indexed transitions, we also sort by the length of the merge_events vector.

Since this needs to be a full rewrite, I'm opening this as a draft until we add back all the features. If we don't complete the rewrite within this PR, it at least serves as documentation of the fundamental concept, which seems to work.

Feature list:

• Basic indexed operators and commutation relations
• Summation
• Conversion to MTK and numeric solution
• Scaling (replacing sums by appropriate factors)
• Correlation functions and spectra together with indexing
• Measurement back action with indexing

FYI, @ChristophHotter and @j-moser

[david-pl]

Little update on the status here: derivation of equations, averaging and the cumulant expansion all work so far. I haven't checked the results in detail, but the equations look okay. I'll check that by comparing to a "brute-force" approach later, when I can actually get a numerical solution out. The next step would be to implement complete.

A couple of things I've learned so far:

• The basic idea here, i.e. implementing commutation relations for indexed operators directly via multiplication (not bothering with a concept of "not equal indices" on sums etc.) seems to work. It makes the code way shorter and easier to read. It's still messy right now (leftover TODOs and a bunch of old commented out code), but it should be easy to clean that up.
• Indices are always contained in SymbolicUtils.arguments(ex) now, except for IndexedOperators. This makes finding and swapping indices so much easier.
• Indices are now c-numbers, so you can do some math on them. This is useful, because it allows to express $\delta_{ij}$ simply as i == j. Similarly, instead of $1 - \delta_{ij}$ you can just write i != j (realizing I could do that took me way longer than I'd like to admit).
• We still need a nice usable way to construct a product of transition operators such as $\sigma_i^{eg}\sigma_j^{ge}$ where $i \neq j$. Right now that's not easily possible, but you may want that e.g. when adding such a product to the list of operators you want to derive an equation for. Otherwise, you always end up with an addition of the two different cases where $i = j$ and where $i \neq j$.
• I'm still not entirely certain how to deal with "scaling". But I think it should be possible to do this by just substituting the index on all occurring sums. For example, if I turn $\Sigma_i g_i \sigma_i$ into $\Sigma_j g_i\sigma_i$, then the latter is just $\Sigma_j g_i\sigma_i = N g_i \sigma_i$. The nice thing about this would be that this would just be naturally implemented as there's a short-hand when creating a sum that checks whether the summation index actually occurs in the expression you're summing over. If not, then you just return the range of the index times the expression, so N * x.

Some other notes, that can be addressed in a follow-up PR:

• Why does an Index need a Hilbert space? It's just a number, I don't understand and I think removing it has little effect. I essentially ignored the index Hilbert space throughout the implementation.
• I'm also not sure whether defining the range of an index should be done on the index itself. Usually, you define that on the object. At the very least, we should have an upper and lower bound on an index. Right now, you always assume it starts at 1 and there's no real reason for that.
• I'm not sure how far the concept of "doing math on indices" works (not far probably), but it would be interesting to see. It would be nice to e.g. be able to write something such as $\Sigma_i \sigma_i \sigma_{i+1}$ and have it just work.

@ChristophHotter anything you'd like to add here?

[ChristophHotter]

That sounds very good, nice!
I have a few comments:

1. Regarding the product of transition operators: If the product (i≠j)*(i=j) = 0 is defined one could maybe write (i≠j) s_i*s_j
2. I think the scaling idea is very nice. In order to account for the correct terms (e.g. N-1) a possible solution would probably be to set all indices on the LHS of the eqs to a specific value (e.g. j=1) but also the corresponding ones on the RHS. Then substitute the indices of the sums. At the end substitute the remaining indices (on the RHS).
3. For the function complete() the number of operators that need to be derived can be reduced if the equations get scaled afterward since there are more "redundant" terms - but I guess this is not so important at the moment.
4. I agree that the index probably does not need a HilbertSpace and the range should be defined later.

[david-pl]

Okay, completing equations now kinda works. However, simplification of sum terms is an issue. They currently don't get fully simplified which leads to "leftover" terms that should actually be 0. Therefore, we add some unnecessary equations when completing a set that has some sums. For documentation purposes, here are my thoughts:

• So long as we're dealing with QNumbers, we can use the same approach as always and bake some simplification into the methods for +,*,-. However, I'm not sure whether it's best for this purpose to write a Sum(::QAdd, index) as a QAdd containing Sums or not.
• The implementation for CNumber summation needs to be redone and properly cleaned up. Simplification should ideally be done with a set of specific @rules via SymbolicUtils. However, that can be tricky since you'd need to change indices at some points in order to check for equality, e.g. Sum(x[j], j] - Sum(x[i], i) should be 0 so long as i and j have the same range.

[david-pl]

Alright, significant update here: the following example works, and actually agrees with the results when doing this in a "brute-force" approach, where we create all the Hilbert spaces explicitly and derive all the equations.

using QuantumCumulants
using OrdinaryDiffEq, ModelingToolkit

# Hilbertspace
hc = FockSpace(:cavity)
ha = NLevelSpace(:atom,2)
h = hc ⊗ ha

# operators
@qnumbers a::Destroy(h)
σ(α,β,i) = IndexedOperator(Transition(h, :σ, α, β),i)


@cnumbers N Δ κ Γ R ν
g(i) = IndexedVariable(:g, i)

i = Index(h,:i,N,ha)
j = Index(h,:j,N,ha)
k = Index(h,:k,N,ha)

# Hamiltonian
H = -Δ*a'a + Σ(g(i)*( a'*σ(1,2,i) + a*σ(2,1,i) ),i)

# Jump operators with corresponding rates
J = [a, σ(1,2,i), σ(2,1,i), σ(2,2,i)]
rates = [κ, Γ, R, ν]

# Derive equations
ops = [a'*a, σ(2,2,j)]
eqs = meanfield(ops,H,J;rates=rates,order=2)

# custom filter function
φ(x::Average) = φ(x.arguments[1])
φ(::Destroy) = -1
φ(::Create) =1
φ(x::QTerm) = sum(map(φ, x.args_nc))
φ(x::Transition) = x.i - x.j
φ(x::IndexedOperator) = x.op.i - x.op.j
φ(x::Sum) = φ(x.term)
phase_invariant(x) = iszero(φ(x))

# Complete equations
eqs_c = complete(eqs; filter_func=phase_invariant)

N0 = 2
eqs_eval = evaluate(eqs_c; limits=Dict(N => N0))

@named sys = ODESystem(eqs_eval)

u0 = zeros(ComplexF64, length(eqs_eval))
g_ = g(i).arguments[1]  # TODO: this is awkward; need a clean way to define the Array that is g without any indices
ps = [
    Δ => 0.0;
    g_ => [0.5 for i=1:N0];
    κ => 1.0;
    Γ => 0.1;
    R => 0.9;
    ν => 0.0;
]

prob = ODEProblem(sys, u0, (0.0, 10.0), ps)
sol = solve(prob, Tsit5())

using PyPlot; pygui(true)
plot(sol.t, sol[a'*a])

[david-pl]

While the basic idea works, and I'm still in favor of using this approach, there are still some open topics. Especially towards the end, I got a bit sloppy and the code isn't super nice or easy to debug.

• One fundamental issue comes with the basic idea of using the term $ (i != j) \sigma_i \sigma_j$ when merging transitions and additionally storing this in a merge_events field: when deriving equations, find_missing will only pick up the average/operator, not the i != j factor. The resulting equations, however, actually need to be multiplied by i != j as they don't make sense when i == j. Similarly, this needs to be respected when e.g. using change_index in order to actually insert integers in the end. We can still infer that indices must not be equal by checking whether two transition operators have been merged previously, but I dislike this as it's not always clear when you need to check for that and debugging this is a pain. So we might need to come up with a cleaner way of doing this.
• Summation, especially when working with CNumbers is tricky. It kinda works right now though.
• Setting the shape for parameters with indexing is quite ugly, since we need to somehow set this in the metadata when inserting actual numbers for N. Otherwise, ModelingToolkit will complain. Setting the metadata is a bit awkward (when compared to just manipulating the expression tree) and the current implementation is quite ugly.
• We need a good way to easily distinguish whenever we are dealing with indexed expressions in meanfield. I'm not a fan of having a separate type for that, I think it can just be a boolean field on the normal MeanfieldEquations struct.

Things that still need to be done:

• Correlations, scaling, etc. Basically anything besides the standard time dynamics.
• Generating sum(x for i=1:N) instead of expanding them literally.
• Everything needs at least one more round of cleaning up.
• I did not care for performance at all so far.

I'll leave it at that for now, however.

@ChristophHotter I think now is the time to really have a look at the suggested implementation here. You can also use the above example to test things a bit. I left a whole bunch of TODO comments throughout, with some thoughts for possible issues or different approaches.

Also, basically all the old files are outdated. I didn't delete them, but they are not actually included anymore. The files you need to look at are:

• indexing.jl
• indexing_sums.jl
• indexing_utils.jl

And, of course, check the changes to other exisiting files such as meanfield.jl