††thanks: Corresponding Author

# An Ultrasonically Actuated Fine-Needle Creates Cavitation in Bovine Liver

Emanuele Perra<EMAIL_ADDRESS>Medical Ultrasonics Laboratory
(MEDUSA), Department of Neuroscience and Biomedical Engineering, Aalto
University, Espoo, 02150, Finland  Nick Hayward  Medical Ultrasonics
Laboratory (MEDUSA), Department of Neuroscience and Biomedical Engineering,
Aalto University, Espoo, 02150, Finland  Kenneth P.H. Pritzker  Department of
Laboratory Medicine and Pathobiology, University of Toronto, Toronto, M5S 1A8,
Canada  Department of Pathology and Laboratory Medicine, Mount Sinai Hospital,
Toronto, M5G 1X5, Canada  Heikki J. Nieminen  Medical Ultrasonics Laboratory
(MEDUSA), Department of Neuroscience and Biomedical Engineering, Aalto
University, Espoo, 02150, Finland

(August 27, 2024)

###### Abstract

Ultrasonic cavitation is being used in medical applications as a way to
influence matter, such as tissue or drug vehicles, on a micro-scale.
Oscillating or collapsing cavitation bubbles provide transient mechanical
force fields, which can, e.g., fractionate soft tissue or even disintegrate
solid objects such as calculi. Our recent study demonstrates that an
ultrasonically actuated medical needle can create cavitation phenomena inside
water. However, the presence and behavior of cavitation and related bioeffects
in diagnostic and therapeutic applications with ultrasonically actuated
needles are not known. Using simulations, we demonstrate numerically and
experimentally the cavitation phenomena near ultrasonically actuated needles.
We define the cavitation onset within a liver tissue model with different
total acoustic power levels. We directly visualize and quantitatively
characterize cavitation events generated by the ultrasonic needle in thin
fresh bovine liver sections enabled by high speed imaging. On a qualitative
basis, the numerical and experimental results show a close resemblance in
threshold and spatial distribution of cavitation. These findings are crucial
for developing new methods and technologies employing ultrasonically actuated
fine-needles such as ultrasound-enhanced fine-needle biopsy, drug delivery and
histotripsy.

††preprint: Emanuele Perra, JASA

##  I Introduction

In recent years, ultrasonic cavitation has emerged in various medical
applications as a way to influence matter non-invasively. Cavitation is a
phenomenon that can be described as the interaction between small spherical
gas bubbles and pressure perturbations taking place in a medium. When the peak
rarefactional pressure amplitude (PRPA) of an ultrasound field is low enough,
gas bubbles can undergo stable oscillations about their equilibrium radius,
which is usually referred to as stable cavitationNeppiras1980 . However, at
elevated PRPAs, if certain threshold conditions are metHolland1990 , gas
bubbles can collapse giving rise to transient cavitationApfel1981 . The
collapse of a cavitation bubble may generate different nonlinear acoustic
phenomena in the surrounding medium, such as generation of rapid liquid
microjets, acoustic emission in the form of shock waves and formation of high
stress fields. These physical effects have been widely investigated and
employed in different medical applications with the intent to, e.g., ablate
tumorsKennedy2005 , fractionate calculiMatsumoto2005  or tissueMaxwell2009a
and enhance the permeability of cells for drug delivery
applicationsLentacker2014 .

In our recent study, we have demonstrated that cavitation events can be
generated in water by an ultrasonically actuated medical needle Perra2021 .
Moreover, it has been shown that at $\sim$$30\text{\,}\mathrm{kHz}$,
ultrasound-enhanced fine-needle aspiration biopsy (USeFNAB) enhances the yield
of a biopsy by 3-5$\times$ in liver compared to when a fine-needle aspiration
biopsy (FNAB) procedure is conducted using a similar needle. These results
suggested that the nonlinear acoustic phenomena generated at the needle tip,
including cavitation, might play an important role in the tissue cutting
mechanism in the context of biopsy applications and beyond. However, the
potential presence of cavitation and related bioeffects in diagnostic and
therapeutic applications with ultrasonically actuated needles require a more
thorough understanding.

Actuation of medical needles by ultrasound is not a new concept, and a number
of studies related to the topic can be found in the literatureGui2014a ;
Kuang2016 ; Sadiq2014a ; Liao2014 ; Liao2013 . However, the applications have
been limited to improve the needle visibility in ultrasound-guided regional
anaesthesia and tissue biopsyGui2014a ; Kuang2016 ; Sadiq2014a  or to reduce
the penetration force of a standard needleLiao2014 ; Liao2013 . So far, no
research seems to have been conducted on studying the potential generation of
nonlinear acoustic phenomena and their interaction with soft tissue.

In this study, we aim at studying the influence of cavitation on soft tissue
under the action of an ultrasonically actuated needle. Numerical modeling is
first used to simulate the time-dependent acoustic field generated by the
ultrasonic needle and the cavitation bubble dynamics in a liver tissue model.
The numerical results provided fundamental understanding on the cavitation
nucleation threshold, spatial distribution and maximum size of the cavitation
bubbles, and their influence on the tissue, according to a cavitation/tissue
interaction model proposed by Mancia et al.Mancia2019 . Experimentally, we
developed a method to visualize cavitation bubbles in thin portions of fresh
liver tissue, involving high-speed (HS) imaging using light transmission. Such
understanding is crucial for optimising the safety and efficacy of clinical
interventional procedures, including many for diagnostics and cancer
treatmentsBlana2004a .

##  II Methods

###  II.1 Numerical Simulation

The computational software COMSOL Multiphysics v5.5 Multiphysics2020  was used
to solve the different equations governing the cavitation bubble dynamics
taking place in soft tissue. We assumed a scenario where a 21G $\times$
$80\text{\,}\mathrm{m}\mathrm{m}$ hypodermic needle is used, since it
represents a common medical needle normally employed in FNAB applications. The
needle is partially placed into a $10\text{\,}\mathrm{m}\mathrm{m}$ $\times$
$12\text{\,}\mathrm{m}\mathrm{m}$ cylinder representing a liver tissue sample
and actuated at the ultrasonic frequency of $33\text{\,}\mathrm{kHz}$ and
total acoustic power (TAP) of 0.2, 0.5 and $0.8\text{\,}\mathrm{W}$.
Ultrasonic flexural standing waves are enabled in the needle shaft via an
ultrasonic device with similar geometry and material properties to the one
employed in the actual experiments (Fig. 1).

The displacement field in the needle and the acoustic pressure field inside
the tissue domain are first calculated in the time domain (time step =
$\frac{1}{100f}=$0.3\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, simulated time =
$2/f$ = $60.60\text{\,}\mathrm{\SIUnitSymbolMicro s}$). In order to study the
bubble dynamics, the Keller–Miksis equationKeller1980  (Eq. (4)) is then
solved in the liver domain as a global ODE, by having the surroundings of the
needle seeded with cavitation nuclei periodically spaced by
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The initial size $R_{0}$ of the
nuclei was assumed to be $500\text{\,}\mathrm{nm}$, according to previous
studies on cavitation behaviour in soft tissueVlaisavljevich2015 . The
acoustic forcing pressure $p_{f}$ is given by the pressure field calculated in
a separate study step, hence, the bubble motion is assumed not to contribute
to the total pressure field. For simplicity, the inter-bubble interaction is
neglected.

The three-dimensional model was meshed with free tetrahedral elements,
considering at least 20 nodes per wavelength, which was considered an
appropriate number for minimizing the local approximation errors Thompson2006
. A detailed list of the model parameters is given in Table 1.

####  II.1.1 Acoustic Wave Propagation in Tissue

The acoustic wave propagation in soft tissue has been modeled by adopting the
Westervelt equationhamilton1998nonlinear :

$\displaystyle\nabla^{2}p-\frac{1}{c_{\infty}^{2}}\frac{\partial^{2}p}{\partial
t^{2}}+\frac{\delta}{c_{\infty}^{4}}\frac{\partial^{3}p}{\partial
t^{3}}+\frac{\beta}{\rho_{\infty}c_{\infty}^{4}}\frac{\partial^{2}p}{\partial
t^{2}}=0,$ (1)

where $p$ is the pressure, $c_{\infty}$ and $\rho_{\infty}$ are the speed of
sound and the density of the medium, respectively. The first two terms of Eq.
(1) describe the linear lossless propagation of sound in a medium, while the
third term is associated with viscous losses, and the last one accounts for
the nonlinear propagation in soft tissue. The sound diffusivity $\delta$ is
defined as Solovchuk2013 ; Demi2014 :

$\displaystyle\delta=\frac{2c_{\infty}^{3}\alpha}{\omega^{2}},$ (2)

where $\omega$ is the angular frequency, $\alpha$ is the acoustic absorption
coefficient. The term $\beta$ represents the nonlinearity coefficient,
expressed as:

$\displaystyle\beta=1+\frac{B}{2A},$ (3)

with $B/A$ being the nonlinearity parameter given for a specific material.

####  II.1.2 Cavitation Model

The Keller–Miksis equationKeller1980  was used to describe the bubble dynamics
in soft tissue:

$\displaystyle\left(1-\frac{\dot{R}}{c_{\infty}}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{\dot{R}}{3c_{\infty}}\right)\dot{R}^{2}=\frac{1}{\rho_{\infty}}\Biggl{(}1+\frac{\dot{R}}{c_{\infty}}+\frac{R}{c_{\infty}}\frac{d}{dt}\Biggr{)}\times\left[p_{B}-(p_{\infty}+p_{f}(t))-\frac{2S}{R}+J\right].$
(4)

In the above, $R$, $\dot{R}$ and $\ddot{R}$ denote the radial displacement,
velocity and acceleration of the cavitation bubble wall, respectively, and the
constants $c_{\infty}$ and $\rho_{\infty}$ denote the speed of sound and the
density of the medium. The driving pressure is expressed by $p_{f}(t)$, while
the pressure at the air-liquid interface of the bubble is defined asGaitan2005
:

$\displaystyle p_{B}=p_{0}\left(\frac{R_{0}}{R}\right)^{3\kappa},$ (5)

where $R_{0}$ the bubble radius at rest and $\kappa$ is the polytropic
exponent. The term $p_{0}$ represents the internal pressure of the bubble when
the bubble is at equilibrium, expressed as:

$\displaystyle p_{0}=p_{\infty}+2\frac{S}{R_{0}},$ (6)

where $p_{\infty}$ indicates the ambient pressure and $S$ the surface tension
of the bubble. Eq. (4) is combined with the Kelvin–Voigt modelGaudron2015 ,
which leads the integral of the deviatoric stress $J$, accounting for the
viscoelastic behaviour of soft tissue, to be expressed as followsMancia2017 :

$\displaystyle
J=2\int_{R}^{\infty}{\frac{\tau_{rr}-\tau_{\theta\theta}}{r}dr}=\frac{-4\mu\dot{R}}{R}-\frac{G}{2}\left[5-4\left(\frac{R_{0}}{R}\right)-\left(\frac{R_{0}}{R}\right)^{4}\right],$
(7)

where $\mu$ is the tissue viscosity and $G$ is the tissue shear modulus.
$\tau_{rr}$ and $\tau_{\theta\theta}$ represent the the radial and tangential
stresses, respectively, due to the bubble deformation. They are related as
follows:

$\displaystyle\tau_{rr}=-4\mu\frac{R^{2}\dot{R}}{r^{3}}+2G\left[\Biggl{(}\frac{r_{0}}{r}\Biggr{)}^{4}-\Biggl{(}\frac{r}{r_{0}}\Biggr{)}^{2}\right]=-2\tau_{\theta\theta},$
(8)

where $r$ is the radial coordinate and $r_{0}=\sqrt[3]{r^{3}-R^{3}-R_{0}^{3}}$
relates the coordinate $r$ to its initial position in the undeformed
configuration of the surrounding tissue. The strain field in the surrounding
tissue is defined asMancia2019 ; Estrada2017 :

$\displaystyle
E_{rr}=-2\text{ln}\left(\frac{r}{r_{0}}\right)=-2E_{\theta\theta},$ (9)

being $E_{rr}$ and $E_{\theta\theta}$ the radial and tangential strain,
respectively.

The model for cavitation/tissue interaction proposed by Mancia et al. (2019)
was adopted to estimate the amount of tissue volume influenced by the
cavitation activity. Specifically, this is evaluated by identifying the
regions where the von Mises strain (Eq. (10)) exceeds the ultimate fractional
strain measured for liver ($0.38\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$)Bamber1981 .

$\displaystyle
E_{mises}=\sqrt{\frac{2}{3}\left[E_{rr}^{2}+2+\left(-\frac{1}{2}E_{rr}\right)^{2}\right]}=|E_{rr}|.$
(10)

###  II.2 Experiments in ex vivo Tissue

####  II.2.1 Experimental Arrangement

A custom-built ultrasonic device Perra2021  was used to excite a flexural
vibration mode ($f=$33\text{\,}\mathrm{kHz}$$) in a 21G hypodermic needle
(length = $80\text{\,}\mathrm{mm}$) (model: 4665465, 100 Sterican, B Braun,
Melsungen, Germany) (Fig. 1a). The needle was coupled to an S-shaped
3D-printed aluminium waveguide (3D Step Oy, Ylöjärvi, Finland) that acts as a
mode converter, translating the longitudinal motion provided by the Langevin
transducer into a flexural motion of the needle (Fig. 1b). The ultrasonic
device was driven by an RF amplifier (model: AG 1012LF, Amplifier/Generator,
T&C Power Conversion, Inc., Rochester, NY, United States) controlled by a
function generator (model: Analog Discovery 2, Digilent, Inc., Henley Court
Pullman, WA, United States). The spatial coordinates of the needle were
controlled by using a motorized three-axis translation stage (model:
8MT50-100BS1-XYZ, Motorized Translation Stage, Standa, Vilnius, Lithuania).
The cavitation events induced by the needle action were filmed in liver tissue
using a HS camera (model: Phantom V1612, Vision Research, Wayne, NJ, United
States) in combination with a macro lens (model: Canon MP-E 65 mm f / 2.8 1-5x
Macro Photo, Canon Inc., Ōta, Tokyo, Japan). A collimated beam of light
(model: OSL2COL, Collimation Package for OSL2IR, Thorlabs, Inc., Newton, NJ,
United States), generated by a halogen fiber optic illuminator (model: OSL2IR,
High-Intensity Fiber-Coupled Illuminator, Thorlabs, Inc., Newton, NJ, United
States), was used to produce back-lit shadowgraph footages of the needle
actuation inside tissue.

####  II.2.2 Sample Preparation

The liver specimen (from a 26 months old female cow) was retrieved from the
slaughterhouse (Vainion Teurastamo Oy, Orimattila, Finland) within 2 h post
mortem and experiments were performed within 6 h post mortem at room
temperature (22–$24\text{\,}\mathrm{\SIUnitSymbolCelsius}$). The specimen was
first rinsed with 1×PBS (BP399-4, Phosphate Buffered Saline, 10× Solution,
Fisher BioReagents, Fisher Scientific, Hampton, NH, United States) to wash
away any excess of blood from its surface. Thin slices, approximately
$1\text{\,}\mathrm{m}\mathrm{m}$ thick, were carefully extracted from the
specimen by using a pair of microtome blades (12101840, Epredia Ultra
Disposable Microtome Blades, Epredia, Portsmouth, NH, United States) fixed to
a spacing of $1\text{\,}\mathrm{m}\mathrm{m}$ from each other. The liver
slices were further washed in 1×PBS, cut into
$2\text{\,}\mathrm{c}\mathrm{m}$$\times$$1\text{\,}\mathrm{c}\mathrm{m}$
portions and inserted into a custom made glass sample holder (Fig. 1c). The
sample holder was created by cutting a
$2\text{\,}\mathrm{c}\mathrm{m}$$\times$$1\text{\,}\mathrm{c}\mathrm{m}$
portion of glass from the upper part of a
$51\text{\,}\mathrm{mm}$$\times$$75\text{\,}\mathrm{mm}$ microscope slide
(J1800BMNZ,Epredia SuperFrost Plus Adhesion slides, Special Size, Epredia,
Portsmouth, NH, United States), which was placed between two intact microscope
slides, in order to form a pocket for the tissue sample (Fig. 1c,d).

####  II.2.3 Data Acquisition

Since the penetration depth $\delta$ of light into bovine liver is estimated
to be, for example, $1.44\text{\,}\mathrm{mm}$ for a wavelength of
$635\text{\,}\mathrm{nm}$Arslan2018 , the thickness of the sample and the
light source spectrum were considered appropriate to ensure a good visibility
of the needle inside tissue during the HS recordings. During the experiments,
the needle was first carefully inserted into the specimen at a depth of
$5\text{\,}\mathrm{mm}$ and penetration speed of
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$.
Ultrasound waves ($33\text{\,}\mathrm{kHz}$, pulse repetition frequency (PRF)
= $55\text{\,}\mathrm{Hz}$, duty cycle (DC) = $50\text{\,}\%$) at 3 different
TAP levels (0.2 W (n = 5), 0.5 W (n = 5) and 0.8 W (n = 5)) were then applied
to the device, while the needle movement inside tissue was recorded with the
following settings: sample rate =
$130\,000\text{\,}\mathrm{f}\mathrm{p}\mathrm{s}$, exposure =
$7.1\text{\,}\mathrm{\SIUnitSymbolMicro s}$, resolution = 256 pixels × 256
pixels, lens aperture = 2.8, spatial resolution =
$5.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{/}\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}$.

####  II.2.4 Data Analysis

The HS frames were analysed in MATLAB (R2020b) Mathworks2016  to quantify the
projected area of cavitation and the needle displacement using a similar
method presented in our recent publicationPerra2021 . A cross-correlation
based image registration was performed along the x-axis between the reference
frame $I_{1}$ and the the i-th frame $I_{i}$, in order to estimate the needle
displacement $\Delta x_{i}$ from its reference position. The image $I_{1}$ was
then rigidly translated by $\Delta x_{i}$ and thresholded with the Otsu
methodOtsu1996 , while an Otsu thresholding followed by a morphological
closing operation (circular structuring element, diameter = 7 pixels) was
applied to the image $I_{i}$. The segmented image $I_{cav,i}$ showing only the
cavitation activity is obtained by subtracting the binarized reference image
$I_{bw,1}$ from the closed image $I_{bw,i}$. Since the needle shape in
$I_{bw,1}$ does not perfectly match the one in $I_{bw,i}$, a final
morphological opening operation (circular structuring element, diameter = 3
pixels) was applied to the output $I_{cav,i}$ in order remove any pixels that
may have remained after the subtraction operation and that are not
representative of the cavitation activity.

Probability maps (Eq. (11)), showing the probability of cavitation manifesting
around the needle tip, and projected areas of cavitation activity over time
(Eq. (12)) were finally calculated as follows:

$\displaystyle P_{cav}=\frac{100}{N}\sum_{i=1}^{N}{I_{cav,i}},$ (11)
$\displaystyle A_{cav,i}=\int\int{I_{cav,i}dxdy},$ (12)

where $N=25000$ is the total number of frames.

Velocity maps, shear and strain rate maps were generated using the PIVlab
toolbox Thielicke2021 .

##  III Results

###  III.1 Simulation of Cavitation in Liver

Based on the simulations, the ultrasonic action of the needle induced
expansion of bubbles in the proximity of the needle tip, when the needle was
driven by ultrasonic waves at the frequency of 33 kHz (Fig. 3a.1). The TAP
levels were 0.2, 0.5 and $0.8\text{\,}\mathrm{W}$, which helped generate
driving pressure waves with peak amplitudes of 90, 500 and
$900\text{\,}\mathrm{kPa}$ evaluated at the location B1 (Fig. 3a.2). The
cavitation bubble dynamics evaluated at location B1 exhibited a TAP dependent
behaviour. At the lowest TAP level, the bubble oscillated around its initial
radius $R_{0}$ with velocities lower than
$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. At increased TAP levels, the
bubble radius expanded up to 45 and $100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$, reaching a maximum velocity of 30 and
$100\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ during the collapse (Fig.
3a.2) at employed TAP levels of 0.5 and $0.8\text{\,}\mathrm{W}$,
respectively.

The probability maps presented in Fig. 3b.1 show the chance for cavitation
events to occur around the needle tip, defined as when the radius $R$ of a
cavitation nucleus becomes as large as $2R_{0}$Gaitan2005 . According to this
criterion, no cavitation activity was detected at $0.2\text{\,}\mathrm{W}$,
while the probability of cavitation occurrence became higher at increased TAP
levels, being up to 50% at $0.8\text{\,}\mathrm{W}$.

The predicted volume of tissue influenced by cavitation activity, calculated
by identifying the regions where $E_{mis}>0.38$, was almost zero at low TAP,
suggesting that no important deformations were induced in the surrounding
tissue (Fig. 3b.2). However, higher strains can be generated in tissue at
higher TAP levels due to the elevated cavitation bubble activity, leading to
an increase in the influenced volume of tissue by cavitation up to
$3.2\text{\,}{\mathrm{mm}}^{3}$ after 2 acoustic cycles at
$0.8\text{\,}\mathrm{W}$.

###  III.2 Observation of Cavitation Events in ex vivo Liver

Thin slices of liver tissue were sonicated at different TAP levels (0.2 W (n =
5), 0.5 W (n = 5) and 0.8 W (n = 5)) with the ultrasonically actuated needle.
Fig. 4a shows some exemplary frames acquired with the HS camera, when a
halogen fiber optic light source was used to produce shadowgraph images of the
needle movement inside tissue. When the delivered TAP was
$0.2\text{\,}\mathrm{W}$, no cavitation activity was detected. At the TAP
level of $0.5\text{\,}\mathrm{W}$, the needle motion induced the formation of
cavitation bubbles, which mostly took place at the distal end of the needle.
However, when the delivered TAP was increased to $0.8\text{\,}\mathrm{W}$,
multiple cavitation bubbles can be noticed along the needle tip, extending to
a few hundreds of $\text{\,}\mathrm{\SIUnitSymbolMicro m}$ from the needle
boundaries along the directions parallel to the needle motion. Fig. 4b shows
the cavitation probability maps calculated across the entire duration of the
HS footages, which suggest that no cavitation events were observed at the
lowest TAP $0.2\text{\,}\mathrm{W}$ employed. By increasing the TAP, the
probability of seeing cavitation bubbles was up to 10% in the region within
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ from the needle tip along the
positive x-axis and $300\text{\,}\mathrm{\SIUnitSymbolMicro m}$ along the
negative z-axis, while this region became considerably greater in area and
uniformly distributed around the needle tip, when the highest TAP of
$0.8\text{\,}\mathrm{W}$ was employed.

Fig. 5b represents the time evolution of the needle tip peak displacement,
obtained by computing the moving maximum of the raw data and using a window
with a size of approximately 2 acoustic cycles
($60\text{\,}\mathrm{\SIUnitSymbolMicro s}$). In all experiments, the peak
displacement reached its maximum value within the first burst, being $\sim$ 9,
45 and $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at the TAP levels of 0.2,
0.5 and $0.8\text{\,}\mathrm{W}$, respectively. Fig. 5c shows the projected
area of the cavitation activity (filtered with a moving average, window size
$\sim$10 acoustic cycles) as a function of time. It can be noted that no
cavitation activity was present at $0.2\text{\,}\mathrm{W}$, while some
activity was detected at $0.5\text{\,}\mathrm{W}$, and, at
$0.8\text{\,}\mathrm{W}$, the measured cavitation activity was relatively
elevated.

Measurements of needle displacement and cavitation activity exhibited high
repeatability within the same power groups, as shown in Fig. 5d. In the first
burst, the needle tip peak displacements were $9.5\pm
6.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (average $\pm$ standard deviation,
n = 5), $34.3\pm 1.0\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $61.6\pm
1.9\text{\,}\mathrm{\SIUnitSymbolMicro m}$, when the employed TAP levels were
0.2, 0.5 and $0.8\text{\,}\mathrm{W}$, respectively. The needle tip peak
displacement stabilised within 3 bursts, reaching the values of $10.8\pm
8.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ ($0.2\text{\,}\mathrm{W}$),
$41.5\pm 1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ ($0.5\text{\,}\mathrm{W}$)
and $91.9\pm 2.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$
($0.8\text{\,}\mathrm{W}$) after 10 bursts. Fig. 5e shows the time integral of
cavitation activity calculated for each individual burst. In all experiments,
no cavitation activity was recorded for TAP = $0.2\text{\,}\mathrm{W}$, while
it slowly built-up over time at TAP = $0.5\text{\,}\mathrm{W}$, being
$0.014\pm 0.002\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$ after the first
burst, and reaching the value of $0.055\pm
0.012\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$ during the 10th burst. At
the highest TAP employed, the cavitation level observed during the first burst
was $0.177\pm 0.060\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$, and
reached its maximum intensity in the last burst ($0.386\pm
0.016\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$). Overall, the temporally
local peak displacement of the needle tip measured for each TAP level were
$11.3\pm 6.8\text{\,}\mathrm{\SIUnitSymbolMicro m}$, $44.4\pm
1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $97.2\pm
1.8\text{\,}\mathrm{\SIUnitSymbolMicro m}$, which led to total cavitation
activity values of $\sim$ $0\pm
0.\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$, $0.459\pm
0.032\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$ and $3.44\pm
0.78\text{\,}{\mathrm{mm}}^{2}\text{\,}\mathrm{ms}$ for TAP = 0.2, 0.5,
$0.8\text{\,}\mathrm{W}$, respectively.

Velocity maps were generated out of 2 consecutive frames of the HS videos in
order to estimate the velocity field of the tissue at the moment of a
cavitation bubble collapse. Fig. 6a shows the velocity vector field
distribution overlapped to a HS frame showing a cavitation event, when the
highest TAP is employed. The velocities are the highest at the very tip of the
needle, being approximately $3\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
in this region (Fig. 6b). Importantly, according to the simulation, the
velocity of the tissue-air -interface can be remarkably greater, i.e. up to
$100\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. However, the limited frame
rate adopted during the recordings (130000 fps) did not permit capture of the
very moment of the cavitation collapse, which resulted in underestimation of
its maximum velocity. Fig. 6c shows the shear rate distribution around the
needle, being the highest in magnitude ($20\text{\,}{\mathrm{ms}}^{-1}$) at
the proximity of the cavitation bubble boundary. This is reasonable since the
cavitation bubble deformation is known to exert considerably high shears and
stresses in the surrounding medium. In Fig. 6d the strain rates assume
negative values, which denote a compression state, on the left hand side of
the needle and positive values in the proximity of the cavitation bubble. The
needle movement is in the direction of the negative x-axis, which causes the
adjacent portion of tissue to be compressed on the left-hand side of the
needle, and to be stretched on the right-hand side.

##  IV Discussion

The results indicate that cavitation events can be triggered by actuating a
standard medical needle with ultrasonic flexural waves in liver tissue at the
frequency of $33\text{\,}\mathrm{kHz}$. The numerical results suggested that
the cavitation activity mostly took place at the needle tip, which was
optically confirmed with HS photography. This is explained by the flexural
vibration mode induced in the needle, which makes the needle oscillate with
its highest displacements at its tip, thus enabling higher pressure amplitudes
in this region. Since cavitation is a strictly related threshold phenomenon,
cavitation events are most likely to appear at the needle tip location, where
most of the acoustic intensity is concentrated. Moreover, due to the geometric
spreading of the acoustic wavefront, directed outwards from the needle shaft,
the acoustic intensity decays rapidly further away from the needle, hence
limiting the cavitation effects to the proximity of the needle tip. However,
by looking at a cross section of the numerical model, we noticed that a
greater number of cavitation nuclei seemed to interact with the acoustic field
within the needle cannula. This has to do with the concave shape of the inner
walls, which tends to focus the acoustic energy towards the center axis of the
needle. Also, the acoustic wave is reflected within the inner walls of the
needle, which might contribute to raising the magnitude of the acoustic
intensity in this region. This justifies the elevated cavitation activity in
this region, which causes the portion of tissue located at the needle opening
to be influenced the most (Fig. 3b.2). However, the volume of tissue
influenced by the cavitation activity seems to decrease drastically at lower
TAP levels.

The experimental results showed that the probability of triggering cavitation
events in soft tissue is a function of TAP, suggesting the existence of a
threshold ($$0.2\text{\,}\mathrm{W}$<$ TAP $<$0.5\text{\,}\mathrm{W}$$ ) for
enabling cavitation, when a standard medical needle is actuated in soft
tissue. Even though the cavitation threshold value was not confirmed through a
direct measurement, based on the simulation results, this threshold seems to
be approximately $-500\text{\,}\mathrm{kPa}$. Assuming an initial size of the
cavitation nuclei of $500\text{\,}\mathrm{nm}$, the natural frequency of a
bubble with this radius would fall in the MHz rangeAtchley1989 . Since the
excitation frequency $f_{0}$ = $33\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$
used in this study is far below the resonance frequency $f_{n}$ of the bubble,
the cavitation threshold criterion is governed by the Blake
pressureAtchley1989 , which determines the critical negative pressure below
which a cavitation event will occur:

$\displaystyle
P_{B}=P_{\infty}+\frac{8\sigma}{9}\sqrt{\frac{3\sigma}{2R_{B}^{3}(P_{\infty}+(2\sigma/R_{B}))}},$
(13)

where $P_{B}$ is the Blake pressure, $\sigma$ is the surface tension and
$R_{B}$ is the Blake bubble radius. Under these assumptions, the Blake
threshold for a bubble of $R_{0}=$500\text{\,}\mathrm{nm}$$ is $\sim$
$-200\text{\,}\mathrm{kPa}$, which is in line with the numerical results, in
which the estimated threshold was around $-500\text{\,}\mathrm{kPa}$.

The interpretation of the results here presented is of fundamental importance
in the context of different medical applications. In our recent study, the
influence of the ultrasonic action of a medical needle was exemplified in
liver tissue by comparing the yield mass collected with the USeFNAB technique
to the one obtained with the conventional FNAB approach. The major finding was
that, by increasing the TAP level, the yield of a liver biopsy was increased
up to 5$\times$as compared to when a standard FNA was performed, without
inducing major alterations to the sample quality up to TAP of 0.8 W. More
importantly, a TAP of $0.2\text{\,}\mathrm{W}$ was enough to increase the
biopsy yield by almost 2$\times$. Based on the findings of the present study
(Fig. 4 and Fig. 5), it seems that this TAP level is unlikely to generate
cavitation events in liver, indicating that the tissue yield increase observed
at this TAP can be in part associated with the tissue cutting mechanisms
arising from shear and hydrodynamic effects promoted by the ultrasonic
vibration of the needle tip, rather than being induced by cavitation. However,
higher TAP levels allowed us to obtain even larger tissue sample
massesPerra2021  as well as more frequent cavitation activity. Although a
clear correlation between the cavitation activity and the tissue yield is yet
to be proven, these observations do not exclude the possibility for cavitation
to be contributing to the enhancement of tissue collection. In fact, the high
strain rates generated in the proximity the gas/tissue interface can
potentially induce different viscoelastic mechanical responses (namely
stiffening, softening, hardening and tissue failure observed in porcine liver
under high strain rate compression testingsChen2009a ) that might facilitate
the tissue cutting mechanisms that yield an increase in sample extraction.

Regarding the safety aspects of the biopsy application in relation to the
potential cavitation-induced effects in tissue, one should consider the
mechanical index MI = $P_{r}/\sqrt{f_{c}}$, where $P_{r}$ is the peak negative
pressure (MPa) and $f_{c}$ is the excitation frequency (MHz). According to our
simulations, at the lowest TAP employed MI $\sim$ 0.4, which would ensure a
cavitation-free biopsy procedure, since cavitation is unlikely to take place
at MI $<$ 0.5Apfel1991 . At $0.5\text{\,}\mathrm{W}$ the MI is approximately
2.7, which will most likely induce the formation of cavitation bubbles; this
might impact on the safety, as at MI values greater than 1.9 potential
bioeffects might be induced in the tissueSen2015 . These bioeffects may
include cell lysis and extravasation of bloodHallow2006a ; Church2008a ;
Haar2010 . The highest TAP employed should be avoided for biopsy applications,
as the high MI value (4.9) suggests that bioeffects and tissue damage due to
the bubble collapses are likely to appear.

If uncontrolled, cavitation events can lead to deleterious effects in soft
tissue. However, this could be turned into a therapeutic advantage if one aims
to treat unhealthy tissue, such as tumors. At high levels (TAP $>$
$0.8\text{\,}\mathrm{W}$), the needle vibration is anticipated to cause the
formation of large clouds of cavitation bubbles and elevated tissue heating,
which may arise from the viscous friction forces that can appear at the bubble
surface. These effects can be used in medical applications such as tumor
ablation Zhou2011 , histotripsy or lithotripsy Yoshizawa2009 , where the
medical intent is to achieve a complete or partial destruction of the target
by mechanical and thermal means. Since a fine hypodermic needle is employed to
bring the acoustic energy directly into the target, one may be able to easily
access different locations inside the body to provide minimally invasive
treatment of solid organ cancers of the prostate Chen2016 , thyroid
Pacella2013  and pancreas Larghi2021 , and other lesions tooPiccioni2019 . In
addition, the cavitation phenomena generated with this technology could
potentially find use in other applications as a way to ultrasonically activate
sonosensitive carriers for the release of drugsGeorge2019 ; Husseini2005 ;
Zhou2014a , mediate drug or gene delivery into cellsDachs1997 , or to improve
the permeation of tissue allowing the entry of therapeutic agents (e.g. as
with ultrasonically mediated blood-brain barrier opening)Aryal2014 .

The limitations of this study include the inability to replicate a numerical
model representative of the real-world scenario. The equations adopted in the
simulations are highly parameter-dependent, and since some of the viscoelastic
and acoustic properties of tissue are largely unknown, some assumptions had to
be made, for example, on the tissue viscosity, surface tension and initial
radius of cavitation nuclei. For simplicity, the interactions between
individual bubbles were neglected as well as potential temperature related
effects, which, conversely, are known to play an important role in the context
of ultrasonic cavitation in soft tissue. Moreover, since the use of thin
portions of liver tissue were necessary to visualize the needle and cavitation
activity during the experiments, the setup might not replicate the same
acoustic and mechanical conditions as in an FNAB procedure.

Nevertheless, the presented results offer an understanding of the cavitation
phenomena in liver tissue near the ultrasonically actuated medical needle.
Such findings could serve as a starting point for designing and developing an
ultrasonic biopsy device in compliance with the safety standards for clinical
applications, and for exploring its potential in other medical applications
involving pathological destruction of tissue.

##  V Conclusion

To conclude, we have studied numerically the dynamics of cavitation bubbles
generated in liver tissue near the tip of an ultrasonically actuated needle.
Experimentally, we have developed a method to capture and quantify the
cavitation activity within thin slices of fresh bovine liver. The main finding
was that cavitation exhibited a TAP dependent behavior, manifesting at TAP $>$
$0.2\text{\,}\mathrm{W}$ and with intensity proportional to the TAP level.
Based on a qualitative comparison, the numerical and the experimental results
presented similarities concerning the cavitation threshold and the spatial
probability of cavitation occurrence around the needle tip. The results are
important since they broaden the understanding of the onset and spatio-
temporal behavior of cavitation near ultrasonically actuated medical needles.
This is especially relevant for ensuring appropriate safety in clinical
scenarios, but also for employing the information in the development of
USeFNAB and other new applications of ultrasonically actuated medical needles.

###### Acknowledgements.

We thank all members of the Medical Ultrasonics Laboratory (MEDUSA) at Aalto
University and Prof. Edward Hæggström at University of Helsinki for
constructive discussions related to the topic. The Academy of Finland is
acknowledged for financial support (grants 314286, 311586 and 335799).

## Author Contributions Statement

All authors contributed to the design of the study, writing or reviewing the
manuscript and have approved the final version of the manuscript. Emanuele
Perra produced all data and conducted the data analyses.

## Conflict of Interest

Heikki J. Nieminen, Kenneth P.H. Pritzker have stock ownership in Swan
Cytologics Inc., Toronto, ON, Canada and are inventors within the patent
application WO2018000102A1. Emanuele Perra and Nick Hayward do not have any
competing interests in relation to this work.

## Data Availability

The datasets are available upon request.

## Code Availability

The codes are available upon request.

## References

  * (1) E. A. Neppiras, “Acoustic cavitation,” Phys. Rep. 61(3), 159–251 (1980) doi: 10.1016/0370-1573(80)90115-5.
  * (2) C. K. Holland and R. E. Apfel, “Thresholds for transient cavitation produced by pulsed ultrasound in a controlled nuclei environment,” J. Acoust. Soc. Am. 88(5), 2059–2069 (1990) doi: 10.1121/1.400102.
  * (3) R. E. Apfel, “Acoustic cavitation,” Methods Exp. Phys. 19(C), 355–411 (1981) doi: 10.1016/S0076-695X(08)60338-5.
  * (4) J. E. Kennedy, “High-intensity focused ultrasound in the treatment of solid tumours,” Nat. Rev. Cancer 5(4), 321–327 (2005) doi: 10.1038/nrc1591.
  * (5) Y. Matsumoto, J. S. Allen, S. Yoshizawa, T. Ikeda, and Y. Kaneko, “Medical ultrasound with microbubbles,” Exp. Therm. Fluid Sci. 29(3), 255–265 (2005) doi: 10.1016/j.expthermflusci.2004.05.008.
  * (6) A. D. Maxwell, C. A. Cain, A. P. Duryea, L. Yuan, H. S. Gurm, and Z. Xu, “Noninvasive Thrombolysis Using Pulsed Ultrasound Cavitation Therapy - Histotripsy,” Ultrasound Med. Biol. 35(12), 1982–1994 (2009) doi: 10.1016/j.ultrasmedbio.2009.07.001.
  * (7) I. Lentacker, I. De Cock, R. Deckers, S. C. De Smedt, and C. T. Moonen, “Understanding ultrasound induced sonoporation: Definitions and underlying mechanisms,” Adv. Drug Deliv. Rev. 72, 49–64 (2014) doi: 10.1016/j.addr.2013.11.008.
  * (8) E. Perra, E. Lampsijärvi, G. Barreto, M. Arif, T. Puranen, E. Hæggström, K. P. H. Pritzker, and H. J. Nieminen, “Ultrasonic actuation of a fine-needle improves biopsy yield,” Sci. Rep. 11(1), 8234 (2021) doi: 10.1038/s41598-021-87303-x.
  * (9) L. Gui, R. Xiao, X. Liao, Y. Kuang, M. Sadiq, Z. Chunnian, S. Cochran, and Z. Huang, “Automatic frequency tracking system for needle actuating device,” in _IEEE Int. Ultrason. Symp. IUS_ , IL, USA (2014), pp. 815–818, doi: 10.1109/ULTSYM.2014.0200.
  * (10) Y. Kuang, A. Hilgers, M. Sadiq, S. Cochran, G. Corner, and Z. Huang, “Modelling and characterisation of a ultrasound-actuated needle for improved visibility in ultrasound-guided regional anaesthesia and tissue biopsy,” Ultrasonics 69, 38–46 (2016) doi: 10.1016/j.ultras.2016.02.018.
  * (11) M. R. Sadiq, S. Cochran, X. Liao, and Z. Huang, “Enhanced US-guided needle intervention through ultrasound actuation of a standard needle,” in _IEEE Int. Ultrason. Symp. IUS_ , IEEE Computer Society (2014), pp. 827–830, doi: 10.1109/ULTSYM.2014.0203.
  * (12) X. Liao, M. Sadiq, Y. Kuang, G. Corner, S. Cochran, and Z. Huang, “Performance optimization of ultrasonic needle actuating device for insertion operation into tissue mimics,” in _IEEE Int. Ultrason. Symp. IUS_ , IEEE, IL, USA (2014), pp. 823–826, doi: 10.1109/ULTSYM.2014.0202.
  * (13) X. Liao, M. Sadiq, G. Corner, S. Cochran, and Z. Huang, “Reduced penetration force through ultrasound activation of a standard needle: An experimental and computational study,” in _IEEE Int. Ultrason. Symp. IUS_ , Prague, Czech Republic (2013), Vol. 1, pp. 1436–1439, doi: 10.1109/ULTSYM.2013.0364.
  * (14) L. Mancia, E. Vlaisavljevich, N. Yousefi, M. Rodriguez, T. J. Ziemlewicz, F. T. Lee, D. Henann, C. Franck, Z. Xu, and E. Johnsen, “Modeling tissue-selective cavitation damage,” Phys. Med. Biol. 64(22), 225001 (2019) doi: 10.1088/1361-6560/ab5010.
  * (15) A. Blana, B. Walter, S. Rogenhofer, and W. F. Wieland, “High-intensity focused ultrasound for the treatment of localized prostate cancer: 5-Year experience,” Urology 63(2), 297–300 (2004) doi: 10.1016/j.urology.2003.09.020.
  * (16) “COMSOL Multiphysics,” Available at www.comsol.com (2014).
  * (17) J. B. Keller and M. Miksis, “Bubble oscillations of large amplitude,” J. Acoust. Soc. Am. 68(2), 628–633 (1980) doi: 10.1121/1.384720.
  * (18) E. Vlaisavljevich, “Investigations of the Cavitation and Damage Thresholds of Histotripsy and Applications in Targeted Tissue Ablation,” Ph.D. thesis, 2015.
  * (19) L. L. Thompson, “A review of finite-element methods for time-harmonic acoustics,” J. Acoust. Soc. Am. 119(3), 1315–1330 (2006) doi: 10.1121/1.2164987.
  * (20) M. F. Hamilton and D. T. Blackstock, _Nonlinear acoustics_ , Vol. 237 (Academic press, San Diego, 1998).
  * (21) M. A. Solovchuk, T. W. Sheu, M. Thiriet, and W. L. Lin, “On a computational study for investigating acoustic streaming and heating during focused ultrasound ablation of liver tumor,” Appl. Therm. Eng. 56(1-2), 62–76 (2013) doi: 10.1016/j.applthermaleng.2013.02.041.
  * (22) R. T. Beyer, “Nonlinear Acoustics,” in _Compr. Biomed. Phys._ , Vol. 2 (Elsevier B.V., 1965), pp. 231–264, doi: 10.1016/B978-0-12-395662-0.50014-X.
  * (23) D. elipe Gaitan, L. A. Crum, C. C. Church, and R. A. Roy, “Sonoluminescence and bubble dynamics for a single, stable, cavitation bubble,” J. Acoust. Soc. Am. 91(6), 3166–3183 (1992) doi: 10.1121/1.402855.
  * (24) R. Gaudron, M. T. Warnez, and E. Johnsen, “Bubble dynamics in a viscoelastic medium with nonlinear elasticity,” J. Fluid Mech. 766(1), 54–75 (2015) https://doi.org/10.1017/jfm.2015.7 doi: 10.1017/jfm.2015.7.
  * (25) L. Mancia, E. Vlaisavljevich, Z. Xu, and E. Johnsen, “Predicting Tissue Susceptibility to Mechanical Cavitation Damage in Therapeutic Ultrasound,” Ultrasound Med. Biol. 43(7), 1421–1440 (2017) doi: 10.1016/j.ultrasmedbio.2017.02.020.
  * (26) J. B. Estrada, M. T. Scimone, H. C. Cramer, L. Mancia, E. Johnsen, and C. Franck, “Microcavitation as a Neuronal Damage Mechanism in an In Vitro Model of Blast Traumatic Brain Injury,” Biophys. J. 112(3), 159a (2017) doi: 10.1016/j.bpj.2016.11.877.
  * (27) J. C. Bamber and C. R. Hill, “Acoustic properties of normal and cancerous human liver-I. Dependence on pathological condition,” Ultrasound Med. Biol. 7(2), 121–133 (1981) doi: 10.1016/0301-5629(81)90001-6.
  * (28) H. Arslan and Y. B. Dolukan, “Determination of the optical properties of bovine liver tissue using integrating sphere technique,” Acad. Perspect. Procedia 1(1), 248–252 (2018) doi: 10.33793/acperpro.01.01.46.
  * (29) Mathworks, “MATLAB,” Available at: www.mathworks.com (2016).
  * (30) Otsu and N., “A threshold selection method from gray-level histograms,” IEEE Trans. Syst. Man Cybern. 9(1), 62–66 (1996).
  * (31) W. Thielicke and R. Sonntag, “Particle Image Velocimetry for MATLAB: Accuracy and enhanced algorithms in PIVlab,” J. Open Res. Softw. 9(1) (2021) doi: 10.5334/jors.334.
  * (32) A. A. Atchley, “The blake threshold of a cavitation nucleus having a radius-dependent surface tension,” J. Acoust. Soc. Am. 85(1), 152–157 (1989) doi: 10.1121/1.397724.
  * (33) J. Chen, L. B. Priddy, R. Prabhu, E. Marin, L. Williams, M. Horstemeyer, and J. Liao, “Traumatic injury: Mechanical response of porcine liver tissue under high strain rate compression testing,” in _Proc. ASME Summer Bioeng. Conf. 2009, SBC2009_ , PART A (2009), pp. 601–602, doi: 10.1115/SBC2009-206532.
  * (34) R. E. Apfel and C. K. Holland, “Gauging the likelihood of cavitation from short-pulse, low-duty cycle diagnostic ultrasound,” Ultrasound Med. Biol. 17(2), 179–185 (1991) doi: 10.1016/0301-5629(91)90125-G.
  * (35) T. Şen, O. Tüfekçioğlu, and Y. Koza, “Mechanical index,” Anadolu Kardiyol. Derg. 15(4), 334–336 (2015) doi: 10.5152/akd.2015.6061.
  * (36) D. M. Hallow, A. D. Mahajan, T. E. McCutchen, and M. R. Prausnitz, “Measurement and correlation of acoustic cavitation with cellular bioeffects,” Ultrasound Med. Biol. 32(7), 1111–1122 (2006) doi: 10.1016/j.ultrasmedbio.2006.03.008.
  * (37) C. C. Church, E. L. Carstensen, W. L. Nyborg, P. L. Carson, L. A. Frizzell, and M. R. Bailey, “The risk of exposure to diagnostic ultrasound in postnatal subjects: Nonthermal mechanisms,” in _J. Ultrasound Med._ , John Wiley and Sons Ltd (2008), Vol. 27, pp. 565–592, doi: 10.7863/jum.2008.27.4.565.
  * (38) G. T. Haar, “Ultrasound bioeffects and safety,” Proc. Inst. Mech. Eng. Part H J. Eng. Med. 224(2), 363–373 (2010) doi: 10.1243/09544119JEIM613.
  * (39) Y. Zhou, “High intensity focused ultrasound in clinical tumor ablation,” World J. Clin. Oncol. 2(1), 8 (2011) doi: 10.5306/wjco.v2.i1.8.
  * (40) S. Yoshizawa, T. Ikeda, A. Ito, R. Ota, S. Takagi, and Y. Matsumoto, “High intensity focused ultrasound lithotripsy with cavitating microbubbles,” Med. Biol. Eng. Comput. 47(8), 851–860 (2009) doi: 10.1007/s11517-009-0471-y.
  * (41) F. K. Chen, A. L. De Castro Abreu, and S. L. Palmer, “Utility of ultrasound in the diagnosis, treatment, and follow-up of prostate cancer: State of the art,” J. Nucl. Med. 57(Suppl 3), 13S–18S (2016) doi: 10.2967/jnumed.116.177196.
  * (42) C. M. Pacella and E. Papini, “Image-guided percutaneous ablation therapies for local recurrences of thyroid tumors,” J. Endocrinol. Invest. 36(1), 61–70 (2013) doi: 10.1007/BF03346744.
  * (43) A. Larghi, M. Rimbaș, G. Rizzatti, C. Carbone, A. Gasbarrini, G. Costamagna, S. Alfieri, and G. Tortora, “Endoscopic ultrasound-guided therapies for pancreatic solid tumors: An overview,” Semin. Oncol. (2021) doi: 10.1053/j.seminoncol.2021.01.004.
  * (44) F. Piccioni, A. Poli, L. C. Templeton, T. W. Templeton, M. Rispoli, L. Vetrugno, D. Santonastaso, and F. Valenza, “Anesthesia for percutaneous radiofrequency tumor ablation (PRFA): A review of current practice and techniques,” Local Reg. Anesth. 12, 127–137 (2019) doi: 10.2147/LRA.S185765.
  * (45) B. George, M. Fink, H. Ermert, S. J. Rupitsch, P. T. Hiltl, and G. Lee, “Investigation of the Inertial Cavitation Activity of Sonosensitive and Biocompatible Nanoparticles for Drug Delivery Applications Employing High Intensity Focused Ultrasound,” Curr. Dir. Biomed. Eng. 5(1), 585–588 (2019) doi: 10.1515/cdbme-2019-0147.
  * (46) G. A. Husseini, M. A. Diaz De La Rosa, E. S. Richardson, D. A. Christensen, and W. G. Pitt, “The role of cavitation in acoustically activated drug delivery,” J. Control. Release 107(2), 253–261 (2005) doi: 10.1016/j.jconrel.2005.06.015.
  * (47) Q. L. Zhou, Z. Y. Chen, Y. X. Wang, F. Yang, Y. Lin, and Y. Y. Liao, “Ultrasound-mediated local drug and gene delivery using nanocarriers,” (2014), doi: 10.1155/2014/963891.
  * (48) G. U. Dachs, G. J. Dougherty, I. J. Stratford, and D. J. Chaplin, “Targeting gene therapy to cancer: A review,” Oncol. Res. 9(6-7), 313–325 (1997).
  * (49) M. Aryal, C. D. Arvanitis, P. M. Alexander, and N. McDannold, “Ultrasound-mediated blood-brain barrier disruption for targeted drug delivery in the central nervous system,” Adv. Drug Deliv. Rev. 72, 94–109 (2014) doi: 10.1016/j.addr.2014.01.008.
  * (50) V. R. Amin, “Ultrasonic attenuation estimation for tissue characterization,” (1989), doi: 10.31274/rtd-180813-8099.
  * (51) K. J. Pahk, M. O. de Andrade, P. Gélat, H. Kim, and N. Saffari, “Mechanical damage induced by the appearance of rectified bubble growth in a viscoelastic medium during boiling histotripsy exposure,” Ultrason. Sonochem. 53, 164–177 (2019) doi: 10.1016/j.ultsonch.2019.01.001.
  * (52) M. Warnez, E. Vlaisavljevich, Z. Xu, and E. Johnsen, “Damage mechanisms for ultrasound-induced cavitation in tissue,” in _AIP Conf. Proc._ (2017), Vol. 1821, doi: 10.1063/1.4977634.
  * (53) M. L. Palmeri, M. H. Wang, J. J. Dahl, K. D. Frinkley, and K. R. Nightingale, “Quantifying Hepatic Shear Modulus In Vivo Using Acoustic Radiation Force,” Ultrasound Med. Biol. 34(4), 546–558 (2008) doi: 10.1016/j.ultrasmedbio.2007.10.009.

Table 1: A table of the general parameters of the numerical model at the ambient temperature of $25\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Properties | Values | References
---|---|---
Ultrasound frequency, $f$ | $33\text{\,}\mathrm{kHz}$ | -
Speed of sound, $c_{\infty}$ | $1575\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ | Amin1989
Density, $\rho_{\infty}$ | $1060\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ | Pahk2019
Nonlinearity parameter, $B/A$ | 7.14 | Bamber1981
Attenuation coefficient, $\alpha$ | $1.39\text{\,}\mathrm{dB}\text{\,}{\mathrm{m}}^{-1}$ | Amin1989
Ambient pressure, $p_{\infty}$ | $101.325\text{\,}\mathrm{kPa}$ | Mancia2019
Polytropic index, $\kappa$ | 1 | Mancia2017
Surface tension, $S$ | $56\text{\,}\mathrm{mN}$ | Warnez2017
Viscosity, $\mu$ | $30\text{\,}\mathrm{mPa}\text{\,}\mathrm{s}$ | Warnez2017
Shear modulus, $G$ | $1.8\text{\,}\mathrm{kPa}$ | Palmeri2008
Initial bubble radius, $R_{0}$ | $500\text{\,}\mathrm{nm}$ | Vlaisavljevich2015
Figure 1:  (a,b) Schematics representing the experimental setup. (c) A custom-
made glass sample holder is used as housing for a thin liver tissue slice.
(d,e) The ultrasonic needle is made to vibrate sideways inside the sample in
the direction of the positive x-axis, while the generated cavitation events
are recorded with a HS camera using a collimated beam of light to produce
back-lit shadowgraph footages. Figure 2:  General description of the
segmentation process of the HS video frames. A cross-correlation based image
registration is applied along the x-axis between (a) the reference frame and
the (b) i-th frame, in order to estimate the displacement $\Delta x_{i}$ of
the needle tip from its reference position. (c) The reference frame is then
translated horizontally by $\Delta x_{i}$, while (d) the input frame is
thresholded with the Otsu method. (e) The binary image showing the translated
reference frame is subtracted to (f) the thresholded input frame, which was
previously closed with a morphological closing operation. (g) The result of
the subtraction is finally filtered with a morphological opening operation, in
order to produce (h) the binary mask for the cavitation activity. Figure 3:
(a1) Simulation of cavitation around the tip of an ultrasonic hypodermic
needle embedded in liver tissue, where the selected time points represent
fractions of an acoustic cycle of duration $T=1/f$. The surrounding medium is
seeded with cavitation nuclei and the needle is actuated at TAP levels of 0.2,
0.3 and $0.8\text{\,}\mathrm{W}$ at a frequency of $33\text{\,}\mathrm{kHz}$.
(a2) Evolution of the driving pressure (MPa), radius (mm) and velocity (m/s)
of a bubble evaluated at position B1. (b1) The probability maps show the
cavitation occurrence probability around the needle tip, when TAP levels of
0.2, 0.3 and $0.8\text{\,}\mathrm{W}$ are delivered to the needle. (b2) The
amount of tissue volume influenced by cavitation activity is estimated by
identifying the regions where the von Mises strain exceeds the ultimate
fractional strain measured for liver (0.38 $\upmu$m/$\upmu$m). The numerical
results suggest that measurable cavitation events are triggered at TAP $>$
$0.2\text{\,}\mathrm{W}$, with greater influence on tissue at higher TAP
levels than at lower TAP levels. Figure 4: (a) Exemplary HS video frames
showing cavitation events taking place at the needle tip at 3 different TAP
levels. (b) The probability maps demonstrate that no cavitation activity is
recorded at $0.2\text{\,}\mathrm{W}$, while it is more frequent at higher TAP
levels. Figure 5:  (a) Ultrasonic sequence adopted during the experiments ($f$
= $33\text{\,}\mathrm{kHz}$, pulse repetition frequency (PRF) =
$55\text{\,}\mathrm{Hz}$, duty cycle (DC) = $50\text{\,}\%$). (b) Time
evolution of the needle-tip displacements, obtained by computing the moving
maximum of the raw data (window size $\sim$ 2 acoustic cycles), and (c)
projected area of the cavitation activity (filtered with moving average,
window size $\sim$10 acoustic cycles). d and e represent the peak-tip
displacement and time integral of cavitation activity calculated for each
individual burst of the ultrasonic sequence, while f represents the same
information evaluated across a time window of $180\text{\,}\mathrm{ms}$. In
the bar charts, the bar height represents the mean of the data set and the
error bar indicates the standard deviation. Figure 6:  (a) Velocity vector
field distribution overlapped to an exemplary HS frame showing a cavitation
event and (b) the velocity magnitude map, when the highest TAP
($0.8\text{\,}\mathrm{W}$) is employed. c and d show the shear and strain rate
distribution in the tissue surrounding the needle tip.


# Kinetic numerical modelling of a magnetopause current sheet with different
energy distribution of particles

A. A. Nechaev<EMAIL_ADDRESS>M. A. Garasev  Vl. V. Kocharovsky
Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod,
603950, Russia

(August 28, 2024)

###### Abstract

We carried out particle-in-cell numerical modelling of various multicomponent
current sheets of the family of analytically found Vlasov–Maxwell equilibria
that allow for an arbitrary particle energy distribution. Depending on the
parameters, the development of small-scale Weibel-type instability inside the
sheet is either possible or excluded. Both cases were studied and compared for
two variants of the particle energy distribution — Maxwellian and Kappa. We
discuss the applicability of the developed model of current sheets and the
results on their small-scale instability to describe physical phenomena in the
vicinity of the magnetopauses of planets and late-type stars.

##  I Introduction

It is well-known that magnetopauses and associated boundary current sheets are
one of the main large-scale structures in the outer magnetosphere of planets
and late-type stars Lotekar  _et al._ (2022); Burch and Phan (2016); Nakamura
_et al._ (2018); Norgren  _et al._ (2018); Louarn  _et al._ (2004); Dunlop
_et al._ (2008); Haaland  _et al._ (2019); Burlaga  _et al._ (2019); Krimigis
_et al._ (2019); Richardson  _et al._ (2019); Shuster  _et al._ (2021);
Allanson  _et al._ (2017); Neukirch  _et al._ (2020). Until now, the
connection between these quasi-static large-scale current structures and
plasma dynamics of small-scale turbulence in the magnetopause region remains
poorly understood. In particular, it is unclear how the possibility of the
generation and existence of quasi-magnetostatic turbulence of the Weibel type
inside it depends on the structure and parameters of the current sheet.
Especially, if the latter has an extended region with a relatively weak
magnetic field and a sufficiently strong anisotropy of the particle velocity
distribution.

In this work, we attempt to answer this question on the basis of particular,
but quite representative examples of two- and four-component current sheets,
found analytically for an arbitrary energy distribution of particles Nechaev
_et al._ (2023); Kocharovsky, Kocharovsky, and Nechaev (2021); Kocharovsky
_et al._ (2019, 2016). For definiteness, we limit ourselves to two such
distributions, Maxwellian and Kappa, and use the resulting analytical self-
consistent stationary solution of the Vlasov–Maxwell equations as the initial
condition in the numerical simulation of the problem of current sheet
stability. The particle-in-cell simulations was performed with the code EPOCH
Arber  _et al._ (2015). Previously such analysis of the development of the
Weibel turbulence in the magnetopause was missing.

The purpose of this work is to elucidate the structure and parameters of
current sheets in which the Weibel-type instability is possible despite the
presence of a self-consistent magnetic field, a low level of anisotropy of the
particle velocity distribution and a limited region of sufficiently dense
nonequilibrium plasma. This paper presents some illustrative examples of
numerical simulation of the evolution of the proposed current sheets both in
the absence and in the presence of the Weibel-type instability.

##  II Analytical estimates

Let us first recall some basic equations used in our models of boundary
current sheets following our work Nechaev  _et al._ (2023).

The distribution function for each plasma component is chosen as an isotropic
function of the absolute value of the momentum $p=|\vec{p}|$, i.e. the
particle energy, multiplied by the Heaviside step function of one of the
projections of the generalized momentum:

$\displaystyle f_{\mathrm{e}}(x,p,p_{y})$ $\displaystyle=$ $\displaystyle
N_{\mathrm{e}1}F_{\mathrm{e}1}(p)\,H\\!\left(-p_{y}+\frac{e}{c}[A_{y}(x)-A_{\mathrm{e}1}]\right)+N_{\mathrm{e}2}F_{\mathrm{e}2}(p)\,H\\!\left(p_{y}-\frac{e}{c}[A_{y}(x)-A_{\mathrm{e}2}]\right)$
(1) $\displaystyle+n_{\mathrm{e}0}(x)F_{\mathrm{e}0}(p),$ $\displaystyle
f_{\mathrm{i}}(x,p,p_{y})$ $\displaystyle=$ $\displaystyle
N_{\mathrm{i}1}F_{\mathrm{i}1}(p)\,H\\!\left(p_{y}+\frac{e}{c}[A_{y}(x)-A_{\mathrm{i}1}]\right)+N_{\mathrm{i}2}F_{\mathrm{i}2}(p)\,H\\!\left(-p_{y}-\frac{e}{c}[A_{y}(x)-A_{\mathrm{i}2}]\right)$
(2) $\displaystyle+n_{\mathrm{i}0}(x)F_{\mathrm{i}0}(p).$

Here the plasma components are denoted by the indices
$\alpha=\mathrm{e},\mathrm{i}$, corresponding to electrons or ions
respectively, and $s=1,2$, which serves to distinguish countercurrents of
particles: $s=1$ determines the current directed along the Cartesian
coordinate axis $y$, $s=2$ is against it. For each component, the constant
$N_{\alpha s}$ and the isotropic factor $F_{\alpha s}(p)$ characterize the
number density and energy distribution of particles on one side of the sheet,
in the region $(-1)^{s+1}x\to+\infty$, where this component is isotropic.
$H(\xi)$ is the Heaviside function such that $H(\xi>0)=1$ and $H(\xi<=0)=0$.
The vector potential $\vec{A}$ defining the magnetic field $B_{z}(x)$ of the
one-dimensional current sheet has a single nontrivial component $A_{y}$, so
that $B_{z}=\mathrm{d}A_{y}/\mathrm{d}x$. Constant shifts $A_{\alpha
s}=\mathrm{const}$ can be chosen different for different plasma components.
For simplicity, ions are assumed to be single-charged, $e$ is the elementary
charge, $c$ is the speed of light in vacuum. Finally, the velocity isotropic
terms with indices $0$ in each function do not contribute to the current and
are introduced in order to guarantee the electrical neutrality of the current
sheet.

An analysis of the general properties of current sheets specified by the
distribution functions (1) and (2) with arbitrary $F_{\alpha s}(p)$ was
carried out in the article Nechaev  _et al._ (2023).

In the numerical simulations presented in Section III we are using only
Maxwellian and Kappa distributions, so let us provide specific formulas for
these two cases.

For the current sheets with Maxwellian distribution of the form

$F_{\alpha s}(p)=\left(2\pi m_{\alpha}T_{\alpha
s}\right)^{-3/2}\exp\left(-\frac{p^{2}}{p_{\alpha s}^{2}}\right)$ (3)

it is easy to obtain the following expressions for the profiles of the current
density $j_{y}(x)$ and the number density $n_{\alpha s}(x)$ of particles by
simply integrating the functions (1) and (2) over the momenta:

$j_{y}(A_{y})=\sum\limits_{\alpha,s}j_{\alpha
s}=\sum\limits_{\alpha,s}\frac{(-1)^{s+1}}{2\pi^{1/2}}eN_{\alpha
s}\frac{p_{\alpha s}}{m_{\alpha}}\exp\left(-a_{\alpha s}^{2}\right),$ (4)
$n_{\alpha}(A_{y})=\sum\limits_{\alpha,s}n_{\alpha
s}=\sum\limits_{s}\frac{1}{2}N_{\alpha
s}\left[1+(-1)^{s+1}\textrm{erf}\left(a_{\alpha s}\right)\right],$ (5)

where we denote the dimensionless vector potential $a_{\alpha
s}(x)=[A_{y}(x)-A_{\alpha s}]\,e/(p_{\alpha s}c)$, effective thermal particle
momentum $p_{\alpha s}=\left(2m_{\alpha}T_{\alpha s}\right)^{1/2}$, and
introduce the effective temperature $T_{\alpha s}=\langle p^{2}\rangle_{\alpha
s}/(3m_{\alpha})$ of the $\alpha s$-component of particles in the region of
its isotropy, at $(-1)^{s+1}A_{y}\to+\infty$. Here
$\langle\dots\rangle_{\alpha s}=\int\\!(\dots)f_{\alpha
s}\,d^{3}\vec{p}\,/n_{\alpha s}$ and
$\textrm{erf}\left(\xi\right)=2\pi^{-1/2}\int_{0}^{\xi}\exp\left(-t^{2}\right)dt$.
As we see, the current density of each component, $j_{\alpha s}$, has the same
sign for all values of $A_{y}$, has a single extremum at $a_{\alpha s}=0$ and
tends to $0$ at $A_{y}\ to\pm\infty$. Number densities $n_{\alpha s}$ are
monotonic functions of their argument.

The parameter $\tau_{\alpha}=1-\left(\langle p_{y}^{2}\rangle-\langle
p_{y}\rangle^{2}\right)/\langle p_{x}^{2}\rangle$, characterizing the
anisotropy of particles of the $\alpha$ variety, in the Maxwellian case is
also easily found Kocharovsky  _et al._ (2019):

$\tau_{\alpha}(A_{y})=\frac{e^{-1}m_{\alpha}n_{\alpha}^{-1}j_{\alpha}^{2}+\sum\limits_{s}p_{\alpha
s}a_{\alpha s}j_{\alpha s}}{e\sum\limits_{s}(-1)^{s+1}N_{\alpha s}T_{\alpha
s}\,\textrm{erf}\left(a_{\alpha s}\right)}\,.$ (6)

For each individual component $\tau_{\alpha s}>0$, at $(-1)^{s+1}a_{\alpha
s}\to+\infty$ its distribution function is isotropic, and at
$(-1)^{s+1}a_{\alpha s}\rightarrow-\infty$, where the particle number density
vanishes, $n_{\alpha s}\to 0$, anisotropy reaches its maximum,
$\tau_{\alpha}\to 1$, since there are only particles with large $y$-momentum
projections.

The magnetic field of the sheet is

$\displaystyle
B_{y}=2\pi^{1/2}\left[2P_{0}-\sum\limits_{\alpha,s}(-1)^{s+1}N_{\alpha
s}T_{\alpha s}\left(\textrm{erf}\,a_{\alpha s}-\textrm{erf}\left.a_{\alpha
s}\right|_{x=0}\right)\right]^{1/2}.$ (7)

To construct current sheets, we will choose the parameter $P_{0}$ such that
the expression in square brackets is positive. Then the magnetic field does
not change sign throughout the current sheet, changing from one constant far
on the left side of the sheet to another constant on the right side. For
convenience, the direction of the $y$ axis is chosen so that this sign is
positive.

Finally, the transition from the vector potential to the coordinate is
determined by the expression

$x(A_{y})=\int\limits_{0}^{A_{y}}\\!\frac{dA^{\prime}}{B_{z}(A^{\prime})},$
(8)

where the limits of integration are chosen so that at the conditional center
of the sheet, at $x=0$, the vector potential $A_{y}=0$. Since the magnetic
field is positive everywhere, the function $x(A_{y})$ is monotonic and
invertible, and from $x\to\pm\infty$ it follows $A_{y}\to\pm\infty$.

The expressions (4), (5), (7) and (8) define a model of an equilibrium current
sheet in collisionless plasma corresponding to the magnetopause, i.e. the
transition between the plasma half-spaces with different number densities and
temperatures of particles, and different values of magnetic field.

Similarly for a current sheet with Kappa energy distributions

$F_{\alpha s}(p)=\frac{M_{\kappa}}{\pi^{3/2}p_{\alpha
s}^{3}}\left(1+\frac{p^{2}}{\left(\kappa-3/2\right)p_{\alpha
s}^{2}}\right)^{-\kappa-1},$ (9)

where $\kappa>3/2$, as before $p_{\alpha s}=\left(2m_{\alpha}T_{\alpha
s}\right)^{1/2}$, and the normalization factor is introduced

$M_{\kappa}=\left(\kappa-\frac{3}{2}\right)^{-3/2}\frac{\Gamma(\kappa+1)}{\Gamma(\kappa-1/2)},$

it is easy to obtain expressions Kocharovsky, Kocharovsky, and Nechaev (2021)
for the current and number density of particles:

$\displaystyle j_{\alpha s}(A_{y})$ $\displaystyle=$
$\displaystyle\frac{(-1)^{s+1}}{2\pi^{1/2}}eN_{\alpha s}\frac{p_{\alpha
s}}{m_{\alpha}}M_{\kappa}\frac{\left(\kappa-3/2\right)^{2}}{\kappa(\kappa-1)}\left(1+\frac{a_{\alpha
s}^{2}}{\kappa-3/2}\right)^{-\kappa+1},$ (10) $n_{\alpha
s}(A_{y})=\frac{N_{\alpha s}}{2}+\frac{(-1)^{s+1}}{\pi^{1/2}}N_{\alpha
s}M_{\kappa}\frac{\kappa-3/2}{\kappa}\,a_{\alpha
s}\times{}_{2}F_{1}\\!\left(\frac{1}{2},\kappa;\frac{3}{2};-\frac{a_{\alpha
s}^{2}}{\kappa-3/2}\right),$ (11)

where ${}_{2}F_{1}\\!\left(\dots\right)$ is the Gaussian hypergeometric
function. The explicit formula for the magnetic field $B_{z}(A_{y})$ is
cumbersome and can be found in the paper Kocharovsky, Kocharovsky, and Nechaev
(2021). In practice, to numerically construct the current sheet, we obtain the
field from the expression

$B_{z}(A_{y})=(8\pi)^{1/2}\left[P_{0}-P_{xx}(A_{y})\right]^{1/2},$ (12)

where, in turn, the $xx$-component of the plasma pressure tensor $P_{xx}$ can
be found from the integral

$P_{xx}(A_{y})=c^{-1}\int\limits_{0}^{A_{y}}\\!j_{y}(A^{\prime})\,\mathrm{d}A^{\prime}.$
(13)

In general, the profiles of the magnetic field, current density and number
density for the Kappa and Maxwellian distributions in the case of plasma
components with the same effective temperatures $T_{\alpha s}$ and number
densities $N_{\alpha s}$ are qualitatively similar, but can differ
significantly quantitatively, primarily in terms of their localization scales
Kocharovsky, Kocharovsky, and Nechaev (2021).

Let us now proceed to the analysis of the properties of the electron Weibel
instability (the growth rate of the ion instability is very small due to their
large mass) using the well-known linear theory for the bi-Maxwellian particle
distribution, which we will use for qualitative order-of-magnitude estimates
for both the case of truncated Maxwellian distributions (1)–(2) with
substitution (3), and for the case of Kappa distributions (9). (Note that in
general case, the properties of instability can significantly depend on the
type of distribution function, see, e.g., Kocharovsky  _et al._ (2016); Silva,
Afeyan, and Silva (2021); Kuznetsov  _et al._ (2022).)

Consider a simple situation where a single anisotropic electron component,
$s=1$, is present. As already mentioned, the quantity $\tau_{\mathrm{e}}$
given by (6) is positive everywhere, and the distribution of electrons in
velocity space has the form of a disk, since
$\tilde{T}_{x}=\tilde{T}_{z}=T\geq\tilde{T}_{y}(x)$, where
$\tilde{T}_{j}=(\langle p_{j}^{2}\rangle-\langle
p_{j}\rangle^{2})/m_{\mathrm{e}}$ and the index $\alpha=\mathrm{e}$ is omitted
for brevity:

$\tilde{T}_{x}=\tilde{T}_{z}\geq\tilde{T}_{y},\quad\vec{B}(x)\parallel Oz$
(14)

As was said, for rough estimates we will assume the electron distribution to
be bi-Maxwellian and use the detailed linear theory Weibel (1959); Davidson
(1989); Vagin and Uryupin (2014); Borodachev  _et al._ (2017), ignoring the
presence of the magnetic field of the current sheet. As is known, for a bi-
Maxwellian distribution, two aperiodic modes are unstable (see, e.g., Vagin
and Uryupin (2014)), and perturbations with wave vectors oriented in the
direction of the smallest dispersion of particle velocities have the highest
growth rates. From those, only the perturbations with a wavelength less than
$\lambda_{\mathrm{e},\min}$ can grow. In the approximation of a small
anisotropy parameter, the perturbations on the wavelength of
$\lambda_{\mathrm{e}}=\sqrt{3}\lambda_{\mathrm{e},\min}$ have the highest
growth rate of $\Gamma$, where

$\lambda_{\mathrm{e},\min}\approx\frac{2\pi
c}{\omega_{\mathrm{p}}}\theta_{\mathrm{e}}^{-1/2},$ (15)
$\Gamma\approx\omega_{\mathrm{p}}\frac{p_{\mathrm{e}}}{{m_{\mathrm{e}}}c}\left(\frac{1+\mathrm{erf}\,a_{\mathrm{e}}}{2}\right)^{1/2}\frac{2}{\pi^{1/2}}\left(\frac{\tau_{\mathrm{e}}}{3}\right)^{3/2},$
(16)

$\omega_{\mathrm{p}}=(4\pi e^{2}N_{\mathrm{e}}/{m_{\mathrm{e}}})^{1/2}$ is the
plasma frequency, and
$\theta_{\mathrm{e}}=(n_{\mathrm{e}}/N_{\mathrm{e}})\,\tau_{\mathrm{e}}(1-\tau_{\mathrm{e}})^{-1}$
is the effective anisotropy that determines the instability properties.

In the simplest case of a single fraction of particles, as seen from Eqs. (5),
(6), the quantity $\theta_{\mathrm{e}}$ is positive only near
$a_{\mathrm{e}}=0$ and almost vanishes already at $|a_{\mathrm{e}}|>2$, being
equivalent to $\propto|a_{\mathrm{e}}|\exp(-a_{\mathrm{e}}^{2})$, as
$a_{\mathrm{e}}\rightarrow\pm\infty$. Its maximum is at
$a_{\mathrm{e}}\approx-0.1$, so we can estimate:
$\max\theta_{\mathrm{e}}\approx\theta_{\mathrm{e}}(0)=(\pi-2)^{-1}$. (Since
$\theta_{\mathrm{e}}<1$, the approximation of plasma dispersion function used
to write Eqs. (15), (16) is adequate.)

Let us now compare the electron gyroradius
$r_{\mathrm{L}\alpha}=(2{m_{\mathrm{e}}}\tilde{T}_{x})^{1/2}c(eB_{z})^{-1}$,
where $\tilde{T}_{x}(x)\equiv T_{\mathrm{e}}$ for the Maxwellian energy
distribution (3), with the half of the minimum unstable wavelength
$\lambda_{\mathrm{e},\min}/2$. Using Eq. (15), we obtain

$\displaystyle\varkappa_{\mathrm{e}}\equiv\left(\frac{2\,r_{\mathrm{L}\alpha}}{\lambda_{\mathrm{e},\min}}\right)^{2}=\frac{\theta_{\mathrm{e}}}{\pi^{2}}\,\frac{8\pi
N_{\mathrm{e}}T_{\mathrm{e}}}{B_{z}^{2}}.$ (17)

At $a_{\mathrm{e}}=0$, where anisotropy $\theta_{\mathrm{e}}$ is nearly
maximal, $B_{z}^{2}(0)=8\pi P_{0}$ and the parameter
$\varkappa_{\mathrm{e}}(0)\leq 2\pi^{-2}\theta_{\mathrm{e}}(0)<0.2$, which
follows from the condition $2P_{0}\geq\max
P_{xx}(A_{y})=N_{\mathrm{e}}T_{\mathrm{e}}$ assumed in Eq. (12) for all
constructed models. Note, that even in the limiting case,
$2P_{0}=N_{\mathrm{e}}T_{\mathrm{e}}$, we have
$\varkappa_{\mathrm{e}}(|a_{\mathrm{e}}|\leq 2)\lesssim 0.9$, although outside
of this region, for $a_{\mathrm{e}}>2$, the field $B_{z}\rightarrow 0$ and
formally $\varkappa_{\mathrm{e}}\rightarrow\infty$. However, even the maximum
growth rate, given by Eq. (16), rapidly goes to zero here, $\Gamma\rightarrow
0$ (and $\lambda_{\mathrm{e}}\rightarrow\infty$), so there is no instability.
Thus, in a simplest model particles are always magnetized in the sheet’s field
regardless of the $N_{\mathrm{e}}$ and $T_{\mathrm{e}}$ values, so the Weibel
instability is inhibited even for the minimal wavelengths.

In the case when there are multiple fractions of particles, but all of them
form currents in the same direction (e.g., for electrons the distribution
$f_{\mathrm{e}}=\sum_{s}N_{\mathrm{e}s}F_{\mathrm{e}s}(p)\,H(-p_{y}+\tilde{A}_{\mathrm{e}s})$
means that all $j_{\mathrm{e}s}>0$), the function $P_{xx}(A_{y})$ is again
monotonic, so we can presume that $\varkappa_{\alpha}\lesssim
2\pi^{-2}\theta_{\alpha}$ in the finite region of instability. If all
$A_{\alpha s}=0$ (or equal to each other), we have Kocharovsky  _et al._
(2019) $\tau_{\alpha}(0)\leq 2/\pi$ leading to
$\theta_{\alpha}(0)\leq(\pi-2)^{-1}$, although it is not clear to the authors,
how close $\max\theta_{\alpha}$ is to $\theta_{\alpha}(0)$ and whether it can
be proven that $\max\varkappa_{\mathrm{e}}<1$ for all sets of plasma
parameters.

When counter-currents are present, the condition
$2P_{0}\geq\max\sum_{s}(-1)^{s+1}N_{\alpha s}T_{\alpha
s}\,\textrm{erf}\,a_{\alpha s}$ does not infer any restrictions on
$\varkappa_{\mathrm{e}}$. This possibility is used in Section III to construct
an unstable current sheet. Moreover, a certain choice of plasma parameters can
provide $\tau_{\alpha}<0$ in some regions of $A_{y}$. Nevertheless, by
increasing the value $P_{0}$ one can assure that $\varkappa_{\mathrm{e}}<1$
for any current sheet of the presented class.

The effect of the Weibel instability inhibiting by a sufficiently strong
magnetic field parallel to the axis of highest temperature, $\vec{B}\parallel
Oz$, can also be obtained from the linear theory. Indeed, as follows from the
dispersion equation for transverse modes in a homogeneous bi-Maxwellian plasma
with temperatures $\tilde{T}_{z}>\tilde{T}_{x}=\tilde{T}_{y}$ in the presence
of an external magnetic field, the inequality $\varkappa_{\mathrm{e}}<1$
approximately coincides with the condition of the collapse of the interval of
unstable wave numbers Emelyanov and Kocharovsky (2023).

##  III Numerical simulations

To verify the analytical estimates, we carried out kinetic numerical particle-
in-cell simulations using the code EPOCH Arber  _et al._ (2015). Calculations
were performed for two models of current sheets to demonstrate both the
stability of simple models without countercurrents and the possibility of
constructing a sheet that is unstable to the Weibel perturbations. In all
simulations, a model ion mass of ${m_{\mathrm{i}}}=18\,{m_{\mathrm{e}}}$ was
used, which reduced the total width of the sheet, proportional to ion
gyroradius, and hence the required size of the simulation domain. In addition,
light ions respond more quickly to the force exerted by the electric field,
which arises due to incomplete compensation of the charge in the central part
of the sheet and disrupts its stationarity. In current sheets with ions of
realistic mass, the influence of this factor will be weaker.

###  III.1 Two-component stable current sheet

Let us consider a sheet formed by one electron and one ion components with
Maxwellian energy distributions. In the isotropic region, the temperatures of
electrons and ions were $T_{\mathrm{{e}}{1}}=150$ eV and
$T_{\mathrm{{i}}{1}}=300$ eV, respectively, the number densities were equal to
$N_{\mathrm{{e1,i1}}{}}=10^{6}$ cm-3. The parameter ${P}_{0}$ was chosen equal
to $1.01\left(\sum_{\alpha,j}N_{\mathrm{{}}{\alpha j}}T_{\mathrm{{}}{\alpha
j}}/2\right)$, so the magnetic field (7) on the right side of the sheet
approached zero. Background electrons and ions, compensating the charge in the
central part of the sheet, had Maxwellian isotropic velocity distributions
with the temperature $T_{\mathrm{{e0,i0}}{}}=10$ eV and number density
profiles $n_{\mathrm{e0}}(x)=\max\\{0,n_{\mathrm{i1}}-n_{\mathrm{e1}}\\}$,
$n_{\mathrm{i0}}(x)=\max\\{0,n_{\mathrm{e1}}-n_{\mathrm{i1}}\\}$. Profiles of
the sheet’s magnetic field, current densities, and particle number densities
at the initial time are shown in Fig. 1.

Figure 1:  Structure of the two-component current sheet at the initial time
moment: (a) black — initial profile of the magnetic field $B_{z0}$ in the
simulation domain $-30<x/r_{\mathrm{L0}}<90$, normalized to its maximum
$B_{z0}(-\infty)$; red — ion current density $j_{\mathrm{i1}}$ in projection
on the axis $y$, normalized to its maximum
$j_{\mathrm{i1}}(0)=(2\pi^{1/2})^{-1}eN_{\mathrm{i1}}p_{\mathrm{i1}}/{m_{\mathrm{i}}}$;
(b) the same projection of electron current density $j_{\mathrm{e1}}$,
normalized to $j_{\mathrm{i1}}(0)$; (c) red — normalized number density
profile of current-carrying ions $n_{\mathrm{i1}}/N_{\mathrm{e1}}$, blue —
same for electrons $n_{\mathrm{e1}}/N_{\mathrm{e1}}$, gree — their difference
$(n_{\mathrm{i1}}-n_{\mathrm{e1}})/N_{\mathrm{e1}}$, proportional to their
charge density; (d) black — effective anisotropy of current-carrying electrons
$\theta_{\mathrm{e}}$, red — effective temperature $\tilde{T}_{y}$ of all
electrons including background, blue — their effective anisotropy, green —
their mutual magnetization parameter $\varkappa_{\mathrm{e}}$.

The gyroradius of the thermal electron in the center of the sheet is equal to
$r_{\mathrm{L0}}\approx 430$ cm, and the width of the sheet does not exceed
$70r_{\mathrm{L0}}$. The effective anisotropy parameter $\theta_{\mathrm{e}}$
(see Fig. 1) has a maximum at $x\approx-0.1\,r_{\mathrm{L0}}$ of the order of
$\theta_{\mathrm{e}}\approx 0.9$, which, according to 15, in the absence of a
magnetic field would correspond to the wavelengths of the fastest growing
Weibel perturbations $\sqrt{3}\lambda_{\mathrm{e},\min}\approx
12\,c/\omega_{\mathrm{p}}$. (Background electrons with the selected density
profile $n_{\mathrm{{e}}{0}}(x)$ reduce the overall anisotropy of electrons on
the left slope of the sheet, where they dominate in number, but do not
significantly change the maximum of the parameter $\theta_{\mathrm{e}}$.)
Taking these estimates into account, the size of the two-dimensional
simulation region were chosen to be equal to
$L_{x}=100\,c/\omega_{\mathrm{p}}\approx 120\,r_{\mathrm{L0}}$ and
$L_{y}=100\,c/\omega_{\mathrm{p}}\approx 15\,\lambda_{\mathrm{e},\min}$ along
the axes $x$ and $y$ respectively. The boundary conditions at $x=-L_{x}/4$ and
$x=3L_{x}/4$ were reflective for particles and absorbing for fields, the
boundaries $y=0$ and $y=L_{y}$ were periodic. The region was divided into
$1250\times 1250$ cells, the current-carrying components of the plasma
comprised approximately $75\times 10^{6}$ macroparticles, and the background
components for $25\times 10^{6}$.

According to analytical estimates from Section II, in a sheet without
countercurrents one should not expect the development of the Weibel
instability. Indeed, as an exact calculation shows, the maximum magnetization
parameter (17) of all electrons in the sheet, $\varkappa_{\mathrm{e}}(x\approx
0.2\,r_{\mathrm{L0}})\approx 0.25$, is significantly less than $1$. If, for
the same plasma parameters, the electrons were not magnetized by the current
sheet’s field, the growth time of the Weibel perturbations (their inverse
maximum growth rate $\Gamma^{-1}$) would be of the order of
$370\,\omega_{\mathrm{p}}^{-1}$. (The value is obtained from the relation (16)
substituting the effective electron temperature $\tilde{T}_{y}\approx
0.4\,T_{\mathrm{e1}}$, see Fig. 1.)

The simulation was carried out until $t=17,500\,\omega_{\mathrm{p}}^{-1}$.
Already at times of the order of $800\,\omega_{\mathrm{p}}^{-1}$ to the left
of the sheet’s center in the distribution of the magnetic field $B_{z}$ and
the number densities of all electrons $n_{\mathrm{e1}}+n_{\mathrm{e0}}$ a
modulation appears along the $y$ axis, reaching maximum at the moment
$t\approx 4900\,\omega_{\mathrm{p}}^{-1}$ shown in Fig. 2. However, this
instability, which develops at the plasma-vacuum boundary, is not associated
with the excitation of currents and, apparently, is not of a Weibel type.
Indeed, in the region $x/r_{\mathrm{L0}}<-1$, where the magnetic field is
high, the number density of particles is insignificant, and the anisotropy of
electrons vanishes, the Weibel instability should be inhibited. Moreover, the
current density components $j_{x}$, $j_{z}$ contain only numerical noise, and
the main component $j_{y}$, which forms the sheet, is not perturbed (Fig.
2(b)). To the right of the boundary with vacuum, the sheet’s magnetic field
$B_{z}$ and the plasma number density remain uniform along the $y$ axis. The
profiles of their dependence on the $x$ coordinate only slowly flatten out,
probably due to influence of the electric field in the central part of the
sheet arising due to thermal movement of the background particles and
violation of charge compensation. By the end of the simulation, the magnetic
field profile $B_{z}(x)$ averaged along the $y$ axis remains monotonic and the
magnitude of its relative change, $|B_{z}(t)/B_{z0}-1|$, where
$B_{z0}=B_{z}(t=0,x)$, does not exceed 30% (Fig. 3(a)). It is noteworthy that
the number density profiles of current-carrying plasma components (Fig.
3(b),(c)) become equal almost everywhere, so that the charge density
associated with them decreases drastically compared to the initial one.

Figure 2:  Instability of plasma-vacuum boundary at the time moment
$t=4900\,\omega_{\mathrm{p}}^{-1}$: (a) two-dimensional distribution of the
magnetic field perturbation $B_{z}(x,y)-B_{z0}(x,y)$, normalized to $\max
B_{z0}$; (b) distribution of the $y$-component of the current density of all
particles $j_{y}(x,y)$, normalized to the maximum of the initial electron
current density
$j_{\mathrm{e1}}(t=0,x=0)=(2\pi^{1/2})^{-1}eN_{\mathrm{e1}}p_{\mathrm{e1}}/{m_{\mathrm{e}}}$;
(c) distribution of the number density of all electrons
$(n_{\mathrm{e1}}+n_{\mathrm{e0}})/N_{\mathrm{e1}}$.  Figure 3:  Structure of
the two-component current sheet at the time moment
$t=17500\,\omega_{\mathrm{p}}^{-1}$: (a) black solid line — magnetic field
profile $B_{z}(x)/B_{z0}(-\infty)$, black dashes — its initial profile
$B_{z0}(x)/B_{z0}(-\infty)$; blue dash-dotted line — their relative difference
$B_{z}/B_{z0}-1$; (b) red solid line — normalized number density profile of
current-carrying ions $n_{\mathrm{i1}}(x)/N_{\mathrm{e1}}$ at the moment
$t=17500\,\omega_{\mathrm{p}}^{-1}$, blue solid — the same for electrons
$n_{\mathrm{e1}}(x)/N_{\mathrm{e1}}$, dashes — their initial profiles; (c)
green — the difference $(n_{\mathrm{i1}}-n_{\mathrm{e1}})/N_{\mathrm{e1}}$,
proportional to the charge density of current-carrying particles; red and blue
— normalized number densities of the background ions,
$n_{\mathrm{i0}}(x)/N_{\mathrm{e1}}$, and electrons,
$n_{\mathrm{e0}}(x)/N_{\mathrm{e1}}$, respectively.

###  III.2 Four-component unstable current sheet

To demonstrate the possibility of the Weibel instability development, we chose
a sheet in which two pairs of electron and ion components with Maxwellian
energy distributions form countercurrents. Electrons and ions forming currents
in the direction of the $y$ axis, in the region of their isotropy to the right
of the sheet, had temperatures of $T_{\mathrm{{e}}{1}}=1500$ eV and
$T_{\mathrm{{i}}{1}}=300$ eV and equal number densities
$N_{\mathrm{{e1,i1}}{}}=10^{6}$ cm-3; those forming countercurrents —
temperatures $T_{\mathrm{{e}}{2}}=4500$ eV and $T_{\mathrm{{i}}{2}}=300$ eV
and equal number densities $N_{\mathrm{{e2,i2}}{}}=N_{\mathrm{{e}}{1}}/3$ to
the left of the sheet. The countercurrents were shifted in space by adding
constant shifts for the second pair of components in the distributions
(2)–(2): $A_{\mathrm{e2}}=2(2{m_{\mathrm{e}}}T_{\mathrm{{e}}{1}})^{1/2}$ and
$A_{\mathrm{i2}}=2(2{m_{\mathrm{i}}}T_{\mathrm{{i}}{1}})^{1/2}$. To ensure the
existence of a region with a weak magnetic field, the parameter $P_{0}$ was
chosen equal to $0.282\left(\sum_{\alpha,j}N_{\mathrm{{}}{\alpha
j}}T_{\mathrm{{}}{\alpha j}}/2\right)$. In order to save computational
resources, background particles were not included in the calculation, i.e. it
was assumed that they had an infinitely large mass and were motionless. The
structure of the sheet at the initial moment is shown in Fig. 4 (a part of the
simulation domain is shown).

Figure 4:  Structure of the four-component current sheet at the initial
moment: (a) black line — magnetic field profile $B_{z}$, normalized to its
maximum $B_{z}(-\infty)$; red and blue — current densities of ions and
electrons respectively, normalized to their maxima; (b) red — normalized
number density of current-carrying ions,
$(n_{\mathrm{i1}}+n_{\mathrm{i2}})/N_{\mathrm{e1}}$, blues — the same for
current-carrying electrons
$(n_{\mathrm{e1}}+n_{\mathrm{e2}})/N_{\mathrm{e1}}$, green — their difference,
proportional to their charge density; (c) blue — effective anisotropy of all
electrons including the background, $\theta_{\mathrm{e}}$, red — their
effective temperature $\tilde{T}_{y}$, green — their mutual magnetization
parameter $\varkappa_{\mathrm{e}}/10$.

The gyroradius of a thermal electron at the center of the sheet is
$r_{\mathrm{L0}}\approx 940$ cm. In the region $1<x/r_{\mathrm{L0}}<7$, where
the magnetic field is suppressed and the magnetization parameter noticeably
exceeds $1$, reaching a maximum of $\varkappa_{\mathrm{e}}(x\approx
4.5\,r_{\mathrm{L0}})\approx 9$, the effective anisotropy is approximately
$\theta_{\mathrm{e}}\approx 0.4$. (Taking into account the background
electrons does not significantly change these values.) Consequently, in the
constructed sheet the Weibel perturbations with characteristic wavelengths
$\sqrt{3}\lambda_{\mathrm{e},\min}\approx 20\,c/\omega_{\mathrm{p}}$ can grow
at times of the order of $\Gamma^{-1}\sim 330\,\ wplinv$. Due to limited
computational resources, a small two-dimensional region with dimensions
$L_{x}=60\,c/\omega_{\mathrm{p}}\approx 34\,r_{\mathrm{L0}}$ and
$L_{y}=96\,c/\omega_{\mathrm{p}}\approx 9\,\lambda_{\mathrm{e},\min}$ was
used, divided by $2500\times 1000$ cells. The boundary conditions at
$x=-0.275\,L_{x}$ and $x=0.725\,L_{x}$ were reflective for particles and
absorbing for fields, the boundaries $y=0$ and $y=L_{y}$ were periodic. The
four current-carrying plasma components comprised approximately $8\times
10^{8}$ macroparticles, and the two background components comprised
approximately $3\times 10^{8}$.

The results of numerical simulation of the sheet’s evolution are shown in Fig.
5. Starting from the moment $t\approx 1500\,\omega_{\mathrm{p}}^{-1}$ in the
two-dimensional distributions of the magnetic field $B_{z}$ and the main
current density component $j_{y}$, modulation appears in the direction of the
$y$ axis, and the current density components $j_{x}$, $j_{z}$ also begin to
grow. In their spatial spectrum along the $y$ axis, a stage of exponential
growth of $6$ lower harmonics is visible, continuing until the moment
$t\approx 3400\,\omega_{\mathrm{p}}^{-1}$, after which their growth becomes
power-law. The calculation continued until the moment $t\approx
7500\,\omega_{\mathrm{p}}^{-1}$.

Figure 5:  Current density of all particles, normalized to
$j_{\mathrm{e1}}(t=0,x=0)=(2\pi^{1/2})^{-1}eN_{\mathrm{e1}}p_{\mathrm{e1}}/{m_{\mathrm{e}}}$,
at the moments of time (a)–(c) $t\approx 3400\,\omega_{\mathrm{p}}^{-1}$,
(d)–(f) $t\approx 7500\,\omega_{\mathrm{p}}^{-1}$. Different projections are
shown: (a), (d) $j_{y}$; (b), (e) $j_{z}$; (c), (f) $j_{x}$;

Although the small transverse size of the region, $L_{y}$, and the short
duration of the simulation did not allow us to see the entire spectrum of
unstable modes and trace their evolution to the decay stage of the Weibel
magnetic field, even this limited simulation indicates:

  * •

the possibility of the Weibel-type instability within the considered class of
multicomponent highly inhomogeneous current sheets;

  * •

the difference between the spectrum of the emerging turbulence and its
nonlinear evolution towards longer wavelengths, and similar spectrum and
evolution in a spatially homogeneous plasma with the same average anisotropy
of the particle velocity distribution;

  * •

the presence of a long stage of power-law, rather than exponential, growth of
individual modes and the root-mean-square level of turbulence until the moment
of its saturation;

  * •

and the possibility of a significant change of the type of the magnetic field
from regular to turbulent in the inner part of the current sheet.

###  III.3 Four-component current sheet with a Kappa energy distribution of
particles

In order to clarify the dependence of the obtained conclusions on the type of
particle energy distribution, we chose a sheet with the same plasma parameters
and parameter ${P}_{0}$ as in subsection III.2, but with a Kappa particle
energy distribution $F_{\alpha s}(p)$, which was sampled using the algorithm
from the work Zenitani and Nakano (2022). The kappa parameter was equal to
$\kappa=3$ for all plasma components. Background electrons and ions had a
temperature of $T_{\mathrm{e0,i0}}=1$ eV and number density profiles
$n_{\mathrm{e0}}(x)=\max\\{0,n_{\mathrm{i1}}-n_{\mathrm{e1}}\\}$,
$n_{\mathrm{i0}}(x)=\max\\{0,n_{\mathrm{e1}}-n_{\mathrm{i1}}\\}$; their motion
was taken into account completely. To speed up the calculations, the number of
macroparticles for all plasma components was taken to be $16$ times less than
in subsection III.2, other simulation parameters were the same. The structure
of the sheet at the initial moment qualitatively resembles the structure of
the Maxwellian sheet (Fig. 4), but electric currents and density gradients
have a larger spatial localization scale, which is due to the longer tail of
the Kappa distribution compared to Maxwellian Kocharovsky, Kocharovsky, and
Nechaev (2021). The calculation continued until
$t=7500\,\omega_{\mathrm{p}}^{-1}$.

On the whole, the simulation of the considered class of current sheets with
the Kappa energy distribution of particles confirms the conclusions 1–4 made
in subsection III.2 and shows that the type of particle energy distribution,
provided it is sufficiently smooth and qualitatively similar to a Maxwellian,
does not significantly change, in comparison with the latter, either the
global large-scale structure of the current sheet, or the general properties
of the Weibel-type instability in its inner part with a weak regular magnetic
field. Therefore, the possibility of long-term existence of small-scale quasi-
magnetostatic turbulence in distributed multicomponent current sheets seems to
be quite universal.

##  IV Conclusion

In the presented work we investigated the development of small-scale
instability of the Weibel type of two- and four-component current sheets,
allowing an analytical description of a distributed magnetopause, which embeds
a considerably wide region, where magnetic field is weak and particles have
significant velocity anisotropy. Based on numerical particle-in-cell
simulations of the evolution of these current sheets with the EPOCH code, it
is shown that, depending on the parameters of the sheet, the Weibel-type
instability can be inhibited or effectively develop, eliminating or ensuring,
respectively, the presence of quasi-magnetostatic turbulence inside the
current sheet. The properties of the Weibel-type instability inside the
current sheet also depends on the type of particle energy distribution, which
we verified by comparing the results of modelling similar sheets with
Maxwellian and Kappa distributions.

The results of numerical simulations made it possible to identify the features
of the emergence and nonlinear evolution of quasi-magnetostatic Weibel
turbulence in the presence of countercurrents of particles of the same sort
and intersecting currents of electrons and ions, as well as in the presence of
a non-uniform self-consistent magnetic field and a non-uniform distribution of
nonequilibrium plasma with an anisotropic particle velocity distribution. It
was established that the presence of gradients of local properties of
magnetoactive plasma and its inhomogeneity in general significantly change the
spectrum of emerging turbulence and the nature of its gradual evolution from
short to long wavelengths in comparison with those previously studied for the
case of spatially homogeneous magnetoactive plasma. It is shown that in the
inner part of the current sheet, a growing small-scale turbulent field can
significantly change the regular field of the current sheet, but cannot
significantly affect the large-scale structure of the entire current
configuration.

The proposed current sheet models and the established properties of their
small-scale Weibel-type instability, leading to the formation of quasi-
magnetostatic turbulence, are of interest for describing physical phenomena in
the vicinity of the outer parts of planetary magnetospheres exposed to stellar
winds and the coronas of late-type stars.

###### Acknowledgements.

The work was supported by the Russian Science Foundation grant No.
20-12-00268. The computer time for simulations were provided by the Joint
Supercomputer Center of the Russian Academy of Sciences.

## References

  * Lotekar  _et al._ (2022) A. B. Lotekar, I. Y. Vasko, T. Phan, S. D. Bale, T. A. Bowen, J. Halekas, A. V. Artemyev, Y. V. Khotyaintsev, and F. S. Mozer, “Kinetic-scale current sheets in near-Sun solar wind: Properties, scale-dependent features and reconnection onset,” The Astrophysical Journal 929, 58 (2022).
  * Burch and Phan (2016) J. L. Burch and T. D. Phan, “Magnetic reconnection at the dayside magnetopause: Advances with MMS,” Geophys. Res. Lett. 43, 8327––8338 (2016).
  * Nakamura  _et al._ (2018) R. Nakamura, A. Varsani, K. J. Genestreti, O. L. Contel, T. Nakamura, W. Baumjohann, T. Nagai, A. Artemyev, J. Birn, V. A. Sergeev, S. Apatenkov, R. E. Ergun, S. A. Fuselier, D. J. Gershman, B. J. Giles, Y. V. Khotyaintsev, P.-A. Lindqvist, W. Magnes, B. Mauk, A. Petrukovich, C. T. Russell, J. Stawarz, R. J. Strangeway, B. Anderson, J. L. Burch, K. R. Bromund, I. Cohen, D. Fischer, A. Jaynes, L. Kepko, G. Le, F. Plaschke, G. Reeves, H. J. Singer, J. A. Slavin, R. B. Torbert, and D. L. Turner, “Multiscale currents observed by MMS in the flow braking region,” Journal of Geophysical Research: Space Physics 123, 1260–1278 (2018).
  * Norgren  _et al._ (2018) C. Norgren, D. B. Graham, Y. V. Khotyaintsev, M. André, A. Vaivads, M. Hesse, E. Eriksson, P.-A. Lindqvist, B. Lavraud, J. Burch, S. Fuselier, W. Magnes, D. J. Gershman, and C. T. Russell, “Electron reconnection in the magnetopause current layer,” Journal of Geophysical Research: Space Physics 123, 9222–9238 (2018).
  * Louarn  _et al._ (2004) P. Louarn, A. Fedorov, E. Budnik, G. Fruit, J. A. Sauvaud, C. C. Harvey, I. Dandouras, H. Rème, M. C. Dunlop, and A. Balogh, “Cluster observations of complex 3d magnetic structures at the magnetopause,” Geophysical Research Letters 31, L19805 (2004).
  * Dunlop  _et al._ (2008) M. W. Dunlop, M. G. G. T. Taylor, Y. V. Bogdanova, C. Shen, F. Pitout, Z. Pu, J. A. Davies, Q.-H. Zhang, J. Wang, B. Lavraud, A. N. Fazakerley, A. Walsh, C. J. Owen, H. Laakso, Q.-G. Zong, Z.-X. Liu, C. P. Escoubet, C. M. Carr, and H. Rème, “Electron structure of the magnetopause boundary layer: Cluster/Double Star observations,” Journal of Geophysical Research: Space Physics 113, A07S19 (2008).
  * Haaland  _et al._ (2019) S. Haaland, A. Runov, A. Artemyev, and V. Angelopoulos, “Characteristics of the flank magnetopause: THEMIS observations,” Journal of Geophysical Research: Space Physics 124, 3421–3435 (2019).
  * Burlaga  _et al._ (2019) L. F. Burlaga, N. F. Ness, D. B. Berdichevsky, J. Park, L. K. Jian, A. Szabo, E. C. Stone, and J. D. Richardson, “Magnetic field and particle measurements made by Voyager 2 at and near the heliopause,” Nature Astronomy 3, 1007–1012 (2019).
  * Krimigis  _et al._ (2019) S. M. Krimigis, R. B. Decker, E. C. Roelof, M. E. Hill, C. O. Bostrom, K. Dialynas, G. Gloeckler, D. C. Hamilton, E. P. Keath, and L. J. Lanzerotti, “Energetic charged particle measurements from Voyager 2 at the heliopause and beyond,” Nature Astronomy 3, 997–1006 (2019).
  * Richardson  _et al._ (2019) J. D. Richardson, J. W. Belcher, P. Garcia-Galindo, and L. F. Burlaga, “Voyager 2 plasma observations of the heliopause and interstellar medium,” Nature Astronomy 3, 1019–1023 (2019).
  * Shuster  _et al._ (2021) J. R. Shuster, D. J. Gershman, J. C. Dorelli, B. L. Giles, S. Wang, N. Bessho, L.-J. Chen, P. A. Cassak, S. J. Schwartz, R. E. Denton, V. M. Uritsky, W. R. Paterson, C. Schiff, A. F. Viñas, J. Ng, L. A. Avanov, D. E. da Silva, and R. B. Torbert, “Structures in the terms of the Vlasov equation observed at Earth’s magnetopause,” Nature Physics 17, 1056–1065 (2021).
  * Allanson  _et al._ (2017) O. Allanson, F. Wilson, T. Neukirch, Y.-H. Liu, and J. D. B. Hodgson, “Exact Vlasov-Maxwell equilibria for asymmetric current sheets,” Geophysical Research Letters 44, 8685–8695 (2017).
  * Neukirch  _et al._ (2020) T. Neukirch, I. Y. Vasko, A. V. Artemyev, and O. Allanson, “Kinetic models of tangential discontinuities in the solar wind,” The Astrophysical Journal 891, 86 (2020).
  * Nechaev  _et al._ (2023) A. A. Nechaev, V. V. Kocharovsky, V. V. Kocharovsky, and M. A. Garasev, “Multicomponent current sheet of the magnetopause with an arbitrary energy distribution of particles,” JETP Letters 117, 214–221 (2023).
  * Kocharovsky, Kocharovsky, and Nechaev (2021) V. V. Kocharovsky, V. V. Kocharovsky, and A. A. Nechaev, “Analytical model of a magnetopause in a multicomponent collisionless plasma with a Kappa energy distribution of particles,” Doklady Physics 66, 9–13 (2021).
  * Kocharovsky  _et al._ (2019) V. V. Kocharovsky, V. V. Kocharovsky, V. Y. Martyanov, and A. A. Nechaev, “An analytical model for the current structure of the magnetosheath boundary in a collisionless plasma,” Astronomy Letters 45, 551–564 (2019).
  * Kocharovsky  _et al._ (2016) V. V. Kocharovsky, V. V. Kocharovsky, V. Y. Martyanov, and S. V. Tarasov, “Analytical theory of self-consistent current structures in a collisionless plasma,” Physics-Uspekhi 59, 1165–1210 (2016).
  * Arber  _et al._ (2015) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Physics and Controlled Fusion 57, 113001 (2015).
  * Silva, Afeyan, and Silva (2021) T. Silva, B. Afeyan, and L. O. Silva, “Weibel instability beyond bi-Maxwellian anisotropy,” Physical Review E 104, 035201 (2021).
  * Kuznetsov  _et al._ (2022) A. A. Kuznetsov, V. V. Kocharovskii, V. V. Kocharovskii, A. A. Nechaev, and M. A. Garasev, “Saturating magnetic field of Weibel instability in plasmas with bi-Maxwellian and bikappa particle distributions,” Plasma Phys. Rep. 48, 973–982 (2022).
  * Weibel (1959) E. S. Weibel, “Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution,” Physical Review Letters 2, 83–84 (1959).
  * Davidson (1989) R. C. Davidson, “Kinetic waves and instabilities in a uniform plasma,” in  _Basic Plasma Physics: Selected chapters from the Handbook of Plasma Physics, vol. 1 and 2_ , edited by A. A. Galeev and R. N. Sudan (North-Holland Publishing Company, Amsterdam, 1989) p. 229.
  * Vagin and Uryupin (2014) K. Y. Vagin and S. A. Uryupin, “On the growth rate of aperiodic instability in plasma with an anisotropic bi-Maxwellian electron velocity distribution,” Plasma Physics Reports 40, 393–403 (2014).
  * Borodachev  _et al._ (2017) L. V. Borodachev, M. A. Garasev, D. O. Kolomiets, V. V. Kocharovsky, V. Y. Martyanov, and A. A. Nechaev, “Dynamics of a self-consistent magnetic field and diffusive scattering of ions in a plasma with strong thermal anisotropy,” Radiophysics and Quantum Electronics 59, 991–999 (2017).
  * Emelyanov and Kocharovsky (2023) N. A. Emelyanov and V. V. Kocharovsky, “The Weibel instability in an external magnetic field: Analytical estimates,” Radiophys. Quantum El. 66 (2023), (In press.).
  * Zenitani and Nakano (2022) S. Zenitani and S. Nakano, “Loading a relativistic Kappa distribution in particle simulations,” Physics of Plasmas 29, 113904 (2022).