Lattice Boltzmann models are simplifications of the continuum Boltzmann equation obtained by discretizing in both physical space and velocity space. The discrete velocities retained typically correspond to lattice vectors of the discrete spatial lattice. That is, each lattice vertex is linked to a finite number of neighboring vertices by lattice vectors. A particle distribution is therefore parametrized by its components in each of the discrete velocities, the vertex of the spatial lattice, and the discrete time. A time step of a classical lattice Boltzmann scheme then contains two steps: a relaxation step where distribution at each vertex is locally modified into a new distribution, and an advection step based on the method of characteristics as an exact time-integration operator. We show in a first part of our presentation that a single particle distribution for the D2Q13 lattice Boltzmann scheme can simulate coupled effects involving advection and diffusion of velocity and temperature. We consider various test cases: non-linear waves with periodic boundary conditions, a test case with buoyancy, propagation of transverse waves, Couette and Poiseuille flows. We test various boundary conditions and propose to mix bounce- back and anti-bounce-back numerical boundary conditions to take into account velocity and temperature Dirichlet conditions. We present also results for the de Vahl Davis heated cavity. Our results are compared with the coupled D2Q9-D2Q5 lattice Boltzmann approach for the Boussinesq system and with an elementary finite differences solver for the compressible Navier-Stokes equations. In a second time, we study the convergence of the scalar D2Q9 lattice Boltzmann scheme with multiple relaxation times when the time step is proportional to the space step and tends to zero. The classical formal analysis when all the relaxation parameters are fixed and the time step tends to zero shows that the numerical solution converges to solutions of the heat equation. If the diffusivity is fixed and the space step tends to zero, the relaxation parameter for the momentum is very small, causing a discrepency between the previous analysis and the numerical results. A new asymptotic partial differential equation, the damped acoustic system, is emergent as a result of this formal analysis. Complementary numerical experi- ments establish the convergence of the scalar D2Q9 lattice Boltzmann scheme with multiple relaxation times and acoustic scaling in this specific case of evanescent relaxation towards the numerical solution of the damped acoustic system.