  We study an interquark QQ^bar potential for the charmonium system, that is
determined from the the equal-time and Coulomb gauge QQ^bar Bethe-Salpeter (BS)
wavefunction through the effective Schr\"odinger equation. This novel approach
enables us to evaluate a kinetic heavy quark mass m_Q and a proper interquark
potential at finite quark mass m_Q, which receives all orders of 1/m_Q
corrections on the static QQ^bar potential from Wilson loops, simultaneously.
Precise information of the interquark potential for both charmonium and
bottomonium states directly from lattice QCD provides us a chance to improve
quark potential models, where the spin-independent interquark potential is
phenomenologically described by the Cornell potential and the spin-dependent
parts are deduced within the framework of perturbative QCD, from
first-principles calculations. In this study, calculations are carried out in
both quenched and dynamical fermion simulations. We first demonstrate that the
interquark potential at finite quark mass calculated by the BS amplitude method
smoothly approaches the conventional static heavy quark potential from Wilson
loops in the infinitely heavy quark limit within quenched lattice QCD
simulations. Secondly, we determine both spin-independent and -dependent parts
of the interquark potential for the charmonium system in 2+1 flavor dynamical
lattice QCD using the PACS-CS gauge configurations at the lightest pion mass,
M_\pi=156 MeV.
