## 2 Numerical Model
### Formulations
#### 2.1.3 Time-Domain  
Quasi-stationary time-domain magnetic problems are solved in FLARE with a MVP formulation under the quasi-magnetostatic approximation (\(\partial\bar{D}/\partial t=0\), \(\partial\bar{B}/\partial t\
eq 0\)):  
\[-\
abla\cdot\left(\frac{1}{\mu_{r}}\
abla A_{z}\right)+\sigma\mu_{0}\frac{ \partial A_{z}}{\partial t}=\mu_{0}(-\sigma\
abla\varphi). \tag{3}\]  
Note that the previous expression--which is the time-domain version of (2)--reduces to (1) when \(\partial/\partial t=0\). The discretized problem is  
\[[K][A_{z}]+[S]\left[\frac{dA_{z}}{dt}\right]=[T], \tag{4}\]  
where \([A_{z}]\), \(\left[\frac{dA_{z}}{dt}\right]\) and \([T]\) are arrays containing the nodal values of the MVP, its derivative, and the right-hand side of (3), respectively. \([K]\) and \([S]\) are assembled from the respective local matrices. Details on the employed vectorized assembly routine are provided in Section 2.2. As anticipated, two different strategies can be adopted by the user to obtain the time-domain solution. Time-Stepping Method  
The first strategy is to use a time-marching strategy, based on the well-known Crank-Nicolson implicit method. The idea is to discretize the time-derivative in (4) by means of a centered finite difference scheme, as opposed to the forward or backward finite difference that would be employed for an explicit or implicit Euler scheme, respectively. In this way, one has for the given node \(i\):  
\[\frac{d}{dt}\left(A_{z,i}^{(k+1/2)}\right)\approx\frac{A_{z,i}^{(k+1)}-A_{z,i }^{(k)}}{\Delta t}, \tag{5}\]  
where  
\[A_{z,i}^{(k+1/2)}\approx\frac{1}{2}\left(A_{z,i}^{(k+1)}+A_{z,i}^{(k)}\right) \tag{6}\]  
Substituting (5) and (6) in (4), one finds:  
\[\frac{1}{2}[K]\left([A_{z}]^{(k+1)}+[A_{z}]^{(k)}\right)+\frac{1}{\Delta t}[S ]\left([A_{z}]^{(k+1)}-[A_{z}]^{(k)}\right)=[T]^{(k+1/2)}. \tag{7}\]  
Then, factoring out the terms at time instants \((k+1)\) and \((k)\) one obtains the final expression:  
\[[M_{1}][A_{z}]^{(k+1)}=[T]^{(k+1/2)}+[M_{2}][A_{z}]^{(k)}, \tag{8}\]  
where  
\[[M_{1}]=\left[\frac{1}{2}[K]+\frac{1}{\Delta t}[S]\right];\quad[M_{2}]=\left[ \frac{1}{\Delta t}[S]-\frac{1}{2}[K]\right] \tag{9}\]  
The Crank-Nicolson scheme is unconditionally stable and second-order accurate in time [38]. The code can be run using either constant or variable time step lengths. When a constant time step is selected the code uses the decompose MATLAB(r) function on the matrix \([M_{1}]\) in (8). decompose performs an automatic selection between several matrix decomposition algorithms, based on the properties of the provided matrix. We refer the reader to the MATLAB(r) guide for further details on the adopted procedure [39]. While the decomposition procedure can be time-intensive for matrices with large numbers of non-zero entries, we found it to be quite efficient for the sizes considered in this work. For example, decompose took \(\approx\)0.41 s for the sparse matrix that will be considered in Section 3, having rank \(\approx\)7 \(\times\) 104 and \(\approx\)5 \(\times\) 105 non-zero entries. In any case, the decomposition of thematrix is performed only once per simulation, and used at each step of the time-loop to compute \([A_{z}]^{(k+1)}\) without repeating the decomposition phase of the (direct) linear system solving procedure. This, at least for the use-cases explored by the authors in this work, allows for an execution time reduction of approximately one order of magnitude when using decompose. Conversely, when a variable time step is selected, the \(\Delta t\) in the left-hand side of the time-discrete version of (8) changes at each iteration. Hence, decompose cannot be used in this case. Nevertheless, using an adaptive time step can be extremely advantageous for exponentially varying sources.