  We study renormalization group multicritical fixed points in the
$\epsilon$-expansion of scalar field theories characterized by the symmetry of
the (hyper)cubic point group $H_N$. After reviewing the algebra of
$H_N$-invariant polynomials and arguing that there can be an entire family of
multicritical (hyper)cubic solutions with $\phi^{2n}$ interactions in
$d=\frac{2n}{n-1}-\epsilon$ dimensions, we use the general multicomponent beta
functionals formalism to study the special cases $d = 3-\epsilon$ and $d
=\frac{8}{3}-\epsilon$, deriving explicitly the beta functions describing the
flow of three- and four-critical (hyper)cubic models. We perform a study of
their fixed points, critical exponents and quadratic deformations for various
values of $N$, including the limit $N=0$, that was reported in another paper in
relation to the randomly diluted single-spin models, and an analysis of the
large $N$ limit, which turns out to be particularly interesting since it
depends on the specific multicriticality. We see that, in general, the
continuation in $N$ of the random solutions is different from the continuation
coming from large-$N$, and only the latter interpolates with the physically
interesting cases of low-$N$ such as $N=3$. Finally, we also include an
analysis of a theory with quintic interactions in $d =\frac{10}{3}-\epsilon$
and, for completeness, the NNLO computations in $d=4-\epsilon$.
