The building up of observed stellar scaling relations of massive galaxies and the connection to black hole growth in the TNG50 simulation

[abridged] We study how mock-observed stellar morphological and structural properties of massive galaxies are built up between $z=0.5$ and $z=3$ in the TNG50 cosmological simulation. We generate mock images with the properties of the CANDELS survey and derive Sersic parameters and optical rest-frame morphologies as usually done in the observations. Overall, the simulation reproduces the observed evolution of the abundances of different galaxy morphological types of star-forming and quiescent galaxies. The $\log{M_*}-\log R_e$ and $\log{M_*}-\log\Sigma_1$ relations of the simulated star-forming and quenched galaxies also match the observed slopes and zeropoints to within 1-$\sigma$. In the simulation, galaxies increase their observed central stellar mass density ($\Sigma_1$) and transform in morphology from irregular/clumpy systems to normal Hubble-type systems in the Star Formation Main Sequence at a characteristic stellar mass of $\sim 10^{10.5}~M_\odot$. This morphological transformation is connected to the activity of the central Super Massive Black Holes (SMBHs). At low stellar masses ($10^9$ < $M_*/M_\odot$ < $10^{10}$) SMBHs grow rapidly, while at higher mass SMBHs switch into the kinetic feedback mode and grow more slowly. During this low-accretion phase, SMBH feedback leads to the quenching of star-formation, along with a simultaneous growth in $\Sigma_1$. More compact massive galaxies grow their SMBHs faster than extended ones of the same mass and end up quenching earlier. In the TNG50 simulation, SMBHs predominantly grow via gas accretion before galaxies quench, and $\Sigma_1$ increases substantially after SMBH growth slows down. The simulation predicts therefore that quiescent galaxies have higher $\Sigma_1$ values than star-forming galaxies for the same SMBH mass, which disagrees with alternative models, and may potentially be in tension with some observations.