Mixing of a passive scalar in statistically homogeneous, isotropic, and stationary turbulence with a mean scalar gradient is investigated via direct numerical simulation, for Taylor‐scale Reynolds numbers,R&lgr;, from 28 to 185. Multiple independent simulations are performed to get confidence intervals, and local regression smoothing is used to further reduce statistical fluctuations. The scalar fluctuation field, &fgr;(x,t), is initially zero, and develops to a statistically stationary state after about four eddy turnover times. Quantities investigated include the dissipation of scalar flux, which is found to be significant; probability density functions (pdfs) and joint‐pdfs of the scalar, its derivatives, scalar dissipation, and mechanical dissipation; and conditional expectations of scalar mixing, ∇2&fgr;. A linear model for scalar mixing jointly conditioned on the scalar andv‐velocity is developed, and reproduces the data quite well. Also considered is scalar mixing jointly conditioned on the scalar and scalar dissipation. Terms appearing in the balance equation for the pdf of &fgr; are examined. From a solution of the scalar pdf equation two sufficient conditions arise for the scalar pdf to be Gaussian. These are shown to be well satisfied for moderate values of the scalar, and approximately so for large fluctuations. Many correlations are also presented, including &rgr;(v,&fgr;), which changes during the evolution of the scalar from a value of unity when initialized to the stationary value of 0.5–0.6. ©1996 American Institute of Physics.