Context : Aims : We constrain the space density and properties of massive galaxy candidates at redshifts of z \geq 3.5 in the Great Observatories Origin Deep Survey North ( GOODS-N ) field . By selecting sources in the Spitzer + IRAC bands , a highly stellar mass-complete sample is assembled , including massive galaxies which are very faint in the optical/near-IR bands that would be missed by samples selected at shorter wavelengths . Methods : The z \geq 3.5 sample was selected down to m _ { AB } = 23 mag at 4.5 ~ { } \mu m using photometric redshifts that have been obtained by fitting the galaxies spectral energy distribution at optical , near-IR bands and IRAC bands . We also require that the brightest band ( in AB scale ) in which candidates are detected is the IRAC 8.0 ~ { } \mu m band in order to ensure that the near-IR 1.6 ~ { } \mu m ( rest-frame ) peak is falling in or beyond this band . Results : We found 53 z~ { } \geq 3.5 candidates , with masses in the range of M _ { \star } \sim 10 ^ { 10 } -10 ^ { 11 } M _ { \odot } . At least \sim 81 % of these galaxies are missed by traditional Lyman Break selection methods based on ultraviolet light . Spitzer + MIPS emission is detected for 60 % of the sample of z \geq 3.5 galaxy candidates . Although in some cases this might suggest a residual contamination from lower redshift star-forming galaxies or Active Galactic Nuclei , 37 % of these objects are also detected in the sub-mm/mm bands in recent SCUBA , AzTEC and MAMBO surveys , and have properties fully consistent with vigorous starburst galaxies at z \geq 3.5 . The comoving number density of galaxies with stellar masses \geq 5 \times 10 ^ { 10 } M _ { \odot } ( a reasonable stellar mass completeness limit for our sample ) is 2.6 \times 10 ^ { -5 } Mpc ^ { -3 } ( using the volume within 3.5 < z < 5 ) , and the corresponding stellar mass density is \sim ( 2.9 \pm 1.5 ) \times 10 ^ { 6 } M _ { \odot } Mpc ^ { -3 } , or about 3 % of the local density above the same stellar mass limit . For the sub-sample of MIPS-undetected galaxies , we find a number density of \sim 0.97 \times 10 ^ { -5 } Mpc ^ { -3 } and a stellar mass density of \sim ( 1.15 \pm 0.7 ) \times 10 ^ { 6 } M _ { \odot } Mpc ^ { -3 } . Even in the unlikely case that these are all truly quiescent galaxies , this would imply an increase in the space density of passive galaxies by a factor of \sim 15 from z \sim 4 to z = 2 , and by \sim 100 to z = 0 . Conclusions :