Viscoelastic active diffusion governed by nonequilibrium fractional Langevin equations: underdamped dynamics and ergodicity breaking

Sungmin Joo Department of Physics, Pohang University of Science and Technology (POSTECH), Pohang 37673, Republic of Korea Jae-Hyung Jeon jeonjh@gmail.com Department of Physics, Pohang University of Science and Technology (POSTECH), Pohang 37673, Republic of Korea Asia Pacific Center for Theoretical Physics (APCTP), Pohang 37673, Republic of Korea

Abstract

In this work, we investigate the active dynamics and ergodicity breaking of a nonequilibrium fractional Langevin equation (FLE) with a power-law memory kernel of the form $K(t)\sim t^{-(2-2H)}$ , where $1/2<H<1$ represents the Hurst exponent. The system is subjected to two distinct noises: a thermal noise satisfying the fluctuation-dissipation theorem and an active noise characterized by an active Ornstein-Uhlenbeck process with a propulsion memory time $\tau_{\mathrm{A}}$ . We provide analytic solutions for the underdamped active fractional Langevin equation, performing both analytical and computational investigations of dynamic observables such as velocity autocorrelation, the two-time position correlation, ensemble- and time-averaged mean-squared displacements (MSDs), and ergodicity-breaking parameters. Our results reveal that the interplay between the active noise and long-time viscoelastic memory effect leads to unusual and complex nonequilibrium dynamics in the active FLE systems. Furthermore, the active FLE displays a new type of discrepancy between ensemble- and time-averaged observables. The active component of the system exhibits ultraweak ergodicity breaking where both ensemble- and time-averaged MSDs have the same functional form with unequal amplitudes. However, the combined dynamics of the active and thermal components of the active FLE system are eventually ergodic in the infinite-time limit. Intriguingly, the system has a long-standing ergodicity-breaking state before recovering the ergodicity. This apparent ergodicity-breaking state becomes exceptionally long-lived as $H\to 1$ , making it difficult to observe ergodicity within practical measurement times. Our findings provide insight into related problems, such as the transport dynamics for self-propelled particles in crowded or polymeric media.

Active diffusion, self-propelled particle, viscoelastic material, Fractional Langevin equation

I Introduction

Active particles can be widely found in nature, appearing in various examples such as molecular motors-macromolecule complexes within living cell Brangwynne et al. (2009); Chen et al. (2015), self-propelled Janus colloids Palacci et al. (2010); Rao et al. (2018), and biological swimmers Lauga and Goldstein (2012); Chamolly et al. (2017). Active particles exhibit self-propelled motion with directional and/or rotational memory, leading to dynamic behaviors that are far from equilibrium. Its nonequilibrium motions are modeled by several stochastic models, such as run-and-tumble model Patteson et al. (2015); Bechinger et al. (2016), active Brownian particle Romanczuk et al. (2012); Shaebani et al. (2020), active Langevin particle Löwen (2020), and active Ornstein-Uhlenbeck particle (AOUP) Fodor et al. (2016); Bonilla (2019); Nguyen et al. (2021a). The understanding the active diffusion has broader implications for active viscoelastic systems, where the nonequilibrium motion and viscoelastic effect of crowded environments coexist. Examples include the transport of passive tracer in active baths Wu and Libchaber (2000a) or in living cells Gal and Weihs (2010); Caspi et al. (2000), chromosomal dynamics Bronstein et al. (2009); Newman et al. (2022), lateral diffusion of membrane proteins Jeon et al. (2011); Lizana et al. (2016), and the tracer diffusion in dense colloidal systems Narinder et al. (2019); Ando and Skolnick (2010). Figure 1 illustrates a schematic representation of several active viscoelastic systems. Figure 1(a) describes a typical model of intracellular environments, i.e., nonequilibrium biopolymer networks. Active particles (green sphere) diffuse in polymeric environments (gray sphere), such as acto-myosin network Wong et al. (2004); Pollard and Cooper (2009), endoplasmic reticulum network Lin et al. (2014); Speckner et al. (2018), and microtubules Vale and Hotani (1988); Caspi et al. (2000), or ATP-consuming macromolecules strongly bound to polymer strands, such as DNA Sanchez et al. (2012) and chromosomes Bronstein et al. (2009); Bronshtein et al. (2015). The presence of embedded polymers leads to viscoelastic feedback to active particles within the network. In Fig. 1(b), thermal tracers (gray spheres) are placed in an active bath. Dense active particles (semi-red spheres) provide athermal noises to the tracers as well as lead to viscoelastic feedbacks Granek et al. (2022); Lau et al. (2003); Goychuk and Hänggi (2002).

To describe the diffusion dynamics in these active complex systems, fractional Langevin equations (FLEs) have emerged as a valuable framework. These equations take the form:

\int_{0}^{t}K(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}=\zeta_{\mathrm{th}}(t)+\zeta_{\mathrm{ac}}(t).

(1)

Here, the memory kernel $K(t)$ accounts for a viscoelastic relaxation of a given embedding medium. Among a variety of viscoelastic systems, we here consider the systems that, at the timescales of interest where the many-body collective dynamics emerge, $K(t)$ is described by a power-law function

K(t)\sim t^{-\alpha}

(2)

where $0<\alpha<1$ . For example, a single monomer in a flexible polymer is governed by $K(t)\sim t^{-1/2}$ Lizana et al. (2010); Saito and Sakaue (2015); Vandebroek and Vanderzande (2017a); Joo et al. (2020) and $K(t)\sim t^{-2\nu/(1+2\nu)}$ if the excluded volume interactions are present ( $\nu$ : Flory exponent) De Gennes (1976); Panja (2010); Kappler et al. (2019). In a semi-flexible filament, the memory kernel behaves as $K(t)\sim t^{-3/4}$ Granek (1997); Caspi et al. (2000); Han et al. (2023). Moreover, $K(t)\sim t^{-2/3}$ for the height undulation dynamics of a planar membrane Granek (1997); Sung (2018) as well as a single monomer dynamics in the Zimm model Panja (2010); Sung (2018). In the presence of polymer condensation and entanglement, the memory kernel is given by $K(t)\sim t^{-1/4}$ for a single monomer in a reptating polymer Grest (2016). The $\zeta_{\mathrm{th}}(t)$ is a thermal noise whose auto-covariance is related to the memory kernel via the fluctuation-dissipation theorem (FDT), i.e., $K(t-t^{\prime})=k_{B}\mathcal{T}\langle\zeta_{\mathrm{th}}(t)\zeta_{\mathrm{th}}(t^{\prime})\rangle$ . The $\zeta_{\mathrm{ac}}(t)$ represents an active noise that violates FDT and usually has a finite memory time. It is responsible for the self-propelled movement of the active particle Romanczuk et al. (2012); Buttinoni et al. (2013); Shaebani et al. (2020) or mimics a nonequilibrium ambient noise from the active bath Wu and Libchaber (2000a); Santra (2023); Shaebani et al. (2020). The special case of $\zeta_{\mathrm{ac}}(t)=0$ is referred to as the equilibrium FLE. This equilibrium model has been employed to study the diffusion dynamics in thermal viscoelastic media Qian (2000); Chung and Hassanabadi (2017); Bonfanti et al. (2020); Höfling and Franosch (2013) induced by macromolecular crowding and/or hydrodynamics interactions Kupferman (2004); Kupferman et al. (2002); Lee et al. (2019); Fodor et al. (2015). The resulting diffusion dynamics are typically viscoelastic subdiffusion where the mean-squared displacement (MSD) scales as

\langle x^{2}(t)\rangle\sim t^{\mu}~{}\hbox{with}~{}\mu=\alpha.

(3)

Examples are the tracer diffusion in crowded dextran matrices Szymanski and Weiss (2009), micellar solutions Szymański et al. (2006); Jeon et al. (2013), and lipid membranes Flenner et al. (2009); Jeon et al. (2012); Javanainen et al. (2013).

For the nonequilibrium case of the FLE (1) with $\zeta_{\mathrm{ac}}(t)\neq 0$ , previous studies reported the following properties. When the active noise has an exponentially decaying correlation having $\langle\zeta_{\mathrm{ac}}(t)\zeta_{\mathrm{ac}}(0)\rangle\sim\exp(-t/\tau_{\mathrm{A}})$ , the viscoelastic active diffusion exhibits $\langle x^{2}(t)\rangle\sim t^{2\alpha}$ for $t<\tau_{\mathrm{A}}$ and $\langle x^{2}(t)\rangle\sim t^{2\alpha-1}$ for $t>\tau_{\mathrm{A}}$ Porrà et al. (1996); Joo et al. (2020); Han et al. (2023); Caspi et al. (2000); Lü and Bao (2005); Vandebroek and Vanderzande (2017b); Grimm and Dolgushev (2018); Sakaue and Saito (2017). For another type of active noises having a long-range correlation of a power-law form $\langle\zeta_{\mathrm{ac}}(t)\zeta_{\mathrm{ac}}(0)\rangle\sim t^{-\beta}$ , the overdamped long-time dynamics follow the MSD of $\langle x^{2}(t)\rangle\sim t^{2\alpha-\beta}$ Muralidhar et al. (1991); Porrà et al. (1996); Caspi et al. (2002); Fa (2006); Despósito and Viñales (2008); Dechant et al. (2014). Compared to the equilibrium viscoelastic subdiffusion with $\alpha=\beta$ , the active viscoelastic diffusion can exhibit more diverse dynamic patterns from superdiffusion, Fickian, subdiffusion, ultraslow diffusion ( $\sim\ln t$ ) to confined motion.

Refer to caption — Figure 1: A cartoon for active viscoelastic systems. (a) Active particles (green spheres) in a dense polymeric network represented by the gray spheres. The viscoelastic feedback originates from the interaction between the active tracers and background networks. (b) Passive tracer particles immersed in an active medium. The semi-red spheres depict self-propelled Janus particles. The passive tracer particles experience random active forces due to collisions with the self-propelled particles. These active particles collectively create a densely packed viscoelastic environment surrounding the tracers.

Extending the FLE (1), here we investigate the active viscoelastic diffusion governed by an underdamped nonequilibrium FLE:

m\ddot{x}(t)+\int_{0}^{t}K(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}=\zeta_{\mathrm{th}}(t)+\zeta_{\mathrm{ac}}(t).

(4)

In the absence of $\zeta_{\mathrm{ac}}(t)$ , the above FLE is the equilibrium underdamped FLE Porrà et al. (1996); Kursawe et al. (2013); Cherstvy et al. (2021); Bao (2017); Burov and Barkai (2008). This equilibrium model undergoes ballistic dynamics $\langle x^{2}(t)\rangle\sim t^{2}$ regardless of $\alpha$ (when $v_{0}\neq 0$ ) in the underdamped regime and cross-overs to subdiffusion with $\langle x^{2}(t)\rangle\sim t^{\alpha}$ in the overdamped regime Fricks et al. (2009); McKinley et al. (2009); Grebenkov and Vahabi (2014). As an exceptional case, when $v_{0}=0$ , the system exhibits hyperdiffusion with $\langle x^{2}(t)\rangle\sim t^{4-\alpha}$ in the underdamped regime. It was demonstrated that micron-sized particles in an actin filament network and membrane proteins in a lipid bilayer perform the viscoelastic thermal diffusion from the ballistic to a subdiffusion Grebenkov et al. (2013); Di Cairano et al. (2021).

In the case of active FLE systems ( $\zeta_{\mathrm{ac}}(t)\neq 0$ ), little is known about its full-time dynamics and ergodic properties. Here, we analytically and numerically solve the active FLE (4) and systematically investigate the dynamic properties of this non-equilibrium process. It turns out that compared to the equilibrium counterpart the active FLE system reveals unexpectedly intricate dynamic patterns as consequences of the combined effects of the non-Markovian active noise, the underdamped inertia dynamics, and the long-time viscoelastic feedback. In contrast to the equilibrium FLE, the active FLE exhibits multi-scaling underdamped dynamics involving a combination of ballistic and hyperdiffusion ( $\sim t^{4}$ , $\sim t^{3}$ ). After the underdamped dynamics, the system undergoes one of the two distinct overdamped dynamics, depending on the relative magnitude of the momentum relaxation time and the noise memory time ( $\tau_{\mathrm{A}}$ ). If the former is larger than the latter, the overdamped dynamics is simply the active subdiffusion ( $\sim t^{2\alpha-1}$ ) followed by the thermal subdiffusion, Eq. (3), beyond a cross-over time $\tau^{*}$ . For the other case, the system has the three-scaling overdamped dynamics consisting of the active superdiffusion ( $\sim t^{2\alpha}$ ), active subdiffusion ( $\sim t^{2\alpha-1}$ ), and finally thermal subdiffusion ( $\sim t^{\alpha}$ ).

In this study, additionally, we are interested in time-averaged observables and the ergodicity of the active FLE system. To our best knowledge, our study is the first analytic work to evaluate the time-averaged MSD and the ergodicity-breaking parameter. It turns out that the time-averaged MSDs significantly differ from the ensemble-averaged MSDs. In the underdamped regime, the activity-induced hyperdiffusion is absent. Instead, it simply manifests the ballistic dynamics. In the overdamped regime, the active dynamics exhibit ultraweak ergodicity breaking in which both ensemble- and time-averaged MSDs have the same scaling function but differ in their amplitudes. We also find that the activity-induced ergodicity-breaking phenomena are apparently long-standing before the active FLE system finally enters the ergodic regime for $t\gg\tau^{*}$ . The cross-over time $\tau^{*}$ has a non-monotonic dependence on the memory kernel exponent $\alpha$ . Importantly, $\tau^{*}$ rapidly increases as $\alpha$ approaches zero. Within practical observation times, accordingly, the active FLE system cannot reach the ergodic state, residing in the apparent ergodicity-breaking state.

The paper is organized in the following. In Sec. II, we provide a comprehensive explanation of our active FLE model, including the exact solutions for velocity and position variables. Additionally, we encapsulate the simulation method to solve the active FLE. In Sec. III, we present the main results of our work. This encompasses analytic and simulation studies of velocity autocorrelation, ensemble- and time-averaged MSDs, as well as an examination of the ergodicity-breaking parameters. Lastly, Sec. IV summarizes our primary findings while discussing the practical application of our theories in diverse active viscoelastic systems.

II The Model & Method

II.1 The active fractional Langevin equation model

Within the framework of the generalized Langevin equation (4), we introduce a class of nonequilibrium dynamic models referred to as the active fractional Langevin equation (FLE). This model incorporates a viscoelastic memory effect with a power-law kernel $K(t)$ in the following form:

m\ddot{x}(t)+\gamma_{0}\int_{0}^{t}(t-t^{\prime})^{2H-2}\dot{x}(t^{\prime})dt^{\prime}=\eta\xi_{H}(t)+\eta_{\mathrm{A}}\xi_{\mathrm{A}}(t).

(5)

Here, the Hurst exponent $H$ is restricted to the range of $1/2<H<1$ such that the memory term can be expressed as

\gamma_{0}\int_{0}^{t}(t-t^{\prime})^{2H-2}\dot{x}(t^{\prime})dt^{\prime}=\gamma_{0}\Gamma(2H-1)\frac{d^{2-2H}}{dt^{2-2H}}x(t).

(6)

In this expression, $\frac{d^{2-2H}}{dt^{2-2H}}f(t)=\Gamma(2H-1)^{-1}\int_{0}^{t}(t-t^{\prime})^{2H-2}\frac{d}{dt^{\prime}}[f(t^{\prime})]dt^{\prime}$ represents the Caputo fractional derivative of order $2-2H$ (with $0<2-2H<1$ ) de Oliveira and Tenreiro Machado (2014). The $\gamma_{0}(>0)$ is a generalized frictional coefficient with the physical dimension of [ $\mathrm{kg}\cdot\mathrm{s}^{-2H}$ ]. In RHS of Eq. (5), the first term describes the thermal equilibrium noise, which is connected to the memory kernel via FDT, thus,

\eta^{2}\langle\xi_{H}(t)\xi_{H}(t^{\prime})\rangle=\gamma_{0}k_{B}\mathcal{T}|t-t^{\prime}|^{2H-2}.

(7)

The thermal noise, which follows the power-law auto-covariance given by Eq. (7), is realized by fractional Gaussian noise (fGn) Mandelbrot and Ness (1968). fGn is a stationary incremental Gaussian process with zero mean $\langle\xi_{H}(t)\rangle=0$ and auto-covariance $\langle\xi_{H}(t)\xi_{H}(t^{\prime})\rangle=2H(2H-1)|t-t^{\prime}|^{2H-2}$ with $H\neq 1/2$ . The coupling constant $\eta$ is then determined to be $\eta=\sqrt{\frac{k_{B}\mathcal{T}\gamma_{0}}{2H(2H-1)}}$ . The FLE (5) without the active noise $\xi_{\mathrm{A}}$ is commonly referred to as the equilibrium FLE, which describes viscoelastic subdiffusion under equilibrium conditions Lutz (2001); Burov and Barkai (2008); Jeon and Metzler (2010); Kursawe et al. (2013); Grebenkov et al. (2013); McKinley et al. (2009); Zhu et al. (2023); Szymanski and Weiss (2009).

The $\xi_{\mathrm{A}}$ with the amplitude $\eta_{\mathrm{A}}=\sqrt{k_{B}\mathcal{T}\gamma_{0}}$ in Eq. (5) represents the active noise responsible for the self-propulsion motion of an active particle or a nonequilibrium noise from an active heat bath. Following Refs. Wu and Libchaber (2000a); Bonilla (2019); Joo et al. (2020); Kim et al. (2022); Sarkar et al. (2023), we model the active noise with a symmetric Gaussian process with a finite memory time described by an exponentially decaying auto-covariance

\langle\xi_{\mathrm{A}}(t)\xi_{\mathrm{A}}(t^{\prime})\rangle=D_{\mathrm{A}}e^{-\frac{|t-t^{\prime}|}{\tau_{\mathrm{A}}}}.

(8)

Here, $\tau_{\mathrm{A}}$ is the propulsion memory time and $D_{\mathrm{A}}$ is the (dimensionless) noise strength. Such Gaussian active noises are known as active Ornstein-Uhlenbeck (OU) noise, which is governed by the following Langevin equation Uhlenbeck and Ornstein (1930)

\frac{d\xi_{\mathrm{A}}(t)}{dt}=-\frac{1}{\tau_{\mathrm{A}}}\xi_{\mathrm{A}}(t)+\sqrt{\frac{2D_{\mathrm{A}}}{\tau_{\mathrm{A}}}}\varkappa(t)

(9)

where $\varkappa(t)$ is a delta-correlated Gaussian noise of $\langle\varkappa(t)\rangle=0$ and $\langle\varkappa(t)\varkappa(t^{\prime})\rangle=\delta(t-t^{\prime})$ . The active particle driven by the active OU noise is referred to as active Ornstein-Uhlenbeck particle (AOUP) Fodor et al. (2016); Bonilla (2019); Nguyen et al. (2021a). In this sense, the active FLE (5) describes the nonequilibrium diffusion dynamics of AOUPs embedded in a viscoelastic environment.

It is worthwhile to introduce two limiting models of the active FLE (5). When $H\to 1/2$ , it becomes the ordinary underdamped AOUP model Nguyen et al. (2021a)

m\ddot{x}(t)+\gamma_{0}\dot{x}(t)=\eta\xi_{H}(t)+\eta_{\mathrm{A}}\xi_{\mathrm{A}}(t).

(10)

This model explains the diffusion dynamics of AOUPs in a viscous fluid Fodor et al. (2016); Bonilla (2019). When $H\to 1$ , the memory kernel term reduces to $\gamma_{0}[x(t)-x(0)]$ , so the corresponding Langevin equation is responsible for AOUPs in an elastic material Szamel (2023). For $H\in(1/2,~{}1)$ , the system is viscoelastic.

The active FLE (5) can be analytically solved in the Laplace space ( $\tilde{f}(s)\equiv\mathcal{L}[f(t)]=\int_{0}^{\infty}e^{-st}f(t)dt$ ). Considering the initial conditions $\dot{x}(0)=v_{0}$ and $x(0)=0$ , the Laplace-transformed velocity is given by

\tilde{v}(s)=\frac{v_{0}+\eta m^{-1}\widetilde{\xi_{H}}(s)+\eta_{\mathrm{A}}m^{-1}\widetilde{\xi_{\mathrm{A}}}(s)}{s+\tau_{H}^{-2H}s^{1-2H}}

(11)

where $\tau_{H}=(m/[\gamma_{0}\Gamma(2H-1)])^{1/2H}$ is a characteristic time of the system. Using the generalized Mittag-Leffler function Bateman and Erdélyi (1955); Haubold et al. (2011)

E_{a,b}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(ak+b)}

(12)

where $a,b\in\mathbb{R}$ and $z\in\mathbb{C}$ (see Appendix B for further information), the analytic expression for the velocity $v(t)$ is obtained as follows:

$\displaystyle v(t)=$	$\displaystyle v_{0}E_{2H,1}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]$	(13)
	$\displaystyle+\frac{\eta}{m}\int_{0}^{t}\xi_{H}(t^{\prime})E_{2H,1}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}$
	$\displaystyle+\frac{\eta_{\mathrm{A}}}{m}\int_{0}^{t}\xi_{\mathrm{A}}(t^{\prime})E_{2H,1}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}.$

In the limit of $H\to 1/2$ (viscous fluids), the Mittag-Leffler function is simplified to $E_{1,1}(x)=\exp(-x)$ and Eq. (13) is the solution to the ordinary Langevin equation (10). In this case, the system’s relaxation dynamics is described by an exponential function with the momentum relaxation time $\tau_{1/2}=m/\gamma_{0}$ . For $1/2<H<1$ , the relaxation function $E_{2H,1}[-(t/\tau_{H})^{2H}]\approx\exp[-(t/\tau_{H})^{2H}/\Gamma(2H+1)]$ for $t\ll\tau_{H}$ and $\approx\frac{\tau_{H}^{2H}}{\Gamma(1-2H)}t^{-{2H}}$ for $t\gg\tau_{H}$ . Accordingly, the $\tau_{H}$ can be viewed as the momentum relaxation time for the underdamped dynamics of a viscoelastic system to relax.

The velocity, Eq. (13), can be viewed as the sum of the two contributions:

v(t)=v_{\mathrm{th}}(t)+v_{\mathrm{ac}}(t).

(14)

The former is the solution for the equilibrium FLE without the active noise Kursawe et al. (2013); Cherstvy et al. (2021); Bao (2017), i.e.,

	$\displaystyle v_{\mathrm{th}}(t)=$	$\displaystyle v_{0}E_{2H,1}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]$		(15)
		$\displaystyle+\frac{\eta}{m}\int_{0}^{t}\xi_{H}(t^{\prime})E_{2H,1}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}.$		(15)

The latter is the contribution of the active noise

v_{\mathrm{ac}}(t)=\frac{\eta_{\mathrm{A}}}{m}\int_{0}^{t}\xi_{\mathrm{A}}(t^{\prime})E_{2H,1}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}.

(16)

The position $x(t)$ can be obtained analytically using the similar Laplace transform technique. The position is also written as the sum of the thermal and active parts, i.e., $x(t)=x_{\mathrm{th}}(t)+x_{\mathrm{ac}}(t)$ where

	$\displaystyle x_{\mathrm{th}}(t)=$	$\displaystyle v_{0}tE_{2H,2}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]$		(17)
		$\displaystyle+\frac{\eta}{m}\int_{0}^{t}\xi_{H}(t^{\prime})t^{\prime}E_{2H,2}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}$		(17)

is the solution for the equilibrium FLE Kursawe et al. (2013); Cherstvy et al. (2021); Bao (2017) and

x_{\mathrm{ac}}(t)=\frac{\eta_{\mathrm{A}}}{m}\int_{0}^{t}\xi_{\mathrm{A}}(t^{\prime})t^{\prime}E_{2H,2}\left[-\left(\frac{t-t^{\prime}}{\tau_{H}}\right)^{2H}\right]dt^{\prime}.

(18)

II.2 Numerical method

To complement our analytic study, we numerically solve our FLE model for various parameter conditions and perform the Langevin dynamics simulation. The numerical technique to solve the active FLE is explained in detail in Appendix A. Here we briefly outline the simulation scheme. The same protocol is used in Refs. Deng and Barkai (2009); Jeon and Metzler (2010); Guo et al. (2013). In simulating the active FLE (5), we set $\gamma_{0}=1$ or $100$ , $m=1$ , $k_{B}\mathcal{T}=1$ , and $H\in\{5/8,~{}3/4,~{}7/8\}$ . The propulsion time of the active OU noise was set to be $\tau_{\mathrm{A}}=0.01$ , $0.1$ , $1$ , and $10$ . The noise strength was $D_{\mathrm{A}}=0$ (thermal), $25$ , $100$ , and $400$ . In the simulation, the initial position was always fixed to $x_{0}=0$ . On the other hand, we studied the $v_{0}$ -dependence on the transport dynamics: For the case of a fixed initial velocity, it was put to be $v_{0}=0$ . For a Gaussian-distributed random initial velocity, it was chosen to have the variance $\langle v_{0}^{2}\rangle_{\mathrm{th}}=k_{B}\mathcal{T}/m$ for the FDT-satisfying equilibrium systems and $\langle v_{0}^{2}\rangle=\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ (i.e., the stationary value) for the FDT-violating systems.

We conducted simulations with a total observation time of $T=Nh$ , where the total time step was set to $N=15\times 10^{3}$ – $10^{8}$ and the integration time was $h=10^{-4}$ , $10^{-3}$ , or $10^{-2}$ . Each parameter set was simulated using a range of $500$ – $20000$ runs to ensure statistically significant ensemble-averaged results. Prior to the main numerical simulations, we rigorously validated our simulation code by applying it to solve equilibrium FLE systems reported in Ref. Kursawe et al. (2013) for various initial conditions. The simulation results were in excellent agreement with the existing analytical theory. Further details about this validation process can be found in Appendix A & Fig. A1.

III Results: Active diffusion in viscoelastic media

III.1 Velocity autocorrelation function (VACF)

The two-time VACF $\langle v(t)v(t^{\prime})\rangle$ can be formally obtained using the solution $v(t)$ [Eq. (13)]. Based on the fact that the thermal and active noises are independent of each other, we can write the VACF as the sum of the thermal and active VACFs: $\langle v(t)v(t^{\prime})\rangle=\langle v_{\mathrm{th}}(t)v_{\mathrm{th}}(t^{\prime})\rangle+\langle v_{\mathrm{ac}}(t)v_{\mathrm{ac}}(t^{\prime})\rangle$ . The thermal VACF is evaluated in accordance with Ref. Bao (2017) as:

\langle v_{\mathrm{th}}(t)v_{\mathrm{th}}(t^{\prime})\rangle=\langle v^{2}\rangle_{\mathrm{th}}F_{1}(\tau)+\left(\langle v_{0}^{2}\rangle-\langle v^{2}\rangle_{\mathrm{th}}\right)F_{1}(t)F_{1}(t^{\prime}),

(19)

where $\tau=t^{\prime}-t$ is the time lag between two observation times ( $t<t^{\prime}$ ), $\langle v^{2}\rangle_{\mathrm{th}}=k_{B}\mathcal{T}/m$ , and $F_{1}(t)=E_{2H,1}\left[-(t/\tau_{H})^{2H}\right]$ . The VACF for the thermal velocity $v_{\mathrm{th}}(t)$ satisfies the stationary condition at thermal equilibrium. The product term of $F_{1}(t)$ and $F_{1}(t^{\prime})$ is ignored at the equilibrium initial condition of $\langle v_{0}^{2}\rangle=\langle v^{2}\rangle_{\mathrm{th}}$ or in the infinite-time limit of $t,t^{\prime}\to\infty$ . From Eq. (16), the active part of the VACF can be expressed as

\langle v_{\mathrm{ac}}(t)v_{\mathrm{ac}}(t^{\prime})\rangle=\frac{\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}}{m^{2}}\int_{0}^{t}\{F_{1}(x+\tau)F_{2}(x)+F_{1}(x)F_{2}(x+\tau)\}dx

(20)

where $F_{2}(t)=\mathcal{L}^{-1}[\tau_{\mathrm{A}}(1+s\tau_{\mathrm{A}})^{-1}(s+\tau_{H}^{-2H}s^{1-2H})^{-1}](t)$ . While $v_{\mathrm{ac}}(t)$ does not reach the equilibrium state, its VACF can satisfy the time-translation symmetry $\langle v_{\mathrm{ac}}(t)v_{\mathrm{ac}}(t^{\prime})\rangle=\langle v_{\mathrm{ac}}(\tau)v_{\mathrm{ac}}(0)\rangle$ in the limit of $t,t^{\prime}\to\infty$ due to the stationary nature of the active OU noise. Thus, we can introduce the stationary VACF as follows:

C_{v}(t)\equiv\lim_{t_{0}\to\infty}\langle v(t_{0}+t)v(t_{0})\rangle=C_{v,\mathrm{th}}(t)+C_{v,\mathrm{ac}}(t)

(21)

where

C_{v,\mathrm{th}}(t)=\langle v^{2}\rangle_{\mathrm{th}}E_{2H,1}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]

(22)

and

C_{v,\mathrm{ac}}(t)=\frac{\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}}{m^{2}}\int_{0}^{\infty}\{F_{1}(t^{\prime}+t)F_{2}(t^{\prime})dt^{\prime}+F_{1}(t^{\prime})F_{2}(t^{\prime}+t)\}dt^{\prime}.

(23)

In these expressions, $t$ acts as time lag. Now let us examine the behaviors of $C_{v}(t)$ upon varying the active noise.

First, the stationary mean-squared velocity, i.e., the kinetic energy of the active FLE system, is obtained by evaluating $C_{v}(t\to 0)$ , yielding

\langle\overline{v^{2}}\rangle_{\mathrm{st}}=\langle v^{2}\rangle_{\mathrm{th}}\left(1+2\frac{\gamma_{0}D_{\mathrm{A}}}{m}\int_{0}^{\infty}F_{1}(t^{\prime})F_{2}(t^{\prime})dt^{\prime}\right).

(24)

Here, the second term explains the contribution of the active component of the velocity. This term dominates over the thermal counterpart when the active noise is sufficiently stronger than the thermal noise. The mathematical expression suggests that the mean-squared velocity monotonically increases with the strength of the active noise $D_{\mathrm{A}}$ . On the other hand, it has a non-monotonic dependence on $\tau_{\mathrm{A}}$ . Figure 2(a) shows $\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ as a function of $\tau_{\mathrm{A}}$ where the simulation result (symbol) is plotted together with the theoretical expression (24) (solid line). Notably, $\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ has a uni-modal profile with the maximum occurring at $\tau_{\mathrm{A}}\simeq\tau_{H}~{}(\approx 0.03)$ . This indicates that the energy injection into the system through the active OU noise reaches its maximal efficiency when $\tau_{\mathrm{A}}\approx\tau_{H}$ . This behavior can be understood as follows. In the limit of $\tau_{\mathrm{A}}\approx 0$ , the active noise changes too rapidly compared to the viscoelastic relaxation, so the active noise hardly injects energy into the system; in this limit, $\langle\overline{v^{2}}\rangle_{\mathrm{st}}=\langle v^{2}\rangle_{\mathrm{th}}\left(1+\frac{\gamma_{0}D_{\mathrm{A}}\tau_{\mathrm{A}}}{m}\int_{0}^{\infty}F_{1}(t^{\prime})^{2}dt^{\prime}\right)\approx\langle v^{2}\rangle_{\mathrm{th}}(1+\mathcal{O}(\tau_{\mathrm{A}}))$ . In the regime where $0\ll\tau_{\mathrm{A}}<\tau_{H}$ , a longer propulsion memory time allows the system to move coherently with the viscoelastic directional movement, resulting in an increase in kinetic energy. However, in the opposite regime where $\tau_{\mathrm{A}}>\tau_{H}$ , the active noise hampers the system’s directional motion, leading to the dissipation of the active energy into the environment. The energy dissipation becomes more significant as $\tau_{\mathrm{A}}$ increases. In the limit of $\tau_{\mathrm{A}}\to\infty$ (i.e., the constant-force limit), it is found that $\int_{0}^{\infty}F_{1}(t^{\prime})F_{2}(t^{\prime})dt^{\prime}=F_{2}(t)^{2}|^{\infty}_{0}=0$ , and thus $\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ converges to $\langle v^{2}\rangle_{\mathrm{th}}$ .

In Fig. 2(b) we plot the profile of $C_{v,\mathrm{th}}(t)$ for several values of $H$ (symbols: simulation, lines: Eq. (22)). The thermal component of the VACF generally shows a negative correlation. As the system has a more elastic component ( $H\to 1$ ), the negative correlation becomes stronger, and the thermal VACF exhibits oscillatory behaviors. Our numerical analysis reveals that the oscillation starts at $H\approx 0.7111$ , and these oscillations become more pronounced as $H$ approaches 1.

In Fig. 2(c) we examine the behaviors of the active VACF, $C_{v,\mathrm{ac}}(t)$ , for increasing $\tau_{\mathrm{A}}$ at fixed $H$ . It exhibits distinct characteristics depending on $\tau_{\mathrm{A}}$ . When $\tau_{\mathrm{A}}\ll\tau_{H}$ , $C_{v,\mathrm{ac}}(t)$ has a very similar profile with $C_{v,\mathrm{th}}(t)$ . For example, compare the active VACF for $\tau_{\mathrm{A}}=0.01$ with $C_{v,\mathrm{th}}(t)$ (plotted in the dashed line). Both curves are almost the same over the entire time window. Such agreement is consistently observed for other values of $H$ (not shown). Given the fact that $C_{v,\mathrm{th}}(t)\sim F_{1}(t)$ is the viscoelastic relaxation function of the medium (see Eqs. (15) and (16)), driving the system with a delta-correlated-like athermal noise is expected to give dynamic information on the medium’s inherent viscoelastic response ¹¹1Rigorously, in the limit of $\tau_{\mathrm{A}}\to 0$ for a fixed $D_{\mathrm{A}}$ , $C_{\mathrm{v,ac}}(t)\approx 2A\int_{0}^{\infty}F_{1}(x+t)F_{1}(x)dx$ while $C_{v,\mathrm{th}}(t)\sim F_{1}(t)$ . Although $C_{\mathrm{v,ac}}(t)$ does not precisely converges to $C_{v,\mathrm{th}}(t)$ as $\tau_{\mathrm{A}}\to 0$ , two functions are numerically very similar to each other.. Evidently, $C_{v,\mathrm{ac}}(t)$ deviates from the thermal VACF with increasing $\tau_{\mathrm{A}}$ . It has a pronounced positive correlation regime for large values of $\tau_{\mathrm{A}}$ , which lasts roughly up to $t\approx\tau_{\mathrm{A}}$ . In this regime, the system is driven by a nearly constant active force, causing its velocity to be in the same direction as the force. In the other regime of $t\gg\tau_{\mathrm{A}}$ , the active VACF eventually becomes negative due to the effect of the viscoelastic medium, decaying as $C_{v,\mathrm{ac}}(t)\sim t^{1-4H}$ irrespective of $\tau_{\mathrm{A}}$ .

The profile of $C_{v}(t)$ is given by the sum of the thermal and active VACFs. Figure 2(d) shows the simulated $C_{v}(t)$ with the theoretical expression (21) for $D_{\mathrm{A}}=100$ and $400$ . For comparison, we also plot the theoretical lines for $C_{v,\mathrm{th}}(t)$ (dashed, black) and $C_{v,\mathrm{ac}}(t)$ (dashed, red). When the active noise is weak, it captures the essential features of $C_{v,\mathrm{th}}(t)$ ; it exhibits a negative correlation and oscillatory behaviors for large values of $H$ . As the active noise is stronger, $C_{v}(t)$ tends to follow $C_{v,\mathrm{ac}}(t)$ . Instead of the negative correlation, it has a positive correlation for $t<\tau_{H}$ and a weak negative correlation for $t>\tau_{H}$ .

III.2 Mean-squared displacements (MSDs)

Following the same line of procedure in the previous section, we first evaluate the two-time position correlator $\langle x(t)x(t^{\prime})\rangle$ . According to Eqs. (17) and (18), $\langle x(t)x(t^{\prime})\rangle$ is the combination of both the thermal and active contributions, which can be expressed as follows:

	$\displaystyle\langle x_{\mathrm{th}}(t)x_{\mathrm{th}}(t^{\prime})\rangle=$	$\displaystyle(\langle v_{0}^{2}\rangle-\langle v^{2}\rangle_{\mathrm{th}})G_{1}(t)G_{1}(t^{\prime})$		(25)
		$\displaystyle+\langle v^{2}\rangle_{\mathrm{th}}\left\{G_{2}(t)+G_{2}(t^{\prime})-G_{2}(\tau)\right\}$		(25)

and

	$\displaystyle\langle x_{\mathrm{ac}}(t)x_{\mathrm{ac}}(t^{\prime})\rangle=\frac{\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}}{m^{2}}\int_{0}^{t}$	$\displaystyle\left\{G_{1}(t^{\prime}+\tau)G_{3}(t^{\prime})\right.$		(26)
		$\displaystyle+\left.G_{1}(t^{\prime})G_{3}(t^{\prime}+\tau)\right\}dt^{\prime}$		(26)

where $G_{1}$ , $G_{2}$ , and $G_{3}$ are the Mittag-Leffler variation functions defined as:


$\displaystyle G_{1}(t)$	$\displaystyle=tE_{2H,2}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right],$	(27a)
$\displaystyle G_{2}(t)$	$\displaystyle=\int_{0}^{t}G_{1}(t^{\prime})dt^{\prime}=t^{2}E_{2H,3}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right],$	(27b)
$\displaystyle G_{3}(t)$	$\displaystyle=\mathcal{L}^{-1}\left[\frac{\tau_{\mathrm{A}}}{(s\tau_{\mathrm{A}}+1)(s^{2}+\tau_{H}^{-2H}s^{2-2H})}\right](t).$	(27c)

With the given system parameters, $\langle v^{2}\rangle_{\mathrm{th}}$ is set to be unity. Starting from these analytically derived correlator expressions, we proceed to formulate the general equations for the ensemble- and time-averaged MSDs and analyze their inherent characteristics.

III.2.1 The ensemble-averaged MSD

From the above expressions, the ensemble-averaged MSD is obtained as

\langle x^{2}(t)\rangle=\langle x_{\mathrm{th}}^{2}(t)\rangle+\langle x_{\mathrm{ac}}^{2}(t)\rangle

(28)

where the thermal part of MSD reads

\langle x_{\mathrm{th}}^{2}(t)\rangle=2\langle v^{2}\rangle_{\mathrm{th}}G_{2}(t)+\left(\langle v_{0}^{2}\rangle-\langle v^{2}\rangle_{\mathrm{th}}\right)G_{1}(t)^{2}

(29)

and the active part is expressed as

\langle x_{\mathrm{ac}}^{2}(t)\rangle=2\frac{\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}}{m^{2}}\int_{0}^{t}G_{1}(t^{\prime})G_{3}(t^{\prime})dt^{\prime}

(30)

where $G_{i}~{}(i=1,2,3)$ are the aforementioned Mittag-Leffler functions shown in Eq. (27). At the thermal (equilibrium) initial condition, the thermal MSD (29) increases as $\propto G_{2}(t)$ . Thus, the thermal MSD grows with the ballistic diffusion exponent $\mu=2$ at short times ( $t\ll\tau_{H}$ ) for all $H$ values and then cross-overs to the subdiffusion regime of $\sim t^{2-2H}$ in the overdamped limit ( $t\gg\tau_{H}$ ). As a special case, if the system starts with $v_{0}=0$ , the system initially performs a hyperdiffusion as $\sim t^{2+2H}$ instead of the ballistic dynamics Deng and Barkai (2009); Cherstvy et al. (2021); Kursawe et al. (2013); Siegle et al. (2010a, b). Refer to the plots of $\langle x_{\mathrm{th}}^{2}(t)\rangle$ in Fig. A1(a) (Appendix A).

The active part of MSD, $\langle x_{\mathrm{ac}}^{2}(t)\rangle$ , has three distinct scaling regimes over time. Plugging $G_{3}(t)\approx G_{2}(t)$ for $t\ll\tau_{\mathrm{A}}$ and $\approx\tau_{\mathrm{A}}G_{1}(t)$ for $t\gg\tau_{\mathrm{A}}$ into Eq. (30), we find that MSDs follow power-law scalings if $H\in(1/2,~{}3/4)$ . In Case 1 where the self-propulsion time $\tau_{\mathrm{A}}$ is greater than the momentum relaxation time $\tau_{H}$ —a typical condition for active colloidal systems—the MSD behaves as

\langle x_{\mathrm{ac}}^{2}(t)\rangle\sim\left\{\begin{array}[]{lc}t^{4},&~{}0\leq t\ll\tau_{H}\\ t^{4-4H},&~{}\tau_{H}\ll t\ll\tau_{\mathrm{A}}\\ t^{3-4H},&~{}t\gg\tau_{\mathrm{A}}\end{array}\right..

(31)

This indicates that the active displacement exhibits a hyperdiffusion of $\sim t^{4}$ in the underdamped regime and, subsequently, displays a sub-ballistic superdiffusion ( $\sim t^{4-4H}$ ) and subdiffusion ( $\sim t^{3-4H}$ ). In Case $2$ where $\tau_{H}>\tau_{\mathrm{A}}$ , the active MSD follows the scaling

\langle x_{\mathrm{ac}}^{2}(t)\rangle\sim\left\{\begin{array}[]{lc}t^{4},&~{}0\leq t\ll\tau_{\mathrm{A}}\\ t^{3},&~{}\tau_{\mathrm{A}}\ll t\ll\tau_{H}\\ t^{3-4H},&~{}t\gg\tau_{H}\end{array}\right..

(32)

In this case, the system has a complicated underdamped dynamics such that the $t^{4}$ hyperdiffusion is followed by a $t^{3}$ hyperdiffusion when $t$ is greater than $\tau_{\mathrm{A}}$ . After these underdamped dynamics, the system suffers the subdiffusion of $\sim t^{3-4H}$ in the overdamped regime. The three distinct dynamics are separated by two important time constants of the systems, i.e., $\tau_{\mathrm{A}}$ and $\tau_{H}$ . Moreover, depending on their relative magnitude, the underdamped dynamics display two distinct dynamic patterns (see the first two scalings in Eqs. (31) & (32)). We cross-check these analytic results with the simulation. Figure 3(a) presents the plots of $\langle x_{\mathrm{ac}}^{2}(t)\rangle$ with $(\tau_{H},~{}\tau_{\mathrm{A}})=(0.357,~{}0.01)$ and $(0.009,~{}1)$ . Indeed, three scaling regimes, predicted by the above theory, are clearly visible for both cases. Numerically we further find that if $\tau_{\mathrm{A}}\approx\tau_{H}$ , the second regime ( $t^{4-4H}$ in Case $1$ or $t^{3}$ in Case $2$ ) vanishes and the MSD grows from $\sim t^{4}$ to $\sim t^{3-4H}$ (not shown). Also note that if $H\geq 3/4$ the active MSD does not exhibit the power-law scaling of $\sim t^{3-4H}$ in the overdamped regime (while the underdamped power-laws remain valid), which will be explained below Eq. (34).

We have a few comments on the underdamped dynamics in the active part of MSD. (1) When $\tau_{\mathrm{A}}<\tau_{H}$ , the system has two-fold hyperballistic underdamped dynamics, $\sim t^{4}$ and $\sim t^{3}$ , separated by $\tau_{\mathrm{A}}$ (see Eq. (32)). We notice that these active hyperdiffusion can be understood in the context of the hyperdiffusion ( $\sim t^{2+2H}$ ) shown in the equilibrium FLE starting with the initial condition of $v_{0}=0$ . At the timescale of $0<t\ll\tau_{\mathrm{A}}$ , the active noise does not lose its directional correlation and acts as fGn with $H=1$ in the equilibrium FLE, so leading to the hyperdiffusion of $\sim t^{4}$ . For $\tau_{\mathrm{A}}\ll t\ll\tau_{H}$ , the active noise completely loses the directionality, which can be treated as a $\delta$ -correlated noise, i.e., fGn with $H=1/2$ . Thus, a hyperdiffusion of $\sim t^{3}$ emerges at this timescale. (2) When $\tau_{\mathrm{A}}>\tau_{H}$ , the system has a single underdamped dynamics of $\sim t^{4}$ for $t\ll\tau_{H}$ by the same reasoning explained in (1). A very interesting overdamped regime occurs for this case at the timescale of $\tau_{H}\ll t\ll\tau_{\mathrm{A}}$ . Here the system’s momentum relaxation is sufficiently reached while the active noise is persistently imposed on the system like a constant force. The resultant active dynamics turns out to be the sub-ballistic superdiffusion with the MSD of $\sim t^{4-4H}$ , capturing the viscoelastic property of the environment. If the medium is viscous ( $H=1/2$ ), the system shows ballistic dynamics due to the very persistent active force. If the medium is elastic, the active force is fully stored into the elastic energy without any kinetic energy ( $\sim t^{0}$ ). In a viscoelastic medium between these two limits, the system displays superdiffusive ( $1/2<H<3/4$ ), Fickian ( $H=3/4$ ), and subdiffusive ( $3/4<H<1$ ) dynamics depending on the Hurst exponent.

In the long-time regime where $t\gg\tau_{H}$ and $\tau_{\mathrm{A}}$ , Eq. (30) yields

\langle x_{\mathrm{ac}}^{2}(t)\rangle\approx\frac{2k_{B}\mathcal{T}D_{\mathrm{A}}\tau_{\mathrm{A}}\sin{(\pi(2-2H))}}{\gamma_{0}\pi\Gamma(2H-1)\Gamma(3-2H)}\int_{0}^{t}y^{2-4H}dy.

(33)

This shows that the long-time active dynamics are classified into the following three cases:

\langle x_{\mathrm{ac}}^{2}(t)\rangle\approx A\times\left\{\begin{array}[]{lc}\frac{1}{3-4H}t^{3-4H},&\frac{1}{2}<H<\frac{3}{4}\\ \ln{t},&H=\frac{3}{4}\\ \mathrm{const}.,&\frac{3}{4}<H<1\end{array}\right.

(34)

with $A=\frac{2k_{B}\mathcal{T}D_{\mathrm{A}}\tau_{\mathrm{A}}\sin{(\pi(2-2H))}}{\gamma_{0}\pi\Gamma(2H-1)\Gamma(3-2H)}$ . We note that Eq. (34) is consistent with the results reported in some previous studies on FDT-violating overdamped systems Porrà et al. (1996); Wang and Tokuyama (1999); Fa (2006). For example, an overdamped nonequilibrium FLE having a power-law memory kernel ( $\sim t^{-\alpha}$ ) conjugated with a fGn of autocovariance decaying as $\sim t^{-\beta}$ exhibits $\langle x^{2}(t)\rangle\sim t^{2\alpha-\beta}$ for $\alpha>\beta/2$ , $\ln t$ at $\alpha=\beta/2$ , and a $\mathrm{constant}$ for $\alpha<\beta/2$ Porrà et al. (1996); Wang and Tokuyama (1999); Fa (2006). In Refs. Porrà et al. (1996); Joo et al. (2020); Han et al. (2023); Caspi et al. (2000); Lü and Bao (2005); Vandebroek and Vanderzande (2017b); Grimm and Dolgushev (2018); Sakaue and Saito (2017), a nonequilibrium FLE with an FDT-violating noise having a finite correlation time (e.g., an exponential function) was investigated. The study demonstrated that $\langle x^{2}(t)\rangle\sim t^{2\alpha-1}$ for $1/2<\alpha<1$ , $\ln t$ at $\alpha=1/2$ , and $\mathrm{a~{}constant}$ for $0<\alpha<1/2$ .

In Fig. 3(b), we plot the total MSD, $\langle x^{2}(t)\rangle=\langle x_{\mathrm{th}}^{2}(t)\rangle+\langle x_{\mathrm{ac}}^{2}(t)\rangle$ , for varying $D_{\mathrm{A}}$ and $v_{0}$ . In the absence of the active noise ( $D_{\mathrm{A}}=0$ ), the system’s dynamics follows that of the equilibrium FLE. Here, the total MSD corresponds to $\langle x_{\mathrm{th}}^{2}(t)\rangle$ , where the short-time dynamics are either ballistic if $v_{0}\neq 0$ or hyperdiffusive ( $\sim t^{2+2H}$ ) if $v_{0}=0$ . At the overdamped timescale, they crossover to subdiffusive dynamics ( $\sim t^{2-2H}$ ). When a strong noise is introduced ( $D_{\mathrm{A}}=400$ ), the scaling properties of $\langle x^{2}(t)\rangle$ mostly follow those of $\langle x_{\mathrm{ac}}^{2}(t)\rangle$ . However, the short-time regime is altered. If $v_{0}\neq 0$ , the ballistic scaling ( $\sim t^{2}$ ) from the thermal MSD dominates over the hyperdiffusion of the active component ( $\sim t^{4}$ ). Accordingly, the full MSD initially increases with $\sim t^{2}$ , then has a transient accelerating dynamics ( $\sim t^{4}$ ), and enters the scaling regimes of the superdiffusion $\sim t^{3/2}$ and the active subdiffusion $\sim t^{1/2}$ as time progresses. This corresponds to the scenario described by Eq. (31) with $\tau_{\mathrm{A}}>\tau_{H}$ . When $v_{0}=0$ , the MSD grows with $\sim t^{2+2H}$ , accompanied with a short cross-over ( $\sim t^{4}$ ), and exhibits the regimes of $\sim t^{3/2}$ and $\sim t^{1/2}$ with increasing time. We note that the active subdiffusion persists until a cross-over time

\tau^{*}\approx O(D_{\mathrm{A}}^{1/(2H-1)}).

(35)

For example, when $H=3/4$ and $D_{\mathrm{A}}\sim O(10^{2})$ , $\tau^{*}$ is estimated to be $O(10^{6})$ with other constant factors. For $t\gg\tau^{*}$ , the thermal motion eventually dominates in EA MSD, and the MSD exhibits the thermal subdiffusion ( $\sim t^{2-2H}$ ).

III.2.2 The time-averaged MSD

Instead of the ensemble-averaged MSD above, we can alternatively define an MSD from a single trajectory $x(t)$ in the time-averaging sense:

\overline{\delta^{2}(t,T)}=\frac{1}{T-t}\int_{0}^{T-t}[x(t^{\prime}+t)-x(t^{\prime})]^{2}dt^{\prime}.

(36)

In this expression, $t$ acts as time lag and $T$ is the total measurement time. Given the fact that the active FLE (5) is a stationary process, here we define and study the time-averaged (TA) MSD $\overline{\delta^{2}(t)}=\lim_{T\to\infty}\overline{\delta(t,T)}$ for sufficiently long trajectories. Further, we define the ensemble-average of $\overline{\delta^{2}(t)}$ over many trajectories, which is expressed as

\langle\overline{\delta^{2}(t)}\rangle=\lim_{T\to\infty}\frac{2}{T-t}\int_{0}^{T-t}[\langle x^{2}(t^{\prime})\rangle-\langle x(t^{\prime}+t)x(t^{\prime})\rangle]dt^{\prime}.

(37)

In this integral, we utilize the stationary parts of $\langle x^{2}(t)\rangle$ and $\langle x(t^{\prime}+t)x(t^{\prime})\rangle$ because the effect from the $v_{0}$ -dependent terms in Eq. (37) becomes negligible if $T$ is sufficiently long.

The TA MSD is also expressed as the superposition of the equilibrium and active components, i.e, $\langle\overline{\delta^{2}(t)}\rangle=\langle\overline{\delta_{\mathrm{th}}^{2}(t)}\rangle+\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle$ . The thermal part of TA MSD is given by Deng and Barkai (2009); Jeon and Metzler (2010)

\langle\overline{\delta_{\mathrm{th}}^{2}(t)}\rangle=2\langle v^{2}\rangle_{\mathrm{th}}G_{2}(t),

(38)

which precisely corresponds to the thermal part of EA MSD (7) under the equilibrium condition. The active part of TA MSD is approximated as:

	$\displaystyle\frac{\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle}{2A^{\prime}}\approx$	$\displaystyle\int_{t}^{\infty}G_{1}(t^{\prime}-t)G_{3}(t^{\prime}-t)dt^{\prime}+\int_{0}^{\infty}G_{1}(t^{\prime})G_{3}(t^{\prime})dt^{\prime}$		(39)
		$\displaystyle-\int_{0}^{\infty}\left[G_{1}(t^{\prime}+t)G_{3}(t^{\prime})+G_{1}(t^{\prime})G_{3}(t^{\prime}+t)\right]dt^{\prime}$		(39)

where $A^{\prime}=\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}/m^{2}$ , and $G_{1}(t)$ and $G_{3}(t)$ are the Mittag-Leffler variation functions, Eqs. (27a) & (27c), respectively. By utilizing the Taylor expansion of Eq. (39), we deduce the scaling properties of the TA MSD at $t\to 0$ and $t\to\infty$ . In the former case, $\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle$ behaves as

	$\displaystyle\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle$	$\displaystyle\approx 2\frac{\eta_{\mathrm{A}}^{2}D_{\mathrm{A}}}{m^{2}}\left(\int_{0}^{\infty}F_{1}(t^{\prime})F_{2}(t^{\prime})dt^{\prime}\right)t^{2}$		(40)
		$\displaystyle=\left(\langle\overline{v^{2}}\rangle_{\mathrm{st}}-\langle v^{2}\rangle_{\mathrm{th}}\right)t^{2},$		(40)

indicating that the active part of TA MSD always exhibits a ballistic scaling, regardless of the values of $H$ . In contrast, the long-time overdamped asymptotic for $t\gg\tau_{\mathrm{A}},\tau_{H}$ is highly dependent on the values of $H$ . (1) For the range $1/2<H<3/4$ , we obtain

\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle\approx-2A\frac{\Gamma(4H-3)\Gamma(2-2H)}{\Gamma(2H-1)}t^{3-4H}.

(41)

Here $A=\frac{2k_{B}\mathcal{T}D_{\mathrm{A}}\tau_{\mathrm{A}}\sin{(\pi(2-2H))}}{\gamma_{0}\pi\Gamma(2H-1)\Gamma(3-2H)}$ is the same prefactor that also appears in Eq. (34) and $\Gamma(4H-3)<0$ . It is noted that the active part of TA MSD exhibits the same subdiffusive power-law exponent $3-4H$ as the corresponding EA MSDs (31) & (32). However, the amplitudes of the EA and TA MSDs are not identical. (2) In the case of $H=3/4$ , the TA MSD is obtained to be

\frac{\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle}{A}\approx 2\gamma_{\mathrm{E}}+2\ln{t}+\frac{\Gamma^{\prime}(2-2H)}{\Gamma(2-2H)}+\frac{\Gamma^{\prime}(2H-1)}{\Gamma(2H-1)}

(42)

where $\gamma_{\mathrm{E}}=0.5772$ is the Euler’s constant and $\Gamma^{\prime}(z)=d\Gamma(z)/dz$ . It is important to note that the active TA MSD differs from $\langle x_{\mathrm{ac}}^{2}(t)\rangle$ by a factor of $2$ in front of the $\ln t$ term and also by other extra constants. (3) For the range $3/4<H<1$ , the active TA MSD follows that of a confined diffusion, where the long-time asymptotic value of the TA MSD is twice the value of the counterpart EA MSD. The disagreement between the EA and TA MSDs is further analyzed in detail in Sec. III.3.

Figure 3(c) illustrates the simulated $\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle$ for $H=5/8$ , considering three cases: $\tau_{H}>\tau_{\mathrm{A}}$ , $\tau_{H}\approx\tau_{\mathrm{A}}$ , and $\tau_{H}<\tau_{\mathrm{A}}$ . In all three cases, the active part of TA MSD exhibits the expected two distinct power-law regimes ( $\sim t^{2}$ and $\sim t^{1/2}$ ) over time. This is in contrast with the complicated scaling behaviors observed in $\langle x_{\mathrm{ac}}^{2}(t)\rangle$ . Notably, the initial hyperdiffusion ( $\sim t^{4}$ ) observed in the EA counterpart is replaced by ballistic diffusion, and the accompanied transient scaling regimes ( $\sim t^{3}$ or $\sim t^{4-4H}$ ) vanish. Furthermore, the simulation confirms that the long-time dynamics occur for $t\gtrsim\mathrm{max}(\tau_{H},~{}\tau_{\mathrm{A}})$ .

In Fig. 3(d), we present the total TA MSD $\langle\overline{\delta^{2}(t)}\rangle$ for increasing $D_{\mathrm{A}}$ while keeping $\tau_{\mathrm{A}}$ and $\tau_{H}$ fixed. At $D_{\mathrm{A}}=0$ , the TA MSD exhibits only two scaling regimes: an initial ballistic dynamics for $t\lesssim\mathrm{max}(\tau_{H},\tau_{\mathrm{A}})$ , followed by the thermal subdiffusion characterized by a power-law of $\sim t^{2-2H}$ beyond the cross-over time. However, in the presence of the active noise ( $D_{\mathrm{A}}\neq 0$ ), the scaling behaviors become more complex. After the ballistic regime, a transient regime of superdiffusion emerges in the time window of $\tau_{H}<t<\tau_{\mathrm{A}}$ . In this plot, we observe an apparent anomalous exponent of approximately $1.7$ , which arises from the superposition of $\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle\sim t^{2}$ and $\langle\overline{\delta_{\mathrm{th}}^{2}(t)}\rangle\sim t^{2-2H}$ . After this regime, the TA MSD exhibits active subdiffusion with a scaling relation of $\sim t^{3-4H}$ (e.g., see $t^{1/2}$ at $D_{\mathrm{A}}=400$ ). For $t\gg\tau^{*}$ , the TA MSD is expected to enter the final regime of the thermal subdiffusion ( $\sim t^{2-2H}$ ).

Figure 4 compares the EA and TA MSDs for different values of the Hurst exponent: (a) $H=5/8$ , (b) $3/4$ , and (c) $7/8$ for Case $1$ ( $\tau_{\mathrm{A}}>\tau_{H}$ ). In each panel, the EA and TA MSDs are simulated with $D_{\mathrm{A}}=0$ and $400$ with an initial condition of $\langle v_{0}^{2}\rangle=1$ . As expected, the active-free system ( $D_{\mathrm{A}}=0$ ) has identical EA and TA MSDs, demonstrating the ergodicity of the equilibrium FLE process. On the contrary, the EA and TA MSDs deviate from each other for the active FLEs ( $D_{\mathrm{A}}=400$ ). These discrepancies signify that the active FLE (5) is non-ergodic within our observation time window; the discrepancies become more pronounced as $H$ increases. In Fig. A2, we present the EA and TA MSDs for Case $2$ ( $\tau_{H}>\tau_{\mathrm{A}}$ ), with the same values of $H$ . Analogously to Fig. 4, the two MSDs exhibit increasing discrepancies with larger values of $H$ . However, a difference is that the TA MSDs are always greater than the EA MSDs across all times. This observation will be further examined below in Sec. III.3.

We emphasize that the observed non-ergodic behavior arises from the combined influence of the viscoelastic memory and underdamped dynamics in the system. The conventional AOUP model, which neglects the memory term in the active FLE (5), is ergodic regardless of the active motion. In the active FLE model where the viscoelastic effect is incorporated, the active component reveals ultraweak ergodicity breaking in the sense that the EA and TA MSDs exhibit the same scaling form but differ only in their amplitudes, as Eqs. (34), (41), & (42) show above.

III.3 Ergodicity-breaking parameter

The analysis of MSDs in the active FLE model reveals that the system exhibits a non-ergodic behavior at certain timescales, and the propensity of non-ergodicity strongly depends on the elapsed time $t$ and Hurst exponent $H$ . To quantify the extent of ergodicity breaking, we introduce the ergodicity-breaking (EB) parameter, given by the expression Godec and Metzler (2013):

\mathcal{EB}(t)=\frac{\langle\overline{\delta^{2}(t)}\rangle}{\langle x^{2}(t)\rangle}.

(43)

Using the EB parameter, we begin by investigating the non-ergodic behavior of the active component dynamics in the active FLE system. In Fig. 5, we present $\mathcal{EB}_{\mathrm{ac}}(t)=\langle\overline{\delta_{\mathrm{ac}}^{2}(t)}\rangle/\langle x_{\mathrm{ac}}^{2}(t)\rangle$ for $H=5/8$ (a, b) and $H=7/8$ (c, d). For each $H$ value, we compare the case of $\tau_{\mathrm{A}}>\tau_{H}$ (Left panels) with the other one of $\tau_{H}>\tau_{\mathrm{A}}$ (Right panels). The EB parameters evidently show that the active viscoelastic dynamics are non-ergodic, and the degree of ergodicity breaking significantly changes with time. In Case $1$ (Left panels), the system exhibits a strong disparity between the two MSDs in the underdamped regime where $\mathcal{EB}_{\mathrm{ac}}(t)$ decays as $\sim t^{-2}$ for $t\lesssim\tau_{H}$ (see the inset plots therein). The EB parameter decreases below unity in the intermediate timescale of $\tau_{H}\lesssim t\lesssim\tau_{\mathrm{A}}$ . In the final regime of $t\gg\tau_{\mathrm{A}}$ , $\mathcal{EB}_{\mathrm{ac}}(t)$ reaches a stationary value, $\mathcal{EB}_{\mathrm{ac,\infty}}(\neq 1)$ and the active component displays ultraweak ergodicity-breaking (UWEB), as aforementioned. We analytically obtain $\mathcal{EB}_{\mathrm{ac,\infty}}(\neq 1)$ from the formal expressions for the two MSDs, which reads

\mathcal{EB}_{\mathrm{ac},\infty}\approx\left\{\begin{array}[]{lc}1&H=\frac{1}{2}\\ 2\frac{\Gamma(4H-2)\Gamma(2-2H)}{\Gamma(2H-1)}&\frac{1}{2}<H<\frac{3}{4}\\ 2&H=\frac{3}{4}\\ 2&\frac{3}{4}<H<1\end{array}\right..

(44)

Here, the case of $H=1/2$ refers to the active LE, i.e., the AOUP in a viscous medium. For the active viscoelastic system with $H\in(1/2,~{}1)$ , the EB parameter satisfies $\mathcal{EB}_{\mathrm{ac},\infty}\in(1,~{}2]$ and increases monotonically with $H$ (see Fig. A3). Figure 5 indeed shows that the theoretical prediction [Eq. (44), the solid line] explains successfully the simulated $\mathcal{EB}_{\mathrm{ac},\infty}$ . Surprisingly, the results suggest that the extent of ergodicity breaking only depends on the viscoelastic memory parameter $H$ , but the noise strength does not contribute to the ergodicity breaking. In Case $2$ (Right panels), the system exhibits time-dependent ergodicity breaking as well, however, of which pattern is clearly distinguished from Case $1$ . In this case, the ensemble-averaged MSD has two distinct underdamped scaling laws [see Eq. (32)], which makes the two power-law scaling of $\sim t^{-2}$ and $\sim t^{-1}$ in $\mathcal{EB}_{\mathrm{ac}}(t)$ for $t<\tau_{H}$ ; see the insets. In the overdamped regime, the EB parameter is eventually saturated to $\mathcal{EB}_{\mathrm{ac},\infty}$ predicted by Eq. (44). Intriguingly, $\mathcal{EB}_{\mathrm{ac}}$ approaches the stationary value with oscillation when $H$ is close to $1$ .

In Fig. 6, we examine the EB parameter for the total dynamics of the active FLE (5), i.e., the sum of the thermal and active components. Before investigating the active FLE system, we examine the active LE model, the so-called conventional AOUP model that neglects the viscoelastic memory Wu and Libchaber (2000b); Löwen (2020); Nguyen et al. (2021b); Martin et al. (2021); Bonilla (2019); Lisin et al. (2021), corresponding to the active FLE (5) with $H\to 1/2$ . In Fig. 6(a, upper), we evaluate $\mathcal{EB}(t)$ for this model with initial velocities taken from the stationary velocity distribution ( $\langle v_{0}^{2}\rangle=\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ ). The simulation results show that the EB parameter, starting from unity at short times, deviates from it in the timescale of $\tau_{1/2}\lesssim t\lesssim\tau_{\mathrm{A}}$ (here $\tau_{1/2}=m/\gamma_{0}$ ), and eventually recovers unity. This pattern remains the same regardless of the noise strength. These results demonstrate that the memory-free AOUP system is ergodic Cherstvy et al. (2018); Wang et al. (2020), even though it violates FDT. A noteworthy feature is that the system suffers transient ergodicity breaking on the intermediate timescale, even under the stationary initial condition. This behavior is clearly distinguished from the ordinary underdamped LE (the case at $D_{\mathrm{A}}=0$ in Fig. 6(a)), where $\mathcal{EB}(t)=1$ for all times. This transient ergodicity breaking is attributed to the fact that the AOUP model is subject to two independent noise sources of which stationary properties are distinguished. This indicates that either the thermal or active components or both of them will start initially in a non-stationary state. The thermal and active parts are then separately relaxed to their stationary states with the timescales of $\tau_{H}$ (thermal) and $\tau_{\mathrm{A}}$ (active). In the lower panel of Fig. 6(a), we plot $\mathcal{EB}(t)$ with a non-stationary initial condition $\langle v_{0}^{2}\rangle=1$ for the same AOUP model. As will be explained further below, in this case, the EB parameter starts from the short-time asymptotic value $\mathcal{EB}_{0}\approx\langle\overline{v^{2}}\rangle_{\mathrm{st}}/\langle v_{0}^{2}\rangle$ . Because $\langle\overline{v^{2}}\rangle_{\mathrm{st}}$ increases proportionally to $D_{\mathrm{A}}$ , the short-time EB, $\mathcal{EB}_{0}$ , can be arbitrarily large. As the AOUP system is ergodic, the EB parameters then monotonically decay to $\mathcal{EB}=1$ after $t\gtrsim\tau_{\mathrm{A}}$ .

Now we investigate $\mathcal{EB}(t)$ for the active viscoelastic system where the system has a power-law decaying memory kernel. Figures 6(b) and 6(c) shows the variation of the EB parameters for $H=5/8$ (b) and $7/8$ (c) with the stationary initial condition and $\tau_{\mathrm{A}}>\tau_{H}$ . For the active FLE system, importantly, $\mathcal{EB}(t)$ experiences four distinct regimes over time. The first regime is the initial-condition dominating regime in the underdamped domain ( $t\ll\tau_{H}$ , $\tau_{\mathrm{A}}$ ), where $\mathcal{EB}(t)$ has the following asymptotic behavior:

\mathcal{EB}_{0}=\left\{\begin{array}[]{ll}\langle\overline{v^{2}}\rangle_{\mathrm{st}}/\langle v_{0}^{2}\rangle,&v_{0}\neq 0\\ t^{-2H},&v_{0}=0\end{array}\right..

(45)

In Figs. 6(b) and 6(c), $\mathcal{EB}_{0}$ is unity under the stationary initial condition. Eq. (45) indicates that at a non-stationary condition $\mathcal{EB}_{0}$ can be arbitrary and its value strongly depends on the initial velocity, $D_{\mathrm{A}}$ , and $\tau_{\mathrm{A}}$ . As a special case, when $v_{0}=0$ , $\mathcal{EB}_{0}$ has no asymptotic value but decays algebraically with an exponent of $-2H$ .

After this regime, $\mathcal{EB}(t)$ enters the transient ergodicity-breaking regime in the time window of $\tau_{H}\lesssim t\lesssim\tau_{\mathrm{A}}$ , where the EB parameter significantly varies over time. Beyond this regime, the system approaches the active dynamics-dominating ergodicity-breaking regime. Here, $\mathcal{EB}(t)$ is saturated to a time-independent constant that deviates from unity. If the active noise is sufficiently stronger than the thermal, the constant is approximated to $\mathcal{EB}_{\mathrm{ac},\infty}$ explained in Eq. (44). Then, the property of $\mathcal{EB}(t)$ in this regime follows that of the active component of the active FLE studied above. This ergodicity-breaking regime persists up to $t\approx\tau^{*}$ , at which the thermal MSD becomes comparable with the active counterpart [see Eq. (35)]. For $t>\tau^{*}$ , the thermal dynamics dominates over the active part and, thus, $\mathcal{EB}(t)$ will converge to unity as $t\to\infty$ . Therefore, the active FLE (5) is, in principle, ergodic in the infinite-time limit.

It is noteworthy that despite the ergodic nature of the FLE model, the presence of viscoelasticity gives rise to a remarkably long active dynamics-dominating ergodicity-breaking state in the time window of $\tau_{\mathrm{A}}\lesssim t\lesssim\tau^{*}$ . To quantify this behavior, we plot $\tau^{*}(H)$ in Fig. 6(d) by numerically solving $\langle x_{\mathrm{th}}^{2}(\tau^{*})\rangle=\langle x_{\mathrm{ac}}^{2}(\tau^{*})\rangle$ . The numerical data (circles) are explained by two analytic approximations (depicted as dashed & dotted lines). In the regime of $1/2<H<3/4$ , the approximate expression for $\tau^{*}(H)$ can be obtained by equating Eqs. (8) and (34), yielding:

\tau^{*}(H)=\left(\frac{D_{\mathrm{A}}\tau_{\mathrm{A}}\sin(\pi(2-2H))}{\pi(3-4H)}\right)^{1/(2H-1)}~{}\hbox{for}~{}H\in(1/2,~{}3/4).

(46)

This expression is depicted by the dotted lines. In the regime of $3/4<H<1$ , the active MSD follows that of a confined diffusion [see Eq. (34)]. By approximating the active MSD as a constant $C$ independent of $H$ (see further information in the Caption), we obtain the following approximation:

\tau^{*}(H)=\left(\frac{C\Gamma(2H-1)\Gamma(3-2H)}{2\gamma_{0}}\right)^{1/(2-2H)}~{}\hbox{for}~{}H\in(3/4,~{}1).

(47)

The theoretical curve (47) is depicted by the dashed lines in the plot. Both analytic curves are shown to explain the numerical data successfully. The plot reveals a non-monotonic dependence of $\tau^{*}$ on $H$ . For $1/2<H<3/4$ , $\tau^{*}$ exhibits a local maximum around $H\sim 1/2$ , and its peak value increases significantly with increasing $H$ . Beyond this point, $\tau^{*}$ monotonically decreases and reaches a local minimum point at $H=3/4$ . In the domain of $3/4<H<1$ , $\tau^{*}$ rapidly increases with $H$ and diverges in the limit as $H\to 1$ . The same trend is observed at a higher $D_{\mathrm{A}}$ while the magnitude of $\tau^{*}$ is sensitive to $D_{\mathrm{A}}$ .

The behaviors of $\tau^{*}(H)$ suggest that observing the final ergodic regime where $\mathcal{EB}\to 1$ is nontrivial for the active FLE system. This effect becomes particularly prominent for systems with $H>3/4$ , where $\tau^{*}$ diverges as $H\to 1$ . As an illustration, Fig. 6(c) displays the plot of $\mathcal{EB}(t)$ for $H=7/8$ up to $t=10^{4}$ . Once the system enters the active dynamics-dominating ergodicity-breaking state, characterized by $\mathcal{EB}_{\mathrm{ac},\infty}\approx 2$ , the EB parameter does not decay within the observed time window. By estimating $\tau^{*}$ to be approximately $10^{8}$ for $D_{\mathrm{A}}=100$ , we can infer that this ergodicity-breaking state will persist at least up to $t\sim 10^{8}$ . Therefore, it is worth noting that while the active FLE system with $1/2<H<1$ can exhibit ergodic behavior in the infinite measurement time, the presence of viscoelastic long-time memory introduces effectively ergodicity-breaking phenomena within the finite measurement window. This indicates that the interplay between active noise and viscoelasticity plays a crucial role in shaping the dynamics and the observed ergodic properties of the system.

IV Concluding remarks

In this study, we have performed a comprehensive investigation into the under- and overdamped dynamics and ergodicity-breaking phenomena for the active FLE systems modeled by Eq. (5). A representative example within this class of models is the nonequilibrium diffusion of active Ornstein-Uhlenbeck particles (AOUPs) embedded in a viscoelastic medium. Our approach involved analytic calculations of various dynamic quantities, including the velocity and position autocorrelation functions, ensemble- and time-averaged MSDs, as well as the ergodicity-breaking parameters.

We found that the active FLE systems are decomposed into two independent dynamic components, i.e., the thermal and active parts. On the one hand, the former follows an equilibrium FLE Lutz (2001); Jeon and Metzler (2010); Kursawe et al. (2013); Grebenkov et al. (2013); McKinley et al. (2009), where the system displays ballistic dynamics in the underdamped regime ( $t\ll\tau_{H}$ ) and subdiffusion with an anomalous exponent of $\mu=2+2H$ in the overdamped regime ( $t\gg\tau_{H}$ ). On the other hand, the active component demonstrates multi-scale, complicated diffusion patterns that are strongly affected by the viscoelastic memory strength (quantified by $H$ ), as well as the values of $\tau_{H}$ and $\tau_{\mathrm{A}}$ . When $1/2<H<3/4$ , the active diffusion reveals power-law scalings, leading to two distinct MSD patterns depending on the relationship between $\tau_{\mathrm{A}}$ and $\tau_{H}$ . In Case $1$ (where $\tau_{\mathrm{A}}>\tau_{H}$ ), the active component exhibits hyperdiffusion with $\mu=4$ in the underdamped regime, followed by the multiple overdamped motion, i.e., sub-ballistic active superdiffusion with $\mu=4-4H$ for $\tau_{H}\ll t\ll\tau_{\mathrm{A}}$ and active subdiffusion with $\mu=3-4H$ for $t\gg\tau_{\mathrm{A}}$ [Fig. 3(a)]. In Case $2$ (where $\tau_{H}>\tau_{\mathrm{A}}$ ), the underdamped dynamics became more complicated. It starts with hyperdiffusion with $\mu=4$ for $t\ll\tau_{\mathrm{A}}$ , followed by another hyperdiffusion with $\mu=3$ for $\tau_{\mathrm{A}}\ll t\ll\tau_{H}$ . In the overdamped limit, the active dynamics was simply subdiffusion with $\mu=3-4H$ [Fig. 3(a)]. In the case of $H>3/4$ , the overdamped active dynamics deviate from power-law scalings. In both Cases $1$ & $2$ , the overdamped dynamics transformed into logarithmic ultraslow diffusion ( $\sim\ln t$ ) at $H=3/4$ or confined diffusion for $3/4<H<1$ . It is worth noting that the thermal dynamics dominate the active FLE systems in the limits of $t\to 0$ and $t\gg\tau^{*}$ , while the active dynamics usually prevail in between the two limits. The cross-over times were sensitively determined by the strength of active noise ( $D_{\mathrm{A}}$ ).

Notably, our analytic and numerical studies demonstrated that the active FLE system exhibits simpler dynamic characteristics when considering time-averaged observables. In the underdamped regime, the active component of TA MSD always displays ballistic diffusion for $t\ll\tau_{H}~{}(\&~{}\tau_{\mathrm{A}})$ , regardless of the initial conditions and the values of $\tau_{\mathrm{A}}$ and $\tau_{H}$ [Fig. 3(c)]. This behavior is in stark contrast to the hyperdiffusions ( $\langle x^{2}(t)\rangle\sim t^{4}$ , $t^{3}$ ) observed in the EA MSDs. The ballistic dynamics then smoothly cross-overs to the overdamped dynamics, which are characterized by subdiffusion with a scaling of $t^{3-4H}$ for $1/2<H<3/4$ , $\sim\ln t$ at $H=3/4$ , and confined diffusion for $3/4<H<1$ . Importantly, there is no sub-ballistic superdiffusion of $\sim t^{4-4H}$ observed in the overdamped counterpart of the active EA MSD when $\tau_{\mathrm{A}}>\tau_{H}$ . Furthermore, in the overdamped regime for $t\gg\tau_{H}~{}\&~{}\tau_{\mathrm{A}}$ , the TA MSD shares the same scaling behaviors with the EA counterpart. However, the amplitudes of both MSDs are not identical, indicating that the active component exhibits (ultra)weak ergodicity breaking (UWEB). We note that the observed UWEB in this work is distinguished from that reported in superdiffusive Lévy walk Froemberg and Barkai (2013); Godec and Metzler (2013). In the latter, the UWEB arises due to insufficient self-averaging in finite measurements. Similarly, in the subdiffusive continuous-time random walk (CTRW) model, the lack of time scale in the waiting times results in WEB even at infinite measurement times He et al. (2008); Albers and Radons (2013a); Barkai et al. (2012); Albers and Radons (2013b). In our active FLE system, the viscoelastic memory combined with a nonequilibrium noise plays a key role in inducing ergodicity-breaking.

We studied the ergodic properties of the active FLE systems using the ergodicity-breaking parameter, $\mathcal{EB}$ , defined by Eq. (43). Our findings revealed that long-time viscoelastic memory plays a crucial role in determining the ergodic behaviors of these systems. (1) In the absence of viscoelastic memory (corresponding to the limit $H\to 1/2$ ), the active LE–also known as the AOUP model Nguyen et al. (2021a); Bonilla (2019)–exhibits transient ergodicity breaking at intermediate timescales in between $\tau_{H}$ and $\tau_{\mathrm{A}}$ . Beyond this timescale, the viscoelasticity-free active system rapidly recovers ergodicity as time increases. (2) When the viscoelastic effect comes into play, the active FLE with $H\in(1/2,1)$ exhibits $H$ -value sensitive ergodic dynamics. Despite the violation of ergodicity observed in the active dynamics [Fig. 5], the overall dynamics of the system, which combines the nonergodic active component with the thermal ergodic component, is always ergodic in the long run in the sense that $\mathcal{EB}\to 1$ as $t$ increases beyond the cross-over time $\tau^{*}$ [Eq. (35)]. The observed ergodicity in the infinite-time regime simply reflects the ergodic property of the thermal component that dominates over the active component at large times of $t\gg\tau^{*}$ . However, before entering this ergodic regime, the active FLE system undergoes a nontrivial ergodicity-breaking state in the time window of $\tau_{\mathrm{A}}\lesssim t\lesssim\tau^{*}(H)$ . Importantly, the magnitude of $\tau^{*}(H)$ is much larger than $\tau_{\mathrm{A}}$ even at moderate strengths of $D_{\mathrm{A}}$ . As a result, the active FLE system is apparently non-ergodic within a typical finite measurement time window (as observed in Figs. 6(b) & 6(c)). Particularly in the domain of $3/4<H<1$ , it becomes almost infeasible to observe the infinite-time ergodic regime. The cross-over time $\tau^{*}(H)$ diverges as $H\to 1$ , as the analytic expressions (46) & (47) reveal. This divergence indicates that the apparent ergodicity-breaking state becomes exceptionally long-lived, making it challenging to observe ergodicity within practical measurement times.

The apparent WEB in our active FLE system is a novel phenomenon that has not been reported in other systems. A well-known WEB is the discrepancy of the scaling relations between EA MSD ( $\sim t^{\alpha}$ ) and TA MSD ( $\sim t/T^{1-\alpha}$ ) in the CTRW or heterogeneous diffusion models He et al. (2008); Albers and Radons (2013a); Barkai et al. (2012); Albers and Radons (2013b); Cherstvy et al. (2013); Massignan et al. (2014). In these models, the presence of long-lasting particle trapping events leads to non-stationary solution and breaks ergodicity. Regarding the aforementioned UWEB in, e.g., the superdiffusive Lévy walk, EA and TA MSDs have the same scaling relation but different amplitudes. In contrast, our active FLE system displays a unique combination of these features. The active component exhibits UWEB, which is combined with the ergodic thermal component in the presence of a power-law viscoelastic memory. As a result, the violation of ergodicity persists during a remarkably long time window, with the duration expected to become infinite as $H$ approaches $1$ .

The theories of the active FLE system developed in this work can be applied to quantitatively understand transport phenomena in various active viscoelastic systems. For example, it can be used to analyze the diffusion of AOUPs strongly bound (or cross-linked) to flexible polymers. In a computational investigation reported in Ref. Joo et al. (2020), it was found that, counterintuitively, the AOUPs interacting with flexible chains become more subdiffusive as their active mobility increases. Their analytic study demonstrated that the viscoelastic feedback from the flexible polymers induces an ultraslow subdiffusion with $\langle x_{\mathrm{ac}}^{2}(t)\rangle\sim\ln t$ Joo et al. (2020). This type of AOUP systems can be described by our active FLE model with $H=3/4$ corresponding to the memory kernel for flexible polymers. In agreement with Ref. Joo et al. (2020), the active component of our active FLE system has EA and TA MSDs of $\sim\ln t$ . Furthermore, our study predicts that the corresponding system exhibits significant ergodicity-breaking, characterized by $\mathcal{EB}_{\mathrm{ac}}=2$ , indicating a visible discrepancy between the two MSDs in the intermediate timescale. However, the system recovers ergodicity at timescales shorter than those for other $H$ cases (see Fig. 6(d)).

The diffusion of AOUPs connected to semiflexible filaments is evidently different from those in flexible polymers. Previous studies have reported that the overdamped diffusion dynamics of AOUPs in such a situation reveal scaling behaviors of $\langle x^{2}(t)\rangle\sim t^{3/2}$ for $t<\tau_{\mathrm{A}}$ and $\langle x^{2}(t)\rangle\sim t^{1/2}$ for $t>\tau_{\mathrm{A}}$ Han et al. (2023); Caspi et al. (2000). The active FLE with $H=5/8$ (explaining the kernel for semiflexible filaments) with Case $1$ ( $\tau_{\mathrm{A}}>\tau_{H}$ ) precisely explains these scaling behaviors [Eq. (31)]. Our study gives further information about this system. For example, in the underdamped limit, the AOUP is expected to display ballistic movement ( $\sim t^{2}$ ) as well as hyperdiffusion ( $\sim t^{4}$ ). However, if its diffusion is observed with TA MSD, it will simply display the ballistic motion in the underdamped regime and a subdiffusion of $\sim t^{1/2}$ [Fig. 3]. The underdamped hyperdiffusion and the overdamped superdiffusion of $\sim t^{3/2}$ observed in the EA MSD will be masked.

Our theoretical framework can be applicable to more complex active viscoelastic systems where multiple distinct viscoelastic memory effects are combined. For example, consider the diffusion of flexible polymers in an active viscoelastic bath, as reported in Refs. Vandebroek and Vanderzande (2015); Grimm and Dolgushev (2018). In these simulations, the center-of-mass (COM) of the polymer was observed to perform anomalous diffusion, with $\langle x_{\mathrm{CM}}^{2}(t)\rangle$ scaling as $\sim t^{4-4H}$ for $t<\tau_{\mathrm{A}}$ and $\sim t^{3-4H}$ for $t>\tau_{\mathrm{A}}$ , where $\tau_{\mathrm{A}}$ is the persistence memory time of the active Ornstein-Uhlenbeck noise of the surrounding medium. The observed active diffusion of the COM can be understood as the overdamped diffusion of an AOUP in our active FLE model with Case $1$ ( $\tau_{\mathrm{A}}>\tau_{H}$ ); see Eq. (31). Now consider the diffusion of a single monomer comprising the polymer. In this case, the dynamics of the single monomer are affected by two viscoelastic memory effects; one from the flexible polymer itself and the other from the surrounding viscoelastic medium. These two memory effects were found to induce complex diffusion dynamics of the single monomer, with $\langle x_{\mathrm{m}}^{2}(t)\rangle$ scaling as $\sim t^{\frac{3}{2}(2-2H)}$ for $t<\tau_{\mathrm{A}}$ and $\sim t^{\frac{3}{2}(2-2H)-1}$ for $t>\tau_{\mathrm{A}}$ . The observed scaling relations can be effectively explained by considering the overdamped active diffusion of Case $1$ in our active FLE model, with the memory kernel modified as $K(t)\sim t^{\frac{3}{4}(2H-2)}$ to account for the combined effects of the two viscoelastic memory components.

Our active FLE model serves as a theoretical framework for quantitatively describing activity-induced transport dynamics in viscoelastic media, e.g., synthesized hydrogels, chromosomes, cytoskeletal networks, and the endoplasmic reticulum network. Expanding the current formalism to incorporate confining or random potentials could yield insights into various trapped-and-hopping phenomena within active biological systems, offering a fundamental understanding of the active viscoelastic transport in complex environments.

Acknowledgements

This work was supported by the National Research Foundation (NRF) of Korea, Grant No. 2021R1A6A1A10042944 & No. RS-2023-00218927.

Appendix A Validation of numerical method &
FLE with thermal noise

Let us introduce a numerical algorithm to solve the active FLE (5), using the method introduced in Deng and Barkai (2009); Jeon and Metzler (2010); Guo et al. (2013). Firstly, by integrating Eq. (5) from $0$ to $t$ , we obtain the stochastic integral equation for velocity $v(t)=\dot{x}(t)$ with the thermal noise $\eta\xi_{H}(t)$ and active noise $\eta_{\mathrm{A}}\xi_{\mathrm{A}}(t)$ , given by

	$\displaystyle v(t)=$	$\displaystyle v_{0}-\frac{\gamma_{0}}{m}\int_{0}^{t}\frac{(t-t^{\prime})^{2H-1}}{2H-1}v(t^{\prime})dt^{\prime}$		(1)
		$\displaystyle+\frac{\eta}{m}\int_{0}^{t}\xi_{H}(t^{\prime})dt^{\prime}+\frac{\eta_{\mathrm{A}}}{m}\int_{0}^{t}\xi_{\mathrm{A}}(t^{\prime})dt^{\prime}.$		(1)

To numerically calculate the stochastic integral equation (1), we employ the product trapezoidal quadrature formula on discrete time points $t_{n}$ ( $n=0,~{}1,~{}\dots,~{}N$ ) within the time domain $[0,~{}T]$ Diethelm et al. (2002):

	$\displaystyle v(t_{n+1})=$	$\displaystyle v_{0}-\frac{\gamma_{0}}{m}\frac{h^{2H}}{(2H+1)!}\sum_{j=0}^{n+1}a_{j,n+1}v(t_{j})$		(2)
		$\displaystyle+\frac{\eta}{m}\sum_{j=0}^{n}\xi_{H}(t_{j})+\frac{\eta_{\mathrm{A}}}{m}\sum_{j=0}^{n}\xi_{\mathrm{A}}(t_{j}).$		(2)

The Eq. (2) provides the discrete form of the velocity, where $h(=t_{n+1}-t_{n})$ represents the time increment, and the coefficients $a_{j,n+1}$ are $n^{2H+1}-(n-2H)(n+1)^{2H}$ for $j=0$ , $(n-j+2)^{2H+1}+(n-j)^{2H-1}-2(n-j+1)^{2H+1}$ for $1\leq j\leq n$ , and $1$ for $j=n+1$ . Rearranging Eq. (2), we find the expression for the velocity at $t_{n+1}$ in the following:

$\displaystyle v(t_{n+1})=$	$\displaystyle\left(\frac{m(2H+1)!}{m(2H+1)!+h^{2H}\gamma_{0}}\right)$	(3)
	$\displaystyle\times\left(v_{0}-\frac{\gamma_{0}}{m}\frac{h^{2H}}{(2H+1)!}\sum_{j=0}^{n}a_{j,n+1}v(t_{j})\right.$
	$\displaystyle\qquad+\left.\frac{\eta}{m}\sum_{j=0}^{n}\xi_{H}(t_{j})+\frac{\eta_{\mathrm{A}}}{m}\sum_{j=0}^{n}\xi_{\mathrm{A}}(t_{j})\right).$

Once $v(t_{n+1})$ is obtained, the trajectory on $t_{n+1}$ is evaluated using the second-order Euler method

x(t_{n+1})=x(t_{n})+\frac{h}{2}(v(t_{n})+v(t_{n+1})).

(4)

In our active FLE model, we have two noise terms, i.e., $\eta\xi_{H}(t)$ and $\eta_{\mathrm{A}}\xi_{\mathrm{A}}(t)$ . The former noise is fGn explained in Sec. II, and the summation of fGn, i.e., $\sum_{j=0}^{n}\xi_{H}(t_{j})$ in Eq. (3), produces fractional Brownian motion (fBm) $x_{H}(t)\equiv\int_{0}^{t}\xi_{H}(t^{\prime})dt^{\prime}$ , which has zero mean of $\langle x_{H}(t)\rangle=0$ and the covariance $\langle x_{H}(t_{1})x_{H}(t_{2})\rangle=D_{H}\left(|t_{1}|^{2H}+|t_{2}|^{2H}-|t_{1}-t_{2}|^{2H}\right)$ . In our simulation, fBm is generated using the FFGN function implemented in Matlab Davies and Harte (1987); Beran (1998); Bardet (2002); Bardet and Bertrand (2007) for a given trajectory length $N$ , the Hurst exponent $H$ , and $D_{H}$ . In our simulation, we set $D_{H}=1/2$ . Note that the simulation of the active FLE does not depend on the value of $D_{H}$ because the prefactor $\eta\propto\sqrt{1/D_{H}}$ cancels out the $D_{H}$ -dependency. The active OU noise $\xi_{\mathrm{A}}(t_{n})$ in Eq. (3) is obtained by solving the stochastic equation (9) in a discrete version, which is

\xi_{\mathrm{A}}(t_{n+1})=\left(1-\frac{h}{\tau_{\mathrm{A}}}\right)\xi_{\mathrm{A}}(t_{n})+\sqrt{\frac{2hD_{\mathrm{A}}}{\tau_{\mathrm{A}}}}\omega

(5)

where $\omega$ is a random number chosen from a normal distribution $\mathcal{N}(0,1)$ .

We validate our code by simulating the equilibrium FLE model (neglecting the active noise part by setting $D_{\mathrm{A}}=0$ in the active FLE (5)) with various initial conditions. We measure the ensemble-averaged MSD and compare it with theoretical expectations. The MSD is defined as $\langle x^{2}(t)\rangle\equiv\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\langle v(t_{1})v(t_{2})\rangle$ , and for a stationary process we have $\langle x^{2}(t)\rangle=2\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}d\tau\langle v(\tau)v(0)\rangle$ with $\tau=|t_{1}-t_{2}|$ . The Laplace transform of the MSD follows a simple relation

\mathcal{L}[\langle x^{2}(t)\rangle](s)=\frac{2}{s^{2}}\mathcal{L}[\langle v(t)v(0)\rangle](s)=\frac{2\langle v_{0}^{2}\rangle}{s^{3}+\tau_{H}^{-2H}s^{3-2H}},

(6)

and its inverse transform leads to the the exact solution expressed in terms of a generalized Mittag-Leffler function

\langle x^{2}(t)\rangle=2\langle v^{2}\rangle_{\mathrm{th}}t^{2}E_{2H,3}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right].

(7)

By the definition of the generalized Mittag-Leffler function (12), we have $\langle x^{2}(t)\rangle\approx\langle v_{0}^{2}\rangle t^{2}$ for the short-time limit. The long-time asymptotic expression for the MSD is found to be $\langle x^{2}(t)\rangle\approx 2v_{0}^{2}t^{2-2H}/\Gamma(3-2H)$ from the series expansion (12). Thus, the MSD shows a cross-over from the ballistic motion at short times to the subdiffusive motion at long times, which has the form Deng and Barkai (2009); Jeon and Metzler (2010)

\langle x^{2}(t)\rangle\approx\frac{2k_{B}\mathcal{T}}{m}\left\{\begin{array}[]{ll}\frac{1}{\Gamma(3)}t^{2},&t\to 0\\ \frac{\tau_{H}^{2H}}{\Gamma(3-2H)}t^{2-2H},&t\to\infty.\end{array}\right.

(8)

The exact analytic expression of the MSD is reported in some previous work Kursawe et al. (2013); Bao (2017). When the initial condition is not equilibrium, i.e., $\langle v_{0}^{2}\rangle\neq k_{B}\mathcal{T}/m$ , the MSD is written as

	$\displaystyle\langle x^{2}(t)\rangle=$	$\displaystyle 2\langle v^{2}\rangle_{\mathrm{th}}t^{2}E_{2H,3}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]$		(9)
		$\displaystyle+\left(\langle v_{0}^{2}\rangle-\langle v^{2}\rangle_{\mathrm{th}}\right)t^{2}E_{2H,2}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right]^{2}$		(9)

with $\tau=|t-t^{\prime}|$ . If $\langle v_{0}^{2}\rangle=\langle v^{2}\rangle_{\mathrm{th}}$ or $t\to\infty$ , Eq. (9) yields Eq. (7). For short times $t\ll\tau_{H}$ , $\langle x^{2}(t)\rangle\sim\langle v_{0}^{2}\rangle t^{2}$ for $v_{0}\neq 0$ , while, for $v_{0}=0$ , hyperdiffusion is expected with $\langle x^{2}(t)\rangle\sim t^{2+2H}$ .

For both equilibrium and nonequilibrium initial conditions, the time-averaged MSD is obtained to be Deng and Barkai (2009); Jeon and Metzler (2010)

\langle\overline{\delta^{2}(t)}\rangle=2\langle v^{2}\rangle_{\mathrm{th}}t^{2}E_{2H,3}\left[-\left(\frac{t}{\tau_{H}}\right)^{2H}\right],

(10)

which is the same functional form as Eq. (7). In Fig. A1, we simulate the thermal FLE systems with several different initial conditions and compare them with the exact theoretical expressions above. As shown, the simulation data show excellent agreement with the theory for all the three cases.

Appendix B Generalized Mittag-Leffler functions

Here we introduce the generalized Mittag-Leffler function (MLF), $E_{a,b}(z)$ , and some important properties of this function. The generalized MLF is defined by the series expansion as follows Bateman and Erdélyi (1955); Haubold et al. (2011); Jia et al. (2019):

E_{a,b}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(ak+b)}

(11)

where $a,b\in\mathbb{C}$ , $\mathcal{R}(a)>0$ , $\mathcal{R}(b)>0$ , $z\in\mathbb{C}$ with $\mathbb{C}$ being the set of complex numbers. For large- $z$ expansion, the generalized MLF is expressed as

E_{a,b}(z)=-\sum_{k=1}^{\infty}\frac{z^{-k}}{\Gamma(b-ak)}.

(12)

When $a=1$ and $b=1$ , the MLF reduces to the exponential function $E_{1,1}(z)=e^{z}$ . Using the above series expansion formulae, we find that $E_{a,b}[-ct^{a}]\approx\exp[-ct^{a}/\Gamma(a+b)]$ for $t\ll 1$ and $\approx t^{-a}/\Gamma(b-a)$ for $t\gg 1$ .

For the Mittag-Leffler function $E_{a,b}[-ct^{a}]$ , a useful Laplace transform formula is the following:

\mathcal{L}\left[t^{b-1}E_{a,b}(\pm ct^{a})\right](s)=\frac{s^{-b}}{\left(1\mp cs^{-a}\right)}.

(13)

We employ this formula to compute the Laplace transform of the convolution forms (such as Eqs. (15), (16), (17) & (18)) that involve the generalized MLF. Another useful property of the MLF is its relation to the derivation of MLF:

\frac{d^{n}}{dt^{n}}\left[t^{b-1}E_{a,b}(ct^{a})\right]=t^{b-n-1}E_{a,b-n}(ct^{a})

(14)

where $b\geq n+1$ .

Appendix C Supplementary figures

In this section, we plot supplementary figures supporting our main results. In Fig. A2, we present the ensemble- and time-averaged MSDs for Case $2$ ( $\tau_{H}>\tau_{\mathrm{A}}$ ) with the Hurst exponent $H=5/8$ (a), $3/4$ (b), and $7/8$ (c). In three panels, the EA and TA MSDs are simulated with $D_{\mathrm{A}}=0$ and $400$ with a non-stationary initial condition of $\langle v_{0}^{2}\rangle=1$ . In the thermal FLE process ( $D_{\mathrm{A}}=0$ ), both MSDs are identical. For the active FLE process ( $D_{\mathrm{A}}=400$ ), the EA and TA MSDs deviate from each other, and the discrepancies increase as $H$ increases. Note that TA MSDs are always greater tan the EA MSDs across all times.

In Fig. A3, we plot $\mathcal{EB}_{\mathrm{ac},\infty}$ as a function of $H$ (see Eq. (44)). At $H=1/2$ , the system becomes the conventional AOUP model in viscous media, which is ergodic, thus, resulting in $\mathcal{EB}_{\mathrm{ac},\infty}=1$ . As $H$ is in the range $1/2<H<3/4$ , we observe a monotonic increase in $\mathcal{EB}_{\mathrm{ac},\infty}$ with increasing $H$ , eventually reaching a value of $2$ as $H\to 3/4$ . At $H=3/4$ , $\mathcal{EB}_{\mathrm{ac},\infty}=2$ originates from the logarithmic growths of both MSDs (Eqs. (34) & (42)). If $H>3/4$ , the confined motion of MSD maintains $\mathcal{EB}_{\mathrm{ac},\infty}=2$ (Eq. (44)).

References

Brangwynne et al. (2009) C. P. Brangwynne, G. H. Koenderink, F. C. MacKintosh, and D. A. Weitz, Trends Cell Biol. 19, 423 (2009).
Chen et al. (2015) K. Chen, B. Wang, and S. Granick, Nat. Mater. 14, 589 (2015).
Palacci et al. (2010) J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 105, 088304 (2010).
Rao et al. (2018) D. V. Rao, N. Reddy, J. Fransaer, and C. Clasen, J. Phys. D: Appl. Phys. 52, 014002 (2018).
Lauga and Goldstein (2012) E. Lauga and R. Goldstein, Phys. Today 65, 30 (2012).
Chamolly et al. (2017) A. Chamolly, T. Ishikawa, and E. Lauga, New J. Phys. 19, 115001 (2017).
Patteson et al. (2015) A. E. Patteson, A. Gopinath, M. Goulian, and P. E. Arratia, Sci. Rep. 5, 15761 (2015).
Bechinger et al. (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
Romanczuk et al. (2012) P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur. Phys. J.: Spec. Top. 202, 1 (2012).
Shaebani et al. (2020) M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Nat. Rev. Phys. 2, 181 (2020).
Löwen (2020) H. Löwen, J. Chem. Phys. 152, 040901 (2020).
Fodor et al. (2016) E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, Phys. Rev. Lett. 117, 038103 (2016).
Bonilla (2019) L. L. Bonilla, Phys. Rev. E 100, 022601 (2019).
Nguyen et al. (2021a) G. H. P. Nguyen, R. Wittmann, and H. Löwen, J. Phys. Condens. Matter 34, 035101 (2021a).
Wu and Libchaber (2000a) X.-L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000a).
Gal and Weihs (2010) N. Gal and D. Weihs, Phys. Rev. E 81, 020903 (2010).
Caspi et al. (2000) A. Caspi, R. Granek, and M. Elbaum, Phys. Rev. Lett. 85, 5655 (2000).
Bronstein et al. (2009) I. Bronstein, Y. Israel, E. Kepten, S. Mai, Y. Shav-Tal, E. Barkai, and Y. Garini, Phys. Rev. Lett. 103, 018102 (2009).
Newman et al. (2022) T. A. C. Newman, B. Beltran, J. M. McGehee, D. Elnatan, C. K. Cahoon, M. R. Paddy, D. B. Chu, A. J. Spakowitz, and S. M. Burgess, Proc. Natl. Acad. Sci. U.S.A 119, e2115883119 (2022).
Jeon et al. (2011) J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys. Rev. Lett. 106, 048103 (2011).
Lizana et al. (2016) L. Lizana, T. Ambjörnsson, and M. A. Lomholt, J. Phys. A Math. Theor. 50, 034001 (2016).
Narinder et al. (2019) N. Narinder, J. R. Gomez-Solano, and C. Bechinger, New J. Phys. 21, 093058 (2019).
Ando and Skolnick (2010) T. Ando and J. Skolnick, Proc. Natl. Acad. Sci. U.S.A 107, 18457 (2010).
Wong et al. (2004) I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch, and D. A. Weitz, Phys. Rev. Lett. 92, 178101 (2004).
Pollard and Cooper (2009) T. D. Pollard and J. A. Cooper, Science 326, 1208 (2009).
Lin et al. (2014) C. Lin, Y. Zhang, I. Sparkes, and P. Ashwin, Biophys. J. 107, 763 (2014).
Speckner et al. (2018) K. Speckner, L. Stadler, and M. Weiss, Phys. Rev. E 98, 012406 (2018).
Vale and Hotani (1988) R. D. Vale and H. Hotani, J. Cell Biol. 107, 2233 (1988).
Sanchez et al. (2012) T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
Bronshtein et al. (2015) I. Bronshtein, E. Kepten, I. Kanter, S. Berezin, M. Lindner, A. B. Redwood, S. Mai, S. Gonzalo, R. Foisner, Y. Shav-Tal, and Y. Garini, Nat. Commun. 6, 8044 (2015).
Granek et al. (2022) O. Granek, Y. Kafri, and J. Tailleur, Phys. Rev. Lett. 129, 038001 (2022).
Lau et al. (2003) A. W. C. Lau, B. D. Hoffman, A. Davies, J. C. Crocker, and T. C. Lubensky, Phys. Rev. Lett. 91, 198101 (2003).
Goychuk and Hänggi (2002) I. Goychuk and P. Hänggi, Proc. Natl. Acad. Sci. U.S.A 99, 3552 (2002).
Lizana et al. (2010) L. Lizana, T. Ambjörnsson, A. Taloni, E. Barkai, and M. A. Lomholt, Phys. Rev. E 81, 051118 (2010).
Saito and Sakaue (2015) T. Saito and T. Sakaue, Phys. Rev. E 92, 012601 (2015).
Vandebroek and Vanderzande (2017a) H. Vandebroek and C. Vanderzande, J. Stat. Phys. 167, 14 (2017a).
Joo et al. (2020) S. Joo, X. Durang, O.-c. Lee, and J.-H. Jeon, Soft Matter 16, 9188 (2020).
De Gennes (1976) P. G. De Gennes, Macromolecules 9, 587 (1976).
Panja (2010) D. Panja, J. Stat. Mech. Theory Exp. 2010, L02001 (2010).
Kappler et al. (2019) J. Kappler, F. Noé, and R. R. Netz, Phys. Rev. Lett. 122, 067801 (2019).
Granek (1997) R. Granek, Journal de Physique II 7, 1761 (1997).
Han et al. (2023) H.-T. Han, S. Joo, T. Sakaue, and J.-H. Jeon, J. Chem. Phys. 159, 024901 (2023).
Sung (2018) W. Sung, “Dynamics of polymers and membranes in fluids,” in Statistical Physics for Biological Matter (Springer Netherlands, Dordrecht, 2018) pp. 391–421.
Grest (2016) G. S. Grest, J. Chem. Phys. 145, 141101 (2016).
Buttinoni et al. (2013) I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Physical review letters 110, 238301 (2013).
Santra (2023) I. Santra, Journal of Physics: Complexity 4, 015013 (2023).
Qian (2000) H. Qian, Biophys. J 79, 137 (2000).
Chung and Hassanabadi (2017) W. S. Chung and H. Hassanabadi, Int. J. Theor. Phys. 56, 851 (2017).
Bonfanti et al. (2020) A. Bonfanti, J. L. Kaplan, G. Charras, and A. Kabla, Soft Matter 16, 6002 (2020).
Höfling and Franosch (2013) F. Höfling and T. Franosch, Rep. Prog. Phys. 76, 046602 (2013).
Kupferman (2004) R. Kupferman, J. Stat. Phys. 114, 291 (2004).
Kupferman et al. (2002) R. Kupferman, A. M. Stuart, J. R. Terry, and P. F. Tupper, Stoch. Dyn. 02, 533 (2002).
Lee et al. (2019) H. S. Lee, S.-H. Ahn, and E. F. Darve, J. Chem. Phys. 150, 174113 (2019).
Fodor et al. (2015) E. Fodor, D. S. Grebenkov, P. Visco, and F. van Wijland, Physica A 422, 107 (2015).
Szymanski and Weiss (2009) J. Szymanski and M. Weiss, Phys. Rev. Lett. 103, 038102 (2009).
Szymański et al. (2006) J. Szymański, A. Patkowski, J. Gapiński, A. Wilk, and R. Hołyst, J. Phys. Chem. B 110, 7367 (2006).
Jeon et al. (2013) J.-H. Jeon, N. Leijnse, L. B. Oddershede, and R. Metzler, New Journal of Physics 15, 045011 (2013).
Flenner et al. (2009) E. Flenner, J. Das, M. C. Rheinstädter, and I. Kosztin, Phys. Rev. E 79, 011907 (2009).
Jeon et al. (2012) J.-H. Jeon, H. M.-S. Monne, M. Javanainen, and R. Metzler, Phys. Rev. Lett. 109, 188103 (2012).
Javanainen et al. (2013) M. Javanainen, H. Hammaren, L. Monticelli, J.-H. Jeon, M. S. Miettinen, H. Martinez-Seara, R. Metzler, and I. Vattulainen, Faraday Discuss. 161, 397 (2013).
Porrà et al. (1996) J. M. Porrà, K.-G. Wang, and J. Masoliver, Phys. Rev. E 53, 5872 (1996).
Lü and Bao (2005) K. Lü and J.-D. Bao, Phys. Rev. E 72, 067701 (2005).
Vandebroek and Vanderzande (2017b) H. Vandebroek and C. Vanderzande, Soft Matter 13, 2181 (2017b).
Grimm and Dolgushev (2018) J. Grimm and M. Dolgushev, Soft Matter 14, 1171 (2018).
Sakaue and Saito (2017) T. Sakaue and T. Saito, Soft Matter 13, 81 (2017).
Muralidhar et al. (1991) R. Muralidhar, D. J. Jacobs, D. Ramkrishna, and H. Nakanishi, Phys. Rev. A 43, 6503 (1991).
Caspi et al. (2002) A. Caspi, R. Granek, and M. Elbaum, Phys. Rev. E 66, 011916 (2002).
Fa (2006) K. S. Fa, Phys. Rev. E 73, 061104 (2006).
Despósito and Viñales (2008) M. A. Despósito and A. D. Viñales, Phys. Rev. E 77, 031123 (2008).
Dechant et al. (2014) A. Dechant, E. Lutz, D. A. Kessler, and E. Barkai, Phys. Rev. X 4, 011022 (2014).
Kursawe et al. (2013) J. Kursawe, J. Schulz, and R. Metzler, Phys. Rev. E 88, 062124 (2013).
Cherstvy et al. (2021) A. G. Cherstvy, W. Wang, R. Metzler, and I. M. Sokolov, Phys. Rev. E 104, 024115 (2021).
Bao (2017) J.-D. Bao, J. Stat. Phys. 168, 561 (2017).
Burov and Barkai (2008) S. Burov and E. Barkai, Phys. Rev. E 78, 031112 (2008).
Fricks et al. (2009) J. Fricks, L. Yao, T. C. Elston, and M. G. Forest, SIAM J. Appl. Math. 69, 1277 (2009).
McKinley et al. (2009) S. A. McKinley, L. Yao, and M. G. Forest, J. Rheol. 53, 1487 (2009).
Grebenkov and Vahabi (2014) D. S. Grebenkov and M. Vahabi, Phys. Rev. E 89, 012130 (2014).
Grebenkov et al. (2013) D. S. Grebenkov, M. Vahabi, E. Bertseva, L. Forró, and S. Jeney, Phys. Rev. E 88, 040701 (2013).
Di Cairano et al. (2021) L. Di Cairano, B. Stamm, and V. Calandrini, Biophysical Journal 120, 4722 (2021).
de Oliveira and Tenreiro Machado (2014) E. C. de Oliveira and J. A. Tenreiro Machado, Math. Probl. Eng. 2014, 238459 (2014).
Mandelbrot and Ness (1968) B. B. Mandelbrot and J. W. V. Ness, SIAM Review 10, 422 (1968).
Lutz (2001) E. Lutz, Phys. Rev. E 64, 051106 (2001).
Jeon and Metzler (2010) J.-H. Jeon and R. Metzler, Phys. Rev. E 81, 021103 (2010).
Zhu et al. (2023) K. Zhu, S. Naserifar, W. A. Goddard, and H. Su, Phys. Chem. Chem. Phys. 25, 10353 (2023).
Kim et al. (2022) Y. Kim, S. Joo, W. K. Kim, and J.-H. Jeon, Macromolecules 55, 7136 (2022).
Sarkar et al. (2023) R. Sarkar, I. Santra, and U. Basu, Phys. Rev. E 107, 014123 (2023).
Uhlenbeck and Ornstein (1930) G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 (1930).
Szamel (2023) G. Szamel, Phys. Rev. E 107, 054602 (2023).
Bateman and Erdélyi (1955) H. Bateman and A. Erdélyi, Higher transcendental functions, California Institute of technology. Bateman Manuscript project (McGraw-Hill, New York, NY, 1955).
Haubold et al. (2011) H. J. Haubold, A. M. Mathai, and R. K. Saxena, J. Appl. Math. Comput. 2011, 298628 (2011).
Deng and Barkai (2009) W. Deng and E. Barkai, Phys. Rev. E 79, 011112 (2009).
Guo et al. (2013) P. Guo, C. Zeng, C. Li, and Y. Chen, Fract. Calc. Appl. Anal. 16, 123 (2013).
Note (1) Rigorously, in the limit of $\tau_{\mathrm{A}}\to 0$ for a fixed $D_{\mathrm{A}}$ , $C_{\mathrm{v,ac}}(t)\approx 2A\intop\ilimits@_{0}^{\infty}F_{1}(x+t)F_{1}(x)dx$ while $C_{v,\mathrm{th}}(t)\sim F_{1}(t)$ . Although $C_{\mathrm{v,ac}}(t)$ does not precisely converges to $C_{v,\mathrm{th}}(t)$ as $\tau_{\mathrm{A}}\to 0$ , two functions are numerically very similar to each other.
Siegle et al. (2010a) P. Siegle, I. Goychuk, and P. Hänggi, Phys. Rev. Lett. 105, 100602 (2010a).
Siegle et al. (2010b) P. Siegle, I. Goychuk, P. Talkner, and P. Hänggi, Phys. Rev. E 81, 011136 (2010b).
Wang and Tokuyama (1999) K. Wang and M. Tokuyama, Phys. A: Stat. 265, 341 (1999).
Godec and Metzler (2013) A. c. v. Godec and R. Metzler, Phys. Rev. Lett. 110, 020603 (2013).
Wu and Libchaber (2000b) X.-L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000b).
Nguyen et al. (2021b) G. H. P. Nguyen, R. Wittmann, and H. Löwen, J. Phys. Condens. Matter 34, 035101 (2021b).
Martin et al. (2021) D. Martin, J. O’Byrne, M. E. Cates, E. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, Phys. Rev. E 103, 032607 (2021).
Lisin et al. (2021) E. A. Lisin, O. S. Vaulina, I. I. Lisina, and O. F. Petrov, Phys. Chem. Chem. Phys. 23, 16248 (2021).
Cherstvy et al. (2018) A. G. Cherstvy, S. Thapa, Y. Mardoukhi, A. V. Chechkin, and R. Metzler, Phys. Rev. E 98, 022134 (2018).
Wang et al. (2020) W. Wang, A. G. Cherstvy, A. V. Chechkin, S. Thapa, F. Seno, X. Liu, and R. Metzler, J. Phys. A Math. Theor. 53, 474001 (2020).
Froemberg and Barkai (2013) D. Froemberg and E. Barkai, Phys. Rev. E 87, 030104 (2013).
He et al. (2008) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
Albers and Radons (2013a) T. Albers and G. Radons, Europhys. Lett. 102, 40006 (2013a).
Barkai et al. (2012) E. Barkai, Y. Garini, and R. Metzler, Phys. Today 65, 29 (2012).
Albers and Radons (2013b) T. Albers and G. Radons, Europhys. Lett. 102, 40006 (2013b).
Cherstvy et al. (2013) A. G. Cherstvy, A. V. Chechkin, and R. Metzler, New J. Phys. 15, 083039 (2013).
Massignan et al. (2014) P. Massignan, C. Manzo, J. A. Torreno-Pina, M. F. García-Parajo, M. Lewenstein, and G. J. Lapeyre, Phys. Rev. Lett. 112, 150603 (2014).
Vandebroek and Vanderzande (2015) H. Vandebroek and C. Vanderzande, Phys. Rev. E 92, 060601 (2015).
Diethelm et al. (2002) K. Diethelm, N. J. Ford, and A. D. Freed, Nonlinear Dyn. 29, 3 (2002).
Davies and Harte (1987) R. B. Davies and D. S. Harte, Biometrika 74, 95 (1987).
Beran (1998) J. Beran, Statistics for long-memory processes (Chapman & Hall, 1998).
Bardet (2002) J.-M. Bardet, IEEE Trans. Inf. 48, 991 (2002).
Bardet and Bertrand (2007) J.-M. Bardet and P. Bertrand, J. Time Ser. Anal. 28, 1 (2007).
Jia et al. (2019) J. Jia, Z. Wang, X. Huang, and Y. Wei, J. Funct. Spaces 2019, 6091602 (2019).

Viscoelastic active diffusion governed by nonequilibrium fractional Langevin equations: underdamped dynamics and ergodicity breaking

Abstract

I Introduction

II The Model & Method

II.1 The active fractional Langevin equation model

II.2 Numerical method

III Results: Active diffusion in viscoelastic media

III.1 Velocity autocorrelation function (VACF)

III.2 Mean-squared displacements (MSDs)

III.2.1 The ensemble-averaged MSD

III.2.2 The time-averaged MSD

III.3 Ergodicity-breaking parameter

IV Concluding remarks

Acknowledgements

Appendix A Validation of numerical method & FLE with thermal noise

Appendix B Generalized Mittag-Leffler functions

Appendix C Supplementary figures

References

References

Appendix A Validation of numerical method &
FLE with thermal noise