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

HAL
open science

# Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu

► To cite this version:

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu. Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method. 2024. hal-04530917

HAL Id: hal-04530917

https://hal.science/hal-04530917

Preprint submitted on 3 Apr 2024

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.
# Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán,

Sijing Liu

► To cite this version:

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu. Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method. Numerical Algorithms, 2024, 10.1007/s11075-024-01902-z, http://dx.doi.org/10.1007/s11075-024-01902-z

HAL Id: hal-04530917

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

Submitted on 3 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
ESTIMATES OF DISCRETE TIME DERIVATIVES FOR THE
PARABOLIC-PARABOLIC ROBIN-ROBIN COUPLING METHOD

ERIK BURMAN, REBECCA DURST, MIGUEL A. FERNÁNDEZ, JOHNNY GUZMÁN, AND SIJING LIU

ABSTRACT. We consider a loosely coupled, non-iterative Robin-Robin coupling method proposed and analyzed in [J. Numer. Math., 31(1):59-77, 2023] for a parabolic-parabolic interface problem and prove estimates for the discrete time derivatives of the scalar field in different norms. When the interface is flat and perpendicular to two of the edges of the domain we prove error estimates in the $H^2$-norm. Such estimates are key ingredients to analyze a defect correction method for the parabolic-parabolic interface problem. Numerical results are shown to support our findings.

# 1. INTRODUCTION

Time splitting methods are popular for fluid-structure interaction (FSI) problems (see, e.g., [2, 16, 4, 17, 9, 10, 24, 25, 30, 28, 8]). One of the first stable splitting methods was proposed by Burman and Fernández [16]. Later the same authors [17] developed a related method which was coined the genuine Robin-Robin splitting method. Both these methods however suffered from a coupling between the space and time discretization parameters that reduced the accuracy. This constraint was lifted by eliminating the mesh dependence of the Robin parameter in [15]. The resulting method has been analyzed in [15, 13, 11]. A very similar method was developed and analyzed by Bukac and Seboldt [29]. In [13] it was proved for the FSI problems that the method converges as $O(\sqrt{\Delta t})$ where $\Delta t$ is the time step. Numerical evidence suggested that those estimates were not sharp, and nearly first-order accuracy was proved (mod possibly a logarithmic factor) in [11] for the analogue method applied to the parabolic-parabolic and wave-parabolic problems. The analysis was extended to the FSI problem in [22, 14]. For parabolic-parabolic couplings, there is a rich literature on splitting schemes motivated by models of ocean-atmosphere interaction. In these models, friction forces on the interface render the physical coupling dissipative through a Robin-type coupling condition, as discussed in [26]. This aspect has been successfully exploited in the design of splitting methods [20, 21, 19, 33, 32, 34, 31, 27]. In our case, the coupling conditions consist of continuity of both the primal variables and the fluxes across the interface. This coupling is conservative, and hence the approach suggested in the above references fails. Instead, the splitting method uses a Robin condition for the computational coupling, which turns out to lead to an unconditionally stable algorithm. An approach for conservative fluid-fluid coupling problems was proposed in [23], using Nitsche or Robin type couplings similar to those introduced in [17]. In a domain decomposition framework, a splitting method based on subcycling, i.e., iterative solution, was proposed and analyzed in [6, 5]. In a similar spirit, but focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection-diffusion was introduced in [18], with numerical investigation of the stability.

It is well known that discrete time differences for time stepping methods (e.g., backward Euler method) superconverge. In particular, when the backward Euler method is applied to a parabolic problem, the first-time difference of the errors converge with order $O((\Delta t)^2)$. In this paper, we prove similar results for the splitting method [11] applied to an interface problem. In particular, we will prove second-order convergence for the scalar fields living on the two sub-domains in the $L^2$-norm. Moreover, in special configurations (i.e., the interface is flat and perpendicular to two sides of the domain), we will be able to prove second-order convergence in the $H^2$-norm. Numerically, the $H^2$ second-order convergence rates seem to hold on more general configurations. This appears to be the
ESTIMATES OF DISCRETE TIME DERIVATIVES FOR THE
PARABOLIC-PARABOLIC ROBIN-ROBIN COUPLING METHOD

ERIK BURMAN, REBECCA DURST, MIGUEL A. FERNÁNDEZ, JOHNNY GUZMÁN, AND SIJING LIU

ABSTRACT. We consider a loosely coupled, non-iterative Robin-Robin coupling method proposed and analyzed in [J. Numer. Math., 31(1):59-77, 2023] for a parabolic-parabolic interface problem and prove estimates for the discrete time derivatives of the scalar field in different norms. When the interface is flat and perpendicular to two of the edges of the domain we prove error estimates in the $H^2$-norm. Such estimates are key ingredients to analyze a defect correction method for the parabolic-parabolic interface problem. Numerical results are shown to support our findings.

# 1. INTRODUCTION

Time splitting methods are popular for fluid-structure interaction (FSI) problems (see, e.g., [2, 16, 4, 17, 9, 10, 24, 25, 30, 28, 8]). One of the first stable splitting methods was proposed by Burman and Fernández [16]. Later the same authors [17] developed a related method which was coined the genuine Robin-Robin splitting method. Both these methods however suffered from a coupling between the space and time discretization parameters that reduced the accuracy. This constraint was lifted by eliminating the mesh dependence of the Robin parameter in [15]. The resulting method has been analyzed in [15, 13, 11]. A very similar method was developed and analyzed by Bukac and Seboldt [29]. In [13] it was proved for the FSI problems that the method converges as $O(\sqrt{\Delta t})$ where $\Delta t$ is the time step. Numerical evidence suggested that those estimates were not sharp, and nearly first-order accuracy was proved (mod possibly a logarithmic factor) in [11] for the analogue method applied to the parabolic-parabolic and wave-parabolic problems. The analysis was extended to the FSI problem in [22, 14]. For parabolic-parabolic couplings, there is a rich literature on splitting schemes motivated by models of ocean-atmosphere interaction. In these models, friction forces on the interface render the physical coupling dissipative through a Robin-type coupling condition, as discussed in [26]. This aspect has been successfully exploited in the design of splitting methods [20, 21, 19, 33, 32, 34, 31, 27]. In our case, the coupling conditions consist of continuity of both the primal variables and the fluxes across the interface. This coupling is conservative, and hence the approach suggested in the above references fails. Instead, the splitting method uses a Robin condition for the computational coupling, which turns out to lead to an unconditionally stable algorithm. An approach for conservative fluid-fluid coupling problems was proposed in [23], using Nitsche or Robin type couplings similar to those introduced in [17]. In a domain decomposition framework, a splitting method based on subcycling, i.e., iterative solution, was proposed and analyzed in [6, 5]. In a similar spirit, but focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection-diffusion was introduced in [18], with numerical investigation of the stability.

It is well known that discrete time differences for time stepping methods (e.g., backward Euler method) superconverge. In particular, when the backward Euler method is applied to a parabolic problem, the first-time difference of the errors converge with order $O((\Delta t)^2)$. In this paper, we prove similar results for the splitting method [11] applied to an interface problem. In particular, we will prove second-order convergence for the scalar fields living on the two sub-domains in the $L^2$-norm. Moreover, in special configurations (i.e., the interface is flat and perpendicular to two sides of the domain), we will be able to prove second-order convergence in the $H^2$-norm. Numerically, the $H^2$ second-order convergence rates seem to hold on more general configurations. This appears to be the
first work where error estimates in stronger norms for splitting methods applied to interface problems have been considered.

Error estimates of derivatives of the solution are, of course, of interest in their own right in many applications. However, our main motivation for this work is the application of these estimates in the analysis of a prediction-correction method [12] that we are concurrently developing. The aim is to improve the first order convergence of [11] to second order convergence in time through a defect-correction procedure [7]. The method in [12] uses a prediction step, which is exactly the splitting method we analyze here. The second-order accuracy of the correction step depends on the second-order accuracy of time differences of the prediction step, which is precisely the subject of this paper.

As for a single parabolic problem, the idea to prove higher convergence for time differences is to use that the time differences satisfy a similar discrete equation with new right-hand sides that have time differences themselves. Then, one uses the error analysis for the original method to proceed. In our case, we will use the analysis provided in [11] to do this. The main difference here is that we consider a problem with Neumann boundary conditions on two of the sides of a square instead of pure Dirichlet boundary conditions, which were considered in [11]. This, in fact, simplifies the analysis slightly, and additionally, one can remove the logarithmic factor that appears in [11]. It should be mentioned that we only consider the time discrete case. The fully discrete case is more involved.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding Robin-Robin coupling method. In Section 3, we present stability results for a Robin-Robin method. Section 4 is devoted to the error estimates. Finally, we provide some numerical results in Section 5 and end with some concluding remarks in Section 6.

## 2. THE PARABOLIC-PARABOLIC INTERFACE PROBLEM AND ROBIN-ROBIN COUPLING METHOD

Let $\Omega = (0, 1)^2$ and suppose that $\Omega = \Omega_f \cup \Omega_s \cup \Sigma$. The interface $\Sigma$ is assumed to be a line segment that intersects $\Omega$ on the two side edges; see Figure 1. We let $\Gamma_N$ denote the two side edges of $\Omega$ and we let $\Gamma_D$ be the bottom and top edges of $\Omega$. We let $\Gamma_N^i = \Gamma_N^i \cap \partial\Omega_i$ for $i = s, f$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The domains $\Omega_f$ and $\Omega_s$ with interface $\Sigma$ and Neumann boundaries.</figcaption></figure>

### 2.1. The Parabolic-parabolic problem. We consider the interface problem

$$
\begin{align}
(2.1a) \qquad \partial_t u - \nu_f \Delta u &= 0, && \text{in } [0, T] \times \Omega_f, \nonumber \\
&u(0, x) = u_0(x), && \text{on } \Omega_f, \nonumber \\
&u = 0, && \text{on } [0, T] \times \Gamma_D^f, \nonumber \\
&\partial_n u = 0, && \text{on } [0, T] \times \Gamma_N^f, \nonumber
\end{align}
$$
first work where error estimates in stronger norms for splitting methods applied to interface problems have been considered.

Error estimates of derivatives of the solution are, of course, of interest in their own right in many applications. However, our main motivation for this work is the application of these estimates in the analysis of a prediction-correction method [12] that we are concurrently developing. The aim is to improve the first order convergence of [11] to second order convergence in time through a defect-correction procedure [7]. The method in [12] uses a prediction step, which is exactly the splitting method we analyze here. The second-order accuracy of the correction step depends on the second-order accuracy of time differences of the prediction step, which is precisely the subject of this paper.

As for a single parabolic problem, the idea to prove higher convergence for time differences is to use that the time differences satisfy a similar discrete equation with new right-hand sides that have time differences themselves. Then, one uses the error analysis for the original method to proceed. In our case, we will use the analysis provided in [11] to do this. The main difference here is that we consider a problem with Neumann boundary conditions on two of the sides of a square instead of pure Dirichlet boundary conditions, which were considered in [11]. This, in fact, simplifies the analysis slightly, and additionally, one can remove the logarithmic factor that appears in [11]. It should be mentioned that we only consider the time discrete case. The fully discrete case is more involved.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding Robin-Robin coupling method. In Section 3, we present stability results for a Robin-Robin method. Section 4 is devoted to the error estimates. Finally, we provide some numerical results in Section 5 and end with some concluding remarks in Section 6.

## 2. THE PARABOLIC-PARABOLIC INTERFACE PROBLEM AND ROBIN-ROBIN COUPLING METHOD

Let $\Omega = (0, 1)^2$ and suppose that $\Omega = \Omega_f \cup \Omega_s \cup \Sigma$. The interface $\Sigma$ is assumed to be a line segment that intersects $\Omega$ on the two side edges; see Figure 1. We let $\Gamma_N$ denote the two side edges of $\Omega$ and we let $\Gamma_D$ be the bottom and top edges of $\Omega$. We let $\Gamma_N^i = \Gamma_N^i \cap \partial\Omega_i$ for $i = s, f$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The domains $\Omega_f$ and $\Omega_s$ with interface $\Sigma$ and Neumann boundaries.</figcaption></figure>

### 2.1. The Parabolic-parabolic problem. We consider the interface problem

$$
\begin{align*}
(2.1a) \qquad \partial_t u - \nu_f \Delta u &= 0, && \text{in } [0, T] \times \Omega_f, \\
&& u(0, x) &= u_0(x), && \text{on } \Omega_f, \\
&& u &= 0, && \text{on } [0, T] \times \Gamma_D^f, \\
&& \partial_n u &= 0, && \text{on } [0, T] \times \Gamma_N^f,
\end{align*}
$$
$$
\begin{align}
& \partial_t \mathbf{w} - \nu_s \Delta \mathbf{w} = 0, && \text{in } [0, T] \times \Omega_s, \nonumber \\
(2.1\text{b}) & \mathbf{w}(0, x) = \mathbf{w}_0(x), && \text{on } \Omega_s, \nonumber \\
& \mathbf{w} = 0, && \text{on } [0, T] \times \Gamma_D^s, \nonumber \\
& \partial_n \mathbf{w} = 0, && \text{on } [0, T] \times \Gamma_N^s, \nonumber
\end{align}
$$

$$
(2.1c) \qquad \mathbf{w} - \mathbf{u} = 0, \quad \text{in } [0, T] \times \Sigma,
$$

$$
(2.1\text{d}) \qquad \nu_s \nabla \mathbf{w} \cdot \mathbf{n}_s + \nu_f \nabla \mathbf{u} \cdot \mathbf{n}_f = 0, \quad \text{in } [0, T] \times \Sigma,
$$

where $\mathbf{n}_f$ and $\mathbf{n}_s$ are the outward facing normal vectors for $\Omega_f$ and $\Omega_s$, respectively. We assume that the initial data is smooth and that $\mathbf{u}$ and $\mathbf{w}$ are smooth on $\Omega_f$ and $\Omega_s$, respectively.

2.2. **Variational form.** Let $(\cdot, \cdot)_i$ be the $L^2$-inner product on $\Omega_i$ for $i = f, s$. Moreover, let $\langle \cdot, \cdot \rangle$ be the $L^2$-inner product on $\Sigma$. Let $N > 0$ be an integer, and define $\Delta t := \frac{T}{N}$, and let $\mathbf{u}^n := \mathbf{u}(t_n, \cdot)$, where $t_n := n\Delta t$ for $n \in \{0, 1, 2, \dots, N\}$. We consider the spaces

$$
(2.2a) \qquad V_f := \{v \in H^1(\Omega_f) : v = 0 \text{ on } \Gamma_D^f\},
$$

$$
(2.2b) \qquad V_s := \{v \in H^1(\Omega_s) : v = 0 \text{ on } \Gamma_D^s\},
$$

$$
(2.2c) \qquad V_g := L^2(\Sigma).
$$

By setting $|n+1 := \nu_f \partial_{n_f} \mathbf{u}^{n+1}$ and assuming that $|n+1 \in L^2(\Sigma)$ for all $n$, the solution to (2.1) also
satisfies the following variational formulation, for $n = 0, \dots, N-1$:

$$
(2.3a) \qquad (\partial_t \mathbf{w}^{n+1}, z)_s + \nu_s (\nabla \mathbf{w}^{n+1}, \nabla z)_s + \langle |^{n+1}, z \rangle = 0, \quad z \in V_s,
$$

$$
(2.3b) \qquad (\partial_t \mathbf{u}^{n+1}, v)_f + \nu_f (\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle = 0, \quad v \in V_f,
$$

$$
(2.3c) \qquad \langle \mathbf{w}^{n+1} - \mathbf{u}^{n+1}, \mu \rangle = 0, \quad \mu \in V_g.
$$

2.3. Robin-Robin coupling: time discrete method. We define the discrete time derivatives:

$$
\partial_{\Delta t} v^{n+1} = \frac{v^{n+1} - v^n}{\Delta t},
$$

and

$$
\partial_{\Delta t}^2 v^{n+1} = \frac{v^{n+1} - 2v^n + v^{n-1}}{(\Delta t)^2}.
$$

The Robin-Robin method solves sequentially:

$$
(2.4a) \qquad \partial_{\Delta t} w^{n+1} - \nu_s \Delta w^{n+1} = 0, \qquad \text{on } \Omega_s,
$$

$$
(2.4b) \qquad w^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(2.4c) \qquad \partial_n w^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(2.4d) \qquad \alpha w^{n+1} + \nu_s \partial_{n_s} w^{n+1} = \alpha u^n - \nu_f \partial_{n_f} u^n \quad \text{on } \Sigma.
$$

$$
(2.5a) \qquad \partial_{\Delta t} u^{n+1} - \nu_f \Delta u^{n+1} = 0, \qquad \text{in } \Omega_f,
$$

$$
(2.5b) \qquad u^{n+1} = 0, \qquad \text{on } \Gamma_D^f,
$$

$$
(2.5c) \qquad \partial_n u^{n+1} = 0, \qquad \text{on } \Gamma_N^f,
$$

$$
(2.5d) \qquad \alpha u^{n+1} + \nu_f \partial_{n_f} u^{n+1} = \alpha w^{n+1} + \nu_f \partial_{n_f} u^n \quad \text{on } \Sigma.
$$
$$
\begin{align}
& \partial_t w - \nu_s \Delta w = 0, && \text{in } [0, T] \times \Omega_s, \nonumber \\
(2.1\text{b}) & w(0, x) = w_0(x), && \text{on } \Omega_s, \nonumber \\
& w = 0, && \text{on } [0, T] \times \Gamma_D^s, \nonumber \\
& \partial_n w = 0, && \text{on } [0, T] \times \Gamma_N^s, \nonumber
\end{align}
$$

$$
(2.1c) \qquad w - u = 0, \quad \text{in } [0, T] \times \Sigma,
$$

$$
(2.1\text{d}) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0, \quad \text{in } [0, T] \times \Sigma,
$$

where $\mathbf{n}_f$ and $\mathbf{n}_s$ are the outward facing normal vectors for $\Omega_f$ and $\Omega_s$, respectively. We assume that the initial data is smooth and that $\mathbf{u}$ and $\mathbf{w}$ are smooth on $\Omega_f$ and $\Omega_s$, respectively.

2.2. **Variational form.** Let $(\cdot, \cdot)_i$ be the $L^2$-inner product on $\Omega_i$ for $i = f, s$. Moreover, let $\langle \cdot, \cdot \rangle$ be the $L^2$-inner product on $\Sigma$. Let $N > 0$ be an integer, and define $\Delta t := \frac{T}{N}$, and let $\mathbf{u}^n := \mathbf{u}(t_n, \cdot)$, where $t_n := n\Delta t$ for $n \in \{0, 1, 2, \dots, N\}$. We consider the spaces

$$
(2.2a) \qquad V_f := \{v \in H^1(\Omega_f) : v = 0 \text{ on } \Gamma_D^f\},
$$

$$
(2.2b) \qquad V_s := \{v \in H^1(\Omega_s) : v = 0 \text{ on } \Gamma_D^s\},
$$

$$
(2.2c) \qquad V_g := L^2(\Sigma).
$$

By setting $|n^{n+1}| := \nu_f \partial_{n_f} \mathbf{u}^{n+1}$ and assuming that $|n^{n+1} \in L^2(\Sigma)$ for all $n$, the solution to (2.1) also satisfies the following variational formulation, for $n = 0, \dots, N-1$:

$$
(2.3a) \qquad (\partial_t \mathbf{w}^{n+1}, z)_s + \nu_s (\nabla \mathbf{w}^{n+1}, \nabla z)_s + \langle |^{n+1}, z \rangle = 0, \quad z \in V_s,
$$

$$
(2.3b) \qquad (\partial_t \mathbf{u}^{n+1}, v)_f + \nu_f (\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle = 0, \quad v \in V_f,
$$

$$
(2.3c) \qquad \langle \mathbf{w}^{n+1} - \mathbf{u}^{n+1}, \mu \rangle = 0, \quad \mu \in V_g.
$$

2.3. Robin-Robin coupling: time discrete method. We define the discrete time derivatives:

$$
\partial_{\Delta t} v^{n+1} = \frac{v^{n+1} - v^n}{\Delta t},
$$

and

$$
\partial_{\Delta t}^2 v^{n+1} = \frac{v^{n+1} - 2v^n + v^{n-1}}{(\Delta t)^2}.
$$

The Robin-Robin method solves sequentially:

$$
(2.4a) \qquad \partial_{\Delta t} w^{n+1} - \nu_s \Delta w^{n+1} = 0, \qquad \text{on } \Omega_s,
$$

$$
(2.4b) \qquad w^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(2.4c) \qquad \partial_n w^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(2.4d) \qquad \alpha w^{n+1} + \nu_s \partial_{n_s} w^{n+1} = \alpha u^n - \nu_f \partial_{n_f} u^n \quad \text{on } \Sigma.
$$

$$
(2.5a) \qquad \partial_{\Delta t} u^{n+1} - \nu_f \Delta u^{n+1} = 0, \qquad \text{in } \Omega_f,
$$

$$
(2.5b) \qquad u^{n+1} = 0, \qquad \text{on } \Gamma_D^f,
$$

$$
(2.5c) \qquad \partial_n u^{n+1} = 0, \qquad \text{on } \Gamma_N^f,
$$

$$
(2.5d) \qquad \alpha u^{n+1} + \nu_f \partial_{n_f} u^{n+1} = \alpha w^{n+1} + \nu_f \partial_{n_f} u^n \quad \text{on } \Sigma.
$$
We let $\lambda^{n+1} = \nu_f \partial_{n_f} u^{n+1}$. Then, the time semi-discrete solution solves the following: Find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, and $\lambda^{n+1} \in V_g$ such that, for $n \ge 0$,

$$ (2.6a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = 0, \quad z \in V_s, $$

$$ (2.6b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = 0, \quad v \in V_f, $$

$$ (2.6c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = 0, \quad \mu \in V_g, $$

with $u^0 = u_0(x)$ and $w^0 = w_0(x)$.

One can re-write (2.6c) as

$$ (2.7) \qquad \alpha(u^{n+1} - w^{n+1}) = \lambda^n - \lambda^{n+1}, \quad \text{on } \Sigma. $$

### 3. STABILITY RESULT

In this section we will prove stability results for the Robin-Robin method with a more general right-hand side: Find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, such that, for $n \ge 0$,

$$ (3.1a) \qquad \partial_{\Delta t} w^{n+1} - \nu_s \Delta w^{n+1} = b_1^{n+1}, \qquad \text{on } \Omega_s, $$

$$ (3.1b) \qquad w^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.1c) \qquad \partial_n w^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.1d) \qquad \alpha w^{n+1} + \nu_s \partial_{n_s} w^{n+1} = \alpha u^n - \nu_f \partial_{n_f} u^n + \varepsilon_1^{n+1} \qquad \text{on } \Sigma. $$

$$ (3.2a) \qquad \partial_{\Delta t} u^{n+1} - \nu_f \Delta u^{n+1} = b_2^{n+1}, \qquad \text{in } \Omega_f, $$

$$ (3.2b) \qquad u^{n+1} = 0, \qquad \text{on } \Gamma_D^f, $$

$$ (3.2c) \qquad \partial_n u^{n+1} = 0, \qquad \text{on } \Gamma_N^f, $$

$$ (3.2d) \qquad \alpha u^{n+1} + \nu_f \partial_{n_f} u^{n+1} = \alpha w^{n+1} + \nu_f \partial_{n_f} u^n + \varepsilon_2^{n+1} \qquad \text{on } \Sigma, $$

with

$$ (3.3) \qquad u^0 = w^0 = b_i^0 = \varepsilon_i^0 = 0 \quad i = 1, 2. $$

We assumed zero initial conditions for simplicity and this will be enough for the error analysis below. When analyzing the error of the Robin-Robin method the terms $\{b_i^n\}$, $\{\varepsilon_i^n\}$ will be the residual terms.

**Remark 3.1 (Regularity).** We assume that $u^{n+1} \in H^2(\Omega_f)$ and $w^{n+1} \in H^2(\Omega_s)$ for the remainder of this paper.

**Remark 3.2 (Notations).** For any negative superscripts, we set the values of the terms on the right-hand side to be zero, i.e.,

$$ (3.4) \qquad u^j = w^j = b_i^j = \varepsilon_i^j = 0 \quad j < 0, \quad i = 1, 2. $$

We let $\lambda^{n+1} = \nu_f \partial_{n_f} u^{n+1}$ and we see that the above solution satisfies, for $n \ge 0$,

$$ (3.5a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = L_1(z), \quad z \in V_s, $$

$$ (3.5b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = L_2(v), \quad v \in V_f, $$

$$ (3.5c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = L_3(\mu), \quad \mu \in V_g, $$
We let $\lambda^{n+1} = \nu_f \partial_{n_f} u^{n+1}$. Then, the time semi-discrete solution solves the following: Find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, and $\lambda^{n+1} \in V_g$ such that, for $n \ge 0$,

$$ (2.6a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = 0, \quad z \in V_s, $$

$$ (2.6b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = 0, \quad v \in V_f, $$

$$ (2.6c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = 0, \quad \mu \in V_g, $$

with $u^0 = u_0(x)$ and $w^0 = w_0(x)$.

One can re-write (2.6c) as

$$ (2.7) \qquad \alpha(u^{n+1} - w^{n+1}) = \lambda^n - \lambda^{n+1}, \quad \text{on } \Sigma. $$

### 3. STABILITY RESULT

In this section we will prove stability results for the Robin-Robin method with a more general right-hand side: Find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, such that, for $n \ge 0$,

$$ (3.1a) \qquad \partial_{\Delta t} w^{n+1} - \nu_s \Delta w^{n+1} = b_1^{n+1}, \qquad \text{on } \Omega_s, $$

$$ (3.1b) \qquad w^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.1c) \qquad \partial_n w^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.1d) \qquad \alpha w^{n+1} + \nu_s \partial_{n_s} w^{n+1} = \alpha u^n - \nu_f \partial_{n_f} u^n + \epsilon_1^{n+1} \qquad \text{on } \Sigma. $$

$$ (3.2a) \qquad \partial_{\Delta t} u^{n+1} - \nu_f \Delta u^{n+1} = b_2^{n+1}, \qquad \text{in } \Omega_f, $$

$$ (3.2b) \qquad u^{n+1} = 0, \qquad \text{on } \Gamma_D^f, $$

$$ (3.2c) \qquad \partial_n u^{n+1} = 0, \qquad \text{on } \Gamma_N^f, $$

$$ (3.2d) \qquad \alpha u^{n+1} + \nu_f \partial_{n_f} u^{n+1} = \alpha w^{n+1} + \nu_f \partial_{n_f} u^n + \epsilon_2^{n+1} \qquad \text{on } \Sigma, $$

with

$$ (3.3) \qquad u^0 = w^0 = b_i^0 = \epsilon_i^0 = 0 \quad i = 1, 2. $$

We assumed zero initial conditions for simplicity and this will be enough for the error analysis below. When analyzing the error of the Robin-Robin method the terms $\{b_i^n\}$, $\{\epsilon_i^n\}$ will be the residual terms.

**Remark 3.1 (Regularity).** We assume that $u^{n+1} \in H^2(\Omega_f)$ and $w^{n+1} \in H^2(\Omega_s)$ for the remainder of this paper.

**Remark 3.2 (Notations).** For any negative superscripts, we set the values of the terms on the right-hand side to be zero, i.e.,

$$ (3.4) \qquad u^j = w^j = b_i^j = \epsilon_i^j = 0 \quad j < 0, \quad i = 1, 2. $$

We let $\lambda^{n+1} = \nu_f \partial_{n_f} u^{n+1}$ and we see that the above solution satisfies, for $n \ge 0$,

$$ (3.5a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = L_1(z), \quad z \in V_s, $$

$$ (3.5b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = L_2(v), \quad v \in V_f, $$

$$ (3.5c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = L_3(\mu), \quad \mu \in V_g, $$
where

$$L_1(z) := (b_1^{n+1}, z)_s + \langle \epsilon_1^{n+1}, z \rangle,$$

$$L_2(v) := (b_2^{n+1}, v)_f,$$

$$L_3(\mu) := \langle \epsilon_2^{n+1}, \mu \rangle.$$

For convenience we rewrite (3.5c) as:

$$(3.6) \qquad \alpha(u^{n+1} - w^{n+1}) = \lambda^n - \lambda^{n+1} + \epsilon_2^{n+1}, \quad \text{on } \Sigma.$$

We will need to define the following quantities for the stability estimates.

$$Z^{n+1}(\psi, \phi, \theta) := \frac{1}{2} \|\phi^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|\psi^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{\Delta t \alpha}{2} \|\phi^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\theta^{n+1}\|_{L^2(\Sigma)}^2,$$

$$S^{n+1}(\psi, \phi, \theta) := \Delta t (\nu_f \|\nabla \phi^{n+1}\|_{L^2(\Omega_f)}^2 + \nu_s \|\nabla \psi^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|\psi^{n+1} - \psi^n\|_{L^2(\Omega_s)}^2 + \|\phi^{n+1} - \phi^n\|_{L^2(\Omega_f)}^2) \\ + \frac{\alpha \Delta t}{2} \|\phi^{n+1} - \phi^n + \frac{1}{\alpha}(\theta^{n+1} - \theta^n)\|_{L^2(\Sigma)}^2.$$

We first state a preliminary result.

**Lemma 3.3.** Let $w, u$ solve (3.1) and (3.2) then the following identity holds.

$$(3.7) \qquad Z^{n+1}(w, u, \lambda) + S^{n+1}(w, u, \lambda) = Z^n(w, u, \lambda) + \Delta t F^{n+1}(w, u) + \frac{\Delta t}{\alpha} \langle \epsilon_2^{n+1}, \lambda^{n+1} \rangle,$$

where

$$F^{n+1}(w, u) := (b_1^{n+1}, w^{n+1})_s + (b_2^{n+1}, u^{n+1})_f + \langle \epsilon_1^{n+1} + \epsilon_2^{n+1}, w^{n+1} \rangle + \langle u^{n+1} - u^n, \epsilon_2^{n+1} \rangle.$$

*Proof.* To begin, we set $z = \Delta t w^{n+1}$ in (3.5a) and $v = \Delta t u^{n+1}$ in (3.5b) to get

$$(3.8) \qquad \frac{1}{2} \|w^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|w^{n+1} - w^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^{n+1} - u^n\|_{L^2(\Omega_f)}^2 \\ + \nu_s \Delta t \|\nabla w^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \Delta t \|\nabla u^{n+1}\|_{L^2(\Omega_f)}^2 \\ = \frac{1}{2} \|w^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^n\|_{L^2(\Omega_f)}^2 + \Delta t J^{n+1},$$

where

$$(3.9) \qquad J^{n+1} := -\alpha \langle w^{n+1} - u^n, w^{n+1} \rangle - \langle \lambda^n, w^{n+1} \rangle + \langle \lambda^{n+1}, u^{n+1} \rangle + L_1(w^{n+1}) + L_2(u^{n+1}).$$

Manipulating the first three terms in (3.9) and using (3.6), we obtain

$$(3.10) \qquad -\alpha \langle w^{n+1} - u^n, w^{n+1} \rangle - \langle \lambda^n, w^{n+1} \rangle + \langle \lambda^{n+1}, u^{n+1} \rangle \\ = J^{n+1} + \frac{1}{\alpha} \langle \epsilon_2^{n+1}, \lambda^{n+1} \rangle + \langle \epsilon_2^{n+1}, w^{n+1} \rangle - \langle u^n - u^{n+1}, \epsilon_2^{n+1} \rangle,$$

with

$$J^{n+1} := \alpha \langle u^n - u^{n+1}, u^{n+1} \rangle + \frac{1}{\alpha} \langle \lambda^n - \lambda^{n+1}, \lambda^{n+1} \rangle - \langle u^n - u^{n+1}, \lambda^n - \lambda^{n+1} \rangle.$$

One can easily show that

$$J^{n+1} = \frac{\alpha}{2} (\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\ - \frac{\alpha}{2} \|(u^n - u^{n+1}) + \frac{1}{\alpha}(\lambda^n - \lambda^{n+1})\|_{L^2(\Sigma)}^2.$$
where

$$L_1(z) := (b_1^{n+1}, z)_s + \langle \epsilon_1^{n+1}, z \rangle,$$

$$L_2(v) := (b_2^{n+1}, v)_f,$$

$$L_3(\mu) := \langle \epsilon_2^{n+1}, \mu \rangle.$$

For convenience we rewrite (3.5c) as:

$$(3.6) \qquad \alpha(u^{n+1} - w^{n+1}) = \lambda^n - \lambda^{n+1} + \epsilon_2^{n+1}, \quad \text{on } \Sigma.$$

We will need to define the following quantities for the stability estimates.

$$Z^{n+1}(\psi, \phi, \theta) := \frac{1}{2} \|\phi^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|\psi^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{\Delta t \alpha}{2} \|\phi^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\theta^{n+1}\|_{L^2(\Sigma)}^2,$$

$$S^{n+1}(\psi, \phi, \theta) := \Delta t (\nu_f \|\nabla \phi^{n+1}\|_{L^2(\Omega_f)}^2 + \nu_s \|\nabla \psi^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|\psi^{n+1} - \psi^n\|_{L^2(\Omega_s)}^2 + \|\phi^{n+1} - \phi^n\|_{L^2(\Omega_f)}^2) \\ + \frac{\alpha \Delta t}{2} \|\phi^{n+1} - \phi^n + \frac{1}{\alpha}(\theta^{n+1} - \theta^n)\|_{L^2(\Sigma)}^2.$$

We first state a preliminary result.

**Lemma 3.3.** Let $w, u$ solve (3.1) and (3.2) then the following identity holds.

$$(3.7) \qquad Z^{n+1}(w, u, \lambda) + S^{n+1}(w, u, \lambda) = Z^n(w, u, \lambda) + \Delta t F^{n+1}(w, u) + \frac{\Delta t}{\alpha} \langle \epsilon_2^{n+1}, \lambda^{n+1} \rangle,$$

where

$$F^{n+1}(w, u) := (b_1^{n+1}, w^{n+1})_s + (b_2^{n+1}, u^{n+1})_f + \langle \epsilon_1^{n+1} + \epsilon_2^{n+1}, w^{n+1} \rangle + \langle u^{n+1} - u^n, \epsilon_2^{n+1} \rangle.$$

*Proof.* To begin, we set $z = \Delta t w^{n+1}$ in (3.5a) and $v = \Delta t u^{n+1}$ in (3.5b) to get

$$(3.8) \qquad \frac{1}{2} \|w^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|w^{n+1} - w^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^{n+1} - u^n\|_{L^2(\Omega_f)}^2 \\ + \nu_s \Delta t \|\nabla w^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \Delta t \|\nabla u^{n+1}\|_{L^2(\Omega_f)}^2 \\ = \frac{1}{2} \|w^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|u^n\|_{L^2(\Omega_f)}^2 + \Delta t J^{n+1},$$

where

$$(3.9) \qquad J^{n+1} := -\alpha \langle w^{n+1} - u^n, w^{n+1} \rangle - \langle \lambda^n, w^{n+1} \rangle + \langle \lambda^{n+1}, u^{n+1} \rangle + L_1(w^{n+1}) + L_2(u^{n+1}).$$

Manipulating the first three terms in (3.9) and using (3.6), we obtain

$$(3.10) \qquad -\alpha \langle w^{n+1} - u^n, w^{n+1} \rangle - \langle \lambda^n, w^{n+1} \rangle + \langle \lambda^{n+1}, u^{n+1} \rangle \\ = J^{n+1} + \frac{1}{\alpha} \langle \epsilon_2^{n+1}, \lambda^{n+1} \rangle + \langle \epsilon_2^{n+1}, w^{n+1} \rangle - \langle u^n - u^{n+1}, \epsilon_2^{n+1} \rangle,$$

with

$$J^{n+1} := \alpha \langle u^n - u^{n+1}, u^{n+1} \rangle + \frac{1}{\alpha} \langle \lambda^n - \lambda^{n+1}, \lambda^{n+1} \rangle - \langle u^n - u^{n+1}, \lambda^n - \lambda^{n+1} \rangle.$$

One can easily show that

$$J^{n+1} = \frac{\alpha}{2} (\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\ - \frac{\alpha}{2} \|(u^n - u^{n+1}) + \frac{1}{\alpha}(\lambda^n - \lambda^{n+1})\|_{L^2(\Sigma)}^2.$$
By combining this identity with (3.9) and (3.10), we arrive at

$$
\begin{align*}
J^{n+1} &= \frac{\alpha}{2} (\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\
&\quad - \frac{\alpha}{2} \|(u^n - u^{n+1}) + \frac{1}{\alpha}(\lambda^n - \lambda^{n+1})\|_{L^2(\Sigma)}^2 \\
&\quad + L_1(w^{n+1}) + L_2(u^{n+1}) + \frac{1}{\alpha} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle + \langle \varepsilon_2^{n+1}, w^{n+1} \rangle + \langle \varepsilon_2^{n+1}, u^{n+1} - u^n \rangle.
\end{align*}
$$

If we plug in these results to (3.8) we arrive at the identity.

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

We can now state an identity for the last term in (3.7).

**Lemma 3.4.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\varepsilon_2^{n+1} \in V_f$ then the following identity holds

$$
\frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = - \frac{\Delta t}{\alpha} \sum_{n=1}^{N-1} (u^n, \partial_{\Delta t} \varepsilon_2^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \left( \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - (b_2^{n+1}, \varepsilon_2^{n+1})_f \right) \\
+ \frac{1}{\alpha} (u^N, \varepsilon_2^N)_f.
$$

*Proof.* We first take $v = \Delta t \varepsilon_2^{n+1}$ in (3.5b) to obtain

$$
\Delta t \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = (u^{n+1} - u^n, \varepsilon_2^{n+1})_f + \Delta t \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - \Delta t (b_2^{n+1}, \varepsilon_2^{n+1})_f.
$$

If we take the sum over $n = 0, \dots, N-1$ and use summation by parts, we get

$$
\Delta t \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = - \Delta t \sum_{n=1}^{N-1} (u^n, \partial_{\Delta t} \varepsilon_2^{n+1})_f + \Delta t \sum_{n=0}^{N-1} \left( \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - (b_2^{n+1}, \varepsilon_2^{n+1})_f \right) \\
+ (u^N, \varepsilon_2^N)_f - (u^0, \varepsilon_2^1).
$$

We conclude the proof by using (3.3).

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

To state the stability estimate we need the next definition.

$$
\begin{align*}
\Xi^N(m_1, m_2, s_1, s_2) :={}& \Delta t \sum_{n=0}^{N-1} \left[ \frac{1}{\nu_s} \|m_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|m_2^{n+1}\|_{L^2(\Omega_f)}^2 \right] \\
& + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|s_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
& + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|s_1^{n+1} + s_2^{n+1}\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|s_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|s_2^N\|_{L^2(\Omega_f)}^2.
\end{align*}
$$

**Theorem 3.5.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\varepsilon_2^{n+1} \in V_f$ then the following estimate
holds

$$
Z^N(w, u, \lambda) + \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \leq C \Xi^N(b_1, b_2, \varepsilon_1, \varepsilon_2).
$$

*Proof.* Using Lemma 3.3 and taking the sum we get

$$
Z^N(w, u, \lambda) + \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) = Z^0(w, u, \lambda) + \Delta t \sum_{n=0}^{N-1} F^{n+1}(w, u) + \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle.
$$
By combining this identity with (3.9) and (3.10), we arrive at

$$
\begin{aligned}
J^{n+1} = & \frac{\alpha}{2} (\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\
& - \frac{\alpha}{2} \|(u^n - u^{n+1}) + \frac{1}{\alpha}(\lambda^n - \lambda^{n+1})\|_{L^2(\Sigma)}^2 \\
& + L_1(w^{n+1}) + L_2(u^{n+1}) + \frac{1}{\alpha} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle + \langle \varepsilon_2^{n+1}, w^{n+1} \rangle + \langle \varepsilon_2^{n+1}, u^{n+1} - u^n \rangle.
\end{aligned}
$$

If we plug in these results to (3.8) we arrive at the identity.

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

We can now state an identity for the last term in (3.7).

**Lemma 3.4.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\varepsilon_2^{n+1} \in V_f$ then the following identity holds

$$
\frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = - \frac{\Delta t}{\alpha} \sum_{n=1}^{N-1} (u^n, \partial_{\Delta t} \varepsilon_2^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \left( \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - (b_2^{n+1}, \varepsilon_2^{n+1})_f \right) \\
+ \frac{1}{\alpha} (u^N, \varepsilon_2^N)_f.
$$

*Proof.* We first take $v = \Delta t \varepsilon_2^{n+1}$ in (3.5b) to obtain

$$
\Delta t \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = (u^{n+1} - u^n, \varepsilon_2^{n+1})_f + \Delta t \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - \Delta t (b_2^{n+1}, \varepsilon_2^{n+1})_f.
$$

If we take the sum over $n = 0, \dots, N-1$ and use summation by parts, we get

$$
\Delta t \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle = - \Delta t \sum_{n=1}^{N-1} (u^n, \partial_{\Delta t} \varepsilon_2^{n+1})_f + \Delta t \sum_{n=0}^{N-1} \left( \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f - (b_2^{n+1}, \varepsilon_2^{n+1})_f \right) \\
+ (u^N, \varepsilon_2^N)_f - (u^0, \varepsilon_2^1).
$$

We conclude the proof by using (3.3).

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

To state the stability estimate we need the next definition.

$$
\begin{aligned}
\Xi^N(m_1, m_2, s_1, s_2) := & \Delta t \sum_{n=0}^{N-1} \left[ \frac{1}{\nu_s} \|m_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|m_2^{n+1}\|_{L^2(\Omega_f)}^2 \right] \\
& + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|s_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
& + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|s_1^{n+1} + s_2^{n+1}\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|s_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|s_2^N\|_{L^2(\Omega_f)}^2.
\end{aligned}
$$

**Theorem 3.5.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\varepsilon_2^{n+1} \in V_f$ then the following estimate holds

$$
Z^N(w, u, \lambda) + \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \leq C \Xi^N(b_1, b_2, \varepsilon_1, \varepsilon_2).
$$

*Proof.* Using Lemma 3.3 and taking the sum we get

$$
Z^N(w, u, \lambda) + \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) = Z^0(w, u, \lambda) + \Delta t \sum_{n=0}^{N-1} F^{n+1}(w, u) + \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle.
$$
Using the Poincaré and trace inequalities, we easily obtain the following bound

$$ \Delta t \sum_{n=0}^{N-1} F^{n+1}(w, u) \le \frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|b_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{\nu_f} \|b_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\varepsilon_1^{n+1} + \varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2 \right). $$

Using Lemma 3.4 and the Poincaré inequality we have

$$ \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle \le \frac{1}{4} \|u^N\|_{L^2(\Omega_f)}^2 + \frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \\ + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ + C \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \|b_2^{n+1}\|_{L^2(\Omega_f)}^2 + C \frac{1}{\alpha^2} \|\varepsilon_2^N\|_{L^2(\Omega_f)}^2. $$

It follows from (3.3) and the definition of $Z^0$ that $Z^0(w, u, \lambda) = 0$. We finish the proof by combining the above estimates. □

Due to Remark 3.2, the discrete time derivative of $u$ and $w$ solves (3.1) and (3.2) with discrete time derivatives of the data as the right-hand sides, we immediately get the following.

**Corollary 3.6.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_{\Delta t} \varepsilon_2^{n+1} \in V_f$ then

$$ Z^N(\partial_{\Delta t} w, \partial_{\Delta t} u, \partial_{\Delta t} \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} w, \partial_{\Delta t} u, \partial_{\Delta t} \lambda) \le C \Xi^N(\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) $$

Similarly, the second-order discrete time derivative of $u$ and $w$ solves (3.1) and (3.2) with the second-order discrete time derivatives of the data as the right-hand sides, we also have the following.

**Corollary 3.7.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_{\Delta t}^2 \varepsilon_2^{n+1} \in V_f$ then

$$ Z^N(\partial_{\Delta t}^2 w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2 \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2 \lambda) \le C \Xi^N(\partial_{\Delta t}^2 b_1, \partial_{\Delta t}^2 b_2, \partial_{\Delta t}^2 \varepsilon_1, \partial_{\Delta t}^2 \varepsilon_2) $$

**Remark 3.8.** We can also analyze the problem (3.1)-(3.2) but replacing the homogeneous Neumann boundary conditions with homogeneous Dirichlet boundary conditions. In other words, the problem with pure Dirichlet boundary conditions that is $\Gamma_D^i = \partial \Omega_i \setminus \Sigma$ for $i = s, f$. In this case exactly the same estimates hold. Now of course, we assume $\varepsilon_2^n \in V_f$ where $V_f$ has zero Dirichlet boundary conditions on $\Gamma_D^f = \partial \Omega_f \setminus \Sigma$.

**3.1. H**$^2$ stability in a special case. In this section we prove $H^2$ estimates for $u$ in a special configuration. Again, we assume that $\Omega = (0, 1)^2$, and now assume $\Sigma$ is parallel to the x-axis (see Figure 2). We take advantage of the fact that the sides composing $\Gamma_N$ are perpendicular to the x-axis and, moreover, we also notice that $\Sigma$ is parallel to the x-axis.

Then, in this particular case we have, for $n \ge 0$,

$$ (3.11a) \qquad \partial_{\Delta t} \partial_x w^{n+1} - \nu_s \Delta \partial_x w^{n+1} = \partial_x b_1^{n+1}, \qquad \text{on } \Omega_s, $$

$$ (3.11b) \qquad \partial_x w^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.11c) \qquad \partial_x w^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.11d) \qquad \alpha \partial_x w^{n+1} + \nu_s \partial_{n_f} \partial_x w^{n+1} = \alpha \partial_x u^n - \nu_f \partial_{n_f} \partial_x u^n + \partial_x \varepsilon_1^{n+1} \qquad \text{on } \Sigma. $$
Using the Poincaré and trace inequalities, we easily obtain the following bound

$$ \Delta t \sum_{n=0}^{N-1} F^{n+1}(w, u) \le \frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|b_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{\nu_f} \|b_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\varepsilon_1^{n+1} + \varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2 \right). $$

Using Lemma 3.4 and the Poincaré inequality we have

$$ \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \langle \varepsilon_2^{n+1}, \lambda^{n+1} \rangle \le \frac{1}{4} \|u^N\|_{L^2(\Omega_f)}^2 + \frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \\ + C \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ + C \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \|b_2^{n+1}\|_{L^2(\Omega_f)}^2 + C \frac{1}{\alpha^2} \|\varepsilon_2^N\|_{L^2(\Omega_f)}^2. $$

It follows from (3.3) and the definition of $Z^0$ that $Z^0(w, u, \lambda) = 0$. We finish the proof by combining the above estimates. □

Due to Remark 3.2, the discrete time derivative of $u$ and $w$ solves (3.1) and (3.2) with discrete time derivatives of the data as the right-hand sides, we immediately get the following.

**Corollary 3.6.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_{\Delta t} \varepsilon_2^{n+1} \in V_f$ then

$$ Z^N(\partial_{\Delta t} w, \partial_{\Delta t} u, \partial_{\Delta t} \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} w, \partial_{\Delta t} u, \partial_{\Delta t} \lambda) \le C \Xi^N(\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) $$

Similarly, the second-order discrete time derivative of $u$ and $w$ solves (3.1) and (3.2) with the second-order discrete time derivatives of the data as the right-hand sides, we also have the following.

**Corollary 3.7.** Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_{\Delta t}^2 \varepsilon_2^{n+1} \in V_f$ then

$$ Z^N(\partial_{\Delta t}^2 w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2 \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2 \lambda) \le C \Xi^N(\partial_{\Delta t}^2 b_1, \partial_{\Delta t}^2 b_2, \partial_{\Delta t}^2 \varepsilon_1, \partial_{\Delta t}^2 \varepsilon_2) $$

*Remark 3.8.* We can also analyze the problem (3.1)-(3.2) but replacing the homogeneous Neumann boundary conditions with homogeneous Dirichlet boundary conditions. In other words, the problem with pure Dirichlet boundary conditions that is $\Gamma_D^i = \partial \Omega_i \setminus \Sigma$ for $i = s, f$. In this case exactly the same estimates hold. Now of course, we assume $\varepsilon_2^n \in V_f$ where $V_f$ has zero Dirichlet boundary conditions on $\Gamma_D^f = \partial \Omega_f \setminus \Sigma$.

### 3.1. $H^2$ stability in a special case.
In this section we prove $H^2$ estimates for $u$ in a special configuration. Again, we assume that $\Omega = (0, 1)^2$, and now assume $\Sigma$ is parallel to the $x$-axis (see Figure 2). We take advantage of the fact that the sides composing $\Gamma_N$ are perpendicular to the $x$-axis and, moreover, we also notice that $\Sigma$ is parallel to the $x$-axis.

Then, in this particular case we have, for $n \ge 0$,

$$ (3.11a) \qquad \partial_{\Delta t} \partial_x w^{n+1} - \nu_s \Delta \partial_x w^{n+1} = \partial_x b_1^{n+1}, \qquad \text{on } \Omega_s, $$

$$ (3.11b) \qquad \partial_x w^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.11c) \qquad \partial_x w^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.11d) \qquad \alpha \partial_x w^{n+1} + \nu_s \partial_{n_f} \partial_x w^{n+1} = \alpha \partial_x u^n - \nu_f \partial_{n_f} \partial_x u^n + \partial_x \varepsilon_1^{n+1} \qquad \text{on } \Sigma. $$
<figure><img src="image_5.png" /><figcaption>FIGURE 2. The domains $\Omega_f$ and $\Omega_s$ with horizontal interface $\Sigma$.</figcaption></figure>

$$ (3.12a) \qquad \partial_{\Delta t} \partial_x u^{n+1} - \nu_f \Delta \partial_x u^{n+1} = \partial_x b_2^{n+1}, \qquad \text{in } \Omega_f, $$

$$ (3.12b) \qquad \partial_x u^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.12c) \qquad \partial_x u^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.12d) \qquad \alpha \partial_x u^{n+1} + \nu_f \partial_{n_f} \partial_x u^{n+1} = \alpha \partial_x w^{n+1} + \nu_f \partial_{n_f} \partial_x u^n + \partial_x \varepsilon_2^{n+1} \quad \text{on } \Sigma, $$

with

$$ (3.13) \qquad \partial_x w^0 = \partial_x u^0 = \partial_x b_i^0 = \partial_x \varepsilon_i^0 = 0 \quad i = 1, 2, $$

according to (3.3).

We then get an immediate Corollary from Remark 3.8 and Theorem 3.5.

**Corollary 3.9.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_x \varepsilon_2^n \in \{v \in H^1(\Omega_f) : v = 0 \text{ on } \partial\Omega_f \setminus \Sigma\}$ on $\Gamma_N^f$ then the following estimate holds

$$ Z^N(\partial_x w, \partial_x u, \partial_x \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_x w, \partial_x u, \partial_x \lambda) \leq C \Xi^N(\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2). $$

**Corollary 3.10.** Under the hypothesis of Corollaries 3.9 and 3.6 we have

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2 u^{n+1}\|_{L^2(\Omega_f)}^2 \leq C \left( \Xi^N(\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) + \Xi^N(\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2) \right) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2. $$

*Proof.* From Corollary 3.9 we get

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \| \nabla \partial_x u^{n+1} \|_{L^2(\Omega_f)}^2 \leq C \Xi^N (\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2). $$
<figure><img src="image_5.png" /><figcaption>FIGURE 2. The domains $\Omega_f$ and $\Omega_s$ with horizontal interface $\Sigma$.</figcaption></figure>

$$ (3.12a) \qquad \partial_{\Delta t} \partial_x u^{n+1} - \nu_f \Delta \partial_x u^{n+1} = \partial_x b_2^{n+1}, \qquad \text{in } \Omega_f, $$

$$ (3.12b) \qquad \partial_x u^{n+1} = 0, \qquad \text{on } \Gamma_D^s, $$

$$ (3.12c) \qquad \partial_x u^{n+1} = 0, \qquad \text{on } \Gamma_N^s, $$

$$ (3.12d) \qquad \alpha \partial_x u^{n+1} + \nu_f \partial_{n_f} \partial_x u^{n+1} = \alpha \partial_x w^{n+1} + \nu_f \partial_{n_f} \partial_x u^n + \partial_x \varepsilon_2^{n+1} \quad \text{on } \Sigma, $$

with

$$ (3.13) \qquad \partial_x w^0 = \partial_x u^0 = \partial_x b_i^0 = \partial_x \varepsilon_i^0 = 0 \quad i = 1, 2, $$

according to (3.3).

We then get an immediate Corollary from Remark 3.8 and Theorem 3.5.

**Corollary 3.9.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $w, u$ solve (3.1) and (3.2) and assuming that $\partial_x \varepsilon_2^n \in \{v \in H^1(\Omega_f) : v = 0 \text{ on } \partial\Omega_f \setminus \Sigma\}$ on $\Gamma_N^f$ then the following estimate holds

$$ Z^N(\partial_x w, \partial_x u, \partial_x \lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_x w, \partial_x u, \partial_x \lambda) \le C \Xi^N(\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2). $$

**Corollary 3.10.** Under the hypothesis of Corollaries 3.9 and 3.6 we have

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2 u^{n+1}\|_{L^2(\Omega_f)}^2 \le C \left( \Xi^N(\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) + \Xi^N(\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2) \right) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2. $$

*Proof.* From Corollary 3.9 we get

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|\nabla \partial_x u^{n+1}\|_{L^2(\Omega_f)}^2 \le C \Xi^N(\partial_x b_1, \partial_x b_2, \partial_x \varepsilon_1, \partial_x \varepsilon_2). $$
Moreover, using (3.2a) and Corollary 3.6, we have

$$
\begin{align*}
\Delta t \sum_{n=0}^{N-1} \nu_f \|\Delta u^{n+1}\|_{L^2(\Omega_f)}^2 &= \Delta t \sum_{n=0}^{N-1} \nu_f \|\partial_{\Delta t} u^{n+1} - b_1^{n+1}\|_{L^2(\Omega_f)}^2 \\
&\le C \Xi^N (\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) + 2 \Delta t \sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2.
\end{align*}
$$

Finally, using the following estimate,

$$
\Delta t \sum_{n=0}^{N-1} \nu_f \| \partial_y^2 u^{n+1} \|_{L^2(\Omega_f)}^2 \le 2 \Delta t \sum_{n=0}^{N-1} \nu_f \left( \| \partial_x^2 u^{n+1} \|_{L^2(\Omega_f)}^2 + \| \Delta u^{n+1} \|_{L^2(\Omega_f)}^2 \right),
$$

and combining the above estimates we obtain the result.

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

4. ERROR ESTIMATES OF THE ROBIN-ROBIN METHOD

In this section we apply the stability results of the previous sections to obtain error estimates of
the Robin-Robin splitting method (2.4)-(2.5) applied to (2.1). We use the following notation for the
errors :

$$
U^n := u^n - u^n, \quad W^n := w^n - w^n, \quad \Lambda^n := |n| - \lambda^n.
$$

We use the convention $u^j = u^j = u_0$ and $w^j = w^j = w_0$ for $j < 0$.

Then the error equations read, for $n \ge 0$,

$$
(4.1a) \qquad \partial_{\Delta t} W^{n+1} - \nu_s \Delta W^{n+1} = -h_1^{n+1}, \qquad \text{on } \Omega_s,
$$

$$
(4.1b) \qquad W^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(4.1c) \qquad \partial_n W^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(4.1d) \qquad \alpha W^{n+1} + \nu_s \partial_{n_s} W^{n+1} = \alpha U^n - \nu_f \partial_{n_f} U^n + \alpha g_1^{n+1} - g_2^{n+1} \quad \text{on } \Sigma.
$$

$$
(4.2a) \qquad \partial_{\Delta t} U^{n+1} - \nu_f \Delta U^{n+1} = -h_2^{n+1}, \qquad \text{in } \Omega_f,
$$

$$
(4.2b) \qquad U^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(4.2c) \qquad \partial_n U^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(4.2d) \qquad \alpha U^{n+1} + \nu_f \partial_{n_f} U^{n+1} = \alpha W^{n+1} + \nu_f \partial_{n_f} U^n + g_2^{n+1} \quad \text{on } \Sigma.
$$

where

$$
h_1^{n+1} := \partial_t w^{n+1} - \partial_{\Delta t} w^{n+1}, \quad g_1^{n+1} := u^{n+1} - u^n,
$$

$$
h_2^{n+1} := \partial_t u^{n+1} - \partial_{\Delta t} u^{n+1}, \quad g_2^{n+1} := |n+1| - |n|.
$$

The equations (4.1)-(4.2) are well-defined and we also have

$$
(4.3) \qquad W^0 = U^0 = h_i^0 = g_i^0 = 0, \quad i = 1, 2,
$$

$$
(4.4) \qquad \partial_{\Delta t} W^0 = \partial_{\Delta t} U^0 = \partial_{\Delta t} h_i^0 = \partial_{\Delta t} g_i^0 = 0, \quad i = 1, 2,
$$

$$
(4.5) \qquad \partial_{\Delta t}^2 W^0 = \partial_{\Delta t}^2 U^0 = \partial_{\Delta t}^2 h_i^0 = \partial_{\Delta t}^2 g_i^0 = 0, \quad i = 1, 2.
$$

We then extend $|\cdot|$ to $\Omega_f$ in a natural way. We let $\tilde{|} = \phi \nu_f \nabla \cdot \mathbf{n}_f$ where $\phi$ is a function that is one on $\Sigma$ and vanishes on $\Gamma_D^f$. Then, we define $\tilde{g}_2^{n+1} = |\tilde{n}^{n+1} - |\tilde{n}|^n$. By construction $\tilde{g}_2^{n+1} \in V_f$ and $\tilde{g}_2^n = g_2^n$ on $\Sigma$. We immediately get the following result if we apply Theorem 3.5 and Corollary 3.6.
Moreover, using (3.2a) and Corollary 3.6, we have

$$
\begin{align*}
\Delta t \sum_{n=0}^{N-1} \nu_f \|\Delta u^{n+1}\|_{L^2(\Omega_f)}^2 &= \Delta t \sum_{n=0}^{N-1} \nu_f \|\partial_{\Delta t} u^{n+1} - b_1^{n+1}\|_{L^2(\Omega_f)}^2 \\
&\le C \Xi^N (\partial_{\Delta t} b_1, \partial_{\Delta t} b_2, \partial_{\Delta t} \varepsilon_1, \partial_{\Delta t} \varepsilon_2) + 2 \Delta t \sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2.
\end{align*}
$$

Finally, using the following estimate,

$$
\Delta t \sum_{n=0}^{N-1} \nu_f \| \partial_y^2 u^{n+1} \|_{L^2(\Omega_f)}^2 \le 2 \Delta t \sum_{n=0}^{N-1} \nu_f ( \| \partial_x^2 u^{n+1} \|_{L^2(\Omega_f)}^2 + \| \Delta u^{n+1} \|_{L^2(\Omega_f)}^2 ),
$$

and combining the above estimates we obtain the result.

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

4. ERROR ESTIMATES OF THE ROBIN-ROBIN METHOD

In this section we apply the stability results of the previous sections to obtain error estimates of
the Robin-Robin splitting method (2.4)-(2.5) applied to (2.1). We use the following notation for the
errors :

$$
U^n := u^n - u^n, \quad W^n := w^n - w^n, \quad \Lambda^n := |n| - \lambda^n.
$$

We use the convention $u^j = u^j = u_0$ and $w^j = w^j = w_0$ for $j < 0$.

Then the error equations read, for $n \ge 0$,

$$
(4.1a) \qquad \partial_{\Delta t} W^{n+1} - \nu_s \Delta W^{n+1} = -h_1^{n+1}, \qquad \text{on } \Omega_s,
$$

$$
(4.1b) \qquad W^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(4.1c) \qquad \partial_n W^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(4.1d) \qquad \alpha W^{n+1} + \nu_s \partial_{n_s} W^{n+1} = \alpha U^n - \nu_f \partial_{n_f} U^n + \alpha g_1^{n+1} - g_2^{n+1} \quad \text{on } \Sigma.
$$

$$
(4.2a) \qquad \partial_{\Delta t} U^{n+1} - \nu_f \Delta U^{n+1} = -h_2^{n+1}, \qquad \text{in } \Omega_f,
$$

$$
(4.2b) \qquad U^{n+1} = 0, \qquad \text{on } \Gamma_D^s,
$$

$$
(4.2c) \qquad \partial_n U^{n+1} = 0, \qquad \text{on } \Gamma_N^s,
$$

$$
(4.2d) \qquad \alpha U^{n+1} + \nu_f \partial_{n_f} U^{n+1} = \alpha W^{n+1} + \nu_f \partial_{n_f} U^n + g_2^{n+1} \quad \text{on } \Sigma.
$$

where

$$
h_1^{n+1} := \partial_t w^{n+1} - \partial_{\Delta t} w^{n+1}, \quad g_1^{n+1} := u^{n+1} - u^n,
$$

$$
h_2^{n+1} := \partial_t u^{n+1} - \partial_{\Delta t} u^{n+1}, \quad g_2^{n+1} := |n+1| - |n|.
$$

The equations (4.1)-(4.2) are well-defined and we also have

$$
(4.3) \qquad W^0 = U^0 = h_i^0 = g_i^0 = 0, \quad i = 1, 2,
$$

$$
(4.4) \qquad \partial_{\Delta t} W^0 = \partial_{\Delta t} U^0 = \partial_{\Delta t} h_i^0 = \partial_{\Delta t} g_i^0 = 0, \quad i = 1, 2,
$$

$$
(4.5) \qquad \partial_{\Delta t}^2 W^0 = \partial_{\Delta t}^2 U^0 = \partial_{\Delta t}^2 h_i^0 = \partial_{\Delta t}^2 g_i^0 = 0, \quad i = 1, 2.
$$

We then extend $|\cdot|$ to $\Omega_f$ in a natural way. We let $\tilde{|} = \phi \nu_f \nabla \cdot n_f$ where $\phi$ is a function that is one on $\Sigma$ and vanishes on $\Gamma_D^f$. Then, we define $\tilde{g}_2^{n+1} = |\tilde{n}^{n+1} - |\tilde{n}|^n$. By construction $\tilde{g}_2^{n+1} \in V_f$ and $\tilde{g}_2^n = g_2^n$ on $\Sigma$. We immediately get the following result if we apply Theorem 3.5 and Corollary 3.6.
**Corollary 4.1.** Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ (4.6) \qquad Z^N(W, U, \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \leq C \Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2), $$

$$ (4.7) \qquad \begin{aligned} & Z^N(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) \\ & \leq C \Xi^N(\partial_{\Delta t} h_1, \partial_{\Delta t} h_2, \alpha \partial_{\Delta t} g_1 - \partial_{\Delta t} g_2, \partial_{\Delta t} \tilde{g}_2), \end{aligned} $$

and

$$ (4.8) \qquad \begin{aligned} & Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \\ & \leq C \Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \alpha \partial_{\Delta t}^2 g_1 - \partial_{\Delta t}^2 g_2, \partial_{\Delta t}^2 \tilde{g}_2). \end{aligned} $$

Then, it is quite straightforward to get a convergence rate by estimating the right-hand sides. The proof of the following Corollary can be found in Appendix A.

**Corollary 4.2.** Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ (4.9) \qquad Z^N(W, U, \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \leq C(\Delta t)^2 Y, $$

$$ (4.10) \qquad Z^N(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) \leq C(\Delta t)^2 Y, $$

and

$$ (4.11) \qquad Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \leq C(\Delta t)^2 Z, $$

where $Y$ is defined as

$$ Y := \frac{1}{\nu_s} \|\partial_t^2 w\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^2 u\|_{L^2(0,T;L^2(\Omega_f))}^2 $$

$$ (4.12) \qquad \begin{aligned} & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t u\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t u\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t u\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t l\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t u\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$

$Y$ is defined as

$$ Y := \frac{1}{\nu_s} \|\partial_t^3 w\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^3 u\|_{L^2(0,T;L^2(\Omega_f))}^2 $$

$$ (4.13) \qquad \begin{aligned} & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^2 u\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^2 u\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^2 u\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^2 l\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^2 u\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$
**Corollary 4.1.** Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ (4.6) \qquad Z^N(W, U, \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \leq C \Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2), $$

$$ (4.7) \qquad \begin{aligned} & Z^N(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) \\ & \leq C \Xi^N(\partial_{\Delta t} h_1, \partial_{\Delta t} h_2, \alpha \partial_{\Delta t} g_1 - \partial_{\Delta t} g_2, \partial_{\Delta t} \tilde{g}_2), \end{aligned} $$

and

$$ (4.8) \qquad \begin{aligned} & Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \\ & \leq C \Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \alpha \partial_{\Delta t}^2 g_1 - \partial_{\Delta t}^2 g_2, \partial_{\Delta t}^2 \tilde{g}_2). \end{aligned} $$

Then, it is quite straightforward to get a convergence rate by estimating the right-hand sides. The proof of the following Corollary can be found in Appendix A.

**Corollary 4.2.** Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ (4.9) \qquad Z^N(W, U, \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \leq C(\Delta t)^2 Y, $$

$$ (4.10) \qquad Z^N(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t} W, \partial_{\Delta t} U, \partial_{\Delta t} \Lambda) \leq C(\Delta t)^2 Y, $$

and

$$ (4.11) \qquad Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) + \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \leq C(\Delta t)^2 \mathfrak{Y}, $$

where $Y$ is defined as

$$ Y := \frac{1}{\nu_s} \|\partial_t^2 w\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^2 u\|_{L^2(0,T;L^2(\Omega_f))}^2 $$

$$ (4.12) \qquad \begin{aligned} & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t u\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t u\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t u\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t l\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t u\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$

$Y$ is defined as

$$ Y := \frac{1}{\nu_s} \|\partial_t^3 w\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^3 u\|_{L^2(0,T;L^2(\Omega_f))}^2 $$

$$ (4.13) \qquad \begin{aligned} & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^2 u\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^2 u\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^2 u\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^2 l\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^2 u\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$
and $\mathfrak{Y}$ is defined as

$$ (4.14) \qquad \begin{aligned} \mathfrak{Y} := & \frac{1}{\nu_s} \|\partial_t^4 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^4 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^3 \mathbf{I}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{aligned} $$

4.1. $H^2$ error estimates in a special case. Here we assume that we have the configuration as in Figure 1. Then, we see that $\partial_x \tilde{g}_2^{n+1}$ vanishes on $\Gamma_N^f$ and $\Gamma_D^f$ and hence belongs to $V_f$. Hence, Corollary 3.10 gives the following corollary.

**Corollary 4.3.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(U^{n+1})\|_{L^2(\Omega_f)}^2 \le C \Xi^N (\partial_{\Delta t} h_1, \partial_{\Delta t} h_2, \partial_{\Delta t}(\alpha g_1 - g_2), \partial_{\Delta t} \tilde{g}_2) \\ + C \Xi^N (\partial_x h_1, \partial_x h_2, \partial_x(\alpha g_1 - g_2), \partial_{\Delta t} \partial_x \tilde{g}_2) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|h_1^{n+1}\|_{L^2(\Omega_f)}^2. $$

Since $\partial_{\Delta t} U$ and $\partial_{\Delta t} W$ satisfy the same equations as $U, W$ with time difference right-hand sides, we have

**Corollary 4.4.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5), then

$$ (4.15) \qquad \begin{aligned} \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(\partial_{\Delta t} U^{n+1})\|_{L^2(\Omega_f)}^2 \le C \Xi^N (\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2(\alpha g_1 - g_2), \partial_{\Delta t}^2 \tilde{g}_2) \\ + C \Xi^N (\partial_{\Delta t} \partial_x h_1, \partial_{\Delta t} \partial_x h_2, \partial_{\Delta t} \partial_x(\alpha g_1 - g_2), \partial_{\Delta t} \partial_x \tilde{g}_2) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|\partial_{\Delta t} h_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{aligned} $$

Then, it is straight-forward to obtain the convergence rate by estimate the right-hand sides. See Appendix B for a proof for the following Corollary.

**Corollary 4.5.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(\partial_{\Delta t} U^{n+1})\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^2 (\mathfrak{Y} + \mathbb{Y} + \nu_f \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2) $$

where $\mathfrak{Y}$ is defined in Corollary 4.2 and $\mathbb{Y}$ is defined as,

$$ (4.16) \qquad \begin{aligned} \mathbb{Y} := & \frac{1}{\nu_s} \|\partial_x \partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t^2 \mathbf{I}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{aligned} $$
and $\mathfrak{Y}$ is defined as

$$ (4.14) \qquad \begin{aligned} \mathfrak{Y} := & \frac{1}{\nu_s} \|\partial_t^4 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^4 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^3 \mathbf{I}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{aligned} $$

4.1. $H^2$ error estimates in a special case. Here we assume that we have the configuration as in Figure 1. Then, we see that $\partial_x \tilde{g}_2^{n+1}$ vanishes on $\Gamma_N^f$ and $\Gamma_D^f$ and hence belongs to $V_f$. Hence, Corollary 3.10 gives the following corollary.

**Corollary 4.3.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(U^{n+1})\|_{L^2(\Omega_f)}^2 \le C \Xi^N (\partial_{\Delta t} h_1, \partial_{\Delta t} h_2, \partial_{\Delta t}(\alpha g_1 - g_2), \partial_{\Delta t} \tilde{g}_2) \\ + C \Xi^N (\partial_x h_1, \partial_x h_2, \partial_x(\alpha g_1 - g_2), \partial_{\Delta t} \partial_x \tilde{g}_2) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|h_1^{n+1}\|_{L^2(\Omega_f)}^2. $$

Since $\partial_{\Delta t} U$ and $\partial_{\Delta t} W$ satisfy the same equations as $U, W$ with time difference right-hand sides, we have

**Corollary 4.4.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5), then

$$ (4.15) \qquad \begin{aligned} \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(\partial_{\Delta t} U^{n+1})\|_{L^2(\Omega_f)}^2 \le C \Xi^N (\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2(\alpha g_1 - g_2), \partial_{\Delta t}^2 \tilde{g}_2) \\ + C \Xi^N (\partial_{\Delta t} \partial_x h_1, \partial_{\Delta t} \partial_x h_2, \partial_{\Delta t} \partial_x(\alpha g_1 - g_2), \partial_{\Delta t} \partial_x \tilde{g}_2) \\ + C \Delta t \sum_{n=0}^{N-1} \nu_f \|\partial_{\Delta t} h_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{aligned} $$

Then, it is straight-forward to obtain the convergence rate by estimate the right-hand sides. See Appendix B for a proof for the following Corollary.

**Corollary 4.5.** Suppose that $\Sigma$ is perpendicular to the two sides of $\Gamma_N$ as in Figure 2. Let $u, w$ solve (2.1) and $w, u$ solve (2.4) and (2.5) then

$$ \Delta t \sum_{n=0}^{N-1} \nu_f \|D^2(\partial_{\Delta t} U^{n+1})\|_{L^2(\Omega_f)}^2 \le C (\Delta t)^2 (\mathfrak{Y} + \mathbb{Y} + \nu_f \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2) $$

where $\mathfrak{Y}$ is defined in Corollary 4.2 and $\mathbb{Y}$ is defined as,

$$ (4.16) \qquad \begin{aligned} \mathbb{Y} := & \frac{1}{\nu_s} \|\partial_x \partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t^2 \mathbf{I}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{aligned} $$
## 5. NUMERICAL EXPERIMENTS

In this section we provide numerical experiments that agree with our theoretical results. We let

$$e_u = \|U^N\|_{L^2(\Omega_f)}, \quad e_{1,u} = \|U^N - U^{N-1}\|_{L^2(\Omega_f)},$$

$$e_{2,u} = \|(U^N - U^{N-1}) - (U^{N-1} - U^{N-2})\|_{L^2(\Omega_f)},$$

$$e_{\lambda} = \|\Lambda^N\|_{L^2(\Sigma)}, \quad e_{1,\lambda} = \|\Lambda^N - \Lambda^{N-1}\|_{L^2(\Sigma)}$$

$$e_{1,u,2} = \|U^N - U^{N-1}\|_{H^2(\Omega_f)}.$$

All numerical experiments are performed using FEniCS and multiphenics [1, 3]. Although we only analyze the semi-discrete method, here we present the results for a fully discrete method where we use the piecewise linear finite element method for the spatial discretization, except for the computation of $e_{1,u,2}$, where we use the piecewise quadratic finite element method because piecewise linear function do not approximation a function well in $H^2$. In addition, we also present convergence rates for the Lagrange multiplier.

*Example 5.1.* We consider the domain $\Omega = (0,1)^2$, $\Omega_f = (0,1) \times (0, .75)$ and $\Omega_s = (0,1) \times (.75, 1)$. See Figure 2 for an illustration. We take $\nu_f = 1 = \nu_s$ and take the solution of (2.3) to be

$$w = u = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

We take $h = \Delta t$, $T = 0.25$ and $\alpha = 4$ where $h$ is the mesh size of the triangulation.

<table><caption>TABLE 1. Errors and convergence rates of $U^N$ for Example 5.1</caption><thead><tr><th>Δt</th><th>e<sub>u</sub></th><th>rates</th><th>e<sub>1,u</sub></th><th>rates</th><th>e<sub>2,u</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>7.65e-02</td><td>-</td><td>7.73e-02</td><td>-</td><td>6.57e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>3.72e-02</td><td>1.04</td><td>5.87e-02</td><td>0.40</td><td>1.55e-01</td><td>8.73</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.74e-02</td><td>1.10</td><td>1.47e-02</td><td>1.99</td><td>7.91e-03</td><td>4.29</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.95e-03</td><td>1.13</td><td>3.58e-03</td><td>2.04</td><td>1.27e-03</td><td>2.64</td></tr><tr><td>(1/2)<sup>6</sup></td><td>3.52e-03</td><td>1.17</td><td>8.41e-04</td><td>2.09</td><td>1.75e-04</td><td>2.85</td></tr><tr><td>(1/2)<sup>7</sup></td><td>1.62e-03</td><td>1.12</td><td>1.96e-04</td><td>2.10</td><td>2.15e-05</td><td>3.03</td></tr><tr><td>(1/2)<sup>8</sup></td><td>7.70e-04</td><td>1.07</td><td>4.69e-05</td><td>2.06</td><td>2.62e-06</td><td>3.04</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.75e-04</td><td>1.04</td><td>1.14e-05</td><td>2.04</td><td>3.22e-07</td><td>3.02</td></tr></tbody></table>

<table><caption>TABLE 2. Error and convergence rates of $\Lambda^N$ for Example 5.1</caption><thead><tr><th>Δt</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>1,u,2</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.55e-01</td><td>-</td><td>2.55e-01</td><td>-</td><td>1.33e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>9.73e-02</td><td>1.39</td><td>2.41e-01</td><td>0.08</td><td>1.11e+01</td><td>0.26</td></tr><tr><td>(1/2)<sup>4</sup></td><td>3.06e-02</td><td>1.67</td><td>4.41e-02</td><td>2.45</td><td>2.61e-01</td><td>2.08</td></tr><tr><td>(1/2)<sup>5</sup></td><td>1.62e-02</td><td>0.92</td><td>6.69e-03</td><td>2.72</td><td>5.76e-02</td><td>2.18</td></tr><tr><td>(1/2)<sup>6</sup></td><td>9.20e-03</td><td>0.82</td><td>2.06e-03</td><td>1.70</td><td>1.40e-02</td><td>2.04</td></tr><tr><td>(1/2)<sup>7</sup></td><td>4.55e-03</td><td>1.01</td><td>5.28e-04</td><td>1.97</td><td>3.32e-03</td><td>2.07</td></tr><tr><td>(1/2)<sup>8</sup></td><td>2.23e-03</td><td>1.03</td><td>1.31e-04</td><td>2.01</td><td>8.03e-04</td><td>2.05</td></tr><tr><td>(1/2)<sup>9</sup></td><td>1.10e-03</td><td>1.02</td><td>3.26e-05</td><td>2.01</td><td>1.97e-04</td><td>2.03</td></tr></tbody></table>

As we can see from Tables 1-2, the $L^2$ error at the final step, $e_u$ is of order $(\Delta t)$ whereas the difference of two consecutive errors, $e_{1,u}$ is $(\Delta t)^2$ and the second difference, $e_{2,u}$ is $(\Delta t)^3$. The $L^2$ error of the Lagrange multiplier at the final time, $e_\lambda$ is of order $\Delta t$ and the difference $e_{1,\lambda}$, is of order $(\Delta t)^2$. It is also clear that the $H^2$ error of the difference of $U^N$, $e_{1,u,2}$ is of order $(\Delta t)^2$ as we proved in Corollary 4.5.
## 5. NUMERICAL EXPERIMENTS

In this section we provide numerical experiments that agree with our theoretical results. We let

$$e_u = \|U^N\|_{L^2(\Omega_f)}, \quad e_{1,u} = \|U^N - U^{N-1}\|_{L^2(\Omega_f)},$$

$$e_{2,u} = \|(U^N - U^{N-1}) - (U^{N-1} - U^{N-2})\|_{L^2(\Omega_f)},$$

$$e_{\lambda} = \|\Lambda^N\|_{L^2(\Sigma)}, \quad e_{1,\lambda} = \|\Lambda^N - \Lambda^{N-1}\|_{L^2(\Sigma)}$$

$$e_{1,u,2} = \|U^N - U^{N-1}\|_{H^2(\Omega_f)}.$$

All numerical experiments are performed using FEniCS and multiphenics [1, 3]. Although we only analyze the semi-discrete method, here we present the results for a fully discrete method where we use the piecewise linear finite element method for the spatial discretization, except for the computation of $e_{1,u,2}$, where we use the piecewise quadratic finite element method because piecewise linear function do not approximation a function well in $H^2$. In addition, we also present convergence rates for the Lagrange multiplier.

*Example 5.1.* We consider the domain $\Omega = (0, 1)^2$, $\Omega_f = (0, 1) \times (0, .75)$ and $\Omega_s = (0, 1) \times (.75, 1)$. See Figure 2 for an illustration. We take $\nu_f = 1 = \nu_s$ and take the solution of (2.3) to be

$$w = u = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

We take $h = \Delta t$, $T = 0.25$ and $\alpha = 4$ where $h$ is the mesh size of the triangulation.

<table><caption>TABLE 1. Errors and convergence rates of $U^N$ for Example 5.1</caption><thead><tr><th>Δt</th><th>e<sub>u</sub></th><th>rates</th><th>e<sub>1,u</sub></th><th>rates</th><th>e<sub>2,u</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>7.65e-02</td><td>-</td><td>7.73e-02</td><td>-</td><td>6.57e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>3.72e-02</td><td>1.04</td><td>5.87e-02</td><td>0.40</td><td>1.55e-01</td><td>8.73</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.74e-02</td><td>1.10</td><td>1.47e-02</td><td>1.99</td><td>7.91e-03</td><td>4.29</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.95e-03</td><td>1.13</td><td>3.58e-03</td><td>2.04</td><td>1.27e-03</td><td>2.64</td></tr><tr><td>(1/2)<sup>6</sup></td><td>3.52e-03</td><td>1.17</td><td>8.41e-04</td><td>2.09</td><td>1.75e-04</td><td>2.85</td></tr><tr><td>(1/2)<sup>7</sup></td><td>1.62e-03</td><td>1.12</td><td>1.96e-04</td><td>2.10</td><td>2.15e-05</td><td>3.03</td></tr><tr><td>(1/2)<sup>8</sup></td><td>7.70e-04</td><td>1.07</td><td>4.69e-05</td><td>2.06</td><td>2.62e-06</td><td>3.04</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.75e-04</td><td>1.04</td><td>1.14e-05</td><td>2.04</td><td>3.22e-07</td><td>3.02</td></tr></tbody></table>

<table><caption>TABLE 2. Error and convergence rates of $\Lambda^N$ for Example 5.1</caption><thead><tr><th>Δt</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>1,u,2</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.55e-01</td><td>-</td><td>2.55e-01</td><td>-</td><td>1.33e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>9.73e-02</td><td>1.39</td><td>2.41e-01</td><td>0.08</td><td>1.11e+01</td><td>0.26</td></tr><tr><td>(1/2)<sup>4</sup></td><td>3.06e-02</td><td>1.67</td><td>4.41e-02</td><td>2.45</td><td>2.61e-01</td><td>2.08</td></tr><tr><td>(1/2)<sup>5</sup></td><td>1.62e-02</td><td>0.92</td><td>6.69e-03</td><td>2.72</td><td>5.76e-02</td><td>2.18</td></tr><tr><td>(1/2)<sup>6</sup></td><td>9.20e-03</td><td>0.82</td><td>2.06e-03</td><td>1.70</td><td>1.40e-02</td><td>2.04</td></tr><tr><td>(1/2)<sup>7</sup></td><td>4.55e-03</td><td>1.01</td><td>5.28e-04</td><td>1.97</td><td>3.32e-03</td><td>2.07</td></tr><tr><td>(1/2)<sup>8</sup></td><td>2.23e-03</td><td>1.03</td><td>1.31e-04</td><td>2.01</td><td>8.03e-04</td><td>2.05</td></tr><tr><td>(1/2)<sup>9</sup></td><td>1.10e-03</td><td>1.02</td><td>3.26e-05</td><td>2.01</td><td>1.97e-04</td><td>2.03</td></tr></tbody></table>

As we can see from Tables 1-2, the $L^2$ error at the final step, $e_u$ is of order $(\Delta t)$ whereas the difference of two consecutive errors, $e_{1,u}$ is $(\Delta t)^2$ and the second difference, $e_{2,u}$ is $(\Delta t)^3$. The $L^2$ error of the Lagrange multiplier at the final time, $e_\lambda$ is of order $\Delta t$ and the difference $e_{1,\lambda}$, is of order $(\Delta t)^2$. It is also clear that the $H^2$ error of the difference of $U^N$, $e_{1,u,2}$ is of order $(\Delta t)^2$ as we proved in Corollary 4.5.
*Example 5.2.* In this example, we test our algorithm for a non-horizontal interface problem. We consider the domain $\Omega = (0, 1)^2$ and we let $\Sigma$ be defined as the straight line connecting $(0, 0.25)$ and $(1, 0.75)$. We then define $\Omega_s$ as the region above $\Sigma$ and $\Omega_f$ as the region below $\Sigma$. We take $\nu_f = 1 = \nu_s$ and take the solution of (2.3) to be

$$w = u = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

Other parameters are identical to those of Example 5.1.

We report the convergence results in Tables 3-4. We again observe expected convergence rates for both $U^N$ and $\Lambda^N$. It indicates that our methods also work for a more general interface problem.

<table><caption>TABLE 3. Errors and convergence rates of $U^N$ for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>u</sub></th><th>rates</th><th>e<sub>1,u</sub></th><th>rates</th><th>e<sub>2,u</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>8.05e-02</td><td>-</td><td>9.64e-02</td><td>-</td><td>4.87e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.10e-02</td><td>0.66</td><td>4.23e-02</td><td>1.19</td><td>1.36e-01</td><td>8.48</td></tr><tr><td>(1/2)<sup>4</sup></td><td>2.46e-02</td><td>1.05</td><td>1.57e-02</td><td>1.43</td><td>4.94e-03</td><td>4.78</td></tr><tr><td>(1/2)<sup>5</sup></td><td>9.20e-03</td><td>1.42</td><td>4.07e-03</td><td>1.95</td><td>1.47e-03</td><td>1.74</td></tr><tr><td>(1/2)<sup>6</sup></td><td>3.52e-03</td><td>1.38</td><td>8.46e-04</td><td>2.27</td><td>1.82e-04</td><td>3.02</td></tr><tr><td>(1/2)<sup>7</sup></td><td>1.52e-03</td><td>1.22</td><td>1.86e-04</td><td>2.18</td><td>2.11e-05</td><td>3.11</td></tr><tr><td>(1/2)<sup>8</sup></td><td>7.00e-04</td><td>1.11</td><td>4.33e-05</td><td>2.10</td><td>2.49e-06</td><td>3.08</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.36e-04</td><td>1.06</td><td>1.04e-05</td><td>2.06</td><td>3.02e-07</td><td>3.04</td></tr></tbody></table>

<table><caption>TABLE 4. Error and convergence rates of $\Lambda^N$ for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>1,u,2</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>8.51e-01</td><td>-</td><td>8.54e-01</td><td>-</td><td>2.04e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.01e-01</td><td>0.76</td><td>3.81e-01</td><td>1.16</td><td>9.67e-01</td><td>1.08</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.94e-01</td><td>1.37</td><td>1.45e-01</td><td>1.39</td><td>3.66e-01</td><td>1.40</td></tr><tr><td>(1/2)<sup>5</sup></td><td>5.34e-02</td><td>1.86</td><td>2.44e-02</td><td>2.58</td><td>9.08e-02</td><td>2.01</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.84e-02</td><td>1.54</td><td>4.38e-03</td><td>2.48</td><td>1.86e-02</td><td>2.29</td></tr><tr><td>(1/2)<sup>7</sup></td><td>7.49e-03</td><td>1.30</td><td>9.10e-04</td><td>2.27</td><td>4.04e-03</td><td>2.20</td></tr><tr><td>(1/2)<sup>8</sup></td><td>3.39e-03</td><td>1.14</td><td>2.07e-04</td><td>2.13</td><td>9.31e-04</td><td>2.11</td></tr><tr><td>(1/2)<sup>9</sup></td><td>1.61e-03</td><td>1.07</td><td>4.95e-05</td><td>2.07</td><td>2.23e-04</td><td>2.06</td></tr></tbody></table>

## 6. CONCLUDING REMARKS

We analyzed the Robin-Robin coupling methods [11] for parabolic\parabolic interface problems and proved higher convergence rates in time for the first-order and second-order discrete time derivatives. We also prove $H^2$ estimates of the discrete time derivatives in a special case. All the estimates in this work are key ingredients in proving that a prediction correction method [12] produces a $O((\Delta t)^2)$ convergence rate.

## ACKNOWLEDGEMENT

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.
*Example 5.2.* In this example, we test our algorithm for a non-horizontal interface problem. We consider the domain $\Omega = (0, 1)^2$ and we let $\Sigma$ be defined as the straight line connecting $(0, 0.25)$ and $(1, 0.75)$. We then define $\Omega_s$ as the region above $\Sigma$ and $\Omega_f$ as the region below $\Sigma$. We take $\nu_f = 1 = \nu_s$ and take the solution of (2.3) to be

$$w = u = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

Other parameters are identical to those of Example 5.1.

We report the convergence results in Tables 3-4. We again observe expected convergence rates for both $U^N$ and $\Lambda^N$. It indicates that our methods also work for a more general interface problem.

<table><caption>TABLE 3. Errors and convergence rates of $U^N$ for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>u</sub></th><th>rates</th><th>e<sub>1,u</sub></th><th>rates</th><th>e<sub>2,u</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>8.05e-02</td><td>-</td><td>9.64e-02</td><td>-</td><td>4.87e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.10e-02</td><td>0.66</td><td>4.23e-02</td><td>1.19</td><td>1.36e-01</td><td>8.48</td></tr><tr><td>(1/2)<sup>4</sup></td><td>2.46e-02</td><td>1.05</td><td>1.57e-02</td><td>1.43</td><td>4.94e-03</td><td>4.78</td></tr><tr><td>(1/2)<sup>5</sup></td><td>9.20e-03</td><td>1.42</td><td>4.07e-03</td><td>1.95</td><td>1.47e-03</td><td>1.74</td></tr><tr><td>(1/2)<sup>6</sup></td><td>3.52e-03</td><td>1.38</td><td>8.46e-04</td><td>2.27</td><td>1.82e-04</td><td>3.02</td></tr><tr><td>(1/2)<sup>7</sup></td><td>1.52e-03</td><td>1.22</td><td>1.86e-04</td><td>2.18</td><td>2.11e-05</td><td>3.11</td></tr><tr><td>(1/2)<sup>8</sup></td><td>7.00e-04</td><td>1.11</td><td>4.33e-05</td><td>2.10</td><td>2.49e-06</td><td>3.08</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.36e-04</td><td>1.06</td><td>1.04e-05</td><td>2.06</td><td>3.02e-07</td><td>3.04</td></tr></tbody></table>

<table><caption>TABLE 4. Error and convergence rates of $\Lambda^N$ for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>1,u,2</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>8.51e-01</td><td>-</td><td>8.54e-01</td><td>-</td><td>2.04e+01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.01e-01</td><td>0.76</td><td>3.81e-01</td><td>1.16</td><td>9.67e-01</td><td>1.08</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.94e-01</td><td>1.37</td><td>1.45e-01</td><td>1.39</td><td>3.66e-01</td><td>1.40</td></tr><tr><td>(1/2)<sup>5</sup></td><td>5.34e-02</td><td>1.86</td><td>2.44e-02</td><td>2.58</td><td>9.08e-02</td><td>2.01</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.84e-02</td><td>1.54</td><td>4.38e-03</td><td>2.48</td><td>1.86e-02</td><td>2.29</td></tr><tr><td>(1/2)<sup>7</sup></td><td>7.49e-03</td><td>1.30</td><td>9.10e-04</td><td>2.27</td><td>4.04e-03</td><td>2.20</td></tr><tr><td>(1/2)<sup>8</sup></td><td>3.39e-03</td><td>1.14</td><td>2.07e-04</td><td>2.13</td><td>9.31e-04</td><td>2.11</td></tr><tr><td>(1/2)<sup>9</sup></td><td>1.61e-03</td><td>1.07</td><td>4.95e-05</td><td>2.07</td><td>2.23e-04</td><td>2.06</td></tr></tbody></table>

## 6. CONCLUDING REMARKS

We analyzed the Robin-Robin coupling methods [11] for parabolic\parabolic interface problems and proved higher convergence rates in time for the first-order and second-order discrete time derivatives. We also prove $H^2$ estimates of the discrete time derivatives in a special case. All the estimates in this work are key ingredients in proving that a prediction correction method [12] produces a $O((\Delta t)^2)$ convergence rate.

## ACKNOWLEDGEMENT

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.
## APPENDIX A. SKETCH OF PROOF: COROLLARY 4.2

*Proof of (4.9) in Corollary 4.2.* To prove (4.9), it is sufficient to bound the term $\Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2)$. Note that since $g_2^{n+1} = \tilde{g}_2^{n+1}$ on $\Sigma$, we have

$$ \begin{aligned} \Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2) := & \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|h_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|h_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\alpha g_1\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\tilde{g}_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|\tilde{g}_2^N\|_{L^2(\Omega_f)}^2 \\ = & T_1 + T_2 + \dots + T_8. \end{aligned} $$

All the terms in $\Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2)$ can be easily estimated by (A.2b) except the $T_3$ and $T_4$. For $T_3$, it follows from (A.5) that,

$$ \begin{aligned} \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 = & (\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{I}^{n+1}\|_{L^2(\Omega_f)}^2 \\ \le & C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^2 \tilde{I}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ \le & C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^2 u\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{aligned} \quad (A.1) $$

The term $T_4$ can be bounded as follows,

$$ \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C \Delta t \frac{(\nu_f)^3}{\alpha^2} \sum_{n=0}^{N-1} \left( \|\nabla(u^{n+1} - u^n)\|_{L^2(\Omega_f)}^2 + \|D^2(u^{n+1} - u^n)\|_{L^2(\Omega_f)}^2 \right). $$

Here $D^2u$ denotes the Hessian of $u$. Using (A.2a) we get

$$ \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C (\Delta t)^2 \frac{(\nu_f)^3}{\alpha^2} \|\partial_t u\|_{L^2(0,T;H^2(\Omega_f))}^2. $$

The estimates above imply (4.9). $\square$

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

Before we prove the estimate (4.10), we need the following preliminary results. We use the Bochner norms $\|v\|_{L^2(a,b;X)} = \left( \int_a^b \|v(\cdot, s)\|_X^2 ds \right)^{1/2}$ and $\|v\|_{L^\infty(a,b;X)} = \sup_{a \le s \le b} \|v(\cdot, s)\|_X$. For a Sobolev space $X$, it is well known that

$$ (A.2a) \qquad \|v^{n+1} - v^n\|_X^2 \le C \Delta t \int_{t_n}^{t_{n+1}} \|\partial_t v(\cdot, s)\|_X^2 ds, $$

$$ (A.2b) \qquad \|\partial_{\Delta t} v^{n+1} - \partial_t v^{n+1}\|_X^2 \le C \Delta t \int_{t_n}^{t_{n+1}} \|\partial_t^2 v(\cdot, s)\|_X^2 ds, $$

$$ (A.2c) \qquad \int_a^b \|v(\cdot, s)\|_X^2 ds \le (b-a) \|v\|_{L^\infty(a,b;X)}^2. $$

The following identities can easily be shown

$$ (A.3) \qquad \partial_{\Delta t}^2 v^n = \frac{1}{(\Delta t)^2} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \partial_t^2 v(\cdot, t_n + s) ds, $$
## APPENDIX A. SKETCH OF PROOF: COROLLARY 4.2

*Proof of (4.9) in Corollary 4.2.* To prove (4.9), it is sufficient to bound the term $\Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2)$. Note that since $g_2^{n+1} = \tilde{g}_2^{n+1}$ on $\Sigma$, we have

$$ \begin{aligned} \Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2) := & \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|h_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|h_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\alpha g_1\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\tilde{g}_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|\tilde{g}_2^N\|_{L^2(\Omega_f)}^2 \\ = & T_1 + T_2 + \dots + T_8. \end{aligned} $$

All the terms in $\Xi^N(h_1, h_2, \alpha g_1 - g_2, \tilde{g}_2)$ can be easily estimated by (A.2b) except the $T_3$ and $T_4$. For $T_3$, it follows from (A.5) that,

$$ \begin{aligned} \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 = & (\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{I}^{n+1}\|_{L^2(\Omega_f)}^2 \\ \le & C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^2 \tilde{I}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ \le & C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^2 u\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{aligned} \quad (A.1) $$

The term $T_4$ can be bounded as follows,

$$ \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C \Delta t \frac{(\nu_f)^3}{\alpha^2} \sum_{n=0}^{N-1} \left( \|\nabla(u^{n+1} - u^n)\|_{L^2(\Omega_f)}^2 + \|D^2(u^{n+1} - u^n)\|_{L^2(\Omega_f)}^2 \right). $$

Here $D^2u$ denotes the Hessian of $u$. Using (A.2a) we get

$$ \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1} \|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C (\Delta t)^2 \frac{(\nu_f)^3}{\alpha^2} \|\partial_t u\|_{L^2(0,T;H^2(\Omega_f))}^2. $$

The estimates above imply (4.9). □

Before we prove the estimate (4.10), we need the following preliminary results. We use the Bochner norms $\|v\|_{L^2(a,b;X)} = \left( \int_a^b \|v(\cdot, s)\|_X^2 ds \right)^{1/2}$ and $\|v\|_{L^\infty(a,b;X)} = \sup_{a \le s \le b} \|v(\cdot, s)\|_X$. For a Sobolev space $X$, it is well known that

$$ (A.2a) \qquad \|v^{n+1} - v^n\|_X^2 \le C \Delta t \int_{t_n}^{t_{n+1}} \|\partial_t v(\cdot, s)\|_X^2 ds, $$

$$ (A.2b) \qquad \|\partial_{\Delta t} v^{n+1} - \partial_t v^{n+1}\|_X^2 \le C \Delta t \int_{t_n}^{t_{n+1}} \|\partial_t^2 v(\cdot, s)\|_X^2 ds, $$

$$ (A.2c) \qquad \int_a^b \|v(\cdot, s)\|_X^2 ds \le (b-a) \|v\|_{L^\infty(a,b;X)}^2. $$

The following identities can easily be shown

$$ (A.3) \qquad \partial_{\Delta t}^2 v^n = \frac{1}{(\Delta t)^2} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \partial_t^2 v(\cdot, t_n + s) ds, $$
and

$$
\begin{align}
\frac{\partial^2 v^n - \partial_{\Delta t}^2 v^{n-1}}{\Delta t} &= \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|)(\partial_t^2 v(\cdot, t_n + s) - \partial_t^2 v(\cdot, t_{n-1} + s))ds \nonumber \\
&= \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-1}+s}^{t_{n}+s} \partial_t^3 v(\cdot, r) dr ds. \tag{A.4}
\end{align}
$$

From these we can show that

$$
(A.5) \qquad \| \partial_{\Delta t}^2 v^n \|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} \| \partial_t^2 v(\cdot, s) \|_X^2 ds,
$$

and

$$
(A.6) \qquad \left\| \frac{\partial_{\Delta t}^2 v^n - \partial_{\Delta t}^2 v^{n-1}}{\Delta t} \right\|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-2}}^{t_{n+1}} \| \partial_t^3 v(\cdot, r) \|_X^2 dr.
$$

Proof of (4.10) in Corollary 4.2. It is similar to prove (4.10). Indeed, let us define, for $j = 1, 2$

$$
G_j^{n+1} = \partial_{\Delta t} g_j^{n+1}, \quad \tilde{G}_2^{n+1} = \partial_{\Delta t} \tilde{g}_2^{n+1}, \quad H_j^{n+1} = \partial_{\Delta t} h_j^{n+1}.
$$

We then need to bound the term

$$
\begin{align*}
\Xi^N(H_1, H_2, \alpha G_1 - G_2, \tilde{G}_2)
&= \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|H_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|H_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
&\quad + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
&\quad + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\alpha G_1\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\tilde{G}_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|\tilde{G}_2^N\|_{L^2(\Omega_f)}^2 \\
&= R_1 + R_2 + \dots + R_8.
\end{align*}
$$

For $R_1$, it follows from the definitions of $H_1^{n+1}$ and $h_1^{n+1}$ that,

$$
\begin{align*}
\| H_1^{n+1} \|_{L^2(\Omega_s)}^2 &= \int_{\Omega_s} \left( \frac{h_1^{n+1} - h_1^n}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_s} \left( \frac{\partial_t w^{n+1} - \partial_{\Delta t} w^{n+1} - \partial_t w^n + \partial_{\Delta t} w^n}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_s} (\partial_{\Delta t} (\partial_t w)^{n+1} - \partial_{\Delta t}^2 w^n)^2 dx \\
(A.7) \qquad \qquad &= \int_{\Omega_s} (\partial_{\Delta t} (\partial_t w)^{n+1} - \partial_t^2 w^n + \partial_t^2 w^n - \partial_{\Delta t}^2 w^{n+1})^2 dx \\
&\le C \int_{\Omega_s} \left( \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} (\partial_t^2 w(t_n) - \partial_t^2 w(s)) ds \right)^2 dx \\
&\quad + C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^2} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s-t_n|) (\partial_t^2 w(s) - \partial_t^2 w(t_n)) ds \right)^2 dx \\
&\le C \Delta t \| \partial_t^3 w \|_{L^2((t_{n-1}, t_{n+1}), L^2(\Omega_s))}^2.
\end{align*}
$$
and

$$
\begin{align}
\frac{\partial^2 v^n - \partial_{\Delta t}^2 v^{n-1}}{\Delta t} &= \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|)(\partial_t^2 v(\cdot, t_n + s) - \partial_t^2 v(\cdot, t_{n-1} + s))ds \nonumber \\
&= \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-1}+s}^{t_{n}+s} \partial_t^3 v(\cdot, r) dr ds. \tag{A.4}
\end{align}
$$

From these we can show that

$$
(A.5) \qquad \| \partial_{\Delta t}^2 v^n \|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} \| \partial_t^2 v(\cdot, s) \|_X^2 ds,
$$

and

$$
(A.6) \qquad \left\| \frac{\partial_{\Delta t}^2 v^n - \partial_{\Delta t}^2 v^{n-1}}{\Delta t} \right\|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-2}}^{t_{n+1}} \| \partial_t^3 v(\cdot, r) \|_X^2 dr.
$$

*Proof of (4.10) in Corollary 4.2.* It is similar to prove (4.10). Indeed, let us define, for $j = 1, 2$

$$
G_j^{n+1} = \partial_{\Delta t} g_j^{n+1}, \quad \tilde{G}_2^{n+1} = \partial_{\Delta t} \tilde{g}_2^{n+1}, \quad H_j^{n+1} = \partial_{\Delta t} h_j^{n+1}.
$$

We then need to bound the term

$$
\begin{align*}
\Xi^N(H_1, H_2, \alpha G_1 - G_2, \tilde{G}_2)
&= \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|H_1^{n+1}\|_{L^2(\Omega_s)}^2 + \left(\frac{1}{\nu_f} + 1\right) \|H_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
&\quad + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha} \|\tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\
&\quad + \Delta t \sum_{n=0}^{N-1} \left( \frac{1}{\nu_s} \|\alpha G_1\|_{L^2(\Sigma)}^2 + \frac{1}{\nu_f} \|\tilde{G}_2^{n+1}\|_{L^2(\Sigma)}^2 \right) + \frac{1}{\alpha^2} \|\tilde{G}_2^N\|_{L^2(\Omega_f)}^2 \\
&= R_1 + R_2 + \dots + R_8.
\end{align*}
$$

For $R_1$, it follows from the definitions of $H_1^{n+1}$ and $h_1^{n+1}$ that,

$$
\begin{align*}
\| H_1^{n+1} \|_{L^2(\Omega_s)}^2 &= \int_{\Omega_s} \left( \frac{h_1^{n+1} - h_1^n}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_s} \left( \frac{\partial_t w^{n+1} - \partial_{\Delta t} w^{n+1} - \partial_t w^n + \partial_{\Delta t} w^n}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_s} (\partial_{\Delta t} (\partial_t w)^{n+1} - \partial_{\Delta t}^2 w^n)^2 dx \\
(A.7) \qquad \qquad \qquad \qquad &= \int_{\Omega_s} (\partial_{\Delta t} (\partial_t w)^{n+1} - \partial_t^2 w^n + \partial_t^2 w^n - \partial_{\Delta t}^2 w^{n+1})^2 dx \\
&\le C \int_{\Omega_s} \left( \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} (\partial_t^2 w(t_n) - \partial_t^2 w(s)) ds \right)^2 dx \\
&\quad + C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^2} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s-t_n|) (\partial_t^2 w(s) - \partial_t^2 w(t_n)) ds \right)^2 dx \\
&\le C \Delta t \| \partial_t^3 w \|_{L^2((t_{n-1}, t_{n+1}), L^2(\Omega_s))}^2.
\end{align*}
$$
Therefore, we obtain

$$
(A.8) \qquad \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_s} \|H_1^{n+1}\|_{L^2(\Omega_s)}^2 \le C(\Delta t)^2 \frac{1}{\nu_s} \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2
$$

The estimate of $R_2$ is similar to that of $R_1$. The term $R_3$ can be estimated as follows by (A.6),

$$
\begin{align*}
\Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= (\Delta t)^3 \sum_{n=2}^{N-1} \frac{1}{\nu_f \alpha^2} \left\| \frac{\partial_{\Delta t}^2 \tilde{l}^n - \partial_{\Delta t}^2 \tilde{l}^{n-1}}{\Delta t} \right\|_{L^2(\Omega_f)}^2 \\
(A.9) \qquad &\le C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^3 \tilde{l}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\
&\le C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2.
\end{align*}
$$

For $R_4$, it follows from the definition of $\tilde{l}$ and the fact $\|\nabla\phi\|_{L^2(\Omega_f)} \le C$ that,

$$
\begin{align*}
\|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= \int_{\Omega_f} \left( \frac{\nabla^{n+1} - 2\nabla^{n} + \nabla^{n-1}}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_f} \left( \frac{1}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) \partial_t^2 \nabla \tilde{l}(s) ds \right)^2 dx \\
(A.10) \qquad &\le C \int_{\Omega_f} \left( \int_{t_{n-1}}^{t_{n+1}} |\partial_t^2 \nabla \tilde{l}(s)| ds \right)^2 dx \\
&\le C \Delta t \int_{\Omega_f} \int_{t_{n-1}}^{t_{n+1}} (\partial_t^2 \nabla \tilde{l}(s))^2 ds dx \\
&\le C \nu_f^2 \Delta t \|\partial_t^2 \mathbf{u}\|_{L^2(t_{n-1}, t_{n+1}; H^2(\Omega_f))}^2,
\end{align*}
$$

and hence,

$$
(A.11) \quad \Delta t \sum_{n=0}^{N-1} \frac{\nu_f}{\alpha^2} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^2 \frac{(\nu_f^3)}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2.
$$

The remaining terms $R_5$ to $R_8$ can be estimated similarly, we give the estimate of $R_6$ here:

$$
\begin{align*}
\|G_1^{n+1}\|_{L^2(\Sigma)}^2 &= \int_{\Sigma} \left( \frac{g_1^{n+1} - g_1^n}{\Delta t} \right)^2 d\sigma \\
&= \int_{\Sigma} \left( \frac{u^{n+1} - 2u^n + u^{n-1}}{\Delta t} \right)^2 d\sigma \\
(A.12) \qquad &= \int_{\Sigma} \left( \frac{1}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) \partial_t^2 u(s) ds \right)^2 d\sigma \\
&\le C \Delta t \|\partial_t^2 u\|_{L^2((t_{n-1}, t_{n+1}), L^2(\Sigma))}^2,
\end{align*}
$$

and therefore, we obtain

$$
(A.13) \quad \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_s} \|\alpha G_1\|_{L^2(\Sigma)}^2 \le C \frac{\alpha^2}{\nu_s} \|\partial_t^2 \mathbf{u}\|_{L^2((0,T), L^2(\Sigma))}^2.
$$

We finish the proof by combining all the estimates above.

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

In order to prove (4.11), we define the following quantities for $j = 1, 2$,

$$
(A.14) \qquad G_j^{n+1} = \partial_{\Delta t} G_j^{n+1}, \quad \mathcal{H}_j^{n+1} = \partial_{\Delta t} H_j^{n+1}, \quad \tilde{G}_2^{n+1} = \partial_{\Delta t} \tilde{G}_2^{n+1}.
$$
Therefore, we obtain

$$
(A.8) \qquad \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_s} \|H_1^{n+1}\|_{L^2(\Omega_s)}^2 \le C(\Delta t)^2 \frac{1}{\nu_s} \|\partial_t^3 w\|_{L^2((0,T),L^2(\Omega_s))}^2
$$

The estimate of $R_2$ is similar to that of $R_1$. The term $R_3$ can be estimated as follows by (A.6),

$$
\begin{align*}
\Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= (\Delta t)^3 \sum_{n=2}^{N-1} \frac{1}{\nu_f \alpha^2} \left\| \frac{\partial_{\Delta t}^2 \tilde{l}^n - \partial_{\Delta t}^2 \tilde{l}^{n-1}}{\Delta t} \right\|_{L^2(\Omega_f)}^2 \\
(A.9) \qquad &\le C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^3 \tilde{l}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\
&\le C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2.
\end{align*}
$$

For $R_4$, it follows from the definition of $\tilde{l}$ and the fact $\|\nabla\phi\|_{L^2(\Omega_f)} \le C$ that,

$$
\begin{align*}
\|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= \int_{\Omega_f} \left( \frac{\nabla^{n+1} - 2\nabla^{n} + \nabla^{n-1}}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_f} \left( \frac{1}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) \partial_t^2 \nabla \tilde{l}(s) ds \right)^2 dx \\
(A.10) \qquad &\le C \int_{\Omega_f} \left( \int_{t_{n-1}}^{t_{n+1}} |\partial_t^2 \nabla \tilde{l}(s)| ds \right)^2 dx \\
&\le C \Delta t \int_{\Omega_f} \int_{t_{n-1}}^{t_{n+1}} (\partial_t^2 \nabla \tilde{l}(s))^2 ds dx \\
&\le C \nu_f^2 \Delta t \|\partial_t^2 \mathbf{u}\|_{L^2(t_{n-1}, t_{n+1}; H^2(\Omega_f))}^2,
\end{align*}
$$

and hence,

$$
(A.11) \quad \Delta t \sum_{n=0}^{N-1} \frac{\nu_f}{\alpha^2} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C (\Delta t)^2 \frac{(\nu_f^3)}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2.
$$

The remaining terms $R_5$ to $R_8$ can be estimated similarly, we give the estimate of $R_6$ here:

$$
\begin{align*}
\|G_1^{n+1}\|_{L^2(\Sigma)}^2 &= \int_{\Sigma} \left( \frac{g_1^{n+1} - g_1^n}{\Delta t} \right)^2 d\sigma \\
&= \int_{\Sigma} \left( \frac{u^{n+1} - 2u^n + u^{n-1}}{\Delta t} \right)^2 d\sigma \\
(A.12) \qquad &= \int_{\Sigma} \left( \frac{1}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) \partial_t^2 u(s) ds \right)^2 d\sigma \\
&\le C \Delta t \|\partial_t^2 u\|_{L^2((t_{n-1}, t_{n+1}), L^2(\Sigma))}^2,
\end{align*}
$$

and therefore, we obtain

$$
(A.13) \quad \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_s} \|\alpha G_1\|_{L^2(\Sigma)}^2 \le C \frac{\alpha^2}{\nu_s} \|\partial_t^2 u\|_{L^2((0,T),L^2(\Sigma))}^2.
$$

We finish the proof by combining all the estimates above.

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

In order to prove (4.11), we define the following quantities for $j = 1, 2$,

$$
(A.14) \qquad G_j^{n+1} = \partial_{\Delta t} G_j^{n+1}, \quad \mathcal{H}_j^{n+1} = \partial_{\Delta t} H_j^{n+1}, \quad \tilde{G}_2^{n+1} = \partial_{\Delta t} \tilde{G}_2^{n+1}.
$$
Note that we also have the following

$$
(A.15) \qquad \tilde{G}_2^{n+1} = \frac{\tilde{g}_2^{n+1} - 2\tilde{g}_2^n + \tilde{g}_2^{n-1}}{(\Delta t)^2} = \frac{\tilde{I}^{n+1} - 3\tilde{I}^n + 3\tilde{I}^{n-1} - \tilde{I}^{n-2}}{(\Delta t)^2}.
$$

Denote

$$
(A.16) \qquad \partial_{\Delta t}^3 v^{n+1} = \frac{v^{n+1} - 3v^n + 3v^{n-1} - v^{n-2}}{(\Delta t)^3},
$$

thus we see that

$$
(A.17) \qquad \frac{\tilde{G}_2^{n+1} - \tilde{G}_2^n}{(\Delta t)^2} = \frac{\partial_{\Delta t}^3 \tilde{I}^{n+1} - \partial_{\Delta t}^3 \tilde{I}^n}{\Delta t}.
$$

Moreover, we see that

$$
\begin{align*}
\frac{\partial_{\Delta t}^3 v^n - \partial_{\Delta t}^3 v^{n-1}}{\Delta t} &= \frac{1}{(\Delta t)^4} \left( \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_n + s) - \partial_t^2 v(\cdot, t_{n-1} + s)) ds \right. \\
&\qquad \left. - \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_{n-1} + s) - \partial_t^2 v(\cdot, t_{n-2} + s)) ds \right) \\
&= \frac{1}{(\Delta t)^4} \left( \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_n + s) - 2\partial_t^2 v(\cdot, t_{n-1} + s) + \partial_t^2 v(\cdot, t_{n-2} + s)) ds \right) \\
&= \frac{1}{(\Delta t)^4} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-2}+s}^{t_{n+}s} (\Delta t - |r - t_{n-1}|) \partial_t^4 v(\cdot, r) dr ds.
\end{align*}
$$

Therefore, we obtain the following estimate

$$
(A.18) \qquad \left\| \frac{\partial_{\Delta t}^3 v^n - \partial_{\Delta t}^3 v^{n-1}}{\Delta t} \right\|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-3}}^{t_{n+1}} \left\| \partial_t^4 v(\cdot, r) \right\|_X^2 dr.
$$

Proof of (4.11) in Corollary 4.2. We now bound the term $\Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2 (\alpha g_1 - g_2), \partial_{\Delta t}^2 \tilde{g}_2)$ which is $\Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \mathcal{G}_1 - \mathcal{G}_2), \tilde{\mathcal{G}}_2)$ according to (A.14). The techniques are similar to those mentioned above. We present the key estimates involving $\mathcal{H}_1$ and $\tilde{\mathcal{G}}_2$ first. It follows from the definition of $\mathcal{H}_1^{n+1}$ that,

$$
\begin{align*}
\| \mathcal{H}_1^{n+1} \|_{L^2(\Omega_s)}^2 &= \int_{\Omega_s} \left( \frac{h_1^{n+1} - 2h_1^n + h_1^{n-1}}{(\Delta t)^2} \right)^2 dx \\
&= \int_{\Omega_s} \left( \frac{\partial_{\Delta t} (\partial_t w)^{n+1} - \partial_{\Delta t}^2 w^{n+1} - (\partial_{\Delta t} (\partial_t w)^n - \partial_{\Delta t}^2 w^n)}{\Delta t} \right)^2 dx \\
&= \int_{\Omega_s} (\partial_{\Delta t}^2 (\partial_t w)^{n+1} - \partial_{\Delta t}^3 w^{n+1})^2 dx \\
(A.19) \qquad &= \int_{\Omega_s} (\partial_{\Delta t}^2 (\partial_t w)^{n+1} - \partial_t^3 w^n + \partial_t^3 w^n - \partial_{\Delta t}^3 w^{n+1})^2 dx \\
&\le C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^2} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) (\partial_t^3 w(s) - \partial_t^3 w(t_n)) ds \right)^2 dx \\
&\quad + C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-1}+s}^{t_{n+}s} \partial_t^3 w(t_n) - \partial_t^3 w(\cdot, r) ds \right)^2 dx \\
&\le C \Delta t \| \partial_t^4 w \|_{L^2((t_{n-2}, t_{n+1}), L^2(\Omega_s))}^2.
\end{align*}
$$
Note that we also have the following

$$ (A.15) \qquad \tilde{G}_2^{n+1} = \frac{\tilde{g}_2^{n+1} - 2\tilde{g}_2^n + \tilde{g}_2^{n-1}}{(\Delta t)^2} = \frac{\tilde{I}^{n+1} - 3\tilde{I}^n + 3\tilde{I}^{n-1} - \tilde{I}^{n-2}}{(\Delta t)^2}. $$

Denote

$$ (A.16) \qquad \partial_{\Delta t}^3 v^{n+1} = \frac{v^{n+1} - 3v^n + 3v^{n-1} - v^{n-2}}{(\Delta t)^3}, $$

thus we see that

$$ (A.17) \qquad \frac{\tilde{G}_2^{n+1} - \tilde{G}_2^n}{(\Delta t)^2} = \frac{\partial_{\Delta t}^3 \tilde{I}^{n+1} - \partial_{\Delta t}^3 \tilde{I}^n}{\Delta t}. $$

Moreover, we see that

$$ \begin{aligned} \frac{\partial_{\Delta t}^3 v^n - \partial_{\Delta t}^3 v^{n-1}}{\Delta t} &= \frac{1}{(\Delta t)^4} \left( \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_n + s) - \partial_t^2 v(\cdot, t_{n-1} + s)) ds \right. \\ & \qquad \left. - \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_{n-1} + s) - \partial_t^2 v(\cdot, t_{n-2} + s)) ds \right) \\ &= \frac{1}{(\Delta t)^4} \left( \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) (\partial_t^2 v(\cdot, t_n + s) - 2\partial_t^2 v(\cdot, t_{n-1} + s) + \partial_t^2 v(\cdot, t_{n-2} + s)) ds \right) \\ &= \frac{1}{(\Delta t)^4} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-2}+s}^{t_{n}+s} (\Delta t - |r - t_{n-1}|) \partial_t^4 v(\cdot, r) dr ds. \end{aligned} $$

Therefore, we obtain the following estimate

$$ (A.18) \qquad \left\| \frac{\partial_{\Delta t}^3 v^n - \partial_{\Delta t}^3 v^{n-1}}{\Delta t} \right\|_{X}^2 \le \frac{C}{\Delta t} \int_{t_{n-3}}^{t_{n+1}} \left\| \partial_t^4 v(\cdot, r) \right\|_{X}^2 dr. $$

*Proof of (4.11) in Corollary 4.2.* We now bound the term $\Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2 (\alpha g_1 - g_2), \partial_{\Delta t}^2 \tilde{g}_2)$ which is $\Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \tilde{G}_1 - \tilde{G}_2), \tilde{G}_2)$ according to (A.14). The techniques are similar to those mentioned above. We present the key estimates involving $\mathcal{H}_1$ and $\tilde{G}_2$ first. It follows from the definition of $\mathcal{H}_1^{n+1}$ that,

$$ \begin{aligned} \|\mathcal{H}_1^{n+1}\|_{L^2(\Omega_s)}^2 &= \int_{\Omega_s} \left( \frac{h_1^{n+1} - 2h_1^n + h_1^{n-1}}{(\Delta t)^2} \right)^2 dx \\ &= \int_{\Omega_s} \left( \frac{\partial_{\Delta t}(\partial_t w)^{n+1} - \partial_{\Delta t}^2 w^{n+1} - (\partial_{\Delta t}(\partial_t w)^n - \partial_{\Delta t}^2 w^n)}{\Delta t} \right)^2 dx \\ &= \int_{\Omega_s} (\partial_{\Delta t}^2 (\partial_t w)^{n+1} - \partial_{\Delta t}^3 w^{n+1})^2 dx \\ &= \int_{\Omega_s} (\partial_{\Delta t}^2 (\partial_t w)^{n+1} - \partial_t^3 w^n + \partial_t^3 w^n - \partial_{\Delta t}^3 w^{n+1})^2 dx \\ &\le C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^2} \int_{t_{n-1}}^{t_{n+1}} (\Delta t - |s - t_n|) (\partial_t^3 w(s) - \partial_t^3 w(t_n)) ds \right)^2 dx \\ &\quad + C \int_{\Omega_s} \left( \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|) \int_{t_{n-1}+s}^{t_{n}+s} \partial_t^3 w(t_n) - \partial_t^3 w(\cdot, r) ds \right)^2 dx \\ &\le C \Delta t \|\partial_t^4 w\|_{L^2((t_{n-2}, t_{n+1}), L^2(\Omega_s))}^2. \end{aligned} $$
We also have, according to the definition of $\tilde{l}$ and the fact $\|\nabla\phi\|_{L^2(\Omega_f)} \le C$ that,

$$ (A.20) \qquad \begin{aligned} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= \int_{\Omega_f} \left( \frac{\nabla^{\tilde{l}^{n+1}-2\nabla^{\tilde{l}^n}+\nabla^{\tilde{l}^{n-1}}}}{\Delta t} - \frac{\nabla^{\tilde{l}^n-2\nabla^{\tilde{l}^{n-1}}+\nabla^{\tilde{l}^{n-2}}}}{\Delta t} \right)^2 dx \\ &= \int_{\Omega_f} \left( \frac{1}{(\Delta t)^2} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|)(\partial_t^2 \nabla \tilde{l}(\cdot, t_n + s) - \partial_t^2 \nabla \tilde{l}(\cdot, t_{n-1} + s)) ds \right)^2 dx \\ &\le C \int_{\Omega_f} \left( \frac{1}{\Delta t} \int_{-\Delta t}^{\Delta t} \int_{t_{n-1}+s}^{t_{n+}s} \partial_t^3 \nabla \tilde{l}(\cdot, r) dr ds \right)^2 dx \\ &\le C \Delta t \int_{\Omega_f} \int_{t_{n-2}}^{t_{n+1}} (\partial_t^3 \nabla \tilde{l}(\cdot, r))^2 dr dx \\ &\le C \nu_f^2 \Delta t \|\partial_t^3 u\|_{L^2(t_{n-2}, t_{n+1}; H^2(\Omega_f))}^2. \end{aligned} $$

It follows from (A.18) that,

$$ (A.21) \qquad \begin{aligned} \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= (\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \left\| \frac{\partial_{\Delta t}^3 \tilde{l}^n - \partial_{\Delta t}^3 \tilde{l}^{n-1}}{\Delta t} \right\|_{L^2(\Omega_f)}^2 \\ &\le C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^4 \tilde{l}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ &\le C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^4 u\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{aligned} $$

Then we follow the same idea as before and obtain the following bound:

$$ (A.22) \qquad \Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \mathcal{G}_1 - \mathcal{G}_2), \tilde{\mathcal{G}}_2) \le C(\Delta t)^2 \mathfrak{Z}, $$

where $\mathfrak{Z}$ is defined in (4.14). $\square$

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

## APPENDIX B. SKETCH OF PROOF: COROLLARY 4.5

*Proof.* According to Corollary 4.4, we need to bound the terms on the right-hand side of the inequality (4.15). The first term in (4.15) is bounded in Corollary 4.2 while the last term in (4.15) can be easily bounded according to (A.7). At last, the analysis to bound the term $\Xi^N(\partial_{\Delta t}\partial_x h_1, \partial_{\Delta t}\partial_x h_2, \partial_{\Delta t}\partial_x(\alpha g_1 - g_2), \partial_{\Delta t}\partial_x \tilde{g}_2)$ is almost identical to that of the term $\Xi^N(H_1, H_2, \alpha G_1 - G_2, \tilde{G}_2)$ except the additional partial derivative which does not affect the techniques. Therefore, we obtain the bound

$$ (B.1) \qquad \Xi^N(\partial_{\Delta t}\partial_x h_1, \partial_{\Delta t}\partial_x h_2, \partial_{\Delta t}\partial_x(\alpha g_1 - g_2), \partial_{\Delta t}\partial_x \tilde{g}_2) \le C(\Delta t)^2 \mathfrak{Y}, $$

where $\mathfrak{Y}$ is defined in (4.16). This finishes the proof. $\square$

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

## REFERENCES

[1] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The fenics project version 1.5. *Archive of numerical software*, 3(100), 2015.

[2] S. Badia, A. Quaini, and A. Quarteroni. Splitting methods based on algebraic factorization for fluid-structure interaction. *SIAM Journal on Scientific Computing*, 30(4):1778–1805, 2008.

[3] F. Ballarin and G. Rozza. https://multiphenics.github.io/, 2016.

[4] J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *J. Comput. Phys.*, 269:108–137, 2014.

[5] M. Beneš. Convergence and stability analysis of heterogeneous time step coupling schemes for parabolic problems. *Appl. Numer. Math.*, 121:198–222, 2017.

[6] M. Beneš, A. Nekvinda, and M. K. Yadav. Multi-time-step domain decomposition method with non-matching grids for parabolic problems. *Appl. Math. Comput.*, 267:571–582, 2015.

[7] K. Böhmer, P. W. Hemker, and H. J. Stetter. The defect correction approach. *Defect Correction Methods: Theory and Applications*, pages 1–32, 1984.

[8] M. Bucelli, M. Geraint Gabriel, A. Quarteroni, G. Gigante, and C. Vergara. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. *Journal of Computational Physics*, 490:112326, 2023.
We also have, according to the definition of $\tilde{l}$ and the fact $\|\nabla\phi\|_{L^2(\Omega_f)} \le C$ that,

$$ (A.20) \qquad \begin{aligned} \|\nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= \int_{\Omega_f} \left( \frac{\nabla^{\tilde{l}^{n+1}-2\nabla^{\tilde{l}^n}+\nabla^{\tilde{l}^{n-1}}}}{\Delta t} - \frac{\nabla^{\tilde{l}^n-2\nabla^{\tilde{l}^{n-1}}+\nabla^{\tilde{l}^{n-2}}}}{\Delta t} \right)^2 dx \\ &= \int_{\Omega_f} \left( \frac{1}{(\Delta t)^2} \int_{-\Delta t}^{\Delta t} (\Delta t - |s|)(\partial_t^2 \nabla \tilde{l}(\cdot, t_n + s) - \partial_t^2 \nabla \tilde{l}(\cdot, t_{n-1} + s)) ds \right)^2 dx \\ &\le C \int_{\Omega_f} \left( \frac{1}{\Delta t} \int_{-\Delta t}^{\Delta t} \int_{t_{n-1}+s}^{t_{n+}s} \partial_t^3 \nabla \tilde{l}(\cdot, r) dr ds \right)^2 dx \\ &\le C \Delta t \int_{\Omega_f} \int_{t_{n-2}}^{t_{n+1}} (\partial_t^3 \nabla \tilde{l}(\cdot, r))^2 dr dx \\ &\le C \nu_f^2 \Delta t \|\partial_t^3 u\|_{L^2(t_{n-2}, t_{n+1}; H^2(\Omega_f))}^2. \end{aligned} $$

It follows from (A.18) that,

$$ (A.21) \qquad \begin{aligned} \Delta t \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \|\partial_{\Delta t} \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 &= (\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f \alpha^2} \left\| \frac{\partial_{\Delta t}^3 \tilde{l}^n - \partial_{\Delta t}^3 \tilde{l}^{n-1}}{\Delta t} \right\|_{L^2(\Omega_f)}^2 \\ &\le C \frac{(\Delta t)^2}{\nu_f \alpha^2} \|\partial_t^4 \tilde{l}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ &\le C (\Delta t)^2 \frac{\nu_f}{\alpha^2} \|\partial_t^4 u\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{aligned} $$

Then we follow the same idea as before and obtain the following bound:

$$ (A.22) \qquad \Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \mathcal{G}_1 - \mathcal{G}_2), \tilde{\mathcal{G}}_2) \le C(\Delta t)^2 \mathfrak{Z}, $$

where $\mathfrak{Z}$ is defined in (4.14).

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

## APPENDIX B. SKETCH OF PROOF: COROLLARY 4.5

*Proof.* According to Corollary 4.4, we need to bound the terms on the right-hand side of the inequality (4.15). The first term in (4.15) is bounded in Corollary 4.2 while the last term in (4.15) can be easily bounded according to (A.7). At last, the analysis to bound the term $\Xi^N(\partial_{\Delta t}\partial_x h_1, \partial_{\Delta t}\partial_x h_2, \partial_{\Delta t}\partial_x(\alpha g_1 - g_2), \partial_{\Delta t}\partial_x \tilde{g}_2)$ is almost identical to that of the term $\Xi^N(H_1, H_2, \alpha G_1 - G_2, \tilde{G}_2)$ except the additional partial derivative which does not affect the techniques. Therefore, we obtain the bound

$$ (B.1) \qquad \Xi^N(\partial_{\Delta t}\partial_x h_1, \partial_{\Delta t}\partial_x h_2, \partial_{\Delta t}\partial_x(\alpha g_1 - g_2), \partial_{\Delta t}\partial_x \tilde{g}_2) \le C(\Delta t)^2 \mathfrak{Y}, $$

where $\mathfrak{Y}$ is defined in (4.16). This finishes the proof.

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

## REFERENCES

[1] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The fenics project version 1.5. *Archive of numerical software*, 3(100), 2015.

[2] S. Badia, A. Quaini, and A. Quarteroni. Splitting methods based on algebraic factorization for fluid-structure interaction. *SIAM Journal on Scientific Computing*, 30(4):1778–1805, 2008.

[3] F. Ballarin and G. Rozza. https://multiphenics.github.io/, 2016.

[4] J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *J. Comput. Phys.*, 269:108–137, 2014.

[5] M. Beneš. Convergence and stability analysis of heterogeneous time step coupling schemes for parabolic problems. *Appl. Numer. Math.*, 121:198–222, 2017.

[6] M. Beneš, A. Nekvinda, and M. K. Yadav. Multi-time-step domain decomposition method with non-matching grids for parabolic problems. *Appl. Math. Comput.*, 267:571–582, 2015.

[7] K. Böhmer, P. W. Hemker, and H. J. Stetter. The defect correction approach. *Defect Correction Methods: Theory and Applications*, pages 1–32, 1984.

[8] M. Bucelli, M. Geraint Gabriel, A. Quarteroni, G. Gigante, and C. Vergara. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. *Journal of Computational Physics*, 490:112326, 2023.
[9] M. Bukač, S. Čanić, R. Glowinski, B. Muha, and A. Quaini. A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. *Internat. J. Numer. Methods Fluids*, 74(8):577–604, 2014.

[10] M. Bukač. An extension of explicit coupling for fluid-structure interaction problems. *Mathematics*, 9(15), 2021.

[11] E. Burman, R. Durst, M. Fernández, and J. Guzmán. Loosely coupled, non-iterative time-splitting scheme based on Robin–Robin coupling: Unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. *Journal of Numerical Mathematics*, 31(1):59–77, 2023.

[12] E. Burman, R. Durst, M. Fernández, J. Guzmán, and S. Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. *Preprint*, 2024.

[13] E. Burman, R. Durst, M. A. Fernández, and J. Guzmán. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. *Numerische Mathematik*, 151(4):807–840, 2022.

[14] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and O. Ruz. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. 2023. https://inria.hal.science/hal-04258861.

[15] E. Burman, R. Durst, and J. Guzmán. Stability and error analysis of a splitting method using Robin–Robin coupling applied to a fluid–structure interaction problem. *Numerical Methods for Partial Differential Equations*, 38(5):1396–1406, 2022.

[16] E. Burman and M. A. Fernández. Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. *Computer Methods in Applied Mechanics and Engineering*, 198(5):766–784, 2009.

[17] E. Burman and M. A. Fernández. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin–Robin coupling. *International Journal for Numerical Methods in Engineering*, 97(10):739–758, 2014.

[18] C. Canuto and A. Lo Giudice. A multi-timestep Robin-Robin domain decomposition method for time dependent advection-diffusion problems. *Appl. Math. Comput.*, 363:124596, 14, 2019.

[19] J. M. Connors and J. S. Howell. A fluid-fluid interaction method using decoupled subproblems and differing time steps. *Numer. Methods Partial Differential Equations*, 28(4):1283–1308, 2012.

[20] J. M. Connors, J. S. Howell, and W. J. Layton. Partitioned time stepping for a parabolic two domain problem. *SIAM Journal on Numerical Analysis*, 47(5):3526–3549, 2009.

[21] J. M. Connors, J. S. Howell, and W. J. Layton. Decoupled time stepping methods for fluid-fluid interaction. *SIAM J. Numer. Anal.*, 50(3):1297–1319, 2012.

[22] R. Durst. *Recent Advances in Splitting Methods Based on Robin-Robin Coupling Conditions*. PhD thesis, Brown University, 2022.

[23] M. A. Fernández, J.-F. Gerbeau, and S. Smaldone. Explicit coupling schemes for a fluid-fluid interaction problem arising in hemodynamics. *SIAM J. Sci. Comput.*, 36(6):A2557–A2583, 2014.

[24] G. Gigante and C. Vergara. On the choice of interface parameters in Robin–Robin loosely coupled schemes for fluid–structure interaction. *Fluids*, 6(6), 2021.

[25] G. Gigante and C. Vergara. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Comput. Math. Appl.*, 96:109–119, 2021.

[26] W. Layton and A. Takhirov. Energy stability of a first order partitioned method for systems with general coupling. *Int. J. Numer. Anal. Model. Ser. B*, 4(3):203–214, 2013.

[27] W. Li, P. Huang, and Y. He. An unconditionally energy stable finite element scheme for a nonlinear fluid-fluid interaction model. *IMA J. Numer. Anal.*, 44(1):157–191, 2024.

[28] P. Minev and R. Usubov. Splitting schemes for the stress formulation of fluid-structure interaction problems. *Applications in Engineering Science*, 9:100082, 2022.

[29] A. Seboldt and M. Bukač. A non-iterative domain decomposition method for the interaction between a fluid and a thick structure. *Numerical Methods for Partial Differential Equations*, 37(4):2803–2832, 2021.

[30] D. A. Serino, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. *SIAM Journal on Scientific Computing*, 41(4):A2464–A2484, 2019.

[31] K. C. Sockwell, P. Bochev, K. Peterson, and P. Kuberry. Interface flux recovery framework for constructing partitioned heterogeneous time-integration methods. *Numer. Methods Partial Differential Equations*, 39(5):3572–3593, 2023.

[32] H. Zhang, Z. Liu, E. Constantinescu, and R. Jacob. Stability analysis of interface conditions for ocean-atmosphere coupling. *J. Sci. Comput.*, 84(3):Paper No. 44, 25, 2020.

[33] Y. Zhang, Y. Hou, and L. Shan. Error estimates of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 333:266–291, 2018.

[34] Y. Zhang, L. Shan, and Y. Hou. New approach to prove the stability of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 371:112695, 19, 2020.
[9] M. Bukač, S. Čanić, R. Glowinski, B. Muha, and A. Quaini. A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. *Internat. J. Numer. Methods Fluids*, 74(8):577–604, 2014.

[10] M. Bukač. An extension of explicit coupling for fluid-structure interaction problems. *Mathematics*, 9(15), 2021.

[11] E. Burman, R. Durst, M. Fernández, and J. Guzmán. Loosely coupled, non-iterative time-splitting scheme based on Robin–Robin coupling: Unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. *Journal of Numerical Mathematics*, 31(1):59–77, 2023.

[12] E. Burman, R. Durst, M. Fernández, J. Guzmán, and S. Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. *Preprint*, 2024.

[13] E. Burman, R. Durst, M. A. Fernández, and J. Guzmán. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. *Numerische Mathematik*, 151(4):807–840, 2022.

[14] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and O. Ruz. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. 2023. https://inria.hal.science/hal-04258861.

[15] E. Burman, R. Durst, and J. Guzmán. Stability and error analysis of a splitting method using Robin–Robin coupling applied to a fluid–structure interaction problem. *Numerical Methods for Partial Differential Equations*, 38(5):1396–1406, 2022.

[16] E. Burman and M. A. Fernández. Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. *Computer Methods in Applied Mechanics and Engineering*, 198(5):766–784, 2009.

[17] E. Burman and M. A. Fernández. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin–Robin coupling. *International Journal for Numerical Methods in Engineering*, 97(10):739–758, 2014.

[18] C. Canuto and A. Lo Giudice. A multi-timestep Robin-Robin domain decomposition method for time dependent advection-diffusion problems. *Appl. Math. Comput.*, 363:124596, 14, 2019.

[19] J. M. Connors and J. S. Howell. A fluid-fluid interaction method using decoupled subproblems and differing time steps. *Numer. Methods Partial Differential Equations*, 28(4):1283–1308, 2012.

[20] J. M. Connors, J. S. Howell, and W. J. Layton. Partitioned time stepping for a parabolic two domain problem. *SIAM Journal on Numerical Analysis*, 47(5):3526–3549, 2009.

[21] J. M. Connors, J. S. Howell, and W. J. Layton. Decoupled time stepping methods for fluid-fluid interaction. *SIAM J. Numer. Anal.*, 50(3):1297–1319, 2012.

[22] R. Durst. *Recent Advances in Splitting Methods Based on Robin-Robin Coupling Conditions*. PhD thesis, Brown University, 2022.

[23] M. A. Fernández, J.-F. Gerbeau, and S. Smaldone. Explicit coupling schemes for a fluid-fluid interaction problem arising in hemodynamics. *SIAM J. Sci. Comput.*, 36(6):A2557–A2583, 2014.

[24] G. Gigante and C. Vergara. On the choice of interface parameters in Robin–Robin loosely coupled schemes for fluid–structure interaction. *Fluids*, 6(6), 2021.

[25] G. Gigante and C. Vergara. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Comput. Math. Appl.*, 96:109–119, 2021.

[26] W. Layton and A. Takhirov. Energy stability of a first order partitioned method for systems with general coupling. *Int. J. Numer. Anal. Model. Ser. B*, 4(3):203–214, 2013.

[27] W. Li, P. Huang, and Y. He. An unconditionally energy stable finite element scheme for a nonlinear fluid-fluid interaction model. *IMA J. Numer. Anal.*, 44(1):157–191, 2024.

[28] P. Minev and R. Usubov. Splitting schemes for the stress formulation of fluid-structure interaction problems. *Applications in Engineering Science*, 9:100082, 2022.

[29] A. Seboldt and M. Bukač. A non-iterative domain decomposition method for the interaction between a fluid and a thick structure. *Numerical Methods for Partial Differential Equations*, 37(4):2803–2832, 2021.

[30] D. A. Serino, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. *SIAM Journal on Scientific Computing*, 41(4):A2464–A2484, 2019.

[31] K. C. Sockwell, P. Bochev, K. Peterson, and P. Kuberry. Interface flux recovery framework for constructing partitioned heterogeneous time-integration methods. *Numer. Methods Partial Differential Equations*, 39(5):3572–3593, 2023.

[32] H. Zhang, Z. Liu, E. Constantinescu, and R. Jacob. Stability analysis of interface conditions for ocean-atmosphere coupling. *J. Sci. Comput.*, 84(3):Paper No. 44, 25, 2020.

[33] Y. Zhang, Y. Hou, and L. Shan. Error estimates of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 333:266–291, 2018.

[34] Y. Zhang, L. Shan, and Y. Hou. New approach to prove the stability of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 371:112695, 19, 2020.
<footnote n="1">¹DEPARTMENT OF MATHEMATICS, UNIVERSITY COLLEGE LONDON, LONDON, UK-WC1E 6BT, UNITED KINGDOM
*Email address: e.burman@ucl.ac.uk*</footnote>

<footnote n="2">²DEPARTMENT OF MATHEMATICS, UNIVERSITY OF PITTSBURGH, PITTSBURGH, PENNSYLVANIA 15261, USA
*Email address: rebecca_durst@alumni.brown.edu*</footnote>

<footnote n="3">³SORBONNE UNIVERSITÉ & CNRS, UMR 7598 LJLL, 75005 PARIS, FRANCE - INRIA, 75012 PARIS, FRANCE
*Email address: miguel.fernandez@inria.fr*</footnote>

<footnote n="4">⁴DIVISION OF APPLIED MATHEMATICS, BROWN UNIVERSITY, 182 GEORGE STREET, PROVIDENCE, RI 02912, USA
*Email address: johnny_guzman@brown.edu*</footnote>

<footnote n="5">⁵THE INSTITUTE FOR COMPUTATIONAL AND EXPERIMENTAL RESEARCH IN MATHEMATICS, BROWN UNIVERSITY, 121 SOUTH MAIN STREET, PROVIDENCE, RI 02912, USA
*Email address: sijing LIU@brown.edu*</footnote>
<footnote n="1">¹DEPARTMENT OF MATHEMATICS, UNIVERSITY COLLEGE LONDON, LONDON, UK-WC1E 6BT, UNITED KINGDOM
*Email address: e.burman@ucl.ac.uk*</footnote>

<footnote n="2">²DEPARTMENT OF MATHEMATICS, UNIVERSITY OF PITTSBURGH, PITTSBURGH, PENNSYLVANIA 15261, USA
*Email address: rebecca_durst@alumni.brown.edu*</footnote>

<footnote n="3">³SORBONNE UNIVERSITÉ & CNRS, UMR 7598 LJLL, 75005 PARIS, FRANCE - INRIA, 75012 PARIS, FRANCE
*Email address: miguel.fernandez@inria.fr*</footnote>

<footnote n="4">⁴DIVISION OF APPLIED MATHEMATICS, BROWN UNIVERSITY, 182 GEORGE STREET, PROVIDENCE, RI 02912, USA
*Email address: johnny_guzman@brown.edu*</footnote>

<footnote n="5">⁵THE INSTITUTE FOR COMPUTATIONAL AND EXPERIMENTAL RESEARCH IN MATHEMATICS, BROWN UNIVERSITY, 121 SOUTH MAIN STREET, PROVIDENCE, RI 02912, USA
*Email address: sijing LIU@brown.edu*</footnote>