How to distinguish the $p_x$, $p_y$, $p_z$ orbitals when one can choose the $x-$, $y-$, and $z-$ axes arbitrary? How to distinguish the $p_x$, $p_y$, $p_z$ orbitals since one could choose the $x-$, $y-$ and $z-$axes arbitrary?
How could the can the wavefunctions that describe the probability of the presence of electrons in the associated orbitals have a physical meaning? One could decide to make any arbitrary tiny 3D rotations so that the resulting orbitals would no more privilege some specific directions?
 A: It looks like you think that $p_x$, $p_y$ and $p_z$ have some special absolute meaning. They don't. A single physical state $|\psi\rangle$ could be described as $|p_x\rangle$ in one reference frame, as $|p_z\rangle$ in another, and as $\frac1{\sqrt2}(|p_x\rangle+|p_z\rangle)$ in yet another frame, with all these frames being related by rotations.
These states become special when external electric field is applied along a coordinate axis. Let $\vec E$ be aligned with $\vec{e}_x$. Then, since the electric field splits energy levels of these $p$ states (Stark effect), you'll be able to distinguish $|p_x\rangle$ from $|p_y\rangle$ or from $|p_z\rangle$ by noting the difference in e.g. frequencies of electromagnetic emissions.
A: The eigenstates of angular momentum, in particular $L^2$ and an arbitrary $L_z$, are the spherical harmonics:
$$ Y_l^m(\theta, \phi)$$
For fixed $l$ (the $L^2$ quantum number), the $2l + 1$ functions are eigenstates of $L_z$ with magnetic quantum number $-l \le m \le l$. They constitute a representation of the rotation group SO(3). That means they are closed under rotation.
If you have $p_z$ and rotate your coordinate system, it's now a mixture of $p_{z'}, p_{x'}, p_{y'}$, but there are no terms of the form $l\ne 1$. Of course, the energies are all degenerate, so without an external reference (say a magnetic or electric field), the choice of coordinates is completely arbitrary.
A: $\newcommand{\ket}[1]{\left|#1\right>}$
If you have two coordinate systems, $xyz$ and $x'y'z'$, you can construct the associated wave functions $\ket{p_x}$, $\ket{p_y}$, $\ket{p_z}$, and $\ket{p_{x'}}$, $\ket{p_{y'}}$, $\ket{p_{z'}}$. Then the latter set can be constructed as (generally complex) linear combinations of the former set.
So the arbitrary choice of coordinate system decides which orbitals you'll get to cover the full space of $\ell=1$ wave functions.
