%~Mouliné par MaN_auto v.0.36.0 (4a84613e) 2025-10-29 16:00:45
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Simon_20250035.bib}

\usepackage{subcaption}
\usepackage{siunitx}

\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\I}{\mathcal{I}}
\renewcommand{\O}{\mathcal{O}}
\newcommand{\F}{\mathcal{F}}
\renewcommand{\L}{\mathcal{L}}
\newcommand{\Dt}{\Delta t}
\renewcommand{\k}{\underline{Y}}

\newcommand{\eg}{{e.g.}}
\newcommand{\ie}{{i.e.}}

\newcommand{\testeqdeuxdeux}{$2\times2$ coupled test equation}
\newcommand{\matdeuxdeux}{\mathcal{A}}
\newcommand{\cstezerostab}{K}
\newcommand{\onestepoperator}{\varphi}
\newcommand{\phiimpdemo}{\Phi}
\newcommand{\termeun}{e_1}
\newcommand{\termedeux}{e_2}

\newcommand{\exa}{\mathrm{exa}}
\newcommand{\intern}{\mathrm{int}}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\DeclareMathOperator{\Real}{Re}

\makeatletter
\renewcommand\thesubfigure{(\alph{subfigure})}
\DeclareCaptionLabelFormat{mysublabelfmt}{(\alph{subfigure})}
\makeatother
\captionsetup[sub]{labelfont=bf,labelsep=space,labelformat=mysublabelfmt}

%%%---------------------------------------------------------------------------

%%% ARROWS

%%% \to
\renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}}

%%% \mapsto
\let\oldmapsto\mapsto
\renewcommand*{\mapsto}{\mathchoice{\longmapsto}{\oldmapsto}{\oldmapsto}{\oldmapsto}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

%%%---------------------------------------------------------------------------

\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
\begingroup\hsize=12.5cm\advance\hsize\abstractmarginL\raggedright
\noindent
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}

\let\cdr@preabstracthook\@empty
\makeatother

%%%---------------------------------------------------------------------------

\title{High-order multistep coupling: convergence, stability and PDE application}
\alttitle{Couplage multipas d'ordre élevé\tralicstex{}{~}: convergence, stabilité et application aux EDPs}

\author{\firstname{Antoine} \middlename{E.} \lastname{Simon}\IsCorresp}
\address{ONERA, 6 Chemin de la Vauve aux Granges, 91123 Palaiseau, France}
\address{CMAP, CNRS --- École polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91120~Palaiseau~Cedex, France}
\email{antoine.simon@polytechnique.edu}

\author{\firstname{Laurent} \lastname{François}}
\address[1]{ONERA, 6 Chemin de la Vauve aux Granges, 91123 Palaiseau, France}
\email{laurent.francois@onera.fr}

\author{\firstname{Marc} \lastname{Massot}}
\address[2]{CMAP, CNRS --- École polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91120 Palaiseau Cedex, France}
\email{marc.massot@polytechnique.edu}

\keywords{\kwd{Multiphysics simulation} \kwd{high-order} \kwd{multistep coupling} \kwd{stability} \kwd{conjugate heat transfer}}
\altkeywords{\kwd{Simulation multi-physique} \kwd{ordre élevé} \kwd{couplage multipas} \kwd{stabilité} \kwd{transfert de chaleur conjugué}}

\thanks{The main support of ONERA and the Exa-MA project (Methods and Algorithms for Exascale, part of the PEPR NumPEx, Grant/Award Number: ANR-22-EXNU-0002 --- \url{https://numpex.org/}) is gratefully acknowledged. The support of Fondation de l'\'Ecole polytechnique through Initiative HPC@Maths is also gratefully acknowledged}
\CDRGrant[ONERA]{ANR-22-EXNU-0002}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\begin{abstract}
Designing coupling schemes for specialized advanced mono-physics solvers in order to conduct accurate and efficient multiphysics simulations is a key issue that has recently received a lot of attention. A novel high-order adaptive multistep coupling strategy has shown potential to improve the efficiency and accuracy of such simulations, but requires further analysis. The purpose of the present contribution is to conduct the numerical analysis of convergence of the explicit and implicit variants of the method and to provide a first analysis of its absolute stability. A simplified coupled problem is constructed to assess the stability of the method along the lines of the Dahlquist's test equation for ODEs. We propose a connection with the stability analysis of other methods such as splitting and ImEx schemes. A stability analysis on a representative conjugate heat transfer case is also presented. This work constitutes a first building block to an a priori analysis of the stability of coupled PDEs.
\end{abstract}

\begin{altabstract}
Concevoir des schémas de couplage de solveurs mono-physiques spécialisés pour mener des simulations multi-physiques précises et efficaces est un problème clé qui a été l'objet de beaucoup de travaux récemment. Un nouveau schéma de couplage multipas d'ordre élevé et à pas de temps adaptatif a montré un potentiel certain pour l'amélioration de l'efficacité et la précision de telles simulations, mais requiert toutefois une analyse plus approfondie. L'objectif de cette contribution est de mener une analyse numérique de la convergence des variantes explicite et implicite de cette méthode et d'introduire une première analyse de sa stabilité absolue. Un problème de couplage simplifié est construit pour évaluer la stabilité de la méthode, à l'instar de l'équation test de Dahlquist pour l'analyse des équations différentielles ordinaires. Nous comparons aussi la stabilité d'autres méthodes, comme le splitting d'opérateurs et des schémas ImEx. Une analyse de stabilité sur un cas test représentatif de transfert de chaleur conjugué est aussi présentée.
\end{altabstract}

\dateposted{2025-11-18}
\begin{document}

\maketitle

\section{Introduction}\label{sec:introduction}

The study of complex systems requires simulating the interactions between various physical phenomena involving multiple scales in time and space. The simulation of such multiphysics systems seldom relies on implementing a new solver for the coupled problem since each physics already possesses optimized and specialized solvers with an important and long-term legacy. As a consequence, multiphysics simulations are often built by coupling multiple existing solvers, which resolve separately each submodel and regularly exchange some quantities (volume or surface coupling terms) that we will refer to as \emph{coupling variables}. A wide range of applications resort to this approach, \eg\ fluid-structure interactions~\cite{farhat_two_2000}, fluid-plasma interactions~\cite{bourlet2024}, conjugate heat transfer~\cite{kazemi2013high,boulet2018modeling}, fluid-fluid coupling for acoustics~\cite{langenais2019accurate}.

So far, the key effort has been put on developing high-performance-computing (HPC) coupling libraries specialised in data transfer and spatial interpolation from one mesh to the other, \eg\ \textsc{Cwipi}~\cite{refloch2011cedre,cwipi_github}, and has been applied in an industrial context~\cite{Moureau_2024}. The time integration is however often basic, with the codes being coupled at a fixed rate (constant coupling time step), and with the exchanged quantities being frozen over each coupling step~\cite{keyes2013multiphysics}. This induces a first-order overall global error in the coupled dynamics. Besides, in most cases the coupling scheme is explicit since an implicit coupling is more difficult to implement and generally requires an iterative procedure. Consequently, for each simulation case, a fair amount of work is needed to determine a coupling time step that leads to a stable and precise coupled computation, which may be even more difficult if the dynamics and time scales strongly evolve with time. Thus, there is a need for a more efficient temporal coupling scheme, ensuring high-order accuracy in time, automatic adaptation of the coupling time step, and the ability to perform an implicit coupling. Solutions mimicking operator splitting~\cite{farhat_two_2000}, although easy to implement and fairly stable despite their explicit nature, cannot go beyond second-order accuracy. Implicit-Explicit (ImEx) schemes~\cite{kennedy_additive_2003} have been originally designed for the case where multiple operators, \eg\ convection and diffusion, act on a single physical system and spatial domain. They can however be advantageously applied as coupling schemes~\cite{kazemi2013high}, solving implicitly the uncoupled dynamics of each coupled model, and handling the coupling terms explicitly. Still, ImEx schemes require the use of the same implicit time integration scheme for all solvers and prevent subcycling. Moreover, the coupling remains explicit and thus prone to instability for strongly-coupled systems.

Multirate time schemes have a long standing history~\cite{gunther_multirate_2001} since their introduction in~\cite{gear_1980,gear_multirate_1984}. They rely on a separation of time scales between different terms and introduce a fixed time step ratio that prescribes the amount of subcycling performed for one term compared to the others. The approach has been used both for coupling~\cite{kang2023multirate} and for treating several coupled operators within the same domain simulation~\cite{wensch2009multirate,loffeld_2024} in mono-physics solvers. This approach is however more suited to mono-physics solvers as they require some important modifications to the time integration procedure. Moreover, the amount of subcycling cannot be changed, thus making it difficult to cope with configurations where two coupled models may have similar or evolving time scales, and/or have a varying coupling strength as time evolves.

An interesting alternative approach has been developed in the coupling library \textsc{Precice}~\cite{ruth2021quasi}, in which the coupling variables are interpolated within a coupling time step, based on their values at the substeps of the different solvers. This enables a high-order implicit coupling in time with arbitrary subcycling. However, an explicit variant relying on extrapolation from the last coupling step is not available, and the coupling time step cannot be dynamically adapted.

A new generic coupling strategy has been recently introduced to remedy these problems~\cite{francois_multistep_2023,Francois-2022,ali_2025}, where the coupling variables are approximated by high-order polynomials of time during a coupling time step, enabling each solver to have accurate exchange terms during their internal substepping. These polynomials are built by an interpolation procedure involving the values of these variables at the previous coupling times. The polynomial representation is similar to that used in some co-simulation algorithms~\cite{kubler2000two,busch2012}, which however are more focused on mechanical systems and assume a different structure for the coupling variables. The multistep coupling approach allows each solver of each subsystem to use any time integration scheme, potentially involving subcycling. The polynomial representation provides high orders of convergence and allows for the derivation of error estimates to dynamically adapt the coupling time step. Both explicit and implicit formulations of the coupling are available. Owing to its similarity with multistep schemes for ordinary differential equations (ODEs), this scheme is named \emph{multistep coupling},
even if our method does not fall into the class of classical linear multistep methods.

A natural question, while examining the literature, is to investigate if the multistep coupling can be thought of as a domain decomposition method in space-time, where the space decomposition is prescribed by the physical domains associated with each model. The time decomposition could consist of (potentially adaptive) discretisation of the overall time interval into multiple non-overlapping coupling time windows, solved in a sequential manner, each time with a high-order approximation of the coupling conditions (transmission conditions). However, the features of our scheme (non-monolithic method with black-box solvers for each subdomain involving different time integrators and substeps, different models and non-overlapping domains) imply that it does not fall into any of the traditional frameworks of domain decomposition, which are focused on the efficient parallelisation of monolithic-like implicit schemes or steady-state problems, such as domain overlap, optimised transmission conditions and coarse spaces (see~\cite{toselli2004domain,gander_landmarks} and references therein). Our method thus requires a dedicated numerical analysis investigation.

This multistep approach is promising, all the more if integrated within existing coupling libraries such as \textsc{Cwipi}. However, while for ODE integration schemes, be it one-step or multistep methods, the notion of absolute stability is well-defined through Dahlquist's test equation~\cite{hairer_book2} linked to an eigenvalue analysis of the Jacobian matrix, no such framework exists for coupling schemes. In general, a much more involved study is conducted, where stability properties are determined numerically through a large parametric study on simple test cases~\cite{busch2012}.

The purpose of this contribution is to conduct the numerical analysis of the multistep coupling scheme in terms of order and convergence, and to propose a framework for the stability analysis of coupling schemes, in the same spirit as Dahlquist's test equation for ODE solvers. It provides a first building block for a connection with coupled partial differential equations (PDEs), as encountered in practical applications, and paves the way for a complete a priori analysis of coupling stability.

The contribution is organized as follows. In a second section after the introduction, the high-order multistep coupling is presented and the consistency/order and convergence analysis conducted. In a third section, we present a simple coupled problem in order to numerically illustrate the convergence properties. It turns out that this simple problem is an interesting generalisation to coupled systems of the well-known Dahlquist's test equation. Its use for the construction of stability diagrams is discussed in a fourth section, were a comprehensive stability analysis is conducted. We also conduct a comparison between the multistep coupling scheme and some splitting~\cite{strang_construction_1968} and ImEx Additive Runge--Kutta (ARK) schemes~\cite{kennedy_additive_2003, kennedy_higher-order_2019} from the literature. The last section presents a numerical stability analysis of a conjugate heat transfer problem, and discussions on the link between the test equation and the conjugate heat model are presented in a conclusion.

\section{The multistep coupling scheme}\label{sec:multistep}

\subsection{Explicit and implicit formulations of the multistep coupling scheme}\label{sec:multistep_desc}

The multistep coupling scheme is a strategy to integrate a set of coupled ODEs, in particular those arising from the spatial semi-discretisation of coupled PDEs. For such systems, the global state vector $y \colon [0, T] \to \R^m$ with $T\in\R_+^*$ and $m\in\N^*$ can be formally decomposed into components as ${y = \parens{y_1^t \dotsm y_M^t}^t}$ ($M\in\N^*$). The component $y_i \colon [0, T] \to \R^{m_i}$ is the state vector of the $i$-th subsystem (associated with one particular solver in the multistep coupling scheme) with $m_i\in\N^*$. In the framework of coupled equations arising from a multiphysics problem, the dynamics of each subsystem $i=1,\dotsc,M$ depends both on the state vector $y_i$ and on the \emph{coupling variables} $u_i \colon [0, T] \to \R^{l_i}$ ($l_i\in\N^*$) which may depend on all subsystems $j=1,\dotsc,M$. We assume the form $u_i = h_i(y_1,\dotsc, y_M)$, with $h_i$ a Lipschitz function. In general, the coupling involves only a lower-dimensional interface between spatial domains so that in most cases $l_i \ll m_i$, \eg\ coupling a 3D fluid flow to a 3D wing model by pressure and displacement terms defined on the 2D wing surface. The whole coupled system's dynamics is described by:
\begin{equation}\label{eq:coupled_pb}
\begin{cases}
\dif_t y_i(t) = \F_i\parens[\big]{y_i(t), u_i(t), t} \\
u_i = h_i(y_1,\dotsc, y_M) \\
y_i \in \mathcal{C}^1\parens[\big]{[0, T], \R^{m_i}} \\
y_i(0) = y_{i,0} \in \R^{m_i}
\end{cases}
\qquad \text{for $i=1,\dotsc,M$}
\end{equation}
where, to ensure well-posedness of this equation, $\F_i \colon \R^m \times [0, T] \to \R^m$ are continuous functions, Lipschitz in terms of $y_i$ and $u_i$, and $h_i$ are Lipschitz functions. For each $i$, the problem described in equation~\eqref{eq:coupled_pb} is called a mono-physics problem, or subsystem.

The multistep coupling scheme consists in building, for the $n$-th coupling step $[t_n, t_{n+1}]$, polynomial predictions $\hat{u}_i^n \colon [t_n, t_{n+1}] \to \R^{l_i}$ of $u_i(t)$, $i=1,\dotsc,M$. These predictors are obtained, in the explicit case, by an extrapolation of the past states $u_{n-j}$, $j \in \lbrace 0,\dotsc,k \rbrace$, $k\in\N$. The prediction operators are denoted $\Psi_i$. In the explicit case, they are defined by:
\begin{equation}\label{eq:extrap_operator}
\Psi_i \colon
\left\lbrace
\begin{aligned}
	\parens{\R^{l_i}}^{k+1} \times [0, T]^{k+1} & \to \mathcal{C}^{\infty}\parens[\big]{[t_{n}, t_{n+1}], \R^{l_i}}
\\	\parens{u_i^{n-k},\dotsc, u_i^{n}, t_{n-k}, \dotsc, t_n} & \mapsto \hat{u}_i^n
\end{aligned}\right.
\end{equation}
where
\begin{equation}\label{eq:extrap_poly}
\hat{u}_i^n(t) = \Psi_i\parens{u_i^{n-k},\dotsc, u_i^{n}, t_{n-k}, \dotsc, t_n}(t) = \sum_{j=0}^{k} u_i^{n-j} \prod_{\substack{l=0 \\ l \neq j}}^k \dfrac{t-t_{n-l}}{t_{n-j}-t_{n-l}}.
\end{equation}
These predictors allow to integrate all subsystems independently over a coupling time step of size $\Dt_n$:
\begin{equation}\label{eq:coupled_pb_interp}
\dif_t y_i(t) = \F_i\parens*{y_i(t), \hat{u}_i^n(t), t}
\end{equation}
for all $i=1,\dotsc,M$, using $M$ mono-physics solvers. These solvers can use any reasonable time integration technique, potentially involving subcycling.

As mentioned in the introduction, the simplest and ubiquitous coupling scheme consists in setting $\hat{u}_i^n(t) = u_i^n$. This corresponds to a constant extrapolation ($k=0$) in equation~\eqref{eq:extrap_poly}. This strategy results in a first-order global error. We will show in Theorem~\ref{thm:convergence} that the multistep coupling technique has a global error of order $k+1$, with $k$ the polynomial degree of the predictions. Thus, it provides a relatively simple way to improve the time accuracy of a coupled simulation.

The definition of $\hat{u}_i$ in equation~\eqref{eq:extrap_poly} leads to an explicit scheme, since data is produced by extrapolating from the past. An implicit formulation can be obtained by using the future (and yet unknown) value $u_i^{n+1}$ as an additional node to define $\hat{u}_i^n$, so that $\Psi_i$ becomes an interpolant of~$u_i$:
\[
\hat{u}_i^n(t) = \Psi_i\parens{u_i^{n-k+1}, \dotsc, u_i^{n+1}, t_{n-k+1}, \dotsc, t_{n+1}}(t).
\]
Applying this to all subsystems, a nonlinear system on $u_{i}^{n+1}$, $i=1, \dotsc, M$, is obtained, which can be solved by Picard iterations or Newton-like schemes~\cite{degroote2013partitioned}. Some model-specific techniques stemming from the domain decomposition community~\cite{toselli2004domain} may also be applied to precondition the fixed-point operator, for instance by adopting optimized Robin-type coupling conditions if applicable. The solution method for the fixed-point problem in the implicit case is however beyond the scope of the present paper.

In the following, it will be useful to recast the multistep coupling scheme as a one-step method by writing it in vector form. For $n = 0, \dotsc, N-1$, the overall state vector at time $t_n$ is denoted $y^n = {\parens*{(y_1^n)^t,\dotsc, (y_M^n)^t}^t}$ and the concatenation of the last $k+1$ values of $y$ at times $t_n, \dotsc, t_{n-k}$ is:
\begin{equation}\label{eq:y_concat}
Y_n = {\parens[\big]{\parens{y^n}^t,\dotsc, \parens{y^{n-k}}^t}^t} \in E
\end{equation}
where $E = (\R^m)^{k+1}$. We denote as $\Phi_i$, $i=1,\dotsc,M$, the mono-physics integrators associated with each mono-physics problem, which we assume to be zero-stable as will be defined in Section~\ref{sec:zero_stability}. The integrator $\Phi_i$ corresponds to the application of the $i$-th solver for the time integration of equation~\eqref{eq:coupled_pb_interp} for subsystem $i$ over a coupling time step:
\begin{equation}
\Phi_i \colon
\left\lbrace
\begin{aligned}
	[0, T] \times \R^{m_i} \times \mathcal{C}^{0}\parens[\big]{[t_{n}, t_{n+1}], \R^{l_i}} \times]0, T] & \to \R^{m_i}
\\	(t, y_i^n, \hat{u}_i^n, \Dt_n) & \mapsto \Phi_i(t, y_i^n, \hat{u}_i^n, \Dt_n)
\end{aligned}\right.
\end{equation}
such that, for a time step $\Dt_n$, the multistep coupling procedure is:
\begin{equation}\label{eq:integration_operator}
y_i^{n+1} = y_i^n + \Dt_n\Phi_i\parens{t_n, y_i^n, \hat{u}_i^n, \Dt_n}.
\end{equation}
This last equation allows to rewrite the whole scheme as follows:
\begin{equation}\label{eq:scheme_expr_standard}
\begin{cases}
	Y_{n+1} = Y_n + \Dt_n\Phi\parens*{t_{n-k}, \dotsc, t_n, Y_n, \Dt_n}		& \text{(explicit)}
\\	Y_{n+1} = Y_n + \Dt_n\Phi\parens*{t_{n-k+1}, \dotsc, t_{n+1}, Y_{n+1}, \Dt_n}	& \text{(implicit)}
\end{cases}
\end{equation}
where $\Phi$ denotes the concatenation of all $\Phi_i$ operators composed with the $\Psi_i$ operators.
Note that the initial value $y_n$ is contained within $Y_{n+1}$, thus $\Phi$ need not depend on $Y_n$ for the implicit case. Assuming a solution to the implicit formulation exists (see Lemma~\ref{lem:imp_well_posed}), the implicit scheme can be written in the usual one-step form:
\begin{equation}\label{eq:implicit:onestep}
Y_{n+1} = Y_n + \Dt_n \tilde{\Phi}\parens{t_{n-k}, \dotsc, t_{n}, Y_{n}, \Dt_n}
\end{equation}
where $\tilde{\Phi}$ involves the corresponding solution operator.

Note that the multistep coupling scheme is not self-starting for $k>0$ in the explicit case or $k>1$ in the implicit case. Just as for multistep ODE integrators, multiple solutions can be considered. A first one is to perform the first $k$ coupling steps with another high-order coupling scheme. Another solution is to simply start at $k=0$ with a very small time step, and increase $k$ (and potentially the time step) at each subsequent step, until the targeted value of $k$ is reached. A final solution for the implicit scheme is to perform the first $k$ steps in a strongly coupled manner, \ie\ solving a more complex interpolation problem involving the first $k$ values of the coupling variables, used to define a single high-order interpolant for these first steps.

\subsection{Convergence of the multistep coupling scheme}\label{sec:multistep_proofs}

We now demonstrate that the multistep coupling scheme, either in explicit or implicit form, is zero-stable and consistent as long as each mono-physics integrator $\Phi_i$, $i=1,\dotsc,M$, satisfies the usual Lipschitz properties leading to their zero-stability. Convergence of the scheme can then be obtained. The result is valid for exact or approximate subsystem integration. In this section, we assume that the coupling time step is constant, set to $\Dt$. The overall simulation time is $T= N\Dt$, where $N \in \N^*$. Using a constant coupling time step allows to simplify the dependency of operators defined above, $\Psi_i$ then only depends on $\parens{u_i^{n-k},\dotsc, u_i^{n}, t_n, \Dt}$, and identically, $\Phi$ only depends on $\parens*{t_n, Y_n, \Dt}$ and $\tilde{\Phi}$ on $\parens*{t_n, Y_{n+1}, \Dt}$.

\subsubsection{Zero-stability}\label{sec:zero_stability}

A scheme is said to be \emph{zero-stable} if the amplification across multiple steps of perturbations is bounded according to the following definition.
\begin{defi}[Zero-stability]\label{def:zero_stab}
Let $\norm{\,\cdot\,}$ be any norm on $E$, and $\onestepoperator$ be an integration scheme expressed as $Y_{n+1} = Y_n + \Dt\onestepoperator(t_n, Y_n, \Dt)$. This scheme is said to be \emph{zero-stable} if there exists a constant $\cstezerostab > 0$ independent of $\Dt$ such that $\forall Y_0, Z_0 \in E$ and $\forall  {(\varepsilon_n)}_{n=0, \dotsc, N-1} \in E^N$, the sequences:
\[
\begin{cases}
	Y_{n+1} = Y_n + \Dt\onestepoperator(t_n, Y_n, \Dt)
\\	Z_{n+1} = Z_n + \Dt\onestepoperator(t_n, Z_n, \Dt) + \varepsilon_n
\end{cases}
\]
satisfy the following inequality:
\begin{equation}\label{eq:def_zero_stab}
\max_{0 \leqslant n \leqslant N} \norm{Y_n - Z_n} \leqslant \cstezerostab \bracks[\Bigg]{\norm{Y_0 - Z_0} + \sum_{n=0}^{N-1} \norm{\varepsilon_n}}.
\end{equation}
\end{defi}

Note that this definition is classical for one-step schemes~\cite[Section~3.2]{ascher1998computer} and may sometimes be referred to simply as \emph{stability}~\cite[equation~(16.1.7)]{schatzman2002numerical}. When applied on the one-step reformulation of a linear multistep ODE scheme, this definition is equivalent to the traditional root condition from the linear multistep literature~\cite[Section~(17.3)]{schatzman2002numerical}. Since our scheme does not belong to the category of linear multistep methods, but is a generalisation thereof, the more general Definition~\ref{def:zero_stab} should be used. This definition can also be used more directly to study the convergence of the method (see Theorem~\ref{thm:convergence}.).

Let us stress that zero-stability is a prerequisite for an integrator to be convergent. As already mentioned, all subsystem integrators $\Phi_i$ are assumed to be Lipschitz operators and thus correspond to zero-stable time integrators.

To establish the zero-stability of the multistep coupling scheme, we first need to show that the multistep coupling integrator $\Phi$ defined by equation~\eqref{eq:scheme_expr_standard} constitutes a Lipschitz operator.

\begin{lemm}[$\Phi$ is Lipschitz]\label{lem:lipsch}
Provided that the integrators $\Phi_i$ are Lipschitz operators, the multistep coupling operator $\Phi$ from equation~\eqref{eq:scheme_expr_standard} is continuous and Lipschitz in $Y$.
\end{lemm}

\begin{proof}
We first show that the extrapolation/interpolation operators $\Psi_i$ are all Lipschitz in terms of their $k+1$ first variables. Let us consider $u_i^{n-k},\dotsc, u_i^{n}, \tilde{u}_i^{n-k},\dotsc, \tilde{u}_i^{n} \in \R^{l_i}$ and denote $u_i = \parens[\big]{(u_i^{n-k})^t,\dotsc, (u_i^{n})^t}^t$ and $\tilde{u}_i = \parens[\big]{(\tilde{u}_i^{n-k})^t,\dotsc, (\tilde{u}_i^{n})^t}^t$. Consider the infinity norm $\norm{\,\cdot\,}_\infty$ on the space $\mathcal{C}^{0}\parens[\big]{[t_{n}, t_{n+1}], \R^{l_i}}$, the infinity norm $\norm{\,\cdot\,}_{\infty, \I}$ on the space $\mathcal{C}^{0}(\I, \R)$ where $\I$ is a closed interval of $\R$, and the infinity norm $\abs{\,\cdot\,}_\infty$ on $(\R^{l_i})^{k+1}$. Using the dimensionless time $\tau=\frac{t-t_n}{\Dt}$, we denote $\L_j, j=0,\dotsc,k$ the basis of Lagrange interpolation polynomials at points $\tau_0=0$, $\tau_1=-1$, \dots, $\tau_k=-k$. The explicit predictors are expressed as follows:
\[
\hat{u}_i(t) = \Psi_i(u_i^{n},\dotsc, u_i^{n-k}, t_n, \Dt)(t) = \sum_{j=0}^{k}u_i^{n-j}\L_j\parens*{\tau}
\]
with a similar expression for $\hat{\widetilde{u}}_i(t)$. The following bound can be obtained directly:
\[
\norm{\hat{u}_i - \hat{\widetilde{u}}_i}_\infty \leqslant \abs{u_i - \tilde{u}_i}_\infty \cdot \sum_{j=0}^{k} \norm{\L_j}_{\infty, [0, 1]}.
\]
Hence, all the $\Psi_i$ are Lipschitz in terms of their $k+1$ first variables. Besides, the Lipschitz constant $\sum_{j=0}^{k} \norm{\L_j}_{\infty, [0, 1]}$ {does not depend on $\Dt$.}

A similar result can be obtained for the case of interpolation, using $\tilde{\L}_j$, $j=0,\dotsc,k$ the basis of Lagrange interpolation polynomials at points $\tau_0=1$, $\tau_1=0$, \dots, $\tau_k=-k+1$.

Since all $h_i$ and all integrators $\Phi_i$ are Lipschitz, then, $\Phi$ is also Lipschitz in terms of $Y_n$ by composition of Lipschitz functions.
\end{proof}

The zero-stability of the explicit multistep coupling scheme can now be stated in the following proposition.

\begin{prop}[Zero-stability of the explicit multistep coupling scheme]\label{thm:zero_stab_exp}
The explicit multistep coupling scheme given by equation~\eqref{eq:scheme_expr_standard} is zero-stable.
\end{prop}

\begin{proof}
As already stated, the subsystem integrators $\Phi_i$ are Lipschitz with a Lipschitz constant independent of $\Dt$ for all $0 < \Dt \leqslant T$. Lemma~\ref{lem:lipsch} ensures that there exists a Lipschitz constant $L_\Phi > 0$ of the overall integrator $\Phi$ independent of $\Dt$. Consider $Y_0, Z_0 \in E$, ${(\varepsilon_n)}_{n=0, \dots, N-1} \in E^N$ and the sequences:
\[
\begin{cases}
	Y_{n+1} = Y_n + \Dt\Phi(t_n, Y_n, \Dt)
\\	Z_{n+1} = Z_n + \Dt\Phi(t_n, Z_n, \Dt) + \varepsilon_n.
\end{cases}
\]
Then, for any $n \leqslant N-1$:
\[
\begin{split}
\abs{Y_{n+1} - Z_{n+1}}_\infty
	& \leqslant \abs[\Big]{Y_{n} - Z_{n} + \Dt \bracks[\big]{\Phi(t_n, Y_n, \Dt) - \Phi(t_n, Z_n, \Dt)} - \varepsilon_n}_\infty
\\	& \leqslant (1 + L_\Phi\Dt) \abs{Y_{n} - Z_{n}}_\infty + \abs{\varepsilon_n}_\infty
\\	& \leqslant (1 + L_\Phi\Dt)^{N} \abs{Y_{0} - Z_{0}}_\infty + \sum_{i=0}^{N-1} (1 + L_\Phi\Dt)^{N}\abs{\varepsilon_{i}}_\infty
\\	& \leqslant \underbrace{e^{L_\Phi T}}_{\cstezerostab}\abs{Y_{0} - Z_{0}}_\infty + \underbrace{e^{L_\Phi T}}_{\cstezerostab}\sum_{i=0}^{N-1} \abs{\varepsilon_{i}}_\infty. \qedhere
\end{split}
\]
\end{proof}

To prove the same property for the implicit scheme, we first need to show that the implicit system in equation~\eqref{eq:scheme_expr_standard} is well-posed.

\begin{lemm}[Well-posedness of the implicit multistep coupling scheme]\label{lem:imp_well_posed}
There exists $\Dt^* > 0$ such that, for all $0 < \Dt \leq \Dt^*$ and for all $Y_n\in E$, there exists a unique solution $Y_{n+1}$ to the non-linear equation $Y_{n+1} = Y_n + \Dt \Phi(t_n, Y_{n+1}, \Dt)$.
\end{lemm}

\begin{proof}
Let $Y_n$ be an element of $E$ and $t_n\in[0, T\mathclose{[}$. Lemma~\ref{lem:lipsch} ensures that there exists $L_{\phiimpdemo} > 0$ such that for all $\Dt > 0$ and for all $Y, Z \in E$: $\norm[\big]{\phiimpdemo(t_n, Y, \Dt) - \phiimpdemo(t_n, Z, \Dt)} \leq L_{\phiimpdemo} \norm*{Y - Z}$.

Then, for all $\Dt > 0$ and for all $Y, Z \in E$:
\[
\norm[\Big]{\parens[\big]{Y_n + \Dt\phiimpdemo(t_n, Y, \Dt)} - \parens[\big]{Y_n + \Dt\phiimpdemo(t_n, Z, \Dt)}} \leq \Dt L_{\phiimpdemo} \norm{Y - Z}.
\]
Let $\Dt^*$ be a positive real number such that $\Dt^* L_{\phiimpdemo} < 1$, $\Dt^*$ does not depend on $Y_n$. For all $0 < \Dt \leq \Dt^*$, the function $Y \mapsto Y_n + \Dt\phiimpdemo(t_n, Y, \Dt)$ is a contraction, therefore, it admits exactly one fixed point that we can denote $Y_{n+1}$.
\end{proof}

This proof also ensures that the previously introduced implicit integration operator $\tilde{\Phi}$ is well defined. Note that the value of $\Dt^*$ and the constant $K$ are related to the Lipschitz constant~$L_{\phiimpdemo}$, which is directly affected by the evaluation of the Lipschitz constant of the one-step ODE integrators $\Phi_i$. For $\Dt^*$ and $K$ to be close to optimality, the proper framework is that of the one-sided Lipschitz condition~\cite[Section~IV.14]{hairer_book2}, which is much more favourable in particular for contracting dynamics, the category into which dissipative systems fall.

The zero-stability of the implicit multistep coupling scheme can now be established.

\begin{prop}[Zero-stability of the implicit multistep coupling schemes]\label{thm:zero_stab_imp}
The implicit multistep coupling scheme given by equation~\eqref{eq:scheme_expr_standard} is zero-stable for time steps ensuring its well-posedness.
\end{prop}

\begin{proof}
We consider $Y_0, Z_0 \in E$, ${(\varepsilon_n)}_{n=0, \dotsc N-1} \in E^N$ and the sequences defined below:
\[
\begin{cases}
	Y_{n+1} = Y_n + \Dt\Phi(t_n, Y_{n+1}, \Dt)
\\	Z_{n+1} = Z_n + \Dt\Phi(t_n, Z_{n+1}, \Dt) + \varepsilon_n
\end{cases}
\]
which are well-defined for all $0 < \Dt < \Dt^*$ according to Lemma~\ref{lem:imp_well_posed}. First, for any $0 < \Dt < \Dt^*$:
\[
\frac{1}{1 - L_\Phi\Dt} = 1 + \frac{\Dt L_\Phi}{1 - L_\Phi\Dt} \leq 1 + \Dt K^*
\]
where $K^* \equiv \frac{L_\Phi}{1 - L_\Phi\Dt^*}$. Then, for any $n \leqslant N-1$:
\[
\begin{split}
\abs{Y_{n+1} - Z_{n+1}}_\infty
	& \leqslant \abs[\Big]{Y_{n} - Z_{n} + \Dt \bracks[\big]{\Phi(t_n, Y_{n+1}, \Dt) - \Phi(t_n, Z_{n+1}, \Dt)} - \varepsilon_n}_\infty
\\	& \leqslant \frac{1}{1 - L_\Phi\Dt} \abs{Y_{n} - Z_{n}}_\infty + \frac{1}{1 - L_\Phi\Dt} \abs{\varepsilon_n}_\infty 
\\	& \leqslant \parens[\Bigg]{\frac{1}{1 - L_\Phi\Dt}}^{N} \abs{Y_{0} - Z_{0}}_\infty + \sum_{i=0}^{N-1} \parens[\Bigg]{\frac{1}{1 - L_\Phi\Dt}}^{N}\abs{\varepsilon_{i}}_\infty
\\	& \leqslant \parens{1 + \Dt K^*}^{N} \abs{Y_{0} - Z_{0}}_\infty + \sum_{i=0}^{N-1} \parens{1 + \Dt K^*}^{N}\abs{\varepsilon_{i}}_\infty
\\	& \leqslant \underbrace{e^{K^* T}}_{\cstezerostab}\abs{Y_{0} - Z_{0}}_\infty + \underbrace{e^{K^* T}}_{\cstezerostab}\sum_{i=0}^{N-1} \abs{\varepsilon_{i}}_\infty. \qedhere
\end{split}
\]
\end{proof}

\subsubsection{Consistency and convergence}

The accuracy of the solution computed with the multistep coupling scheme can now be assessed by considering its local and global errors, defined below. In this section, the consistency of the scheme is proved, which, combined with the previous results, establishes its convergence.

\begin{defi}[Global and local error --- convergence]\label{def:error}
Let us consider the integration scheme equation~\eqref{eq:scheme_expr_standard}, $y$ a solution of equation~\eqref{eq:coupled_pb} and $Y(t) = \parens[\big]{y(t)^t, y(t-\Dt)^t,\dotsc, y(t-k\Dt)^t}^t$. Let $\norm{\,\cdot\,}$ be any norm on $E$. The \emph{global error} of the scheme is:
\[
e_n = Y(t_n) - Y_n
\]
and the scheme is \emph{convergent} if:
\[
\lim_{\Dt \to 0} \; \max_{n=0, \dotsc, N} \norm{e_n} = 0.
\]
We call \emph{local error} or \emph{local truncation error} the quantity:
\[
\mathcal{E}_n = Y(t_{n+1}) - Y(t_n) - \Dt\Phi\parens[\big]{t_n, Y(t_n), \Dt}.
\]
\end{defi}

\begin{defi}[Consistency]
With the same hypothesis as in Definition~\ref{def:error} the scheme is said to be \emph{consistent} if the local truncation error verifies:
\[
\abs{\mathcal{E}_n} \leqslant \Dt \epsilon(\Dt)
\]
with $\epsilon$ a function such that $\epsilon(\Dt) \to 0$ when $\Dt \to 0$.

Besides, if $\epsilon(\Dt) = \O(\Dt^p)$, the scheme is said to be consistent at order $p$.
\end{defi}

\begin{prop}[Consistency of the multistep coupling scheme]\label{thm:consistency}
Assuming that the subsystem integrators\/ $\Phi_i$ are consistent at order $p_i \geq k+1$, the multistep coupling scheme given in equation~\eqref{eq:scheme_expr_standard} is consistent and the order of consistency is $k+1$ where $k$ is the degree of the polynomial predictor.
\end{prop}

The consistency has been proven to be of order $k+1$ prediction polynomials of degree $k$, assuming exact integration of the subsystems~\cite{francois_multistep_2023}. In the present contribution, we propose a more general proof considering non-exact mono-physics integrators $\Phi_i$ instead.

\begin{proof}
All $\F_i(y_i, u_i, t)$ are Lipschitz both in terms of $y_i$ and $u_i$. Let us denote the corresponding Lipschitz constants $L_{i, y}$ and $L_{i, u}$. Consider $\widetilde{\widetilde{y}}_i^n \in \R^{m_i}$ the solution at time $t_n$. One step of the multistep scheme yields for the $i$-th subsystem:
\[
\widetilde{\widetilde{y}}_i^{n+1} = \widetilde{\widetilde{y}}_i^n + \Dt \Phi_i(t_n, \widetilde{\widetilde{y}}_i^n, \hat{u}_i, \Dt).
\]
An exact integration of equation~\eqref{eq:coupled_pb_interp} with the same initial condition gives a solution $\tilde{y}_i$:
\[
\tilde{y}_i(t) = \widetilde{\widetilde{y}}_i^n + \int_{t_n}^{t} \F_i\parens[\big]{\tilde{y}_i(s), \hat{u}_i(s), s} \dif s
\]
where $\hat{u}_i(s)$ is the polynomial prediction of the coupling variables. Finally, an exact integration of equation~\eqref{eq:coupled_pb} with the same initial condition as above gives a solution $y_i$ of the original problem such that:
\[
y_i(t) = \widetilde{\widetilde{y}}_i^n + \int_{t_n}^{t} \F_i\parens[\big]{y_i(s), u_i(s), s} \dif s
\]
where $u_i(s)$ is the exact evolution of the coupling variables.

Let us notice that $y_i(t_n) = \tilde{y}_i(t_n) = \widetilde{\widetilde{y}}_i^n$.

The absolute local truncation error $\abs{\mathcal{E}_n} = \abs[\big]{\widetilde{\widetilde{y}}_i^{n+1} - y_i(t_{n+1})}$ then can be written:
\[
\begin{split}
\abs[\Big]{y_i(t_{n+1}) - y_i(t_{n}) - \Dt \Phi_i\parens[\big]{t_n, y_i(t_n), \hat{u}_i, \Dt}}
	& = \abs[\big]{y_i(t_{n+1}) - \widetilde{\widetilde{y}}_i^{n+1}}
\\	& = \abs[\big]{y_i(t_{n+1}) - \tilde{y}_i(t_{n+1}) + \tilde{y}_i(t_{n+1}) - \widetilde{\widetilde{y}}_i^{n+1}}
\\	& \leq \underbrace{\abs[\big]{y_i(t_{n+1}) - \tilde{y}_i(t_{n+1})}}_{\termeun} + \underbrace{\abs[\big]{\tilde{y}_i(t_{n+1}) - \widetilde{\widetilde{y}}_i^{n+1}}}_{\termedeux}
\end{split}
\]
where $\termeun$ represents the error arising from the polynomial approximation of the coupling variables, and $\termedeux$ represents the error produced by the non-exact subsystem integration.

First, let us consider $\termeun$:
\[
\begin{split}
\termeun
	& = \abs[\big]{y_i(t_{n+1}) - \tilde{y}_i(t_{n+1})}
\\	& = \abs*{\widetilde{\widetilde{y}}_i^{n} + \int_{t_n}^{t_{n+1}} \F_i\parens[\big]{y_i(s), u_i(s), s} \dif s - \widetilde{\widetilde{y}}_i^{n} - \int_{t_n}^{t_{n+1}} \F_i\parens[\big]{\tilde{y}_i(s), \hat{u}_i(s), s} \dif s}
\\	& \leq \int_{t_n}^{t_{n+1}} \abs[\Big]{\F_i\parens[\big]{y_i(s), u_i(s), s} - \F_i\parens[\big]{\tilde{y}_i(s), \hat{u}_i(s), s}} \dif s
\\	& \leq L_{i, u}\Dt \max_{[t_n, t_{n+1}]}\abs{u_i - \hat{u}_i} + \int_{t_n}^{t_{n+1}}L_{i, y}\abs[\big]{\tilde{y}_i(s) - y_i(s)} \dif s.
\end{split}
\]
Since $t \mapsto L_{i, u}\Dt \max_{[t_n, t]}\abs{u_i - \hat{u}_i}$ is a non-decreasing function of $t$, Grönwall's inequality ensures that:
\[
\termeun \leq L_{i, u}\Dt \max_{[t_n, t_{n+1}]}\abs{u_i - \hat{u}_i} \exp(L_{i, y}\Dt).
\]
Besides, the interpolation approximation error $\max_{[t_n, t]}\abs{u_i - \hat{u}_i}$ is of order $\O(\Dt^{k+1})$ for an interpolation polynomial of degree $k$, hence:
\[
\termeun = L_{i, u} \Dt e^{L_{i, y}\Dt} \O(\Dt^{k+1}) = \O(\Dt^{k+2}).
\]
Then, consider $\termedeux$:
\[
\termedeux = \abs[\big]{\tilde{y}_i(t_{n+1}) - \widetilde{\widetilde{y}}_i^{n+1}} = \abs*{\widetilde{\widetilde{y}}_i^n + \int_{t_n}^{t_{n+1}} \F_i\parens[\big]{\tilde{y}_i(s), \hat{u}_i(s), s} \dif s - \widetilde{\widetilde{y}}_i^{n} - \Dt \Phi_i(t_n, \widetilde{\widetilde{y}}_i^n, \hat{u}_i, \Dt)}.
\]
Since the integrators $\Phi_i$ are consistent at order $p_i \geq k+1$, $\termedeux$ satisfies:
\[
\termedeux = \O(\Dt^{p_i+1}) = \O(\Dt^{k+2}).
\]
Finally, $\abs[\big]{y_i(t_{n+1}) - y_i(t_{n}) - \Dt \Phi_i\parens[\big]{t_n, y_i(t_n), \hat{u}_i, \Dt}} = \O(\Dt^{k+2})$ and the scheme is consistent at order $k+1$.
\end{proof}

Thus, the order of consistency is the degree $k$ of the prediction polynomials, and the subsystem integrators should be at least of order $k$ to preserve the order of the coupling scheme.
A formula to compute the error constant is derived in Appendix~\ref{appendix:error}. Much like multistep ODE integrators, increasing the order does not necessarily induce an increase in computational cost, but only an increase in storage for keeping past values of $u_i$. This is a strong advantage of the multistep coupling method.

We may now prove that the multistep coupling scheme satisfies Definition~\ref{def:error} and therefore is convergent.

\begin{theo}[Convergence of the multistep coupling scheme]\label{thm:convergence}
The multistep coupling scheme is convergent and converges at order $k+1$ for any $k \geq 0$.
\end{theo}

\begin{proof}
When plugging into Definition~\ref{def:zero_stab} of zero-stability $Y_0 = Z_0$ and $\varepsilon_n = \mathcal{E}_n$ (local truncation error from Definition~\ref{def:error}), it can be seen that the sequence $Z_n$ is the sequence $Z_n = Y(t_n)$ of the values of the exact solution of equation~\eqref{eq:coupled_pb}, while $Y_n$ is the sequence of values given by the multistep coupling scheme. Then the results of Theorems~\ref{thm:zero_stab_exp}, \ref{thm:zero_stab_imp} and~\ref{thm:consistency} give:
\[
\abs[\big]{Y_n - Y(t_n)} \leq \cstezerostab \underbrace{\abs[\big]{Y_0 - Y(0)}}_{= 0} + \cstezerostab \underbrace{\sum_{i=0}^{N-1}\underbrace{\abs{\mathcal{E}_i}}_{= \mathcal{O}(\Dt^{k+2})}}_{= \mathcal{O}(\Dt^{k+1})} = \mathcal{O}(\Dt^{k+1}).
\]
Therefore, per Definition~\ref{def:error}, we can conclude the proof.
\end{proof}

\section{A simple model for the study of coupled equations}\label{sec:model_equation}

Now that we have demonstrated the previous essential theoretical properties of the scheme, we search for a simple coupled model to test the method, in particular to investigate its stability. In the literature, couplings are generally referred to as \emph{weak} when stability is not a particular issue~\cite{kazemi2013high}. On the opposite, couplings inducing strong stability limits on explicit schemes are referred to as \emph{strong} couplings~\cite{kazemi_kamyab_analysis_2014, birken2015fast}. However, no generic criterion has been proposed to predict how strong a coupling is.

In the context of co-simulation, the stability of coupling strategies is analyzed using chains of mass-spring systems~\cite{busch2012}. Such test cases are already too complex, so that their stability can only be studied by a systematic numerical scan of the parameter space to locate stable set of parameter values and time step. This makes the analysis cumbersome and time-consuming, while also restricting the thorough understanding of coupling stability that can be obtained from~it.

Having a complete theoretical picture of the coupling stability for generic coupled systems seems however out of reach. In the framework of ODEs, the Dahlquist's test equation $\dif_t y = \lambda y$ is a well-established choice for the study of absolute stability and is built as follows. Any ODE $\dif_t y = f(t, y)$ can be linearised locally along the solution path to obtain $\dif_t y = \parens{\partial_y f} y$. Diagonalisation of the Jacobian of $f$ yields a system of scalar equations $\dif_t x = \mu_i x$, where $\{\mu_i\}_i$ is the spectrum of the Jacobian. Stability is then ensured if the spectrum multiplied by the time step is contained within the stability domain of the chosen ODE integrator. The use of Dahlquist's test equation relies on a diagonalisation procedure of the linearised ODE, which is not relevant for coupled problems, since a global diagonalisation of the problem would destroy the partitioned nature of the approach. Hence a different test model should be sought.

In this section, we introduce a simple equation that aims at reproducing a coupled problem in its simplest form. Convergence results on this equation are presented to verify the convergence orders from Theorem~\ref{thm:convergence}. With certain parameter choice, stability seems to be an issue with the multistep explicit scheme, hence a deeper analysis of the scheme stability is proposed.

\subsection{The \texorpdfstring{$2\times2$}{2 x 2} test equation}\label{sec:gen_dahlquist_and_2x2}

The simplest possible coupled problem is a coupling between two scalar equations as given in equation~\eqref{eq:two_scalar_coupled}:
\begin{equation}\label{eq:two_scalar_coupled}
\begin{cases}
	\dif_t y_1 = a_1y_1 + a_{12}y_2
\\	\dif_t y_2 = a_{21}y_1 + a_2y_2
\end{cases}
\quad \Longleftrightarrow \quad
\dif_t \underbrace{
\begin{pmatrix}
	y_1
\\	y_2
\end{pmatrix}}_{\equiv y} = \underbrace{
\begin{pmatrix}
	a_{1}	& a_{12}
\\	a_{21} 	& a_{2}
\end{pmatrix}}_{\equiv \matdeuxdeux}
\begin{pmatrix}
	y_1
\\	y_2
\end{pmatrix}
\end{equation}
where $a_i\in\mathbb {C}$. We call this system the \emph{\testeqdeuxdeux{}}. It corresponds to equation~\eqref{eq:coupled_pb} with $M = 2$, $u_1 = h_1(y_1, y_2) = y_2$, $u_2 = h_2(y_1, y_2) = y_1$, $\F_1 \colon (y_1, u_1) \mapsto a_1y_1 + a_{12} u_1$ and $\F_2 \colon (y_2, u_2) \mapsto a_{21} u_2 + a_2y_2$. Variants of this model will be presented in Section~\ref{sec:stability_analysis}.

When the multistep coupling scheme is applied to this model, the two scalar ODEs $\dif_t y_1 = a_1 y_1 + a_{12} \hat{u}_1$ and $\dif_t y_2 = a_{21} \hat{u}_2 + a_2 y_2$ are integrated in parallel, which is sometimes referred to as a Jacobi coupling~\cite{busch2012}. At each coupling time step both predictors $\hat{u}_1$ and $\hat{u}_2$ are updated.

\subsection{Numerical convergence results}\label{sec:conv_num}

We now use the previous model to check that the convergence rates established in Theorem~\ref{thm:convergence} are obtained in practice. Both the implicit and explicit formulations are applied to the simple \testeqdeuxdeux{}, with parameters $a_{12}=40$, $a_{21}=40$, $a_2=-80$ and $y(0) = (5, -2)^t$. Multiple fixed-time-step simulations are performed over the time interval $[0,1]$ with different values of the time step. As discussed in Section~\ref{sec:multistep_desc}, the scheme is not self-starting. Therefore, we initialise the prediction polynomials by using exact values of the coupling variables in negative time ($t_{-j} < 0$). Even though the system is dissipative, this does not produce any numerical artifacts for this simple model. The error is computed as the $L^2$-norm in time of the difference between these solutions and the analytical solution $\exp(At) y(0)$. The obtained convergence graphs are presented in Figure~\ref{fig:conv_multistep}, where the numerically computed orders are indicated.

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{conv_graph_explicit_a1=-80.0_a12=40.0_a21=40.0_a2=-80.0_notitle.png}
\caption{Convergence of the explicit scheme}\label{fig:conv_multistep:err_exp}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{conv_graph_implicit_a1=-80.0_a12=40.0_a21=40.0_a2=-80.0_notitle.png}
\caption{Convergence of the implicit scheme}\label{fig:conv_multistep:err_imp}
\end{subfigure}
\caption{Convergence of the multistep scheme applied to the \testeqdeuxdeux{} with $a_1=-80$, $a_{12}=40$, $a_{21}=40$, $a_2=-80$ and $y(0) = (5, -2)^t$ integrated over the time interval $[0, 1]$.}\label{fig:conv_multistep}
\end{figure}

All schemes reach their design order of convergence. The higher-order formulations provide an important decrease of the error, as expected. It may be noted that, for a given order $k>1$, the implicit formulation possesses a lower error constant and is more precise.
This is coherent with the error constants derived theoretically in Appendix~\ref{appendix:error}.

It should however be observed that the explicit multistep scheme, especially at higher orders, has an error which diverges when the time step is too large. It is actually related to a blow-up of the coupled solution. This instability phenomenon is investigated in the following section.

\section{Stability analysis}\label{sec:model_stability}

It was observed in Section~\ref{sec:conv_num} that the multistep coupling scheme may suffer from instabilities. It would be interesting to be able to predict them from an a priori spectral analysis of the linearised coupled system, as done for ODEs with Dahlquist's test equation.
In this section, we first show that the \testeqdeuxdeux{}~\eqref{eq:two_scalar_coupled} is a relevant case for the study of coupled ODEs.
We then recast the multistep scheme as a recurrence formula involving a matrix whose properties dictate the stability of the coupling, and we draw stability diagrams, which are compared with those from common splitting schemes and state-of-the-art ImEx ARK methods.

\subsection{On the characterisation of the stability domain}\label{sec:stability_analysis}

One can notice that, assuming complex values for $\matdeuxdeux$, the space of model parameters to analyze is of dimension~$8$. This makes it difficult to characterize intuitively and efficiently the complete stability domain of a coupling scheme. Therefore, two simplifying approaches are proposed to reduce the number of dimensions.

\subsubsection{Real asymmetrical case}\label{sec:asymmetrical_model}

A first approach is to consider a general asymmetrical matrix $\matdeuxdeux$ with real coefficients to reduce the number of dimensions to 4. Then, cross-sections of the parameter space can be made with a specific parametrization of $\matdeuxdeux$, already used for the analysis of multirate schemes~\cite{kvaerno_stability_2000}. We introduce the following parameters:
\begin{equation}
\delta = \dfrac{a_2}{a_1}
\qquad \text{and} \qquad
\gamma = \dfrac{a_{12}a_{21}}{a_1a_2}.
\end{equation}
The parameter $\delta$ represents the ratio of the internal time scales of both scalar subsystems, while $\gamma$ can be interpreted as a coupling strength, which measures how much the eigenvalues of the coupled problem depart from those of the scalar subsystems (\ie\ from the diagonal of $\matdeuxdeux$). These settings are relevant for several reasons. First, the stability of the continuous problem is guaranteed if and only if $a_1 + a_2 < 0$ and $\gamma < 1$. Besides, the stability of the multistep coupling scheme for a given $\Dt$ only depends on $(a_1\Dt, \delta, \gamma)$.\footnote{This was confirmed for orders up to 15 using symbolic computation of the coefficients of the characteristic polynomial of the amplification matrix defined later in Section~\ref{sec:rec_mac_multi}.} Thus, the use of $\gamma$ reduces by 1 the number of dimensions, and two-dimensional cuts can be performed by specifying the value of $\delta$. This approach is analysed in Section~\ref{sec:kvaerno_stability}.

\subsubsection{Complex symmetrical case}\label{sec:symmetrical_model}

A second approach is to assume a symmetrical structure for $\matdeuxdeux$ with equal diagonal terms. Considering complex coefficients, these assumptions on $\matdeuxdeux$ reduce to 4 the number of dimensions for the parameter space. A cut can be performed (\eg\ fixing real and imaginary parts of the diagonal term) to reduce it to a two-dimensional space in the spirit of ImEx schemes and thus allows to characterize it and use it more easily. This approach is developed in Section~\ref{sec:symmetrical_model_eq_stability}. The symmetrical coupled system can be linked to the generalised Dahlquist's test equation used in the study of ImEx schemes:
\begin{equation}\label{eq:dahlquist_gen}
\dif_t x = (a_1 + a_{12}) x.
\end{equation}
Indeed, assuming initial conditions $y_1(0) = y_2(0) = x(0)$, the solutions to both equations verify $y_1(t) = y_2(t) = x(t)$, $\forall t \geq 0$. The generalised Dahlquist's test equation actually corresponds to the solution of equation~\eqref{eq:two_scalar_coupled} along the eigenvector $(1, 1)^t$ of $\matdeuxdeux$ which is associated with the eigenvalue $a_1 + a_{12}$. Thus, stability domains obtained on both equations are similar, except near the stability limit of the continuous system, since the second eigenvector $(1,-1)^t$ of $\matdeuxdeux$ (with eigenvalue $a_1 - a_{12}$) is not taken into account by the generalised Dahlquist's test equation.

Both real asymmetrical and complex symmetrical cases are relevant for our study of stability. On the one hand, the asymmetrical case allows to study the coupling of subsystems possessing different internal time scales, and to see the effect of the off-diagonal coupling terms through a single parameter $\gamma$. On the other hand, the complex symmetrical case, limited to subsystems with identical time scales, allows to study the effect of imaginary parts in the coefficients. Combining both approaches will allow for a more complete picture of stability to be obtained.

\subsection{Recurrence matrix for the multistep coupling scheme}\label{sec:rec_mac_multi}

In this section, we recast the multistep coupling method in one-step vectorial form. For a more generic formulation we temporarily consider the coupling of two ODE systems, of size $m_1$ and~$m_2$. The case of the previous \testeqdeuxdeux{} can be easily recovered from the results presented in this section.

Let $m = m_1+m_2$ and let $A_1 \in \C^{m_1\times m_1}$, $A_2 \in \C^{m_2\times m_2}$, $A_{12} \in \C^{m_1\times m_2}$ and $A_{21} \in \C^{m_2\times m_1}$ be complex-valued matrices. We consider the following linear coupled system:
\begin{equation}\label{eq:lin_coupled_lin}
\dif_t y = Ay, \qquad A \equiv
\begin{pmatrix} A_1 & A_{12} \\ A_{21} & A_2
\end{pmatrix}
\end{equation}
where the first $m_1$ components of $y$ constitute the state vector $y_1$ of the first subsystem, and the $m_2$ remaining ones constitute those of the state vector $y_2$ from the second subsystem. Let $\Dt > 0$ be the time step used for our coupling scheme and $Z_i \equiv A_i \Dt$. Then, for any degree $k \geq 0$, there exists a matrix $R(Z_1, Z_{12}, Z_{21}, Z_2, k)$ called the \emph{recurrence matrix} or \emph{amplification matrix} such that one step of the multistep coupling scheme can be expressed as:
\begin{equation}
Y_{n+1} = RY_n.
\end{equation}
The coupling scheme is stable if the spectral radius $\rho(R)$ is lower than~1. In the case of the \testeqdeuxdeux{}, the same expressions of amplification matrices hold with $z_i = a_i \Dt$.

The stability analysis will be conducted assuming an exact integration of the subsystems, so as to avoid making particular assumptions on the methods used for the subsystems integration. Although not presented here, numerical and theoretical stability analyses with non-exact integration have shown that using a non-subcycled low-order implicit scheme for the subsystem integration may increase the stability of the coupling scheme at the expense of accuracy.

Introducing $\tau = \frac{t-t_n}{\Dt}$, the exact integration of equation~\eqref{eq:coupled_pb_interp} for the first subsystem yields:
\begin{equation}
y_1^{n+1} = e^{Z_1}y_1^n + \int_{0}^{1} e^{Z_1(1-\tau)} Z_{12}\hat{u}_1^n(t_n + \tau \Dt)\dif \tau.
\end{equation}
For polynomial predictions of degree $k$, the amplification matrix of the explicit multistep scheme is:
\begin{equation}\label{eq:exp_rec_mat}
R_{\mathrm{MS},\mathrm{exp}} =
\begin{pmatrix}
e^{Z_1} & S^0_1 & 0 & S^1_1 & \cdots & 0 & S^k_1 \\
S^0_2 & e^{Z_2} & S^1_2 & 0 & \cdots & S^k_1 & 0 \\
I_{m_1} & 0 & \cdots & \cdots & \cdots & 0 & 0 \\
0 & I_{m_2} & \ddots & & & \vdots & \vdots \\
\vdots & \ddots & \ddots & \ddots & & \vdots & \vdots \\
\vdots & & \ddots & I_{m_1} & \ddots & \vdots & \vdots \\
0 & \cdots & \cdots & 0 & I_{m_2} & 0 & 0
\end{pmatrix} \in \mathbf{C}^{m(k+1)\times m(k+1)}.
\end{equation}
The basis of Lagrange interpolation polynomial at points $\tau_0=-k$, $\tau_{1}=-k+1$, \dots, $\tau_k=0$ is denoted $\{\L_0,\dotsc, \L_k\}$ so that the values of $S^j_1$ and $S^j_2$ are given by:
\begin{equation}
S^j_1 = \int_{0}^{1} e^{Z_1(1-\tau)} \L_j(\tau) \dif \tau \cdot Z_{12}
\qquad \text{and} \qquad
S^j_2 = \int_{0}^{1} e^{Z_2(1-\tau)} \L_j(\tau) \dif \tau \cdot Z_{21}.
\end{equation}
The amplification matrix of the implicit multistep scheme is:
\begin{equation}
R_{\mathrm{MS},\mathrm{imp}} =
\begin{pmatrix}
I_{m_1} & S^0_1 & 0 & \cdots & 0 \\
S^0_2 & I_{m_2} & 0 & \cdots & 0 \\
0 & 0 & \ddots & \ddots & \vdots \\
\vdots & & \ddots & I_{m_1} & 0 \\
0 & \cdots & \cdots & 0 & I_{m_2}
\end{pmatrix}^{-1} \cdot
\begin{pmatrix}
e^{Z_1} & S^1_1 & 0 & S^2_1 & \cdots & 0 & S^k_1 \\
S^1_2 & e^{Z_2} & S^2_2 & 0 & \cdots & S^k_1 & 0 \\
I_{m_1} & 0 & \cdots & \cdots & \cdots & 0 & 0 \\
0 & I_{m_2} & \ddots & & & \vdots & \vdots \\
\vdots & \ddots & \ddots & \ddots & & \vdots & \vdots \\
\vdots & & \ddots & I_{m_1} & \ddots & \vdots & \vdots \\
0 & \cdots & \cdots & 0 & I_{m_2} & 0 & 0
\end{pmatrix} \in \mathbf{C}^{mk\times mk}
\end{equation}
where $S_1^j$ and $S_2^j$ are defined as previously, but with the Lagrange basis polynomial at points $\tau_0=-k+1$, $\tau_{1}=-k+2$, \dots, $\tau_k=1$.

In Figures~\ref{fig:conv_multistep:sprad_exp} and~\ref{fig:conv_multistep:sprad_imp}, we show as an example the evolution of the spectral radius of $R_{\mathrm{MS},\mathrm{exp}}$ and $R_{\mathrm{MS},\mathrm{imp}}$ as a function of $\Dt$ for the coupled system from Section~\ref{sec:conv_num}. It is clear that the explicit scheme may be unstable, and the onset of stability matches the divergence of the error in Figure~\ref{fig:conv_multistep:err_exp}. The implicit scheme does not suffer from instabilities in the range of time steps considered.

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{spec_rad_explicit}
\caption{Spectral radius of the explicit scheme}\label{fig:conv_multistep:sprad_exp}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{spec_rad_implicit_a1=-80.0_a12=40.0_a21=40.0_a2=-80.0.png}
\caption{Spectral radius of the explicit scheme}\label{fig:conv_multistep:sprad_imp}
\end{subfigure}
\caption{Spectral radii of the recurrence matrices for the multistep scheme applied to the \testeqdeuxdeux{} with the same parameters as in Figure~\ref{fig:conv_multistep}.}\label{fig:conv_multistep:sprad}
\end{figure}

\subsection{Other schemes from the literature and their amplification matrices}\label{sec:other_schemes}

It is instructive to compare the stability of our multistep coupling scheme to that of other coupling schemes found in the literature. In particular, we focus on splitting schemes and ImEx ARK methods already discussed in the introduction, and derive the associated recurrence matrices.

\subsubsection{Splitting schemes}

In numerous works, the coupling is implemented in a pragmatic manner, where each solver is run sequentially. This is sometimes referred to as a Gauss--Seidel or sequential splitting, or Conventional Serial Staggered (CSS) scheme for instance in the fluid-structure community~\cite{farhat_two_2000}. This actually corresponds to a first-order Lie splitting of the original coupled problem, where equation~\eqref{eq:lin_coupled_lin} is decomposed as:
\begin{equation}
\dif_t y = \underbrace{
\begin{pmatrix} A_1 & A_{12} \\ 0 & 0
\end{pmatrix}}_{\equiv B} y + \underbrace{
\begin{pmatrix} 0 & 0 \\ A_{21} & A_2
\end{pmatrix}}_{\equiv C} y.
\end{equation}
As for the multistep coupling scheme, the study of stability is conducted considering an exact integration of each subsystem, \ie\ only splitting errors are considered. We will use the first-order (Lie) and second-order (Strang) splittings~\cite{strang1968construction}, which still are extensively used in a wide variety of applications~\cite{keyes2013multiphysics}. The corresponding amplification matrices are:
\begin{align}
	R_{\mathrm{Lie}} & = \exp\parens{\Dt C} \exp\parens{\Dt B},
\\	R_{\mathrm{Strang}} & = \exp\parens[\big]{\tfrac{1}{2} \Dt C} \exp\parens{\Dt B} \exp\parens[\big]{\tfrac{1}{2} \Dt C}.
\end{align}
There also exists another version of each splitting, in which the operator $C$ is integrated before~$B$. However, we observed that these schemes behaved identically in terms of stability. Higher-order splitting schemes can only be obtained if complex or negative time steps are used, which is usually impossible with existing solvers, and even more so when applied to dissipative problems. Note that the Strang splitting has been independently ``reinvented'' in the fluid-structure as the Improved Serial Staggered (ISS) scheme~\cite{farhat_two_2000}.

It should be stressed that dynamic adaptation of the coupling time step (splitting time step) can only be done with a large additional cost~\cite{descombes2011adaptive}.

\subsubsection{ImEx ARK schemes} In the context of ImEx ARK schemes, the operator $A$ from equation~\eqref{eq:lin_coupled_lin} is decomposed as follows:
\[
\dif_t y = \underbrace{
\begin{pmatrix} A_1 & 0 \\ 0 & A_2
\end{pmatrix}}_{\equiv M} y + \underbrace{
\begin{pmatrix} 0 & A_{12} \\ A_{21} & 0
\end{pmatrix}}_{\equiv N} y
\]
\ie\ the internal dynamics are integrated implicitly, while the coupling terms are treated explicitly. Codes that already provide an optimised algorithm (\eg\ Newton with preconditioned linear solver) for an implicit scheme can thus be coupled via an ImEx scheme and reuse this algorithm to solve each stage of the implicit part of the ImEx scheme. Such a scheme may also offer time adaptation capabilities~\cite{kennedy_higher-order_2019} with low-cost error estimates.

In the most general case, such methods are built to solve the following type of problem~\cite{sandu_generalized-structure_2015}: $\dif_t y = \sum_{\nu=1}^N F^{(\nu)}(y)$ where $y\in\R^m$, $m\in\N^*$. At the $n$-th time step, for any $y_n\in\R^m$, an $s$-stage ARK scheme is defined by the Butcher tableaus $A^{(\nu)} = \parens{a_{ij}^{(\nu)}}_{1\leq i,j \leq s}$, $b^{(\nu)} = \parens{b_i^{(\nu)}}_{1\leq i \leq s}$ and \smash{$c^{(\nu)} = \parens{c_i^{(\nu)}}_{1\leq i \leq s}$}. A single step is given by:
\begin{equation}\label{eq:ark_def}
\begin{cases}
	\k_i = y_n + \Dt \sum_{\nu=1}^{N} \sum_{j=1}^{s} a_{ij}^{(\nu)} F^{(\nu)}(\k_j)
\\	y_{n+1} = y_n + \Dt \sum_{\nu=1}^{N}\sum_{i=1}^{s} b_i^{(\nu)} F^{(\nu)}(\k_i).
\end{cases}
\end{equation}
The comparisons conducted in this paper are restricted to the case $N=2$ operators. These schemes are called ARK\textsubscript{2} schemes. Therefore, operators $M$ and $N$ given above correspond to $F^{(1)}(y) = My$ and $F^{(2)}(y) = Ny$. The concatenation of the vectors $\k_i$ is denoted $\k = (\k_1^t, \dotsc, \k_s^t)^t$.
For clarity we use the Kronecker product $\otimes$ and the column matrix $\mathbb{1}_s\in\R^{s\times 1}$ filled with ones. Following equation~\eqref{eq:ark_def}:
\[
\k = \underbrace{\Dt \parens[\big]{A^{(1)} \otimes M + A^{(2)} \otimes N}}_{\text{dimension $ms\times ms$}} \k + \mathbb{1}_s \otimes y_n.
\]
This leads to the following result:
\begin{equation}\label{eq:ark_k_vals}
\k = \bracks[\Big]{I_{ms} - \Dt \parens[\big]{A^{(1)} \otimes M + A^{(2)} \otimes N}}^{-1} \parens{\mathbb{1}_s \otimes I_m} y_n.
\end{equation}
Equation~\eqref{eq:ark_def} can be rewritten as $y_{n+1} = y_n + \Dt \parens[\big]{b^{(1)} \otimes M + b^{(2)} \otimes N} \k$, giving the expression of the recurrence matrix for ARK\textsubscript{2} schemes:
\begin{equation}\label{eq:ark2_rec_mat}
y_{n+1} = \underbrace{\braces[\Big]{I_m + \Dt \parens[\big]{b^{(1)} \otimes M + b^{(2)} \otimes N} \bracks[\Big]{I_{ms} - \Dt \parens[\big]{A^{(1)} \otimes M + A^{(2)} \otimes N}}^{-1} \parens{\mathbb{1}_s \otimes I_m}}}_{\equiv R_{\mathrm{ARK}_2}} y_n.
\end{equation}
The multistep coupling scheme will be compared to some state-of-the-art ImEx ARK\textsubscript{2} schemes. The chosen ARK\textsubscript{2} schemes are of orders~3 to~5. They integrate the block-diagonal operator~$M$ (representing the subsystems internal dynamics) with an L-stable, stiffly-accurate, singly diagonally implicit Runge--Kutta method with an explicit first stage (ESDIRK), while they integrate the coupling operator $N$ with an explicit Runge--Kutta method (ERK). The third-order scheme (ARK3(2)4L[2]SA) is taken from~\cite{kennedy_additive_2003} and the fourth- and fifth-order schemes (ARK4(3)7L[2]SA\textsubscript{1} and ARK5(4)8L[2]SA\textsubscript{2}) are found in~\cite{kennedy_higher-order_2019}. A simple first-order ImEx method combining the implicit and explicit Euler schemes is also considered.

\subsection{Stability analysis for the asymmetrical case}\label{sec:kvaerno_stability}

In this section, we study the stability of the multistep coupling scheme on the \testeqdeuxdeux{}~\eqref{eq:two_scalar_coupled} using the first approach suggested in Section~\ref{sec:model_stability}.

The parameter $\gamma = \frac{a_{12}a_{21}}{a_1a_2}$ is considered to be a measure of the strength of the coupling~\cite{kvaerno_stability_2000}, while $\delta$ is the ratio of the internal characteristic times of both subsystems. Figure~\ref{fig:stab_comp_kvaerno} shows the stability limits for the various methods. Each curve denotes the lower limit of the stability domain. The continuous system is only stable for $\gamma<1$, which is used as the upper limit for the vertical axis. We do not consider the zone $z_1=a_1\Dt > 0$ which corresponds to an unstable internal dynamics for subsystem~1. The stability limits of the implicit multistep coupling schemes of order~1 to~4 do not appear in these figures, since they are much wider than those of the other methods.

For a given scheme and a given set $\{z_1, \delta\}$, there exist $1 \geqslant \gamma_+(z_1, \delta) > 0 > \gamma_-(z_1, \delta)$ such that the scheme is stable in the range $\gamma \in [\gamma_-, \gamma_+]$. The most stable schemes have $\gamma^+=1$. The case $\gamma = 0$, which corresponds to a one-way coupling or to the absence of coupling, is stable for all schemes for any $\delta$ and $z_1$. It can also be seen that the width $\gamma_+ - \gamma_-$ of the stability domain tends to $\infty$ as $z_1$ tends to 0.

The implicit multistep scheme has a large stability domain even for high orders. The first-order explicit multistep scheme is comparable to the splittings and the first-order ImEx scheme, however higher-order variants are much less stable. Both splitting schemes have the same stability domain which is competitive with the other explicit or ImEx methods. ImEx schemes show good stability properties when the $\delta$ is small and when the internal dynamics is nonstiff ($z_1$ small). When the internal stiffness increases as $\abs{z_1}$ is increased, the stable range of $\gamma$ is reduced. For a large value of $\delta$, $z_2=\delta z_1$ becomes a very stiff term and causes a large diminution of the stability domain for ImEx schemes. Note that, in Figure~\ref{fig:stab_comp_kvaerno_delta_10} for large values of $\gamma$, \eg\ $-5$, and low values of $\abs{z_1}$, the stability of the ImEx ARK schemes can be comparable or better than the implicit multistep scheme, and much more favourable than the explicit variant.

\begin{figure}[p!]
\begin{subfigure}{0.7\textwidth}
\centering
\includegraphics[trim={0 0 0 1cm},clip,width=\textwidth]{stab_comp_delta_01.png}
\caption{$\delta = 0.1$}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.7\textwidth}
\centering
\includegraphics[trim={0 0 0 1cm},clip,width=\textwidth]{stab_comp_delta=10.0.png}
\caption{$\delta = 10$}\label{fig:stab_comp_kvaerno_delta_10}
\end{subfigure}
\caption{Stability comparison between the implicit and explicit multistep methods, the Lie and Strang splittings and some ImEx schemes (asymmetrical case).}\label{fig:stab_comp_kvaerno}
\end{figure}

\subsection{Stability analysis for the symmetrical case}\label{sec:symmetrical_model_eq_stability}

We now turn to the symmetrical model from Section~\ref{sec:symmetrical_model}, which enables us to study the effect of imaginary components in the coupling terms.

Figure~\ref{fig:stab_comp_sym} shows the stability limit of the studied coupling schemes for two values of the diagonal term $z_1 = a_1 \Dt$. The stability limit is thus a function of the real and imaginary parts of the coupling terms $z_{12} = z_{21}$, and is here normalised by $\abs{z_1}$ to help with the visual comparison. The interior of the stability domain is simply found by noting that it contains the origin $z_{12}=0$. The stability limit is symmetrical with respect to both the imaginary and real parts of $z_{12}$, thus only the positive quarter-plane is shown. The red dotted line $\Real(z_{12})=\abs{{z_{1}}}$ indicates the stability limit of the continuous system.

We may first observe that the stability domains tend to shrink as the order increases, with some exception for the ImEx schemes. The implicit formulation is more stable than the explicit one and some of its stability domains are unbounded for the lowest orders. For both tested values of $z_1$, the implicit stability domains cover the whole stability domain of the continuous system along the real axis at all orders, while the explicit formulation only verifies this property for all $z_1$ at order 1, or at higher-orders but for smaller values of $z_1$ not shown here. This behaviour for purely real coefficients is coherent with the results from the previous section, where it corresponds to the case $\gamma \in [0,1]$.

In the case of the explicit multistep and splitting schemes, the stiffer the internal dynamics is (large negative $z_1$), the smaller the stability domain is in the figure, \ie\ relative to $\abs{z_1}$, which is coherent with our observations in the asymmetrical real case in Section~\ref{sec:kvaerno_stability}. However, in the absolute space (without the scaling $1/\abs{z_1}$) the stability domains are actually growing. The explicit multistep coupling scheme is comparable to the splitting schemes in stability only at order 1, and its stability domain is close to that of the first-order ImEx method.

At fixed order, the implicit stability domains grow when the internal dynamics of the subsystems becomes stiffer (larger $\abs{z_1}$). In Figure~\ref{fig:stab_sym_stiff}, it can even be seen that the implicit stability domains of order~3 to~5 have become unbounded and cover the whole continuous stability domain $\Real(z_{12})<\abs{z_1}$, except for a small portion near the point $\mathrm{Im}(z_{12})=0$, $\Real(z_{12}) = \abs{z_1}$ for higher-order variants. On the opposite, ImEx, explicit multistep coupling and splitting schemes always have bounded stability domains. This observation suggests that a notion of A-stability, adapted from that for ODE integrators, could be defined, requiring that a scheme be stable over the whole stability domain of the continuous system.

\begin{figure}[p!]
\begin{subfigure}{0.7\textwidth}
\centering
\includegraphics[trim={0 0 0 1cm},clip,width=\textwidth]{stab_comp_z1=-1.0.png}
\caption{$z_1 = -1$}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.7\textwidth}
\centering
\includegraphics[trim={0 0 0 1cm},clip,width=\textwidth]{stab_comp_z1=-10.0.png}
\caption{$z_1 = -10$}\label{fig:stab_sym_stiff}
\end{subfigure}
\caption{Stability comparison between the implicit and explicit multistep methods, the Lie and Strang splittings and some ImEx schemes (symmetrical case). The red dotted line $\Real(z_{12})=\abs{{z_{1}}}$ indicates the stability limit of the continuous system.}\label{fig:stab_comp_sym}
\end{figure}

\subsection{Discussion}

Both previous sections have shown that the stability domains may vary strongly between the different schemes and also depending on the choice of parameters. The choice of a coupling scheme should be driven by multiple aspects: stability, efficiency, ease of use, ease of implementation. In this section, we discuss each scheme and indicate situations where the multistep coupling scheme is a valuable choice.

The Lie splitting scheme is the simplest scheme available, as it is the easiest to implement, requiring only a sequential call to each solver in turn (Gauss--Seidel coupling). It is slightly more stable than the first-order explicit multistep scheme, where the solvers can however be called in parallel (Jacobi coupling), which may greatly improve the computational efficiency. The Strang splitting is an extension to order 2 which has no overcost, since it is essentially a shifted variant of the Lie splitting. It retains the same stability domain while improving accuracy, and should thus be a scheme of choice for applications requiring simple coupling algorithms.

Higher-order explicit multistep schemes can be seen as high-order generalisation of parallel (Jacobi) coupling. The stability domain is reduced, however the accuracy is greatly increased. Note that the explicit multistep coupling schemes provide natural error estimates to guarantee a certain level of accuracy on the approximation of the coupling variables~\cite{Francois-2022,francois_multistep_2023}. This is a strong advantage compared to splitting schemes, which do not offer time adaptation natively, or only at a large additional cost~\cite{descombes2011adaptive}. Indeed, this allows a user to directly specify an error tolerance, without having to try different coupling time steps to assess the temporal convergence of a simulation.

The chosen ImEx schemes perform well in terms of stability, and have a stability up to~5 times larger than the explicit multistep coupling scheme for the higher-order variants. The next section will somewhat mitigate this observation for certain configurations commonly encountered in practice. ImEx schemes are however closer in their principle to monolithic methods. Indeed, every subsystem must evolve with a single time step, \ie\ no substepping is possible, which makes it difficult to accurately and efficiently couple systems with very different time scales. Furthermore, they require that each solver possesses a nonlinear system solution algorithm (Newton or similar) for the solution of the implicit part of the scheme. Thus, solvers only using explicit time integrators internally cannot be used, and solvers possessing an implicit scheme must however be heavily modified to permit a synchronised temporal advancement of all solvers at each stage of the ImEx scheme. Moreover, the implicit Runge--Kutta method used in the ImEx scheme may not be suited to all physics considered in a given multiphysics problem.

The implicit multistep scheme generally offers the best stability of all studied scheme, and is the only one to offer unbounded stability domains in certain configurations. This however comes at the cost of having to solve a potentially nonlinear system on the coupling variables at each step. For nonstiff couplings, a simple relaxation procedure can suffice. If the coupling is stronger, Newton-like methods may be needed. For surface-coupled systems, interface quasi-Newton method have emerged recently~\cite{degroote2013partitioned}, where a Jacobian of the coupling residuals is iteratively approximated from initial relaxation iterations. These methods have proved valuable for many applications, in particular fluid-structure interaction. However, there is still much research to be conducted on the topic of efficient solution algorithms for the implicit coupling.

The multistep coupling scheme overall has acceptable stability limits in its explicit variant, and much better ones in its implicit form. It seems to be a valuable choice for many applications, due to its ease of implementation. Indeed, solvers can still use their existing physics-optimised time integrator and are free to subcycle. Modifications to be made to the codes are minor, mostly requiring the ability to update the coupling terms via the polynomial predictors $\hat{u}_i^n(t)$ during the solver's subcycles. Another advantage is that the computational cost remains the same at all orders, similarly to multistep ODE integrators.

The implementation of the high-order coupling procedure in an HPC context is being tackled in the coupling library \textsc{Cwipi}~\cite{cwipi_github}, with the aim of providing an easy-to-use interface for large scale and industrial simulations.

\section{A conjugate heat transfer test case}\label{sec:cht}

The previous section has proposed an analysis of the coupling stability for the simple \testeqdeuxdeux{}~\eqref{sec:model_equation}. In this section, we consider a more complex coupled model to determine how the stability of each scheme is affected. We study the application of the multistep coupling strategy to a conjugate heat transfer (CHT) test case, with two solid slabs of finite length exchanging heat through an interface. This case is recognized as a representative configuration for multidimensional fluid-solid thermal interactions, in particular with respect to the stability restrictions on the time step. Indeed, the stability limit is usually controlled by the diffusive processes near the coupling surface~\cite{giles1997stability}, whose stiffness is linked to the mesh refinement in the direction perpendicular to the coupling surface.

This classical problem~\cite{francois_multistep_2023, Francois-2022} is briefly recalled for completeness. Figure~\ref{fig:schema_cht1d} provides a sketch of the configuration. For simplicity, both slabs have the same length, same space discretisation (with a second-order finite-volume approach), and same thermal properties: $\lambda_1=\lambda_2=\qty{1}{W/m/K}$, $c_1=c_2=\qty{1}{J/kg/K}$. They are only different in their densities $\rho_1$ and $\rho_2$. We set $\rho_1=\qty{1}{kg/m}$. Each slab has a length $L=\qty{1}{m}$ and is discretised uniformly with 50~cells.

\begin{figure}[h]
\includegraphics[width=0.8\textwidth]{cht_1d.png}
\caption{Sketch of the 1D conjugate heat transfer case.}\label{fig:schema_cht1d}
\end{figure}

The coupling is of the Neumann--Dirichlet type, with the first slab receiving a heat flux prescribed by the second solver (extrapolated to the interface by a one-sided second-order scheme), and the second slab being imposed an interface temperature prescribed by the first solver in a similar fashion. Taking each slab as a standalone system with a dedicated solver, we obtain a coupled system representative of the CHT between a fluid solver and a heat diffusion solver. Following the same approach as in Section~\ref{sec:rec_mac_multi}, we construct for each coupling scheme the recurrence matrix for the CHT problem, with the submatrices $A_{ij}$ corresponding to the terms arising from the finite-volume discretisation~\cite{Francois-2022}. By varying the time step through a dichotomy process, we can determine the corresponding maximum stable time step of a coupling scheme applied to this configuration.

Existing studies on the stability of this problem show for the simple first-order coupling scheme that the ratio of effusivities is the main parameter impacting stability~\cite{giles1997stability}, which in our case amounts to the ratio of densities. Thus, we propose to determine the maximum stable time step for each coupling method as a function of the ratio $\rho_1/\rho_2 \in [10^{-2}, 10^{3}]$, using a dichotomy process for each ratio and method to find $\Dt$ such that the spectral radius of the associated recurrence matrix is 1.

The results are shown in Figure~\ref{fig:dt_lim_cht_exp} for the explicit multistep coupling scheme, and in Figure~\ref{fig:dt_lim_cht_imp} for the implicit scheme. Both figures include the stability limits of the ImEx schemes and splitting schemes from Section~\ref{sec:other_schemes}. We also include the stability limit of the explicit Euler scheme, which indicates how a standard monolithic explicit scheme would perform.

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{dt_lim_stab_comparison_explicit_}
\caption{Explicit formulation}\label{fig:dt_lim_cht_exp}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{dt_lim_stab_comparison_implicit_.png}
\caption{Implicit formulation}\label{fig:dt_lim_cht_imp}
\end{subfigure}
\caption{Maximum stable time steps for the 1D CHT case.}\label{fig:dt_lim_cht}
\end{figure}

It can be noted that the stability limits obtained for the ImEx schemes and the explicit multistep coupling schemes closely reproduce those obtained in~\cite{francois_multistep_2023}, where the Python demonstrator \textsc{Rhapsopy}~\cite{rhapsopy_github} of the multistep coupling has been used to perform several fixed-time-step integrations, and stability limits have been found based on a divergence criterion for the obtained solutions. For the implicit multistep scheme however, high-order variants have obtained here a lower stability limit. This may be attributed to the low-order initialisation (first steps) of the simulations from the previous work, which may improve the perceived numerical stability of the method.

Both splittings have the same stability limit, which is virtually infinite for density ratios above~$0.5$. We observe that the explicit multistep coupling scheme is competitive for large values of $\rho_1/\rho_2$, and that an increase in order yields a decrease in stability, as seen with the test equation from Section~\ref{sec:symmetrical_model_eq_stability}. ImEx schemes have a relatively low stability limit which is only advantageous for intermediate density ratios. This could indicate that the explicit multistep scheme has a lower stiffness leakage~\cite{kennedy_additive_2003} than the chosen ImEx schemes, \ie\ stiffness from the internal dynamics of the subsystems may leak more into the explicit part under certain conditions, causing stronger stability restrictions.

The implicit scheme has, for all orders, a much larger stability limit, even virtually unlimited for schemes of order up to 2, or at all orders for the higher density ratios. This clearly demonstrate the advantage of the implicit formula. Note that the interface boundary conditions should be reversed when $\rho_1<\rho_2$~\cite{giles1997stability}, hence the implicit scheme can be seen to remain robust even when the boundary conditions choice is inadequate and favours instability.

For most gas-thermal coupling, the density ratio is above~100, hence the explicit multistep scheme or the Strang splitting would be a good choice, owing to their simplicity of use and implementation. For stronger coupling where the density ratio is closer to~1, \eg\ liquid-solid thermal coupling, the multistep implicit coupling scheme would be preferable. The splitting schemes may also be advantageous due to their explicit nature. A study on the achieved precision should be conducted to discriminate between these methods.

This case, which is representative of a category of realistic multiphysics applications where diffusion dominates the coupling process, shows that the multistep coupling scheme in its implicit variant has a very good stability. The explicit variant can be competitive for intermediate density ratios, and very stable for large ratios. The chosen ImEx scheme, although they are stable for larger time steps than a standard monolithic explicit scheme, are only relevant for intermediate values of the density ratio.

Note that, in the context of adaptive time stepping based on error estimates, the coupling time step will remain within or at the border of the stability region, thus preventing instability. However, for a poorly stable scheme, the time step may end up being dictated by the stability limit instead of the accuracy requirement, greatly reducing the computational efficiency~\cite[Section IV.2]{hairer_book2}. Therefore, stability considerations may be helpful in diagnosing such an issue and in choosing a well-suited coupling scheme.

The conjugate heat transfer highlights the strong influence of physical parameters on the stability of the schemes. It would be interesting to search for a link between this complex coupled system of PDEs and the simple \testeqdeuxdeux{}\ (see Section~\ref{sec:gen_dahlquist_and_2x2}). We are currently investigating whether a set of $2\times2$ test equations can be generated from the complete CHT system, in such a way that each of these simpler models could be located within the stability domains explored in Section~\ref{sec:model_stability} to assess the stability of the original CHT system. Such an approach would have the capability of predicting stability a priori, offering a very efficient way to choose the proper time step.

\section{Conclusion}\label{sec:conclusion}

In this work, a thorough analysis of the essential convergence properties of the multistep coupling scheme has been presented, extending results initiated in~\cite{francois_multistep_2023}. Inspired by Dahlquist's test equation for the stability analysis of ODEs, a \testeqdeuxdeux{} has been proposed. Relying on a recurrence matrix formulation of the multistep coupling scheme, the stability limit can be studied systematically by simple considerations on the spectral radius of the matrix. The analysis shows that the multistep scheme stability is good compared to other schemes, but is reduced as the polynomial order is increased. Implicit coupling provides much greater stability than the explicit variant. State-of-the-art ARK ImEx methods show good stability limits, surpassing the stability of the multistep coupling scheme at higher orders, however at the cost of being closer to a monolithic integration scheme, thus more demanding in terms of implementation, while prohibiting subcycling. Low-order Lie and Strang splittings show a very good performance in terms of stability, thus providing interesting coupling schemes for applications less affected by coupling accuracy. Application on a 1D conjugate heat transfer test case shows that the tendencies obtained with the \testeqdeuxdeux{} convey into the analysis of a system representative of actual multiphysics applications.

Overall, the multistep coupling scheme is a valuable coupling scheme. It forms a natural extension of the basic parallel coupling schemes often used in practice, only requiring minor modifications in the codes to be used. The order of convergence in time can be arbitrary, and an increase in order does not result in an increase in computational cost, as it only requires storing more past values of the coupling variables. Although not discussed in this article, this scheme allows for error estimates to be simply derived from the polynomials, enabling a dynamic adaptation of the coupling time step to guarantee accuracy and stability, while also optimising the computational cost.

The multistep coupling method is the subject of ongoing work. Regarding the stability analysis, a quantitative link should be made between the \testeqdeuxdeux{}~\eqref{eq:lin_coupled_lin} and actual coupled partial differential equations, so as to allow an a priori analysis of the coupling stability. The multistep scheme is being implemented in the open-source HPC library \textsc{Cwipi}~\cite{refloch2011cedre,cwipi_github} with the aim of being used in a wide range of applications. Such a library will provide an easy-to-use package to easily incorporate a high-order adaptive multistep coupling scheme into new or existing coupled applications. Finally, the optimisation of the computational cost of the implicit variant will be a subject of specific work, in particular regarding the application of specific Newton-like algorithms for the solution of the nonlinear system obtained at each step.

\appendix

\section{Local truncation error constants}\label{appendix:error}

Let us consider the linear problem from equation~\eqref{eq:lin_coupled_lin}. It can be reformulated as:
\begin{equation}
\dif_t y = A_{\intern} y + A_{\mathrm{cpl}} y
\end{equation}
with the internal dynamics operator $A_{\intern} \equiv
\parens[\big]{\begin{smallmatrix} A_1 & 0 \\ 0 & A_2 \end{smallmatrix}}$ and the coupling operator $A_{\mathrm{cpl}} \equiv
\parens[\big]{\begin{smallmatrix} 0 & A_{12} \\ A_{21} & 0 \end{smallmatrix}}$.
This corresponds to the case $y=(y_1^t, y_2^t)^t$, $u_1=y_2$ and $u_2=y_1$.

Over a coupling time step $[t_n, t_{n+1}]$, the multistep scheme corresponds to the integration of the following modified equation:
\begin{equation}
\dif_t y = A_{\intern} y + A_{\mathrm{cpl}} \hat{y}_n
\end{equation}
where $\hat{y}_n$ is the polynomial approximation of $y$:
\begin{equation}
\hat{y}_n(t) = \sum_{j=\delta k}^{k+\delta k} y_{n-j} \prod_{\substack{l=\delta k \\ l \neq j}}^{k+\delta k} \dfrac{t-t_{n-l}}{t_{n-j}-t_{n-l}}
\end{equation}
with $k$ the degree of the interpolation polynomial, and $\delta k=0$ in the explicit case, -1 in the implicit case (so that $t_{n+1}$ is a sampling point of the interpolant).

The exact solution $y_{\exa}$ and the multistep solution $y_{n+1}$ read:
\begin{equation}
\begin{aligned}
	y_{\exa}(t_{n+1}) & = e^{A_\intern\Dt}y_{\exa}(t_n) + \int_{0}^{1} e^{A_{\intern}\Dt(1-\eta)} A_{\mathrm{cpl}}\Dt y_{\exa}(t_n +\eta \Dt) \dif \eta
\\	y_{n+1} & = e^{A_{\intern}\Dt}y_n + \int_{0}^{1} e^{A_{\intern}\Dt(1-\eta)} A_{\mathrm{cpl}}\Dt\hat{y}^n(t_n +\eta \Dt) \dif \eta
\end{aligned}
\end{equation}
with $\Dt=t_{n+1}-t_{n}$.

The local truncation error can be obtained by applying the multistep coupling scheme to the exact solution, \ie\ using the exact solution to generate the interpolation polynomial $\hat{y}_n$ and to set the initial condition $y_n=y_{\exa}(t_n)$. Subtracting both solutions yields the local truncation error $\delta y_{n+1} \equiv y_{n+1} - y_{\exa}(t_{n+1})$ for a single coupling step:
\begin{equation}
\delta y_{n+1} = \int_{0}^{1} e^{A_{\intern}\Dt(1-\eta)} A_{\mathrm{cpl}}\Dt \underbrace{\bracks{y_{\exa} - \hat{y}_n}(t_n +\eta \Dt)}_{\equiv \epsilon} \dif \eta.
\end{equation}
Since $\hat{y}_n$ interpolates the exact solution, its interpolation error can be expressed as:
\begin{equation}
\epsilon(t_n + \eta \Dt) = \dfrac{\Dt^{k+1} y_{\exa}^{(k+1)}(t_n)}{(k+1)!} \prod_{j=\delta k}^{k+\delta k} (\eta + j) + \mathcal{O}(\Dt^{k+2}).
\end{equation}
Then, the local truncation error reads:
\begin{equation}
\delta y_{n+1} = \int_{0}^{1} e^{A_{\intern}\Dt(1-\eta)} A_{\mathrm{cpl}} \bracks*{{\dfrac{\Dt^{k+2} y_{\exa}^{(k+1)}(t_n)}{(k+1)!} \prod_{j=\delta k}^{k+\delta k} (\eta + j)} + \mathcal{O}(\Dt^{k+3})} \dif \eta.
\end{equation}
Since we are interested in the asymptotic convergence regime where $\Dt$ is small, we can expand the exponential as $e^{A_{\intern}\Dt(1-\eta)} = I + A_{\intern}\Dt(1-\eta) + \mathcal{O}(\Dt^2)$. The first-order term $A_{\intern}\Dt(1-\eta)$ multiplied by the polynomial of $\eta$ produces a contribution of order $k+3$ that is absorbed by the $\mathcal{O}(\Dt^{k+3})$ term, which allows for the following simplified form:
\begin{equation}
\delta y_{n+1} = \Dt^{k+2} \underbrace{\cdot \dfrac{1}{(k+1)!} \cdot {\int_{0}^{1} \prod_{j=\delta k}^{k+\delta k} (\eta + j) \dif \eta}}_{\text{error constant}} \cdot A_{\mathrm{cpl}} \cdot y_{\exa}^{(k+1)}(t_n) + \mathcal{O}(\Dt^{k+3}).
\end{equation}
This simplification is coherent with the fact that the internal dynamics is assumed to be integrated exactly, so that its accuracy can only be indirectly worsened due to the polynomial approximation of the coupling term. In particular, when there is no coupling, $A_{\mathrm{cpl}}=0$ and the error is 0.

The error constants exactly correspond to those of the Adams--Bashforth multistep schemes in the explicit coupling case, and to those of the Adams--Moulton methods in the implicit case~\cite[Section~5.1.1]{ascher1998computer}. This is coherent, since the simplification following the series expansion has the same effect as setting $A_\intern=0$, in which case the multistep coupling scheme reduces to an Adams multistep scheme. Note that, in this particular setting, $y_{\exa}'= A_{\mathrm{cpl}} y$, therefore the error can be reformulated in the more traditional linear multistep form~\cite[equation~(5.2)]{ascher1998computer}:
\begin{equation}
\delta y^{n+1} = \dfrac{1}{(k+1)!} \cdot \int_{0}^{1} \prod_{j=\delta k}^{k+\delta k} (\eta + j) \dif \eta \cdot \Dt^{k+2} \cdot y_{\exa}^{(k+2)}(t_n) + \mathcal{O}(\Dt^{k+3})
\end{equation}
where the operators prescribing the system dynamics no longer appear.

The predicted ratios of error constants between the implicit and explicit multistep coupling schemes of the same order are listed in Table~\ref{table:ratio_err}, alongside the numerically evaluated ratios based on fits of the asymptotic parts of the convergence curves from Figure~\ref{fig:conv_multistep}. The ratios indicate that the implicit coupling scheme is increasingly more accurate than its explicit counterpart as the order is increased. This can be understood by the behaviour of the error polynomial (main term of $\epsilon$), which diverges monotonically on $[t_n, t_{n+1}]$ in the extrapolation case, while it is bounded on that interval in the interpolation case (implicit coupling).

\begin{table}[h]
\caption{Ratio of the implicit error constant divided by the explicit error constant, as predicted by the theory, and as evaluated numerically from Figure~\ref{fig:conv_multistep}.}
\begin{tabular}{|c|c|c|}
\hline
\\[-2.7ex]
\hline
\textbf{Order} & \textbf{Theory} & \textbf{Figure~\ref{fig:conv_multistep}}
\\
\hline
\\[-2.7ex]
\hline
1 & 1 & 1.0 \\\hline 2 & 5 & 5.2 \\\hline 3 & 9 & 9.4 \\\hline 4 & 13.2 & 13.5 \\\hline 5 & 17.6 & 17.7 \\\hline 6 & 22.1 & 22.8 \\\hline
\end{tabular}\label{table:ratio_err}
\end{table}

\section*{Acknowledgments}

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

The authors would like to thank Josselin Massot and Nicole Spillane (CMAP) for the reading of the draft of the paper and fruitful discussions.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}