Colm-cille Caulfield University of Cambridge We consider the nonlinear optimisation of the mixing of a passive scalar, initially arranged in two layers, in a two-dimensional plane Poiseuille flow at finite Reynolds and P\’eclet numbers, below the linear instability threshold. We use a nonlinear-adjoint looping approach to identify optimal perturbations leading to maximum time-averaged energy as well as maximum mixing in a freely evolving flow, measured through the minimisation of either the passive scalar variance or the so-called mix-norm, as defined by Mathew, Mezi\’c \& Petzold (Physica D, vol. 211, 2005, pp. 23–46). We show that energy optimisation appears to lead to very weak mixing of the scalar field whereas the optimal mixing initial perturbations, despite being less energetic, are able to homogenise the scalar field very effectively. For sufficiently long time horizons, minimising the mix-norm identifies optimal initial perturbations which are very similar to those which minimise scalar variance, demonstrating that minimisation of the mix-norm is an excellent proxy for effective mixing in this finite-Péclet-number bounded flow. By analysing the time evolution from initial perturbations of several optimal mixing solutions, we demonstrate that our optimisation method can identify the dominant underlying mixing mechanism, which appears to be classical Taylor dispersion, i.e. shear-augmented diffusion. The optimal mixing proceeds in three stages. First, the optimal mixing perturbation, energised through transient amplitude growth, transports the scalar field across the channel width. In a second stage, the mean flow shear acts to disperse the scalar distribution leading to enhanced diffusion. In a final third stage, linear relaxation diffusion is observed. We also demonstrate the usefulness of the developed variational framework in a more realistic control case: mixing optimisation by prescribed streamwise velocity boundary conditions.