  The mass spectra and wave functions of both $1S$ and $2S$ state heavy
pseudoscalar ($P$) and vector ($V$) mesons are analyzed within the light-front
quark model. Important empirical constraints employed in our analysis of the
mass spectra and wave functions are the experimental mass-gap relation, $\Delta
M_P > \Delta M_V$, where $\Delta M_{P(V)}=M^{2S}_{P(V)}-M^{1S}_{P(V)}$ and the
hierarchy of the decay constants, $f_{1S}>f_{2S}$, between $1S$ and $2S$ meson
states. We maintain the orthogonality of the trial wave functions of the $1S$
and $2S$ states in our variational calculation of the Hamiltonian with the
Coulomb plus confining potentials and treat the hyperfine interaction
perturbatively for the heavy-heavy and heavy-light $P$ and $V$ mesons due to
the nature of the heavy quark symmetry. Realizing that the empirical
constraints cannot be satisfied without mixing of the $1S$ and $2S$ states, we
find the lower bound of the mixing angle $\theta$ between $1S$ and $2S$ states
as $\theta_c = {\rm cot}^{-1}(2\sqrt{6})/2\simeq 6^{\circ}$ and obtain the
optimum value of the mixing angle around $12^\circ$ to cover both the charm and
bottom flavors of the heavy quark. The mixing effects are found to be more
significant to the $2S$ state mesons than to the $1S$ state mesons. The
properties of $1S$ and $2S$ state mesons including the mass spectra, decay
constants, twist-2 distribution amplitudes, and electromagnetic form factors
are computed. Our results are found to be in a good agreement with the
available data and lattice simulations. In particular, the $2S$ state
pseudoscalar $D_s$ meson is predicted to have a mass of $2600$ MeV, which is
very close to the mass of the newly discovered $D_{s0}(2590)^+$ meson by the
LHCb Collaboration. This supports the interpretation of the observed state as a
radial excitation of the $D^+_s$ meson.
