Buffon's Brownian needles: harnessing thermal motion for stochastic sampling

Abstract

We demonstrate a physical implementation of Monte Carlo sampling using the Brownian motion of microscopic rods, applied to the classical Buffon's needle experiment. In this way, a problem in geometric probability is mapped onto a Monte Carlo method, with a physical system performing key aspects of the computation. The experiment's parameters are embedded directly: the rods length encodes the probability integral, while their thermal motion supplies the sampling. Although only a toy-model system, this approach illustrates how embedding probabilistic structure into soft matter can provide a low-energy pathway for stochastic computation that exploits freely available thermal noise.

---

1 Introduction

According to the International Energy Agency (IEA) the global energy consumption of computing, including data centers, networks and individual devices make up a significant amount of the total energy consumption. Data centers alone are assumed to consume about 1–1.5% of the global electricity, trends strongly increasing.^1,2 To ensure energy security, not only new sustainable energy sources, also new, more efficient computing strategies are needed on all levels, from hardware and architecture, to software based approaches, but also new, emergent forms of computing are developing.³

Especially quantum-based computing has attracted scientific interest in general,^4,5 been awarded a very recent Nobel Prize and received a high volume of investment capital.^6,7 Alternative forms of computing based on biological^8,9 or chemical properties¹⁰ intrinsic to the systems have drawn attention over time. Colloidal and active matter systems are increasingly attracting interest for computing applications: colloidal solutions have been used in a number of publications as a physical reservoir for reservoir computing.^11–13 The first realisation of a physical reservoir computer using self-propelled active microparticles was based on a nonlinear dynamical system due to time delays in retarded interactions.¹⁴ Immense potential is also ascribed to the use of inherent material properties to embody intelligence.¹⁵

Randomness plays a central role in modern computing. From secure cryptographic protocols and gambling systems to simulations and stochastic algorithms, the need for high-quality random number generators (RNGs) has only grown with the increasing complexity of computational tasks.^16–19 However, generating randomness is not without cost. Whether it is achieved via deterministic, classical computer-based pseudo-RNGs; or specialized, hardware-based true-RNGs, producing randomness in digital systems incurs significant energy costs – through both computational cycles and memory operations.^20,21

A notable application of RNGs is found in Monte Carlo methods – a class of stochastic sampling techniques widely used to solve problems in numerical integration, physics, finance, and machine-learning.^22–24 Monte Carlo methods rely, fundamentally, on the generation of random sampling – the quality and quantity of which can significantly impact the accuracy and efficiency of the computation. As digital technology continues to progress, especially in the field of machine-learning, the energy cost of computational tasks is becoming highly scrutinized and alternative paradigms of computing are now under exploration.^25,26 Within this scope, alternative, low-energy forms of random number generation demand investigation.

Thermal noise – the microscopic energy fluctuations caused by the thermal motion of particles – is a ubiquitous and naturally occurring form of randomness. It is a direct consequence of the statistical nature of many-bodied, thermodynamic systems in which the kinetic energies of bodies are distributed following the Boltzmann distribution.²⁷ In classical computing, where results should be deterministic and non-probabilistic, thermal fluctuations are suppressed by averaging results out. Conversely, in quantum computing they are suppressed by holding systems at low temperature in order to maintain quantum coherence.²⁸ Given that thermal noise is freely available at ambient temperature, a question arises: instead of suppression, can thermal fluctuations be leveraged as a computational resource? This question has motivated research into exploiting thermal noise for stochastic computing.^29–31

An overtly physical manifestation of ambient thermal fluctuations is Brownian motion. First described by Robert Brown in 1827 and explained by Einstein in 1905,† Brownian motion describes the erratic movement of microscopic particles suspended in a fluid, whose rheological properties and density of the surrounding fluid influence the observed phenomenon. The erratic motion of the micronscale objects is driven by collisions with the surrounding molecules, which in turn are driven by ambient thermal noise.^33,34 Importantly, while the underlying physics is deterministic, the macroscopic behaviour is effectively unpredictable and thus serves as a potent physical source of randomness.³⁵

In this work, we present an experimental system of microscopic rods in an aqueous medium which exhibit Brownian motion and diffuse over a surface digitally patterned with evenly spaced lines, creating a dynamic realisation of the classical Buffon needle experiment. In the experiment, the probability that a randomly dropped needle intersects a set of parallel lines depends on the geometry of the system and the value of π (Fig. 1).^36,37 Buffon demonstrated mathematically that by knowing both the rods length and the separation of the parallel lines, π could be estimated by simply counting the number of scattered needles and the number of needle-line crossings. Specifically, by setting l = 2d, then the probability of crossing is 1/π. Thus, π can be estimated by counting the total number of needles and dividing by the number of crossings. The accuracy of the estimation can be increased by re-scattering the needles and increasing the sample count, N. First posed as a question in geometric probability, it embodies the same statistical principles of Monte Carlo methods: drawing repeated random samples from a well-defined probability space, recording the occurrence of a specific event, and using the observed event frequency to estimate an underlying constant or integral.^38,39 In our system, re-scattering of the needles, i.e., resampling, is realised entirely by thermal fluctuations in a fluidic environment, driving the rods into new, randomised positions and angles. Through this conceptual toy model, we hope to demonstrate Brownian motion as a naturally occurring, zero-energy-input source of stochasticity required for solving probability integrals.

Fig. 1 (A) Schematic of Buffon's needle experiment. Rods of uniform length are randomly scattered on a surface of parallel lines with uniform separation. (B) Probability of rod crossing a line depends on rod length, l, line separation, d, and π. Thus, π can be estimated by counting the number of rods and crossings. In (A), the estimate would be [small pi, Greek, circumflex] = 9/3 = 3.

Similar uses of thermodynamic noise to generate ‘high-quality’ random numbers have been demonstrated previously,^40–42 and Brownian motion of colloids for RNG has been demonstrated by others.^35,43 We do not aim to compete with state-of-the-art physical RNGs for the quality or frequency of random numbers generated. Rather, we demonstrate a conceptual proof-of-principle by using inherent properties of the physical systems to directly solve Monte Carlo integrals. Specifically, the rod-shape geometry of the particle naturally embodies variables of the Buffon needle experiment.

2 Results and discussion

2.1 Rod selection and experimental set-up

To reproduce Buffon's needle experiment on the microscale, the first step is the selection of rod-like particles which exhibit both translational and angular Brownian motion in order to re-scatter themselves. These motions are the engine of the sampling within the Monte Carlo framework, continuously randomising without any external input. Three species of rod-like particle, spanning a range of sizes, were measured. These were: 1 μm long, 0.2 μm diameter SiO₂ rods; 3 μm long, 0.5 μm diameter ZnO rods; and 15 μm long, 2.5 μm diameter SiO₂ rods. A measure for the degree of Brownian motion can be obtained from the diffusion coefficient of a particle, which in turn is strongly dependent on both, the colloidal properties, but also the surrounding fluid. As size decreases, both the angular and translational diffusion coefficients of the rods increases according to the Stokes–Einstein equation as well as equations describing rods derived by Broersma.^44,45 With physical stochastic systems, the timescales of the processes are inherent to the components of the system. Thus, to achieve faster ‘computation’ rates, smaller systems which display faster kinetics are desirable. On the other hand, larger particles can be more easily resolved in optical microscopy and so can provide higher resolution data.

Using adapted Broersma relations^46,47 we estimate decorrelation times τ from the rotational diffusion coefficients D_rot for rods with length, l, and diameter, w.

τ_1,0.2 ∼ 0.1 s,

τ_3,0.5 ∼ 2 s,

τ_15,2.5 ∼ 300 s

These estimations indicate clearly the importance of particle size on resampling rate, with decorrelation times spanning 4 orders of magnitude. They indicate that larger SiO₂ rods would be unsuitable for this experimental work, let alone a reasonable random-number generator.

After dispersing the microrods in a dilute surfactant solution on a plasma-cleaned glass slide to minimise aggregation and sticking to the surface, the rods behaviour was observed by optical microscopy. It was found that both the large and small SiO₂ rods were unsuitable for analysis. The large rods, as predicted, were measured showing limited diffusion (D ∼ 10⁻¹⁴ m² s⁻¹). Conversely, the small SiO₂ rods were highly dynamic, translating more than a body length on the order of 10 s of milliseconds. However, they had a sufficiently low buoyant mass to show off-plane angular diffusion, resulting in them coming in and out of focus and with a changing projected shape (SI). This meant they could not be reliably tracked and were unsuitable for the strictly two-dimensional Buffon needle experiment. Balancing the need for resolvable planar motion with a sufficient reorientation rate, we selected 3 μm × 0.5 μm ZnO rods as the optimal compromise for further experiments. The ZnO rods showed clear Brownian motion both translationally and rotationally, and with all rotation being planar they were suitable for performing Buffon's test.

Measurements of their angle over time display a decorrelation time (autocorrelation function (ACF) < 1/e) of 0.81 s with a deviation of ±0.34 s (Fig. 2C) – slightly lower than the Broersma-estimated value which could be due to electrostatic effects or low-level photocatalytic activity of the rods, both of which are neglected in the Broersma estimation. Measurements of the mean-square-displacements of the rods provide a measure for the translational diffusion coefficient: D_trans = 16.9 μm² s⁻¹. This means within 1 s each rod can reposition a body length (3 μm) away from its initial position. Thus, with every 2 s it can be assumed that the angle and position of the rods has been sufficiently randomised, or re-scattered as with the analogy to Buffon's needle.

Fig. 2 (A) Schematic showing reorientation (Δθ) and translation (Δs) of dynamic rods after a time Δt. (B) Optical and electron microscopy images of ZnO rods sedimented on the surface. (C) Plot of the mean autocorrelation function of the angles (ACF) of rods over time (blue) and angles of rods over time (red). Decorrelation time is measured as the lag time at which the ACF drops below 1/e, denoted with the dashed black line and the upper standard deviation denoted with the grey dashed line.

ZnO rods were dispersed as described and images were recorded at 2 s intervals for just over 17 h (31 058 frames total). The resulting video was binarised and processed using a custom OpenCV platform (SI). A set of parallel lines, separated by a distance 2× the median rod length, was digitally imposed on the images and in each frame every rod as well as every rod-line crossing was counted. For each frame, a π estimate was measured from just rods and crossings in the frame as well as a π estimate from the cumulative rods and crossings.

2.2 From Brownian dynamics to Pi

Fig. 3A shows the cumulative estimate of π, defined as the ratio of the running-total numbers of rods and crossings. Superimposed on this are the single-frame estimates (plotted as red crosses) which exhibit significant variability. These individual estimates range from as low as 1.863 to as high as 7.0, reflecting the stochastic nature of the process on short timescales (Fig. 3B). Nonetheless, the cumulative estimate converges towards the true value of π, displaying clear law of large numbers behaviour with the absolute error scaling by (Fig. 3C and D). This directly mirrors the statistical efficiency of digital Monte Carlo sampling, confirming that passive Brownian sampling can, in principle, achieve the same asymptotic accuracy without consuming energy in re-sampling cycles. To quantify uncertainty in the running estimate of π, we computed 95% Bayesian credible intervals from the cumulative rod and crossing counts, treating crossings as a Poisson-distributed variable. The intervals narrow steadily over time, reflecting the increased statistical confidence as more data accumulates.

Fig. 3 Results of π estimation by Brownian rods. (A) π estimate across the full set of frames with cumulative estimate (blue line) being calculated from the running total numbers of rods and crossings. Individual frame estimates are shown as red crosses. The 95% credible interval is highlighted in blue. (B) Same dataset as (A) but showing the full range of individual frame π estimates. The standard deviation of the individual estimates is highlighted in orange. (C) Final 1000 frames of the dataset showing the estimation becoming stable within 3 significant figures of π. (D) log–log plot of the absolute error (ε) versus frames with a fitted linear regression line (coefficient ∼ −0.5).

To explore whether the convergence behaviour depended on the temporal extent of the experiment, we re-analyzed the same video starting from various frame indices (5000, 10 000, 15 000, 20 000 and 25 000) (SI). In all cases, the cumulative estimates converge within π ± 0.02 suggesting the long-term convergence is robust. For example, over the original dataset, the estimate follows an increasing trend which could be interpreted as a physical effect – such as rod length gradually decreasing over time due to dissolution. However, when the analysis is restricted to a subset of the data beginning at later frames, the trend is predominantly decreasing. This suggests that trends in the cumulative estimate may result from statistical fluctuations, rather than underlying changes in the physical system.

At the conclusion of the experiment, a total of 2 820 336 rods and 898 483 crossings had been detected, yielding an estimate,

and thus an absolute error, ε, of

ε = |[small pi, Greek, circumflex] − π| = 2.6 × 10⁻³.

Addressing the statistical uncertainty in our system, we obtain

[small pi, Greek, circumflex] = 3.1390 ± 0.0053(stat.)

where the statistical error is obtained from the 95% Bayesian credible intervals. A systematic error exists for the estimate which comes from setting d = 2 × 〈l〉. This is measured to be,

〈l〉 = 3.11 μm ± 0.20 μm,

resulting in a systematic error of,

δπ(syst.) = 0.21,

which renders the final π estimate as,

[small pi, Greek, circumflex] = 3.14 ± 0.21.

This reflects the limitations of a physical system, specifically, the polydispersity of the colloidal rods. Nonetheless, the much lower statistical uncertainty, which decreases with increasing sample size indicates that our system behaves as a Monte-Carlo solver. A more experimentally precise system (e.g., by using a higher resolution camera) would exhibit the same statistical behaviour.

At this total rod count, the number of crossings required to produce an estimate closest to π would be 897 741, resulting in

ε = 5.94 × 10⁻⁷.

This means at the conclusion of the measurement, there was a deviation of 742 crossings from the ideal behaviour. As would be expected, the magnitude of the deviation rate (in this case 742/898 483 ≈ 8 × 10⁻⁴) is inversely proportional to the magnitude of the loss of accuracy (in this case 10⁻³/10⁻⁷).

The closest approach to the true value of π occurred at frame 28 388, at which point 2 580 938 rods and 821 538 crossings had been recorded. The resulting estimate,

This was the closest possible value at this rod count and corresponds to an absolute error of

ε = 3.09 × 10⁻⁷.

Interestingly, within the dataset we identified 13 separate instances of four-frame sequences that recorded 355 rods and 113 crossings. These combinations yield

with an error of just

ε = 2.67 × 10⁻⁷.

If the experiment had consisted of only such brief, 8 second sequences, the estimates would have been extraordinarily accurate. This would, of course, be cheating the system by using approximations of π known since antiquity.⁴⁸

As a comparison, we implemented an OpenCV video generator that created binary videos mimicking the videos of our Brownian rods after pre-processing. The difference here was that, instead of exhibiting Brownian motion, between every frame the rods would be randomly distributed with a new angle and position using the Numpy PCG64 pseudo-RNG. After running the Buffon experiment on this new video, we compared the quality of π estimation finding that the Numpy based approach yielded much stronger convergence towards π (SI). The cumulative estimation for π remains stable within 3 significant figures after only 5000 frames. Nonetheless, this method requires energy consumption on the order of 4000 J min⁻¹ whereas the Brownian motion comes free.⁴⁹

Notably, bands form in the scatter plots of frame-specific π estimates, clustering around integers and common fractions (e.g., 3, 3.5, 4). This pattern is a direct consequence of measuring a discrete, physical system – especially when the sample size per frame, N, is relatively small and the mathematical operations are minimal. Because [small pi, Greek, circumflex] is computed as a ratio of two integers, it is constrained to rational values. In essence, the inverse of our estimates form a subset of Farey sequences F_N for N ∈ [58, 120].⁵⁰ As with Farey sequences, values formed from integers with many common divisors (e.g., 2, 3, 4, 5) appear with higher frequency. For example, [small pi, Greek, circumflex] = 3 can be formed by many rod-crossing pairs ((60, 20), (63, 21), …, (117, 39), and (120, 40)). In contrast, less common values like [small pi, Greek, circumflex] = 4.1 only occur in two specific combinations: (62, 20) or (93, 30). A more detailed description and visualisation can be found in SI.

This feature reflects a fundamental aspect of measuring physical systems: such systems rely on discrete observations, which inherently produce discretised outputs. As a result, certain values are simply inaccessible. While increasing the measurement resolution (e.g., by observing more events) or performing more mathematical operations (i.e., increasing permutations of event counts) can reduce this effect, a fundamental limitation remains – tied to the discrete nature of counting and observation. Notably, this limit is shared with electronic computational methods using floating points but which must still be embodied in a physical system, i.e., n-bit bytes are discretised to 2ⁿ values. Interestingly, this limitation parallels that of electronic computation, where floating-point numbers are constrained by finite bit-depth – n-bit registers can only represent 2ⁿ discrete states.

3 Conclusions

This microscale Buffon's needle system is a proof-of-principle that Brownian motion can perform unbiased Monte Carlo integration without digital computation or external energy input. Here, the physical properties of the rods – length, orientation, planar motion – directly map onto the variables of the Buffon problem. The colloidal architecture of the system is inherently parallel in which each unit randomises autonomously and independently.

Currently, the set-up relies heavily on digital systems to do both the measurement and mathematical operations. The next fundamental step would be to encode processing power into the physical system, for example through electrochemical signals with capacitive sensing, or microfluidics with on-chip counting.

While the present system is a minimal demonstration, it establishes a conceptual foundation for more complex physically implemented Monte Carlo algorithms. By re-framing passive Brownian motion as a computational engine, we link the physics of soft matter to the mathematics of sampling, opening new directions in alternative computing paradigms. Naturally, we do not envision such systems as being advanced π-estimators. But by careful selection of a colloidal system and a geometric probability problem that can be mapped to stochastic-computations; they may offer a zero-energy-input resource if an energy-friendly readout method is devised in the future.

4 Methods

4.1 Reagents

Zinc nitrate hexahydrate (ZnO(NO₃)₂·6H₂O), hexamethylenetetramine (HMTA), polyvinylpyrrolidone (M_W 40 000) (PVP40), sodium citrate (NaCit), 1-pentanol, ethanol, NH₃, and Tween-20 surfactant were all purchased from Sigma-Aldrich. All chemicals were analytical grade and used as purchased without any further purification.

4.2 SiO₂ microrod synthesis

SiO₂ microrods (1 μm long, 0.2 μm diameter) were synthesised according to procedure.⁵¹ Briefly, 2.5 g PVP40 was dispersed in 25 mL 1-pentanol by ultrasonication in a 50 mL Falcon tube. 2.5 mL ethanol, 0.7 mL H₂O and 0.25 mL 0.18 M NaCit were added. After shaking, 0.5 mL of a fresh NH₃ solution (25% in H₂O) was added, followed by shaking and the addition of 0.25 mL TEOS. The solution was allowed to rest for 14 h. The obtained rods were separated by centrifugation for 20 min at 5000 rpm and washed three times with ethanol. Larger SiO₂ glass micro rods with a diameter of 3 μm were purchased from Nippon Electric Glass Co., Ltd.

4.3 ZnO microrod synthesis

ZnO microrods (3 μm long, 0.5 μm diameter) were synthesised according to procedure.^52,53 Briefly, 0.594 g of ZnO(NO₃)·6H₂O and 0.056 g of HMTA were dissolved in 150 mL of DI water in a round-bottomed flask. The flask was closed with a rubber stopper and syringe then heated under stirring in a water bath at 90 °C for 30 minutes. After this the products were separated by centrifugation for 5 in at 5000 rpm before and washed 3 times with DI water. The rods were characterised using scanning electron microscopy (ZEISS Gemini SEM 300) and optical microscopy (Zeiss Axioscope 5).

4.4 Rod measurement

ZnO rods were dispersed in 0.15% Tween-20 surfactant solution and 20 μL were deposited into a well on a plasma-treated glass slide, and the well was covered with a glass slip to prevent evaporation. After sedimentation, the resulting number density on the surface was approximately 3000 rods per mm². The rods were visualised in microscopy using a Zeiss Axioscope 5 upright microscope equipped with a 50×/0.55 HD DIC lens. Images were recorded in 8-bit grayscale (2048 px × 1536 px = 212 μm × 160 μm) using an Imaging Source DFK 33UX265 Camera. To measure the angular autocorrelation function, the rods were recorded diffusing at 120 fps and the autocorrelation was measured with lags of 1/120 s. The autocorrelation of the rods was measured using,
where x is the variable (angle or position), τ is the lag time, σ² is the variance, n is the total number of time points, and μ is the mean of the time series, calculated as:

To perform the Buffon's needle experiment rods were recorded at intervals of 2 s.

4.5 Rod tracking

All image processing and rod tracking was done using an OpenCV-based python code.⁵⁴ In summary, the 8-bit grayscale videos were first binarised by thresholding to a fixed value such that only rods were visualised as contours. OpenCV was then used to measure the area of individual contours and perform principal component analysis on each. The median contour area was measured and if the area of a detected contour exceeded 1.5× the median, then it was assumed to be the measurement of two overlapping rods and was discounted from analysis. The length of the rods was calculated as the furthest distance between two points along the principal axis and the angle was derived from these two points. A line connecting the two points was generated to represent the 2D rod objects as 1D lines. A series of vertical grid lines was imposed on the images, separated by a fixed value equal to 2 times the median length of the rods in the sample. This was performed based on the initial frame of the measurement and remained fixed throughout. Intersections between rod-lines and grid-lines were detected by performing a counter-clockwise detection method on each rod with the grid-lines. For every frame, the total number of rods in the frame and the total number of intersections is counted and added to a cumulative total for each. The cumulative π estimate was calculated as the total number of rods divided by the total number of crossings.

Author contributions

C. Maslen: conceptualisation, data curation, formal analysis, investigation, software, writing – original Draft. L. Nicholson: investigation, visualisation, validation. J. Simmchen: conceptualisation, funding acquisition, project administration, supervision, writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data used in this manuscript is available in the supplementary information (SI) or from the corresponding author. The code for rod tracking and performing the Buffon experiment is available on Github at https://github.com/SimmchenGroup/Buffon_Needle.

Supplementary information including supporting videos, diffusion behaviour of differently sized rods, an example of image tresholding, pi-estimates with different starting frames, a comparison with simulated rods and information on band formations. See DOI: https://doi.org/10.1039/d5sm00844a.

Acknowledgements

We thank Martin Wittman and Maximilian Voigtmann for providing rods for our studies. This project is funded in part by the Advanced Research + Invention Agency (ARIA). JS acknowledges the DFG-ANR project Rodrolls - 490954343, ANR-21-CE30-0058.

Notes and references

1. I. E. Agency, Energy and AI: Energy Demand from AI, https://www.iea.org/reports/energy-and-ai/energy-demand-from-ai, 2025, Accessed: 2025-10-08.
2. E. García-Martín, C. F. Rodrigues, G. Riley and H. Grahn, J. Parallel Distrib. Comput., 2019, 134, 75–88.
3. A. Adamatzky, Advances in unconventional computing: Volume 1: Theory, Springer, 2016, vol. 22.
4. A. Steane, Rep. Prog. Phys., 1998, 61, 117.
5. Quantum Computing: Progress and Prospects, ed. M. Horowitz, E. Grumbling, Quantum Computing: Progress and Prospects, 2019 DOI:10.17226/25196, https://nap.nationalacademies.org/catalog/25196/.
6. P.-N. Nguyen, Digital Finance, 2025, 1–40.
7. M. AbuGhanem, J. Supercomput., 2025, 81, 687.
8. L. M. Adleman, Sci. Am., 1998, 279, 54–61.
9. A. Lopiccolo, B. Shirt-Ediss, E. Torelli, A. F. A. Olulana, M. Castronovo, H. Fellermann and N. Krasnogor, Nat. Commun., 2021, 12, 4861.
10. J. Gorecki, K. Gizynski, J. Guzowski, J. N. Gorecka, P. Garstecki, G. Gruenert and P. Dittrich, Philos. Trans. R. Soc., A, 2015, 373, 20140219.
11. Y. Viero, D. Guérin, A. Vladyka, F. Alibart, S. Lenfant, M. Calame and D. Vuillaume, Adv. Funct. Mater., 2018, 28, 1801506.
12. M. Crepaldi, C. Mohan, E. Garofalo, A. Adamatzky, K. Szaciłowski and A. Chiolerio, Adv. Mater., 2023, 35, 2211406.
13. R. Fortulan, N. Raeisi Kheirabadi, P. Mougkogiannis, A. Chiolerio and A. Adamatzky, ACS Omega, 2024, 9, 42127–42136.
14. X. Wang and F. Cichos, Nat. Commun., 2024, 15, 774.
15. V. A. Baulin, A. Giacometti, D. A. Fedosov, S. Ebbens, N. R. Varela-Rosales, N. Feliu, M. Chowdhury, M. Hu, R. Füchslin and M. Dijkstra, et al. , Soft Matter, 2025, 21, 4129–4145.
16. J. E. Gentle, Random number generation and Monte Carlo methods, Springer, 2003, pp. 1–40.
17. T. E. Hull and A. R. Dobell, SIAM Rev., 1962, 4, 230–254.
18. S.-I. Kim, H.-J. You, M.-S. Kim, U.-S. An, M.-S. Kim, D. H. Lee, S.-T. Ryu and Y.-K. Choi, Sci. Adv., 2024, 10, eadk6042.
19. M. Stipčević and Ç. K. Koç, Open problems in mathematics and computational science, Springer, 2014, pp. 275–315.
20. D. Vodenicarevic, N. Locatelli, A. Mizrahi, J. S. Friedman, A. F. Vincent, M. Romera, A. Fukushima, K. Yakushiji, H. Kubota and S. Yuasa, et al. , Phys. Rev. Appl., 2017, 8, 054045.
21. S. Misra, L. C. Bland, S. G. Cardwell, J. A. C. Incorvia, C. D. James, A. D. Kent, C. D. Schuman, J. D. Smith and J. B. Aimone, Adv. Mater., 2023, 35, 2204569.
22. D. P. Kroese, T. Brereton, T. Taimre and Z. I. Botev, Wiley Interdiscip. Rev.: Comput. Stat., 2014, 6, 386–392.
23. E. Dickinson and S. R. Euston, Adv. Colloid Interface Sci., 1992, 42, 89–148.
24. S. Mohamed, M. Rosca, M. Figurnov and A. Mnih, J. Mach. Learn. Res., 2020, 21, 1–62.
25. A. Mehonic, D. Ielmini, K. Roy, O. Mutlu, S. Kvatinsky, T. Serrano-Gotarredona, B. Linares-Barranco, S. Spiga, S. Savel’ev and A. G. Balanov, et al. , APL Mater., 2024, 12, 109201.
26. Y. Zhang, B. Zou, X. Jin, Y. Luo, M. Song, Y. Ye, Q. Hu, Q. Chen and A. C. Zambroni, Appl. Energy, 2025, 377, 124697.
27. F. Reif, Fundamentals of statistical and thermal physics, Waveland Press, 2009.
28. T. Hylton, MDPI Proc., 2020, 47, 23.
29. A. Agrawal, I. Chakraborty, D. Roy, U. Saxena, S. Sharmin, M. Koo, Y. Shim, G. Srinivasan, C. Liyanagedera and A. Sengupta, et al. , IEEE Trans. Very Large Scale Integr. Syst., 2020, 28, 2481–2494.
30. L. B. Kish, Appl. Phys. Lett., 2006, 89, 144104.
31. M. Aifer, S. Duffield, K. Donatella, D. Melanson, P. Klett, Z. Belateche, G. Crooks, A. J. Martinez and P. J. Coles, 2024 IEEE International Conference on Rebooting Computing (ICRC), 2024, pp. 1–20.
32. L. Bachelier, Annales scientifiques de l’École normale supérieure, 1900, pp. 21–86.
33. A. Einstein, Ann. Phys., 1906, 19, 371–381.
34. P. Hänggi and F. Marchesoni, Chaos, 2005, 15, 26101.
35. A. Scholz, P. Arya and J. Aghassi-Hagmann, AIP Adv., 2024, 14, 105114.
36. J. F. Ramaley, Am. Math. Monthly, 1969, 76, 916–918.
37. G. L. L. de Buffon, Euvres Philosophiques, 1777, 4, 46–123.
38. A. F. Bielajew, Monte Carlo techniques in radiation therapy, CRC Press, 2021, pp. 3–15.
39. K. Kamat, P. M. Naullage, V. Molinero and B. Peters, J. Chem. Phys., 2022, 157, 214113.
40. C. S. Petrie and J. A. Connelly, IEEE Trans. Circ. Syst. I: Fundam. Theory Appl., 2000, 47, 615–621.
41. N. C. Laurenciu and S. D. Cotofana, 2015 IEEE international symposium on circuits and systems (ISCAS), 2015, pp. 2724–2727.
42. S. Matsuoka, S. Ichikawa and N. Fujieda, Int. J. Circ. Theory Appl., 2021, 49, 3354–3367.
43. V. Reyes-Alejo and J. R. Guzman-Sepulveda, Opt. Lasers Eng., 2025, 194, 109176.
44. S. Broersma, J. Chem. Phys., 1960, 32, 1626–1631.
45. C. C. Miller, Proc. R. Soc. London, Ser. A, 1924, 106, 724–749.
46. M. M. Tirado and J. G. De La Torre, J. Chem. Phys., 1979, 71, 2581–2587.
47. M. M. Tirado, C. L. Martnez and J. G. De La Torre, J. Chem. Phys., 1984, 81, 2047–2052.
48. R. A. Nowlan, Masters of Mathematics: The Problems They Solved, Why These Are Important, and What You Should Know about Them, Springer, 2017, pp. 191–202.
49. B. Antunes and D. R. Hill, J. Data Sci. Intell. Syst., 2024, 53, 100655.
50. G. H. Hardy and E. M. Wright, An introduction to the theory of numbers, Oxford University Press, 1979.
51. A. Kuijk, A. Van Blaaderen and A. Imhof, J. Am. Chem. Soc., 2011, 133, 2346–2349.
52. Y. Bao, C. Wang and J.-Z. Ma, Ceram. Int., 2016, 42, 10289–10296.
53. L. Wang, M. Borrelli and J. Simmchen, ChemPhotoChem, 2021, 5, 933–939.
54. G. Bradski, Dr. Dobb's Journal of Software Tools, 2000.

---

Footnote

† The mathematical formalism for Brownian motion predates its physical explanation. In 1900, Louis Bachelier described similar stochastic processes in his attempts to model the Paris stock exchange – underscoring the deep connection between mathematical stochasticity and physical thermodynamic systems.³²

---