# Breakthrough-induced loop formation in evolving transport networks

Stanislaw Żukowski, Annemiek Johanna Maria Cornelissen, Florian Osselin, Stéphane Douady, Piotr Szymczak

► To cite this version:

Stanislaw Żukowski, Annemiek Johanna Maria Cornelissen, Florian Osselin, Stéphane Douady, Piotr Szymczak. Breakthrough-induced loop formation in evolving transport networks. Proceedings of the National Academy of Sciences of the United States of America, 2024, 121, 10.1073/pnas.2401200121. insu-04676495

HAL Id: insu-04676495

https://insu.hal.science/insu-04676495v1

Submitted on 15 Jul 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Breakthrough-Induced Loop Formation in Evolving Transport Networks

Stanisław Żukowski,¹,²,* Annemiek Johanna Maria Cornelissen,²
Florian Osselin,³ Stéphane Douady,² and Piotr Szymczak¹

¹Institute of Theoretical Physics, Faculty of Physics,
University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

²Laboratoire Matière et Systèmes Complexes (MSC), UMR 7057,
CNRS & Université Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France

³Institut des Sciences de la Terre d'Orléans, UMR 7327,
CNRS & BRGM & Université d'Orléans, rue de la Ferollerie, 45100 Orléans, France

(Dated: July 31, 2024)

Transport networks, such as vasculature or river networks, provide key functions in organisms and the environment. They usually contain loops whose significance for the stability and robustness of the network is well documented. However, the dynamics of their formation is usually not considered. Such structures often grow in response to the gradient of an external field. During evolution, extending branches compete for the available flux of the field, which leads to effective repulsion between them and screening of the shorter ones. Yet, in remarkably diverse processes, from unstable fluid flows to the canal system of jellyfish, loops suddenly form near the breakthrough when the longest branch reaches the boundary of the system. We provide a physical explanation for this universal behavior. Using a 1D model, we explain that the appearance of effective attractive forces results from the field drop inside the leading finger as it approaches the outlet. Furthermore, we numerically study the interactions between two fingers, including screening in the system and its disappearance near the breakthrough. Finally, we perform simulations of the temporal evolution of the fingers to show how revival and attraction to the longest finger leads to dynamic loop formation. We compare the simulations to the experiments and find that the dynamics of the shorter finger are well reproduced. Our results demonstrate that reconnection is a prevalent phenomenon in systems driven by diffusive fluxes, occurring both when the ratio of the mobility inside the growing structure to the mobility outside is low and near the breakthrough.

Keywords: nonlinear physics | unstable growth processes | transport networks

## I. SIGNIFICANCE STATEMENT

Loops are ubiquitous in animate and inanimate transport networks, from leaf venation to river deltas. Yet, current physical models fail to reproduce them accurately. We report a way of loop formation that manifests in remarkably diverse processes – from unstable fluid flows to the canal system of jellyfish. Near the breakthrough, when one of the branches reaches the outlet of the system, shorter branches revive and accelerate towards the longest one, ultimately leading to a reconnection event. This results in a hierarchical structure of interconnected fingers. Our explanation of the shared underlying mechanism behind the phenomenon significantly advances the understanding of how looping transport networks dynamically emerge.

## II. INTRODUCTION

Nature offers us a wide spectrum of spatial, transport networks, which provide key functions in living organisms and the surrounding environment. Examples span from leaf venation [2], blood vessels [3], jellyfish gastrovascular

canal system [1], to river networks [4, 5], deltas [6, 7], or cave conduits [8, 9]. Topologically, these networks can take the form of either branched, tree-like structures or looping patterns. The latter are more robust against damage [10], hence they are often favored by biological evolution. For example, it is thought that leaf venation is branched in ancestral plants, but highly reticulated in more evolutionary recent ones [11–13].

While the question of why looping networks might have been chosen by evolution and what quantity, if any, they optimize, was extensively studied [6, 10, 14–18], the question of how branches in the network interact to form a looping network remains unclear. A natural candidate for an underlying mechanism leading to reconnections is a tensorial stress field [19, 20]. However, many networks in nature develop as a result of unstable growth processes in a scalar field [21, 22].

The growth in such a case is driven by the gradient of an external field, such as electric potential or concentration. Growing parts of the structure compete for the available flux of the field, effectively interacting with each other. The branches in the network try to go away from each other, to maximize the flux coming to their tips and avoid being shadowed [21–23]. This mechanism, however, cannot explain the formation of loops, for which an effective attraction between the growing branches have to be present.

Surprisingly, loops suddenly appear across numerous

<footnote>* s.zukowski@uw.edu.pl</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Systems near breakthrough in nature. (a) Each system consists of two phases: an invading phase ($\mathcal{F}$) with mobility $\lambda_1$, and a displaced phase ($\mathcal{O}$) with mobility $\lambda_2$, separated by an interface. The growth direction is from left to right. Red arrows mark places where the fingers are near breakthrough or have already broken through, and blue arrows mark reconnections. (b-e) Examples of reconnection near breakthrough in various systems: (b) a fracture dissolution experiment in a Hele-Shaw cell, (c) viscous fingers in the Saffman-Taylor experiment, (d) streamer channels in air (photo: Sander Nijdam, Eindhoven University of Technology, by permission), (e) an octant of the jellyfish *Aurelia*, showing the gastrovascular canal network in dark gray [1].</figcaption></figure>

systems as the longest branch reaches the outlet of the system. Here we report that near the breakthrough, the shorter branch grows toward the longest one and reconnects to it (Fig. 1a). This results in a characteristic hierarchical pattern found in a large variety of systems (Fig. 1b-e): dissolving fractures [24, 25], viscous fingering [26], discharge patterns [27], and even the growth of a gastrovascular canal network of the jellyfish *Aurelia* [1]; see *Materials and Methods*, section VA for a detailed description of the examples and experiments. The ubiquity of this process suggests the existence of a shared underlying mechanism, which we elucidate here.

### A. Laplacian growth

In many unstable growth processes, the boundary – the interface between two phases – moves due to the external forcing, such as pressure gradient between the inlet and outlet of the system (Fig. 1a). Important parameters in such models are the mobilities of the two phases, e.g. hydraulic permeabilities for pressure driven growth. Whenever the mobility of the invading phase ($\lambda_1$) is larger than the mobility of the displaced phase ($\lambda_2$) the flux concentrates on small protrusions of the interface and the front can break into fingers [22, 28, 29]. Because of the flux concentration the fingers tend to grow more in length than in diameter. The width of the fingers is then set by

surface tension, or other short-scale regularization mechanisms. Additional effects, such as tip splitting, can give rise to a highly ramified, hierarchical tree-like structure.

A paradigm for such growth processes is Laplacian growth, where the fingers extend with velocity proportional to the flux of a diffusive field ($\phi_i$) given by: $\vec{J}_i = -\lambda_i \nabla \phi_i$, for $i = 1, 2$ depending on the phase. The conservation of the flux results in the Laplace equation for the field in both the invading ($\mathcal{F}$) and the displaced phase ($\mathcal{O}$):

$$ \Delta\phi_1(\vec{x}) = 0 \quad \vec{x} \in \mathcal{F} \quad \text{and} \quad \Delta\phi_2(\vec{x}) = 0 \quad \vec{x} \in \mathcal{O}. \quad (1) $$

These equations are supplemented with the continuity condition for the field and its flux at the interface ($\mathcal{T}$):

$$ \phi_1(\vec{x}) = \phi_2(\vec{x}) \quad \vec{x} \in \mathcal{T}, \\ \lambda_1(\nabla\phi_1(\vec{x}))_n = \lambda_2(\nabla\phi_2(\vec{x}))_n \quad \vec{x} \in \mathcal{T}, \quad (2) $$

where $n$ denotes the normal to the interface. Additionally, the Dirichlet boundary condition is imposed on the inlet ($x = 0$) and outlet ($x = 1$) of the system:

$$ \phi_1(x = 0) = 1 \quad \text{and} \quad \phi_2(x = 1) = 0. \quad (3) $$

Note that the coordinates here are rescaled by the system length, and the field is rescaled by the value at the inlet.
## B. Effect of mobility ratio

However, reticulated networks are usually not obtained in such models [21–23]. This is due to the simplifying assumption that neglects the mobility of the invading phase, which is equivalent to taking the limit of mobility ratio going to infinity, $M = \lambda_1/\lambda_2 \to \infty$. This leads to the omission of the field drop within the fingers. The Laplace equation is then solved only in the displaced domain with a constant value of the field directly on the moving boundary. The Dirichlet boundary condition on the fingers results in: (i) long fingers screening the shorter ones and hindering their growth; (ii) two parallel fingers growing away from each other as they get more flux from the sides.

If one takes into account a finite mobility ratio the effective repulsion and screening diminishes. The field inside the invading phase is then no longer constant and the resulting field gradients can make the fingers attract each other and create loops [30, 31]. As shown by Budek et al. [31], this can occur only for a specific range of mobility and finger length ratios.

However, as presented in Fig. 1, the breakthrough reconnections are not limited to these specific scenarios. Notably, they can occur even in the general case of high mobility ratio, which was previously thought to be impossible due to screening effects. Nonetheless, near the breakthrough, screening diminishes, enabling the revival and growth of shorter branches towards the longest one, ultimately leading to a reconnection event. We demonstrate that contrary to previous studies reconnection is a prevalent phenomenon in Laplacian growth, occurring both in low mobility ratio cases and near the breakthrough.

# III. RESULTS AND DISCUSSION

## A. Drop of potential along a single finger

To understand the generality of the breakthrough reconnections we begin with a 1D case. Here, $\mathcal{F} = \{x \in [0, x_t]\}$, $\mathcal{O} = \{x \in [x_t, 1]\}$ and $\mathcal{T} = \{x_t\}$ is just the fingertip. The solution of the Eqs. (1)-(3) are then two piecewise linear functions stitched at the fingertip (Fig. 2a)

$$ \phi(x) = \begin{cases} 1 - \frac{x}{M - x_t(M-1)}, & \text{if } x \in [0, x_t] \\ \frac{1-x}{1-x_t(1-1/M)}, & \text{if } x \in [x_t, 1]. \end{cases} \quad (4) $$

A few remarks can be made based on this simple example.

First, for a finite mobility ratio ($M = 100$, Fig. 2a) when the fingertip is far from the outlet ($x_t = l_s \ll 1$), approximation of a constant field on the moving boundary works well, and the field inside the short finger can be treated as constant. However, if we take a longer finger ($x_t = l_1 \approx 1$), there is a significant potential drop inside it. By comparing the field values at $x = l_s$ we see that

<figure><img src="image_3.png" /><figcaption>FIG. 2. Solutions of Eqs. (1)-(3) for a 1D system. (a) Potential in the system with mobility ratio $M = 100$ for a short finger of length $l_s$ and a longer finger of length $l_1 > l_s$. The dark line marks the field inside the finger ($\phi_i$), the lighter line marks the field outside ($\phi_o$), and the black dot marks the position of the fingertip. (b) A profile of the potential at the fingertip ($\phi_t$) as a function of the tip position ($x_t$) for different mobility ratios. Such a profile for $M = 100$ is also plotted in (a) with a thin gray line.</figcaption></figure>

there is a difference of potential between the two fingers:

$$ \delta\phi \approx 1 - \phi(x = l_s, x_t = l_1) = \frac{l_s}{1 + d(M - 1)}, \quad (5) $$

where $d = 1 - l_1$. If such fingers were placed at a distance $\delta y$ next to each other, there would be a gradient of potential $\delta\phi/\delta y$, between the shorter finger tip and the longer finger, provided that the fingers do not influence each other. As a result, the shorter finger would be attracted toward the longer one.

Second, let us focus on the potential at the fingertip as a function of the tip position $\phi_t(x_t) = \phi(x = x_t)$ (Fig. 2a, light gray line and Fig. 2b), and analyze how it changes with the mobility ratio. The higher the $M$, the steeper the profile of the potential at the fingertip, and the longest finger must be closer to the outlet to feel its impact. As can be seen in (5), $\delta\phi$ is inversely proportional to the product $Md$. This suggest that the critical distance, $d_c$, at which the pressure inside the longer finger begins to drop and the attraction between two fingers would appear, scales as $d_c \sim M^{-1}$.

In the limiting case $M \to \infty$ the functional dependence $\phi_t(x_t)$ becomes a step function (Fig. 2b). The fingers have a constant potential along their length, no matter how close they are to the outlet, and there is no difference in potential between them. However, when one of the fingers breaks through, the potential inside it takes the form $\phi(x) = 1 - x$. This instantaneous pressure drop inside the longer finger induces a pressure difference with respect to the shorter finger, $\delta\phi = l_s$. Consequently, at the moment of breakthrough, we observe a sudden transition from no interaction between the fingers to attraction of the shorter finger to the longer one.
## B. Finger interactions in two dimensions—screening and revival

Although insightful, the 1D model treats two fingers independently of each other, and does not take into account effects of finger interactions such as screening. To investigate these effects we conduct numerical simulations in a 2D geometry. Two fingers are placed in a cell of length $L = 1$ and width $W = L/3$ with periodic boundary conditions on the bottom ($y = 0$) and top ($y = W$) wall. The long finger of length $l_1$ is positioned at $y = 0$, and the short finger of length $l_s$ is at $y = W/3$, as shown in Figs. 3a-b. The fingers have a shape of thin rectangles of width $W/15$ with semicircular caps. We solve the equations for the field with the finite element method implemented in the FreeFEM++ software [32] (Materials and Methods, section V B).

Figs. 3a-b represent isolines of the field (the same set of values) in two cases: $M = 10^6$ and $M = 10^2$. We observe that for high mobility ratio there is a negligible potential drop inside the fingers, and the longer finger attracts almost all of the flux. As a result, it screens the shorter one and would suppress its growth. For lower mobility ratio, due to the potential drop inside the longer finger, some flux can reach the shorter one, so it can still grow.

We consider the field values along the center line of the long and short fingers for the two mobility ratios (Figs. 3c-d). The profiles along the longer finger are piecewise linear, similar to the 1D solutions presented in Fig. 2a. However, the field values outside the shorter fingers ($\phi_o^s$) are strongly influenced by the longer ones. In fact, the difference between the field inside the longer finger and the field outside the shorter one, $\phi_i^1 - \phi_o^s$, depends exponentially from the perspective of the longer finger tip, as shown by the black dashed lines in the insets in Figs. 3c-d. As a result, the shorter finger is strongly screened and almost no flux reaches its tip.

The situation changes drastically as the longer finger approaches the outlet. In Fig. 4a we analyze the total flux $Q$ through the tip of the shorter finger as a function of gap $d$ – the distance from the longer finger tip to the outlet of the system. The length ratio is kept constant $l_s/l_1 = 1/3$ and $Q$ is calculated by integrating the flux at the shorter finger tip for $\theta \in [-\pi/2, \pi/2]$, see inset in Fig. 3a.

For relatively low mobility ratio, the shorter finger is weakly screened and receives some flux, even when the longer finger is far from the outlet (large $d$). The higher the $M$, the stronger the screening effect, hence for large $d$ the total flux in the shorter finger is close to zero, and it only starts to increase as the gap gets smaller. Note that the critical distance from the outlet, at which the screening starts to diminish, decreases as $M$ gets larger. This is a further manifestation of the fact that the critical gap at which screening vanishes scales as $M^{-1}$, as already described with the 1D model. The above confirms that for sufficiently small $d$, or after the breakthrough, short

<figure><img src="image_4.png" /><figcaption>FIG. 3. Numerical solutions of Eqs. (1)-(3) for a 2D system. (a)-(b) A schematic of the 2D setup with two fingers (of lengths $l_s = 0.3$ and $l_1 = 0.9$) in a periodic cell of length $L = 1$ and width $W = L/3$ with the isolines of the field for: (a) $M = 10^2$ and (b) $M = 10^6$. In the inset we mark the attractive ($Q_a$) and repulsive ($Q_r$) fluxes accepted by the right and the left side of the fingertip, respectively. (c)-(d) Field values along the longer and shorter finger (cross sections marked with dashed lines in (a) and (b)) for the two mobility ratios. Insets in (c)-(d): the difference between the field inside the longer finger and outside the shorter one: $\phi_i^1 - \phi_o^s$. The black dashed line is a fit to the linear part of the plot in inset (c), highlighting the exponential decay of $\phi_i^1 - \phi_o^s$ (from the longer finger tip perspective).</figcaption></figure>

fingers that were strongly screened would revive and start growing again, as also noted in Ref. [25, 26].

Interestingly, for large $M$ the $Q(d)$ dependence approaches a universal function that simply shifts in $d$ as $M$ increases. This suggests that the flow profiles should collapse onto a master curve when plotted as a function of $Md$, which is indeed the case (Fig. 4b). Moreover, if we define rescaled $d$: $d' = Md$, plug it into (5), and take the limit of $M \to \infty$ we see that the master curve will be of the form $1/(1 + d')$. This can be adapted for a 2D system and written as

$$Q(d') = q_0 + (q_1 - q_0)/(1 + d'/d'_c). \quad (6)$$

The three parameters $q_0, q_1, d'_c$ can be easily extracted from the simulation results and have a clear physical interpretation: $q_0$ and $q_1$ are the flux reaching the shorter
<figure><img src="image_5.png" /><figcaption>FIG. 4. Revival and interactions between the fingers when the gap changes. (a) Total flux intercepted by the shorter finger tip, $Q = Q_a + Q_r$ (as defined in inset in Fig. 3a), as a function of $d$. The length ratio is kept constant, $l_s/l_1 = 1/3$. (b) Total flux plotted as a function of $M_d$. The master curve, red dashed line, corresponds to (6) plotted with $q_0 = 0.05$, $q_1 = 9.76$, $d_c = 0.23$. (c) Maps of interactions between the fingers, $Q_a - Q_r$, as a function of two parameters: mobility ratio ($M$) and finger length ratio ($l_s/l_1$). Maps are plotted for four gaps: $d = 1/2, 10^{-2}, 10^{-4}, 0$. Red colors indicate attraction, gray repulsion, black dashed line separates the two.</figcaption></figure>

finger tip during screening phase and after revival, respectively; $d'_c$, the inflection point of the sigmoid, is the rescaled critical gap at which the shorter finger in the two-dimensional system revives. As can be seen in Fig. 4b, dashed line, (6) perfectly captures revival in our system.

### C. Finger interactions in two dimensions—repulsion and attraction

Having described the revival of the shorter finger, let us discuss the change in their growth direction near the breakthrough. To quantify this, we calculate the attractive, $Q_a$, and repulsive flux, $Q_r$, by integrating the flux on the shorter finger tip for $\theta \in [0, \pi/2]$ and $\theta \in [-\pi/2, 0]$, respectively (see inset in Fig. 3a). We then take the difference of the two fluxes, $Q_a - Q_r$. The positive value of this quantity will be interpreted as an effective attraction, and should result in the shorter finger growing toward the longer one. Conversely, if the value is negative the fingers repel. Fig. 4c represents repulsion-attraction maps as a function of two parameters: mobility ratio ($M$) and finger length ratio ($l_s/l_1$). The maps are shown for four gaps: $d = 1/2, 10^{-2}, 10^{-4}, 0$; the latter corresponds to the breakthrough.

For large gap (Fig. 4c, $d = 1/2$), we observe an island of attraction in the region of lower mobility ratio, $M \in (10^0 - 10^2)$, and the length ratio, $l_s/l_1 \in (0, 0.6)$. This quantitatively agrees with the results of the resistor model presented by Budek et al. [31], where the interactions between the fingers distant from the outlet were analyzed. As described there and as can be seen in Fig. 4c, $d = 1/2$, for the systems of low mobility ratio the screening is weak and the fingers can attract each other to form loops. Note that in the other maps in Fig. 4c, regardless of the gap size $d$, the value of $Q_a - Q_r$ remains almost the same in the region of low mobility and length ratio. This suggests that in such systems the breakthrough will not drastically affect the dynamics and interactions between

the fingers.

The impact of the gap size on the interactions is more apparent in the region of higher mobility ratio, $M \in (10^3 - 10^6)$, and length ratio, $l_s/l_1 \in (0.6, 0.9)$. Here, the interactions change from slight to strong repulsion when the leading finger is distant from the outlet, $d = 1/2, 10^{-2}$, and finally, after the breakthrough, $d = 0$, we observe a transition from strong repulsion to strong attraction. Hence, for systems of higher mobility ratio the impact of the breakthrough should be more striking in the dynamics of the fingers. In particular, even for very high mobility ratio, the longer finger will always attract the shorter ones after the breakthrough.

### D. Temporal evolution of the fingers

Finally, we perform dynamic simulations of finger growth. During growth, we do not change the finger shape (constant width and semicircular tips) and extend it only in the direction from which the highest flux is coming (*Materials and Methods*, section V B). Here we analyze the evolution of two fingers with slightly different initial lengths, $l_s = 0.04$ vs $l_1 = 0.05$, and mobility ratio in the system $M = 1000$. With such initial conditions we can observe all previously described interactions between the fingers: competition and repulsion, screening, revival and attraction. The dynamics of the fingers are presented in a supplementary video *Movie S1*.

In Fig. 5 we show how the velocity of the shorter finger changes over time, along with the snapshots from the simulation. First, as the fingers start with similar lengths, we observe competition and repulsion between them – the fingers accelerate and grow with similar rate, until the longer finger wins and starts to screen the shorter one. At this moment ($t = t_1$), the shorter finger starts to slow down. As the longer finger approaches the outlet, the field inside it begins to decrease (as visible in *Movie S1*). Consequently, screening disappears and the overall flux within the shorter finger increases. The shorter
<figure><img src="image_6.png" /></figure>

FIG. 5. Velocity of the shorter finger in the dynamical simulation of growth with mobility ratio $M = 1000$. Black dashed line marks the total velocity of the finger, gray its $x$-component (pointing towards the outlet), red the $y$-component (repulsion from the longer finger) and green the absolute value of $v_y$ (attraction to the longer finger). Gray vertical lines mark: initial time ($t_0$); start of screening ($t_1$), when the acceleration of the finger becomes negative; start of revival ($t_2$), when the acceleration changes its sign again; breakthrough (BT, $t_3$) and reconnection ($t_4$). Snapshots from the simulation were taken at the corresponding moments of the evolution (lighter color on the fingers marks the parts that grew between successive moments).

finger revives and accelerates again ($t = t_2$).

Shortly after $t_2$ the attraction towards the longer finger appears: the positive $v_y$, away from the longer finger (depicted by the red line in Fig. 5) transitions to negative, toward the longer finger ($|v_y|$, marked with the green line in Fig. 5). This stage corresponds to the expansion of the attraction region on the interaction maps prior to the breakthrough (Figs. 4a). Then, just before the breakthrough ($t = t_3$), there is a sudden jump in velocity. This last boost significantly expedites the loop formation process. As the distance between the two fingers diminishes, the attraction intensifies. This cumulative, snowball-like effect, triggered by the breakthrough, culminates in the eventual reconnection at $t = t_4$.

## E. Comparison of the simulations to experiments

We additionally present two simulations that qualitatively reproduce the behavior of the fingers in the fracture dissolution and Saffman-Taylor experiments seen in Figs. 1b-c, with approximate mobility ratios of $M = 14$ and $M = 100$, respectively (see Materials and Methods, section VA 1 and VA 2 for details on how the mobility ratio in the experiments is calculated). In these simulations, the fingers were initiated with lengths of $l_s = 0.002$ and $l_1 = 0.18$. In Supplemental Material we include videos from the simulations (Movie S2 and Movie S3) and videos from the experiments (Movie S4 and Movie S5). In Fig. 6 we show the plots of the velocity of the shorter finger, both in the simulations and in the experiments, in the case of fracture dissolution ($M = 14$) and viscous fingering ($M = 100$).

As predicted in the previous section, in the case of the lower mobility ratio, the shorter finger is attracted

towards the longer one almost from the very beginning of the evolution, both in the simulation (Fig. 6a, Movie S2) and the experiment (Fig. 6b, Movie S4). There is also no effect of screening, and flux is non-negligible even when the longer finger is far away from the outlet. As a result, the finger grows with an almost constant total velocity, which remains unchanged even during the breakthrough.

In contrast, for higher mobility ratio, the shorter finger in the simulation (Fig. 6c, Movie S3) and in the experiment (Fig. 6d, Movie S5) is initially screened by the longer one and grows relatively slowly. Only after the breakthrough does it revive and eventually reconnect forming a loop.

## IV. CONCLUSIONS

As previously described in the literature, attractive interactions between the fingers can appear in the systems of low mobility ratio, leading to reconnections. For high mobility ratio, only screening and repulsion have been observed, resulting in branched loopless structures. However, when studying systems such as viscous fingering or fracture dissolution, experiments are often terminated when the invading phase reaches the border of the system. An unexpected behavior arises in a broad class of unstable growth processes when the proximity of the outlet is considered and the evolving structure breaks through.

We have shown that a striking transition in the system dynamics occurs when the breakthrough is reached, especially in the case of infinite mobility ratio, where a singular limit emerges. Prior to the breakthrough, the field along the longest finger remains constant regardless of its proximity to the outlet, and the finger screens the rest of the system. At the moment of breakthrough, however, the field in the longer finger drops dramatically, allowing the flux to reach the shorter neighboring fingers. This affects the dynamics of the shorter fingers, causing their revival and a strong attraction to the longer finger. The interplay of these two factors leads to reconnection and loop formation.

To explain breakthrough reconnections, we used a simplified model of growing Laplacian fingers: (i) we neglected the jump in the field across the interface associated with the regularization mechanism; (ii) we assumed that the change of mobility across the interface is discontinuous, rather than smooth as in, for instance, fracture dissolution; (iii) we kept the shape of the fingers unchanged during their evolution, whereas in some systems, such as viscous fingering, when the fingers receive more flux, their tips begin to grow in width, which impacts their velocity [33, 34]. Regardless, for high mobility ratio, boundary condition on the fingers transitions from a constant field in the initial state to a linear gradient near the breakthrough. This field change near the breakthrough is much larger than any jump in the field resulting from regularization mechanisms (Materials and
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of simulations to experiments. (a)-(b) Velocities of the shorter finger in the fracture dissolution simulation and experiment, respectively. The mobilities here are related to the aperture values in the dissolved and undissolved regions (*Materials and Methods*, section V A 1). The mobility ratio is $M = 14$. (c-d) Velocity of the shorter finger in the viscous fingering simulation and experiment. Here, the mobility ratio is the inverse of the viscosity ratio and is $M = 100$. The line colors and time points in the plots indicate the same as in Fig. 5. Snapshots were taken at the corresponding moments of the evolution. The lighter color on the fingers from the simulations marks the parts that grew between successive moments. The black arrow in (d) marks the shorter finger.</figcaption></figure>

*Methods*, section V A 2). The transition between the two boundary conditions might be more gradual due to (ii) and (iii), but it still significantly impacts the finger dynamics. This is evidenced by sudden loop formation observed in physical systems such as viscous fingering or discharge patterns, and also in living organisms such as the jellyfish *Aurelia* (Figs. 1b-e).

The breakthrough reconnections are expected to occur in virtually any system driven by diffusive fluxes. Observing it in a system of a yet unknown growth mechanism, such as the gastrovascular canal network of the jellyfish (Fig. 1e), is a strong indication that the system dynamics are controlled by the effective diffusion of a morphogenic parameter. This sheds new light on the possibility of dynamical loop formation in many systems.

# V. MATERIALS AND METHODS

## A. Experiments

### 1. Fracture dissolution

The fracture dissolution experiment presented in Fig. 1b was performed in a microfluidic setup described in detail in Ref. [35]. It consists of two polycarbonate disks. The bottom one contains a rectangular indentation (3.3 cm × 3.8 cm × 100 µm), which is initially filled with a soluble material, in our case plaster (Plaster of Paris,

Blik Modelarski Alabastrowy). The top plate contains a hierarchical system of inlet and outlet channels connected to large inlet/outlet reservoirs (4.5 cm × 5 mm × 2 mm). Such a design helps to maintain uniform pressure across the width of the plaster. The aperture above the plaster is created by gluing the plates together with an ultrathin PET-based double-coated tape 70 µm thick with a rectangular hole the size of the plaster block. The cast was prepared with a 60% (w/w) ratio of water to plaster. This yields an average porosity of the block of $\phi = 50\%$ (measured porosity to water) and a permeability of 45 mD (measured by injection of isopropanol). Pure water is injected into the system with a syringe pump (Harvard Apparatus PHD2000) at a rate of $q = 0.5$ ml/h. We recorded the experiment with a UI 1550LE-C-HQ CCD camera (IOS, Germany), acquiring photographic images of the system every 100 s. In order to ensure homogeneous light intensity over the system, we used a circular fluorescent illuminator.

For experiments conducted in a Hele-Shaw cell, the mobility can be expressed as $\lambda = h^3/12\mu$, where $h$ is the aperture of the Hele-Shaw cell available to the fluid and $\mu$ is the fluid viscosity. Thus, increasing the aperture in the dissolved part of the system (black area in Fig. 1b) effectively introduces two phases with different mobilities. The mobility ratio in this experiment was approximately $M = (h_1/h_2)^3 = (170/70)^3 \approx 14$.
## 2. Saffman-Taylor experiment

The Saffman-Taylor experiment presented in Fig. 1c was performed in a circular Hele-Shaw cell with a 1 mm separation between the bottom and top plates. The Hele-Shaw cell was initially filled with oil (dyed with paprika to increase contrast in the images). Water was injected into the system through an inlet located in the center of the top plate. The experiment was recorded with a Nikon D3000 camera. A LED panel was used below the Hele-Shaw cell to ensure uniform light intensity throughout the system. the mobility ratio is the inverse of the viscosity ratio and was of the order of $M = \mu_2/\mu_1 \approx 100$.

The pressure jump in a Hele-Shaw is estimated using the Young-Laplace equation: $\Delta p^{\text{cap}} = \gamma(2/b + 1/r)$, where $\gamma = 49 \text{ mN/m}$ is the surface tension between oil and water, $b = 1$ mm is the spacing between the plates, $r \sim 1$ cm is a typical radius of curvature for the fingers. This yields $\Delta p^{\text{cap}} \sim 10^{-3}$ Pa. When the breakthrough occurs, the pressure inside the longest finger drops from the inlet value to a linear gradient between the inlet and outlet. The typical pressure difference between inlet and outlet in our Saffman-Taylor experiment was of the order of $10^2$ Pa. Thus, the change in $\phi$ induced by the breakthrough is 5 orders of magnitude larger than the pressure jump due to surface tension, and it can be neglected.

## 3. Jellyfish

The breakthrough reconnections can also be encountered in biological systems. The gastrovascular system of the jellyfish *Aurelia* is composed of canals in which seawater flows, carrying nutrients and oxygen to the surrounding tissues. New canals (sprouts) appear on the circulatory canal at the rim of the jellyfish sub-umbrella and grow toward four stomachs in the center of the jellyfish. We observed that the smaller sprouts reconnect to the long canal that just connected to the gastrocircular groove around the stomach [1] (Fig. 1e).

Jellyfish were reared at room temperature (22°C) in artificial seawater prepared by diluting 28 g of synthetic sea salt (Instant Ocean; Spectrum Brands, Madison, WI) per liter of osmosis water (osmolarity 1100 mOsm). In the laboratory, polyps of the Roscoff strain [36] are used (obtained by courtesy of Konstantin Khalturin from the Marine Genomics Unit, Okinawa Institute of Science and Technology Graduate University, Onna, Okinawa, Japan). A more detailed description of jellyfish rearing can be found in Ref. [1].

To observe the jellyfish gastrovascular system, the jellyfish are caught from the aquarium approximately 3 hours after feeding them with artemia. In this way, the gastrovascular canals are colored orange by the digested artemia. They are anesthetized with magnesium chloride dissolved in seawater. Then, they are placed in a Petri dish in shallow seawater with the subumbrella facing up. The images are taken by transillumination using

a Leica macro zoom (MACROFLUO LEICA Z16 APO S/No: 5763648) and a Photron Fastcam SA3 camera. The images were stitched using Adobe Photoshop. The canals in the images were highlighted with gray color.

## 4. Streamer discharge

Streamer discharges are fundamental building blocks of sparks and lightnings. An experimental photo presented in Fig. 1d was obtained by courtesy of Sander Nijdam [37]. Here the mobilities are the conductivities inside the discharge channel and the air. Breakthrough reconnections in this context have been observed in many experiments [38–40]. In Ref. [27] stereo photography was used to fully reconstruct the three-dimensional structure of the reconnection events.

## 5. Velocity of the fingers

To estimate the velocity of the fingers in the experiments we binarized the images and prepared a mask that covered the area of the image in which the finger was growing. Next, we calculated the area of the shorter finger (within the mask) in successive frames. The approximate velocity of the finger was defined as the rate of change of this area. Finally, the velocity was smoothed using the Savitzky-Golay filter (function savgol_filter from the Scipy package).

## B. Numerical simulations

The numerical calculations and simulations were performed using the codes available on the GitHub repository (https://github.com/stzukowski/reticuler) [41–44]. In the temporal simulations, the fingers have a constant width and end with semi-circular tips, as in the static numerical calculations. In each time step of the evolution we solve the equations for the field (Eqs. (1)-(3)) with the finite element method implemented in the FreeFEM++ software [32]. To calculate the field in the system, we decompose it into two domains: the fingers and the outside. The mobility field is declared, $M(x, y)$, which is equal to $\lambda_1$ inside the fingers and $\lambda_2$ outside. Our script solves the equation for the field in the form: $\nabla \cdot (M(x, y)\nabla\phi(x, y)) = 0$.

To determine the velocity of the finger, we integrate the flux entering its tip. Following the principle of local symmetry [41, 45], we extend the finger in the direction from which the highest flux comes (while keeping its width constant). The dynamical simulation in Fig. 5 was initiated with branches of length $l_s = 0.04$ and $l_l = 0.05$, while the simulations in Fig. 6 started with branches of lengths: $l_s = 0.002$ and $l_l = 0.18$. The velocity in y direction in Figs. 5 and 6 was treated with the Savitzky-Golay filter (savgol_filter function from the Scipy package) to
smooth the numerical noise caused by the small order of magnitude of this velocity.

ACKNOWLEDGMENTS

This project was partially supported by the Na-
tional Science Center (Poland) under research grant
2020/02/Y/ST3/00121 and by the Mission for Transver-

sal and Interdisciplinary Initiatives (MITI) of the French
National Center of Scientific Research (CNRS) in the
program Auto-Organisation 2021 and 2022 under re-
search grants 212680, 284366, 238661. We also acknowl-
edge the funding from the Polish National Agency for
Academic Exchange and Campus France under grant
BPN/BFR/2022/1/00033 in the program PHC Polo-
nium. We also thank Sélène Jeammet for her help with
the Saffman-Taylor experiment and Solène Song for the
jellyfish timelapse images.

[1] S. Song, S. Żukowski, C. Gambini, P. Dantan, B. Mau-
roy, S. Douady, and A. J. M. Cornelissen, Morphogenesis
of the gastrovascular canal network in Aurelia jellyfish:
Variability and possible mechanisms, Frontiers in Physics
**10**, 10.3389/fphy.2022.966327 (2023).

[2] G. J. Mitchison, A model for vein formation in higher plants, Proceedings of the Royal Society B **207**, 79 (1980).

[3] T.-H. Nguyen, A. Eichmann, F. Le Noble, and V. Fleury, Dynamics of vascular branching morphogenesis: The effect of blood and tissue flow, Physical Review E **73**, 061907 (2006).

[4] A. P. Petroff, O. Devauchelle, D. M. Abrams, A. E. Lobkovsky, A. Kudrolli, and D. H. Rothman, Geometry of valley growth, Journal of Fluid Mechanics **673**, 245 (2011).

[5] O. Devauchelle, A. P. Petroff, H. J. Seybold, and D. H. Rothman, Ramification of stream networks, Proceedings of the National Academy of Sciences **109**, 20832 (2012).

[6] A. Konkol, J. Schwenk, E. Katifori, and J. B. Shaw, Interplay of River and Tidal Forcings Promotes Loops in Coastal Channel Networks, Geophysical Research Letters **49**, e2022GL098284 (2022).

[7] W.-T. Ke, J. B. Shaw, R. C. Mahon, and C. A. Cathcart, Distributary Channel Networks as Moving Boundaries: Causes and Morphodynamic Effects, Journal of Geophysical Research **124**, 1878 (2019).

[8] P. Szymczak and A. J. C. Ladd, The initial stages of cave formation: Beyond the one-dimensional paradigm, Earth and Planetary Science Letters **301**, 424 (2011).

[9] I. Ghani, D. Koehn, R. Toussaint, and C. W. Passchier, Dynamic development of hydrofracture, Pure and Applied Geophysics **170**, 1685 (2013).

[10] E. Katifori, G. J. Szöllősi, and M. O. Magnasco, Dam-
age and fluctuations induce loops in optimal transport
networks, Physical Review Letters **104**, 048704 (2010).

[11] C. K. Boyce and A. H. Knoll, Evolution of developmental potential and the multiple independent origins of leaves in Paleozoic vascular plants, Paleobiology, 70 (2002).

[12] C. K. Boyce, The evolutionary history of roots and leaves, in *Vascular Transport in Plants*, Physiological Ecology, edited by N. M. Holbrook and M. A. Zwieniecki (Academic Press, Burlington, 2005) pp. 479–499.

[13] S. Douady, C. Lagesse, M. Atashinbar, P. Bonnin, R. Pousse, and P. Valcke, A work on reticulated patterns, Comptes Rendus. Mécanique **348**, 659 (2020).

[14] J. R. Banavar, F. Colaiori, A. Flammini, A. Maritan, and A. Rinaldo, Topology of the Fittest Transportation Network, Physical Review Letters **84**, 4745 (2000).

[15] S. Bohn and M. O. Magnasco, Structure, Scaling, and Phase Transition in the Optimal Transport Network, Physical Review Letters **98**, 088702 (2007).

[16] M. Durand, Structure of Optimal Transport Networks Subject to a Global Constraint, Physical Review Letters **98**, 088701 (2007).

[17] H. Ronellenfitsch and E. Katifori, Global optimization, local adaptation, and the role of growth in distribution networks, Physical Review Letters **117**, 138301 (2016).

[18] F. Kaiser, H. Ronellenfitsch, and D. Witthaut, Discon-
tinuous transition to loop formation in optimal supply
networks, Nature Communications **11**, 5796 (2020).

[19] Y. Couder, L. Pauchard, C. Allain, M. Adda-Bedia, and S. Douady, The leaf venation as formed in a tensorial field, The European Physical Journal B **28**, 135 (2002).

[20] S. Bohn, L. Pauchard, and Y. Couder, Hierarchical crack pattern as formed by successive domain divisions., Physical Review E **71**, 046214 (2005).

[21] V. Fleury, J.-F. Gouyet, and M. Leonetti, *Branching in Nature: Dynamics and Morphogenesis of Branching Structures, from Cell to River Networks*, Vol. 14 (Springer Science & Business Media, 2001).

[22] P. Ball, *Branches: Nature's Patterns: A Tapestry in Three Parts* (OUP Oxford, 2009).

[23] P. Meakin, *Fractals, Scaling and Growth Far from Equilibrium*, Vol. 5 (Cambridge University Press, 1998).

[24] R. L. Detwiler, R. J. Glass, and W. L. Bourcier, Experimental observations of fracture dissolution: The role of Peclet number on evolving aperture variability, Geophysical Research Letters **30** (2003).

[25] V. Starchenko, C. J. Marra, and A. J. C. Ladd, Three-dimensional simulations of fracture dissolution, Journal of Geophysical Research **121**, 6421 (2016).

[26] S. Yang, G. Kong, Z. Cao, Z. Wu, and H. Li, Hydrody-
namics of gas-liquid displacement in porous media: Fin-
gering pattern evolution at the breakthrough moment
and the stable state, Advances in Water Resources ,
104324 (2022).

[27] S. Nijdam, C. G. C. Geurts, E. M. van Veldhuizen, and U. Ebert, Reconnection and merging of positive streamers in air, Journal of Physics D **42**, 045201 (2009).

[28] M. Muskat, The flow of fluids through porous media, Journal of Applied Physics **8**, 274 (1937).

[29] P. S. Stevens, *Patterns in Nature*, 1st ed. (Little, Brown and Company, 1974).

[30] A. Luque and U. Ebert, Growing discharge trees with self-consistent charge transport: The collective dynamics of streamers, New Journal of Physics **16**, 013039 (2014).

[31] A. Budek, K. Kwiatkowski, and P. Szymczak, Effect of
mobility ratio on interaction between the fingers in unstable growth processes, Physical Review E 96, 042218 (2017).

[32] F. Hecht, New development in FreeFem++, Journal of Numerical Mathematics 20, 251 (2012).

[33] I. Bischofberger, R. Ramachandran, and S. R. Nagel, An island of stability in a sea of fingers: Emergent global features of the viscous-flow instability, Soft Matter 11, 7428 (2015).

[34] O. Zik and E. Moses, Fingering instability in combustion: An extended view, Physical Review E 60, 518 (1999).

[35] F. Osselin, P. Kondratiuk, A. Budek, O. Cybulski, P. Garstecki, and P. Szymczak, Microfluidic observation of the onset of reactive-infiltration instability in an analog fracture, Geophysical Research Letters 43, 6907 (2016).

[36] K. Khalturin, C. Shinzato, M. Khalturina, M. Hamada, M. Fujie, R. Koyanagi, M. Kanda, H. Goto, F. Anton-Erxleben, M. Toyokawa, S. Toshino, and N. Satoh, Medu-sozoan genomes inform the evolution of the jellyfish body plan, Nature Ecology & Evolution 3, 811 (2019).

[37] S. Nijdam, J. Teunissen, and U. Ebert, The physics of streamer discharge phenomena, Plasma Sources Science and Technology 29, 103001 (2020).

[38] E. M. van Veldhuizen and W. R. Rutgers, Pulsed positive corona streamer propagation and branching, Journal of Physics D: Applied Physics 35, 2169 (2002).

[39] G. J. J. Winands, Z. Liu, A. J. M. Pemen, E. J. M. van Heesch, K. Yan, and E. M. van Veldhuizen, Temporal development and chemical efficiency of positive streamers in a large scale wire-plate reactor as a function of voltage waveform parameters, Journal of Physics D: Applied Physics 39, 3010 (2006).

[40] T. M. P. Briels, J. Kos, E. M. van Veldhuizen, and U. Ebert, Circuit dependence of the diameter of pulsed positive streamers in air, Journal of Physics D: Applied Physics 39, 5201 (2006).

[41] S. Żukowski, P. Morawiecki, H. J. Seybold, and P. Szymczak, Through history to growth dynamics: Deciphering the evolution of spatial networks, Scientific Reports 12, 20407 (2022).

[42] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Shepard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).

[43] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Papat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental algorithms for scientific computing in python, Nature Methods 17, 261 (2020).

[44] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing in Science & Engineering 9, 90 (2007).

[45] Y. Cohen, O. Devauchelle, H. J. Seybold, R. S. Yi, P. Szymczak, and D. H. Rothman, Path selection in the growth of rivers, Proceedings of the National Academy of Sciences 112, 14132 (2015).