diff --git a/src/ChaosTools.jl b/src/ChaosTools.jl
index 8012283e..7e7e61c9 100644
--- a/src/ChaosTools.jl
+++ b/src/ChaosTools.jl
@@ -26,6 +26,7 @@ include("periodicity/davidchacklai.jl")
 include("rareevents/mean_return_times/mrt_api.jl")
 
 include("chaosdetection/lyapunovs/lyapunov.jl")
+include("chaosdetection/lyapunovs/EAPD.jl")
 include("chaosdetection/lyapunovs/lyapunov_from_data.jl")
 include("chaosdetection/lyapunovs/lyapunovspectrum.jl")
 include("chaosdetection/lyapunovs/local_growth_rates.jl")
diff --git a/src/chaosdetection/lyapunovs/EAPD.jl b/src/chaosdetection/lyapunovs/EAPD.jl
new file mode 100644
index 00000000..e652afd8
--- /dev/null
+++ b/src/chaosdetection/lyapunovs/EAPD.jl
@@ -0,0 +1,155 @@
+export ensemble_averaged_pairwise_distance,lyapunov_instant
+
+"""
+    lyapunov_instant(ρ,times;interval=1:length(times)) -> λ(t)
+
+Convenience function that calculates the instantaneous Lyapunov exponent by taking the slope of 
+the ensemble-averaged pairwise distance function `ρ` wrt. to the saved time points `times` in `interval`.
+"""
+function lyapunov_instant(ρ,times;interval=1:length(times))
+    _,s = linreg(times[interval], ρ[interval]) #return estimated slope  
+    return s
+end
+
+"""
+    ensemble_averaged_pairwise_distance(ds, init_states::StateSpaceSet, T, pidx;kwargs...) -> ρ,t
+
+Calculate the ensemble-averaged pairwise distance function `ρ` for non-autonomous dynamical systems 
+with a time-dependent parameter, using the metod described by [^Jánosi, Tél]. Time-dependence is assumed to be a linear drift. The rate of change
+of the parameter needs to be stored in the parameter container of the system `p = current_parameters(ds)`,
+at the index `pidx`. In case of autonomous systems (with no drift), `pidx` can be set to any index as a dummy.
+To every member of the ensemble `init_states`, a perturbed initial condition is assigned.
+`ρ(t)` is the natural log of phase space distance between the original and perturbed states averaged 
+over all pairs, calculated for all time steps up to `T`. 
+
+
+## Keyword arguments
+
+* `initial_params = deepcopy(current_parameters(ds))`: initial parameters 
+* `Ttr = 0`: transient time used to evolve initial states to reach 
+    initial autonomous attractor (without drift)
+* `perturbation = perturbation_normal`: if given, it should be a function `perturbation(ds,ϵ)`,
+   which outputs perturbed state vector of `ds` (preferrably `SVector`). If not given, a normally distributed 
+   random perturbation with norm `ϵ` is added.
+*  `Δt = 1`: step size 
+*  `ϵ = sqrt(dimension(ds))*1e-10`: initial distance between pairs of original and perturbed initial conditions 
+
+
+## Description
+In non-autonomous systems with parameter drift, long time averages are less useful to assess chaoticity.
+Thus, quantities using time averages are rather calculated using ensemble averages. Here, a new 
+quantity called the Ensemble-averaged pairwise distance (EAPD) is used to measure chaoticity of 
+the snapshot attractor/ phase space object traced out by the ensemble [^Jánosi, Tél].
+
+To any member of the original ensemble (`init_states`) a close neighbour (test) is added at an initial distance `ϵ`. Quantity `d(t)` is the 
+dimensionless phase space distance between a test particle and an ensemble member at time t .
+If `init_states`` are randomly initialized (far from the attractor at the initial parameter), and there's no transient, 
+the first few time steps cannot be used to calculate any reliable averages.
+The function of the EAPD `ρ(t)` is defined as the average logarithmic distance between original and 
+perturbed initial conditions at every time step: `ρ(t) = ⟨ln d(t)⟩`
+
+An analog of the classical largest Lyapunov exponent can be extracted from the 
+EAPD function `ρ`: the local slope can be considered an instantaneous Lyapunov exponent.
+
+## Example
+Example of a rate-dependent (linearly drifting parameter) system
+
+```julia
+#r parameter is replaced by r(n) = r0 + R*n drift term
+function drifting_logistic(x,p,n)
+    r0,R = p
+    return SVector( (r0 + R*n)*x[1]*(1-x[1]) )
+end
+
+r0 = 3.8 #inital parameter
+R = 0.001 #rate parameter
+p = [r0,R] # pidx = 2
+
+init_states = StateSpaceSet(rand(1000)) #initial states of the ensemble 
+ds = DeterministicIteratedMap(drifting_logistic, [0.1], p)
+ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,2;Ttr=1000)
+```
+
+[^Jánosi, Tél]: Dániel Jánosi, Tamás Tél, Physics Reports **1092**, pp 1-64 (2024)
+
+"""
+function ensemble_averaged_pairwise_distance(ds,init_states::StateSpaceSet,T,pidx;
+    initial_params = deepcopy(current_parameters(ds)),Ttr=0,perturbation=perturbation_normal,Δt = 1,ϵ=sqrt(dimension(ds))*1e-10)
+
+	set_parameters!(ds,initial_params)
+    original_rate = current_parameter(ds, pidx)
+    N = length(init_states)
+    d = dimension(ds)
+    dimension(ds) != d && throw(AssertionError("Dimension of `ds` doesn't match dimension of states in init_states!"))
+    
+    nt = length(0:Δt:T) #number of time steps
+    ρ = zeros(nt) #store ρ(t)
+    times = zeros(nt) #store t
+    
+    #duplicate every state
+    #(add test particle to every ensemble member)
+    init_states_plus_copies = StateSpaceSet(vcat(init_states,init_states))
+
+    #create a pds for the ensemble
+    #pds is a ParallelDynamicalSystem
+    pds = ParallelDynamicalSystem(ds,init_states_plus_copies)
+
+	#set to non-drifting for initial ensemble
+    set_parameter!(pds,pidx,0.0) 
+
+    #step system pds to reach attractor(non-drifting)
+    #system starts to drift at t0=0.0
+    for _ in 0:Δt:Ttr
+        step!(pds,Δt,true)
+    end
+    
+    #rescale test states 
+    #add perturbation to test states
+    for i in 1:N 
+        state_i = current_state(pds,i)
+        perturbed_state_i = state_i .+ perturbation(ds,ϵ)
+        #set_state!(pds.systems[N+i],perturbed_state_i)
+        set_state!(pds,perturbed_state_i,N+i)
+    end
+
+	#set to drifting for initial ensemble
+    set_parameter!(pds,pidx,original_rate) 
+
+    #set back time to t0 = 0
+    reinit!(pds,current_states(pds))
+
+    #calculate EAPD for each time step
+    ensemble_averaged_pairwise_distance!(ρ,times,pds,T,Δt)
+    return ρ,times
+
+end
+
+#calc distance for every time step until T
+function ensemble_averaged_pairwise_distance!(ρ,times,pds,T,Δt)
+    for (i,t) in enumerate(0:Δt:T)
+        ρ[i] = ensemble_averaged_pairwise_distance(pds)
+        times[i] = current_time(pds)
+        step!(pds,Δt,true)
+    end
+end
+
+#calc distance for current states of pds
+function ensemble_averaged_pairwise_distance(pds)
+
+    states = current_states(pds)
+    N = Int(length(states)/2)
+
+    #calculate distance averages
+    ρ = 0.0
+    for i in 1:N
+        ρ += log.(norm(states[i] - states[N+i]))
+    end
+    return ρ/N 
+
+end
+
+function perturbation_normal(ds,ϵ)
+    D, T = dimension(ds), eltype(ds)
+    p0 = randn(SVector{D, T})
+    p0 = ϵ * p0 / norm(p0)  
+end
diff --git a/test/chaosdetection/EAPD.jl b/test/chaosdetection/EAPD.jl
new file mode 100644
index 00000000..dca45284
--- /dev/null
+++ b/test/chaosdetection/EAPD.jl
@@ -0,0 +1,59 @@
+using ChaosTools, Test
+using LinearAlgebra
+
+henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+henon() = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3, 0.0]) # add third dummy parameter
+
+#test if ensemble averaging gives the same as 
+#the usual lyapunov exponent for autonomous system
+@testset "time averaged and ensemble averaged lyapunov exponent" begin
+    ds = henon()
+
+    #eapd slope
+    init_states = StateSpaceSet(0.2 .* rand(1000,2))
+    pidx = 3 # set to dummy, not used anywhere (no drift)
+    ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,pidx;Ttr=5000)
+    lyap_instant = lyapunov_instant(ρ,times;interval=20:30)
+
+    #lyapunov exponent
+    λ = lyapunov(ds,1000;Ttr=5000)
+    @test isapprox(lyap_instant,λ;atol=0.01)
+end
+
+#test sliding Duffing map 
+#-------------------------duffing stuff-----------------------
+#https://doi.org/10.1016/j.physrep.2024.09.003
+
+function duffing_drift(u0 = [0.1, 0.25]; ω = 1.0, β = 0.2, ε0 = 0.4, α=0.00045)
+    return CoupledODEs(duffing_drift_rule, u0, [ω, β, ε0, α])
+end
+
+@inbounds function duffing_drift_rule(x, p, t)
+    ω, β, ε0, α = p
+    dx1 = x[2]
+    dx2 = (ε0+α*t)*cos(ω*t) + x[1] - x[1]^3 - 2β * x[2]
+    return SVector(dx1, dx2)
+end
+
+@testset "Duffing map" begin
+    #----------------------------------hamiltonian case--------------------------------------
+    duffing = duffing_drift(;β = 0.0,α=0.0,ε0=0.08) #no dissipation -> Hamiltonian case
+    duffing_map = StroboscopicMap(duffing,2π)
+    init_states_auto,_ = trajectory(duffing_map,5000,[-0.85,0.0];Ttr=0) #initial condition for a snapshot torus
+    #set system to sliding
+    set_parameter!(duffing_map,4,0.0005)
+    pidx=4
+
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,init_states_auto,100,pidx;Ttr=0)
+    lyap_instant = lyapunov_instant(ρ,times;interval=50:60)
+    @test isapprox(lyap_instant,0.87;atol=0.01) #0.87 approximate value from article
+
+    #-----------------------------------dissipative case------------------------------------
+    duffing = duffing_drift() #no dissipation -> Hamiltonian case
+    duffing_map = StroboscopicMap(duffing,2π)
+    init_states = randn(5000,2) 
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,StateSpaceSet(init_states),100,pidx;Ttr=20)
+    lyap_instant = lyapunov_instant(ρ,times;interval=2:20)
+    @test isapprox(lyap_instant,0.61;atol=0.01) #0.61 approximate value from article
+
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f4fad1ef..7596c8d0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,6 +25,7 @@ testfile("chaosdetection/gali.jl")
 testfile("chaosdetection/partially_predictable.jl")
 testfile("chaosdetection/01test.jl")
 testfile("chaosdetection/expansionentropy.jl")
+testfile("chaosdetection/EAPD.jl")
 
 testfile("stability/fixedpoints.jl")
 testfile("periodicity/periodicorbits.jl")