Coupled Boltzmann computation of mixed axion neutralino dark matter in
  the SUSY DFSZ axion model
  The supersymmetrized DFSZ axion model is highly motivated not only because it
offers solutions to both the gauge hierarchy and strong CP problems, but also
because it provides a solution to the SUSY mu problem which naturally allows
for a Little Hierarchy. We compute the expected mixed axion-neutralino dark
matter abundance for the SUSY DFSZ axion model in two benchmark cases-- a
natural SUSY model with a standard neutralino underabundance (SUA) and an
mSUGRA/CMSSM model with a standard overabundance (SOA). Our computation
implements coupled Boltzmann equations which track the radiation density along
with neutralino, axion (produced thermally (TH) and via coherent oscillations
(CO)), saxion (TH- and CO-produced), axino and gravitino densities. In the SUSY
DFSZ model, axions, axinos and saxions go through the process of freeze-in-- in
contrast to freeze-out or out-of-equilibrium production as in the SUSY KSVZ
model-- resulting in thermal yields which are largely independent of the
re-heat temperature. We find the SUA case with suppressed saxion-axion
couplings (\xi=0) only admits solutions for PQ breaking scale f_a~< 5\times
10^{12} GeV where the bulk of parameter space tends to be axion-dominated. For
SUA with allowed saxion-axion couplings (\xi =1), then f_a values up to ~
2\times 10^{14} GeV are allowed. For the SOA case, almost all of SUSY DFSZ
parameter space is disallowed by a combination of overproduction of dark
matter, overproduction of dark radiation or violation of BBN constraints. An
exception occurs at very large f_a~ 10^{15}-10^{16} GeV where large entropy
dilution from CO-produced saxions leads to allowed models.
