  In this article, we review the HAL QCD method to investigate baryon-baryon
interactions such as nuclear forces in lattice QCD. We first explain our
strategy in detail to investigate baryon-baryon interactions by defining
potentials in field theories such as QCD. We introduce the Nambu-Bethe-Salpeter
(NBS) wave functions in QCD for two baryons below the inelastic threshold. We
then define the potential from NBS wave functions in terms of the derivative
expansion, which is shown to reproduce the scattering phase shifts correctly
below the inelastic threshold. Using this definition, we formulate a method to
extract the potential in lattice QCD. Secondly, we discuss pros and cons of the
HAL QCD method, by comparing it with the conventional method, where one
directly extracts the scattering phase shifts from the finite volume energies
through the L\"uscher's formula. We give several theoretical and numerical
evidences that the conventional method combined with the naive plateau fitting
for the finite volume energies in the literature so far fails to work on
baryon-baryon interactions due to contaminations of elastic excited states. On
the other hand, we show that such a serious problem can be avoided in the HAL
QCD method by defining the potential in an energy-independent way. We also
discuss systematics of the HAL QCD method, in particular errors associated with
a truncation of the derivative expansion. Thirdly, we present several results
obtained from the HAL QCD method, which include (central) nuclear force, tensor
force, spin-orbital force, and three nucleon force. We finally show the latest
results calculated at the nearly physical pion mass, $m_\pi \simeq 146$ MeV,
including hyperon forces which lead to form $\Omega\Omega$ and $N\Omega$
dibaryons.
