Lancelot Barthe

Martinus H V Werts
email: martinus.werts@ens-rennes.fr

Sedimentation of colloidal nanoparticles in fluids: efficient and robust numerical evaluation of analytic solutions of the Mason-Weaver equation

The Mason-Weaver equation (MWE) is an advection-diffusion equation that describes the timeevolution of the concentration profile of a solution of nanoparticles undergoing sedimentation and Brownian motion from an initially homogeneous state towards sedimentation equilibrium. In spite of the availability of analytic solutions, recent work has used numerical schemes to obtain practical solutions for the MWE. Here, the numerical evaluation of analytic solutions of the MWE is investigated using standard floating-point computations. It was found that the numerical evaluation is not always straightforward as this involves summing over an infinite series of terms which is sometimes automatically truncated due to limitations in floating-point computation. By combining several analytic expressions, each having its own range of validity in the MWE parameter space, robust and computationally efficient numerical evaluation of the solution is finally achieved. The expressions and the numerical procedure have been coded into a computer program enabling practical calculation of nanoparticle sedimentation profiles.

I. INTRODUCTION

The development of synthetic functional nanoparticles for applications in medicine, solar energy and so forth has brought renewed interest in the sedimentation behaviour of such small particles dispersed in fluids (i.e. colloidal solutions [START_REF] Shimizu | Thermodynamic Stability Condition Can Judge Whether a Nanoparticle Dispersion Can Be Considered a Solution in a Single Phase[END_REF]). The sedimentation behaviour of sufficiently small colloidal particles is governed both by sedimentation (which forces the particles in a specific direction, typically towards the bottom of the container) and diffusion (which tends to disperse the particles through the entire available volume). The basic model describing this behaviour was established almost a century ago in the form of a partial differential equation now known as the Mason-Weaver equation (MWE) [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF].

In their original publication, Mason and Weaver (M&W) gave an analytic solution for their equation. This solution takes the form of a periodic infinite series (i.e. a Fourier series) that converges slowly to the true solution, in particular at short times, requiring summation of many terms of the series to overcome the oscillatory behaviour known as the Gibbs phenomenon. More anecdotally, there are typographical errors in some of the formulae in the original paper, which does not facilitate the task of coding the expressions into a working computer program. It may be for these reasons that recent work has resorted to using finite-difference schemes to calculate numerical solutions for the MWE [START_REF] Alexander | Gravitational Sedimentation of Gold Nanoparticles[END_REF] [START_REF] Midelet | The Sedimentation of Colloidal Nanoparticles in Solution and Its Study Using Quantitative Digital Photography[END_REF]. These finitedifference schemes are computationally inexpensive for short times, and may also be adapted to cases that deviate from the classic MWE (e.g. different initial conditions, non-constant diffusion and sedimentation coefficients). However, finite-difference schemes do not generally conserve the total mass in the system, may display unphysical numerical divergence at very long times, and become expensive for calculation of the solution at long times.

Therefore, it is desirable to have an efficient method for evaluating the analytic solution that exists. The analytic solution holds for all times, will not display deviations from mass conservation, and its computation is likely to be inexpensive even for longer times. Very recently, Holt et al. used the analytic solution in their work on the rôle of nanoparticle sedimentation in cytotoxicity studies. [START_REF] Holt | The Blanket Effect: How Turning the World Upside Down Reveals the Nature of Graphene Oxide Cytocompatibility[END_REF] They reported that they needed to numerically sum a great number of terms (over 10 5 ) in the expression by M&W to obtain satisfactory results, but did not give further details. Here, the numerical evaluation of the expressions derived by Mason & Weaver is further investigated over a wide range of parameters.

In addition to the infinite-series solution (ISS) requiring computation of many terms, M&W derived an integral solution (IntS) that works for short times and does not involve summation over an infinite series. On basis of the original work, further expressions can be derived for the case where sedimentation begins to dominate over the Brownian diffusion term, i.e. at high Péclet numbers (hiPe). We found that, with current 'double precision' floating-point arithmetic and scientific software libraries, each of these three solutions has its own range of parameters where it produces reliable results. The standard numerical floating-point precision of 64 bits was used. In certain critical cases, 80 bits 'extended precision' was tested, but this did not substantially improve the numerical results. In contrast, switching to 32 bits 'single precision' floating point proved highly detrimental to numerical evaluation of the expressions.

A computer code was developed for numerical evaluation the different analytic solutions of the MWE in its dimensionless form. The program also provides a function that automatically selects evaluation of the relevant expression on basis of the input parameters, and a function that handles the interface between the 'real world' parameters and the dimensionless expressions. The latter can readily and transparently be used to calculate concentration profiles for Brownian nanoparticles undergoing sedimentation in solution.

II. THE MASON-WEAVER EQUATION

With concentration c, position z and time t, the timeevolution of the nanoparticle concentration profile c = c(z, t) of a colloidal solution undergoing sedimentation is illustrated in Figure 1, from an initially homogeneous distribution of nanoparticles to an equilibrium nanoparticle gradient. This evolution is governed by the MWE, Eqn. [START_REF] Shimizu | Thermodynamic Stability Condition Can Judge Whether a Nanoparticle Dispersion Can Be Considered a Solution in a Single Phase[END_REF].

∂c ∂t = D ∂ 2 c ∂z 2 + sg ∂c ∂z (1)
which has boundary conditions

D ∂c ∂z + sgc = 0 @(z = 0, z = z max ) (2)
and initial condition

c(z, t = 0) = c 0 (3)
D, s and g are, respectively, the diffusion coefficient of the particles in the fluid, their sedimentation coefficient and the gravitational acceleration. The z-axis corresponds to the height position in the cell, with z = 0 being the bottom and z = z max the top.

A dimensionless form of the Mason-Weaver equation is obtained through application of the changes of variable expressed in Eqn. [START_REF] Midelet | The Sedimentation of Colloidal Nanoparticles in Solution and Its Study Using Quantitative Digital Photography[END_REF].

z 0 = D sg α = z 0 z max t sed = z max sg τ = t t sed y = 1 - z z max (4)
The dimensionless form of the MWE for c = c(y, τ ) is then given by Eqns. ( 5)- [START_REF] Barthe | MasonWeaver-analytic[END_REF].

∂c ∂τ = α ∂ 2 c ∂y 2 - ∂c ∂y (5)
α ∂c ∂y -c = 0 @(y = 0, y = 1) (6) c(y, τ = 0) = c 0 (7)

III. INFINITE-SERIES SOLUTION (ISS)

M&W showed that Eqn. ( 8) is a solution to the dimensionless form of the MWE, Eqn. ( 5). [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF] c(y, τ )

c 0 ISS = exp y α α exp 1 α -1 + 16α 2 π exp 2y -τ 4α ∞ m=1 exp (-αm 2 π 2 τ )m (1 + 4π 2 m 2 α 2 ) 2 1 -(-1) m exp -1 2α sin (mπy) + 2πmα cos (mπy) (8)
The first term of this 'infinite series' solution (ISS) corresponds to the steady-state solution, c ss (y), while the second term (i.e. the infinite series) corresponds to the transient solution, c t (y, τ ). This is expressed as Eqn. [START_REF] Harris | Array Programming with NumPy[END_REF].

c(y, τ ) c 0 = c ss (y) + c t (y, τ ) (9)
Taking the limit of the complete solution when τ tends towards infinity yields the steady-state solution, since the transient solution vanishes at long times (Eqn. [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]).

lim τ →∞ c(y, τ ) c 0 = c ss (y) (10)
The transient term of the solution has the following form (Eqn. ( 11)).

c t (y, τ ) = s f (y, τ ) ∞ m=1 s m (y, τ ) ( 11
)
where s m (y, τ ) and s f (y, τ ) are, respectively, the terms of the series and a constant pre-factor with respect to m that can be taken out of the summation.

As pointed out earlier, [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF][5] [START_REF] Van Holde | Rapid attainment of Sedimentation Equilibrium[END_REF] the series converges rapidly for large values of {τ, α}, but slowly for small values. A large number of terms for the series (e.g., > 100000) might be evaluated and added to the final result in the hope of reaching convergence, as was done by Holt et al. [START_REF] Holt | The Blanket Effect: How Turning the World Upside Down Reveals the Nature of Graphene Oxide Cytocompatibility[END_REF]. However, it was found in the present work using standard double precision floating-point arithmetic that the factor exp (-αm 2 π 2 τ )m

(1+4π 2 m 2 α 2 ) 2
in the series summation of (8), numerically evaluates as strictly zero already for modest values of m. Any terms evaluated beyond that point do not contribute numerically to the solution, leading to a false convergence. Switching from the standard 64-bit to 80-bit floating-point precision did not substantially alleviate this problem.

This numerical artefact limits the number of terms that can be usefully evaluated for the ISS. It also provides a clear condition for stopping evaluation of the series: the loop evaluating terms for increasing m can be ended once exp (-αm 2 π 2 τ )m (1+4π 2 m 2 α 2 ) 2 numerically evaluates to zero. The result obtained in this manner is the best result numerically attainable, which does not necessarily mean that that result is sufficiently close to the correct solution.

Satisfactory results with the numerical evaluation of the ISS are obtained systematically only for α ⩾ 0.02. For α < 0.02, significant high-frequency numerical oscillations are often present in the 'best' numerically attainable result, in particular at low values for τ (short times). It is numerically impossible to add further terms. Low values of α correspond to high Péclet numbers, where sedimentation dominates diffusion giving rise to a sharp sedimentation front, requiring a large number high fre-quency terms to be included.

The limits for α and τ for correct evaluation of the ISS were investigated in more detail, and it was found possible to work at values for α slightly lower than 0.02 (say, 0.01), but only for values for larger values of τ (say, 0.5). It is possible to automatically detect numerical artefacts in evaluated ISS concentration profiles, by testing if the result is monotonous (within a certain tolerance) and by detecting deviations is the mass conservation. This may identify satisfactory numerical evaluations of the ISS at α < 0.02. However, the safe choice was made finally to constrain validity of the numerically evaluated ISS to α ⩾ 0.02. At α < 0.02, other analytic expressions exist that have better numerical behaviour, not needing any summing over a series (see the "high Péclet" solution below).

Under the constraint that α ⩾ 0.02 and using the aforementioned numerical vanishing of the factor

exp (-αm 2 π 2 τ )m (1+4π 2 m 2 α 2 ) 2
as a stopping condition, the number of terms of the series that is included in the numerical evaluation of the ISS is plotted as a function of τ for various values of α in Figure 2. This graph clearly illustrates the strong increase in number of terms for low α, τ . This makes evaluation of the ISS in those cases computationally expensive. The number of terms reaches almost 2000 for α = 0.02 and τ = 0.001, but drops rapidly as diffusion starts dominating sedimentation and when solutions close to the steady-state are evaluated.

IV. INTEGRAL SOLUTION (INTS) FOR SMALL

{α, τ } At small {α, τ }, where a large number of terms need to be included for evaluation of the ISS, M&W suggested [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF] an alternative solution that they obtained using the image method. This approximate analytic solution will be referred to as the 'integral' solution (IntS, Eqn. ( 12)). It is highly accurate at very small values of {α, τ }, as is clear from the very low RMS error between IntS and ISS.

c(y, τ ) c 0 IntS = 1 2 exp y -1 α 1 + y + τ -1 α θ 2 -y -τ √ 4ατ -θ 1 -y -τ √ 4ατ - τ απ exp y -1 α   exp - 2 -y -τ √ 4ατ 2 -exp - 1 -y -τ √ 4ατ 2   + 1 2 exp y α 1 + y + τ α θ 1 + y + τ √ 4ατ -θ y + τ √ 4ατ + τ απ exp y α   exp - 1 + y + τ √ 4ατ 2 -exp - τ + y √ 4ατ 2   + 1 2 θ 1 -y + τ √ 4ατ -θ τ -y √ 4ατ (12)
Using the same symbol as M&W, θ(x) = (2/ √ π)

x 0 exp (-t 2 ) dt is the error function, erf(x). IntS becomes less accurate at larger values, and fails badly when used outside of its domain of validity. For instance, the IntS exhibits convergence issues when α < 0.02 due to the limited stiffness of the error function at the junction (i.e anti-symmetry points) in the construction process by mirror image. Thus, both ISS and IntS cannot be used for α < 0.02.

In practice, the principal advantage of the IntS is the reduction of the computational effort compared to ISS. The range of values {α, τ } where the IntS is favourably applied instead of the ISS can be found by considering the mean square error (MSE), calculated between the IntS and the ISS, and identifying where this MSE is very small. Subsequently, a function τ switch (α) should be found in the aim of giving the time τ value below which the IntS will be used instead of the ISS, while remaining below a given MSE. The MSE of IntS relative to ISS is calculated over a grid spanning a range of α and τ . Visualized as a 'heat map' (Figure 3) it is clearly seen that there is a region at low {α, τ } where the MSE is very small, and where IntS can be used instead of ISS in order to limit the computation effort.

Using curve fitting, a parametric function τ switch (α) was obtained that defines a boundary for the region where IntS has low MSE (typically, less than 10 -8 ). The function is given by Eqn. (13), and is traced in red in Figure 3). Below the red trace, the IntS is used. Above the trace, the ISS is used. The range of α for which this τ switch is considered has been chosen to go from 0.02 to 20.

τ switch (α) = A 1 α k1 if α < 0.04625 A 2 α k2 else (13)
where A i , k i are the constants given in Eqn. (14).

A 1 = 8.134 k 1 = 0.9666

A 2 = 0.02934 k 2 = -0.8634 (14)
This parametrization was made on a specific computer system (Intel Core i7 (12th gen.) processor, Python 3.9.13 with numpy 1.22.3, scipy 1.9.1), and might depend somewhat on the numerical properties of the specific system, e.g. the software implementation for the evaluation of the error function θ(x). The optimal parameters for defining the switch between IntS and ISS are not expected to change significantly for other computer systems and languages (Matlab, FORTRAN), since these use the same basic numerical libraries and floating-point operations.

As illustrated in Figure 4, IntS is roughly between 50 to 10000 times faster than ISS, while having a MSE relative to ISS of less than 10 -8 . This is a useful gain in performance, especially in cases where multiple evaluations of profiles need to be made, such as in iterative curve fitting. The Figure also shows the switching between IntS and ISS at τ switch , after which only the ISS is considered (MSE becomes trivially 0). Taking the decision on which solution (IntS or ISS) to use generates negligible overhead.

V. APPROXIMATE "HIGH P ÉCLET" SOLUTIONS FOR VERY SMALL α

For very small α (α < 0.02, as explained above), neither the IntS nor the ISS can be reliably evaluated numerically. This is illustrated in Figure 5 which presents a typical case where both IntS and ISS display numerical artefacts. The case shown is for {α = 0.01, τ = 0.23}, but many combinations of parameters yield this kind of artefacts when α < 0.02, indicating the need for an alternative solution for small α. The situation where α is small corresponds to a high Péclet number, i.e. the contribution of nanoparticle diffusion becomes small compared to the sedimentation. The particles sediment with a clear sedimentation front that shows only little diffusive broadening. The advection due to sedimentation and the diffusive behaviour can be decoupled which leads to an approximate solution in the form of Eqn. (15), referred as "high Péclet" (hiPe) solution.

c(y, τ ) c 0 hiPe = 1 2 1 + θ y -τ √ 4ατ + K(τ ) exp y α α exp 1 α -1 (15)
In Eqn. (15), K(τ ) corresponds to a mixing between the steady state and the propagating diffusive front solutions. K(τ ) was obtained using mass conservation con-sideration (i.e.

K(τ ) = 1- 4ατ π   exp -(τ -1) 2 4ατ -exp -τ 4α   + (τ -1)θ 1 -τ √ 4ατ + τ θ τ √ 4ατ (16)
This "hiPe" solution was found to work well, numerically, for 0.0015 ⩽ α < 0.02.

For α < 0.0015, Eqn. ( 15) can be modified keeping only the sedimentation front and excluding the steadystate part by limiting the range of y to 0.99. This "hiPe2" solution, Eqn (17), does not conserve mass, and basically only describes the moving sedimentation front, which undergoes a minor broadening due to diffusion. For these very small values of α, the MWE is less relevant, since these systems are almost purely sedimentation systems, with a sharp sedimentation front. We include the "hiPe2" solution so that the entire range of α can be handled.

c(y, τ ) c 0 hiPe2 = 1 2 1 + θ y -τ √ 4ατ , 0 < y < 0.99 (17)

VI. COMBINED, PIECE-WISE ANALYTIC SOLUTION (PWS) FOR SUCCESSFUL NUMERICAL EVALUATION AT ANY α, τ

Considering the different regions in the {α, τ } plane where each of the analytic solutions to the MWE (ISS, IntS, hiPe, hiPe2) can be successfully evaluated numerically, a piece-wise analytic solution (PWS) can be constructed that covers the entire parameter space for α > 0 and τ > 0. Figure 6 contains a plot of the domains of application for the different analytic solutions discussed.

The PWS can be formulated according to Eqn. (18), with the understanding that α > 0 and τ > 0. Eqn. 18 forms the point of entry to the computer program module masonweaver analytic [START_REF] Barthe | MasonWeaver-analytic[END_REF] implementing the computational strategy for evaluation of the analytic solution of the MWE explained in this work. It was coded in the Python programming language [START_REF] Oliphant | Python for Scientific Computing[END_REF] with use of the numpy [START_REF] Harris | Array Programming with NumPy[END_REF] and scipy [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] scientific libraries, and may be readily ported to other languages. The dimensionless Eqn. 18 is implemented as the function routine MW adim(y, τ , α) which returns a vector1 containing the numerical values of the nanoparticle concentration profile for the (dimensionless) positions given by the vector y.

c(y, τ ) c 0 =                                    c hiPe2 (y, τ, α) if α < 0.
Finally, a 'real-world' (dimensional) function routine MW c profile(t, z max , D, sg, . . . ) is provided. It returns a pair of vectors (z, c) giving the numerical values for the (dimensional) positions and corresponding concentrations reproducing the full nanoparticle concentration profile using the analytic solutions of the MWE.

VII. COMPARISON WITH CURVES PUBLISHED BY MASON & WEAVER AND EXPERIMENTAL DATA

In addition to extensive testing of the PWS against the numerical solutions generated by the previously published finite-difference program, [START_REF] Midelet | The Sedimentation of Colloidal Nanoparticles in Solution and Its Study Using Quantitative Digital Photography[END_REF] the computational accuracy of the combined piece-wise solution was demonstrated by comparing its results to curves from the original publication by M&W. [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF]. To this end, these curves were digitized using the Engauge Digitizer software package, [START_REF] Mitchell | engauge-digitizer[END_REF] converting them into sets of data points representing c(y, τ )/c 0 as a function of y for different values of τ . These curves are shown in Figure 7, together with curves obtained using the numerical evaluation of the PWS by our computer program.

It is evident from Figure 7 that the curves evaluated by our PWS procedure are in accordance with the curves from the original publication. No precise error estimation is proposed here, since the digitization process by itself is a source of significant uncertainty. Interestingly, the PWS result for τ = 1, does not match the original M&W curves. This is a minor mistake in the original publication.

The numerically evaluated analytic solution also fits previously published experimental sedimentation profiles obtained using quantitative digital photography of gold nanoparticles in aqueous solution undergoing sedimentation. [START_REF] Midelet | The Sedimentation of Colloidal Nanoparticles in Solution and Its Study Using Quantitative Digital Photography[END_REF] This is shown in Figure 8 which revisits Figure 2 from the original publication, using the analytic solution instead of the numerical finite-difference solution. There is no obvious difference between the two sets of theoretical curves (analytic vs finite-difference). The analytic function is calculated more quickly, especially for longer times in the sedimentation process. It instantly generates the corresponding concentration profile, whereas the finite-difference solver needs to start at the beginning and iterate over time steps. The analytic solution perfectly preserves the total mass in the system, whereas the finite-difference scheme causes a small increase in the total mass. [START_REF] Midelet | The Sedimentation of Colloidal Nanoparticles in Solution and Its Study Using Quantitative Digital Photography[END_REF] VIII. CONCLUSION It is not an entirely trivial task to evaluate numerically the analytic solution to the Mason-Weaver equation. Straightforward evaluation of the ISS by simply including a great number of terms does not always lead to satisfactory solutions, due to higher terms being evaluated as strictly zero by the floating-point arithmetic. In the original paper by M&W, [START_REF] Mason | The Settling of Small Particles in a Fluid[END_REF] alternative approxi- mate solutions were suggested that can be applied to circumvent this problem (hiPe) or significantly alleviate the computational burden (IntS). Implemented in a Python computer program, [START_REF] Barthe | MasonWeaver-analytic[END_REF] the combination of these different solutions, each being operational in a specific range of parameters, leads to a computationally robust and efficient method for calculating nanoparticle sedimentation profiles.

Figure 1 .

 1 Figure 1. Sedimentation behaviour of Brownian particles as described by the Mason-Weaver equation (MWE). The particles are initially homogeneously distributed. Top: concentration profiles as a function of height at (a) t = 0, (b) intermediate t, and (c) sedimentation equilibrium (t → ∞, steadystate solution of the MWE). Bottom: visualization of the corresponding vertical gradients such as can be observed, e.g. with gold nanoparticles in solution.[4]. The z coordinates refer to the 'real-world' dimensional version of the MWE, where z = 0 is the bottom of the cell. The y coordinates in the dimensionless form move in the opposite direction, y = 0 being the top of the cell.

Figure 2 .

 2 Figure 2. The number of terms, Nterms, used in evaluation of the ISS as a function of τ for different values of α (⩾ 0.02). Inclusion of further terms is not possible due to very small factors numerically evaluating as zero in floating-point arithmetic (see Text).

Figure 3 .

 3 Figure 3. Logarithm of the MSE between the IntS and ISS (color map) as a function of α and τ (plotted as their logarithms). The red trace is the boundary defined by τ switch in Eqn. (13).

Figure 4 .

 4 Figure 4. Comparison of pure ISS and mixed IntS/ISS for the calculation of MW concentration profiles at α = 0.02 as a function of τ . Top: Calculation time for a single profile. Bottom: Mean-square error between IntS/ISS and ISS solutions, dropping to zero when ISS is used in both calculations. Profiles were calculated over y from 0 to 1 in 200 points. The calculation time was measured for the Python program executed on a single core of an Intel Core i7 (12th gen.) 64-bit processor running at 2.4 GHz. The green dotted vertical line indicates where the switch between IntS and ISS is made.

Figure 5 .

 5 Figure 5. Example of the failure of both ISS and IntS, when α < 0.02. Three analytic solutions to the Mason-Weaver were evaluated for {α = 0.01, τ = 0.23}. Only the hiPe solution gives a completely artefact-free concentration profile.

  τ ) c0 dy = 1) as proposed in Eqn. (16).

Figure 6 .

 6 Figure 6. Graph representing the range of application of each numerically evaluated analytic solution to the MWE as a function of {τ, α} in the combined PWS. Evaluation using standard floating-point arithmetic with 64 bits precision.

Figure 7 .

 7 Figure 7. Comparison of the curves from first Figure in the paper by M&W [2] (dotted lines) and the numerically evaluated analytic solution from the present work (continuous lines) for α = 0.025.

Figure 8 .

 8 Figure 8. Evolution of the vertical particle density gradient in an aqueous solution of 40 nm diameter gold nanospheres toward sedimentation equilibrium, observed using quantitative digital photography at 277 K. Photos were taken at t = 0, 1 h, 21 h, 27 h, 44 h, 52 h, 7 d, 9 d, 11 d, 14 d, 21 d, 23 d, 28 d, 35 d, and 39 d. Top: Experimental optical density profiles obtained from digital photographs.[4] The arrow indicates the direction of time. Bottom: theoretical optical density (concentration) profiles evaluated using the analytic solution of the Mason-Weaver equation, with D = 5.1 × 10 -12 m 2 s -1 and s = 7.9 × 10 -10 s.

In this text, a 'vector' is a one-dimensional array of floating-point numbers.

IX. ACKNOWLEDGEMENTS

We thank Prof. Steven Abbott for his enthusiastic and eloquent communications about colloid science. Our correspondence with Prof. Abbott about solving the Mason-Weaver equation was a source of motivation for the present work.

We thank Dr Rozenn Texier-Picard (département mathématiques, ENS Rennes) for insightful and stimulating discussions.


Raja Chatila

Rachid Alami

Thierry Siméon

Julien Pettré

Léonard Jaillet

Safe, Reliable and Friendly Intera ction Between Humans And Humanoids

Humanoid robots include a mobile base (walking or rolling) and manipulation devices, usually two arms. By the number of actuators and accessible workspace, they are certainly the most complex robotics systems designed. Unlike classical robots in industry, they have to interact directly with humans, and unlike service robots which have usually only a mobility capacity or limited manipulation, this interaction is a very close one. It is therefore mandatory that the safety and human-friendly aspects of this interaction be addressed as such, considering the problem in its entirety. This paper discusses these issues considering several constraints that must be respected by the robot motion planning system, related to the human-robot interaction distance. It introduces the concepts of safety distance, social distance, manipulation zone, and friendly posture.

Introduction

Research on humanoid robots has achieved several results on the issues related to their design and control for walking, navigation, manipulation, etc. On the other hand several issues related to human-friendly robots have also been addressed including the study of human behavior and modeling of human environment. Finally, the dependability issues of robots in terms of the reliability of achieving their tasks and the safety aspects related to their interaction with humans are also subject of investigation. This paper addresses new issues related all these aspects and specifically to the close interactions between humans and humanoid robots from the standpoint of the motion decisions that must be taken by the robot in order to ensure a :

• A safe interaction, i.e., that cannot harm the human, • A reliable interaction, i.e, that achieves the task adequately considering the motion capacities of the robot, and • A user friendly interaction, i.e, that takes into account a motion model and a communication model of the human. We discuss next these notions that are based on the identification of several models for interaction.

Safe Interaction

In human environments, safety issues are of concern. The basic functions for a safe sensory human/robot interaction include for instance compliant motions. But a central point is that robot must be able to take into account the safety of the human body. In [START_REF] Kuffner | Self-Collision Detection and Prevention for Humanoid Robots[END_REF], the authors address a similar question concerning the robot own body. In this paper we address the interaction with humans. In any of its motions the robot should not be able to harm the human, and therefore it must identify a safety distance to the human. However, the robot must be able to approach the human because of the necessity of their interaction (for example to handle some object to the human). Hence, this distance is not uniform and fixed, but depends on the interaction. Figure 1. depicts a safety distance which is represented by bubbles surrounding the human.

Human safety distance

The human has an articulated multi-body structure. We can identify a "safety suit" which is attached to the human body, composed of bubbles at a predefined distance from its skin and which represents a forbidden region in the robot configuration space. By permanently computing the minimal distance between its own body and the human safety suit, the robot will produce motions that are safe guarded and will never collide with the suit, and therefore keep the human body out of robot reach (figure 2). Even if mechanical safety devices such as touch sensors are built in the robot for safety, it is desirable that the robot plans motion that are safe. The motion planning system of the humanoid robot, considering its whole body, takes into account the geometrical model of the human body and includes a collision checker. The model may include a large number of geometric primitives, and it is well known that the most time consuming operation in motion planning are interference detection and collision-checking. The motion planning system Move3D that we developed [START_REF] Siméon | Move3D: a generic platform for path planning[END_REF][START_REF] Siméon | Cortès Computer aided motion: Move 3D within MOLOG[END_REF] is based on probabilistic motion planning techniques [START_REF] Kavraki | Probabilistic Roadmaps for Path Planning in High-Dimensional Configuration Spaces[END_REF] and includes an Interferencedetector [START_REF] Van Geem | KCD: a collision detector for path planning in factory models[END_REF] based on algorithms combining techniques proposed in [START_REF] Lin | Collision detection: Algorithms and applications[END_REF] while allowing to process non convex polyhedra together with other volumic primitives (eg. spheres, tubes, torus) as in [START_REF] Van Den Bergen | A fast and robust GJK implementation for collision detection of convex objects[END_REF]. Computation is performed in two stages: selection of possibly colliding pairs based on simple bounding volumes that approximate the geometry of the objects, and a precise interference detection limited to the pairs selected at the first stage. A hierarchical structure based on OBB-trees is constructed on top of the convex primitives (or facets), instead of the triangular decomposition required by [START_REF] Gootschalk | OBB-Tree: A hierarchical structure for rapid interference detection[END_REF]. This allows to reduce the size of the data structures in large models where most of the objects are modeled by simple primitives. The Collision-checker developed for determining whether a given path is collisionfree or not, is currently performed by multiple calls to the interference detection algorithms, using a dichotomic sampling of the path and a non uniform step computed from the distance to the obstacles.

Human interaction distances

In usual interactions between humans, some non written rules are respected that determine the distance between two persons (see the controversial proxemic theory of E. T. Hall [START_REF] Hall | The Hidden Dimension[END_REF]). Hall has proposed that four interpersonal distances were respected in society : Intimate distance (0 to 1.5 feet) for very close relationships, Personal distance (1.5 to 4 feet) for most usual friendly interactions and conversations, Social distance (4 to 12 feet) for more formal social interactions, and finally Public distance (over 12 feet) with minimal interactions. Of course these figures are approximate and related to culture and several other factors. However in human-robot interaction distances will also be defined because humans are used to respect such distances and would be surprised or would not feel at ease if a humanoid robot does not respect them. In addition, if the interaction with the humanoid robot is similar to the interaction with other humans, i.e., through speech, vision and touch, then such distance will be very natural to the human. Our purpose is to make them also natural for the robot. Hence, when a robot must listen to a human or speak to him or her, or when it moves along the human or when a human approaches the robots (figure 3) etc., it must be always keep the adequate distance. Again, distance checking in the robot motion planning and control system plays a central role.

Reliable interaction

In some situations, the robot and the human must exchange objects that are manipulated and handed from one to the other. In this case, the safety distance must be adapted to this kind of interaction, as well as the proximity of the robot and human bodies. But in order to actually be human friendly and reliable, the robot must understand what is the adequate position or motion to adopt so that the object manipulation task is feasible by the human in a manner that is the most comfortable for them. For this a human motion model is required, as well as a human-robot manipulation model for object exchange. Here again, the motion planning system plays a central role in checking the adequate constraints. Figure 4. shows how a specific zone can be defined for object exchange, which enables to violate the constraint of the safety distance. However, a closer safety distance is then defined to forbid direct contact (figure 5). For object exchange, the safety region can be defined closer to the human body. We can note that is such exchanges, the robot motions are more local and slower.

Friendly Interaction

In order to feel safe, the human must be able to see the robot and to expect its motions. This in turn means that the robot must have a model of human communication behavior . For example, when the robot speaks to the human it is more adequate to face him or her at, or when it must hand something to him/her, it is again better that the robot comes to face the human, even if this requires it to avoid obstacles or to plan a complex motion [START_REF] Alami | On the influence of sensor capacities and environment dynamics on collision-free motion plans[END_REF]. The robot should be endowed with a human communication model defining its capacities and adapted to the human which the robot interacts with (e.g, handicapped or not). Figure 6 shows an example of robot reaction when understanding the human "calling" gesture requiring a service from the robot. If the human requires some action from the robot, the robot motion must be adapted to the human behavior. Figure 7 shows the robot pouring water in the glass handed out by the human. Here the motion planning system must reason on the human reachable workspace and on the obstacles in the environment to compute a path of the whole robot body so that the pouring operation is possible and seems natural for the human.

Conclusion

In this paper we have introduced the requirements on the humanoid robot motion system for producing safe, reliable and friendly motions during its interaction with humans. These requirements translate into an efficient motion planning system with permanent collision checking of the robot whole body with the human body using the notion of safe distance and models of human interactions related to interaction distances, object manipulation, human motion and human communication. Work is in progress for experimenting these concepts in the motion planning system MOVE 3D.

Figure 1 :

 1 Figure 1: safety distances. Virtual "bubbles" surround the human.

Figure 2 :

 2 Figure 2: The human safety distance is a forbidden region frontier in the configuration space of the humanoid robot.

Figure 3 :

 3 Figure 3: The forbidden distance represents the human intimate distance. The robot must move away when in contact with it.

Figure 4 :

 4 Figure 4: An exchange zone related to the human manipulation model is defined for direct interaction.

Figure 5 :

 5 Figure 5: The forbidden region during object exchange is closer to human body.

Figure 6 :

 6 Figure 6: Even if the robot has gesture recognition capacities the glass is shown for service), it must produce also understandable motions in cluttered environments.

Figure 7 :

 7 Figure 7: The motion planner must plan a position of the robot and of the pouring gesture that look natural and are safe for the human.


Matt Eames

Zachary Larrabee

Arik Hananel

Frederic Padilla

Jean-Francois Aubry
email: jean-francois.aubry@espci.fr

Low-Cost Thermochromic Quality Assurance Phantom for Therapeutic Ultrasound Devices: a proof of concept

Keywords: Therapeutics, System & device design, General Physical Acoustics, Medical Transducers

High intensity focused ultrasound (HIFU) transducer acoustic output can vary over time due to an inconsistent power supply, damage to the transducer, or deterioration over time. Therefore, easy implementation of a daily quality assurance (DQA) method is of great importance for preclinical research and clinical applications.

We present here a thermochromic material-based phantom, validated by thermal simulations, and shown to provide repeatable visual power output assessments in fewer than 15 seconds that are accurate to within 10%.

Whereas current available methods such as radiation force balance measurements provide an estimate of the total acoustic power, we demonstrate here that the thermochromic phantom is sensitive to the shape of the acoustic field at focus, by changing the aperture of a multielement transducer with a fixed acoustic power. The proposed phantom allows the end-user to visually assess the transducer's functionality without resorting to expensive, timeconsuming hydrophone measurements or image analysis.

INTRODUCTION

With the increasing use of focused ultrasound (FUS) for research and development applications and for the treatment of medical conditions in the clinical arena, there is a growing need for a simple, expeditious, and inexpensive method for assessing constant and safe device output. Literature reports that, within a sample of 64 physiotherapy therapeutic ultrasound systems, the measured output power of 59% of the systems falls outside a +/-20% deviation from prescribed power output [START_REF] Schabrun | How Accurate are Therapeutic Ultrasound Machines?[END_REF]. Patient safety may be compromised by an unexpected higher dosage. On the contrary, the efficacy of the treatment would be affected by a lower dosage.

Existing methods for conducting quality assurance on ultrasound systems require complex set-ups and specific expertise. The most comprehensive method is the measurement of the 3D pressure field generated by the transducer with a hydrophone. Due to the intensities involved in therapeutic ultrasound, an appropriate hydrophone should be used, such as fiber-optic displacement sensors, commercial HIFU needle hydrophone, or membrane hydrophone with a protective coating against cavitation effects [START_REF] Haller | A comparative evaluation of three hydrophones and a numerical model in high intensity focused ultrasound fields[END_REF]. Three-dimensional measurements are time consuming. The 3D pressure field can be inferred from pre-focal 2D planar field acquisition using back-and forward-projection [START_REF] Clement | Field characterization of therapeutic ultrasound phased arrays through forward and backward planar projection[END_REF][START_REF] Sapozhnikov | Acoustic holography as a metrological tool for characterizing medical ultrasound sources and fields[END_REF][START_REF] Schafer | Transducer characterization using the angular spectrum method[END_REF]. The 2D projection of the HIFU beam can be quickly visualized by Schlieren imaging [START_REF] Kudo | A simple Schlieren system for visualizing a sound field of pulsed ultrasound[END_REF]). This technique takes advantage of acoustic-induced changes in optical refractive index. Nevertheless, quantitative characterization of the 3D pressure field requires several projections at varying angles coupled with tomographic backprojection, and is therefore not straightforward to extract [START_REF] Civale | Quality assurance for clinical high intensity focused ultrasound fields[END_REF]. Infrared cameras have been shown efficient for mapping the 3D heating distribution induced by the transducer in an absorbing material [START_REF] Khokhlova | Infrared mapping of ultrasound fields generated by medical transducers: Feasibility of determining absolute intensity levels[END_REF][START_REF] Shaw | The feasibility of an infrared system for real-time visualization and mapping of ultrasound fields[END_REF]. It provides excellent spatial and thermal resolutions but the main disadvantages are the cost of the camera and the inability of infrared radiation to travel through water.

Radiation force balance measurements can be used to quantify the total acoustic power emitted by the transducer (Payne et al. 2021). In these measurements, the transducer is aimed at a target (an ultrasound absorber) suspended in water tank from precision balance in a draft-free environment [START_REF] Aubry | High-intensity therapeutic ultrasound: metrological requirements versus clinical usage[END_REF]. The momentum carried by the ultrasonic wave is measured as a weight shift on the balance after reflection or absorption on the target [START_REF] Gélat | Relationship between acoustic power and acoustic radiation force on absorbing and reflecting targets for spherically focusing radiators[END_REF][START_REF] Maruvada | Acoustic power calibration of high-intensity focused ultrasound transducers using a radiation force technique[END_REF].

Non-invasive quantitative thermometry can also be achieved with ultrasound, based on detection and estimation of minute temperature-induced echo shifts, which can be processed to generate the 2D maps of temperature change. Important challenges remain however, requiring development of dedicated software architecture and signal processing, to compensate artifacts arising from tissues or samples deformation leading to echo decorrelation and ultrasound beam distortion [START_REF] Ebbini | Real-Time Ultrasound Thermography and Thermometry[END_REF]. For Magnetic Resonance guided Focused Ultrasound Treatment (MRgFUS), calibration can be achieved with the sole use of a tissue-mimicking hydrogel in conjunction with the MR. The thermal rise induced by the ultrasound beam inside the hydrogel can indeed be quantified and mapped with MR thermometry [START_REF] Aubry | High-intensity therapeutic ultrasound: metrological requirements versus clinical usage[END_REF][START_REF] Mcdannold | Quality assurance and system stability of a clinical MRI-guided focused ultrasound system: Four-year experience[END_REF]. Permanent thermochromic phantoms have been proposed to counteract the need for MR to image the thermal rise: the color of the 3D phantom changes permanently when the temperature reaches a threshold. The cross section can be assessed by cutting the phantom open [START_REF] Negussie | Thermochromic tissuemimicking phantom for optimisation of thermal tumour ablation[END_REF]. Thermochromic paper has been proposed as the most rapid quality assurance technique [START_REF] Butterworth | Exploiting Thermochromic Materials for the Rapid Quality Assurance of Physiotherapy Ultrasound Treatment Heads[END_REF]. Nevertheless, it relies on the analysis of images captured by a digital camera and on automatic fitting of the beam profile. Good and consistent lighting conditions are critical. Thermochromic paper approaches thus concentrated on improving the precision of the estimation of the spatial dimension of the color change. We propose here to use a thermochromic liquid crystal affixed to a highly absorbing material to record the time for inducing a color change, after the heat radiates through the absorbing material. We hypothesize that the time needed for the heat to diffuse from the transducer side to the thermochromic side is a function of power and that this measurement can be used to assess the stability of the transducer. We additionally investigate the sensitivity of the thermochromic phantom to the shape of the acoustic field at focus, by changing the aperture of a multielement transducer with a fixed acoustic power.

MATERIALS AND METHODS

Phantom

The phantom is a 2.8 mm thick highly attenuating ultrasound absorber (AptFlex F28, Precision Acoustics Ltd, Dorchester, UK) (Fig. 1 and Table 1), with a surface layer of thermochromic liquid crystals (TLCs) (R35C5B, LCR Hallcrest LLC, Glenview, IL, USA). The TLC material is black in color until it transitions to red at 35 °C followed by green at 36 °C and blue at 40 °C, allowing an observer to visually detect when temperature reaches and exceed these thresholds (Fig 1). The time to see color change, denoted as changing time, is nonzero because the absorbed heat must diffuse from the non-TLC side of the absorber to the TLC surface. Thus, the changing time for TLC of a given thickness is a function of baseline temperature, acoustic power and center frequency.

The phantom can be designed to be positioned either at the surface or at the acoustic focus of the transducer of interest. Fig. 2 depicts the two transducers and phantom geometries tested: one with a flat 11 mm aperture transducer (Model XDR058, Sonic Concepts, Bothel, WA) to which the phantom was designed to be coupled with acoustic gel directly to the face of the transducer; and another the Insightec Neuro device, a 30 cm-diameter hemispherical transducer, to which the phantom was mounted at the acoustic focus using a 3D printed holder. In both cases, the phantom measures 2.8 mm in thickness.

Numerical simulations

A time-domain finite difference simulation of ultrasound-induced thermal rise in the acoustic absorbing material was performed.

The pressure field was first measured experimentally. The pressure was then used as an input to the bioheat equation.

Pressure Field measurements to characterize the acoustical output of the transducers

The pressure field of a high-power 1.1 MHz transducer (Model XDR058, Sonic Concepts, Bothell, WA, USA) was measured in water with a 400 µm needle hydrophone (HNA-0400, Onda Corp., Sunnyvale, CA, USA) mounted on a three-axis motor gantry (AIMS III, Onda Corp., Sunnyvale, CA, USA). The active surface of the transducer is an 11 mm diameter disc. An x-z plane scan was performed, z being the axisymmetric axis of the transducer. The transducer was mounted on a mechanical holder to position it along the z-gantry axis. The scan was 30 mm along the x axis and 40 mm along the z axis with a 0.25 mm step for each axis. Assuming symmetry in the x and y directions, this scan was then radially projected to obtain a three dimensional pressure field, P(x,y,z).

The 3D pressure field pAptflex within the 2.8 mm-thick absorbing layer was deduced from the measured field in water pWater by adding an exponential attenuation as a function of distance:

𝑝 𝐴𝑝𝑡𝑓𝑙𝑒𝑥 (x, y, z) = 𝑝 𝑊𝑎𝑡𝑒𝑟 (x, y, z) • e -αz ( 1
)
where α is the attenuation coefficient in Np/m at the (x,y,z) coordinate. The pressure field in the Aptflex was then used as an input for the thermal simulations (see next section).

Table 1 provides material properties used in the simulation.

The pressure field of the Neuro system (Insightec, Tirat Carmel, Israel) was also measured with the same setup.

The transducer was operated at 670 kHz, with the central elements activated only (256 active elements out of the 1024-elements). The z axis corresponds to the axisymmetric axis of the transducer. Two planes were acquired:

the transverse XY plane and the longitudinal XZ plane. Each plane was 10x10 mm 2 with a 0.25 mm step in each direction.

Estimation of the thermal rise and time to color change by solving the bioheat equation

The physical heat transfer is modeled by a Matlab (Mathworks, Natick, MA, USA) simulation of the bio-heat equation throughout the experiment volume. The equation is solved with a finite-difference scheme in which the derivatives of the equation are replaced by second order polynomial approximation:

ρc p ∂𝑇 ∂𝑡 = k∇ 2 𝑇 + 𝛼𝑃 2 ρc (2)
Where T is temperature, 𝑘 is thermal conductivity, 𝛼 is the absorption coefficient, 𝑐 is speed of sound, 𝜌 is density, 𝑐𝑝 is specific heat, and P is the acoustic pressure amplitude field. The pressure profile from one scan is assumed to be constant for all simulations. This pressure is normalized and then modulated by the square root of the electrical power to model thermal change for escalating power values. The simulations are performed in a 61x61x44 grid with an isotropic resolution (0.3215 mm) and a 0.05 s time step. All voxels corresponding to water are set to a fixed temperature (17 °C) at each time step to account for its cooling effect. Dirichlet conditions are applied at each boundary planes except for the air interface where a Neumann condition is applied to account for the air thermal insulation, as compared to water.

The minimum time required to achieve a peak temperature of 35 °C in a single voxel at the distal plane of the absorbing medium is recorded, and the simulation is re-run across several acoustic output intensities.

Time to color change measurements

Measurements of the time until first-observed TLC material color change (T=35 °C) were recorded using a stopwatch on the two different setups.

Single element transducer

Measurements were performed across a [0.4-12] W range of electric power. All power levels were measured 8 times within six months.

Two DQA phantoms were assembled from different sections of single sheets of AptFlex and TLC material.

Coupling gel (Aquasonic 100, Parker Laboratories) was placed between the transducer and the phantom. Technology, Cambridgeshire, UK) connected to a PC. An example of the longest cool-down period (for the highest power) can be seen in Fig. 3.

Multielement transducer

Measurements were performed using the central 256 elements of the 1024-element hemi-spherical, 30 cm diameter, 15 cm focal length Exablate Neuro high frequency transducer (Insightec Inc., Tirat Carmel, Israel), across a [0.2-1.5] W range of electric power. Over the course of 6 months, all power levels have been measured

6 times with 4 repetitions. The system was operated at 670 kHz The dedicated Insightec user interface and hardware were used to power the transducer in continuous-wave mode for up to 45 s in duration, as necessary, to measure time to first TLC material color change. The Insightec logfile measures and records the actual total electrical energy delivered by the system during the sonication. This value was divided by sonication time to estimate the delivered average power.

Two DQA phantoms were assembled to suit the multielement transducer, and the same 3D-printed holder was used to ensure reliable repositioning of the phantom at the geometric focus across all measurements. A three minute cool down period was used between measurements, and the phantom, which was large enough to permit it to be placed 2 cm off-center in the plane of the focal spot, was rotated 90  between sonications such that more than 12 minutes elapsed between sonications at the same region of the phantom. Four measurements (one at each rotation angle) were recorded per power setting, per experiment session.

RESULTS

Pressure field measurements

The peak negative pressure is displayed, in decibels, in Fig. 4 in the transverse and longitudinal planes for both the single-element (a) and Insightec Neuro (b) transducers.

Thermochromic phantom measurements

The time to first color change is plotted for the single-element transducers (Fig. 5). During the course of testing, one transducer used for these measurements was damaged, as is evident in the collected data.

The time to first color change for the Insightec Neuro transducer is plotted in (Fig. 6). In order to reproduce the effect of reduced acoustic output, data were collected with the central 256 elements (labeled "Apodized") and with the full array, which had the effect of distributing the requested power across more elements by energizing the outer ring of 768 elements in addition to the central 256 (labeled "Unapodized"). This served to produce a shorter focal spot in the propagation direction. As with the single element case, the effect of this change is evident in the data. Of note, one set of measurements (blue, Fig 6) were supposed to be performed without apodization.

As this set of measurement seems to align with the apodized data, we double checked the logs of the Exablate

Neuro system and noticed that the time of acquisition corresponded to apodized sonications. The have been labeled

as "Apodized unintent.".

Single-element data and Insightec Neuro data are plotted against the associated bioheat simulation results in the respective Figures.

DISCUSSION

The single-element simulation reproduces the general trend of the experimental results (Fig. 5) but at high electric powers the simulation indicates a time-to-color-change that is 2 seconds faster than experimental results.

This could be due to both latency in the TLC material and observer reaction time.

Two different behaviors, both for simulations and experiments, were observed for each of the two apodization scenario tested with the Insightec Neuro (Fig. 6): the time to see a color change was much shorter without apodization (using only 256 central elements out of the 1024) that with apodization (using the full 1024 elements) due to the longer focal spot in the propagation direction without apodization. Nevertheless, at lower electric powers, the simulation indicates a time-to-color-change that is faster than experimental measurements. This is likely due primarily to observer reaction time at longer times to first color change, in which the color change is more gradual than those that occur at higher powers.

For optimal sensitivity, it is recommended to select an electrical power of 1.5 -2.0 W in the single-element case and an electrical power of 0.3 -0.5 W in the case of the Insightec Neuro. At these powers, DQA assessments are complete within 20 seconds, the color change is relatively sudden, and the typical user (with reaction time of 0.5 s) can detect a 10 % deviation in acoustic power output.

It must be mentioned here that at the starting point of the study, the idea was to measure the diameter of the color change (as can be seen on Fig 1) as a function of power. A color change can indeed be seen when TLC material is placed in front of a transducer, even without any absorbing material. Nevertheless, the first experiments showed that the diameter of the color change depends primarily on the size of the pressure field generated by the transducer at the surface of the TLC. For this reason, we added an absorbing material between the TLC and the transducer in order to take advantage of the time needed for the temperature to build up in the absorbing material.

The time until first visible TLC color change was found to depend more significantly on the acoustic power output of the transducer, as presented in this paper.

It is noted that one of the single-element transducers used in this study was damaged during the course of the project, and its reduced performance was easily detected by the thermochromic DQA. The thermochromic phantom thus appears to be an effective way to check the proper functioning of focused ultrasound transducers.

To explore the phantom's sensitivity to transducer deficiencies during measurements of the Insightec Neuro transducer, sonications were performed with the central 256 elements ("Apodized" data in Fig. 6) and will the full 1024 elements in the hemispherical array ("Unapodized" data in Fig. 6). The Insightec software is designed to distribute the requested power across the available elements in the array. By enabling the outer 768 array elements (and keeping power fixed), the acoustic focus becomes shorter in the axis of propagation and therefore extends a shorter distance into the absorber, having the effect of increasing the distance of thermal diffusion until the thermochromic liquid crystal is heated. The data are in broad agreement with bioheat simulation of this scenario, both indicating a perceived reduction in the thermal output in the lengthened time to first color change.

Interestingly, our measurements highlighted a set of odd data. A set of measurements (Fig 6, in blue) were supposed to be performed without apodisation but appeared to align with the apodised data. The thermochromic setup appeared to be sensitive enough to the shape of the ultrasound beam to allow us to detect an error in the emission settings. This error was confirmed by looking at the logs of the Exablate neuro that automatically record all the parameters that have been used at each sonication.

The thermochromic phantom is inexpensive and straightforward to establish and set up in a lab. The cost of the materials are summarized in Table 2.

The 3mm sheet of AptFlex F28 was cut with scissors into smaller pieces for use in the thermochromic phantoms.

It is reasonable to assume that a piece with length and width of 2.5 cm would suffice (for transducers with diameter or focal spot diameter less than 2.5 cm), so the 300x300x3 mm sheet of F28 could be cut into greater than 100 pieces of 30x30x3 mm absorbers. Since the sheets of Hallcrest TLC are slightly larger in length and width than the AptFlex F28, the same cutting idea was applied. The cost per phantom is estimated to be seven dollars.

Furthermore, the TLC material used here changes color reversibly and can be re-used. In this paper we used two phantoms for the single element measurements: one for each element. Alternatively, F28 could be replaced by a cheaper rubber. This would bring the cost of the phantom down drastically. However, the acoustic and thermal properties of the rubber may need to be measured in order to use the phantom with confidence.

For larger apertures, a 3D printed or other custom holder is recommended in order to limit the amount of error a user can produce from assembling the setup. A holder ensures the phantom is at a fixed distance from the transducer so that the pressure at the absorber face closest to the transducer can assumed to be constant for each trial.

The thermochromic phantom calibration method is straightforward and relatively rapid. Once the user has the phantom, ultrasound gel, and holding mechanism (either a 3D printed holder or some equivalent) calibration can be done in minutes. If calibration with one sonication at high power is sufficient for an experiment, then the calibration can be done in about a minute -the user places the transducer in the holder, puts ultrasound gel at interfaces, turns on the transducer, then waits to see a color change. The color-change time for a high power sonication would range from 2-10 seconds. If more data points are necessary, the user has multiple options. If the user has more than one phantom, they can switch out the gel pad and phantom and sonicate again without waiting between trials. This would ensure that the initial state of the system is the same between sonications. If the user only has one absorber, they can either run the phantom under cold water or wait 5-10 minutes between sonications to roughly ensure that the state of the system is the same between sonications. The temperature in the absorber could be estimated by an infrared thermometer or by a thermocouple if necessary. Regardless, the calibration time frame is in the order of minutes and not hours.

Overall, this work is a preliminary proof of concept and further developments are needed to assess the sensitivity of the approach. Different potential sources of errors could affect the approach.

Air bubbles are a possible source of error. They can form inside the gel pad, in between interfaces coupled by ultrasound gel, and in between the TLC sheet and the absorber. Air bubbles in the gel pad near the transducer surface can potentially damage the transducer. Air bubbles in the gel pad or near an interface may interfere with the ultrasound propagation and alter the lesion time recorded by the user. Air bubbles in-between the TLC sheet and absorber may alter heat transfer at the air interface and alter the lesion time recorded by the user. Fortunately, air bubbles in the gel pad are easy for the user to see and seldom form. Air bubbles in-between the TLC and absorber are harder to see.

Ambient temperature

The ambient temperature of the room may change over time. Since the TLC shows absolute temperature change, a change in the ambient temperature would change the lesion time recorded by the user.

Heterogeneity in Materials.

It is possible that the properties of the TLC sheet or the acoustic absorber are heterogeneous. Since each phantom is created by cutting away one section from each material, it is possible that each phantom could have different thermal or acoustic properties due to heterogeneities in the materials.

Bonding Difference.

It is possible that the bond between the TLC sheet and the acoustic absorber is different for different phantoms.

If so, this would mean that heat would conduct from the absorber to the TLC sheet at different rates for different phantoms. This implies that the color-change time would be different depending on which phantom is used.

Prior work on thermochromic phantoms was in the form of a traditional, volumetric ultrasound phantom [START_REF] Negussie | Thermochromic tissuemimicking phantom for optimisation of thermal tumour ablation[END_REF]. While this phantom proved useful for determining the volume of tissue mimic that experienced a peak temperature of 60 °C, it was a single-use phantom because of the permanent nature of the change in color of the thermochromic ink used. While the phantom we present is not volumetric, the ability to reuse the phantom dozens of times contributes to the consistency / longevity of the measurement accuracy and the low relative cost of the phantom. Such an approach could benefit to recent low cost devices developed for central nervous system applications [START_REF] Biase | Transcranial Focused Ultrasound (tFUS) and Transcranial Unfocused Ultrasound (tUS) Neuromodulation: From Theoretical Principles to Stimulation Practices[END_REF]Chen et al. 2020b;[START_REF] Chen | Neuronavigation-guided focused ultrasound for transcranial blood-brain barrier opening and immunostimulation in brain tumors[END_REF][START_REF] Ferri | On the Evaluation of the Suitability of the Materials Used to 3D Print Holographic Acoustic Lenses to Correct Transcranial Focused Ultrasound Aberrations[END_REF]Maimbourg et al. 2018a;Maimbourg et al. 2018b;Maimbourg et al. 2018c;[START_REF] Pouget | Neuronavigated Repetitive Transcranial Ultrasound Stimulation Induces Long-Lasting and Reversible Effects on Oculomotor Performance in Non-human Primates[END_REF][START_REF] Rabut | Ultrasound Technologies for Imaging and Modulating Neural Activity[END_REF][START_REF] Schafer | Design, Development, and Operation of a Low-Intensity Focused Ultrasound Pulsation (LIFUP) System for Clinical Use[END_REF][START_REF] Wu | Efficient Blood-Brain Barrier Opening in Primates with Neuronavigation-Guided Ultrasound and Real-Time Acoustic Mapping[END_REF][START_REF] Zou | Ultrasound Neuromodulation Inhibits Seizures in Acute Epileptic Monkeys[END_REF], low back pain treatment [START_REF] Aginsky | Tolerability and Feasibility of X-ray Guided Non-Invasive Ablation of the Medial Branch Nerve with Focused Ultrasound: Preliminary Proof of Concept in a Pre-clinical Model[END_REF]Chen et al. 2020a;[START_REF] Perez | Fluoroscopy-Guided High-Intensity Focused Ultrasound Neurotomy of the Lumbar Zygapophyseal Joints: A Clinical Pilot Study[END_REF] or varicose veins (Barnat et al. 2020a;Barnat et al. 2020b;[START_REF] Obermayer | Extracorporeal Treatment with High Intensity Focused Ultrasound of an Incompetent Perforating Vein in a Patient with Active Venous Ulcers[END_REF].

CONCLUSION

We demonstrated on two transducers and two DQA phantom geometries that the thermochromic phantom approach is an effective way to assess the stability of focused ultrasound transducers. Next steps involve refining the simulation methods used to better validate the DQA design and deploying custom-mounted DQA phantoms to end users. For example, FUS mobile (Alpharetta, GA, USA) started using such a thermochromic phantom for targeting confirmation and for validation of the full functionalities of the Neurolyser platform.

  Two identical SonicConcepts transducers (Model XDR058, Sonic Concepts, Bothel, WA) with 1.1 MHz acoustic working frequency were used to sonicate the phantoms.A waveform generator (DG1022, Rigol, Beaverton, OR, USA) supplied a 1.1 MHz sinusoid to a power amplifier (100A100, Amplifier Research Corp., Souderton, PA, USA). The amplified signal was monitored by an oscilloscope (TDS2001C, Tektronix, Beaverton, OR, USA) and connected to the input of a Sonic Concepts matching network and ultrasound transducer. The impedance of the transducer with matching network was measured with a network analyzer (4194A Impedance/Gain-Phase Analyzer, Yokogawa Hewlett Packard, Tokyo, Japan) and was used to calculate the electric power applied to the transducer based on the measured voltage.The time until first visible TLC color change was recorded for electrical power in the range 0.4-12 W.Temperature was monitored by Copper-Constantan thermocouples (Type T, 0.010 in diameter, Omega, Stamford, CT, USA) attached to the side of the transducers housing in order to ensure the temperature returned to baseline before performing a new measurement. Recording was achieved by an acquisition board (TC-08, Pico

FIGURESFigure 1 :

 1 FIGURES

Figure 2 .Figure 3 :

 23 Figure 2. CAD drawing of the geometry of the phantom. X-, Y-, and Z-axes are red, green, and blue,

Figure 4 :

 4 Figure 4: (a)Pressure field generated by the XDR058 single element transducer, measured by a needle

  20

Figure 5 :

 5 Figure 5: Time to see a color change with the single element transducer. Error bars indicate standard deviation

Figure 6 :

 6 Figure 6: Time to see a color change with the Exablate Neuro system operated at 670 kHz.

ACKNOWLEDGEMENTS

This work was supported by the Focused Ultrasound Foundation.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare. 1) 7500 ( 5) 997 (8) 1.21 (9) Speed of sound (m/s) 1500 ( 1) 4350 ( 5) 1523 ( 8) 345 ( 4) Attenuation (Np/m/MHz) 330 (1) 2.5x10 -3 0 (7) 18.9 (3) Specific heat (J/kg*K) 1800 ( 2) 420 ( 6) 4310 ( 7) 1006 (9) Thermal conductivity (W/m*K) 0.165 (2) 1.45 (6) 0.603 (7) 0.025 (9) (1) Provided by manufacturer in the datasheet: Precision Acoustics Ltd (Dorchester, Dorset, UK)