Massive galaxies at higher redshifts ($\emph{z}$ $>$ 2) show different characteristics from their local counterparts: They are compact and most likely have a disk. In this study, we trace the evolution of local massive galaxies by performing a detailed morphological analysis, namely, fitting single S\'{e}rsic profiles and performing bulge+disk decompositions. We analyze $\sim$ 250 massive galaxies selected from all CANDELS fields (COSMOS, UDS, EGS, GOODS-South and GOODS-North). We confirm that both star-forming and quiescent galaxies increase their sizes significantly from $\emph{z}$ $\approx$ 2.5 to the present day. The global S\'{e}rsic index of quiescent galaxies increases over time (from $n$ $\approx$ 2.5 to $n$ $>$ 4), while that of star-forming galaxies remains roughly constant ($n$ $\approx$ 2.5). By decomposing galaxy profiles into bulge+disk components, we find that massive galaxies at high redshift have prominent stellar disks, which are also evident from visual inspection of the images. By $z$ $\approx$ 0.5, the majority of the disks disappear and massive quiescent galaxies begin to resemble the local elliptical galaxies. Star-forming galaxies have lower bulge-to-total ratios ($B/T$) than their quiescent counterparts at each redshift bin. The bulges of star-forming and quiescent galaxies follow different evolutionary histories, while their disks evolve similarly. We conclude that major mergers, along with minor mergers, have played a crucial role in the significant size increase of high-\emph{z} galaxies and the destruction of their massive and large-scale disks.

A considerable fraction of the massive quiescent galaxies at \emph{z} $\approx$ 2, which are known to be much more compact than galaxies of comparable mass today, appear to have a disk. How well can we measure the bulge and disk properties of these systems? We simulate two-component model galaxies in order to systematically quantify the effects of non-homology in structures and the methods employed. We employ empirical scaling relations to produce realistic-looking local galaxies with a uniform and wide range of bulge-to-total ratios ($B/T$), and then rescale them to mimic the signal-to-noise ratios and sizes of observed galaxies at \emph{z} $\approx$ 2. This provides the most complete set of simulations to date for which we can examine the robustness of two-component decomposition of compact disk galaxies at different $B/T$. We confirm that the size of these massive, compact galaxies can be measured robustly using a single S\'{e}rsic fit. We can measure $B/T$ accurately without imposing any constraints on the light profile shape of the bulge, but, due to the small angular sizes of bulges at high redshift, their detailed properties can only be recovered for galaxies with $B/T$ \gax\ 0.2. The disk component, by contrast, can be measured with little difficulty.

Massive quiescent galaxies at $z \approx 2$ are apparently much more compact than galaxies of comparable mass today. How robust are these size measurements? We perform comprehensive simulations to determine possible biases and uncertainties in fitting single-component light distributions to real galaxies. In particular, we examine the robustness of the measurements of the luminosity, size, and other structural parameters. We devise simulations with increasing realism to systematically disentangle effects due to the technique (specifically using GALFIT) and the intrinsic structures of the galaxies. By accurately capturing the detailed substructures of nearby elliptical galaxies and then rescaling their sizes and signal-to-noise to mimic galaxies at different redshifts, we confirm that the massive quiescent galaxies at $z \approx 2$ are significantly more compact intrinsically than their local counterparts. Their observed compactness is not a result of missing faint outer light due to systematic errors in modeling. In fact, we find that fitting multi-component galaxies with a single S\'ersic profile, the procedure most commonly adopted in the literature, biases the inferred sizes higher by up to 10% - 20%, which accentuates the amount of size evolution required. If the sky estimation has been done robustly and the model for the point-spread function is fairly accurate, GALFIT can retrieve the properties of single-component galaxies over a wide range of signal-to-noise ratios without introducing any systematic errors.