# Gradient descent optimization of acoustic holograms for transcranial focused
ultrasound

Ahmed Sallam  Ceren Cengiz  Mihir Pewekar  Eric Hoffmann  Department of
Mechanical Engineering, Virginia Tech, Blacksburg, VA, $24061$, USA  Wynn
Legon  Fralin Biomedical Research Institute, Virginia Tech, Roanoke, VA,
$24016$, USA  Eli Vlaisavljevich  Department of Biomedical Engineering and
Mechanics, Virginia Tech, Blacksburg, VA, $24061$, USA  Shima Shahab
<EMAIL_ADDRESS>Department of Mechanical Engineering, Virginia Tech,
Blacksburg, VA, $24061$, USA

###### Abstract

Acoustic holographic lenses, also known as acoustic holograms, can change the
phase of a transmitted wavefront in order to shape and construct complex
ultrasound pressure fields, often for focusing the acoustic energy on a target
region. These lenses have been proposed for transcranial focused ultrasound
(tFUS) to create diffraction-limited focal zones that target specific brain
regions while compensating for skull aberration. Holograms for tFUS are
currently designed using time-reversal approaches in full-wave time-domain
numerical simulations. However, such simulations need time-consuming
computations, which severely limits the adoption of iterative optimization
strategies. Furthermore, in the time-reversal method, the number and
distribution of virtual sources can significantly influence the final sound
field. Because of the computational constraints, predicting these effects and
determining the optimal arrangement is challenging. This study introduces an
efficient method for designing acoustic holograms using a volumetric
holographic technique to generate focused fields inside the skull. The
proposed method combines a modified mixed-domain method for ultrasonic
propagation with a gradient descent iterative optimization algorithm. This
approach enables substantially faster holographic computation than previously
reported techniques. The iterative process uses explicitly defined loss
functions to bias the ultrasound field’s optimization parameters to specific
desired characteristics, such as axial resolution, transversal resolution,
coverage, and focal region uniformity, while eliminating the uncertainty
associated with virtual sources in time-reversal techniques. Numerical studies
are conducted on four brain structures: the anterior insula, hippocampus,
caudate nucleus, and amygdala. The findings are further validated in
underwater experiments with a realistic 3D-printed skull phantom.

††preprint: AIP/123-QED

##  I Introduction

Transcranial focused ultrasound (tFUS) is a transformative technology that
employs concentrated ultrasound energy to treat various neurological
conditions. Meng, Hynynen, and Lipsman (2021); Darmani  _et al._ (2022). Non-
ionizing acoustic beams can penetrate through the skull and precisely target
specific structures deep within the brain. The energy of the ultrasound beam
elicits thermal and mechanical effects in brain tissue. Those effects are
utilized for a variety of therapeutic treatments and procedures. At lower
intensities, tFUS can accurately transmit energy to specific regions in the
brain and induce local modulation of neural activity Legon  _et al._ (2014);
Lee  _et al._ (2016); Mueller  _et al._ (2014); Legon  _et al._ (2018a).
Furthermore, tFUS has been demonstrated to locally disrupt the blood-brain
barrier (BBB) Deffieux and Konofagou (2010); Legon  _et al._ (2018b, c);
Strohman  _et al._ (2024), allowing for improved delivery of therapeutic
drugs. High-intensity ultrasound is also leveraged in tFUS to treat various
diseases and disorders Meng, Hynynen, and Lipsman (2021), for which ultrasound
energy can be used to thermally ablate diseased tissue Lipsman  _et al._
(2013), and high pressures can be used for histotripsy Xu  _et al._ (2021);
Kim  _et al._ (2014); Sukovich  _et al._ (2018), in which targeted tissues are
mechanically fractionated.

Focusing the ultrasound beam in tFUS is challenging due to the highly stiff
skull. The high acoustic impedance, irregular geometry, and nonuniform
thickness of the skull bones strongly scatter and attenuate the propagating
ultrasound beam. This phenomenon is known as skull-induced aberrations
Kyriakou  _et al._ (2014); Almquist, Parker, and Christensen (2016), and it
causes severe defocusing and phase distortion of the beam, resulting in
inefficient and out-of-target focusing. To address this issue, 3D-printed
acoustic holographic lenses, known as acoustic holograms, have been recently
proposed for tFUS and compensation of skull aberrations Maimbourg  _et al._
(2018); Kook  _et al._ (2023); Hu  _et al._ (2022); Jiménez-Gambín  _et al._
(2021); He  _et al._ (2021); Jiménez-Gambín  _et al._ (2019); Jiménez-Gambín,
Jiménez, and Camarena (2020). Acoustic holograms act like phase plates,
modulating the phase of the wavefront with unique thickness maps. The
thickness map enables the reconstruction of complicated acoustic fields with a
desired spatial distribution of pressure or phase using a single-element
acoustic transducer. These lenses have been shown to modulate reflected and
transmitted wavefronts Sallam  _et al._ (2021); Sallam, Meesala, and Shahab
(2021); Sallam and Shahab (2022); Melde  _et al._ (2016), forming multi-focal
spots or arbitrary complex patterns for applications such as contactless
ultrasound energy transfer Bakhtiari-Nejad  _et al._ (2018), particle
manipulation Melde  _et al._ (2016, 2023), controlling cavitation zones Kim
_et al._ (2021), and hyperthermia Andrés  _et al._ (2022a). These intricate
modulation and manipulation capabilities make them ideal for tFUS
applications. The acoustic lens compensates for skull aberrations and creates
highly localized pressure fields that direct the ultrasound beam to specific
brain areas within the skull. Acoustic holograms provide an order of magnitude
higher resolution and acoustic pixel density than phased array transducers at
a substantially lower cost and simpler fabrication procedure Melde  _et al._
(2016). They have been tested successfully in bench-top setups using
3D-printed phantoms and ex-vivo human skulls Jiménez-Gambín  _et al._ (2019);
Jiménez-Gambín, Jiménez, and Camarena (2020), as well as in-vivo mouse models
Jiménez-Gambín  _et al._ (2021); He  _et al._ (2021) for neuromodulation and
BBB opening procedures.

Time-reversal numerical simulations and direct error minimization techniques
are currently employed to generate an appropriate thickness map that accounts
for skull aberrations and directs the ultrasonic field to the target structure
in the brain. In the numerical domain, monopole virtual sources are introduced
at the target regions of the skull. The ultrasound is then transmitted from
the simulated sources to the hologram plane outside of the skull, where the
acoustic signal is recorded. The extracted signal is subsequently processed to
determine the phase and amplitude distributions in the hologram plane. The
phase and amplitude information are then used to construct the lens thickness
map. Ultrasound propagation in strongly heterogeneous media is often carried
out using full-wave numerical simulations, namely the pseudospectral time-
domain method implemented in the open-source software k-Wave Treeby and Cox
(2010). However, full-wave numerical simulations are computationally
expensive, lasting up to 20 hours when using CPU parallel processing Jiménez-
Gambín  _et al._ (2019). The long computation time also prevents the use of
iterative optimization techniques, which could help improve focusing quality
and allow the implementation of explicitly specified loss functions to bias
the ultrasound field’s optimization parameters to specific desired
characteristics, such as axial resolution, transversal resolution, coverage,
and focal region uniformity. Furthermore, iterative algorithms take into
account the transducer’s experimental amplitude output, improving the accuracy
of the thickness map and eliminating the need for error diffusion schemes.
Acoustic holography is a technique that allows for the reconstruction of 3-D
acoustic fields by projecting measured or simulated 2-D data from a
measurement plane to the remainder of the volume Sallam and Shahab (2023).
Because of its computing efficiency, holography is used in iterative
optimization techniques such as the Gerchberg-Saxton (GS) and iterative
angular spectrum algorithms to compute and design the thickness map of
holographic lenses Melde  _et al._ (2016, 2023). In those techniques, the
wavefront is iteratively propagated forward and backward between the hologram
and target planes until a satisfying boundary condition is obtained. However,
such techniques are intended for freely propagating media and optimizing
planar pressure patterns. This paper introduces an analogous iterative
optimization approach that takes into account strongly heterogeneous media and
aberrating layers, including the effects of reflection, transmission,
scattering, and attenuation. Furthermore, the proposed method optimizes for
volumetric targets, such as directing sound to specific volumetric components
in the brain. The modified mixed domain method (MMDM) Gu and Jing (2020, 2018,
2021) is used to simulate the forward and backward propagation inside the
strongly heterogeneous media. In addition, a non-convex gradient descent
optimization algorithm is used to find the optimal phase for targeting a
specific volume. The proposed technique is substantially more efficient,
allowing for the creation of a hologram in less than 50 minutes while also
allowing for the explicit defining of loss functions. The proposed algorithm
is evaluated numerically by considering four test cases that target different
brain structures: the anterior insula, hippocampus, caudate nucleus, and
amygdala. Each target resides at a different depth in the brain and has a
distinct complex volume and geometric structure in 3D space. The holographic
lenses obtained for each target are then 3D printed and tested in underwater
trials and on a 3D-printed skull phantom.

##  II Methods

###  II.1 Acoustic propagation

The MMDM Gu and Jing (2021, 2020, 2018) is implemented to execute
computationally efficient forward and backward propagations through the
strongly heterogeneous domain, which includes the skull and water. This
methodology is based on the generalized Westervelt equation and can be thought
of as an expanded form of the angular spectrum method. The acoustic field is
calculated by marching in the spatial domain along the axis that is normal to
the source plane.

Using the generalized Westervelt equation as a starting point:

$\rho\nabla\cdot\left(\frac{1}{\rho}\nabla
p\right)-\frac{1}{c^{2}}\frac{\partial^{2}p}{\partial
t^{2}}+\frac{\delta}{c^{4}}\frac{\partial^{3}p}{\partial
t^{3}}+\frac{\beta}{\rho c^{4}}\frac{\partial^{2}p^{2}}{\partial
t^{2}}=\gamma\frac{\partial p}{\partial t}$ (1)

where $p$ is the acoustic pressure, $\rho$ is the ambient density, $c$ is the
speed of sound, $\delta$ is the sound diffusivity and $\beta$ is the
nonlinearity coefficient. The frequency-independent absorption term
$\gamma\partial p/\partial t$ accounts for the absorption layer close to the
computational boundaries. Applying the normalized wavefield $f=p/\sqrt{\rho}$,
assuming homogeneous density, and performing the Fourier transform with
respect to $X$, $Y$, and $t$ yields:

$\displaystyle\frac{\partial^{2}}{\partial Z^{2}}\tilde{P}+K^{2}\tilde{P}=$
(2) $\displaystyle
F_{XY}\left\\{\left[-\frac{\omega^{2}}{c_{0}^{2}}\left(\frac{c_{0}^{2}}{c^{2}}-1\right)+\frac{i\delta\omega^{3}}{c^{4}}+i\omega\gamma\right]F_{t}(p)\right\\}$
$\displaystyle\quad+F_{XY}\left(\frac{\beta\omega^{2}}{\rho
c^{4}}F_{t}\left(p^{2}\right)\right)$

where $\tilde{P}$ is the Fourier transform (FT) of $p,F_{XY}$ is the FT
operator in $X$ \- and $Y$-dimensions, $F_{t}$ is the FT operator in the time
domain, $c_{0}$ is the background sound speed and
$K^{2}=\omega^{2}/c_{0}^{2}-k_{X}^{2}-k_{Y}^{2}$, $k_{X}$ and $k_{Y}$ are the
wavenumbers in $X$ \- and $Y$-dimensions, respectively. An implicit, one-way
propagation can be derived from the 1-D Green’s function in the form of an
integral equation, such that

$\tilde{P}(Z)=\tilde{P}(0)e^{iKZ}+\frac{e^{iKZ}}{2iK}\int_{0}^{Z}e^{-iKZ^{\prime}}M\left(p\left(Z^{\prime}\right)\right)dZ^{\prime}$
(3)

where,

$\displaystyle M(p)=$ (4) $\displaystyle
F_{XY}\left\\{\left[-\frac{\omega^{2}}{c_{0}^{2}}\left(\frac{c_{0}^{2}}{c^{2}}-1\right)+\frac{i\delta\omega^{3}}{c^{4}}+i\omega\gamma\right]F_{t}(p)\right\\}$
$\displaystyle+F_{XY}\left(\frac{\beta\omega^{2}}{\rho
c^{4}}F_{t}\left(p^{2}\right)\right)$

Equation (3) is solved by using a Simpson-like rule.

The pressure distribution on the source plane or any arbitrary plane is
reconstructed using the backward projection. By changing $Z$ in (3) to $-Z$,
the backward projection is obtained as:

$\tilde{P}(-Z)=\tilde{P}(0)e^{-iKZ}-\frac{e^{-iKZ}}{2iK}\int_{0}^{-Z}e^{iKZ^{\prime}}M\left(p\left(-Z^{\prime}\right)\right)dZ^{\prime}$
(5)

The Kramers-Kronig dispersion relation is applied directly by replacing the
speed of sound $c$ with $c_{p}$ and $c_{p}=\left(1/\hat{c}+\alpha_{0}\tan(\pi
y/2)\omega^{y-1}\right)^{-1}$, where $\hat{c}$ is the sound speed at zero
frequency, $y$ is the power-law exponent, $\alpha_{0}$ is the absorption in
nepers per megahertz-y per meter and $\alpha_{o}=\alpha_{NP}\omega^{-y}$.
Using this model, wave effects such as full-wave diffraction, attenuation, and
dispersion are considered. However, this model is only accurate for weakly
heterogeneous media Gu and Jing (2018, 2020). To have a more general model
that could be applied to strongly heterogeneous media, phase correction and
transmission compensation are introduced.

Assuming 1D propagation through an inhomogeneous media with a varied speed of
sound distribution, the following phase correction term is added:

$\displaystyle\left(P_{z+\Delta z}\right)_{\text{corr}}=$
$\displaystyle\left(P_{z+\Delta z}\right)$ (6)
$\displaystyle+P_{z}e^{iK^{\prime}\Delta
z}\left(1-e^{i\left(K-K^{\prime}\right)\Delta
z}\frac{1+\frac{M_{z}}{4iK}\Delta z}{1-\frac{M_{z+\Delta z}}{4iK}\Delta
z}\right)$

Where $P_{z}$ is the wave pressure at $z$, $(P_{z+\Delta z})_{\text{corr}}$ is
the phase corrected pressure, $K^{\prime}$ is the wave number, and
$K^{\prime}=\omega/c_{2}$ with a frequency of $\omega$.

Moreover, the original formulation does not take into account the transmission
coefficient caused by variations in sound speed. In order to handle the
density and speed of sound heterogeneities, amplitude correction is employed.
The compensation term:

$\displaystyle T(x,y,z)$ (7) $\displaystyle\quad=\frac{2\rho(x,y,z+\Delta
z)c(x,y,z+\Delta z)}{\rho(x,y,z)c(x,y,z)+\rho(x,y,z+\Delta z)c(x,y,z+\Delta
z)}$

where $c(x,y,z)$ and $\rho(x,y,z)$ are the speed of sound and density at plane
$z$, respectively, and $c(x,y,z+\Delta z)$ and $\rho(x,y,z+\Delta z)$ are the
speed of sound and density at plane $z+\Delta z$, respectively.

Finally, reflections are further added to the formulation by using the
following equation:

$P_{\text{reflection }}=P_{\text{incident }}(T-1)\text{, }$ (8)

Here, $P_{\text{incident }}$ is the incident wave used for calculating the
reflected wave. The second-order reflection field is formed by propagating the
resulting $P_{\text{reflection }}$ in the forward direction. This process is
repeated until the intended maximum order of reflection is achieved. This
makes it possible to model scattering from multilayer interfaces and
distributed scattering events. Superposing each solution yields the final wave
field.

###  II.2 Holograms computation

Computer-generated holography (CGH) can be implemented to construct desired 3D
pressure patterns by altering the wavefront of a coherent acoustic beam Melde
_et al._ (2023). CGH determines the spatial phase and amplitude distribution
of an incident wavefront, which, upon diffraction, forms the target volume.
Computing the 2D phase mask that yields certain intensity patterns is
generally a non-trivial and ill-posed task Zhang  _et al._ (2017). Most three-
dimensional intensity patterns are not physically viable. The wave field must
diffract through space in a way that is consistent with Helmholtz’s equations
and maintains intensity conservation at each depth plane. Furthermore,
resolution is limited due to the finite numerical aperture (NA). The phase-
only hologram constraint also compounds these limitations further. As a
result, any required 3D pressure pattern cannot be precisely reproduced. The
optimal thickness map is determined by optimizing the output acoustic field
under experimental limitations such as incident beam profile, source NA, and
wave shaping modality. Following prior efforts in acoustics and optics, we use
a gradient descent optimization approach in the present investigation Zhang
_et al._ (2017); Melde  _et al._ (2023).

The defined loss function greatly influences the outcomes of 3D CGH, where
there is no exact solution but numerous approximate solutions. Our algorithm
searches for optimal holograms using custom cost functions that users can
select to best suit a certain application. The optimization aims to obtain a
phase distribution $\phi$ at the hologram plane that minimizes a specified
error function $L[P_{3D}(\phi),T]$ between the generated acoustic field
$P_{3D}(\phi)$ in the volumetric domain and the volumetric target pattern T,
with the constraint that the pressure amplitude at the transducer plane must
always equal to the experimental pressure $|p0|=A_{0}$. The numerical solution
of this optimization problem necessitates calculating the gradient of the loss
function with respect to the decision variables. By implementing Wirtinger
subdifferential calculus and following the derivation in Melde  _et al._
(2023), the gradient of the loss function $L$ is calculated as:

$\displaystyle\nabla_{p}L=2A^{*}\left(\left.\frac{A\bm{p}}{|A\bm{p}|}\circ\nabla_{x}\ell(\bm{x},\bm{T})\right|_{x=|A\bm{p}|}\right)$
(9)

Here $A$ is a linear parameter representing the propagation of the acoustic
field at the transducer plan to the rest of the volume, such that
$A\bm{p}=P_{3D}$. Correspondingly, $A^{*}$ represents the inverse propagation
from each plane in the volume to the transducer plane such that
$A^{*}P_{3D}=\bm{p}$. Note that the propagation is implemented using the MMDM.
Here $\circ$ denotes the element-wise product. Assuming square loss

$L(\bm{x},\bm{T})=\frac{1}{2}\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\left(x_{ij}-T_{ij}\right)^{2}$
(10)

Substituting equation (10) into (9) yields:

$\nabla_{p}L=2A^{*}\left(A\bm{p}-\bm{g}\circ\frac{A\bm{p}}{|A\bm{p}|}\right)$
(11)

Thus the iterative gradient descent algorithm update rule becomes

$p^{k+1}=\mathcal{P}_{U}\left[A^{*}\left(\bm{T}\circ\frac{A\bm{p}^{k}}{\left|A\bm{p}^{k}\right|}\right)\right]$
(12)

where $\mathcal{P}_{U}(\bm{p}):=\bm{u}\circ\bm{p}/|\bm{p}|$ is the Euclidean
projection onto the constrain set $U$.

This algorithm is analogous to the GS algorithm and follows similar steps.
However, with the aid of the MMDM, this algorithm correctly takes into account
the heterogeneity of the media and uses adaptive step sizes with the aid of
the Armijo line-search parameters $\beta$ and $\sigma$ as shown in Algorithm 1
below following Melde  _et al._ (2023). In words, the algorithm’s steps are as
follows: (1) propagate the pressure from the transducer plane to the rest of
the volume, (2) calculate the gradient of the loss function and back propagate
each plane in the volume back to the hologram plane and obtain the new phase
distribution at the transducer plane, (3) impose the new phase distribution
and the initial transducer amplitude at the transducer plane, and (4) repeat
steps 1-3 until the desired convergence criterion is met.

The algorithm output provides a phase map for the holographic element, which
we translate it into a surface map for fabrication. The resulting thickness
map is then calculated asMelde  _et al._ (2016): The holographic plate starts
with an initial thickness $T_{0}$. Removing material at pixel position $(X,Y)$
causes a relative phase change

$\Delta\phi(X,Y)=\left(k_{\mathrm{m}}-k_{\mathrm{h}}\right)\Delta T(X,X)$

Where $T(X,Y)=T_{0}-\Delta T(X,Y)$ is the thickness of the pixel positioned at
coordinates $(X,Y)$ in the hologram plane and $k_{\mathrm{h}},k_{\mathrm{m}}$
are the wave numbers in the hologram plate and its surrounding medium,
respectively.

Result: Phase pattern $\phi=\arg(p)$

Initial field configuration $p=p_{0}$;

Initial stepsize $\tau=1.0$;

Armijo line-search parameters $\beta=0.5,\sigma=0.1$;

while not converged do

$\partial
p=2A^{*}\left(\left.\frac{Ap}{\left|Ap\right|}\circ\nabla_{x}\ell(x,T)\right|_{x=\left|Ap\right|}\right)$

$p_{\text{new}}=\mathcal{P}_{U}(p-\tau\partial p)$

while $\ell\left(p_{\text{new}}\right)>\ell(p)-\sigma\frac{\tau}{2}\|\partial
p\|^{2}$ do

$\tau=\beta\tau$

$p_{\text{new}}=\mathcal{P}_{U}(p-\tau\partial p)$

end while

$p=p_{\text{new}}$

end while

Gradient descent with line search for computer-generated volumetric holography

##  III Numerical results

Figure 1: Summary of the numerical findings for targeting the desired volumes
within the skull phantom. Columns (a) and (b) depict axial and transverse
sections of the normalized pressure field, respectively. The dashed line
represents the desired volume, whereas the solid white line represents the
skull phantom. Columns (c) and (d) depict the axial and transverse sections of
the MRI scans, respectively. The target volumes’ geometries are emphasized in
distinctive colors.

In this section, we numerically generate four acoustic holograms using the
aforementioned optimization techniques to target four distinct human brain
areas via the temporal window. A flat single-element transducer combined with
an acoustic hologram can correct for skull aberrations and create a focal
point tailored to a certain target location and volume. The geometry of the
skull is derived from a computed tomography (CT) scan and processed to a
resolution of the numerical grid. The skull geometry is sliced along a plane
parallel to the sagittal plane to match the 3D-printed skull in the
experimental setting. The skull section is integrated into the numerical
simulations with a resolution of $0.5$ mm, which corresponds to $\lambda/6.67$
for the transducer drive frequency in water. The numerical grid’s resolution
is uniform in all directions ($dX=dY=dZ$). The geometry and placement of the
therapeutic targets chosen for this investigation have been retrieved from the
Magnetic resonance imaging (MRI) scans and fitted and mapped into their proper
locations in the cranial cavity. The flat circular transducer and hologram are
situated at the temporal bone window, at a distance of $15\mathrm{~{}mm}$
between the hologram’s surface and the skull’s surface. The transducer’s
output is experimentally evaluated in a water tank, with hydrophone
measurements taken at an excitation frequency of $444$ kHz. The numerical
algorithm is constrained by the transducer’s measured experimental amplitude,
and the thickness map is designed using the initial experimental phase, as
detailed in Algorithm 1. In all simulations, the average amplitude at the
transducer’s surface is normalized to $1\mathrm{Pa}$. Water is assumed
everywhere else in the numerical domain since it is a good approximation of
soft tissue and to match the water-submerged experiments. The domain size is
$384X212X199$ grid points in the X, Y, and Z directions, respectively. The
Armijo line-search parameters in Algorithm 1 are set to $\beta=0.9$ and
$\sigma=0.002$ with an initial step size of $\tau=1$. The target volume $T$ in
Algorithm 1 is a binary volumetric matrix with nonzero values at the targeted
brain structure. The numerical simulations of acoustic propagation are
performed using the MMDM and implemented in the open-source toolbox mSOUND on
MATLAB R2019b Gu and Jing (2021). A computational node equipped with AMD EPYC
7702 2.5 GHz (total of 18 cores) and 36 GB of RAM is used for all simulations.
The backward propagation step in Algorithm 1 is the most computationally time-
consuming. In Algorithm 1, the backward propagation stage consumes the most
processing time. The pressure at each plane in the entire volume must be
propagated back to the transducer plane. Because each propagation is
independent, a parallel backpropagation technique was created using MATLAB’s
$parfor$ function, which considerably accelerated the procedure. In the
current study, each iteration lasted an average of 380 seconds, and each
simulation produced satisfactory results after seven iterations on average.
The average overall simulation time was 44 minutes, which is much faster than
previously reported time-reversal techniques.

To assess the proposed algorithm’s performance, we investigate four separate
test cases involving various targets inside the brain. We created four
distinct holograms that focus the ultrasound field on the anterior insula,
hippocampus, caudate nucleus, and amygdala. Each of these targets has distinct
properties in terms of shape, volume, orientation, and spatial placement
within the brain structure. This diversity allows for a more comprehensive and
general assessment of the holograms’ performance capabilities. We shall look
at the axial and transversal distribution of the ultrasonic field within the
brain. In addition, we look into the sonication volume and peak focal pressure
gain inside and outside the target regions. Sonication volumes are calculated
as the region with a pressure gain of at least $-3\mathrm{~{}dB}$ with respect
to the peak pressure in the brain for each test instance. Normalized peak
pressure is calculated as the pressure gain relative to the average pressure
at the transducer surface. The results of the numerical simulations for all
the test cases are summarised in Fig. 1 and quantitatively in Table 1.

###  III.1 Anterior Insula

The anterior insula is a brain area that processes taste, emotions, and
interoceptive awareness, including hunger and pain. It also plays roles in
autonomic functioning and sensory processing, including auditory and
vestibular functions, and is essential for empathy and social cognition Craig
(2009); Uddin (2015). In the present simulation, the insula’s volume is
$0.875$ cm3, and it is positioned at a depth of $8.32$ cm from the hologram’s
surface along the axial direction, making it the deepest target among all the
test cases. The dimensions of the insula extend to lengths of $0.9$ cm, $1.9$
cm, and $1$ cm along the X, Y, and Z axes, respectively. The normalized cross-
section of the maximum pressure field in the YZ and XY planes are shown in the
first row of Fig. 1 (a) and (b), respectively. Each cross-section intersects
the centroid of the target volume. The cartesian axes are aligned with the
axes of the source. The segmented MR image with the highlighted geometry of
the anterior insula is illustrated in 1 (c) and (d). As illustrated in the
axial and transversal cross-sections, the pressure is concentrated on the
Insula volume, which is marked by the white dashed line. Good coverage is
obtained, with the main focal encompassing the majority of the volume.
However, while the energy is focused on the Insula, significant acoustic
pressure occurs both before and after the focal zone. The maximum normalized
pressures inside and outside the target volume are close, at a value of $4.55$
and $4.67$, respectively. A total of $0.45$ $\mathrm{cm}^{3}$ of the Insula
was sonicated, representing $51.4\%$ of the total volume. While $2.25$
$\mathrm{cm}^{3}$ was sonicated outside the target. It is also crucial to note
the low-pressure amplitudes seen in the skull (which is outlined by solid
white lines). This is crucial since it reduces the possibility of skull
heating effects, which improves the therapeutic technique’s safety profile.

###  III.2 Hippocampus

The hippocampus, located in the temporal lobe of the brain, is responsible for
memory formation and storage, spatial navigation, and stress and mood
regulation. It also affects the learning process. Dysfunction in this area can
lead to memory issues, spatial disorientation, and emotional problems
Eichenbaum, Otto, and Cohen (1992). The volume of the segmented Hippocampus is
$3.85$ cm3, and it is located at a depth of $7.55$ cm from the hologram’s
surface along the axial direction. The dimensions of the Hippocampus extend to
lengths of $4.45$ cm, $2.55$ cm, and $1.55$ cm along the X, Y, and Z axes,
respectively. The normalized cross-section of the maximum pressure field in
the YZ and XY planes are shown in the second row of Fig. 1 (a) and (b),
respectively. MRI scans, presented in Fig. 1 (c) and (d), show that the
hippocampus extends more prominently in the transversal plane along the
X-axis. At 4.45 cm, the hippocampus is the largest in terms of length in the X
direction among all the target volumes investigated. The axial and transverse
cross-sections clearly show the pressure focus within the hippocampus volume.
The transverse cross-section shows that pressure is effectively contained
inside the intended area. However, axially, substantial pressure exists
outside the target zone. The maximum normalized pressure inside and outside
the target volume are close, at a value of 2.97 and 3.15, respectively. Using
the hologram $1.61$ $\mathrm{cm}^{3}$ or $41.8\%$ of the Hippocampus is
sonicated. However, a significant portion of the sonicated volume, amounting
to $8.1$ $\mathrm{cm}^{3}$, occurs outside the targeted region, with $83.4\%$
of the total sonicated volume falling outside the intended area. Furthermore,
high-pressure amplitudes in the skull increase the risk of causing unwanted
skull heating effects when performing therapeutic operations.

###  III.3 Caudate Nucleus

The caudate nucleus in the brain’s basal ganglia is crucial for motor control,
procedural learning, reward processing, decision-making, and emotional
processing. Its dysfunction is associated with various neurological and
psychiatric disorders, including Parkinson’s and Huntington’s diseases Grahn,
Parkinson, and Owen (2008). The volume of the segmented caudate nucleus is
$4.73$ cm3, and it is positioned at a depth of $6.2$ cm from the hologram’s
surface along the axial direction, making it the largest target and the closet
to the hologram surface among all the targets considered in the present
investigation. The dimensions of the caudate extend to lengths of $4.1$ cm,
$2.3$ cm, and $1.75$ cm along the X, Y, and Z axes, respectively. The
normalized cross-section of the maximum pressure field in the YZ and XY planes
are shown in the third row of Fig. 1 (a) and (b), respectively. MRI scans, as
presented in Fig. 1 (c) and (d), show that the caudate also extends more
prominently in the transversal plane along the X-axis. The axial and
transverse cross-sections show the pressure concentration inside the caudate’s
volume, however there is less effective coverage of the target area than in
earlier tests. This is emphasized by the relatively smaller sonicated volume
of the caudate, amounting to $1.2$ $\mathrm{cm}^{3}$ or $25.4\%$ of its total
volume. This level of sonication is notably the lowest among all the
considered test cases. The sonicated volume outside of the target volume is
$8.81$ $\mathrm{cm}^{3}$ which is $86.3\%$ of the total sonicated volume. The
maximum normalized pressure inside and outside the target volume are close, at
a value of $3.32$ and $3.38$, respectively.

###  III.4 Amygdala

The amygdala, a tiny region in the brain, is critical in emotional processing,
notably fear and aggression. It is essential for developing emotional memories
and influencing decision-making. The amygdala is also involved in stress and
anxiety reactions, as well as social interactions through the interpretation
of social cues. Amygdala dysfunction can cause emotional disorders, memory
problems, and difficulties with social behaviorJanak and Tye (2015). The
volume of the segmented Amygdala is $0.491$ cm3, and it is positioned at a
depth of $7.49$ cm from the hologram’s surface along the axial direction. The
dimensions of the amygdala extend to lengths of $0.8$ cm, $1.2$ cm, and $1.15$
cm along the X, Y, and Z axes, respectively. The amygdala is the smallest and
most contained target volume among all the test cases, as highlighted by the
MRI scans in Fig. 1 (c) and (d). The normalized cross-section of the maximum
pressure field in the YZ and XY planes are shown in the fourth row of Fig. 1
(a) and (b), respectively. The axial and transversal cross-sections
demonstrate adequate coverage and containment of the focus zone in the target
region. This is highlighted by the high sonication rate within the amygdala
with a total sonication volume of $0.32$ $\mathrm{cm}^{3}$ or $65.3\%$, which
is the highest sonication percentage among all the test cases. The sonicated
volume outside of the target volume is $1.93$ $\mathrm{cm}^{3}$, which
represents $85.8\%$ of the total sonicated volume. The maximum normalized
pressure inside and outside the target volume is the same, at 5.49, which is
the largest gain achieved among all test cases. Low-pressure amplitudes in the
skull were also observed when targeting the amygdala.

Table 1: Summary of the Numerical Results |  Volume [$\mathrm{cm}^{3}$] |  Total sonicated volume [$\mathrm{cm}^{3}$] |  Sonication inside target [$\mathrm{cm}^{3}$] |  Percent of target sonication [%] |  Gain inside target |  Gain outside target
---|---|---|---|---|---|---
Insula |  0.875 |  2.7 |  0.45 |  51.4 |  4.55 |  4.67
Hippocampus |  3.85 |  9.9 |  1.61 |  41.8 |  2.97 |  3.15
Caudate |  4.73 |  8.81 |  1.2 |  25.4 |  3.32 |  3.38
Amygdala |  0.491 |  2.22 |  0.321 |  65.3 |  5.49 |  5.49

###  III.5 Discussion

The pressure fields created by acoustic holograms are inherently constrained
by the diffraction limit and must adhere to the Helmholtz equation and the
principle of energy conservation. This implies that sound fields that do not
conform to these fundamental principles are not physically viable. As a
result, there will inevitably be differences between the theoretically desired
target field and the actual field that is physically realized due to these
inherent limitations and principles governing acoustic propagation. The sound
field can be significantly improved and optimized by carefully considering key
parameters such as the numerical aperture, operating frequency, and the
location and orientation of both the transducer and the hologram in relation
to the size, position, orientation, and geometry of the target volume. This
rigorous calibration guarantees that the acoustic field is fine-tuned to meet
the appropriate targeting characteristics, maximizing the treatment process’s
efficacy, precision, and safety.

The numerical results clearly show that acoustic holograms are more effective
at restricting the focal region within the transversal portion than in the
axial section, especially for bigger targets such as the hippocampus and
caudate. The wider the target’s extent and cross-sectional area in the
transversal plane, the weaker the field’s focus is, resulting in a lower gain
and a focused zone that is more elongated and distributed in the axial
direction, causing greater sonication outside the target volume. This impact
can be minimized by increasing the numerical aperture (transducer size) or
adjusting the operating frequency, albeit both of these modifications have
limitations. Another technique is to orient the transducer so that its axis is
parallel to the target’s longest dimension, resulting in an elongated focal
area that is substantially contained inside the target volume. However, this
strategy has limitations; for example, targeting the hippocampus will
necessitate passing through the parietal bone, which is thicker and more
curved, which could lead to new complications and issues such as skull heating
and mode convergence.

These constraints are inherent in acoustic holographic lenses, not the
optimization algorithm itself. However, the optimization methodology provided
here, with its explicitly defined loss function, provides a method for
somewhat overcoming these limits. Different loss and weighting functions can
be used to direct the optimization algorithm to prioritize specific ultrasonic
parameters. For example, to promote axial resolution by giving axial errors
more weight than transversal errors. This approach could provide better axial
confinement at the expense of transversal containment. Other parameters, such
as pressure uniformity and coverage, can also be optimized to fit specific
therapeutic demands and objectives, resulting in a more specialized
optimization process that improves the efficacy of ultrasonic applications.
The loss functions can be customized further by designating regions with low
pressure, directing the algorithm to avoid sonication in these areas. This
adjustment provides more precise control over pressure distribution,
protecting vital or sensitive areas from unwanted ultrasonic wave exposure.
Furthermore, the proposed algorithm effectively eliminates the uncertainties
inherent in time-reversal techniques, where the impact of the number,
distribution, and locations of virtual sources on focusing parameters is
unclear. Determining the optimal distribution of virtual sources for specific
desired focusing parameters can be challenging due to the non-iterative nature
of these techniques and the computational limitations associated with them
Andrés  _et al._ (2022b).

##  IV Experimental verification

###  IV.1 Experimental Setup

The experimental setup is submerged in a water tank, as shown in Fig. 2. The
walls of the tank are partially covered with Aptflex F28 acoustic absorbing
sheets manufactured by Precision Acoustics, Ltd. The water tank is filled with
degassed and deionized water in a temperature-controlled environment. A
custom-made holder is fabricated to hold the skull, the ultrasound transducer,
and the holographic lens. The position of the transducer can be moved manually
relative to the skull in the X,Y, and Z directions. The transducer is custom-
made and composed of a single circular piezoelectric disc with a diameter of
$60$ mm and an operating frequency of $444$ kHz. The transducer is excited
with a sinusoidal pulse signal (20 cycles) by a signal generator (Keysight
33500B) and an amplifier (Electronics and Innovation A075). The pressure field
on a specified target plane is mapped with the ONDA HNR–0500 needle hydrophone
attached to a custom-built in-house 3D positioning system. The measured signal
from the hydrophone is sampled and digitized by a Tektronix TBS2104
oscilloscope and averaged 64 times to increase the signal-to-noise ratio. To
eliminate the effect of any DC bias and higher frequency noises, the signal is
post-processed using MATLAB by performing an FFT on the time series. A MATLAB
program is used to control the positioning system and triggering of the
function generator. The obtained voltage amplitude is then converted into
acoustic pressure according to the calibration factor provided by the
manufacturer (-250.13 dB). The geometry of the skull is obtained using a
segmented CT scan. A section of the skull is then 3D-printed with
Polycarbonate (PC) FDM by STRATASYS (F270, Eden Prairie, MN) with a smooth
hand finish. The acoustic properties of PC have been previously investigated
underwater Antoniou  _et al._ (2021) with $c=2380$ m/s, $\rho=1153$ kg/m3, and
$\alpha=12.23$ dB/cm at $1.1$ MHz. A Formlabs (Formlabs, Somerville, MA, US)
Form 3 grey resin was used for printing the acoustic lenses with a lateral and
axial resolution of 25 microns. The acoustic properties of the resin were
previously reported in an underwater experimental study Bakaric  _et al._
(2021) with $c=2591$ m/s, $\rho=1178$ kg/m3, and $\alpha=2.922$ dB/cm at $1$
MHz.

Figure 2: (a) shows the experimental setup, including the acoustic source,
3D-printed acoustic hologram, 3D-printed skull phantom, and the hydrophone.
The 3D-printed holograms for targetting the insula (b.1), hippocampus (b.2),
caudate (b.3), and the amygdala (b.4). (c) closeup of the 3D-printed hologram.
Figure 3: (a) The geometry of each target evaluated in the current study. (b)
Schematic of the experimental setup, with the targets positioned in their
proper placements within the cranial cavity of the skull phantom. (c), (d),
(e), and (f) show the experimental and numerical results for the test cases of
the insula, hippocampus, caudate, and amygdala.

###  IV.2 Experimental results

Four holograms are designed based on the numerical results to validate the
proposed algorithm. Underwater pressure scans in both the transversal and
axial directions are carried out for each target using a hydrophone attached
to a positioning system. Due to the hydrophone’s limited mobility range within
the skull, only smaller portions surrounding the target location are scanned.
The experimental measurements are then compared to numerical simulations.
Figure 3 (a) and (b), respectively, show the geometry of the four target
volumes and their corresponding location inside the 3D-printed skull phantom.
Figure 3 (c), (d), (e), and (f) respectively compare the experimental
measurements and their corresponding numerical results for the insula,
hippocampus, caudate, and amygdala. Good overall agreement is observed between
the experimental measurements and the numerical results. The maximum
experimental pressures measured in the transversal plane were $85.2$ kPa for
the insula, $98.5$ kPa for the Hippocampus, $127$ kPa for the Caudate, and
$98.2$ kPa for the Amygdala with an excitation input of $200$ Vpp.

For a quantitative comparison between experimental and numerical results, we
calculated the focusing efficiency Sallam  _et al._ (2021) in the transversal
section. By comparing the acoustic power within the focal region to the total
acoustic power in the entire plane. The calculations assume planar wave
propagation. The focusing efficiency is given as:

$\eta=\frac{\Pi_{focal}}{\Pi_{total}}=\frac{\sum
P_{focal}^{2}(X,Y,Z_{0})}{\sum P^{2}(X,Y,Z_{0})}$ (13)

Where $P_{focal}(X,Y,Z_{0})$ and $\Pi_{focal}$ are the acoustic pressure
distribution and the acoustic power in the focal regions, respectively.
$\Pi_{total}$ is the total acoustic power in the transversal plane. Examining
Fig. 3 (c), which presents insula results, we notice more pressure leakage
outside of the focal region in the experimental data. This is evident from the
disparity in focusing efficiency, with the simulations calculated at $76.5\%$
compared to $46.1\%$ in the experimental measurements. A similar pattern is
observed in other test cases: the hippocampus has a numerical efficiency of
$84.2\%$ compared to an experimental efficiency of $49.2\%$, and the amygdala
with a numerical efficiency of $70.1\%$ against an experimental efficiency of
$41.6\%$. Interestingly, the caudate is the only test case where the
experimental results outperform the simulations, with an efficiency of
$64.2\%$ versus $57.1\%$, potentially due to a more pronounced top lobe in the
experimental data as shown in Fig. 3(e). The differences between experimental
and numerical results could be attributable to factors such as measurement
mistakes, alignment challenges, variances in material qualities, and flaws in
the produced skull phantom and hologram.

##  V Conclusions

We introduced an efficient volumetric holographic technique for creating
3D-printed acoustic holograms specifically for low-intensity non-thermal tFUS
applications. The fundamental goal of these acoustic lenses is to generate
precise FUS fields within the brain while successfully compensating for
aberrations induced by the skull’s structure. The proposed optimization
strategy uses the modified mixed domain method to perform efficient ultrasonic
propagation in strongly heterogeneous media, which is then incorporated into a
gradient descent iterative volumetric holography optimization algorithm. The
holographic approach outperforms the previously reported time-reversal
strategies used in full-wave time-domain numerical simulations. Using the
proposed approach, a hologram can be created in less than 50 minutes.
Furthermore, due to the more efficient propagation, iterative optimization
techniques can be used with explicitly assigned loss functions. The
optimization technique can be fine-tuned to prioritize individual ultrasound
characteristics by using different loss and weighting functions. For example,
increasing axial resolution entails giving greater weight to axial errors than
transversal ones, which improves axial containment but may reduce transversal
containment. Other factors like pressure uniformity and coverage can also be
adjusted to suit particular therapeutic requirements, enabling a more
customized optimization process. The loss functions can also be specifically
tailored to designate regions where minimal pressure is crucial, thereby
directing the algorithm to prevent sonication in those sensitive areas,
ensuring a more precise and effective ultrasound treatment. We applied the
proposed algorithm to target four brain structures: the anterior insula,
hippocampus, caudate nucleus, and amygdala. Verification was done through
underwater experiments using realistic 3D-printed skull phantoms. Good
agreement was observed between the numerical results and the measured
experimental data. Our results demonstrate that the acoustic holograms
achieved effective coverage and confinement in the transversal plane. However,
in the axial direction, sonication leakage varied depending on the target
geometry, orientation, and volume, with lower axial containment observed in
larger and transversely elongated targets like the hippocampus and caudate.
The proposed approach enables more rapid computation of acoustic holograms,
enabling quicker administration of therapeutic procedures and allowing for
more flexibility in tailoring ultrasound fields for specific therapeutic
requirements.

## ACKNOWLEDGEMENT

This work was supported by the U.S. National Science Foundation (NSF) under
CAREER grant, Award No. CMMI $2143788$, which is gratefully acknowledged.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## REFERENCES

## References

  * Meng, Hynynen, and Lipsman (2021) Y. Meng, K. Hynynen, and N. Lipsman, Nature Reviews Neurology 17, 7 (2021).
  * Darmani  _et al._ (2022) G. Darmani, T. Bergmann, K. B. Pauly, C. Caskey, L. De Lecea, A. Fomenko, E. Fouragnan, W. Legon, K. Murphy, T. Nandi, _et al._ , Clinical Neurophysiology 135, 51 (2022).
  * Legon  _et al._ (2014) W. Legon, T. F. Sato, A. Opitz, J. Mueller, A. Barbour, A. Williams, and W. J. Tyler, Nature neuroscience 17, 322 (2014).
  * Lee  _et al._ (2016) W. Lee, H.-C. Kim, Y. Jung, Y. A. Chung, I.-U. Song, J.-H. Lee, and S.-S. Yoo, Scientific reports 6, 34026 (2016).
  * Mueller  _et al._ (2014) J. Mueller, W. Legon, A. Opitz, T. F. Sato, and W. J. Tyler, Brain stimulation 7, 900 (2014).
  * Legon  _et al._ (2018a) W. Legon, L. Ai, P. Bansal, and J. K. Mueller, Human brain mapping 39, 1995 (2018a).
  * Deffieux and Konofagou (2010) T. Deffieux and E. E. Konofagou, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 57, 2637 (2010).
  * Legon  _et al._ (2018b) W. Legon, L. Ai, P. Bansal, and J. K. Mueller, Human brain mapping 39, 1995 (2018b).
  * Legon  _et al._ (2018c) W. Legon, P. Bansal, R. Tyshynsky, L. Ai, and J. K. Mueller, Scientific reports 8, 10007 (2018c).
  * Strohman  _et al._ (2024) A. Strohman, B. Payne, A. In, K. Stebbins, and W. Legon, Journal of Neuroscience (2024).
  * Lipsman  _et al._ (2013) N. Lipsman, M. L. Schwartz, Y. Huang, L. Lee, T. Sankar, M. Chapman, K. Hynynen, and A. M. Lozano, The Lancet Neurology 12, 462 (2013).
  * Xu  _et al._ (2021) Z. Xu, T. L. Hall, E. Vlaisavljevich, and F. T. Lee Jr, International Journal of Hyperthermia 38, 561 (2021).
  * Kim  _et al._ (2014) Y. Kim, T. L. Hall, Z. Xu, and C. A. Cain, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 61, 582 (2014).
  * Sukovich  _et al._ (2018) J. R. Sukovich, C. A. Cain, A. S. Pandey, N. Chaudhary, S. Camelo-Piragua, S. P. Allen, T. L. Hall, J. Snell, Z. Xu, J. M. Cannata, _et al._ , Journal of neurosurgery 131, 1331 (2018).
  * Kyriakou  _et al._ (2014) A. Kyriakou, E. Neufeld, B. Werner, M. M. Paulides, G. Szekely, and N. Kuster, International journal of hyperthermia 30, 36 (2014).
  * Almquist, Parker, and Christensen (2016) S. Almquist, D. L. Parker, and D. A. Christensen, Journal of therapeutic ultrasound 4, 1 (2016).
  * Maimbourg  _et al._ (2018) G. Maimbourg, A. Houdouin, T. Deffieux, M. Tanter, and J.-F. Aubry, Physics in Medicine & Biology 63, 025026 (2018).
  * Kook  _et al._ (2023) G. Kook, Y. Jo, C. Oh, X. Liang, J. Kim, S.-M. Lee, S. Kim, J.-W. Choi, and H. J. Lee, Microsystems & Nanoengineering 9, 45 (2023).
  * Hu  _et al._ (2022) Z. Hu, Y. Yang, L. Xu, Y. Hao, and H. Chen, Frontiers in Neuroscience 16, 984953 (2022).
  * Jiménez-Gambín  _et al._ (2021) S. Jiménez-Gambín, N. Jiménez, A. N. Pouliopoulos, J. M. Benlloch, E. E. Konofagou, and F. Camarena, IEEE Transactions on Biomedical Engineering 69, 1359 (2021).
  * He  _et al._ (2021) J. He, J. Wu, Y. Zhu, Y. Chen, M. Yuan, L. Zeng, and X. Ji, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 69, 662 (2021).
  * Jiménez-Gambín  _et al._ (2019) S. Jiménez-Gambín, N. Jiménez, J. M. Benlloch, and F. Camarena, Physical Review Applied 12, 014016 (2019).
  * Jiménez-Gambín, Jiménez, and Camarena (2020) S. Jiménez-Gambín, N. Jiménez, and F. Camarena, Physical Review Applied 14, 054070 (2020).
  * Sallam  _et al._ (2021) A. Sallam, V. C. Meesala, M. R. Hajj, and S. Shahab, Applied Physics Letters 119 (2021).
  * Sallam, Meesala, and Shahab (2021) A. Sallam, V. C. Meesala, and S. Shahab, in  _Active and Passive Smart Structures and Integrated Systems XV_ , Vol. 11588 (SPIE, 2021) p. 1158807.
  * Sallam and Shahab (2022) A. Sallam and S. Shahab, Applied Physics Letters 121 (2022).
  * Melde  _et al._ (2016) K. Melde, A. G. Mark, T. Qiu, and P. Fischer, Nature 537, 518 (2016).
  * Bakhtiari-Nejad  _et al._ (2018) M. Bakhtiari-Nejad, A. Elnahhas, M. R. Hajj, and S. Shahab, Journal of Applied Physics 124 (2018).
  * Melde  _et al._ (2023) K. Melde, H. Kremer, M. Shi, S. Seneca, C. Frey, I. Platzman, C. Degel, D. Schmitt, B. Schölkopf, and P. Fischer, Science Advances 9, eadf6182 (2023).
  * Kim  _et al._ (2021) J. Kim, S. Kasoji, P. G. Durham, and P. A. Dayton, Applied Physics Letters 118 (2021).
  * Andrés  _et al._ (2022a) D. Andrés, J. Vappou, N. Jiménez, and F. Camarena, Applied Physics Letters 120 (2022a).
  * Treeby and Cox (2010) B. E. Treeby and B. T. Cox, Journal of biomedical optics 15, 021314 (2010).
  * Sallam and Shahab (2023) A. Sallam and S. Shahab, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control (2023).
  * Gu and Jing (2020) J. Gu and Y. Jing, The Journal of the Acoustical Society of America 147, 4055 (2020).
  * Gu and Jing (2018) J. Gu and Y. Jing, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 65, 1258 (2018).
  * Gu and Jing (2021) J. Gu and Y. Jing, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 68, 1476 (2021).
  * Zhang  _et al._ (2017) J. Zhang, N. Pégard, J. Zhong, H. Adesnik, and L. Waller, Optica 4, 1306 (2017).
  * Craig (2009) A. D. Craig, Nature reviews neuroscience 10, 59 (2009).
  * Uddin (2015) L. Q. Uddin, Nature reviews neuroscience 16, 55 (2015).
  * Eichenbaum, Otto, and Cohen (1992) H. Eichenbaum, T. Otto, and N. J. Cohen, Behavioral and neural biology 57, 2 (1992).
  * Grahn, Parkinson, and Owen (2008) J. A. Grahn, J. A. Parkinson, and A. M. Owen, Progress in neurobiology 86, 141 (2008).
  * Janak and Tye (2015) P. H. Janak and K. M. Tye, Nature 517, 284 (2015).
  * Andrés  _et al._ (2022b) D. Andrés, N. Jiménez, J. M. Benlloch, and F. Camarena, Ultrasound in Medicine & Biology 48, 872 (2022b).
  * Antoniou  _et al._ (2021) A. Antoniou, N. Evripidou, M. Giannakou, G. Constantinides, and C. Damianou, The Journal of the Acoustical Society of America 149, 2854 (2021).
  * Bakaric  _et al._ (2021) M. Bakaric, P. Miloro, A. Javaherian, B. T. Cox, B. E. Treeby, and M. D. Brown, The Journal of the Acoustical Society of America 150, 2798 (2021).


#  A hyperbolic proof of Pascal’s Theorem

Miguel Acosta  and  Jean-Marc Schlenker  Department of mathematics, Université
du Luxembourg, Maison du Nombre, 6, Avenue de la Fonte, L-4364 Esch-sur-
Alzette, Luxembourg<EMAIL_ADDRESS>, jean-
<EMAIL_ADDRESS>
###### Abstract.

We provide a simple proof of Pascal’s Theorem on cyclic hexagons, as well as a
generalization by Möbius, using hyperbolic geometry.

Both authors were partially supported by FNR grant AGoLoM OPEN/16/11405402

##  1\. Pascal’s Theorem

Blaise Pascal (1623–1662) is a towering intellectual figure of the XVIIth
century. He is credited with inventing and building the first mechanical
calculator, the _Pascaline_ , and with laying the foundations of probability
theory, in particular in his correspondence with Fermat – he came up for
instance with Pascal’s triangle. He is also known for his work on hydrostatics
and Pascal’s law (as well as the invention of the syringe), and for
discovering the variation of air pressure with altitude – the SI unit of
pressure is called the Pascal. However, he was most influential in his time as
a philosopher and theologian, and well-known for “Pascal’s wager”.

Pascal was raised and educated by his father, Etienne Pascal, who had a strong
interest in the intellectual developments of his time and was an active member
of a group of scientists meeting around Martin Mersenne, including Desargues,
Descartes and others. According to contemporary sources [7, p. 176], Blaise
was extraordinarily precocious, and so passionate in studying mathematics
that, when he was 11, his father forbid him to read any mathematics book
before he turned 15 and knew latin and greek. Blaise therefore continued
studying geometry by himself and in secret and, at 16, published his first
treatise on projective geometry, the _Traité sur les coniques_ , which
contained the following remarkable statement.

######  Theorem A (Pascal’s Theorem).

Let $ABCDEF$ be a cyclic hexagon. Let $X$ be the intersection point of $AB$
and $DE$, $Y$ the intersection point of $BC$ and $EF$ and $Z$ the intersection
point of $CD$ and $FA$. Then $X$, $Y$ and $Z$ are aligned.

The main goal of this note is to provide a simple proof, based on hyperbolic
geometry, of Pascal’s Theorem. The statement has a natural setting in the
projective plane. Instead of considering a cyclic hexagon, one then considers
a hexagon with vertices on a conic. Any non-degenerate conic is projectively
equivalent to a circle, while the statements for degenerate conics can be
obtained by a limiting argument where a hexagon with vertices on a degenerate
conic is obtained as a limit of hexagons with vertices on non-degenerate
conics.

By considering different permutations of the set of vertices, the
configuration of Pascal’s theorem leads to a fascinating set of $95$ lines and
$95$ points with a beautiful incidence structure. We recommend the articles of
Conway and Ryba [1] and [2] describing the full configuration, that they call
_mysticum hexagrammaticum_. On the one hand, this name describes their
approach to label the lines and points by permutations of the set
$\\{A,B,C,D,E,F\\}$, giving a natural way to state incidence relations and be
able to identify symmetries. On the other hand, it recalls the original name
of A, given by Pascal in $1639$: _hexagrammum mysticum_. The original work
from Pascal did not survive, but there are many proofs of A, using a wide
variety of tools. A proof using algebraic geometry is sketched by Conway and
Ryba in [1]. Other proofs involve cross-ratios and symmetries of the
projective plane, as in [6] or Euclidean lengths and Menalaus’s theorem like
in [3, p.77].

The proof given here, based on hyperbolic geometry, is not really novel (it
can be deduced easily from [6, page 436]) but it brings to light a striking
link between Pascal’s theorem and elementary hyperbolic geometry. A recent and
similar link is exhibited by Drach and Schwarz in [4], where they revisit the
seven-circles theorem, in Euclidean geometry, in terms of hyperbolic geometry.
Even if those links between hyperbolic geometry and results on projective or
Euclidean geometry do not lead to novel results, they give a beautiful
perspective to them.

##  2\. A hyperbolic statement of Pascal’s Theorem

$A$$B$$C$$D$$E$$F$$Z$$X$$Y$ Figure 1. The configuration of Pascal’s theorem
(A) and B.

We are going to use the Klein model of the hyperbolic plane. Consider an open
disk $\Delta$ on the projective plane, bounded by a circle $\Gamma$. The
_hyperbolic plane_ is defined as $\Delta$, endowed with the _Hilbert distance_
, defined as follows

######  Definition.

Let $P,Q\in\Delta$, and $A,B$ be the intersection points of the line $PQ$ with
$\Gamma$. Then

$d(P,Q)=\frac{1}{2}\log\left(\frac{BP}{BQ}\cdot\frac{AQ}{AP}\right).$

This distance induces a complete Riemannian metric on $\Delta$, and thus
notions of angles and lengths. A key property of this Hilbert distance is that
it is invariant under projective transformations that leave $\Delta$
invariant. The geodesics are precisely the straight lines, but the angles are
not the Euclidean ones. The circle $\Gamma$ is then the _boundary at infinity
of the hyperbolic plane_ , and its points are called _ideal points_. The model
is very rich, but we focus here only on a few points related to the polarity
in the projective plane and orthogonality in the hyperbolic plane.

Given a circle in the projective plane, one can define a polarity relation
between points and lines: there is a _polar line_ for each point, and a _pole_
for each line. By definition, given points $P\in\Delta$ and
$Q\not\in\bar{\Delta}$, $Q$ is in the polar line of $P$, and conversely, if
and only if

$\frac{BP}{BQ}\cdot\frac{AQ}{AP}=-1~{},$

where $A$ and $B$ are again the intersections of the line $PQ$ with $\Gamma$.
We will denote the polar line of $P$ by $P^{*}$. The polarity relation is also
invariant under projective transformations that leave $\Delta$ invariant.

It follows from the definition that if $P,Q\in\Delta$, then $P$ and $Q$ are
both in the polar line of $P^{*}\cap Q^{*}$, the intersection point of the
polar lines $P^{*}$ and $Q^{*}$. It follows that the line $PQ$ is the polar
line of $P^{*}\cap Q^{*}$. As a consequence, the polarity relation preserves
the incidence: three points are aligned if and only if their polar lines are
concurrent.

This last fact allows to have _dual_ statements. For example, the dual
statement of Pascal’s theorem is Brianchon’s theorem: if a conic is inscribed
in a hexagon, then the three diagonals joining the opposite vertices of the
hexagon are concurrent.

Back to hyperbolic geometry, the link between polarity and the Klein model
that we are going to use is the following:

######  Proposition.

Two lines $l_{1}$ and $l_{2}$ in the hyperbolic plane are orthogonal if and
only if the pole of $l_{1}$ is contained in the extension of $l_{2}$ to the
projective plane.

In particular, if $l_{1}$ and $l_{2}$ are two lines that do not intersect in
the hyperbolic plane, they have a unique common perpendicular, that is the
polar line of their intersection point in the projective plane.

The proposition follows from the projective invariance of the hyperbolic
metric and of the polarity relation under projective transformations leaving
$\Delta$ invariant, since one can always find such a projective transformation
leaving $\Delta$ invariant and bringing the intersection of $l_{1}$ and
$l_{2}$ to the center of $\Delta$ – in this case the proposition is easy to
check.

We can now restate Pascal’s theorem in terms of hyperbolic geometry.
Considering the polar lines $l_{1}$, $l_{2}$ and $l_{3}$ of the points $X$,
$Y$ and $Z$ of the statement of A, we obtain the following equivalent
statement.

######  Proposition B.

Let $ABCDEF$ be an ideal hyperbolic hexagon. Let $l_{1}$ be the common
perpendicular to $AB$ and $DE$, $l_{2}$ the common perpendicular to $BC$ and
$EF$ and $l_{3}$ the common perpendicular to $CD$ and $FA$. Then $l_{1}$,
$l_{2}$ and $l_{3}$ are concurrent.

##  3\. Proof of Proposition B

The proof of Proposition B is based on the hyperbolic version of an elementary
statement on triangles.

######  Theorem C.

Let $PQR$ be a triangle. Then the angle bisectors of $PQR$ are concurrent.

######  Remark.

C holds for Euclidean, spherical and hyperbolic triangles. The proof is, in
all three cases, elementary. A point of a triangle is in the bisector of an
angle if and only if it is at equal distance from the two corresponding edges.
Therefore, the intersection point of two angle bisectors is at equal distance
from all three edges, and is therefore contained in the third bisector.

######  Lemma D.

Let $ABDE$ be an ideal hyperbolic quadrilateral. Then, the common orthogonal
to $AB$ and $DE$ is the angle bisector of the lines $AD$ and $BE$.

###### Proof.

Let $l$ be the angle bisector of the lines $AD$ and $BE$. The hyperbolic
reflection on $l$ exchanges $A$ and $B$, so the line $AB$ is preserved by the
reflection. Thus, the angle bisector of the two lines is orthogonal to $AB$.
Similarly, $l$ is also perpendicular to $DE$, so it is precisely the common
perpendicular to $AB$ and $DE$. ∎

$A$$B$$C$$D$$E$$F$$P$$Q$$R$ Figure 2. The triangle $PQR$ and its hyperbolic
angle bisectors.

###### Proof of Proposition B.

If the diagonals $AD$, $BE$ and $CF$ are concurrent at a point $P$, then D
implies that $l_{1}$, $l_{2}$ and $l_{3}$ contain $P$, and are therefore
concurrent.

Now, suppose that the diagonals $AD$, $BE$ and $CF$ are not concurrent. We
call $P$ the intersection point of $BE$ and $CF$, $Q$ the intersection point
of $AD$ and $CF$, and $R$ the intersection point of $AD$ and $BE$.

It follows from Lemma D that the angle bisector of $PQR$ at $P$ is the common
perpendicular $l_{1}$ to $AB$ and $DE$, while the angle bisector of $PQR$ at
$Q$ is the common perpendicular $l_{2}$ to $BC$ and $EF$, and the angle
bisector of $PQR$ at $R$ is the common perpendicular $l_{3}$ to $CD$ and $FA$.

By Theorem C applied to the triangle $PQR$, the lines $l_{1}$, $l_{2}$ and
$l_{3}$ are concurrent, and the result follows. ∎

##  4\. The Möbius generalization

In 1847, Möbius proved a generalization of Pascal’s Theorem for $(4n+2)$-gons
[5].

######  Theorem E (Möbius, 1847).

Let $A_{1}A_{2}\cdots A_{4n+2}$ be a cyclic $(4n+2)$-gon. Let
$X_{1},\dots,X_{2n+1}$ be the intersection points of the pairs of opposite
sides of $A_{1}A_{2}\cdots A_{4n+2}$. If $X_{1},\dots,X_{2n}$ are aligned,
then $X_{2n+1}$ lies in the same line as $X_{1},\dots,X_{2n}$.

By considering the polar lines $l_{1},\dots,l_{2n+1}$ of
$X_{1},\dots,X_{2n+1}$, we obtain the corresponding hyperbolic statement, for
which the proof of B extends easily.

######  Proposition F.

Let $A_{1}A_{2}\cdots A_{4n+2}$ be a hyperbolic ideal $(4n+2)$-gon. Let
$l_{1},\dots,l_{2n+1}$ be the common perpendicular to the pairs of opposite
sides of $A_{1}A_{2}\cdots A_{4n+2}$. If $l_{1},\dots,l_{2n}$ are concurrent,
then the common intersection point belongs also to $l_{2n+1}$.

For each $i\in\\{1,\dots,2n+1\\}$, let $m_{i}$ be the line joining $A_{i}$ and
$A_{i+2n+1}$. Let $R_{i}$ be the be the union of the two quarters defined by
$m_{i}$ and $m_{i+1}$ that contain the sides $A_{i}A_{i+1}$ and
$A_{i+2n+1}A_{i+2n+2}$, as in Figure 3, where the indices are taken modulo
$2n+1$. Observe that, by D, the line $l_{i}$ is the set of points of $R_{i}$
that is at the same distance from $m_{i}$ and $m_{i+1}$. If
$l_{1},\dots,l_{2n}$ are concurrent at a point $P$, then $P$ is at the same
distance from all the lines $m_{i}$, so, in particular, it is at the same
distance from $m_{2n+1}$ and $m_{1}$. The only remaining point to complete the
proof (and the reason the statement is false for $4n$-gons) is given by the
following lemma.

$A_{i}$$A_{i+1}$$A_{i+2n+1}$$A_{i+2n+2}$$m_{i+1}$$m_{i}$$l_{i}$ Figure 3. The
region $R_{i}$.

######  Lemma G.

If $P\in R_{1}\cap\cdots\cap R_{2n}$, then $P\in R_{1}\cap\cdots\cap R_{2n+1}$

###### Proof.

Consider Cartesian equations for the lines $m_{i}$, that we still denote by
$m_{i}$, so $m_{i}(A_{i})=m_{i}(A_{i+2n+1})=0$. Up to changing signs, we can
suppose that $m_{i}(A_{i+1})>0$. Thus, for each $i\in\\{1,\dots,2n\\}$ the
region $R_{i}$ is defined by the inequality $m_{i}m_{i+1}<0$, but the region
$R_{2n+1}$, bounded by $m_{2n+1}$ and $m_{1}$, is defined by
$m_{2n+1}m_{1}>0$. Now, if $P\in R_{1}\cap\cdots\cap R_{2n}$, then
$m_{1}(P)m_{2}(P)<0$, $\dots$ , $m_{2n}(P)m_{2n+1}(P)<0$. By multiplying this
even number of inequalities, we obtain $m_{2n+1}(P)m_{1}(P)>0$, so $P\in
R_{2n+1}$. ∎

### Acknowledgement

The authors are grateful to Arseniy Akopyan for pointing out reference [6].

## References

  * [1] John Conway and Alex Ryba.  The Pascal Mysticum Demystified.  The Mathematical Intelligencer, 34(3):4–8, September 2012.
  * [2] John Conway and Alex Ryba.  Extending the Pascal Mysticum.  The Mathematical Intelligencer, 35(2):44–51, June 2013.
  * [3] Harold S. M. Coxeter and S. L. Greitzer.  Geometry Revisited.  American Mathematical Society, Washington, UNITED STATES, 1967.
  * [4] Kostiantyn Drach and Richard Evan Schwartz.  A Hyperbolic View of the Seven Circles Theorem.  The Mathematical Intelligencer, December 2019.
  * [5] A. F. Möbius.  Verallgemeinerung des Pascalschen Theorems, das in einen Kegelschnitt beschriebene Sechseck betreffend.  Journal für die reine und angewandte Mathematik, 1848(36):216–220, January 1848.
  * [6] Jürgen Richter-Gebert.  Perspectives on Projective Geometry: A Guided Tour Through Real and Complex Geometry.  Springer Science & Business Media, February 2011.
  * [7] Gédéon Tallemant des Réaux.  Les Historiettes de Tallemant des Réaux: Mémoires pour servir à l’histoire du XVIIe siècle.  Delloye, Paris, France, 1840.