Since the pioneering work of Einstein [1], Brownian motion has become the paradigm for the description and understanding of a large variety of diffusion processes that are present in numerous physical, biological, and chemical systems. In recent years, the dynamics of macromolecules and nanoparticles on surfaces or curved spaces has been the subject of intensive investigations, especially because particle diffusion shows a richer dynamical behavior at different time scales [2, 3] than its counterpart in open and flat geometries. In particular, diffusion plays a key role in the dynamics of molecular motors moving along heterogeneous substrates [4], in the transport of biomacromolecules in the cell due to crowding [5, 6], and in the lateral diffusion of proteins on fluctuating membranes [7, 8]. Membranes can be simplified as two-dimensional differentiable manifolds based on their typical length scales. This approach allows us to describe many of its characteristics, such as the shape and dynamics, as if they were continuous pieces of material [9–11]. As a result, particles that move along the membrane due to certain underlying interactions can be explained as the random motion of particles restricted to the surface. Although this perspective might ignore relevant molecular details, most of the diffusion properties of tiny particles confined ¹ on surfaces depend strongly on the generic features of the surface or, strictly speaking, on the surface geometry [12]. Typically, particle dynamics is not only influenced by geometrical features but also by local and thermodynamic properties that experience the effects of the geometry of the manifold where the particles are embedded [13–15]. In a more general setting, the relevance of studying Brownian motion in differential manifolds is due to its ability to describe the system of Brownian particles with holonomic constraints that can model their interactions with the medium that houses them, as explicitly discussed in [16, 17] and recently in [18].

A great effort for understanding Brownian motion on surfaces can be found in colloidal soft matter, where the dynamics of colloidal particles on quasi-two-dimensional geometries have been both experimentally and theoretically investigated by using optical techniques such as digital videomicroscopy, computer simulations, and theoretical approximations [19, 20]. Nonetheless, such investigations deal basically with (almost) flat surfaces, i.e., without including curvature effects. The interest in the use of colloids resides in the fact that they are small (nanometer to micrometer-sized) particles and typically are considered model systems because of other interesting features [21]. Their characteristic time and length scales are experimentally accessible, which allows us to follow the colloidal dynamics and transport processes in real-time [21]. Furthermore, since the colloidal interactions are relatively weak, colloids are highly susceptible to external forces, and hence, their static and dynamical properties can be controlled through the application of external fields or by imposing geometrical restrictions. In other words, it is assumed that the molecular forces holding particles to the manifold are stronger than the forces between colloids (see, e.g., Ref. [21] and references therein). Then, colloids represent an ideal model system to account for the effects of geometry on the nature and dynamics of many-body systems.

In particular, it has already been demonstrated and experimentally corroborated that the curvature dependence of a fluctuating membrane affects the diffusion processes of molecules on the membrane surface [22–24]. These geometrical effects, although important, are still difficult to interpret. The lack of a precise interpretation resides in the fact that, unfortunately, there is not a unique way to define diffusion observables on a curved surface (see, for instance [22, 23, 25]). The description of colloidal dynamics in curved spaces is a non-trivial task; it represents a formidable physical and mathematical challenge. Recently, one of us proposed the generalization of the Smoluchowski equation on curved spaces [26]. Furthermore, Castro–Villarreal also put forward different geometrical observables to quantify the displacement of a single colloidal particle [25]. Within this approach, it was shown that the geodesic mean-squared displacement captures the intrinsic elements of the manifold, whereas the Euclidean displacement provides extrinsic information from the surface. An interesting extension of the theory now provides the description of the motion of active Brownian particles [27], where the mean-squared geodesic displacement captures the relationship between the curvature and the activity of the active colloid. This theoretical framework provided evidence that an active Brownian particle experiences a dynamical transition in any compact surface from a monotonic to an oscillating behavior observed in the mean-squared geodesic displacement [27]; a theoretical prediction of a dynamic transition of this type can be established using a run-and-tumble active particle confined on circle S ¹. This prediction was recently corroborated in experiments using a non-vibrating magnetic granular system (see, e.g., Ref. [28]). However, we still face challenges in colloidal dynamics on curved manifolds, for example, the generalization of this approach to the situation where the colloids interact not only with other macromolecules, i.e., direct forces, but also with the inclusion of all those geometrical mechanisms originating from the curvature and to situations when passive and active colloids are not just restricted geometrically but also are immersed in a viscoelastic fluid [29]. Furthermore, the transitions occurring in compact manifolds previously are examples of a non-trivial recurrence. It would be interesting to approach the problem from the perspective of dynamical system theory [30].

The aforementioned theoretical formalism has also allowed us to determine the equation of motion of interacting colloids in curved spaces; a generalized Ermack–McCammon algorithm has been developed to study a broader class of transport phenomena in curved manifolds [31]. Interestingly, the predictions of the particle transport in non-Euclidean spaces have been partially corroborated in a series of experiments [3, 24]; superparamagnetic colloids embedded in a circle and subjected to external magnetic fields [3] and polystyrene nanoparticles diffusing on highly curved water–silicone oil interfaces [24]. However, further experimental, computational, and theoretical studies are needed to better understand the rich diffusion mechanisms, particle distribution, and thermodynamic properties that emerge in colloidal dispersions when the curvature of the space plays an important role.

From a rigorous mathematical perspective, the study of stochastic processes on differential manifolds has a long history, dating back almost a century ago. A. Kolmogoroff made a seminal contribution to the covariant formulation of the Fokker–Planck (FP) equation (32), which was further developed in [33, 34] by proving the existence of solutions to the covariant FP equation and stochastic differential equations in manifolds, respectively. In addition, through different limiting cases reported in [35], when t → 0, asymptotic expressions were observed for the probability density function ρ(ξ, ξ′, t) in the cases of nearby and distant points [35]. These covariant versions of the Fokker–Planck equation are also introduced in [36], where one can identify a Riemannian geometry defined by a metric tensor given by the diffusion tensor [37], while a recent formulation starting from a coordinate covariance and gauge invariance in [38] was able to show the existence of a Langevin equation with a space-dependent diffusion matrix, whose corresponding Fokker–Planck equation is also given in a covariant form. The covariant formulation of stochastic equations has been further developed intensely during the last decade (e.g., see, [18, 39]). In this work, we take a little step forward to develop a covariant approach for a system of interacting colloidal particles in a sub-manifold of the Euclidean space M ⊂ R d + 1 with dimension d without the explicit inclusion of hydrodynamic interactions. The physical scenario consists of N interacting colloidal particles moving on M ; particularly, we choose the case when M is a curved surface with d = 2, but the formulation derives the steps to generalize to other sub-manifolds. There are four physical assumptions that we need to take into account. To begin with, the shape of the sub-manifold must remain constant over time. In addition, the interaction between colloids and the molecules that constitute the sub-manifold (or the external field that defines the sub-manifold) is stronger than the interaction between colloids. The collective effects of the solvent molecules on the colloid are described by a stochastic force that is defined in the Euclidean space R d + 1 . Finally, the size of the colloid is significantly smaller than the usual lengths of the sub-manifold.

The article is structured as follows. Section 2 introduces the covariant approach of the colloidal dynamics from a stochastic differential equation in the Stratonovich sense. In particular, we can simplify the equation for the N-particle system moving on M to the equation of a one-particle system moving in a hyper-dimensional manifold M N . Section 3 uses the covariant approach to analyze the short-time behavior of the joint probability density and the mean-squared displacement of a tagged particle from the interacting system. Furthermore, we present two examples of a hard-sphere system and the soft interaction to estimate the order of the curvature when particles are confined to the sphere. In particular, it presents an explicit expression for the mean-squared displacement when the interaction is central. Lastly, Section 4 concludes with our remarks, challenges, and perspectives on the future of soft condensed matter explicitly applying our covariant approach.

As discussed previously, one of the main challenges in understanding the effects of geometry on the dynamics of colloids embedded in a curved space is to develop experimental tools and theoretical frameworks that account for the transport properties that occur on the manifold. In the following section, we then provide the first preliminary steps to build a covariant theoretical formulation of the dynamics of an interacting colloidal system based on the many-body Langevin equation in the so-called overdamped limit [31], which allows us to deduce a Smoluchowski equation [40] for the interacting system on the manifold. We should mention that, without loss of generality, within this approach, hydrodynamic interactions are disregarded, but their explicit inclusion will be reported elsewhere. Before starting with the covariant formulation, let us introduce our notation. Let us consider the colloidal system confined on a d − dimensional manifold M embedded in a d + 1 − dimensional Euclidean space R d + 1 and described with the parameterization X : U ⊂ R d → R d + 1 , where a particular point in M is given by X(x), being x ≡ x 1 , x 2 , … , x d ∈ U local coordinates of the neighborhood U. Using the embedding function X(x), one can define a Riemannian metric tensor by g _αβ = e _α ⋅e _β , where e α = ∂ ∂ x α X x , with α = 1, … , d. Further notions like normal vector, extrinsic curvature tensor, and Weingarten–Gauss equations are introduced in Appendix A from [25]. Typically, spatial dimensions of interest are d = 1 and d = 2.

As we have pointed out previously, our starting point to describe the dynamics of colloids confined in a curved manifold is based on a previous contribution [31], where the many-body Langevin stochastic equations are posed in the overdamped regime, i.e., the diffusive time scale, in local coordinates, is written as x ̇ i α = 1 ζ e α x i ⋅ f i t + ∑ j ≠ i F i j x i , x j , (1) where ζ is the friction coefficient and with x i α being the i − th particle position with i = 1, … , N and x ̇ i α ≡ d x i α d t . The quantity f i t represents the collective effects of the solvent molecules on the colloid, and it is expressed by a stochastic force over the ith-particle, which satisfies the fluctuation–dissipation theorem in the Euclidean space R d + 1 , that is, f i t = 0 and f i ( t ) f j ( τ ) = 2 ζ k B T 1 δ i j δ t − τ , where k _B T is the thermal energy with T being the temperature and k _B the Boltzmann’s constant. The term F i j x i , x j is the force that the ith-particle experiences at the position x _i and is due to the interaction with the jth-particle located at the position x _j . In Eq. 1, the tangent vector e _α ≡ ∂ _α X projects the dynamics on the tangent space T X M since the dynamics is occurring intrinsically on the manifold. Note that rising and lowering indices are carried out by the standard fashion using the metric tensor to lowering indices and inverse metric tensor g ^αβ for rising indices, for instance, v ^α = g ^αβ v _β for an arbitrary vector v.

In the present exposition, we adopt the consideration that Eq. 1 is a set of N stochastic differential equations in Stratonovich’s sense [41], d x i α = 1 ζ F i α d t + 2 D 0 e i , a α d W i , a t , (2) where F i α ≡ ∑ j ≠ i F i j α , with F i j α as the tangent projection of the interacting term F _ij , and D ₀ = k _B T/ζ is the collective-diffusion coefficient. In addition, there is an implicit sum over the indices a = 1, … , d + 1 to take into account the tangent projection with the stochastic term in Eq. 1, which has been identified with a Wiener process for each particle d W _i (t) = (dW _i,1(t), dW _i,2(t), …, dW _i,d+1(t)), so that the total Wiener process d W (t) is such that dim[d W (t)] = (d + 1)N. Since the dynamics occurs in the curved space, the Wiener process should also be projected on it. Therefore, we introduce a block diagonal projection operator P ̂ = diag ( e 1 , a α , e 2 , a α , … , e N , a α ) , with e i , a α ≡ e α x i a , where the blocks are individual operators for each particle given by the tensorial product on the basis of the tangent space and the basis of the Euclidean space.

It is a well-known fact that given a set of differential stochastic equation in the Stratonovich form, such as Eq. 2, one can find its associate Chapman–Kolmogorov differential equation for the joint probability density function p : M N × R → R [41]. For this, we only have to identify the components of the drift vector and the diffusion matrix, which, in this case, are A i α = F i α / ζ and B i , a α = 2 D 0 e i , a α , respectively. Then, we obtain the following expression: ∂ t p = − 1 ζ ∑ i = 1 N ∂ α F i α p + D 0 ∑ i = 1 N ∂ α e i , a α ∂ β e i , a β p . (3)

In this equation, let us note that the partial derivation ∂ α = ∂ ∂ x i α depends on the index i, which is associated with the particle label. Although this last equation has information on the geometry of the surface through the tangent vectors, it is not written in a covariant form yet. To this end, we define the probability density appropriately normalized with the volume element d V = ∏ i = 1 N d v g i , where d v g i is the Riemannian volume element defined by d v g i ≡ d d x i g ( x i ) for each particle. Thus, it is convenient to define a covariant joint probability density function ρ x 1 , … , x N ; t such as p x 1 , … , x N ; t = ∏ i = 1 N g ( x i ) ρ x 1 , … , x N ; t , where g(x _i ) is the determinant of the metric tensor g α β x i . After this change and using the Weingarten–Gauss equation mentioned previously, Eq. 3 takes the following mathematical form: ∂ t ρ = − 1 ζ ∑ i = 1 N ∇ α , i F i α ρ + D 0 ∑ i = 1 N 1 g ∂ α g α β g ∂ β ρ + g α β ρ ∂ β g − g Γ ν β ν , where the covariant derivative acting on a vector field v ^α is ∇ α , i v α = 1 g ∂ α g v α using the coordinates of the i − th particle, and x i α and Γ ν β α are the Christoffel symbols [42]. Additionally, applying the identity Γ ν β ν = ∂ β ⁡ log g and identifying that the Laplace–Beltrami operator acts on the scalars Δ g , i = g − 1 ∂ α ( g α β g ∂ β ) (also using local coordinates, x i α ), it is straightforward to obtain the desired covariant expression: ∂ t ρ = D 0 ∑ i = 1 N Δ g , i ρ − 1 ζ ∑ i = 1 N ∇ α , i F i α ρ . (4)

Equation 4 represents the covariant formulation of the Smoluchowski equation of a colloidal system of interacting particles constrained to a curved space M , where all the geometrical features are included in the Laplace–Beltrami operator and the covariant derivative. This equation is reduced to the standard Smoluchowski equation when the manifold M is the open Euclidean space R d , where the metric tensor is g _αβ = δ _αβ . Notice that hydrodynamic interactions have not been included in the derivation of Eq. 4.

Furthermore, one can write down Eq. 4 in a more compact form that allows us to prove that both systems shown in Figure 1, that is, the system of N interacting particles confined to a d-dimensional manifold M , and the system of a single particle in an external force confined to a D -dimensional manifold M represent equivalent systems. For this purpose, let us define a hyper-dimensional Riemannian geometry by N cartesian products of the manifold M , that is, M = M × M × ⋯ × M ≡ M N of dimension D = N d , where a local patch is described with the local coordinate ξ A = x i α , where the Greek superscript, α, describes the local coordinate component, while the Latin subscript, i, describes the particle and A = 1 , … , D . Now, this manifold M is equipped with a Riemannian metric defined through the following line element: d s 2 = ∑ i = 1 N g α β x i d x i α d x i β , (5) in terms of the metric tensor g _αβ of the coordinates of each particle. Thus, the metric tensor associated with the line element (5) for the manifold M is given by the block diagonal matrix G A B = diag g α β x 1 , … , g μ ν ( x N ) . It is not difficult to see that the covariant derivative compatible with the metric G _AB for the manifold M can be written as ∇ A = ∇ α , 1 , ∇ β , 2 , … , ∇ μ , N , (6) and the corresponding Laplace–Beltrami operator acting on scalars is simply the sum of each Laplace–Beltrami operator: Δ G = ∇ A ∇ A = ∑ i = 1 N Δ g , i . (7)

Now, defining F A = F 1 α , F 2 β , … , F N μ as the components of a vector field at the point ξ ∈ M , it is straightforward to write down the Smoluchowski equation for the full N − particle colloidal system confined on the curved space (4) as ∂ t ρ = D 0 Δ G ρ − 1 ζ ∇ A F A ρ . (8)

By expressing the Smoluchowski equation in this compact manner, it is now clear in what sense one can interpret the problem of the interacting colloidal system as the Brownian motion of a single particle in an external field F but in a hyper-dimensional space M . This identification was already implemented in a previous contribution [43], where an interacting colloidal system was studied on the line. Moreover, if we suppose that the interaction forces encoded in F A can be written as F A = − ∇ A Φ , where Φ is a certain interacting potential, one can see that the expected equilibrium distribution is satisfied at long times, namely, ρ ξ , t = 1 Z e − β Φ ξ , where Z is the partition function for the particle system confined to the curved manifold: Z = ∫ ∏ i = 1 N d v g i e − β Φ ξ , (9) where β = 1/(ζD ₀) = 1/(k _B T). Let us note that the expression of this partition function can also be obtained by integrating the momentum p i α variables from the Boltzmann weight using the Hamiltonian H = ∑ i = 1 N 1 2 m p i α g α β ( x i ) p i β + Φ ξ . Usually, the potential Φ ξ is considered a pairwise additive; thus, one can carry out the usual cluster diagrammatic expansion for the colloidal system in the curved space in equilibrium conditions [12].

Consequently, Eqs (4) and (8) represent the starting point of a covariant description, without the inclusion of hydrodynamic forces that allows us to study in detail the colloid dynamics in curved spaces. In the following paragraphs, we will discuss some applications of this formulation and highlight some challenges and future perspectives that can be tackled within this approach.

In this section, we study an application of the advantage of deriving the Smoluchowski equation in curved spaces in a covariant formulation (8). This consists in providing a general behavior of the joint probability density function at the short-time regime, or equivalently, in a neighborhood around a point of the manifold M . This calculation allows us to determine the short-time expansion of the self-diffusion coefficient; this is carried out in the following section by calculating the mean-squared displacement at the short-time regime. Since M is a Riemannian manifold with the metric tensor G _AB , one can explore the curvature effects on the colloidal interacting system using the Riemann normal coordinates (RNC) (see, e.g., Refs. [42, 44, 45]) in the neighborhood of a point p ∈ M in an entirely analog manner as it has been performed for a single particle [26]. To derive an approximate expression for the joint probability density function (PDF) at a short time, it is common to write the Smoluchowski Equation 8 as a heat-kernel equation: ∂ t + O ̂ ρ ξ , ξ ′ , t = 1 G δ ξ − ξ ′ δ t , (10) where ρ ξ , ξ ′ , t d V is the probability to find a diffusing particle in the hyper-volume element dV centered in ξ, at time t, when the particle started in ξ′ at t = 0. The operator O ̂ is defined as O ̂ = − D 0 Δ G + 1 ζ ∇ A F A ⋅ . At the initial condition, t → 0, the PDF acquires the form of a Dirac delta: ρ ( ξ , ξ ′ , t → 0 ) = 1 G δ ξ − ξ ′ . This initial condition establishes that the system is at the configuration ξ′ at the starting time. Then, by performing a Fourier transform on the time parameter, the aforementioned equation can be written as i E + O ̂ ρ ξ , ξ ′ , E = δ ξ − ξ ′ / G , where E is the time-conjugate Fourier variable. We should remark that ρ ξ , ξ ′ , t is equal to ρ(ξ, t) with the initial condition that ξ = ξ′ at time t = 0.

In the following section, we use the De Witt procedure [46], that is, we first separate the points to write the term G in front of the Dirac delta as the expression G → G 1 4 ( ξ ) G 1 4 ξ ′ . Now, we redefine the PDF as ρ ̄ ξ , ξ ′ , t = G 1 4 ξ ρ ξ , ξ ′ , t G 1 4 ξ ′ . Thus, after some algebraic rearrangements, Eq. 10 can be rewritten as i E + H ̂ ρ ̄ ξ , ξ ′ , E = δ ξ − ξ ′ , (11) where H ̂ = G 1 4 O ̂ G − 1 4 , or explicitly, this operator has the following mathematical form: H ̂ = − D 0 ∂ A G A B ∂ B + G − 1 4 ∂ A G 1 2 G A B ∂ B G − 1 4 − β ∂ A F A ⋅ + 1 4 G − 1 ∂ A G F A . (12)

Next, we take the advantage of the fact that the Smoluchowski equation is invariant under a general coordinate transformation. Thus, we choose Riemann normal coordinates (RNC) y ^A in a local neighborhood N ξ ′ ∈ M centered at ξ′. In RNC, the neighborhood N _ξ′ looks like Euclidean space, so we choose ξ′ to be the origin of this Euclidean space ² . The advantage of these coordinates is that one can express the metric tensor as G A B = δ A B + 1 3 R A C D B y C y D + ⋯ , where R A C D B is the Riemann curvature tensor of M evaluated at ξ′. In addition, we express the interaction terms in a Taylor expansion around the origin of the neighborhood F A ξ = F A ξ ′ + ∇ B F A ξ ′ y B + 1 2 ∇ B ∇ C F A ξ ′ y B y C + ⋯ , where the coefficients are evaluated at point ξ′.

In the subsequent section, we have all the pieces to split the operator (12) as H ̂ = H ̂ 0 + H ̂ I , where H ̂ 0 = D 0 p ̂ 2 − D 0 6 R + D 0 β ∇ A F A , (13) is a free “Hamiltonian” and H ̂ I = D 0 β ∇ B ∇ A F A − 1 6 R B A F A y B − D 0 β 6 R B A ∇ C F A y B y C + i D 0 β F A p ̂ A + i D 0 β ∇ B F A y B p ̂ A + i D 0 β 2 ∇ B ∇ C F A y B y C p ̂ A − D 0 3 R C A B D p ̂ A y C y D p ̂ B , (14) an interacting “Hamiltonian,” where we have defined a “momentum operator” as p ̂ A = − i ∂ A in an analogy with quantum mechanics. Now, the solution for the probability density function can be obtained by identifying δ ( ξ − ξ ′ ) = ξ ξ ′ and solving Eq. 11 as follows: ρ ̄ ξ , ξ ′ , E = ξ K ̂ ξ ′ , where K ̂ = 1 / ( i E + H ̂ ) is the resolvent operator. Next, we carry out a standard perturbation theory at the first order again in an entire analogy with quantum mechanics: Thus, the approximation of the resolvent operator through the perturbation theory is K ̂ = K ̂ 0 + K ̂ 0 H ̂ I K ̂ 0 + ⋯ . At this approximation, there are just six terms to evaluate, corresponding to the quantities of the form I i ξ , ξ ′ = ξ K ̂ 0 O ̂ K ̂ 0 ξ ′ , with i = 1, … , 6, where O ̂ i is one of the six terms: y ^B , y ^B y ^C , p ̂ A , y B p ̂ A , y B y C p ̂ A , and p ̂ A y C y D p ̂ B , respectively. Since K ̂ 0 depends just on the “momentum operator” p ̂ , it is convenient to introduce two completeness relations using the momentum basis { p } to compute the contributions from the interacting Hamiltonian. Hence, one can write I i ξ , ξ ′ = ∫ d D p 2 π D ∫ d D q K 0 p , α * e i ξ ⋅ p p O ̂ i q K 0 q , α * e − i ξ ′ ⋅ q , (15) where α * = − D 0 6 R + D 0 β ∇ A F A and K ₀(p, α) = 1/(iE., + D ₀ p ² + α) are simply functions of the value of the “momentum” p = p A p A and energy E. In addition, we have used the transformation from the position to the momentum basis as usual ξ | p = e i ξ ⋅ p / ( 2 π ) D 2 . We should recall that we have chosen ξ′ = 0 as the origin of the neighborhood N _ξ′; this allows us to simplify the calculation of the integrals I i ξ , ξ ′ . In Supplementary Appendix S4.1, we explicitly explain the procedure implemented to evaluate these integrals. After a straightforward calculation, the short-time approximation for the probability density function ρ(ξ, 0, t) of the full interacting system can be written as G ρ ξ , 0 , t = 1 4 π D 0 t D / 2 e − ξ 2 4 D 0 t 1 + τ 0 + τ B 1 ξ B + τ B C 2 ξ B ξ C + ⋯ , (16) where the terms τ ⁽⁰⁾, τ B ( 1 ) , and τ B C ( 2 ) are tensors given by τ 0 = D 0 t 1 6 R − 1 2 β ∇ A F A , (17) τ B 1 = β 2 G B A 1 + D 0 t R 6 − β ∇ C F C + D 0 t 6 R B A + G B A Δ G − 16 ∇ B ∇ A F A , (18) τ B C 2 = β 4 1 + D 0 t R 6 − β ∇ A F A ∇ B F C − 2 D 0 t 9 R B A ∇ C F A − 1 12 1 + D 0 t R 6 − 1 2 β ∇ C F C R B C . (19)

Equation 16 represents the probability distribution function of the interacting particle system at the short-time regime ³ ; it can be appreciated that the leading term, ρ 0 ( ξ , 0 , t ) ≡ exp − ξ 2 / ( 4 D 0 t ) / 4 π D 0 t D / 2 , is given by the Gaussian probability density valid for a very dilute system, while the sub-leading terms capture the corrections due to the curvature effects and interactions. One should notice that the joint probability distribution function has the same structure as the one developed in the asymptotic limit as t approaches zero for near points [35].

The expectation values of the observables can be calculated using the standard definition O ( ξ ) = ∫ M d D ξ G ρ ( ξ , 0 , t ) O ( ξ ) . Within the approximation given by Eq. 16, the expectation values can be estimated in the short-time regime using expectation values O ( ξ ) 0 with the leading term ρ ₀(ξ, 0, t); in other words, O ( ξ ) = O ( ξ ) 0 1 + τ ( 0 ) + τ B ( 1 ) ξ B O ( ξ ) 0 + τ B C ( 2 ) ξ B ξ C O ( ξ ) 0 + ⋯ . Expectation values of polynomial observables are particularly easy to compute due to the Gaussian structure of ρ 0 ξ , 0 , t .

We are interested in the calculation of the mean-squared geodesic displacement s 2 , where s = δ A B ξ A ξ B is the geodesic displacement in RNC. In addition, it is interesting to calculate the expectation value of the coordinate itself ξ ^B . For these expectation values, it is not a very difficult task to show by means of the standard calculation of the moments of a Brownian motion in a D -dimensional space that 1 0 = 1 , which is consistent with the normalization of the leading distribution ρ ₀(ξ, 0, t) also, it can be shown that the odd products are ξ A 1 ξ A 2 ⋯ ξ A 2 k − 1 0 = 0 , for any positive integer k, and for even products ξ B 0 = ξ B ξ 2 0 = ξ B ξ C ξ A 0 = 0 , ξ A ξ B 0 = 2 D 0 t G A B , and ξ A ξ B ξ 2 0 = 4 D + 2 D 0 t 2 G A B , where G ^BC is evaluated at ξ′. Since the previous approximation neglects the quadratic curvature effects that correspond to pre-factors of order D 0 t 3 in the mean-squared displacement [26], we only present the result up to order ( D 0 t ) 2 ; an alternative method for the calculation of the moments is given in Supplementary Appendix 4.2. This means that we basically neglect the linear terms of D ₀ t in τ B C ( 2 ) . Thus, the mean-squared displacement for the full N-particle system is given by s 2 = 2 D D 0 t − 2 3 R − 2 β ∇ A F A D 0 t 2 + ⋯ . (20)

One can notice in Eq. 20 that in the absence of the interaction term, that is, when F A = 0 , the mean squared displacement reduces to the previous result proportional to ( D 0 t ) 2 [26]. In addition, it is not difficult to elucidate that the subsequent correction of order ( D 0 t ) 3 involves pre-factors where curvature and interactions are coupled, for instance, terms proportional to R ∇ A F A and R B A ∇ C F A from the tensor τ B C ( 2 ) appeared as pre-factors; the cubic correction will be computed elsewhere in a future communication. In addition, note that similar terms appear in the expectation value of ξ _B , ξ B = 2 D 0 t τ B ( 1 ) .

One can notice that in the absence of the curvature, ξ B reduces to the well-known term β D 0 t F B , which establishes, on average, a preferential direction of the Brownian motion. In addition, this equation shows how the curvature is coupled to the interaction term within ( D 0 t ) 2 approximation. Finally, given an interacting force F _ij and specific sub-manifold M ⊂ R d + 1 , one can compute the mean-squared displacement for a tagged particle of the colloidal system by defining M S D ( t ) = 1 N s 2 , which is a quantity that can be easily calculated in dynamic-like simulations [21].

For practical purposes, we will provide an explicit expression for the mean-squared displacement of a tagged particle in the interacting system. We can utilize the fact that the big metric G _AB can be viewed as a block diagonal matrix, which enables us to express the Ricci curvature R of the big manifold M as the sum of the curvatures of the physical manifold M . Additionally, we will apply the explicit expressions for F A and the covariant derivative ∇ A to write ∇ A F A = ∑ i = 1 N ∇ α , i F i α . Moreover, using F i α = e α ( x i ) ⋅ ∑ j ≠ i F i j and the Weingarten–Gauss equation ∇ _α e ^α = −K n, one can calculate the expression straightforwardly as M S D t = 2 d D 0 t − 2 3 N ∑ i = 1 N R g x i − 2 β N ∑ i , j N e α x i ⋅ ∂ α , i F i j − K x i n x i ⋅ F i j D 0 t 2 + ⋯ , (21)

It can be observed that unlike the Brownian motion of a single particle [26], where the curvature effects are solely intrinsic, in an interacting particle system where the interaction happens in Euclidean space, extrinsic curvature terms are introduced by the mean curvature K in the geodesic mean-squared displacement. As an additional observation, we should comment that the aforementioned result (21) is consistent with the standard result for flat surfaces with R _g = K = 0 [47].

In the following paragraphs, we present two illustrative examples that will highlight the importance of the geometry on the particle diffusion.

Example 1: We now provide the estimation of the order of the curvature effect for a dilute hard-sphere gas on a spherical surface. According to [40], the equation of motion for the density of hard-sphere-like Brownian particles (that might include smooth inhomogeneities) at low concentrations, ϕ, is given by the effective diffusion equation ∂ρ/∂t = D∇² ρ, where the self-diffusion coefficient of gas of hard-spheres, D gas (at low concentrations), is given by D = D ₀(1 + α _∇ ϕ), where α _∇ ≃ 1.55 is a number obtained by considering hydrodynamic interactions (e.g., see chapter 6 of [40]). If one constrains the particles to diffusive exclusively on a curved surface, one expects, at least at the short-time regime, that one just needs to replace the Laplacian ∇² by the Laplace–Beltrami Δ _g operator associated with the surface. Thus, the geodesic-mean squared displacement is of the form s 2 = 4 D t − 2 3 R g ( D t ) 2 + ⋯ . Then, the short-time self-diffusion coefficient take the mathematical form D e ff = D 1 − 1 6 R g ( D t ) + ⋯ . Now, taking τ as the typical time that a particle diffuses a distance of the order of its own size, such as D ₀ τ ≃ σ ², where σ is the diameter of the particle, and a spherical surface of radius R, i.e., its curvature is R _g = 2/R ², one can immediately estimate the curvature effects as follows: 1 6 R g D τ ≃ σ 2 1 + α ∇ ϕ 3 R 2 . (22)

Let us take particularly the values R/σ = 10, α _∇ = 1.55 [40] and a packing fraction of ϕ ≃ 0.1; thus, one has 1 6 R g ( D τ ) ≃ 1 0 − 3 . In addition, if the diffusing distance ( ∼ D 0 τ ) increases, the curvature effects will be notorious.

Example 2: We now provide an illustrative example for interacting particles confined on a sphere S ². Let us first consider a generic force F _ij = f(r _ij )r _ij , where r _ij is the distance between the particles in the Euclidean space R 3 and r _ij is a vector from particle j to particle i; consequently F _ij is the force of particle j acting on particle i. Since particles are confined to the sphere, the vector r _ij = X _i − X _j , where X _i is the embedding function of the i − th particle. Since the manifold is the sphere, we have the following advantages X _i = R n(x _i ), where x _i = (θ _i , φ _i ) with θ _i and φ _i being the usual spherical coordinates for the i − th particle, and R is the radius of the sphere. Notably, the distance between the particles can be written as r i j = 2 R 1 − cos γ i j 1 2 , where cos γ i j = cos θ i ⁡ cos θ j + sin θ i ⁡ sin θ j ⁡ cos φ i − φ j . In addition, the curvatures R g ( x i ) = 2 R 2 and K ( x i ) = 2 R are independent from the coordinates; thus, the expression for the mean-squared displacement can be written as M S D t = 4 D 0 t − 4 3 D 0 t 2 R 2 − 2 β N ∑ i , j N 1 2 f ′ r i j r i j 1 + cos γ i j + 2 f r i j cos γ i j D 0 t 2 + ⋯ . (23)

Notably, β f ( r ) = − 1 r d d r β u ( r ) has units of inverse of square length, where u(r) is the pair potential. To be concrete, we use the soft potential u(r) defined by β u ( r ) = σ e / r 6 to compare with the previous results [31]. For this case, it is not difficult to see that f′(r)r = −8f. Then, one obtains the expression: M S D t = 4 D 0 t − 4 3 1 + 9 8 N σ e R 6 I N ; γ i j D 0 t 2 R 2 + ⋯ , (24) where I N ; γ i j is defined by I N ; γ i j ≔ ∑ i , j = 1 N 2 + cos γ i j 1 − cos γ i j 4 . (25)

In addition, note that the value of the ratio σ _e /R depends on the strength of the interaction potential. Now, we carry out a numerical bound of the finite sum (25) as follows: note that first I N ; γ i j is a positive number, second, we used cos γ _ij < 1, and third if σ/R is the minimal angle γ _ij between two colloids on the sphere, thus 1 − cos γ i j − 4 < 1 − cos σ R − 4 . Then, one should follow the following condition: I N ; γ i j < 3 N 2 1 − cos σ R − 4 . Thus, the contribution obtained from interactions obeys 9 8 N σ e R 6 I N , γ i j < 27 N 8 σ e R 6 1 1 − cos σ R 4 . (26)

Now, using the value N = 10, weak interaction σ e R = 1 0 − 2 , and R = 10σ (see, e.g., Ref. [31]) one obtains 9 8 N ( σ e R ) 6 I ( N , γ i j ) < 1 0 − 2 . Therefore, in this case, the curvature effects are visibly more pronounced.

The covariant form of Smoluchowski Eq. 4 opens up the possibility of developing a theoretical framework to study different interesting phenomena that cannot be understood with the standard statistical mechanics approximations based on a Euclidean formulation. For example, one of the topics that can be tackled with this approach is the initiation of the spinodal separation of particles interacting with short-ranged attractive forces and constrained to curved space in an analogy with the procedure presented by Jan Dhont in the case of Euclidean space [40]. Following these ideas, we need to convert Eq. 4 into an expression for the probability density of one particle instead of the joint probability of all the particles. To this end, it is necessary to perform a hierarchy of equations that allows us to marginalize the joint probability density function. Once the reduced Smoluchowski equation is obtained, it is necessary to take advantage of the short-range interactions to relate out-of-equilibrium phenomena with their counterparts in equilibrium. The connection between both cases, as usual, is made through approximations concerning the equilibrium values; at this point, there exists a wide range of ways to proceed. For instance, a perturbation approach can be combined using Riemann normal coordinate formalism, Monge’s parameterization, or covariant Fourier series to calculate all the relevant observables. On the other hand, a covariant Taylor expansion [48] approach can also be performed to compare the results with their flat counterparts [40].

In addition, the covariant formalism provided by Eq. 4 can be straightforwardly used to highlight the role of the geometry on the equilibrium equation of the state of colloidal dispersions embedded in a curved space, to elucidate the geometrical contributions during the onset of non-equilibrium states, such as gels and glasses, to study the dynamics of either passive or active colloidal particles on manifolds, and to investigate the curvature effects on the structural, kinetic, and phase transitions of attractive colloids, to mention a few examples of interest in the colloidal soft matter domain. As mere speculation and motivated by the recent contribution presented in [28], the formalism here presented can also be considered to study the dynamics of granular matter in curved manifolds.

Two aspects that need to be considered toward the extension of the covariant description of the Smoluchowski equation are the generalization of the fluctuation–dissipation theorem and the inclusion of the hydrodynamic interactions in the manifold. Here, we considered that the former is satisfied in Euclidean space and that the latter are completely neglected. However, both aspects are definitely crucial to account for the dynamical properties on the manifold since they are also intimately related to the onset of non-equilibrium states. Furthermore, this theoretical framework can be extended to include rotational motion in the manifold to consider cases where the particles are anisotropic in nature, i.e., anistropic particles, such as rods, and not all the orientations are independent in the case of board-like colloids or even in those physical situations where anisotropy is associated with the particle interaction, like in patchy colloids.

Finally, we should mention that the covariant compact form of Smoluchowski Eq. 8 allowed us to obtain an expression for the joint probability density function for the full system in the short-time regime. The method implemented can be extended to capture corrections of the order of ( D 0 t ) 3 . The short-time expression of the PDF can be used to give the curvature effects in the mean-squared displacement and the search role of the coupling between the curvature and the interactions; for instance, using this procedure, we can choose a specific interaction force and specific manifold M and give an estimation of the mean-squared geodesic displacement of a tagged particle of the colloidal system at short times. Moreover, the short-time expression of the PDF (16) can also serve to define a computational scheme to study the behavior of the full system using a Monte Carlo dynamics simulation that considers curvature effects. Additionally, the covariant compact form (8) allows us to formulate the N-particle system using a Feynman path integral representation, following the steps already implemented in [43]. Last but not the least, the study of some limiting cases of Eqs 4, 8 will also serve as a benchmark to computational or molecular simulation schemes adapted to study the behavior of colloids in non-Euclidean spaces.