  We propose a minimal extension of the Standard Model by incorporating sterile
neutrinos and a QCD axion to account for the mass and mixing hierarchies of
quarks and leptons and to solve the strong CP problem, and by introducing
$G_{\rm SM}\times \Gamma_N\times U(1)_X$ symmetry. We demonstrate that the
K{\"a}hler transformation corrects the weight of modular forms in the
superpotential and show that the model is consistent with the modular and
$U(1)_X$ anomaly-free conditions. This enables a simple construction of a
modular-independent superpotential for scalar potential. Using minimal
supermultiplets, we demonstrate a level 3 modular form-induced superpotential.
Sterile neutrinos explain small active neutrino masses via the seesaw mechanism
and provide a well-motivated $U(1)_X$ breaking scale, whereas gauge singlet
scalar fields play crucial roles in generating the QCD axion, heavy neutrino
mass, and fermion mass hierarchy. The model predicts a range for the $U(1)_X$
breaking scale from $10^{13}$ GeV to $10^{15}$ GeV for $1\,\mbox{TeV}<
m_{3/2}<10^6\,\mbox{TeV}$. In the supersymmetric limit, all Yukawa coefficients
in the superpotential are given by complex numbers with an absolute value of
unity, implying a democratic distribution. Performing numerical analysis, we
study how model parameters are constrained by current experimental results. In
particular, the model predicts that the value of the quark Dirac CP phase falls
between $38^\circ$ to $87^\circ$, which is consistent with experimental data,
and the favored value of the neutrino Dirac CP phase is around $250^\circ$.
Furthermore, the model can be tested by ongoing and future experiments on axion
searches, neutrino oscillations, and $0\nu\beta\beta$-decay.
