Lagrangian data for velocity, scalars, and energy and scalar dissipation from direct numerical simulations are used to validate Lagrangian mixing models for inert passive scalars in stationary isotropic turbulence. The scalar fluctuations are nearly Gaussian, and, as a result of production by uniform mean gradients, statistically stationary. Comparisons are made for Taylor-scale Reynolds numbers in the range 38 to about 240 and Schmidt numbers in the range 1/8 to 1. Model predictions for one-point, one-time Eulerian statistics ~Eulerian correspondence! and one-particle, two-time Lagrangian statistics ~Lagrangian correspondence! are examined. Two scalar mixing models, namely the Lagrangian Fokker–Planck model and the Lagrangian colored-noise ~LCN! model, are proposed and written in terms of stochastic differential equations ~SDE! with specified drift and diffusion terms. Both of these models rely on statistics of the scalar field conditioned upon the energy dissipation, as provided by the Lagrangian spectral relaxation ~LSR! model. With the exception of the scalar dissipation, the models are shown to capture the Reynolds and Schmidt-number dependence of the Lagrangian integral time scales. However, the LCN model provides a more realistic description of the Lagrangian scalar fluctuations as differentiable time series having the correct form of the scalar autocorrelation function. Further extensions of the new mixing models to non-Gaussian scalars are conceptually straightforward, but require a closure for the scalar-conditioned scalar dissipation rate matrix. Likewise, accurate prediction of joint statistics for differential diffusion between different scalars with unequal molecular diffusivities will require the formulation of a multiscale SDE similar to the LSR model.