Brownian diffusion in non-harmonic potentials

Abstract

Brownian motion in confined systems is widespread in soft matter physics, biophysics, statistical physics and related fields. In most of these systems, a Brownian particle cannot freely diffuse in the space but is confined by a potential well in a limited range of positions. When performing data analysis, typically the harmonic assumption is made, assuming that in the regions explored by the particle during its dynamics, the confining potential is fairly well described by a harmonic potential. This is however not valid a priori. In this work, it is shown how the diffusion coefficient and the potential width obtained through standard analysis underlying a harmonic approximation are affected by increasing errors when moving away from the conditions under which harmonic approximation is legitimate. These observations motivate the research of a more general method for properly obtaining the diffusion coefficient for a particle diffusing in a generic potential well. Here, a method is proposed that allows retrieving the correct diffusion coefficient by comparing the original data and ad hoc simulations without any a priori knowledge of the potential.

---

Introduction

The measurement of the mean squared displacement (MSD) of a Brownian particle is a useful tool to address the particle diffusivity and the viscous properties of a fluid at thermal energy scales. For a free particle, indeed, the slope of the MSD versus the lag-time τ gives directly the diffusion coefficient of the particle and ultimately the viscosity of the fluid in which the particle is moving. Such a technique has been largely used to find the viscosity in simple and complex fluids¹ and to probe the flow boundary conditions close to confining walls.² In the case of anisotropic particles, Brownian motion also allows to simply address the coupling between orientational and translational degrees of freedom.³ Unfortunately in most relevant situations, this simple method cannot be easily employed, as the particle is not free but has a dynamics confined by external fields. In this case, extracting the linear regime proportional to the diffusion coefficient is not always straightforward and lean on the experimentally accessible timescales. The examples of such situations include the movement of particles confined in an optical trap,⁴ by steric walls,⁵ in the magnetic field⁶ and by DLVO interactions.^7–9 Caged dynamics can also be observed in out-of-equilibrium crowded complex systems, as in the case of foams^10,11 and confluent cell monolayers.^12–14 In all these systems, the MSD increases at short time scales, while for longer timescales, it plateaus. The plateau is the signature of the presence of a confinement. In this situation, extracting the particle diffusivity† by the linear best fit of the MSD at short time scales can be tricky, as the linear part of the MSD is not always accessible for typical experimental sampling rates.^15–22 The typical solution for this problem is to fit the full MSD with an analytical fit function considering both the effect of the potential and the thermal diffusivity. Such an expression is available for harmonic potentials,^23,24 but is missing for a generic potential. These last cases are typically treated with an effective harmonic potential by assuming that the generic potential does not deviate too much from the harmonic one in the range of positions explored by the Brownian particle. Consequently, the MSD expression for a harmonic potential is typically used also for non-harmonic confinements.²⁵ However, to what extent this assumption can be considered valid has not been investigated.

In the present work, we propose a method to disentangle the viscous dominated behavior at short time scales from the potential dominated one at large time scales for a Brownian particle in a generic non-harmonic potential. We focus on two examples of potentials whose non-harmonic terms are tuned by control parameters. We computionally generate particle Brownian displacements and MSD in a fluid with given viscosity and external potential. The diffusion coefficient and plateau values, retrieved from the best fit of the full MSD assuming a harmonic potential, are then compared with the nominal ones to quantify the error underlying the harmonic assumption during the analysis. Building up on these results, a general method for obtaining the appropriate diffusion coefficient independent of the specific expression of the confining potential is proposed.

It must be noted that in some real systems – especially in biological ones – anomalous diffusion is present.^{18,19,26–28} The present work focuses on the simpler case of normal diffusion, but the proposed method can be extended to a wider plethora of phenomena, upon a proper adaptation of the simulation.

Methods

In typical diffusion experiments, the information on the dynamics is accessed through the MSD analysis. For a given trajectory and considering a given lag time τ, the displacement occurring along the direction x at time t during the interval τ is Δx(t,τ) = x(t + τ) − x(t). The MSD at lag time τ is then given by the time average 〈Δx(t,τ)²〉_t.

Noting U(x) as the conservative potential acting on the particle and as the stochastic force, where W(t) is the white noise, the Langevin equation is written as:

where m is the mass of the particle, v is its velocity along x, F(x) = dU(x)/dx and is the drag coefficient, with D being the diffusion coefficient. In the overdamped limit, typical in soft matter systems, eqn (1) results in

At first, the case of a harmonic potential U(x) = mω₀²(x − x₀)²/2 is considered, where x₀ and ω₀ are the equilibrium position and the characteristic frequency of the harmonic potential, respectively. In this case, an analytical expression for the MSD can be found^23,29 as:

where ω_μ = ξ/2m and .

When ω₀ ≪ ω_μ, the MSD results in the simpler expression:

Eqn (4) is commonly employed as the fitting function of the MSD in order to obtain the drag coefficient of particles confined in harmonic potentials, as for example in an optical trap.²⁹ Please note that in the limit of τ → 0, eqn (4) reduces to 〈Δx(t,τ)²〉_t = 2Dτ, as for a free Brownian motion. Meanwhile at large time lag τ → ∞, 〈Δx(t,τ)²〉_t reaches a plateau value given by the equipartion theorem: .

For general potentials, there is usually no analytical expressions for the MSD and numerical simulations are required. In order to fit an MSD in a generic potential, we first need to numerically generate the MSD data. To this scope, we generalized the procedure of Volpe and Volpe³⁰ developed for harmonic potentials. From eqn (2) a finite differential equation can be obtained in the form of:

where the discrete index i runs from 1 to the trajectory length N and δt is the time interval used as the simulation step. w_i is a Gaussian random number with zero mean and unit variance.‡ The underlying assumption is that the force F acting on the colloid does not change significantly over the time ranging from discrete time points i − 1 and i. In order to satisfy this condition, the simulation time step δt must be chosen sufficiently small. Eqn (5) is then used to simulate a trajectory of N points. From the obtained particle positions x_i at times iδt, the discrete squared displacement at time τ = nδt is calculated as Δx_i,n² = (x_i+n − x_i)² and then averaged over all the N − n points to obtain the MSD 〈Δx_i,n²〉_N−n.

In order to be as general as possible, the physical quantities have been made non-dimensional by normalizing over the timestep δt and the lengthscale [small script l] as follows: U′ = U/k_BT, x′ = x/[small script l], D′ = (Dδt)/[small script l]², and . For notational simplicity from here on, we redefine the non-dimensioned quantities without the prime: U′ → U, x′ → x, D′ → D and .

The force F(x_i) in eqn (5) is given by the x-derivative of a conservative potential U. In the present work, two sets of non-harmonic potentials have been chosen – one symmetric and the other non-symmetric with respect to their minimum at x = x₀. A useful parameter to systematically study the effect of non-harmonicity upon the simulated MSD is the half width of the potential h_w at a value U_d = 1 (k_BT in dimensioned units).

By defining

x₋ = x(U = U_d, x < 0) (6)

x₊ = x(U = U_d, x > 0) (7)

the half-width at U = U_d is given by:

At first, it is considered the case of a symmetric potential. To this end, a polynomial with even exponents is taken having the form:

where ρ is an non-harmonicity parameter ranging from 0 to 1 that establishes the relative weight of the non-harmonic terms compared to the harmonic one: for ρ = 0 the potential is harmonic, while for ρ = 1 the harmonic term is missing. Also, in eqn (9), n_p is an integer number indicating the number of even exponents larger than 2 considered in the polynomial series. In the present work, n_p is systematically chosen equal to 9 (corresponding to a polynomial degree 20), as no significant changes in the results have been observed for higher values of n_p. The spatial dependence is chosen so that x₀ = 0, U(0) = 0 and the half potential width at U = U_d is equal to h_w.

In Fig. 1a are shown examples of the considered non-harmonic potentials for different values of ρ.

Fig. 1 Non-harmonic (a) and asymmetric (b) potentials as defined in the text for different values of ρ and α, respectively, as a function of the non-dimensioned space coordinate x.

Although considering higher order symmetric terms is the first logical step for exploring the effect of non-harmonic potentials, in most real cases, the potential is also non-symmetric. For such types of potential, one can define an asymmetric parameter α in order to quantify the degree of asymmetry:

α is zero for symmetric potential wells and increases for increasing asymmetry between the left and the right sides of the potential with respect to the equilibrium position.

A typical example of a non-symmetric potential in colloidal dynamics experiments is the combination of van der Waals (or DLVO) and gravitational potentials. This is the case of a colloid sedimented near a wall or an interface.⁷ A simplified version of such potential can be obtained adding an hyperbolic repulsive term (van der Waals) to a linear term (gravity), in the form of:

defined for x > x_d. Such potential allows the independent tuning of h_w and α while keeping simple the algebras. By imposing the minimum of U at x = 0 and U(0) = 0, the expression becomes:with . By making explicit h_w and α according to their definition in eqn (8) and (10), respectively, A and B result in:

Examples of the asymmetric potentials obtained with the expression in 11 are reported in Fig. 1b for α ranging from 0.02 to 1.9.

Results

Using the simulation procedure described in the previous section, about 1.6 × 10⁵ Brownian independent trajectories (each one made of 2 × 10⁶·h_w points) have been simulated using the potential functions given in eqn (9) and (12). Simulations have been made for different values of potential well width h_w, non-harmonic parameter ρ and potential asymmetry α. It must be noted that for avoiding divergence at x = x_d in simulations with the asymmetric potential, the left-side of U(x) was approximated for U > 5 with the tangent to the potential at U = 5.

The simulation parameters and other relevant quantities are reported in Table 1.

Table 1 Summary of the simulation parameters and other relevant quantities

D = 1.76	U d = 1
h w = 1	n p = 9
ρ ∈ [0,1]	α ∈ [0,2]
n. points per simulation 2 × 10^6hw

---

For each generated trajectory, the MSD was then computed and fitted with the analytical MSD expression for a harmonic potential (eqn (4)) to obtain the diffusion coefficient D_fit and the plateau value h_fit². For the MSD, only one point every 2000·h_w was considered. Thus, the MSD time step is dt_MSD = 2000·h_wδt. This is made to reduce the computational time and to be closer to the sampling rate of real experiments, where the linear part of the confined MSD typically does not exceed the first 2–3 experimental points. The obtained D_fit and h_fit are then compared with the simulation inputs D and h_w in order to assess the validity of the harmonic approximation in the MSD analysis. For the sake of comparison, the relative errors e_D = (D_fit − D)/D in the diffusion coefficient and e_h = (h_fit − h_w)/h_w in the potential width related to the MSD plateau have then be computed.

Please note that such a comparison is independent of the chosen value of the input diffusion coefficient in the limit of relatively slow diffusion when the force can safely be assumed constant between two successive time steps.

As an output of the analysis, it was also noticed that the discrepancies are also independent of the value of h_w, as it can be seen in Fig. 2. Simulations have been performed for different values of h_w in the range 0.25–25. In Fig. 2, e_D and e_h are reported as a function of h_w for different values of ρ for simulations using the symmetric non-harmonic potential. Points are the averages over different simulations, while the error bars indicate the standard deviation and the standard deviation of the mean. Within a given potential and for a given control parameter ρ, they are all equal within the statistical incertitude. Similar results are obtained for the asymmetric potential for each considered value of α. Consequently, in the following discussion, we focus on the dependence of e_D and e_h on ρ and α.

Fig. 2 Values of e_h (a and b) and e_D (c and d) as a function of h_w obtained by fitting trajectories simulated from eqn (9). In the shown data, ρ is equal to 10⁻⁵ (a and c) and to 1 (b and d). Points are averages over 35 different simulations, while the error bars indicate the standard deviation and the standard deviation of the mean.

In Fig. 3a–c are shown the main results relative to the simulations within the non-harmonic symmetric potential described by eqn (9). In Fig. 3a are shown the MSDs obtained for h_w = 1 for increasing values of the non-harmonic parameter ρ. As expected, for low values of ρ the MSD is close to the one expected for a particle moving within a harmonic potential (black crosses). As ρ increases, both the plateau values and the diffusion coefficients decrease. This trend can be better visualized in Fig. 3b and c, where e_h and e_D are respectively reported as a function of ρ. The black lines are averages over 35 simulations, while the gray shadow region represents the corresponding standard deviation of the mean. The MSD fit undervalues the diffusion coefficient up to 8% as ρ increases from 0 to 1. It may be surprising that the slope of the MSD changes with ρ if the input D is the same for all simulations, but it must be pointed out that in the present data, the linear regime of the MSD holds only for very low values of τ/dt_MSD. For larger time scales, the effect of the potential is strong enough to affect the slope of the MSD (i.e., the dynamics is already subdiffusive). This is the case for the intermediate region of the curve shown in the inset of Fig. 3a. The deviation of D_fit from D is therefore originated by the wrong assumption of a harmonic potential underlying eqn (4) used for the fit.

Fig. 3 (a and d) Normalized mean squared displacement obtained from different simulations using the unharmonic (a) and asymmetric (d) potentials. Different colors represent different values of ρ (a) and α (d) according to the legend. On the x-axis it is reported that the delay time scaled over dt_MSD = 2000·h_wdt (here, h_w = 1). Values of e_h (b and e) and e_D (c and f) as a function of ρ (b and c) and α (e and f) obtained from simulations using the non-harmonic and asymmetric potentials, respectively. Each black curve is the average over 35 different simulations. The shadow region represent the standard deviation of the mean.

Concerning the potential width, it can be seen in Fig. 3b that the MSD plateaus at lower values when ρ increases, in spite of the fact that all the potentials have been built with the same half-width h_w at k_BT. For the upper limit of ρ, this discrepancy exceeds 20%.

The MSD's plateau value represents the square of the maximum displacement the particle explores on average inside the potential well. In order to rationalize why for the same potential width the presence of non-harmonic symmetric terms reduces h_w, it is possible to resort to a simple qualitative argument. As can be seen in eqn (5), the instantaneous displacement of the particle depends on the resultant force given by the sum of the stochastic force plus the conservative one. If the simulated particle at a given time step is located in a position x, in the following simulation step, it can move further away from the equilibrium position only if the stochastic noise is larger than the absolute value of the conservative force in that point. In this view, the limit of the displacements is therefore given by the comparison between F(x) and the width of the stochastic term distribution. In other words, for a given stochastic noise, the maximum displacement from the minimum depends on the local slope of the potential. This has been qualitatively tested by computing for a non-harmonic potential the distance Δx between the position of the potential minimum and the point where the slope of the potential is the same as one of the harmonic potential at k_BT (Δx = h_w in the harmonic case). In Fig. 4a, the relative difference between Δx and h_w has been plotted versus ρ. There it can be seen that indeed the trend with ρ is similar to the one of e_h. For increasing ρ, therefore, the distance from the equilibrium that the particle can reach before reaching an overwhelming recall force decreases, thus resulting in a lower plateau value in the MSD, as observed in Fig. 3b.

Fig. 4 Relative difference between Δx and h_w computed as a function of ρ (a) and α (b) using eqn (9) and (11), respectively. Here Δx is defined as the average distance (left and right) from the potential minimum to the coordinate where the slope of the potential is equal to the slope of the harmonic potential at k_BT.

The slow increase that can be seen in e_h for ρ approaching 1 can be understood considering the limit case of a box potential. There, the potential is zero in the range (−h_w, h_w) and becomes infinite for x < −h_w or x > h_w. The particle is therefore expected to diffuse as a free particle in the range −h_w < x < h_w and the MSD should plateau at hw², as in the harmonic case. Consequently, we expect that by approaching the box potential at the largest values of ρ, the difference e_h with respect to the harmonic case has to reduce as observed in Fig. 3b.

The results of the simulations with the asymmetric potential can be seen in Fig. 3d–f. As for the non-harmonic potential, in Fig. 3d, some examples of the MSD are reported, in this case for increasing α. Good agreement is found between the MSD simulated at low α and the one for the harmonic potential (black crosses). As α increases, both the fitted diffusion coefficient and plateau value deviate from those of the harmonic potential, but with some differences in the trend. Contrary to the case of the non-harmonic potential, now h_fit is systematically larger than the potential width (Fig. 3e), signaling that by increasing the asymmetry, the particle thermally explores a larger region of the potential. As for the non-harmonic potential, the sign of the deviation is the same as that of the relative difference between Δx and h_w as a function of α, as shown in Fig. 4b, thus highlighting a similar mechanism. In this case, Δx is always larger than h_w. As α increases, indeed, the coordinate on the right-side of the potential where the slope reaches that of the parabolic potential at k_BT moves away from the equilibrium position faster than how the corresponding coordinate on the left-side approaches the equilibrium.

Concerning the diffusion coefficient, e_D decreases as the potential deviates from the harmonic case for low values of α. In this region, the trend is therefore the same as for the non-harmonic potential: an underestimation of the real diffusion coefficient. The amount of the deviation is greater but comparable to the one of the non-harmonic case.

In general, the deviation of the diffusion coefficient in the case of a non harmonic potential may therefore be significant and potential-dependent, thus making questionable the use of a harmonic potential to fit the MSD data. This raises the question as how to properly extract a diffusion coefficient from experimental data of Browinian diffusion of a confined particle in the most general case, when the analytical expression for the MSD is missing.

In the following discussion, a method is proposed to overcome such difficulties which is based on the combination of experimental data analysis and ad hoc simulations, named the iterated simulation (IS) method.

The iterated simulation (IS) method for retrieving diffusivity in a generic confinement

To illustrate the proposed method, here a simulated trajectory x(t) has been chosen with known potential U_k(x) and diffusion coefficient D_k. The chosen potential is asymmetric. In the following discussion, such simulation will be referred to as the original dataset, playing the role of experimental data from which we want to extract the diffusion coefficient. Consequently, both the potential and the diffusion coefficient have to be recovered only knowing the trajectories of the original dataset and then compared with U_k(x) and D_k to test the efficiency of the proposed method.

From the particle positions x, a probability distribution p(x) is calculated, as shown in Fig. 5a. From the probability distribution, the confining potential is retrieved using the Boltzmann equation:^7,32

where x₀ is the equilibrium position that corresponds to the x − coordinate of the p(x) maximum and U₀ is the value of the potential in x₀. After a proper smoothing of the potential, performed with a moving average, the corresponding force can be evaluated as F(x) = (U_x+dx − U_x−dx)/2dt, where dx is the bin size.

Fig. 5 Simulated datasets used for testing the IS method. In the simulations, α (from red to blue): 0.002, 0.8, 1.2, 1.7. h_w = 1 nm (0.5). D = 1.76. (a) Normalized histogram of the positions for each trajectory. (b–e) Potentials obtained from (a) using the Boltzmann relation (continuous lines) compared with the potential used for the dataset simulation (green dashed lines) for increasing asymmetry.

In Fig. 5b–e are shown potential profiles obtained from the probability distributions shown in Fig. 5a. Green dashed lines represent the analytical function of the potential given as an input to numerically calculate the trajectories used in the p(x), as shown in Fig. 5a. A good agreement can be observed between the input potential and the one obtained from the statistical analysis of the original datasets, except for the higher values of x due to the scarce statistics.

The force obtained from the probability distribution is then used to replace the analytical expression for the force in eqn (5) and used to simulate a series j of simulations having the same timestep of x(t) and different diffusion coefficients D_in varying within a reasonable range of values.

In order to retrieve the input diffusion coefficient D_k, we fit the MSD of the original dataset x(t) with eqn (4), thus obtaining D_fit, as in standard analyses. Similarly, each MSD obtained from the set of simulations for different D_in is fitted with the same expression, obtaining a different value of for each D_in.

The correct diffusion coefficient is found when both fits return the same value of the diffusion coefficient in the range of the error bars. Practically, because of the randomly generated white noise and the finite number of points per simulation, simulations performed with identical parameters yield different MSDs. To gain in stability, a number of simulations are thus generated for a given D_in and fitted with eqn (4) until the standard deviation of the obtained values of the diffusion coefficient stabilizes at a given value σ_D. The average for a given D_in is then compared with D_fit obtained by fitting the original dataset. When , the corresponding D_in is taken as the effective diffusion coefficient D_out resulting from the IS method. A summary of the IS method is shown in Fig. 6.

Fig. 6 Sketch summarizing the procedure to measure the diffusion coefficient for a particle diffusing in a generic confinement potential, as described in the text.

In order to test the efficiency of the IS method, some examples are reported in Fig. 7 for different values of α. Each plot refers to one of the datasets shown in Fig. 5. Among all the possible simulations, for each α the original dataset has been chosen so that they have a discrepancy between D_k and D_fit corresponding to the average one as depicted in Fig. 3. Each blue point is the average made over 500 independent simulations with a given D_in. The continuous and the dashed blue lines are respectively the linear fit of the different and the input D_in. Once again, the discrepancy can be seen between the input and the fitted diffusion coefficients, negligible for small α but increasing for larger values of the asymmetry parameter. Similarly, red continuous and dashed lines represent the fitted D_fit and the diffusion coefficients D_k used to build the original dataset. In the same figures, it is also illustrated how D_out is defined: the coordinate where the linear fit of intersects D_fit. For helping with the comparison between D_k, D_fit, and D_out, lines transfer the intersect coordinate on the y-axis (green line). It can be seen that for low values of α, both D_fit and D_out coincide with D_k. For increasing values of the asymmetry, D_fit underestimates more and more D_k, while D_out shows a better estimation of D_k. The corresponding numerical values are reported in Table 2: for the reported cases, the maximum deviation of the obtained diffusion coefficient from the value used to build the trajectory is of the order of 2%, while the error made with the standard fit of the MSDs is systematic and reaches up to 10%. The IS method used here on simulated data sets can thus be reliably employed for measuring the diffusion coefficient of experimental data for particles diffusing within a generic potential.

Fig. 7 Results of the process for obtaining D_out using the simulation procedure described in the text. The original datasets are simulations made using the potential in eqn (11) with α equal to 0.002 (a), 0.8 (b), 1.2 (c) and 1.7 (d). Blue points represent the obtained averaging the diffusion coefficient fitted over 500 simulations. Error bars represent the standard deviation σ_D. The continuous and the dashed blue lines are the linear fit of the scattered points and D_in, respectively. The continuous and dashed red lines refer to the original dataset and represent D_fit and D_k, respectively. The green continuous lines represent D_out.

Table 2 Comparison of the diffusion coefficients obtained using MSD fitting and IS methods with the ones used to build the data for different values of α. The numbers correspond to the data shown in Fig. 7. In the columns, D_fit and D_out are reported as the obtained diffusion coefficients and their percentage deviations from D_k

α	D k	D fit	D out
0.002	1.76	1.76, 0%	1.76, 0%
0.8	1.76	1.66, −6%	1.73, −2%
1.2	1.76	1.62, −8%	1.74, −1%
1.7	1.76	1.59, −10%	1.78, 1%

---

Although this method is applied here to the one-dimensional case, a generalization to the multidimensional case – 2D and 3D, also considering rotational degrees of freedom – is possible. A thorough examination of these cases is beyond the scope of this paper. In this section, we will provide a concise overview of the primary features of the transition to the multi-dimensional case. The extension to a multidimensional case results in an increase in the number of differential equations and the introduction of coupling terms between them. The coupling terms may arise from the non-separability of the potential and from the roto-translational coupling. In order to generalize the IS method, the potential must first be obtained as for the 1D case, from the multidimensional histogram of positions. It is then possible to recover the different components of the force by means of space derivatives of the given potential. Subsequently, the set of coupled differential equations (one for each degree of freedom) can be obtained by utilizing the recovered force. The coupling between different degrees of freedom (e.g. roto-translation) should be eventually taken into account. The differential equations can then be used for implementing numerical simulations. Finally, the MSDs can be obtained from each degree of freedom from the original and simulated datasets and compared, as previously described, in order to get the diffusivity value for each degree of freedom.

Conclusions

The present work has pointed out the limitation of the MSD fitting method for retrieving the potential well width and the diffusion coefficient for Brownian trajectories confined in non-harmonic potentials. Two types of potentials have been considered: a symmetric but non-harmonic one and an asymmetric one. By comparing the parameters used to build an input numerical trajectory with the best-fit outputs of the MSD obtained from simulations, it has been shown how the commonly assumed equivalence between the MSD plateau value and the potential square width at U = k_BT breaks down for nonharmonic potentials. In addition, an incorrect value of the diffusion coefficient is also found. Building up on this, the simulation framework has been used to devise a method able to correctly evaluate the diffusion coefficient without any a priori knowledge of the confining potential. The method is based on the comparison of the MSD of the relevant dataset with ad hoc simulations. This method can be used to check if the effect of non-harmonicity is relevant in specific cases and to increase the precision of the obtained diffusion coefficient by removing the systematic error introduced by non-harmonicity. This approach is particularly valued for situations where a precise measurement of the diffusion coefficient is important.² Moreover, it can in principle be generalized to the 2D or 3D case, provided that attention is paid in adding the proper coupling terms to the differential equations that will replace eqn (5). The present work focused on tracking-based methods for quantifying the dynamics. However, other techniques can also be used to estimate the MSD of a single particle or an ensemble of particles (e.g., FCS,³³ DWS,³⁴ and DDM³⁵). In principle, the proposed IS methodology could be applied in conjunction with experimental techniques that do not rely on single-particle trajectory reconstruction. However, with tracking-free methods, the direct access to the potential measurement is lost. Alternative ways for accessing (or knowing) the potential would therefore be needed. Provided this, and for symmetric potentials (radial in the 2D case), the presented analysis can be extended to tracking-free methods.

Even though the relevant cases of anomalous diffusion have not been considered in the present work, the same working principle can per se be adapted to such rich systems. For this, simulations should, however, be adapted case by case for modeling super- or sub-diffusive behaviors. The present work only considers the passive case. The possibility of generalizing the same arguments and method to the active case depends on whether it is allowed to consider the activity as an effective temperature. This can be verified when the active motion is isotropic, sufficiently homogeneous in time and space, and with a persistent length smaller than the size of the well.³⁶

Recently, the IS method has been applied to the study of the viscous drag of spherical²¹ and ellipsoidal²² colloids in the vicinity of an air–water interface, where the superposition of a gravity potential and DLVO interactions results in an asymmetric non-harmonic potential.⁷ There, the method has made it possible to verify the validity of the harmonic approximation in the studied cases. We believe that the proposed methodology is not only limited to the colloidal field but can potentially be employed in any confined Brownian dynamics when a precise evaluation of diffusivity is required.

Author contributions

MN and SV conceived the project. SV worked on the simulations, the data analysis and the original draft of the paper. All authors contributed to the writing (review and editing).

Conflicts of interest

There are no conflicts to declare.

Data availability

The code implementing the IS method for retrieving diffusivity in a generic confinement can be found at https://github.com/stevilla1/ISmethod with DOI identifier https://doi.org/10.5281/zenodo.15373268. The version of the code employed for this study is version v1.0.0.

Acknowledgements

Open Access funding provided by the Max Planck Society.

Notes and references

1. T. Mason, K. Ganesan, J. H. J. H. van Zanten, D. Wirtz and S. Kuo, Phys. Rev. Lett., 1997, 79, 3282.
2. S. Villa, G. Boniello, A. Stocco and M. Nobili, Adv. Colloid Interface Sci., 2020, 284, 102262.
3. Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky and A. G. Yodh, Science, 2006, 314, 626–630.
4. K. Berg-Sørensen and H. Flyvbjerg, Rev. Sci. Instrum., 2004, 75, 594–612.
5. S. Ghosh, D. Wijnperlé, F. Mugele and M. Duits, Soft Matter, 2016, 12, 1621–1630.
6. R. Bubeck, C. Bechinger, S. Neser and P. Leiderer, Phys. Rev. Lett., 1999, 82, 3364.
7. S. Villa, A. Stocco, C. Blanc and M. Nobili, Soft Matter, 2020, 16, 960–969.
8. J. A. Rivera-Morán, Y. Liu, S. Monter, C.-P. Hsu, P. Ruckdeschel, M. Retsch, M. Lisicki and P. R. Lang, Soft Matter, 2021, 17, 10301–10311.
9. S. L. Eichmann, S. G. Anekal and M. A. Bevan, Langmuir, 2008, 24, 714–721.
10. A. Gittings and D. J. Durian, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 066313.
11. F. Giavazzi, V. Trappe and R. Cerbino, J. Phys.: Condens. Matter, 2020, 33, 024002.
12. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
13. R. Cerbino, S. Villa, A. Palamidessi, E. Frittoli, G. Scita and F. Giavazzi, Soft Matter, 2021, 17, 3550–3559.
14. S. Villa, A. Palamidessi, E. Frittoli, G. Scita, R. Cerbino and F. Giavazzi, Eur. Phys. J. E: Soft Matter Biol. Phys., 2022, 45, 50.
15. G. Boniello, J. Malinge, C. Tribet, E. Marie and D. Zanchi, Colloids Surf., A, 2017, 532, 510–515.
16. R. Sarfati, C. P. Calderon and D. K. Schwartz, ACS Nano, 2021, 15, 7392–7398.
17. R. J. Oetama and J. Y. Walz, J. Colloid Interface Sci., 2005, 284, 323–331.
18. 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., 2004, 92, 178101.
19. M. Levin, G. Bel and Y. Roichman, J. Chem. Phys., 2021, 154, 14.
20. A. Birjiniuk, N. Billings, E. Nance, J. Hanes, K. Ribbeck and P. S. Doyle, New J. Phys., 2014, 16, 085014.
21. S. Villa, C. Blanc, A. Daddi-Moussa-Ider, A. Stocco and M. Nobili, J. Colloid Interface Sci., 2023, 629, 917–927.
22. S. Villa, D. Larobina, A. Stocco, C. Blanc, M. M. Villone, G. d'Avino and M. Nobili, Soft Matter, 2023, 19, 2646–2653.
23. M. C. Wang and G. E. Uhlenbeck, Rev. Mod. Phys., 1945, 17, 323.
24. P. N. Pusey and W. Van Megen, Phys. A, 1989, 157, 705–741.
25. G. Nisato, P. Hebraud, J.-P. Munch and S. Candau, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 61, 2879.
26. M. Schoen, J. H. Cushman and D. J. Diestler, Mol. Phys., 1994, 81, 475–490.
27. N. Gal and D. Weihs, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 020903.
28. I. Abelenda-Núñez, F. Ortega, R. G. Rubio and E. Guzmán, Small, 2023, 19, 2302115.
29. T. Li, S. Kheifets, D. Medellin and M. G. Raizen, Science, 2010, 328, 1673–1675.
30. G. Volpe and G. Volpe, Am. J. Phys., 2013, 81, 224–230.
31. B. Øksendal, Stochastic Differential Equations, Springer, 2003.
32. B. M. Alexander and D. C. Prieve, Langmuir, 1987, 3, 788–795.
33. J. Kubečka, F. Uhlk and P. Košovan, Soft Matter, 2016, 12, 3760–3769.
34. D. J. Pine, D. A. Weitz, P. M. Chaikin and E. Herbolzheimer, Photon Correlation Techniques and Applications, 1988, p. QONC35.
35. P. Edera, D. Bergamini, V. Trappe, F. Giavazzi and R. Cerbino, Phys. Rev. Mater., 2017, 1, 073804.
36. G. Volpe, S. Gigan and G. Volpe, Am. J. Phys., 2014, 82, 659–664.

---

Footnotes

† Here and in the following we only consider thermal passive systems, where particle diffusivity is properly defined. This can be in principle generalized to out-of-equilibrium systems where an effective diffusivity can be defined, but a case-by-case assessment must be carried out.

‡ For the reason why W(t) can be discretized as , see Volpe and Volpe³⁰ and Øksendal³¹ dissertations.

---