Simulating airborne ultrasound vibrations in human skin for haptic applications

The ability to focus high amplitude ultrasound onto human skin and elicit a tactile sensation in mid-air poses a number of scientific opportunities and challenges. For example, a common issue with mid-air haptics is the relatively low forces generated on the user’s skin, thus limiting the dynamic range of vibrotactile sensations that can be induced. To that end, we develop a viscoelastic rheological finite element model and apply it to show that different acoustic radiation patterns can generate varying forces, displacements, and oscillatory patterns on the targeted skin and its nearby vicinity. Our framework allows for the high resolution computational exploration of a wide range of mid-air stimulus sensations saving time and resources spent on running expensive experiments while.


I. INTRODUCTION
Mid-air haptic technology often refers to the use of focusing airborne ultrasound waves to remotely induce vibrotactile sensations on the human skin [1], [2]. The applications of this technology are wide-ranging and include augmented and virtual reality (AR/VR), robotic teleoperation, humancar interfaces, digital signage and interactive kiosks (see recent review article [3]). Moreover, many of the advancements relating to this emerging field relate to hardware and software improvements in the form of new device prototypes capable of producing higher acoustic pressures at improved spatial and temporal resolutions [4]- [7]. On the other end of the spectrum, psychophysical studies have advanced our understanding of how different mid-air haptic stimuli are perceived by users to convey shapes, textures and even human emotion [8]- [10]. The bridge between these two ends of the spectrum is supported by a somewhat limited physical understanding of how focused ultrasound waves interact with the human skin, a topic that has recently been studied using laser Doppler vibrometry (LDV) experiments, both in vivo and in silico using silicone phantoms [11]- [13].
This paper aims to strengthen this bridge through the use of experimentally calibrated computer simulation tools.
Interest in skin mechanics and wave propagation, both on and through the skin has steadily proliferated in recent years, mostly due to advancements in robotics and humancomputer interaction (HCI) applications. In the former, roboticists can gain intuition about how to produce secure grasping mechanisms for dexterous manipulation by understanding the frictional forces applied by human fingers [14], [15]. In the latter, haptic engineers can design and program novel devices that can produce tactile experiences, finding application in a variety of sectors including automotive, education, health, VR and AR [16].
Meanwhile, recent studies have shown that tactile information is not restricted to the points of contact of the tactile stimulation, i.e., the grasping area of the fingertips, but rather that tactile information propagates through the whole hand at surface wave speeds of 5-7 m/s and that the vibration signals intrinsically encode sufficient information about the interaction type and the touched object's properties [17]. The running hypothesis is that these signals are then decoded by a distribution of receptors embedded in our skin [18] and then spatiotemporally integrated and frequency-filtered by our somatosensory system to inform and ascribe meaning about the world around us [19].
To that end, a major challenge in haptic display engineering is to find methods for efficiently stimulating this large continuous medium, the human skin, with a variety of haptic devices. Moreover, to do so effectively, one needs to take into account the viscoelastic properties of the skin which display both viscous and elastic responses [20]. Recall that viscosity is related to the delayed recovery from deformation, like a buffer that counteracts velocity of the external force resulting in a delay or damping of movement. Equally, elasticity is related to the rebounding and quick recovery from a deformation, like a spring that produces a force proportional to the deformation.
While numerous haptic interfaces and devices have been developed over the years to deliver a variety of tactile stimuli [21], these are mostly either wearable [22] or embedded in other devices (e.g., controllers, smartphones, touchpads and touchscreens) and require direct physical contact. In contrast, non-contact (touchless) haptic interfaces can deliver vibrotactile sensations remotely through the use of focused ultrasound waves, lasers or air jets. The latter suffers from time delays due to relatively low wind speed, and is diffusive or even possibly turbulent and thus cannot project high-resolution haptics at range. In contrast, mid-air haptic devices such as those provided by Ultraleap operate at 40 kHz airborne ultrasound and use 256 piezoelectric transducers (usually Murata MA40S4S) that can be electronically controlled to focus acoustic energy to a point and generate RMS pressures that can reach up to 160 dB SPL (relative to 20 µPa), i.e., 2000 Pascals, at 20 cm away from the device centre. The focal volume has an ellipsoid shape with a −6 dB cross-sectional diameter of full width half maximum FWHM = 14.5 mm, exerting a localised force onto an acoustically reflecting surface (e.g., skin or skinmimicking silicone) of about 10.8 mN [23] and generating a peak to peak displacement of about 1 µm [12]. Multiple focal points can be created and distributed in space to form shapes and then amplitude modulated (AM), i.e., pulsed on and off at some frequency f AM as to stimulate the Pacinian corpuscle (PC) mechanoreceptors that respond to vibratory stimuli in the frequency range of 20-1000 Hz, however are most sensitive at f AM = 200 − 250 Hz [24]. Sophisticated algorithms can spatiotemporally modulate (STM) those focal points and rapidly update their locations (e.g., at 40 kHz update rate) along open or closed paths that repeat at some frequency f STM in order to stimulate a response from the receptors that the STM path overlaps with [25], [26].
Despite all these efforts, experiments are unable to sample the huge stimulus space available to ultrasonic mid-air haptic displays that consist of multi-point, multi-path, multifrequency dimensions. Physical experiments are difficult, expensive, and time-consuming to perform.
In this paper, we advance our understanding of mid-air haptics by developing a computer simulation model that predicts the skin-response to any mid-air haptic stimulus. We experimentally calibrate our model through LDV data, and then use the model to explore a variety of stimuli and response properties. The three main contributions of the paper are: 1) we develop an experimentally validated viscoelastic rheological skin model, 2) we port this model into a finite element computer simulator, and 3) we explore and analyse different mid-air haptic stimuli effects using the model. The paper is structured as follows: In Sec. II we describe how a phased array of ultrasonic transducers is used to generate a focusing acoustic field, and how it is mathematically modelled. In Sec. III we describe how a high intensity acoustic focus impinges a force onto an absorbing or reflecting material and also validate our model through experimental measurements using a precision weighing scale. In Sec. IV we detail the dynamic rheological model that we have constructed to describe the behaviour of viscoelastic materials such as skin and skin-mimicking silicone. In Sec. V we describe how to bring all the previous components together in a COMSOL simulation environment, i.e., the acoustic load and the solid mechanics, whiles also describing the data collection methods for post-processing. In Sec. VI we simulate three mid-air haptic stimuli that are validated against experimental data from previous LDV measurements [12] and also describe how these relate to standard rheological tests. In Sec. VII we describe a number of parameter investigations and explorations enabled by our constructed multiphysics simulation model of mid-air haptic stimuli and their effects on viscoelastic surfaces. Finally, in Sec. VIII we summarise the significance of our results and discuss future work.

II. ACOUSTIC FIELD MODEL
In basic terms we can think of the haptic system as follows: 1) a phased array of transducers emits a time varying ultrasonic pulse at a carrier frequency of 40 kHz, 2) as these waves propagate they constructively interfere and create a high amplitude focal point at z ∈ R 3 , 3) this focal point is usually targeted onto the surface of a human user's skin (e.g., fingertip or palm), 4) the acoustic pressure at z is high enough to deform the skin and elicit a tactile sensation. A full physical model capturing a mid-air haptic device and the acoustically induced physical effects on the skin is a significant challenge and therefore linear acoustic models are used as 'good enough' approximations. Advanced acoustic simulation tools leveraging the finite element method (FEM) or pseudo-spectral methods such as k-Wave can be used to model the acoustic fields generated by transducer arrays, and can even capture non-linear acoustic effects which may become non-negligible at very high levels (e.g., greater than 145 dB SPL). Time domain simulations are even more complex and computationally expensive and can become prohibitive for large domain sizes as they need very fine meshes and many degrees of freedom. To simplify these requirements and make progress we will model the acoustic output of an array of transducers as baffled pistons and apply Huygens principle of superposition to construct the crosssectional pressure profile of a focal point. We shall then convert this into acoustic radiation pressure, and later use that to calculate the applied force on a viscoelastic surface.
We start by describing the acoustic field at some sampled location x produced by N electronically controlled ultrasonic transducers that are arranged on a flat surface centred at the coordinate origin facing up towards the z axis. To do so, one can start by calculating the complex pressure P i (x) at a point x due to the ith emitter located at point y i on the array surface using where j = √ −1 and P 0 is a constant that is defined by the transducer amplitude, d(x, y) is the Euclidean distance between points x and y, the transducer directivity function is defined by 2 J 1 (ka sin θ ) ka sin θ , where J 1 (x) is the Bessel function of the first kind, k = 2π λ is the wave-number, λ is the wavelength of the ultrasonic signal, a is the transducer radius, θ is the polar angle between points x and y i , and φ i is the phase delay applied to transducer i. To calculate the total peak pressure P T (x) in Pascals generated by N transducers at (x), one must compute the summation of the contribution of each transducer and take the absolute value to obtain To achieve a single focusing at some location z, the individual phase delays φ i of each transducer must be chosen such that they counterbalance the difference in distances from each transducer to the focal point position z and thus always constructively interfering at that point Using this formulation one can simulate and plot the total pressure P T for a whole region in space. Fig. 1a) shows the acoustic field near a focal point targeted at z = (0, 0, 20) cm simulated using the HandyBeam open-source software [29]. It is observed that the total acoustic pressure reaches a maximum of 2 kPa (160 dB SPL), and drops away rapidly within ±10 mm. Note that this is the RMS pressure, not the peak pressure P T (z) which can be obtained by multiplying 2 kPa by √ 2. It is known that the tightness of the focus depends on several aspects including the aperture of the array, the inter-transducer spacings, and the distance of the focus from the array. For the purposes of simulating the acoustic radiation pressure exerted onto a viscoelastic surface, we first approximate the pressure distribution of the focal point by a Gaussian function with peak pressure P T (z) Pascals and width w meters where r is the radial distance on the x − y slice of the focal point location z, and w = FWHM 2 √ 2 ln 2 . Thus, the acoustic pressure distribution shown in Fig. 1a) is very well represented by P RMS (r, z) with w = 6.12 mm at z = {0, 0, 20} cm. We note that non-linear acoustic effects such as waveform distortion will likely come into play at high enough pressures however these are beyond the scope of our present model and investigations. Further, we will ignore any additional fringes of acoustic pressures surrounding the focal point such as grating lobes, thus enabling us to model multiple focal points as a collection of Gaussian profiles added together.

III. ACOUSTIC RADIATION FORCE
When a focused acoustic wave encounters another medium, e.g., from air to a viscoelastic surface, some of the wave VOLUME 4, 2016 will reflect and some of it will be transmitted. The reflection coefficient R ∈ [0, 1] can be determined by their respective acoustic impedances. Importantly, for mid-air haptic applications, due to the large impedance mismatch between the air and skin interface we have that R = 99.9%, implying that the effective acoustic pressure experienced by the skin when targeted by a focal point is actually almost double that of free space, especially if the wave propagation direction is perpendicular to the skin. During this interaction, momentum must be conserved, and therefore a force F is exerted onto the surface. To calculate this force we first calculate the Langevin radiation pressure P R given by [30] where c 0 = 343 m/s is the speed of sound in air, ρ 0 = 1.225kg/m 3 is the density of air, and α = 1 + R 2 . Note that since P RMS is a Gaussian, then P R will also be Gaussian shaped, but narrower due to the squaring in (5). Finally, to obtain the acoustic radiation force we simply integrate P R over the target surface where r 0 is the radius of the surface area being irradiated by the focal point. Note that when r 0 is large then F R ≈ α 2 P T (z) 2 4ρ 0 c 2 0 πw 2 Newtons. To validate equation (6) we have conducted a series of experimental measurements using a precision force balance (Kern-PCB 2500-2) with an accuracy of 0.01 g. The experimental setup is shown in Fig. 2. An Ultraleap USX device is held upside down by a stage and is programmed to focus at z = (0, 0, 20) cm. The precision scale reads the force applied by the focal point onto a Perspex cylindrical pillar of diameter 2.1 cm which sits on its plate. The plate is shielded from non-focused ultrasound waves by a foam board that has a circle cutout to fit the pillar and stands on the table therefore is independent of the scale. The pillar has a high impedance mismatch resulting in α ≈ 2 that perfectly reflects any incident acoustic waves.
Through the Ultraleap Application Programming Interface (API), the device is programmed to create a focus at different intensities from 0 to 1, at 0.1 increments, with intensity 1 corresponding to a simulated RMS value of 2 kPa at the focus. Note that the focus is not modulated in any way. The focus is created for 20 seconds and then turned off for 40 seconds (cool-off period) before moving on to the next intensity level. To enable the force balance to reach a steady state, a reading is taken at approximately 10 seconds, i.e., the midpoint of the focus time. Each sample point is repeated 5 times before an average is taken. The results FIGURE 2: Force balance measurement apparatus. a) An Ultraleap USX device is held upside down by a stage. The array is programmed to focus at z = (0, 0, 20) cm where we also place a Perspex cylinder of diameter 2.1 cm that is standing on the force balance plate. A foam board with a 2.2 cm diameter hole cut-out is positioned as to shelter the balance from non-focused waves therefore recording only the force applied due to the focus. b) Shows the experimental results (blue points), and the analytical prediction of equation (6) (orange points). The device API intensity (x-axis) linearly maps from 0 to 2 kPa of RMS acoustic pressure. are plotted in Fig. 2 for a range of intensities reporting a maximum force of 10.89 mN measured when the focal point pressure output from the device is expected to be 2 kPa RMS, or equivalently 160 dB (SPL). This force value is appropriate and similar to ones previously reported in the literature.
The analytical prediction of the acoustic radiation force (6) (orange points in Fig. 2) is underestimating the experimental results, reaching a maximum of 6.25 mN. We attribute this discrepancy to non-linear acoustic effects such as acoustic streaming which can generate airflow and therefore a downward force onto the plate. Pittera et al. [31] have recently reported airflow speeds of s = 2.9 m/s generated by a 2 kPa focal point. We can use this information to calculate the corresponding acoustic streaming force F S applied to the Perspex pillar and plate by the air jet using the formula F S = ρπr 2 0 s 2 . We find that for jet speeds of s = 2.9 m/s the F S = 3.82 mN, which is almost exactly the difference between the measured force 10.89 mN and the calculated F R = 6.25 mN. We therefore highlight in Fig. 2 the contributions to the total measured force. By comparing these regions we observe that F S also grows quadratically with pressure but becomes noticeable only after an API intensity of 0.3 (i.e., a focal point of 600 Pa RMS) after which the ratio F R /F S ≈ 1.45. Having described and experimentally validated the acoustic radiation force applied by a focused ultrasonic field onto a reflective surface, we next proceed to investigate the response of a viscoelastic skin-mimicking material by constructing a rheological model that will capture the timedependent behaviour of the material under the influence of different localized and time varying stresses.

IV. RHEOLOGICAL MODEL
The term viscoelastic is derived from the words "viscous" and "elastic"; a viscoelastic material presents both viscous and elastic behaviour and can be described through rheological models [32]. The linear elastic behaviour is modelled by a spring that follows Hooke's law, relating the stress applied σ to the resulting strain ε through σ = Gε, where G is the elasticity modulus. The viscous behaviour is modelled by a damper that follows the constitutive law σ = ηε whereε = dε dt is the rate of change of strain, and η is the material viscosity. Combining these components in series or in parallel makes it possible to describe more complicated viscoelastic behaviours. For example, the Maxwell model consists of a spring and a damper connected in series, while the Kelvin-Voigt model has them in parallel, both of which are commonly used in the literature. However, the Maxwell model cannot describe creep or recovery, and the Kelvin-Voigt model cannot describe stress relaxation [33]. Both creep and relaxation mechanical tests are important and present in vibrotactile simulations of the skin. We therefore resort to using a Standard Linear Solid (SLS) rheological model in its Maxwell representation (SLSM) consisting of a spring and a Maxwell body connected in parallel [34] schematically shown in Fig. 3. Note that the SLSM arrangement is the simplest model that predicts both creep and recovery phenomena and has previously been a good fit for skin-mimicking viscoelastic materials [35].
The viscoelastic components relevant of the SLSM model are characterised by the low-frequency modulus G 1 , the high-frequency modulus G 2 , and the viscosity η 3 . Moreover, its constitutive relations of the individual components are from which we can derive the constitutive equation of the system in Fig. 3a) as follows Rearranging, we arrive at the constitutive equation of the SLSM system with a relaxation time constant of τ = η 3 G 1 . We note however that the elastic and the viscous properties of all solid materials like the one described in (9) are to some extent dependent on the frequency of the dynamic load being applied, i.e., the high-frequency modulus G 2 value increases with the stimulus frequency which for mid-air haptics is the carrier frequency of the ultrasound f = 40 kHz and can also be directly related to the rate of change of the strain f =ε ε . Therefore, we substitute σ = G d ε andσ = G dε into (9) to get where G d represents a frequency dependent elastic (dynamic) modulus that acts as a transfer function between the applied stress σ and the resulting strain ε. Dividing both sides by ε and identifying f =ε ε we finally arrive at which can be rearranged to get a formula for calculating Note that equation (12) essentially increases from G d = G 2 at f = 0 (static case), and asymptotically approaches is the big O notation describing the limiting behaviour of G d for large values of f . The modified constitutive equation of our rheological model with a dynamic elastic modulus therefore reads and is implemented in COMSOL in the next section. In our experimental and numerical simulation studies that follow we used an Ecoflex phantom sample material (Ecoflex-0010) as it is simpler than human skin and its mechanical behaviours of interest are very similar to human skin [36], [37]. The parameters used are taken from [35] and equation (12)

V. FINITE ELEMENT SIMULATION IMPLEMENTATION
A three-dimensional FEM model was constructed in COM-SOL Multiphysics 5.5 capturing all previous sub-models described in this paper. Namely, the acoustic radiation pressure due to a high intensity ultrasound focus was modelled as a Gaussian shaped indenter P R (5) and the rheological model describing the viscoelastic properties of a skin-mimicking phantom was modelled through the constitutive equation (13). The corresponding simulations were then performed on a high-performance computer. The implementation of the FEM has three distinct aspects: 1) the acoustic load stimulus, 2) the material representation (constitutive laws) and the boundary conditions (loading and restraints), and 3) the data collections and analysis. Therefore, we have structured this section along these three main points: acoustic load implementation, solid mechanics implementation, and data collection.

A. ACOUSTIC LOAD IMPLEMENTATION
The boundary load condition is the application of a force or pressure towards a specific domain. The load amplitude, direction, shape, and dynamics are all different parameters that needed to be formulated in the COMSOL environment.
To that end, we transferred the acoustic radiation pressure distribution P R derived in Sec. II into a dynamic domain that would modulate both the intensity and location of the focus such that the applied pressure on a target surface at location x and time t due to a focal point at z(t) is given by where A(t) ∈ [0, 1] modulates the acoustic pressure of the focus from 0 to 1, while z(t) captures dynamics of the focus position at time t ≥ 0 as it moves along a specific path. Basically, equation (14) describes the applied pressure P onto a target surface at each point x and time t, due to an acoustic load stimulus (the focal point) that is amplitude modulated (AM) by A(t), and spatially modulated by changing its focus location z(t). Equation (14) can be used to define a variety of stimuli and can be directly implemented into COMSOL. Three stimuli that we have used as benchmarks for our investigations in the following sections are In the third case (STM), we use A(t) = 1 and z(t) = (r 0 cos(2π f STM t), r 0 sin(2π f STM t), 20) resulting in a non-localised vibrotactile point that is moving round a circular path of radius r 0 cm at a fundamental frequency of f STM times a second, with harmonics at integer multiples n × f STM . The surface speed of the focal point can be calculated through u = 2πr 0 f STM cm/s and is an important parameter since displacement appears to be maximal when u equals the speed of surface wave propagation of the viscoelastic target [11]. More complex stimuli can be defined using the generative equation (14) having non-zero f AM and f STM .
In COMSOL, we have implemented the dynamic load using a Boundary Load condition to apply the desired force distributed on the silicone model. This condition permits to choose between force or pressure and direction. Adding more boundary load conditions permits the creation of multi-focal point stimuli, however this requires some care since one needs to consider power limitations of the physical devices which will have to distribute their total acoustic pressure output capacity between multiple focal points. Also, overlapping focal points in the simulation should not amplify one another.

B. SILICONE SOLID MECHANICS IMPLEMENTATION
A three-dimensional FEM was built in COMSOL using the Structural Mechanics Module to simulate the excitations and vibrations resulting from the acoustic load stimulus. A homogeneous viscoelastic material was defined in the shape of a thin cylinder and an SLSM rheological model was created like in Fig. 3 with one purely elastic branch and one viscoelastic using the material parameters described in Sec. IV and Tab. 1 as to mimic Ecoflex which resembles skin behaviour well enough. Appropriate boundary conditions and mesh densities were also implemented and the geometries were taken from that used in the LDV experiment described in [12], namely 120 mm in diameter and 20 mm thickness.
Suitably fine meshing is a requirement for the accuracy of the results in the FEM model but also for computational efficiency. To ensure a suitably fine representation of our model, a mesh sensitivity study was performed resulting in a recommended maximum element size of 1 mm in the tangential dimension (x − y plane), and 2 mm through the thickness of the sample (z direction). Similarly, to ensure accurate time stepping we conducted time step convergence studies and used a time step size of at least 1/ 40 max( f AM , f STM ) seconds. So for an AM stimulus of 200 Hz, the simulation time step size is 0.125 ms.
Finally, modelling the propagation of surface waves on a large vibrating structure can be a challenging task. It requires balancing the reduction of the computational domain's size with the decrease of reflection at surface boundaries due to attenuation. Therefore, we use the "Low-Reflecting Boundary" condition along the perimeter of the cylinder to remove any reflections from the domain edges and reducing the computational domain that needed to be simulated. All simulations were run on a high performance desktop PC.

C. DATA COLLECTION AND PROCESSING
We employ two approaches towards data collection and processing. In the first instance, we store the full dynamic mesh excitation matrix of the viscoelastic surface, thereby allowing us to produce 3D animations and displacement heat-maps like the one shown in Fig. 3c). In the second instance, a single probe is inserted into the surface of the model at some location x = z(t = t 0 ) from which we record the surface displacement time-series during a certain time window t ∈ [t 0 ,t 1 ]. Taking the Fourier transform of the displacement time-series we can then obtain a frequency representation of the surface excitations at x. This helps identify how input frequencies like f AM and f STM transfer over onto the viscoelastic medium. Generally, it is believed that a stronger tactile response will be perceived when the frequency spectrum of the viscoelastic medium is matched with that of the embedded mechanoreceptors [24].
Another key observable is the peak-to-peak (P2P) displacement at x, allowing for a Peak-and-Trough analysis highlighting differences between the various stimuli applied (UM, AM, STM, etc.). Generally, it is believed that a larger P2P will lead to a stronger tactile response. Many other results can be extracted from the suggested model like the velocity of the surface waves and the strain in the tangential directions (in-plane). Finally, we remark that the displacement measured by the probe at x is due to two physical effects, 1) the AM or STM of the focal point when x = z(t), and 2) the surface wave interference or propagation effects generated by the focal point modulation when x ̸ = z(t). The former is due to the localized tactile stimulus, while the latter is a non-local tactile effect and is much harder to predict.

VI. MODEL VALIDATION
Chilles et al. [12] conducted LDV experiments on an Ecoflex sample while being excited by a focal point generated by an Ultraleap device. The focal point was either stationary and UM, stationary and AM with f AM = 200 Hz, or an STM circle with f STM = 70 Hz and r 0 = 1 cm, as described in Sec. V-A. These three different stimuli are replicated in our COMSOL simulation environment and are used here to validate our FEM model through a comparison with the experimental LDV data. In each instance, a probe is placed on the surface of the simulated model, either at the centre of the focal point position, or along the STM circle path. The probe records the time series vibrations of surface displacement which are compared to the corresponding LDV measurements. A Fourier transform is also performed to compare the excitation spectra. These results are reported in Figs. 4, 6, and 9. We will discuss each of these cases separately with further analysis in the following subsections.
Tab. 2 summarizes the response of the silicone sample to the different stimuli and presents a comparison of the P2P and RMS surface displacements. The table also compares the experimental and the simulation results with previous attempts at a FEM simulation [12] denoted as the static model. It is noted that the dynamic theological model introduced in this paper matches well with the P2P and RMS excitations measured in LDV experiments. In contrast, the static model displays a noticeable error between the simulation and experimentally measured displacements and this is because it did not benefit from the dynamic modulus derivation presented in equation (12) and its inclusion into the modified constitutive equation (13).

A. UNMODULATED SCENARIO
This case describes an unmodulated focal point described by A(t) = 1 and z(t) = (0, 0, 20) cm. The comparison results of the displacement time series are shown in Fig. 4.a) and are observed to be quantitatively similar. Note that the time series displacement at the centre of the stimulus resembles that of an underdamped system to a step function. Namely, the displacement is characterised by an initial overshoot, before settling into a steady state value. Even if the initial overshoot peak is about 20% higher in the LDV measured data than in the COMSOL simulation, the steady-state value matches very well. Fig. 4.b) reports on the Fourier transform of the displacement generated over time.
As expected, there is no clear excitation frequency since the stimulus is unmodulated. In addition, a 40 kHz peak is only just noticeable in the LDV measurements, but not in the simulated data. This is expected since the silicone material and our skin shows almost no response to the 40 kHz component, which is also the main reason why we do not simulate such a high-frequency component in COMSOL, the other reason being that it would be computationally very expensive. It is worth noting here that the unmodulated scenario is very similar to a creep test. The creep test consists of measuring the time-dependent strain ε(t) resulting from the application of steady uniaxial stress σ 0 . The creep test helps understand two situations: how long would it take until a sample starts deforming at a constant rate, and once the applied stress is removed, how long does it take for the material to recover? The LDV and simulated data from the UM scenario can help address the first part of this question.   retarded elastic deformation (creep strain), and c) steadystate viscous flow. These three regions are also present in the experimental data shown in Fig. 4.
In the SLSM model, the creep test equation can be obtained through a Laplace transform and solved in the case when the applied stress is a simple step function. The UM stimuli however is not a step function but rather a smooth but rapid build up of a Gaussian shaped radiation pressure. Therefore, we notice some oscillations before the steadystate is reached. Despite this, the characteristics of a creep test can be deduced from Fig. 5. Firstly, the instantaneous elastic response or the initial strain is of 2.4 × 10 −5 and duration of 0.1 ms. Then, the viscoelastic region where the strain starts at 2.4 × 10 −5 to 7.5 × 10 −5 with a duration of 0.4 ms. This could be interpreted as a fundamental timedelay limitation imposed here by the physics of the skin onto mid-air haptic technology. Finally, at the steady-state region the strain settles at 5.25 × 10 −5 after 25 ms from the beginning of the test.

B. AMPLITUDE MODULATION SCENARIO
The AM stimulus resembles a dynamic indentation test. Here, the pressure applied to the surface alternates in a sinusoidal fashion due to the acoustic pressure which is modulated through A(t) = 1 2 1 − sin(2π f AM t) while keeping the focal point stationary at z(t) = (0, 0, 20) cm resulting in a localised static vibrotactile stimulus with a fundamental stimulus frequency of f AM and a harmonic at 2 f AM . Fig. 6 demonstrates the displacement over time for a focal point that is AM at a frequency of f AM = 200 Hz. The corresponding frequency spectrum is also shown, both of which are in very good accordance with the experimental LDV data, both capturing the fundamental frequency and the first harmonic at 2 f AM , thus indicating a strong coupling between driving frequency and surface response. The main discrepancy between the simulated and experimental plots is the under estimation of displacement rebound which does not rebound out of the surface and peaks at 0 µm after a few cycles. In contrast, the experimental data rebounds out of the surface height at rest to reach +0.25 µm in between of the pressure indentations down to −0.95 µm. This is however a relatively small discrepancy with both time-series and spectra being quantitatively similar.
Interpreting the AM stimuli as a dynamic indentation test enables some further viscoelastic investigations. Namely, by expressing the time-dependent stress as (15) and the resulting time-dependent strain as we can see that the strain ε is oscillating with the same frequency content as the stress but lags behind by a phase angle δ ∈ [0, 2π] due to the response time of the viscoelastic material. This angle is often referred to as the loss angle of  the material [33]. A loss angle of δ = 0 is characteristic of elastic behaviour where stress and strain are in-phase to each other, while δ = π/2 is characteristic of viscous behaviour where they said to be are out-of-phase. Moreover, tan(δ ) is known as the mechanical loss or the loss factor and can be considered as a measure of the material's damping effect, i.e., its ability in restraining vibratory motion. Fig 7 uses our COMSOL model to simulate and plot the applied stress and resulting strain during the AM condition, thus demonstrating a small phase off-set observed by the material.
To measure the effect of this phase off-set, we first consider equations 15 and 16 as a parametrization in time t. The resulting stress-strain curve forms an elliptic hysteresis loop as seen in Fig. 8a), the enclosed area of which can be interpreted as the energy dissipated during each cycle, through internal friction or heat. The constructed model thus allows us to investigate all these aspects during a sweep of different dynamic indentation frequencies f AM and identify any preferential stimulus frequencies.
The normalized energy loss after such a frequency sweep is plotted as a function of f AM in Fig. 8b). We observe a monotonically increasing curve indicating that higher AM frequencies introduce a growing phase-offset and thus a larger hysteresis loop indicating that more and more energy is dissipated during each cycle due to damping. Also in Fig. 8b), we plot the simulated maximum dis-placement into the surface generated by different f AM . We observe a clear peak at 200 Hz reaching slightly more than 0.9 µm indicating a kind of resonance phenomenon, i.e., where the driver frequency of the applied force matches that of the natural frequency response of the material dictated by the constitutive relationships (13) and the corresponding silicone properties used in Tab. 1. Note that at f AM = 0 the stimulus is essentially the same as the UM scenario and thus the maximum displacement matches that of the simulated steady state in Fig. 4.

C. SPATIO-TEMPORAL MODULATION SCENARIO
The STM stimulus resembles a moving load scenario where the force applied changes the point at which it is applied over time. Here, the pressure is constant at A(t) = 1 and the focus location follows a closed circular loop described by z(t) = (r 0 cos(2π f STM t), r 0 sin(2π f STM t), 20) cm with a circle radius is r 0 cm resulting in a non-localised vibrotactile stimulus with a rotation frequency of f STM Hz. The surface speed of the focal point can be calculated through u = 2πr 0 f STM cm/s.  As previously described in the Sec. V-C, a stationary probe is simulated and placed along the trajectory of the moving focus which we sample in time thereby observing the oscillatory vibrations observed in Fig. 9a). Every dip in the time series represents the applied pressure of a moving focal point as it passes over the probe inducing a surface indentation followed by a relaxation back to its equilibrium. The time-width of each dip can be approximately calculated if we ignore any surface wave effects by setting it equal VOLUME 4, 2016 to the width of the focal point radiation pressure divided by its speed u. In this instance, the circular trajectory repeats at 70 Hz, therefore the focal point is moving at a speed of u = 4.4 m/s. If the width of the focal point pressure is approximately given by the full width at tenth of maximum (FWTM) of the Gaussian fit in equation (4), then the radiation pressure width is less by a factor of √ 2 due to the squaring of the pressure in equation (5). Using that FWTM = ln 10 ln 2 FWHM we obtain that the width of the focal point radiation pressure is 18.5 mm, and the time-width of the dips in Fig. 9a) are approximately 4.2 ms. It therefore follows that the periodic time series observed in Fig. 9a) can be represented by a Fourier series x(t) = ∑ n c n e jn f STM t . Taking the Fourier transform of x(t) we obtain X( f ) = 2π ∑ n c n δ ( f − n f STM ) demonstrating why we observe harmonics at integer multiples of f STM of decreasing amplitude in Fig. 9b).

VII. PARAMETER AND TACTILE PATTERN INVESTIGATIONS
Having developed and then validated our SLSM viscoelastic simulation model, in this section we investigate the behaviour of different viscoelastic materials during ultrasound mid-air haptic stimulation. Namely, we investigate the effect of different viscoelastic parameters (viscosity η 3 and elastic modulus G d ) and different STM tactile patterns (square, line, and circle). The first investigation can provide insights on how different skin types and age will respond to mid-air haptic stimuli, while the second can help us understand how the skin responds to different stimuli.

A. VISCOELASTIC PROPERTIES INVESTIGATION
It is well known that the viscoelastic properties differ greatly between individuals, change over time, and that they also differ on average between men and women. Namely, depending on the measurement, anatomical site, skin hydration level, age, individual person, and theoretical model applied, elastic moduli of human skin have been reported to vary between 4.4 kPa to 57 MPa [38]. While our simple SLSM model is built upon silicone properties frequently used as a hand skin mimicking phantom in medical engineering research [35] (see also Tab. 1), in this subsection we shall try to extrapolate our model and make predictions on how efficiently acoustic energy generates vibrations for different skin characteristics. Specifically, ageing generally has the effect of decreasing the skin thickness, elasticity (i.e., increase in elastic modulus G d ), and viscosity η 3 . Meanwhile, skin hydration can generally reduces the elasticity of human skin by an order of magnitude. Fig. 10a) presents the simulated surface displacement with an AM stimuli at f AM = 200 Hz. The results report a displacement heatmap of different viscoelastic properties. On the y−axis we vary the elastic modulus G d with values ranging from 3 kPa to 300 kPa. On the x−axis we vary the viscosity η 3 taking values ranging from 0.342 Pa · s to 34.2 Pa · s. This large parameter range search of viscoelastic properties permits to study the displacement induced by mid-air stimuli for very different materials (very soft to medium). A zoomed in version is also displayed in Fig. 10b) for parameter values closer to the ones used to calibrate our model (see Tab. 1).
We observe changes in peak-to-peak displacement ranging from 0.1 µm to 5 µm. It is also observed that the generated displacement is more affected by the elasticity parameter G d than viscosity η 3 . In particular, we observe that a higher elastic modulus skin characteristic of older and dryer skin will result in less displacement. Despite this, however, we note that the minimum simulated displacement within our exploration range is still higher than the reported human detection threshold [39].

B. DIFFERENT TACTILE PATTERN INVESTIGATION
In Sec. VI we validated our SLSM silicone model using three specific types of mid-air haptic stimuli: UM, AM, and STM. These were in some sense analogous to standard rheological tests, such as a creep test, a dynamic indentation test, and a moving load scenario. In this section, we study the spatial and temporal vibrations produced on our viscoelastic material due to different STM stimuli.
We study three different shapes that are commonly found and used in mid-air haptic literature and demonstrators: circle, line, and square, defined and implemented into COMSOL using the notation introduced in equation (14) with A(t) = 1 in all three cases as shown in Fig. 11. The circle STM stimulus is defined as in the previous sections using z(t) = (r 0 cos(2π f STM t), r 0 sin(2π f STM t), 20) cm while using r 0 = 1 cm and f STM = 70 Hz. The line STM stimulus is a modification of the circle stimulus in that we limit the focal point's motion along the x direction on the surface such that z(t) = (r 0 cos(2π f STM t), 0, 20) cm while using r 0 = 3 cm and f STM = 150 Hz. Finally, the square STM is also a modification of the circle stimulus in that we use the Fernández-Guasti squircle formulation [40] that allows for a κ ∈ [0, 1] parametrization that smoothly interpolates between a circle and square such that z(t) = (ρ(t) cos(2π f STM t), ρ(t) sin(2π f STM t), 20) cm using where κ = 0 produces a circle of radius r 0 , and κ = 1 a square of side 2r 0 . The square implemented in Fig. 11 is using r 0 = 1 cm, f STM = 150 Hz, and κ = 0.95 to avoid sharp corners. Fig. 11 shows the resulting time-averaged peak to peak (P2P) displacement of the viscoelastic surface along and near the stimulus path in the form of a heatmap (top row), the displacement time-series at different probe locations along the stimulus path (middle row), and the corresponding Fourier frequency spectrum (bottom row). The methods of collecting and processing the data displayed in Fig. 11 were described in Sec. V-C. The surface displacements generated are summarised in Tab. 3. The results from the circle STM stimulus are very much what we expected and in line with the LDV validation study presented earlier in this paper. Namely, the displacement heatmap is that of a circle with a P2P displacements of about 1.23 µm, the displacement time series consists of a sequence of spikes that repeat every 1 f STM = 14.28 ms, and the frequency spectrum displays a main peak at the fundamental frequency of f STM = 70 Hz followed by harmonics at integer multiples. It is worth noting that there is a kind of doublehump rebound effect that follows each spike indentation. Presumably this is a wave train effect whereby the main indentation caused by the moving focus leaves behind it a trail of vibration waves that self interfere constructively or destructively.
The results from the line STM stimulus are more interesting. Namely, the displacement heatmap resembles a dashed line with clearly noticeable gaps. The P2P displacement at the edges and at the mid-point of the line are about 1.42 µm but only 0.85 µm at the two gaps seen at x = ±18 mm. Also, the displacement time series and corresponding frequency spectrum displays very different vibrational signatures depending on where the probe is positioned. Specifically, while in both probe locations the driving frequency of the stimulus is f STM = 150 Hz, this frequency is lost (not observed) at the mid-line probe because of a doubling effect when the focus passes over the mid-point of the line twice per cycle. Therefore, the spectrum at the mid-point of the line is composed of a fundamental frequency of 2 f STM and harmonics at integer multiples of that, i.e., at 2n f STM and n ∈ N + . In contrast, f STM is clearly visible at the edge probe however is not the dominant one. The reason for this can be extracted from the time series where we can clearly VOLUME 4, 2016 see oscillations taking place at a frequency of 2 f STM even though the peak displacement spike caused by the moving focus repeats at f STM .
A closer inspection and comparison of the line stimulus time series and that of the circle suggests that the doublehump rebound effect that follows each peak displacement spike is intensified here, presumably through constructive interference of different wave trains travelling in opposite directions caused by the focus that decelerates as it approaches the edge, stops, and then accelerates in the reverse direction. Recall that the linear speed of the focal point is given by u(t) = |ż(t)| = 2π f STM | sin(2π f STM t)|, where we use the dot notationż to indicate a derivative with respect to time dz dt . A similar phenomenon could explain the reduced P2P displacement at x = ±18 mm, namely that the applied pressure by the moving focus is opposed by a rebounding wave train propagating in the opposite direction.
While further analysis is needed to fully understand these surface wave non-linear interference effects, one thing which is very interesting is that the frequency response of the skin, and therefore the activation of the embedded mechanoreceptors, varies greatly along the line stimulus. One could therefore argue that the line stimulus could potentially be approximated by three AM focal points sparsely positioned along a straight line with the middle focus driven at f AM = 2 f STM and the two outer ones at f AM = f STM . It is worth mentioning here that the line stimulus is very similar to the Lateral Modulation (LM) studied by Takahashi et al. [41].
The square with smoothed corners STM stimulus has a more predictable behaviour, similar to the circle. Namely, the displacement heatmap is that of a square with smoothed corners with a P2P displacements of about 1.19 µm, the displacement time series consists of a sequence of spikes that repeat every 1 f STM = 6.6 ms, and the frequency spectrum displays a main peak at the fundamental frequency of f STM = 150 Hz followed by harmonics at integer multiples. The probe at the corner and that at the edge of the pattern are similar, but not completely identical. The corner one has a more pronounced double-hump rebound while the edge does not. Perhaps this could be explained by looking at the non-constant speed of the focal point that travels faster round the corners than on the straights. Recall that the linear speed of the focal point round the squircle can be calculated by u(t) = |ż(t)| = 4 f 2 STM π 2 ρ(t) 2 +ρ(t) 2 . The insights provided by the above discussions and Fig.  11 suggest that different tactile patterns and the way they are implemented (rendering algorithm) can generate very rich and very different vibrations on the skin which themselves then are likely to produce very different haptic sensations, perceptions, and experiences.

VIII. CONCLUSION AND DISCUSSION
FEM is a powerful engineering tool that can model a complex material and its response to different forces. In this paper, we have applied this tool to the field of ultrasound mid-air haptics [1] in order to better understand and make predictions about how skin-mimicking viscoelastic materials react to different mid-air haptic stimuli produced by high intensity focused acoustic forces. Namely, we were able to study for the first time ever how varying acoustic pressures of the order of 160 dB SPL can produce localised indentations and dynamic surface vibrations of the order of a micrometer, and how these then induce damped surface waves that propagate and interfere constructively or destructively thereby modifying the original tactile input stimulus. Never before could we get such accurate and high-resolution data about these interactions at an acoustic-viscoelastic interface. Previous investigations were limited to costly and timeconsuming in-silico and in-vivo LDV experiments [11]- [13]. Our FEM approach can therefore accelerate such explorations while providing high resolution rheological data in a highly flexible computer environment. Importantly, we believe that the engineering and mechanical insights unveiled about what is happening on the skin during mid-air haptic stimulation can help one design and optimize better stimuli that efficiently activate our skin receptors to produce rich somatosensory responses and touch perceptions.
First, we mathematically modelled and experimentally validating the acoustic forces applied by a focusing phased array of transducers (see equation (6)). These include forces from both the acoustic radiation pressure and acoustic streaming. For the first time, we were able to separate their respective contributions towards the total applied force (see Fig. 2). We then constructed a dynamic SLSM rheological model of a viscoelastic skin-mimicking material with a frequency dependent elastic modulus (see (12)) and derive its constitutive equation (see (13)). The two models were then ported into a COMSOL environment which we used to simulate different interface vibrational responses and validated them against LDV data reported in [12]. Our simulation results were in very good quantitative and qualitative agreement with the LDV data used as a benchmark (see Tab. 2). Moreover, the physical understanding of each scenario studied was compared to well-understood physical rheological tests providing us with a deeper understanding about the different phases of the viscoelastic material's strain response (see Fig. 5) and the efficiency of different stimuli in inducing a resonant or high displacement response (see Fig. 8). In the second half of the paper we performed a number of computer explorations and analysed the response of the simulated Ecoflex phantom material to different midair haptic stimuli used for rendering shapes like circles, lines and squares (see Tab. 3). We also studied how the response depends on the viscoelastic properties (viscosity and elastic modulus) of the material itself, thus emulating age and hydration effects and how these may affect the induced surface oscillations and displacements (see Fig.  10). Importantly, we were able to obtain spatial, temporal, and frequency-rich data showing how different tactile input patterns (e.g., circle, line, and square STM curves) produce different vibrational signatures at different parts of the targeted surface (see Fig. 11). Notably, we observed strong surface wave interactions that can produce wave-trains of vibrations, a doubling of the frequency spectrum, and a skinrebound (outdent) of up to 0.45 µm.
Human skin is a very complex material and varies greatly within the population. Despite the useful insights that we have obtained above, the models and simulations presented in this paper have significantly simplified the problem and are therefore in need of further experimental validations and model enhancements. Namely, we have not modelled the geometry and structure of the skin, and particularly that of the hand which is where mid-air haptics are usually applied to. Fingers, palms, nails and skin creases will surely affect how surface wave vibrations propagate, dampen and interfere. Also, the various layers of the skin (epidermis, dermis, etc.) will display different viscoelastic properties. Therefore, future studies should aim to upgrade our model to include a more complex hand geometry composed of a multi-layer rheological model. Finally, it would be very interesting to connect and integrate our mechanical model to neurocognitive ones that model how tactile afferents respond to skin indentations [42] and how these then translate into tactile sensations [43] and experiences [44].