# Non-unique self-similar blowups in shell models: insights from dynamical systems and machine-learning

Ciro Campolina, Eric Simonnet, Simon Thalabard

► To cite this version:

Ciro Campolina, Eric Simonnet, Simon Thalabard. Non-unique self-similar blowups in shell models: insights from dynamical systems and machine-learning. Journal of Physics A: Mathematical and Theoretical, 2025, 58 (17), pp.175701. 10.1088/1751-8121/adc773. hal-05407713

HAL Id: hal-05407713

https://hal.science/hal-05407713v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
**PAPER • OPEN ACCESS**

# Non-unique self-similar blowups in shell models: insights from dynamical systems and machine-learning

To cite this article: Ciro Campolina et al 2025 J. Phys. A: Math. Theor. 58 175701

View the [article online](#) for updates and enhancements.

## You may also like

- *Sub-grid effects of the Voigt viscoelastic regularization of a singular dyadic model of turbulence*
Rafael Borges and Fabio Ramos

- *RG approach to the inviscid limit for shell models of turbulence*
Alexei A Mailybaev

- *A Consistent Set of Empirical Scaling Relations for Spiral Galaxies: The $(v_{\text{max}}, M_{\text{BH}})-(M_{\text{BH}})$ Relations*
Benjamin L. Davis, Alister W. Graham and Françoise Combes
# Non-unique self-similar blowups in shell models: insights from dynamical systems and machine-learning

Ciro Campolina<sup>1,2</sup> and Eric Simonnet<sup>2</sup> and Simon Thalabard<sup>2,*</sup>

<sup>1</sup> Calisto, Inria, Université Côte d'Azur, 2004 route des Lucioles, 06902 Sophia Antipolis, France

<sup>2</sup> Institut de Physique de Nice, Université Côte d'Azur CNRS—UMR 7010, 17 rue Julien Lauprêtre, 06200 Nice, France

E-mail: simon.thalabard@univ-cotedazur.fr

Received 13 January 2025; revised 10 March 2025

Accepted for publication 31 March 2025

Published 23 April 2025

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Strong numerical hints exist in favor of a universal blowup scenario in the Sabra shell model, a popular cascade model of 3D turbulence, which features complex velocity variables on a geometric progression of scales $l_n \propto \lambda^{-n}$. The blowup is thought to be of self-similar type and characterized by the finite-time convergence towards a universal profile with non-Kolmogorov (anomalous) small-scale scaling $\propto l_n^x$. Solving the underlying nonlinear eigenvalue problem has however proven challenging, and prior insights mainly used the Dombre–Gilson renormalization scheme, transforming self-similar solutions into solitons propagating over infinite rescaled time horizon. Here, we further characterize Sabra blowups by implementing two strategies targeting the eigenvalue problem. The first involves formal expansion in terms of the bookkeeping parameter $\delta = (1-x)\log\lambda$, and interpretes the self-similar solution as a (degenerate) homoclinic bifurcation. Using standard bifurcation toolkits, we show that the homoclinic bifurcations identified under finite-truncation of the series converge to the observed Sabra solution. The second strategy uses machine-learning optimization to solve directly for the Sabra eigenvalue. It reveals an intricate phase space, with the presence of a continuous family of non-universal blowup profiles, characterized by various number $N$ of pulses and exponents $x_N \ge x$.

Keywords: blow-up, asymptotic self-similarity, cascade models, turbulence

<figure><img src="image_3.png" /></figure>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
# 1. Introduction

The blowup problem in mathematical fluid dynamics addresses the settlement of non-Lipschitz singularities within a finite-time horizon, out of sufficiently smooth initial data. The presence of blowups and their nature are believed to be key for the classification of a variety of commonly observed statistical features of turbulent systems, from power laws and multiscaling [1, 2] to stochastic predictability [3, 4]. With the 3D Euler as a paradigmatic case, the direct mathematical evidence of blowup from generic initial conditions is however a notoriously hard problem [5, 6], with even direct numerical observations leading to conflicting insights [7–12].

Presumably, one of the simplest blowup scenarios is that of asymptotic self-similarity, which describes the finite-time algebraic convergence towards a universal self-similar profile, with specific scaling exponent $x$ prescribed as a non-linear eigenvalue problem [6, 13]. In the 3D Euler setting, an explicit self-similar axially symmetric solution was recently constructed through machine-learning optimization techniques, pointing towards the relevance of the self-similar scenario at least for certain specific flow configurations [14]. Self-similar blowups were also identified in the simplified setting of cascade models, including in particular diffusion-approximation (Leith) [15, 16] and shell model representations of fluid dynamics [17–23]. In Leith models, self-similar attractors were explicitly identified in terms of homoclinic bifurcations of low-order but likely intricate dynamical systems [24–26], generically prescribing non-Kolmogorov (anomalous) scaling exponent. In shell models, following the seminal work by Gilson et al [1], the identification of self-similar blowups relied on numerical simulations exploiting formulations in terms of similarity variables [6, 27], which in essence map self-similar blowups onto travelling waves propagating within an infinite time horizon, whose speed is prescribed by the eigenvalue $x$. A later work [28] approached the eigenvalue problem by reformulating it as a variation procedure for approximated self-similar solitary peaks (solitons). Studying local minima of the associated cost functionals, the contribution of trajectories near the self-similar soliton (called quasisolitons) was considered for the asymptotic multiscaling nature of turbulent solutions.

The purpose of the present paper is to revisit the existence of self-similar blowup solutions in shell models [17], by implementing two strategies directly addressing the underlying non-linear eigenvalue problem. For the sake of clarity, we here focus on the case of the Sabra shell model [29], known to exhibit self-similar blowups over a wide range of parameters [23]. Still, our results also immediately hold for the GOY model [30, 31], and our techniques can be extended and applied to other shell models. The two distinct strategies provide complementary views on the complexity of the Sabra blowup. The first strategy is perturbative and involves a formal expansion in terms of the bookkeeping parameter $\delta = (1-x)\log\lambda$. Such expansion yields a hierarchy of associated dynamical systems, for which the appearance of a self-similar solution is interpreted as a (degenerate) homoclinic bifurcation. Using standard numerical tools of bifurcation theory, we show that under any finite-truncation of the expansion series, the homoclinic bifurcation can be tracked down from a classical Hopf bifurcation scenario, with fast convergence to the observed Sabra solution. The second strategy is non-perturbative: it uses machine-learning optimization to solve directly the Sabra eigenvalue problem. This approach reveals an intricate phase space, with the presence of a continuous family of non-universal blowup profiles.

The paper is organized as follows. Section 2 recalls the general Sabra dynamics and its Dombre–Gilson formulation, and formulates the eigenvalue problem. Section 3 describes the perturbative approach and interpretation of the similarity profile as a global bifurcation. Section 4 describes the non-perturbative solution using machine learning optimization. Section 5 formulates concluding remarks. Some additional material is left to appendices.
## 2. Sabra dynamics and the eigenvalue problem

### 2.1. Sabra shell model

Shell models of turbulence [17, 32] are infinite dimensional dynamical systems of continuous time $t \ge 0$, whose variables are complex velocities $u_n(t)$, $n \in \mathbb{Z}$, associated to a geometric progression of scales $\ell_n = k_n^{-1} = \lambda^{-n}$, with $\lambda > 1$ the intershell ratio. The specific form of the equations is chosen to mimic the spectral Navier–Stokes equations and to keep some of their main properties, like the symmetries and the inviscid invariants of motion. Many different formulations have been proposed [30, 31, 33, 34]. We refer the reader to [35] for a comprehensive list and deduction of shell models in the context of magnetohydrodynamics. In this paper, we shall consider the inviscid Sabra model [29], which prescribes a cascade dynamics with nearest and next-nearest shell interactions as

$$ \dot{u}_n = ik_n [\lambda u_{n+2} u_{n+1}^* - (1+c)u_{n+1}u_{n-1}^* - c\lambda^{-1}u_{n-1}u_{n-2}], \quad n \in \mathbb{Z}, \quad (1) $$

where i is the imaginary unit, and the upper-script stars stand for complex conjugation. The coefficient c is here prescribed to be negative, and we later write it as $c = -\lambda^{-g}$ with g a positive constant. In this form, Sabra dynamics exhibits the two inviscid conserved quantities $E = \sum_n |u_n|^2$ and $C = \sum_n c^{-n}|u_n|^2 = \sum_n (-1)^n k_n^g |u_n|^2$. The first invariant is the kinetic energy, and the second is often called kinetic helicity for not been a sign-defined quantity. Classical Sabra model is associated to the choice $g = 1$, corresponding to $c = -\lambda^{-1}$. In this case, the second invariant can be rewritten as $C = \sum_n (-1)^n |u_n\omega_n|$, by introducing the vorticities $\omega_n = ik_n u_n$, and thus C has the same dimension as the full Euler helicity.

Forced and dissipated Sabra shares a number of features with the 3D Navier Stokes equations. The constant energy flux solution yields Kolmogorov scaling $u_n \propto l_n^{1/3}$, and exact inertial laws, akin to Kolmogorov’s four-fifths law, also hold [36, 37]. Most importantly, the structure functions feature algebraic decay with anomalous exponents [36, 38, 39], signaling a multifractal behavior similar to the heuristics of Navier–Stokes turbulence [29]. The strong reduction in degrees of freedom however allows for accurate numerical estimations of arbitrarily high-order structure functions [40, 41].

### 2.2. Finite-time blowup

The main feature of Sabra that we study in this paper is the development of singularities in finite time from regular initial data. By regular data, we mean that the values of $u_n$ decay sufficiently fast with respect to $n$, say the norm $\|u\| = (\sum_n k_n^2 |u_n|^2)^{1/2} < \infty$. In this case, local-in-time well-posedness theorems [42] assure that a unique solution exists, and it remains regular over a maximal interval of time $[0, t_*)$. In principle, either the solution keeps the regularity for arbitrarily large times, i.e. $t_* = \infty$, or there is a singularity at $t_* < \infty$, after which the solution cannot be continued in the same regularity class. In the latter case, we necessarily have the blowup of the norm $\|u(t)\| \to \infty$ as $t \to t_*$. More specifically, a Beale–Kato–Majda blowup criterion [42] shows that $\int_0^{t_*} \sup_n |\omega_n| dt = \infty$ for the vorticities $\omega_n = ik_n u_n$ in the case of finite-time singularity. For instance, this implies that if $\sup_n |\omega_n| \sim (t_* - t)^{-\beta}$, we must have $\beta \ge 1$. While the finite-time blowup has been proved for other shell model formulations [21, 43], there are no rigorous statements about whether $t_*$ is finite or not for Sabra. Nevertheless, all previous numerical works [3, 19, 23, 44, 45] attest the existence of a finite-time singularity $t_* < \infty$.

Figure 1 illustrates an archetypical numerical integration of the classical Sabra model. The solution initially restricted to small $n$ extends to ultraviolet scales by developing the power
<figure><img src="image_4.png" /><figcaption>Figure 1. Direct numerical integrations of classical Sabra dynamics (1), $\lambda=2$ and $c = -\lambda^{-1}$. For the computations, we simulate only the points $n \ge 0$, and we use the boundary conditions $u_{-1} = u_{-2} = 0$. The initial data is nonzero only for the two first points $u_0 = -i$ and $u_1 = -0.5i$. Left: shell velocities $|u_n|$ in log scale as function of $n$ for different instants of time approaching the blowup time $t_* \approx 1.243$. The solution develops a power law $|u_n| \propto k_n^{-x}$ with $x \approx 0.281$. Right: the velocity self-similar profile $\mathcal{U}(\eta)$ computed from the numerical integrations using formula (4). Colors correspond to different instants of time. The profile displays the infra-red power law $\mathcal{U} \propto \eta^{-x}$ and the ultra-violet sharp front.</figcaption></figure>

law $|u_n| \propto k_n^{-x}$ in the finite time $t_* \approx 1.243$. The scaling exponent is empirically estimated as $x \approx 0.281$, and it was previously found in literature [23]. Numerical computations also reveal the asymptotics $\|u\| \sim (t_* - t)^{-1}$ and $\sup_n |\omega_n| \sim (t_* - t)^{-1}$ as $t \to t_*$, which are compatible with blowup criteria. Most importantly, one verifies that the shape of the spectra front extending from infrared to ultraviolet is preserved up to rescaling, suggesting self-similar dynamics.

## 2.3. Dombre-Gilson renormalization

Self-similarity is better characterized when considering the renormalization scheme introduced by Dombre and Gilson (DG) [19]. Strictly speaking, the DG scheme introduces the logarithmic time and the rescaled vorticities

$$t_* - t \mapsto s = -\log(t_* - t), \quad u_n \mapsto w_n = ik_n u_n e^{-s}. \qquad (2)$$

Such transformation maps the finite-time blowup to an infinite-horizon essentially-bounded dynamics. Sabra equations (1) are then written in those new variables as

$$\frac{dw_n}{ds} = -w_n - [\lambda^{-2}w_{n+2}w_{n+1}^* - (1+c)w_{n+1}w_{n-1}^* + c\lambda^2w_{n-1}w_{n-2}]. \quad (3)$$

Numerical simulations show that the self-similar blowups in the original Sabra equations correspond to traveling wave solutions in these new renormalized variables [19, 23]. Figure 2 illustrates such traveling waves. The scaling exponent $x$ of the blowup solution depends upon the parameters of the problem, namely the intershell ratio $\lambda$ and the $g$ parameter prescribing the second inviscid invariant. More interestingly, the wave profile may vary for different parameters. Numerical experiments suggest the existence of at least one traveling wave for each pair of values $(\lambda, g)$. Figure 3 compares distinct solutions obtained from three different choices of parameters.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: Sabra dynamics for the Dombre–Gilson renormalized vorticities $|w_n|$ as function of $n$ for different logarithmic times $s$, which correspond exactly to the instants of time $t$ displayed in figure 1. Right: the homogeneous vorticity self-similar profile $W(\tau)$ computed from the direct numerical simulations using formula (8).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Typical blow-up profiles obtained by integrating (3) in log scale and renormalized to have unit $\ell^2$ norm. The DG solutions are shown as a function of $n$ for some fixed time well after transients. All these solutions are indeed travelling waves. Three cases are shown: the red curve corresponds to $(\lambda, g) = (7.389, 0.2)$, the black curve is for $(\lambda, g) = (2, 1)$, and the blue curve is $(\lambda, g) = (1.221, 2)$.</figcaption></figure>

Self-similarity is related to the fact that the waves preserve exactly the profile in their constant speed propagation, while the speed of propagation itself can be related to the scaling exponent $x$ (see [45] for more details). The wave profile is recovered by considering a Poincaré map defined over a set of successive discrete instants of renormalized time $s_n$ for which the center of mass of the wave lies on the nodes $n$. In this case, the self-similar blowup corresponds to a fixed-point attractor of the Poincaré map, which is precisely the wave profile. Using this approach to other shell models reveals rather intricate blowup scenarios, like periodic, quasi-periodic and chaotic blowups [46]. Nevertheless, the blowup of Sabra was verified to be self-similar for a wide range of parameters. It is even conjectured that the blowup in Sabra might be universal, as a manifestation of the renormalization group formalism [23]. In dynamical systems terminology, this would correspond to a fixed point attractor with non-zero measure basin of attraction. These particular statements are in the core of our investigations on the next sections.
Particularly, we shall support the existence of multiple unstable self-similar blowups in Sabra, which have not been yet reported in literature, up to our knowledge. This reveals new insights about the important problem of universality in Sabra finite-time blowup.

## 2.4. Sabra self-similarity eigenvalue problem

While the DG scheme proves a very useful tool to classify the different types of blowups and to track their bifurcations in general shell models, here we focus only on the self-similar scenario in Sabra. In this case, the usual and most effective approach is to deduce governing equations directly for the self-similar profile [47], which we do now. The strategy is to write a general ansatz for the profile together with the proper boundary conditions. Then we undergo some variable transformations to make the system autonomous. Such procedure leads to a nonlinear eigenvalue problem, whose eigenvector is the renormalized self-similar profile, and the eigenvalue is the scaling exponent x. For simplicity, we here focus on the case of purely imaginary variables $u_n \in i\mathbb{R}$. For this reason, our results equally apply to the GOY shell model [30, 31], which differs from Sabra only on the complex conjugations of the nonlinear term. Specifically, we consider nontrivial solutions of the form

$$u_n(t) = -i(t_* - t)^{\frac{x}{1-x}} \mathcal{U}(\eta_n), \quad \text{for } \eta_n = (t_* - t)^{\frac{1}{1-x}} k_n. \qquad (4)$$

The minus sign is conventional. The real-valued function $\mathcal{U}(\eta)$ is the self-similar profile, which couples the infrared power-law scaling

$$\mathcal{U} \propto \eta^{-x} \quad \text{as } \eta \to 0, \qquad (5)$$

to a sharp ultraviolet front

$$\mathcal{U} \to 0 \text{ faster than any power law as } \eta \to \infty. \qquad (6)$$

We draw the self-similar profile $\mathcal{U}(\eta)$ from direct numerical integrations in figure 1. Ansatz (4) and boundary conditions (5) and (6) are motivated by the scaling symmetry of Sabra and the empirical evidence of an infrared power-law scaling $|u_n| \propto k_n^{-x}$ extending up to an ultraviolet sharp front—see figure 1 for the numerical picture. We remark that the scaling exponent is necessarily bounded above, $x < 1$. With this bound, equation (4) describes a front $k_* \propto (t_* - t)^{-1/(1-x)}$ propagating to the UV in finite-time as $t \to t_*$. Under the ansatz (4), Sabra dynamics (1) prescribes the equation for the profile $\mathcal{U}(\eta)$ as

$$x\mathcal{U}(\eta_n) + \eta_n \mathcal{U}'(\eta_n) = (1-x)\eta_n \left[ \lambda\mathcal{U}(\eta_{n+2})\mathcal{U}(\eta_{n+1}) - (1+c)\mathcal{U}(\eta_{n+1})\mathcal{U}(\eta_{n-1}) \right. \\ \left. \qquad + c\lambda^{-1}\mathcal{U}(\eta_{n-1})\mathcal{U}(\eta_{n-2}) \right]. \qquad (7)$$

Together with the boundary conditions (5) and (6), this is the fundamental boundary value problem (BVP) for the self-similar blowup profile $\mathcal{U}(\eta)$.

For further analysis, it proves convenient to recast equation (7) into an autonomous form. We achieve this upon introducing the logarithmic similarity variable and the rescaled vorticity profile

$$\tau_n = \log \eta_n, \quad \mathcal{W}(\tau_n) = e^{\tau_n} \mathcal{U}(e^{\tau_n}). \qquad (8)$$
After rescaling the new variable $\tau_n \mapsto \tau_n/(1-x)$, the equation (7) becomes the autonomous dynamical system

$$W'(\tau_n) = W(\tau_n) + \lambda^{-2}W(\tau_{n+2})W(\tau_{n+1}) - (1+c)W(\tau_{n+1})W(\tau_{n-1}) \\ + c\lambda^2W(\tau_{n-1})W(\tau_{n-2}), \quad (9)$$

subject to the boundary conditions

$$W \underset{-\infty}{\sim} e^{\tau} \to 0 \quad \text{and} \quad W \underset{\infty}{\sim} 0. \quad (10)$$

Observe that the logarithmic transformation (8) maps shifts in $n$ into additions to $\tau_n$ as

$$\tau_{n+j} = \tau_n + j\delta, \quad j \in \mathbb{Z}, \quad (11)$$

for the parameter

$$\delta = (1-x) \log \lambda. \quad (12)$$

Figure 2 plots the profile $W$ computed from direct numerical integrations of Sabra.

In principle, Sabra dynamics prescribes the profile $W$ only at the discrete scales $\tau_n$. This poses difficulties, as the left-hand side of (9) involves evaluating the derivative $W'$. To analyze the Sabra BVP, we therefore extend the problem to the full line, i.e., we consider the profile $W(\tau)$ for the continuous variable $\tau \in \mathbb{R}$. We remark that this is well justified since such continuous limit does not require estimating differential operators. Finally, the Sabra BVP becomes the functional fixed-point problem

$$\dot{W} = \mathcal{F}[W], \mathcal{F}[W] = W + \lambda^{-2} W_{2\delta} W_{\delta} - (1+c) W_{\delta} W_{-\delta} + c\lambda^2 W_{-\delta} W_{-2\delta}, \quad (13)$$

with the shorthand $W_{j\delta}(\tau) = W(\tau + j\delta)$, $j \in \mathbb{Z}$, the parameter $\delta$ given by (12), and the boundary conditions (10). Equation (13) can also be interpreted as a nonlinear eigenvalue problem, whose eigenvector is the renormalized self-similar profile $W$, and the associated eigenvalue is the scaling exponent $x$. This is the main problem we shall analyze in this paper.

Some remarks should be addressed. We point out that the solution of the autonomous BVP gives the profile of the traveling wave obtained from the DG scheme. Indeed, the dynamics featured in equation (13) is recovered from DG variables (3) by considering the travelling wave solution

$$w_n(\tau) = W(\tau_n), \quad \tau_n = (1-x) \log \lambda(n-n_*), \quad n_*(\tau) = \frac{\tau}{(1-x) \log \lambda}, \quad (14)$$

and identifying $w_{n+j}(\tau) = W(\tau_{n+j}) = W_j(\tau_n)$—recalling the shorthand $W_{j\delta}(\tau) = W(\tau + j\delta)$. As such, the DG rescaling scheme, be it obtained from the change of variables (8) or from (2), maps the similarity ansatz onto a traveling wave solution $W(\tau_n)$, which propagates in shell space over the infinite time $\tau$ horizon. In the end, the final BVP problem (13) emerges as a continuous limit (i.e. taking $\tau_n$, $n \in \mathbb{N} \mapsto \tau \in \mathbb{R}$) of either a self-similar ansatz or a Dombre–Gilson traveling wave solution. Figure 4 illustrates the correspondence between the different approaches.

From a topological viewpoint, the origin $W \equiv 0$ is a fixed point of $\mathcal{F}$. The vanishing boundary conditions (10) at $\tau \to \pm\infty$ imply that the similarity profile $W$ is obtained as a homoclinic trajectory of the homogeneous functional dynamics (13). At this point, the difficulty lies in the fact that (13) is not a classical delay differential equation (see e.g. [48]), but rather a general Riccati-like functional equation mixing positive and negative time shifts. We are not aware of any general theory for the analysis of such kind of system.
<figure><img src="image_7.png" /><figcaption>Figure 4. Starting from Sabra shell model (1): correspondence between Dombre-Gilson (3) and our BVP problem (13).</figcaption></figure>

Our main goal in the next sections is to study existence, uniqueness and stability of blowup solutions in Sabra through the autonomous BVP (13). In the absence of general analytical theories, we approach this problem using two distinct strategies. The first strategy consists in Taylor-expanding the shifted functions $W_{j\delta}$ with respect to the bookkeeping parameter $\delta$, and then use dynamical systems tools. This perturbative approach is reported in section 3. The second strategy consists in solving directly the problem without any a priori simplification. For this, we use machine learning tools. We address this non-perturbative approach in section 4.

### 3. Perturbative analysis & Sabra hierarchy

The main difficulty of the similarity equation resides on the field time shifts. To overcome this, we consider a systematic perturbative approach. We undertake Taylor expansions of the shifted profiles $W_{j\delta}$ in equation (13) with respect to the bookkeeping parameter $\delta$. This yields a hierarchy of dynamical systems whose bifurcations describe Sabra blowup solutions. We start by the local analysis of the main fixed point equation.

#### 3.1. Fixed points and Hopf bifurcation

Local analysis of Sabra BVP (13) identifies two elementary fixed points $\mathcal{F}[W] = 0$, namely the constant profiles

$$W_0 = 0 \quad \text{and} \quad W_H = \frac{\lambda^2}{(\lambda^2 - 1)(1 - c\lambda^2)} > 0. \qquad (15)$$

Direct calculations reveal that $W_0$ behaves like a saddle. We are interested in studying the stability of $W_H$ when $x$ varies. For this, we are able to show rigorously the following result, that we state as a theorem.

**Theorem 1.** Assume $c < 0$, and $\lambda > 1$. The fixed point $W_H$ loses its stability through an Hopf bifurcation for the critical exponent $x_{Hopf}$

$$x_{Hopf} = 1 - \frac{(\lambda^2 - 1)(1 - c\lambda^2)}{(1 - c\lambda^4)\log\lambda} \frac{\arccos\rho}{(1 + 2\rho)\sqrt{1 - \rho^2}} < 1, \qquad (16)$$
with frequency

$$ \xi_{\text{Hopf}} = \frac{1 - c\lambda^4}{(\lambda^2 - 1)(1 - c\lambda^2)} (1 + 2\rho) \sqrt{1 - \rho^2}, \quad (17) $$

where $\rho$ is the unique root $|\rho| \le 1$ of $\rho^2 + (\frac{1}{2} - R)\rho + \frac{1}{2}R - 1 = 0$, with $R = \frac{c+1}{c\lambda^2 + \lambda^{-2}}$. The solution is stable when $x > x_{\text{Hopf}}$ and unstable otherwise.

The proof is given in the appendix A. The remarkable aspect of this result is that it applies for all values of $c < 0$ and $\lambda > 1$. Another striking fact is that, in the limit $\lambda \to 1^+$, the Hopf bifurcation parameters $x_{\text{Hopf}}$ have a well-defined finite limit, with the frequencies $\xi_{\text{Hopf}}$ growing to infinity.

System (13) appears to always follow the same bifurcation sequence when $x$ decreases:

$$ \text{stable fixed point} \to \text{Hopf bifurcation} \to \text{stable limit cycle} \to \text{homoclinic orbit.} $$

However, in its continuous functional form, the main difficulty in the analysis of the Sabra BVP (13) comes from the presence of the $x$-dependent time shifts through integer multiples of $\delta = (1-x)\log\lambda$. This constitutes a class of dynamical systems for which (infinite-dimensional) bifurcation theory might be non-standard.

## 3.2. The Sabra hierarchy

The idea of the perturbative approach is to consider expansions of the Sabra BVP in terms on the bookkeeping parameter $\delta = (1-x)\log\lambda$, in order to interpret its solution as a global (homoclinic) bifurcation arising in an infinite-dimensional dynamical system. We show below that this system can be truncated and analyzed with standard dynamical system tools at any finite order $M$.

Explicitly, when truncated at order $O(\delta^M)$, the Sabra BVP (13) maps to the dynamical system

$$ \left\{ \begin{array}{l} X_k'(\theta) = X_0 X_{k+1}, \quad (0 \le k \le M-2), \quad X_{M-1}'(\theta) = \frac{1}{\sigma_{0M}} \mathcal{G}_M[\mathbf{X}], \\ \text{with } \mathcal{G}_M[\mathbf{X}] = -X_0 + \frac{X_1}{\delta} - \frac{1}{2} \sum_{k=0}^{M-1} \sum_{\substack{i+j=k \\ 0 \le i,j \le k}} \sigma_{ij} X_i X_j - \frac{1}{2} \sum_{\substack{i+j=M \\ 1 \le i,j \le M-1}} \sigma_{ij} X_i X_j. \end{array} \right. \quad (18) $$

It is formulated in terms of the renormalized time $\theta$ and the rescaled vorticity derivatives $\mathbf{X}$ as

$$ \mathbf{X} = (X_0 \dots X_{M-1}), \quad X_k(\theta) = \delta^k \mathcal{W}^{(k)}(\tau), \quad \theta = \frac{1}{\delta} \int_0^\tau \frac{d\tau'}{\mathcal{W}(\tau')}, \quad (19) $$

and featuring the symmetric coefficients

$$ \sigma_{ij} = \frac{\lambda^{-2} (2^i + 2^j)}{i!j!} \left( (-1)^{i+j} c\lambda^4 - \frac{(-1)^i + (-1)^j}{2^i + 2^j} (1+c)\lambda^2 + 1 \right). \quad (20) $$

Besides, the boundary conditions (10) prescribe

$$ \mathbf{X} \to 0 \quad \text{as} \quad \theta \to \pm\infty \quad (21) $$

for the similarity solutions. We refer the reader to the appendix B for a detailed derivation of (18)–(20). System (18) defines a hierarchy of dynamical systems of increasing dimension
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: Hopf bifurcation parameters $x_{\text{Hopf}}$ as a function of $\lambda$, for $M$ from 2 to 7. Here, $g=1$ is kept fixed. The value $M = \infty$ is computed from the exact formula (16). Right: Hopf parameters in the plane $(\lambda, g)$ for $M = \infty$ using the exact formula (16). The white curve indicates the reference case $\lambda = 2, g = 1$. Note that the $\lambda$ range is different in the left and right panels.</figcaption></figure>

$M$, which formally retrieves the Sabra BVP (13) when $M \to \infty$, and which we later refer to as the Sabra hierarchy. The low-order systems can be explicitly given as

$$ \begin{cases} M = 2 & X'_{0} = X_{0}X_{1}, \sigma_{02}X'_{1} = -X_{0} + \frac{X_{1}}{\delta} - \frac{1}{2}\sigma_{00}X_{0}^{2} - \sigma_{01}X_{0}X_{1} - \frac{1}{2}\sigma_{11}X_{1}^{2}, \\ M = 3 & X'_{0} = X_{0}X_{1}, X'_{1} = X_{0}X_{2}, \sigma_{03}X'_{2} = G_{2} - \sigma_{02}X_{0}X_{2} - \sigma_{12}X_{1}X_{2} \\ M = 4 & X'_{0} = X_{0}X_{1}, X'_{1} = X_{0}X_{2}, X'_{2} = X_{0}X_{3}, \sigma_{04}X'_{3} = G_{3} - \sigma_{03}X_{0}X_{3} - \sigma_{13}X_{1}X_{3} - \frac{\sigma_{22}}{2}X_{2}^{2} \end{cases} \quad (22) $$

and so on, where we have used the shorthand $G_M = \mathcal{G}_M[\mathbf{X}] = \sigma_{0M}X'_{M-1}$, satisfying

$$ G_{m+1} = G_m - \sigma_{0M}X_0 X_M - \frac{1}{2} \sum_{\substack{i+j=M+1 \\ 1 \le i,j \le M}} \sigma_{ij} X_i X_j. \quad (23) $$

### 3.3. The homoclinic explosion

The hierarchy (18) captures the homoclinic bifurcation relevant for the similarity profiles observed in Sabra direct numerical simulations. Local analysis reveals that all the truncated systems have the same (M-independent) fixed points than (13), namely the origin $\mathbf{X}_0 = (0,0,\dots)$ and $\mathbf{X}_H = (W_H,0,0,\dots)$ — see Expression (15) of $W_H$. However, only the point $\mathbf{X}_H$ is necessarily isolated. The origin is part of a wider set

$$ Z_0^{(M)} = \{ \mathbf{X} = (0, \mathbf{X}_+) \text{, with } \mathbf{X}_+ \in \mathbb{R}^{M-1} \text{ & } \mathcal{G}_M[(0, \mathbf{X}_+)] = 0 \}. \quad (24) $$

This structure gives rise to a classical co-dimension 1 scenario of homoclinic bifurcation, which we identify computationally using standard bifurcation tracking software (PyCont and BifurcationToolkit [49, 50]).

For each $M$, the isolated fixed point $\mathbf{X}_H$ undergoes a Hopf bifurcation at $x_{\text{Hopf}}^{(M)}$. The resulting limit cycle then collides with the invariant set $Z_0^{(M)}$ at $x = x_{\star}^{(M)}$. The critical values $x_{\text{Hopf}}^{(M)}$ and $x_{\star}^{(M)}$ depend upon the order $M$. Nevertheless, in the limit $M \to +\infty$, the value $x_{\text{Hopf}}^{(M)}$ converges towards the theoretical parameter $x_{\text{Hopf}}$ from theorem 1, and $x_{\star}^{(M)}$ converges to the numerically observed exponent $x^*$. We show this in figures 5 and 6. In figure 5, we plot the Hopf bifurcations
<figure><img src="image_9.png" /><figcaption>Figure 6. Left: homoclinic explosion showing period T vs x for various M, and reference parameters λ = 2 and g = 1. Inset shows the limit cycles tracked for M = 7. Right: critical x in Sabra (1) as a function of both λ and g by direct computations. The white curve corresponds to the reference case with g = 1, λ = 2.</figcaption></figure>

for different Sabra parameters $\lambda$ and $g$. Particularly, one observes a very fast convergence of the truncated parameters $x_{\text{Hopf}}^{(M)}$ to the theoretical one $x_{\text{Hopf}}$ with increasing $M$. The curves are visually indistinguishable for $M = 5, 6, 7$. In figure 6, we illustrate the homoclinic explosion, i.e. the continuous evolution from the limit cycle to the homoclinic orbit as we vary $x$ from $x_{\text{Hopf}}^{(M)}$ to $x_{*}^{(M)}$. The homoclinic parameter $x_{*}^{(M)}$ also quickly converges to the numerically observed values $x^{*}$. The fast convergence suggests that, although of infinite dimensionality, the Sabra blowup is accurately described by a low-order model.

### 3.4. Degeneracies

At any order, the Sabra hierarchy (18) exhibits a continuous set of fixed points: this type of degeneracy is inherited from the singular nature of the Sabra BVP, and also appears in Leith models [24, 26, 51]. The degeneracy is removed through the nonlinear time change (19) but signals that the homoclinic explosion might be degenerate. This is better grasped by considering the low-order models of the Sabra hierarchy, prescribed by equation (22). For $M = 2$, the set of fixed points $\mathcal{Z}_{0}^{(2)}$ is the two-point set

$$ \mathcal{Z}_{0}^{(2)} = \{\mathbf{0}\} \cup \{\mathbf{X}_{c}\}, \quad \mathbf{X}_{c}^{(2)} = \left(0, \frac{2}{\sigma_{11}\delta}\right), \quad \sigma_{11} = 4\lambda^{-2}\left(c\lambda^{4} + \frac{1}{2}(1+c)\lambda^{2} + 1\right), \quad (25) $$

which for the particular shell ratio $\lambda = 2$ yields $\mathbf{X}_{c} \simeq (0, -0.78)$. The left panel of figure 7 shows that as $x$ approaches $x_{*}^{(2)} \approx 0.389$, the limit cycle does not collide with the origin, but rather with $\mathbf{X}_{c}^{(2)}$. This suggests that the homoclinic cycle is made of two heteroclinic components, namely $0 \to \mathbf{X}_{c}^{(2)}$ and $\mathbf{X}_{c}^{(2)} \to 0$. The same type of scenario is seen explicitly for the case $M = 3$. The set of fixed points $\mathcal{Z}_{0}^{(3)}$ identifies with the union of two lines, that is

$$ X_0 = 0, X_1 = 0, X_2 \in \mathbb{R} \quad \text{or} \quad X_0 = 0, X_1 \in \mathbb{R}, \quad \sigma_{12}X_2 = -\sigma_{11}X_1 + \frac{1}{\delta} = 0. \quad (26) $$

Figure 7(right) suggests that the limit cycles intersect $\mathcal{Z}_{0}^{(3)}$ at the point $\mathbf{X}_{c}^{(3)} = (0, 0, 1/(\delta\sigma_{12}))$.
<figure><img src="image_10.png" /><figcaption>Figure 7. Degenerate homoclinic orbits for $M=2$ (left) and $M=3$ (right), resulting from the collision of the limit cycles with the continuum set $\mathcal{Z}_0^{(M)}$.</figcaption></figure>

## 4. Non-perturbative approach using machine-learning

In the previous section, we showed that the perturbative approach captures the main numerically observed blowup features. We want to attest that those observations are indeed properties of the original BVP, and not only of the truncated systems. Therefore, we would like to directly solve the functional equation (13) without going through any perturbative approximation. One cannot rely on a time marching approach since both future and past can influence the state at $\tau$. A simple idea is therefore to formulate (13) as an optimization problem. A similar idea was considered in [28], where a variation principle was formulated for self-similar solitary peaks. There, such solitons were approximated by rational trial functions. In the present work, we parametrize the solution by a deep neural network (NN), an approach in the spirit of the so-called physics-informed-NNs method [52, 53]. In the end, we show that our approach enlarges our understanding of the problem: it recovers the self-similar blowups previously discussed, but also uncovers new similarity profiles, which to our knowledge were not reported before.

The problem consists in finding $(x, W)$ such that $\dot{W} = \mathcal{F}[W]$, where $x$ appears in $\mathcal{F}$ only in the delays, but not in the coefficients. Therefore, it is convenient to understand $W$ as a function of both $\tau$ and $x$. Let us consider the following cost functional

$$C[W] = \int_{\mathbb{R}} \int_{x_{\min}}^{x_{\max}} |\dot{W} - \mathcal{F}[W]|^2 dx d\tau. \quad (27)$$

The corresponding optimization problem formulates as finding the argmin solution, called $W^*$,

$$W^*(\tau, x) = \underset{W \in V}{\operatorname{argmin}} C[W] \quad (28)$$

within the functional space

$$V = \left\{ W : (\tau, x) \mapsto W(\tau, x) \mid W(\tau, x) \underset{\tau \to -\infty}{\sim} C_{-\infty} e^{\tau} \text{ and } \lim_{\tau \to +\infty} W(\tau, x) = 0 \right\}, \quad (29)$$

whose functions satisfy the problem boundary conditions. Here, $C_{-\infty}$ is the real prefactor of the infrared asymptotics, which will be important in the computational formulation. In this
approach, we actually obtain a one-parameter family of solutions indexed by x. Then, we compute for which x the fixed point problem is exactly satisfied. The main difficulty is to prescribe the boundary conditions, i.e. to guarantee that $W \in V$. We explain the details next.

## 4.1. NN parametrization and ansatz

First of all, we affirm that we can set the prefactor $C_{-\infty}$ in the definition of V to be a fixed nonzero constant. Indeed, the equations are invariant to time translations. Therefore, fixing a value is equivalent to choosing some time reference frame. Moreover, a nonzero value filters out the trivial identically zero solution. Hereafter we set $C_{-\infty} = 1$. Next, one must be careful in the representation of the whole $\tau$ line. In practice, we choose $\tau \in [-T, T]$ for $T \gg 1$. We then consider the following neural ansatz

$$W_{nn}(\tau, x; \theta) = \chi_{-}(\tau)\chi_{+}(\tau)(1 + \mathcal{N}(\tau, x; \theta))^{2}, \quad \chi_{\pm}(\tau) = (1 + e^{-(\tau/2\pm\sigma)})^{-1}, \quad (30)$$

where $\mathcal{N}$ is a simple fully-connected feedforward (or perceptron) NN, i.e. $\mathcal{N}(\tau, x; \theta) = (\mathcal{N}_{\text{out}} \circ \mathcal{N}_{N} \circ \dots \circ \mathcal{N}_{1} \circ \mathcal{N}_{\text{in}})(\tau, x; \theta)$, $\mathcal{N}_{k}(y; \theta_{k}) = f(W_{k}y + b_{k})$. Here $\theta = (\theta_{k})_{k=0, \dots, N+1}$, $\theta_{k} = (W_{k}, b_{k})$ are the neural parameters to be found, f is called an activation function, which must be nonlinear, and N is the depth of the NN. The size of the matrices $W_{k}$ determines the so-called capacity: the larger it is, the more precise and time-consuming the method becomes.

Let us address important remarks about the choice of ansatz (30). The boundary conditions are actually relaxed: the NN is free to choose the best decays at $\tau \to \pm\infty$ in accordance with the problem (13). For instance imposing a wrong decay, say $W \sim e^{2\tau}$ will have no effect, and the ending result will recover the correct one $W \sim e^{\tau}$. The main purpose in this formulation is to avoid being attracted directly to the trivial zero solution. That said, we observe that the precise details of the ansatz we adopt are unimportant for solving the problem, and many other formulations are possible. Indeed, the solutions obtained do not (and must not) depend on the ansatz chosen. This was checked many times. As we are going to present below, the formulation we chose yields good results, but we do not claim in any way it is optimal.

In this formulation, the problem consists in finding optimal NN parameters. Therefore, we rewrite the cost functional (27) as an optimization of $\theta$ for the NN ansatz (30) in the form

$$C[\theta] = \int_{-T}^{T} \int_{x_{\min}}^{x_{\max}} |\dot{W}_{nn} - \mathcal{F}[W_{nn}]|^2 (\tau, x; \theta) \, dx \, d\tau. \quad (31)$$

We call $\theta^* = \operatorname*{argmin}_{\theta} C[\theta]$. The NN function $W_{nn}(\tau, x; \theta^*)$ represents a parametrized family of solutions from which we identify the best candidate for the eigenvalue problem by solving $\dot{W} - \mathcal{F}[W] = 0$. This technique enables us to avoid minimizing separately many functionals for every choice of x. It works well in practice provided the interval $[x_{\min}, x_{\max}]$ is small enough. The critical $x^*$ one is looking for is therefore

$$x^* = \operatorname*{argmin}_{x} \log \int_{-T}^{T} |\dot{W}_{nn} - \mathcal{F}[W_{nn}]|^2 (\tau, x; \theta^*) \, d\tau. \quad (32)$$

Numerically, we discretize the cost functionals, and we employ the stochastic gradient descent algorithm ADAM [54]. Further numerical and algorithmic details are given in the appendix C.

We show in figure 8 the results obtained after a few iterations. Here we consider classical Sabra, i.e. $\lambda = 2$ and $g = 1$. The algorithmic stochastic descent is shown together with the relative error as a function of x and the converged solution $W_{nn}$. One can observe the typical profile of the traveling wave in DG with its power-law behavior to the left and its fast decay to the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** One-bump self-similar profile for $\lambda = 2$, $g = 1$ obtained using a machine-learning approach. The background figure corresponds to the evolution of the cost functional during the stochastic gradient descent (ADAM). The upper panel is a measure of the relative error as a function of $x$, by plotting the expression inside the argmin of equation (32). Numerical parameters are $T = 15$, $\sigma = 8$. The neural network has 10 hidden layers with 10 neurons/layer and swish activation functions $f$. Batch size is 100 × 50 points for $(\tau, x)$.</figcaption></figure>

right—compare it, for instance, with the profile from direct numerical integrations of figure 2. Moreover, we also observe that the minimum of the cost functional is achieved in the same scaling exponent $x \approx 0.281$ obtained from the direct numerical integrations. The expression inside the argmin in equation (32) is shown in the insert of figure 8. Those first observations attest the efficacy of the method.

## 4.2. N-bumps solutions

We recall that our machine-learning approach described previously relies on stochastic gradient descent. Due to that, depending on the realization of the noise and/or on changing slightly the hyperparameter setup, one can obtain radically different solutions. We show some of them in figure 9. In particular, one is able to obtain solutions exhibiting an arbitrary number $N$ of bumps, each one being very similar to the reference one-bump solution. For example, one may compare the one-bump solution in figure 7 with the one- and two-bumps solutions in figure 9.

A curious observation is the behavior of a *N*-bumps solution with respect to the distance between the many bumps. When the bumps come closer to each other, the shapes and the peaks of the various bumps may differ from each other, and the global picture of the solution becomes more complex—for instance, see the five-bumps solutions in figure 9.
<figure><img src="image_12.png" /><figcaption>Figure 9. N-bumps self-similar solutions for $\lambda = 2, g = 1$, obtained using a machine-learning approach. The phase portraits in the plane $W, \dot{W}$ are displayed showing the increase of the winding number. As the number of bumps increases, so does the critical $x^*$, as long as the bumps do not remain departed from each other.</figcaption></figure>

Another observation that is correlated with the distance between the bumps is the value of the scaling exponents $x_N^*$ associated to the N-bumps solution. They take values rather close to the reference one-bump exponent $x_1^* \approx 0.281$, but more importantly they always seem to be larger than the ground (eigen)value $x_1^* \le x_N^*$, and they approach this lower bound as the distance between the bumps increase. We explain this observation by the following argument. Consider the BVP (13) for the self-similar profile $W(\tau)$, and make the change of variables $\tau = \delta s$ and $H(s) = \delta W(\delta s)$. Then, the governing equation becomes $\dot{H} - \delta H = \lambda^{-2}H_1H_2 - (1+c)H_{-1}H_1 + c\lambda^2H_{-2}H_{-1}$, where we use the notation $H_j(s) = H(s+j)$ for the shifted functions. In this formulation, the shifts are independent of $x$, which appears in the equation through the coefficient $\delta = (1-x)\log \lambda$. One verifies through the energy balance that the scaling exponent $x^*$ satisfies

$$1 - x^* = \frac{W_H^{-1} \int_{\mathbb{R}} H_{-1} H_1 \, ds}{\log \lambda \int_{\mathbb{R}} H^2 \, ds}, \qquad (33)$$

where $W_H$ is a constant given below in (15). Next, we consider a N-bumps solution $H^{(N)}$ approximately composed by the sum of N one-bump solutions sufficiently departed from each
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Numerical simulations of *N*-bump self-similar solutions in the DG representation. Each row displays the time evolution from a different initial condition.</figcaption></figure>

other. We write $\mathcal{H}^{(N)} \approx \sum_{k=1}^{N} H_k$, where $H_k(s) = \mathcal{H}^{(1)}(s+j_k)$ with $j_1 \ll j_2 \ll \cdots \ll j_N$, such that $H_k H_l \approx 0$ if $k \neq l$. In this case, we have $\int_{\mathbb{R}} \mathcal{H}_{-1}^{(N)} \mathcal{H}^{(N)} \mathcal{H}_1^{(N)} ds \approx N \int_{\mathbb{R}} \mathcal{H}_{-1}^{(1)} \mathcal{H}^{(1)} \mathcal{H}_1^{(1)} ds$ and $\int_{\mathbb{R}} (\mathcal{H}^{(N)})^2 ds \approx N \int_{\mathbb{R}} (\mathcal{H}^{(1)})^2 ds$. Therefore, $x_N^* \approx x_1^*$. See, for instance, the two-bumps solution with reasonably departed peaks in figure 9, for which $x_2^* \approx x_1^*$ within numerical estimation. On the other hand, if the bumps become closer to each other, the different bumps have nontrivial contributions to the nonlinearities in expression (33), and this might lead to distinct exponents. This is, for example, the case of the five-bumps solution displayed in figure 9, whose exponent $x_5^* \approx 0.324$, which is slightly bigger than $x_1^*$. We conjecture that the one-bump solution provides the maximal value of (33). We next discuss the stability of those *N*-bumps solutions.

### 4.3. *Numerical simulations*

Our machine-learning approach reveals strong nonuniqueness of self-similar blowups in Sabra. While the classical one-bump blowup has been verified numerically to be rather stable, a natural question would be whether those multiple-bump profiles are stable or not. We investigate this question through direct numerical integrations of DG Sabra dynamics (3).

We consider initial conditions exhibiting two or three bumps of the same shape, but with distinct intensities. The initial shape of each bump is arbitrarily chosen. We show the time evolution of some of the solutions in figure 10. We observe a relaxation from our arbitrarily chosen initial condition to the *N*-bump profiles found with machine learning—for instance, compare the numerical integrations from figure 10 to the machine learning outputs in figure 9. Such relaxation occurs in the first instants of evolution. After some integration time, we observe the suppression of the *N*-bump configuration by the classical one-bump profile. By fine tuning the intensities of the initial bumps, the profile is kept for longer times. In principle, one could carry the *N*-bump profile for arbitrarily large times, as long as the machine precision allows the fine adjustment of the bumps initial intensities.

This whole scenario suggests that the multiple-bump solutions are saddle points of the BVP (13). The fine-tuning of parameters just guarantee closeness of the initial conditions to the stable manifold, whose contraction proportionates the observed initial relaxation. But perturbations (even the numerical integration) may eventually disturb the system on unstable directions, and this causes the suppression of the *N*-bump profile, and a subsequent convergence to the stable one-bump traveling wave. Therefore, strictly speaking, those multiple-bump solutions do not seem to be stable in the sense of Lyapunov. We claim, however, that more systematic studies are needed for a precise characterization of the stability of such solutions.
## 5. Conclusions

We investigated self-similar finite-time blowups in the inviscid Sabra shell model. To this end, we derived the associated similarity equations and formulated them as a BVP parametrized by the scaling exponent $x$. Since we restricted our analysis to purely imaginary solutions, all our results also immediately hold for the GOY model. The resulting equations involve past and future time shifts, which turns out to be difficult treating analytically. To address this, we employed two distinct and complementary strategies.

First, we considered a perturbative approach by Taylor-expanding the shifts in the equations with respect to the bookkeeping parameter $\delta = (1-x)\log\lambda$. This allowed us to map the similarity equations to a hierarchy of dynamical systems of increasing order. There, the blowup scenario is linked to a robust sequence of bifurcations: fixed point, Hopf bifurcation, limit cycle, and finally a homoclinic orbit, which corresponds to the desired blowup solution. Despite the infinite-dimensional nature of the problem, we verified that the main features of the blowup are accurately described by a low-order truncation of the Sabra hierarchy. Indeed, the 7th order expansion already presented results indistinguishable from the theoretical predicted bifurcations and the numerically computed scaling exponents. Similar to blowups observed in differential (Leith) models of turbulence, blowup solutions in Sabra emerge as codimension-1 global (homoclinic) bifurcations. Unlike Leith blowups [51], the high-dimensionality of the BVP revealed by the Sabra hierarchy a priori allows for less generic bifurcations, meaning that the hypothesis of uniqueness of Sabra blowups should be studied in more depth.

Second, we employed a nonperturbative approach to directly solve the similarity equations without approximations. This was possible by reformulating the problem in an optimization framework, and leveraging standard machine learning techniques. With this approach, we discovered infinitely many distinct blowup solutions characterized by $N$-bump profiles. While the classical one-bump profile was found to be stable, the multiple-bump profiles appeared to lack stability, although they could still be accurately tracked using direct numerical simulations.

We also observed that the scaling exponent of the one-bump solutions $x_{\star}^{(1)}$ is likely to be the minimum of all scaling exponents, i.e. $x_{\star}^{(N)} \ge x_{\star}^{(1)}$ for all $N$. We leave the question open as to whether the multiple-bump profiles can be found in the Sabra hierarchy. The $N$-bump blowup solution would then be associated to multiple homoclinic orbits with winding number $N$, but further investigations are required to better understand the role of such structures in blow-up studies.

Together, the two approaches reveal the nontrivial topology of self-similar blowups in the Sabra shell model. While only one (stable) finite-time blowup is generically observed in numerical simulations, there are actually infinitely many other blowups with more complex profiles.

Their absence from prior reports might be due to their inherent instability. The machine learning approach proves to be an excellent tool for tracking these unstable solutions and these findings open interesting perspectives for the study of more complex models. For instance, one might wonder whether scenarios even more intricate than those described here for Sabra could arise in the 3D Euler equations or related models. It is well-known that self-similar finite-time blowups in 3D Euler cannot occur from smooth initial data with sufficiently strong vorticity decay at infinity [55], necessitating more specific initial conditions [56]. Machine learning approaches could serve as promising tools to identify not only stable blowups but also potentially unstable ones in these more complex equations.

Beyond self-similar blowups, the Dombre–Gilson scheme and its variants have revealed that cascade models could exhibit a variety of blowup scenarios [23], including the possibility
of a chaotic blowup [11]. One may wishfully ask whether combining hierarchy and machine learning approaches could elucidate the roles of unstable similarity solutions in the emergence of those more complex roads to blowups. We leave those interesting matters for future work.

## Data availability statement

The data can be reproduced from the text. The data that support the findings of this study are available upon reasonable request from the authors.

## Acknowledgments

The authors wish to thank J Bec and N Valade for their valuable continuing discussions and debates on this topic. C C was supported by the French National Research Agency (ANR Project TILT ANR-20-CE30-0035), and dedicates this work to Notre Dame de Laghet.

## Appendix A. Local analysis

In the following, we discuss bifurcations results for the system (13) written in the dynamical system form $\dot{W} = \mathcal{F}[W]$. It is important to notice that, as such, one deals with a genuine infinite-dimensional dynamical system since the right-hand side involves delays. For convenience, we will refer as 'time' the variable $\tau$. We first briefly discuss the fixed points.

One is looking for solutions satisfying $\mathcal{F}[W] = 0$ for all time. Since $\dot{W} = 0$, only the constant solutions are relevant and they must obey $W(1 + W(\lambda^{-2} - 1 - c + c\lambda^2)) = 0$. Since $c < 0$ (3D case) and $\lambda > 1$, it gives (15).

### A.1. Linearized dynamics and Hopf bifurcation

We consider here the linearisation of (13) around a fixed point $\bar{W}$. The interest is that one can obtain explicit formula. Writing $W = \bar{W} + \phi$, it is straightforward to obtain at first order in the perturbation $\phi$, the linearized dynamics

$$ \dot{\phi} - \phi = \bar{W} (\lambda^{-2} (\phi_{+\delta} + \phi_{+2\delta}) - (1+c)(\phi_{-\delta} + \phi_{+\delta}) + c\lambda^2 (\phi_{-2\delta} + \phi_{-\delta})), $$

where we use again the shorthand $\phi_{k\delta}(\tau) = \phi(\tau + (1-x)k\log\lambda)$ for the shifted functions with respect to integer $j$ multiples of $\delta = (1-x)\log\lambda$. Applying the Fourier transform yields (with the convention $\hat{\phi}(\xi) = \int_{\mathbb{R}} \phi(\tau)e^{-i\xi\tau} d\tau$),

$$ \Gamma(\xi) \hat{\phi}(\xi) = 0, \quad \Gamma(\xi) = (i\xi - 1) - \bar{W} (\lambda^{-2} (e^{i\delta\xi} + e^{2i\delta\xi}) - (1+c)(e^{-i\delta\xi} + e^{+i\delta\xi}) + c\lambda^2 (e^{-i\delta\xi} + e^{-2i\delta\xi})). \quad (A1) $$

The case $\bar{W} = 0$ corresponds to $\phi(\tau) = C\epsilon^\tau$, indicating that $0$ behaves as a saddle. Note that it is spectrally unstable since the Jacobian is the identity. Let us turn now to the other fixed point $\bar{W} = W_H$. Its stability is dictated by the behavior of the zeros of $\Gamma$ in the complex plane. In particular, a transfer of stability occurs when $\Im\xi = 0$ with nonzero parameter velocity (e.g. parameter $x$). We thus denote $\Omega = \delta\xi \in \mathbb{R}$ and compute the real and imaginary parts of $\Gamma$, namely

$$ \Re\Gamma = -1 - \bar{W} ((\lambda^{-2} - 2(1+c) + c\lambda^2)\cos\Omega + (\lambda^{-2} + c\lambda^2)\cos 2\Omega), \quad (A2) $$
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Poles (black curves) of the function $\xi \to \Gamma(\xi)$, $\xi \in \mathbb{C}$ when $x$ varies. The norm $\log|\Gamma|$ in the background is computed at the ending value of $x$. The Hopf bifurcation corresponds to the crossing of the black curves with the real axis.</figcaption></figure>

and

$$ \Im\Gamma = \delta^{-1}\Omega - \overline{W}(\lambda^{-2} - c\lambda^2)(\sin\Omega + \sin 2\Omega). \quad (A3) $$

The real part is quadratic in term of $\cos\Omega$. We now take $\overline{W} = W_H$ (see (15)) and solve $\Re\Gamma = 0$. Let $\rho_{1,2}$ be the roots of

$$ \rho^2 + \left(\frac{1}{2} - R\right)\rho + \frac{R}{2} - 1 = 0, \quad R := \frac{c+1}{c\lambda^2 + \lambda^{-2}}. \quad (A4) $$

We are looking for a root such that $|\rho| \le 1$. The two roots are $2\rho = R - \frac{1}{2} \pm \sqrt{R^2 - 3R + \frac{17}{4}}$, where the expression inside the squareroot is always positive and larger than 2. The condition for having $|\rho| \le 1$ imposes either $R \le 1$ or $R \ge \frac{1}{3}$. As a consequence, there is always a root such that $|\rho| \le 1$. In fact the two roots satisfy this property provided $R \in [\frac{1}{3}, 1]$. It turns out that this scenario never occurs when $c < 0$. Otherwise, it would implies not one, but at least two bifurcations, together with possible codim-2 bifurcations.

One can then proceed with $\Omega = \arccos\rho$, namely $\xi = \gamma^{-1}\arccos\rho$ and the imaginary part $\Im\Gamma$ gives $\xi = W_H(\lambda^{-2} - c\lambda^2)(\sin\Omega + \sin 2\Omega) = W_H(\lambda^{-2} - c\lambda^2)\sqrt{1 - \rho^2(1 + 2\rho)}$. One is interested in the dependence on the power-law exponent $x$. It appears in the term $\delta$ only in the real part of $\Gamma$ (since $W_H$ does not depend on $x$). In other words, the necessary condition for a bifurcation to occur is

$$ 1 - x_{\text{Hopf}} = \frac{(\lambda^2 - 1)(1 - c\lambda^2)}{(1 - c\lambda^4)\log\lambda} \frac{\arccos\rho}{(1 + 2\rho)\sqrt{1 - \rho^2}}, \quad (A5) $$

with frequency

$$ \xi_{\text{Hopf}} = \frac{1 - c\lambda^4}{(\lambda^2 - 1)(1 - c\lambda^2)} (1 + 2\rho) \sqrt{1 - \rho^2}. \quad (A6) $$

### A.2. The limit $\lambda \to 1^+$

Let us write $\lambda = 1 + \epsilon$ and assume that $c = -\lambda^{-g}$. One easily obtains $1 - x \sim \frac{2\arccos\rho}{(1+2\rho)\sqrt{1-\rho^2}}$ and $R \sim \frac{g}{4-g}$ so that $\rho$ is well defined except maybe for $g = 4$. In fact when $g \to 4$, the roots are
$\{-\infty, \frac{1}{2}\}$. We obtain in this particular case $1 - x_{\text{Hopf}} \to 1 - \frac{2\pi}{3\sqrt{3}}$. The frequency term behaves as $\xi \sim \frac{1}{2\epsilon}(1 + 2\rho)\sqrt{1 - \rho^2}$. It is quite remarkable that although the fixed point $W_{\text{Hopf}} \sim \frac{1}{4\epsilon}$ and the frequency are pushed to infinity when $\epsilon \to 0$, a bifurcation still occurs at a well defined value of $x$ which depends on $g$ using the same formula than before but with $R = \frac{g}{4-g}$.

## Appendix B. The Sabra hierarchy

We here provide a technical derivation of the Sabra hierarchy, and comment on the choice of $\delta = (1 - x) \log \lambda$ as a natural bookkeeping parameter. The hierarchy proposed here is justified only in the case $|\delta| < 1$ and converges to the continuous problem modulo an ad-hoc rescaling.

### B.1. Technical derivation

We derive the Sabra hierarchy (18) by Taylor expanding the shifted functions as

$$W_k(\tau) = \sum_{i=0}^{M} \frac{k^i \delta^i}{i!} W^{(i)}(\tau) + O(\delta^{M+1}), \quad (B1)$$

which in turn yields for the $\mathcal{F}$-functional

$$\mathcal{F}[W] = W^{(0)}(\tau) + \frac{1}{2} \sum_{k=0}^{M} \delta^k \sum_{i+j=k} \sigma_{ij} W^{(i)}(\tau) W^{(j)}(\tau) + O(\delta^{M+1}) \quad (B2)$$

in terms of the symmetric coefficients $\sigma_{ij}$ defined by equation (20). Neglecting the terms $O(\delta^{M+1})$ in equation (B2), the Sabra BVP (13) truncates into the $M$th order system

$$W^{(1)}(\tau) = \mathcal{F}_M[W], \quad \mathcal{F}_M[W] = W^{(0)}(\tau) + \frac{1}{2} \sum_{k=0}^{M} \delta^k \sum_{i+j=k} \sigma_{ij} W^{(i)}(\tau) W^{(j)}(\tau). \quad (B3)$$

Besides, the boundary conditions (10) become

$$\forall k, \beta \ge 0 \quad W^{(k)}(\tau) \underset{-\infty}{\propto} e^{\tau} \to 0, \quad W^{(k)}(\tau) \underset{\infty}{\propto} o(e^{-\beta\tau}) \to 0. \quad (B4)$$

Writing $X_k = \delta^k W^{(k)}$ and $\mathbf{X}(\tau) = (X_0, \dots, X_{M-1})$, equation (B3) recasts into the first order system

$$\left\{ \begin{array}{l} \delta \dot{X}_k(\tau) = X_{k+1}, \quad 0 \le k \le M-2, \quad \sigma_{0M} \delta X_0 \dot{X}_{M-1}(\tau) = \mathcal{G}_M[\mathbf{X}], \\ \mathcal{G}_M[\mathbf{X}] = -X_0 + \frac{X_1}{\delta} - \frac{1}{2} \sum_{k=0}^{M-1} \sum_{\substack{i+j=k \\ 0 \le i,j \le k}} \sigma_{ij} X_i X_j - \frac{1}{2} \sum_{\substack{i+j=M \\ 1 \le i,j \le M-1}} \sigma_{ij} X_i X_j, \end{array} \right. \quad (B5)$$

with boundary conditions $\mathbf{X}(\tau) \to 0$. System (B5) features a term $X_0$ on its left-hand side (lhs) and is therefore singular. To desingularize it, we introduce the new time variable $\theta = \int_0^\tau \frac{d\tau'}{\delta X_0(\tau')}$ as prescribed by equation (19)—in short $d\theta = (\delta X_0)^{-1} d\tau$. This yields the Sabra hierarchy (18). For the boundaries, the change of variable $\tau \to \theta$ obviously prescribes $\theta(\tau) \to -\infty - \infty$. Does it obviously prescribe $\theta(\tau) \to \infty$ as $\tau \to \infty$? One must assume that $W > 0$ ($X_0 > 0$) for all $\tau$, if not the change of variable $\theta$ is not well-defined. The boundary condition at $+\infty$ is $X_0 = o(e^{-\beta\tau})$. It is equivalent to write $X_0 e^{\beta\tau} = g(\tau) \to 0$ at $+\infty$ with $g > 0$. Then, there exist $\tilde{\tau}$ and a constant $c$ such that $g(\tau) < 1$ for all $\tau \ge \tilde{\tau}$ and $g(\tau) > c > 0$ for all $\tau \le \tilde{\tau}$ so that $\theta = \int_0^\tau \frac{e^{\beta s}}{g(s)} ds + \int_{\tilde{\tau}}^\tau \frac{e^{\beta s}}{g(s)} ds \ge c s t + \int_{\tilde{\tau}}^\tau e^{\beta s}$, therefore $\theta \to \infty$ when $\tau \to \infty$.
## B.2. The asymptotics $\delta \to 0$

The asymptotics $\delta \to 0$ describes the case $x \to 1^-$ or the formal limit $\lambda \to 1$. In this asymptotics, the Sabra BVP (13) becomes local and the hierarchy (B5) reduces to the case $M = 2$. In this situation, the Sabra BVP and the hierarchy reduce to

$$ \dot{W} = W + \frac{1}{2}\sigma_{00}W^2, \quad \sigma_{00} = -2W_H^{-1}, \quad (+BC) \qquad (B6) $$

yielding the explicit (logistic) solution $W(\tau) = \frac{1}{W_H^{-1} + C^{-1}e^{-\tau}}$. The logistic solution satisfies the vanishing boundary condition only for $\tau \to -\infty$, and hence is not a solution to the Sabra BVP. It converges to the fixed point $W_H > 0$ for $\tau \to \infty$. Heuristically, we however expect the true solution to the BVP to be close to the logistic solution, a feature indeed observed in the non-perturbative optimization approach. This motivates to consider variations in the small parameter $\delta = (1-x) \log \lambda$ to match the ultraviolet boundary condition, and provides a physical rational for the Taylor expansion in $\delta$.

## Appendix C. Machine-learning approach

We provide here more details on the numerical simulations of section 4, which uses a NN parametrisation. First of all, the cost functional (31) must be discretized. We consider a set of $N_\tau \times N_x$ independent and identically distributed (i.i.d.) random points $(\tau, x)$ in the functional domain $D_{\tau,x} = [-T, T] \times [x_{\min}, x_{\max}]$ with uniform distribution $U([-T, T]) \times U(x_{\min}, x_{\max})$. Then, we use a Monte-Carlo integral approximation for the cost functional

$$ C[\theta] \approx \frac{1}{N_{\tau}} \frac{1}{N_x} \sum_{k=1}^{N_x} \sum_{j=1}^{N_{\tau}} \left| \dot{W}_{nn}(\tau_{j,k}, x_{j,k}; \theta) - \mathcal{F}[W_{nn}(\tau_{j,k}, x_{j,k}; \theta)] \right|^2, \qquad (C1) $$

where $W_{nn} = W(\tau, x; \theta)$ establishes a NN with parameters $\theta$ to be found—recall (30). The gradients of (C1) are found by automatic differentiation (backpropagation). The algorithm follows the simple rules:

* **Hyperparameter setup:** choose $T, \sigma$ and build some NNs with input dimension equals to 2 (for $\tau$) and output dimension equals to 1. Choose $N_\tau, N_x$ and a learning rate $\eta$.

* **Main loop:** REPEAT until satisfactory convergence ($C \ll 1$)

(1) Draw $N_\tau \times N_x$ i.i.d. points in $D_{\tau,x}$, denoted by $\mathbf{m} = (\tau_{ij}, x_{ij})_{i,j}$;

(2) Compute $\nabla_{\theta}C[\theta](\mathbf{m})$ from (C1) by automatic differentiation;

(3) Update $\theta: \theta \leftarrow \theta - \eta \nabla_{\theta}C$. We employ the stochastic descent algorithm ADAM [54];

(4) Monitor every few steps $W$ and $C$.

In practice, we use a feedforward NN with 10 hidden layers, with 10 neurons/layer (the weights matrices $W$ are of size $10 \times 10$), and we use a swish activation function $f$ [57]. The training rate is $\eta = 10^{-3}$ (ADAM) decreased to $\eta = 10^{-4}$ at the end of the descent. The typical parameters in the simulations are $15 \le T \le 20$, $\sigma = 8$, $N_\tau = 100$, $N_x = 50$. These parameters are not fine-tuned and can be modified. More powerful setups are possible, but they are not discussed here.

## ORCID iDs

Ciro Campolina https://orcid.org/0009-0000-3282-3675
Eric Simonnet https://orcid.org/0000-0003-4636-4910

Simon Thalabard https://orcid.org/0000-0003-2371-7242

## References

[1] Gilson J, Daumont I and Dombre T 1998 A two-fluid picture of intermittency in shell-models of turbulence *Advances in Turbulence VII* (Springer) pp 219–22

[2] Mailybaev A A 2012 Computation of anomalous scaling exponents of turbulence from self-similar instanton dynamics *Phys. Rev. E* **86** 025301

[3] Mailybaev A A 2016 Spontaneous stochasticity of velocity in turbulence models *Multiscale Model. Simul.* **14** 96

[4] Drivas T and Mailybaev A A 2021 Life after death in ordinary differential equations with a non-Lipschitz singularity *Nonlinearity* **34** 2296

[5] Chae D 2008 Incompressible Euler equations: the blow-up problem and related results *Handbook of Differential Equations: Evolutionary Equations* vol 4 (North-Holland) p 1

[6] Eggers J and Fontelos M 2008 The role of self-similarity in singularities of partial differential equations *Nonlinearity* **22** R1

[7] Frisch U 1997 Is there finite-time blow-up in 3-D Euler flow? *Current Developments in Mathematics* vol 1997 (International Press of Boston) p 193 (available at: https://intlpress.com/site/pub/files/_fulltext/journals/cdm/1997/1997/0001/CDM-1997-1997-0001-a007.pdf)

[8] Gibbon J D 2008 The three-dimensional euler equations: where do we stand? *Physica D* **237** 1894

[9] Hou T and Li R 2008 Blowup or no blowup? The interplay between theory and numerics *Physica D* **237** 1937

[10] Bustamante M and Brachet M 2012 Interplay between the Beale-Kato-Majda theorem and the analyticity-strip method to investigate numerically the incompressible Euler singularity problem *Phys. Rev. E* **86** 066302

[11] Campolina C and Mailybaev A 2018 Chaotic blowup in the 3D incompressible Euler equations on a logarithmic lattice *Phys. Rev. Lett.* **121** 064501

[12] Krstulovic G and Nazarenko S 2024 Initial evolution of three-dimensional turbulence en route to the Kolmogorov state: emergence and transformations of coherent structures, self-similarity and instabilities *Phys. Rev. Fluids* **9** 034610

[13] Barenblatt G and Zel'Dovich Y 1972 Self-similar solutions as intermediate asymptotics *Annu. Rev. Fluid Mech.* **4** 285

[14] Wang Y, Lai C-Y, Gómez-Serrano J and Buckmaster T 2023 Asymptotic self-similar blow-up profile for three-dimensional axisymmetric Euler equations using neural networks *Phys. Rev. Lett.* **130** 244002

[15] Leith C 1967 Diffusion approximation to inertial energy transfer in isotropic turbulence *Phys. Fluids* **10** 1409

[16] Nazarenko S 2011 *Wave Turbulence* vol 825 (Springer)

[17] Biferale L 2003 Shell models of energy cascade in turbulence *Annu. Rev. Fluid Mech.* **35** 441

[18] Nakano T 1988 Determination of a dynamical scaling function in a cascade model of turbulence *Prog. Theor. Phys.* **79** 569

[19] Dombre T and Gilson J-L 1998 Intermittency, chaos and singular fluctuations in the mixed Obukhov–Novikov shell model of turbulence *Physica D* **111** 265

[20] Andersen K, Bohr T, Jensen M, Nielsen J and Olesen P 2000 Pulses in the zero-spacing limit of the GOY model *Physica D* **138** 44

[21] Katz N and Pavlović N 2005 Finite time blow-up for a dyadic model of the Euler equations *Trans. Am. Math. Soc.* **357** 695

[22] Barbato D, Flandoli F and Morandin F 2011 Energy dissipation and self-similar solutions for an unforced inviscid dyadic model *Trans. Am. Math. Soc.* **363** 1925

[23] Mailybaev A A 2012 Renormalization and universality of blowup in hydrodynamic flows *Phys. Rev. E* **85** 066317

[24] Thalabard S, Nazarenko S, Galtier S and Medvedev S 2015 Anomalous spectral laws in differential models of turbulence *J. Phys. A: Math. Theor.* **48** 285501

[25] Galtier S, Nazarenko S, Buchlin E and Thalabard S 2019 Nonlinear diffusion models for gravitational wave turbulence *Physica D* **390** 84
[26] Thalabard S, Medvedev S, Grebenev V and Nazarenko S 2021 Inverse cascade anomalies in fourth-order Leith models *J. Phys. A: Math. Theor.* **55** 015702

[27] Giga Y and Kohn R V 1987 Characterizing blowup using similarity variables *Indiana Univ. Math. J.* **36** 1

[28] L'vov V S 2002 Quasisolitons and asymptotic multiscaling in shell models of turbulence *Phys. Rev. E* **65** 026309

[29] L'vov V S, Podivilov E, Pomyalov A, Procaccia I and Vandembroucq D 1998 Improved shell model of turbulence *Phys. Rev. E* **58** 1811

[30] Gledzer E 1973 System of hydrodynamic type admitting two quadratic integrals of motion *Sov. Phys. Dokl.* **18** 216

[31] Yamada M and Ohkitani K 1987 Lyapunov spectrum of a chaotic model of three-dimensional turbulence *J. Phys. Soc. Japan* **56** 4210

[32] Ditlevsen P 2010 *Turbulence and Shell Models* (Cambridge University Press)

[33] Desnyansky V and Novikov E 1974 The evolution of turbulence spectra to the similarity regime *Izv. Akad. Nauk SSSR Fiz. Atmos. Okeana* **10** 127

[34] Ohkitani K and Yamada M 1989 Temporal intermittency in the energy cascade process and local Lyapunov analysis in fully-developed model turbulence *Prog. Theor. Phys.* **81** 329

[35] Plunian F, Stepanov R and Frick P 2013 Shell models of magnetohydrodynamic turbulence *Phys. Rep.* **523** 1

[36] Pisarenko D, Biferale L, Courvoisier D, Frisch U and Vergassola M 1993 Further results on multifractality in shell models *Phys. Fluids A* **5** 2533

[37] Biferale L, Pierotti D and Toschi F 1998 Helicity transfer in turbulent models *Phys. Rev. E* **57** R2515

[38] Jensen M, Paladin G and Vulpiani A 1991 Intermittency in a cascade model for three-dimensional turbulence *Phys. Rev. A* **43** 798

[39] Angheluta L, Benzi R, Biferale L, Procaccia I and Toschi F 2006 Anomalous scaling exponents in nonlinear models of turbulence *Phys. Rev. Lett.* **97** 160601

[40] de Wit X, Ortali G, Corbetta A, Mailybaev A, Biferale L and Toschi F 2024 Extreme statistics and extreme events in dynamical models of turbulence *Phys. Rev. E* **109** 055106

[41] Creswell J, Mukhanov V and Oz Y 2024 Anomalous scaling in shell model turbulence (arXiv:2402.07813)

[42] Constantin P, Levant B and Titi E 2007 Regularity of inviscid shell models of turbulence *Phys. Rev. E* **75** 016304

[43] Cheskidov A 2008 Blow-up in finite time for the dyadic model of the Navier-Stokes equations *Trans. Am. Math. Soc.* **360** 5101

[44] Mailybaev A A 2013 Blowup as a driving mechanism of turbulence in shell models *Phys. Rev. E* **87** 053011

[45] Mailybaev A A 2016 Spontaneously stochastic solutions in one-dimensional inviscid systems *Nonlinearity* **29** 2238

[46] Mailybaev A A 2013 Bifurcations of blowup in inviscid shell models of convective turbulence *Nonlinearity* **26** 1105

[47] Barenblatt G I 1996 *Scaling, Self-Similarity and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics* vol 14 (Cambridge University Press)

[48] Diekmann O, Van Gils S A, Lunel S M and Walther H-O 2012 *Delay Equations: Functional-, Complex- and Nonlinear Analysis* vol 110 (Springer)

[49] Veltz R 2020 BifurcationKit.jl (available at: https://hal.archives-ouvertes.fr/hal-02902346)

[50] Clewley R, Sherwood W, LaMar M and Guckenheimer J 2007 PyDSTool, a software environment for dynamical systems modeling (available at: http://pydstool.sourceforge.net)

[51] Grebenev V, Nazarenko S, Medvedev S, Schwab I and Chirkunov Y 2013 Self-similar solution in the leith model of turbulence: anomalous power law and asymptotic analysis *J. Phys. A: Math. Theor.* **47** 025501

[52] Lagaris I, Likas A and Fotiadis D 1998 Artificial neural networks for solving ordinary and partial differential equations *IEEE Trans. Neural Netw.* **9** 987

[53] Raissi M, Perdikaris P and Karniadakis G 2019 Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations *J. Comput. Phys.* **378** 686

[54] Kingma D P and Ba J 2015 Adam: a method for stochastic optimization *Proc. 3rd Int. Conf. on Learning Representations (ICLR)*
[55] Chae D 2007 Nonexistence of self-similar singularities for the 3D incompressible Euler equations
*Commun. Math. Phys.* **273** 203

[56] Elgindi T M 2021 Finite-time singularity formation for $C^{1,\alpha}$ solutions to the incompressible Euler equations on $\mathbb{R}^3$ *Ann. Math.* **194** 647

[57] Ramachandran P, Zoph B and Le Q 2017 Searching for activation functions (arXiv:1710.05941)