MvNormal much slower in pymc5 than pymc3?

Almost certainly! For future reference, here is the code that runs extremely slowly:

import numpy as np
import matplotlib.pyplot as plt
import arviz as az
import pytensor.tensor as pt
nt = 100
n = 30

x = np.arange(nt)
#arbitrary good-looking function
true_mu=(50*x**2-.5*x**3-100*x)/1.e4
innovations = np.random.normal(scale=10, size=(n, nt))
true_rho = np.random.uniform()
ensemble = np.zeros((n, nt))
for t in range(1, nt):
    ensemble[:, t] = true_rho * ensemble[:, t-1] + innovations[:, t]+true_mu[t]
with pm.Model() as from_pymc3:

    σ=pm.HalfNormal("σ",10)
    ϕ=pm.Uniform("ϕ",-1,1)
   
    
    #AR(1) covariance matrix
    δ=np.arange(nt)
    imj=np.abs(np.subtract.outer(δ,δ))
    cov=pt.pow(σ,2)/(1-pt.pow(ϕ,2))*pt.pow(ϕ,imj)

    chol=pt.slinalg.cholesky(cov)
    
  
    #Uncertain mean value across all ensemble members
    μ=pm.Normal("μ",0,10,shape=nt)
    
    F=pm.MvNormal("likelihood_hist",mu=μ,chol=chol,\
                             observed= ensemble)
    pymc3_trace=pm.sample()