  The Galactic disk is expected to be spatially, kinematically, and chemically
clustered on many scales due to both star formation and non-axisymmetries in
the Galactic potential. In this work we calculate the spatial and kinematic
two-point correlation functions using a sample of $1.7 \times 10^6$ stars
within 1 kpc of the Sun with 6D phase space information available from
\textit{Gaia} DR2. Clustering is detected on spatial scales of 1-300 pc and
velocity scales of at least 15 km s$^{-1}$. With bound structures included, the
data have a power-law index ($\xi(\Delta r) \propto \Delta r^{\gamma}$) of
$\gamma\approx-2$ at most spatial scales, which is in line with theoretical
predictions. After removing bound structures, the data have a power-law index
of $\gamma\approx-1$ for $<100$ pc and $\gamma \lesssim -1.5$ for $>100$ pc. We
interpret these results with the aid of a novel star-by-star simulation of the
Galaxy in which stars are born in clusters orbiting in a realistic potential
that includes spiral arms, a bar, and GMCs. We find that the simulation largely
agrees with the observations (within a factor of 2-3) at all spatial and
kinematic scales. In detail, the correlation function in the simulation is
shallower than the data at $< 20$ pc scales, and steeper than the data at $>
30$ pc scales. We also find a persistent clustering signal in the kinematic
correlation function for the data at large $\Delta v$ ($>5$ km s$^{-1}$) not
present in the simulations. We speculate that this mismatch between
observations and simulations may be due to two processes not included in the
simulation: hierarchical star formation and transient spiral arms. We also use
the simulations to predict the clustering signal as a function of pair-wise
metallicity and age separations. Ages and metallicities measured with a
precision of $50\%$ and $0.05$ dex are required in order to enhance the
clustering signal.
