# Mathematical Derivation of a Homogenized Bidomain Model for Pulsed Field Ablation

Simon Bihoreau-Duchemin, Annabelle Collin, Michael Leguèbe, Clair

Poignard

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

## To cite this version:

Simon Bihoreau-Duchemin, Annabelle Collin, Michael Leguèbe, Clair Poignard. Mathematical Derivation of a Homogenized Bidomain Model for Pulsed Field Ablation. 2025. hal-05368825

HAL Id: hal-05368825

https://inria.hal.science/hal-05368825v1

Preprint submitted on 17 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical Derivation of a Homogenized Bidomain Model for Pulsed Field Ablation

Simon Bihoreau Duchemin<sup>1,2</sup>, Annabelle Collin<sup>3</sup>, Michael Leguèbe<sup>1,2,4</sup>, and Clair Poignard<sup>5</sup>

<sup>1</sup>Centre Inria de l'Université de Bordeaux, Talence, F-33400, France

<sup>2</sup>Institut de Mathématiques de Bordeaux, Université de Bordeaux, Talence, F-33400, France

<sup>3</sup>Laboratoire de Mathématiques Jean Leray, Nantes Université, Nantes, F-44100, France

<sup>4</sup>IHU Liryc, Fondation Bordeaux Université, Pessac, F-33600, France

<sup>5</sup>Centre Inria de l'Université de Rennes, Rennes, F-35042, France

September 2025

## Abstract

We propose a static multiscale mathematical model for pulsed field ablation (PFA) in cardiac tissue, grounded in a detailed microscopic description of intra- and extracellular electrical activity. Using high-order two-scale asymptotic expansions, we derive a physically meaningful macroscopic bidomain system that captures the nonlinear effects of electroporation through membrane conductivity which depends on homogenized transmembrane voltage. A key feature of our model is the presence of the small parameter $\varepsilon$ within the nonlinear membrane law, which leads to a fundamentally different macroscopic behavior compared to classical electrophysiological models. We then conduct a numerical study to investigate the convergence of the model with respect to $\varepsilon$ and to illustrate qualitatively the behavior of the homogenized system under biologically relevant stimulation protocols. Finally, we prove the convergence of the asymptotic expansion to arbitrary order, providing a full mathematical justification of the model.

## 1 Introduction

Cardiac ablation is a widespread therapy for the treatment of cardiac arrhythmias. It consists in the destruction of cardiac cells in targeted areas responsible for the disruption of the electrophysiological current pattern which drives heart contraction. A thin and flexible tube, called catheter, is inserted inside the heart through a blood vessel. The catheter contains a device, which depends on the technology used, that delivers energy into the region to ablate. The most common ablation technologies are Radio-Frequency Ablation (RFA) [43] relying on heat, cryoablation [10] relying on cold and Pulsed Field Ablation (PFA) relying on a phenomenon called irreversible electroporation.

Recent research focuses on PFA, a promising technology that drastically reduces the risks of damaging adjacent structures as lungs during the ablation while improving the clinical performance [35, 13, 50, 59]. PFA consists in delivering short and intense electric pulses in the cardiac tissue, which induces cell death through irreversible electroporation. Under the effect of the intense electric field, pores appear on the cell membranes which become permeable and
are destroyed as a result of excessive exchanges between intra- and extracellular media. Electroporation has been in use for oncology for approximately 30 years, and its use in cardiology is now increasing [53].

Given the complexity of the electroporation phenomenon and of the cardiac tissue, modeling electrical conduction at tissue scale during PFA - which is the central objective of our work - remains a major challenge for improving both its understanding and its clinical applicability. Several models have already been developed for electrically passive tissues. In most cases, an electrostatic model is used to compute a single electric potential $u$ which satisfies $-\nabla \cdot \sigma \nabla u = f$, where electroporation is modeled through a conductivity $\sigma$ that depends nonlinearly on the electric field $\|\nabla u\|$ [41, 30, 32]. This approach, that we call *standard* in the following, is phenomenological: the non-linearity of $\sigma$ is not based on a detailed description of the physics on the cell membranes.

Several variations of the standard model for electrically passive tissues have been proposed in the literature: either with a dynamical response of the conductivity with the intensity of the electric field [57], or, more commonly, by adding a dependence on temperature for irreversible electroporation (see for instance Supp. Material 3 of [1] for a collection of standard model variations). Concerning cardiac tissue specifically, some studies use directly the standard model as for passive tissue [31, 40, 48], while others focus on the anisotropic property of muscle fibers [14, 36]. Another approach was proposed earlier in [2], where the authors propose a two dimensional bidomain model - which models cardiac electrical activity and results from the periodic homogenization of microscopic bidomain model accounting for the ionic exchanges at cell membranes - with an electric conductivity that depends solely on the transmembrane voltage instead of $\|\nabla u\|$. Simply plugging an electroporation model into cardiac bidomain equations may not be physically justified, for at least two reasons. First, there is a large magnitude difference between the physiological ionic currents and PFA pulses, which should not lead to the classical bidomain equations after the homogenization process. Secondly, in the case of two electrodes at opposite potentials of the same intensity, the transmembrane potential in the middle of the electrodes is zero while the electroporated region clearly includes this point [19].

Here, we derive a physically inspired model specifically for PFA, consistent with a microscopic description that explicitly includes intra- and extracellular media. While our approach is conceptually related to the cardiac bidomain model, our microscopic problem is different and tailored to capture the effects of electroporation. Starting from this microscopic PFA problem, we rigorously carry out the homogenization process with appropriate scaling of the physical quantities involved. Our goal is to derive a macroscopic bidomain model for PFA that is both mathematically justified and physically interpretable, and that accurately captures the effects of electroporation.

The microscopic model from which we start is itself of great interest, as it provides a detailed description of the physical processes at the cellular scale. However, due to the extremely high computational cost of microscopic model - see for instance [15] in a cardiac electrophysiology context - it is impractical for large-scale simulations and would be completely intractable in any realistic clinical setting. Deriving a macroscopic homogenized model is therefore essential for enabling clinically relevant numerical studies.

From a mathematical standpoint, several homogenization techniques can be used to derive a macroscopic model from a microscopic one. In the context of cardiac electrophysiology, notable approaches include periodic unfolding methods [12], $\Gamma$-convergence [46], and two-scale convergence as introduced by Nguetseng [45] and Allaire [5] and applied in various contexts, including cardiac models [17]. In this work, we will carry out the expansion to high orders, since the classical zero-order limit, obtained by two-scale convergence, is not sufficient in our setting.
While such higher-order expansions have not yet been developed for cardiac bidomain models to our knowledge, they have been studied in other contexts; see for instance [26, 25, 56, 7, 16, 51].

In Section 2, we introduce the model which accounts for the cell scale, and its nondimensionalization in the context of PFA. We then provide a mathematical analysis including the proof of well-posedness of the equations and the derivation of energy estimates with respect to the small parameter $\epsilon$ corresponding to the ratio between the cellular and tissue scales. In Section 3, we study the limit problem as $\epsilon$ tends to zero using two-scale convergence, and formally derive the equation satisfied by the next-order asymptotic term. Together, these constitute the proposed tissue-scale model for PFA, which is the main result of this work and is summarized in Section 3.4 along with an effective model for means of comparison with the classical bidomain for electrophysiology. Our homogenized model and then its effective version (24) differs from the bidomain model by involving both the electric field magnitude and the macroscopic transmembrane potential in the nonlinearity, while at the microscopic scale only the transmembrane potential matters. This shows that electroporation is driven by transmembrane potential at the cell scale but by the electric field at the tissue scale. We then present numerical evidence of the convergence of the cell-scale model towards this macroscopic model, and introduce a quantity to assess electroporation at the tissue scale, supported by numerical illustrations. Encouraged by the observed numerical convergence rates, we provide a rigorous justification of the asymptotic expansion to any order $n$ in Section 4.

# 2 Microscopic modeling and analysis

This section is devoted to the microscopic model considered in this work. We first introduce the geometrical setup, then present the model itself and its nondimensionalization, and finally establish existence and uniqueness results along with energy estimates.

## 2.1 Geometrical setup

Let $d = 2$ or 3. Let $Y = [0, 1]^d$ be the unit mathematical cell with boundary $\partial Y$. Consider two subdomains $Y_i$ and $Y_e$ of $Y$, representing the intra- and extracellular media, such that $\bar{Y}_i \cup \bar{Y}_e = Y$ and $\hat{Y}_i \cap \hat{Y}_e = \emptyset$. The interface between these subdomains, corresponding to the cardiac cell membrane, is noted $\Gamma = \partial \bar{Y}_i \cap \partial \bar{Y}_e$ (see Figure 1).

To represent the tissue, let $\Omega := \prod_{i=1}^d ]a_i, b_i[ \subset \mathbb{R}^d$, composed of multiple periodically arranged cardiac cells. We assume that there exists a parameter $\epsilon > 0$ and an index set $N_\epsilon \subset \mathbb{Z}^d$ such that

$$ \bar{\Omega} = \bigcup_{n \in N_{\epsilon}} \epsilon(Y + n). $$

It is worth noting that $\epsilon$ represents the ratio between the characteristic length scale of the cells and that of the macroscopic tissue domain.

We then define the intra- and extracellular subdomains of $\Omega$, and their interface by

$$ \bar{\Omega}_i^\epsilon := \bigcup_{n \in N_\epsilon} \epsilon(\bar{Y}_i + n), \quad \bar{\Omega}_e^\epsilon := \bigcup_{n \in N_\epsilon} \epsilon(\bar{Y}_e + n), \quad \Gamma_m^\epsilon := \bigcup_{n \in N_\epsilon} \epsilon(\Gamma + n), $$

so that they are unions of rescaled and translated unit-cell subdomains and interfaces. We assume that $\Omega_e^\epsilon$ is connected, which is necessary for the Poincaré-Wirtinger inequality in Theorem A.1 to hold, and that $\Gamma_m^\epsilon$ is infinitely smooth.
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of unit cell Y (left) and macroscopic domain Ω (right)</figcaption></figure>

## 2.2 Microscopic bidomain equations

### 2.2.1 Model statement

We consider the following microscopic bidomain model in which the intra- and extra-cellular electric potentials $u_i^\epsilon$ and $u_e^\epsilon$ verify

$$-\nabla \cdot (\sigma_i^\epsilon \nabla u_i^\epsilon) = f^\epsilon|_{\Omega_i^\epsilon}, \qquad \text{in } \Omega_i^\epsilon, \qquad (1a)$$

$$-\nabla \cdot (\sigma_e^\epsilon \nabla u_e^\epsilon) = f^\epsilon|_{\Omega_e^\epsilon}, \qquad \text{in } \Omega_e^\epsilon, \qquad (1b)$$

$$\sigma_i^\epsilon \nabla u_i^\epsilon \cdot n_i = \sigma_e^\epsilon \nabla u_e^\epsilon \cdot n_i, \qquad \text{on } \Gamma_m^\epsilon, \qquad (1c)$$

$$-\sigma_i^\epsilon \nabla u_i^\epsilon \cdot n_i = I_m^\epsilon (u_i^\epsilon - u_e^\epsilon), \qquad \text{on } \Gamma_m^\epsilon, \qquad (1d)$$

where $\sigma_i^\epsilon > 0$ and $\sigma_e^\epsilon > 0$ are the conductivities of the intra- and extra-cellular media (in $S \text{ m}^{-1}$) and $f^\epsilon$ (in $\text{A m}^{-3}$) is a source term representing the electric pulses. System (1) is complemented with a gauge condition $\int_{\Omega_e^\epsilon} u_e^\epsilon = 0$, and with periodic boundary conditions.

In cardiac electrophysiology, the transmembrane current $I_m^\epsilon$ originates mainly from complex ionic exchanges across the cell membrane, driven by the transmembrane potential $v_m^\epsilon := u_i^\epsilon - u_e^\epsilon$ and gating variables, and described by physiological [29, 54] or phenomenological [3, 42] models. However, in the context of pulsed electric fields, the high intensity during the applied short pulses justify neglecting the physiological ionic currents [18, 38, 34]. Instead, the dominant contribution to $I_m^\epsilon$ stems from a resistive current induced by the sharp increase in membrane permeability during electroporation. In this work, we adopt the modeling approach proposed in [38, 34], where the transmembrane current is simplified as a purely resistive flux across the membrane. The membrane conductivity is then modeled by a nonlinear function $S_m^\epsilon$, which depends on the transmembrane voltage:

$$I_m^\epsilon(v) := S_m^\epsilon(v) v. \qquad (2)$$

At rest, the cell membranes have a high resistive effect, which is reduced by electroporation. As consequence, $S_m^\epsilon$, which is the inverse of the surface resistance, always take values that are strictly positive, although very low when the cell is at rest.

### 2.2.2 Nondimensionalization

In order to highlight the relative amplitudes of the different quantities that are involved, in the specific use of electroporation ablation on cardiac tissue, we proceed by nondimensionalizing System (1). A piece of tissue considered in cardiac ablation is smaller than the whole heart (around $10^{-1}$ m), typically its characteristic size is $L_0 = 10^{-2}$ m [49]. The characteristic cardiac
cell size is $10^{-4}$ m [52], therefore the ratio between the two scales at stake is $\epsilon = 10^{-2}$. An electrical pulse for PFA leads to a potential that is around $U_0 = 10^2$ V [49]. At the cell level, electroporation occurs after the transmembrane potential reaches a threshold of characteristic value $V_0^{ep} = 1$ V [55]. A significant increase in membrane conductance is then observed, with the baseline value being on the order of $1$ S m$^{-2}$. According to [30, 37], the exact maximum value of this conductance has little influence on the overall dynamics. We therefore choose to model the characteristic electroporated membrane conductance $S_0$ as substantially above the baseline, and set $S_0 = 10^2$ S m$^{-2}$. As in every cell scale electroporation model, we assume that the increase of $S_m^\epsilon$ is located around $V_0^{ep}$. In consequence, $S_m^\epsilon$ can be rewritten in the electroporation regime as $S_m^\epsilon(\lambda V_0^{ep})$, where $\lambda$ is a unitless variable of order 1. The physiological value of the intra- and extra-cellular conductivities is around $\Sigma_0 = 1$ S m$^{-1}$ [52]. The injected current is of order 1 A ($= \Sigma_0 L_0 U_0$) so the volume current source $f^\epsilon$ has a characteristic value $F_0 = \frac{1}{L_0^3}$ A = $10^6$ A m$^{-3}$.

Using these normalization constants, we define the following unitless quantities

$$ \tilde{u}_{i,e}^{\epsilon} \left( \frac{\cdot}{L_0} \right) = \frac{1}{U_0} u_{i,e}^{\epsilon}(\cdot), $$

$$ \tilde{\sigma}_{i,e} = \frac{1}{\Sigma_0} \sigma_{i,e}^{\epsilon}, $$

$$ \tilde{S}_m \left( \frac{\cdot}{V_0^{ep}} \right) = \frac{1}{S_0} S_m^{\epsilon}(\cdot), $$

$$ \tilde{f} \left( \frac{\cdot}{L_0} \right) = \frac{1}{F_0} f^{\epsilon}(\cdot), $$

where we assumed that the dependence relatively to $\epsilon$ only lies implicitly in the characteristic quantities. By substituting these quantities in System (1), we obtain

$$ -\tilde{\nabla} \cdot (\tilde{\sigma}_i \tilde{\nabla} \tilde{u}_i^\epsilon) = \frac{L_0^2 F_0}{\Sigma_0 U_0} \tilde{f}|_{\tilde{\Omega}_i^\epsilon}, \qquad \tilde{\Omega}_i^\epsilon, \qquad (3a) $$

$$ -\tilde{\nabla} \cdot (\tilde{\sigma}_e \tilde{\nabla} \tilde{u}_e^\epsilon) = \frac{L_0^2 F_0}{\Sigma_0 U_0} \tilde{f}|_{\tilde{\Omega}_e^\epsilon}, \qquad \tilde{\Omega}_e^\epsilon, \qquad (3b) $$

$$ \tilde{\sigma}_i \tilde{\nabla} \tilde{u}_i^\epsilon \cdot n_i = \tilde{\sigma}_e \tilde{\nabla} \tilde{u}_e^\epsilon \cdot n_i, \qquad \tilde{\Gamma}_m^\epsilon, \qquad (3c) $$

$$ -\tilde{\sigma}_i \tilde{\nabla} \tilde{u}_i^\epsilon \cdot n_i = \frac{S_0 L_0}{\Sigma_0} \tilde{S}_m \left( \frac{U_0}{V_0^{ep}} (\tilde{u}_i^\epsilon - \tilde{u}_e^\epsilon) \right) (\tilde{u}_i^\epsilon - \tilde{u}_e^\epsilon), \qquad \tilde{\Gamma}_m^\epsilon, \qquad (3d) $$

where the domains $\tilde{\Omega}_{i,e}^\epsilon$ and $\tilde{\Gamma}_m^\epsilon$ are obtained by rescaling $\Omega_{i,e}^\epsilon$ and $\Gamma_m^\epsilon$ by $L_0$. This leads to the following nondimensionalized system by omitting the tilde notation for better readability and denoting $I_m(v) := S_m(v)v$,

$$ -\nabla \cdot (\sigma_i \nabla u_i^\epsilon) = f|_{\Omega_i^\epsilon}, \qquad \Omega_i^\epsilon, \qquad (4a) $$

$$ -\nabla \cdot (\sigma_e \nabla u_e^\epsilon) = f|_{\Omega_e^\epsilon}, \qquad \Omega_e^\epsilon, \qquad (4b) $$

$$ \sigma_i \nabla u_i^\epsilon \cdot n_i = \sigma_e \nabla u_e^\epsilon \cdot n_i, \qquad \Gamma_m^\epsilon, \qquad (4c) $$

$$ -\sigma_i \nabla u_i^\epsilon \cdot n_i = \varepsilon I_m \left( \frac{u_i^\epsilon - u_e^\epsilon}{\varepsilon} \right), \qquad \Gamma_m^\epsilon, \qquad (4d) $$

which is the starting point of the homogenization process.

*Remark 2.1.* The main difference with the microscopic bidomain model traditionally used in cardiac electrophysiology lies in the last equation, which, in a static framework and after nondimensionalization (see [17]), typically reads

$$ -\sigma_i \nabla u_i^\epsilon \cdot n_i = \varepsilon I_{\text{ion}}(u_i^\epsilon - u_e^\epsilon), \quad \Gamma_m^\epsilon, $$

where $I_{\text{ion}}$ denotes the ionic current model. The presence of $\epsilon$ inside the nonlinearity plays a crucial role and, as we will see later, leads to a different behavior.
## 2.3 Analysis of the microscopic problem

### 2.3.1 Periodic spaces definition.

We introduce several notations for periodic spaces. Let $D = \Omega$ or $Y$ and let $\mathcal{D}$ a subset of $D$, typically $\mathcal{D} = \Omega_i^\varepsilon, \Omega_e^\varepsilon$ or $\Omega$ if $D = \Omega$ and $\mathcal{D} = Y_i, Y_e, Y$ or $\Gamma$ if $D = Y$. We denote by

$$C_{\text{per}}^{\infty}(\mathcal{D}) = \{u|_{\mathcal{D}}, u \in C^{\infty}(\mathbb{R}^d), u \text{ is D-periodic}\},$$

and $H_{\text{per}}^1(\mathcal{D})$ its completion for the $H^1(\mathcal{D})$-norm. We also define, for $\mathcal{D}$ a subset of $Y$,

$$C_{\text{per}}^{\infty}(\Omega \times \mathcal{D}) = \{u|_{\Omega \times \mathcal{D}}, u \in C^{\infty}(\mathbb{R}^d \times \mathbb{R}^d), u \text{ is } \Omega \times Y\text{-periodic}\}.$$

For $E(\mathcal{D})$ being any of these spaces, we denote by

$$\bar{E}(\mathcal{D}) := \{u \in E(\mathcal{D}), \int_{\mathcal{D}} u = 0\}.$$

### 2.3.2 Existence, uniqueness and energy estimates.

In this section, we prove the well-posedness of the microscopic system and the energy estimates verified by the solution. The analysis is done under the following mathematical assumptions.

**Assumption 1** (Source term). $f \in \bar{C}_{\text{per}}^{\infty}(\Omega)$.

**Assumption 2** (Non linearity). $S_m \in C^{\infty}(\mathbb{R})$ is even on $\mathbb{R}$, non-decreasing on $\mathbb{R}^+$ and there exist $S_m^0, S_m^1, L$ and $M$ positive constants such that

$$\begin{aligned}
& \forall \lambda \in \mathbb{R}, \quad 0 < S_m^0 \le S_m(\lambda) \le S_m^1, \\
& \forall (\lambda, \mu) \in \mathbb{R}^2, \quad |S_m(\lambda) - S_m(\mu)| \le L|\lambda - \mu|, \\
& \forall \lambda \in \mathbb{R}, \quad |I_m'(\lambda)| = |(S_m(\cdot)\cdot)'(\lambda)| \le M.
\end{aligned}$$

*Remark 2.2.* We can show that Assumption 2 implies that for all $(\lambda, \mu) \in \mathbb{R}^2$,

$$(S_m(\lambda)\lambda - S_m(\mu)\mu)(\lambda - \mu) \ge S_m^0(\lambda - \mu)^2.$$

or equivalently,

$$(S_m(\cdot)\cdot)'(\lambda) = S_m'(\lambda)\lambda + S_m(\lambda) > S_m^0 > 0.$$

**Theorem 2.3.** Under Assumptions 1 and 2, System (4) has a unique solution $(u_i^\varepsilon, u_e^\varepsilon)$ in $C_{\text{per}}^{\infty}(\Omega_i^\varepsilon) \times C_{\text{per}}^{\infty}(\Omega_e^\varepsilon)$. Moreover, the following estimates hold for all $\varepsilon > 0$ small enough,

$$\|u_e^\varepsilon\|_{H^1(\Omega_e^\varepsilon)} \le C \|f\|_{L^2(\Omega)}$$

$$\|u_i^\varepsilon\|_{H^1(\Omega_i^\varepsilon)} \le C \|f\|_{L^2(\Omega)}$$

$$\|v_m^\varepsilon\|_{L^2(\Gamma_m^\varepsilon)} \le C \|f\|_{L^2(\Omega)},$$

where C does not depend on $\varepsilon$, and $v_m^\varepsilon := u_i^\varepsilon - u_e^\varepsilon$ in the sense of trace operators on the interface $\Gamma_m^\varepsilon$.
Proof. This proof relies on Theorems A.1 and A.2, originally established in [9] and recalled in Appendix A for completeness. We begin the proof by establishing existence and uniqueness. First, we show the existence of a solution to the variational formulation, then we prove its regularity in $C_{\text{per}}^{\infty}$, and finally we conclude with the proof of energy estimates. Existence and uniqueness of the variational formulation. We first show that the following variational formulation has a unique solution: find $(u_i^\varepsilon, u_e^\varepsilon) \in H_{\text{per}}^1(\Omega_i^\varepsilon) \times \overline{H}_{\text{per}}^1(\Omega_e^\varepsilon)$ such that for all $(\phi_i, \phi_e) \in H_{\text{per}}^1(\Omega_i^\varepsilon) \times \overline{H}_{\text{per}}^1(\Omega_e^\varepsilon)$,

$$ \int_{\Omega_i^\varepsilon} \sigma_i \nabla u_i^\varepsilon \cdot \nabla \phi_i + \int_{\Omega_e^\varepsilon} \sigma_e \nabla u_e^\varepsilon \cdot \nabla \phi_e + \int_{\Gamma_m^\varepsilon} \varepsilon I_m \left( \frac{u_i^\varepsilon - u_e^\varepsilon}{\varepsilon} \right) (\phi_i - \phi_e) = \int_{\Omega_i^\varepsilon} f \phi_i + \int_{\Omega_e^\varepsilon} f \phi_e. \quad (5) $$

Let us first prove the unicity. Assume that $(u_i, u_e)$ and $(U_i, U_e)$ are two solutions. The difference $(\nu_i, \nu_e) := (u_i - U_i, u_e - U_e)$ verifies

$$ \sigma_i \|\nabla \nu_i\|_{L^2(\Omega_i^\varepsilon)}^2 + \sigma_e \|\nabla \nu_e\|_{L^2(\Omega_e^\varepsilon)}^2 + \varepsilon \int_{\Gamma_m^\varepsilon} \left( I_m \left( \frac{u_i - u_e}{\varepsilon} \right) - I_m \left( \frac{U_i - U_e}{\varepsilon} \right) \right) (\nu_i - \nu_e) = 0. $$

Using Assumption 2 to bound the third term, we show that $\|\nabla \nu_i\|_{L^2(\Omega_i^\varepsilon)} = \|\nabla \nu_e\|_{L^2(\Omega_e^\varepsilon)} = \|\nu_i - \nu_e\|_{L^2(\Gamma_m^\varepsilon)} = 0$. Thanks to the gauge condition, it implies that $\nu_e = 0$ in $\Omega_e^\varepsilon$. Therefore $0 = \|\nu_i - \nu_e\|_{L^2(\Gamma_m^\varepsilon)} = \|\nu_i\|_{L^2(\Gamma_m^\varepsilon)}$, and $\nu_i = 0$ in $\Omega_i^\varepsilon$ because it is constant.

Let us now prove the existence of the solution. Consider the following functional defined for all $(v_i, v_e) \in H_{\text{per}}^1(\Omega_i^\varepsilon) \times \overline{H}_{\text{per}}^1(\Omega_e^\varepsilon)$ by

$$ \mathcal{F}(v_i, v_e) = \frac{1}{2} \int_{\Omega_i^\varepsilon} \sigma_i |\nabla v_i|^2 + \frac{1}{2} \int_{\Omega_e^\varepsilon} \sigma_e |\nabla v_e|^2 + \varepsilon^2 \int_{\Gamma_m^\varepsilon} \int_0^{\frac{(v_i - v_e)(x)}{\varepsilon}} I_m(s) \, ds \, dx - \int_{\Omega_i^\varepsilon} f v_i - \int_{\Omega_e^\varepsilon} f v_e. $$

$\mathcal{F}$ is $C^1$ and its derivative $\mathcal{F}'$ is

$$ \mathcal{F}'(v_i, v_e) \cdot (\phi_i, \phi_e) = \int_{\Omega_i^\varepsilon} \sigma_i \nabla v_i \cdot \nabla \phi_i + \int_{\Omega_e^\varepsilon} \sigma_e \nabla v_e \cdot \nabla \phi_e + \varepsilon \int_{\Gamma_m^\varepsilon} I_m \left( \frac{v_i - v_e}{\varepsilon} \right) (\phi_i - \phi_e) - \int_{\Omega_i^\varepsilon} f \phi_i - \int_{\Omega_e^\varepsilon} f \phi_e. $$

We show in the following that $\mathcal{F}$ is convex and coercive, thus it admits a critical point (see Corollary 1.4 of Chapter 3 of [33]), which is exactly a solution of the variational formulation.

We first show that $\mathcal{F}$ is convex. Let $(v_i, v_e)$ and $(w_i, w_e)$ in $H_{\text{per}}^1(\Omega_i^\varepsilon) \times \overline{H}_{\text{per}}^1(\Omega_e^\varepsilon)$, we have

$$ (\mathcal{F}'(v_i, v_e) - \mathcal{F}'(w_i, w_e)) \cdot ((v_i, v_e) - (w_i, w_e)) = \int_{\Omega_i^\varepsilon} \sigma_i |\nabla (v_i - w_i)|^2 + \int_{\Omega_e^\varepsilon} \sigma_e |\nabla (v_e - w_e)|^2 + \varepsilon \int_{\Gamma_m^\varepsilon} \left( I_m \left( \frac{v_i - v_e}{\varepsilon} \right) - I_m \left( \frac{w_i - w_e}{\varepsilon} \right) \right) (v_i - v_e - (w_i - w_e)) \ge 0, $$

because every term is positive, thanks to Assumption 2.
We then show the coercivity of $\mathcal{F}$. Let $(v_i, v_e) \in H_{\text{per}}^1(\Omega_i^\varepsilon) \times \overline{H_{\text{per}}^1(\Omega_e^\varepsilon)}$ and $\alpha > 0$, we get

$$
\begin{align*}
\mathcal{F}(v_i, v_e) &\ge \frac{1}{2} \left( \sigma_i \|\nabla v_i\|_{L^2(\Omega_i^\varepsilon)}^2 + \sigma_e \|\nabla v_e\|_{L^2(\Omega_e^\varepsilon)}^2 \right) + \frac{S_m^0}{2} \|v_i - v_e\|_{L^2(\Gamma_m^\varepsilon)}^2 \\
&\quad - \int_{\Omega_i^\varepsilon} f v_i - \int_{\Omega_e^\varepsilon} f v_e, \\
&\ge C \left( \|v_i\|_{H^1(\Omega_i^\varepsilon)}^2 + \|v_e\|_{H^1(\Omega_e^\varepsilon)}^2 \right) - \int_{\Omega_i^\varepsilon} f v_i - \int_{\Omega_e^\varepsilon} f v_e, \\
&\ge C \left( \|v_i\|_{H^1(\Omega_i^\varepsilon)}^2 + \|v_e\|_{H^1(\Omega_e^\varepsilon)}^2 \right) - \frac{\alpha}{2} \|v_i\|_{L^2(\Omega_i^\varepsilon)}^2 - \frac{\alpha}{2} \|v_e\|_{L^2(\Omega_e^\varepsilon)}^2 \\
&\quad - \frac{1}{2\alpha} \|f\|_{L^2(\Omega_i^\varepsilon)}^2 - \frac{1}{2\alpha} \|f\|_{L^2(\Omega_e^\varepsilon)}^2, \\
&\ge C \left( \|v_i\|_{H^1(\Omega_i^\varepsilon)}^2 + \|v_e\|_{H^1(\Omega_e^\varepsilon)}^2 \right) - C' \frac{\alpha}{2} \|v_i\|_{H^1(\Omega_i^\varepsilon)}^2 - C' \frac{\alpha}{2} \|v_e\|_{H^1(\Omega_e^\varepsilon)}^2 \\
&\quad - \frac{1}{2\alpha} \|f\|_{L^2(\Omega_i^\varepsilon)}^2 - \frac{1}{2\alpha} \|f\|_{L^2(\Omega_e^\varepsilon)}^2.
\end{align*}
$$

First line uses Assumption 2, second line uses norm equivalence given in Theorem A.2, third line uses Young inequality and finally fourth line uses Poincaré-Wirtinger inequality given in Theorem A.1. By choosing $\alpha$ such that $C - C'\frac{\alpha}{2} > 0$, we obtain that $\lim_{\varepsilon \to 0} \mathcal{F}(v_i, v_e) = +\infty$ when $\|v_i\|_{H^1(\Omega_i^\varepsilon)} + \|v_e\|_{H^1(\Omega_e^\varepsilon)} \to +\infty$.

**Regularity.** We identify $\Omega$ with a torus, so that the boundaries of $\Omega_i^\varepsilon$ and $\Omega_e^\varepsilon$ consist solely of the membrane interface $\Gamma_m^\varepsilon$, which is assumed to be infinitely smooth. $(u_i^\varepsilon - u_e^\varepsilon)/\varepsilon \in H^{1/2}(\Gamma_m^\varepsilon)$ so there exists $v \in H^1(\Omega_i^\varepsilon)$ such that $v|_{\Gamma_m^\varepsilon} = (u_i^\varepsilon - u_e^\varepsilon)/\varepsilon$. We have

$$ \nabla (I_m(v)) = I'_m(v) \nabla v. $$

By Assumption 2, $I'_m(v) \in L^\infty(\Omega_i^\varepsilon)$ so $I_m(v) \in H^1(\Omega_i^\varepsilon)$ and therefore

$$ \varepsilon I_m \left( \frac{u_i^\varepsilon - u_e^\varepsilon}{\varepsilon} \right) \in H^{1/2}(\Gamma_m^\varepsilon). $$

As $u_i^\varepsilon$ verifies

$$
\begin{align}
-\nabla \cdot (\sigma_i \nabla u_i^\varepsilon) &= f && \in L^2(\Omega_i^\varepsilon), \nonumber \\
-\sigma_i \nabla u_i^\varepsilon \cdot n_i &= \varepsilon I_m \left( \frac{u_i^\varepsilon - u_e^\varepsilon}{\varepsilon} \right) && \in H^{1/2}(\Gamma_m^\varepsilon), \tag{6}
\end{align}
$$

we have $u_i^\varepsilon \in H^2(\Omega_i^\varepsilon)$ (by elliptic regularity, see Proposition 2.5.2.3 of [27]). With the same reasoning, we show that $u_e^\varepsilon \in H^2(\Omega_e^\varepsilon)$.

We prove by induction that $u_i^\varepsilon \in H^k(\Omega_i^\varepsilon)$ and $u_e^\varepsilon \in H^k(\Omega_e^\varepsilon)$ for any $k \ge 2$. At step $k$, we have $v \in H^k(\Omega_i^\varepsilon)$. If we show that $I_m(v) \in H^k(\Omega_i^\varepsilon)$, then $\varepsilon I_m((u_i^\varepsilon - u_e^\varepsilon)/\varepsilon) \in H^{k-1/2}(\Gamma_m^\varepsilon)$ and, as $f \in C^\infty(\Omega_i^\varepsilon)$ and $u_i^\varepsilon$ still verifies Equation (6), we would get $u_i^\varepsilon \in H^{k+1}(\Omega_i^\varepsilon)$, using for instance Theorem 2.5.1.1 of [27] (the reasoning is the same for $u_e^\varepsilon$). To demonstrate that $I_m(v) \in H^k(\Omega_i^\varepsilon)$, let $\alpha \in \mathbb{N}^d$ and $\beta \in \mathbb{N}^d$, such that $|\alpha| = k$, $|\beta| = 2$ (we denote by $\beta_1$ and $\beta_2$ the two non-zero components of $\beta$) and $\beta \le \alpha$, we have

$$ \partial^\alpha (I_m(v)) = \partial^{\alpha-\beta} \partial^\beta (I_m(v)) = \partial^{\alpha-\beta} (I_m''(v) \partial^{\beta_2} v \partial^{\beta_1} v + I_m'(v) \partial^\beta v), $$

which gives, by Leibniz formula,

$$ \sum_{\nu \le \alpha - \beta} \binom{\alpha - \beta}{\nu} \left( \partial^\nu \left( \partial^{\beta_2} v \partial^{\beta_1} v \right) \partial^{\alpha - \beta - \nu} (I_m''(v)) + \partial^{\nu + \beta} v \partial^{\alpha - \beta - \nu} (I_m'(v)) \right). $$
The terms $\partial^{\alpha-\beta-\nu}(I_m''(v))$ and $\partial^{\alpha-\beta-\nu}(I_m'(v))$ contain derivatives of $v$ of order at most $k-2$ so, by Sobolev injection $H^k(\Omega_i^\epsilon) \subset W^{k-2,\infty}(\Omega_i^\epsilon)$ (see for instance [24], Section 5.6, and recall that we work in $\mathbb{R}^d$, with $d=2$ or 3), these terms are in $L^\infty(\Omega_i^\epsilon)$. Moreover, by Sobolev injection $H^k(\Omega_i^\epsilon) \subset W^{k-1,4}(\Omega_i^\epsilon)$, we have $\partial^\nu(\partial^{\beta_2}v\partial^{\beta_1}v) \in L^2(\Omega_i^\epsilon)$. And finally, $\partial^{\nu+\beta}v \in L^2(\Omega_i^\epsilon)$ because $|\nu+\beta| \le k$. This gives $\partial^\alpha(I_m(v)) \in L^2(\Omega_i^\epsilon)$, which concludes the induction. $u_i^\epsilon$ and $u_e^\epsilon$ being of regularity $H^k$ for all $k \in \mathbb{N}$, they belong to $C^\infty$.

**Energy estimates.** To conclude the proof, it remains to establish the energy estimates. We recall that $v_m^\epsilon = u_i^\epsilon - u_e^\epsilon$. Setting $\phi_i = u_i^\epsilon$ and $\phi_e = u_e^\epsilon$ in the variational formulation (5) leads to

$$ \sigma_i \|\nabla u_i^\epsilon\|_{L^2(\Omega_i^\epsilon)}^2 + \sigma_e \|\nabla u_e^\epsilon\|_{L^2(\Omega_e^\epsilon)}^2 + \varepsilon \int_{\Gamma_m^\epsilon} I_m \left(\frac{v_m^\epsilon}{\varepsilon}\right) v_m^\epsilon \\ = (f, u_i^\epsilon)_{L^2(\Omega_i^\epsilon)} + (f, u_e^\epsilon)_{L^2(\Omega_e^\epsilon)} $$

where the right hand side term can be bounded, using Schwarz inequality, by

$$ \left( \|u_i^\epsilon\|_{H^1(\Omega_i^\epsilon)} + \|u_e^\epsilon\|_{H^1(\Omega_e^\epsilon)} \right) \|f\|_{L^2(\Omega)} $$

By assumption, $S_m(\lambda) \ge S_m^0$ for any $\lambda$ and $\varepsilon$ is a small parameter then one can assume that $\varepsilon \le 1$. Hence, by using the norm equivalence given in Theorem A.2 (see Appendix A),

$$ \begin{aligned} \|u_i^\epsilon\|_{H^1(\Omega_i^\epsilon)}^2 + \|u_e^\epsilon\|_{H^1(\Omega_e^\epsilon)}^2 & \le C \left( \sigma_i \|\nabla u_i^\epsilon\|_{L^2(\Omega_i^\epsilon)}^2 + \sigma_e \|\nabla u_e^\epsilon\|_{L^2(\Omega_e^\epsilon)}^2 + \varepsilon S_m^0 \|v_m^\epsilon\|_{L^2(\Gamma_m^\epsilon)}^2 \right) \\ & \le C \left( \sigma_i \|\nabla u_i^\epsilon\|_{L^2(\Omega_i^\epsilon)}^2 + \sigma_e \|\nabla u_e^\epsilon\|_{L^2(\Omega_e^\epsilon)}^2 + S_m^0 \|v_m^\epsilon\|_{L^2(\Gamma_m^\epsilon)}^2 \right) \\ & \le C \left( \sigma_i \|\nabla u_i^\epsilon\|_{L^2(\Omega_i^\epsilon)}^2 + \sigma_e \|\nabla u_e^\epsilon\|_{L^2(\Omega_e^\epsilon)}^2 + \int_{\Gamma_m^\epsilon} S_m \left(\frac{v_m^\epsilon}{\varepsilon}\right) (v_m^\epsilon)^2 \right) \\ & \le C \left( \|u_i^\epsilon\|_{H^1(\Omega_i^\epsilon)} + \|u_e^\epsilon\|_{H^1(\Omega_e^\epsilon)} \right) \|f\|_{L^2(\Omega)} \end{aligned} $$

where $C$ is a positive constant. Therefore $\|u_i^\epsilon\|_{H^1(\Omega_i^\epsilon)}$ and $\|u_e^\epsilon\|_{H^1(\Omega_e^\epsilon)}$ are uniformly bounded (i.e. relatively to $\varepsilon$). Moreover,

$$ \begin{aligned} S_m^0 \|v_m^\epsilon\|_{L^2(\Gamma_m^\epsilon)}^2 & \le \int_{\Gamma_m^\epsilon} S_m \left(\frac{v_m^\epsilon}{\varepsilon}\right) (v_m^\epsilon)^2 \\ & \le \left( \|u_i^\epsilon\|_{H^1(\Omega_i^\epsilon)} + \|u_e^\epsilon\|_{H^1(\Omega_e^\epsilon)} \right) \|f\|_{L^2(\Omega)} \\ & \le C \|f\|_{L^2(\Omega)} \end{aligned} $$

thus the final estimate.

□

**Remark 2.4.** The first two estimates imply $\varepsilon \|v_m^\epsilon\|_{L^2(\Gamma_m^\epsilon)}^2 \le C \|f\|_{L^2(\Omega)}$ (see Theorem A.2 in Appendix A) but we showed that, due to the structure of the equations, a stronger result holds, namely the third estimate $\|v_m^\epsilon\|_{L^2(\Gamma_m^\epsilon)} \le C \|f\|_{L^2(\Omega)}$.

# 3 Tissue scale model

After a context reviewing related contribution, the new proposed model – which forms the main contribution of this work – is introduced. It is derived from an asymptotic analysis of the microscopic system and consists in a limit problem (order 0) and a first-order correction. We also present numerical simulations to illustrate the behavior of the model solution. For the sake of readability, the rigorous justification of the asymptotic derivation of the model is deferred to the next section.
## 3.1 Context and related work

Regarding the asymptotic analysis, the non-linearity of $\lambda \mapsto I_m(\lambda)$ is the first main difficulty. Let us first consider a linear mapping, i.e. $I_m(\lambda) = S_m^0\lambda$. Under this assumption, the homogenization of System (4) falls into a broader class of problems where the interfacial flux takes the form $\varepsilon^\gamma \kappa v_m$, with $\gamma \in \mathbb{R}$ and $\kappa > 0$. Such a setting arises in various contexts, including heat conduction in composite media with interfacial barriers [39], modeling of contact between two elastic bodies [23], electric conduction in biological tissues [20], and diffusion in fissured media [22]. A detailed discussion of the possible values of $\gamma$ can be found in [8, 11], and further generalizations are presented in [47]. For the treatment of nonlinear interfacial resistance, we refer to [21].

However, all these works concentrate on the limit problem - i.e., the homogenized system obtained as $\varepsilon \to 0$ - which in our setting reads, as will be shown later,

$$u_i^0 = u_e^0 := u^0, \quad \text{in } \Omega, \qquad -\nabla \cdot \left( (\Sigma_i + \Sigma_e)\nabla u^0 \right) = f, \quad \text{in } \Omega, \quad \int_{\Omega} u^0 = 0, \qquad (7)$$

with periodic boundary conditions and effective tensors $\Sigma_i$ and $\Sigma_e$ to be defined later. System (7) corresponds to a monodomain model, in the sense that a single macroscopic potential $u^0$ characterizes the system at the limit. Notably, it does not retain any information about the flux transmission across the membrane.

To obtain a more precise approximation than the one given by the limit problem, we need to extend the asymptotic expansion to higher orders - this constitutes the second main difficulty of our work. We follow a procedure already explored in various contexts. In [7], the authors highlight the necessity of a second-order corrector in the homogenization of conductive-radiative heat transfer in perforated domains, particularly when the source term is a periodically oscillating function, in order to achieve an accurate temperature approximation. Similarly, in [51], the full two-scale expansion of the equation $-\nabla \cdot (A(x/\varepsilon)\nabla u^\varepsilon) = f$ is rigorously derived in the framework of elasticity, for a fixed domain with macroscopic periodic boundary conditions (see also [16] for the quasilinear case), leading to what is referred to as an "infinite-order" homogenized solution. More generally, these works propose a formal expansion, whose validity is ensured by proving that the difference between the exact solution and the expansion satisfies a well-posed PDE, which allows for quantitative error estimates. We adopt the same approach in our work - see the discussion in the beginning of Section 4. Such high order development is also investigated in [6], where the authors compare classical asymptotic expansions with Bloch wave expansions for periodic homogenization, and advocate for higher-order corrections, notably suggesting that the error is better captured with a fourth-order term. In [25], the complete two-scale expansion of a perturbed convection equation is justified using two-scale convergence techniques (see also [26] for related results in the context of gyrokinetic approximations). In [56], effective models for linear elastic beams are obtained by retaining asymptotic terms up to order 2. To the best of our knowledge, the rigorous justification of a full two-scale expansion for System (4) together with the derivation of a non-linear physical model for PFA is original and constitutes the core contribution of the present article.

## 3.2 Limit problem

In the following, the $\alpha$ subscript designates either i or e, for the intra- and extracellular quantities, respectively.

We can have a first insight on the behaviour of the microscopic system of equations in the asymptotic regime by studying the limit problem as $\varepsilon \to 0$. Using the estimates from
Theorem 2.3, we apply in this section the theory of two-scale convergence (for an introduction, see [5], [4]) to derive the limit problem in Theorem 3.3. This problem is defined on the macroscopic domain $\Omega$, yet it retains information about the microscopic structure through the solutions of the so-called cell problems introduced in the following proposition and defined on the single cell domain (see Section 2.1).

**Proposition-Definition 3.1.** Let $F \in C_{\text{per}}^{\infty}(Y_{\alpha})$ and $G \in C_{\text{per}}^{\infty}(\Gamma)$ such that $\int_{Y_{\alpha}} F = \int_{\Gamma} G$. We denote by

$$ \mathcal{P}_{\alpha}(F, G) := \chi, $$

the unique solution in $C_{\text{per}}^{\infty}(Y_{\alpha})$ of

$$ -\Delta\chi = F, \text{ in } Y_{\alpha}, \qquad (8a) $$

$$ -\nabla\chi \cdot n_{\alpha} = G, \text{ on } \Gamma, \qquad (8b) $$

$$ \int_{Y_{\alpha}} \chi = 0. \qquad (8c) $$

It verifies

$$ \|\chi\|_{H^1(Y_{\alpha})} \le C \left( \|F\|_{L^2(Y_{\alpha})} + \|G\|_{H^{-\frac{1}{2}}(\Gamma)} \right). $$

*Proof.* The proof is classic and the compatibility condition $\int_{Y_{\alpha}} F = \int_{\Gamma} G$ is used to show the equivalence between the variational formulation and System (8). By identifying $Y$ to a torus and using Theorem 2.5.1.1 of [27], we conclude that $\chi \in C_{\text{per}}^{\infty}(Y_{\alpha})$. $\square$

**Definition 3.2 (Cell problems).** For $1 \le j \le d$, we define

$$ \chi_i^j = \mathcal{P}_i(0, e^j \cdot n_i), \quad \chi_e^j = \mathcal{P}_e(0, e^j \cdot n_e), $$

$$ \chi_i = (\chi_i^1, \dots, \chi_i^d)^T, \quad \chi_e = (\chi_e^1, \dots, \chi_e^d)^T, \quad [\chi] = \chi_i - \chi_e, $$

where $e^j$ is the j-th vector of the canonical base of $\mathbb{R}^d$. For consistency with the development at any order, we also introduce the following notations

$$ (\chi_i^{1,0,1})_j = \chi_i^j, \quad (\chi_e^{1,0,1})_j = \chi_e^j, \quad \chi_i^{0,1,1} = \chi_e^{0,1,1} = 1, \quad \theta_i^{0,0,1} = \theta_e^{0,0,1} = 0. $$

**Proposition 3.3 (Limit problem).** Under Assumptions 1 and 2, the sequence $(u_i^{\varepsilon}, u_e^{\varepsilon})$ of solutions of System (4) in $C_{\text{per}}^{\infty}(\Omega_i^{\varepsilon}) \times \overline{C}_{\text{per}}^{\infty}(\Omega_e^{\varepsilon})$ verifies, for all $\phi \in \mathcal{D}(\Omega; C_{\text{per}}^{\infty}(Y))$,

$$ \lim_{\varepsilon \to 0} \int_{\Omega_i^{\varepsilon}} u_i^{\varepsilon}(x) \phi \left( x, \frac{x}{\varepsilon} \right) = \int_{\Omega} \int_{Y_i} u_i^0(x) \phi(x, y), $$

$$ \lim_{\varepsilon \to 0} \int_{\Omega_e^{\varepsilon}} u_e^{\varepsilon}(x) \phi \left( x, \frac{x}{\varepsilon} \right) = \int_{\Omega} \int_{Y_e} u_e^0(x) \phi(x, y), $$

where

$$ u_i^0 = u_e^0 = u^0, \quad \Omega, \qquad (9) $$

and $u^0$ is the unique solution in $C_{\text{per}}^{\infty}(\Omega)$ of

$$ -\nabla \cdot \left( (\Sigma_i + \Sigma_e) \nabla u^0 \right) = f, \quad \Omega, \qquad (10a) $$

$$ \int_{\Omega} u^0 = 0, \qquad (10b) $$
where

$$ \Sigma_i = \sigma_i \left( |Y_i| \mathrm{Id} + \int_{Y_i} \nabla \chi_i \right), \quad \Sigma_e = \sigma_e \left( |Y_e| \mathrm{Id} + \int_{Y_e} \nabla \chi_e \right). \qquad (11) $$

Remark 3.4. The tensors $\Sigma_i$ and $\Sigma_e$ can be rewritten as follows:

$$ \begin{aligned} (\Sigma_i)_{j,k} &= \sigma_i \int_{Y_i} (e^j + \nabla \chi_i^j) \cdot (e^k + \nabla \chi_i^k), \\ (\Sigma_e)_{j,k} &= \sigma_e \int_{Y_e} (e^j + \nabla \chi_e^j) \cdot (e^k + \nabla \chi_e^k), \end{aligned} $$

or even

$$ \begin{aligned} (\Sigma_i)_{j,k} &= \sigma_i \int_{\partial Y \cap \partial Y_i} y_k (e^j + \nabla \chi_i^j) \cdot n_i, \\ (\Sigma_e)_{j,k} &= \sigma_e \int_{\partial Y \cap \partial Y_e} y_k (e^j + \nabla \chi_e^j) \cdot n_e. \end{aligned} $$

The first expression shows that the tensors $\Sigma_i$ and $\Sigma_e$ are symmetric, while the second shows that if $\partial Y \cap \partial Y_i$ (resp. $\partial Y \cap \partial Y_e$) has a null measure, then $\Sigma_i$ (resp. $\Sigma_e$) is null. In addition, for all $\xi \in \mathbb{R}^d \setminus \{0\}$, we have

$$ \xi^T \cdot \Sigma_i \cdot \xi = \sigma_i \int_{Y_i} |\xi + \nabla \psi_{\xi}|^2 \geq 0, $$

where $\psi_{\xi} = \sum \xi_j \chi_i^j$. If $\xi^T \cdot \Sigma_i \cdot \xi = 0$, then

$$ \nabla \psi_{\xi} = -\xi, $$

and therefore

$$ \psi_{\xi} = -\xi \cdot y + C, $$

for some constant $C$, which is not possible if $\partial Y \cap \partial Y_i$ has a non null measure, because then $\psi_{\xi}$ has periodic constraints on the shared boundary. Therefore $\Sigma_i$ (resp. $\Sigma_e$) is positive definite if and only if $\partial Y \cap \partial Y_i$ (resp. $\partial Y \cap \partial Y_e$) has a non null measure. In this work, we assumed that $\Omega_e^\epsilon$ are connected, so $\Sigma_e$ are symmetric positive definite, ensuring that System (10) is elliptic.

**Proof of Proposition 3.3.** We first establish the well-posedness of System (10), and then apply classical two-scale convergence results to prove the convergence.

**Existence of System (10).** Thanks to the gauge condition, the regularity of $f$, the regularity of $\Sigma_i$ and $\Sigma_e$ (through the regularity of $\chi_i$ and $\chi_e$) and the periodic boundary conditions, the well-posedness in $C_{\mathrm{per}}^{\infty}(\Omega)$ follows directly identifying $\Omega$ to a torus and using for instance Theorem 2.5.1.1 of [27].

**Convergence.** This proof does not contain new difficulties and is a direct application of two-scale convergence. Using estimates of Theorem 2.3, we can apply a classical result of two-scale convergence: there exists

$$ (u_i^0, u_e^0) \in (H^1(\Omega))^2, \quad (u_i^1, u_e^1) \in L^2[\Omega; \overline{H_{\mathrm{per}}^1}(Y_i)] \times L^2[\Omega; \overline{H_{\mathrm{per}}^1}(Y_e)] $$
such that for all $\phi \in \mathcal{D}(\Omega; C_{\text{per}}^{\infty}(Y))$ and $\psi \in (\mathcal{D}(\Omega; C_{\text{per}}^{\infty}(Y)))^d$, the following convergences hold when $\varepsilon \to 0$ and up to a subsequence

$$ \int_{\Omega_i^{\varepsilon}} u_i^{\varepsilon}(x) \phi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{Y_i} u_i^0(x) \phi(x, y), \quad (12a) $$

$$ \int_{\Omega_e^{\varepsilon}} u_e^{\varepsilon}(x) \phi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{Y_e} u_e^0(x) \phi(x, y), \quad (12b) $$

$$ \int_{\Omega_i^{\varepsilon}} \nabla u_i^{\varepsilon}(x) \cdot \psi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{Y_i} \left( \nabla u_i^0(x) + \nabla_y u_i^1(x, y) \right) \cdot \psi(x, y), \quad (12c) $$

$$ \int_{\Omega_e^{\varepsilon}} \nabla u_e^{\varepsilon}(x) \cdot \psi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{Y_e} \left( \nabla u_e^0(x) + \nabla_y u_e^1(x, y) \right) \cdot \psi(x, y), \quad (12d) $$

$$ \varepsilon \int_{\Gamma_m^{\varepsilon}} u_i^{\varepsilon}(x) \psi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{\Gamma} u_i^0(x) \psi(x, y), \quad (12e) $$

$$ \varepsilon \int_{\Gamma_m^{\varepsilon}} u_e^{\varepsilon}(x) \psi \left( x, \frac{x}{\varepsilon} \right) \to \int_{\Omega} \int_{\Gamma} u_e^0(x) \psi(x, y). \quad (12f) $$

Since $\varepsilon \int_{\Gamma_m^{\varepsilon}} (u_i^{\varepsilon}(x) - u_e^{\varepsilon}(x)) \psi(x, \frac{x}{\varepsilon}) \to \int_{\Omega} \int_{\Gamma} (u_i^0(x) - u_e^0(x)) \psi(x, y)$, another classical result of two-scale convergence (Proposition 2.5 in [4]) gives

$$ \int_{\Omega} \int_{\Gamma} (u_i^0(x) - u_e^0(x))^2 \leq \lim_{\varepsilon \to 0} \varepsilon \int_{\Gamma_m^{\varepsilon}} (u_i^{\varepsilon} - u_e^{\varepsilon})^2. $$

We showed in Theorem 2.3 that $\|u_i^{\varepsilon} - u_e^{\varepsilon}\|_{L^2(\Gamma_m^{\varepsilon})} \leq C$, hence $u_i^0 = u_e^0$.

To prove that (10) holds, we start from the variational formulation of System (4). We can show that, for all $(\phi_i, \phi_e) \in (\mathcal{D}(\Omega; C_{\text{per}}^{\infty}(Y)) \cap C_{\text{per}}^{\infty}(\Omega \times Y))^2$, one has

$$ \begin{aligned} & \int_{\Omega_i^{\varepsilon}} \sigma_i \nabla u_i^{\varepsilon}(x) \cdot \nabla \phi_i \left( x, \frac{x}{\varepsilon} \right) dx + \int_{\Omega_e^{\varepsilon}} \sigma_e \nabla u_e^{\varepsilon}(x) \cdot \nabla \phi_e \left( x, \frac{x}{\varepsilon} \right) dx \\ & + \varepsilon \int_{\Gamma_m^{\varepsilon}} I_m \left( \frac{u_i^{\varepsilon}(x) - u_e^{\varepsilon}(x)}{\varepsilon} \right) \left( \phi_i \left( x, \frac{x}{\varepsilon} \right) - \phi_e \left( x, \frac{x}{\varepsilon} \right) \right) \\ & = \int_{\Omega_i^{\varepsilon}} f(x) \phi_i \left( x, \frac{x}{\varepsilon} \right) + \int_{\Omega_e^{\varepsilon}} f(x) \phi_e \left( x, \frac{x}{\varepsilon} \right), \end{aligned} \quad (13) $$

where $\nabla \phi_{i,e}(x, x/\varepsilon) = \nabla_x \phi_{i,e}(x, x/\varepsilon) + \varepsilon^{-1} \nabla_y \phi_{i,e}(x, x/\varepsilon)$. Let us multiply Eq. (13) by $\varepsilon$ and make it tend towards 0. The surface term verifies

$$ \begin{aligned} & \varepsilon^2 \int_{\Gamma_m^{\varepsilon}} I_m \left( \frac{u_i^{\varepsilon}(x) - u_e^{\varepsilon}(x)}{\varepsilon} \right) \left( \phi_i \left( x, \frac{x}{\varepsilon} \right) - \phi_e \left( x, \frac{x}{\varepsilon} \right) \right) \\ & \leq S_m^1 \sqrt{\varepsilon} \|u_i^{\varepsilon} - u_e^{\varepsilon}\|_{L^2(\Gamma_m^{\varepsilon})} \sqrt{\varepsilon} \left\| \phi_i \left( \cdot, \frac{\cdot}{\varepsilon} \right) - \phi_e \left( \cdot, \frac{\cdot}{\varepsilon} \right) \right\|_{L^2(\Gamma_m^{\varepsilon})} \to 0 \end{aligned} $$

where $\sqrt{\varepsilon} \|\phi_i(\cdot, \cdot/\varepsilon) - \phi_e(\cdot, \cdot/\varepsilon)\|_{L^2(\Gamma_m^{\varepsilon})}$ is bounded because $\phi_i$ and $\phi_e$ are smooth enough (see for instance Lemma 2.4 in [4]) and $\sqrt{\varepsilon} \|u_i^{\varepsilon} - u_e^{\varepsilon}\|_{L^2(\Gamma_m^{\varepsilon})} \to 0$ because of the surface estimate in Theorem 2.3. Therefore, using Eq. (12), one has at the limit

$$ \int_{\Omega} \int_{Y_i} \sigma_i \left( \nabla u_e^0 + \nabla_y u_i^1 \right) \cdot \nabla_y \phi_i + \int_{\Omega} \int_{Y_e} \sigma_e \left( \nabla u_e^0 + \nabla_y u_e^1 \right) \cdot \nabla_y \phi_e = 0. $$

Remark that $v_i(x, y) := \sum_{j=1}^d \chi_i^j(y) \frac{\partial u_e^0}{\partial x_j}(x)$ and $v_e(x, y) := \sum_{j=1}^d \chi_e^j(y) \frac{\partial u_e^0}{\partial x_j}(x)$ also verify

$$ \int_{\Omega} \int_{Y_i} \sigma_i \left( \nabla u_e^0 + \nabla_y v_i \right) \cdot \nabla_y \phi_i + \int_{\Omega} \int_{Y_e} \sigma_e \left( \nabla u_e^0 + \nabla_y v_e \right) \cdot \nabla_y \phi_e = 0. $$
It implies that for all $(\phi_i, \phi_e) \in (\mathcal{D}(\Omega; \mathcal{C}_{\text{per}}^{\infty}(Y)) \cap \mathcal{C}_{\text{per}}^{\infty}(\Omega \times Y))^2,$

$$ \int_{\Omega} \int_{Y_i} \sigma_i \nabla_y (u_i^1 - v_i) \cdot \nabla_y \phi_i + \int_{\Omega} \int_{Y_e} \sigma_e \nabla_y (u_e^1 - v_e) \cdot \nabla_y \phi_e = 0. \quad (14) $$

Now, we use $\phi_i = \phi_e := \phi \in \mathcal{D}(\Omega) \cap \mathcal{C}_{\text{per}}^{\infty}(\Omega)$ in (13) and make $\varepsilon$ tend towards 0. Using (12), one has

$$ \int_{\Omega} \int_{Y_i} \sigma_i (\nabla u_e^0 + \nabla_y u_i^1) \cdot \nabla \phi + \int_{\Omega} \int_{Y_e} \sigma_e (\nabla u_e^0 + \nabla_y u_e^1) \cdot \nabla \phi = |Y| \int_{\Omega} f \phi. $$

Using (14) and $|Y| = 1$, one finally has

$$ \int_{\Omega} \sum_{j=1}^{d} \left( \sigma_i \left( |Y_i| e^j + \int_{Y_i} \nabla_y \chi_i^j \right) + \sigma_e \left( |Y_e| e^j + \int_{Y_e} \nabla_y \chi_e^j \right) \right) \frac{\partial u_e^0}{\partial x_j} \cdot \nabla \phi = \int_{\Omega} f \phi, $$

which rewrites into

$$ \int_{\Omega} (\Sigma_i + \Sigma_e) \nabla u_e^0 \cdot \nabla \phi = \int_{\Omega} f \phi. $$

Moreover, $\int_{\Omega} u_e^0 = |Y_e|^{-1} \lim_{\varepsilon \to 0} \int_{\Omega_e^{\varepsilon}} u_e^{\varepsilon} = 0.$

Therefore, since System (10) has a unique solution $u_e^0 = u^0$ and the whole sequence $(u_i^\varepsilon, u_e^\varepsilon)$ converges. $\square$

### 3.3 Ansatz and first order formal derivation

The leading-order (limit) problem (9)-(10) does not involve the membrane conductance $S_m$ at all. This is a structural limitation of the order 0 model: the essential nonlinear effects due to electroporation are completely absent from the limit equations. As a result, describing the impact of $S_m$ - and thus of electroporation - requires analyzing higher-order terms in the asymptotic expansion. This observation is consistent with both our numerical findings and with the fact that $\varepsilon = 10^{-2}$ is not particularly small. In contrast, in classical electrophysiology without PFA, $\varepsilon$ is typically of the order $10^{-3}$, and the order 0 model already provides a rich and accurate description of the dynamics.

We apply the method of two-scale asymptotic expansion to formally obtain the first-order development of the solution $(u_i^\varepsilon, u_e^\varepsilon)$ of System (4). This method consists in introducing a macroscopic variable $x$ and a microscopic variable $y = x/\varepsilon$, thereby explicitly separating the scales involved in the model, and in assuming the following ansatz

$$ u_i^\varepsilon(x) = u_i^0 \left( x, \frac{x}{\varepsilon} \right) + \sum_{n=1}^{\infty} \varepsilon^n u_i^n \left( x, \frac{x}{\varepsilon} \right), $$

$$ u_e^\varepsilon(x) = u_e^0 \left( x, \frac{x}{\varepsilon} \right) + \sum_{n=1}^{\infty} \varepsilon^n u_e^n \left( x, \frac{x}{\varepsilon} \right), $$

where the functions $y \mapsto u_i^n(x, y)$ and $y \mapsto u_e^n(x, y)$ are Y-periodic for all $n \in \mathbb{N}$. This ansatz will be justified in Section 4.3 where we derive rigorous error estimates.

Plugging the ansatz in System (4) and identifying the powers of $\varepsilon$ results in a cascade of equations that uniquely define all $u_i^n$ and $u_e^n$. It involves the characterization of the macroscopic parts (which only depend on $x$), denoted by a superscript $\sharp$:

$$ u_i^{n, \sharp}(x) := \frac{1}{|Y_i|} \int_{Y_i} u_i^n(x, y), \quad u_e^{n, \sharp}(x) := \frac{1}{|Y_e|} \int_{Y_e} u_e^n(x, y). $$
We then denote

$$v_m^n := u_i^n - u_e^n, \quad v_m^{n, \sharp} = u_i^{n, \sharp} - u_e^{n, \sharp}.$$

The first terms of each ansatz, $u_i^0$ and $u_e^0$, are the two-scale limits of $u_i^\xi$ and $u_e^\xi$. We showed in the previous section that they only depend on $x$ and are equal to $u^0$ (as defined in Proposition 3.3).

To handle the non linearity, we remark that $v_m^0 = 0$ and use Taylor expansion around $\varepsilon v_m^1$: plugging the ansatz into $I_m$ gives, at first order,

$$I_m \left( \frac{v_m^\varepsilon}{\varepsilon} \right) = I_m(v_m^1) + O(\varepsilon).$$

This expansion of the non linearity will be extended and proved later in Section 4.

Let us inject the ansatz into the microscopic equations and identify the powers of $\varepsilon$. The first system of equations of the cascade is

$$-\nabla_y \cdot (\sigma_i (\nabla_y u_i^1 + \nabla_x u_i^0)) = 0, \qquad \Omega \times Y_i, \qquad (15a)$$

$$-\nabla_y \cdot (\sigma_e (\nabla_y u_e^1 + \nabla_x u_e^0)) = 0, \qquad \Omega \times Y_e, \qquad (15b)$$

$$\sigma_i (\nabla_y u_i^1 + \nabla_x u_i^0) \cdot n_i = \sigma_e (\nabla_y u_e^1 + \nabla_x u_e^0) \cdot n_i, \qquad \Omega \times \Gamma, \qquad (15c)$$

$$-\sigma_i (\nabla_y u_i^1 + \nabla_x u_i^0) \cdot n_i = 0, \qquad \Omega \times \Gamma, \qquad (15d)$$

$$+ x \text{ and } y \text{ periodic conditions.} \qquad (15e)$$

Considering $x$ as a parameter and $(y \mapsto u_i^1(x, y), y \mapsto u_e^1(x, y))$ as the unknown, this system has a unique solution, up to a constant which depends on $x$. Note, as in the proof of Theorem 3.3, that $(\chi_i \cdot \nabla u^0, \chi_e \cdot \nabla u^0)$ is a solution. Therefore, we can write

$$u_i^1(x, y) = \chi_i(y) \cdot \nabla u^0(x) + u_i^{1, \sharp}(x), \qquad (16a)$$

$$u_e^1(x, y) = \chi_e(y) \cdot \nabla u^0(x) + u_e^{1, \sharp}(x), \qquad (16b)$$

where $u_i^{1, \sharp}, u_e^{1, \sharp}$ are periodic functions depending only on $x$. The next system reads

$$-\nabla_y \cdot (\sigma_i (\nabla_y u_i^2 + \nabla_x u_i^1)) = \nabla_x \cdot (\sigma_i (\nabla_y u_i^1 + \nabla_x u_i^0)) + f, \qquad \Omega \times Y_i, \qquad (17a)$$

$$-\nabla_y \cdot (\sigma_e (\nabla_y u_e^2 + \nabla_x u_e^1)) = \nabla_x \cdot (\sigma_e (\nabla_y u_e^1 + \nabla_x u_e^0)) + f, \qquad \Omega \times Y_e, \qquad (17b)$$

$$\sigma_i (\nabla_y u_i^2 + \nabla_x u_i^1) \cdot n_i = \sigma_e (\nabla_y u_e^2 + \nabla_x u_e^1) \cdot n_i, \qquad \Omega \times \Gamma, \qquad (17c)$$

$$-\sigma_i (\nabla_y u_i^2 + \nabla_x u_i^1) \cdot n_i = I_m(v_m^1), \qquad \Omega \times \Gamma, \qquad (17d)$$

$$+ x \text{ and } y \text{ periodic conditions.} \qquad (17e)$$

Again, considering $x$ as a parameter and $(y \mapsto u_i^2(x, y), y \mapsto u_e^2(x, y))$ as the unknown, this system has a unique solution, up to a constant which depends on $x$, if and only if the following compatibility condition is verified, for all $x \in \Omega$,

$$\int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y u_i^1 + \nabla_x u_i^0)) + |Y_i| f = \int_{\Gamma} I_m(v_m^1), \\ \int_{Y_e} \nabla_x \cdot (\sigma_e (\nabla_y u_e^1 + \nabla_x u_e^0)) + |Y_e| f = - \int_{\Gamma} I_m(v_m^1).$$
By using (16) and recalling that

$$ \Sigma_i = \sigma_i \left( |Y_i| \mathrm{Id} + \int_{Y_i} \nabla \chi_i \right), \quad \Sigma_e = \sigma_e \left( |Y_e| \mathrm{Id} + \int_{Y_e} \nabla \chi_e \right), $$

it rewrites

$$ \nabla \cdot (\Sigma_i \nabla u^0) + |Y_i| f = \int_{\Gamma} I_m((\chi_i - \chi_e) \cdot \nabla u^0 + u_i^{1,\#} - u_e^{1,\#}), $$

$$ \nabla \cdot (\Sigma_e \nabla u^0) + |Y_e| f = - \int_{\Gamma} I_m((\chi_i - \chi_e) \cdot \nabla u^0 + u_i^{1,\#} - u_e^{1,\#}). $$

Summing both equations, we get back to the limit problem 10,

$$ -\nabla \cdot (\Sigma_i + \Sigma_e) \nabla u^0 = f, \quad \Omega, \qquad (18) $$

and it remains the equation verified by $v_m^{1,\#}$,

$$ \nabla \cdot (\Sigma_i \nabla u^0) + |Y_i| f = \int_{\Gamma} I_m((\chi_i - \chi_e) \cdot \nabla u^0 + v_m^{1,\#}), \quad \Omega, \qquad (19) $$

which uniquely defines $v_m^{1,\#}$ as it will be shown later in Proposition 4.1 applied to $B(x, y) = [\chi](y) \cdot \nabla u^0(x)$ and $F(x) = \nabla \cdot \Sigma_i \nabla u^0(x) + |Y_i| f(x)$. These two equations allow to entirely define $u^0$ and $v_m^1$ (defined as the difference between the 2 equations of (16)) such that the compatibility condition of System (17) is verified. We can then express the solutions of this system thanks to new cell problems.

We recall that in Definition 3.2, we denote $\mathcal{P}_\alpha(F, G)$ as the unique solution in $C_{\mathrm{per}}^\infty(Y_\alpha)$ of System (8) with source terms $F$ and $G$.

**Definition 3.5.** We define for $1 \le j \le d$ and $1 \le k \le d$,

$$ (\chi_i^{1,0,2})_j = \mathcal{P}_i \left( \frac{1}{|Y_i|} \int_{\Gamma} S_m(v_m^1) \chi_i^j, S_m(v_m^1) \chi_i^j \right), $$

$$ (\theta_i^{1,0,2})_j = \mathcal{P}_i \left( \frac{1}{|Y_i|} \int_{\Gamma} S_m(v_m^1) \chi_e^j, S_m(v_m^1) \chi_e^j \right), $$

$$ \chi_i^{0,1,2} = \mathcal{P}_i \left( \frac{1}{|Y_i|} \int_{\Gamma} S_m(v_m^1), S_m(v_m^1) \right), $$

$$ (\chi_i^{2,0,2})_{j,k} = \mathcal{P}_i \left( 2 \frac{\partial \chi_i^j}{\partial y_k} - \frac{1}{|Y_i|} \int_{Y_i} \frac{\partial \chi_i^j}{\partial y_k}, \chi_i^j n_{i,k} \right), $$

and

$$ (\chi_e^{1,0,2})_j = \mathcal{P}_e \left( -\frac{1}{|Y_e|} \int_{\Gamma} S_m(v_m^1) \chi_e^j, -S_m(v_m^1) \chi_e^j \right), $$

$$ (\theta_e^{1,0,2})_j = \mathcal{P}_e \left( -\frac{1}{|Y_e|} \int_{\Gamma} S_m(v_m^1) \chi_i^j, -S_m(v_m^1) \chi_i^j \right), $$

$$ \chi_e^{0,1,2} = \mathcal{P}_e \left( -\frac{1}{|Y_e|} \int_{\Gamma} S_m(v_m^1), -S_m(v_m^1) \right), $$

$$ (\chi_e^{2,0,2})_{j,k} = \mathcal{P}_e \left( 2 \frac{\partial \chi_e^j}{\partial y_k} - \frac{1}{|Y_e|} \int_{Y_e} \frac{\partial \chi_e^j}{\partial y_k}, \chi_e^j n_{e,k} \right). $$

For consistency with the generalization at any order, we also introduce

$$ \theta_i^{0,1,2} = \chi_i^{0,1,2}, \quad \theta_e^{0,1,2} = \chi_e^{0,1,2}. $$
Thanks to these new correctors, we can explicitly write the solutions of System (17) up to constants $u_i^{2,\#}(x)$ and $u_e^{2,\#}(x)$,

$$u_i^2(x, y) = \chi_i(y) \cdot \nabla u_i^{1,\#}(x) + \chi_i^{2,0,2}(y) \cdot \nabla^2 u^0(x) + \frac{1}{\sigma_i} \chi_i^{0,1,2}(x, y) v_m^{1,\#}(x) \quad (20a)$$

$$+ \frac{1}{\sigma_i} (\chi_i^{1,0,2}(x, y) - \theta_i^{1,0,2}(x, y)) \cdot \nabla u^0(x) + u_i^{2,\#}(x), \quad (20b)$$

$$u_e^2(x, y) = \chi_e(y) \cdot \nabla u_e^{1,\#}(x) + \chi_e^{2,0,2}(y) \cdot \nabla^2 u^0(x) + \frac{1}{\sigma_e} \chi_e^{0,1,2}(x, y) v_m^{1,\#}(x) \quad (20c)$$

$$+ \frac{1}{\sigma_e} (-\chi_e^{1,0,2}(x, y) + \theta_e^{1,0,2}(x, y)) \cdot \nabla u^0(x) + u_e^{2,\#}(x). \quad (20d)$$

Similarly as before, we can show that the next system in the cascade has solutions if and only if some compatibility conditions are verified. The first one is

$$-\int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y u_i^2 + \nabla_x u_i^1)) - \int_{Y_e} \nabla_x \cdot (\sigma_e (\nabla_y u_e^2 + \nabla_x u_e^1)) = 0. \quad (21)$$

Combining it with Eq. (20), we deduce an equation on $u_e^{1,\#}$ which can be rewritten using the following lemma.

**Lemma 3.6.** Let $\chi \in C_{\text{per}}^{\infty}(Y_\alpha)$ and let $1 \le j \le d$, we have

$$\int_{Y_\alpha} \frac{\partial \chi}{\partial y_j} = \int_{Y_\alpha} \Delta \chi \chi_\alpha^j - \int_{\Gamma} \nabla \chi \cdot n_\alpha \chi_\alpha^j.$$

*Proof.*

$$\int_{Y_\alpha} \frac{\partial \chi}{\partial y_j} = \int_{\Gamma} \chi n_{\alpha,j} = - \int_{\Gamma} \chi \nabla \chi_\alpha^j \cdot n_\alpha = - \int_{Y_\alpha} \nabla \chi \cdot \nabla \chi_\alpha^j \\ = \int_{Y_\alpha} \Delta \chi \chi_\alpha^j - \int_{\Gamma} \nabla \chi \cdot n_\alpha \chi_\alpha^j.$$

□

After some computation, Eq. (21) can be rewritten as:

$$-\nabla \cdot (\Sigma_i \nabla u_i^{1,\#}) - \nabla \cdot (\Sigma_e \nabla u_e^{1,\#}) = \nabla \cdot \left( \int_{\Gamma} I_m(v_m^1)[\chi] \right).$$

Replacing $u_i^{1,\#}$ by $u_e^{1,\#} + v_m^{1,\#}$, we get an equation for $u_e^{1,\#}$

$$-\nabla \cdot \left( (\Sigma_i + \Sigma_e) \nabla u_e^{1,\#} \right) = \nabla \cdot \left( \Sigma_i \nabla v_m^{1,\#} \right) + \nabla \cdot \left( \int_{\Gamma} I_m(v_m^1)[\chi] \right), \quad \Omega, \quad (22)$$

that we complement with the gauge condition

$$\int_{\Omega} u_e^{1,\#} = 0.$$

Remark that $f_1 := \nabla \cdot (\Sigma_i \nabla v_m^{1,\#}) + \nabla \cdot (\int_{\Gamma} I_m(v_m^1)[\chi]) \in C_{\text{per}}^{\infty}(\Omega)$ verifies $\int_{\Omega} f_1 = 0$, as the divergence of periodic quantities. Therefore this system has a unique solution, as it will be shown in the next section in Theorem 4.2.

In order to progress in the expansion, the rest of the process is the same. The second compatibility condition defines $v_m^{2,\#} := u_i^{2,\#} - u_e^{2,\#}$, so that $u_i^d$ and $u_e^d$ can be expressed with new cell correctors and injected in the next compatibility condition which defines $v_m^{3,\#} := u_i^{3,\#} - u_e^{3,\#}$ and $u_e^{2,\#}$, and so forth. The reasoning by induction is summed up in Section 4, along with the rigorous justification of the full asymptotic expansion.
## 3.4 Effective model

The first-order model is composed of equations (10), (19), and (22) and can be rewritten as

$$u_e^0 - u_i^0 = 0, \quad \Omega, \qquad (23a)$$

$$-\nabla \cdot \left( \left( \Sigma_i \nabla u_i^0 + \Sigma_e \nabla u_e^0 \right) \right) = f, \quad \Omega, \qquad (23b)$$

$$2 \int_{\Gamma} I_m (u_i^{1,\#} - u_e^{1,\#} + (\chi_i - \chi_e) \cdot \nabla \frac{u_e^0 + u_i^0}{2}) = \nabla \cdot \left( \Sigma_i \nabla u_i^0 - \Sigma_e \nabla u_e^0 \right), \quad \Omega, \qquad (23c) \\ + (|Y_i| - |Y_e|) f$$

$$-\nabla \cdot \left( \left( \Sigma_i \nabla u_i^{1,\#} + \Sigma_e \nabla u_e^{1,\#} \right) \right) = \nabla \cdot \left( \int_{\Gamma} I_m (u_i^{1,\#} - u_e^{1,\#}) \right. \\ \left. + (\chi_i - \chi_e) \cdot \nabla \frac{u_e^0 + u_i^0}{2}) [\chi] \right), \quad \Omega, \qquad (23d)$$

$$\int_{\Omega} u_e^0 = \int_{\Omega} u_e^{1,\#} = 0, \qquad (23e)$$

along with periodic boundary conditions. It is interesting to derive from these 2 first orders an effective model at the order 1, that can be compared for instance to the classical bidomain model for cardiac electrophysiology [17].

Introducing the variables $(U_i^\epsilon, U_e^\epsilon)$ solution to the following problem

$$-\nabla \cdot (\Sigma_i \nabla U_i^\epsilon + \Sigma_e \nabla U_e^\epsilon) = f + \epsilon \nabla \cdot \left( \int_{\Gamma} I_m \left( \frac{U_i^\epsilon - U_e^\epsilon}{\epsilon} + (\chi_i - \chi_e) \cdot \nabla \frac{U_i^\epsilon + U_e^\epsilon}{2} \right) [\chi] \right) \qquad (24a)$$

$$2\epsilon \int_{\Gamma} I_m \left( \frac{U_i^{\epsilon} - U_e^{\epsilon}}{\epsilon} + (\chi_i - \chi_e) \cdot \nabla \frac{U_i^{\epsilon} + U_e^{\epsilon}}{2} \right) dy = \epsilon \left( \nabla \cdot (\Sigma_i \nabla U_i^{\epsilon} - \Sigma_e \nabla U_e^{\epsilon}) + (|Y_i| - |Y_e|)f \right). \qquad (24b)$$

One can show that if the solution $(U_i^\epsilon, U_e^\epsilon)$ of the above system exists, then thanks to an expansion in power of $\epsilon$ of the solution, tedious but simple calculations show that, at least formally, the 2 first orders of the expansion of $(U_i^\epsilon, U_e^\epsilon)$ satisfy the same problem as $(u_e^0, u_i^0)$ and $(u_i^{1,\#}, u_e^{1,\#})$.

This effective model is presented here for informational purposes only and is not used in the following numerical simulations, which deals separately with the 2 first orders. More precisely, the first-order model used in the subsequent numerical simulations and recommended for practical computations reads as follows:

$$-\nabla \cdot \left( (\Sigma_i + \Sigma_e) \nabla u^0 \right) = f, \quad \Omega, \qquad (25a)$$

$$\nabla \cdot (\Sigma_i \nabla u^0) + |Y_i| f = \int_{\Gamma} I_m ((\chi_i - \chi_e) \cdot \nabla u^0 + v_m^{1,\#}), \quad \Omega, \qquad (25b)$$

$$-\nabla \cdot (\Sigma_i + \Sigma_e) \nabla u_e^{1,\#} = \nabla \cdot (\Sigma_i \nabla v_m^{1,\#}) + \nabla \cdot \left( \int_{\Gamma} I_m (v_m^1) [\chi] \right), \quad \Omega, \qquad (25c)$$

$$\int_{\Omega} u^0 = \int_{\Omega} u_e^{1,\#} = 0, \qquad (25d)$$

along with periodic boundary conditions. This model enables the computation of approximations of the electric field, as detailed in the next section (see Eq. (26)), as well as other quantities of interest, as discussed in Section 3.5.3.

## 3.5 Numerical simulations

Once the paper accepted, the code will be publicly available online.
### 3.5.1 Implementation

In this section, we detail our strategy and our choices concerning the numerical illustrations.

**Algorithm** Numerical simulations are performed using the FreeFem++ software [28]. The spatial discretization and fixed-point convergence criteria have been selected finely enough to ensure that numerical errors remain negligible and do not interfere with the observed convergence or the validity of the results. The steps are the following:

* compute the order 1 cell correctors $\chi_\alpha$ by solving the cell problems (3.2), with $P^2$ Lagrangian finite elements and a Lagrange multiplier to force the gauge condition,

* compute $u^0$ by solving the limit problem (10) with $P^2$ Lagrangian finite elements,

* compute $v_m^{1,\#}$ by solving the non linear algebraic equation (19), using a fixed point method on each point of the mesh,

* compute $u_e^{1,\#}$ by solving the problem (22), using $P^2$ Lagrangian finite elements,

* compute the order 2 cell correctors $\chi_\alpha^{2,0,2}$, $\chi_\alpha^{0,1,2}$, $\chi_\alpha^{1,0,2}$, $\theta_\alpha^{1,0,2}$ by solving the cell problems given in Theorem 3.5, with $P^2$ Lagrangian finite elements and a Lagrange multiplier to force the gauge condition.

Then, for a given $\varepsilon$, the approximations of the potentials $u_\alpha^{1,\varepsilon}(x, \frac{x}{\varepsilon})$, and the electric fields $E_\alpha^0(x, \frac{x}{\varepsilon})$ and $E_\alpha^{1,\varepsilon}(x, \frac{x}{\varepsilon})$ are given by

$$u_{\alpha}^{1,\varepsilon}(x, y) := u^{0}(x) + \varepsilon u_{\alpha}^{1}(x, y), \qquad (26a)$$

$$E_{\alpha}^{0}(x, y) := - \left( \nabla u^{0}(x) + \nabla_{y} u_{\alpha}^{1}(x, y) \right), \qquad (26b)$$

$$E_{\alpha}^{1,\varepsilon}(x, y) := E_{\alpha}^{0}(x, y) - \varepsilon \left( \nabla_x u_{\alpha}^{1}(x, y) + \nabla_y u_{\alpha}^{2}(x, y) \right). \qquad (26c)$$

As a consequence of the problem being non-linear, the cell problems for $\chi_\alpha^{0,1,2}$, $\chi_\alpha^{1,0,2}$ and $\theta_\alpha^{1,0,2}$ cannot be computed once and for all. Indeed, their source term depends on

$$I_m \left( \left( \chi_i(y) - \chi_e(y) \right) \cdot \nabla u^0(x) + v_m^{1,\#}(x) \right)$$

and therefore on the macroscopic variable $x$. Consequently, it must be solved on each point of the macroscopic domain. We loop over every point $x_j \in \mathbb{R}^2$ of the macroscopic mesh and compute $\chi_\alpha^{0,1,2}(x_j, \cdot)$, $\chi_\alpha^{1,0,2}(x_j, \cdot)$ and $\theta_\alpha^{1,0,2}(x_j, \cdot)$. In practice, for a given $\varepsilon$, we only store $\chi_\alpha^{0,1,2}(x_j, \frac{x_j}{\varepsilon})$, $\chi_\alpha^{1,0,2}(x_j, \frac{x_j}{\varepsilon})$ and $\theta_\alpha^{1,0,2}(x_j, \frac{x_j}{\varepsilon})$ and not the whole functions.

Although solving these corrector problems at each macroscopic point is computationally demanding, it still requires a much coarser mesh than what would be needed to solve the full multiscale problem. In practice, depending on the quantity of interest, it may not even be necessary to compute all these correctors, which further reduces the computational cost.

**Membrane conductivity** As often in electroporation modeling, either at tissue [58] or cell scale [34], $S_m$ is modeled as a sigmoid function. We choose the following profile:

$$S_m(\lambda) = S_m^0 + \frac{S_m^1 - S_m^0}{1 + \exp\left(-k\left(\frac{\lambda^2}{V_{ep}^2} - 1\right)\right)},$$

with the following physical parameters:

$$\sigma_i = 2, \sigma_e = 1, S_m^0 = 0.01, S_m^1 = 0.5, k = 5, V_{ep} = 0.1.$$
<figure><img src="image_4.png" /><figcaption>(a) Isotropic cell</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Ellipsoid cell</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: Profile of the normalized source term $f/A$. Contour lines are those of $|f/A|$.</figcaption></figure>

**Geometry** The domain is $\Omega = [0, 1]^2$, and, as explained in Section 2.1, is the juxtaposition of $N \in \mathbb{N}$ square cell domains $Y$ in each direction, meaning that $\varepsilon$ takes values $1/N$. We considered two types of cells to study the effect of the cell shape on the convergence and the solution itself, since cardiac muscle cells are very elongated. We considered isotropic cells (Fig. 2a) and ellipses oriented towards the x-axis (Fig. 2b).

**Source term** The source term $f$ is modeled as

$$f(x_1, x_2) = A \left( \exp \left( -\nu \left( (x_1 - 0.5)^2 + (x_2 - 0.35)^2 \right) \right) + \exp \left( -\nu \left( (x_1 - 0.5)^2 + (x_2 - 0.65)^2 \right) \right) \right), \quad (27)$$

where $\nu = 100$ and $A = 10$, which is the adimensioned amplitude of a typical electroporation current $\sim 1$ A. This source term aims to mimic the insertion of two electrodes, the first one centered at point $(0.5, 0.35)$ and the second one at point $(0.5, 0.65)$, see Fig. 3.

### 3.5.2 Convergence study

We conduct a numerical convergence study with respect to $\varepsilon$. For a range of decreasing values of $\varepsilon$, we compute the non linear microscopic problem (1) using a fixed point method and $P^2$ Lagrangian finite elements, and compare its solutions to the homogenized approximations. More precisely, we study the order 0 and order 1 approximations of $(u_i^\varepsilon, u_e^\varepsilon)$, $v_m^\varepsilon$ and $(E_i^\varepsilon, E_e^\varepsilon) := (-\nabla u_i^\varepsilon, -\nabla u_e^\varepsilon)$ given by Equation (26).
Consequently, the following errors will be numerically approximated and monitored:

$$e_{L_2}^{\alpha,0} := \frac{\|u_{\alpha}^{\epsilon}(x) - u^{0}(x)\|_{L^2(\Omega_{\alpha}^{\epsilon})}}{\sqrt{|Y_{\alpha}|} \|u^{0}(x)\|_{L^2(\Omega)}}, \qquad (28a)$$

$$e_{L_2}^{\alpha,1} := \frac{\|u_{\alpha}^{\epsilon}(x) - u^{0}(x) - \epsilon u_{\alpha}^{1}(x, \frac{x}{\epsilon})\|_{L^2(\Omega_{\alpha}^{\epsilon})}}{\sqrt{|Y_{\alpha}|} \|u^{0}(x)\|_{L^2(\Omega)}}, \qquad (28b)$$

$$e_{H^1}^{\alpha,0} := \frac{\|\mathbf{E}_{\alpha}^{\epsilon}(x) - \mathbf{E}_{\alpha}^{0}(x, \frac{x}{\epsilon})\|_{L^2(\Omega_{\alpha}^{\epsilon})}}{\sqrt{|Y_{\alpha}|} \|\nabla u^{0}(x)\|_{L^2(\Omega)}}, \qquad (28c)$$

$$e_{H^1}^{\alpha,1} := \frac{\|\mathbf{E}_{\alpha}^{\epsilon}(x) - \mathbf{E}_{\alpha}^{1,\epsilon}(x, \frac{x}{\epsilon})\|_{L^2(\Omega_{\alpha}^{\epsilon})}}{\sqrt{|Y_{\alpha}|} \|\nabla u^{0}(x)\|_{L^2(\Omega)}}, \qquad (28d)$$

$$e_{L_2}^{v_m} := \left\| (u_i^{\epsilon} - u_e^{\epsilon}) - \epsilon(u_i^{1} - u_e^{1}) \left(x, \frac{x}{\epsilon}\right) \right\|_{L^2(\Gamma_m^{\epsilon})}. \qquad (28e)$$

**Convergence rates.** Numerical convergence rates are shown in Fig. 4. For each cell type, the errors defined in Eq. (28) are computed for four different values of $\epsilon$. A linear regression is performed to estimate the convergence rate. The meshes consist of approximately $10^6$ nodes, with computations for each value of $\epsilon$ taking around half an hour on 32 cores. This highlights the advantage of a homogenized model, which reduces the need for fine meshing, especially in 3D.

The observed convergence rates are as follows: order $\epsilon$ for the zeroth-order approximations of $(u_i^\epsilon, u_e^\epsilon)$ and $(\mathbf{E}_i^\epsilon, \mathbf{E}_e^\epsilon)$, order $\epsilon^{3/2}$ for the first-order approximation of $v_m^\epsilon$ and order $\epsilon^2$ for the first-order approximations of $(u_i^\epsilon, u_e^\epsilon)$ and $(\mathbf{E}_i^\epsilon, \mathbf{E}_e^\epsilon)$. These encouraging results suggest that theoretical convergence estimates can indeed be established. This will be the focus of the next section.

**Comparison between microscopic and homogenized models** Numerical illustrations of the approximate quantities of interest obtained from the microscopic and homogenized models — when the anisotropic cell is aligned with the x-axis and $\epsilon = 1/15$ — are displayed in Figure 5. The corresponding modeling errors for these quantities are reported in Figs. 6 and 7.

### 3.5.3 Characterization of the ablated zone

To conclude this numerical study and as a first step towards potential biological applications, we illustrate the average membrane conductance over the unit cell membrane, defined by

$$S_{m}^{\text{eff},0}(x) = \frac{1}{|\Gamma|} \int_{\Gamma} S_{m}(v_{m}^{1}(x, y)) dy.$$

This quantity may serve as a preliminary marker for characterizing the effectively ablated zone in PFA.

The numerical setup is qualitatively the same as before. The non linearity profile and the physical parameters $\sigma_i$, $\sigma_e$, $S_m^0$, $S_m^1$, $k$ and $V_{ep}$ remain unchanged. The size of the domain and the distance between the electrodes are inspired from [58] to be closer to a realistic experimental setup. Instead of setting A the amplitude of the source term, we set I the total current injected into the domain and choose A accordingly ($\nu = 100$ is set). We compare specific isolines of $S_m^{\text{eff},0}$ for different current values ($I = 0.5$A and $I = 1$A) and for cells with varying privileged directions: aligned along the x-axis, aligned along the y-axis, and isotropic (no privileged direction).
<figure><img src="image_7.png" /><figcaption>Figure 4: Errors defined in Eq. (28), for the order 0 (left) and order 1 (right) homogenized models as a function of $\varepsilon$ in a log-log scale. $m$ is the slope of the linear regression. The values of $\varepsilon$ are (1/6, 1/9, 1/12, 1/15) for the isotropic cell and (1/12, 1/15, 1/18, 1/21) for the anisotropic cell.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: Solution of the microscopic problem with $\epsilon = 1/15$ and an anisotropic cell oriented towards the x-axis. From left to right, the electric potential $u^\epsilon = u_i^\epsilon 1_{\Omega_i^\epsilon} + u_e^\epsilon 1_{\Omega_e^\epsilon}$, the amplitude of the electric field $|E^\epsilon| = |E_i^\epsilon| 1_{\Omega_i^\epsilon} + |E_e^\epsilon| 1_{\Omega_e^\epsilon}$, the transmembrane potential $v_m^\epsilon$ and the membrane conductance $S_m(v_m^\epsilon/\epsilon)$. All the quantities are unitless.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Difference between the electric fields computed with the microscopic model ($E^\epsilon$) and the homogenized models at order 0 ($E^0$) and order 1 ($E^1$) for $\epsilon = 1/15$ and the anisotropic cell. Note the 6-fold difference between the errors at order 0 and order 1.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Approximation of the transmembrane potential and the resulting cell membrane conductance for $\epsilon = 1/15$. All quantities are unitless. Recall that $S_m$ ranges from $S_m^0 = 0.01$ to $S_m^1 = 0.5$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: Average membrane conductance over the unit cell membrane ranging from $S_m^0 = 0.01$ (no cell membrane electroporated) to $S_m^1 = 0.5$ (all the cell membrane is electroporated). Recall that all quantities are unitless. Delivered current is $I = 0.5$A (top), $I = 1$A (bottom). Three cell types are considered, displayed below the results, representing fibers oriented towards $x$ (left), towards $y$ (right) and no privileged direction (middle). Isolines are 0.1 (light pink), 0.2 (blue), 0.3 (light green).</figcaption></figure>
The size of the cells was set to have similar intra/extra volume ratios, and similar membrane length. This is why the number of connected components, which actually represent the cells, may vary within the cellular domain. Figure 8 presents these unit cells along with the results of the comparison.

The orientation of the fibers, as determined by the cell geometry, has a non-negligible impact on both the shape and the area of the ablated zone. In both anisotropic configurations, the area of the ablated region is relatively similar, albeit for different reasons. When the cells are aligned along the x-axis, a large portion of the cell membrane exhibits conductance values above $S_m^0$. In contrast, when the cells are aligned along the y-axis, only a small portion of the membrane reaches that threshold, but the conductance values are significantly higher, resulting in a similar average $S_m^{\text{eff},0}$. This suggests that additional quantities may need to be considered - such as the variance of the membrane conductance across the membrane, or the root mean square of the membrane conductance - in order to better understand and characterize the ablation behavior.

# 4 Rigorous justification of the full asymptotic expansion

The goal of this section is to provide a rigorous mathematical justification of the asymptotic model introduced earlier. While the derivation was formal, the numerical results presented in the previous section strongly support the validity of the expansion: the observed convergence rates align well with the expected orders in $\epsilon$. This motivates a detailed analysis to establish precise error bounds and confirm that the homogenized system, together with its higher-order corrections, accurately approximates the microscopic model in the asymptotic limit.

We now discuss the strategy adopted in this section. In the variational formulation (13), we formally have

$$ \epsilon \int_{\Gamma_m^\epsilon} I_m \left( \frac{u_i^\epsilon - u_e^\epsilon}{\epsilon} \right) (\phi_i - \phi_e) \to \int_{\Omega} \int_{\Gamma} I_m (u_i^1 - u_e^1) (\phi_i - \phi_e). $$

Using $\phi_i = \phi_i(x)$ and $\phi_e = \phi_e(x)$, and making $\epsilon \to 0$ in Eq. (13), we find

$$ \begin{aligned} \nabla \cdot (\Sigma_i \nabla u^0) + |Y_i| f &= \int_{\Gamma} I_m (u_i^1 - u_e^1), \quad \Omega \\ \nabla \cdot (\Sigma_e \nabla u^0) + |Y_e| f &= -\int_{\Gamma} I_m (u_i^1 - u_e^1), \quad \Omega. \end{aligned} $$

These equations are consistent with the limit problem (10) (already rigorously justified) and the first correction equation (19). However, rigorously justifying the above convergence would require proving that $u_i^\epsilon - u_e^\epsilon = O(\epsilon)$ (in a suitable norm), which is stronger than what we already have i.e. $u_i^\epsilon - u_e^\epsilon = o(1)$ in $H^1$ norm. It usually requires that $u_i^1 - u_e^1$ is already defined. In fact, the process of proving these improved estimates, via the construction of $u_i^1$ and $u_e^1$ and the study of the error, already gives the wanted result, i.e. the justification of the first-order term in the ansatz, making it unnecessary to rely on two-scale convergence arguments.

For these reasons, we adopt a constructive approach: we explicitly define the higher-order correctors via the formal multi-scale expansion and then derive quantitative error estimates to rigorously justify the asymptotic approximation. Importantly, we extend this justification to the full asymptotic expansion at all orders, since the process of constructing correctors and proving error bounds is analogous at each step.
We recall the ansatz

$$u_i^\varepsilon(x) = u_i^0(x) + \sum_{n=1}^\infty \varepsilon^n u_i^n \left(x, \frac{x}{\varepsilon}\right),$$

$$u_e^\varepsilon(x) = u_e^0(x) + \sum_{n=1}^\infty \varepsilon^n u_e^n \left(x, \frac{x}{\varepsilon}\right),$$

and before pursuing the formal multi-scale development initiated in Section 3.3, we give in the following sub-section several mathematical tools that will be required. In the final sub-section, we prove rigorous estimates at any order by studying the system verified by the error.

## 4.1 Mathematical tools

In this subsection, we introduce two types of mathematical tools required for our asymptotic analysis: (1) existence and uniqueness results for PDE systems arising in the formal derivation, and (2) the expansion of the nonlinear term.

**Existence and uniqueness** We present existence and uniqueness results for two prototype systems that arise in our asymptotic analysis. The first system appears in the determination of the transmembrane potential at order one, using the potential at the order zero.

**Proposition 4.1.** Let $B \in C_{\text{per}}^{\infty}(\Omega \times \Gamma)$ and $F \in C_{\text{per}}^{\infty}(\Omega)$, there exists a unique $w \in C_{\text{per}}^{\infty}(\Omega)$ such that

$$\int_{\Gamma} I_{\text{m}}(B(x, y) + w(x)) = F(x), \quad \Omega,$$

where $I_{\text{m}} : \lambda \mapsto S_{\text{m}}(\lambda)\lambda$ verifies Assumption 2.

*Proof.* We identify $\Omega$ to a torus. Let $G : (x, \lambda) \in C^{\infty}(\Omega \times \mathbb{R})$ such that

$$G(x, \lambda) = \int_{\Gamma} I_{\text{m}}(B(x, y) + \lambda) - F(x).$$

Let $x_0 \in \Omega$, the function $\lambda \in \mathbb{R} \mapsto G(x_0, \lambda)$ is in $C^{\infty}(\mathbb{R})$, is strictly increasing (see Assumption 2 and Theorem 2.2) and tends towards $\pm\infty$ as $\lambda \to \pm\infty$. Therefore, there exists a unique $\lambda$, denoted by $w(x_0)$ such that $G(x_0, \lambda) = 0$. Moreover,

$$\partial_{\lambda} G(x_0, w(x_0)) \neq 0.$$

So, by the implicit function theorem, there exists a neighborhood $U$ of $(x_0, w(x_0))$, a neighborhood $V$ of $x_0$ and a function $\psi \in C^{\infty}(V, \Omega)$ such that for all $(x, y) \in \Omega \times \mathbb{R}$,

$$((x, y) \in U \text{ and } G(x, y) = 0) \iff (x \in V \text{ and } y = \psi(x)).$$

Therefore, for any $x \in V$, $G(x, \psi(x)) = 0$ and thus, by unicity, $\psi(x) = w(x)$. It is true in the neighborhood of any point $x_0 \in \Omega$, hence $w \in C^{\infty}(\Omega)$. $\square$

The approximation at order $n$ of $u_e^\varepsilon$ includes the solution of a partial differential equation that has the same structure as the limit problem (10), as described in the following proposition.
**Proposition 4.2.** Let $n \in \mathbb{N}$, let $f_n \in C_{\text{per}}^{\infty}(\Omega)$ such that $\int_{\Omega} f_n = 0$. There exists a unique $u_{e}^{n, \sharp} \in C_{\text{per}}^{\infty}(\Omega)$ such that

$$ -\nabla \cdot ((\Sigma_i + \Sigma_e) \nabla u_e^{n, \sharp}) = f_n, \quad \Omega, \qquad (29a) $$

$$ \int_{\Omega} u_{e}^{n, \sharp} = 0, \qquad (29b) $$

where $\Sigma_i$ and $\Sigma_e$ are given by (11). It verifies

$$ \|u_{e}^{n, \sharp}\|_{H^{1}(\Omega)} \leq C \|f_{n}\|_{L^{2}(\Omega)}. $$

*Proof.* Existence and uniqueness of a solution in $H^1(\Omega)$, along with the continuity relatively to the source term, are obtained by classic arguments. To get more smoothness, we identify $\Omega$ to a torus. As $f_n \in C^\infty(\Omega)$, we get $u_e^{n, \sharp} \in C^\infty(\Omega)$, using for instance Theorem 2.5.1.1 of [27]. □

**Expansion of the non-linear term** Because $v_m^0 = 0$, we expand the non-linear term around $\varepsilon v_m^1$,

**Lemma-Definition 4.3.** Let $N \in \mathbb{N}$, let $v_m^{\varepsilon, N} = \sum_{k=1}^{N} \varepsilon^k v_m^k \in \mathbb{R}$, the following series expansion holds

$$ I_m \left( \frac{v_m^{\varepsilon, N}}{\varepsilon} \right) = I_m(v_m^1) + \sum_{p=2}^{N} \left( \varepsilon^{p-1} I_m^{(1)}(v_m^1) v_m^p + F_p \left( \left( v_m^i \right)_{1 \le i \le p-1} \right) \right) + O(\varepsilon^N), $$

where

$$ F_p \left( \left( v_m^i \right)_{1 \le i \le p-1} \right) = \sum_{\substack{k \text{ s.t.} \\ \sum ik_i = p-1, \quad |\vec{k}|>1}} \frac{1}{|\vec{k}|} \left( \left( \vec{k} \right) I_m^{|\vec{k}|} (v_m^1) \prod_{i=1}^{+\infty} \left( v_m^{i+1} \right)^{k_i} \right). $$

*Proof.* It involves a sequence of computations derived from the Taylor-Lagrange expansion. A similar case is developed in [44]. □

**Remark 4.4.** Note that $F_0 = F_1 = F_2 = 0$.

The term $F_n\left(\left(v_m^k\right)_{1 \le k \le n-1}\right)$ appears in the two-scale expansion at order $n$. Handling it explicitly, as a combination of $v_m^k$, adds a lot of complexity in the notations and computations as the $n$-th term of the expansion, defined in Equation (37), would contain combinations of products of derivatives of lower order terms. Instead, we treat it as a source term accounting for the non-linear behaviour of $I_m$. The cumulative effect of these "source terms" $F_n$ appears in the expression of $(u_i^n, u_e^n)$ as a term $(G_i^n, G_e^n)$ (null if $n \le 2$) that is defined by induction in Theorem 4.5 and could be made explicit with the same reasoning developed in this section.

## 4.2 Development at any order

The cascade of equations writes, for $n \ge 2$,

$$ -\nabla_y \cdot \left( \sigma_i \left( \nabla_y u_i^{n+1} + \nabla_x u_i^n \right) \right) = \nabla_x \cdot \left( \sigma_i \left( \nabla_y u_i^n + \nabla_x u_i^{n-1} \right) \right), \quad \Omega \times Y_i, \qquad (30a) $$

$$ -\nabla_y \cdot \left( \sigma_e \left( \nabla_y u_e^{n+1} + \nabla_x u_e^n \right) \right) = \nabla_x \cdot \left( \sigma_e \left( \nabla_y u_e^n + \nabla_x u_e^{n-1} \right) \right), \quad \Omega \times Y_e, \qquad (30b) $$

$$ \sigma_i \left( \nabla_y u_i^{n+1} + \nabla_x u_i^n \right) \cdot n_i = \sigma_e \left( \nabla_y u_e^{n+1} + \nabla_x u_e^n \right) \cdot n_i, \quad \Omega \times \Gamma, \qquad (30c) $$

$$ -\sigma_i \left( \nabla_y u_i^{n+1} + \nabla_x u_i^n \right) \cdot n_i = I'_m(v_m^1) (u_i^n - u_e^n) + F_n, \quad \Omega \times \Gamma, \qquad (30d) $$

$$ + x \text{ and } y \text{ periodic conditions,} \qquad (30e) $$
where, for readability, we write $F_n$ instead of $F_n((v_m^k)_{1\le k\le n-1})$. Once again, seeing it as a system verified by $y \mapsto u_i^{n+1}(x, y)$ and $y \mapsto u_e^{n+1}(x, y)$ for all $x \in \Omega$, it has a unique solution (up to a constant depending on $x$) if and only if the following compatibility conditions are verified

$$ \int_{Y_i} \nabla_x \cdot \sigma_i (\nabla_y u_i^n + \nabla_x u_i^{n-1}) = \int_{\Gamma} (I'_m(v_m^1)(u_i^n - u_e^n) + F_n), \quad \Omega, \quad (31a) $$

$$ \int_{Y_e} \nabla_x \cdot \sigma_e (\nabla_y u_e^n + \nabla_x u_e^{n-1}) = - \int_{\Gamma} (I'_m(v_m^1)(u_i^n - u_e^n) + F_n), \quad \Omega. \quad (31b) $$

As explained before, for the sake of simplicity and readability, the non-linear term is not handled explicitly. Instead, we prove that we can define $((G_i^n, G_e^n))_{n \in \mathbb{N}}$ as in the following proposition.

**Proposition-Definition 4.5.** Given $(v_m^k)_{n \in \mathbb{N}} \in C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(\Gamma))^{\mathbb{N}}$, there exists for all $n \in \mathbb{N}$, $(G_i^n, G_e^n) \in C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_i)) \times \overline{C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_e))}$ such that

$$ -\nabla_y \cdot (\sigma_i (\nabla_y G_i^{n+1} + \nabla_x G_i^n)) = \nabla_x \cdot (\sigma_i (\nabla_y G_i^n + \nabla_x G_i^{n-1})) \quad \Omega \times Y_i \quad (32a) $$

$$ -\nabla_y \cdot (\sigma_e (\nabla_y G_e^{n+1} + \nabla_x G_e^n)) = \nabla_x \cdot (\sigma_e (\nabla_y G_e^n + \nabla_x G_e^{n-1})) \quad \Omega \times Y_e \quad (32b) $$

$$ \sigma_i (\nabla_y G_i^{n+1} + \nabla_x G_i^n) \cdot n_i = \sigma_e (\nabla_y G_e^{n+1} + \nabla_x G_e^n) \cdot n_i \quad \Omega \times \Gamma \quad (32c) $$

$$ -\sigma_i (\nabla_y G_i^{n+1} + \nabla_x G_i^n) \cdot n_i = I'_m(v_m^1)(G_i^n - G_e^n) + F_n \quad \Omega \times \Gamma, \quad (32d) $$

and

$$ \int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y G_i^n + \nabla_x G_i^{n-1})) = \int_{\Gamma} I'_m(v_m^1)(G_i^n - G_e^n) + F_n, \quad \Omega, \quad (33) $$

$$ \int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y G_i^n + \nabla_x G_i^{n-1})) \quad (34a) $$

$$ + \int_{Y_e} \nabla_x \cdot (\sigma_e (\nabla_y G_e^n + \nabla_x G_e^{n-1})) = 0, \quad \Omega, \quad (34b) $$

where $G_i^n = G_e^n = 0$ for $n \le 2$.

**Remark 4.6.** Let $N \in \mathbb{N}$, note that we only need the knowledge of $(v_m^k)_{1\le k\le N-1}$ to define uniquely $G_\alpha^n$ for all $0 \le n \le N$ and $G_\alpha^{N+1}$ up to a constant depending on $x$.

*Proof.* We show by induction, for $N \ge 2$, that:

* for all $0 \le n \le N-1$, there exist $G_i^n \in C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_i))$ and $G_e^n \in \overline{C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_e))}$ verifying Equations (32), (33) and (34),

* there exist $\tilde{G}_i^N \in C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_i))$ and $\tilde{G}_e^N \in \overline{C_{\text{per}}^{\infty}(\Omega; C_{\text{per}}^{\infty}(Y_e))}$ such that $\int_{Y_i} \tilde{G}_i^N = \int_{Y_e} \tilde{G}_e^N = 0$ and Equation (32) holds with $n = N-1$ and $G_\alpha^{n+1}$ replaced by $\tilde{G}_\alpha^N$. (Conceptually, $\tilde{G}_\alpha^N$ is the part of $G_\alpha^N$ with null average over the $y$ variable.)

The induction is initialized with $G_\alpha^0 = G_\alpha^1 = G_\alpha^2 = 0$ because $F_0 = F_1 = F_2 = 0$. We now describe how to go from rank $N$ to rank $N+1$. Let

$$ v_G^N := \frac{1}{\int_{\Gamma} I'_m(v_m^1)} \left( \begin{aligned} & \int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y \tilde{G}_i^N + \nabla_x G_i^{N-1})) - F_{N-1} \\ & - \int_{\Gamma} I'_m(v_m^1) (\tilde{G}_i^N - \tilde{G}_e^N) \end{aligned} \right). $$
Conceptually, $(G_i^N - G_e^N)(x, y) = (\tilde{G}_i^N - \tilde{G}_e^N)(x, y) + v_G^N(x)$. Let $(v_i, v_e)$ the solution in $C_{\text{per}}^{\infty}(\Omega; \bar{C}_{\text{per}}^{\infty}(Y_i)) \times C_{\text{per}}^{\infty}(\Omega; \bar{C}_{\text{per}}^{\infty}(Y_e))$ of

$$ 
\begin{aligned}
-\nabla_y \cdot (\sigma_i (\nabla_y v_i + \nabla_x \tilde{G}_i^N)) &= \nabla_x \cdot (\sigma_i (\nabla_y \tilde{G}_i^N + \nabla_x G_i^{N-1})) & \Omega \times Y_i \\
-\nabla_y \cdot (\sigma_e (\nabla_y v_e + \nabla_x \tilde{G}_e^N)) &= \nabla_x \cdot (\sigma_e (\nabla_y \tilde{G}_e^N + \nabla_x G_e^{N-1})) & \Omega \times Y_e \\
\sigma_i (\nabla_y v_i + \nabla_x \tilde{G}_i^N) \cdot n_i &= \sigma_e (\nabla_y v_e + \nabla_x \tilde{G}_e^N) \cdot n_i & \Omega \times \Gamma \\
-\sigma_i (\nabla_y v_i + \nabla_x \tilde{G}_i^N) \cdot n_i &= I'_m(v_m^1)(\tilde{G}_i^N - \tilde{G}_e^N + v_G^N) + F_n & \Omega \times \Gamma,
\end{aligned}
$$

where the variable $x$ is seen as a parameter. Conceptually, $v_\alpha = \tilde{G}_\alpha^{N+1} - \chi_\alpha \cdot \nabla \overline{G_\alpha}^N$. This system has a solution because its compatibility condition is verified using the definitions of $\tilde{G}_\alpha^N$ and $v_G^N$. The regularity of $(v_i, v_e)$ in the $x$-variable comes from the $x$-regularity of the source terms and the use of error estimates. Let us define $\overline{G}_e^N$, and thus $\overline{G}_i^N := \overline{G}_e^N + v_G^N$, as the solution in $\bar{C}_{\text{per}}^{\infty}(\Omega)$ of

$$ -\nabla \cdot ((\Sigma_i + \Sigma_e) \cdot \nabla \overline{G_e}^N) - \nabla \cdot (\Sigma_i \cdot \nabla v_G^N) = f_N, \quad \Omega, $$

where

$$ f_N = \int_{Y_i} \nabla_x \cdot (\sigma_i (\nabla_y v_i + \nabla_x \tilde{G}_i^n)) + \int_{Y_e} \nabla_x \cdot (\sigma_e (\nabla_y v_e + \nabla_x \tilde{G}_e^n)). $$

The existence of a solution is ensured by $f_N$ and $\nabla \cdot (\Sigma_i \cdot \nabla v_G^N)$ having a null average over $\Omega$ as divergences of $\Omega$-periodic functions. To conclude the induction, let us define

$$ 
\begin{aligned}
G_{\alpha}^{N} &= \tilde{G}_{\alpha}^{N} + \overline{G}_{\alpha}^{N}, \\
\tilde{G}_{\alpha}^{N+1} &= \chi_{\alpha} \cdot \nabla \overline{G}_{\alpha}^{N} + v_{\alpha}.
\end{aligned}
$$

Computations show that indeed these quantities verify the rank $N + 1$ induction step. $\square$

We previously defined the first and second order cell problems. In the same spirit, we define the cell problems at order $n \ge 3$.

**Definition 4.7** (Cell problems at order $n$). Let us denote $\mathcal{P}_\alpha(F, G)$ the unique solution in $C_{\text{per}}^\infty(Y_\alpha)$ of System (8) with source terms $F$ and $G$. The subscript $\beta$ is $e$, resp. $i$, if $\alpha$ is $i$, resp. $e$. Let $n \ge 3$, let $p \in \mathbb{N}$ and $k \in \mathbb{N}$. If $p + k > n$, then $\chi_\alpha^{p,k,n} = 0$. Otherwise, for $(l_1, \dots, l_p) \in [1; d]^p$,

$$ 
\begin{aligned}
(\chi_{\alpha}^{p,k,n})_{(l_1, \dots, l_p)} &= \delta_{p \ge 1} \mathcal{P}_{\alpha} \left( 2 \frac{\partial \chi_1}{\partial y_{l_p}} - \frac{1}{|Y_{\alpha}|} \int_{Y_{\alpha}} \frac{\partial \chi_1}{\partial y_{l_p}}, \chi_1(n_{\alpha})_{l_p} \right) \\
&\quad + \delta_{p \ge 2} \mathcal{P}_{\alpha} \left( \chi_2 - \frac{1}{|Y_{\alpha}|} \int_{Y_{\alpha}} \chi_2, 0 \right) \delta_{l_{p-1}=l_p} \\
&\quad + \frac{1}{\sigma_{\alpha}} \mathcal{P}_{\alpha} \left( \frac{1}{|Y_{\alpha}|} \int_{\Gamma} I'_{m}(v_{m}^{1})\chi_{3}, I'_{m}(v_{m}^{1})\chi_{3} \right) + \delta_{p=0, k=n},
\end{aligned}
$$

where

$$ \chi_1 = (\chi_{\alpha}^{p-1,k,n-1})_{(l_1, \dots, l_{p-1})}, $$

$$ \chi_2 = (\chi_{\alpha}^{p-2,k,n-2})_{(l_1, \dots, l_{p-2})}, $$

$$ \chi_3 = (\chi_{\alpha}^{p,k,n-1} - \theta_{\beta}^{p,k,n-1})_{(l_1, \dots, l_p)}. $$
If $p+k \ge n$, then $\theta_{\alpha}^{p,k,n} = 0$. Otherwise, for $(l_1, \dots, l_p) \in [1; d]^p$,

$$ \begin{aligned} \left(\theta_{\alpha}^{p,k,n}\right)_{(l_1, \dots, l_p)} ={}& \delta_{p \ge 1} \mathcal{P}_{\alpha} \left( 2 \frac{\partial \theta_1}{\partial y_{l_p}} - \frac{1}{|Y_{\alpha}|} \int_{Y_{\alpha}} \frac{\partial \theta_1}{\partial y_{l_p}}, \theta_1 (n_{\alpha})_{l_p} \right) \\ &+ \delta_{p \ge 2} \mathcal{P}_{\alpha} \left( \theta_2 - \frac{1}{|Y_{\alpha}|} \int_{Y_{\alpha}} \theta_2, 0 \right) \delta_{l_{p-1} = l_p} \\ &+ \frac{1}{\sigma_{\alpha}} \mathcal{P}_{\alpha} \left( \frac{1}{|Y_{\alpha}|} \int_{\Gamma} I'_{\mathrm{m}}(v_m^1) \theta_3, I'_{\mathrm{m}}(v_m^1) \theta_3 \right), \end{aligned} $$

where

$$ \theta_1 = (\theta_{\alpha}^{p-1,k,n-1})_{(l_1, \dots, l_{p-1})}, $$

$$ \theta_2 = (\theta_{\alpha}^{p-2,k,n-2})_{(l_1, \dots, l_{p-2})}, $$

$$ \theta_3 = (\theta_{\alpha}^{p,k,n-1} - \chi_{\beta}^{p,k,n-1})_{(l_1, \dots, l_p)}. $$

Remark 4.8. By construction,

$$ \int_{Y_{\alpha}} \chi_{\alpha}^{p,k,n} = \int_{Y_{\alpha}} \theta_{\alpha}^{p,k,n} = 0, $$

except for $p=0$ and $k=n$, where

$$ \int_{Y_{\alpha}} \chi_{\alpha}^{0,n,n} = |Y_{\alpha}|. $$

We can now define by induction the terms $(u_i^n, u_e^n)$. At each step, the macroscopic quantities $u_i^{n,\sharp} - u_e^{n,\sharp}$ and $u_e^{n,\sharp}$ are defined such that the compatibility conditions (31) are verified, then the solution of System (30) is written as a linear combination of derivatives of the non-oscillating terms $(u_i^{k,\sharp}, u_e^{k,\sharp})$ weighted by the correctors that we previously defined. The first terms $(u_i^0, u_e^0)$ and $(u_i^1, u_e^1)$ were defined previously, we deal with the terms $n \ge 2$ by induction in the following Proposition-Definition:

**Proposition-Definition 4.9.** Let $n \ge 2$. Assume that for all $0 \le k \le n-1$, $u_i^k \in C_{\mathrm{per}}^{\infty}(\Omega \times Y_i)$ and $u_e^k \in C_{\mathrm{per}}^{\infty}(\Omega \times Y_e)$ are defined. We proceed in four steps:

1. With $(v^k)_{0 \le k \le n-1} = (v_m^k)_{0 \le k \le n-1}$, we can define uniquely $G_{i,e}^k$ for $0 \le k \le n$ and $G_{i,e}^{n+1}$ up to a constant, depending on $x$, thanks to Theorem 4.5.

2. We define $v_m^{n,\sharp} \in C_{\mathrm{per}}^{\infty}(\Omega)$ such that

$$ \left( \int_{\Gamma} I'_{\mathrm{m}}(v_m^1) \right) v_m^{n,\sharp} = g_n, \qquad (35) $$

where

$$ \begin{aligned} g_n ={}& \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_x u_i^{n-1} - \int_{\Gamma} F_n - \int_{\Gamma} I'_{\mathrm{m}}(v_m^1) (G_i^n - G_e^n) \\ &+ \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_y \left( \chi_i^{p,k,n} \nabla^p u_i^{k,\sharp} + \theta_i^{p,k,n} \nabla^p u_e^{k,\sharp} \right) + \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_y G_i^n \\ &- \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \int_{\Gamma} I'_{\mathrm{m}}(v_m^1) \left( \chi_i^{p,k,n} \nabla^p u_i^{k,\sharp} + \theta_i^{p,k,n} \nabla^p u_e^{k,\sharp} - \chi_e^{p,k,n} \nabla^p u_e^{k,\sharp} - \theta_e^{p,k,n} \nabla^p u_i^{k,\sharp} \right). \end{aligned} $$

It is well defined because $I'_{\mathrm{m}}(v_m^1) > 0$ (see Theorem 2.2) and because $g_n \in C_{\mathrm{per}}^{\infty}(\Omega)$, as a sum of terms belonging to $C_{\mathrm{per}}^{\infty}(\Omega)$.
3. We then define $u_e^{n,\#} \in C_{\text{per}}^{\infty}(\Omega)$ such that

$$ -\nabla \cdot ((\Sigma_i + \Sigma_e) \nabla u_e^{n,\#}) = f_n - \frac{1}{|\Omega|} \int_{\Omega} f_n, \quad \Omega, \quad (36a) $$

$$ \int_{\Omega} u_e^{n,\#} = 0, \quad (36b) $$

where

$$ \begin{aligned} f_n = \nabla_x \cdot \left( \int_{\Gamma} I'_m(v_m^1) (\chi_i - \chi_e) \right) v_m^{n,\#} + \nabla_x \cdot \Sigma_i \nabla_x v_m^{n,\#} + \delta_{n=0} |Y_e| f \\ + \sum_{k=0}^{n-1} \sum_{p=0}^{n-k+1} \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_y \left( \chi_i^{p,k,n+1} \nabla^p u_i^{k,\#} + \theta_i^{p,k,n+1} \nabla^p u_e^{k,\#} \right) + \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_y G_i^{n+1} \\ + \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_x \left( \chi_i^{p,k,n} \nabla^p u_i^{k,\#} + \theta_i^{p,k,n} \nabla^p u_e^{k,\#} \right) + \int_{Y_i} \nabla_x \cdot \sigma_i \nabla_x G_i^n \\ + \sum_{k=0}^{n-1} \sum_{p=0}^{n-k+1} \int_{Y_e} \nabla_x \cdot \sigma_e \nabla_y \left( \chi_e^{p,k,n+1} \nabla^p u_e^{k,\#} + \theta_e^{p,k,n+1} \nabla^p u_i^{k,\#} \right) + \int_{Y_e} \nabla_x \cdot \sigma_e \nabla_y G_e^{n+1} \\ + \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \int_{Y_e} \nabla_x \cdot \sigma_e \nabla_x \left( \chi_e^{p,k,n} \nabla^p u_e^{k,\#} + \theta_e^{p,k,n} \nabla^p u_i^{k,\#} \right) + \int_{Y_e} \nabla_x \cdot \sigma_e \nabla_y G_e^n. \end{aligned} $$

It is well defined by use of Theorem 4.2, and noting that $f_n \in C_{\text{per}}^{\infty}$ and $\int_{\Omega} f_n = 0$ ($f_n$ is a periodic divergence).

4. Finally, we define

$$ \begin{aligned} u_i^n(x, y) = & \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \left( \chi_i^{p,k,n}(x, y) \nabla^p u_i^{k,\#}(x) + \theta_i^{p,k,n}(x, y) \nabla^p u_e^{k,\#}(x) \right) \\ & + G_i^n(x, y) + u_i^{n,\#}(x), \end{aligned} \quad (37) $$

$$ \begin{aligned} u_e^n(x, y) = & \sum_{k=0}^{n-1} \sum_{p=0}^{n-k} \left( \chi_e^{p,k,n}(x, y) \nabla^p u_e^{k,\#}(x) + \theta_e^{p,k,n}(x, y) \nabla^p u_i^{k,\#}(x) \right) \\ & + G_e^n(x, y) + u_e^{n,\#}(x). \end{aligned} \quad (38) $$

Remark 4.10. The definition of $v_m^{n,\#}$ (35) reads

$$ \int_{Y_i} \nabla_x \cdot \sigma_i (\nabla_y u_i^n + \nabla_x u_i^{n-1}) = \int_{\Gamma} (I'_m (u_i^n - u_e^n) + F_n), $$

and the definition of $u_e^{n,\#}$ (36) reads

$$ \int_{Y_i} \nabla_x \cdot \sigma_i (\nabla_y u_i^{n+1} + \nabla_x u_i^n) + \int_{Y_e} \nabla_x \cdot \sigma_e (\nabla_y u_e^{n+1} + \nabla_x u_e^n) = 0, $$

which are the compatibility conditions (31).

The cell correctors are defined such that the following proposition holds. We do not give the proof which consists in heavy computations.

**Proposition 4.11.** For all $n \ge 2$, $u_i^n$ and $u_e^n$, as defined in Definition 4.9, are solutions of System (30).
## 4.3 Error estimates

The formal two-scale asymptotic expansion led us towards a definition of $u_i^n$, $u_e^n$. Now comes the time to rigorously prove that this definition is the right one, in the sense that it is consistent with the ansatz. We therefore prove that the finite development $\sum_{k=0}^n \varepsilon^k u_\alpha^k$ is an approximation of $u_\alpha^\varepsilon$ at order $\varepsilon^{n+1}$.

Note that in the previous sub-section, we approximated the gauge condition $\int_{\Omega_\varepsilon} u_\varepsilon = 0$ by $\int_{\Omega} \int_{Y_e} u_e^n = 0$ for all $n \in \mathbb{N}$. Theorem B.3, given with its proof in Appendix B confirms that it is accurate enough, as it shows that the difference between these two terms can be controlled by $C(n, f)\varepsilon^n$ for all $n \in \mathbb{N}$ where $C(n, f)$ is a constant dependent only on $f$ and $n$.

**Proposition 4.12 (Error estimates).** Under Assumptions 1 and 2, let $(u_i^\varepsilon, u_e^\varepsilon) \in C_{\text{per}}^\infty(\Omega_i^\varepsilon) \times \overline{C_{\text{per}}^\infty}(\Omega_e^\varepsilon)$ the sequence of solutions of System (1) and let $(u_i^k, u_e^k)$ defined as in Theorem 4.9. For all $n \in \mathbb{N}$ and $\varepsilon > 0$ small enough, there exists $C$ independent of $\varepsilon$, but depending on $n$, such that the quantities

$$ 
\begin{aligned}
z_i^{\varepsilon,n}(x) &= u_i^{\varepsilon}(x) - \left( u_i^0(x) + \sum_{k=1}^n \varepsilon^k u_i^k(x, \frac{x}{\varepsilon}) \right), \\
z_e^{\varepsilon,n}(x) &= u_e^{\varepsilon}(x) - \left( u_e^0(x) + \sum_{k=1}^n \varepsilon^k u_e^k(x, \frac{x}{\varepsilon}) \right),
\end{aligned}
$$

verify

$$ 
\begin{aligned}
\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^\varepsilon)} + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^\varepsilon)} &\le C\varepsilon^n, \\
\|z_i^{\varepsilon,n} - z_e^{\varepsilon,n}\|_{L^2(\Gamma_m^\varepsilon)} &\le C\varepsilon^{n+\frac{1}{2}}.
\end{aligned}
$$

*Proof.* This proof uses Theorems A.2 and B.3 given with their proofs in Appendices A and B, respectively.

It relies on firstly finding the system verified by the errors $z_i^{\varepsilon,n}$ and $z_e^{\varepsilon,n}$, which has the same structure as the original system, and then deriving estimates to get a bound on the error in some energy norms. The computation of the error system is straightforward and rely on the fact that the correctors $u_i^n$ and $u_e^n$ are defined to verify the cascade of equations found previously. We find

$$ 
\begin{aligned}
-\nabla \cdot \sigma_i \nabla z_i^{\varepsilon,n} &= -\varepsilon^{n-1} \nabla_y \cdot \sigma_i \nabla_y u_i^{n+1} + \varepsilon^n \nabla_x \cdot \sigma_i \nabla_x u_i^n, && \Omega_i^\varepsilon, \\
-\nabla \cdot \sigma_e \nabla z_e^{\varepsilon,n} &= -\varepsilon^{n-1} \nabla_y \cdot \sigma_e \nabla_y u_e^{n+1} + \varepsilon^n \nabla_x \cdot \sigma_e \nabla_x u_e^n, && \Omega_e^\varepsilon, \\
-\sigma_i \nabla z_i^{\varepsilon,n} \cdot n_i &= \varepsilon I_m \left( \frac{u_i^\varepsilon - u_e^\varepsilon}{\varepsilon} \right) - \varepsilon I_m (v_m^1) && \Gamma_m^\varepsilon, \\
&\quad - \sum_{k=2}^n \varepsilon^k \left( I_m'(v_m^1) (u_i^k - u_e^k) + F_k \right) - \varepsilon^n \sigma_i \nabla_y u_i^{n+1} \cdot n_i, && \Gamma_m^\varepsilon, \\
-\sigma_e \nabla z_e^{\varepsilon,n} \cdot n_e &= -\varepsilon I_m \left( \frac{u_i^\varepsilon - u_e^\varepsilon}{\varepsilon} \right) + \varepsilon I_m (v_m^1) && \Gamma_m^\varepsilon, \\
&\quad + \sum_{k=2}^n \varepsilon^k \left( I_m'(v_m^1) (u_i^k - u_e^k) + F_k \right) - \varepsilon^n \sigma_e \nabla_y u_e^{n+1} \cdot n_e, && \Gamma_m^\varepsilon.
\end{aligned}
$$

In order to find estimates, we multiply the volume equations by $z_i^{\varepsilon,n}$ and $z_e^{\varepsilon,n}$, integrate over $\Omega_i^\varepsilon$ and $\Omega_e^\varepsilon$ and sum. We use the fact that if $(x, y) \mapsto \psi(x, y)$ is $y$-periodic and smooth enough then $\int_{\Omega_\varepsilon} \psi(x, \frac{x}{\varepsilon})^2$ is uniformly bounded because it converges to $\int_{\Omega} \int_{Y_\alpha} \psi(x, y)^2$. Together with
the series development given in Theorem 4.3, we can show that the right hand side is bounded by

$$C \varepsilon^n (\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^{\varepsilon})} + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^{\varepsilon})}).$$

Therefore, using Assumption 2, the estimate becomes

$$\int_{\Omega_i^{\varepsilon}} \sigma_i |\nabla z_i^{\varepsilon,n}|^2 + \int_{\Omega_e^{\varepsilon}} \sigma_e |\nabla z_e^{\varepsilon,n}|^2 + \int_{\Gamma_m^{\varepsilon}} S_m^0 |z_i^{\varepsilon,n} - z_e^{\varepsilon,n}|^2 \\ \leq C \varepsilon^n (\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^{\varepsilon})} + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^{\varepsilon})}).$$

Then, using the norm equivalence given in Theorem A.2 (Appendix A) on $(z_i^{\varepsilon,n} - \frac{1}{|\Omega_e^{\varepsilon}|} \int_{\Omega_e^{\varepsilon}} z_e^{\varepsilon,n}, z_e^{\varepsilon,n} - \frac{1}{|\Omega_e^{\varepsilon}|} \int_{\Omega_e^{\varepsilon}} z_e^{\varepsilon,n}$ we get for $\varepsilon$ small enough,

$$\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^{\varepsilon})}^2 + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^{\varepsilon})}^2 \\ \leq C \left( \varepsilon^n (\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^{\varepsilon})} + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^{\varepsilon})}) + \frac{|\Omega|}{|\Omega_e^{\varepsilon}|^2} \left( \int_{\Omega_e^{\varepsilon}} z_e^{\varepsilon,n} \right)^2 \right).$$

Now remark that $\int_{\Omega_e^{\varepsilon}} z_e^{\varepsilon,n} = - \sum_{k=0}^n \int_{\Omega_e^{\varepsilon}} \varepsilon^k u_e^k(x, \frac{x}{\varepsilon}) \leq C \varepsilon^{n+1}$ by using Theorem B.3 (Appendix B), where $C$ depends on $n$. We obtain finally

$$\|z_i^{\varepsilon,n}\|_{H^1(\Omega_i^{\varepsilon})} + \|z_e^{\varepsilon,n}\|_{H^1(\Omega_e^{\varepsilon})} \leq C \varepsilon^n, \\ \|z_i^{\varepsilon,n} - z_e^{\varepsilon,n}\|_{L^2(\Gamma_m^{\varepsilon})} \leq C \varepsilon^n.$$

To improve the last result, note that

$$\|z_i^{\varepsilon,n} - z_e^{\varepsilon,n}\|_{L^2(\Gamma_m^{\varepsilon})} \leq \|z_i^{\varepsilon,n+1} - z_e^{\varepsilon,n+1}\|_{L^2(\Gamma_m^{\varepsilon})} \\ + \varepsilon^{n+1} \left\| u_i^{n+1}\left(x, \frac{x}{\varepsilon}\right) - u_e^{n+1}\left(x, \frac{x}{\varepsilon}\right) \right\|_{L^2(\Gamma_m^{\varepsilon})} \\ \leq C \varepsilon^{n+\frac{1}{2}},$$

because $\sqrt{\varepsilon} \|u_i^{n+1}(x, \frac{x}{\varepsilon}) - u_e^{n+1}(x, \frac{x}{\varepsilon})\|_{L^2(\Gamma_m^{\varepsilon})}$ is bounded. $\square$

Remark 4.13. We have

$$\|z_{\alpha}^{\varepsilon,n}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} \leq \|z_{\alpha}^{\varepsilon,n+1}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} + \varepsilon^{n+1} \|u_{\alpha}^{n+1}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} \leq C \varepsilon^{n+1}$$

and

$$\|\nabla z_{\alpha}^{\varepsilon,n} + \varepsilon^n \nabla_y u_{\alpha}^{n+1}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} \leq \|\nabla z_{\alpha}^{\varepsilon,n+1}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} + \varepsilon^{n+1} \|\nabla_x u_{\alpha}^{n+1}\|_{L^2(\Omega_{\alpha}^{\varepsilon})} \\ \leq C \varepsilon^{n+1}.$$

The last inequality is useful in practice because it allows us to gain one order of convergence without having to solve a new macroscopic problem, indeed $\nabla_y u_{\alpha}^{n+1}(x, y)$ does not contain the term $u_{\alpha}^{n+1,\#}(x)$, it is a combination of cell problems solutions and of the previous macroscopic quantities $u_i^{k,\#}$ and $u_e^{k,\#}$ for $k \leq n$.

Remark 4.14 (On the regularity). We assumed that $f, S_m$ and $\Gamma$ are infinitely smooth, so that the expansion can be done at any order $n$. These assumptions on regularity can obviously be weakened if the expansion is stopped at a fixed order $n_0$. We would then work in Hilbert spaces of order depending on $n_0$.
# 5 Conclusion

In this work, we have developed and rigorously justified a multiscale model for pulsed field ablation in cardiac tissue, based on high-order two-scale asymptotic expansions. The resulting homogenized bidomain system captures key nonlinear effects of electroporation at the macroscopic level. An interesting feature of our model lies in the fact that the gradient of the potential, thus the electric field plays a role in the nonlinearity. This is the main difference with the usual bidomain model used in electrocardiology, in which the nonlinearity involves only the transmembrane voltage. It is also important to note that the electric field magnitude and the macroscopic transmembrane voltage appear in the nonlinearity of the macroscopic effective model, while the nonlinearity of the microscopic model involved only the microscopic transmembrane voltage. From the modeling view point, our result clarifies and quantifies the fact that at the cell scale, electroporation is driven by the increased of the transmembrane voltage while at the tissue scale, the magnitude of the electric field is the key variable. Our analysis ensures the convergence of the expansion at arbitrary order, providing a solid mathematical foundation for the model. Thanks to this high-order asymptotic approach, we now have a reliable and well-structured framework ready to be confronted with biological data, which constitutes the natural next step of this study.

# Acknowledgements

This work was partly funded by the French Federation of Cardiology through the Dielectric project. It also received support the French National Research Agency through several grants (Liryc ANR-10-IAHU-04, MIRE4VTACH ANR-22-CE45-0014, IMITATE ANR-22-CE51-0043), and from Inserm Plan Cancer MECI project (21CM119_00). Simulations were carried out on the computing facilities of the MCIA (Mésocentre de Calcul Intensif Aquitain).

# A Spaces that depend on $\varepsilon$

**Lemma A.1** (Usual inequalities - dependance on $\varepsilon$). There exists a constant $C$ independent of $\varepsilon$ such that for all $v_e \in H^1(\Omega_e^\varepsilon)$ and $v_i \in H^1(\Omega_i^\varepsilon)$, we have

$$ \left\| v_e - \frac{1}{|\Omega_e^\varepsilon|} \int_{\Omega_e^\varepsilon} v_e \right\|_{L^2(\Omega_e^\varepsilon)} \le C \|\nabla v_e\|_{L^2(\Omega_e^\varepsilon)}, \quad (\text{Poincaré-Wirtinger ineq.}), \\ \|v_\alpha\|_{L^2(\Gamma_m^\varepsilon)}^2 \le C\varepsilon^{-1} \|v_\alpha\|_{L^2(\Omega_\alpha^\varepsilon)}^2 + C\varepsilon \|\nabla v_\alpha\|_{L^2(\Omega_\alpha^\varepsilon)}^2, \quad (\text{Trace ineq.}). $$

*Proof.* See appendices of [9]. □

**Lemma A.2** (Norm equivalence). There exist $c > 0$ and $C > 0$ independent of $\varepsilon$ such that for all $H^1(\Omega_i^\varepsilon) \times \overline{H^1}(\Omega_e^\varepsilon)$,

$$ c \|(v_i, v_e)\|_{H^1(\Omega_i^\varepsilon) \times H^1(\Omega_e^\varepsilon)} \le \left( \|\nabla v_i\|_{L^2(\Omega_i^\varepsilon)}^2 + \|\nabla v_e\|_{L^2(\Omega_e^\varepsilon)}^2 + \varepsilon \|v_i - v_e\|_{L^2(\Gamma_m^\varepsilon)}^2 \right)^{\frac{1}{2}} \\ \le C \|(v_i, v_e)\|_{H^1(\Omega_i^\varepsilon) \times H^1(\Omega_e^\varepsilon)}. $$

*Proof.* See Appendix C of [9]. □
## B Approximation of the gauge condition

In Section 4.2, the quantities $u_e^{k,\#}$ (and thus $u_e^k$) can be defined up to a constant. For $\sum_{k=0}^n \varepsilon^k u_e^k$ to be a good approximation of $u_e^\varepsilon$, the average $\sum_{k=0}^n \varepsilon^k \int_{\Omega_e^\varepsilon} u_e^k(x, x/\varepsilon)$ has to be close to $\int_{\Omega_e^\varepsilon} u_e^\varepsilon = 0$, as seen in the proof of Theorem 4.12.

A classic result tells us that

$$ \int_{\Omega_e^\varepsilon} u_e^k \left( x, \frac{x}{\varepsilon} \right) = \int_{\Omega} \int_{Y_e} u_e^k(x, y) + O(\varepsilon). \quad (40) $$

So, for an approximation at order $\varepsilon$, it suffices to set

$$ \int_{\Omega} u_e^{k,\#}(x) = \int_{\Omega} \int_{Y_e} u_e^k(x, y) = 0. $$

But for approximation at order $\varepsilon^{n+1}$, we need to develop the remainder $O(\varepsilon)$ in (40). This is done in Theorem B.1. The new terms still depend on $x/\varepsilon$, therefore, in Theorem B.2, we recursively apply Theorem B.1 to approximate them. Finally, in Theorem B.3, we get the final approximation result, that is, in fact

$$ \int_{\Omega_e^\varepsilon} u_e^k \left( x, \frac{x}{\varepsilon} \right) = \int_{\Omega} \int_{Y_e} u_e^k(x, y) + O(\varepsilon^{n+1}). $$

This result is a generalization of Lemma 1 in [16] and a variant of Lemma 3.1 in [7] for higher order approximation. It is also close to the Theorem in Appendix C of [51].

**Lemma B.1.** Let $(D, E) = (Y, \Omega)$ or $(Y_e, \Omega_e^\varepsilon)$. Let $f: \Omega \times D \to \mathbb{R}$ such that for all $y \in D$, $x \mapsto f(x, y) \in C_{\text{per}}^\infty(\Omega)$ and for all $x \in \Omega$, $y \mapsto f(x, y) \in L^1(D)$ and is $Y$-periodic. For all $n \in \mathbb{N}$, we have

$$ \int_E f \left( x, \frac{x}{\varepsilon} \right) = \int_{\Omega} \int_D f(x, y) + \sum_{k=1}^n \varepsilon^k \int_{\Omega} S_k^D(f) \left( x, \frac{x}{\varepsilon} \right) + \varepsilon^{n+1} R_n^{\varepsilon,D}(f), $$

where

$$ S_k^D(f)(x, y) = \sum_{|\alpha|=k} \frac{1}{\alpha!} \int_D (z - \{y\})^\alpha \partial_x^\alpha f(x, z) dz, \quad (41) $$

where $\{y\}$ is the fractional part of $y$, and

$$ R_n^{\varepsilon,D}(f) = \sum_{|\alpha|=n+1} \frac{n+1}{\alpha!} \int_{\Omega \times D} \left(y - \left\{\frac{x}{\varepsilon}\right\}\right)^{\alpha} \int_{0}^{1} (1-t)^n (\partial_x^{\alpha} f) \left(x + \varepsilon t \left(y - \left\{\frac{x}{\varepsilon}\right\}\right), y\right) d^3 tyx. \quad (42) $$

*Proof.* Let $\varepsilon > 0$, let $m \in N_\varepsilon$. Let $y \in D$ and $z \in Y$, the Taylor development of $f(\cdot, y)$ gives for $n \in \mathbb{N}$,

$$ f(\varepsilon(y+m), y) = f(\varepsilon(z+m), y) + \sum_{k=1}^{n} \varepsilon^k \sum_{|\alpha|=k} \frac{1}{\alpha!} (y-z)^{\alpha} \partial_x^{\alpha} f(\varepsilon(z+m), y) \\ + \varepsilon^{n+1} \sum_{|\alpha|=n+1} \frac{n+1}{\alpha!} (y-z)^{\alpha} \int_{0}^{1} (1-t)^n \partial_x^{\alpha} f(\varepsilon(z+m+t(y-z)), y). $$
Integrate this equation over $Y$ and $D$, and by change of variable, one gets

$$ egin{aligned} |Y| \int_{\varepsilon(D+m)} f \left( x, \frac{x}{\varepsilon} \right) = & \int_{\varepsilon(Y+m)} \int_{D} f(x, y) \\ & + \sum_{k=1}^{n} \varepsilon^k \sum_{|\alpha|=k} \frac{1}{\alpha!} \int_{\varepsilon(Y+m)} \int_{D} \left( y - \left\{ \frac{x}{\varepsilon} \right\} \right)^{\alpha} \partial_x^{\alpha} f(x, y) \\ & + \varepsilon^{n+1} \sum_{|\alpha|=n+1} \frac{n+1}{\alpha!} \int_{\varepsilon(Y+m)} \int_{D} \left( y - \left\{ \frac{x}{\varepsilon} \right\} \right)^{\alpha} \int_{0}^{1} (1-t)^n \partial_x^{\alpha} f \left( x + \varepsilon t \left( y - \left\{ \frac{x}{\varepsilon} \right\} \right), y \right). \end{aligned} $$

Sum over all $m \in N_{\varepsilon}$ and the proof is over. $\square$

The correction $\int_{\Omega} S_k^D(f)(x, x/\varepsilon)$ depends on $\varepsilon$, therefore it also needs to be approximated.

**Lemma B.2.** Under the same hypothesis as Theorem B.1, for all $n \in \mathbb{N}$, we have

$$ egin{aligned} \int_E f \left( x, \frac{x}{\varepsilon} \right) = & \frac{1}{|Y|} \int_{\Omega} \int_D f(x, y) \\ & + \frac{1}{|Y|} \sum_{k=1}^{n} \varepsilon^k \int_{\Omega} \int_Y \sum_{(\beta_j) \in (\mathbb{N}^*)^J \atop |\beta|=k} (S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f)(x, y) \\ & + \varepsilon^{n+1} \left( \sum_{l=0}^{n-1} \sum_{(\beta_j) \in (\mathbb{N}^*)^J \atop |\beta|=n-l} R_l^{\varepsilon, Y}(S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f) + R_n^{\varepsilon, D}(f) \right), \end{aligned} $$

where $S_k^D$ and $R_k^{\varepsilon, D}$ are defined in (41) and (42).

*Proof.* The proof is done by induction. The case $n = 0$ is true because

$$ \int_E f \left( x, \frac{x}{\varepsilon} \right) = \frac{1}{|Y|} \int_{\Omega} \int_D f(x, y) + \varepsilon R_0^{\varepsilon, D}(f). $$

To go from rank $n$ to rank $n + 1$, remark the induction relation (using Theorem B.1 at rank $n$ and $n + 1$)

$$ R_n^{\varepsilon, D}(f) = \int_{\Omega} S_{n+1}^{D}(f) \left( x, \frac{x}{\varepsilon} \right) + \varepsilon R_{n+1}^{\varepsilon, D}(f). $$

Using Theorem B.1 on $S_{n+1}^D(f)$, we get

$$ R_n^{\varepsilon, D}(f) = \frac{1}{|Y|} \int_{\Omega} \int_{Y} S_{n+1}^{D}(f)(x, y) + \varepsilon(R_0^{\varepsilon, Y}(S_{n+1}(f)) + R_{n+1}^{\varepsilon, D}(f)). $$

Therefore

$$ egin{aligned} R_l^{\varepsilon, Y}(S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f) \\ = & \frac{1}{|Y|} \int_{\Omega} \int_{Y} (S_{l+1}^Y \circ S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f)(x, y) \\ & + \varepsilon(R_0^{\varepsilon, Y}(S_{l+1}^Y \circ S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f) \\ & + R_{l+1}^{\varepsilon, Y}(S_{\beta_1}^Y \circ \cdots \circ S_{\beta_{J-1}}^Y \circ S_{\beta_J}^D \circ f)), \end{aligned} $$

and by change of indexes, the proof is over. $\square$
**Lemma B.3.** Under the same hypothesis as Theorem B.1, for all $n \in \mathbb{N}$, we have

$$ \left| \int_E f \left( x, \frac{x}{\varepsilon} \right) - \frac{1}{|Y|} \int_{\Omega} \int_D f(x, y) \right| \leq C(n, f) \varepsilon^{n+1}, $$

where $C(n, f) > 0$ is a constant depending on $n$ and $f$.

*Proof.* Starting from the result of Theorem B.2, we first show that for any $g: \Omega \times D \to \mathbb{R}$ such that for all $y \in D$, $x \mapsto g(x, y) \in C_{\text{per}}^{\infty}(\Omega)$ and for all $x \in \Omega$, $y \mapsto g(x, y)$ is $Y$-periodic, for any $k \geq 1$,

$$ \int_{\Omega} \int_{Y} S_{k}^{D}(g)(x, y) = 0. $$

It suffices to perform an integration by parts relatively to the $x$-variable and recall that $g$ and all its derivatives are $x$-periodic. Finally we bound the remainder, using the fact that for the same kind of $g$ as before, for $l \in \mathbb{N}$,

$$ |R_l^{\varepsilon,D}(g)| \leq \frac{d^{\frac{3(l+1)}{2}}}{(l+1)!} |\Omega||D| \sup_{(x,y) \in \Omega \times D, |\alpha|=l+1} |\partial_x^\alpha g(x, y)| $$

and

$$ |S_l^D(g)(x, y)| \leq \frac{d^{\frac{3l}{2}}}{l!} |D| \sup_{(x', y') \in \Omega \times D, |\alpha|=l} |\partial^\alpha g(x', y')|. $$

□

## References

[1] Pierre Agnass, Eran van Veldhuisen, Martin J. C. van Gemert, Cees W. M. van der Geld, Krijn P. van Lienden, Thomas M. van Gulik, Martijn R. Meijerink, Marc G. Besselink, H. Petra Kok, and Johannes Crezee. Mathematical modeling of the thermal effects of irreversible electroporation for in vitro, in vivo, and clinical use: a systematic review. *International Journal of Hyperthermia*, 37(1):486–505, 2020. PMID: 32423258.

[2] Felipe Aguel, Katherine A. Debrutn, Wanda Krassowska, and Natalia A. Trayanova. Effects of electroporation on the transmembrane potential distribution in a two-dimensional bidomain model of cardiac tissue. *Journal of Cardiovascular Electrophysiology*, 10(5):701–714, 1999.

[3] R. Aliev and A. Panfilov. A simple two-variable model of cardiac excitation. *Chaos, Solitons and Fractals*, 7(3):293–301, 1996.

[4] G. Allaire, A. Damlamian, and U. Hornung. Two-scale convergence on periodic surfaces and applications. *Mathematical Modelling of Flow through Porous Media*, 09 1995.

[5] Grégoire Allaire. Homogenization and two-scale convergence. *SIAM Journal on Mathematical Analysis*, 23(6):1482–1518, 1992.

[6] Grégoire Allaire, Marc Briane, and Muthusamy Vanninathan. A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures. *SeMA Journal: Boletín de la Sociedad Española de Matemática Aplicada*, 73(3):237–259, 2016.
[7] Grégoire Allaire and Zakaria Habibi. Second order corrector in the homogenization of a conductive-radiative heat transfer problem. *Discrete and Continuous Dynamical Systems* - *B*, 18(1):1–36, 2013.

[8] M. Amar, D. Andreucci, P. Bisegna, and R. Gianni. On a hierarchy of models for electrical conduction in biological tissues. *Mathematical Methods in the Applied Sciences*, 29(7):767–787, 2006.

[9] Habib Ammari, Josselin Garnier, Laure Giovangigli, Wenjia Jing, and Jin Keun Seo. Spectroscopic imaging of a dilute cell suspension. *Journal de Mathématiques Pures et Appliquées*, 105(5):603–661, 2016.

[10] Jason G Andrade. Cryoablation for atrial fibrillation. *Heart Rhythm O2*, 1(1):44–58, 2020.

[11] Jean-Louis Auriault and Horia I. Ene. Macroscopic modelling of heat transfer in composites with interfacial thermal barrier. *International Journal of Heat and Mass Transfer*, 37(18):2885–2892, 1994.

[12] Mostafa Bendahmane, Fatima Mroue, Mazen Saad, and Raafat Talhouk. Unfolding homogenization method applied to physiological and phenomenological bidomain models in electrocardiology. *Nonlinear Analysis: Real World Applications*, 50:413–447, 2019.

[13] Guido Caluori, Eva Odehnalova, Tomasz Jadczyk, Martin Pesl, Iveta Pavlova, Lucia Valikova, Steffen Holzinger, Veronika Novotna, Vladimir Rotrekl, Ales Hampl, et al. Ac pulsed field ablation is feasible and safe in atrial and ventricular settings: a proof-of-concept chronic animal study. *Frontiers in Bioengineering and Biotechnology*, 8:552357, 2020.

[14] Quim Castellvi and Antoni Ivorra. Computational multiscale modeling of pulsed field ablation considering conductivity and damage anisotropy reveals deep lesion morphologies. *International Journal for Numerical Methods in Biomedical Engineering*, 41(8):e70077, 2025. e70077 CNM-Feb-25-0196.R1.

[15] Fatemeh Chegini, Algiane Froehly, Ngoc Mai Monica Huynh, Luca F. Pavarino, Mark Potse, Simone Scacchi, and Martin Weiser. Efficient numerical methods for simulating cardiac electrophysiology with cellular resolution. In *COUPLED 2023 - 10th International Conference on Computational Methods for Coupled Problems in Science and Engineering*, Chania, Greece, June 2023.

[16] Kirill D. Cherednichenko and Valery P. Smyshlyaev. On full two-scale expansion of the solutions of nonlinear periodic rapidly oscillating problems and higher-order homogenised variational problems. *Archive for Rational Mechanics and Analysis*, 174:385–442, 2004.

[17] Annabelle Collin and Sébastien Imperiale. Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model. *Mathematical Models and Methods in Applied Sciences*, 28:1–57, 03 2018.

[18] K. DeBruin and W. Krassowska. Modelling electroporation in a single cell. I. Effects of field strength and rest potential. *Biophysical Journal*, 77(3):1213–1224, 1999.

[19] Manon Deville. *Mathematical Modeling of enhanced drug delivery by mean of Electroporation or Enzymatic treatment*. Theses, Université de Bordeaux ; Università degli studi di Roma "Tor Vergata" (1972-....), November 2017.
[20] Manon Deville. *Mathematical Modeling of enhanced drug delivery by mean of Electroporation or Enzymatic treatment*. PhD thesis, Université de Bordeaux ; Università degli studi di Roma "Tor Vergata", November 2017.

[21] Patrizia Donato and Kim Hang Le Nguyen. Homogenization of diffusion problems with a nonlinear interfacial resistance. *Nonlinear Differential Equations and Applications*, 22(5):1345–1380, October 2015.

[22] H. I. Ene and D. Poliševski. Model of diffusion in partially fissured media. *Zeitschrift für angewandte Mathematik und Physik ZAMP*, 53(6):1052–1059, November 2002.

[23] Horia Ene and Claudia Timofte. Microstructure models for composites with imperfect interface via the periodic unfolding method. *Asymptotic Analysis*, 89(1-2):111–122, 2014.

[24] A. Evans. *Partial Differential Equations*. American Mathematical Society, 1998.

[25] E Frénod, P.-A Raviart, and E Sonnendrücker. Two-scale expansion of a singularly perturbed convection equation. *Journal de Mathématiques Pures et Appliquées*, 80(8):815–843, 2001.

[26] Emmanuel Frénod. Application of the averaging method to the gyrokinetic plasma. *Asymptotic Analysis*, 46, 03 2007.

[27] Pierre Grisvard. *Elliptic Problems in Nonsmooth Domains*. Society for Industrial and Applied Mathematics, 2011.

[28] F. Hecht. New development in freefem++. *Journal of Numerical Mathematics*, 20(3-4):251–265, 2012.

[29] A. Hodgkin and A. F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. *Journal of Physiology*, 117(4):500–544, 1952.

[30] A. Ivorra, J. Villemejane, and LM. Mir. Electrical modeling of the influence of medium conductivity on electroporation. *Physical Chemistry Chemical Physics*, 12(34), 2010.

[31] Edward J. Jacobs, Pedro P. Santos, and Rafael V. Davalos. Comprehensive characterization of waveform-dependent cardiac tissue electroporation for pulsed field ablation. *Biosensors and Bioelectronics*, 289:117890, 2025.

[32] Gaspard Jankowiak, Cécile Taing, Clair Poignard, and Annabelle Collin. Comparison and calibration of different electroporation models - Application to rabbit livers experiments. *ESAIM: Proceedings and Surveys*, 67, June 2020.

[33] Otared Kavian. *Introduction à la théorie des points critiques et applications aux problèmes elliptiques*. Springer Berlin, Heidelberg, 1991.

[34] Otared Kavian, Michael Leguèbe, Clair Poignard, and Lisl Weynans. "Classical" Electroporationmeabilization Modeling at the Cell Scale. *Journal of Mathematical Biology*, 68(1-2):30, January 2014.

[35] Jacob Koruth, Kenji Kuroki, Jin Iwasawa, Yoshinari Enomoto, Raju Viswanathan, Richard Brose, Eric D Buck, Molly Speltz, Srinivas R Dukkipati, and Vivek Y Reddy. Preclinical evaluation of pulsed field ablation: electrophysiological and histological assessment of thoracic vein isolation. *Circulation: Arrhythmia and Electrophysiology*, 12(12):e007781, 2019.
[36] Nicholas Labarbera and Corina Drapaca. Anistropically varying conductivity in irreversible electroporation simulations. *Theoretical Biology and Medical Modelling*, 14(1), nov 1 2017.

[37] Michael Leguebe. *Modélisation de l'électroperméabilisation à l'échelle cellulaire*. Theses, Université de Bordeaux, September 2014.

[38] M. Leguèbe, A. Silve, L.M. Mir, and C. Poignard. Conducting and permeable states of cell membrane submitted to high voltage pulses: Mathematical and numerical studies validated by the experiments. *Journal of Theoretical Biology*, 360:83-94, 2014.

[39] Robert Lipton. Heat conduction in fine scale mixtures with interfacial contact resistance. *SIAM Journal on Applied Mathematics*, 58(1):55-72, 1998.

[40] Daniel Meckes, Mehrdad Emami, Ian Fong, Dennis H. Lau, and Prashanthan Sanders. Pulsed-field ablation: Computational modeling of electric fields for lesion depth analysis. *Heart Rhythm O2*, 3(4):433-440, 2022.

[41] Damijan Miklavcic and Leila Towhidi. Numerical study of the electroporation pulse shape effect on molecular uptake of biological cells. *Radiology and Oncology*, 44(1):34-41, 2010.

[42] C. C. Mitchell and D. G. Schaeffer. A two-current model for the dynamics of cardiac membrane. *Bulletin of mathematical biology*, 65(5):767-793, 2003.

[43] Sunil Nath, John P DiMarco, and David E Haines. Basic aspects of radiofrequency catheter ablation. *Journal of Cardiovascular Electrophysiology*, 5(10):863-876, 1994.

[44] Simone Nati Poltri. *Mathematical modeling of cardiac tissue response after Pulsed Field Ablation*. PhD thesis, Université de Bordeaux, December 2024.

[45] Gabriel Nguetseng. A general convergence result for a functional related to the theory of homogenization. *SIAM Journal on Mathematical Analysis*, 20(3):608-623, 1989.

[46] Micol Pennacchio, Giuseppe Savaré, and Piero Colli Franzone. Multiscale modeling for the bioelectric activity of the heart. *SIAM Journal on Mathematical Analysis*, 37:1333-1370, 2005.

[47] Malte A. Peter and Michael Böhm. Different choices of scaling in homogenization of diffusion and interfacial exchange in a porous medium. *Mathematical Methods in the Applied Sciences*, 31(11):1257-1282, 2008.

[48] Argyrios Petras, Gerard Amoros Figueras, Zoraida Moreno Weidmann, Tomás García-Sánchez, David Viladés Medel, Antoni Ivorra, Jose M. Guerra, and Luca Gerardo-Giorda. Is a single lethal electric field threshold sufficient to characterize the lesion size in computational modeling of cardiac pulsed-field ablation? *Heart Rhythm O2*, 6(5):671-677, 2025.

[49] Argyrios Petras, Gerard Amoros Figueras, Zoraida Moreno Weidmann, Tomás García-Sánchez, David Viladés Medel, Antoni Ivorra, Jose M. Guerra, and Luca Gerardo-Giorda. Is a single lethal electric field threshold sufficient to characterize the lesion size in computational modeling of cardiac pulsed-field ablation? *Heart Rhythm O2*, 6(5):671-677, 2025.
[50] Vivek Y Reddy, Srinivas R Dukkipati, Petr Neuzil, Ante Anic, Jan Petru, Moritoshi Furasako, Hubert Cochet, Kentaro Minami, Toni Breskovic, Ivan Sikiric, et al. Pulsed field ablation of paroxysmal atrial fibrillation: 1-year outcomes of impulse, pefcat, and pefcat ii. *Clinical Electrophysiology*, 7(5):614-627, 2021.

[51] V. P. Smyshlyaev and K. D. Cherednichenko. On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 48(6):1325–1357, June 2000.

[52] Jeroen G. Stinstra, Bruce Hopenfeld, and Rob S. MacLeod. On the passive cardiac conductivity. *Annals of Biomedical Engineering*, 33:1743–1751, 2005.

[53] Chadi Tabaja, Arwa Younis, Ayman A. Hussein, Tyler L. Taigen, Hiroshi Nakagawa, Walid I. Saliba, Jakub Sroubek, Pasquale Santangeli, and Oussama M. Wazni. Catheter-based electroporation: A novel technique for catheter ablation of cardiac arrhythmias. *JACC: Clinical Electrophysiology*, 9(9):2008–2023, 2023.

[54] K. ten Tusscher and A. Panfilov. Alternans and spiral breakup in a human ventricular tissue model. *American Journal of Physiology-Heart and Circulatory Physiology*, 291(3), 2006.

[55] O. Tovar and L. Tung. Electroporation and recovery of cardiac cell membrane with rectangular voltage pulses. *American Journal of Physiology*, 263(4), 1992.

[56] L. Trabucho and J. M. Viaño. Existence and characterization of higher-order terms in an asymptotic expansion method for linearized elastic beams. *Asymptotic Analysis*, 2(3):223–255, January 1989.

[57] Damien Voyer, Aude Silve, Lluis M. Mir, Riccardo Scorretti, and Clair Poignard. Dynamical modeling of tissue electroporation. *Bioelectrochemistry*, 119:98–110, 2018.

[58] D. Šel, D. Cukjati, D. Batiuskaite, T. Slivnik, L.M. Mir, and D. Miklavčič. Sequential finite element model of tissue electropermeabilization. *IEEE Transactions on Biomedical Engineering*, 52(5):816–827, 2005.

[59] Fred HM Wittkampf and Hiroshi Nakagawa. RF catheter ablation: lessons on lesions. *Pacing and Clinical Electrophysiology*, 29(11):1285–1297, 2006.