  The nature of the thermal phase transition of two flavor QCD in the chiral
limit has an important implication for the QCD phase diagram. We carry out
lattice QCD simulations in an attempt to address this problem. Simulations are
conducted with a Symanzik-improved gauge action and the HISQ fermion action.
Within the imaginary chemical potential formulation, five different quark
masses, $am=0.020,\, 0.018, \, 0.015, \, 0.013,\, 0.010$, and four different
lattice volumes $N_s=8, \, 12,\, 16, \, 20$ with temporal extent $N_t=4$ are
used to explore the scaling behavior. At each of the quark masses, the Binder
cumulants of the chiral condensate on different lattice volumes approximately
intersect at one point. We find that at the intersection point, the Binder
cumulant $B_4(am,a\mu_c) $ is around $3$ which deviates from the $Z(2)$
universality class value 1.604. However, based on the expectations of $Z(2)$
criticality, the fitting result only with the data from the largest lattice
volume $N_s=20$ agrees well with earlier result [ Phys. Rev., D90, 074030(2014)
]\cite{Bonati:2014kpa}. This fact implies that, although the finite cut-off
effects could be reduced with HISQ fermions even on $N_t=4$ lattices, larger
lattices with spatial extent $N_s>=20$ for such studies are needed to control
finite volume effects.
