## 2 Materials and Methods
### Solution  
To solve the system of SIDEs (2)-(7), we used the first-order integration method, generalizing the Euler method for stochastic differential equations (finite difference approach). For memory integral integration in the arbitrary viscous torque evaluation of the acceleration torque, an order one quadrature scheme, similar to those presented in [50], is generally used. However, in contrast, we have used a different kernel function (6) arising from the definition of viscous torque equation (5). The integral occurring in this equation can be transcribed using integration by parts:  
\[\int\limits_{t_{0}}^{t}K(t-\tau)\frac{\mathrm{d}\omega(\tau)}{\mathrm{d}\tau} \,\mathrm{d}\tau+K(t-t_{0})\omega(t_{0})=\frac{\mathrm{d}}{\mathrm{d}t}\int \limits_{t_{0}}^{t}K(t-\tau)\omega(\tau)\mathrm{d}\tau. \tag{8}\]  
Now, we can divide the time span by the sequence of \(n-1\) constant timestep intervals \(h=\tau_{i+1}-\tau_{i}\) for \(i=1\) to \(n-1\), where \(\tau_{1}=t_{0}\) and \(\tau_{n}=t\), which gives:  
\[\frac{\mathrm{d}}{\mathrm{d}t}\int\limits_{t_{0}}^{t}K(t-\tau)\omega(\tau) \mathrm{d}\tau=\frac{\mathrm{d}}{\mathrm{d}t}\sum\limits_{i=1}^{n-1}\int \limits_{\tau_{i}}^{\tau_{i+1}}K(t-\tau)\omega(\tau)\mathrm{d}\tau\equiv \frac{\mathrm{d}}{\mathrm{d}t}\sum\limits_{i=1}^{n-1}I_{i}(t). \tag{9}\]  
If we now examine the simplest case, a linear approximation, the calculation leads to an order one quadrature scheme. By approximating \(\omega(\tau)\) linearly in the interval \(\tau\in[\tau_{i},\,\tau_{i+1}]\):  
\[\omega(\tau)=\omega(\tau_{i})+\frac{\omega(\tau_{i+1})-\omega(\tau_{i})}{h}( \tau-\tau_{i})+\mathcal{O}(h^{2}), \tag{10}\]  
we obtain:  
\[I_{i}(t)\equiv\int\limits_{\tau_{i}}^{\tau_{i+1}}K(t-\tau)\omega(\tau)\, \mathrm{d}\tau=\left[\omega(\tau_{i})+\mathcal{O}(h^{2})\right]\int\limits_{0 }^{h}K(t-\tau_{i}-\tau)\,\mathrm{d}\tau+\frac{\omega(\tau_{i+1})-\omega(\tau_ {i})}{h}\int\limits_{0}^{h}\tau K(t-\tau_{i}-\tau)\,\mathrm{d}\tau. \tag{11}\]  
The integrals in Equation (11), with kernel function (6), can be computed analytically to yield:\[I_{i}(t) \approx\omega(\tau_{i})\sqrt{\pi}\frac{R^{2}}{\
u}\left[\exp\frac{ \
u(t-\tau_{i}-\tau)}{R^{2}}\operatorname{erfc}\frac{\sqrt{\
u(t-\tau_{i}-\tau)} }{R}\right]_{\tau=0}^{h}+\] \[+\frac{\omega(\tau_{i+1})-\omega(\tau_{i})}{h}\left[\frac{2R^{3}}{ \
u^{3/2}}\sqrt{t-\tau_{i}-\tau}+\sqrt{\pi}\bigg{(}\frac{R^{4}}{\
u^{2}}+\frac {R^{2}}{\
u}\tau\bigg{)}\exp\frac{\
u(t-\tau_{i}-\tau)}{R^{2}}\operatorname{ erfc}\frac{\sqrt{\
u(t-\tau_{i}-\tau)}}{R}\right]_{\tau=0}^{h}, \tag{12}\]  
if higher orders of \(h\), than the first, are omitted. A method that considers the stochastic term (7) in the solution of stochastic differential equations (SDEs), using a finite difference approach was shown in [42] and also used in [33]. The approach utilizes a discrete sequence of random numbers \(W_{i}\) that mimics the properties of \(W(t)\) and is stationary with zero mean, as well as having to fulfill \(\langle W(t)^{2}\rangle=1\) for each value of \(t\), which in a discrete sequence sense means that \(\langle(W_{i}\Delta t)^{2}\rangle/\Delta t=1\), i.e., the \(W_{i}\) has variance \(1/\Delta t\). Furthermore, because \(W(t)\) is uncorrelated, we assume \(W_{i}\) and \(W_{j}\) to be independent for \(i\
eq j\), i.e., we use a sequence of uncorrelated random numbers with zero mean and variance \(1/\Delta t\). The realization of a such sequence is as simple as:  
\[W_{i}=\frac{w_{i}}{\sqrt{\Delta t}}\equiv\frac{w_{i}}{\sqrt{h}}, \tag{13}\]  
where \(w_{i}\) is a Gaussian random number with zero mean and unit variance. For the timestep \(\Delta t\), we used the standard notation \(h\), as used above.