Toward 1966-1967, M. S. Baouendi, [1], has considered the boundary value problem and hypoellipticity on essentially the operators P — Dl + yDx, k=l, 2, -. Shortly later, V. V. Grushin introduced in the papers [8] and [10] a wide class of degenerate elliptic operators including the above operator P with three conditions for them to be hypoelliptic. After then there have been investigated the problem of analytic and non-analytic hypoellipticity of the Grushin operators, (cf. [19], [2], [9], [27], [28], [30], [31], [11]). The aim of this paper is to give a nearly complete answer to this problem by expanding the idea developed in [10] to utilize the operator-valued pseudodifferential operators. Such a method might be called a kind of separation of variables in x and y. The Gevrey index of each Grushin operator may be determined depending on the quasihomogeneity condition of the symbol, (cf. Condition 1.1, Condition 6.1 and Condition 7.1 and Theorem 1.2, Theorem 6.1 and Theorem 7.1 respectively). We treat the Grushin operators dividing into three groups. The operators in the first group introduced in § 1 including the above operator P are analytic hypoelliptic in the space of hyperfunctions % as mentioned in Theorem 1.2. The operators in the second and the third group are Gevrey hypoelliptic in the corresponding ultradistribution spaces ®, (0> 1) and treated similarly as those of the first group, (cf. Theorem 6.1 and Theorem 7.1). In § 2~§ 4, we prepare necessary steps to prove these three theorems. In § 2, it will be shown that any eigenfunction of Grushin operators of the form (2.1) satisfying ellipticity condition (2.2) belongs to the jJ-space of Gel'fand-Shilov, (in fact ^i/d+^). This fact (especially l=^+j^) plays an