From a molecular point of view the differences among phases lie in different average values of several quantities, for example distance between molecules, orientation, momentum. The kinds of quantities and the precise average values depend on the system under consideration (for example, orientation is very important in liquid crystals but unimportant in monatomic gases).
For a number of molecules of the order of the Avogadro constant, the bulk averages assume very precise values for statistical reasons, and are constant in time at equilibrium and also quite uniform when computed for slightly smaller parts of the system.
For a small number of molecules the "bulk" averages can have large fluctuations in time and across smaller subgroups of molecules, so the situation is not clear-cut. If we don't see any time trends in the time average of these "bulk" averages, we can use the time averages as indicators of a solid phase, liquid phase, etc. In a way we are actually extending the definition of "solid" etc in such situations. It makes indeed sense to indicate the spatial and temporal scale on which our "solid", "liquid", etc denomination applies.
A similar situation also exists for bulk matter when we consider the homogeneity of the averages between very small groups of molecules (spatial homogeneity), or the constancy for very short times. Even if we have a bulk system in, say, a liquid phase, we could find very small group of molecules for which the averages are, for some short time, more characteristic of the solid or gaseous phase, for example. This is indeed what happens in phase transitions: we can identify (in the extended sense above) different phases in different, small spatial regions; and the bulk is a mixture of those.
From the point of view of statistical mechanics each phase is characterized by a probability distribution of the relevant physical quantities, and this distribution can be made to depend on the number of particles. So we have to measure the frequency distribution of the quantities in our small system and compare it with the relevant probability distribution. The frequency distribution could be over time, or over space, or over repeated ex-novo preparations of the small system. We are basically doing an inverse-inference problem to infer the phase of the system.
This also leads to an answer to your question: even with only two molecules we can give a probability – once we specify a temporal scale and we give the measured values of the relevant quantities – that they are in a solid, liquid, gaseous, or other phase. As you increase the number of molecules the probability becomes more and more peaked on a specific phase. So there is not a clear threshold: you could say "I want to be sure at a 90% level" or "at a 60% level" depending on your application. It can also happen that the distribution remains peaked on more than one phase, indicating a phase transition; see the references below on this possibility.
I'd recommend to read the following work, because it studies this problem in detail and experimentally:
Also this:
And this review for the statistical mechanics of small systems in general: