  We present the first detailed analysis of the connection between galaxies and
their dark matter halos for the unWISE galaxy catalog -- a full-sky,
infrared-selected sample built from WISE data, containing over 500 million
galaxies. Using unWISE galaxy-galaxy auto-correlation and Planck CMB
lensing-galaxy cross-correlation measurements down to 10 arcmin angular scales,
we constrain the halo occupation distribution (HOD), a model describing how
central and satellite galaxies are distributed within dark matter halos, for
three unWISE} galaxy samples at mean redshifts $\bar{z} \approx 0.6$, $1.1$,
and $1.5$. We constrain the characteristic minimum halo mass to host a central
galaxy, $M_\mathrm{min}^\mathrm{HOD} = 1.83^{+0.41}_{-1.63} \times 10^{12}
M_\odot/h$, $5.22^{+0.34}_{-4.80} \times 10^{12} M_\odot/h$, $6.60
^{+0.30}_{-1.11} \times 10^{13} M_\odot/h$ for the unWISE samples at
$\bar{z}\approx 0.6$, $1.1$, and $1.5$, respectively. We find that all three
samples are dominated by central galaxies, rather than satellites. Using our
constrained HOD models, we infer the effective linear galaxy bias for each
unWISE sample, and find that it does not evolve as steeply with redshift as
found in previous perturbation-theory-based analyses of these galaxies. We
discuss possible sources of systematic uncertainty in our results, the most
significant of which is the uncertainty on the galaxy redshift distribution.
Our HOD constraints provide a detailed, quantitative understanding of how the
unWISE galaxies populate the underlying dark matter halo distribution. These
constraints will have a direct impact on future studies employing the unWISE
galaxies as a cosmological and astrophysical probe, including measurements of
ionized gas thermodynamics and dark matter profiles via Sunyaev-Zel'dovich and
lensing cross-correlations.
