  We explore the phase space spanned by the temperature and the chemical
potential for 4-flavor lattice QCD using the Wilson-clover quark action. In
order to determine the order of the phase transition, we apply finite size
scaling analyses to gluonic and quark observables including plaquette, Polyakov
loop and quark number density, and examine their susceptibility, skewness,
kurtosis and Challa-Landau-Binder cumulant. Simulations were carried out on
lattices of a temporal size fixed at $N_{\text{t}}=4$ and spatial sizes chosen
from $6^3$ up to $10^3$. Configurations were generated using the phase
reweighting approach, while the value of the phase of the quark determinant
were carefully monitored. The $\mu$-parameter reweighting technique is employed
to precisely locate the point of the phase transition. Among various
approximation schemes for calculating the ratio of quark determinants needed
for $\mu$-reweighting, we found the Taylor expansion of the logarithm of the
quark determinant to be the most reliable. Our finite-size analyses show that
the transition is first order at $(\beta, \kappa, \mu/T)=(1.58, 0.1385,
0.584\pm 0.008)$ where $(m_\pi/m_\rho, T/m_\rho)=(0.822, 0.154)$. It weakens
considerably at $(\beta, \kappa, \mu/T)=(1.60, 0.1371, 0.821\pm 0.008)$ where
$(m_\pi/m_\rho, T/m_\rho)=(0.839, 0.150)$, and a crossover rather than a first
order phase transition cannot be ruled out.
