text:
# Non-linear effects and shock formation in the focusing of a spherical acoustic wave

## 1 Introduction

Recent experiments have shown that acoustic waves can be used to study the nucleation of phase transitions far from equilibrium under very clean conditions lambare98 ; caupin\_b01 ; chavanne\_b\_c01 . Thanks to hemispherical piezo-electric transducers, we have focused 1MHz acoustic waves in liquid helium and produced large pressure and density oscillations. These waves are quasi-spherical and, at the acoustic focus (the center), their amplitude can be very large. We used an optical method to detect the nucleation of bubbles by the negative swings of the waves lambare98 ; caupin\_b01 . This nucleation occurs beyond a certain threshold in the sound amplitude which needs to be determined as accurately as possible, in order to compare with independent theoretical predictions. We later obtained evidence for the nucleation of crystals by the positive swings chavanne\_b\_c01  and had the same need.

In the absence of non-linear effects, the measurement of the nucleation threshold would be simple to do. For example, one could study the nucleation as a function of the static pressure in the experimental cell, and then use a linear extrapolation caupin\_b01 . However, we expect non-linear effects to occur, especially in cavitation studies. Indeed, the homogeneous nucleation of bubbles occurs near the “spinodal limit” where the compressibility diverges and the sound velocity vanishes. When an acoustic wave is produced in a fluid, with an amplitude such that during the negative swings the sound velocity approaches zero, it is clear that the wave must be highly distorted. Non-linear effects have been already noticed by several authors sirotyuk ; roy\_m\_a90 ; nemirovskii90 .

We thus try to calculate the non-linear focusing of the acoustic waves. We start with the spherical geometry, because in a first approximation, everything depends only on the radial distance $`r`$ from the center. As we shall see (section 2), this calculation still appears difficult because the focusing of acoustic waves leads to the formation of shocks at all amplitudes in a spherical geometry, and their treatment is not quite straightforward. We first obtain this result and the associated shock length from an analytic calculation which uses the methods of characteristics (section 2.2). Our calculation extends the former work of Nemirovskii nemirovskii90  to the spherical case, except that we neglect the coupling with heat modes. It is done in the spirit of Greenspan and Nadim’s work greenspan\_n93 , though in our case it is slightly more tricky due to the shape of the equation of state. We make it quantitative by using the equation of state of liquid helium maris91  which is well established. For weak oscillation amplitudes of the transducer, shocks can be ignored and the pressure calculated at the focal point by simulating the Euler equations using a finite difference method (section 2.3). Indeed, shocks form with an infinitesimal amplitude at a distance r from the center which is much less than our mesh size, so that one can neglect them. At larger amplitude, shocks cannot be ignored. In order to treat the shocks, we have adapted to the case of helium a code devoted to shock simulation, based on a WENO scheme (section 2.4). In the end, we obtain the amplitude of the density oscillation at the focus as a function of an important parameter, the amplitude of the displacement of the transducer surface where waves are generated. This transducer is a piezo-electric ceramic. We choose 1MHz for the frequency of the waves, in order to compare with the experiments. In parallel, we have built an experiment to measure the focusing in a quasi-spherical geometry. As explained in section 3, the results of this experiment allow a precise comparison with our theoretical and numerical work chavanne02 ; appert01a . We find that the shape of the acoustic wave is indeed distorted at high amplitude and very well described by our calculations, thus validating our theoretical method. The final comparison with our calculation allows us to calibrate the efficiency of the ceramics. As described in our conclusion, this work should now be extended to different geometries. One of them is the hemispherical geometry where, according to other experimental results caupin\_b01 , non-linear effects are apparently less important, an observation which needs to be understood and compared with future calculations.

## 2 Theory

### 2.1 Description of the model

Throughout this paper, we consider a spherical geometry. We take it as one-dimensional since the pressure and density fields only depend on the radial distance $`r`$. We neglected dissipation since our main goal was to compare with experiments in superfluid helium 4 which has zero viscosity and where the attenuation of sound vanishes in the low temperature limit.

In the case of liquid helium 4 at zero temperature, the equation of state has been obtained by three different methods (sound velocity extrapolations, density functional calculations, and Monte Carlo simulations) with similar results. Maris maris91  uses the simple form to relate the pressure ($`P`$) to the density ($`\rho `$):

$$PP_{\mathrm{sp}}=\frac{b^2}{27}(\rho \rho _{\mathrm{sp}})^3$$
(1)
with

$`P_{\mathrm{sp}}`$ $`=`$ $`9.6435\mathrm{bar}`$
$`\rho _{\mathrm{sp}}`$ $`=`$ $`94.18\mathrm{kg}\mathrm{m}^3`$
$`b`$ $`=`$ $`14.030\mathrm{m}^4\mathrm{s}^1\mathrm{kg}^1.`$
$`P_{\mathrm{sp}}`$ is the spinodal limit where the compressibility diverges, the sound velocity vanishes and the liquid becomes totally unstable against the formation of the vapor. At $`P=0`$, the density is $`\rho _0=145.13\mathrm{kg}\mathrm{m}^3`$. The sound speed is then $`c_{\mathrm{s0}}=238.3\mathrm{m}\mathrm{s}^1`$ and the wavelength $`\lambda _0=c_{\mathrm{s0}}T=0.238\mathrm{mm}`$ for 1 MHz waves whose period is $`T=1/f=1\mu \mathrm{s}`$.

We also considered helium 3, a lighter liquid which is not superfluid except at very low temperature - i.e. below the achievable temperature in our experiment. The same form is used for the equation of state, now with

$`P_{\mathrm{sp}}`$ $`=`$ $`3.1534\mathrm{bar}`$
$`\rho _{\mathrm{sp}}`$ $`=`$ $`53.50\mathrm{kg}\mathrm{m}^3`$
$`b`$ $`=`$ $`19.262\mathrm{m}^4\mathrm{s}^1\mathrm{kg}^1.`$
The value of $`P_{\mathrm{sp}}`$ is less negative, which means that the inner cohesion of liquid helium 3 is weaker than for helium 4. At $`P=0`$, the density of liquid helium 3 is $`\rho _0=81.916\mathrm{kg}\mathrm{m}^3`$. The sound speed is then $`c_{\mathrm{s0}}=182.5\mathrm{m}\mathrm{s}^1`$ and the wavelength $`\lambda _0=c_{\mathrm{s0}}T=0.182\mathrm{mm}`$.

At the temperatures considered, the viscosity is very weak and can be neglected though, in case of helium 3, neglecting dissipation is an approximation which would need to be better justified. Therefore, the numerical approximation considers the Euler equations. In order to use a dimensionless form of the equations, we have chosen as a time scale the period of the wave $`T`$, as a length scale the wavelength at zero pressure $`\lambda _0=c_{\mathrm{s0}}T`$ and as a density scale the spinodal density $`\rho _{\mathrm{sp}}`$. If we now consider $`\rho `$ as the dimensionless density and $`u`$ as the dimensionless velocity of the fluid, the Euler equations are written as follows:

$`_\mathrm{t}\rho +u_r\rho +\rho _ru`$ $`=`$ $`{\displaystyle \frac{\theta \rho u}{r}}`$
$`_\mathrm{t}u+u_ru`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}_rP,`$ (2)
or, by using the conservative variables $`\rho `$ and $`j=\rho u`$,

$`_\mathrm{t}\rho +_rj`$ $`=`$ $`{\displaystyle \frac{\theta j}{r}}`$
$`_\mathrm{t}j+_r[j^2/\rho +P(\rho )]`$ $`=`$ $`{\displaystyle \frac{\theta j^2}{\rho r}}`$ (3)
where $`\theta `$ is respectively $`0`$, $`1`$, or $`2`$ in planar, axisymetric cylindrical, or spherical geometry. The equation of state, governing the dimensionless pressure variations, reads:

$$P\frac{P_{\mathrm{sp}}}{\rho _{\mathrm{sp}}c_{\mathrm{s0}}^2}=\frac{C_\mathrm{o}^2}{3}(\rho 1)^3$$
(4)
with $`C_\mathrm{o}=1.848`$ for helium 4 and $`C_\mathrm{o}=1.883`$ for helium 3. Note that, in order to have simple notations, we use the same names for physical and reduced variables. In general, calculations will be performed with reduced variables, while numerical results and experimental parameters will be given as physical quantities.

Boundary conditions are imposed at the center ($`r=0`$)

$$u=0,$$
(5)
and on the transducer surface ($`r=L_0`$)

$$u(t)=\omega \mathrm{\Delta }x_0\mathrm{sin}(\omega t)\left[1\mathrm{exp}(t/1.5)\right].$$
(6)
Indeed, the motion of the transducer surface is much smaller than a mesh size in all our simulations, and it is sufficient to impose an oscillating velocity on the ceramics. The exponential term represents the response time of the transducer. In the experiment the response time is rather equal to $`8`$ $`\mu `$s. Here we took it shorter to have a more rapid convergence of the calculation to the steady regime. This has no effect on the final result, as our simulations will always be used in the stationary regime.

In the whole paper, all the calculations are based on Euler equations, for which shock waves do occur. We are aware that, in superfluids as helium 4, there is no shock wave strictly speaking. Actually, when a steep gradient appears, it is regularized by dispersion instead of dissipation. This means that modes propagate with a different velocity depending on their frequency, and thus highest frequency modes are expelled from the steep region. However, we assume that, as shocks form only in the focal region, and during a limited time, not too much momentum is lost locally due to dispersion, and that in a first approximation, a viscous or dispersive regularization of the shocks is equivalent (note that any numerical scheme able to handle shocks will always introduce either a dissipation or a dispersion term in order to make shocks regular, even if Euler equations do not contain any viscous term). Besides, as we shall see (section 2.3.2), we do not need to describe the shock structure exactly, as we are above all interested in the relaxation part of the wave. We have found from our simulations that the negative pressure swing only depends very weakly on the numerical viscous regularization we use - and thus on the shock amplitude near the focal point. Still, it would of course be interesting to be able to quantify further the effects of dispersion, and this could be the object of a further work.

Throughout the paper, we shall consider helium 4 unless it is explicitly specified that it is helium 3. In the next sections, we show that shock waves can occur in this system, and we compute the radius (denoted as shock length in the followings) at which a shock wave occurs for various oscillation amplitudes $`\mathrm{\Delta }x_0`$ of the transducer surface.

### 2.2 The method of characteristics

In the study of compressible fluids, the method of characteristics is a standard one characteristics . It has already been used for helium in a planar geometry (see for example nemirovskii90 ).

Here we are interested in the shock length in spherical geometry, for the equation of state (1). In this paper we define the shock length as the distance from the center where the shock forms. We shall use the method of characteristics to predict a lower bound for the shock length.

#### 2.2.1 Rieman invariants and characteristics

Let us first recall the principle of the method. Solving the Euler equations means knowing the density $`\rho `$ and the velocity $`u`$ everywhere within a certain domain of the (r,t) plane. In our case, the “characteristics” are two families of curves $`\left\{𝒞_i^{(+)}\right\}_{iR}`$ and $`\left\{𝒞_i^{()}\right\}_{iR}`$ parameterized by $`i`$, where the parameter $`i`$ could for example be defined as the $`r`$-coordinate at which the characteristic cuts the axis $`t=0`$ (see Fig. 1 for an illustration). Thus, $`𝒞_i^{(+)}`$ (resp. $`𝒞_i^{()}`$) refers to a single curve in the (r,t) plane, which belongs to the $`(+)`$ family (resp. $`()`$), and has in our case a positive slope (resp. negative). Each family completely covers the (r,t) plane as $`i`$ is varied. The two families of curves intersect each other. Then, instead of locating a point in the plane by its coordinates (r,t), it is equivalent to indicate which particular characteristics $`𝒞_i^{(+)}`$ and $`𝒞_j^{()}`$ intersect each other at this point. The point (r,t) will then equivalently be referred to as point $`(i,j)`$.

These families are chosen so that, for each family $`\left\{𝒞_i^{(k)}\right\}`$ with $`k=\pm `$, there exists a quantity $`I_k`$ called “Rieman invariant”, which is a function of $`\rho `$ and $`u`$, and obeys a simple evolution equation along any characteristic $`𝒞_i^{(k)}`$ of the family. If the value of $`I_k`$ is known at one point of a characteristic $`𝒞_i^{(k)}`$ (for example at the initial time), then it is easy to compute it on the whole curve.

As each point of the (r,t) plane is the intersection of two characteristics $`𝒞_i^{(+)}`$ and $`𝒞_j^{()}`$, we know the values of $`I_+(\rho ,u)`$ and $`I_{}(\rho ,u)`$ at this point. Then the density $`\rho `$ and the velocity $`u`$ are entirely determined.

Now the precise form of the Rieman invariants and corresponding characteristics have to be derived from the Euler equations (2) and the equation of state for helium (4). The detailed calculations are given in Appendix A. Here we just give the results concerning the shape of the characteristics and the Rieman invariants. By definition, at any point, the slope of a characteristic $`𝒞_i^{(+)}`$ with equation $`r=x(t)`$ is $`dx/dt=u+c_\mathrm{s}`$ where $`u`$ and $`c_\mathrm{s}`$ are taken in $`(x(t),t)`$. For a characteristic $`𝒞_i^{()}`$, it is $`dx/dt=uc_\mathrm{s}`$. The sound speed $`c_\mathrm{s}`$ is given from the equation of state (4) by

$$c_\mathrm{s}=C_\mathrm{o}(\rho 1).$$
(7)
Then the derivative along the characteristic reads

$$\frac{d}{dt}=_\mathrm{t}+(u\pm c_\mathrm{s})_r.$$

The Rieman invariants are found to be

$$I_\pm C_\mathrm{o}(\rho \mathrm{ln}\rho )\pm u$$
(8)
The equations verified by the “invariants” $`I_k`$ with $`k=\pm `$ read

$$\frac{d}{dt}\left[I_\pm \right]+\theta \frac{C_\mathrm{o}(\rho 1)u}{r}=0$$
(9)
where the derivative $`\frac{d}{dt}`$ is taken along a characteristic $`𝒞_i^{(k)}`$. Again, $`\theta `$ is respectively $`0`$, $`1`$, or $`2`$ in planar, cylindrical, or spherical geometry.

In the case of planar geometry ($`\theta =0`$), $`I_k`$ is a true invariant, since it is constant along the corresponding characteristic, hence its name. In spherical or cylindrical geometries, it is not constant, due to the source term in Eq.(9), though it is still called an “invariant”.

All this is valid at least as long as the characteristics belonging to the same family do not intersect each other. When they do, the corresponding $`I_k`$ becomes multivalued. This is an indication that a shock has formed, and beyond the corresponding time, the description used in this section breaks down.

#### 2.2.2 Lower bound for the shock length - Analytic calculation

We shall now show that such an intersection does occur in the system, and calculate the corresponding shock length. We consider a spherical domain bounded by a spherical piston. At $`t=0`$, the fluid is at rest with a density $`\rho _{\mathrm{st}}`$ (label “st” standing for static) and the piston surface is a sphere of radius $`L_0`$. As long as the fluid is at rest, all the characteristics are straight lines (see Fig. 2). The respective slopes of the characteristics $`\left\{𝒞_i^{(+)}\right\}_{iR}`$ and $`\left\{𝒞_i^{()}\right\}_{iR}`$ are $`+c_{\mathrm{st}}`$ and $`c_{\mathrm{st}}`$, $`c_{\mathrm{st}}`$ being the sound velocity for the initial density $`\rho _{\mathrm{st}}`$. The piston starts to move at $`t=0`$ with a velocity $`v_p(t)=\mathrm{\Delta }v_0\mathrm{sin}(\omega t)`$.

We denote $`𝒞_0^{()}`$ the characteristic which originates from $`r=L_0`$ when t=0, with slope $`c_{\mathrm{st}}`$. The domain to the left of $`𝒞_0^{()}`$ is unperturbed unless some characteristic $`𝒞_i^{()}`$ crosses $`𝒞_0^{()}`$, i.e. unless there is a shock.

The aim of this calculation is to find an upper bound for the time necessary to form a first shock in the system. As the piston moves, it emits some characteristics which will cut $`𝒞_0^{()}`$ after a while, leading to a shock. We only study how the characteristics emitted at early times, and almost parallel to $`𝒞_0^{()}`$, will eventually cross it. Of course, some other characteristics emitted later could cross $`𝒞_0^{()}`$ earlier, or some shocks could occur somewhere else at earlier times. That is why our calculation only gives an upper bound for the shock time.

The details of the calculation are given in Appendix B. We find that an upper bound for the time of shock formation is

$$t_{\mathrm{shock}}\frac{L_0}{c_{\mathrm{st}}}\left\{1\mathrm{exp}\left[\frac{c_{\mathrm{st}}^2}{2L_0\omega \mathrm{\Delta }v_0}\frac{\rho _{\mathrm{st}}1}{\rho _{\mathrm{st}}\frac{1}{2}}\right]\right\}\frac{L_0}{c_{\mathrm{st}}}.$$
(10)
As the corresponding shock length $`r_{\mathrm{shock}}`$ is measured from the center of the sphere, a lower bound for $`r_{\mathrm{shock}}`$ is

$$L_0r_{\mathrm{shock}}L_0\mathrm{exp}\left[\frac{c_{\mathrm{st}}^2}{2L_0\omega \mathrm{\Delta }v_0}\frac{\rho _{\mathrm{st}}1}{\rho _{\mathrm{st}}\frac{1}{2}}\right]>0.$$
(11)
From this study we conclude that in spherical geometry, there is always formation of a shock, whatever the velocity of the piston is, as long as it has a nonzero acceleration towards the center of the sphere. However, when the oscillation amplitude $`\mathrm{\Delta }x_0=\mathrm{\Delta }v_0/\omega `$ goes to zero, the shock forms very near the focal point. As shocks are formed by the intersection of tangential characteristics, their initial amplitude is zero. If they are formed very near the center, their amplitude does not have time to grow much. This is true especially because of the existence of a cut-off: the notion of shock becomes meaningless when the width of the shock becomes of the same size as the shock length itself. So, as the oscillation amplitude tends to zero, the jump in density and velocity at the shock also vanishes, and their is no contradiction with the fact that the solution is expected to approach the linear solution landau\_l6 .

The above calculation is analytic and gives only a lower bound for the shock length because we have assumed that the shock occurs on the first characteristic $`𝒞_0^{()}`$.

#### 2.2.3 Numerical calculation of the shock length

As we shall see now, a numerical simulation shows that our calculation of the lower bound is in fact a good estimate of the shock length itself. Let us now solve Euler’s equations numerically by computing a network of characteristics $`\left\{𝒞_i^{()}\right\}_{iZ}`$ and $`\left\{𝒞_i^{(+)}\right\}_{iZ}`$. The values of the Rieman invariants, and thus of the density $`\rho `$ and velocity $`u`$, will be defined on the intersections of the left $`𝒞_j^{()}`$ and right $`𝒞_i^{(+)}`$ characteristics. As the fluid is at rest on the left of $`𝒞_0^{()}`$, we restrict our calculation to the $`(r,t)`$ domain located between $`𝒞_0^{()}`$ and the piston trajectory.

Several parameterization choices are possible. For the $`(+)`$ family, we have chosen to take $`i=t_i/\delta t`$ for each characteristic $`𝒞_i^{(+)}`$, where $`t_i`$ is the time at which $`𝒞_i^{(+)}`$ and $`𝒞_0^{()}`$ intersect, and $`\delta t`$ is an arbitrary fixed time step. We discretize the system by restricting $`i`$ to integer values. This means that the subset of characteristics that will be considered for the simulation are initialized at regular time intervals $`i\delta t`$ when they cross the first characteristic $`𝒞_0^{()}`$ emitted by the piston (see Fig. 2).

When $`𝒞_i^{(+)}`$ meets the piston, a new characteristic $`𝒞_i^{()}`$ (with the same label $`i`$) is emitted from the piston at the same time. This defines the parameterization of the $`()`$ family. The intersections of these two families of characteristics form a network of points whose locations will be determined in the following. The intersection of the characteristics $`𝒞_i^{(+)}`$ and $`𝒞_j^{()}`$ is referred to as point $`(i,j)`$.

During the $`n^{th}`$ step, we compute all the points $`(i,j)`$ such that $`i+j=n`$ with $`i`$ and $`j`$ integers. Let us describe now how a point $`(i,j)`$ can be computed from the points of the previous step (see Fig. 3). All the information comes from the two sites $`(i,j1)`$ and $`(i1,j)`$. We must extrapolate each characteristic $`𝒞_i^{(+)}`$ and $`𝒞_j^{()}`$ up to the next intersection $`(i,j)`$. We compute the local slopes $`u\pm c_\mathrm{s}`$ of the characteristics $`𝒞_i^{(+)}`$ and $`𝒞_j^{()}`$ in sites $`(i,j1)`$ and $`(i1,j)`$ respectively. Then $`(i,j)`$ is the intersection of the two straight lines which respectively go through $`(i,j1)`$ and $`(i1,j)`$ with these slopes.

The values of the Rieman invariants $`I_+`$ and $`I_{}`$ at $`(i,j)`$ are found by numerical integration of the equations (8,9) along the two involved characteristics. From these values, both $`\rho `$ and $`j=\rho u`$ can be obtained.

The program stops whenever two characteristics of the same family cross each other, as shown on Fig. 2. Then a shock occurs and the calculation based on characteristics breaks down (there would be multivalued points).

Let us now compare analytic and numerical results. Both calculations were done by taking $`L_0`$ equal to the experimental cell radius, i.e. 8 mm = 33.6 $`\lambda `$. In Fig. 4, we plot the shock distance $`r_{\mathrm{shock}}`$ measured from the center of the sphere as a function of the amplitude of the oscillation of the piston. We compare it with the analytic result. The agreement is excellent. This shows that our analytic calculation not only gives a lower bound but in fact a good estimate for the shock distance itself.

We have also performed the analytic calculation in the case of helium 3, neglecting the role of viscosity. It is expected that shocks form at smaller wave amplitude in helium 3, because, in this lighter liquid, the spinodal pressure is less negative than in helium 4.

### 2.3 Numerical simulations of Euler’s equation

#### 2.3.1 Numerical method

As we shall see now, a simple finite difference numerical scheme is sufficient to simulate our system with moderate amplitude, at least as long as the shock length is smaller than the spatial discretisation step. Our aim is to calculate numerically the pressure and density oscillation at the center. We chose to have two lattices, one for mass and the other for momentum. They are staggered (Fig. 5) and allow us to enforce exactly the conservation of mass:

$$\rho _k^{t+\delta t/2}=\rho _k^{t\delta t/2}+\frac{\delta t}{\delta r}\frac{\left[\left(k\frac{1}{2}\right)^2j_{k1}^t\left(k+\frac{1}{2}\right)^2j_k^t\right]}{k^2+\frac{1}{12}}.$$
(12)
We have taken the momentum equation in the form:

$`j_k^{t+\delta t}`$ $`=`$ $`j_k^t{\displaystyle \frac{\delta t}{\delta r}}\left[{\displaystyle \frac{(j_k^t+j_{k+1}^t)^2}{4\rho _{k+1}^{t+\delta t/2}}}{\displaystyle \frac{(j_k^t+j_{k1}^t)^2}{4\rho _k^{t+\delta t/2}}}\right]`$ (13)
$``$ $`C_0^2\left({\displaystyle \frac{\rho _k^{t+\delta t/2}+\rho _{k+1}^{t+\delta t/2}}{2}}1\right)^2{\displaystyle \frac{\rho _{k+1}^{t+\delta t/2}\rho _k^{t+\delta t/2}}{\delta r}}\delta t`$
$``$ $`4{\displaystyle \frac{\delta t}{\delta r}}{\displaystyle \frac{(j_k^t)^2}{\left(k+\frac{1}{2}\right)\left[\rho _k^{t+\delta t/2}+\rho _{k+1}^{t+\delta t/2}\right]}}.`$
As we use a staggered lattice, we only have to specify the boundary conditions for the momentum. It is vanishing at the center of the sphere (r=0), and thus the symmetry with respect to the center imposes $`j_0=j_1`$ (see Fig. 5). On the other hand, the motion of the piston is implemented by

$$j_K(t)=\rho _K\omega \mathrm{\Delta }x_0\mathrm{sin}(\omega t)\left[1\mathrm{exp}(t/1.5)\right].$$

#### 2.3.2 Focal pressure

In Fig. 6, the focal pressure is represented as a function of time. It is calculated from the average density in the central cell. The results of Fig. 6 were obtained for a cell length equal to the experimental one, i.e. $`L_0=33.6\lambda `$, an oscillation amplitude $`\mathrm{\Delta }x_0=0.7`$ nm, and a zero static pressure. Then the reduced sound velocity is equal to 1, and the wave needs 33.6 time units to reach the center of the cell. We are interested in the steady regime, which is established around $`t/T>45`$. For $`t/T>33.6`$, the wave reflected on the focal point propagates towards the ceramic, and then back to the focal point. It reaches the latter at $`t/T=100.8=3\times 33.6`$. Then we stop our measurements, to be consistent with the experiment, which uses short bursts.

As long as the oscillation amplitude of the ceramic is not too large, non-linear effects are negligible far from the focal point. If we observe the focal pressure during the steady regime, we find it also nearly sinusoidal, and the positive swings are only slightly larger than the negative ones (Fig. 7).

We have checked that the amplitude of the density is well represented by the function $`\mathrm{sin}(kr)/(kr)`$, as predicted from the linear theory landau\_l6  (see Fig. 8).

From our calculation of section 2.2, we expect a shock to occur near the focal point. However, for the case of Fig. 8, it happens on a region around $`r=0`$ that cannot be seen because it is much smaller than the mesh size. Indeed, the predicted reduced shock length would be $`r_c/\lambda =4.2\times 10^{28}`$, to be compared with the mesh size $`\delta r/\lambda =0.01`$.

When $`\mathrm{\Delta }x_0`$ is increased to $`5`$ nm, non-linear effects become more important. One sees the formation of fronts (Fig. 9). The reduced shock length (from the center) is equal to $`0.003`$ and is thus still smaller than the spatial step $`\delta r/\lambda =10^2`$. At the center, the positive swings of the pressure are now much larger than the negative ones (Fig. 10). There are several sources for non linearities: the equation of state; the inertial term in Euler equations; both enhanced by the spherical geometry. The relative importance of these factors will be discussed now. We have performed various simulations where we suppressed one term or the other. Results are summarized in Table 1.

If we take a constant sound speed $`c_\mathrm{s}=c_{\mathrm{st}}`$, i.e. a linear equation of state, then non-linear effects are strongly reduced and the maximal and minimal pressure excursions are almost symmetrical.

If we rather suppress the inertial term $`uu`$ from the Euler equation for the velocity, non-linear effects are also reduced, but to a lesser extent.

Thus all non-linear terms reinforce each other, but the dominant non-linear effect comes from the equation of state (the geometry also plays a crucial role of course, which cannot be separated from the other effects).

This is confirmed if we do the analytic calculation of section 2.2.2 again, now with a constant sound speed. We find then that

$$r_{\mathrm{shock}}L_0\mathrm{exp}\left[\frac{c_{\mathrm{st}}^2}{L_0\omega \mathrm{\Delta }v_0}\right].$$
(14)
On Fig. 11, we compare with the result (11) of section 2.2.2 obtained with the full equation of state. It turns out that, for realistic oscillation amplitudes, the shock appears almost at the center if the nonlinearities of the equation of state are not taken into account.

### 2.4 Shock formation at large amplitudes. The WENO scheme

Eventually, when the amplitude $`\mathrm{\Delta }x_0`$ is further increased, the formation of shocks can be observed in our simulations. The numerical scheme described in the previous section becomes unstable when $`r_{\mathrm{shock}}\delta r`$, i.e. for a vanishing static pressure and a reduced $`\delta r=10^2`$, when $`\mathrm{\Delta }x_06\mathrm{nm}`$. Then a new numerical scheme has to be used.

For performing fine analysis of the Euler flow dynamics, the numerical scheme must recover low dissipative property. Non-dissipative high-order accurate schemes (like spectral or Padé schemes) have been identified as suitable tools as far as regular numerical solutions are searched. Nevertheless, when dealing with compressible flows involving discontinuities, non-dissipative high-order schemes introduce spurious oscillations in the vicinity of the discontinuity and one must use a numerical scheme which can both represent the smooth regions of the solution with the minimum of numerical dissipation, and capture the discontinuities by using an ad hoc scheme with a robust discontinuity-capturing features. Therefore, as shock waves may occur in the computational domain for the present calculations, the numerical method we use is based on a high-order intrinsically-dissipative scheme originally designed to capture discontinuities. The method, called WENO (Weighted Essentially Non-Oscillatory) jiang\_s96 , is presented in the next section, and some more technical details about WENO and ENO schemes can be found in the literature jiang\_s96 ; tenaud\_g\_s00 ; shu\_o89 ; shu\_o88 ; roe81

As in many schemes devoted to computations involving shocks, the WENO scheme uses the Riemann invariants as variables. Computing the evolution of these variables requires to know their values not only at integer space coordinates, but also at some intermediate locations. An extrapolation from integer positions is thus necessary and this is where the fundamental idea of WENO schemes comes in. The time evolution scheme (a third order Runge-Kutta) and the change of variables towards the characteristic plane are more classical, but we shall still recall them for self-consistency.

#### 2.4.1 The WENO method in more details

For simplicity, the governing equations (3) are recast in the following abridged form:

$$\frac{Q}{t}=(r,Q)$$
(15)
where $`Q=(\rho ,\rho u)^t`$ is the vector of the conservative variables and

$$(r,Q)=\frac{F\left(Q\right)}{r}+𝒮(r,Q)$$

stands for a spatial operator, applied on $`Q`$, based on both the Euler flux vector $`F\left(Q\right)=(\rho u,\rho u^2+P)^t`$ and the source term $`𝒮(r,Q)={\displaystyle \frac{\theta }{r}}(\rho u,\rho u^2)^t`$.

In view of the discretization of the Euler equations (15), we will denote by $`\delta t`$ and $`\delta r`$ the time step and cell width respectively. $`Q_i^n`$ will denote the numerical vector solution at a time $`t=t_0+n\delta t`$ and at a position $`r=i\delta r`$. For simplicity, the integration of these equations has been performed by means of a decoupled time and space algorithm.

a) Time integration

The time integration is then performed by means of a third-order accurate Runge-Kutta method, proposed by Shu and Osher shu\_o89 , chosen because this high order accurate scheme does not increase the Total Variation of the right-hand-side of the equations ($`(r,Q)`$). When dealing with discontinuities, this Total Variation Diminishing (TVD) property is important since it ensures that no local extremum can be created during the time integration, meaning that non oscillation may occur in the shock wave vicinity due to the time scheme. At each point of the computational grid, the time integration is then obtained via a multi-step algorithm as follows:

$`Q^{(0)}`$ $`=`$ $`Q_i^n`$
$`Q^{(1)}`$ $`=`$ $`Q^{(0)}+\delta t\left(Q^{(0)}\right)`$
$`Q^{(2)}`$ $`=`$ $`{\displaystyle \frac{3}{4}}Q^{(0)}+{\displaystyle \frac{1}{4}}Q^{(1)}+{\displaystyle \frac{1}{4}}\delta t\left(Q^{(1)}\right)`$ (16)
$`Q^{(3)}`$ $`=`$ $`{\displaystyle \frac{1}{3}}Q^{(0)}+{\displaystyle \frac{2}{3}}Q^{(2)}+{\displaystyle \frac{2}{3}}\delta t\left(Q^{(1)}\right)`$
$`Q_i^{n+1}`$ $`=`$ $`Q^{(3)}`$
The explicit Runge-Kutta scheme exhibits a stability condition based on the $`CFL`$ number: $`CFL=\mu \delta t/\delta r`$, where $`\mu `$ is the maximal value of the eigenvalues $`\mu _i^k`$, to be defined below. All the calculations presented herein have been obtained by considering $`CFL=0.5`$, which corresponds to a nearly optimal value for the considered Runge-Kutta scheme. Let us mention that this value ensures a good representation of all the time scales of the flow considered.

b) Spatial integration

The spatial discretization of the right-hand-side term $`(r,Q)`$ of equations (15) is obtained by means of a high-order finite difference scheme:

$$(r_i,Q_i^n)=\frac{1}{\delta r}\left[\overline{F}_{i+1/2}^n\overline{F}_{i1/2}^n\right]+𝒮(r_i,Q_i^n)$$
(17)
where $`\overline{F}_{i+1/2}^n`$ is the numerical flux evaluated at the cell interface ($`r_{i+1/2}`$). To reconstruct the numerical flux at the cell interfaces, a scheme with a discontinuity-capturing feature must be employed to prevent oscillation in the vicinity of the shock wave. Following a previous study tenaud\_g\_s00  on the capability of some recent high-order shock capturing schemes to recover basic fluid mechanic phenomena, the numerical flux has been evaluated by means of Essentially Non-Oscillatory (ENO) family scheme jiang\_s96 ; shu\_o89 ; shu\_o88 . The numerical flux is approximated by means of polynomial reconstruction over several grid points (the set of these points is named “stencil”) around the cell interface considered. We shall now describe this reconstruction of the fluxes in details.

$``$ Change of variable: For simplicity and accuracy purposes, the discretization of the Euler flux is based on a polynomial reconstruction applied on the local characteristic variables (Riemann invariants, see § 2.2.1) since the equations recover the scalar form. Then, the propagation directions can easily be followed in the characteristic plane. In order to perform this change of variables, the method is the same as in Appendix A. First, one linearizes the Euler equations and finds the eigenvectors of the Euler flux jacobian evaluated at the cell interface ($`\left(F/Q\right)_{i+1/2}`$). Note that, as we are now using the Euler equations in conservative form (3) instead of (2), the jacobian differs from the one given in appendix A (matrix $`A`$), and reads now

$$F/Q=\left(\begin{array}{ccc}0& & 1\\ \frac{j^2}{\rho ^2}+c_\mathrm{s}^2& & 2\frac{j}{\rho }\end{array}\right)$$
(18)
In order to compute $`\left(F/Q\right)_{i+1/2}`$. \- and then the eigenvalues ($`\mu _{i+1/2}^k)_{k=\pm }`$, and the left ($`𝐥_{i+1/2}^k`$) and right ($`𝐫_{i+1/2}^k`$) eigenvectors -, the conservative variables $`Q_i^n`$ must be evaluated at the cell interface $`r_{i+1/2}`$ (note that a bold $`𝐫_{i+1/2}^k`$ refers to the eigenvector, while $`r_{i+1/2}`$ stands for the scalar spatial coordinate). As these variables do not vary linearly in the cell if a shock is present, one cannot use a simple arithmetic or geometric average, but rather a Roe average, whose description can be found in Ref. roe81 . This ensures the consistency of the scheme, i.e. that $`\left(F/Q\right)_{i+1/2}`$ converges towards $`\left(F/Q\right)_i`$ when $`\delta r`$ tends to zero.

The numerical Euler flux is then projected onto the left eigenvector matrix ($`𝐥_{i+1/2}^{(+)},𝐥_{i+1/2}^{()}`$). The scalar ENO reconstruction procedure is applied to the projected fluxes shu\_o89 . In the physical domain, the numerical Euler flux are then obtained by a projection onto the right eigenvectors ($`𝐫_{i+1/2}=𝐥_{i+1/2}^1`$) and read:

$$\overline{F}_{i+1/2}=\underset{k=\pm }{}\left[f_{}^{ENO}{}_{i+1/2}{}^{k}𝐫_{i+1/2}^k\right]$$
(19)
where $`f_{}^{ENO}{}_{i+1/2}{}^{k}`$ stands for the scalar ENO reconstruction, which will be defined below.

$``$ Extrapolation of variables at non-integer locations - the core of (W)ENO schemes:

First, we shall present the strategy used in ENO (Essentially Non–Oscillatory) schemes shu\_o89 . For a given non integer location, there are several ways to perform the extrapolation, depending on from how many integer points it will be performed, and how these points will be located with respect to the non-integer one. The set of these points is called a stencil. Fig. 12 illustrates the different possible choices for a stencil of length $`p=3`$. Simple finite difference schemes use a stencil defined once for all. Here, only the length $`p`$ of the stencil is fixed for a given simulation. All the $`p+1`$ possible locations are considered as candidates, provided that there is at least one point adjacent to the non-integer point (see again Fig. 12).

If we denote by $`p`$ the order of the reconstruction, the ENO procedure shu\_o89  chooses the most regular stencil among the $`p+1`$ stencil candidates. As an exemple for $`p=3`$, we can see all the stencil candidates on Fig. 12.

A first selection among the $`p+1`$ stencil candidates is performed according to the sign of the two eigenvalues ($`\mu _{i+1/2}^k`$): one keeps the $`p`$ most left stencils for the positive or null eigenvalues and the $`p`$ most right otherwise (Fig. 12). Indeed, the sign of the eigenvalue gives the propagation direction of the associated characteristics, and thus of the relevant information.

The regularity of the function on each of the $`p`$ remaining stencils is measured by the undivided difference Table shu\_o89  evaluated on each stencil and the most regular stencil is chosen among all the $`p`$ stencil candidates. Of course, this stencil may be different at each time step, for each location, and for each eigenvalue $`(\mu _{i+1/2}^k)_{k=\pm }`$.

The scalar ENO reconstruction is then applied on this specific stencil by means of the following polynomial development:

$$f_{}^{ENO}{}_{i+1/2}{}^{k,m}=\underset{j=0}{\overset{p1}{}}\zeta _j^{p,k}𝐥_{i+1/2}^kF\left(Q_{i+m+j}\right)$$
(20)
The integer $`m`$ refers to the index of the most left point of the chosen stencil. The sum goes over all the points of this stencil. One recognizes the projection of the fluxes onto the left eigenvectors, that allows to go from physical variables to characteristic variables. The constant coefficients of the polynomial ($`\zeta _j^{p,k}`$) are calculated in order to recover a scheme of order $`p`$ in regular regions. The values of $`\zeta _j^p`$ can be found in Table 2 up to $`p=5`$.

Using this generic ENO scheme (19, 20), the Euler flux derivative is estimated with a $`p^{th}`$-order of accuracy at best (in regular regions). However, when the stencil used at the cell interface $`r_{i+1/2}`$ is different from the one at $`r_{i1/2}`$ (which is the case in strong gradients or shock regions), the order of accuracy decreases.

One of the drawbacks in the generic ENO scheme is the necessity to check and choose between $`p`$ stencil candidate, which is quite CPU consuming. To overcome this disadvantage, we preferred using a Weighted ENO scheme (WENO) since it improves the order of accuracy of the generic ENO scheme by using a weighted combination of the $`p`$ possible stencils. The weights jiang\_s96  depend on the degree of regularity of the solution. In regular regions, they can be computed to achieve $`(2p1)^{th}`$-order of accuracy whereas in regions with discontinuities they are set to zero, leading to a standard ENO scheme. The WENO flux is estimated by:

$$f_{}^{WENO}{}_{i+1/2}{}^{k}=\underset{n=0}{\overset{p1}{}}\omega _nf_{}^{ENO}{}_{i+1/2}{}^{k,np+1}$$
(21)
where $`f_{}^{ENO}{}_{i+1/2}{}^{k,np+1}`$ is the generic ENO reconstruction given by equation (20) and $`\omega _n`$ are the weights defined as follows:

$`{\displaystyle \underset{n=0}{\overset{p1}{}}}\omega _n=1,`$ (22)
$`\omega _n={\displaystyle \frac{\beta _n}{{\displaystyle \underset{l=0}{\overset{p1}{}}}\beta _l}}`$ with $`\beta _n={\displaystyle \frac{C_n^p}{\left(\epsilon +IS_n\right)^2}}`$ (23)
$`\epsilon `$ is a small positive number to avoid denominator to be zero (hereafter we set $`\epsilon =10^6`$) and $`IS_n`$ is a measure of the flux function regularity for the $`n^{th}`$ ENO stencil candidates. The evaluation of the smoothness measurement ($`IS_n`$) is based on the undivided-differences and the $`IS_n`$ formulation can be found in jiang\_s96 . It is such that a more regular curve gives a smallest $`IS_n`$, and thus a largest weight $`\beta _n`$. The $`C_n^p`$ coefficients are reported in Table 3 up to the order $`p=5`$, for positive eigenvalues at the cell interface ($`\mu _{i+1/2}^\pm >0`$). For the negative eigenvalue case, the WENO coefficients can be obtained by symmetry with respect to the considered cell interface ($`r_{i+1/2}`$).

Using this WENO scheme (21-23), the Euler flux derivative is estimated with a $`(2p1)^{th}`$-order of accuracy at best (in regular regions). Moreover, let us underline that, if the solution is regular enough, the WENO procedure recovers a high-order centered scheme, which is of course non-dissipative.

If the length of the stencil is $`p=3`$ \- which is our case, the order of the scheme in regular regions is then 5. However, as the order of the scheme drops to 1 in shock regions - as for any other scheme - it is not interesting to increase too much the order in regular regions. As the scheme is less dissipative with $`p=4`$ in the non regular regions, this could even make the results worse (indeed, in our case, $`p=4`$ gave less satisfying results), since spurious oscillations may occur.

c) Boundary conditions

To solve the system of equations we need boundary conditions at the center of the sphere ($`r=0`$) and at the ceramic surface ($`r=L(t)`$). At the center of the sphere, since the velocity of the fluid is an antisymmetric quantity, the fluid is at rest ($`u|_{r=0}=0`$). Moreover, as the pressure and the density are symmetric quantities, their radial derivatives are set to zero ($`{\displaystyle \frac{}{r}}(P,\rho )|_{r=0}=0`$). The singular behavior of the source terms involves a specific treatment at the sphere center. While the velocity tends to zero when the radius tends to zero, the ratio $`{\displaystyle \frac{\rho u}{r}}`$ tends to a limit which has to be determined though $`{\displaystyle \frac{\rho u^2}{r}}`$ tends to zero. As $`{\displaystyle \frac{\rho u}{r}}`$ is a symmetric quantity, its limit (reached at $`r=0`$) is calculated by writing that its radial derivative vanishes at the sphere center by using a second order upwind difference:

$`{\displaystyle \frac{}{r}}\left({\displaystyle \frac{\rho u}{r}}\right)|_{r=0}`$ $`=`$ $`{\displaystyle \frac{1}{2r}}\left[3\left({\displaystyle \frac{\rho u}{r}}\right)_{i=0}+4\left({\displaystyle \frac{\rho u}{r}}\right)_{i=1}\left({\displaystyle \frac{\rho u}{r}}\right)_{i=2}\right]`$ (24)
$`=`$ $`0.`$
At the ceramic surface, one has to impose a condition corresponding to the motion of the ceramics. In our case, as the displacement of the sphere does not reach a sonic velocity, the two eigenvalues $`\mu ^+`$ and $`\mu ^{}`$ are of opposite sign. This means that the two informations necessary to determine the two unknowns $`\rho `$ and $`u`$ come from opposite directions. In particular, at the boundary, one comes from the interior of the domain and the other from the outside. Thus, it is possible to prescribe one of the variable, and the other one will be determined by an upwind scheme, i.e. a scheme asymmetric towards the inner domain. In the present problem, it is much more natural to prescribe the velocity at the ceramic surface:

$$u(L(t),t)=\omega \mathrm{\Delta }x_0\mathrm{sin}(\omega t)\left[1\mathrm{exp}(t/1.5)\right].$$
(25)
$`\mathrm{\Delta }x_0`$ is the amplitude of the displacement of the ceramic and $`\omega `$ is the ceramic pulsation. These two parameters are, of course, prescribed. The exponential term represents the response time of the transducer (see Sec. 2.1).

#### 2.4.2 Numerical results

Of course, all the results presented in section 2.3.2 can be reproduced with our WENO scheme. Besides, one can increase further the oscillation amplitude. For $`\mathrm{\Delta }x_0=7\mathrm{nm}`$, the reduced shock length (from the center) is equal to $`0.043`$, which is now larger than the spatial step $`\delta r/\lambda =10^2`$. One clearly sees the formation of sharp fronts (Fig. 13). At the center, the positive swings of the pressure are now very sharp compared to the negative ones (Fig. 14). The shape of the negative swings is also clearly asymmetric in time.

It would be tempting, in order to reduce the computational effort, to assume that non-linearities play a role only in the last wavelengths. Then, one could simulate a reduced box with radius $`L_{\mathrm{red}}`$, and take as an input condition :

$$\mathrm{\Delta }x_{\mathrm{red}}=\mathrm{\Delta }x_{\mathrm{exp}}\frac{L_{\mathrm{exp}}}{L_{\mathrm{red}}}$$
(26)
where $`\mathrm{\Delta }x_{\mathrm{exp}}`$ is the experimental oscillation amplitude of the transducer, and $`L_{\mathrm{exp}}`$ its radius. Actually, this is fine as long as non-linearities do not play a role at all, i.e. as long as the focal pressure is sinusoidal. In all other cases, as this is shown on Fig. 15, non-linearities are built through the whole propagation process, and one cannot neglect them even far from the center, without modifying the focal pressure. This could also be seen in the shock length analytic expression (11), which is directly proportional to the cell radius. In Fig. 15, we have performed simulations with various cell radii, corresponding to the same experimental oscillation amplitude on the ceramic. As the cell radius increases, both the positive and negative pressure swings decrease (in absolute values). Non-linearities make it more difficult to reach extreme values. Thus in all the simulations presented in this paper, we have simulated the whole experimental cell, with radius $`8`$ mm.

We have seen on Fig. 14 that positive pressure peaks can reach tremendously high values. Of course one may wonder whether this is physical. The obvious answer is no - but it is worth to discuss this point in some details.

The solution of the problem as we defined it in Section 2.1 becomes singular when the shock reaches the center of the sphere, and positive peaks actually diverge in the simulations as the spatial discretisation step $`\delta r`$ is decreased (see Fig. 16). For each fixed $`\delta r`$, one still finds a finite focal pressure, as it is defined as an average over the central cell - with radius $`\delta r/2`$. We check on Fig. 17 that the minimal pressure does converge for a decreasing $`\delta r`$.

A first remark is that this singularity involves only a very small region around $`r=0`$ (ultimately, it is singular only in $`r=0`$). In the experiment, one cannot measure the pressure exactly in $`r=0`$, but rather over the whole region reached by the laser beam (see section 3). The intensity of the laser beam through its cross section is Gaussian, with a waist (half width) equal to 7 microns. In order to take this averaging effect into account in the simulation, we have also computed a spatial average value of the pressure weighted by a Gaussian of waist $`7\mu m`$ (Fig. 18). Positive pressure peaks are lowered, while the negative swing is not very sensitive to the averaging process, as gradients are much weaker than during the positive swing. But above all, one now obtains extrema which are converging as $`\delta r0`$ (Fig. 16).

But it is not enough to rule out the singularity by an averaging process. Actually, the singularity is not expected to hold as such in a more realistic description. First, the third Euler equation for energy should also be taken into account in this regime, as well as regularization mechanisms (dispersion, dissipation). Secondly, in the experiment, we expect diffraction, and the fact that actually the flux is not zero at the focal point, to break the geometrical symmetry. Making quantitative estimates of these effects is not simple, and is postponed to future work. Finally, even with averaging, positive pressure peaks in Fig. 18 are still so high that they are far above the solidification pressure (25.3 bar at T=0). As we shall see in the last section, our recent experiments show that indeed acoustic waves can trigger crystallization, not only cavitation chavanne\_b\_c01 . In our simulations, we are not able to treat the possible crystallization.

As a conclusion, the value of the positive pressure swings is not expected to be reliable for high amplitudes of the transducer surface. But we checked, using different types of numerical regularization of the shocks (an example is given in the next section), that even if it affects the maximal pressure, it has no effect on the negative swings, and thus it is still possible to use our simulations to draw conclusions for negative pressures.

In Figs. 19 and 20, we show the density profile and focal pressure obtained for $`\mathrm{\Delta }x_0=30`$ nm, a value which can be reached in cavitation experiments. Fig. 19 illustrates how shocks are formed when the wave arrives near the focal point. The amplitude of the shock increases when the shock itself arrives at the focal point, leading to tremendous pressure maxima in the simulations (here about 12000 bars). This is of course unphysical, as we just discussed it, but the important point is that the minimum (negative) pressure does not depend on the value of the maximal pressure. Thus it is still possible to calculate a minimum pressure value and compare with experiments.

Fig. 21 shows the amplitude of the transducer oscillation which is necessary to reach a given minimal pressure at the center, starting from a static pressure $`P_{\mathrm{st}}`$. Actually, the variable we plot is rather $`P_{\mathrm{st}}`$ as a function of $`\rho _{\mathrm{st}}.\mathrm{\Delta }x_0`$, in order to see the departure from the linear theory

$$P_{\mathrm{min}}=P_{\mathrm{st}}+\omega ^2L_{\mathrm{exp}}\rho _{\mathrm{st}}\mathrm{\Delta }x_0.$$
(27)
Our interest in such a curve came from a first version of the experiment, for which it was not possible to measure the focal pressure. Still, we were able to measure the oscillation amplitude necessary to obtain cavitation for various initial static pressures, i.e. we could plot a curve such as those of Fig. 21, for the special case $`P_{\mathrm{min}}=P_{\mathrm{cav}}`$, the cavitation pressure. However this could be done only for positive static pressures, while we were interested in the zero amplitude value of the curve, for which static pressure and cavitation pressure should be the same. Thus simulations are useful to determine which kind of extrapolation should be used for negative static pressure values.

It is interesting to see that the effect of non-linearities is to bend such curves in a concave way (Fig. 21). In Ref. caupin\_b01 , the sign of this curvature was used to show that a linear extrapolation provides an upper bound of the cavitation pressure. Whether the shock formation affects the nucleation mechanism is an open question. If we try to plot the same oscillation amplitude as a function of the static density in the cell (Fig. 22), instead of the static pressure (Fig. 21), non-linearities are even more pronounced, as would have been expected from the equation of state (a concave function of a concave function is still more concave).

More important for the validation of our theoretical methods is that we succeeded later in measuring at the focal point the temporal signal itself. As explained below, this was done in a quasi-spherical geometry and very good agreement was found between theory and experiments.

## 3 Experiments

A hemispherical piezoelectric transducer is held against a clean glass plate. In a first approximation, the glass reflects the sound wave so that this is equivalent to a full spherical geometry. The main interest of the glass plate is that it allowed us to measure the instantaneous density at the center from the reflection of light at the glass/helium interface. Indeed, the reflectance depends on the refraction index of liquid helium which depends on its density as is well known from the Clausius-Mossoti relation.

### 3.1 Experimental method

Our experimental method is described in full details elsewhere chavanne\_b\_c01c ; chavanne\_b\_c02 . Let us only summarize it here. The transducer radius is $`8\mathrm{mm}`$ and its thickness is $`2\mathrm{mm}`$. It resonates in a thickness mode at $`f`$ = 1.019 MHz. At this frequency, it has a minimum impedance $`Z`$ = 22 Ohm. Its quality factor is $`Q`$ = 50 $`\pm `$ 5 when immersed in liquid helium at 25 bar. We usually pulse it with bursts of 6 oscillations. The $`300\mathrm{cm}^3`$ experimental cell is full of liquid helium and attached to the mixing chamber of a dilution refrigerator. We can run the experiment between 30 mK and 1.5 K, at static pressures from 0 to 25 bar. Fig. 23 shows our optical setup: a brass piece holds a lens whose focal length is 21.8 mm in liquid helium and a wedged glass plate (20 mm in diameter, $``$ 2 mm thick). The transducer is pressed against the plate. All the space inside is filled with liquid He. Thanks to the lens, the radius of the laser waist is reduced from 320 $`\mu \mathrm{m}`$ to 7 $`\mu \mathrm{m}`$. This means that the spatial resolution is about 14 $`\mu \mathrm{m}`$, the diameter of the optical focal region. This is small compared to the size of the acoustic focal region which is set by the acoustic wavelength at 1 MHz: from 240 $`\mu \mathrm{m}`$ at 0 bar to 360 $`\mu \mathrm{m}`$ at 25 bar. The distance of the lens to the glass plate has been carefully adjusted to have the laser focused at the glass/helium interface, on the transducer side. This was checked from the parallelism of the reflected beam. The 2 degrees wedge of the glass plate avoids interferences with reflections on its front face. A 1.7 mm hole in the transducer allows the transmitted light to be analyzed on the other side of the cryostat (see Fig. 24).

We use a single mode $`\mathrm{Ar}^+`$ laser. It is operated around 10 mW where its stability is best. In front of the laser is an absorber which reduces the light entering the cell (Fig. 24). We use the reflection on the front face of the absorber to monitor possible drifts of the laser power, but its stability is better than 0.5 % per day. For low temperature measurements, we reduced the dissipation in the cell with a small electromechanical shutter; it was synchronized by the pulse generator so that the cell was illuminated during 10 msec only, around the arrival time of each acoustic pulse on the glass plate. Two moving mirrors allow us to translate the laser beam vertically and horizontally. The last mirror is mounted on a rotating plate, so that the angle of incidence of the beam can be adjusted as well. These rotations correspond to translations of the optical focus in the focal plane of the lens. By successive operations, the laser spot was brought to the center of the acoustic focal region. The final adjustment was obtained by maximizing the modulation of the reflected signal by the acoustic wave.

The transmitted light is collected by a photomultiplier and used to detect nucleation events one by one. The light which is reflected at the glass/helium interface is separated from the incident beam by means of a semi-transparent plate and directed towards a photodiode. We used either a Hamamatsu C5331-03 avalanche photodiode (the“APD”), for the detection of the ac-modulation of the reflected beam, or a Hamamatsu S1406 silicon photodiode (the“SPD”) for the detection of the dc-component. The output from the photodiodes is digitized with a LeCroy 9344 CM oscilloscope at 1 GS/s with 8-bit resolution (6.5 effective bits due to the clock jitter at 1 GS/s).

The ac-component of the signal is related to the modulation of the helium density by the acoustic wave, and it is at most a few percent of the dc part which is related to the static density. The acoustic transmission into the glass being very small, we have neglected its effect on the reflected light. In order to achieve a 1% accuracy on the acoustic wave amplitude, we reduced the noise by averaging on 10000 sound bursts with a repetition rate of 1 to 10 Hz. There are two main sources of noise. The first one is photon noise in the reflected light. For a 5 $`\mu `$W power, this quantum noise is typically 10 nW, much more than the resolution we need. The other noise source is the oscilloscope jitter and it has a comparable amplitude. After averaging, the signal is well enough extracted from the noise as shown on Fig. 25.

### 3.2 Calibration

The intensity of the reflected light is proportional to the normal reflectance R at the glass/helium interface:

$$R=\left(\frac{n_\mathrm{g}n}{n_\mathrm{g}+n}\right)^2$$
(28)
where $`n_\mathrm{g}=1.5205`$ is the refractive index of glass for 514.5 nm green light. As for the index n of helium it is given by the Clausius Mossoti relation:

$$\frac{n^21}{n^2+2}=\frac{4\pi \rho \alpha _M}{3M}$$
(29)
where M = 4.0026 g, $`\rho `$ is the helium density, and $`\alpha _M`$ = 0.1245 $`\mathrm{cm}^3\mathrm{mol}^1`$ is the molar polarizability for the same green light. Note that $`\alpha _M`$ increases slightly as a function of frequency from its zero frequency value $`\alpha _{M0}=0.1233\mathrm{cm}^3\mathrm{mol}^1`$, as explained by successive authors polar .

We proceed as follows. We first measure the static pressure in the cell. Knowing the equation of state $`P(\rho )`$, we obtain the static density maris91 . From the static density, we calculate the normal reflectance in the absence of modulation by the wave. This is our reference. We then measure the ratio of the ac- to the dc-component of the reflected light, and obtain the amplitude of the density modulation in the acoustic wave. Unfortunately, we cannot do this with a single diode, so that we have to calibrate the ratio of the respective gains of our two photodiodes. This is achieved in the overlap of their bandwidths, with the help of an acousto-optic modulator operated at 200 kHz. We also have to check the linearity of the APD. Its gain is found constant up to incoming powers of 6 $`\mu W`$ where it starts decreasing. Most of the time, we use the APD in its linear regime; otherwise, a small correction is applied.

We finally have the following sources of uncertainties: the static density can be known within $`2`$ to $`\mathrm{4\hspace{0.17em}10}^2\mathrm{kg}\mathrm{m}^3`$. About the same uncertainty comes from the determination of the base line of the APD signal. Uncertainties in the gain ratio and in the APD measurement lead to a 1% uncertainty in the wave amplitude. We could finally check our calibration by studying heterogeneous nucleation of bubbles on the glass plate. We observe various nucleation mechanisms. One of them occurs at saturated vapor pressure ($`P_{\mathrm{sv}}=0\mathrm{bar}`$ in the low temperature limit) where the liquid density is $`145.13\mathrm{kg}\mathrm{m}^3`$. In a series of measurement at a static pressure $`P_{\mathrm{st}}=4.30\mathrm{bar}`$, we found cavitation at $`145.15\mathrm{kg}\mathrm{m}^3`$; in another series of measurements at $`P_{\mathrm{st}}=2.95\mathrm{bar}`$, we found cavitation at $`145.12\mathrm{kg}\mathrm{m}^3`$. This illustrates the final uncertainty in our measurements.

### 3.3 Comparison with calculations

Fig. 25 shows two recordings obtained at 0.1K with respective excitation amplitudes of 9.05 and 20.4 V on the transducer. In the cell, the static pressure is 9.80 bar, corresponding to a static density $`\rho =158.51\mathrm{kg}\mathrm{m}^3`$. The density oscillation is found asymmetric at large amplitude: negative swings are broader with a smaller amplitude than positive swings. Moreover, the negative swings are not symmetric in time. We have chosen this recording at intermediate pressure and moderate amplitude because the signal shape is not modified by any nucleation of crystals or bubbles.

We can compare the experimental recordings with the numerical calculations described in sections 2.3 and 2.4. Simulations are performed with the same static pressure as in the experiment. Then there is only one free parameter in the simulation (the oscillation amplitude $`\mathrm{\Delta }x_0`$) in order to adjust both the amplitude and the shape of the signal. The adjustment with the experimental signal is made only on the central oscillation of the latter. Indeed, in the experiment, the transducer is excited with an electrical burst of six oscillations. Since the transducer has a finite quality factor $`Q`$ 50, the amplitude of the sound wave increases during six periods and slowly decreases afterwards. The numerical result is used in the steady regime, after the initial transient. As can be seen, we find a very good agreement for numerical oscillation amplitudes $`\mathrm{\Delta }x_0=3`$ and $`6`$ nm: both the asymmetry with respect to the horizontal axis and the asymmetry in time are well reproduced by the calculation.

We make similar adjustments with several recordings, for different oscillation amplitudes and different static pressures. Our results are summarized on Fig. 26. For each adjustment, the numerical amplitude obtained by fitting the central oscillation is associated to the experimental voltage applied to the transducer. On Fig. 26, we also indicate as a second x-axis the estimate given by (37) for the oscillation amplitude $`\mathrm{\Delta }x_0`$. As this estimate is obtained independently from the simulations, as this will be detailed below, we refer to it as the experimental oscillation amplitude in the figure. If the estimate and the numerical adjustment were in perfect agreement, one would expect a slope equal to one in Fig. 26. Actually, we find that both methods give different results, and this will be discussed below. The adjustment between experimental signals and numerical simulations yields the following calibration :

$$\frac{\mathrm{\Delta }x_0}{V}=0.30\pm 0.02\mathrm{nm}\mathrm{V}^1.$$
(30)
It is interesting to compare this value with an estimate from the measurement of the electrical characteristics of the ceramic. Let us summarize the derivation of this estimate, which is given in Ref. caupin\_b01 .

Indeed, the quality factor is simply related to the ratio of the acoustic energy $`E_{\mathrm{ac}}`$ which is stored during one period to the average dissipated power $`V^2/(2Z)`$:

$$Q=4\pi f\frac{E_{\mathrm{ac}}Z}{V^2},$$
(31)
and the acoustic energy can be evaluated as follows. Let us call $`R`$ the mean radius of the transducer and $`2e`$ its thickness. For a resonance in a thickness mode, one can assume that the sound wave inside the transducer is a spherical wave with an amplitude proportional to $`1/r\mathrm{sin}[(k(rR)]\mathrm{sin}(\omega t)`$ where the wavevector $`k=\pi /2e`$ and $`\omega =2\pi f`$. For a spherical wave, the elastic displacement inside the transducer is

$$u=u_0\frac{R}{r}\mathrm{sin}[k(rR)]\mathrm{sin}(\omega t).$$
(32)
Since $`E_{\mathrm{ac}}`$ is twice the average kinetic energy in the transducer caupin\_b01  and the local velocity is simply the time derivative of $`u`$, one can integrate over the thickness and write:

$`E_{\mathrm{ac}}=\pi \rho _\mathrm{t}\omega ^2u_0^2R^2e`$ (33)
where $`\rho _\mathrm{t}`$ is the transducer density. After expressing $`E_{\mathrm{ac}}`$ in terms of the maximum displacement $`\mathrm{\Delta }x_0=u_0R/(Re)`$ of the inner surface, we obtain

$$E_{\mathrm{ac}}=\frac{M}{4}\omega ^2\mathrm{\Delta }x_{0}^{}{}_{}{}^{2}\frac{(Re)^2}{R^2+\frac{e^2}{3}}$$
(34)
where $`M=7.510^3\mathrm{kg}`$ is the mass of our transducer. Note that the last factor was forgotten in Ref. caupin\_b01 . Finally, we can express the displacement $`\mathrm{\Delta }x_0`$ as a function of the applied voltage:

$$\mathrm{\Delta }x_0=V\left(\frac{2Q}{\omega ^3MZ}\right)^{1/2}\frac{Re}{R}\sqrt{1+\frac{e^2}{3R^2}}$$
(35)
The above equation leads to $`\mathrm{\Delta }x_0/V=1.7\mathrm{nm}\mathrm{V}^1`$. However, this value would correspond to an excitation with long bursts, when the stored energy saturates. Since we excite it with bursts of six oscillations only, and since the phase is such that they start with positive swings, the maximum pressure is reached (6 + 1/4) periods after time zero, and the displacement to be considered in our case is

$$\mathrm{\Delta }x_0=\mathrm{\Delta }x_{0}^{}{}_{\mathrm{}}{}^{}\left[1\mathrm{exp}(12\pi /Q)\right]\left[\mathrm{exp}(\pi /2Q)\right]$$
(36)
We find $`\mathrm{\Delta }x_0=0.513\mathrm{\Delta }x_{0}^{}{}_{\mathrm{}}{}^{}`$ with $`\mathrm{\Delta }x_{0}^{}{}_{\mathrm{}}{}^{}`$ given by Eq. 35, and our final prediction is

$$\frac{\mathrm{\Delta }x_0}{V}=0.87\mathrm{nm}\mathrm{V}^1$$
(37)
The above value has the right order of magnitude but it is three times more than given by the fit of our numerical calculations (Eq.30). There are several assumptions in the above analysis which can be claimed as responsible for this discrepancy. We list them starting with those that we expect to be the most relevant:
a - the resonance in a thickness mode may be coupled to flexion modes, in which case the efficiency of the transducer can easily be reduced.
b - there is a small hole in the center of the transducer which allows the transmitted light to be analyzed.
c - the sound wave in the transducer cannot be strictly spherical, since there must be edge effects near its free equator.
d - the reflexion by the glass plate is not perfect so that our closed hemispherical geometry is not strictly equivalent to a full spherical geometry. Once more this should reduce the efficiency of the transducer.
e - some of the emitted energy is lost in the various pieces which hold it in the cell.
f - there is also some uncertainty of order 10 $`\%`$ in the measurement of $`Z`$ and $`Q`$.

As a result, we consider the value $`0.3\mathrm{nm}\mathrm{V}^1`$ as a very useful calibration of the efficiency of our transducer, in qualitative agreement with a simple estimate.

## 4 Conclusions and perspectives

In this article, we have presented analytic and numerical calculations of the focusing of a spherical acoustic wave. We have shown that shocks are generated in this geometry and we have obtained an analytic estimate for the shock length based on the characteristics method. Then, in order to perform full numerical simulations of the focussing process, we have used a WENO scheme to treat shocks. We then showed that our method is validated by a comparison with experimental measurements in a quasi-spherical geometry. We have measured a wave distortion which is well reproduced by our calculation and the analysis of its dependence on the excitation amplitude has led us to a very useful calibration of the efficiency of our transducers.

We consider this work as a first step only, and we plan to extend it to a hemispherical geometry for two important reasons. Indeed, in order to study the homogeneous nucleation of bubbles in stretched fluids or that of crystals in pressurized fluids, we need to eliminate the effect of walls. This is achieved by using hemispherical transducers which focus acoustic waves away from any walls lambare98 ; caupin\_b01 . In such experiments, since we have no probe in the acoustic focal region where nucleation takes place, there is a difficult problem of calibration of the sound amplitude, for which any reliable calculations would be very useful. Of course, the calculation in a hemispherical geometry is much more difficult because it is two-dimensional (it depends on both the radial distance and the polar angle). Now that the method is known for the treatment of shocks, the 2-D calculations should be tried. Furthermore, we have estimated the amplitude of non-linear effects in the hemispherical geometry caupin\_b01 , and found them much smaller than in the spherical geometry, though lower pressures seem to be reached. This is interesting in itself and should be tested numerically. One physical explanation could be that the local condition at the center is different: by symmetry, the spherical geometry imposes that the center is a node for the fluid velocity. In the hemispherical geometry, there is no reason why it should be so. On the contrary, the sound wave could even create a flow at the center with non vanishing averaged value. This phenomenon is known in the literature as acoustic streaming. This symmetry difference might lead to a different amplitude for the non-linear effects. It would be very interesting to study this phenomenon numerically.

Another direction of research would deal with the interaction between shocks and nucleation. Until now, all theories predicting the nucleation threshold completely ignore the presence of very steep gradients. This is not necessarily justified.

Acknowledgments: We are grateful to H. Lambaré for his contribution to the calculations in their early stage. C.A. would like to thank Frederic Coquel for interesting discussions.

Appendix : theoretical prediction of the shock length by the characteristics method

Appendix A: Characteristic equations

The system of Euler equations (2) is of the form

$$_\mathrm{t}v_i+A_{ij}_rv_j=b_i,$$
(38)
with $`𝐯(\rho ,u)`$ and $`𝐛=(\theta \rho u/r,0)`$. The matrix

$$A=\left(\begin{array}{cc}u& \rho \\ \frac{c_\mathrm{s}^2}{\rho }& u\end{array}\right)$$

has two eigenvalues $`\mu _+=u+c_\mathrm{s}`$ and $`\mu _{}=uc_\mathrm{s}`$ associated with the left eigenvectors

$`𝐥_+`$ $`=`$ $`\left(\begin{array}{c}c_\mathrm{s}/\rho \\ 1\end{array}\right),`$
$`𝐥_{}`$ $`=`$ $`\left(\begin{array}{c}c_\mathrm{s}/\rho \\ 1\end{array}\right).`$ (39)
If we apply $`𝐥_k`$ on the left of equation (38), we obtain

$$𝐥_k.\left[\frac{d𝐯}{dt}𝐛\right]=0,$$
(40)
where

$$\frac{d}{dt}=_\mathrm{t}+\mu _k_r.$$

Thus the derivative $`d/dt`$ is taken along a curve $`r=x(t)`$ with slope $`dx/dt=\mu _k`$ everywhere ($`\mu _k`$ itself being a function of $`r`$ and $`t`$ via $`v`$). By definition, this curve is called a k-th characteristic and is denoted by $`𝒞^{(k)}`$. There is a whole family of $`𝒞^{(k)}`$ characteristics, covering the whole space, each curve being determined, for example, by the initial conditions.

The leftmost term of equation (40) can be integrated and the equation becomes

$$\frac{d}{dt}\left[C_\mathrm{o}(\rho \mathrm{ln}\rho )\pm u\right]+\theta \frac{C_\mathrm{o}(\rho 1)u}{r}=0,$$
(41)
where again the time-derivative is taken along a characteristic curve.

We call Rieman invariants the quantities

$$I_\pm C_\mathrm{o}(\rho \mathrm{ln}\rho )\pm u$$
(42)
appearing in eq. (41).

Appendix B: Lower bound for the shock length

Considering the geometry and notations of Sec. 2.2.2, we calculate at what time the first characteristic $`𝒞_0^{()}`$ emitted by the piston at $`t=0`$ is cut by another characteristic $`𝒞_i^{()}`$. This is the signature of a shock forming. The calculation assumes that it is the first shock ever formed in the cell.

As $`𝒞_0^{()}`$ is first cut by characteristics $`𝒞^{()}`$ emitted at early times, and almost parallel to $`𝒞_0^{()}`$, we perform the following change of variable. Instead of $`r`$, the location of any point will be given by its distance $`\eta `$ to a point moving on the characteristic $`𝒞_0^{()}`$, taken at the same time :

$$\eta =L_0c_{\mathrm{st}}tr.$$
(43)
The characteristic equations (41) have to be written in the new coordinates $`(\eta ,t)`$. Besides, we would like to eliminate $`u`$ and $`\rho `$ from the equations so that the only remaining unknowns would be the Rieman invariants $`I_+`$ and $`I_{}`$. Following greenspan\_n93 , we have $`d\eta /dt=(u\pm c_\mathrm{s}+c_{\mathrm{st}})`$ and $`u=(I_+I_{})/2`$. But unlike Ref. greenspan\_n93 , the density $`\rho `$ -and thus the sound velocity- cannot be expressed in a simple way in terms of $`I_+`$ and $`I_{}`$. Indeed, we would need to invert the relation

$$I_++I_{}=2C_\mathrm{o}(\rho \mathrm{ln}\rho )2I(\rho ).$$
(44)
At this stage, we will only assume that the inversion can be performed and write

$$\rho =^1\left(\frac{I_++I_{}}{2}\right).$$
(45)
The characteristic equations now read

$`_\mathrm{t}I_+\left\{{\displaystyle \frac{I_+I_{}}{2}}+C_\mathrm{o}\left[^1\left({\displaystyle \frac{I_++I_{}}{2}}\right)1\right]+c_{\mathrm{st}}\right\}_\eta I_+`$
$`+{\displaystyle \frac{2}{Lc_{\mathrm{st}}t\eta }}\left({\displaystyle \frac{I_+I_{}}{2}}\right)C_\mathrm{o}\left[^1\left({\displaystyle \frac{I_++I_{}}{2}}\right)1\right]=0`$
(46)
$`_\mathrm{t}I_{}\left\{{\displaystyle \frac{I_+I_{}}{2}}C_\mathrm{o}\left[^1\left({\displaystyle \frac{I_++I_{}}{2}}\right)1\right]+c_{\mathrm{st}}\right\}_\eta I_{}`$
$`+{\displaystyle \frac{2}{Lc_{\mathrm{st}}t\eta }}\left({\displaystyle \frac{I_+I_{}}{2}}\right)C_\mathrm{o}\left[^1\left({\displaystyle \frac{I_++I_{}}{2}}\right)1\right]=0`$
(47)
We search for a solution under the form

$`I_+(\eta ,t)`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}I_+^{(m)}(t)\eta ^m,`$
$`I_{}(\eta ,t)`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}I_{}^{(m)}(t)\eta ^m.`$ (48)
We choose the lowest order terms $`I_+^{(0)}`$ and $`I_{}^{(0)}`$ equal to their value in the fluid at rest $`I_+^{(0)}=I_{}^{(0)}=C_\mathrm{o}(\rho _{\mathrm{st}}\mathrm{ln}\rho _{\mathrm{st}})`$. Then if we write equations (46-47) at lowest order, all the terms in the second equation vanish. In the first equation, both the first and last terms disappear. The remaining term leads to

$$I_+^{(1)}=0.$$
(49)
This is not very surprising, as $`I_+`$ corresponds to the characteristics moving from the fluid at rest into the perturbed region. Taking (47) to the next order, we obtain an equation for $`I_{}^{(1)}`$

$$\frac{dI_{}^{(1)}}{dt}+\frac{1}{tL_0/c_{\mathrm{st}}}I_{}^{(1)}+𝒦\left[I_{}^{(1)}\right]^2=0,$$
(50)
where

$$𝒦\frac{1}{2}\left(\frac{2\rho _{\mathrm{st}}1}{\rho _{\mathrm{st}}1}\right).$$
(51)
This equation can be solved using a change of variables $`\overline{I_{}}1/I_{}^{(1)}`$. The solution yields

$$I_{}^{(1)}(\tau )=\frac{I_{}^{(1)}(0)}{(1\tau )\left[1𝒦\frac{L_0}{c_{\mathrm{st}}}I_{}^{(1)}(0)\mathrm{ln}(1\tau )\right]},$$
(52)
with

$$\tau \frac{c_{\mathrm{st}}}{L_0}t.$$
(53)
We have now to determine the initial value $`I_{}^{(1)}(0)`$. It refers to small $`t`$, rather than $`t`$ exactly equal to zero.

To determine its value, we calculate for small time $`t`$ the variation of $`I_{}`$ between point $`A=(r,t)=(L_0c_{\mathrm{st}}t,t)`$ which sits on the first characteristic $`𝒞_0^{()}`$, and point $`B=(r,t)=(r_p(t),t)`$ where $`r_p(t)`$ is the location of the piston at time $`t`$.

At $`A`$, the fluid is at rest. We have $`I_{}=I_{}^{(0)}`$ and $`\eta =0`$.

At $`B`$, for $`t`$ small, $`rL_0`$ and thus $`\eta c_{\mathrm{st}}t`$. Besides, $`I_{}=C_\mathrm{o}(\rho \mathrm{ln}\rho )u`$ where $`u`$ is equal to the velocity of the piston $`v_p(t)=\mathrm{\Delta }v_0\mathrm{sin}(\omega t)\mathrm{\Delta }v_0\omega t`$. If we expand $`\rho (r_p(t),t)=\rho _{\mathrm{st}}+\alpha t`$ (where $`\alpha `$ is unknown), then replacing into $`I_{}`$ and expanding in $`t`$ yields

$`I_{}`$ $`=`$ $`C_\mathrm{o}(\rho _{\mathrm{st}}\mathrm{ln}\rho _{\mathrm{st}})+C_\mathrm{o}\alpha tC_\mathrm{o}{\displaystyle \frac{\alpha }{\rho _{\mathrm{st}}}}t+\mathrm{\Delta }v_0\omega t`$
$`=`$ $`I_{}^{(0)}+\left({\displaystyle \frac{c_{\mathrm{st}}}{\rho _{\mathrm{st}}}}\alpha +\mathrm{\Delta }v_0\omega \right)t.`$
Comparing this relation with the expansion

$$I_{}=I_{}^{(0)}+\eta I_{}^{(1)}I_{}^{(0)}c_{\mathrm{st}}I_{}^{(1)}t,$$
(54)
gives

$$I_{}^{(1)}=\frac{\alpha }{\rho _{\mathrm{st}}}\frac{\mathrm{\Delta }v_0}{c_{\mathrm{st}}}\omega .$$
(55)
On the other hand, an expansion of (45) in powers of $`\eta `$ gives

$$\alpha t=\rho _{\mathrm{st}}\rho _{\mathrm{st}}=\frac{c_{\mathrm{st}}}{\rho _{\mathrm{st}}}\left(\eta \frac{I_{}^{(1)}}{2}\right)t$$
(56)
The elimination of $`\alpha `$ between the above two equations yields

$$I_{}^{(1)}(t=0)=2\frac{\mathrm{\Delta }v_0\omega }{c_{\mathrm{st}}}.$$
(57)
The time at which $`I_{}^{(1)}`$ becomes infinite (see equation 52) gives an upper bound $`t_{\mathrm{shock}}`$ for shock formation. It is only an upper bound because some other terms of the expansion (48) in $`\eta `$ may explode before $`I_{}^{(1)}`$. We find

$$t_{\mathrm{shock}}\frac{L_0}{c_{\mathrm{st}}}\left\{1\mathrm{exp}\left[\frac{c_{\mathrm{st}}^2}{2L_0\omega \mathrm{\Delta }v_0}\frac{\rho _{\mathrm{st}}1}{\rho _{\mathrm{st}}\frac{1}{2}}\right]\right\}\frac{L_0}{c_{\mathrm{st}}}.$$
(58)
As the corresponding shock distance $`r_{\mathrm{shock}}`$ is measured from the center of the sphere, a lower bound for $`r_{\mathrm{shock}}`$ is

$$L_0r_{\mathrm{shock}}L_0\mathrm{exp}\left[\frac{c_{\mathrm{st}}^2}{2L_0\omega \mathrm{\Delta }v_0}\frac{\rho _{\mathrm{st}}1}{\rho _{\mathrm{st}}\frac{1}{2}}\right]>0.$$
(59)
All the above calculations are valid for small $`\eta `$, i.e. only for characteristics not too far from $`𝒞_0^{()}`$. These are the characteristics emitted by the initial motion of the piston.


format: text

text:
# Vortex-antivortex configurations and their stability in mesoscopic superconducting square.

## I  Introduction

Recently it was shown that the influence of boundaries can lead to stabilization of the vortex-antivortex molecules in mesoscopic samples. Analysis of the linearized Ginzburg-Landau equation (GLE) has shown that such molecules appear at particular values of the external magnetic field depending on the sample shape and size. The solution of the GLE in the limit of the extreme type-II superconductor shows that such molecules have a very shallow minimum in the free energy and are very sensitive to the change of the sample shape.

In a square mesoscopic thin film with the total vorticity $`m=3`$ the symmetric solution with four vortices and one antivortex is the solution of the linearized GLE with the lowest free energy. According to ref. , away from the $`H_{c2}`$ line the giant vortex with vorticity $`m=3`$ is stable and has the lowest free energy. This implies that a topological phase transition without change of the vorticity and without a reduction of the symmetry should take place with the change of the external field or/and the coherence length away from the critical-field line.

It was proposed that in the limit $`\kappa 1/\sqrt{2}`$, where the vortex-antivortex interaction changes the sign, vortex-antivortex complexes should be more stable for the cylindric sample shape. For a thin film with $`\lambda _{eff}=\lambda ^2/d`$ and $`\kappa _{eff}=\lambda _{eff}/\xi `$, where $`d`$ is the thickness of the film, $`\lambda `$ is the London penetration depth, and $`\xi `$ is the superconducting coherence length, this mechanism should be less effective due to a smaller contribution of the ’magnetic energy’ to the total free energy.

We performed an extensive study of the region of the phase diagram where the vortex-antivortex phase was previously reported for the sample of the square shape. We focused to the region $`4<a/\xi <8`$ and $`\kappa _{eff}>0.05`$. We found that the antivortex phase is stable in a broad range of parameters. The region of stability of the phase does not depend strongly on the value of the parameter $`\kappa _{eff}`$. The energy gain due to the antivortex formation is much smaller then the energy difference between two phases with different vorticities. The giant vortex with $`m=3`$ is unstable for any field, $`\xi /a`$ and $`\kappa _{eff}0.1`$. Phase transition to the phase with three separated vortices takes place when $`\xi /a`$ is driven away from the critical field line. The reduction of $`\kappa _{eff}`$ does not stabilize the antivortex phase for the thin film sample in the contrast to the case of the cylindric sample ref..

## II Formalism and Solution

GLE for the normalized complex order parameter $`\psi =\mathrm{\Psi }/\mathrm{\Psi }_0`$, $`\mathrm{\Psi }_0=\sqrt{\beta /\left|\alpha \right|}`$ has the following form:

$$\xi ^2(i+\frac{2\pi 𝐀}{\mathrm{\Phi }_0})^2\psi \psi +\psi |\psi |^2=0$$
(1)
here $`\xi =\frac{\mathrm{}^2}{4m|\alpha |}`$, $`\alpha `$ and $`\beta `$ are the temperature dependent parameters of the Ginzburg-Landau expansion for the free energy, $`\mathrm{\Phi }_0`$ is the flux quantum, $`𝐀`$ is the vector potential and $`𝐇=\times 𝐀`$ the magnetic field. The second GLE equation for the vector potential reads:

$$\times \times 𝐀=i\frac{\mathrm{\Phi }_0}{4\pi \lambda ^2}(\psi ^{}\psi \psi \psi ^{})\frac{|\psi |^2𝐀}{\lambda ^2}.$$
(2)
In addition to Eq.(1) we assume the boundary condition for the superconductor-insulator junction on the sample edges:

$$(i+\frac{2\pi 𝐀}{\mathrm{\Phi }_0})𝐧\psi =0,$$
(3)
where $`𝐧`$ is the vector normal to the surface of the sample.

As it was described in ref. we introduce $`N\times N`$ discrete points on the square and rewrite Eq.(1) in the form of the nonlinear discrete Schrödinger equation:

$$\underset{𝐥}{}t_{𝐢+𝐥,𝐢}\psi _{𝐢+𝐥}ϵ(𝐢)t_{𝐢,𝐢}\psi _𝐢\psi _𝐢+\psi _𝐢|\psi _𝐢|^2=0,$$
(4)
where the summation index $`𝐥=(\pm 1,0)`$, $`(0,\pm 1)`$ points toward the nearest neighbours and $`t_{𝐢_\mathrm{𝟏},𝐢}=(\xi N/a)^2exp(i\varphi _{𝐢_\mathrm{𝟏},𝐢})`$ and $`\varphi _{𝐢_\mathrm{𝟏},𝐢}=\frac{2\pi }{\mathrm{\Phi }_0}_{𝐫_𝐢}^{𝐫_{𝐢_\mathrm{𝟏}}}𝐀(𝐫)𝑑𝐫`$. The boundary conditions are included in the discrete nonlinear Schr¨odinger equation as in ref. where $`\psi _𝐢=0`$ if $`𝐢`$ is outside of the sample and $`ϵ(𝐢)=4\delta _{i_x,1}\delta _{i_x,N}\delta _{i_y,1}\delta _{i_y,N}`$ where $`𝐢=(i_x=1,\mathrm{},N,i_y=1,\mathrm{},N)`$.

After discretization of Eq.(2) we can obtain the exact expression for the vector potential:

$$A_𝐢^v=\underset{𝐧}{}K(𝐢𝐧)J_𝐧^v$$
(5)
where

$`J_𝐢^v`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_0a}{4\pi \lambda _{eff}N}}\mathrm{}(exp(i\varphi _{𝐢+𝐥_𝐯,𝐢})\psi _𝐢^{}\psi _{𝐢+𝐥_𝐯}`$ (7)
$`exp(i\varphi _{𝐢𝐥_𝐯,𝐢})\psi _𝐢^{}\psi _{𝐢𝐥_𝐯})`$
where $`v\{x,y\}`$ and $`𝐥_𝐱=(1,0)`$, $`𝐥_𝐲=(0,1)`$

$$K(𝐧)=\frac{N}{2\pi ^2a}_0^\pi 𝑑x𝑑y\frac{cos(n_xx)cos(n_yy)}{\sqrt{42cos(x)2cos(y)}}$$
(8)
The numerical self consistent solution of the problem is obtained by iterating the solution of the nonlinear equation for the order parameter Eq.(4) and calculations of the current and the vector potential Eqs.(5,6). We used two ways of solving Eq.(4). The first is similar to that reported in ref. and corresponds to the iterative solution of the linearized Eq.(4). The second relies on the fact that Eq.(4) represents the Euler equation for the free-energy functional whith included boundary conditions. Eq.(4) was therfore solved by the direct minimization of the corresponding functional using the conjugate-gradient method. Both techniques gave identical results.

## III Results

The main goal of the paper is to investigate the phase diagram in the region $`4.5<\mathrm{\Phi }/\mathrm{\Phi }_0<6.5`$ and $`(a/\xi )^2<60`$ where the solution with one antivortex and four vortices (Fig.1) has been reported. We found that the region of the phase diagram where the symmetry induced antivortex solution has the lowest energy is broader than expected from the solution of the linearized GLE. As it is shown in Fig. 2 for $`\kappa _{eff}=\mathrm{}`$ the antivortex phase is stable up to $`(a/\xi )^255`$, depending on $`\mathrm{\Phi }/\mathrm{\Phi }_0`$. For a finite $`\kappa _{eff}`$ this region shifts to the higher field as $`(a/\xi )^2`$ increases (see Fig.2).

The interesting behaviour is observed when the external field is fixed and $`(a/\xi )^2`$ increases. Close to the $`H_{c2}`$ the lowest minimum of the free energy corresponds to the solution with the vorticity $`m=41`$ with the antivortex in the center of the square. Present calculations do not confirm the existence of the giant-vortex solution with $`m=3`$ in this region of the phase diagram as reported previously . The difference is due to increase of the number of discrete points $`N`$ enabling detection of the antivortex. With increase of $`(a/\xi )^2`$ away from the $`H_{c2}`$ line the phase transition to the multivortex state with the same vorticity ($`m=3`$) and a lower symmetry takes place (see Fig. 2). In general, the free energy depends on the vorticity $`m=n_+n_{}`$ and the total number of vortices in the system $`n=n_++n_{}`$. The transition at $`(a/\xi )^255`$ and $`\mathrm{\Phi }/\mathrm{\Phi }_05.5`$ takes place at the constant vorticity $`m=3`$ with the change of $`n`$ from 5 to 3. The transition is therfore not only characterized by an order parameter, but also by the change of the number of vortices at the constant total vorticity $`m`$, suggesting that the transition is close to the first order. This statement is confirmed by the observation that above the transition point, $`(a/\xi )>(a/\xi )_{crit}`$, both solutions with $`m=3`$ and $`m=41`$ coexists (see Fig. 4). Since near the transition the free-energy difference between the phases with the same vorticity $`m`$ and different $`n`$ is small it is difficult to determine the phase boundary between phases with $`m=41`$ and $`m=3`$ accurately. The transition could be easier observed by calculating the two component order parameter $`\eta _x=x|\psi (x,y)|^2𝑑x𝑑y`$, $`\eta _y=y|\psi (x,y)|^2𝑑x𝑑y`$ shown in Fig.3. The transition point is given by the point where the order parameter goes to 0. At the same point the calculated sample magnetization changes slope as it is clearly seen from Fig.4.

Close to the $`H_{c2}`$ line the repulsion of vortices from the boundaries and attraction of the 4 vortices to the antivortex stabilizes the phase with $`m=41`$ and small vortex-vortex distances. At smaller value of $`\xi `$ the repulsion from the boundaries decreases and one vortex annihilates with the antivortex. As a result, the repulsion between the remaining vortices increases leading to an increase of the order parameter with a further decrease of $`\xi `$.

Increasing the field up to $`\mathrm{\Phi }/\mathrm{\Phi }_0=10.6`$ leads to the stabilization of the phase with total vorticity $`m=7`$. Near the $`H_{c2}`$ line similar to the phase with $`m=41`$ the solution with $`m=81`$ is realized. When $`(a/\xi )^2`$ increases in a complete analogy to the case with $`m=41`$ the second order phase transition to the phase with $`m=7`$ and a similar order parameter takes place. Here also, both solutions with $`m=81`$ and $`m=7`$ coexists above $`(a/\xi )_{crit}`$ indicating that the transition is close to the first order. We believe that the situation is quite general for the case of arbitrary $`m=4l1`$ for $`l=1,2,3,\mathrm{}`$.

At the end we would like to discuss the dependence of the stability of the antivortex phase at small $`\kappa `$. According to the arguments of Ref., at small $`\kappa `$ the vortex-vortex interaction changes the sign making the antivortex phase more stable. As a result, the average distance between vortices in the middle of the square increases as well. In order to verify this conjecture for the thin film sample we plot in Fig.5 the vortex-antivortex distance $`r_0`$ as a function of $`1/\kappa _{eff}`$. The distance decreases with the decreasing $`\kappa _{eff}`$. For $`\kappa _{eff}<0.1`$ the distance is smaller than the grid spacing $`a/N`$ so we can not resolve separate vortices. We find that $`r_0\mathrm{exp}(\mathrm{\Lambda }/\lambda _{eff})`$ with $`\mathrm{\Lambda }a`$. The situation is just opposite to that reported in ref.. We believe that in the case of the thin film of the square shape the reduction of $`\kappa `$ does not stabilize the phase with the antivortex.

It is interesting to note differences between samples of different shapes. For the cylindric shape the giant vortex phase with any vorticity is always stable close to the $`H_{c2}`$ line. According to the Ref. for the mesoscopic triangle the giant vortex state with $`m=2`$ is metastable and the solution with the antivortex ($`m=31)`$ is stable. For the case of the square shape the giant vortex solution with $`m=3`$ is never stable for $`\kappa _{eff}0.1`$. For $`\kappa _{eff}<0.1`$ the limited grid spatial resolution prevented us to distinguish the solution with the antivotrex from the possibly (meta)stable giant-vortex solution.

Finally, let us discuss the possibility to detect the state with the antivortex experimentally. Calculation of the magnetic field in the sample shows that the magnetic field has a local minimum in the center of the sample also for the giant vortex solution with $`m=2`$. The local minimum observed for the antivortex state with $`m=41`$ is therfore not due to the antivortex formation (Fig.6) but due to a particular distribution of the current in the sample. Therefore, imaging of the magnetic field distribution cannot provide an evidence for the antivortex. The magnetic field for the multivortex solution with $`m=3`$ has 3 well separated maxima that break the four fold rotational symmetry of the sample allowing a direct imaging of vortices. Since the antivortex state cannot be detected directly the observation of a hysteresis in the vicinity of the transition line from the $`m=41`$ antivortex state to the $`m=3`$ multivortex state could suggest that the symmetric phase is indeed the phase with the antivortex.

Figure captions

Fig. 1 (a) The magnitude of the superconductivity order parameter for $`(a/\xi )^2=35`$, $`\kappa _{eff}=\mathrm{}`$ and $`\mathrm{\Phi }=5.9\mathrm{\Phi }_0`$. (b) The central region where the vortices are located in an expanded scale. The position of the antivortex is indicated by the symbol $``$.

Fig. 2 The calculated phase diagram. Different phases are marked with icons schematically indicating the vortex pattern where the full dot represents a vortex, the open dot reperesents an antivortex and the larger full dot represents a double vortex. The full symbols and continous lines represent the phase boundaries for $`\kappa _{eff}=\mathrm{}`$ while the open symbols and dotted lines reperesent the phase boundaries for $`\kappa _{eff}=1`$. In the later case only the phase boundaries of the region with the total vorticity 3 are shown.

Fig. 3 The magnitude of the order parameter around the transition where one vortex anihilates with the antivortex as a function of $`(a/\xi )^2`$ at the constant magnetic field.

Fig. 4 The magnetic moment of the sample as a function of $`(a/\xi )^2`$ at the constant magnetic field.

Fig. 5 The vortex-antivortex distance (a) and the magnetic moment of the sample (b) as functions of the parameter $`1/\kappa _{eff}`$. In (a) error bars represent the grid spacing and the solid line the exponential fit discussed in the text.

Fig.6 The magnetic field in the film in the case of (a) giant vortex with $`m=2`$, (b) the antivortex solution with $`m=41`$ and (c) three separate vortices with $`m=3`$. Here $`H_0`$ is the external magnetic field.