ModelingToolkit complex system of equations

A few years ago I had a similar problem regarding a system of complex equations and I had found a workaround: ModelingToolkit matrix form system of equations and complex values.
As far as I can tell complex symbolics are still not supported in ModelingToolkit, but now unfortunately my workaround is no longer functioning either. Some of the ModelingToolkit syntax has changed in the almost 3 years since, so here is an updated MWE:

@parameters icomplex Γ32 δl δμ Ωl Ωμ
@variables t (ρ(t))[1:3,1:3]
D = Differential(t)

H = [
    -δl-δμ -Ωμ/2 0
    -Ωμ/2 -δl -Ωl/2
    0  -Ωl/2 0
]

L = [
    0 0 -Γ32*ρ[1,3]/2
    0 Γ32*ρ[3,3] -Γ32*ρ[2,3]/2
    -Γ32*ρ[3,1]/2 -Γ32*ρ[3,2]/2 -Γ32*ρ[3,3]
]

eqns = D(ρ) ~ -icomplex*(H*ρ-ρ*H) + L

@mtkbuild bloch = ODESystem(Symbolics.scalarize(eqns), t)

ρᵢ = zeros(ComplexF64,3,3)
ρᵢ[2,2] = 1
u0 = [ρ[idx,idy] => ρᵢ[idx,idy] for idx in 1:3 for idy in 1:3]

p = [Ωl => 1, Ωμ=> 0, δl => 0., δμ => 0., Γ32 => 1., icomplex => 1.0im]

prob = ODEProblem(bloch, u0, (0., 20.), p, jac = true)
sol = solve(prob, Tsit5())

The issue occurs when constructing the ODEProblem, throwing the an InexactError: Float64(0.0 + 1.0im). It seems like internally it tries to convert icomplex to a Float64.

Is there a different workaround possible?

@variables t (ρ(t))[1:3,1:3]::Complex? I haven’t tried complex arrays but complex variables should work if the type is appropriately declared.

I tried this but unfortunately it fails with the following error, which happens in the Symbolics.scalarize(eqns) step

TypeError: in Complex, in T, expected T<:Real, got Type{Complex{Real}}

It still isn’t functioning, but I’ve narrowed down the issue, but also found a new one.

The dissipator term L has explicit ρ terms, when these are removed from the equation the TypeError dissapears. Constructing the problem then throws the following error: ERROR: Initialization expression (ρ(t))[1, 1] is currently not supported.