Within ZFC, the von Neumann hierarchy consists of sets $V_\alpha$ indexed by ordinals, subject to the following conditions:

• $V_0=\varnothing$.
• $V_{\alpha+1}=\mathcal P(V_\alpha)$.
• $V_\lambda=\bigcup_{\beta<\lambda}V_\beta$ for limit $\lambda$.

My question is: what is the formal justification for the last step?

The axiom of union would allow us to construct $V_\lambda$ if we already had a set $\{V_\beta:\beta<\lambda \}$. However, the existence of $V_\beta$ for $\beta<\lambda$ does not obviously guarantee the existence of this set: $V_{\omega\cdot 2}$ in ZF famously models itself minus the axiom of replacement. This of course suggests that with replacement, this set could be constructed as the image of $\lambda$ (which under the von Neumann ordinal construction is the set of all lower ordinals) under the class function $V$.

But then $V$ would have to be a class function somehow defined through transfinite induction. Since predicates can't refer to themselves, and since we can't just assert the existence and uniqueness of a class function as a theorem the way we can with normal functions, the way this works eludes me.