2. Problem statement

We consider an incompressible Newtonian fluid flow through a square duct with walls located at cross-axial positions $x = \pm h$ and $y = \pm h$ , and periodic in the axial $z$ -direction with period $L_z$ . Unless stated otherwise, $L_z=8\pi\! h$ , which is close to the minimal streamwise length suitable for supporting the previously reported localised solution (Okino Reference Okino2014; Gepner et al. Reference Gepner, Okino and Kawahara2025). The flow, represented by the velocity vector $\boldsymbol {u} = [u, v, w]$ , is driven by a constant-in-time pressure gradient producing a laminar flow $\boldsymbol {u}_l$ with a steady, quasi-parabolic velocity profile. We use the half-width $h$ of the duct as the characteristic length scale, and the laminar centreline velocity $W$ for the given constant-in-time streamwise pressure gradient as the characteristic velocity scale to non-dimensionalise physical quantities. Taking density $\rho$ as unity and kinematic viscosity $\nu$ , the Reynolds number is defined as $\mathit{Re} = \textit{Wh}/\nu$ and the bulk Reynolds number is $\mathit{Re}_b(t) = w_b(t)h/\nu$ , with $w_b(t) = ({1}/{V})\iiint _V w(t) \, \textrm {d}V$ where $V$ is the domain volume. Due to the applied pressure gradient constraint, $\mathit{Re}_b$ is not a fixed quantity – see table 1.

Table 1. Time-averaged flow quantities and growth rate $\sigma$ of the leading eigenmode (if applicable) for different types of solutions at $\mathit{Re}=4000$ . Here, $\langle \boldsymbol{\cdot }\rangle$ indicates a time-averaged value.

The system is studied via the SEM implemented in Nektar $\scriptscriptstyle{++}$ (Cantwell et al. Reference Cantwell2015; Moxey et al. Reference Moxey2020). The numerical resolution matches that used by Gepner et al. (Reference Gepner, Okino and Kawahara2025), comprising four quadrilateral elements in the $(x, y)$ -plane, each with polynomial order $p=24$ and $q=36$ Gauss–Lobatto–Legendre quadrature points per direction, satisfying the $q = 3/2 p$ integration rule for exact nonlinear term quadrature (Karniadakis et al. Reference Karniadakis and Sherwin2005). In the streamwise $z$ -direction, a Fourier discretisation is employed, truncated to $M = \pm 128$ modes prior to applying the $2/3$ dealiasing rule (Patterson Jr & Orszag Reference Patterson and Orszag1971). The resulting spatial resolution corresponds to approximately $1.8 \times 10^6$ degrees of freedom (DoF) for velocity components. Robustness is confirmed by increasing the resolution to $p=30$ , $q=45$ and $M = \pm 144$ ( $\sim 3.1 \times 10^6$ DoF), which yields consistent invariant solution characteristics and stability properties.

3. Edge tracking

Perturbing the laminar flow with a short (approximately $2$ time units) low-intensity white-noise body force (volume forcing) applied to the momentum equations induces a non-laminar dynamics. Depending on the Reynolds number and the disturbance amplitude, this non-laminar response manifests itself either as a chaotic transient, eventually relaminarising, or as sustained turbulence. At $\mathit{Re} = 4000$ , we observe a persistent non-laminar state; throughout long-time simulations exceeding $10^5$ time units, no return to laminar flow is detected. Conversely, below approximately $\mathit{Re} = 3600$ , only transient non-laminar states arise, with relaminarisation typically occurring within $10^4$ time units. This suggests the presence of a critical Reynolds number for sustained turbulence in the square duct, akin to previous observations in pipe flow (Avila et al. Reference Avila, Moxey, De Lozar, Avila, Barkley and Hof2011). Below this threshold, turbulence is transient (a chaotic saddle dynamics), and the notion of the edge should be understood as the codimension-one boundary separating initial conditions that relaminarise directly from those that first undergo chaotic transients (Chantry et al. Reference Chantry, Willis and Kerswell2014). For subsequent edge-tracking analysis, we fix the Reynolds number at $\mathit{Re} = 4000$ , corresponding to a friction Reynolds number $\mathit{Re}_{\tau } \approx 82$ , and the laminar and time-averaged turbulent bulk Reynolds numbers are equal to $\mathit{Re}_b=1908$ and $\langle \mathit{Re}_b \rangle =1136$ , respectively, with $\langle \boldsymbol{\cdot }\rangle$ indicating a time-averaged value (see table 1).

We initiate edge tracking with a snapshot $\boldsymbol {u}_s$ drawn from fully developed turbulent flow. This snapshot is taken sufficiently far from the initial perturbation ( $2000$ time units later or even longer) to ensure it originates from the turbulent attractor rather than the transient phase. No additional criteria are imposed in selecting the initial condition for edge tracking. We then apply the classical bisection algorithm (Itano & Toh Reference Itano and Toh2001; Skufca et al. Reference Skufca, Yorke and Eckhardt2006) to find initial states $ \boldsymbol {u}_i = \boldsymbol {u}_l + \alpha (\boldsymbol {u}_s - \boldsymbol {u}_l),$ varying $\alpha \in (0,1)$ at time $t=0$ such that the flow neither relaminarises nor transitions fully to turbulence but remains on the laminar–turbulent boundary, the edge, separating the two attraction basins. During the bisection process, both the perturbation energy $ E_{3D} = ({1}/{2V}) \sum _{k=1}^M \iiint _V \boldsymbol {u}_{-k} \boldsymbol{\cdot }\boldsymbol {u}_k \, \textrm {d}V$ and the friction factor $ f = 8 ({\mathit{Re}_\tau }/{\mathit{Re}_b} )^2$ are monitored and gradually decrease from their turbulent values, while exhibiting recurrent fluctuations.

Figure 1. Snapshots of the velocity field depicting (a,c) the identified edge state – LB solution – and (b,d) the upper branch (UB) solution at $\mathit{Re}=4000$ . Panels (a,b) show isosurfaces of the second invariant of the velocity gradient tensor $Q$ at $0.01$ , and panels (c,d) display contours of the difference in the streamwise-velocity component $w$ relative to the laminar solution $w_l$ , with negative values indicated by dashed lines; all slices are taken at $z=0$ . Velocity fields are adjusted such that the maximum velocity perturbation aligns with $z=0$ .

This bisection process has been implemented starting from various initial snapshots, with some leading the procedure to converge, after an initial transient, to a regular state with a non-zero asymptote of perturbation energy. The velocity field obtained via bisection was then used as an initial guess for a custom implementation of the Newton–Krylov solver (Viswanath Reference Viswanath2007, Reference Viswanath2009) built around the Nektar $\scriptscriptstyle{++}$ framework. The solver converged to a simple travelling-wave solution. During independent bisections, we observe convergence to the same regular travelling-wave solution, up to a rotation by $\pi /2$ about the duct axis. This is consistent with the $\pi /2$ -rotational equivariance of the square-duct geometry, implying that symmetry-related copies can be located at any wall.

This solution features a streamwise-localised vortex pair situated next to one wall of the duct (here shown at $x=1$ ) that travels downstream at a constant phase speed $w_p=0.695$ . Figure 1(a,c) illustrates the invariant solutions via $ Q = ({1}/{2}) ( \|\varOmega \|^2 - \|S\|^2 )$ contours (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988), with $S$ and $\varOmega$ the symmetric and antisymmetric parts of the velocity gradient tensor, respectively, plotted at $Q = 1.0 \times 10^{-2}$ , and streamwise-velocity perturbations relative to the laminar quasi-parabolic profile on the $z=0$ plane. While the computations are performed in the full phase space without imposing any symmetry constraints, the converged travelling wave exhibits an autonomous reflection symmetry with respect to one of the wall bisectors – here $y=0$ (see figure 1), illustrating that symmetries arise autonomously even in the full, unconstrained phase space.

Figure 2. Variations of (a) the cross-flow energy $E_{\perp }$ and (b) streamwise velocity at the duct centre $w_c$ of the LB solution at $\mathit{Re}=4000$ with streamwise coordinate $z$ . Velocity fields are shifted in the periodic $z$ -direction such that the maximum cross-flow energy is located at $z=0$ .

To investigate streamwise localisation, we continued the solution by increasing the domain length $L_z$ , maintaining $\mathit{Re}=4000$ and adjusting the streamwise resolution to preserve approximately 10 Fourier modes per unit length. Figure 2(a) shows the streamwise variation of the cross-flow energy $ E_{\perp }(z) = ({1}/{4}) \int _{-1}^1 \int _{-1}^1 ({1}/{2}) (u^2 + v^2) \, \textrm {d}x \textrm {d}y$ , and figure 2(b) the centreline streamwise-velocity $w_c$ , for increasing periodic length $L_z$ , illustrating that the cross-flow energy remains confined around a compact core as $L_z$ is increased, whereas the streamwise-velocity perturbation exhibits weaker, non-exponential decay. In this sense, streamwise localisation refers to the energy core (based on $E_{\perp }(z)$ ) rather than implying strict exponential localisation in all velocity components. This trend suggests that, with further increase of the domain length, the solution should persist and remain available for the case of a fully extended system, making it meaningful for experimental investigation (De Lozar et al. Reference De Lozar, Mellibovsky, Avila and Hof2012).

Using the arc-length continuation method (Keller Reference Keller1977; Dijkstra et al. Reference Dijkstra2014), we continued the solution in Reynolds number while keeping the domain length fixed at $L_z = 8\pi$ . This continuation revealed that the identified edge state corresponds to the LB of a saddle-node bifurcation occurring near $\mathit{Re} = 2794.5$ . We further identified the corresponding UB solution of this pair and continued it with increasing Reynolds number. Figure 1(b,d) illustrates the UB topology at $\mathit{Re} = 4000$ . Similarly, by fixing $\mathit{Re} = 4000$ and varying the duct length $L_z$ , we observed an analogous saddle-node bifurcation near $L_z \approx 4.64\pi$ , with an associated UB solution. At $L_z = 8\pi$ , the UB solution matches that shown in figure 1(b,d). These results confirm that the identified solution belongs to a continuous solution branch parametrised by both Reynolds number and domain length.

Stability analysis of the identified solution pair was conducted using Arnoldi iterations (Viswanath Reference Viswanath2009). At $\mathit{Re} = 4000$ and $L_z = 8\pi$ , the LB solution possesses a single purely real unstable eigenvalue (see table 1) in the full phase space, confirming its status as an edge state. Accordingly, this travelling wave represents a full-phase-space edge state whose codimension-one stable manifold locally partitions the phase space into two distinct basins (Duguet et al. Reference Duguet, Willis and Kerswell2008; Gepner et al. Reference Gepner, Okino and Kawahara2025), thus forming the laminar–turbulent edge. Unlike the previously reported symmetric-subspace edge states in square-duct flow (Okino Reference Okino2014; Gepner et al. Reference Gepner, Okino and Kawahara2025) and pipe flow (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013) – both appearing embedded within the turbulent attractor (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013; Gepner et al. Reference Gepner, Okino and Kawahara2025) – the present solution is necessarily separated from the turbulent trajectory, a prerequisite for persistent turbulence (Lustro et al. Reference Lustro, Kawahara, van Veen, Shimizu and Kokubu2019). Extending the stability analysis to longer domains reveals that the identified travelling wave remains an edge state for streamwise lengths of at least $L_z = 32\pi$ ( ${\gt } 100$ ), the longest domain length considered here, with only small variation in the growth rate of the single unstable eigenvalue, suggesting that it should be an edge state in a fully extended system.

4. Edge-state continuation and heteroclinic connection

The performed stability analysis shows that the LB solution remains an edge state, possessing a single unstable eigenvalue across all tested Reynolds numbers up to the bifurcation point at $\mathit{Re} = 2794.5$ . The coexisting UB solution has multiple unstable directions. Past the bifurcation point, for Reynolds numbers in the range $\mathit{Re} \in (2794.5, 2800)$ , the UB solutions feature precisely two real unstable eigenvalues. Beyond $\mathit{Re} = 2800$ , up to at least $\mathit{Re} = 2823$ , the number of real unstable eigenvalues increases to three. Prior to $\mathit{Re} = 2847$ , two of these real unstable eigenvalues transition into a pair of complex conjugate unstable eigenvalues. Additional unstable eigenvalues continue to appear on the UB as the Reynolds number is further increased.

Figure 3. State portrait of the dynamics represented on the $(E_{3D}, f)$ -plane. The green solid curve corresponds to a persistent turbulent trajectory at $\mathit{Re}=4000$ , observed over $10^5$ time units. The grey solid (dashed) line represents the LB (UB) solutions with varying Reynolds number. Solid red (blue) curves start from edge states at different $\mathit{Re}$ and depict excursions toward the turbulent attractor (if one exists) – solid dark red – or the turbulent transient – solid red – (laminarisation) generated from initial conditions formed by perturbing the LB and UB within their respective unstable manifolds. Red and blue curves shown in the inset depict the bisection process started from the UB solution, following the heteroclinic orbit (marked with a thick arrow) connecting the UB to the LB solution at $\mathit{Re}=2796$ .

The single purely real unstable eigenvalue of the LB solution defines a direction necessarily transverse to the edge manifold (Duguet et al. Reference Duguet, Willis and Kerswell2008) at the LB. In contrast, for Reynolds numbers below approximately $2800$ , the UB solution possesses two purely real unstable eigenvalues, which locally define a two-dimensional unstable manifold. We investigate unstable directions of both branches by perturbing the LB solution along its single unstable eigenvector and the UB solutions in a direction that is a combination of its two unstable eigenvectors at $\mathit{Re} = 2796$ and $L_z = 8\pi$ . The resulting trajectories are shown in figure 3 in the $(E_{3D}, f)$ -plane. In each case, the separation of initial conditions is clear between those that relaminarise and those that transition to turbulence. By performing bisection along unstable directions, we track edge trajectories that start near these solutions. For the LB, initial conditions lie increasingly longer near the solution before departing toward either laminar or turbulent states, confirming its role as an edge state. A similar pattern occurs for the UB when perturbed along its leading unstable direction, suggesting this direction remains closely transverse to the edge. At the same time, perturbations that are closely aligned with the secondary unstable eigenvector of the UB, with only a small component along the leading unstable eigenvector (see the inset in figure 3) initially depart but remain confined to the edge. As bisection progresses, trajectories starting at UB steadily approach the LB, as depicted by the a thick arrow in the inset in figure 3. This bisection exhibits recurrent escapes to laminar or turbulent states, consistent with trajectories embedded within the edge manifold. Ultimately, the bisection converges to the LB solution, exhibiting monotonic decay in perturbation energy and friction factor. A final Newton iteration confirms that the approached state is the LB solution and providing sound evidence for a heteroclinic connection embedded in the edge, which links the UB and LB solutions. Consequently, the UB solution must also reside within this edge (but is not an edge state itself).