jats:pWe consider turbulence driven by a large-scale horizontal shear in Kolmogorov flow (i.e. with sinusoidal body forcing) and a background linear stable stratification with buoyancy frequency jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline1" />jats:tex-math$NB2$</jats:tex-math></jats:alternatives></jats:inline-formula> imposed in the third, vertical direction in a fluid with kinematic viscosity jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline2" />jats:tex-math$𝜈$</jats:tex-math></jats:alternatives></jats:inline-formula>. This flow is known to be organised into layers by nonlinear unstable steady states, which incline the background shear in the vertical and can be demonstrated to be the finite-amplitude saturation of a sequence of instabilities, originally from the laminar state. Here, we investigate the next order of motions in this system, i.e. the time-dependent mechanisms by which the density field is irreversibly mixed. This investigation is achieved using ‘recurrent flow analysis’. We identify (unstable) periodic orbits, which are embedded in the turbulent attractor, and use these orbits as proxies for the chaotic flow. We find that the time average of an appropriate measure of the ‘mixing efficiency’ of the flow jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline3" />jats:tex-math$E=𝜒/(𝜒+D)$</jats:tex-math></jats:alternatives></jats:inline-formula> (where jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline4" />jats:tex-math$D$</jats:tex-math></jats:alternatives></jats:inline-formula> is the volume-averaged kinetic energy dissipation rate and jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline5" />jats:tex-math$𝜒$</jats:tex-math></jats:alternatives></jats:inline-formula> is the volume-averaged density variance dissipation rate) varies non-monotonically with the time-averaged buoyancy Reynolds numbers jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline6" />jats:tex-math$Re―B=D―/(𝜈NB2)$</jats:tex-math></jats:alternatives></jats:inline-formula>, and is bounded above by jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline7" />jats:tex-math$1/6$</jats:tex-math></jats:alternatives></jats:inline-formula>, consistently with the classical model of Osborn (jats:italicJ. Phys. Oceanogr.</jats:italic>, vol. 10 (1), 1980, pp. 83–89). There are qualitatively different physical properties between the unstable orbits that have lower irreversible mixing efficiency at low jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline8" />jats:tex-math$Re―B\sim O(1)$</jats:tex-math></jats:alternatives></jats:inline-formula> and those with nearly optimal jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline9" />jats:tex-math$E\lesssim 1/6$</jats:tex-math></jats:alternatives></jats:inline-formula> at intermediate jats:inline-formulajats:alternatives<jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" mime-subtype="gif" xlink:type="simple" xlink:href="S0022112017006954_inline10" />jats:tex-math$Re―B\sim 10$</jats:tex-math></jats:alternatives></jats:inline-formula>. The weaker orbits, inevitably embedded in more strongly stratified flow, are characterised by straining or ‘scouring’ motions, while the more efficient orbits have clear overturning dynamics in more weakly stratified, and apparently shear-unstable flow.</jats:p>