The Brightest Galaxies in the Dark Ages: Galaxies' Dust Continuum
  Emission During the Reionization Era
  Though half of cosmic starlight is absorbed by dust and reradiated at long
wavelengths (3$\mu$m-3mm), constraints on the infrared through millimeter
galaxy luminosity function (the `IRLF') are poor in comparison to the
rest-frame ultraviolet and optical galaxy luminosity function, particularly at
z>2.5. Here we present a backward evolution model for interpreting number
counts, redshift distributions, and cross-band flux density correlations in the
infrared and submillimeter sky, from 70$\mu$m-2mm, using a model for the IRLF
out to the epoch of reionization. Mock submillimeter maps are generated by
injecting sources according to the prescribed IRLF and flux densities drawn
from model spectral energy distributions that mirror the distribution of SEDs
observed in $0<z<5$ dusty star-forming galaxies (DSFGs). We explore two extreme
hypothetical case-studies: a dust-poor early Universe model, where DSFGs
contribute negligibly ($<$10%) to the integrated star-formation rate density at
$z>4$, and an alternate dust-rich early Universe model, where DSFGs dominate
$\sim$90% of $z>4$ star-formation. We find that current submm/mm datasets do
not clearly rule out either of these extreme models. We suggest that future
surveys at 2mm will be crucial to measuring the IRLF beyond $z\sim4$. The model
framework developed in this paper serves as a unique tool for the
interpretation of multiwavelength IR/submm extragalactic datasets and will
enable more refined constraints on the IRLF than can be made from direct
measurements of individual galaxies' integrated dust emission.
