**1. Introduction**

In this paper, we are concerned with random delay differential equations, defined as classical delay differential equations whose inputs (coefficients, forcing term, initial condition, ...) are considered as random variables or regular stochastic processes on an underlying complete probability space (Ω, <sup>F</sup>, <sup>P</sup>), which may take a wide variety of probability distributions, such as Binomial, Poisson, Gamma, Gaussian, etc.

Equations of this kind should not be confused with stochastic differential equations of Itô type, forced by an irregular error term called White noise process (formal derivative of Brownian motion). In contrast to random differential equations, the solutions to stochastic differential equations exhibit nondifferentiable sample-paths. See [1] (pp. 96–98) for a detailed explanation of the difference between random and stochastic differential equations. See [1–6], for instance, for applications of random differential equations in engineering, physics, biology, etc. Thus, random differential equations require their own treatment and study: they model smooth random phenomena, with any type of input probability distributions.

From a theoretical viewpoint, random differential equations may be studied in two senses: the sample-path sense or the L*p*-sense. The former case considers the trajectories of the stochastic processes involved, so that the realizations of the random system correspond to deterministic versions of the problem. The latter case works with the topology of the Lebesgue space (L*p*, ·*p*) of random variables with finite absolute *<sup>p</sup>*-th moment, where the norm ·*<sup>p</sup>* is defined as: *U<sup>p</sup>* <sup>=</sup> <sup>E</sup>[|*U*<sup>|</sup> *p*] 1/*p* for 1 <sup>≤</sup> *<sup>p</sup>* <sup>&</sup>lt; <sup>∞</sup> (<sup>E</sup> denotes the expectation operator), and *U*<sup>∞</sup> <sup>=</sup> inf{*<sup>C</sup>* <sup>≥</sup> <sup>0</sup> : <sup>|</sup>*U*| ≤ *<sup>C</sup>* almost surely} (essential supremum), *<sup>U</sup>* : <sup>Ω</sup> <sup>→</sup> <sup>R</sup> being any random variable. The Lebesgue space (L*p*, ·*p*) has the structure of a Banach space. Continuity, differentiability, Riemann integrability, etc., can be considered in the aforementioned space L*p*, which gives rise to the random L*p*-calculus.

In order to fix concepts, given a stochastic process *<sup>x</sup>*(*t*) <sup>≡</sup> *<sup>x</sup>*(*t*, *<sup>ω</sup>*) on *<sup>I</sup>* <sup>×</sup> <sup>Ω</sup>, where *<sup>I</sup>* <sup>⊆</sup> <sup>R</sup> is an interval (notice that as usual the *ω*-sample notation might be hidden), we say that *x* is <sup>L</sup>*p*-continuous at *<sup>t</sup>*<sup>0</sup> <sup>∈</sup> *<sup>I</sup>* if lim*h*→<sup>0</sup> *x*(*t*<sup>0</sup> <sup>+</sup> *<sup>h</sup>*) <sup>−</sup> *<sup>x</sup>*(*t*0)*<sup>p</sup>* <sup>=</sup> 0. We say that *<sup>x</sup>* is L*p*-differentiable at *<sup>t</sup>*<sup>0</sup> <sup>∈</sup> *<sup>I</sup>* if lim*h*→<sup>0</sup> *<sup>x</sup>*(*t*0+*h*)−*x*(*t*0) *<sup>h</sup>* − *x* (*t*0)*<sup>p</sup>* = 0, for certain random variable *x* (*t*0) (called the derivative of *x* at *t*0). Finally, if *I* = [*a*, *b*], we say that *x* is L*p*-Riemann integrable on [*a*, *b*] if there exists a sequence of partitions {*Pn*}<sup>∞</sup> *<sup>n</sup>*=<sup>1</sup> with mesh tending to 0, *Pn* = {*a* = *t n* <sup>0</sup> < *t n* <sup>1</sup> < ... < *t n rn* = *b*}, such that, for any choice of points *s<sup>n</sup> <sup>i</sup>* ∈ [*t n <sup>i</sup>*−1, *<sup>t</sup> n <sup>i</sup>* ], *<sup>i</sup>* <sup>=</sup> 1, ... ,*rn*, the limit lim*n*→<sup>∞</sup> <sup>∑</sup>*rn <sup>i</sup>*=<sup>1</sup> *<sup>x</sup>*(*s<sup>n</sup> <sup>i</sup>* )(*t n <sup>i</sup>* − *t n <sup>i</sup>*−1) exists in L*p*. In this case, these Riemann sums have the same limit, which is a random variable and is denoted by *b <sup>a</sup> x*(*t*) d*t*.

This L*p*-approach has been widely used in the context of random differential equations with no delay, especially the case *p* = 2 which corresponds to the Hilbert space L<sup>2</sup> and yields the so-called mean square calculus; see [5,7–15]. Only recently, a theoretical probabilistic analysis of random differential equations with discrete constant delay has been addressed in [16–18]. In [16], general random delay differential equations in L*<sup>p</sup>* were analyzed, with the goal of extending some of the existing results on random differential equations with no delay from the celebrated book [5]. In [17], we started our study on random delay differential equations with the basic autonomous-homogeneous linear equation, by proving the existence and uniqueness of L*p*-solution under certain conditions. In [18], the authors studied the same autonomous-homogeneous random linear differential equation with discrete delay as [17], but considered the solution in the sample-path sense and computed its probability density function via the random variable transformation technique, for certain forms of the initial condition process. Other recent contributions for random delay differential equations, but focusing on numerical methods instead, are [19–21].

There is still a lack of theoretical analysis for important random delay differential equations. Motivated by this issue, the aim of this contribution is to advance further in the theoretical analysis of relevant random differential equations with discrete delay. In particular, in this paper we extend the recent study performed in [17] for the basic linear equation by adding a stochastic forcing term:

$$\begin{cases} \mathbf{x}'(t,\omega) &= \mathbf{a}(\omega)\mathbf{x}(t,\omega) + \mathbf{b}(\omega)\mathbf{x}(t-\tau,\omega) + f(t,\omega), \; t \ge 0, \; \omega \in \Omega, \\\ \mathbf{x}(t,\omega) &= \mathbf{g}(t,\omega), \; -\tau \le t \le 0, \; \omega \in \Omega. \end{cases} \tag{1}$$

The delay *τ* > 0 is constant. The coefficients *a* and *b* are random variables. The forcing term *f*(*t*) and the initial condition *g*(*t*) are stochastic processes on [0, ∞) and [−*τ*, 0] respectively, which depend on the outcome *ω* ∈ Ω of a random experiment which might be sometimes omitted in notation. The term *x*(*t*) represents the differentiable solution stochastic process in a certain probabilistic sense. Formally, according to the deterministic theory [22], we may express the solution process as

$$\begin{split} \mathbf{x}(t,\omega) &= \mathbf{e}^{a(\omega)(t+\tau)} \mathbf{e}^{b\_{1}(\omega),t}\_{\tau} g(-\tau,\omega) \\ &+ \int\_{-\tau}^{0} \mathbf{e}^{a(\omega)(t-s)} \mathbf{e}^{b\_{1}(\omega),t-\tau-s}\_{\tau} (g'(s,\omega) - a(\omega)g(s,\omega)) \, \mathrm{d}s \\ &+ \int\_{0}^{t} \mathbf{e}^{a(\omega)(t-s)} \mathbf{e}^{b\_{1}(\omega),t-\tau-s}\_{\tau} f(s,\omega) \, \mathrm{d}s, \end{split} \tag{2}$$

where *b*<sup>1</sup> = e−*aτb* and

$$\mathbf{e}\_{\tau}^{\varsigma,t} = \begin{cases} 0, & -\infty < t < -\tau, \\ 1, & -\tau \le t < 0, \\ 1 + c\frac{t}{1!}, & 0 \le t < \tau, \\ 1 + c\frac{t}{1!} + c^2 \frac{(t-\tau)^2}{2!}, & \tau \le t < 2\tau, \\ \vdots & \vdots \\ \sum\_{k=0}^{n} c^k \frac{(t - (k-1)\tau)^k}{k!}, & (n-1)\tau \le t < n\tau, \end{cases}$$

is the delayed exponential function [22] (Definition 1), *<sup>c</sup>*, *<sup>t</sup>* <sup>∈</sup> <sup>R</sup>, and *<sup>n</sup>* <sup>=</sup> *<sup>t</sup>*/*<sup>τ</sup>* <sup>+</sup> 1 (here · denotes the integer part defined by the so-called floor function). This formal solution is obtained with the method of steps and the method of variation of constants.

The primary objective of this paper is to set probabilistic conditions under which *x*(*t*) is an L*p*-solution to (1). We decompose the original problem (1) as

$$\begin{cases} \ y'(t,\omega) &= \left. a(\omega)y(t,\omega) + b(\omega)y(t-\tau,\omega) \right\}, \ t \ge 0, \\\ y(t,\omega) &= \left. g(t,\omega) \right\vert\_{\tau} - \tau \le t \le 0, \end{cases} \tag{3}$$

and

$$\begin{cases} \ z'(t,\omega) &= a(\omega)z(t,\omega) + b(\omega)z(t-\tau,\omega) + f(t,\omega), \; t \ge 0, \\\ z(t,\omega) &= 0, \; -\tau \le t \le 0. \end{cases} \tag{4}$$

System (3) does not possess a stochastic forcing term, and it was deeply studied in the recent contribution [17]. Under certain assumptions, its L*p*-solution is expressed as

$$\begin{split} y(t,\omega) &= \mathbf{e}^{a(\omega)(t+\tau)} \mathbf{e}^{b\_1(\omega),t}\_{\tau} g(-\tau,\omega) \\ &\quad + \int\_{-\tau}^{0} \mathbf{e}^{a(\omega)(t-s)} \mathbf{e}^{b\_1(\omega),t-\tau-s}\_{\tau} (\mathbf{y}'(s,\omega) - a(\omega)\mathbf{y}(s,\omega)) \, \mathrm{d}s,\tag{5} \end{split} \tag{5}$$

as a generalization of the deterministic solution to (3) obtained via the method of steps [22] (Theorem 1). Problem (4) is new and requires an analysis in the L*p*-sense, in order to solve the initial problem (1). Its formal solution is given by

$$z(t,\omega) = \int\_0^t \mathbf{e}^{a(\omega)(t-s)} \mathbf{e}^{b\_1(\omega), t-\tau-s} f(s,\omega) \, \mathrm{d}s,\tag{6}$$

see [22] (Theorem 2). In order to differentiate (6) in the L*p*-sense, one requires the extension of the deterministic Leibniz's integral rule for differentiation to the random scenario. This extension is an important piece of this paper.

In Section 2, we show preliminary results on L*p*-calculus that are used through the exposition, which correspond to those preliminary results from [17] and the new random Leibniz's rule for L*p*-Riemann integration. Auxiliary but novel results to demonstrate the random Leibniz's integral rule are Fubini's theorem and a chain rule theorem. In Section 3, we prove in detail that *x*(*t*) defined by (2) is the unique L*p*-solution to (1), by analyzing problem (4). We also find closed-form expressions for some statistics (expectation and variance) of *x*(*t*) related to its moments. Section 4 deals with the L*p*-convergence of *x*(*t*) as the delay *τ* tends to 0. We then show a numerical example that illustrates the theoretical findings. Finally, Section 5 draws the main conclusions.

In order to complete a fair overview of the existing literature, it must be pointed out that, apart for random delay differential equations (which is the context of this paper), other complementary approaches are stochastic delay differential equations and fuzzy delay differential equations. Stochastic delay differential equations are those in which uncertainty appears due to stochastic processes with irregular sample-paths: the Brownian motion process, Wiener process, Poisson process, etc. A new tool is required to tackle equations of this type, called Itô calculus [23]. Studies on stochastic delay differential equations can be read in [24–28], for example. On the other hand, in fuzzy delay differential equations, uncertainty is driven by fuzzy processes; see [29] for instance. In any of these approaches, the delay might even be considered random; see [30,31].
