Dynamical large deviations of reflected diffusions

We consider an ergodic reflected diffusion $X_{t}$, as defined in the previous section, and a dynamical observable $S_{T}$ of this process having the form (1). We are interested in finding the rate function of $S_{T}$, defined by the limit

where $P_{T}(s)$ is the probability distribution of $S_{T}$. In practice, it is very difficult to obtain this distribution exactly, which is why the rate function is sought instead. Indeed, for a large class of Markov processes and observables, it can be shown that $P_{T}(s)\approx e^{-TI(s)}$ with sub-exponential corrections in $T$, so that $I(s)$ effectively describes the shape of $P_{T}(s)$ at leading order in $T$ Touchette (2009). This holds in the limit of large integration times $T$, typically much longer than the relaxation time-scale of $X_{t}$.

To calculate the rate function, we use the fact that it is dual to another large deviation function, called the scaled cumulant generating function (SCGF) Dembo and Zeitouni (1998) and defined as

Using the Feynman–Kac formula, it can be shown Touchette (2018) that this function coincides with the dominant (real) eigenvalue of a linear operator, called the tilted generator, having the form

where $\mathcal{L}$ is the generator of $X_{t}$ and $f$ is the function appearing in the definition (1) of the observable $S_{T}$. Having the SCGF, we then obtain $I(s)$ by taking a Legendre transform, so that

where $k(s)$ is the unique root of $\lambda^{\prime}(k)=s$. This essentially holds if $\lambda(k)$ is differentiable and strictly convex.

The calculation of the rate function thus reduces to solving the spectral problem

where $r_{k}$ is the eigenfunction of $\mathcal{L}_{k}$ corresponding to the dominant eigenvalue $\lambda(k)$. Given that $\mathcal{L}_{k}$ is not in general Hermitian, this spectral problem must be solved in conjunction with its dual

where $l_{k}$ is the eigenfunction of $\mathcal{L}_{k}^{\dagger}=\mathcal{L}^{\dagger}+kf$ with the same dominant eigenvalue $\lambda(k)$. The boundary conditions that we must impose on these two spectral equations to obtain $\lambda(k)$ are discussed next. Independently of these conditions, $r_{k}$ and $l_{k}$ are dual functions with respect to the standard inner product (15), so they must satisfy $\langle l_{k},r_{k}\rangle<\infty$. They are also positive functions, since they are the dominant modes of $\mathcal{L}_{k}$ and $\mathcal{L}_{k}^{\dagger}$, respectively. In the literature Chetrite and Touchette (2015a), it is common to normalize them in such a way that

and

The latter integral is consistent, as we will see, with the fact that $\mathcal{L}_{k}^{\dagger}$ is related to the Fokker–Planck operator acting on normalized densities.

The boundary conditions that must be imposed on $r_{k}$ and $l_{k}$ to solve the spectral equations (22) and (23) are determined by considering, as done before, the boundary term that arises when transforming $\mathcal{L}_{k}$ to its dual $\mathcal{L}_{k}^{\dagger}$. Since these operators differ from $\mathcal{L}$ and $\mathcal{L}^{\dagger}$, respectively, only by the constant term $kf$, which produces no boundary terms of its own, the result of (16) can readily be used to obtain

where

is the probability current associated with $l_{k}$.

For diffusions evolving in $\mathbb{R}$ without boundaries, we see that the boundary term in (26) vanishes for any $r_{k}(x)$ and $l_{k}(x)$ if the latter eigenfunction decays to 0 sufficiently fast as $|x|\rightarrow\infty$. This is consistent with the normalization integral in (25) and implies with $l_{k}>0$ that both $l^{\prime}_{k}(x)$ and $J_{F,l_{k}}(x)$ vanish as $|x|\rightarrow\infty$. There is no condition on $r_{k}$ alone, as such, since the Markov generator $\mathcal{L}$, and by extension $\mathcal{L}_{k}$, have no natural boundary conditions, except for the fact that both act on functions that have finite inner product, which for $\mathcal{L}_{k}$ translates into the integral condition in (24), normalized to 1.

For reflected diffusions, the normalization integrals (24) and (25) are no longer sufficient on their own to define boundary conditions on $r_{k}$ and $l_{k}$. Instead, we must impose

and

in order for the boundary term in (26) to vanish and, importantly, for $\mathcal{L}_{k}^{\dagger}$ and $\mathcal{L}_{k}$ to have boundary conditions that are consistent with those imposed on $\mathcal{L}^{\dagger}=\mathcal{L}^{\dagger}_{k=0}$ and $\mathcal{L}=\mathcal{L}_{k=0}$, respectively, when $k=0$. The same conditions also follow by noticing again that the boundary term in the duality of $\mathcal{L}_{k}^{\dagger}$ and $\mathcal{L}_{k}$ does not depend on $kf$, so that (28) simply extends the zero-current condition of $\mathcal{L}^{\dagger}$ to $\mathcal{L}_{k}^{\dagger}$, while (29) extends the Neumann boundary condition of $\mathcal{L}$ to $\mathcal{L}_{k}$. As a result, $\mathcal{D}(\mathcal{L}_{k}^{\dagger})=\mathcal{D}(\mathcal{L}^{\dagger})$ and $\mathcal{D}(\mathcal{L}_{k})=\mathcal{D}(\mathcal{L})$. These boundary conditions must be imposed, incidentally, not just on the eigenfunctions related to the dominant eigenvalue, but to all eigenfunctions, thereby determining the full spectrum of $\mathcal{L}_{k}$ which is conjugate to the spectrum of $\mathcal{L}_{k}^{\dagger}$. For large deviation calculations, however, we only need the dominant eigenvalue, which is real, and the corresponding eigenfunctions, which are positive.

While the SCGF and the rate function describe the fluctuations of $S_{T}$, they do not provide any information about how these fluctuations are created in time. Recently, it has been shown Chetrite and Touchette (2013, 2015a, 2015b) that this information is provided by a modified Markov process $\hat{X}_{t}$, called the effective or driven process, which represents in some approximate way the original diffusion $X_{t}$ conditioned on the fluctuation $S_{T}=s$ and so describes the paths of $X_{t}$ that lead to or create that fluctuation.

The details of this process can be found in many studies Chetrite and Touchette (2013, 2015a, 2015b), which provide a full derivation of its properties and interpretation. In the context of ergodic diffusions, the driven process satisfies the new SDE

where the driven force $F_{k}$ is a modification of the force $F$ given in terms of the dominant eigenfunction $r_{k}$ by

The parameter $k$ is the same as that entering in the SCGF: its value is set for a given fluctuation $S_{T}=s$ according to the Legendre transform (21) as the root of $\lambda^{\prime}(k)=s$ or, equivalently, as $k=I^{\prime}(s)$ ³³3These two relations, which follow from the Legendre transform (21), hold when $I(a)$ is strictly convex Chetrite and Touchette (2015a).. With this choice, the stationary density of the driven process, which is known to be given (Chetrite and Touchette, 2015a, Sec. 5.4) by

is such that

This shows that we can also interpret the driven process as a change of process that transforms the fluctuation $S_{T}=s$ into a typical value realized by $\hat{X}_{t}$ in the ergodic limit. The change of drift and density also modifies the stationary current to

Note that for $k=0$, $r_{k=0}=1$ while $l_{k=0}=\rho^{*}$, leading to $F_{k=0}=F$, $\rho_{k=0}^{*}=\rho^{*}$ and $J_{F_{k=0},\rho^{*}_{k=0}}=J_{F,\rho^{*}}$.

For an ergodic diffusion $X_{t}$ evolving on $[a,b]$ with reflecting boundaries, the driven process $\hat{X}_{t}$ also evolves on $[a,b]$, since it represents a conditioning of $X_{t}$, and inherits for the same reason a zero-current boundary conditions at $x=a^{+}$ and $x=b^{-}$, given in terms of the driven current by

This can be verified, in fact, independently of the interpretation of $\hat{X}_{t}$ by applying the boundary conditions on $r_{k}$ and $l_{k}$ discussed previously to the driven process. First, note that the Neumann boundary conditions (29) on $r_{k}$ implies with the definition of the driven force (31) that the latter is constrained to satisfy

Hence, the drift at the boundaries is not modified at the level of the driven process, although it is modified in the interior of $[a,b]$, as we will see in the next section with specific examples. This is a significant result of our work, which also implies that any density $\rho$ satisfying a zero-current condition at the boundaries with respect to the original drift must also satisfy a zero-current condition at the boundaries with respect to the driven force. In other words, $J_{F,\rho}(a^{+})=0$ is equivalent to $J_{F_{k},\rho}(a^{+})=0$ and similarly for $x=b^{-}$.

Next, we note that the boundary conditions (28) and (29) can be combined as

and

We recognize these as the boundary terms in (26), which can be combined to obtain

Upon expanding the expression of the current associated with $l_{k}$ and combining terms, we then obtain

which is a zero-current condition for the product $r_{k}l_{k}$. From this result, we finally recover (35) using (36) and the fact that $r_{k}l_{k}$ is the stationary density of the driven process, as shown in (32).

This confirms in a more direct way that the driven process also evolves in $[a,b]$ with perfect reflections at the boundaries. Of course, since we are dealing with one-dimensional diffusions, the stationary current of the driven process must vanish not just at the boundaries but over the whole of $[a,b]$, similarly to the original diffusion. This is a known result in the theory of the driven process: if the original Markov process is reversible, then so is the driven process in the case where the dynamical observable has the form of (1) Chetrite and Touchette (2015a). This does not mean that $\rho^{*}$ and $\rho^{*}_{k}$ are the same in $[a,b]$ – we will see illustrations of this point in the next section. Moreover, note that requiring $J_{F_{k},\rho_{k}}=0$ at the boundaries is not sufficient to define boundary conditions for $r_{k}$ and $l_{k}$ since the driven current mixes both $F_{k}$ and $\rho_{k}$. These conditions are determined again by studying the boundary term arising in the duality between $\mathcal{L}_{k}$ and $\mathcal{L}_{k}^{\dagger}$.

The spectral calculation of the SCGF is complicated by the fact that the tilted generator $\mathcal{L}_{k}$ is not Hermitian in general. A significant simplification is possible when the spectrum of $\mathcal{L}_{k}$ is real, as is the case, for example, when $X_{t}$ is an ergodic, gradient diffusion characterized by the Gibbs stationary distribution (8) and $S_{T}$ is defined as in (1). Then it is known Touchette (2018) that $\mathcal{L}_{k}$ can be transformed in a unitary way, therefore preserving its spectrum, to the following Hermitian operator:

which involves a quantum-like potential given by

where $U$ is the potential of the gradient force and $f$ the function defining the observable $S_{T}$.

Since the new operator $\mathcal{H}_{k}$ has the same spectrum as $\mathcal{L}_{k}$, the spectral problem associated with the SCGF becomes

where $\psi_{k}$ is the corresponding eigenfunction related to $r_{k}$ and $l_{k}$ by

From (24), $\psi_{k}$ thus satisfies the normalization condition

similarly to that found in quantum mechanics, but for a real eigenfunction. Moreover, using $\psi_{k}^{2}=r_{k}l_{k}=\rho^{*}_{k}$ we find with (40) that $\psi_{k}$ must satisfy for reflected diffusions the zero-current boundary conditions

This can be expressed more explicitly as

with a similar expression holding for $x=b^{-}$.

The first example that we consider is the reflected Ornstein–Uhlenbeck process (ROUP) satisfying the SDE

with $\gamma>0$, $\sigma>0$, $X_{t}\in[0,\infty)$, and perfect reflection at $x=0$. This diffusion is used in engineering as a continuous-space model of queues in the heavy-traffic regime Giorno et al. (1986); Ward and Glynn (2003); Linetsky (2005) and represents, more physically, the dynamics of a Brownian particle attracted to a reflecting wall by a linear force induced, for example, by laser tweezers Ashkin (1997) or an ac trap Cohen and Moerner (2005). For this process, we consider the observable

which for laser tweezers is related to the mechanical power expended by the lasers on the particle van Zon and Cohen (2003).

The tilted generator $\mathcal{L}_{k}$ associated with this process and observable is given from (12) and (20) by

and is clearly non-Hermitian. However, because the ROUP is gradient and ergodic on $[0,\infty)$, we can apply the symmetrization transform described before by substituting the potential $U(x)=\gamma x^{2}/2$ in (42) to obtain from (41)

This defines with (43) the spectral problem that we need to solve in order to obtain the SCGF. The boundary condition (47) reduces in this case to the simple Neumann condition

For $x=\infty$, the normalization (45) requires that $\psi_{k}(x)^{2}\rightarrow 0$ as $x\rightarrow\infty$ and, therefore, $\psi_{k}(x)\rightarrow 0$ as $x\rightarrow\infty$.

The spectral problem (43) defines for $\mathcal{H}_{k}$ above a differential equation whose solutions are taken from the class of parabolic cylinder functions $D_{\nu}(z)$, the dominant solution $\psi_{k}$ having the form

up to a normalization constant, where we have defined

This solution decays to 0 at infinity. By imposing the boundary condition (52) and using well-known properties of the parabolic cylinder functions, we find that $\lambda(k)$ is determined implicitly by the transcendental equation

To be more precise, $\lambda(k)$ is the largest root of this equation; the other roots, which are all real, give the rest of the spectrum of $\mathcal{H}_{k}$ and therefore of $\mathcal{L}_{k}$.

We show in Fig. 2 the plot of $\lambda(k)$ obtained by solving the transcendental equation numerically for a given set of parameters $\gamma$ and $\sigma$ and different values of $k$. We can see that $\lambda(0)=0$, as follows from the definition of the SCGF, and that $\lambda(k)$ appears to be differentiable and strictly convex. As a result, we can take the Legendre transform (21) to obtain the rate function $I(s)$, shown in Fig. 3. There we see that $I(s)$ has two very different branches on either side of the minimum and zero $s^{*}$, corresponding to the most probable value of $S_{T}$ at which $P_{T}(s)$ concentrates exponentially as $T\rightarrow\infty$. The left branch, related to the $k<0$ branch of the SCGF, is steep and therefore indicates that small fluctuations of $S_{T}$ close to $0$, produced by paths that stay close to the reflecting boundary, are very unlikely. The right branch, on the other hand, is less steep and has overall the shape of a parabola, signalling that the large fluctuations of $S_{T}$ are Gaussian-distributed.

This is confirmed by comparing $I(s)$ with the rate function of $S_{T}$ for the normal Ornstein–Uhlenbeck process (OUP) evolving on the whole of $\mathbb{R}$, which is known to be

This rate function is shown with the dashed line in Fig. 3 and closely approximates, as can be seen, the upper tail of the rate function obtained with reflection. This is explained by noting that large fluctuations of $S_{T}$ arise from paths that venture far from the reflecting boundary, and so do not “feel” its influence. The zero at $s^{*}$ is itself a product of the boundary, since $s^{*}=0$ in the normal OUP. We can determine its value by noting from the ergodic theorem that the most probably value of $S_{T}$ is the stationary expectation

Here $\rho^{*}(x)$ is a truncated Gaussian density with potential $U(x)=\gamma x^{2}/2$ restricted to $x\in[0,\infty)$, leading in the integral above to $s^{*}=\sigma/\sqrt{\pi\gamma}$, which gives $s^{*}=1/\sqrt{\pi}\approx 0.564$ for the parameters used in Fig. 3.

The large Gaussian fluctuations of $S_{T}$ are also confirmed by analyzing the driven force $F_{k}$, which we find from (31) using the expression (53) for $\psi_{k}$ and the relation (44), yielding

where

To obtain this expression, we have also used some identities of the parabolic cylinder functions du Buisson (2020). The result is plotted in Fig. 4(a) for different values of $k$, together with the properly normalized $\rho_{k}^{*}(x)=\psi_{k}(x)^{2}$ in Fig. 4(b).

Comparing the two plots, we see that $F_{k}(x)$ has two zeros for $k>0$ and, therefore, two critical points: one at $x=0$, which gives rise to a local minimum in $\rho_{k}^{*}(x)$, and another at some value $x>0$, which is responsible for the maximum of $\rho_{k}^{*}(x)$. Around the latter critical point, $F_{k}(x)$ is approximately linear with slope $-\gamma$, implying that $\rho_{k}^{*}(x)$ is approximately Gaussian around its maximum. Moreover, as $k$ is increased, we see that the maximum of $\rho_{k}^{*}(x)$ moves away from $x=0$, showing that the driven process is repelled from the boundary, thus creating larger typical values of $S_{T}$, similarly to the driven process of the normal OUP Chetrite and Touchette (2015a). The difference between the two processes is that the driven force of the ROUP is “bent” near the boundary so as to have $F_{k}(0)=F(0)$, consistently with (36), whereas that of the OUP is always linear Chetrite and Touchette (2015a).

For $k<0$, the picture is different. The driven force $F_{k}(x)$ only has a single critical point at $x=0$, creating the maximum of $\rho_{k}^{*}(x)$ seen in Fig. 4(b). As $k\rightarrow-\infty$, $\rho_{k}^{*}(x)$ gets more concentrated at $x=0$, as the driven process is attracted to the reflecting boundaries, creating smaller fluctuations of $S_{T}$. Such a behavior is very unlikely in the ROUP, which is why the rate function is steep close to $s=0$. In fact, it is steeper than a parabola because $F_{k}(x)$ is not linear away from $x=0$: its curvature becomes more pronounced for large negative values of $k$, leading $\rho_{k}^{*}(x)$ to be non-Gaussian. Note in all cases that $F_{k}(0)=F(0)=0$, consistently again with (36).

This analysis of $F_{k}(x)$ is useful not only in understanding the different stochastic mechanisms underlying or “producing” different fluctuations of $S_{T}$, but also in deriving accurate approximations of the rate function by following three steps Chetrite and Touchette (2015b). First, approximate $F_{k}$ by some function, say $\tilde{F}_{\theta}$, where $\theta$ denotes a set of parameters. Second, calculate the stationary distribution $\tilde{\rho}^{*}_{\theta}$ that results from this approximation. Third and finally, calculate the ergodic expectation of $S_{T}$ with respect to $\tilde{\rho}^{*}_{\theta}$:

as well as the integral

Then

with equality if and only if $\tilde{F}_{\theta}=F_{k}$ at $s_{\theta}$ Chetrite and Touchette (2015b).

It is beyond this paper to prove this result (see Chetrite and Touchette (2015b)). At this point we only want to use it to find useful upper bounds on $I(s)$, beginning with the left branch of this function, related as mentioned before to the $k<0$ branch of $\lambda(k)$. In this case, we know that $F_{k}(x)$ has a single critical point at $x=0$, so we approximate it in a simple way as

Only $\theta\geq\gamma$ need be considered, since it is clear from Fig. 4(b) that $F^{\prime}_{k}(0)\leq-\gamma$ for $k\leq 0$. This approximation retains the linear form of $F(x)$, which means that $\tilde{\rho}^{*}_{\theta}(x)$ is the same truncated Gaussian density as the ROUP but with $\gamma$ replaced by $\theta$. As a result, we have $s_{\theta}=\sigma/\sqrt{\pi\theta}$ and obtain

by direct integration of (61). Changing the $\theta$ variable to $s$ with $s_{\theta}=s$, we then find

as our approximation of $I(s)$ for $s\in(0,s^{*}]$.

We do not compare this result with the exact rate function obtained from the spectral calculation, as the two are nearly indistinguishable du Buisson (2020). They agree exactly at $s^{*}$, since $\theta=\gamma$ recovers the drift of the ROUP, and start to differ only close to $s=0$ because $F_{k}$ is curved there, as we noticed, whereas $\tilde{F}_{\theta}$ is not. Note that the divergence of $I(s)$ near $s=0$ predicted by the approximation above is $\tilde{I}(s)\sim\sigma^{2}/(4\pi s^{2})$ as $s\rightarrow 0$, which is independent of $\gamma$.

Similar calculations can be carried out for $k>0$ to approximate $I(s)$ for $s>s^{*}$. In that case, the form of $F_{k}$ suggests that we use

as an approximate drift parameterized by $\theta>0$. The associated stationary density $\tilde{\rho}_{\theta}^{*}$ can also be obtained in closed form and yields, after solving a number of Gaussian integrals du Buisson (2020),

and

The presence of the error function in $s_{\theta}$ prevents us from expressing $K_{\theta}$ in closed form as a function of $s$, as in (65). However, the result clearly shows that the rate function becomes a parabola as $\theta\rightarrow\infty$, for in that limit, $s_{\theta}\sim\theta/\gamma$, leading to $\tilde{I}(s)=K_{\theta(s)}=I_{\textrm{OUP}}(s)$ as $s\rightarrow\infty$.

The second example we consider is the reflected Brownian motion with drift (RBMD) governed by the SDE

where $\mu>0$, $\sigma>0$, and $X_{t}\in[0,\infty)$, with reflection imposed at $x=0$ Linetsky (2005). This process was studied by Fatalov Fatalov (2017) and models, similarly to the ROUP, the dynamics of a particle attracted by a force to a reflecting wall. The force is now constant and can be viewed, for instance, as the gravity pulling vertically on a Brownian particle in a container. As for the ROUP, we consider the linear observable defined in (49).

The tilted generator associated with the large deviations of $S_{T}$ for this process is given by (12) with $F(x)=-\mu$ and leads, after symmetrization with the corresponding potential $U(x)=\mu x$, to the Hermitian operator

This defines with (43) the spectral problem that gives the SCGF, which needs to be solved with the Robin boundary condition (47)

at $x=0$ and $\psi_{k}(x)\rightarrow 0$ as $x\rightarrow\infty$.

The eigenfunctions of $\mathcal{H}_{k}$ satisfying these conditions are now expressed in terms of Airy functions of the first kind. The dominant eigenfunction is

while $\lambda(k)$ is given by the largest root $\lambda$ of the following transcendental equation:

As before, we can solve this equation numerically for given values of $\mu$ and $\sigma$ as well as various values of $k$ so as to build an interpolation of $\lambda(k)$, from which we obtain the rate function by computing the Legendre transform (21). The resulting functions are shown in Figs. 5. Unlike the ROUP, $\lambda(k)$ is now defined only for $k\leq 0$ because the potential

which is related to the classic quantum triangular well, is confining only for $k<0$, and so has bound states only for this range of parameters. This implies that the Legendre transform of $\lambda(k)$ gives $I(s)$ only for $s\in(0,s^{*}]$, as shown in Fig. 5, where $s^{*}$ is again the typical value of $S_{T}$, found here from (8) and (57) to be $s^{*}=\sigma^{2}/(2\mu)$.

Above this value, $S_{T}$ does have fluctuations, but its large deviations are not covered by the spectral calculation, which is a sign generally that $P_{T}(s)$ scales weaker than exponentially in $T$. An example of such a scaling was discussed recently for the OUP Nickelsen and Touchette (2018) using path integral techniques that predict a stretched exponential scaling in $T$, although the exact rate function cannot be found. The application of these techniques is beyond the scope of this paper, so we leave the study of the fluctuation region $S_{T}>s^{*}$ as an open problem ⁴⁴4This region is not considered by Fatalov Fatalov (2017)..

Note, incidentally, that for $\mu=0$, $S_{T}$ has no large deviations at the scale $T$ because Brownian motion is non-ergodic. The confining force produced by the negative drift along with the reflecting boundary makes that motion ergodic and creates fluctuations $S_{T}<s^{*}$ that are exponentially unlikely with $T$, though it is not strong enough, somehow, to constrain the fluctuations $S_{T}>s^{*}$ in the same exponential way. A similar effect is seen for the Brownian motion with reset den Hollander et al. (2019).

To understand how the fluctuations of $S_{T}$ arise below its mean $s^{*}$, where $P_{T}(s)$ scales exponentially, let us analyze the driven process ⁵⁵5Defining the driven process for $S_{T}>s^{*}$ is also an open problem, related again to the fact that the large deviations in that region are not given by the spectral elements of $\mathcal{L}_{k}$ Nickelsen and Touchette (2018).. The explicit expression of $F_{k}(x)$ for the RBMD is too long to show here, so we provide only the formula

which follows from (31) and (44), and which leads with (72) to a ratio of the derivative of the Airy function and the Airy function. The result, plotted in Fig. 6 for various values of $k$, is similar to the driven force found for the ROUP when $k<0$. Its shape shows that small fluctuations of $S_{T}$ are created, as expected, by squeezing the process close to the reflecting boundary by two forces: the negative constant drift $-\mu$ of the original RBMD and an added nonlinear force, which can be approximated near $x=0$ as a linear force with varying “friction” coefficient $\gamma_{k}=-F^{\prime}_{k}(0^{+})$. The action of these two forces creates a stationary density $\rho_{k}^{*}$ (not shown) which is approximately a shifted Gaussian density truncated at the boundary and increasingly concentrated there as $k\rightarrow-\infty$ or, equivalently, $S_{T}\rightarrow 0$ du Buisson (2020).