Thanks it did!
I tried tt.stack and tt.concatenate but to no avail… The shape is so confusing to me!
I will look into pm.LKJCholeskyCov.
Edit: now looking back, I see it’s exactly the same expression in the data generating process y = np.stack((μ1, μ2)).T + ϵ, no idea why it never occurred to me to try it in the model…