  We 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 $N_B^2$ imposed in the third, vertical
direction in a fluid with kinematic viscosity $\nu$. 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 $\mathscr{E}= \chi/(\chi+\mathcal{D})$ ($\mathcal{D}$ is the
volume-averaged kinetic energy dissipation rate and $\chi$ is the
volume-averaged density variance dissipation rate) varies non-monotonically
with the time-averaged buoyancy Reynolds numbers $\overline{Re}_B=
\overline{\mathcal{D}}/(\nu N_B^2)$, and is bounded above by $1/6$,
consistently with the classical model of Osborn (1980). There are qualitatively
different physical properties between the unstable orbits that have lower
irreversible mixing efficiency at low $\overline{Re}_B \sim O(1)$ and those
with nearly optimal $\mathscr{E} \lesssim 1/6$ at intermediate $\overline{Re}_B
\sim 10$. 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.
