I think this is an xarray question. If you do elementwise operations (like subtraction) between two xarray DataArrays, it will automatically broadcast across dimensions that share the same name. So something like this will do what you want:
phi_post = fit.posterior.phi # (chain, draw, n_circuits)
tau_post = fit.posterior.tau # (chain, draw, n_teams)
difference = phi_post - tau_post # (chain, draw, n_circuits, n_teams)
This is not a correlation matrix, it’s just a matrix of pairwise posterior distances between the parameters. Most obviously, there’s nothing that forces the differences to be between -1 and 1. If you want to estimation correlations between teams and courses, you should explicitly include it in your model. I don’t know anything about STAN, so I can’t help you there 