\section*{Introduction}

    Time irreversibility (TI) is a fundamental property of non-equilibrium systems, which are typically dissipative~\cite{evans1993probability, lebowitz1993boltzmann}. In steady state, TI appears as an asymmetry of the statistical properties of a signal when the time direction is reversed~\cite{cox1991long, Lawrance1991DirectionalityAR, lebowitz1993boltzmann, bertini2001fluctuations, giberti2007temporal}. TI plays an important role in revealing key features of nonlinear dynamical systems, such as long-range nonlinearity and non-Gaussianity, and is inherently associated with entropy production in statistical mechanics and thermodynamics~\cite{cox1991long, porporato2007irreversibility}. Owing to its importance in non-equilibrium systems, TI has been investigated through the lens of time-series analysis in many scientific areas~\cite{costa2005broken, porporato2007irreversibility, zorzetto2018extremes, skinner2021estimating, o2022time}.
    
    Turbulent flows represent a paradigmatic example of dissipative, highly-nonlinear, and far-from-equilibrium systems~\cite{frisch1995turbulence, pope2000turbulent}. A distinctive feature of turbulent flows is the presence of a broad range of scales -- from the largest, inertial, scales to the smallest, dissipative, scales -- across which energy is redistributed through a cascade process~\cite{pope2000turbulent}. While viscous dissipation formally leads to an \textit{intrinsic irreversibility} related to time-irreversible (Navier-Stokes) equations, high-Reynolds-number turbulence (where the large-to-small scale range is wide) is mainly driven by inertial forces that make dissipative effects less evident, leading to a (nearly-)reversible dynamics~\cite{vela2021entropy, vela2022subgrid}. In spite of this, turbulent flows at high Reynolds numbers -- owing to their high complexity -- manifest significant \textit{statistical irreversibility} that is typically associated with the energy cascade process~\cite{frisch1995turbulence, xu2014flight, xu2016lagrangian}. Specifically, it is possible to quantify the energy flux driving the cascade process without accounting for the temporal evolution of the flow, by exploiting Kolmogorov’s 4/5-law that involves the (fixed time) third-order moment of velocity spatial increments~\cite{frisch1995turbulence}. Therefore, a quite paradoxical situation arises, where the energy cascade (entailing an asymmetry of spatial statistics) can be characterized without consideration of time, although the energy flux is itself the cause of statistical TI~\cite{xu2016lagrangian, josserand2017turbulence, Lawrance1991DirectionalityAR}.
    

    \begin{figure*}[t]
        \centering
        \includegraphics[width=0.99\linewidth]{fig1_H-01.eps}
        \caption{\textbf{Summary sketch of wall turbulence and HVGs.} \textbf{a} Schematic of the wall-normal structure of a turbulent boundary layer: (left) vertical profiles of the mean velocity $U^+(y^+)$ (red line, with $U_\infty^+=31.89$) and the root-mean-square velocity $\sigma^+_u(y^+)$ (green line), as a function of the normalized wall-normal coordinate $y^+$ (in a log-scale); (middle) the main vertical layers delimited by conventional $y^+$ limits~\cite{pope2000turbulent}; (right) a qualitative sketch of scale arrangement in the flow. \textbf{b} Sketch of a time-series $u(t_i)$ (black vertical lines), and HVG links (colored arrows) for a node $t_j$. Red-blue number pairs indicate backward and forward degree values, $k_b$ and $k_f$ respectively, for each node illustrated as a black-filled circle on top of each vertical line.}
        \label{fig:fig1}
    \end{figure*}
    
    
    With the aim to shed light on such paradoxical situation and understand the connection between energy cascade-induced spatial asymmetry and time irreversibility in turbulent flows, several studies have been carried out, with a particular focus on homogeneous-isotropic turbulence~\cite{xu2014flight, jucha2014time, buaria2015characteristics, bragg2016forward, cencini2017time} and two-dimensional (2D) turbulence~\cite{rondoni1999fluctuations, gallavotti2004lyapunov, bragg2018irreversibility, porporato2020fluctuation}. In particular, the Lagrangian viewpoint has been adopted as the main framework to investigate TI, by looking at the asymmetry in the statistics of backward and forward dispersion of tracer particles~\cite{xu2014flight, jucha2014time, buaria2015characteristics, bragg2016forward, xu2016lagrangian, cencini2017time, polanco2018relative, zhang2021time, cheminet2022eulerian}. In contrast, the investigation of TI in wall-bounded turbulent flows have received much less attention~\cite{polanco2018relative, zorzetto2018extremes}, although wall turbulence plays a crucial role in several key engineering and geophysical problems~\cite{jimenez2013near, li2019large}. Wall-bounded turbulence includes an additional level of complexity with respect to homogeneous isotropic turbulence due to the inhomogeneity and anisotropy of the flow along the wall-normal direction. This additional peculiarity of wall turbulence implies that statistical irreversibility in wall turbulence can arise, not only via the energy cascade process, but also due to the complex spatio-temporal development of wall-induced coherent motion, which constitutes the backbone of turbulent flows~\cite{jimenez2018coherent}.
    
    
    In order to fill the knowledge gap concerning the connection between TI and coherent motion, in this work we investigate statistical irreversibility in wall-bounded turbulence focusing on how different flow scales contribute to TI. Characteristic flow scales, in fact, are associated with coherent structures, and display a complex dynamics, not only in terms of energy cascade, but also through inter-scale interactions such as modulation phenomena~\cite{baars2015wavelet, iacobello2021large}. Since coherent structures exhibit different features along the wall-normal direction (e.g., in terms of geometry, dynamics, as well as energetic content~\cite{smits2011high, jimenez2018coherent}), we study the effect of the wall-normal coordinate on TI levels, exploring the whole range of wall-normal distances from the near-wall region (where wall-induced effects are dominant) to the outer flow region (where wall-induced effects are negligible). By so doing, we are able to focus on the relation between TI and flow scales while also accounting for the inhomogeneity and anisotropy of wall-bounded turbulence, which affects the spatio-temporal development of coherent structures. To this aim, experimental and numerical data of the streamwise velocity $u$ (which retains fundamental features of wall-bounded turbulence mechanisms~\cite{jimenez2018coherent}) are exploited from both external (turbulent boundary layer) and internal (turbulent channel) flows at high Reynolds number. 
    
    
    In contrast to previous works focusing on the Lagrangian viewpoint, in the present study we adopt an Eulerian, multiscale framework to investigate TI. In the Lagrangian framework, characteristic scales are identified through the average distance between particles, which is a function of time because particle separation increases (on average) with time~\cite{xu2016lagrangian}. It follows that particle motion senses (Eulerian) turbulent scales of varying sizes at different times due to particle relative dispersion in the flow, thus making difficult to perform a scale-dependent study of TI. The Lagrangian multiscale analysis is even more challenging in the case of wall turbulence, where the flow dynamics (hence, scales' features) is strongly dependent on the distance from the wall, and Lagrangian dispersion significantly depends on particle wall-normal position~\cite{polanco2018relative}. The Eulerian framework adopted here, instead, allows us to carry out a scale decomposition with ease by exploiting the Fourier transform, as it is typically done in the study of turbulent signals~\cite{smits2011high, jimenez2018coherent, pope2000turbulent}. Then, we take advantage from tools of nonlinear time-series analysis to quantify TI in the Eulerian viewpoint, by exploiting 1-dimensional velocity signals measured at various vertical (i.e., wall-normal) coordinates. Specifically, a parameter-free methodology based on horizontal visibility graphs (HVGs) is employed as the reference approach to capture temporal asymmetry in the signals~\cite{lacasa2012time}.
    
    
    \begin{figure*}
        \centering
        \includegraphics[width=0.99\textwidth]{fig2-01.eps}
        \caption{\textbf{Full-signal TI results for the turbulent boundary layer.} \textbf{a} Wall-normal behavior of $\widetilde{I}_k$. The maximum and minimum values of the irreversibility ratio $I_{k,r}=(I_k-\mu_{k,r})/\sigma_{k,r}$ are also highlighted, where $\mu_{k,r}$ and $\sigma_{k,r}$ are computed from an ensemble of 100 randomly-shuffled $u(t_i)$ signals. Conditional analyses of $\widetilde{I}_k$ with respect to $y^+$ and, \textbf{b} the local signal slope $|\mathrm{d}u^+/\mathrm{d}t^+|$ or, \textbf{c} $u^\prime/\sigma_u$. Intervals of $|\mathrm{d}u^+/\mathrm{d}t^+|$ and $u^\prime/\sigma_u$ are binned as $2:2:100$ percentiles.}
        \label{fig:fig2}
    \end{figure*}  
    
    
    Visibility networks represent a widely-used tool for nonlinear time-series analysis~\cite{zou2019complex}, including turbulent and vortical flows~\cite{iacobello2021review}, which have recently been exploited for TI analysis of both steady and unsteady phenomena~\cite{donges2013testing, suyal2014visibility, schleussner2015indications, lacasa2015time, xie2016time, flanagan2016irreversibility, gonzalez2020arrow}. The choice to adopt a HVG-based metric of TI is dictated by ease of implementation, and by results' robustness as arbitrarily-defined parameters are not required~\cite{lacasa2012time, lacasa2015time}. In fact, one of the main obstacles in quantifying TI is providing robust estimates, which is rarely achievable via traditional time-series symbolization~\cite{lacasa2012time, kennel2004testing}. However, like other measures of irreversibility, HVG-based metrics provide partial information on the system~\cite{lacasa2012time}. This is inevitable, since a single parameter can not fully represent a complex phenomenon. For this reason, aiming to strengthen our results' reliability, we extend the analysis to two alternative TI metrics: a higher-order (\textit{lag-reversibility}) correlation coefficient~\cite{Lawrance1991DirectionalityAR}, and a measure based on the fluctuation theorem~\cite{porporato2007irreversibility}. The results based on such alternative TI metrics (collected in the Supplementary Information) corroborate the HVG-based analysis thus advocating its applications.
    
    
      
        
\section*{Results}

    \noindent\textbf{Full-signal TI analysis of turbulent boundary layer.} A laboratory zero-pressure-gradient turbulent boundary layer at friction Reynolds number $R_\tau\approx 14{,}750$ is chosen as the representative (high Reynolds number) test case for inner flows (see a schematic in Fig.~\ref{fig:fig1}a)~\cite{baars2015wavelet}. Streamwise velocity time-series, recorded at varying wall-normal distances $y^+$ (the $+$ superscript indicates wall-units normalization by the friction velocity, $u_\tau$, and viscosity, $\nu$), are mapped into directed graphs following the horizontal visibility algorithm~\cite{luque2009horizontal}, as illustrated in Fig.~\ref{fig:fig1}b (see Methods). TI is then quantified as the Kullback-Leibler divergence, $I_k$, of the backward- ($k_b$) and forward-degree ($k_f$) probability distributions, $p(k_b)$ and $p(k_f)$, where the node degree $k$ quantifies the number of links adjacent to each node (e.g., see Fig.~\ref{fig:fig1}b)~\cite{lacasa2012time, lacasa2015time}. Time reversible signals imply $I_k=0$ (where zero is exactly reached for infinitely-long signals), while growing values $I_k>0$ indicate stronger levels of TI~\cite{lacasa2012time}.
        
    
        
    A full-signal (i.e., not scale-dependent) analysis is carried out first. Figure~\ref{fig:fig2}a shows the behavior of $\widetilde{I}_k=I_k/I_{k,max}$ (where $I_{k,max}$ is the maximum $I_k$ value along $y^+$), computed from full-length $u(t_i)$ signals. While a slightly-inclined plateau is observed in the log-layer with a drop in the external region, a $\widetilde{I}_k$ peak is observed in the near-wall region. The emergence of higher TI levels in the buffer layer (which is known to be a very active flow region) is a peculiar result, as this (near-wall) region of wall turbulence is characterized by the development of coherent motion and bursting events~\cite{jimenez2013near, jimenez2018coherent}, as well as modulation mechanisms~\cite{iacobello2021large}, whose complex dynamics can contribute to the generation of statistical TI. 
    
    In order to quantify the degree of reliability of $\widetilde{I}_k$ values, a reliability ratio, $I_{k,r}$, is also evaluated following Gonz{\'a}lez-Espinoza et al.~\cite{gonzalez2020arrow}: $I_{k,r}$ corresponds to the Z-score computed with respect to the mean and standard deviation of $I_{k}$ values from an ensemble of randomly-shuffled $u(t_i)$ signals (which are time-reversible). The maximum and minimum values of the reliability ratio (highlighted in Fig.~\ref{fig:fig2}a) are $I_{k,r}\gg1$, thus allowing one to ascertain with extreme confidence that streamwise velocity signals are irreversible. A parametric analysis on the impact of decreasing the time-series length on the wall-normal behavior of $I_k$ is discussed in Supplementary Section I, while the whole behavior of $I_{k,r}$ as a function of $y^+$ is reported in Supplementary Section II.
    

      
        
    
    \begin{figure*}
        \centering
        \includegraphics[width=0.99\linewidth]{fig3-01.eps}
        \caption{\textbf{Scale-dependent results for the turbulent boundary layer.} \textbf{a} Time series of $u^\prime$ (black) at $y^+=10.5$ and $y^+=464$, and their components $u^\prime_\uparrow$ (green) and $u^\prime_\downarrow$ (blue) obtained for $T_c^+\approx 170.3$. \textbf{b} Ratio $r_{\uparrow\downarrow}$ as a function of $y^+$ and $T_c^+$. \textbf{c} Energy spectra of $u^\prime$, $\phi_{uu}(T^+,y^+)$, pre-multiplied by frequency $f$ and normalized through $U_\tau^2$ as $\phi_{uu}^+=f\phi_{uu}/U_\tau^2$ (contour level range $0.2-2$, level step $0.2$). $T_{tr}^+(y^+)$ is shown as a red line, while the two vertical dotted lines refer to $y_{STS}^+\approx 13$ and $y_{LTS}^+=3.9\sqrt{R_\tau}$.}
        \label{fig:fig3}
    \end{figure*}
    
        
   
    To shed more light on the origin of TI in full-length $u(t)$ signals, a conditional analysis is performed with respect to: (i) the local slope in the time series, $|\mathrm{d}u^+/\mathrm{d}t^+|$, and (ii) the velocity fluctuations $u^{\prime}/\sigma_u$, where $u^{\prime}(t)=u(t)-U$ ($U$ is the mean velocity), while $\sigma_u$ is the root-mean-square velocity (see Fig.~\ref{fig:fig1}a). Values of $\widetilde{I}_k$ conditioned to the local slope $|\mathrm{d}u^+/\mathrm{d}t^+|$ for different $y^+$ coordinates are reported in Fig.~\ref{fig:fig2}b. The highest levels of TI are detected for intense temporal variations in the buffer layer, which are reminiscent of near-wall bursting events~\cite{jimenez2018coherent}. This link between high TI levels and bursting events is confirmed by the conditional analysis on turbulent fluctuations $u^{\prime}/\sigma_u$ (Fig.~\ref{fig:fig2}c), which points out larger $\widetilde{I}_k$ values residing in the range $-1 < u^{\prime}/\sigma_u < 0$ in the buffer layer. In fact, near-wall bursting events are typically detected as intervals of $u^\prime(t_i)$ starting at $u^{\prime}/\sigma_u=-1$ and ending at $u^{\prime}/\sigma_u<-0.25$~\cite{bogard1986burst, vinuesa2015documentation, tang2016bursting}, in very good agreement with intervals of large $\widetilde{I}_k$ in Fig.~\ref{fig:fig2}c. 
    \\Results from Fig.~\ref{fig:fig2} are also in qualitative accordance with previous studies where irreversibility was related to (i) the velocity-gradient tensor perceived by dispersing Lagrangian particles~\cite{jucha2014time}, and (ii) an asymmetry between growth and decay of signal fluctuations (i.e., where local signal slopes are large) in non-equilibrium systems~\cite{luchinsky1997irreversibility, bertini2001fluctuations, giberti2007temporal}. 
   
   
    
        

   
    \vspace{10pt}
    \noindent\textbf{Multiscale analysis of TI.} A scale-dependent investigation is carried out here to highlight how different turbulent scales contribute to TI of the streamwise velocity. Accordingly, $u(t_i)$ signals are Fourier-filtered at various cut-off periods $T_c^+=1/f_c^+$ such that $u(t_i)= u_\uparrow(t_i)+ u_\downarrow(t_i)$, where $u_\uparrow$ and $u_\downarrow$ indicate the low-pass and high-pass filtered velocity components of $u$, respectively (see Fig.~\ref{fig:fig3}a). The HVG-algorithm is then applied to the $u_\uparrow$ and $u_\downarrow$ velocity components, thus obtaining the corresponding TI measures $I_{k,\uparrow}$ and $I_{k,\downarrow}$, respectively. Figure~\ref{fig:fig3}b shows the ratio 
        
    \begin{equation}
        r_{\downarrow\uparrow} = \log_{10}{\frac{I_{k,\downarrow}}{I_{k,\uparrow}}},    
    \end{equation}
        
    \noindent as a function of $y^+$ and $T_c^+$: blue-shaded areas ($r_{\downarrow\uparrow}<0, I_{k,\downarrow}<I_{k,\uparrow}$) indicate a greater contribution to TI coming from temporal scales of $u(t)$ with larger periods than $T_c^+$; green-shaded areas ($r_{\downarrow\uparrow}>0, I_{k,\downarrow}>I_{k,\uparrow}$) instead indicate a greater contribution to TI coming from temporal scales of $u(t)$ with smaller periods than $T_c^+$. It should also be pointed out that, in the limiting cases of $T_c^+\gg 0$ (bounded by Nyquist frequency) and $T_c^+\rightarrow 0$, all temporal scales will contribute either to $I_{k,\downarrow}$ (hence $r_{\downarrow\uparrow}\gg 0$) or $I_{k,\uparrow}$ (hence $r_{\downarrow\uparrow}\ll 0$).
        
    
    The transitional line $T^+_{tr}(y^+)$ for which $r_{\downarrow\uparrow}\approx 0$ (i.e., $I_{k,\downarrow}\approx I_{k,\uparrow}$, see red contour in Fig.~\ref{fig:fig3}b) discriminates between the region where the cut-off filter is $T_c^+<T^+_{tr}$ (such that $r_{\downarrow\uparrow}<0$), and the region where the cut-off filter is $T_c^+>T^+_{tr}$ (such that $r_{\downarrow\uparrow}>0$). The $T^+_{tr}(y^+)$ boundary displays a peculiar behavior, made up of two bumps emerging far from the wall. The first bump arises at $y^+\approx 450$ in the log-layer, while an additional bump is present at $y^+\approx 14{,}500$ in the intermittency region (see sketch in Fig.~\ref{fig:fig1}a), which is the flow layer where the entrainment process occurs owing to the proximity with the external (non-turbulent) flow~\cite{pope2000turbulent}. The occurrence of bumps in the $T^+_{tr}(y^+)$ behavior suggests the presence of relevant phenomena in the flow, leading to higher TI levels at larger temporal scales (since higher cut-off periods $T_c^+=T^+_{tr}$ are needed to have the balance $I_{k,\uparrow}= I_{k,\downarrow}$, i.e., $r_{\downarrow\uparrow}=0$).
        
        
   
    \vspace{10pt}
    \noindent\textbf{Relating TI levels and turbulent flow scales.} With the aim to connect the $r_{\downarrow\uparrow}$ behavior (specifically, the presence of bumps in $T^+_{tr}(y^+)$) with the underlying structure of the turbulent boundary layer, we show in Fig.~\ref{fig:fig3}c the pre-multiplied energy spectrum of $u$, $\phi_{uu}^+$. Two peaks are distinguishable in Fig.~\ref{fig:fig3}c around $y_{STS}^+\approx 13$ (buffer layer) and $y_{LTS}^+\approx 470$ (log-layer): they are associated with the development of organized coherent flow structures commonly referred to as \textit{small turbulent scales} (STS) and \textit{large turbulent scales} (LTS), respectively, the latter emerging at high Reynolds numbers~\cite{smits2011high, baars2015wavelet, iacobello2021large}.
        
    \begin{figure*}
        \centering
        \includegraphics[width=0.99\linewidth]{fig4-01.eps}
        \caption{\textbf{Full-signal TI results for the turbulent channel flow.} \textbf{a} Wall-normal behavior of $\widetilde{I}_k$. The maximum and minimum values of the irreversibility ratio $I_{k,r}$ are also highlighted. Conditional analysis of $\widetilde{I}_k$ with respect to $y^+$ and \textbf{b} the local signal slope $|\mathrm{d}u^+/\mathrm{d}t^+|= U^+|\mathrm{d}u^+/\mathrm{d}x^+|$, and \textbf{c} $u^\prime/\sigma_u$. Intervals of $|\mathrm{d}u^+/\mathrm{d}t^+|$ and $u^\prime/\sigma_u$ are binned as $2:2:100$ percentiles.}
        \label{fig:fig4}
    \end{figure*} 
    
    
    The $T^+_{tr}(y^+)$ line (i.e., where $I_{k,\uparrow}=I_{k,\downarrow}$) is also displayed (in red) in Fig.~\ref{fig:fig3}c, showing that both STS and LTS spectral peaks stand above $T_{tr}^+(y^+)$. Hence, for a cut-off filter $T_c^+$ such that $T_c^+\equiv T_{tr}^+$, the contribution to TI coming from STS and LTS is -- in both cases -- enclosed in $I_{k,\uparrow}$ since, by definition, $I_{k,\uparrow}$ accounts for the contribution to TI coming from flow scales larger than $T_c^+= T_{tr}^+$. By increasing the value of the cut-off period, i.e. $T_c^+>T_{tr}^+$, a growing contribution to TI coming from flow scales associated with STS and LTS is also enclosed in $I_{k,\downarrow}$, leading to $I_{k,\downarrow}>I_{k,\uparrow}$ (i.e., $r_{\downarrow\uparrow}>0$). Therefore, the presence of energetic coherent structures in the flow (STS and LTS) significantly affect the levels of TI of $u(t)$ at varying wall-normal locations, as captured by the relative intensity of $I_{k,\downarrow}$ and $I_{k,\uparrow}$.
    \\Furthermore, we note that a larger extent of dark-green areas (i.e., $r_{\downarrow\uparrow}\geq 1$ or $I_{k,\downarrow}/I_{k,\uparrow}\geq 10$) emerges in Fig.~\ref{fig:fig3}b closer to the $T^+_{tr}(y^+)$ boundary at $y_{STS}^+$ than at $y_{LTS}^+$. The different contribution to TI at various $y^+$ can hence be related to the different intensities of energetic peaks, as the inner peak at $y_{STS}^+$ is stronger than the outer peak at $y_{LTS}^+$. This resonates with the previous results of Fig.~\ref{fig:fig2} where larger levels of TI where observed in the buffer layer from a full-signal perspective. This outcome corroborates the key role of the buffer layer as region where complex dynamical processes are at play which lead to greater TI.

    
    
    Concerning the intermittency region (see Fig.~\ref{fig:fig1}a), previous works~\cite{chauhan2014turbulent} found that a characteristic (large) length scale of entrainment is $\Lambda_x^+\approx 1.7\,\delta^+$. The $T^+_{tr}$ bump in this region occurs at $T^+\approx 600$ and corresponds to a length scale (using Taylor hypothesis) of $600\,U_\infty^+\approx 1.3\,\delta^+$, which is -- similarly to STS and LTS -- close to the characteristic flow scale $1.7\,\delta^+$ although slightly smaller. The present result reveals that, despite very low TI levels are detected from a full-signal analysis in the intermittency region (see the right tail of Fig.~\ref{fig:fig2}a), a significant link between TI and the entrainment process can be established.
    
    
   
    \vspace{10pt}
    \noindent\textbf{TI analysis of turbulent channel flow.} The close relation between TI and the arrangement of coherent flow motion in the turbulent boundary layer is further validated here by using a second test case, consisting of a turbulent channel flow at $R_\tau=2{,}003$~\cite{hoyas2006scaling}. Streamwise velocity signals $u(x)$ are extracted from a direct numerical simulation extensively validated by previous studies~\cite{hoyas2006scaling, jimenez2008turbulent}. The common assumption of Taylor's hypothesis is adopted to investigate $u(x)$ signals as equivalent time series $u(t)$ via HVG~\cite{smits2011high, iacobello2021large}.
    
    Figure~\ref{fig:fig4} reports the same quantities as shown in Fig.~\ref{fig:fig2} for the turbulent boundary layer, highlighting that results for the turbulent channel flow are in agreement with those for the turbulent boundary layer. In particular, in Fig.~\ref{fig:fig4}a, we observe a peak of $\widetilde{I}_k$ in the buffer layer (as for the turbulent boundary layer) with values of the confidence ratio $I_{k,r}>1$ (highlighted by red boxes). However, due to the lack of an intermittency region for the channel flow, we do not observe any drop of $\widetilde{I}_k$ for $y^+\rightarrow y^+_{max}=2{,}003$ in Fig.~\ref{fig:fig4}a. Moreover, results from the conditional analysis in the channel flow (Fig.~\ref{fig:fig4}b-c) match those for the boundary layer flow (Fig.~\ref{fig:fig2}b-c), with a peak of $\widetilde{I}_k$ for strong local slopes $|\mathrm{d}u^+/\mathrm{d}t^+|$ (\figurename~\ref{fig:fig4}b) and for $u^\prime/\sigma_u\approx -0.5$ (\figurename~\ref{fig:fig4}c) in the buffer layer.
    
    
    \begin{figure}
        \centering
        \includegraphics[width=0.5\linewidth]{fig5-01.eps}
        \caption{\textbf{Scale-dependent results for the turbulent channel flow.} Colored areas corresponds to the ratio $r_{\uparrow\downarrow}$, while the black contours refers to the pre-multiplied energy spectra $\phi_{uu}^+$ (contour level range $0.25-2$, level step $0.25$). The two yellow dashed lines refer to the $\lambda_x^+\sim{y^+}^{3/7}$ scaling~\cite{monty2009comparison}, in terms of time scales $T^+=\lambda_x^+/U^+$.}
        \label{fig:fig5}
    \end{figure}
    
    
    Figure~\ref{fig:fig5} shows the energy spectra (black contours) together with the ratio $r_{\downarrow\uparrow}$ (colored area), where $I_{k,\downarrow}$ and $I_{k,\uparrow}$ are evaluated for different cut-off periods $T_c^+=\lambda_x^+/U^+$ (where $\lambda_x$ is the wavelength along $x$). Figure~\ref{fig:fig5} is in excellent agreement with outcomes from Fig.~\ref{fig:fig3}b-c, highlighting larger green-shaded areas in correspondence of the spectral peak in the buffer layer at $y^+_{STS}$. The notable difference is the absence of a bumps at larger $y^+$, due to two main reasons. First, the intermittency region is absent in a channel flow, hence the rightmost bump of Fig.~\ref{fig:fig3}b does not appear in Fig.~\ref{fig:fig5}. Second, the bump in correspondence of $y^+_{LTS}$ does not emerge because very large-scale flow motions extend much further towards higher $y^+$ in channel flows than in boundary layers, although their spectral energy similarly decreases for $y^+\rightarrow y^+_{max}$~\cite{monty2009comparison}. This statement is supported by scaling arguments of $r_{\downarrow\uparrow}$ in the log-layer. In fact, the transitional (red) line $T_{tr}^+$ in Fig.~\ref{fig:fig5} follows the same ${y^+}^{3/7}$ scaling law (see dashed yellow lines in Fig.~\ref{fig:fig5}) found for the size of very large-scale motion in internal wall-bounded turbulence~\cite{monty2009comparison}. This outcome further corroborates the strong connection between time irreversibility and the spatio-temporal development of coherent structures in the flow.
        
        

   
\section*{Discussion}

    The present study provides a novel perspective into the statistical irreversibility of high-Reynolds-number turbulent flows, that complements previous works looking at broken temporal symmetry in the statistics of Lagrangian tracers. In particular, this work represents a first effort to characterize TI from an Eulerian viewpoint in wall-bounded turbulent flows, which have been much less investigated than other flow configurations from the TI point of view. Exploiting tools of nonlinear time-series analysis, we are able to quantify TI from Eulerian data by explicitly highlighting the contributions to TI coming from various flow scales, as well as the effect of the wall-normal (spatial) coordinate, revealing non-trivial TI patterns. 
    
    
    Our findings -- relying on the visibility graph-based approach and corroborated through additional methodologies, as reported in Supplementary Section III -- point out that broken temporal symmetry in the streamwise velocity is significantly linked to the underlying (space-dependent) organized flow structure of wall turbulence. In fact, we show that the scale-dependent TI levels follow non-monotonic trends at all wall-normal coordinates, highlighting a non-trivial contribution to TI from smaller and larger scales. We find that TI in the proximity of the wall is dominated by small turbulent scales associated with the inner spectral-energy peak, thus suggesting a connection between TI and the dynamical process related to the near-wall (regeneration) cycle~\cite{jimenez2018coherent}. This claim is partially supported by the large TI levels detected in this work in correspondence of burst events in the buffer layer (see Fig.s~\ref{fig:fig2}b-c and \ref{fig:fig4}b-c). Moving away from the wall, the largest contribution to TI shifts from small to large turbulent scales, that appear at high Reynolds number and make a significant contribution to the turbulent kinetic energy and Reynolds stress production~\cite{smits2011high, jimenez2018coherent}. Moreover, significant levels of time irreversibility in the intermittency region of turbulent boundary layers are associated with the characteristic scale of entrainment process. 
    \\In this picture, TI patterns in the near-wall region are in agreement between turbulent boundary layer (external flow) and turbulent channel (internal flow), while dissimilarities emerge far from the wall due to intrinsic differences in the flow features of external and internal wall turbulence (e.g., the absence of an intermittency region in internal flows). Overall, the outcomes for both turbulent channel and turbulent boundary layer flows indicate that TI patterns are not limited to a particular flow case, but are distinctive feature of wall-bounded turbulent flows at high Reynolds number.
    
    It is worth noting here that the full-signal time irreversibility behavior shown in Fig.s~\ref{fig:fig2}a and \ref{fig:fig4}a displays a wall-normal trend similar to the skewness of the spatial derivative $\partial u/\partial x$ reported in previous works~\cite{vreman2014statistics, djenidi2017skewness}, with important implications. Such similarity suggests that $\partial u/\partial x$ can be reinterpreted as a surrogate measure of TI, although $\partial u/\partial x$ is computed at fixed time. However, the TI reinterpretation of $\partial u/\partial x$ is conditional to the applicability of Taylor's hypothesis, which allows the conversion of spatial signals $u(x)$ into equivalent time series $u(t)$, and \textit{vice versa}. This point is crucial because an unconsidered use of spatial asymmetry metrics to quantify TI may lead to inconsistent outcomes, especially for non-canonical flow setups~\cite{zorzetto2018extremes}. For the flow cases considered in this study, the application of classical Taylor's hypothesis -- namely, when the convection velocity equals the mean velocity -- led to consistent outcomes from time- and spatial-series, as discussed for the turbulent channel flow case, but caution is still needed. In this sense, refined formulations of the Taylor's hypothesis -- e.g., accounting for modulation mechanisms in wall turbulence~\cite{iacobello2021large} -- could help to shed more light on the duality between temporal and spatial asymmetries, and will be explored in future works. 
    
    Furthermore, it should be pointed out that different metrics need to be examined and compared in order to provide reliable TI outcomes (as done in this study; see Supplementary Section III). Since TI can be quantified using different tools~\cite{Lawrance1991DirectionalityAR, lacasa2012time}, the choice of the methodologies can be problem dependent but should rest upon the method interpretability and robustness (e.g., avoiding subjective binning or symbolization procedures whenever possible). In light of this, visibility graphs represent a versatile and robust approach, which deserves further consideration for the investigation of turbulence signals~\cite{iacobello2021review}. 
    
    
    In conclusion, the present Eulerian analysis identified scale-dependent TI patterns in wall turbulence originating in correspondence of the development of small and large energetic structures at various, characteristic, wall-normal distances. Present outcomes can hence open up new research avenues aiming to gain further insights into the interplay between spatial and temporal asymmetries, stimulate new developments in Eulerian and Lagrangian wall-turbulence modeling~\cite{fang2012time, xia2013lagrangian, bae2022scientific}, as well as foster flow control strategies, owing to the key role played by coherent structures in drag-reduction actuation~\cite{marusic2021energy}.
    
   


    
    

\section*{Methods}    
    
    \noindent\textbf{Horizontal visibility graphs and TI quantification.} Given a time series $u(t_i)$ evaluated at discrete times $t_i$, a horizontal visibility graph representation of $u(t_i)$ is obtained, (i) assigning each time $t_i$ to a network node, and (ii) linking two nodes $t_i$ and $t_j$ when the conditions
        \begin{equation}
            u(t_i)>u(t_l) \wedge u(t_j)>u(t_l),\label{eq:hvg_def}
        \end{equation}
        
        \noindent are satisfied for all $t_l$ such that $t_i<t_l<t_j$~\cite{luque2009horizontal}. Consecutive data points, $u(t_i)$ and $u(t_{i+1})$, are always linked by construction (i.e., the index $l$ can be null). The resulting graph is stored as a binary adjacency matrix $A_{ij}$, whose entries are $A_{ij}=1$ if and only if nodes $i$ (corresponding to datum $u(t_i)$) and $j$ (corresponding to datum $u(t_j)$) are linked with each other as per Eq.~(\ref{eq:hvg_def}). Horizontal visibility networks are typically constructed as \textit{undirected} graphs, namely the link direction is not taken into account and $A_{ij}=A_{ji}$ for any node pair $(i,j)$. In this work, however, we consider \textit{directed} graphs, namely explicitly differentiating between forward-in-time links $A_{ij}=1$ (i.e., when $t_j>t_i$) and backward-in-time links $A_{ji}=1$ (i.e., when $t_j<t_i$). Figure~\ref{fig:fig1}b shows an exemplifying discrete signal, where HVG links for node $t_j$ are represented by colored arrows. 
        
        The amount of links of each node is referred to as degree centrality $k(t_i)$. Here, we distinguish between forward-degree $k_f$ and backward-degree $k_b$ (commonly referred to as out- and in-degree~\cite{zou2019complex}) as the number of links pointing towards increasing and decreasing time, respectively (Fig.~\ref{fig:fig1}b). \citet{lacasa2012time} proposed to quantify TI in a time series as the Kullback-Leibler divergence of the backward- and forward-degree distributions:
        
        \begin{equation}
            I_k=\sum_{k_b,k_f}{p(k_b)\log{\frac{p(k_b)}{p(k_f)}}}, \label{eq:Ik_def}
        \end{equation}
        
        \noindent where $p(k_b)$ and $p(k_f)$ are the marginal probability distributions of $k_b$ and $k_f$, respectively. The degree of reliability of $I_k$ can be assessed evaluating the \textit{irreversibility ratio} 
        \begin{equation}
            I_{k,r}=\frac{I_k-\mu_{k,r}}{\sigma_{k,r}},    
        \end{equation}
        \noindent where $\mu_{k,r}$ and $\sigma_{k,r}$ are the mean and standard deviation, respectively, of $I_k$ values calculated from an ensemble of signals obtained through a random (null) model~\cite{gonzalez2020arrow}. As random signals by definition privilege no direction, $I_{k,r}\gg 1$ signifies that a time series is TI with extreme confidence~\cite{gonzalez2020arrow}.
        
        
        It should be noted that, in general, the Kullback-Leibler divergence is not a symmetric measure, namely $\sum{ p(k_b)\log\left[p(k_b)/p(k_f)\right]}\neq \sum{p(k_f)\log\left[p(k_f)/p(k_b)\right]}$. However, switching the position of $p(k_b)$ and $p(k_f)$ in Eq.~(\ref{eq:Ik_def}) does not significantly change the outcomes and hence our conclusions (as also reported in previous studies~\cite{lacasa2012time,lacasa2015time}). Moreover, $I_k$ is not an additive measure due to the non-linearity of the Kullback-Leibler divergence, thereby the sum $I_{k,\uparrow}+I_{k,\downarrow}$ (obtained through high-pass and low-pass filtering of the signal) is in general not equal to $I_k$ (obtained from the full signal).
        
        \vspace{10pt}
        \noindent\textbf{Turbulent boundary layer experiments.} The boundary layer was experimentally obtained in the wind-tunnel facility of the University of Melbourne~\cite{MARUSIC2020dataset}. The friction Reynolds number is $R_\tau= \delta U_\tau/\nu \approx 14{,}750$, where $\delta= 0.361$\,m and $U_\tau = 0.626\,\mathrm{m\,s^{-1}}$ are the boundary layer thickness and the friction velocity, respectively, while $\nu= 1.532\times 10^{-5}\, \mathrm{m^2s^{-1}}$ is the kinematic viscosity of air. The value of the Reynolds number of the experiment is large enough to ensure a wide range of temporal scales, thereby allowing very-large scale motions to develop~\cite{baars2015wavelet}. Further details of the experiments can be found in Baars et al.~\cite{baars2015wavelet}.
        
        The time series of $u(t_i)=u^\prime(t_i)+U$ -- where $U(y)$ is the local (time-averaged) mean velocity, while $u^\prime$ are turbulent fluctuations (Fig.~\ref{fig:fig1}a) -- were recorded at $41$ wall-normal coordinates $y$, while fixing the streamwise and spanwise coordinates, $x$ and $z$, respectively. Wall-units normalization is indicated through the $+$ superscript such that $u^+=u/U_\tau$ and $y^+=y U_\tau/\nu$. Time series were collected for three acquisition cycles of $120$\,s at sampling frequency $f_s=20$\,kHz, thereby results are intended to be averaged over the three acquisition cycles.
        
        
        \vspace{10pt}
        \noindent\textbf{Turbulent channel numerical simulation.} Data from a direct numerical simulation of a turbulent channel flow at $R_\tau= 2{,}003$ are used as a representative case of internal flow. The resulting dataset is available online on the webpage of the Fluid Dynamics Group of Universidad Politecnica de Madrid (\url{https://torroja.dmt.upm.es/turbdata})~\cite{hoyas2006scaling}. The simulation is run on a smooth-wall channel setup with periodic boundary conditions in the streamwise ($x$) and spanwise ($z$) directions. The domain size is $8\pi \delta \times 2\delta \times 3\pi \delta$ in the streamwise, vertical and spanwise directions respectively, where $\delta$ is the half-channel height. The numerical grid consists of $6{,}144$ and $4{,}608$ uniformly-spaced grid points in the streamwise and spanwise direction, respectively, while a non-uniform grid with $633$ points is used in the wall-normal direction. Further details can be found in Hoyas and Jim{\'e}nez~\cite{hoyas2006scaling}.
       
       
       Spatial series $u(x)$ are used in this work to quantify TI, where the Taylor's hypothesis $\Delta x=-U(y)\Delta t$ is used to transform the spatial signal into time series~\cite{pope2000turbulent, jimenez2018coherent}. In this regard, the spatial resolution $\Delta x^+=8.2$ is an equivalent to the sampling time step of the turbulent boundary layer. It should be noted here that, the application of Taylor's hypothesis does not affect the computation of $I_k$ from spatial signals $U(x)$ in the channel flow, because the horizontal visibility algorithm is insensitive to constant re-scaling of the horizontal axis (more generally, it is insensitive to affine transformations)~\cite{luque2009horizontal}. One time snapshot from the database is used, thereby ensemble-averages of the results are only performed along the spanwise (homogeneous) direction $z$, over $288$ uniformly-spaced locations. Further increases of the averaging sample size do not substantially change the discussed outcomes.    


\section*{Data availability}
All the data analyzed in this paper are openly available at locations referenced herein.

     


\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%

\section*{References}}
\begin{document}


\title{Supplementary Information for:\\``Coherent structures at the origin of time irreversibility in wall turbulence''

\author{G. Iacobello}
    \email{g.iacobello@surrey.ac.uk}
    \affiliation{School of Mechanical Engineering Sciences, University of Surrey, Guildford, GU2 7XH, UK}

\author{S. Chowdhuri}
    \affiliation{Department of Civil and Environmental Engineering, University of California, Irvine, CA 92697, USA}

\author{L. Ridolfi}
    \affiliation{Department of Environmental, Land and Infrastructure Engineering, Politecnico di Torino, Turin, 10129, Italy}

\author{L. Rondoni}
    \affiliation{Department of Mathematical Sciences, Politecnico di Torino, Turin, 10129, Italy}
    \affiliation{Istituto Nazionale di Fisica Nucleare, INFN, Sezione di Torino, Turin, 10125, Italy}
   
\author{S. Scarsoglio}
    \affiliation{Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Turin, 10129, Italy}




            




\begin{abstract}
    In this Supplementary Information we provide evidence of the robustness of the visibility-based measure of time irreversibility in terms of time-series length (Section~\ref{SM_sec:param_leng}) and deviations from random signals (Section~\ref{SM_sec:irr_ratio}). To compare and corroborate visibility-based results discussed in the main text, additional results are provided using alternative methodologies (Section~\ref{SM_sec:res_alt_meth}).
\end{abstract}


\maketitle


\section{Impact of varying signal length on the $I_k(y^+)$ behavior}\label{SM_sec:param_leng}

In this section, we assess the impact on the visibility-based TI measure $I_k$ (Eq.~(1) in the main text) of smaller temporal lengths $T$ of the streamwise velocity signals $u(t_i)$, with respect to the full temporal window $T_f=120\,$s (as used in the main analysis). With this aim, $N$ sub-samples of length $T$ with a $50\%$ overlap between consecutive samples are extracted from the $u$ signals. $I_k$ values are then computed and average over all $N$ sub-samples to obtain $\overline{I_k}(T)$. To quantify the effect of varying $T/T_f$, we compute the relative ratio

\begin{equation}
    \epsilon_{I_k}(T)=\frac{\overline{I_k}(T)-I_k(T_f)}{I_k(T_f)},
\end{equation}

\noindent which is shown in Fig.~\ref{SM_fig1}, where different colored lines correspond to decreasing temporal lengths $T/T_f$. Relative deviations from the $I_k(T_f)$ values are observed to be smaller than $10\%$ if the velocity time series are shortened up to $20\%$ of $T_f$. For the vast majority of the flow domain, relative variations below $5\%$ are observed, even considering $1\%$ of the time-series length (black line in Fig.~\ref{SM_fig1}). Larger deviations are only obtained if signals are significantly shortened (cyan and black lines in Fig.~\ref{SM_fig1}), although such variations are limited to $25\%$ of $I_k(T_f)$ and are confined to the intermittency region ($y^+>6000$~\cite{baars2015wavelet}).  We can conclude that $I_k$ is a robust measure of time irreversibility with respect to the signal length, and hence the values $I_k(T_f)$ used in the main analysis are not significantly affected by the signal length.


    \begin{figure}
        \centering
        \includegraphics[width=0.9\linewidth]{SM_fig1-01.eps}
        \caption{\textbf{Effect of time-series length on $I_k$}. Wall-normal behavior of the $\epsilon_{I_k}(T)$ ratio (in percentage) for varying time-series lengths $T/T_f$ (such that $N=2(T_f/T)-1$) of the streamwise velocity signals. Different colors refer to different $T/T_f$ values.}
        \label{SM_fig1}
    \end{figure}


\section{Wall-normal behavior of the irreversibility ratio $I_{k,r}$}\label{SM_sec:irr_ratio}

In order to assess the degree of confidence that non-zero $I_k$ values represent significant time irreversibility levels, we computed the \textit{irreversibility ratio} defined as

\begin{equation}
    I_{k,r}=\frac{I_k-\mu_{k,r}}{\sigma_{k,r}}, \label{SM_eq:irrev_ratio}
\end{equation}

\noindent where $\mu_{k,r}$ and $\sigma_{k,r}$ are the mean and standard deviation, respectively, of $I_k$ values calculated from an ensemble of signals from a random (null) model~\cite{gonzalez2020arrow}. In fact, random signals are expected to be time reversible, such that $I_k\rightarrow 0$ for $T\rightarrow \infty$. In Fig.~2a of the main text, we have highlighted the maximum and minimum values of $I_{k,r}$ computed with respect to an ensemble of $100$ randomly-shuffled signals. Here, the whole wall-normal behavior of $I_{k,r}$ is shown in Fig.~\ref{SM_fig2}a for the sake of completeness, highlighting that $I_{k,r}\gg 1$ at any $y^+$, thereby concluding that $I_k$ is a reliable measure of time irreversibility.



    \begin{figure}[t]
        \centering
        \includegraphics[width=0.99\linewidth]{SM_fig2-01.eps}
        \caption{\textbf{Confidence analysis on $I_k$.} Irreversibility ratio $I_{k,r}$ as a function of $y^+$ computed from: \textbf{a} the full velocity signals $u$; \textbf{b} filtered velocity signals, $u_\downarrow$ and $u_\uparrow$. The two red $\blacktriangle$ markers in \textbf{a} highlight the maximum and minimum $I_{k,r}$ values, as reported in Fig.~2a of the main text.}
        \label{SM_fig2}
    \end{figure}

To assess the degree of confidence of time irreversibility for scale-filtered signals, we also show in Fig.~\ref{SM_fig2}b the irreversibility ratio computed by replacing $I_k$ with $I_{k,\downarrow}$ or $I_{k,\uparrow}$ in Eq.~(\ref{SM_eq:irrev_ratio}) for three representative cut-off filters $T^+_c=\lbrace 10, 100, 1000 \rbrace$. In all cases, the irreversibility ratio is much larger than unity, thus allowing us to ascertain that full streamwise velocity signals (see Fig.~\ref{SM_fig2}a) and their filtered components (in the range of filters considered, see Fig.~\ref{SM_fig2}b) are irreversible with extreme confidence. 

\newpage
\section{Time irreversibility analysis via alternative methodologies}\label{SM_sec:res_alt_meth}

The aim of this section is to provide evidence that alternative methodologies used to quantify time irreversibility (TI) lead to results qualitatively similar to those obtained through the visibility-based approach, thus corroborating the analysis presented in the main text. Two alternative approaches are here adopted: a measure based on the fluctuation theorem~\cite{porporato2007irreversibility}, and a higher-order (\textit{lag-reversibility}) correlation coefficient~\cite{Lawrance1991DirectionalityAR}. The different methodologies are compared with each other, highlighting their advantages and disadvantages in providing reliable and robust results of TI.


\subsection{Fluctuation theorem-based measure of TI}\label{SM_subsec:FT}

    The first alternative approach rests upon insights deriving from the fluctuation theorem~\cite{evans1993probability, porporato2007irreversibility} to obtain a measure of TI from discrete time series. In particular, Porporato et al.~\cite{porporato2007irreversibility} proposed a measure of asymmetry involving conditional probabilities of a variable and its first time derivative, which is written here for the streamwise velocity $u(t_i)$ as:

    \begin{equation}
        \langle Z_2^*\rangle=\sum_{u}{\sum_{\Delta u}{p(u,\Delta u)\log{\frac{p(u,\Delta u)}{p(u,-\Delta u)}}}},\label{SM_eq:Z2star}
    \end{equation}

    \noindent where $\Delta u(\theta)=\left[u(t_i+\theta)-u(t_i)\right]/\theta$ is the numerical (temporal) derivative of $u$, $\theta=n\Delta t_s$ is a time lag proportional to the sampling time step $\Delta t_s$, and $p(\bullet,\bullet)$ indicates the joint probability distribution. Figure~\ref{SM_fig3} shows the wall-normal behavior of $\langle Z_2^*\rangle$ for $\theta=\Delta t_s$ (i.e., the minimum available time step) and for different number of bins $N_b$ (used to sample the $u$ and $\Delta u$ signals). The wall-normal behavior of $\langle Z_2^*\rangle$ qualitatively resembles the result obtained via the visibility-based approach, i.e. $I_k(y^+)$, illustrated in Fig.~2a of the main text. In fact, a peak of $\langle Z_2^*\rangle$ is detected in the buffer layer with lower values elsewhere, thus corroborating the $I_k$ results. However, significant variations of $\langle Z_2^*\rangle$ values are observed for different number of bins, which is an arbitrary parameter to be set \textit{a priori}. In this regard, the visibility based approach displays a striking advantage as it does not require any \textit{a priori} parameter. 
    

    \begin{figure}
        \centering
        \includegraphics[width=0.99\linewidth]{SM_fig3-01.eps}
        \caption{\textbf{Binning effect on fluctuation theorem-based results at minimum lag.} Wall-normal behavior of $\langle Z_2^*\rangle$ for $\theta=\Delta t_s$. The legend reports the number of bins, $N_b\times N_b$, used to sample the $u$ and $\Delta u$ signals.}
        \label{SM_fig3}
    \end{figure}
    
    
    
    \begin{figure}[t]
        \centering
        \includegraphics[width=0.9\linewidth]{SM_fig4-01.eps}
        \caption{\textbf{Comparison between visibility-based and fluctuation theorem-based results at different scales.} \textbf{a} Two-dimensional plot of $\widetilde{I}_{k,\uparrow}=I_{k,\uparrow}/I_{k,\uparrow,max}$ as a function of the wall-normal coordinate $y^+$ and the cut-off period $T^+_c$. The white, dotted, horizontal line is $T_{c,min}^+=1/f_{Nyquist}^+=2/f_s^+$, being $f_s$ the sampling frequency. The $\widetilde{I}_{k,\uparrow}$ values for $T^+_c<T_{c,min}^+$ have then been inserted as $\widetilde{I}_{k,\uparrow}=\widetilde{I}_{k}$ owing to the fact that for very small cut-off periods $u_\uparrow(t_i)\approx u(t_i)$. (b-d) $\widetilde{\langle Z_2^*\rangle} = \langle Z_2^*\rangle/\langle Z_2^*\rangle_{max}$ as a function of $y^+$ and the time lag $\theta^+$, for different number of bins: \textbf{b} $N_b=100$; \textbf{c} $N_b=500$; \textbf{d} $N_b=1000$. In all plots, the gray line refers to $T_{tr}^+(y^+)$ corresponding to the red line in Fig.~3b-c of the main text.}
        \label{SM_fig4}
    \end{figure}
    
    
    A scale-dependent analysis is also carried out by increasing the time lag $\theta$ in Eq.~(\ref{SM_eq:Z2star})~\cite{porporato2007irreversibility}. It should be noted that temporal scales are here associated with time lags $\theta$, while in the main text temporal scales are associated with Fourier-based cut-off periods $T_c$. For a given $\theta$ value, therefore, temporal scales greater than or equal to $\theta$ are retained in the computation of $\Delta u$ (i.e., the signals are implicitly under-sampled), thus miming the $I_{k,\uparrow}$ behavior. A noteworthy methodological difference here is that $I_{k,\downarrow}$ is also computed using visibility networks (which can be compared with $I_{k,\uparrow}$, hence providing a richer picture of scale-dependent TI), while the $\langle Z_2^*\rangle$-based approach has not such additional information. Figure~\ref{SM_fig4}a shows $\widetilde{I}_{k,\uparrow}={I}_{k,\uparrow}/{I}_{k,\uparrow,max}$ as a function of $y^+$ and $T_c^+$, and, for comparison, $\widetilde{\langle Z_2^*\rangle} = \langle Z_2^*\rangle/\langle Z_2^*\rangle_{max}$ is shown in Fig.~\ref{SM_fig4}b-d as a function of $y^+$ and $\theta^+$ for different number of bins. The transitional line $T_{tr}^+(y^+)$ (red line in Fig.~3b-c of the main text) is also overlapped in Fig.~\ref{SM_fig4} as a gray line. As expected, $\widetilde{\langle Z_2^*\rangle}$ mimics $\widetilde{I}_{k,\uparrow}$; in particular: 
    
    
    \begin{enumerate}[i)]
        \item high-values of $\widetilde{\langle Z_2^*\rangle}$ are detected for small time lags $\theta^+$ in the buffer layer;
        
        \item mild values of $\widetilde{\langle Z_2^*\rangle}$ are also detected for larger $\theta^+$ in the intermittency region and log-layer, leading to bulge-like pattern appearing for $N_b=100$ (Fig.~\ref{SM_fig4}b) and $N_b=500$ (Fig.~\ref{SM_fig4}c). This pattern is not evident for $N_b=1000$ (Fig.~\ref{SM_fig4}d), so we can argue that $N_b=1000$ is supposedly too large to properly sample the $u$ and $\Delta u$ signals. In this regard, the lack of \textit{a priori} parameters (such as $N_b$) is a significant advantage of visibility networks.
    \end{enumerate}
    
    
    In summary, results from Fig.~\ref{SM_fig4} are in qualitative accordance with the scale-dependent analysis shown in Fig.~3b-c of the main text for the visibility-based approach, in spite of the different nature of the two methodologies employed.
    
    

\subsection{Lag-reversibility-based measure of TI}\label{SM_subsec:LR}

    The second alternative approach considered here is based on a less general definition of reversibility referred to as \textit{lag-reversibility}, which involves bi-variate joint probability distributions instead of the more complex multivariate distributions~\cite{Lawrance1991DirectionalityAR, zorzetto2018extremes}. Specifically, lag-reversibility implies that a stationary signal $u(t_i)$ is time reversible if the quantity

    \begin{equation}
      \rho_c(\tau)=\left|\rho\left[u(t)^2,u(t+\tau)\right]-\rho\left[u(t),u(t+\tau)^2\right]\right| \label{SM_eq:lag_rev}
    \end{equation}
    
    \noindent is $\rho_c(\tau)=0$ for any time lag $\tau$~\cite{Lawrance1991DirectionalityAR}, where $\rho\left[\bullet,\bullet\right]$ is the Pearson correlation coefficient. In Eq.~(\ref{SM_eq:lag_rev}) we introduced the absolute value although it is not strictly required for $\rho_c(\tau)=0$ to be satisfied~\cite{Lawrance1991DirectionalityAR}, because we are interested in testing time irreversibility that is $\rho_c\neq 0$. It should be noted that $\theta$ (Section~\ref{SM_subsec:FT}) and $\tau$ represent two different types of time intervals. In fact, while $\theta$ represents the time interval used to calculate the time derivative $\Delta u$, hence acting as a parameter that under-samples the signal for increasing $\theta$ values, $\tau$ is a time shift between two signals in the evaluation of the correlation coefficient (hence, without any under-sampling effect). Moreover, a consequence of Eq.~(\ref{SM_eq:lag_rev}) is that $\rho_c=0$ for $\tau=0$. 
    
    
    \begin{figure}[t]
        \centering
        \includegraphics[width=0.9\linewidth]{SM_fig5-01.eps}
        \caption{\textbf{Lag-reversibility-based results}. \textbf{a} Two-dimensional plot of $\widetilde{\rho_c}=\rho_c/\rho_{c,max}$ as a function of the wall-normal coordinate $y^+$ and the time shift $\tau^+$. \textbf{b} A zoomed-in region of panel \textbf{a} -- highlighted as a red dashed box -- but with a re-scaled colorbar.}
        \label{SM_fig5}
    \end{figure}
    

    Figure~\ref{SM_fig5}a shows $\widetilde{\rho_c}=\rho_c/\rho_{c,max}$ as a function of the wall-normal coordinate $y^+$ and the time shift $\tau^+$. Similar to the previous analyses (e.g., Fig.~\ref{SM_fig4}), larger deviations from $\rho_c=0$ are found in the buffer layer at smaller $\tau^+$, as well as along a bulge-like pattern in the log-layer, thus corroborating the main outcomes of the present analysis. However, the lag-reversibility measure of TI, $\widetilde{\rho_c}$, fails to clearly identify larger TI levels in the intermittency region (dashed red region in Fig.~\ref{SM_fig5}a), which are only detectable by significantly reducing the upper limit of the colorbar as shown in the inset in Fig.~\ref{SM_fig5}b. Therefore, the time-irreversibility patterns found via the lag-reversibility approach (i.e., $\rho_c$) further support the investigation carried out in the main text though visibility networks (i.e., $I_k$), as well the fluctuation theorem-based analysis in Section~\ref{SM_subsec:FT} (i.e., $\langle Z_2^*\rangle$). However, differently from $I_k$ and $\langle Z_2^*\rangle$ which are based on a nonlinear measure (i.e., the Kullback-Leibler divergence), $\rho_c$ relies on Pearson correlation which is a measure of linear (temporal) dependence in the signal. This might explain why high irreversibility levels are not detected in the intermittency region (Fig.~\ref{SM_fig5}a). 
    
    
    In conclusion, results described in the main text via the visibility-based approach are here corroborated, as all three methodologies provide concordant outcomes in terms of TI levels at different scales and wall-normal positions, in spite of the different nature and sensitivity to TI of each methodology.
 
    

 
