* Jolliffe (2002) Jolliffe, I. T. 2002, Principal Component Analysis (New York: Springer New York, NY)
  * Kauffmann et al. (2017) Kauffmann, J., Goldsmith, P. F., Melnick, G., et al. 2017, A&A, 605, L5, doi: 10.1051/0004-6361/201731123
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Krieger et al. (2019) Krieger, N., Bolatto, A. D., Walter, F., et al. 2019, ApJ, 881, 43, doi: 10.3847/1538-4357/ab2d9c
  * Krieger et al. (2019) Krieger, N., Bolatto, A. D., Walter, F., et al. 2019, The Astrophysical Journal, 881, 43, doi: 10.3847/1538-4357/ab2d9c
  * Krieger et al. (2020) Krieger, N., Bolatto, A. D., Leroy, A. K., et al. 2020, ApJ, 897, 176, doi: 10.3847/1538-4357/ab9c23
  * Leroy et al. (2015) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, The Astrophysical Journal, 801, doi: 10.1088/0004-637X/801/1/25
  * Leroy et al. (2018) —. 2018, The Astrophysical Journal, 869, 126, doi: 10.3847/1538-4357/aaecd1
  * Levy et al. (2021) Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2021, ApJ, 912, 4, doi: 10.3847/1538-4357/abec84
  * Levy et al. (2022) —. 2022, ApJ, 935, 19, doi: 10.3847/1538-4357/ac7b7a
  * Liszt & Turner (1978) Liszt, H. S., & Turner, B. E. 1978, ApJ, 224, L73, doi: 10.1086/182762
  * Lupton et al. (2004) Lupton, R., Blanton, M. R., Fekete, G., et al. 2004, PASP, 116, 133, doi: 10.1086/382245
  * MacQueen (1967) MacQueen, J. 1967, in Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, 281–297
  * Martín et al. (2019) Martín, S., Martín-Pintado, J., Blanco-Sánchez, C., et al. 2019, A&A, 631, A159, doi: 10.1051/0004-6361/201936144
  * Martín et al. (2009) Martín, S., Martín-Pintado, J., & Viti, S. 2009, The Astrophysical Journal, 706, 1323, doi: 10.1088/0004-637X/706/2/1323
  * Martín et al. (2003) Martín, S., Mauersberger, R., Martín-Pintado, J., García-Burillo, S., & Henkel, C. 2003, A&A, 411, L465, doi: 10.1051/0004-6361:20031442
  * Martín et al. (2006) Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, ApJS, 164, 450, doi: 10.1086/503297
  * Martín et al. (2011) Martín, S., Krips, M., Martín-Pintado, J., et al. 2011, A&A, 527, A36, doi: 10.1051/0004-6361/201015855
  * Martín et al. (2015) Martín, S., Kohno, K., Izumi, T., et al. 2015, A&A, 573, A116, doi: 10.1051/0004-6361/201425105
  * Martín et al. (2021) Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46, doi: 10.1051/0004-6361/202141567
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Meier & Turner (2005) Meier, D. S., & Turner, J. L. 2005, ApJ, 618, 259, doi: 10.1086/426499
  * Meier & Turner (2012) —. 2012, ApJ, 755, 104, doi: 10.1088/0004-637X/755/2/104
  * Meier et al. (2015) Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ, 801, 63, doi: 10.1088/0004-637X/801/1/63
  * Mills et al. (2021) Mills, E. A. C., Gorski, M., Emig, K. L., et al. 2021, ApJ, 919, 105, doi: 10.3847/1538-4357/ac0fe8
  * Müller et al. (2005) Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, Journal of Molecular Structure, 742, 215, doi: 10.1016/j.molstruc.2005.01.027
  * Müller et al. (2001) Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49, doi: 10.1051/0004-6361:20010367
  * Muller et al. (2013) Muller, S., Beelen, A., Black, J. H., et al. 2013, A&A, 551, A109, doi: 10.1051/0004-6361/201220613
  * Nakajima et al. (2018) Nakajima, T., Takano, S., Kohno, K., Harada, N., & Herbst, E. 2018, PASJ, 70, 7, doi: 10.1093/pasj/psx153
  * Nishimura et al. (2017) Nishimura, Y., Watanabe, Y., Harada, N., et al. 2017, ApJ, 848, 17, doi: 10.3847/1538-4357/aa89e5
  * Pearson (1901) Pearson, K. 1901, Philosophical Magazine, 2, 559â572, doi: 10.1080/14786440109462720
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pety et al. (2012) Pety, J., Gratier, P., Guzmán, V., et al. 2012, A&A, 548, A68, doi: 10.1051/0004-6361/201220062
  * Pety et al. (2017) Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98, doi: 10.1051/0004-6361/201629862
  * Quan et al. (2010) Quan, D., Herbst, E., Osamura, Y., & Roueff, E. 2010, ApJ, 725, 2101, doi: 10.1088/0004-637X/725/2/2101
  * Rekola et al. (2005) Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330, doi: 10.1111/j.1365-2966.2005.09166.x
  * Rico-Villas et al. (2020) Rico-Villas, F., Martín-Pintado, J., González-Alfonso, E., Martín, S., & Rivilla, V. M. 2020, MNRAS, 491, 4573, doi: 10.1093/mnras/stz3347
  * Saito et al. (2022) Saito, T., Takano, S., Harada, N., et al. 2022, ApJ, 935, 155, doi: 10.3847/1538-4357/ac80ff
  * Sakamoto et al. (2011) Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, The Astrophysical Journal, 735, doi: 10.1088/0004-637X/735/1/19
  * Sakamoto et al. (2021) Sakamoto, K., Martín, S., Wilner, D. J., et al. 2021, ApJ, 923, 240, doi: 10.3847/1538-4357/ac29bf
  * Schinnerer et al. (2019) Schinnerer, E., Hughes, A., Leroy, A., et al. 2019, ApJ, 887, 49, doi: 10.3847/1538-4357/ab50c2
  * Shirley (2015) Shirley, Y. L. 2015, PASP, 127, 299, doi: 10.1086/680342
  * Sorai et al. (2000) Sorai, K., Nakai, N., Kuno, N., Nishiyama, K., & Hasegawa, T. 2000, PASJ, 52, 785, doi: 10.1093/pasj/52.5.785
  * Suutarinen et al. (2014) Suutarinen, A. N., Kristensen, L. E., Mottram, J. C., Fraser, H. J., & van Dishoeck, E. F. 2014, MNRAS, 440, 1844, doi: 10.1093/mnras/stu406
  * Tafalla et al. (2021) Tafalla, M., Usero, A., & Hacar, A. 2021, A&A, 646, A97, doi: 10.1051/0004-6361/202038727
  * Takano et al. (2019) Takano, S., Nakajima, T., & Kohno, K. 2019, PASJ, 71, S20, doi: 10.1093/pasj/psz020
  * Takano et al. (2014) Takano, S., Nakajima, T., Kohno, K., et al. 2014, PASJ, 66, 75, doi: 10.1093/pasj/psu052
  * Takeuchi et al. (2022) Takeuchi, T. T., Yata, K., Aoshima, M., et al. 2022, arXiv e-prints, arXiv:2203.04535, doi: 10.48550/arXiv.2203.04535
  * Tanaka et al. (2018) Tanaka, K., Nagai, M., Kamegai, K., Iino, T., & Sakai, T. 2018, ApJS, 236, 40, doi: 10.3847/1538-4365/aab9a5
  * Tanaka et al. (2023) Tanaka, K., Mangum, J., Viti, S., et al. 2023, submitted to ApJ
  * Taniguchi et al. (2022) Taniguchi, K., Tanaka, K. E. I., Zhang, Y., et al. 2022, ApJ, 931, 99, doi: 10.3847/1538-4357/ac69d1
  * Tercero et al. (2020) Tercero, B., Cernicharo, J., Cuadrado, S., de Vicente, P., & Guélin, M. 2020, A&A, 636, L7, doi: 10.1051/0004-6361/202037837
  * Thaddeus et al. (1981) Thaddeus, P., Guelin, M., & Linke, R. A. 1981, ApJ, 246, L41, doi: 10.1086/183549
  * Turner (1985) Turner, B. E. 1985, ApJ, 299, 312, doi: 10.1086/163702
  * Turner & Ho (1985) Turner, J. L., & Ho, P. T. P. 1985, The Astrophysical Journal, 299, L77, doi: 10.1086/184584
  * Ulvestad & Antonucci (1997) Ulvestad, J. S., & Antonucci, R. R. J. 1997, The Astrophysical Journal, 488, 621, doi: 10.1086/304739
  * Ungerechts et al. (1997) Ungerechts, H., Bergin, E. A., Goldsmith, P. F., et al. 1997, ApJ, 482, 245, doi: 10.1086/304110
  * Veilleux et al. (2020) Veilleux, S., Maiolino, R., Bolatto, A. D., & Aalto, S. 2020, A&A Rev., 28, 2, doi: 10.1007/s00159-019-0121-9
  * Walter et al. (2017) Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265, doi: 10.3847/1538-4357/835/2/265
  * Watanabe et al. (2017) Watanabe, Y., Nishimura, Y., Harada, N., et al. 2017, ApJ, 845, 116, doi: 10.3847/1538-4357/aa7ece
  * Wernli et al. (2007) Wernli, M., Wiesenfeld, L., Faure, A., & Valiron, P. 2007, A&A, 464, 1147, doi: 10.1051/0004-6361:20066112
  * Wu et al. (2010) Wu, J., Evans, Neal J., I., Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313, doi: 10.1088/0067-0049/188/2/313


# Role of kidney stones in renal pelvis flow

C. R. Constante-Amores  Mathematical Institute
University of Oxford
Oxford OX2 6GG, UK
Email<EMAIL_ADDRESS>L. Kahouadji
Department of Chemical Engineering
Imperial College London
London SW7 2AZ, UK
J. G. Williams
Department of Urology and Pelvic Health
Boston Scientific Corporation
Marlborough
Massachusetts  B. W. Turney  Nuffield Department of Surgical Sciences
University of Oxford
John Radcliffe Hospital
Oxford OX2 6GG, UK
S. Shin  Dept. of Mech. and System Design Engineering
Hongik University
Seoul 04066, Republic of Korea
J. Chergui  Université Paris Saclay
LISN
CNRS
91400 Orsay, France  D. Juric  Université Paris Saclay
LISN
CNRS
91400 Orsay, France;
DAMTP
University of Cambridge
Cambridge CB3 0WA, UK
D. E. Moulton  Mathematical Institute
University of Oxford
Oxford OX2 6GG, UK
S. L. Waters  Mathematical Institute
University of Oxford
Oxford OX2 6GG, UK

###### Abstract

Ureteroscopy is a commonly performed medical procedure to treat stones in the
kidney and ureter using a ureteroscope. Throughout the procedure, saline is
irrigated through the scope to aid visibility and washout debris from stone
fragmentation. The key challenge that this research addresses is to build a
fundamental understanding of the interaction between the kidney stones/stone
fragments and the flow dynamics in the renal pelvis flow.  We examine the
time-dependent flow dynamics inside an idealised renal pelvis in the context
of a surgical procedure for kidney stone removal. Here, we examine the time-
dependent evolution of these vortical flow structures in three-dimensions, and
incorporate the presence of rigid kidney stones. We perform direct numerical
simulations, solving the transient Navier-Stokes equations in a spherical
domain. Our numerical predictions for the flow dynamics in the absence of
stones are validated with available experimental and numerical data, and the
governing parameters and flow regimes are chosen carefully in order to satisfy
several clinical constraints. The results shed light on the crucial role of
flow circulation in the renal cavity and its effect on the trajectories of
rigid stones. We demonstrate that stones can either be washed out of the
cavity along with the fluid, or be trapped in the cavity via their interaction
with vortical flow structures. Additionally, we study the effect of multiple
stones in the flow field within the cavity in terms of the kinetic energy,
entrapped fluid volume, and the clearance rate of a passive tracer modelled
via an advection–diffusion equation. We demonstrate that the flow in the
presence of stones features a higher vorticity production within the cavity
compared with the stone-free cases.

##  1 Introduction

The renal pelvis is a funnel-like cavity inside the kidney which connects
directly to the ureter, leading to the urinary bladder. The kidney’s primary
function – to remove waste and excess fluid from the body – may be disrupted
by the presence of kidney stones in the renal pelvis, resulting in potentially
life-threatening conditions. Kidney stones affect up to $10\%$ of the global
population, and thus constitute a significant healthcare concern [1, 2, 3].
Chemically, the majority ($80-85\%$) of kidney stones are composed of calcium
oxalate and are often due to poor diet/lifestyle and/or dehydration from low
fluid intake. There are many other types of kidney stones including calcium
phosphate, urate, infection stones and metabolic stones (e.g. cystine). Very
small stones, the size of a grain of sand, may form and pass without causing
any symptoms but larger stones which grow to more than $5$ mm [4] can become
stuck in the renal pelvis. Some stones grow to fill the kidney (staghorn
stones) and can measure several centimeters across. Enlarged stones may cause
blockages within the urinary system, often resulting in severe patient
discomfort [5, 6]. Thus, there is a clear necessity for a clinical procedure
to efficiently remove stones from the renal pelvis.

One method for stone removal, ureteroscopy, uses a laser to ablate and
fragment the stone to reduce its size and facilitate natural evacuation
through the ureter. The laser is inserted through the central lumen of the
ureteroscope. To further assist with stone clearance, a flow of saline is
continuously introduced to the renal pelvis through this central lumen, in a
process called irrigation. The fluid exits the renal cavity via a narrow gap
between the ureteroscope and surrounding ureter. This return flow is often
facilitated by the presence of an access sheath, a hollow tube surrounding the
scope shaft, which increases the gap – now between the scope shaft and access
sheath wall – via which saline (and stone fragments) can exit. Figure 1a
provides a schematic representation of the placement of the ureteroscope and
access sheath during ureteroscopy.

Prior to laser-induced fragmentation, kidney stone sizes range between $5$ and
$10$ mm ($66\%$), and $10$ and $15$ mm ($28\%$) [7]. After ablation, the stone
is reduced to a large number of fragments, ranging in radii from $R_{s}=1$
$\mu$m to $R_{s}=1$ mm. The stone particle sizes therefore vary by several
orders of magnitude, and particle-fluid interaction plays a non-trivial role
in the renal-flow dynamics and associated kidney stone removal. Different
particle-fluid interaction behaviours are anticipated depending on the size of
the stone fragments: (i) for small stone sizes ($R_{s}\lesssim 200~{}\mu$m),
stones behave as a passive tracer without affecting the flow dynamics of the
renal cavity, (ii) for large stone sizes ($R_{s}\gtrsim 5000~{}\mu$m), stones
will settle in the cavity under gravity, and (iii) for intermediate stone
sizes ($200~{}\mu\rm m\lesssim R_{s}\lesssim 5000~{}\mu$m), a fully-coupled
interaction between the rigid solid and the flow is expected [8]. We aim to
develop a holistic understanding of the different fluid-structure interaction
mechanisms at play, including the multi-scale effects of varying solid
particle sizes and their interactions with a viscous fluid. An ability to
model stone fragment behaviour – and its dependence on the number, size, and
arrangement of stone particles in the renal pelvis – has the potential to
advise clinical procedure and device design to optimise the efficiency of
ureteroscopic stone removal.

Previous work by [9] and [10] have proposed theoretical models of ureteral
fluid mechanics, considering the ureteroscope as a single conduit. However,
during the ureteroscopic procedure, there is both the forward flow into the
cavity from the scope lumen, as well as the return flow from the kidney to the
ureter through the access sheath. [11] incorporated the effect of the return
flow through the access sheath with a lumped-parameter model to relate flow
rate, kidney pressure, and ureteroscope and access sheath geometry, modelling
the kidney as a linearly compliant material with constant stiffness. [12] also
considered a lumped-parameter model incorporating a more biologically accurate
exponential constitutive law for the kidney compliance, and studied the role
of auxiliary ‘working tools’ – e.g. laser fibres – which are passed through
the central lumen of the ureteroscope and affect resistance to saline flow.

These lumped parameter models were able to predict the renal pressure due to
irrigation and its dependence on flow rate and scope/access sheath geometry,
but the precise nature of the flow within the renal pelvis was not considered.
More recently, [13] modelled the renal pelvis in an idealised two-dimensional,
rectangular geometry, and studied in detail the steady flow patterns and their
subsequent effect on the clearance time of a passive tracer, using a
combination of numerical methods and high-speed imaging techniques. Williams
et al revealed the intricate vortex structures in a two-dimensional cavity,
and demonstrated the connection between the wash-out time (defined as the time
required for $90\%$ of the initial tracer to leave the cavity) and the vortex
characteristics within the cavity, concluding that large vortices combined
with low tracer diffusivity lead to prolonged wash-out times. In a follow-up
study, [14] studied the role of the inflow/outflow channel geometries on the
wash-out time of a passive tracer. Using shape optimisation techniques with an
objective function based on properties of the steady flow field, Williams et
al demonstrated that changing the ureteroscope shape results in a reduction of
the size of vortical flow structures, which in turn leads to a decrease in
wash-out times. Although these studies provide valuable insights into the
relationship between the underlying flow dynamics and the washout time of a
passive tracer, this knowledge is only transferable to ureteroscopy in the
regime where stone fragments are extremely small, and thus unable to influence
flow characteristics. In reality, larger stone pieces are present during
ureteroscopy; thus it is necessary to consider the two-way coupling of flow
and stone dynamics for an accurate understanding of stone removal.
Additionally, there is a need to study full three-dimensional flows in
cavities to unravel the flow physics and the associated stones trajectories.

In this study, for the first time, we interrogate the two-way coupling between
the fluid flow within the kidney and the transport of kidney stones, and
determine the behaviour of stone wash-out during ureteroscopy. Considering the
stones as rigid, non-porous solid objects inside a three-dimensional idealised
spherical cavity representing the renal pelvis that is connected to coaxial
cylinders representing the ureteroscope and sheath, we perform time-dependent
three-dimensional direct numerical simulations using a fictitious domain
method with direct forcing approach to resolve the interaction between the
rigid solids and the flow. To orientate our study, the values of the governing
parameters are chosen to be consistent with clinically-realisable values.

The paper is organised as follows: Section 2 presents the governing equations,
numerical set-up and the validation of the numerical method. Section 3
presents  the results which are focused on fluid-structure interaction in both
the absence and the presence of stones. Finally, a discussion of the results ,
and concluding remarks and future perspective work are summarised in Section 4
and 5, respectively.

##  2 Methods

Figure 1: (a) Schematic representation of the regions of fluid flow during the
ureteroscopic procedure; we highlight the ureteroscope via which the flow
enters the renal pelvis, and a ureteral access sheath in which the irrigation
flow exits the cavity. The flow direction inside of the renal pelvis is shown
with red arrows. (b) Three-dimensional representation of the idealised
geometry of the renal cavity considered in this study with a single rigid
spherical stone at center. The image shows a transversal cut to aid
visualisation of the spherical stone, it also highlights the nondimensional
parameters of the cavity geometry in a three-dimensional Cartesian domain,
$\textbf{x}=(x,y,z)$ which is discretised for our numerical computation by a
$(256)^{3}$ uniform mesh. The center of the spherical cavity corresponds to
$\textbf{x}=(0,0,0)$. Red and blue arrows show the inlet and outlet channels,
respectively. (c) Three-dimensional representation of the $xy$-, $xz$\- and
$zy$-planes that will be used throughout the discussion of the results.

With the aim of studying the flow dynamics of kidney stone removal, we perform
direct numerical simulations of the transient, incompressible Navier-Stokes
equations in a three-dimensional domain focusing, for this present study, on
the case of the absence of gravity. We assume spherical rigid non-porous
stones which move freely in the computational domain owing to their
interaction with the cavity flow. Figure 1b shows a schematic representation
of the cavity considered in this study. To simplify the tremendous complexity
of the renal cavity, we consider a stiff, non-deformable sphere of radius
$R_{c}$, which is connected to a ureteroscope and a ureteral access sheath.
The ureteroscope and sheath are considered to be coaxial cylinders, in which
flow enters the cavity via the inner scope channel of radius $a$ and exits the
cavity through the gap (width $h$) between the access sheath and scope shaft
of radius $b$. Renal stones are characterised by an irregular shape, but for
the sake of simplicity, we have considered them as perfect spherical non-
porous rigid solids of radius $R_{s}$. They are placed initially in the centre
of the cavity, if not explicitly stated otherwise.

Following the immersed boundary approach [15], we use a one-fluid formulation
to couple the fluid-structure interaction (FSI) problem. Within this
formulation, parameters $\rho_{l}$ and $\mu_{l}$ are the density and viscosity
of the fluid, respectively, and $\rho_{s}$ is the density of the solid.

We consider dimensionless variables

$\quad\tilde{\mathbf{x}}=\frac{\mathbf{x}}{a},\quad\tilde{t}=\frac{t}{a/U},\quad\tilde{\textbf{u}}=\frac{\textbf{u}}{U},\quad\tilde{p}=\frac{p}{\mu_{l}U/a},$
(1)

where t, u, and $p$ represent time, velocity, and pressure, respectively, and
tildes denote dimensionless quantities. We nondimensionalise spatial
coordinates $\mathbf{x}=(x,y,z)$ with respect to the radius of the scope lumen
$a$, $\mathbf{u}=(u,v,w)$ with respect to the average irrigation flow velocity
$U$, and time with respect to the timescale of the flow, $a/U$. The
dimensionless time-dependent, incompressible Navier-Stokes equations for a
Newtonian viscous fluid are thus

$\nabla\cdot\tilde{\textbf{u}}=0,$ (2)

$\tilde{\rho}\left(\frac{\partial\tilde{\textbf{u}}}{\partial\tilde{t}}+\tilde{\textbf{u}}\cdot\nabla\tilde{\textbf{u}}\right)+\nabla\tilde{p}=\frac{1}{Re}\nabla^{2}\tilde{\textbf{u}},$
(3)

where
$\tilde{\rho}=\rho_{s}/\rho_{l}+\left(1-\rho_{s}/\rho_{l}\right)\mathcal{H}\left(\tilde{\textbf{x}},\tilde{t}\right)$
where $\mathcal{H}\left(\tilde{\textbf{x}},\tilde{t}\right)$ is zero in the
solid stones and unity in the fluid. The dimensionless Reynolds number,
$Re=\rho_{l}Ua/\mu_{l}$, in equation 3 relates the inertial to viscous forces.
All the variables appearing in the equations and boundary conditions are
rendered dimensionless using the aforementioned scalings, unless stated
otherwise.

A no-slip condition for the immersed solid-fluid boundary, $\Gamma$, is
enforced, taking into account rigid body motion due to translation and
rotation of the stones

${\tilde{\textbf{u}}}_{\Gamma}=\tilde{\textbf{V}}_{s}+\tilde{\Omega}\times{\tilde{\textbf{r}}}_{c}$
(4)

where ${\tilde{\bf V}}_{s}$ represents the translation velocity of the solid
(non-dimensionalised with respect to $U$), ${\tilde{\Omega}}$ stands for its
angular velocity vector (non-dimensionalised with respect to $U/a$) and
${\tilde{\bf r}}_{c}$ the position vector of the centroid of the solid (non-
dimensionalised with respect to $a$). ${\tilde{\bf V}}_{s}$ and
${\tilde{\Omega}}$ can be computed by averaging over the solid region.  For
the collision modeling we use the the well-known impulse-response model (for
further information we refer the reader to [16]).

###  2.1 Initial and domain boundary conditions

Henceforth we drop tildes, and all quantities are now dimensionless. The
simulations are initialised with fluid and stones at rest in the absence of
gravity. Solutions are sought subject to Poiseuille inlet flow. Thus a
parabolic profile for the inflow channel is specified (i.e., $u=2(1-r^{2})$),
where $r=\sqrt{\left(y-y_{o}\right)^{2}+\left(z-z_{o}\right)^{2}}$, where
$y_{o}$ and $z_{o}$ are the (dimensionless) coordinates of the centre of the
inflow channel. As outflow conditions, we impose Neumann boundary conditions
for the velocity in the streamwise direction (i.e., perpendicular to the
outflow face), and the spanwise velocity components are zero. A Dirichlet
boundary condition is imposed for the pressure. This condition is a result of
imposition of zero normal stress at the outlet together with the assumption of
fully-developed streamwise flow. Additionally, we impose no-slip boundary
conditions on the cavity walls.

###  2.2 Numerical method and validation

The Navier-Stokes equations are solved using classic finite volume techniques
applied on a uniform staggered grid [17]. A multigrid iterative method is used
for solving the elliptic pressure Poisson equation that arises in the
projection method [18] when enforcing the incompressibility condition
(Equation 2). With respect to the spatial derivatives, standard centered-
difference discretisations are used, except for the nonlinear term, which
makes use of a second-order essentially nonoscillatory (ENO) scheme [19, 20].
The Lagrangian motion of the solid centroid (i.e., $\rm d\textbf{r}_{c}/\rm
d\textbf{t}=\textbf{V}_{s}$) is computed using a second-order in time Runga-
Kutta algorithm. The boundary condition Eq. 4 is enforced using the Direct
Forcing approach of [21] where the desired value of velocity is imposed
directly on the solid/fluid boundary. In addition to the extensive analysis
and validation provided in [21] we validate our own implementation of the
Direct Forcing method with results shown in the Appendix. In this section we
have provided only a brief synopsis of the numerical approach, however a
detailed description of the immersed boundary, fictitious domain and direct
forcing methods implemented here can be found in [22, 23, 21, 24].  The code
is wholly written by the authors in Fortran 2008 and uses a domain
decomposition strategy for parallelization with MPI.

The validity of the numerical method was benchmarked against the experimental
and numerical work of [13] for a two-dimensional cavity with excellent
agreement (more information in the Appendix). In terms of the mesh
characteristics, the computational domain is composed of a cubic uniform
Cartesian grid, and the inlet extends in the $x-$direction. With respect to
the resolution, we have ensured that the presented results are mesh-
independent, and therefore for a uniform mesh resolution of $(256)^{3}$ (i.e.,
$a/\Delta{\bf x}=6.4$, where $\Delta{\bf x}$ stands for the cell-size), the
results do not significantly change with decreasing cell size (see Appendix
for more details). Additionally, extensive mesh studies for turbulent two-
phase jets and surface-tension-driven phenomena using the same numerical
method can be found in [25, 26].

###  2.3 Physical parameters

The chosen flow parameters for this study are consistent with values used
during ureteroscopy. In the clinic, irrigation flow rates range from $50$
mL/min to $200$ mL/min, although introduction of a laser fibre (for stone
fragmentation) to the scope channel may reduce the flow rate by about half.
The injected fluid is considered to be water, with constant physical
properties (i.e., $\rho_{l}=1000$ kg/m3 and $\mu_{l}=10^{-3}$ Pa.s).
Therefore, the Reynolds number lies in the range $220<\text{Re}<2000$.

In the ureteroscopy procedure, the stone is first fragmented by the laser to
enable expulsion through the ureter. We consider fragmented stone sizes
ranging between $200\mu\hbox{m}<R_{s}<1000\mu$ m where there will be a fully-
coupled interaction between the flow and stone trajectory. We take the stone
density as $\rho_{s}=1900$ kg/m3 [27]. Gravity will play a minor role in the
flow dynamics, as indicated by the typical Froude number $Fr=U/\sqrt{ga}\sim
O(10^{2})$ (where $g$ stands for the gravitational constant); thus, the effect
of gravity is not considered in this study.

The values for the cavity dimensions are identical to the previous work
performed by [14], and consequently, $a=6\times 10^{-4}$ m, $b=2.5701\times
10^{-3}$ m, and $h=1.1598\times 10^{-3}$ m. The size of the cavity was chosen
to be $R_{c}=10^{-2}$ m, which could be seen as a representation of a small
region of the renal pelvis. The non-dimensional length of the inlet channel is
set to $5$.  We note that the values for $b$ and $h$ are slightly larger than
the typical range of scope and access sheath dimensions for ureteroscopy
procedures. These were chosen to agree with the values used in previous
physical and numerical experiments of flow in rectangular domains, simulating
the dynamics of ureteroscopic flow in the renal pelvis [13, 14]. This
agreement allows us to investigate the effect of scaling up from two to three
dimensions, and although exact flow solutions will depend on the values of $b$
and $h$, we anticipate reported trends in flow structure and the influence of
flow on stone dynamics will be similar for smaller values of $b$ and $h$.
Values for $b$ and $h$ more representative of ureteroscope and access sheath
dimensions were chosen for the simulations in the Discussion (see figure 7),
demonstrating the potential for higher fidelity ureteroscopy simulations using
the techniques outlined in this manuscript.

##  3 Results

In this section we first consider the effect of Reynolds number on three-
dimensional flow in the absence of stones, before subsequently analysing the
effects of both single and multiple stones. It is worth noting that after
$t=5000$, for the lowest Reynolds number ($\text{Re}=50$), no further changes
in the shape of the streamlines or coherent structures were observed; however,
for the higher Reynolds numbers ($\text{Re}=250$, $500$), flow evolution
continued and a steady-state was not reached.

###  3.1 Three-dimensional flow dynamics in the absence of rigid stones

Figure 2:  Time-dependent flow dynamics inside of the idealised renal cavity
for $\text{Re}=(50,250,500)$. Panels (a)-(c) correspond to early flow dynamics
for $Re=500$ at $\text{t}=(27,138,194)$ through two-dimensional
representations of the predicted streamlines in a half of the $xy-$plane.
Panels (d)-(f) correspond to two-dimensional representations of the predicted
streamlines in a $xy-$plane together with the azimuthal vorticity
$\omega_{\theta}$ at $t=5000$ for $\text{Re}=(50,250,500)$ in columns one to
three, respectively; Panel (g) shows the spatio-temporal evolution of the
three-dimensional coherent structures by representing the entrapped fluid
volume when $Q>0$, $V_{v}$ normalised by the volume of the cavity, $V_{c}$,
for $\text{Re}=(50,250,500)$. Additionally, we refer the reader to
Supplementary Material ‘Animation-Fig2.avi’ to see a video from panel (e).

First, we focus on the transient dynamics of the jet. Panels (a)-(c) of Figure
2 show the instantaneous streamlines at early times for $\text{Re}=500$. For
visualization purposes, we display a two-dimensional streamlines onto the
$xy$-plane (for sake of space only the half of the plane is shown). At early
times of the injection, the jet dynamics are characterised by an axisymmetric
behaviour. The homogeneous jet develops a ‘primary-vortex’ (labelled ‘V-1’ in
figure 2a) as a result of the detachment of the velocity boundary layer from
the inlet nozzle, and its consequent roll-up in the quiescent medium. A
similar physical mechanism for the formation of the leading vortex has been
reported previously by [28, 29, 26]. The primary vortex grows over time owing
to the shear-driven interaction between the high-velocity discharge of the
fluid into the quiescent medium. The primary-vortex elongates as it moves
downstream until it impacts against the surface of the cavity wall resulting
in the outward radial spread of the head-vortex (see figure 2b). The formation
of a stagnation point (labelled ‘SP-1’ in figure 2b) near the cavity wall
resulting in flow separation is observed as a result of the impact of the jet
against the wall. For the high Reynolds numbers cases, we observe that the
radial expansion of the jet moving outward leads to an adverse pressure
gradient until a point in which the jet-induced momentum can no longer
overcome the frictional effects of the cavity wall resulting in the formation
of a stagnation point (labelled ‘SP-2’ in figure 2c) and flow separation. This
phenomenon leads to the detachment of a secondary-vortex, labelled ‘V-2’ in
figure 2c (in agreement with [30, 31, 32, 33]). As the secondary vortex grows
over time, a mutual-induction of V-1 and V-2 is predicted (see panel (ii) of
figure 2g): a velocity induction by the primary-vortex over the secondary-
vortex causes V-2 to go through V-1’s centre. This mutual-induction mechanism
is observed during the entire flow dynamics at high $Re$ numbers.

Panels (d)-(f) of figure 2 highlight the instantaneous streamlines at $t=5000$
for $\text{Re}=50$, $\text{Re}=250$, and $\text{Re}=500$, respectively. For
visualization purposes, we have also shown a two-dimensional streamlines onto
the $xy$-plane. The azimuthal component of the vorticity field
$\boldsymbol{\omega}=\triangledown\times\textbf{u}$ (e.g. the circumferential
direction around the axis of the jet), is shown as a contour field underneath
the streamlines. For the smallest Reynolds number, the jet ejection results in
an axisymmetric flow with no formation of secondary vortices. A closer
inspection of the same figure shows that the areas of closed streamlines
(i.e., the main vortical structure) are surrounded by streamlines which show a
direct path between the inflow and outflow channels with no recirculation
regions; these outcomes are in agreement with [13, 14]. Increasing the
Reynolds number results in the breaking of the symmetry of the flow patterns,
and subsequently, enhances the formation of complex flow patterns inside of
the cavity. The primary vortex is no longer aligned axisymmetrically in the
cavity as displayed in figures 2(e)-(f). At higher $Re$, there is a reduction
in the number of direct paths from the inlet to the outlet channels owing to
the increase of the inertia within the cavity. Additionally, by inspection of
the azimuthal vorticity in the $xy$-plane, we observe that the highest
vorticity generation coincides with the velocity boundary layer between the
injected-jet and the surrounding fluid, which leads to strong tangential flow
with respect to the initial quiescent medium. Attention is now turned to the
value of vorticity (red/blue contours in Figures 2d-f) which increases as the
Reynolds number increases. For $Re=50$, the azimuthal vorticity component
dominates over its streamwise counterpart (not shown) because of the low
inertia; however, as the Reynolds number increases, the streamwise component
becomes more prominent, becoming of the same order of magnitude with respect
to its azimuthal counterpart, explaining the disruption of the axisymmetric
behaviour of the coherent structures as explained below.

Figure 2g displays the time-dependent flow dynamics by plotting the volume
ratio $V_{v}/V_{c}$ occupied by the vortical structures (i.e., recirculation
regions $V_{v}$), predicted by the Q-criterion inside of the cavity (with
volume $V_{c}$). The Q-criterion measures the dominance of vorticity $\omega$
over that of strain $\bf{s}$, i.e.,
$Q=1/2(|\boldsymbol{\omega}|^{2}-|\bf{s}|^{2})$ [34]. The flow regions with a
positive Q-criterion value are defined as the recirculation zone (i.e.,
vortex-dominated regions). Thus, the volume representation of figure 2g
denotes flow regions in which $Q>0$. A similar approach was used by [13, 14]
to quantify the vortex-regions for a two-dimensional system (i.e., det
$\triangledown\textbf{u}>0$).

For $Re=50$, the axisymmetric primary-vortex defines the global flow dynamics
and the vortex-dominated region reaches a steady state resulting in a constant
entrapped recirculating fluid within the cavity (see figure 2g). At higher
$Re$, the Q-criterion predicts larger vortex-dominated regions inside of the
cavity with a linear increase of the entrapped-fluid volume at early times of
the simulation (see figure 2g). The increase in the rich dynamics of the flow
structures is evident by observing the spatial development of the coherent
structures (i.e., regions of dominant vorticity). At short times, the jet
impacts the cavity wall (see panel (i) of figure 2g), and eventually results
in the formation of secondary vortices which are advected towards the centre
of the primary vortex (see panel (ii)) while the primary vortex grows over
time. Once the primary vortex reaches its maximum size, the entire flow almost
behaves as a steady state (see panels (iii)-(iv)).

To see the impact of the vortical structures on the trajectories of rigid
stones, in the next section we turn our attention to the effect of flow
dynamics on the motion of stones within the cavity.

###  3.2 Effect of flow dynamics on rigid kidney stones

Figure 3: Trajectories of single rigid stones up to $t=5000$. Panels (a) and
(b) correspond to a three-dimensional representation of the trajectories for
$\text{Re}=250$ and $\text{Re}=500$, respectively. Panels (c)-(h) correspond
to the projection of the trajectory onto the $xy$-, $xz$\- and $yz$-plane for
$\text{Re}=250$ (middle panels) and $\text{Re}=500$ (bottom panels),
respectively. The red circles correspond to the limits of the cavity, where
the center of the spherical cavity corresponds to $\textbf{x}=(0,0,0)$.
Figure 4: Trajectories of multiple rigid stones (i.e., ten) of size
$R_{s}=0.833$ up to $t=5000$. Panels (a)-(c) show the initial location of the
stationary stones, the location of the stones at $t=5000$ for $\text{Re}=250$,
and $\text{Re}=500$, respectively. Additionally, a two-dimensional
representation of the velocity streamlines in the xy-plane have been added to
visualise the flow field. Some stones have left the renal cavity and are found
in the ureteral access sheath, meanwhile others are trapped in the primary
vortex. Panels (d)-(i) correspond to the projection of the trajectory for four
selected stones onto the $xz$-, $xy$\- and $zy$-plane for $\text{Re}=250$
(middle panels) and $\text{Re}=500$ (bottom panels), respectively. The colour
of the trajectories correspond to the coloured stones of panels (a)-(c). See
Supplementary Material ‘Animation-Fig4.avi’.

In this section, we consider the effect of renal pelvis flows on the dynamics
of a single rigid stone, of varying size, placed initially in the centre of
the cavity (e.g., as displayed in Figure 1b). We vary the stone size between
$0.416<R_{s}<1.66$ (nondimensional values). Additionally, we only report
results for $\text{Re}=250$ and $\text{Re}=500$ as these correspond to larger
values of inertia prompting richer flow dynamics within the cavity.

Figure 3 shows the trajectories of the stones over time by tracking their
centroid as both three-dimensional trajectories (Figure 3a,b) and projections
onto the $xy$-, $xz$\- and $zy$-planes (Figures 3c-h). At early stages of the
simulation, the rigid stones are dragged towards the back of the cavity – as a
result of their interaction with the head-vortex – and eventually collide
against the cavity wall. The stones subsequently settle into a circulatory
motion owing to their interaction with the local vorticity field.

For $\text{Re}=250$, the stones become trapped in the primary vortex, leading
to their circulation in the vortex-dominated region (see Figure 3a).
Inspection of the stones’ movement in each projection plane (Figure 3c-d)
shows that the stones are trapped in the primary vortex (we refer the reader
to Supplementary material ‘Animation-Fig3.avi’ where we show the transient
dynamics of the stone of size $R_{s}=0.833$ together with the Q-criterion).
Interestingly, the radii of the trajectories grow over time, so it is expected
that eventually the stones would be displaced to a non vortex-dominated region
(e.g., near the cavity exit), where they will washout together with the fluid.
Nonetheless, for the considered non-dimensional time of $t<5000$, the rigid
stones remain trapped in the primary vortex during the entire simulation.

Still considering $Re=250$, now we compare the effect of stone size on the
flow pattern and particle trajectory. For the smallest stone size (e.g.,
$R_{s}=0.416$), the stone follows the primary vortex, and its circulation
radius grows over time. The simulations predict smaller circulation radii for
the largest stone size and stronger entrapment in the vortical structure
because the stones start stationary, so a larger stone would require a larger
force to attain the same velocity.

For $\text{Re}=500$, richer flow dynamics are predicted owing to the higher
fluid inertia, which in turn leads to more complex stone trajectories (see
Figures 3b,f-h). At early time in the simulations, the rigid stones similarly
enter into a circulatory motion due to their interaction with the primary
vortex; however, the higher fluid inertia induces the stones to leave the
primary vortical structure, which eventually results in the smallest stone,
i.e., $R_{s}=0.416$, being flushed out of the cavity by the flow (at
$t=3507$). Therefore, for the first time (to the best of the authors’
knowledge), we have provided a simulation of stones being washed out of the
cavity along with the fluid during the ureteroscopy procedure.

Above we have considered the effect of flow dynamics on single rigid stones
placed in the centre of the cavity. To investigate how the initial position of
the rigid stones affects the fate of the stone, we consider ten stones within
the cavity, initially placed on random grid points as indicated in figure 4a.
Panels of figure 4(b)-(i) show snapshots of the flow pattern and trajectory of
stones over time via projections onto the $xz$-, $xy$\- and $zy$-planes. We
found that stones initially positioned near the centre of the cavity are
pushed towards the outflow channel as the primary vortex grows in size, and
eventually advect directly out of the cavity. Conversely, stones initially
placed further downstream of the injection point are entrapped in a
circulatory motion owing to their interaction with the primary-vortex, and
subsequently, the stone’s dynamics are similar to those described above for a
single stone that does not leave the domain. Therefore, some stones leave the
cavity into the ‘ureteral access sheath’, whereas others are still entrapped
in the main vortical structure (see Supplementary Material ‘Animation-
Fig4.avi’).

We acknowledge that to better understand the interaction between kidney stones
and the flow cavity, it would be interesting to use Lagrangian coherent
structures (LCS) based on particle trajectories to naturally find borderlines
in the flow that partition different regions with different trajectory
behavior. This will connect the particle initial and final locations. However,
this is out of the scope of the manuscript, and should be addressed in futures
studies.

###  3.3 Effect of rigid kidney stones on flow dynamics

In the previous section, we considered the effect of the flow field on the
trajectories of small stones which are not likely to significantly affect the
flow dynamics. In this section, we will consider the influence of larger
kidney stones on the flow dynamics by exploring the fluid kinetic energy, the
clearance rate of a passive tracer and the entrapped fluid volume when $Q>0$,
depending on the number of stones, $N_{p}$. On this basis, we consider a
passive tracer of concentration $C(\textbf{x},t)$ within the cavity which is
passively advected by the flow field. The nondimensional advection-diffusion
equation for the tracer concentration is expressed as

$\frac{\partial C}{\partial t}+\textbf{u}\cdot\nabla
C=\frac{1}{Pe}\nabla^{2}C,$ (5)

where $Pe=Ua/\mathcal{D}$ represents the ratio of convective to diffusive
time-scales (e.g., $\mathcal{D}$ stands for the tracer diffusion coefficient).
Following [14], we assume that the initial tracer concentration is evenly
distributed inside of the cavity, and zero elsewhere. We have selected
$Pe=500$, guided by the previous work from [13]. We assume no-flux boundary
conditions on the walls of the cavity (i,e., $\partial
C/\partial\mathbf{n}=0$), and Neumann boundary conditions for $C$ on the
outlet channels which allow the passive tracer to leave the computational
domain. To quantify the effect of the rigid stones on the flow dynamics, an
approach proposed by [13, 14] was employed in this work which measures the
clearance rate at which the passive tracer leaves the cavity. Thus, the
nondimensional loss at time $t$ is expressed as

$C_{\text{l}oss}=\frac{\int_{\mathcal{V}}\left[C(\mathbf{x},0)-C(\mathbf{x},t)\right]d\mathcal{V}}{\int_{\mathcal{V}}\left[C(\mathbf{x},0)\right]d\mathcal{V}}$
(6)

where $\mathcal{V}$ refers to the reduced volume of the cavity based on the
number of stones.

In Figure 5, we show the influence of stones on the fluid kinetic energy,
defined as $E_{k}=\int_{\mathcal{V}}(\mathbf{u}^{2}/2)d\mathcal{V}$, the
entrapped fluid volume for $R_{s}=2$ 111$R_{s}=2$ was chosen as it is larger
than the size of the outflow channel, $h$, and therefore the stones will
remain within the cavity and $Pe=500$, and the clearance rate of the passive
tracer. Inspection of the kinetic energy $E_{k}$ in figure 5a,d reveals a
complex relationship with the number of stones, $N_{p}$. For all values of
$N_{p}$, at early times ($t<35$), $E_{k}$ grows independently of the number of
stones. Then the kinetic energy plateaus, but for $N_{p}>0$, we observe
irregular oscillatory behaviour with increased amplitude for higher $N_{p}$.
The fluid must do work to accelerate the stones, taking some of the kinetic
energy of the fluid, but since the motion of the stones is quite irregular,
they experience numerous accelerations and decelerations. Hence, the more
stones there are, the more noisy the kinetic energy.

Panels (b), (e) of figure 5 plot the ratio of entrapped fluid volume (i.e.
with $Q>0$) to total fluid volume. At early times, a linear growth of the
entrapped fluid volume is observed. As mentioned above there is a direct link
between the entrapped volume and vorticity production. [35] showed that
vorticity is produced at boundaries; thus, the addition of stones in the
cavity implies the addition of regions for production of vorticity. However,
the stones may also disrupt vortical structures, depending on their location
within the cavity. This is demonstrated in figure 6, which shows a snapshot of
the instantaneous vorticity field (in a fixed plane) for three scenarios: (a)
the absence of stones, (b) stones with position held fixed, and (c) freely
moving stones. Comparing (b) and (c) to (a), we see the two-sided impact of
stones, both in producing partially detached vortices, and in disrupting
vortical structures; similar findings were reported by [36]. Due to this
complex interplay, a cavity with more stones may have increased entrapped
volume at some times, and decreased entrapped volume at other times, a noisy
and unpredictable relationship as demonstrated in figure 5 (b), (e).

We now turn our attention to the clearance rate plots (panels (c), (f) of
figure 5), which show the percentage of tracer that has exited the cavity over
time. It has been previously demonstrated [14] that the reduction of vortex
structures within the cavity results in a reduction in the time needed for the
clearance of the tracer, while clearance time depends non-monotonically on the
fluid kinetic energy. Given the complex relationships between stone presence
and both kinetic energy and vorticity, it is perhaps not surprising then that
the clearance rate does not exhibit a clear trend. The highest number of
stones does perform the worst in terms of tracer clear out, it is by a small
margin, and the smaller number of stones performs slightly better than the
base case of no stones.

Figure 5:  The effect of rigid stones on the flow field as a function of the
number, $N_{p}$, for $R_{s}=2$ and $Pe=500$. Kinetic energy, entrapped fluid
volume when $Q>0$, and the clearance rate of a passive tracer are shown in
columns one to three, respectively. The top and bottom panels correspond to
$\text{Re}=250$ and $\text{Re}=500$, respectively.  Figure 6:  Two-dimensional
representation of the vorticity field in the $xz$-plane ($y=2$) in the absence
of stones, fixed stone location and stones which move freely in the cavity,
corresponding to panels $(a-c)$, respectively. The conditions for the
simulations are $Re=250$ and $R_{s}=2$ and $N_{p}=18$ (for cases b and c). The
colour represents the vorticity magnitude $|\omega|$, where appropriate scales
are shown in the panel. Irregular vortex production can be seen around the
rigid stones which move freely.

##  4 Discussion

A fundamental challenge in optimising protocols in ureteroscopy procedures is
understanding the flow of irrigation fluid within the kidney, the movement of
kidney stone particles and/or dust, and how larger stone fragments may impact
fluid flow. Stone destruction via laser lithotripsy can create a range of
sizes of stone fragments or dust [37], which, if not removed efficiently via
irrigation fluid, can obscure the operating clinician’s field-of-view [38,
39]. Previous theoretical studies have uncovered the interplay between fluid
structure, in particular the presence of vortical regions, and dust washout.
These studies uncovered qualitative features and enabled large parameter
exploration with minimal computational cost, but were restricted to a highly
idealised two-dimensional geometry. To bridge the gap towards clinical
relevance, in this paper we have extended this framework by incorporating
discrete stone particles of finite size that both move with and modulate the
flow, solved in a more realistic three-dimensional setting. The chosen
parameters for the flow dynamics and the stone sizes are consistent with
clinically-realisable systems, and the numerical framework was validated
against the experimental and numerical work presented by [13] for the two-
dimensional geometry.

In contrast to the aforementioned 2D studies, the vortex structure we have
uncovered in the present study are more complex both in terms of spatial
structure, interaction, and time-dependence. Specifically, we observed that
the injection of a fluid jet via a modelled ureteroscope nozzle results in the
formation of a primary vortex owing to the detachment of the velocity boundary
layer from the nozzle. The primary vortex grows over time, and its
axisymmetric shape becomes more disrupted with increasing Reynolds number. In
the absence of stones, we have quantified the temporal entrapped fluid volume
by identifying the recirculation zones through the Q-criterion, i.e. $Q>0$,
integrated over the cavity domain. At low inertia, a steady state is reached,
and subsequently a constant volume for the entrapped-fluid is predicted. At
increasing inertia, a steady state is not reached owing to the formation of a
secondary vortex from the cavity walls. Mutual induction of the primary and
secondary vortices characterise the complex flow dynamics observed in the
system. The increase of fluid inertia at larger Reynolds numbers also results
in an increase in the streamwise vorticity component with respect to its
azimuthal vorticity counterpart. Thus, it is responsible for the loss of the
axisymmetric behaviour observed in the small Reynolds number case.

With rigid stones introduced in the cavity flow, two different behaviours were
observed: either they become trapped in a circulatory motion in the primary
vortical structure or they are flushed out of the cavity together with the
fluid through the ureteral access sheath. Both regimes were found to coexist
in the studied cases and the outcome of the stones depended on the initial
location of the stationary stone, and the importance of inertia. More fluid
inertia allows for higher probability of a stone to leave the cavity within
the same temporal frame. We have also studied the effect of multiple stones in
the flow dynamics by exploring the fluid kinetic energy, entrapped fluid
volume and the clearance rate of a passive tracer. We have shown that the
solid boundaries, due to the the presence of stones, cause an increase in the
vorticity production, which results in richer flow dynamics within the cavity.

We have simulated for the first time (to the best of our knowledge) the two
possible outcomes for the kidney stone trajectories during ureteroscopy –
trapping or wash-out – and this can be considered to be a significant step
forward in the understanding of the dynamics of kidney stones in the renal
cavity flow during the surgical procedure. The complex behaviour outlined
above highlights the need for computational fluid mechanics tools in
ureteroscopy, and demonstrates a strong potential for optimising driving flow
conditions to promote stone removal. In practice, this may best be achieved
with time-dependent driving flow, a complication we have not considered here.

As well as incorporating time-dependent inlet conditions, the approach we
described in this manuscript may naturally be extended to a systematic follow-
up analysis predicting the flow characteristics and stone motions in an actual
renal cavity geometry, potentially incorporating full three-dimensional time-
dependent simulations and coupling of moving solids, buoyancy effects, and
heat-transfer considerations. We have restricted attention to an idealised
spherical cavity, and an actual renal pelvis would provide a more
representative interaction of stones with the fluid. Figure 7 shows a
preliminary simulation from an actual renal pelvis extracted from imaging in a
patient, in which we have selected a branch of the complex renal pelvis and
introduced the ureteroscope  together with rigid stones; results at long times
show that the stones are trapped in one of the branches of the renal cavity.
This avenue of research could lead to the next generation of in-silico models
that are patient-specific, and translate our understanding ‘from bench-to-
bedside’.

The present work is limited by the assumption of rigid cavity walls. This
assumption is realistic for the inflow and outflow channels; however, the
walls of an actual renal cavity are characterised by deformability that is
inherent to any true physiological structure. Deformability of the renal
pelvis would certainly affect the flow structure inside of the cavity,
potentially resulting in more complex flow dynamics, and constituting a
fruitful area of future research. The numerical method used in this study may
naturally be extended to include deformability of the renal cavity and
inclusion of these effects constitutes a realistic short-term follow-up to the
present study. We have limited our study to perfectly spherical kidney stones,
but their inherent complex three-dimensional shape will play a major role in
the fluid-structure interaction. We have only studied equal-sized rigid stones
in the cavity, but in surgical procedures, stone destruction via laser
lithotripsy leads to obliteration of the stone into a large range of sizes,
and subsequently, future research should consider the potential effects of
uneven-sized stones. Particle-particle interaction should be taken into
account in future studies as it would become more important as the particles
become larger and more numerous. Additionally, future studies would also
quantify the enstrophy within the renal cavity as a function of the time, as a
measure of the rotational energy of the fluid flow.

Figure 7:  The left panel shows a three-dimensional renal pelvis from a
patient which has been reconstructed using CT scans. The right panel shows a
snapshot at long times for a preliminary simulation to demonstrate the
capabilities of the methodology used in this research. We have extracted one
of the branches of the renal pelvis (highlighted in the left panel), and have
introduced a ureteroscope. The close interaction of flow cavity and the stiff
stones result in the stones being trapped in one of the renal pelvis branches.
Computed velocity is shown to visualise the flow dynamics inside of the
cavity.

##  5 Conclusions

This work has focused on understanding the fluid mechanics of ureteroscope
irrigation through three-dimensional numerical simulations. This study is
limited to stiff, non-deformable solids for both the stones and the renal
pelvis cavity. The numerical predictions show the close interplay between the
local vorticity and the stones, and their interaction decides the outcome of
the stones: either they become trapped in a circulatory motion in the primary
vortical structure or they are flushed out of the cavity together with the
fluid through the ureteral access sheath. We have also studied the efficacy of
debris clearance influenced by the presence of different numbers of kidney
stones.

## Acknowledgements

We acknowledge HPC facilities and computing resources provided by the Research
Computing Service (RCS) of Imperial College London and support through
computing time at the Institut du Developpement et des Ressources en
Informatique Scientifique (IDRIS) of the Centre National de la Recherche
Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul
Intensif) Grant No. 2022A0122B06721. L.K. acknowledges the financial support
from the Engineering and Physical Sciences Research Council, United Kingdom,
through the EPSRC PREMIERE (EP/T000414/1) Programme Grant. The numerical
simulations were performed with code BLUE ([40]) and the visualisations have
been generated using ParaView.

The authors declare no conflict of interest.

## References

  * [1] Romero, V., Akpinar, H., and Assimos, D. G., 2010.  “Kidney stones: a global picture of prevalence, incidence, and associated risk factors”.  Reviews in urology, 12(2-3), pp. e86–e96.
  * [2] Cortes, J. A. S., Motamedinia, P., and Gupta, M., 2011.  “Update on technological and selection factors influencing shockwave lithotripsy of renal stones in adults and children”.  Current Opinion in Urology, 21(2).
  * [3] Turney, B. W., Reynard, J. M., Noble, J. G., and Keoghane, S. R., 2012.  “Trends in urological stone disease”.  BJU International, 109(7), pp. 1082–1087.
  * [4] Khan, S. R., Pearle, M. S., Robertson, W. G., G. Gambaro, B. K. C., Doizi, S., Traxer, O., and Tiselius, H. G., 2016.  “Kidney stones”.  Nature Reviews Disease Primers, 2(16008).
  * [5] Segura, J. W., Patterson, D. E., LeRoy, A. J., Williams, H. J., Barrett, D. M., Benson, R. C., May, G. R., and Bender, C. E., 1985.  “Percutaneous removal of kidney stones: Review of 1,000 cases”.  Journal of Urology, 134(6), pp. 1077–1081.
  * [6] Miller, N. L., and Lingeman, J. E., 2007.  “Management of kidney stones”.  BMJ, 334(7591), pp. 468–472.
  * [7] Cui, H., Thomee, E., Noble, J. G., Reynard, J. M., and Turney, B. W., 2013.  “Efficacy of the lithotripsy in treating lower pole renal stones”.  Urolithiasis, 41(3), pp. 231–234.
  * [8] Keller, E., De Coninck, V., Doizi, S., Daudon, M., and Traxer, O., 2021.  “What is the exact definition of stone dust? an in vitro evaluation”.  World J. Urol., 39, 01.
  * [9] Lykoudis, P. S., and Roos, R., 1970.  “The fluid mechanics of the ureter from a lubrication theory point of view”.  J. Fluid Mech., 43(4), p. 661–674.
  * [10] Yang, P. J., Pham, J., Choo, J., and Hu, D. L., 2014.  “Duration of urination does not change with body size”.  PNAS, 111(33), pp. 11932–11937.
  * [11] Oratis, A. T., Subasic, J. J., Hernandez, N., Bird, J. C., and Eisner, B. H., 2018\.  “A simple fluid dynamic model of renal pelvis pressures during ureteroscopic kidney stone treatment”.  PLOS ONE, 13(11), 11, pp. 1–13.
  * [12] Williams, J. G., Rouse, L., Turney, B. W., Waters, S. L., and Moulton, D. E., 2020\.  “A lumped-parameter model for kidney pressure during stone removal”.  IMA J Appl Math, 85(5), 07, pp. 703–723.
  * [13] Williams, J. G., Castrejón-Pita, A. A., Turney, B. W., Farrell, P. E., Tavener, S. J., Moulton, D. E., and Waters, S. L., 2020.  “Cavity flow characteristics and applications to kidney stone removal”.  J. Fluid Mech., 902, p. A16.
  * [14] Williams, J., Wechsung, F., Turney, B., Waters, S., and Moulton, D., 2021.  “Shape optimisation for faster washout in recirculating flows”.  J. Fluid Mech., 914, p. A37.
  * [15] Peskin, C. S., 1977.  “Numerical analysis of blood flow in the heart”.  J. Comput. Phys, 25(3), pp. 220 – 252.
  * [16] Baraff, D., 1997.  “An introduction to physically based modeling: rigid body simulation ii—nonpenetration constraints”.  SIGGRAPH course notes, pp. D31–D68.
  * [17] Harlow, F., and Welch, J., 1965.  “Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface”.  Phys. Fluids, 8, pp. 1958–1988.
  * [18] Chorin, A. J., 1968.  “Numerical solution of the Navier-Stokes equations”.  Maths. Comp., 22(104), pp. 745–745.
  * [19] Shu, C., and Osher, S., 1989.  “Efficient implementation of essentially non-oscillatory shock-capturing schemes, ii”.  J. Comput. Phys., 83(1), pp. 32–78.
  * [20] Sussman, M., Smereka, P., and Osher, S., 1994.  “A level set approach for computing solutions to incompressible two-phase flow”.  J. Comp. Phys., 114(1), pp. 146 – 159.
  * [21] Fadlun, E., Verzicco, R., Orlandi, P., and Mohd-Yusof, J., 2000.  “Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations”.  J. Comput. Phys., 161(1), pp. 35–60.
  * [22] Shin, S., and Juric, D., 2009.  “A hybrid interface method for three-dimensional multiphase flows based on front tracking and level set techniques”.  Int. J. Num. Meth. Fluids, 60, pp. 753–778.
  * [23] Shin, S., Chergui, J., and Juric, D., 2020.  “Interaction of a deformable solid with two-phase flows: An eulerian-based numerical model for fluid-structure interaction using the level contour reconstruction method”.  Int J Numer Methods Fluids, 92(11), pp. 1478–1505.
  * [24] Glowinski, R., Pan, T.-W., Hesla, T., and Joseph, D., 1999.  “A distributed lagrange multiplier/fictitious domain method for particulate flows”.  Int. J. Multiph. Flow, 25(5), pp. 755–794.
  * [25] Constante-Amores, C. R., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D., and Matar, O. K., 2020.  “Rico and the jets: Direct numerical simulations of turbulent liquid jets”.  Phys. Rev. Fluids, 5, Nov, p. 110501.
  * [26] Constante-Amores, C., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D., and Matar, O., 2021.  “Direct numerical simulations of transient turbulent jets: vortex-interface interactions”.  J. Fluid Mech., 922, p. A6.
  * [27] Zhong, P., Chuong, C., and Preminger, G., 1993.  “Characterization of fracture toughness of renal calculi using a microindentation technique”.  J Mater Sci Lett, 12, p. 1460–1462.
  * [28] Gharib, M., Rambod, E., and Shariff, K., 1998.  “A universal time scale for vortex ring formation”.  J. Fluid Mech., 360, p. 121–140.
  * [29] Marugán-Cruz, C., Rodríguez-Rodríguez, J., and Martínez-Bazán, C., 2013.  “Formation regimes of vortex rings in negatively buoyant starting jets”.  J. Fluid Mech., 716, p. 470–486.
  * [30] Garimella, S. V., and Rice, R. A., 1995.  “Confined and Submerged Liquid Jet Impingement Heat Transfer”.  J. Heat Transfer, 117(4), 11, pp. 871–877.
  * [31] Sexton, A., Punch, J., Stafford, J., and Jeffers, N., 2018.  “The thermal and hydrodynamic behaviour of confined, normally impinging laminar slot jets”.  Int. J. Heat Mass Transf., 123, pp. 40–53.
  * [32] Sivasamy, A., Selladurai, V., and Rajesh Kanna, P., 2007.  “Numerical simulation of two-dimensional laminar slot-jet impingement flows confined by a parallel wall”.  Int. J. Numer. Methods Fluids, 55(10), pp. 965–983.
  * [33] Lee, H., Yoon, H., and Ha, M., 2008.  “A numerical investigation on the fluid flow and heat transfer in the confined impinging slot jet in the low reynolds number region for different channel heights”.  Int. J. Heat Mass Transf., 51(15), pp. 4055–4068.
  * [34] Hunt, J., Wray, A., and Moin, P., 1988.  “Eddies, streams, and convergence zones in turbulent flows”.  Studying Turbulence Using Numerical Simulation Databases, 1, 11, pp. 193–208.
  * [35] Batchelor, G. K., 1967.  An Introduction to Fluid Dynamics.  Cambridge Mathematical Library. Cambridge University Press.
  * [36] Essmann, E., Shui, P., Popinet, S., Zaleski, S., Valluri, P., and Govindarajan, R., 2020.  “Chaotic orbits of tumbling ellipsoids”.  J. Fluid Mecha., 903, p. A10.
  * [37] Aldoukhi, A. H., Roberts, W. W., Hall, T. L., and Ghani, K. R., 2017.  “Holmium laser lithotripsy in the new stone age: dust or bust?”.  Frontiers in surgery, 4, p. 57.
  * [38] Moore, R. G., Bishoff, J. T., Loenig, S., and Docimo, S. G., 2005.  Minimally Invasive Urological Surgery.  CRC Press.
  * [39] Smith, A. D., 2007.  Smith’s textbook of endourology.  PMPH-USA.
  * [40] Shin, S., Chergui, J., and Juric, D., 2017.  “A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows”.  J. of Mech. Sci. and Tech., 31, pp. 1739–1751.
  * [41] Constante-Amores, C. R., 2021.  “Three-dimensional computational fluid dynamics simulations of complex multiphase flows with surfactants”.  Imperial College London, PhD Thesis.
  * [42] Mordant, N., and Pinton, J. F., 2000.  “Velocity measurement of a settling sphere”.  The European Physical Journal B, 18, 11, pp. 343–352.

## Appendix: Validation of the numerical method

This section provides some validation studies for the immersed solid solver,
flow-solver and mesh studies to provide conclusive evidence of the accuracy of
our numerical predictions.

In order to assess the accuracy of our flow solver, we have validated our
predictions against the experimental and numerical data from [13]. Figure 8
highlights qualitative validation of our numerical frame-work in terms of
numerical streamlines for $\text{Re}=7$ and $\text{Re}=34$. It is clear that
the transient solutions of the Navier Stokes equations provided by our
numerical method are capable of predicting the rich dynamics observed by the
flow-visualisation experiments, and previous numerical simulations based on
the steady solutions of the flow equations. Additionally, we want to highlight
that the flow-solver has also been successfully validated for two-phase
turbulent jets in previous works, such as [25, 26, 41].

For the validation of the immersed solid solver using the Direct Forcing
approach of [21], we have considered the transient and terminal velocities of
a solid sphere settling under gravity in a quiescent fluid. We compare our
numerical predictions against the experiments of [42], where solid spheres are
released in water. We have considered the case of a particle with density
ratio of $\rho_{s}/\rho_{l}=2.56$ and $Re=d_{s}\rho_{s}U_{t}/\mu_{l}=41$,
where$d_{s}$ and $U_{t}$ correspond to the particle diameter and terminal
velocity of the particle (e.g. computations are assumed to reach steady-state
when the settling velocity is below $0.1\%$ for one period of dimensionless
time). The experiments were done under the assumption of an unbounded liquid
bath, therefore our simulations are performed in a sufficiently large domain
to avoid effects from the boundaries. Figure 9a shows the predicted settling
velocity for two different levels of refinement (e.g., $6$ and $10$ cells per
sphere radius, which corresponds to the low (‘LR’) and high resolution (‘HR’)
simulations, respectively). As shown in figure 9a, the predictions from the
immersed solid solver present an excellent agreement with respect to
experimental measures for the transient and terminal velocity (e.g., the
latter with error of less than $0.1\%$). Additionally, we acknowledge that the
same numerical method has been previously validated with respect to the
interaction of a moving solid with the flow, we refer the reader to [23], who
described extensively the numerical formulation and provided extensive
benchmark tests of the fluid-structure interaction of solids (either as rigid
or deformable structures) with multiphase flows using a combination of
immersed boundary and fictitious domain-direct forcing methods.

Finally, we aim to provide conclusive evidence that our numerical results are
mesh-independent. To this end, the kinetic energy for $Re=250$ in the absence
of stones are tested for different mesh resolutions. Figure 9b shows the
temporal evolution of $E_{k}$ for two types of refinements (e.g., ‘LR’ and
‘HR’ refer to levels of refinement characterised by $a/\Delta{\bf x}=6.4$ and
$a/\Delta{\bf x}=12.8$, respectively). As shown in figure 9b, both level of
refinements are capable of predicting accurately the dynamics at steady state.
Thus, we conclude that a mesh characterised with $a/\Delta{\bf x}=6.4$ is
sufficiently refined to ensure mesh-independent results while providing a good
compromise with the computational cost of the simulation. Therefore, we have
proved that the ‘LR’ mesh is capable of predicting the complex dynamics of the
phenomena, and consequently detailed analysis of the vortical structures is
performed using a ‘LR’ mesh-type (unless stated otherwise).

Figure 8:  Direct comparison of the streamlines from our numerical
predictions, third column, with both experimental and simulation results
reported by [13], first and second columns, for a two-dimensional cavity
characterised by $b=1.933$, $h=4.2835$, and cavity-length $\alpha=21.7$ (here,
$\alpha$ stands for the nondimensional length of the rectangular cavity).
Panels (a)-(c) and (d)-(f) represent the flow dynamics with $\text{Re}=7$ and
$\text{Re}=34$, respectively.  Figure 9:  (a) Settling velocity of a solid
sphere at $Re=41$ in relation to the experimental results reported by [42],
here the velocity has been normalised with respect to the experimental
observation at steady state, $U_{t}*$. (b) Mesh study for $Re=250$ in the
absence of stones highlighting the temporal evolution of the kinetic energy
$E_{k}$.


Fourier analysis and (b) convolution matrix generated by re-arranging
(circulant matrix) Fourier coefficients.

The functions for convolution matrix generation are located in
‘convolution_matrix.py’ file for each backend. This part transforms the
structure from the real space to the Fourier space and returns a convolution
matrix (also called Toeplitz matrix) of the Fourier coefficients to apply
convolution operation with the E and H fields. Figure 19 shows the Fourier
coefficients matrix and convolution matrix made from the coefficient matrix.
Code 8 is the definition of ‘conv_solve()’ method and shows how the
convolution matrix generation is integrated inside. As shown in the code,
meent offers 3 different methods to get convolution matrix since each method
has different input type and implementation. This can be chosen by the
argument ‘fft_type’: 0 is for raster modeling with DFS, 1 for raster with CFS
and 2 for vector with CFS.

⬇

1 def conv_solve(self, **kwargs):

2 [setattr(self, k, v) for k, v in kwargs.items()]

3 # needed for optimization

4

5 if self.fft_type == 0: # raster with DFS

6 E_conv_all, o_E_conv_all = to_conv_mat_raster_discrete(self.ucell,
self.fourier_order[0], self.fourier_order[1], device=self.device,
type_complex=self.type_complex, improve_dft=self.improve_dft)

7

8 elif self.fft_type == 1: # raster with CFS

9 E_conv_all, o_E_conv_all = to_conv_mat_raster_continuous(self.ucell,
self.fourier_order[0], self.fourier_order[1], device=self.device,
type_complex=self.type_complex)

10

11 elif self.fft_type == 2: # vector with CFS

12 E_conv_all, o_E_conv_all = to_conv_mat_vector(self.ucell_info_list,
self.fourier_order[0], self.fourier_order[1], type_complex=self.type_complex)

13

14 else:

15 raise ValueError

16

17 de_ri, de_ti, layer_info_list, T1, kx_vector = self._solve(self.wavelength,
E_conv_all, o_E_conv_all)

18

19 self.layer_info_list = layer_info_list

20 self.T1 = T1

21 self.kx_vector = kx_vector

22

23 return de_ri, de_ti

Code 8: ‘conv_solve()’

####  H.3.2 Maxwell’s equations computation

After generating the convolution matrix, meent solves Maxwell’s equations and
returns diffraction efficiencies with the method ‘solve()’. As in the Code 9,
it is a wrapper of ‘_solve()’ method that actually does the calculations and
returns the diffraction efficiencies with other information that is necessary
for the field calculation.

⬇

1 def solve(self, wavelength, e_conv_all, o_e_conv_all):

2 de_ri, de_ti, layer_info_list, T1, kx_vector = self._solve(wavelength,
e_conv_all, o_e_conv_all)

3

4 # internal info. for the field calculation

5 self.layer_info_list = layer_info_list

6 self.T1 = T1

7 self.kx_vector = kx_vector

8

9 return de_ri, de_ti

Code 9: ‘solve()’

Input parameters:

wavelength

: float
The wavelength of the incident light in vacuum.

e_conv_all

: array of {float or complex}
A stack of convolution matrices of the permittivity array; this is
$\llbracket\varepsilon_{r,g}\rrbracket$ in Chapter G. The order of the axes is
the same as of ucell ($Z$ $Y$ $X$).

o_e_conv_all

: array of {float or complex}
A stack of convolution matrices of the one-over-permittivity array; this is
$\llbracket\varepsilon_{r,g}^{-1}\rrbracket$ in Chapter G. The order of the
axes is the same as of ucell ($Z$ $Y$ $X$).

The diffraction efficiencies are 1D array for 1D and 1D-conical grating and 2D
for 2D grating.

####  H.3.3 Field calculation

The ‘calculate_field()’ method in Code 10 calculates the field distribution
inside the structure. Note that the ‘solve()’ method must be preceded. This
function returns 4 dimensional array that the length of the last axis varies
depending on the grating type as shown in Code 11. 1D TE and TM has 3 elements
(TE has Ey, Hx and Hz in order and TM has Hy, Ex and Ez) while the others have
6 elements (Ex, Ey, Ez, Hx, Hy and Hz) as in Figure 20.

⬇

1 field_cell = mee.calculate_field(res_x=100, res_y=100, res_z=100)

Code 10: ‘calculate_field()’

Input parameters:

res_x

: integer
The field resolution in X direction (number of split which the period of x is
divided by).

res_y

: integer
The field resolution in Y direction (number of split which the period of y is
divided by).

res_z

: integer
The field resolution in Z direction (number of split in thickness of each
layer).

field_algo

: integer
The level of vectorization for the field calculation. Default is 2 which is
fully vectorized for fast calculation while 1 is half-vectorized and 0 is
none. Option 0 and 1 are remained for debugging or future development (such as
parallelization).

  * •

0: Non-vectorized

  * •

1: Semi-vectorized: in X and Y direction

  * •

2: Vectorized: in X, Y and Z direction

⬇

1 # 1D TE and TM case

2 field_cell = torch.zeros((res_z * len(layer_info_list), res_y, res_x, 3),
dtype=type_complex)

3

4 # 1D conincal and 2D case

5 field_cell = torch.zeros((res_z * len(layer_info_list), res_y, res_x, 6),
dtype=type_complex)

Code 11: the shape of returned array from ‘calculate_field()’

(a) Ex

(b) Ey

(c) Ez

(d) Hx

(e) Hy

(f) Hz

Figure 20: Field distribution on XY plane from 2D grating structure. (a)-(c):
absolute value of the electric field in each direction, (d)-(f): absolute
value of the magnetic field in each direction.

##  Appendix I Benchmark

In this section, we will address the 1D metasurface problem covered in our
previous work [16] with meent so that we can benchmark and analyze its
capability and functionality.

###  I.1 case application

Figure 21: The image of 1D diffraction metagrating on a silicon dioxide
substrate.

This metagrating deflector is composed of silicon pillars placed on a silica
substrate. The device period is divided into 64 cells, and each cell can be
filled with either air or silicon. The Figure of Merit for this optimization
is set to the deflection efficiency of the +1st order transmitted wave when TM
polarized wave is normally incident from the silica substrate as in Figure 21.

###  I.2 Fourier series implementations

When the sampling frequency of permittivity distribution is not enough,
Fourier coefficients from DFS is aliased. It can be resolved by increasing the
sampling rate that is implemented in the way of duplicating the elements so
the array is extended to have identical distribution but larger array size. We
will call this Enhanced DFS, and it’s implemented in meent as a default
option.

(a) Convergence test

(b) Histogram of the difference compared to Reticolo

Figure 22: Evaluation of 4 different RCWA implementations: Reticolo, meent
DFS, meent CFS and meent Enhanced DFS.  (a) shows diffraction efficiency at
$+1^{st}$ order by FTO sweep of a particular structure. DFS behaves
differently while the other results seem similar. (b) is the histogram of 600k
simulation result (deflection efficiency) difference. Here Reticolo is the
reference and other 3 implementations in meent are benchmarked.

Figure 22(a) illustrates the convergence tests of a particular structure with
four different RCWA implementations. Considering Reticolo as the reference, we
can see CFS is well-matched but DFS shows different behavior. This is due to
the insufficient sampling rate of permittivity distribution, which can be
resolved by Enhanced DFS. Figure 22(b) is the histogram of the discrepancies
from Reticolo result. About 600k structures were evaluated with 4
implementations and the errors of 3 meent implementations were calculated
based on Reticolo. CFS shows the smallest errors and this is because Reticolo
too uses CFS (CFS algorithms in meent are adopted from Reticolo). Enhanced DFS
decreases the error about three orders of magnitudes (e.g., the median of DFS
is 4.3E-4 and this becomes 1.4E-7).

###  I.3 Computing performance

Table 10: Performance test condition backend | device | bit | alpha server | beta server | gamma server
---|---|---|---|---|---
NumPy | CPU | 64 | (A1) | (B1) | (C1)
NumPy | CPU | 32 | (A2) | (B2) | (C2)
JAX | CPU | 64 | (A3) | (B3) | (C3)
JAX | CPU | 32 | (A4) | (B4) | (C4)
JAX | GPU | 64 | - | (B5) | (C5)
JAX | GPU | 32 | - | (B6) | (C6)
PyTorch | CPU | 64 | (A7) | (B7) | (C7)
PyTorch | CPU | 32 | (A8) | (B8) | (C8)
PyTorch | GPU | 64 | - | (B9) | (C9)
PyTorch | GPU | 32 | - | (B10) | (C10)

In this section, computing options to speed up the calculation - backend,
device (CPU and GPU) and architecture (64bit and 32bit) - will be benchmarked.
Table 1 is the hardware specification of the test server and Table 10 is the
index of each test condition.

Figure 23: Performance test: calculation time with respect to FTO. Top row is
the result from 64bit and bottom is from 32bit. The first column is the result
from the test server alpha and the rest is beta and gamma in order.

The graphs in Figure 23 are calculation time vs FTO with all the data per
machine and architecture. Before look into the details, we will briefly
mention some notice in this figure. (1) JAX can’t afford large FTO regardless
of device. We suspect that this is related to JIT compilation which takes much
time and memory for the compilation at the first run. (2) GPU with JAX and
PyTorch can’t accept large FTO even though GPU memory is more than needed for
array upload. (3) if large amount of calculation is needed, Numpy or PyTorch
on CPU is the option. (4) no golden option exists: it is recommended to find
the best option for the test environment by doing benchmark tests.

We will visit these computing options one by one. The option C9 at FTO 1600
will be excluded in further analyses: this seems an optimization issue in
PyTorch or CUDA.

####  I.3.1 Backend: NumPy, JAX and PyTorch

NumPy, JAX and PyTorch as a backend are benchmarked. NumPy is installed via
PyPI which is compiled with OpenBLAS. There are many types of BLAS libraries
and the most representative ones are OpenBLAS and MKL (Math Kernel Library).
As of now, PyPI provides NumPy with OpenBLAS while conda does one with MKL.
This makes small discrepancy in terms of speed and precision hence pay
attention when doing consistency test between machines.

Figure 24: Performance test: calculation time by FTO sweep. The result is
normalized by NumPy case from the same options to compare the behavior of
other backends. In these plots, black dashed line is $y=1$ and the results of
NumPy cases lie on this line since they are normalized by themselves.

Figure 24 is the relative simulation time per server and architecture
normalized by the time of NumPy case in the same conditions to make comparison
easy. In small FTO regime, all the options were successfully operated and no
champion exists. Hence it is strongly recommended to run benchmark test on
your hardware and pick the most efficient one. In case of X7 (A7, B7 and C7),
Alpha and Gamma show the same behavior - spike in 100 - while beta shows
fluctuation around B1. One possible reason for this is the type of CPU. The
CPUs of Alpha and Gamma belong to ‘Xeon Scalable Processors’ group but Beta is
‘Xeon E Processors’. Currently we don’t know if this actually makes difference
or some other reason (such as the number of threads or BLAS implementation)
does. This result may vary if MKL were used instead of OpenBLAS. In large FTO,
only two options are available: NumPy and PyTorch on CPU in 64 bit. In case of
JAX, the tests were failed: we watched memory occupation surge during the
simulation which seems unrelated to matrix calculation. This might be an issue
of JIT (Just In Time) compilation in JAX. Between NumPy and PyTorch, PyTorch
is about twice faster than NumPy in both architectures at Alpha and Gamma, but
beta shows different behavior. This too, we don’t know the root cause but one
notable difference is the family of CPU type.

####  I.3.2 Device: CPU and GPU

Figure 25: Performance test result. The calculation time of GPU cases are
normalized by CPU cases from the same options to see the efficiency of GPU
utilization. In these plots, black dashed line is $y=1$ where the capability
of both are the same.

Figure 25 shows the relative simulation time of GPU cases normalized by CPU
cases on the same backend and architecture. Note that it is relative time, so
the smaller time does not mean it is a good option for the simulation
experiments: the relative time can be small even if the absolute time of CPU
and GPU are very large compared to other options.

JAX shows good GPU utilization throughout the whole range (except one point in
beta) regardless of the architecture. Considering the architecture, the data
trend in beta is not clear while the gamma clearly shows that GPU utilization
can be more effective in 32bit operation. PyTorch data is a bit noisier than
of JAX, but has the similar behavior per server. The data in beta is hard to
conclude as the JAX cases and the gamma too shows ambiguous trend but we can
consider GPU option is efficient with wide range of FTOs.

Up to date, eigendecomposition for non-hermitian matrix which is the most
expensive step ($O(M^{3}N^{3})$) in RCWA, is not implemented on GPU in JAX and
PyTorch hence the calculations are done on CPU and the results are sent back
to GPU. As a result, we cannot expect great performance enhancement in using
GPUs.

####  I.3.3 Architecture: 64 and 32 bit

Figure 26: Performance test result. The calculation time of 32bit cases are
normalized by 64bit cases from the same options. In these plots, black dashed
line is $y=1$ where the capability of both are the same.

In Figure 26, calculation time of 32bit case is normalized by 64bit case with
the same condition. With some exceptions, most points show that simulation in
32bit is faster than 64bit. Here are some important notes: (1) From our
understanding, the eigendecomposition (Eig) in NumPy operates in 64bit
regardless of the input type - even though the input is 32bit data (float32 or
complex64), the matrix operations inside Eig are done in 64bit but returns the
results in 32bit data type. This is different from JAX and PyTorch - they
provides Eig in 32bit as well as 64bit. Hence the 32bit NumPy cases in the
figure approach to 1 as FTO increases because the calculation time for Eig is
the same and it is the most time-consuming step. (2) Keep in mind that 32bit
data type can handle only 8 digits. This means that 1000 + 0.00001 becomes
1000 without any warnings or error raises. For such a reason, the accuracy of
32bit cases in the figures are not guaranteed - we only consider the
calculation time. (3) Eig in PyTorch shows interesting behavior: as FTO
increases, calculation time in 32bit overtakes 64bit - see A8/A7, B8/B7 and
C8/C7. This is counter-intuitive and we don’t have good explanation but
cautiously guess that this might be related to the accuracy and precision in
Eig or an optimization issue of PyTorch.

##  Appendix J Applications

meent is expected to be useful for solving inverse design and optimization
problems in OCD and metasurface design. In this section, we present some
exemplary cases where meent proves its capabilities. Leveraging the automatic
differentiation function, we successfully carry out optimization for
diffraction gratings and achieve inverse design of the geometric parameters.

###  J.1 Inverse design of 1D diffraction grating

(a) Optimization curve

(b) Electric field

Figure 27: Optimization result of 1D beam deflector. (a) The deflection
efficiencies are calculated for every iterations and this experiment is
repeated 100 times with random starting points (b) Electric field distribution
from the final structure

In this example, we will optimize 1D beam deflector that was used for
benchmark in Chapter I using AD with these options - 256 cells, FTO $=100$,
$\lambda_{0}=900$ nm and the deflection angle $=50^{\circ}$. During the
optimization, each cell can have non-binary refractive index values, leading
to a gray-scale optimization. To obtain the final structure consisting of only
silicon/air binary structures, an additional binary-push process is required.
The initial structure for optimization is randomly generated so the each cell
can have the refractive index value between of air and silicon under uniform
distribution. The Figure of Merit for this optimization process is set to the
$+1^{st}$ order diffraction efficiency, and the gradient is calculated by AD.
The refractive indices are updated over multiple epochs using the ADAM
optimizer [89] with the learning rate of 0.5.

Figure 27(a) shows the deflection efficiency change by iteration. Two solid
lines are averaged value of all the samples at the same iteration step. Shaded
area is marked with $\pm$ standard deviation from the average. The blue line
(Before binarization) is the result of device with any real number between two
refractive indices (silicon and air), which is non-practical, and the orange
line (After binarization) is the final device composed of silicon and air. The
best result we found is 89.4%.

###  J.2 Inverse design of 2D diffraction grating

Figure 28: Optimization result of 2D beam deflector. (a) The schematic of 2D
beam deflector and the final structure after optimization. (b) Convergence
test of the initial structure. (c) Learning curve of structure optimization
for 110 epochs. Spatial blurring and binary push is applied on each epoch (d)
The electric field distribution of the optimized structure in XZ plane.

Here, we demonstrate optimization of a 2D diffraction metagrating as shown in
Figure 28a. Similar to the previous 1D diffraction metagrating, the 2D
diffraction metagrating also consists of silicon pillars located on top of a
silicon dioxide substrate. TM polarized wave with $\lambda=1000$ nm is
normally incident from the bottom of the substrate and the device is designed
to deflect the incident light with deflection angle $\theta=60^{\circ}$ in
$X$-direction. The device has a rectangular unit cell of period
$\lambda/\sin\theta\approx 1150$ nm and $\lambda/2=500nm$ for the x and
y-axis, respectively. Moreover, the unit cell is gridded into $256\times 128$
cells which is either filled by air or silicon. The convergence of RCWA
simulation for different number of Fourier harmonics are plotted in Figure
28b. Considering the trade-off between simulation accuracy and time, we set
$N_{x}=13$ and $N_{y}=10$.

After 110 epochs of optimization, the final structure achieves an efficiency
of 92% and successfully deflects the incoming beam at a 60∘ angle (Figure
28d). The optimized structure and the learning curve are presented in Figure
28a and Figure 28c, respectively.

###  J.3 Inverse design of 1D grating color router

Figure 29: A schematic of a color router. The incoming light is guided to
subpixels of corresponding wavelength.

Until now, we have focused on the problems where the FoM was simply defined.
However, in this example, we aim to demonstrate the optimization process of a
meta color router, which involves a complex FoM.

A meta color router is an optical component designed for next-generation image
sensors. It is designed to route the incoming light to the subpixel region of
corresponding color, as depicted in Figure 29. In this exemplary case, we
consider an RGB meta color router featuring a pixel pitch of $0.5\mu m$ that
consists of vertically stacked 1D binary gratings. The constituent dielectrics
are silicon dioxide and silicon nitride, with fixed refractive indices of 1.5
and 2.0, respectively. The meta device region (width of $1.5\mu m$ and height
of $2\mu m$) is sliced into 8 layers with 64 cells per layer.

The FoM for this meta color router is defined as the average of TE mode
electric field intensity over the corresponding subpixel region, as given by
Equation (74).

$\displaystyle
FoM=\frac{1}{N}\sum_{\lambda_{1}}^{\lambda_{N}}\frac{\int_{x_{1}}^{x_{2}}|\vec{\mathbf{E}}(\lambda)|^{2}dx}{\int_{0}^{P}|\vec{\mathbf{E}}(\lambda)|^{2}dx}\times
T(\lambda)$ (74)

Here, $\mathbf{E}$ represents the electric field within the subpixel region,
while $T$ represents transmittance. The parameter $x\in(x1,x2)$ determines the
desired subpixel region corresponding to the incident beam wavelength. For
simplicity, we define the wavelength ranges for R, G, and B as 600 nm - 700
nm, 500 nm - 600 nm, and 400 nm - 500 nm, respectively. Throughout the
optimization process, optical efficiencies are averaged across 9 wavelength
points to ensure a finely tuned broadband response.

The optimization procedure for the meta color router follows a similar
approach to the previous examples, including random initialization,
optimization via back-propagated gradients with or without binary push. The
optimization curve and the final binarized device structure are shown in
Figure 30.

Figure 30: Optimization result of 1D grating meta color router. (a)
Optimization curve of greyscale device and binary-pushed device at each epoch.
(b) Color sorting efficiency spectrum. (c) The electric field inside the final
color router device.

##  Appendix K Licenses

Development

  * •

Numpy: BSD license

  * •

JAX: Apache License 2.0

  * •

PyTorch: BSD license (BSD-3)

Experiment

  * •

Ray: Apache License 2.0

  * •

SheepRL: Apache License 2.0

  * •

neuraloperator: MIT license