# Chaotic advection in a steady three-dimensional MHD flow

Julien Fontchastagner, Jean-François Scheid, Jean-Régis Angilella, Jean-Pierre Brancher

► To cite this version:

Julien Fontchastagner, Jean-François Scheid, Jean-Régis Angilella, Jean-Pierre Brancher. Chaotic advection in a steady three-dimensional MHD flow. Physical Review Fluids, 2025, 10 (114101), pp.24. 10.1103/6ws1-q767. hal-04875268v5

HAL Id: hal-04875268

https://inria.hal.science/hal-04875268v5

Submitted on 3 Nov 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
# Chaotic advection in a steady three-dimensional MHD flow

Julien Fontchastagner,¹, * Jean-François Scheid,², † Jean-Régis Angilella,³, ‡ and Jean-Pierre Brancher⁴, §

¹Université de Lorraine, GREEN, F-54000 Nancy, France

²Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France

³Normandie Université, UNICAEN, UNIROUEN, ABTE, F-14000 Caen, France

⁴Université de Lorraine, CNRS, IECL, F-54000 Nancy, France

We investigate the 3D stationary flow of a weakly conducting fluid in a cubic cavity, driven by the Lorentz force created by two permanent magnets and a weak constant current. Our goal is to determine the conditions leading to efficient mixing within the cavity. The flow is composed of a large recirculation cell created by one side magnet, superposed to two recirculation cells created by a central magnet perpendicular to the first one. The overall structure of this flow, obtained here by solving the Stokes equations with Lorentz forcing, is similar to the tri-cellular model flow studied by Toussaint et. al. (Phys. Fluids. 7, 1995). Chaotic advection in this flow is analyzed by means of Poincaré sections, Lyapunov exponents and expansion entropies. In addition, we quantify the quality of mixing by computing contamination rates, homogeneity, as well as mixing times. Though individual vortices have poor mixing properties, the superposition of both flows creates chaotic streamlines and efficient mixing.

## I. INTRODUCTION AND MOTIVATION

### A. General considerations

Mixing devices at small Reynolds number have been studied during the past fifty years. These are challenging situations where turbulence cannot be used to produce random motion [1, 2]. The low Reynolds number condition is generally met when viscosity is large, or dimensions are small (e.g. microfluidic devices), or when fluids are fragile and require moderate stirring. Nevertheless, high quality mixing can be achieved even for such fluids, provided the flow exhibits chaotic trajectories corresponding to efficient stretching and folding mechanisms [1, 3, 4].

Chaotic advection has been studied and exploited in two-dimensional time-dependent flows, such as flows between eccentric cylinders [1, 5], flows induced by moving walls in a cavity [6], or electromagnetically driven flows [7]. In three dimensions, the pioneering works of Dombre et al. [8] have shown that the flow does not need to be time-dependent to produce chaotic advection. Indeed, streamlines in steady 3D flows can be chaotic, so that fluid points trajectories, which are identical to streamlines in steady flows, are extremely sensitive to initial positions. This property was exploited in a steady open pipe flow with electrically conducting fluids driven from outside by high frequency helical inductors [9, 10], as well as in steady flows in a cubic periodic box [11, 12]. The authors of the latter reference considered a simple arbitrary analytic flow composed of a steady vortex superposed to a vortex pair perpendicular to the first one. In the present paper, a flow with a similar tri-cellular structure is considered. However, instead of using a simple analytical model, our flow is obtained by means of magnetohydrodynamic (MHD) forcing, involving two magnets to create the system of vortices, and will be calculated by solving the Stokes equation with Lorentz forcing. This is a first step towards an experimental facility.

### B. Qualitative description of the flow: side magnet and central magnet

To create such a tri-cellular flow in a realistic way, we consider a viscous fluid with low electric conductivity (e.g. salted water), confined within a cubic cavity. It is submitted to a magnetic field induced by two magnets located outside the box, and to a uniform horizontal current with density $j_0$ passing through the cavity (FIG. 1(c)). The axis $z$ will be denoted as "vertical" in the following.

In the case where a uniform vertical magnetic field $h_1$ is applied over the right-hand-side of the cavity, and set to zero in the other half (FIG. 1(a)), one can check that the liquid cannot remain at rest, and that a flow will appear [9]. Indeed, if the fluid is at rest, the force balance imposes $f = \text{grad } p$ everywhere, where $f$ is the Lorentz force and

<footnote>* julien fontchastagner@univ-lorraine.fr</footnote>

<footnote>† jean-francois.scheid@univ-lorraine.fr</footnote>

<footnote>‡ jean-regis.angilella@unicaen.fr</footnote>

<footnote>§ jean-pierre brancher@univ-lorraine.fr</footnote>
p is the pressure, so that the curl of the Lorentz force must be zero everywhere. This condition cannot be fulfilled when the magnetic field is present only over a portion of the volume of the cavity: here the curl of the Lorentz force is non-zero along the surface of discontinuity of the magnetic field, the hydrostatic condition is not fulfilled, and motion will appear. Under the combined effect of the Lorentz force and of the walls, a large recirculation cell as the one sketched in FIG. 1 (a) will appear. This situation corresponds to the case where a side magnet like the blue one of FIG. 1 (c) is placed near the box. The corresponding Lorentz force will be denoted $\mathbf{f}_1$ in the following. Similarly, a vertical force $\mathbf{f}_2$ can be created near the mid-plane of the cubic box by means of a magnetic field $\mathbf{h}_2$ centered around this mid plane. It will create a double recirculation cell, as sketched in FIG. 1 (b). This can be achieved by placing a centered magnet like the red one of FIG. 1 (c).

The two forces $\mathbf{f}_1$ and $\mathbf{f}_2$ can be superimposed to generate a complex flow, in view of creating a simple MHD mixing device involving static sources only. As discussed above, the structure of this flow is very close to one investigated in Refs. [11] and [12], which involved a simple analytic velocity field. Note however that our flow has no-slip boundary conditions at the walls of the cubic box, as expected for viscous flows, in contrast with Refs. [11] and [12] where fluid is allowed to slip along the faces of the cubic box.

In this article, we will make use of the two configurations described above, except that the magnetic fields will not be taken piecewise constant, but will be issued from a real physical device and computed from Maxwell's equations. Such a device, based solely on the use of two pairs of permanent magnets, is illustrated in FIG. 1 (c), and described in the next section. The velocity field in the fluid will then be computed numerically with these magnetic forces by means of a vorticity – streamfunction formulation of the Stokes equations.

The equations of the coupled system under study for the magnetic and velocity fields are presented in the next chapter. Then, numerical simulations will be performed for various parameters of interest. The magnetic force that will be imposed to the fluid is a weighted average of the two forces $\mathbf{f}_1$ (the Lorentz force due to the side magnet) and $\mathbf{f}_2$ (due to the central magnet), and the ratio of the two will be the control parameter of the device. It will be shown that chaotic advection is obtained for specific choices of this ratio. The computation of Poincaré sections, Lyapunov exponents and expansion entropies will provide qualitative and quantitative information about the chaotic streamlines we obtain. A brief study of mixing efficiency is presented in section IV C.

## II. MOTION EQUATIONS AND CALCULATION OF THE FLOW

### A. General considerations

A cubic tank of side length $a$ containing an incompressible Newtonian fluid with kinematic viscosity $\nu$ is considered. A uniform current density $\mathbf{j}_0$ (directed along the unit vector $\mathbf{e}_x$) is imposed inside the tank by means of a constant voltage drop. As described in the previous section, the two vortex configurations are obtained by means of two pairs of magnets placed in two orthogonal directions (FIG. 1(c)). The side magnet (magnet #1, in red) is magnetized along $\mathbf{e}_z$ and is used to create the force $\mathbf{f}_1$. The central magnet (magnet #2, in blue) is magnetised along $\mathbf{e}_y$ and produces the force $\mathbf{f}_2$. The length and the depth of the magnets are equal to $a$ whereas their width is equal to $a/2$. As explained in the introduction, the force $\mathbf{f}_1$ alone would produce a single recirculation occupying the whole domain, and the force $\mathbf{f}_2$ alone would produce a double-cell. It is the combination of both forces that will produce a complex flow.

Throughout the paper, the Reynolds number based on the size of the box and on the typical fluid velocity is assumed to be small. The value of the current density is chosen to be sufficiently low compared to the magnetization of the permanent magnets, so that its impact on the total magnetic field, and therefore on the total corresponding force, can be neglected. This assumption was verified by a posteriori simulations with and without current density contribution. The $L_2$ norm of the relative difference was less than 0.1% for the typical values of the source terms that were considered.

### B. Governing equations

Following to the previous assumptions, the flow corresponds to a Stokes flow, governed by

$$ \left\{ \begin{array}{l} \nu \Delta \mathbf{v} - \mathrm{grad} \, \mathbf{p} + \mathbf{f} = \mathbf{0} \\ \mathrm{div} \, \mathbf{v} = \mathbf{0} \end{array} \right. , \quad \text{in } \Omega_f \qquad (1) $$

where $\mathbf{v}$ is the velocity field, $\mathbf{p}$ is the pressure field (including gravity effects if any), $\mathbf{f}$ is the Lorentz force density, and $\Omega_f$ denotes the fluid domain. The fluid being viscous, no-slip conditions are imposed on the tank boundary $\Gamma_f$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Upper graphs: sketches of the flows when a single magnet is present. (a): single-cell in the (x, y) plane induced by a side magnet. (b): double-cell in the (x, z) plane created by a central magnet. Current density is shown in black, force fields are shown in gray. The corresponding magnetic fields **h**₁ (configuration (a)) and **h**₂ (configuration (b)) are shown respectively in red and blue. Velocities **v**₁ and **v**₂ are in dark red and dark blue. Lower graph (c): Setup of the mixing MHD device creating these configurations.</figcaption></figure>
for the velocity field:

$$ \mathbf{v} = \mathbf{0} \text{ on } \Gamma_f. \qquad (2) $$

The contribution of each pair of magnets can be treated separately according to the superposition principle:

$$ \mathbf{f} = \mu_0 \mathbf{j}_0 \times \mathbf{h} = \underbrace{\mu_0 \mathbf{j}_0 \times \mathbf{h}_1}_{\mathbf{f}_1} + \underbrace{\mu_0 \mathbf{j}_0 \times \mathbf{h}_2}_{\mathbf{f}_2} \qquad (3) $$

where $\mathbf{h}_i$, $i = \{1, 2\}$ is the field produced by magnet #i. In this way, due to the linearity of the Stokes equations, the velocity field $\mathbf{v}$ in (1) can be decomposed into

$$ \mathbf{v} = \mathbf{v}_1 + \mathbf{v}_2 \qquad (4) $$

where each part $\mathbf{v}_1$ and $\mathbf{v}_2$ satisfies:

$$ \left\{ \begin{array}{ll} \nu \Delta \mathbf{v}_i - \mathrm{grad} \, p_i + \mu_0 \mathbf{j}_0 \times \mathbf{h}_i = \mathbf{0} & \text{in } \Omega_f \\ \mathrm{div} \, \mathbf{v}_i = 0 & \text{in } \Omega_f, \text{ for } i = \{1, 2\} \\ \mathbf{v}_i = \mathbf{0} & \text{on } \Gamma_f \end{array} \right. \qquad (5) $$

The magnetic field $\mathbf{h}_i$ and flux density $\mathbf{b}_i$ in the overall domain $\Omega (= \mathbb{R}^3)$ are governed by the Maxwell-Ampère and Maxwell-Thomson equations:

$$ \mathrm{curl} \, \mathbf{h}_i = \mathbf{0}, \quad \mathrm{div} \, \mathbf{b}_i = 0, \qquad (6) $$

and both fields vanish far from the magnets ($\|\mathbf{x}\| \to \infty$). The overall domain $\Omega$ corresponds to the air containing the fluid domain (the tank) and the pairs of magnets (see Appendix and FIG.2 for the precise overall domain used for computations).

The source term at the origin of the magnetic flux density is due to the magnetization of the magnets, $\mathbf{m}_1 = m_0 \mathbf{e}_z$ for magnet #1 and $\mathbf{m}_2 = m_0 \mathbf{e}_x$ for magnet #2 (see FIG. 1(c)). Then, magnetization vectors have the same norm, $m_0$. The magnetic behaviour law linking the flux densities and the magnetic fields is then

$$ \left\{ \begin{array}{ll} \mathbf{b}_i = \mu_0 (\mathbf{h}_i + \mathbf{m}_i) & \text{in } \Omega_{mi} \\ \mathbf{b}_i = \mu_0 \mathbf{h}_i & \text{elsewhere} \end{array} \right. , \text{ for } i = \{1, 2\} \qquad (7) $$

Extending magnet magnetizations as piecewise functions on the whole domain:

$$ \mathbf{m}_i = \begin{cases} \mathbf{m}_i & \text{in } \Omega_{mi} \\ \mathbf{0} & \text{in } \Omega \setminus \Omega_{mi} \end{cases}, \qquad (8) $$

the magnetic fields $\mathbf{h}_1$ and $\mathbf{h}_2$ are finally governed by

$$ \left\{ \begin{array}{ll} \mathrm{curl} \, \mathbf{h}_i = \mathbf{0} & \text{in } \Omega \\ \mathrm{div}(\mathbf{h}_i + \mathbf{m}_i) = 0 & \text{in } \Omega, \\ \lim_{\|\mathbf{x}\| \to \infty} \mathbf{h}_i = \mathbf{0} & \end{array} \right. , \text{ for } i = \{1, 2\} \qquad (9) $$

These equations will be set non-dimensional as follows. The side length of the cubic tank $a$ is chosen as the characteristic length of the problem, so that the non-dimensional position vector is defined as: $\hat{\mathbf{x}} = \mathbf{x}/a$. The norm $j_0$ of the current density can be expressed from the total current $I_0$ flowing through the system so that the non-dimensional current density can be defined by

$$ \hat{\mathbf{j}}_0(\hat{\mathbf{x}}) = \frac{\mathbf{j}_0(\mathbf{x})}{j_0} = \frac{a^2}{I_0} \mathbf{j}_0(\mathbf{x}). \qquad (10) $$

In the same way, a non-dimensional magnetization vector $\hat{\mathbf{m}}_i$, and the corresponding $\hat{\mathbf{h}}_i$ field are introduced as follows:

$$ \hat{\mathbf{m}}_i(\hat{\mathbf{x}}) = \frac{\beta}{m_0} \mathbf{m}_i(\mathbf{x}), \quad \hat{\mathbf{h}}_i(\hat{\mathbf{x}}) = \frac{\beta}{m_0} \mathbf{h}_i(\mathbf{x}), \quad \text{for } i = \{1, 2\} \qquad (11) $$
where $\beta$ is an arbitrary scaling factor designed to avoid too low values of magnetic fields (and therefore forces and velocities) when solving the problem. In terms of these non-dimensional quantities, Eq. (9) now reads

$$ egin{cases} \mathrm{curl} \hat{\mathbf{h}}_i = \mathbf{0} & \text{in } \hat{\Omega} \quad (\Omega = a \hat{\Omega}) \\ \mathrm{div} (\hat{\mathbf{h}}_i + \hat{\mathbf{m}}_i) = 0 & \text{in } \hat{\Omega} \\ \lim_{\|\hat{\mathbf{x}}\| \to \infty} \hat{\mathbf{h}}_i = \mathbf{0} & \text{, for } i = \{1, 2\} \end{cases} \quad (12) $$

Velocity and pressure are set non-dimensional as:

$$ \hat{\mathbf{v}}_i(\hat{\mathbf{x}}) = \frac{\mathbf{v}_i(\mathbf{x})}{v_0}, \quad \hat{p}_i(\hat{\mathbf{x}}) = \frac{p_i(\mathbf{x})}{p_0}, \quad \text{for } i = \{1, 2\} \quad (13) $$

where $v_0$ and $p_0$ are characteristic velocity and pressure that will be defined below. The fluid equations becomes

$$ \Delta \hat{\mathbf{v}}_i - \frac{a p_0}{\nu v_0} \mathrm{grad} \hat{p}_i + \frac{\mu_0 I_0 m_0}{\nu v_0 \beta} \hat{\mathbf{j}}_0 \times \hat{\mathbf{h}}_i = \mathbf{0}, \quad \text{for } i = \{1, 2\}. \quad (14) $$

By choosing

$$ v_0 = \frac{\mu_0 I_0 m_0}{\nu \beta} \quad \text{and} \quad p_0 = \frac{\mu_0 I_0 m_0}{a \beta}, \quad (15) $$

the motion equation finally reads:

$$ egin{cases} \Delta \hat{\mathbf{v}}_i - \mathrm{grad} \hat{p}_i + \hat{\mathbf{j}}_0 \times \hat{\mathbf{h}}_i = \mathbf{0} & \text{in } \hat{\Omega}_f \quad (\Omega_f = a \hat{\Omega}_f) \\ \mathrm{div} \hat{\mathbf{v}}_i = 0 & \text{in } \hat{\Omega}_f \\ \hat{\mathbf{v}}_i = \mathbf{0} & \text{on } \hat{\Gamma}_f \end{cases}, \quad \text{for } i = \{1, 2\}. \quad (16) $$

Finally, by sequentially solving Eqs. (12) and (16), for $i=1,2$, we obtain the two velocity fields $\hat{\mathbf{v}}_1$ and $\hat{\mathbf{v}}_2$.

In the remainder of this article, we will be concerned only with non-dimensional quantities, and for convenience the hat symbols ($\hat{\cdot}$) will be omitted.

### III. NUMERICAL COMPUTATIONS OF FIELDS

To avoid projection or interpolation errors due to the use of different meshes or function spaces, we have chosen to solve the governing equations by means of the same computational tool. Mixed finite elements methods, based on both nodal and edge basis functions, provide an efficient and convenient framework for solving problems of low-frequency computational electromagnetism [13–15]. However, they can also be used to correctly approximate the Stokes system [13, 16] as described in section III B.

#### A. Magnetic fields: scalar magnetic potential formulation

First of all, we should point out that the magnetic fields produced by the magnets (and only these fields) can be calculated purely analytically in our particular case. Indeed, analytical formulas exist for rectangular magnets in air, see [17] or more recently [18]. However, the authors preferred to use a numerical model, as this will be necessary in our future work in order to design an experimental test bench. To reduce the total volume of magnets required, ferromagnetic yokes will be added to each pair of magnets to loop their magnetic flux and reinforce the field inside the tank. In this case, an analytical model is no longer possible. However, we compared both approaches on our study case, and the results from the analytical and numerical models were perfectly consistent. No significant difference were observed on the final velocities. This allows us to validate the numerical calculation of the magnetic fields detailed below.

The overall domain of computation $\Omega$ can be reduced to $\Omega = \Omega_f \cup \Omega_{m1} \cup \Omega_{m2} \cup \Omega_a \cup \Omega_\infty$ which is represented in FIG. 2. $\Omega_f$ corresponds to the tank, $\Omega_{m1}$ and $\Omega_{m2}$ to the pairs of magnets, and $\Omega_a$ to the air region around them. Their union forms the domain of interest $\Omega_i = \Omega_f \cup \Omega_{m1} \cup \Omega_{m2} \cup \Omega_a$ in which the magnetic fields are wanted. The external boundary of $\Omega_i$ is that of a sphere of radius $R_i$. This region is surrounded by a spherical finite shell $\Omega_\infty$ with
<figure><img src="image_3.png" /><figcaption>FIG. 2. Representation of the total computational domain, and the corresponding mesh used for the numerical computations (control parameter ne = 15).</figcaption></figure>

inner and outer radii respectively $R_i$ and $R_o$. This will be useful to deal with the open boundary condition on fields as explained below.

The spaces $L^2(\Omega)$ and $L^2(\Omega)$ can be defined by:

$$L^2(\Omega) = \{u : \Omega \to \mathbb{R} \mid \int_{\Omega} u(\mathbf{x})^2 d\mathbf{x} < \infty\}, \quad \text{and} \quad L^2(\Omega) = \{u : \Omega \to \mathbb{R}^3 \mid \int_{\Omega} \|u(\mathbf{x})\|^2 d\mathbf{x} < \infty\} \quad (17)$$

They permit to define the Sobolev space $H(\text{grad}, \Omega)$, more commonly known as $H^1(\Omega)$:

$$ \begin{aligned} H(\text{grad}, \Omega) &= \{u \in L^2(\Omega) \mid \text{grad } u \in L^2(\Omega)\} \\ &= \{u \in L^2(\Omega) \mid \partial_x u, \partial_y u, \partial_z u \in L^2(\Omega)\} = H^1(\Omega) \end{aligned} \quad (18) $$

The Maxwell-Ampère equation $\text{curl} \mathbf{h}_i = 0$ implies that there exists a magnetic scalar potential $\phi_i \in H(\text{grad}, \Omega)$ such that:

$$ \mathbf{h}_i = - \text{grad} \phi_i \quad (19) $$

To deal with the $\mathbf{h}_i \xrightarrow{\|\mathbf{x}\| \to \infty} 0$ condition, a bijective transformation that maps all the region $\mathbb{R}^3 \setminus \Omega_i$ onto $\Omega_\infty$ is constructed thanks to the relation:

$$ \|\mathbf{x}\| = R_i \frac{R_o - R_i}{R_o - \|\mathbf{X}\|} \quad (20) $$
where $\mathbf{X}$ is the position vector in the shell $\Omega_\infty$ [19]. The considered transformation $\mathcal{F}: \mathbb{R}^3\setminus\Omega_i \to \Omega_\infty$ is therefore defined by:

$$ \mathbf{X} = \mathcal{F}(\mathbf{x}) = \frac{\|\mathbf{X}\|}{\|\mathbf{x}\|} \mathbf{x} = \left[ \frac{R_o}{\|\mathbf{x}\|} - \frac{R_i(R_o - R_i)}{\|\mathbf{x}\|^2} \right] \mathbf{x} \quad (21) $$

The open-boundary conditions on $\mathbf{h}_1$ and $\mathbf{h}_2$ can be replaced by imposing homogeneous Dirichlet conditions on $\phi_1$ and $\phi_2$ on $\Gamma_\infty = \partial\Omega$ (the sphere of radius $R_o$). We will also use the transformation $\mathcal{F}$ and its jacobian matrix $\mathbf{J}_\mathcal{F}$ to bring back all integrals posed on $\mathbb{R}^3 \setminus \Omega_i$ to $\Omega_\infty$.

The second equation of $(S_{m,i})$ becomes:

$$ \mathrm{div}(\mathrm{grad} \, \phi_i - \mathbf{m}_i) = 0 \quad \text{in } \Omega \quad (22) $$

Multiplying (22) by a test function $\xi_i \in H_0(\mathrm{grad}, \Omega) = \{u \in H(\mathrm{grad}, \Omega) : u|_{\Gamma_\infty} = 0\}$, then using the transformation $\mathcal{F}$ to write down the integrals terms over $\mathbb{R}^3 \setminus \Omega_i$ to $\Omega_\infty$ and finally integrating by parts, we obtain the weak form $(\Sigma_{m,i})$ of the system $(S_{m,i})$ given by (12):

$$ (\Sigma_{m,i}) \quad \left\{ \begin{array}{ll} \displaystyle \int_{\Omega_i} \mathrm{grad} \, \phi_i \cdot \mathrm{grad} \, \xi_i \, d\mathbf{x} - \int_{\Omega_{m,i}} \mathbf{m}_i \cdot \mathrm{grad} \, \xi_i \, d\mathbf{x} & (23) \\ \displaystyle + \int_{\Omega_\infty} \mathbf{J}_\mathcal{F}^{-1} \mathrm{grad} \, \phi_i \cdot \mathbf{J}_\mathcal{F}^{-1} \mathrm{grad} \, \xi_i \, |\det \mathbf{J}_\mathcal{F}| \, d\mathbf{x} = 0 & \end{array} \right. $$

Once $\phi_i$ is computed from (23), the force density involved in the non-dimensional Stokes flow system in $\Omega_f$ will be:

$$ \mathbf{f}_i = -\mathbf{j}_0 \times \mathrm{grad} \, \phi_i \quad (24) $$

## B. Flows velocities: vorticity – streamfunction formulation

Let us begin by defining the spaces:

$$ \mathbf{H}(\mathrm{curl}, \Omega_f) = \{\mathbf{u} \in \mathbb{L}^2(\Omega_f) \mid \mathrm{curl} \, \mathbf{u} \in \mathbb{L}^2(\Omega_f)\} \quad (25) $$

$$ \mathbf{H}_0(\mathrm{curl}, \Omega_f) = \{\mathbf{u} \in \mathbf{H}(\mathrm{curl}, \Omega_f) \mid \mathbf{u} \times \mathbf{n}|_{\Gamma_f} = \mathbf{0}\} \quad (26) $$

where $\Gamma_f$ is the boundary of $\Omega_f$ ($\Gamma_f = \partial\Omega_f$). It is possible to rewrite the Stokes system (16) in terms of the vorticity fields $\omega_i \in \mathbf{H}(\mathrm{curl}, \Omega_f)$ defined by $\omega_i = \mathrm{curl} \, \mathbf{v}_i$:

$$ \mathrm{curl} \, \omega_i + \mathrm{grad} \, p_i - \mathbf{f}_i = \mathbf{0} \quad (27) $$

And the fact that $\mathbf{v}_i = \mathbf{0}$ on $\Gamma_f$ directly implies:

$$ 0 = \mathrm{curl} \, \mathbf{v}_i \cdot \mathbf{n} = \omega_i \cdot \mathbf{n} \quad \text{on } \Gamma_f \quad (28) $$

Since $\mathbf{v}_i$ is divergence free, there exists a vector potential $\psi_i \in \mathbf{H}_0(\mathrm{curl}, \Omega_f)$ such that [20]:

$$ \mathbf{v}_i = \mathrm{curl} \, \psi_i, \quad \text{and} \quad \mathrm{div} \, \psi_i = 0 \quad (29) $$

By using this potential and applying the curl operator to the first equation of the system $(S_{f,i})$, we obtain:

$$ \left\{ \begin{array}{ll} \mathrm{curl} \, \mathrm{curl} \, \omega_i = \mathrm{curl} \, \mathbf{f}_i, & \mathrm{div} \, \omega_i = 0 \quad \text{in } \Omega_f \\ \omega_i = \mathrm{curl} \, \mathrm{curl} \, \psi_i, & \mathrm{div} \, \psi_i = 0 \quad \text{in } \Omega_f \\ \omega_i \cdot \mathbf{n} = 0, \quad \psi_i \times \mathbf{n} = \mathbf{0}, \quad \mathrm{curl} \, \psi_i \times \mathbf{n} = \mathbf{0} & \text{on } \Gamma_f \end{array} \right. \quad (30) $$

The vorticity $\omega_i$ can be decomposed as $\omega_i = \omega_i^0 + \omega_i^*$ where $\omega_i^0 \in \mathbf{H}_0(\mathrm{curl}, \Omega_f)$ is a regular part and $\omega_i^* \in \mathbf{H}(\mathrm{curl}, \Omega_f)$ is the harmonic contribution [21, 22]. We rewrite the global system in two weakly coupled subsystems as follows:

$$ \left\{ \begin{array}{ll} -\Delta\omega_i^0 = \mathrm{curl} \, \mathrm{curl} \, \omega_i^0 = \mathrm{curl} \, \mathbf{f}_i & \\ \mathrm{div} \, \omega_i^0 = 0 & \\ \omega_i^0 \times \mathbf{n}|_{\Gamma_f} = \mathbf{0} & , \end{array} \right. \quad \left\{ \begin{array}{ll} -\Delta\psi_i = \mathrm{curl} \, \mathrm{curl} \, \psi_i = \omega_i^0 + \omega_i^* & \\ -\Delta\omega_i^* = \mathrm{curl} \, \mathrm{curl} \, \omega_i^* = \mathbf{0} & \\ \mathrm{div} \, \omega_i^* = 0, \quad \mathrm{div} \, \psi_i = 0 & \\ (\omega_i^0 + \omega_i^*) \cdot \mathbf{n}|_{\Gamma_f} = \mathbf{0} & \\ \psi_i \times \mathbf{n}|_{\Gamma_f} = \mathbf{0}, \quad \mathrm{curl} \, \psi_i \times \mathbf{n}|_{\Gamma_f} = \mathbf{0} & \end{array} \right. \quad (31) $$
$$ (\Sigma_{f,i}) \left\{ \begin{array}{l} \text{Find } (\omega_i^0, \lambda_i^0) \in \mathbf{H}_0(\text{curl}, \Omega_f) \times \text{H}(\text{grad}, \Omega_f) \text{ such that:} \\ \left(\Sigma_{f,i}^{(1)}\right) \begin{array}{l} \displaystyle \int_{\Omega_f} \text{curl} \omega_i^0 \cdot \text{curl} \varphi_i \, dx + \int_{\Omega_f} \text{grad} \lambda_i^0 \cdot \varphi_i \, dx = \\ \displaystyle \int_{\Omega_f} \mathbf{f}_i \cdot \text{curl} \varphi_i \, dx, \quad \forall \varphi_i \in \mathbf{H}_0(\text{curl}, \Omega_f) \\ \displaystyle \int_{\Omega_f} \omega_i^0 \cdot \text{grad} \mu_i \, dx = 0, \quad \forall \mu_i \in \mathbf{H}_0(\text{grad}, \Omega_f) \end{array} \\ \\ \left(\Sigma_{f,i}^{(2)}\right) \begin{array}{l} \text{Find } (\omega_i^*, \lambda_i^*) \in \mathbf{H}(\text{curl}, \Omega_f) \times \text{H}(\text{grad}, \Omega_f) \text{ and} \\ \quad (\psi_i, \eta_i) \in \mathbf{H}_0(\text{curl}, \Omega_f) \times \text{H}(\text{grad}, \Omega_f) \text{ such that:} \\ \displaystyle \int_{\Omega_f} \omega_i^* \cdot \varphi_i \, dx - \int_{\Omega_f} \text{curl} \psi_i \cdot \text{curl} \varphi_i \, dx + \int_{\Omega_f} \text{grad} \eta_i \cdot \varphi_i \, dx = \\ \displaystyle - \int_{\Omega_f} \omega_i^0 \cdot \varphi_i \, dx, \quad \forall \varphi_i \in \mathbf{H}(\text{curl}, \Omega_f) \\ \displaystyle \int_{\Omega_f} \text{curl} \omega_i^* \cdot \text{curl} \xi_i \, dx + \int_{\Omega_f} \text{grad} \lambda_i^* \cdot \xi_i \, dx = 0, \quad \forall \xi_i \in \mathbf{H}_0(\text{curl}, \Omega_f) \\ \displaystyle \int_{\Omega_f} \psi_i \cdot \text{grad} \chi_i \, dx = 0, \quad \forall \chi_i \in \mathbf{H}_0(\text{grad}, \Omega_f) \\ \displaystyle \int_{\Omega_f} (\omega_i^0 + \omega_i^*) \cdot \text{grad} \mu_i \, dx = 0, \quad \forall \mu_i \in \mathbf{H}(\text{grad}, \Omega_f) \end{array} \end{array} \right. \quad (32) $$

For each subsystem, a mixed formulation can be deduced by introducing the scalar fields $\lambda_i^0$, $\lambda_i^*$ and $\eta_i$ which permit respectively to weakly impose the divergence free properties of $\omega_i^0$, $\omega_i^*$ and $\psi_i$, by acting as Lagrange multipliers. The systems $(S_{f,i})$ are then equivalent to the following weak form $(\Sigma_{f,i})$, given by (32) (see [21]).

It is important to mention that the condition $(\omega_i^0 + \omega_i^*) \cdot \mathbf{n} = 0$ on $\Gamma_f$ is implicitly imposed by choosing the test functions $\mu_i$ in $\mathbf{H}(\text{grad}, \Omega_f)$; and for the boundary condition $\text{curl} \psi_i \times \mathbf{n} = 0$ on $\Gamma_f$, by choosing in $(\Sigma_{f,i}^{(2)})$ the test functions $\varphi_i$ in $\mathbf{H}(\text{curl}, \Omega_f)$ and not in $\mathbf{H}_0(\text{curl}, \Omega_f)$.

By sequently solving $(\Sigma_{m,i})$, then $(\Sigma_{f,i}^{(1)})$, and lastly $(\Sigma_{f,i}^{(2)})$ for $i = \{1, 2\}$, the desired velocities are finally obtained by $\mathbf{v}_1 = \text{curl} \psi_1$ and $\mathbf{v}_2 = \text{curl} \psi_2$. The practical implementation with high order basis functions and discrete form of the full system solved are detailed in Appendix. The plots of the calculated magnetic fields, the corresponding forces and the resulting velocities (when a single magnet is present) are provided by FIG. 3.

As expected, the magnetic fields created by the magnets remain concentrated in the area between each of them, and so does the electromagnetic forces. The structure of these velocity fields corresponds to the sketches of FIG. 1, with one recirculation cell created by the side magnet #1 and two cells created by the central magnet #2. To make this clearer, a few streamlines of $\mathbf{v}_1$ and $\mathbf{v}_2$ are also displayed in FIG. 4. The motion of tracers in the steady flow obtained by involving the two magnets, i.e. by superposing $\mathbf{v}_1$ and $\mathbf{v}_2$, is studied in the next section.

Once the computation is complete, the two velocity fields are exported on a regular grid composed of 100 × 100 × 100 points covering the whole tank domain $\Omega_f$. Values are computed thanks to our high order basis functions in order to keep the best possible accuracy. These export files constitute the final results of the numerical procedure described in this section, and are the starting point of the study of the chaotic behaviour detailed in the following.

# IV. LAGRANGIAN CHAOS AND MIXING ANALYSIS

## A. Highlighting and quantifying chaos

### 1. Trajectories and sensitivity to initial conditions

In the following, we will assume that it will always be possible to modulate the velocity fields $\mathbf{v}_1$ and $\mathbf{v}_2$ (e.g. by adjusting the magnet length), so that the total velocity $\mathbf{v}$ may be written as :

$$ \mathbf{v} = \alpha \mathbf{v}_1 + (1 - \alpha) \mathbf{v}_2, \quad \text{for } \alpha \in [0, 1] \quad (33) $$
<figure><img src="image_4.png" /></figure>

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

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

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

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

<figure><img src="image_9.png" /><figcaption>FIG. 3. Plots of calculated non-dimensional fields, when a single magnet is present: **h**₁ and **h**₂ are the magnetic fields, **f**₁ and **f**₂ are the corresponding Lorentz forces, **v**₁ and **v**₂ are the resulting velocities (see Appendix for details).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 4. Examples of velocity streamlines when a single magnet is present. Left: $v_1$ (side magnet). Right: $v_2$ (central magnet).</figcaption></figure>

where $\alpha$ is the convex combination coefficient of the two velocity fields. To find out the value of this velocity field at any point $\mathbf{x} = (x, y, z)$ within the tank domain $[-0.5, 0.5]^3$, we simply interpolate the values calculated in the previous section by means of cubic splines. In the following, these interpolated functions will be still denoted $v_1$ and $v_2$.

The key point of this part is to find out the parameter $\alpha$ leading to the best possible mixing. To do this, we will study the nonlinear dynamical system corresponding to the trajectory of a tracer particle with position $\mathbf{x}(t)$, initially released at $\mathbf{x}_0$. In our case, this reads as:

$$ \left\{ \begin{array}{l} \frac{d\mathbf{x}}{dt} = \dot{\mathbf{x}}(t) = \alpha \mathbf{v}_1(\mathbf{x}(t)) + (1-\alpha) \mathbf{v}_2(\mathbf{x}(t)), t \ge 0 \\ \mathbf{x}(0) = \mathbf{x}_0 = (x_0, y_0, z_0). \end{array} \right. \qquad (34) $$

Here also, velocities, lengths and times have been set non-dimensional by means of $v_0$, $a$ and $a/v_0$ respectively.

The ordinary differential system (34) has been implemented with the Julia programming language using the DynamicalSystemsBase.jl package (and its subset ChaosTool.jl) [23]. We choose a solver based on an Adams-Bashforth explicit method: a fixed time step fifth-order Adams-Bashforth-Moulton method, where a 5th order Adams-Bashforth method is used as predictor and an Adams-Moulton 4-steps method is used for the corrector. To calculate starting values, we make use of a Runge-Kutta method of order 4. Throughout the rest of the paper, the non-dimensional time step is fixed to $5 \cdot 10^{-4}$.

To get started, some trajectories are mapped out in order to study the sensitivity to initial conditions. Three very close starting points are chosen: $\mathbf{x}_0 = (0.15, 0.15, 0.15)$, $\mathbf{x}_{01} = \mathbf{x}_0 + \varepsilon$ and $\mathbf{x}_{02} = \mathbf{x}_0 - \varepsilon$, with $\varepsilon = (1, 1, 1) \cdot 10^{-3}$. The early times of the three corresponding trajectories, for six values of $\alpha$ arbitrarily chosen between 0 and 1, are plotted in FIG. 5.

As expected, trajectories correspond to the streamlines of the single-magnet flows (FIG. 4) when $\alpha = 0$ ($\mathbf{v} = \mathbf{v}_2$) and $\alpha = 1$ ($\mathbf{v} = \mathbf{v}_1$). No sensitivity to initial conditions is observed there.

In contrast, other values of $\alpha$ between these two extremes lead to much less regular trajectories (see FIG. 5) and to sensitivity to initial conditions. This property, which characterizes chaos, is particularly pronounced for $\alpha \simeq 0.25$, but is always present for other values of $\alpha$ (except, of course, for $\alpha = 0$ and 1).

It is therefore necessary to study deeper the impact of the coefficient $\alpha$, using other tools such as Poincaré sections and Lyapunov exponents.

## 2. Poincaré sections

For the single starting point $\mathbf{x}_0 = (0.15, 0.15, 0.15)$ we draw the corresponding Poincaré sections in the plane $z = 0$. They are shown in FIG. 6, for eight values of $\alpha$ equally distributed between 0.2 and 0.6. It is important to point out that these plots have been obtained for a unique starting point, and not for a set of several points tightened around the initial condition. In the three cases $\alpha = 0.3, 0.4$ and $0.55$, just a second starting point very close to $\mathbf{x}_0$ has been added to obtain a similar number of impact points for all sections. The high densities of the plots show that the problem is solved with a very good accuracy.
<figure><img src="image_11.png" /><figcaption>FIG. 5. Short-term trajectories for various values of $\alpha$ and for 3 starting points very close to $\mathbf{x}_0 = (0.15, 0.15, 0.15)$. For $\alpha \neq 0$ and 1, trajectories are irregular and sensitive to initial conditions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 6. Influence of the coefficient $\alpha$ on the size and location of tori of periodic trajectories: dead zones represented by a hole in the Poincaré sections ($z = 0, \mathbf{x}_0 = (0.15, 0.15, 0.15)$), the darkest color corresponds to outgoing points and the lightest color to incoming points.</figcaption></figure>
$\alpha = 0.25 / \mathbf{x}_0 = (0.15, 0.15, 0.15)$: orange, $(0.35, 0, 0)$: blue

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

$\alpha = 0.5 / \mathbf{x}_0 = (0.15, 0.15, 0.15)$: orange, $(0.25, 0.25, 0)$: blue, $(-0.25, 0.3, 0)$: green

<figure><img src="image_14.png" /><figcaption>FIG. 7. Trajectories and corresponding $z = 0$ Poincaré sections with highlight of one torus or two tori of periodic trajectories.</figcaption></figure>

First, we observe in FIG. 6 that most of the Poincaré sections are covered, which means that the entire tank is explored by each trajectory. This is reassuring in terms of mixing efficiency. Nonetheless, we note the presence of holes in each case. This seems to indicate the presence of regions of regular trajectories which must take the form of torii of periodic trajectories. Depending on the values of the coefficient $\alpha$, the position of these holes varies, with a $+-$shaped placement for low values of $\alpha$ and a shift towards an $\times$-shape for higher values. The changeover seems to take place around $\alpha = 0.45$.

To better visualize these torii of periodic trajectories, a 3D plot is shown in FIG. 7. For two values of $\alpha$ (0.25 and 0.5), the beginning of trajectories is plotted for the previous starting point ($\mathbf{x}_0 = (0.15, 0.15, 0.15)$) and for points taken in one or two holes. The corresponding Poincaré sections $z = 0$ have been added alongside. The plots clearly show the periodic trajectories inside the orange chaotic trajectory. We can see one torus for the case $\alpha = 0.25$,
<figure><img src="image_15.png" /><figcaption>FIG. 8. Poincaré sections for one single starting point (0.15, 0.15, 0.15) in the case $\alpha = 0.45$. The darkest color corresponds to outgoing points and the lightest color to incoming points.</figcaption></figure>

corresponding to the +-shaped holes. For the case $\alpha = 0.5$, where the holes follow an ×-shape, two tori are present. However, even in the presence of these regular trajectories, the flow may be chaotic from a global point of view. To characterize this, we will calculate some Lyapunov exponents in the next section.

Before that, and to conclude this paragraph, we plot in FIG. 8 a set of Poincaré sections (in three planes in each direction of space), still for the same starting point and in the case $\alpha = 0.45$, which seems to be the one where periodic trajectories are less prominent. Indeed, the plots clearly show that the entire space is almost visited by the particle, suggesting that mixing is efficient.
### 3. Lyapunov exponents

Lyapunov exponents are often used to quantify deterministic chaos. For a given trajectory $\mathbf{x}(t)$ of a chaotic dynamical system, we perturb the initial position $\mathbf{x}(0) = \mathbf{x}_0$ at time $t = 0$. The perturbed trajectory, denoted $\mathbf{y}(t)$ in the following, originates at $\mathbf{y}(0) = \mathbf{y}_0$. As time $t$ goes on, $\mathbf{y}(t)$ moves further away from $\mathbf{x}(t)$ with a distance $\delta(t) \equiv \|\mathbf{x} - \mathbf{y}\| \simeq \delta_0 \exp(\lambda t)$, where $\delta_0 = \|\mathbf{x}_0 - \mathbf{y}_0\|$ and the real number $\lambda$, called Lyapunov exponent, is positive in the case of a chaotic system. Looking for perturbations along the three axes of the frame $(O, x, y, z)$, leads to three Lyapunov exponents $\lambda_1, \lambda_2$ and $\lambda_3$. For the dynamical system (34), we consider a perturbation matrix $Y$ according to the three cardinal axes, which satisfies the linearized dynamical system $\dot{Y}(t) = \nabla\mathbf{v}(\mathbf{x}(t))Y(t)$ where $\mathbf{v}$ is given by (33). The Lyapunov exponents can be computed from a QR-decomposition of the perturbation matrix $Y(t) = Q(t)R(t)$ with $Q$ being an orthogonal matrix and $R$ an upper triangular matrix with positive diagonal coefficients $R_{ii} > 0$ (see [24, Appendix A] and [25]):

$$ \lambda_i = \lim_{t \to +\infty} \frac{1}{t} \ln R_{ii}(t), \quad \text{for } i = 1, 2, 3. \qquad (35) $$

For practical numerical computations of the Lyapunov exponents, we perform an iterated QR-decomposition of $Y$ at discrete times $t_k$ with a time step $\Delta t$. We have $Y(t_k) = Q_k R_k$ with $Y$ satisfying

$$ \left\{ \begin{array}{ll} \dot{Y}((t) = \nabla \mathbf{v}(\mathbf{x}(t))) Y(t), & t \in [t_{k-1}, t_k] \\ Y(t_{k-1}) = Q_{k-1}. & \end{array} \right. \qquad (36) $$

For a final time $T = N\Delta t$, we obtain the iterated QR-decomposition with $Y(T) = Q_N R_N R_{N-1} \dots R_1$ and we compute the Lyapunov exponents according to

$$ \lambda_i \approx \frac{1}{N\Delta t} \sum_{j=1}^{N} \ln (R_j)_{ii}. \qquad (37) $$

The dynamical system (34) will present a chaotic behavior if at least one Lyapunov exponent is positive. Since the velocity $\mathbf{v}$ in (34) is divergence-free, the following property holds:

$$ \lambda_1 + \lambda_2 + \lambda_3 = 0. \qquad (38) $$

We emphasize that the values of the Lyapunov exponents depend on the initial position $\mathbf{x}(0) = \mathbf{x}_0$ of the trajectory. It is known that autonomous dynamical systems always possess a zero Lyapunov exponent provided that the trajectory of the flow does not converge to a steady state. We refer to [12] and [26] §2.5.6 where a sketch of a simple proof is given. In a chaotic region (where the initial position $\mathbf{x}_0$ has been chosen), we then necessarily get from (38)

$$ \lambda_1 > 0, \lambda_2 = 0, \lambda_3 = -\lambda_1 < 0, \qquad (39) $$

while in a smooth or regular region, we expect to have

$$ \lambda_1 = \lambda_2 = \lambda_3 = 0. \qquad (40) $$

This method can easily be called up using the Julia package mentioned previously [23]. We have modified the default implementation to preserve the history of terms calculated during the process. For our starting point $\mathbf{x}_0 = (0.15, 0.15, 0.15)$, and for the two values of $\alpha$ considered in FIG. 7, we plot the Lyapunov exponents $\lambda_i$ versus the number of evolution steps $N$ in FIG. 9.

The long-term shape of these curves agrees with Eq. (39), and the calculation converges fairly quickly, with a final value reached after nearly 7,000 steps. This confirms our choice of the ODE solver and associated time step. In both cases, we find a chaotic behavior with $\lambda_1 > 0$, but it is much more pronounced for $\alpha = 0.25$.

The maximal Lyapunov exponent $\lambda_1$ also brings valuable information for starting points associated with the regular points observed in section IV A 2. It is shown in FIG. 10 as a function of $\alpha$, for the points $x_0 = (0.35, 0, 0)$ (blue dashed line) and $x_0 = (-0.25, 0.3, 0)$ (green dotted line), corresponding respectively to the cases of one torus and two tori in FIG. 7. In addition, the maximal Lyapunov exponent of point $x_0 = (0.15, 0.15, 0.15)$ has been plotted for comparison in FIG. 10 (orange solid line). These three curves confirm the qualitative observations of FIG. 7. Indeed, the blue dashed curve clearly shows a regular periodic trajectory for low values of $\alpha$, becoming chaotic as $\alpha$ increases. The green dotted curve follows the opposite evolution, and the switch between the two behaviors occurs around $\alpha = 0.45$.

To find a more global criterion that reflects the chaotic behavior of the entire flow, we have computed expansion entropies. These are shown in the next section.
<figure><img src="image_16.png" /><figcaption>FIG. 9. Evolution of Lyapunov exponents versus the number of computation steps, for the two previous cases ($\alpha$ = 0.25 and 0.5, $\mathbf{x}_0$ = (0.15, 0.15, 0.15)).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 10. Lyapunov exponents versus $\alpha$, for the trajectories shown in FIG. 7.</figcaption></figure>

## B. Expansion entropy

The behavior of Lyapunov exponents is directly linked to the choice of the initial position of the trajectory. A more satisfactory and comprehensive way to characterize chaos without being too dependent on initial conditions is to consider numerical quantities such as topological entropy or expansion entropy. The topological entropy, with its standard definitions [27], is however quite difficult to compute numerically. This is the reason why we prefer to deal with the expansion entropy whose definition is better suited to numerical evaluation. In particular, expansion entropy does not require to identify a compact invariant set, in contrast with topological entropy. We refer to [28] for a review and a discussion on various chaos indicators and also for the definition and properties of the expansion entropy. An important result states that for autonomous systems on a compact manifold, the topological entropy coincides with the expansion entropy for a regular flow [29].

The expansion entropy is defined as

$$H_0(S) = \lim_{T \to +\infty} \frac{1}{T} \ln E_T(\varphi_T, S) \quad (41)$$
with

$$E_T(\varphi_T, S) = \frac{1}{\mathrm{vol}(S)} \int_{S_T} G(\nabla \varphi_T(\mathbf{x})) \, \mathrm{d}\mathbf{x} \qquad (42)$$

where $\varphi_t$ is the flow associated with dynamical system (34) i.e. $\mathbf{x}(t) = \varphi_t(\mathbf{x}_0)$, $S$ is a subset of the fluid tank $\Omega_f$ and $S_T$ is the set of $\mathbf{x}$ such that $\varphi_t(\mathbf{x}) \in S$ for all $t \in [0, T]$ (the trajectories that remain in $S$ up to time $T$). Finally, $G(A)$ denotes the product of the singular values of a matrix $A$ that are greater than 1 (if none of the singular values are > 1 then we set $G(A) = 1$). We consider that chaos occurs in $S$ if $H_0(S) > 0$. The computations of the expansion entropy are performed in the whole fluid domain, that is with $S = S_T = \Omega_f$.

There is also a link between the Lyapunov exponents and $G(\nabla\varphi_T)$. Indeed, due to the Oseledec's multiplicative ergodic theorem (see [25]) it can be shown [30] that the exponential growth rate of $G(\nabla\varphi_T(\mathbf{x}))$ as $T \to +\infty$ is the sum of the positive Lyapunov exponents of the trajectory starting from $\mathbf{x}$ at time $t = 0$. In a chaotic region, thanks to (39), we have

$$\lim_{T \to \infty} \frac{1}{T} \ln G(\nabla \varphi_T(\mathbf{x})) = \lambda_1. \qquad (43)$$

As a consequence, the limit in (41) can be seen as a space average on $S$ of the sum of the positive Lyapunov exponents of trajectories that remain in $S$.

Numerically, the expansion entropy is computed as follows. First, we randomly generate a uniform distribution of $N$ points $\{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N\}$ in $S = \Omega_f$, then we compute

$$\hat{E}_T(\varphi_T) = \frac{1}{N} \sum_{i=1}^{N} G(\nabla \varphi_T(\mathbf{x}_i)) \qquad (44)$$

which is an estimate of $E_T(\varphi_T, \Omega_f)$ given by (44). Then, we compute the approximate expansion entropy

$$\hat{H}_0 = \frac{1}{T} \ln \hat{E}_T(\varphi_T) \qquad (45)$$

which is an approximation of $H_0(\Omega_f)$ for large $N$ and $T$. In practice, we choose a number $q$ of different batches of $N$ points and we take the average of the approximate expansion entropy computed with each batch. We use the Julia package `ChaosTools.jl` with the function `expansionentropy` which computes $\hat{H}_0$ from $q$ batches.

In FIG. 11, the computed expansion entropy $\hat{H}_0$ is depicted with respect to the parameter $\alpha$. We chose $N = 1000$ points uniformly distributed in $\Omega_f$ and $q = 10$ batches. We see that the expansion entropy is always positive and higher for values of $\alpha$ less than 0.5. A peak with a maximum of $\hat{H}_0$ appears around $\alpha = 0.14$. The chaotic behavior of the system is therefore more pronounced around this value of $\alpha$. We also point out that for some initial conditions, the Lyapunov exponents behave like the expansion entropy. For instance, this is the case for the initial condition $\mathbf{x}_0 = (-0.25, 0.3, 0)$ as can be seen with the green dotted curve in FIG.10.

### C. Chaotic mixing analysis

The Poincaré sections, Lyapunov exponents and expansion entropies studied in the previous sections show that the dynamics of many fluid points is chaotic. The system is therefore expected to mix efficiently. To study the efficiency of mixing in a practical manner, we have computed contamination rates, homogeneities and mixing times. These quantities are widely used in chemical engineering or biology, and are presented in the following lines.

First, we simulate the dispersion of clouds of $N_p = 5000$ particles initially located within a small cube (width $a/100$) near the center of the cube, by solving Eq. (34). The final time of the simulation of the cloud of particles is taken to be of a few convective times

$$t_c = \frac{a}{U} \qquad (46)$$

where $U \approx 0.1$ is an order of magnitude of the fluid velocity $\mathbf{v}$ in the MHD computations. This corresponds to $t_c = 10$ non-dimensional units (as a reminder, $a = 1$ in our case study). Note that, as explained in the previous section, the contribution of each magnet has been parameterized by means of the parameter $\alpha$ (Eq. (33)): $\alpha = 1$ corresponds to stirring by the side magnet alone (magnet #1, blue colors in the following figures), whereas $\alpha = 0$ corresponds to the central magnet alone (magnet #2, red).
<figure><img src="image_18.png" /><figcaption>FIG. 11. Expansion entropy of the entire flow versus $\alpha$.</figcaption></figure>

FIG. 12 shows these clouds at $t = 250 t_c$, for various values of $\alpha$. We observe that the shape of the cloud highly depends on $\alpha$ when $\alpha$ approaches 0 or 1, and mixing appears very poor there. Differences between transport by the approximately poloidal flow ($v_2$ alone, $\alpha = 0$) and the toroidal flow ($v_1$ alone, $\alpha = 1$) are clearly visible. The poloidal flow alone highly deforms the cloud of particles within a quasi-planar domain. In contrast, the sole toroidal flow has little effect on the initial spot of particles, which rotates near the box center without being stretched.

In a second time, we divide the flow domain into $N_c = 125,000$ ($= 50^3$) cubic sub-domains ("cells"), then simulate the dispersion of a much bigger cloud of $N_p = 1,000,000$ particles initially located, as above, within one of these cells near the center of the cube. We then derive two quantities of great interest for practical applications: the contamination rate $C(t)$ and the final homogeneity $H_\infty$. The former is the proportion of cells which were visited by at least one particle during the time interval $[0, t]$. The latter is based on the final particle density $n_l$, that is the number of particles in cell $l \in [1, N_c]$ at the final time of the simulation (see [31] for example). At that time, we calculate the

<figure><img src="image_19.png" /><figcaption>FIG. 12. Clouds of 5000 particles after 250 convective times, for various values of $\alpha$. Particles were initially released in a cube with side 1/100 located at the center of the domain. The poloidal flow $v_2$ ($\alpha = 0$) stretches the cloud but keeps it in a quasi-planar domain. The toroidal flow alone $v_1$ ($\alpha = 1$) has very little effect on the tiny initial cloud.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 13. Upper graph, left: contamination rate $C(t)$ versus time for various values of $\alpha$ from 0 (blue) to 1 (red). Upper graph, right: final values of $C_{\infty}$, at $t = 500 = 50 t_c$. Lower graph: time of mixing (defined by $C(t_{mix}) = 80\%$) versus $\alpha$.</figcaption></figure>

standard deviation $\sigma$ of $n_l$, then the final homogeneity of the distribution defined as

$$H_{\infty} = 1 - \frac{\sigma}{\sigma_{\max}} \qquad (47)$$

where $\sigma_{max}$ is the maximum value of the standard deviation, corresponding to the case where all particles gather in a single cell: $\sigma_{max} = N_p\sqrt{N_c - 1}/N_c \approx N_p/\sqrt{N_c}$. Note that some authors choose slightly different definitions (e.g. $H_{\infty} = 1 - (\sigma/\sigma_{max})^2$) [32] [33]. In any case, poor mixing corresponds to $H_{\infty} = 0$, whereas efficient mixing corresponds to $H_{\infty} \approx 1$.

FIG. 12 suggests that there is a wide range of $\alpha$'s where mixing seems efficient. This is confirmed in FIG. 13, as the contamination rate remains small for the extremal values of $\alpha$ (single magnet), whereas $C(t)$ increases quickly for intermediate values of $\alpha$ (two magnets). FIG. 13 shows the growth of contamination rates $C(t)$ for various $\alpha$'s (upper left graph), as well as the final contamination rate $C_{\infty}$ (upper right graph). For $\alpha$ in the range [0.1, 0.7], mixing is efficient with a contamination rate $C_{\infty} \ge 80\%$. It is slightly more efficient when $0.4 < \alpha < 0.65$; for these values, the largest final contamination rate $C_{\infty} \ge 90\%$ is reached.

The time of mixing is also an important parameter to quantify the quality of the device, and is sensitive to the parameter $\alpha$. It is defined here as the time $t_{mix}$ such that $C(t_{mix}) = 0.8$ (i.e. 80% of the cells have been visited at least once). It is plotted on the lower graph of FIG. 13 which shows that mixing is faster for $0.1 < \alpha < 0.6$ with a mixing time $t_{mix} \le 300$ (i.e. $\le 30 t_c$). The minimal $t_{mix}$ value is reached for $\alpha = 0.14$, which corresponds to the maximum of the expansion entropy $H_0$ plotted in FIG. 11.

The spatial quality of final mixing is quantified by means of the homogeneity coefficient $H_{\infty}$ defined above. We observe that mixing is excellent for $\alpha$ between 0.02 and 0.4 (FIG. 14 left) with $H_{\infty}$ values close to 1, and is a bit
<figure><img src="image_21.png" /><figcaption>FIG. 14. Left: Final homogeneity $H_\infty$ versus $\alpha$, right: product of final contamination and homogeneity ($C_\infty \times H_\infty$) versus $\alpha$.</figcaption></figure>

poorer when $0.5 < \alpha < 0.9$, though it remains close to 80 % there. And for $\alpha > 0.95$, it is clearly insufficient with $H_\infty$ values under 50%.

Finally, by plotting the product of our two quantity, final contamination rate and homogeneity $C_\infty \times H_\infty$ as a function of $\alpha$, we can determine a range of good compromise in terms of contamination rate, homogeneity and time of mixing. According to FIG. 14 (right), mixing is particularly efficient for $\alpha \in [0.08, 0.5]$. These values correspond perfectly with the range where the expansion entropy plotted in FIG. 11 is the highest, which confirms our results.

V. CONCLUDING REMARKS

The steady MHD device presented in this work has been shown to mix efficiently. It generates velocity fields which effectively produce a three-dimensional Stokes flow structure corresponding to a single recirculation cell superposed to an orthogonal double cell. Both elementary flows are created by means of magnets placed in an appropriate manner: the first magnet is located in an asymmetric manner and pushes the flow along the side of the tank, thus creating a large recirculation cell. The second magnet is placed along a mid-plane of the tank, and creates a double recirculation cell. Due to the linearity of the Stokes equation, the flow resulting from the action of both magnets is a linear combination of these elementary flows.

A weak vorticity – streamfunction formulation, based on the mixed finite element method, has been developed and implemented to calculate these velocity fields. Then, by computing Poincaré sections, Lyapunov exponents and expansion entropies, it was shown that most combinations of these two velocity fields produce fluid points that are indeed chaotic. We also observed that, for a wide range of $\alpha$, mixing appears to be particularly efficient in terms of contamination and homogeneity rates. Following this theoretical study, the authors recommend using a value of $\alpha$ between 0.1 and 0.5 for best efficiency.

Lyapunov exponents, as well as expansion entropies, suggest that the structure of the flow is different according to whether the first magnet or the second magnet dominate (i.e. $\alpha$ is nearest to 1 or 0). A detailed analysis of the model flow of Toussaint *et al.* [12] shows that a Shilnikov chaos might be at work there, but the stagnation points involved in this Shilnikov chaos are different whether $\alpha$ is nearest to 1 or 0 [34, 35]. Clearly, this model flow differs from ours in that it has periodic boundary conditions and an explicit solution. However, not too close to the walls, the structure of the MHD flow shares some common features with the flow of [12]. This point will be analyzed in the near future.

The results obtained in the present analysis are based exclusively on simulations. A demonstrator will be set up soon, and an experimental bench will be created. This device will be designed with the help of the present study, and measurements will be performed to validate experimentally our numerical results.
# Appendix: Implementation of numerical fields computation

In this part, we describe the main steps involved in the numerical resolution of the systems described in section III. The first step consists in generating the meshing of the geometry. For this purpose, the free and open-source software Gmsh is used [36]. The mesh size in the whole domain is controlled by the number $n_e$ of elements on the edges of the cubic tank. We therefore impose a structured mesh in $\Omega_f$; and outside, the characteristic length of the elements is then gradually scaled up as the distance from the tank increases, until the boundary $\Gamma_\infty$ is reached. The mesh of the domain with the chosen control parameter value ($n_e = 15$) is displayed on FIG. 2, and some characteristics are provided in TABLE I.

<table><caption>TABLE I. Some statistics of the chosen mesh</caption><thead><tr><th>Domain</th><th>Nodes</th><th>Triangles</th><th>Tetrahedra</th></tr></thead><tbody><tr><td>Ω</td><td>18,350</td><td>6,866</td><td>107,674</td></tr><tr><td>Ω<sub>f</sub></td><td>3,423</td><td>3,256</td><td>15,964</td></tr></tbody></table>

Afterwards, the formulations detailed in the previous section were then implemented in C++ language, using the free and open-source GmshFEM library [37, 38]. Thanks to this library, it is relatively easy to define the approximate subspaces based on our mesh of the global needed function spaces $H_0(\text{grad}, \Omega)$, $H(\text{grad}, \Omega_f)$, $H_0(\text{grad}, \Omega_f)$, $H(\text{curl}, \Omega_f)$, and $H_0(\text{curl}, \Omega_f)$. These finite dimensional function spaces are respectively called $W_0^0(\Omega)$, $W^0(\Omega_f)$, $W_0^0(\Omega_f)$, $W^1(\Omega_f)$ and $W_0^1(\Omega_f)$, and are defined thanks to hierarchical basis functions. The three first spaces are generated using 0-form nodal basis functions which are the classical nodal basis functions of Lagrange elements, their order are noted $p_n$. The two other spaces $W^1(\Omega_f)$ and $W_0^1(\Omega_f)$ are generated by 1-form edge basis functions of order $p_e$ built with the same functions as the previous ones, see [14, 15, 39] for details.

The numerical implementation and runs were performed on a dedicated workstation with 48 threads @3.5 GHz and 192 GiB of RAM. For the best possible accuracy of results, we choose to use high order elements with an interpolation order of $p_n = 4$ for the nodal basis functions of $W_0^0(\Omega)$, $W^0(\Omega_f)$ and $W_0^0(\Omega_f)$, and $p_e = 3$ for the edge basis functions of $W^1(\Omega_f)$ and $W_0^1(\Omega_f)$. These values combined with a mesh parameterized with the chosen value of $n_e$ allow us to make the best use of all our computing resources. The synthesis of the systems we solve (using the MUMPS direct solver [40]) is given by (A.1), and some details upon the numerical runs are reported in TABLE II.

<table><caption>TABLE II. Details of the numerical resolution ($n_e = 15, p_n = 4, p_e = 3$)</caption><thead><tr><th>System</th><th>(Σ<sub>hm,i</sub>)</th><th>(Σ<sup>(1)</sup><sub>hf,i</sub>)</th><th>(Σ<sup>(2)</sup><sub>hf,i</sub>)</th></tr></thead><tbody><tr><td>DOFs</td><td>1,207,089</td><td>527,823</td><td>1,198,914</td></tr><tr><td>Assembling time</td><td>2.47 sec</td><td>1.92 sec</td><td>4.12 sec</td></tr><tr><td>Solving time</td><td>208 sec</td><td>92 sec</td><td>3050 sec</td></tr></tbody></table>

The very low value of assembling time highlights the efficiency of the chosen solver. Moreover, even if the computation time of the last system ($\Sigma_{hf,i}^{(2)}$) seems high compared to the others (due to a non-symmetry of the stiffness matrix), the duration of the whole resolution remains largely acceptable for such a quite big problem.
$$
\begin{equation}
\begin{aligned}
& (\Sigma_{hm,i}) \left\{
\begin{array}{@{}l@{}}
\text{Find } \phi_{hi} \in W_0^0(\Omega) \text{ such that:} \\
\displaystyle \int_{\Omega_i} \mathbf{grad} \phi_{hi} \cdot \mathbf{grad} \xi_{hi} \,dx - \int_{\Omega_{mi}} \mathbf{m}_i \cdot \mathbf{grad} \xi_{hi} \,dx \\
\displaystyle \quad + \int_{\Omega_\infty} \mathbf{J}_F^{-1} \mathbf{grad} \phi_{hi} \cdot \mathbf{J}_F^{-1} \mathbf{grad} \xi_{hi} \,|\det \mathbf{J}_F| \,d\mathbf{X} = 0, \quad \forall \xi_{hi} \in W_0^0(\Omega)
\end{array}
\right\} \\
& \downarrow \phi_{hi} \\
& (\Sigma_{hf,i}^{(1)}) \left\{
\begin{array}{@{}l@{}}
\text{Find } (\omega_{hi}^0, \lambda_{hi}^0) \in W_0^1(\Omega_f) \times W^0(\Omega_f) \text{ such that:} \\
\displaystyle \int_{\Omega_f} \mathbf{curl} \omega_{hi}^0 \cdot \mathbf{curl} \varphi_{hi} \,dx + \int_{\Omega_f} \mathbf{grad} \lambda_{hi}^0 \cdot \varphi_{hi} \,dx \\
\displaystyle \quad + \int_{\Omega_f} (\mathbf{j}_0 \times \mathbf{grad} \phi_{hi}) \cdot \mathbf{curl} \varphi_{hi} \,dx = 0, \quad \forall \varphi_{hi} \in W_0^1(\Omega_f) \\
\displaystyle \int_{\Omega_f} \omega_{hi}^0 \cdot \mathbf{grad} \mu_{hi} \,dx = 0, \quad \forall \mu_{hi} \in W_0^0(\Omega_f)
\end{array}
\right\} \tag{A.1} \\
& \downarrow \omega_{hi}^0 \\
& (\Sigma_{hf,i}^{(2)}) \left\{
\begin{array}{@{}l@{}}
\text{Find } (\omega_{hi}^*, \lambda_{hi}^*, \psi_{hi}, \eta_{hi}) \in W^1(\Omega_f) \times W^0(\Omega_f) \times W_0^1(\Omega_f) \times W^0(\Omega_f) \text{ such that:} \\
\displaystyle \int_{\Omega_f} \omega_{hi}^* \cdot \varphi_{hi} \,dx - \int_{\Omega_f} \mathbf{curl} \psi_{hi} \cdot \mathbf{curl} \varphi_{hi} \,dx + \int_{\Omega_f} \mathbf{grad} \eta_{hi} \cdot \varphi_{hi} \,dx = \\
\displaystyle \quad - \int_{\Omega_f} \omega_{hi}^0 \cdot \varphi_{hi} \,dx, \quad \forall \varphi_{hi} \in W^1(\Omega_f) \\
\displaystyle \int_{\Omega_f} \mathbf{curl} \omega_{hi}^* \cdot \mathbf{curl} \xi_{hi} \,dx + \int_{\Omega_f} \mathbf{grad} \lambda_{hi}^* \cdot \xi_{hi} \,dx = 0, \quad \forall \xi_{hi} \in W_0^1(\Omega_f) \\
\displaystyle \int_{\Omega_f} \psi_{hi} \cdot \mathbf{grad} \chi_{hi} \,dx = 0, \quad \forall \chi_{hi} \in W_0^0(\Omega_f) \\
\displaystyle \int_{\Omega_f} (\omega_{hi}^0 + \omega_{hi}^*) \cdot \mathbf{grad} \mu_{hi} \,dx = 0, \quad \forall \mu_{hi} \in W^0(\Omega_f)
\end{array}
\right\}
\end{aligned}
\end{equation}
$$

ACKNOWLEDGMENTS

All the numerical tools used in this article are free and open-source software. Then, the authors would like to express their gratitude to the developers and contributors of the corresponding projects for providing and sharing their work to the community. We especially thank: Christophe Geuzaine and Jean-François Remacle for the mesh generator Gmsh; Anthony Royer, Eric Béchet and Christophe Geuzaine again for the GmshFEM library; and George Datseris for the DynamicalSystems.jl Julia package.

[1] J. M. Ottino, *The Kinematics of Mixing: Stretching, Chaos and Transport* (Cambridge University Press, Cambridge, 1989).

[2] K. El Omari, E. Younes, T. Burghelea, C. Castelain, Y. Moguen, and Y. Le Guer, Active chaotic mixing in a channel with rotating arc-walls, *Phys. Rev. Fluids* **6**, 024502 (2021).

[3] H. Aref and S. Balachandar, Chaotic advection in a Stokes flow, *Phys. Fluids* **29**, 3515 (1986).

[4] A. Acrivos, IUTAM International Symposium on Fluid Mechanics of Stirring and Mixing, University of California, La Jolla, Aug. 20-24, 1990, *Proceedings, Phys. Fluids A* **3** (1991).

[5] T. J. Kaper and S. Wiggins, An analytical study of transport in stokes flows exhibiting large-scale chaos in the eccentric journal bearing, *J. Fluid Mech.* **253**, 211 (1993).

[6] C. W. Leong and J. M. Ottino, Experiments on mixing due to chaotic advection in a cavity, *J. Fluid Mech.* **209**, 463 (1989).

[7] A. Figueroa, S. Cuevas, and E. Ramos, Lissajous trajectories in electromagnetically driven vortices, *Journal of Fluid Mechanics* **815**, 415–434 (2017).

[8] T. Dombre, U. Frisch, J. M. Greene, M. Hénon, A. Mehr, and A. M. Soward, Chaotic streamlines in the ABC flows, *Journal of Fluid Mechanics* **167**, 353 (1986).

[9] J. C. Leprévost, J. R. Angilella, and J. P. Brancher, Geometrical analysis of chaotic mixing in a low reynolds number magnetohydrodynamic quadripolar flow, *Phys. Rev. E* **63**, 056309 (2001).

[10] J. P. Brancher and J. C. Leprevost, Laminar mixing induced by a twisted quadripolar stokes flow, *Int. J. Heat Mass Transf.* **47**, 1561 (2004).
[11] V. Toussaint, *Etude des phénomènes de mélange dans des écoulements laminaires chaotiques en configuration tridimensionnelle*, Ph.D. thesis, Ecole centrale de Lyon (1997).

[12] V. Toussaint, P. Carrière, and F. Raynal, A numerical Eulerian approach to mixing by chaotic advection, Phys. Fluids **7**, 2587 (1995).

[13] J. C. Nédélec, Mixed finite elements in $R^3$, Numer Math (Heidelb) **35**, 315 (1980).

[14] Z. Ren and N. Ida, High order differential form-based elements for the computation of electromagnetic field, IEEE Trans. Magn. **36**, 1472 (2000).

[15] C. Geuzaine, *High order hybrid finite element schemes for Maxwell's equations taking thin structures and global quantities into account*, Ph.D. thesis, Université de Liège (2001).

[16] J. C. Nédélec, A new family of mixed finite elements in $R^3$, Numer Math (Heidelb) **50**, 57 (1986).

[17] G. Akoun and J.-P. Yonnet, 3D analytical calculation of the forces exerted between two cuboidal magnets, IEEE Transactions on Magnetics **20**, 1962 (1984).

[18] E. P. Furlani, Chapter 4 - permanent magnet applications, in *Permanent Magnet and Electromechanical Devices*, Electro-magnetism, edited by E. P. Furlani (Academic Press, San Diego, 2001) pp. 207-333.

[19] J. F. Imhoff, G. Meunier, X. Brunotte, and J. C. Sabonnadière, An original solution for unbounded electromagnetic 2D- and 3D-problems throughout the finite element method, IEEE Trans. Magn. **26**, 1659 (1990).

[20] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault, Vector potentials in three-dimensional non-smooth domains, Math. Methods Appl. Sci. **21**, 823 (1998).

[21] M. Amara, H. Barucq, and M. Duloué, Une formulation mixte convergente pour le système de stokes tridimensionnel, C. R. Acad. Sci. - Series I - Mathematics **328**, 935 (1999).

[22] M. Amara and C. Bernardi, Convergence of a finite element discretization of the navier-stokes equations in vorticity and stream function formulation, Esaim Math Model Numer Anal **33**, 1033 (1999).

[23] G. Datseris, DynamicalSystems.jl: a Julia software library for chaos and nonlinear dynamics, JOSS **3**, 598 (2018).

[24] G. Datseris and U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code* (Springer Nature, Cham, Switzerland, 2022).

[25] J. P. Eckmann and D. Ruelle, Ergodic theory of chaos and strange attractors, Rev. Mod. Phys. **57**, 617 (1985).

[26] A. Pikovsky and A. Politi, *Lyapunov Exponents: A Tool to Explore Complex Dynamics* (Cambridge University Press, 2016).

[27] R. L. Adler, A. G. Konheim, and M. H. McAndrew, Topological entropy, Transactions of the American Mathematical Society **114**, 309 (1965).

[28] B. R. Hunt and E. Ott, Defining chaos, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **25**, 097618 (2015).

[29] O. S. Kozlovski, An integral formula for topological entropy of $c^\infty$ maps, Ergodic Theory and Dynamical Systems **18**, 405-424 (1998).

[30] A. Katok and B. Hasselblatt, *Introduction to the Modern Theory of Dynamical Systems*, Encyclopedia of Mathematics and its Applications (Cambridge University Press, 1995).

[31] M. Engler, N. Kockmann, T. Kiefer, and P. Woias, Numerical and experimental investigations on liquid mixing in static micromixers, Chemical Engineering Journal **101**, 315 (2004).

[32] P. Danckwerts, The definition and measurement of some characteristics of mixtures, Applied Scientific Research, Section A **3**, 279 (1952).

[33] H. Berthiaux, V. Mosorov, L. Tomczak, C. Gatumel, and J.-F. Demeyre, Principal component analysis for characterising homogeneity in powder mixing using image processing techniques, Chemical Engineering and Processing: Process Intensification **45**, 397 (2006).

[34] D. D. Holm and Y. Kimura, Zero-helicity lagrangian kinematics of three-dimensional advection, Physics of Fluids A: Fluid Dynamics **3**, 1033 (1991).

[35] S. Gonchenko, A. Kazakov, D. Turaev, and A. L. Shilnikov, Leonid shilnikov and mathematical theory of dynamical chaos, *Chaos: An Interdisciplinary Journal of Nonlinear Science* **32** (2022).

[36] C. Geuzaine and J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int J Numer Methods Eng **79**, 1309 (2009).

[37] A. Royer, E. Béchet, and C. Geuzaine, Gmsh-Fem: An efficient finite element library based on Gmsh, in *14th World Congress on Computational Mechanics (WCCM), ECCOMAS Congress 2020* (2021) pp. 1-13.

[38] A. Royer, *Efficient finite element methods for solving high-frequency time-harmonic acoustic wave problems in heterogeneous media*, Ph.D. thesis, ULiège - Université de Liège [Faculté des Sciences Appliquées], Belgique, Belgium (2023).

[39] P. Solin, K. Segeth, and I. Dolezel, *Higher-Order Finite Element Methods*, 1st ed. (Chapman and Hall/CRC, 2003).

[40] P. Amestoy, A. Buttari, J.-Y. L'Excellent, and T. Mary, Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures, ACM Trans Math Softw **45**, 2:1 (2019).