The ceiling is a distribution, not the mean.
The true ceiling is p(y|x), not E[Y|X]. Given that you only observe x, the best any model can possibly do — in expected log-loss, against any proper scoring rule — is the mixture
p(y|x) = ∫ p(y|x, z) p(z|x) dz.
This is the Bayes-optimal predictive distribution for an x-only observer. You're right that no x-conditioned model can describe p(y|x, z) itself: z is marginalized out of your data, so the z-resolved conditional is simply not identifiable from (x, y) pairs. But the best achievable object is this whole mixture, with all its multimodality, not merely its mean. The mean y = f(x) = E[Y|X] is a further reduction that only becomes "the best model" after you additionally commit to point prediction under squared loss. Note by the tower property E[Y|X] = E_z[E[Y|X,Z]], which is exactly the quantity you wrote — so your statement is correct conditional on having chosen L2 point prediction, but that choice threw away the distribution first. Under a different loss the optimal point summary is a different functional of p(y|x) (median under L1, a quantile under pinball loss), and for generative tasks you need the mixture itself.
Your aleatoric intuition is exactly the law of total variance. Take the deterministic mechanism y = g(x, z). Then
Var(Y|X) = E_z[Var(Y|X,Z)] + Var_z(E[Y|X,Z] | X) = 0 + Var_z(g(x, Z) | X).
All the aleatoric spread is the variance of the z-resolved prediction as z fluctuates — precisely "the effect of unobservable sampling from z." (If the mechanism were itself stochastic even given z, the first term would add a residual floor.)
Epistemic is then the gap to the ceiling, at the distribution level. It's the discrepancy between your p̂(y|x) and the mixture p(y|x) — not the gap to E[Y|X,Z], which is a function of z you could never be scored against. In L2 regression this gap manifests as f̂(x) vs E[Y|X], so your formulation is the shadow the general statement casts onto point prediction. The clean decomposition of your expected error at x:
E[(Y − f̂(x))²|X=x] = Var(Y|X=x) + (f̂(x) − E[Y|X=x])²,
irreducible aleatoric floor plus squared epistemic gap. More data shrinks the second term toward zero; nothing shrinks the first — unless you change the conditioning set.
That last clause is the refinement worth internalizing: the aleatoric/epistemic boundary is relative to what you observe, not ontological. If z (or part of it) becomes measurable, spread that was "irreducible" migrates out of the aleatoric bucket — Var(Y|X,Z) < Var(Y|X). What's noise to an x-only observer is signal to an (x,z)-observer. So "aleatoric" should always be read as "aleatoric given this feature set," which is also why feature engineering and UQ are secretly the same conversation.
One final caveat on "no model can describe p(y|x,z)": models with latent variables (mixture density networks, conditional VAEs, diffusion) do posit an internal z and fit y = g(x, z_model). But they are only ever fit to match the mixture p(y|x); the latent they recover need not correspond to the true z, because infinitely many (g, p(z)) pairs induce the same mixture. They represent the right distribution through a possibly wrong mechanism — which is fine for prediction, and fatal if you wanted to interpret z causally.