1. Introduction
Turbulent flows in straight pipes of various cross-sections are of utmost relevance for many industrial processes. For circular pipes, the turbulent dynamics both in the transitional and the fully developed regime is nowadays quite well understood (Smits, McKeon & Marusic Reference Smits, McKeon and Marusic2011; Avila, Barkley & Hof Reference Avila, Barkley and Hof2023). Ducts with a rectangular cross-section, on the other hand, have attained much less attention, despite their relevance in many technical applications such as air ventilation systems. In contrast to the circular case, the rotational symmetry with respect to the centreline is broken by the four surrounding walls, and a two-dimensional mean velocity field with a distinct mean secondary flow of Prandtl’s second kind is the consequence (Bradshaw Reference Bradshaw1987). Although of comparably low amplitude (usually approximately
$\mathcal{O}(1\,\%)$
of the bulk velocity
$u_b$
), these mean secondary currents are of significance for the mass, momentum and heat transfer between the near-wall zones, the corner regions and the duct core (Demuren & Rodi Reference Demuren and Rodi1984; Modesti & Pirozzoli Reference Modesti and Pirozzoli2022). In a streamwise- and time-averaged sense, the turbulence-induced secondary flow is generated by the inhomogeneity and anisotropy of the Reynolds-stress tensor across the duct cross-section (Speziale Reference Speziale1982). Hence, most previous experimental studies (Brundrett & Baines Reference Brundrett and Baines1964; Gessner Reference Gessner1973) and direct numerical simulations (DNS) (Gavrilakis Reference Gavrilakis1992, Reference Gavrilakis2019; Modesti et al. Reference Modesti, Pirozzoli, Orlandi and Grasso2018; Pirozzoli et al. Reference Pirozzoli, Modesti, Orlandi and Grasso2018) laid great emphasis on a detailed analysis of the mean momentum, vorticity and energy budgets (Nikora & Roy Reference Nikora and Roy2012).
While such investigations of the mean flow equations are insightful, for instance, regarding the scaling of the mean secondary flow pattern and its intensity with the Reynolds number, they provide little information about the underlying physical mechanisms and the involved coherent flow structures. Uhlmann et al. (Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007) and Pinelli et al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010) therefore performed detailed coherent structure eduction studies for a range of marginal to moderate Reynolds numbers, assessing the role of individual streaks and quasi-streamwise vortices in the generation of the characteristic mean secondary flow pattern. Uhlmann et al. (Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007) showed that, for sufficiently high Reynolds numbers
${\textit{Re}}_b=u_b H\!/\!\nu$
(where
$H$
is the duct’s half-side length and
$\nu$
is the kinematic viscosity), each of the four bounding walls accommodates one or more streaks with flanking quasi-streamwise vortices that undergo the classical self-sustained near-wall regeneration cycle (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997; Hall & Sherwin Reference Hall and Sherwin2010). However, the quasi-streamwise vortices next to the corner regions are limited in their spatial mobility by the two adjacent walls (Kawahara & Kamada Reference Kawahara and Kamada2000), so that they get essentially locked in their cross-sectional position and leave a clear footprint in the mean streamwise vorticity
${\langle {\omega _{x}} \rangle _{xt}}$
. The resulting preferential organisation of small-scale vortices in certain regions of the cross-section has direct implications for the large-scale secondary flow that is coupled to the small-scale dynamics by a Poisson equation
${\nabla} ^2\langle {\psi } \rangle _{xt}=-\langle {\omega _{x}} \rangle _{xt}$
; with
$\langle {\psi } \rangle _{xt}$
indicating the mean secondary flow streamfunction (Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010). Further theoretical support for this causal relation between coherent structures and the mean secondary flow was provided by Uhlmann, Kawahara & Pinelli (Reference Uhlmann, Kawahara and Pinelli2010), who found a family of travelling wave solutions whose ‘exact coherent structures’ induce a strikingly similar ‘eight-vortex’ mean secondary flow pattern (see also the discussion of invariant solutions below).
At marginal Reynolds numbers
${\textit{Re}}_b$
just large enough to sustain turbulence, on the other hand, the typical diameter of the quasi-streamwise vortices (which naturally scale in wall units) is, at
$\mathcal{O}(0.1 H)$
, so large that independent streak–vortex pairs cannot always be sustained on all four walls simultaneously (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007). Instead, only a subset of the four surrounding walls features a significant turbulent/vortical activity during such episodes, while turbulent fluctuations are weak or even entirely absent near the remaining ones. In the remainder of this study, we will occasionally refer to these two situations as ‘active’ and ‘non-active’ walls, respectively. Figure 1(a,b) present the streamwise-averaged velocity fields for two such flow states extracted from a long-time turbulent trajectory at
${\textit{Re}}_b=1150$
, in which a single and two opposing walls are ‘active’ in the above defined sense. Such localised states are transient in time and a ‘switching’ of the vortical activity to other walls takes place intermittently in intervals of
$\mathcal{O}(100)$
bulk time unit length (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007), although the underlying physical mechanisms are not entirely understood. Only long-time averages over intervals much longer than the lifetime of the individual coherent structures and the characteristic ‘wall-switching’ time scale do recover the characteristic fourfold-symmetric eight-vortex mean secondary flow pattern under such marginal conditions (cf. figure 1
c). Both the ‘wall-switching’ behaviour and the localisation of turbulent fluctuations to individual walls possess analogues in ‘minimal flow units’ of plane channel flow (Jiménez & Moin Reference Jiménez and Moin1991): for extended time periods, only one of the two channel walls hosts a streak with accompanying vortices, while the other one features an essentially laminar dynamics until an intermittent switching event causes a reorganisation of the flow (Neelavara, Duguet & Lusseyran Reference Neelavara, Duguet and Lusseyran2017). Also, there is an equivalence between the discussed temporally intermittent ‘wall-switching’ behaviour and a spatially intermittent ‘wall switching’ observed in DNS of extended domains like the
$L_x=160\pi H$
long duct studied by Sekimoto (Reference Sekimoto2011). For a fixed moment in time, the dynamics along each wall is here characterised by a streamwise alternating pattern of ‘active’ and ‘passive’ regions, all together inducing the characteristic fourfold-symmetric ‘eight-vortex’ mean secondary flow in the instantaneous streamwise average.

Figure 1. (a,b) Selected instantaneous streamwise-averaged velocity fields
$\langle {{\boldsymbol{u}}} \rangle _{x}$
taken from a long-time turbulent trajectory of
$\mathcal{O}(10^4)$
bulk time units length at a bulk Reynolds number
${\textit{Re}}_b=1150$
and a streamwise domain length
$L_x=2\pi \!H$
. (c) Long-time average of the velocity field
$\langle {{\boldsymbol{u}}} \rangle _{xt}$
over the full trajectory. In all panels, solid lines represent isolines of
$\langle {u} \rangle$
between
$0$
and
$1.5$
times the bulk velocity
$u_b$
, with increment
$0.25$
. Intensity and orientation of the secondary flow field
$(\langle {v} \rangle ,\langle {w} \rangle )^T$
is indicated by the vector plot. The arrows in (c) are scaled by a factor of two compared with those in (a,b) for better visibility. For a complete definition of the velocity field and the averaging operators, see § 2.
While the occurrence of spatially localised states and intermittent wall-switching events in low-Reynolds-number square duct turbulence is hence well documented, our understanding of the underlying physical mechanisms remains incomplete. In this study, we analyse both phenomena from a dynamical systems’ perspective of fluid turbulence, in which the time evolution of a turbulent system is understood as a trajectory wandering across the state space associated with the dissipative Navier–Stokes system (Hopf Reference Hopf1948). While the latter is formally of infinite dimension, the turbulent dynamics within a finite spatial domain can be expected to lie on an inertial manifold, and can therefore be expressed with a finite number of degrees of freedom (Foias et al. Reference Foias, Manley, Rosa and Temam2001). Turbulent coherent structures are interpreted as physical incarnations of the least unstable invariant solutions in this state space (Jiménez Reference Jiménez1987; Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2008), which comprises stationary equilibria, travelling waves (TWs, equilibria within a co-moving frame of reference) and periodic orbits (POs). Together with their stable and unstable manifolds, these invariant solutions form the skeleton of the system’s dynamics in state space: chaotic turbulent trajectories visit their neighbourhood transiently and shadow their dynamics during that time, before they get repelled along their unstable manifolds (Waleffe Reference Waleffe2002). In the past decades, this conceptual picture has been extraordinarily fruitful for the understanding of transitional (Kerswell Reference Kerswell2005; Avila et al. Reference Avila, Barkley and Hof2023) and sustained turbulence (Kawahara, Uhlmann & van Veen Reference Kawahara, Uhlmann and van Veen2012; Graham & Floryan Reference Graham and Floryan2021) in a variety of different flow configurations. Most prominently, equilibria, TWs and (relative) POs have been found in plane Couette (Nagata Reference Nagata1990; Kawahara & Kida Reference Kawahara and Kida2001; Gibson et al. Reference Gibson, Halcrow and Cvitanović2008) and plane Poiseuille flow (Park, Shekar & Graham Reference Park, Shekar and Graham2018; Waleffe Reference Waleffe2001, Reference Waleffe2003), asymptotic suction boundary layers (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013) and circular pipe flow (Duguet, Willis & Kerswell Reference Duguet, Willis and Kerswell2008; Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013; Budanur et al. Reference Budanur, Short, Farazmand, Willis and Cvitanović2017). In square duct flows, on the other hand, only a handful of TWs are known to date: the solution documented by Uhlmann et al. (Reference Uhlmann, Kawahara and Pinelli2010) generates an ‘eight-vortex’ mean secondary flow pattern that bears a striking similarity to the long-time average in low-Reynolds-number flows. The ‘four-vortex’ states of Wedin, Bottaro & Nagata (Reference Wedin, Bottaro and Nagata2009) and Okino et al. (Reference Okino, Nagata, Wedin and Bottaro2010) and the ‘two-vortex’ state of Okino & Nagata (Reference Okino and Nagata2012), in turn, are representative of situations in which not all of the four surrounding walls feature a significant vortical activity. In contrast to canonical flows such as circular pipe or plane channel flows, however, no POs have been found so far in the square duct setting to the best of the authors’ knowledge.
For all above discussed flow configurations including square duct flow (Tatsumi & Yoshimura Reference Tatsumi and Yoshimura1990), the laminar flow state is stable with respect to infinitesimal disturbances for Reynolds numbers at which transition to turbulence usually sets in. Hence, finite-amplitude perturbations are required to trigger the transition to turbulence and the laminar flow takes the form of an equilibrium/fixed point in state space, with a distinct basin of attraction (Graham & Floryan Reference Graham and Floryan2021). For transiently turbulent systems, the laminar state is complemented by a chaotic saddle to which nearby trajectories are attracted for a potentially long, but finite time horizon, before they finally experience a relaminarisation. For sustained turbulence, two concurring state space pictures have been discussed: turbulence might be further understood as a chaotic saddle, but with survival times of turbulent perturbations tending to infinity. Alternatively, turbulence is seen as a strange attractor from which trajectories cannot escape (cf. Avila et al. Reference Avila, Barkley and Hof2023 for a recent review).
Irrespective of the chosen viewpoint, the state space of all above considered flows has been seen to feature an invariant manifold of co-dimension one (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983) that is conventionally termed the ‘edge’
${\mathscr{M}}$
: trajectories starting off from a point inside the edge will stay in the latter for all times, neither swinging up to a turbulent state nor decaying directly towards the laminar fixed point. While the existence of the edge is consistent with both discussed state space pictures of turbulence, only in the presence of a turbulent attractor does it become a separatrix of phase space in a strict sense (Avila et al. Reference Avila, Barkley and Hof2023). Analysing the dynamics of such edge trajectories and identifying attracting sets within
${\mathscr{M}}$
(‘edge states’) has provided valuable insights into the transition processes to turbulence (Toh & Itano Reference Toh and Itano2003; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007; Duguet, Schlatter & Henningson Reference Duguet, Schlatter and Henningson2009; Mellibovsky et al. Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009). In this context, the nature of the detected edge states differs significantly among the different flows, ranging from simple equilibria in Couette flow (Schneider et al. Reference Schneider, Gibson, Lagha, De Lillo and Eckhardt2008) and POs in plane channels (Rawat, Cossu & Rincon Reference Rawat, Cossu and Rincon2014; Zammert & Eckhardt Reference Zammert and Eckhardt2014), asymptotic suction boundary layers (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016; Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013) or symmetric subspaces of circular pipe flow (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013) to chaotic edge states in the unconstrained circular pipe (Duguet et al. Reference Duguet, Willis and Kerswell2008).
To the best of the authors’ knowledge, only two groups have made efforts to analyse the edge for square duct flows: Biau & Bottaro (Reference Biau and Bottaro2009) compared a rather short edge trajectory with a flow state that they had found to be optimal in terms of the maximum nonlinear growth, without further discussing the edge dynamics or possible edge states. Okino (Reference Okino2014, in Japanese) and Gepner, Okino & Kawahara (Reference Gepner, Okino and Kawahara2025) searched for edge states in a square duct under certain symmetry restrictions and in domains of streamwise periods
$L_x/H\in \{2\pi ,4\pi ,8\pi \}$
. For the three considered domains, a ‘domain-filling’ TW, a chaotic (relative) attractor and a streamwise-localised TW were found to be edge states of the respective symmetry-reduced systems, respectively. Especially the streamwise-localised solution bears a strong resemblance to localised ‘puff-like’ solutions in the circular pipe (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013), with important implications for the understanding of localised turbulence in rectangular ducts (Takeishi et al. Reference Takeishi, Kawahara, Wakabayashi, Uhlmann and Pinelli2015).
However, a comprehensive study addressing the edge dynamics and how it compares with the characteristics of chaotic low-Reynolds-number turbulence is lacking. In this study, we therefore present a detailed analysis of the edge in square duct flow for both the full state space and symmetric subspaces thereof. We provide evidence that the edge states in the full state space are chaotic attractors within
${\mathscr{M}}$
, whereas their counterparts in the symmetric subspaces turn out to be POs. To the best of the authors’ knowledge, these periodic solutions are the first of their kind and are hence considered to be of relevance for the understanding of transitional and marginally turbulent duct flows. Specific implications for the understanding of transverse turbulence localisation and the process behind the ‘wall-switching’ dynamics are discussed.
The manuscript is organised as follows: in § 2, the considered physical system and the admissible symmetry groups are presented. The essential features of the DNS code and the details of the edge-tracking algorithm are described in § 3. With the aid of these methods, the edge to turbulence in square duct flows is analysed in § 4: the edge dynamics in the full state space is addressed first, before we turn to the edge states in symmetric subspaces. In § 5, the detected edge states are juxtaposed to states in marginally turbulent square duct flow and their properties are discussed in comparison with edge states in other flows. The manuscript ends with a summary of the relevant outcomes and an outlook on future work in § 6.
2. Theoretical framework
2.1. Flow configuration
In the present study, the flow of an incompressible Newtonian fluid in a straight, streamwise-periodic duct of square cross-section with side lengths
$2{{H}}$
is considered. No-slip boundary conditions are imposed along the four surrounding walls. A Cartesian coordinate system with the origin located in the centre of the cross-section is adopted, with basis vectors pointing in the streamwise (
$x$
) and the two cross-stream, wall-normal directions (
$y,z$
). For the sake of clarity, we choose a terminology similar to that in plane Couette or Poiseuille flow and refer to
$y$
and
$z$
as the ‘vertical’ and ‘spanwise’ directions in the remainder of this work. Note that this designation is an arbitrary choice since the two cross-sectional directions are equivalent regarding the rotational symmetry obeyed by the governing equations. In what follows, the velocity field is expressed as
${\boldsymbol{u}}({\boldsymbol{x}},t)=(u,v,w)^T$
, with
$u$
,
$v$
and
$w$
indicating the fluid velocity along the
$x$
-,
$y$
- and
$z$
-directions, respectively. Analogously, the vorticity field is introduced as
${\boldsymbol{\omega }}({\boldsymbol{x}},t)\equiv \boldsymbol{\nabla }\times {\boldsymbol{u}}=({\omega _{x}},{\omega _{y}},{\omega _{z}})^T$
.
All flow variables are periodic in the streamwise direction, with fundamental period
${L_x}=2\pi\! {{H}}/{\alpha }$
and fundamental wavenumber
$\alpha$
. If not otherwise stated,
${\alpha }=1$
and thus
${L_x}=2\pi\! {{H}}$
holds for the simulation results presented in this manuscript. Streamwise, volume and time averaging operators are consequently defined as
\begin{equation} {\langle {\bullet } \rangle _{x}} \equiv \dfrac {{\alpha }}{2\pi }\displaystyle \int \limits _0^{2\pi /{\alpha }} \bullet \;{\textrm{d}}x, \quad {\langle {\bullet } \rangle _{V}} \equiv \dfrac {1}{{\vert {{\varOmega }} \vert }}\displaystyle \int \limits _{\varOmega } \bullet \;{\textrm{d}}{\boldsymbol{x}}, \quad {\langle {\bullet } \rangle _{t}} \equiv \displaystyle \lim _{T\to \infty }\dfrac {1}{T}\displaystyle \int \limits _0^T \bullet \;{\textrm{d}}t, \end{equation}
where
${\varOmega }=[0,2\pi /{\alpha })\times [-{{H}},{{H}}]\times [-{{H}},{{H}}]$
denotes the physical domain. In a similar fashion, the
$L_2$
-norm of a vector field
$\boldsymbol{a}$
is introduced as
\begin{equation} {\left \lVert {{\boldsymbol{a}}} \right \rVert} _{2}^2 \equiv \dfrac {1}{{\vert {{\varOmega }} \vert }}\displaystyle \int \limits _{\varOmega } {\boldsymbol{a}}\boldsymbol{\cdot }{\boldsymbol{a}} \;{\textrm{d}}{\boldsymbol{x}} = {\langle {{\boldsymbol{a}}\boldsymbol{\cdot }{\boldsymbol{a}}} \rangle _{V}}. \end{equation}
A time-dependent streamwise pressure gradient
${\boldsymbol{f}_{\!e}}=({f_e}(t),0,0)^T$
is imposed and adjusted at each time step to drive the flow at a constant mass flow rate
$Q_m$
, guaranteeing a fixed bulk velocity
${u_b}\equiv Q_m/(4{{H}}^2)$
. Unless otherwise stated, time will be measured in bulk time units
${T_b}={{H}}/{u_b}$
in the remainder of this work. Both the perimeter-averaged wall shear stress
\begin{equation} {\overline {\tau }_w}(t)=\dfrac {\rho \nu }{8{{H}}} \left (\, \displaystyle \int \limits _{z=-{{H}}}^{{{H}}} \left .\dfrac {\partial u}{\partial y}\right |_{y\in \{-{{H}},{{H}}\}} {\textrm{d}}z + \displaystyle \int \limits _{y=-{{H}}}^{{{H}}} \left .\dfrac {\partial u}{\partial z}\right |_{z\in \{-{{H}},{{H}}\}} {\textrm{d}}y \,\right ) \end{equation}
and the friction velocity
${u_\tau }(t)=\sqrt {{\overline {\tau }_w}(t)/{\rho }}$
fluctuate in time. Variables non-dimensionalised with the time-averaged friction velocity
$\langle {{u_\tau }} \rangle _{t}$
and
$\nu$
are indicated as
${\bullet }^+$
. Based on this definition, all considered systems feature a domain clearly longer than the critical value required to sustain turbulence (
${L}_x^+|_{\textit{min}}\approx 200$
, cf. Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007 and compare with table 1) and are thus not minimal with respect to the buffer layer structures (Jiménez & Moin Reference Jiménez and Moin1991).
Table 1. Physical properties of the detected edge states. Here,
${\textit{Re}}_b$
is the bulk,
${\textit{Re}}_\tau$
the mean friction Reynolds number and
${L_x}/\!{{H}}$
,
${L}_x^+$
indicate the outer- and inner-scaled streamwise domain periods, respectively. For time-periodic edge states, the period (here given in bulk time units
${T_b}={{H}}/{u_b}$
) is estimated as described in § 4.2 and Appendix A. The last two columns list imposed symmetries and provide information on the shape of the found edge state, respectively (CA: chaotic attractor, PO: stable (relative) PO within
${\mathscr{M}}$
).

For the considered case of a pressure-driven square duct flow, the overall kinetic energy budget simplifies to a balance between the energy input
${I}$
and the dissipation rate
${D}$
per unit mass, viz.
where the individual contributions are defined as
The laminar base flow
${\boldsymbol{u}}_{\textit{lam}}= (u_{\textit{lam}}(y,z),0,0 )^T$
is defined as the solution of the following Poisson equation:
where
${{\nabla}} _{\perp }^2=(\partial ^2/\partial y^2 + \partial ^2/\partial z^2)$
represents a two-dimensional Laplacian acting along the two cross-sectional directions only, and
$f_{e,lam}$
indicates the constant pressure gradient necessary to drive a laminar flow at identical mass flow rate
$Q_m$
. Perturbations of a given velocity field with respect to the laminar equilibrium are henceforth denoted as
${\boldsymbol{u}_{\!p}} \equiv {\boldsymbol{u}} - {\boldsymbol{u}}_{\textit{lam}}$
.
2.2. Symmetric subspaces
Subject to the given boundary conditions, the Navier–Stokes equations admit several continuous and discrete symmetries. Streamwise periodicity induces a continuous translational shift symmetry in
$x$
, viz.
allowing the existence of TWs (relative equilibria) and relative POs, where ‘relative’ refers to invariance within a co-moving frame of reference. Unlike circular pipe flow, the rectangular cross-section cannot feature a continuous rotational symmetry. Instead, the symmetry group reduces to the discrete dihedral group
$D_2=\{{e},{{\boldsymbol{Z}}_y},{{\boldsymbol{Z}}_z},{{\boldsymbol{R}}_\pi }\}$
. Here,
$e$
denotes the identity,
${{\boldsymbol{Z}}_y},{{\boldsymbol{Z}}_z}$
are the two reflectional symmetries along the
$y$
- and
$z$
-directions, respectively, and
${\boldsymbol{R}}_\pi$
is the
$\pi$
-rotational symmetry around the
$x$
-axis, viz.
\begin{equation} \left . \begin{array}{rlcl} {{\boldsymbol{Z}}_y}:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-v,w\right ]\left (x,-y,z\right )\!, \\[1.5ex] {{\boldsymbol{Z}}_z}:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,v,-w\right ]\left (x,y,-z\right )\!, \\[1.5ex] {{\boldsymbol{R}}_\pi }:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-v,-w\right ]\left (x,-y,-z\right )\!. \end{array} \right . \end{equation}
For the special case of a square cross-section, the symmetry group increases to
$D_4=D_2\cup \{{{\boldsymbol{D}}_{yz}},{{\boldsymbol{D}}_{zy}},{{\boldsymbol{R}}_{\pi /2}},{{\boldsymbol{R}}_{3\pi /2}}\}$
. In addition to the elements of
$D_2$
,
$D_4$
hence contains reflectional symmetries with respect to the corner bisectors
$y=z$
and
$y=-z$
(
${{\boldsymbol{D}}_{yz}},{{\boldsymbol{D}}_{zy}}$
) and rotational symmetries by counter-clockwise rotations of
$\pi /2$
and
$3\pi /2$
about the
$x$
-axis (
${{\boldsymbol{R}}_{\pi /2}},{{\boldsymbol{R}}_{3\pi /2}}$
), viz.
\begin{equation} \left . \begin{array}{rlcl} {{\boldsymbol{D}}_{yz}} :& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,w,v\right ]\left (x,z,y\right )\!, \\[1.5ex] {{\boldsymbol{D}}_{zy}} :& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-w,-v\right ]\left (x,-z,-y\right )\!, \\[1.5ex] {{\boldsymbol{R}}_{\pi /2}} :& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-w,v\right ]\left (x,-z,y\right )\!, \\[1.5ex] {{\boldsymbol{R}}_{3\pi /2}} :& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,w,-v\right ]\left (x,z,-y\right )\!. \end{array} \right . \end{equation}
Combining both streamwise and cross-sectional symmetries, it is readily seen that the governing equations are equivariant under the symmetry group
$\varGamma =SO(2)_x\times D_4$
. In the current work, the ‘shift-reflect’ and ‘shift-rotate’ elements of
$\varGamma$
for fixed streamwise shifts
${L_x}/2 = \pi /{\alpha }$
are of specific interest and therefore deserve an explicit definition
\begin{equation} \left . \begin{array}{rlcl} {{\boldsymbol{S}}_y}:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-v,w\right ]\left (x+\dfrac {\pi }{{\alpha }},-y,z\right )\!, \\[1.75ex] {{\boldsymbol{S}}_z}:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,v,-w\right ]\left (x+\dfrac {\pi }{{\alpha }},y,-z\right )\!, \\[1.75ex] {{\boldsymbol{S}}_\pi }:& \left [u,v,w\right ]\left (x,y,z\right ) &\to & \left [u,-v,-w\right ]\left (x+\dfrac {\pi }{{\alpha }},-y,-z\right )\!. \end{array} \right . \end{equation}
3. Numerical method
3.1. Fluid solver
The pseudo-spectral DNS code developed by Pinelli et al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010) is used to integrate a given flow state forward in time. Time stepping is realised using an incremental pressure projection scheme, in the context of which linear and nonlinear contributions are treated with a semi-implicit Crank–Nicolson and a low-storage Runge–Kutta method (Rai & Moin Reference Rai and Moin1991; Verzicco & Orlandi Reference Verzicco and Orlandi1996), respectively. Throughout, simulation time steps are small enough to ensure that the CFL number stays below
$0.3$
. The velocity and pressure fields are expanded in terms of truncated Fourier series in the periodic streamwise and Chebyshev polynomials in the two cross-stream directions. While a grid of equispaced nodes is introduced in the streamwise direction, the standard Gauss–Chebyshev–Lobatto points are chosen as collocation points in the cross-stream directions. Back and forth transformation between physical and Fourier space along the streamwise direction is performed with the aid of fast Fourier transforms and standard de-aliasing according to the
$2/3$
-rule is applied throughout. Fourier expansion of the fields along the streamwise direction allows us to solve linear Poisson and Helmholtz equations for each streamwise wavenumber separately. The resulting two-dimensional systems only depend on the cross-sectional directions and can be efficiently solved with a fast-diagonalisation technique (Haldenwang et al. Reference Haldenwang, Labrosse, Abboudi and Deville1984). The remaining nonlinear contributions are solved in physical space.
For all edge-tracking simulations, a combination of
${N_x}=\{64,128\}$
(for
${L_x}/\!{{H}}=\{2\pi ,4\pi \}$
) Fourier modes and
${N_y}={N_z}=129$
Chebyshev polynomials are chosen in the streamwise and the two cross-stream directions, respectively. The resulting discretisation fulfils the resolution requirements for fully turbulent trajectories at all considered Reynolds numbers, leading to inner-scaled grid spacings of
${\Delta x^+}\leqslant 14$
and
$\max ({\Delta y^+})=\max ({\Delta z^+})\leqslant 3.5$
. Tests in which edge trajectories were integrated with fewer Chebyshev polynomials
${N_y}={N_z}=\{65,97\}$
gave qualitatively similar results in that the trajectories converged to periodic edge states at comparable period. However, the coarser runs failed to resolve the bursting amplitudes of the higher-resolution data with sufficient precision, justifying the relatively high resolution in the two cross-stream directions. A lower resolution of
${N_x}\times {N_y}\times {N_z}=48\times 65\times 65$
is merely used for the reference simulation at a marginal Reynolds number
${{\textit{Re}}_b}=1150$
, in which case this discretisation is sufficient to resolve all relevant flow scales.
3.2. Edge tracking
Tracking trajectories along the edge is commonly done by a bisection between trajectories that approach the laminar fixed point and those which leave the neighbourhood of the edge towards the turbulent attractor/saddle (Itano & Toh Reference Itano and Toh2001; Toh & Itano Reference Toh and Itano2003; Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006; Schneider et al. Reference Schneider, Eckhardt and Yorke2007; Duguet et al. Reference Duguet, Willis and Kerswell2008; Schneider & Eckhardt Reference Schneider and Eckhardt2009). In the remainder, trajectories that return to the laminar state and those experiencing a rise in energy towards a turbulent state are denoted as
$\boldsymbol{u}^L$
and
$\boldsymbol{u}^H$
, respectively. While similar in concept, the different edge-tracking approaches feature some relevant differences concerning the states chosen on either side of the edge to perform the bisection. In this study, we apply a technique similar to that presented by Schneider et al. (Reference Schneider, Eckhardt and Yorke2007) to follow trajectories along the edge. In this approach, bisection is performed along a ray connecting a state
$\boldsymbol{u}^H$
that will tend to the turbulent ‘side’ and the laminar equilibrium
$\boldsymbol{u}_{\textit{lam}}$
. As long as
${\boldsymbol{u}^H} - {\boldsymbol{u}_{\textit{lam}}}$
is not tangential to
${\mathscr{M}}$
, the procedure per se allows us to track edge trajectories for arbitrary time intervals. However, the positive Lyapunov exponents of the chaotic system cause any initially close trajectories to separate exponentially fast. Hence, even if the bisection determines a flow state on the edge up to double precision, a ‘refinement step’ after
${\delta t}=\mathcal{O}(100{T_b})$
is unavoidable to find a new approximation for a state on
${\mathscr{M}}$
(Toh & Itano Reference Toh and Itano2003; Duguet et al. Reference Duguet, Willis and Kerswell2008). The finally obtained solution is therefore – strictly speaking – a piecewise-continuous approximation to a true edge trajectory, with small-amplitude discontinuities occurring at the refinement times
$t_i$
(Itano & Toh Reference Itano and Toh2001).
In the current work, the first bisection is initialised with a flow state
${\boldsymbol{u}}_0({\boldsymbol{x}},t_{0})$
, randomly chosen from a turbulent trajectory at matching
${\textit{Re}}_b$
. The general procedure can then be summarised as follows: the
$i$
th bisection step iteratively seeks a parameter
$\beta _i\in {\mathbb{R}}$
(
$\forall i\in {\mathbb{N}}$
) that rescales the velocity perturbation of a state
${\boldsymbol{u}}^H_{i-1}({\boldsymbol{x}},t_{i-1})$
with respect to the laminar base flow at the beginning of the corresponding time interval
$t_{i-1}$
, viz.
\begin{equation} \left \{ \begin{array}{rcll} {\boldsymbol{u}}_{\beta _i}({\boldsymbol{x}},t_{i-1}) \!\! &\!\!\!\!\!\equiv\!\!\!\!\! &\!\! {\boldsymbol{u}}_{\textit{lam}} + \beta _i \left ({\boldsymbol{u}}^H_{i-1}({\boldsymbol{x}},t_{i-1}) -{\boldsymbol{u}}_{\textit{lam}} \right )\quad & \text{at}\;\, t=t_{i-1}, \\[0.8ex] {\boldsymbol{u}}_{\beta _i}({\boldsymbol{x}},t)\!\! &\!\!\equiv\!\! & \!\!{\boldsymbol{u}}_{\beta _i}({\boldsymbol{x}},t_{i-1}) + \displaystyle \int \limits _{t^{\prime}=t_{i-1}}^{t} \dfrac {\partial }{\partial t} {\boldsymbol{u}} \;{\textrm{d}}t^{\prime} \quad & \forall t\in (t_{i-1},t_{i}). \end{array} \right . \end{equation}
In each of these shooting steps, the rescaled flow state is advanced in time according to (3.1) with the aid of the DNS code outlined in § 3.1, until the flow either relaminarises or becomes turbulent. As a measure for when a trajectory has left the edge, we use the root mean square (r.m.s.) of the kinetic energy of the non-constant streamwise Fourier modes, viz.
with
\begin{equation} [{{E}_{3D}}]_j= \displaystyle \sum \limits _{i=1}^{{N_x}-1} \;\; \displaystyle \iint \limits _{y,z=-{{H}}}^{{{H}}} \; [{\hat {\boldsymbol{u}}}_{i}]_j[{\hat {\boldsymbol{u}}}_{i}^{*}]_j \;{\textrm{d}}z{\textrm{d}}y \quad \forall j\in \{x,y,z\}. \end{equation}
Therein,
${\hat {\boldsymbol{u}}}_{i}$
and
${\hat {\boldsymbol{u}}}_{i}^{*}$
indicate the
$i$
th Fourier coefficient and its complex conjugate, respectively, and
$[\bullet ]_j$
denotes the
$j$
th component of a vector field. The trajectory is considered as deviated from
${\mathscr{M}}$
when
${E}_{\textit{rms}}$
either falls below
${{E}_{\textit{rms}}^l}\approx 2\times10^{-3}{u_b}$
or, following a steep energy rise, exceeds
${{E}_{\textit{rms}}^t}\approx 2\times10^{-1}{u_b}$
. Each bisection is continued until the relative error
${\vert {\beta ^H-\beta ^L} \vert }/{\vert {2\beta _i} \vert }$
falls below
$10^{-12}$
, where
$\beta _i$
represent the current approximation to the parameter
$\beta$
. The bracketing parameter values, eventually leading to a relaminarisation or a turbulent state, are denoted by
$\beta ^L$
and
$\beta ^H$
, respectively.
Since the process is started with a turbulent flow state, the first parameter value
$\beta _1$
is usually clearly lower than unity. For short enough refinement intervals, subsequent bisections start with a much better initial guess for a state within
${\mathscr{M}}$
, so that all
$\beta _i \;\forall i\geqslant 2$
are very close to unity and, thus, the computed piecewise-continuous curve is a good approximation to a true edge trajectory. In order to ensure this quality of the approximation, the trajectories
$\boldsymbol{u}^L$
and
$\boldsymbol{u}^H$
are advanced in time over
${\delta t}=300{T_b}$
(in single cases up to
${\delta t}=600{T_b}$
), as long as it is guaranteed that the relative distance between the two bracketing states
$\boldsymbol{u}^L$
and
$\boldsymbol{u}^H$
\begin{equation} \varepsilon _{\textit{LH}} \equiv \dfrac {{\left \lVert {{\boldsymbol{u}^H}(t_i + {\delta t})-{\boldsymbol{u}^L}(t_i + {\delta t})} \right \rVert _{2}}} {{\left \lVert {{\boldsymbol{u}^H}(t_i + {\delta t})} \right \rVert _{2}} + {\left \lVert {{\boldsymbol{u}^L}(t_i + {\delta t})} \right \rVert _{2}}} \end{equation}
remains below a threshold of
$10^{-4}$
. In case
$\varepsilon _{\textit{LH}}$
exceeds this limit,
$\delta t$
is successively reduced until the criterion is fulfilled (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013).
4. Results
4.1. Edge-state dynamics in the full state space
In a first step, edge trajectories are sought in the full phase space without restriction to any symmetric subspace. The characteristic dynamics will be analysed for trajectories in two domains of different streamwise period,
${L_x}/\!{{H}}\in \{2\pi ,4\pi \}$
, at Reynolds number
${{\textit{Re}}_b}=2000$
(equivalent to a friction Reynolds number
${{\textit{Re}}_\tau } \equiv {{H}}^+ \approx 89$
in the turbulent case). Trajectories at other Reynolds numbers up to
${{\textit{Re}}_b}=3200$
revealed a qualitatively similar behaviour, so that the current findings can be seen as representative of a wider range of Reynolds numbers. In the remainder of this work, the two considered cases will be referred to as
${\textrm{CA}}{2000}_{2\pi }$
and
${\textrm{CA}}{2000}_{4\pi }$
, respectively, and the properties of their edge states are summarised in table 1, together with the periodic edge states discussed in the subsequent section.
The here analysed Reynolds numbers
${{\textit{Re}}_b}\geqslant 2000$
lie well above the critical value for the sustainance of turbulence
${\textit{Re}}_b^{c}\approx 1077$
(equivalent to
${\textit{Re}}_{\tau }^{c}\approx 77$
in the turbulent case), so that all four bounding walls simultaneously feature a pronounced turbulent activity most of the time. A clear separation between the typical perturbation energy level within the edge and along fully turbulent trajectories is for these Reynolds numbers guaranteed. Note that the latter is vital for the edge-tracking protocol, as it allows an unambiguous classification of whether or not a trajectory has left the edge towards the turbulent basin. For much lower Reynolds numbers approaching
${\textit{Re}}_b^{c}$
, on the other hand, states along edge trajectories and those in turbulence are of comparable amplitude and can no longer be properly distinguished from each other (Schneider & Eckhardt Reference Schneider and Eckhardt2009; Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013).

Figure 2. Time evolution of different flow measures for the non-symmetric edge trajectories of cases (a,c,e)
${\textrm{CA}}{2000}_{2\pi }$
and (b,d,f)
${\textrm{CA}}{2000}_{4\pi }$
. (a,b) The r.m.s. energy signal
${{E}_{\textit{rms}}}(t)$
, with solid nodes indicating snapshots that are visualised in subsequent figures. (c,d) Individual contributions
$S_i(t)$
of the triangular sectors to the overall mean streamwise enstrophy. (e,f) Shear stress averaged along each of the four walls separately,
${\tau _{w,i}}(t)$
, normalised by the corresponding laminar value
${\tau _{w,i,lam}}(t)$
. In (c–f), colour coding indicates the respective wall/sub-sector:
$i={S}$
(blue),
$i={E}$
(green),
$i={N}$
(orange),
$i={W}$
(red).
At
${{\textit{Re}}_b}=2000$
, in the absence of imposed symmetries, the edge states in the considered domains are relative chaotic attractors within
${\mathscr{M}}$
: trajectories initiated with appropriately scaled random turbulent snapshots quickly approach a chaotic edge state, as can be seen from the temporal variation of the perturbation energy signals
${E}_{\textit{rms}}$
depicted in figure 2(a,b). A recurrent pattern can nonetheless be found in all analysed chaotic edge trajectories: throughout, the dynamics is characterised by long-lasting, relatively quiescent phases of
$\mathcal{O}(100{T_b})$
length, which are suddenly interrupted by intense bursting events. These latter manifest themselves in form of a rapid rise of the systems’ perturbation energy, followed by a decrease back to the energy level of the quiescent period. In wall-bounded turbulence, such bursting events are usually attributed to the bending and eventual breakup of a streamwise velocity streak due to the action of flanking quasi-streamwise vortices (Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997). In order to deduce near which of the four walls such a bursting event is happening, the enstrophy contained in the streamwise-averaged vorticity (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007)
\begin{equation} S_i(t)=\displaystyle \int \limits _{{\varOmega }_i} {\langle {{\omega _{x}}} \rangle} _{x}^2 \, {\textrm{d}}A \qquad \quad \forall i\in \{W,S,E,N\} \end{equation}
is measured for the four disjoint triangular subsets
of the duct cross-section
${\varOmega _{\perp }}=[-{{H}},{{H}}]\times [-{{H}},{{H}}]$
separately. In order to quantify which of the two pairs of parallel walls features a significant turbulent activity, we define in addition a dimensionless indicator function (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007)
based on the triangular sector contributions defined in (4.1). By definition,
${\textit{Id}}_{\triangle }$
attains values close to plus or minus unity whenever the enstrophy carried by the streamwise-averaged vorticity is predominantly focussed near the walls at
$y\in \{-{{H}},{{H}}\}$
or those at
$z\in \{-{{H}},{{H}}\}$
, respectively.
The temporal variation of the different contributions
$S_i(t)$
presented in figure 2(c,d) clearly indicates that, over a significant portion of the observation interval, the turbulent vortical activity is primarily concentrated near a single wall. As has been shown earlier, this (instantaneous) ‘localisation’ of streaks and vortices to a subset of the four surrounding walls only is not exclusive to edge trajectories, but is visible in marginally turbulent flows as well (cf. figure 1
a). Interestingly, the concentration to a single wall is strongest at the peaks of the bursting phases, with over
$99\,\%$
(
${\textrm{CA}}{2000}_{2\pi }$
) and
$97\,\%$
(
${\textrm{CA}}{2000}_{4\pi }$
) of the enstrophy associated with
$\langle {{\omega _{x}}} \rangle _{x}$
residing in the vicinity of a single wall during these times. Also, figure 2(c,d) reveals that each bursting episode is accompanied by a ‘switching’ of the activity to one of the neighbouring walls. During this switching period, in turn, turbulent activity spreads temporarily over several triangular sectors
${\varOmega }_{i}$
, before focussing again in a single one. The sequence in which the individual contributions
$S_{i}$
dominate the overall vortical activity moreover implies that the ‘switching’ is not necessarily continuously oriented in one direction, but that a back-and-forth switching is equally possible. This is interesting as a qualitatively similar chaotic dynamics is reported for the lateral ‘streak switching’ in the chaotic edge states of the asymptotic suction boundary layer (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013), in this case in form of an irregular left–right switching. In the square duct case, the ‘wall-switching’ dynamics also leaves its footprint in the temporal evolution of the wall shear stress, separately averaged over the four surrounding walls: as can be seen from figure 2(e,f), the wall shear stress
$\tau _{w,i}$
of ‘active walls’ can exceed the values in a corresponding laminar flow by approximately
$70\,\%$
. The vicinity of the remaining walls is then usually quiescent, with an averaged wall shear stress
${\tau _{w,i}}\approx {\tau _{w,i,lam}}$
, underlining the absence of a self-sustaining near-wall regeneration cycle along these walls.
The different stages of a selected bursting event in the interval
$t/{T_b}\in [600,900]$
and the corresponding wall switching are shown in figure 3 in terms of the streamwise-averaged velocity field
${\langle {\boldsymbol{u}} \rangle _{xt}}$
, extracted from the edge trajectory at the instances indicated by solid markers in figure 2(a). Figure 4 complements this view with three-dimensional visualisations of the low-speed streaks (depicted as iso-surfaces of
$u_{\!p}$
) and the flanking quasi-streamwise vortices (measured as regions of high values of the second invariant of the velocity gradient tensor,
$Q$
; see Hunt et al. Reference Hunt, Wray and Moin1988) at the same moment in time. Figure 4(a) reveals that the flow state just before the energetic burst excursion is characterised by a single pronounced low-velocity streak near the centre of the top wall, flanked by a pair of counter-rotating quasi-streamwise vortices in a streamwise staggered arrangement. The vortices are inclined with respect to both the neighbouring wall and the streamwise direction, so that they induce a four-vortex cell pattern in the streamwise-averaged field (cf. figure 3
a). While this configuration promotes the evolution of the low-speed streak in the direct vicinity of the wall, it counteracts it further away – until the low-speed streak eventually turns into a high-speed streak as the burst reaches its peak (figures 3
c and 4
c). Here, the special shape of the duct cross-section comes into play, since the quasi-streamwise vortices responsible for the birth of the high-speed streak at the top wall simultaneously induce new low-speed streaks along the two sidewalls at
$z=\{-{{H}},{{H}}\}$
, while the energy quickly tends back to the pre-burst level (figures 3
d and 4
d). The flow is, however, not entirely symmetric with respect to the
$y$
-axis and, eventually, a pronounced low-speed streak surrounded by a pair of quasi-streamwise vortices is only evolving on one of the sidewalls at
$z=-{{H}}$
(figures 3
e,f and 4
e,f). In contrast to their counterparts at the top wall just before the bursting event, these structures are still in an early stage of their lifetime: the streak is still aligned with the streamwise direction, while the vortices do not lean over the streak yet, as can be seen from the streamwise-averaged secondary flow patterns in figure 3(e,f).

Figure 3. Selected snapshots of the instantaneous streamwise-averaged primary and secondary velocity fields for case
${\textrm{CA}}{2000}_{2\pi }$
. The snapshots are extracted at different times along a selected bursting event
$t/{T_b}\in [600,900]$
, as indicated by the markers in figure 2(a): (a–f)
$t/{T_b}=\{675.4, 712.7, 744.5, 787.4, 825.3, 862.6\}$
. Isovalues are the same as in figure 1.

Figure 4. Three-dimensional visualisations of the velocity field in case
${\textrm{CA}}{2000}_{2\pi }$
in terms of low-speed streaks and quasi-streamwise vortices during the selected bursting event
$t{\kern-0.5pt}/{T_b}\in [600,900]$
, extracted at the same times as those in figure 3 and as indicated by the markers in figure 2(a). Light blue iso-surfaces of the streamwise perturbation velocity
${u_{\!p}}=-0.15{u_b}$
indicate the position of low-speed streaks. Quasi-streamwise vortices are identified in terms of the
$Q$
-criterion of Hunt, Wray & Moin (Reference Hunt, Wray and Moin1988) as regions with
$Q \gt 0.6 \max _{\varOmega }(Q)$
. Clockwise (red) and counter-clockwise rotation (dark blue) is measured by the local sign of
$\omega _{x}$
.
The successive change of a mild low-speed streak to an intense high-speed region in two different bursting events at
$t\approx 750{T_b}$
(top wall
$y={{H}}$
) and
$t\approx 1320{T_b}$
(right sidewall
$z={{H}}$
) can be clearly identified by the footprint these events leave in the wall shear stress profiles in figure 5. Therein, the temporal evolution of the streamwise-averaged wall shear stress (
$\vert {{\tau _{xy}}} \vert$
or
$\vert {{\tau _{xz}}} \vert$
) is shown for all four walls separately, and is contrasted with the time signal of the perimeter-averaged friction Reynolds number. The visualisation underlines that the high shear stress declines quickly in
$\mathcal{O}(10{T_b})$
after the bursting event, but its value at the corresponding wall stays the dominant contribution over a much longer time interval. In between the two bursting events, when the turbulent activity is rather evenly distributed across the different triangular sectors, there is also no pronounced high- or low-speed streak detectable on either of the four walls. The edge state therefore features two characteristic time scales, a slow one with
$\mathcal{O}(100){T_b}$
that is associated with the interval between two bursting events and a much faster related to the structures’ downstream propagation (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013).

Figure 5. Variation of the streamwise-averaged wall shear stress along each wall during two bursting episodes in the interval
$t/{T_b}\in [600,1500]$
of case
${\textrm{CA}}{2000}_{2\pi }$
:
$\langle {{\tau _{xz}}} \rangle _{x}$
along (a) the left (
$z=-{{H}}$
) and (b) right sidewalls (
$z={{H}}$
), as well as
$\langle {{\tau _{xy}}} \rangle _{x}$
along (c) the bottom (
$y=-{{H}}$
) and (d) top walls (
$y={{H}}$
). The wall shear stress is normalised with the perimeter-averaged value in the corresponding laminar state. Isocontours are drawn in the interval
$[0,3.5]$
, with increment
$0.25$
. (e) Temporal variation of the perimeter-averaged friction Reynolds number
${\textit{Re}}_\tau$
over the same time window.
The observed streak–vortex interaction, including the breakup and bursting of the streaks, agrees well with the different phases of the classical self-sustained process/vortex–wave interaction envisaged by Hamilton et al. (Reference Hamilton, Kim and Waleffe1995), Waleffe (Reference Waleffe1997) and Hall & Sherwin (Reference Hall and Sherwin2010). To provide quantitative evidence, the kinetic perturbation energy per unit mass is decomposed into different contributions that can be associated with the different ingredients of the process (Lucas & Kerswell Reference Lucas and Kerswell2017), viz.
\begin{equation} {{E_{\!p}}} = \dfrac {1}{2} {\big\langle {{\boldsymbol{u}}_{\!p}^2} \big\rangle _{V}} = {\underbrace {\dfrac {1}{2}{\big\langle {{\langle {{u_{\!p}}}} \rangle _{x}^2} \big\rangle _{V}}}_{{{{E_{\textit{streak}}}}}}} + {\underbrace {\dfrac {1}{2}{\big\langle {{\langle {{v_{\!p}}} \rangle} _{x}^2+{\langle {{w_{\!p}}} \rangle} _{x}^2} \big\rangle _{V}}}_{{{{E_{\textit{rolls}}}}}}} + {\underbrace {\dfrac {1}{2}{\big\langle {\left ({\boldsymbol{u}_{\!p}}-{\langle {{\boldsymbol{u}_{\!p}}} \rangle _{x}}\right )^2} \big\rangle _{V}}}_{{{{E_{{wave}}}}}}}. \end{equation}
Here,
${E_{\textit{streak}}}$
,
${E_{\textit{rolls}}}$
and
${E_{\textit{wave}}}$
refer to the energy contributions of the streamwise constant streak mode, the quasi-streamwise vortices/rolls and the streamwise-varying wave mode, respectively. In the square duct,
${E_{\textit{rolls}}}$
effectively measures the kinetic energy of the instantaneous secondary flow. Despite being the dominant contribution throughout, figure 6 highlights that the intensity of the streamwise-elongated streaks
${E_{\textit{streak}}}$
regularly drops during the identified bursting intervals. In these phases, the energy is instead transferred into the streamwise-dependent wave field (
${E_{\textit{wave}}}$
), that is, the initially streamwise-aligned streak gets more and more distorted. With a slight delay of several bulk time units, a rise of the roll mode energy
${E_{\textit{rolls}}}$
follows, which indicates the gain in strength of the quasi-streamwise vortices. Owing to the much smaller amplitude of the cross-stream velocity components compared with their streamwise counterpart, the kinetic energy residing in the quasi-streamwise vortices is two orders of magnitude smaller than that of the wave field. The outlined observations are consistent for the shorter and longer domain, with somewhat higher portions of the total perturbation energy
${E_{\!p}}$
being held in
${E_{\textit{wave}}}$
for the longer domain.

Figure 6. Time evolution of the individual contributions to the perturbation energy
${E_{\!p}}$
for cases (a)
${\textrm{CA}}{2000}_{2\pi }$
and (b)
${\textrm{CA}}{2000}_{4\pi }$
. Line styles indicate:
${E_{\textit{streak}}}$
(solid),
${E_{\textit{wave}}}$
(dashed),
$150{{E_{\textit{rolls}}}}$
(dotted). The insets show the evolution of the wave and roll modes’ amplitude enlarged during selected bursting events.
The fact that the edge state in the square duct is a chaotic attractor within
${\mathscr{M}}$
is consistent with circular pipe flows, where the edge state takes the same form (Duguet et al. Reference Duguet, Willis and Kerswell2008; Mellibovsky et al. Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009; Schneider & Eckhardt Reference Schneider and Eckhardt2009; Avila et al. Reference Avila, Barkley and Hof2023). Given that many TW solutions detected in circular pipe flow (Pringle, Duguet & Kerswell Reference Pringle, Duguet and Kerswell2009) possess a strikingly similar counterpart in the square duct (Okino Reference Okino2011), the current findings provide further support for the suspected close relation between the state space structure of both systems. When it comes to a direct comparison of the individual dynamics, however, it turns out that comparable alternating quiescent phases and intermittent bursts have not been reported for the chaotic edge states of circular pipe flows. In other systems including plane channels (Rawat et al. Reference Rawat, Cossu and Rincon2014, Reference Rawat, Cossu and Rincon2016; Toh & Itano Reference Toh and Itano2003; Zammert & Eckhardt Reference Zammert and Eckhardt2014) or asymptotic suction boundary layers (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013), however, the edge states do feature a remarkably similar dynamics. For certain (streamwise and spanwise) domain periods and Reynolds numbers, the edge state in these systems even simplifies to a symmetric time-periodic limit cycle within
${\mathscr{M}}$
, or an unstable PO with respect to the full state space; a more detailed discussion of these solutions and comparison with the edge states in square duct will follow in § 4.2.
4.2. Periodic edge-state dynamics in symmetric subspaces
In circular pipe flow, edge states reduce to simple invariant solutions merely if trajectories are restricted to appropriate symmetric subspaces of the full state space. Examples include TW edge states in
$\pi$
-rotational symmetric flows (Duguet et al. Reference Duguet, Willis and Kerswell2008) and (relative) PO edge states under enforcement of a
$\pi$
-rotational, combined with a reflectional symmetry (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013). Based on these experiences, Gepner et al. (Reference Gepner, Okino and Kawahara2025) (cf. also Okino Reference Okino2014, in Japanese) performed edge-tracking studies in the
${\boldsymbol{R}}_\pi$
-symmetric subspace of square duct flows driven with a constant pressure gradient. For domains of
${L_x}/\!{{H}}=2\pi$
and
${L_x}/\!{{H}}=8\pi$
length, these authors found the edge states to be TWs equivariant under the symmetry groups
$({{\boldsymbol{S}}_y},{{\boldsymbol{S}}_z},{{\boldsymbol{R}}_\pi })$
and
$({{\boldsymbol{Z}}_y},{{\boldsymbol{Z}}_z},{{\boldsymbol{R}}_\pi })$
, respectively. For the intermediate domain length
${L_x}/\!{{H}}=4\pi$
, however, the edge state was chaotic even in the
${\boldsymbol{R}}_\pi$
-symmetric subspace. In the current study, we found a TW edge state with the same characteristics as the one found by Okino (Reference Okino2014) in the
${\boldsymbol{R}}_\pi$
-symmetric subspace (with the solution featuring a
${{\boldsymbol{S}}_y}{{\boldsymbol{S}}_z}{{\boldsymbol{R}}_\pi }$
-symmetry) despite our somewhat different driving, i.e. a time-varying pressure gradient ensuring a constant mass flow rate versus a constant pressure gradient in the study of Okino (Reference Okino2014). Both solutions arguably belong to the same family of solutions, proving the consistency of our results with theirs.
In the course of this work, edge trajectories have been analysed in a number of different subspaces, obeying varying combinations of shift-reflect symmetries (
${{\boldsymbol{S}}_y},{{\boldsymbol{S}}_z}$
), rotational symmetries (
${\boldsymbol{R}}_\pi$
) and full reflectional symmetries (
${{\boldsymbol{Z}}_y},{{\boldsymbol{Z}}_z}$
). In most of the subspaces, the edge states turned out to be TWs known in the community, even though not necessarily in their role as edge states (Wedin et al. Reference Wedin, Bottaro and Nagata2009; Uhlmann et al. Reference Uhlmann, Kawahara and Pinelli2010; Okino Reference Okino2014). In the remainder of this section, we therefore restrict ourselves to the
${\boldsymbol{Z}}_y$
- and
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspaces, for which edge trajectories at various Reynolds numbers in the range
${{\textit{Re}}_b}\in [1600,3200]$
converged to limit cycles with a rich bursting dynamics. At selected Reynolds numbers, edge trajectories have been computed starting from different turbulent snapshots as initial conditions. For almost all parameter points, the different edge trajectories eventually converge to the same periodic edge states modulo discrete or continuous symmetry operations. The only exception is the
${\boldsymbol{Z}}_y$
-symmetric case
${\textrm{PO}}{2000}_y$
, for which one edge trajectory converged to the PO shown in figure 7(
b), while the other approached a TW. The latter has been successfully converged with a standard Newton–Raphson method, using the generalised minimal residual (GMRES) approach to approximately solve the arising linear systems (Chandler & Kerswell Reference Chandler and Kerswell2013). The detected TW lives in a higher-symmetric subspace invariant under
$({{\boldsymbol{Z}}_y},{{\boldsymbol{S}}_z},{{\boldsymbol{S}}_\pi })$
and presumably belongs to the family of solutions found by Wedin et al. (Reference Wedin, Bottaro and Nagata2009), although no effort was made to prove this rigorously via a continuation study. To show that the TW is indeed an alternative edge state rather than a relative saddle within
${\mathscr{M}}$
, we continued tracking the edge trajectory for several thousand bulk time units after it had approached the TW, without any sign of deviation from the latter. State space configurations with multiple co-existing ‘local’ edge states of possibly varying type (e.g. one TW and one PO as observed here) are not exclusive to the square duct, but have been reported for a variety of flows (Duguet et al. Reference Duguet, Willis and Kerswell2008; Suri et al. Reference Suri, Pallantla, Schatz and Grigoriev2019).

Figure 7. Time evolution of the r.m.s. energy signal
${{E}_{\textit{rms}}}(t)$
(upper plots) and individual contributions of the triangular sectors to the overall mean streamwise enstrophy
$S_i(t)$
(lower plots) for periodic edge states in the
${\boldsymbol{Z}}_y$
-symmetric subspace: (a)
${\textrm{PO}}{1600}_y$
, (b)
${\textrm{PO}}{2000}_y$
, (c)
${\textrm{PO}}{2400}_y$
and (d)
${\textrm{PO}}{3200}_y$
. In (a–c), the fundamental period (grey background) is repeated periodically. For the sake of visualisation, in (d), we show only the first half of the POs period (blue background) which represents a pre-PO with period
$T_{\textit{pre}}={T}/2$
(see the discussion in the main text). Colour coding in the lower plots is:
$i={E}$
(green),
$i={W}$
(red),
$i=\{S,N\}$
(orange, identical due to symmetry).
In the remainder, we will restrict our analysis to the periodic edge states. Unfortunately, the limit cycles’ long periods
$T=\mathcal{O}(100{T_b})$
do not allow us to converge these solutions to machine precision with the standard Newton–GMRES approach; a problem that is shared with time-periodic edge states in other systems (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013; Zammert & Eckhardt Reference Zammert and Eckhardt2014). If at all possible, convergence might be achieved with more advanced multipoint-shooting techniques (Budanur et al. Reference Budanur, Short, Farazmand, Willis and Cvitanović2017), or recently developed adjoint-based variational Jacobian-free methods (Ashtari & Schneider Reference Ashtari and Schneider2023). Here, however, time stepping alone provides a fairly good approximation to the PO since the latter is an attracting state in
${\mathscr{M}}$
: the relative distance between two velocity fields separated by a full period
$T$
, viz.
is small – of order
$\mathcal{O}(10^{-5})$
when only the streamwise-dependent modes
${\boldsymbol{u}}-{\langle {{\boldsymbol{u}}} \rangle _{x}}$
are compared in (4.5), and
$\mathcal{O}(10^{-8})$
if the full set of modes is considered. Here, the period of the limit cycles is estimated up to a precision of
$5{\Delta t} \approx 0.07{T_b}$
(which is half the interval in which statistics have been aggregated) based on the peaks of the auto-correlation functions of different measures; see Appendix A for further details and Appendix B for an explanation how the stability of the limit cycles has been verified based on first-return maps of
${{\textit{Id}}_{\triangle }}(t)$
(Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Zammert & Eckhardt Reference Zammert and Eckhardt2014). Note that the minimisation over all possible streamwise shifts
$s_x$
in (4.5) checks for relative POs, which are self-recurrent in a co-moving frame of reference only. It turns out that the identified edge states are relative POs indeed: In case
${\textrm{PO}}{2400}_y$
, for instance,
${\boldsymbol{u}}({\boldsymbol{x}}+{\boldsymbol{s}},t+T)={\boldsymbol{u}}({\boldsymbol{x}},t)$
holds for
${s_x}/\!{{H}}\approx 0.11$
.
4.2.1.
$\mathscr{M}$
under
${\boldsymbol{Z}}_y$
-symmetry
In a first step, we analyse states whose velocity and pressure fields obey a single reflectional symmetry with respect to the wall bisector. As for the chaotic edge states shown in figure 2(
a,b), the edge trajectories in the
${\boldsymbol{Z}}_y$
-symmetric subspace experience regular bursting events that interrupt the otherwise quiescent dynamics (cf. figure 7, upper plots in each panel) for all considered Reynolds numbers. Per periodic cycle, the flow undergoes – depending on the Reynolds number – between two and ten such bursting episodes for
$E_{\textit{rms}}$
with different peak amplitudes. In all cases, the strongest bursts go hand in hand with a rapid increase in energy by approximately a factor three compared with the pre-burst phase.
The temporal variation of the individual sector contributions to the streamwise-averaged enstrophy,
$S_i(t)$
, is presented in figure 7 (lower plots in each panel). For all but case
${\textrm{PO}}{3200}_y$
, the two consecutive peaks of the energy signals come with a subsequent switching from a ‘top-bottom’ (higher peak) to a ‘left-only’ (weaker peak) configuration, or vice versa. The two configurations are seen in the streamwise-averaged velocity fields, extracted along a full cycle of case
${\textrm{PO}}{2400}_y$
(cf. figure 8): the higher energetic peak encodes the simultaneous bursting of a streak at the top wall and its symmetric twin at the bottom wall (figure 8
b), while the lower indicates a bursting event at one of the sidewalls (figure 8
d). The streak–vortex groups at the top and bottom walls arguably undergo the same self-sustaining process as in the chaotic edge state, with the difference that the symmetrical constrains here enforce them to appear in pairs at the two opposite walls. In analogy to the dynamics of the chaotic edge state, each of the quasi-streamwise vortices responsible for the change of the low-speed into a high-speed streak at the bottom and top wall induces a low-speed streak at the neighbouring sidewall – a first one in the lower half (
$y\!/\!{{H}}\lt 0$
) and a second one in the upper half (
$y\!/\!{{H}}\gt 0$
) of the domain (figure 8
c). The two mirror-symmetric streaks then experience a similar instability as their larger counterparts at the top and bottom wall, until they finally break (figure 8
d) and the associated quasi-streamwise vortices give rise to a new generation of streaks at the top and bottom walls (figure 8
f). A comparison of this last state in figure 8
f with the one in figure 8
a once more confirms the exact recurrence of the orbit after a period
${T}=653.3{T_b}$
.

Figure 8. Snapshots of the streamwise-averaged primary and secondary flow fields of case
${\textrm{PO}}{2400}_y$
extracted at different instances along the orbit, as indicated by the markers in figure 7(
c): (a–f)
$t/{T_b}=\{88.0, 166.6, 226.3, 635.4, 653.2, 741.3\}$
. Isovalues are the same as in figure 1.
The switching behaviour between a single streak–vortex family at the bottom and top wall on the one hand and the two narrow streaks with associated vortices at one of the sidewalls leaves a footprint in the temporal evolution of the streamwise-averaged wall shear stress, cf. figure 9. On the other hand, the presented data also indicates that the remaining wall at
$z={{H}}$
never hosts a pronounced streak, again highlighting that the switching occurs only between the top/bottom walls and the sidewall at
$z=-{\!{H}}$
in this case. The observed behaviour is qualitatively similar for the lower-Reynolds-number cases
${\textrm{PO}}{1600}_y$
,
${\textrm{PO}}{2000}_y$
,
${\textrm{PO}}{2200}_y$
(not shown) and
${\textrm{PO}}{2400}_y$
, the only difference being that the period of the orbits monotonically increases with the bulk Reynolds number (cf. discussion of figure 13). Such a trend is consistent with observations for periodic edge states in other flow configurations (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013). These characteristics change for the highest considered Reynolds number (case
${\textrm{PO}}{3200}_y$
), for which the edge state is still periodic, but the dynamics is more complex and the period is with approximately
$3800{T_b}$
one order of magnitude longer than those at the lower Reynolds numbers. In contrast to these latter cases, the edge state in case
${\textrm{PO}}{3200}_y$
switches between three rather than two distinct states in a fairly complicated manner (cf. figure 7
d). The sequence of dominant contributions
$S_i$
over a single cycle reads
i.e. the flow changes between a ‘top-bottom’ state and a flow state with turbulent vortical activity focussing near either of the two sidewalls. It turns out that the two sequences marked as
$\text{pre-PO}_1$
and
$\text{pre-PO}_2$
represent two pre-POs (Cvitanović et al. Reference Cvitanović, Artuso, Mainieri, Tanner, Vattay, Whelan and Wirzba2022, p. 187): after
${T}_{\textit{pre}}={T}/2$
, the flow state repeats itself modulo a reflection such that
${\boldsymbol{u}}({\boldsymbol{x}},t) = {{\boldsymbol{Z}}_z}{\boldsymbol{u}}({\boldsymbol{x}},t + {T}_{\textit{pre}})$
and
${\boldsymbol{u}}({\boldsymbol{x}},t) = {{\boldsymbol{Z}}_z}{{\boldsymbol{Z}}_z}{\boldsymbol{u}}({\boldsymbol{x}},t + {T}) = {\boldsymbol{u}}({\boldsymbol{x}},t + {T})$
holds, ignoring possible streamwise shifts for readability. Note that such sudden changes of the edge-state characteristics due to parameter variations (e.g. changing domain extents or Reynolds number) are reported for edge states in many flow configurations, sometimes even changing from a simple TW or PO to a chaotic state (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013; Zammert & Eckhardt Reference Zammert and Eckhardt2014).

Figure 9. Variation of the streamwise-averaged wall shear stress along each wall over a full period of case
${\textrm{PO}}{2400}_y$
:
$\langle {{\tau _{xz}}} \rangle _{x}$
along (a) the left sidewall (
$z=-{{H}}$
) and (b) right sidewall (
$z={{H}}$
) and (c)
$\langle {{\tau _{xy}}} \rangle _{x}$
along the top/bottom walls (
$y=\pm {{H}}$
). The wall shear stress is normalised with the perimeter-averaged value in the corresponding laminar state. Isocontours are drawn in the interval
$[0,3.5]$
, with increment
$0.25$
. (d) Temporal variation of the perimeter-averaged friction Reynolds number
${\textit{Re}}_\tau$
over the same time window.
4.2.2.
$\mathscr{M}$
under
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetry
In a second step, the flow is further restricted to the twofold mirror-symmetric subspace for which trajectories are invariant under the operations
${\boldsymbol{Z}}_y$
and
${\boldsymbol{Z}}_z$
. It is then readily shown that states in this subspace automatically fulfil a rotational symmetry by
$\pi$
about the
$x$
-axis (
${\boldsymbol{R}}_\pi$
) as well. As for the
${\boldsymbol{Z}}_y$
-symmetric cases, all considered trajectories at various Reynolds numbers converge after an initial transient to a relative PO as edge state. The converged limit cycles can be identified in terms of
${E}_{\textit{rms}}$
in figure 10 for different Reynolds numbers. The general dynamics is similar to that of their counterparts in the
${\boldsymbol{Z}}_y$
-symmetric subspace in that quiescent phases alternate again with intermittent intense high energy and dissipation excursions. Here, however, the systems oscillate periodically between two states in which the overall turbulent activity is almost entirely focussed at the top and bottom walls at
$y=\pm {{H}}$
or the two sidewalls at
$z=\pm {{H}}$
, respectively. All edge states with
${{\textit{Re}}_b}\geqslant 2000$
(cf. cases
${\textrm{PO}}{2000}_{yz}$
and
${\textrm{PO}}{3200}_{yz}$
in figure 10
c,d as examples) undergo two bursting events with identical peak amplitudes per limit cycle, from which one is associated with a streak breakup at the top/bottom walls and the other one with a corresponding event at the sidewalls. Figure 11, in which the temporal evolution of the streamwise-averaged wall shear stress is presented for case
${\textrm{PO}}{2000}_{yz}$
, reveals that not only the peak amplitudes are identical: the limit cycle consists of two pre-POs with
${T}_{\textit{pre}} = {T}/2$
that are identical modulo
${\boldsymbol{R}}_{\pi /2}$
or
${\boldsymbol{R}}_{3\pi /2}$
, i.e. rotations by
$\pi /2$
or
$3\pi /2$
about the
$x$
-axis. Hence, the relative distance between two velocity fields separated by
${T}_{\textit{pre}}$
(following the definition in (4.5) including a correction for the rotation) is at
$\mathcal{O}(10^{-5})$
of the same order as when evaluated for the full period
$T$
. In figure 12, therefore, instantaneous streamwise-averaged velocity fields are shown for the first half of the limit cycle only. A comparison between the states in figure 12(b,f) visually confirms that any two states separated by
${T}_{\textit{pre}}$
are identical except for a rotation by
$\pi /2$
. Apart from that, the snapshots underline that the different phases the two neighbouring streaks undergo during a cycle of the pre-PO on each of the two sidewalls are arguably the same as those observed for the streaks along the sidewall of case
${\textrm{PO}}{2400}_y$
in figure 8 – enforced by the symmetry constraints to appear simultaneously on both opposing walls.

Figure 10. Time evolution of the r.m.s. energy signal
${{E}_{\textit{rms}}}(t)$
(upper plots) and individual contributions of the triangular sectors to the overall mean streamwise enstrophy
$S_i(t)$
(lower plots) for periodic edge states in the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspace: (a)
${\textrm{PO}}{1600}_{yz}$
, (b)
${\textrm{PO}}{1800}_{yz}$
, (c)
${\textrm{PO}}{2000}_{yz}$
and (d)
${\textrm{PO}}{3200}_{yz}$
. Colour coding in the lower plots is:
$i=\{W,E\}$
(red, identical due to symmetry),
$i=\{S,N\}$
(orange, identical due to symmetry).

Figure 11. Variation of the streamwise-averaged wall shear stress over a full period of case
${\textrm{PO}}{2000}_{yz}$
: (a)
$\langle {{\tau _{xz}}} \rangle _{x}$
along the left/right sidewalls (
$z=\pm {{H}}$
) and (b)
$\langle {{\tau _{xy}}} \rangle _{x}$
along the top/bottom walls (
$y=\pm {{H}}$
). (c) Temporal variation of the perimeter-averaged friction Reynolds number
${\textit{Re}}_\tau$
over the same time window. Contour levels and normalisation are identical to those in figure 9.

Figure 12. Snapshots of the streamwise-averaged primary and secondary flow fields of case
${\textrm{PO}}{2000}_{yz}$
, extracted at different instances along the orbit as indicated by the markers in figure 10
c: (a–f)
$t/{T_b}=\{0.5, 37.7, 88.2, 138.2, 178.6, 239.1 \}$
. Isovalues are the same as in figure 1.
It is interesting to see that a reduction of the Reynolds number below
$2000$
breaks the discussed
${\boldsymbol{R}}_{\pi /2}$
-symmetry between the two pre-periodic segments: in case
${\textrm{PO}}{1800}_{yz}$
(cf. figure 10
b), the two energetic peaks associated with the different pairs of opposite walls differ slightly in amplitude, and the dynamics during the decline of the energy after the burst is somewhat different too. Analogue symmetry breaks have been reported for edge states in the asymptotic suction boundary layer by Kreilos et al. (Reference Kreilos, Veble, Schneider and Eckhardt2013), in their case, however, induced by a variation of the domain size. At even lower Reynolds numbers (cf. case
${\textrm{PO}}{1600}_{yz}$
, figure 10
a), the energy signal implies that a period-doubling sequence might occur for
${{\textit{Re}}_b}\lesssim 1800$
: a full cycle period now features eight individual bursting episodes, four associated with the ‘top/bottom-wall states’ and four with the ‘sidewall states’.
To completely clarify whether such period-doubling or period-multiplication sequences take place in the intervals
${{\textit{Re}}_b}\in (2400,3200)$
and
${{\textit{Re}}_b}\in (1600,1800)$
for the
${\boldsymbol{Z}}_y$
- and the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspaces, respectively, many more edge-tracking runs have to be launched to resolve the ‘continuation’ in Reynolds number – an undertaking that is out of the scope of the current study. Nonetheless, the current results indicate that the variation of the edge-state characteristics with the Reynolds number is for both subspaces highly non-uniform across the parameter space: figure 13(a,b) visualises the edge-state periods as a function of the Reynolds number, implying that the period varies essentially linearly in some parameter ranges. In others, however, a similar change in Reynolds number causes abrupt increases of the limit cycle’s period, by factor four and six in the
${\boldsymbol{Z}}_y$
- and the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspaces (as indicated by the open symbols in figure 13), respectively. Such sudden changes of the edge-state characteristics are in qualitative agreement with the findings of Zammert & Eckhardt (Reference Zammert and Eckhardt2014), who reported alternating regimes of periodic and chaotic edge states, separated by intervals of period multiplication sequences, as the Reynolds number increases in plane channel flow. Khapko et al. (Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014) observed a similar variation of the edge-state properties in the asymptotic suction boundary layer when varying the spatial period of their computational domain. To clarify whether the current variations of the edge-state period are of similar nature, however, requires further efforts and goes beyond the scope of the current study.

Figure 13. Edge-state periods (closed squares) as a function of the bulk Reynolds number for (a) the
${\boldsymbol{Z}}_y$
-symmetric and (b) the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspaces. For limit cycles with long periods, additional open symbols indicate the time separation to other local maxima of the auto-correlation function associated with the signal
${{\textit{Id}}_{\triangle }}(t)$
.
5. Discussion
5.1. Comparison with marginal duct turbulence
In figure 14, we compare the perturbation energy
${E}_{\textit{rms}}$
and the mean secondary flow intensity
${u_{\perp }}=\sqrt {2{{E_{\textit{rolls}}}}}$
of the detected edge states with the corresponding values in fully turbulent trajectories and those for the ‘four-’ and ‘eight-vortex state’ TW families of Okino et al. (Reference Okino, Nagata, Wedin and Bottaro2010) and Uhlmann et al. (Reference Uhlmann, Kawahara and Pinelli2010), respectively, in a domain of identical length. The key difference between the two measures is that
${E}_{\textit{rms}}$
quantifies the energy in the streamwise-varying modes of all three velocity components, whereas
$u_{\perp }$
measures the streamwise-constant modes of the cross-stream velocity components only. From figure 14, it can be seen that the period-averaged values for both quantities in the newly detected edge states are comparable to those of the TWs found by Uhlmann et al. (Reference Uhlmann, Kawahara and Pinelli2010) and Okino et al. (Reference Okino, Nagata, Wedin and Bottaro2010). As can be expected from the number of mean secondary flow cells hosted by each of the TW solutions, the mean perturbation energy of the doubly symmetric edge states agrees fairly well with the ‘eight-vortex state’ of Uhlmann et al. (Reference Uhlmann, Kawahara and Pinelli2010), while that of the single-symmetric edge states attains comparable values as the ‘four-vortex state’ found by Okino et al. (Reference Okino, Nagata, Wedin and Bottaro2010). Instantaneously, however, the peak values of both quantities clearly exceed the corresponding period-averaged values and those attained by the TW solutions. While the maxima of the perturbation energy
${E}_{\textit{rms}}$
along a periodic cycle remain clearly below the level of the turbulent long-time averages, the largest values of the mean secondary flow intensity
$u_{\perp }$
are of similar amplitude, sometimes even exceeding the turbulent long-time averages, in particular for the doubly symmetric edge states. The fact that the doubly-symmetric edge state features the highest mean secondary flow strength is not surprising, as the twofold symmetry enforces vortices to appear near two of the four surrounding walls during the entire cycle, while the chaotic edge state features vortical activity near a single wall most of the time.

Figure 14. (a) The r.m.s. energy signal
${{E}_{\textit{rms}}}/{u_b}$
and (b) mean secondary flow intensity
${u_{\perp }}/{u_b}$
of the edge states as a function of the bulk Reynolds number
${\textit{Re}}_b$
. The vertical ‘error bars’ visualise the range of values attained during one cycle of the periodic cases or along the entire chaotic edge trajectory. Time averages over a full cycle (periodic) or the full trajectory (chaotic) are indicated by a short thick horizontal line along the error bar. The black and blue dashed lines indicate the TW families found by Uhlmann et al. (Reference Uhlmann, Kawahara and Pinelli2010) (‘eight-vortex state’) and Okino et al. (Reference Okino, Nagata, Wedin and Bottaro2010) (‘four-vortex state’), respectively, for the same streamwise period
${\alpha }{{H}}=1$
. Black circles indicate long-time averages of turbulent trajectories (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007). Note that data points for
${\boldsymbol{Z}}_y$
- and
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric cases have been shifted to
${{\textit{Re}}_b}\pm 20$
, respectively, to improve their visibility.
A closer look at the temporal variation of energy and dissipation is provided in figure 15, where the edge dynamics is presented as low-dimensional projections onto planes spanned by
${I}$
and
${D}$
as well as by
${\textit{Id}}_{\triangle }$
and
${D}$
. Recalling that edge states, together with their stable and unstable manifolds, act as mediators between the laminar and turbulent basins in phase space (Toh & Itano Reference Toh and Itano2003), they are well known to feature relevant similarities with weakly turbulent states on the edge to relaminarisation concerning both their organisation in physical space and their flow dynamics in general (Mellibovsky et al. Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009). In figure 15, the edge states are therefore contrasted with the joint probability density function (jpdf) for the corresponding quantities in the marginally turbulent reference simulation at
${{\textit{Re}}_b}=1150$
, for which selected snapshots have been presented in figure 1. The Reynolds number considered for the reference trajectory is slightly higher than the critical value
${\textit{Re}}_b^c\approx 1077$
(Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007) and thus just large enough to sustain turbulence over time intervals of
$\mathcal{O}(10^4{T_b})$
length or longer. In contrast to the fully-developed turbulent case, however, not all four walls are able to host a turbulent regeneration cycle simultaneously most of the time, as has been seen in figure 1.

Figure 15. State space projections onto the planes spanned by (a–c) the dissipation
${D}$
and the energy input
${I}$
and (d–f) the indicator function
${\textit{Id}}_{\triangle }$
and the dissipation
${D}$
. Energy input and dissipation are normalised with the respective values in a laminar flow with identical mass flow rate. The coloured curves in (a,d) represent the dynamics along the chaotic edge states of cases
${\textrm{CA}}{2000}_{2\pi }$
(solid) and
${\textrm{CA}}{2000}_{4\pi }$
(dashed). The remaining plots visualise the dynamics of the periodic edge states in (b,e) the
${\boldsymbol{Z}}_y$
-symmetric subspace (
${{\textrm{PO}}{1600}_y}, {{\textrm{PO}}{2400}_y}, {{\textrm{PO}}{3200}_y}$
as solid, dashed, dotted lines) and (c,f) the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspace (
${{\textrm{PO}}{1800}_{yz}}, {{\textrm{PO}}{2600}_{yz}}, {{\textrm{PO}}{3200}_{yz}}$
as solid, dashed, dotted lines). The grey-coloured background maps represent joint probability density functions of the same quantities in the turbulent reference simulation at
${{\textit{Re}}_b}=1150$
, with contours enclosing
$90\,\%$
and
$50\,\%$
of their total masses. The white-filled squares in (a–c) represent the two flow states in figure 1
a,b, with the lower-dissipation state being the ‘two-vortex state’ in figure 1
a.
The energy input-dissipation diagrams depicted in figure 15(a–c) highlight that the edge states reside most of the time in the rather quiescent region of the state space in between the laminar equilibrium
$({{I}_{\textit{lam}}},{{D}_{\textit{lam}}})$
and what is the turbulent regime for
${{\textit{Re}}_b}=1150$
. For all edge states, the flow orbits along the cycles in a clockwise sense, implying that first the dissipation (or equivalently the enstrophy) experiences a steep rise owing to the strengthening of the quasi-streamwise vortices. After a small time lag of
$\mathcal{O}(1{T_b})$
, the former low-speed streak turns into an intense high-speed region, causing the wall shear stress and thus the external energy input
${I}$
to increase abruptly. This is consistent with the different phases of the self-sustaining process discussed with the aid of the energy decomposition in (4.4) (Waleffe Reference Waleffe1997), and is in line with the edge-state dynamics in plane Poiseuille flow (Toh & Itano Reference Toh and Itano2003) and the asymptotic suction boundary layer (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013).
Also, figure 15(c,f) provides further evidence for the symmetry break occurring for
${{\textit{Re}}_b}\leqslant 1800$
within the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric edge-state family: two distinct loops are discernible at the lower Reynolds number that collapse onto each other as the Reynolds number increases, reflecting that the full orbit is composed of two pre-POs that differ by a rotation of
$\pi /2$
about the centreline. During the bursting excursions, the chaotic edge states touch the lower end of the turbulent regime indicating the chaotic saddle, while their periodic counterparts in the symmetric subspaces reach much further into the turbulent ‘regime’. This is expected as the imposed symmetries enforce turbulent activity to appear on two opposing walls for a subinterval of the cycle (in the
${\boldsymbol{Z}}_y$
-symmetric subspace) or even for the full period (in the
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspace). The volume-averaged dissipation and energy input should therefore attain higher values than in the symmetrically unconstrained case, where the vortical activity is – most of the time – of relevant amplitude near a single wall only.
For the sake of comparison, the parameter points associated with the ‘two-vortex’ state and the ‘four-vortex’ state presented in figure 1(a,b) are marked with white squares in figure 15(a–c), illustrating that the energy input and dissipation is in these cases comparable to the values attained by some of the detected edge states. Such ‘two-vortex’ states, characterised by a strong concentration of the turbulent activity to a single wall, are in fact not that rare for flows under marginally turbulent conditions: for
$20\,\%$
of the long-time turbulent trajectory at
${{\textit{Re}}_b}=1150$
, more than
$50\,\%$
of the total streamwise-averaged enstrophy
$\sum _i S_i$
resides in a single triangular sector. The qualitative comparison between a ‘two-vortex’ configuration of the chaotic edge state and one extracted from the chaotic marginally turbulent trajectory in figure 16 highlights that the general flow configuration of a single low-velocity streak, surrounded by several counter-rotating quasi-streamwise vortices, along a single wall is comparable in both situations. The difference lies in the amplitude of the vortices and thus the intensity of the generated secondary flow field, as well as in the streamwise length of the coherent structures: while the edge state features a single pair of streamwise elongated vortices, their counterparts in the turbulent state have a more complex shape and stretch less in the mean flow direction.

Figure 16. Comparison of instantaneous ‘two-vortex’ states in (a) the chaotic edge state
${\textrm{CA}}{2000}_{2\pi }$
and (b) the marginally turbulent state at
${{\textit{Re}}_b}=1150$
. Left column: streamwise-averaged velocity fields
$\langle {{\boldsymbol{u}}} \rangle _{x}$
(iso-levels as in previous plots). The arrows in (a) are scaled by a factor of two compared with those in (b). Right column: three-dimensional visualisation of the velocity field. The green iso-surface marks
$u=0.75{u_b}$
for all
$y\leqslant -z$
. Vortices are identified as (a)
$Q^+ \gt 0.002$
and (b)
$Q^+ \gt 0.004$
, respectively. Clockwise (red) and counter-clockwise rotation (blue) is measured by the local sign of
$\omega _{x}$
. To facilitate the comparison, the flow state in (b) has been rotated by
$\pi /2$
about the
$x$
-axis.
The fact that the turbulent trajectory transiently visits states in which a single wall/a pair of opposing walls dominates
$\sum _i S_i$
is confirmed by the jpdfs presented in figure 15(d–f): the turbulent trajectory stays most of its lifetime at dissipation levels exceeding
$1.6{{D}_{\textit{lam}}}$
, featuring a vortical activity that is with
${\vert {{{\textit{Id}}_{\triangle }}} \vert }\lesssim 0.5$
rather uniformly distributed among the triangular sectors. However, the presented data also indicate that the turbulent trajectory spends a non-negligible portion of its lifetime in a low-dissipation regime too. The latter is associated with a stronger focussing of the turbulent activity to a single wall/a pair of opposing walls, and the detected edge states are seen to oscillate between these two low-dissipation ‘tails’ with indicator values
${\vert {{{\textit{Id}}_{\triangle }}} \vert }\gt 0.5$
. The marginally turbulent trajectory, on the other hand, does not switch directly from a state with
${{\textit{Id}}_{\triangle }}\gtrsim 0.5$
to one with
${{\textit{Id}}_{\triangle }}\lesssim 0.5$
or vice versa, but it usually spends a (potentially long) while in the higher-dissipation region first (cf. also Sekimoto Reference Sekimoto2011).
Our findings thus strongly suggest that the low-dissipation excursions of a weakly turbulent trajectory to localised flow states are transient visits to one of the discussed edge states. The trajectories shadow the edge-state dynamics for a certain time, potentially undergoing a bursting and wall-switching episode. Eventually, however, they get repelled from the edge along its single unstable direction back to the higher-dissipation regime of phase space, where flow states feature a significant vortical activity all over the cross-section.
5.2. The edge state in state space
In their seminal work on the edge state of plane Poiseuille flow, Toh & Itano (Reference Toh and Itano2003) observed a periodic-like dynamics of alternating quiescent phases and energetic bursting episodes along the edge to turbulence, but their limited computational resources did not allow them to track the edge trajectory longer than one and a half cycles. So, Toh & Itano (Reference Toh and Itano2003) were raising the question whether the state space object they had identified was a PO or a heteroclinic cycle. In contrast to a PO, a heteroclinic cycle consists of two saddles that are mutually connected in phase space by two heteroclinic orbits (cf. Guckenheimer & Holmes Reference Guckenheimer and Holmes1983, p. 46 ff and Wiggins Reference Wiggins2003, p. 631 ff).
While Zammert & Eckhardt (Reference Zammert and Eckhardt2014) later re-analysed the edge for plane Poiseuille flow and confirmed that the edge state is indeed a PO, heteroclinic cycles were for a long time interpreted as the state space equivalent of turbulent bursts: Aubry et al. (Reference Aubry, Holmes, Lumley and Stone1988) found such an attracting heteroclinic cycle in a low-dimensional model of the near-wall region of a turbulent boundary layer. For most of their lifetime, trajectories were seen to stay near one of the two weakly unstable fixed points – representing relatively quiescent episodes of the turbulent activity. These episodes are interrupted when the trajectories eventually leave the fixed points along their unstable manifolds, following the heteroclinic orbit to its sibling. In physical space, the dynamics along the heteroclinic orbits represents intense bursting events. In a follow-up study, Armbruster, Guckenheimer & Holmes (Reference Armbruster, Guckenheimer and Holmes1988) showed that it is in fact the
$O(2)$
-symmetry of the system that makes structurally stable heteroclinic cycles possible. In the following years, similar statements were made for other systems obeying a
$O(3)$
-symmetry (Guckenheimer & Holmes Reference Guckenheimer and Holmes1988) and even discrete
$D_n$
-symmetries (Campbell & Holmes Reference Campbell and Holmes1992), underlining the high relevance of the system’s symmetries for the existence of simple invariant sets in state space.
For the square duct, we have provided clear evidence that the symmetric edge states are limit cycles within the edge, and attempts to identify weakly unstable fixed points along the orbit were unsuccessful. In this regard, our results are in good accordance with findings in
$O(2)$
-symmetric flow systems like plane Poiseuille flows (Neelavara et al. Reference Neelavara, Duguet and Lusseyran2017; Rawat et al. Reference Rawat, Cossu and Rincon2014, Reference Rawat, Cossu and Rincon2016; Zammert & Eckhardt Reference Zammert and Eckhardt2014) and the asymptotic suction boundary layer (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013), whose edge states are all POs rather than heteroclinic cycles between two saddles.
Even though simple heteroclinic cycles could not be found within
${\mathscr{M}}$
for the discussed flows, the systems’ symmetries still play a major role for the structure of the edge: topologically more complicated heteroclinic networks between pairs of symmetry-related saddles and stable fixed points (here the edge states) turn out to occur in many systems and can be dynamically linked to periodic edge states. So, Kreilos et al. (Reference Kreilos, Veble, Schneider and Eckhardt2013) showed with a homotopy from plane Couette flow to the asymptotic suction boundary layer – with the suction velocity at the wall as homotopy parameter – that time-periodic edge states arise in a saddle-node infinite-period (SNIPER) bifurcation (Tuckerman & Barkley Reference Tuckerman and Barkley1988; Strogatz Reference Strogatz2018, p. 265f; Tuckerman Reference Tuckerman2020). Under subcritical conditions, the edge features two symmetry-related edge states and two saddles separating them, all mutually connected by heteroclinic orbits. At the critical point, the stable fixed points collide with the saddle points and give rise to a PO whose period is infinite at the bifurcation point, but decreases rapidly when moving away from it. Even for clearly supercritical conditions, the PO still features some properties of the previously present invariant structures in that the dynamics slows down significantly when the trajectory passes by the ‘ghosts’ of the fixed points in phase space. A similar SNIPER bifurcation gives rise to the periodic edge states in plane Poiseuille flow under a variation of the lateral spatial domain size (Rawat, Cossu & Rincon Reference Rawat, Cossu and Rincon2016). In the light of the qualitatively similar dynamics and taking into account the system’s symmetries, we suspect that the periodic edge states in the square duct might also arise in a structurally similar SNIPER bifurcation. In analogy to the situation in plane channel flow (Rawat et al. Reference Rawat, Cossu and Rincon2016), a likely candidate for a parameter that drives the system into a SNIPER bifurcation is the cross-sectional aspect ratio
${A}={L_z}/{L_y}$
. In order to check for the existence of such a bifurcation, however, additional efforts in form of a thorough homotopy study with
$A$
as homotopy parameter are required.
In a more general sense, our findings imply – consistent with the mentioned studies on plane Poiseuille and asymptotic suction boundary layer flows – that the observed bursting dynamics occurs when a system follows connecting orbits between one or more low-dissipation invariant solutions, or the segment of a PO that is dynamically relatable to such a homoclinic or heteroclinic orbit. This is in line with the observations of van Veen & Kawahara (Reference van Veen and Kawahara2011) who argued for Couette flow that the bursting dynamics might represent the physical incarnation of a homoclinic cycle around a single ‘gentle’ PO. The latter is an alternative to earlier considerations by Kawahara & Kida (Reference Kawahara and Kida2001) who proposed a heteroclinic cycle between a ‘gentle’ low-dissipation periodic edge state and a more ‘vigorous’ high-dissipation PO in the turbulent attractor as state space equivalent of bursting events in physical space.
In this context, it is worth mentioning that bursting in all detected edge states (chaotic and time periodic) goes hand in hand with a switching of the turbulent activity to one or both neighbouring walls. Unless the dynamics is further restricted to a higher-symmetric subspace in which states merely differing by a
$n\pi /2$
-fold rotation about the
$x$
-axis (
$n\in {\mathbb{N}}$
) collapse onto each other, the qualitative state space picture is that of two (or more) heteroclinically connected states, or a related PO. A homoclinic orbit, on the other hand, would imply a regeneration cycle that remains attached to the same wall after the burst, not experiencing the ‘wall-switching’ dynamics visible in the here detected edge states. While such dynamics is generally conceivable, we have not observed it in any of our edge-tracking runs (neither in the unconstrained state space nor in its symmetric subspaces).
6. Conclusion
In this study, we have investigated the dynamics of trajectories within the edge
${\mathscr{M}}$
to square duct turbulence for both the full state space and symmetric subspaces of the former. Considering the full state space without symmetry constraints, the edge state is a chaotic attractor within
${\mathscr{M}}$
for the analysed parameter regime. The results are consistent with the chaotic edge states reported for circular pipe flows (Duguet et al. Reference Duguet, Willis and Kerswell2008), which feature many similarities to the dynamics in the square duct. In contrast to their counterparts in the pipe, however, the edge states in the square duct experience phases of strong intermittent bursts that interrupt the otherwise rather quiescent dynamics. In physical space, these bursting episodes are seen to represent the different stages of the well-known self-sustaining near-wall cycle (Waleffe Reference Waleffe1997; Hall & Sherwin Reference Hall and Sherwin2010): the quiescent phases are characterised by a single straight low-speed streak that stretches along one of the four surrounding walls, flanked by counter-rotating quasi-streamwise vortices on either side. The initially straight streak undergoes a linear instability and starts bending laterally, while the vortices intensify and tilt with respect to the streamwise direction. Eventually, the vortices lean over the low-speed streak, tear it apart and induce a pronounced high-speed zone instead. Due to the particular geometry of the square duct, however, the intense vortices are located in the corner regions of the cross-section and interact simultaneously with the near-wall region of both adjacent walls. In other words, the large corner vortices simultaneously push high-momentum fluid towards one wall and induce a low-speed streak at the other. The new low-speed streak at the neighbouring wall slowly starts bending and the bursting cycle starts anew. With a length of
$\mathcal{O}(100){T_b}$
for both chaotic and time-periodic edge states, the characteristic time scale of the bursting events is much longer than that associated with their downstream propagation.
The outlined bursting cycles reveal a strong dynamical similarity to the periodic and chaotic edge states in plane Poiseuille flow (Rawat et al. Reference Rawat, Cossu and Rincon2014; Zammert & Eckhardt Reference Zammert and Eckhardt2014) and the asymptotic suction boundary layer (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013; Khapko et al. Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014). But while the tilted vortices in these ‘planar’ spanwise-periodic flow configurations give birth to a new low-speed streak shifted by (approximately) half the spanwise period
${L_z}/2$
to the left or right, here, the vortices induce a new low-speed streak at a position rotated by
$\pi /2$
or
$3\pi /2$
about the duct centreline. The structural similarity of the edge dynamics in both systems accentuates the role of the systems’ symmetries that make such switching dynamics possible in the first place – a continuous
$O(2)$
-symmetry in Poiseuille flow or the asymptotic suction boundary layer, and a discrete
$D_n$
-symmetry in the square duct.
In the ‘planar’ spanwise-periodic cases, edge trajectories converge for a wide range of Reynolds numbers and domain sizes to periodic edge states that are limit cycles within
${\mathscr{M}}$
. For the square duct, we detect such (relative) periodic edge states merely when restricting trajectories to distinct symmetric subspaces, specifically focussing on states obeying a single or twofold mirror symmetry with respect to the duct’s wall bisectors. To the best of the authors’ knowledge, the reported edge states are the first POs found in a duct with rectangular cross-section. The edge states in both symmetric subspaces feature a very similar dynamics to their chaotic siblings in the unconstrained case – with the difference that the imposed symmetries partly enforce bursting cycles to appear in symmetric pairs on two opposing walls. A systematic variation of the Reynolds number confirmed the robustness of the detected periodic edge-state dynamics over a significant range of Reynolds numbers, but it also revealed that the specific edge-state characteristics change with the governing parameters in a highly non-uniform way – in accordance with observations in ‘planar’ spanwise-periodic shear flows (Zammert & Eckhardt Reference Zammert and Eckhardt2014). While the edge-state period varies almost linearly with
${\textit{Re}}_b$
in some parameter regimes, sudden steep rises of the period in terms of a period doubling or a symmetry break are observed elsewhere in the parameter space.
The current results indicate that the localisation of a single turbulent regeneration cycle to one or two of the surrounding walls is a fundamental feature of the edge-state dynamics in square duct flow. Also, it was revealed that marginally turbulent trajectories regularly attain similar low-dissipation states, in which turbulent perturbations are concentrated near one or two walls only. Together, these observations raise the question whether such ‘two-vortex’ and ‘four-vortex’ episodes along marginally turbulent trajectories can be associated with transient visits to the edge states detected in the course of this study. Within this dynamical systems perspective, the trajectories are expected to approach an edge state via its stable manifold that locally collapses with
${\mathscr{M}}$
, to shadow its dynamics for a certain time interval and to eventually get repelled along the single unstable direction towards the turbulent saddle/attractor. In order to check the validity of this hypothesis, further efforts are necessary to identify such excursions to the vicinity of an edge state along a marginally turbulent trajectory and to verify that these are indeed accompanied by a localisation of the turbulent activity to one or two walls.
More generally, it remains to clarify how the edge-state dynamics changes when the cross-sectional aspect ratio
$A$
is increased towards rectangular ducts of non-square cross-section. Although some preliminary results prove the general existence of spanwise-localised TW solutions in the rectangular case
${A}\gt 1$
(Okino Reference Okino2011), a comprehensive investigation of the edge dynamics that could clarify how the intermittent bursting dynamics and the wall-switching behaviour change with
$A$
is missing. With the observation that periodic edge states can arise in the context of SNIPER bifurcations in mind (Kreilos et al. Reference Kreilos, Veble, Schneider and Eckhardt2013; Rawat et al. Reference Rawat, Cossu and Rincon2016), a comparison of the edge-state dynamics at successively increasing aspect ratios will help verifying whether or not the here identified periodic edge states are born in a similar global bifurcation. Recalling that plane Poiseuille flow represents the asymptotic state of an infinitely wide duct as
${A}\to \infty$
, we expect such a continuation in the duct’s aspect ratio to moreover provide valuable information on how edge states in the rectangular duct are dynamically connected to those in plane channel flow.
Acknowledgements
We want to express our thanks to Professor J. Jiménez and Dr K. Osawa for their kind hospitality during the fifth Madrid turbulence summer workshop. We also thank Professor M. Avila, Dr D. Morón and Professor L. van Veen for very helpful discussions.
Funding
This work was supported in part by the European Research Council under the Caust grant ERC-AdG-101018287. Part of the simulations were performed on the supercomputer bwUniCluster funded by the Ministry of Science, Research and the Arts Baden-Württemberg and the Universities of the State of Baden-Württemberg. The computer resources, technical expertise and assistance provided by the staff are gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Convergence to a periodic edge state
The time signal
${{\textit{Id}}_{\triangle }}(t)$
is used to estimate the fundamental period of the shadowed limit cycle, since it contains more information on the dynamics than the ‘pure’ energy signal
${E}_{\textit{rms}}$
. In particular, it distinguishes different states that are identical modulo finite rotations. To check for convergence, an integral
$L_2$
-error between two consecutive cycles
${\textit{Id}}_{\triangle }^{\;k-1}$
and
${\textit{Id}}_{\triangle }^{\;k}$
of the expected period length
$T$
is introduced as
\begin{equation} \varepsilon = \dfrac {{\big \lVert {{\textit{Id}}_{\triangle }^{\;k}-{\textit{Id}}_{\triangle }^{\;k-1}} \big \rVert _{2,T}}}{{\big \lVert {{\textit{Id}}_{\triangle }^{\;k}} \big \rVert _{2,T}}}, \end{equation}
where the
$L_2$
-norm over a full (guessed) period is defined as
\begin{equation} {\left \lVert {\bullet } \right \rVert} _{2,T}^2 = \displaystyle \int \limits _{t^{\prime}=t}^{t+T}\; \bullet ^{2} \;{\textrm{d}}t^{\prime}. \end{equation}
As period guess
$T$
, we choose the time at which the temporal auto-correlation function of
${{\textit{Id}}_{\triangle }}(t)$
(excluding an initial transient) attains its non-trivial global maximum. An edge trajectory is then considered converged to a PO if the error
$\varepsilon$
defined in (A1) falls below
$5\,\%$
. The only exception is case
${\textrm{PO}}{3200}_y$
, for which the two pre-POs feature the same
${{\textit{Id}}_{\triangle }}(t)$
profile. In this case, we quantify convergence instead based on the signal of the streamwise-averaged enstrophy in the single sectors,
$S_i(t)$
.
Appendix B. Stability of the periodic edge state
That the shadowed POs are indeed limit cycles is shown in figure 17, where a first-return map is presented for the peaks of the signal
${{\textit{Id}}_{\triangle }}(t)$
for selected POs in both considered symmetric subspaces (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Zammert & Eckhardt Reference Zammert and Eckhardt2014). Following classical dynamical systems theory (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983, p. 22), this low-dimensional projection of the full system is strictly periodic if its Poincaré map converges to a fixed point. This is the case if the sequence in the direct neighbourhood of the fixed point follows a straight line with slope
$|\gamma | \lt 1$
. The black lines in figure 17 represent the linear trend for the last few iterations and verify that this condition is indeed fulfilled for the selected cases.

Figure 17. First-return map of the indicator function
${\textit{Id}}_{\triangle }$
in terms of the peak values of the signal for selected POs in the (a)
${\boldsymbol{Z}}_y$
- and (b)
${{\boldsymbol{Z}}_y}{{\boldsymbol{Z}}_z}$
-symmetric subspaces. Thin black lines represent the best fit linear approximation to the map
${{\textit{Id}}}_{max }^{\;k+1} = {\textrm{P}}({{\textit{Id}}}_{max }^{\;k})$
, computed for the last three points in the sequence.





































































































































