\section{Introduction}
GEANT4 is a C++ based open-source toolkit for simulating the passage of particles through matter\cite{doi:10.1016/S0168-9002(03)01368-8}. Often used for high energy physics and accelerator design, GEANT4 also has a range of scientific/engineering applications such as design of radiation shielding, medical equipment, and space hardware\cite{doi:10.1016/j.nuclphysbps.2004.10.083,doi:10.1016/j.nima.2016.06.125}. The transport of particles through aerosols is important both in traditional application such as particle detectors\cite{doi:10.1038/srep41699, doi:10.1016/j.jenvrad.2016.03.025} and in many alternate applications such as environmental studies\cite{isbn:0521807670, doi:10.1007/978-0-387-87664-1_10, doi:10.5772/50248} and astronomical studies\cite{doi:10.3847/0004-637x/818/2/164, doi:10.1088/0004-637x/737/2/105, doi:10.1086/590482, doi:10.3847/1538-4357/ab1b20}, but this transport is complicated\cite{doi:10.1016/j.jaerosci.2018.12.010} since, for some aerosols, their granularity is important\cite{doi:10.3847/0004-637x/818/2/164}.

The simplest aerosol model, a single-volume mixture of the background and droplet materials in the correct mass fractions, is sufficient when this granularity is unimportant, but we experimentally find this only to be true when droplets are small. The failure of this approximation can be understood because a single volume model implicitly posits equal interaction probability throughout the aerosol (uniform distribution of droplet material) but, when aerosols consist of large droplets, there is a heightened probability for interactions in quick succession.

The granularity is most directly respected by modelling the aerosol as a collection of (non-intersecting) droplets, which can be done with GEANT4's `parameterised solid' geometry. A parameterised solid consists of one or more sub-objects that are repeatedly placed in the world according a parameterisation function that allows differences in position, size, rotation, material, and more between copies\cite{bookForApplicationDevelopers}. While such solids are maximally flexible, they experience computational difficulties when the sub-object is populated millions or billions of times, as in the case of many aerosols. To allow for such aerosols with many large droplets to be simulated in GEANT4, we present the `fastAerosol' geometry class which models the aerosol as a collection of droplets but with three following differences from parameterised solids.

First, while both fastAerosol and parameterised solids organize volumes in voxels, fastAerosol utilizes a simpler voxelization algorithm which only considers the average droplet number density and droplet bounding radius in contrast to parameterised solids' algorithm which utilizes the location of the sub-objects in voxelization.

Second, fastAerosol differs from parameterised solids by not needing to load/generate the droplets at the start of the simulation: fastAerosol allows dynamic population of droplets as particles travel through an aerosol (`lazy' population). This allows speedups especially in aerosols where the effective volume (in which particles transport through) is significantly smaller than the total bulk volume or when the number of primaries is low.

Third, the fastAerosol class utilizes the droplet geometry to speed up calculations. The motivation is that, if the distance to the center of a droplet gives knowledge of the distance to the surface $\pm\delta$, then the closest droplet must be in the spherical shell with radial width $R_0\leq R\leq R+2\delta$ for $R_0$ the distance to the closest droplet center.

\section{\texorpdfstring{\MakeLowercase{fast}A\MakeLowercase{erosol}}{fastAerosol} Geometry}
\iffalse
We present a geometry class alternative to parameterised solids for modelling aerosols with roughly spherical droplets. With the assumption on droplet shape, calculations such as finding the distance to the aerosol can be performed by first finding the closest droplet and only then performing the distance calculation with this closest droplet, contrasted with parameterised volumes which must calculate the distance of the query point to every droplet that its checking.
\fi

This report presents an aerosol model, `fastAerosol', implemented in GEANT4 and representing an aerosol as a collection of uniformly shaped droplets. Volumes in GEANT4 are described via a hierarchy\cite{bookForApplicationDevelopers}: at the lowest level, a volume is described by its shape in a `solid volume'. Next, the shape is given a material and physical properties in a `logical volume'. Finally, the shape is given a position and rotation in a `physical volume'. Since the presented fastAerosol class is strictly geometrical, this class is on the solid volume level of the hierarchy.

The model consists of a helper class, `fastAerosol', and a solid volume class, `fastAerosolSolid'. The fastAerosol helper class contains the droplet positions, methods to find nearest droplets, and methods to populate droplets; the fastAerosolSolid class contains the GEANT4 requisite solid functions (e.g., DistanceToIn). We refer to the collection of the two classes as `fastAerosol'.

\subsection{Droplet Distribution}
\label{sec:distributionMethods}
The fastAerosol class was designed to be usable with only knowledge of the aerosol's macroscopic properties. While the aerosol's microstate requires description of each droplet's position and rotation, this information may not be available/relevant to the user. To allow the user to be unconcerned of microstate properties, fastAerosol automatically generates droplets in the aerosol based off of an average droplet number density, a position distribution, and a rotation distribution.

By allowing fastAerosol to control the microstate, computational gains can be achieved by lazily populating droplets. To achieve such population, the aerosol is voxelized and droplets are only generated/loaded in voxels which could contain the closest droplet for a given calculation. The voxelization pitch, $p_{\text{grid}}$, is fixed such that there are an average of $N_{\text{droplets/voxel}}$ droplets per voxel (default$=4$; chosen from experimentation to minimize simulation time):
\begin{equation}
    p_{\text{grid}} = \left( \frac{ N_{\text{droplets/voxel}} }{ \langle n_d \rangle } \right)^{1/3}
\end{equation}
for an aerosol with average droplet number density $\langle n_d \rangle$.

These voxels are contained in a vector of length $N_{\text{tot}} = N_x\, N_y\, N_z$ for $N_x$, $N_y$, and $N_z$ the number of voxels in the $x$, $y$, and $z$-directions respectively (appropriately aligned with the bounding box). A voxel with center $(x,y,z)$ is assigned an index $I(x,y,z)$ in this vector given by
\begin{equation}
    \label{eq:index}
    I(x,y,z) = \frac{(x-x_0)+ N_x \left[(y-y_0) + N_y (z-z_0)\right]}{p_{grid}}
\end{equation}
for $(x_0, y_0, z_0)$ the coordinates of the center of the `minimal' voxel. That is, $(x_0,y_0,z_0)$ are defined so that $x_0\leq x$ and $y_0\leq y$ and $z_0\leq z$ for $(x,y,z)$ the coordinates of the center of any voxel in the bounding box's voxelization. We will call $(x-x_0)/p_{grid}$ the voxel $x$-coordinate and similar for $y$ and $z$.

Two values must be determined for each voxel: a random seed, $s(x,y,z)$, and the expected number of droplets per voxel, $\mu(x,y,z)$. First, $s(x,y,z)$ is chosen to be:
\begin{equation}
    s(x,y,z) = I(x,y,z) + s_{global} N_{\text{tot}},
\end{equation}
for $s_{global}$ a random seed for the entire program. Multiplying $s_{global}$ by $N_{\text{tot}}$ guarantees completely new voxel seeds for different integer global seeds.

Second, $\mu(x,y,z)$ is set based off of the spatial average, $\langle n_d \rangle$, and a distribution, $f$. This is done by first calculating a value, $\mu_0$, proportional to the mean:
\begin{equation}
    \mu_0(x,y,z) = \langle n_d \rangle\, \max\left[0, V(x,y,z)\,f(x,y,z)\right]
\end{equation}
where $V(x,y,z)\in [0,1]$ is the fraction of the volume in the voxel that allows placement of a droplet center (inside the aerosol bulk and at least $r$ away from the edge for droplet radius $r$). If $(x,y,z)$ is at least $r+\sqrt{3}p_\text{grid}/2$ from the aerosol surface, then the voxel is fully contained in the aerosol so $V(x,y,z)=1$. Otherwise, this overlap fraction is manually calculated by placing $100$ points with uniform random distribution in the voxel and setting $V(x,y,z)=N_\text{good}/100$ where $N_\text{good}$ is the number of placed points which allow droplet placement.

Since no restrictions are placed on $f$, this $\mu_0(x,y,z)$ is not guaranteed to give the correct $\langle n_d \rangle$, so $\mu(x,y,z)$ must be scaled from $\mu_0(x,y,z)$:
\begin{equation}
    \mu(x,y,z) = \mu_0(x,y,z) \frac{\langle n_d \rangle\, V_B}{\sum_{\text{voxels}} \mu_0(x',y',z')}
\end{equation}
for $V_B$ the volume of the aerosol bulk. This guarantees that $\sum_{\text{voxels}} \mu(x,y,z) = \langle n_d \rangle\, V_B$, as desired.

When performing distance calculations, voxels are lazily populated: if an unpopulated voxel at $(x,y,z)$ could have the closest droplet, a Poisson random number of droplets (mean $\mu(x,y,z)$ and seed $s(x,y,z)$) are placed with centers uniform randomly distributed in the voxel. If a center is to be placed less than $r$ from aerosol boundary or less than $2r+minSpacing$ from any other droplet center (for some user-specified buffer, $minSpacing$), the position is regenerated. If a single center fails placement a user-specified amount of times (default is $100$ attempts), the simulation skips that placement. In the event that the number of skipped centers exceeds a user-specified maximum (default is $1\%$ of the total droplets expected in the aerosol), the simulation ends with an error.

\subsection{Closest Droplets}
GEANT4 distance calculations are sensitive to the droplet with the closest surface, so discussion is focused on determination of said droplet. The fastAerosol class finds the droplet with the closest surface by searching droplet centers for a collection of droplets guaranteed to contain it and then iterating over this collection to find the one with the closest surface. There are two definitions of distance relevant to GEANT4 calculations: absolute distance and distance along a vector, both of which are covered in the following sections.

To generate this collection of candidates for either definition of distance, the relationship between the distance to a droplet's center and the distance to its surface is required. Define $r_\text{in}$ to be the radius of the largest sphere, centered at a droplet, that can be drawn fully inside the droplet; define $r_{\text{droplet}}$ to be the radius of the smallest sphere, centered at the same droplet, that contains the droplet (also called the `bounding radius'). Since droplets are assumed to be uniform in shape/size, these $r_\text{in}$ and $r_\text{droplet}$ hold for all droplets.

For $R_0$ the distance from a query point to the center of the droplet with the closest center, it is then known that the droplet with the closest surface is centered in the shell (centered at the query point) with radius range $R_0<R<R_0+\sigma$ where $\sigma=r_{\text{droplet}}-r_\text{in}$. This $\sigma$ bound is only strictly correct for finding the absolutely closest droplet, but we also use it as an approximation for the closest droplet along a vector.

\subsubsection{Absolute Closest Droplet}
\label{sec:absoluteClosest}
To calculate the distance from a query point $\vec{p}$ to the fastAerosol aerosol, define $R_\text{bulk}$ to be the distance of $\vec{p}$ from the aerosol bulk's surface. Additionally, let $\mathscr{C}$ be the collection of candidate droplets and $R_0$ be the minimum distance from $\vec{p}$ to the center of any droplet in $\mathscr{C}$ (currently $\mathscr{C} = \{\}$ so we treat $R_0=\infty$).

To fill $\mathscr{C}$, define a `voxel radius' $R_{voxel}=\left \lfloor R_\text{bulk}/p_{grid}\right \rfloor$ and search spherical shells of voxels centered at $\vec{p}$ with radius $R\geq R_{voxel}$ (shells generated by the modified mid-point circle algorithm defined in \cite{doi:10.1016/j.jcp.2013.01.035}). Upon finding any droplet, it is added to $\mathscr{C}$ and $R_0$ is updated. The collection is complete when $R>\lceil 0.25+(R_0+\sigma)/p_{grid} \rceil$ or when $R$ exceeds a user defined maximum. The upper limit on $R$ is chosen so that, at this limit, there are no points in the voxelized shell at radius $R+1$ that are closer than $R_0+\sigma$ (upper limit calculated in $2D$ and used as approximation for $3D$).

The absolute closest droplet is then found by iterating over all droplets in $\mathscr{C}$ that have centers closer than $R_0+\sigma$ and returning the droplet with the closest surface distance (return no droplet if $\mathscr{C}=\{\}$).

\subsubsection{Vector Closest Droplet}
To calculate the closest droplet to a particle at point $\vec{p}$ along normalized vector $\vec{v}$, first, define a voxel search width, $w$, as
\begin{equation}
    w = 2 \left \lceil{\frac{r_{\text{droplet}}}{p_{\text{grid}}}}\right \rceil+1,
\end{equation}
and a temporary position, $\vec{p}\,'$, beginning at $\vec{p}\,' = \vec{p}$ with distance $R_\text{bulk}$ from the aerosol bulk (let $R_\text{bulk}$ automatically update when $\vec{p}\,'$ changes). Additionally, let $\vec{d}$ be the closest droplet found so far and $R_0$ the distance to this droplet's center (treat $R_0=\infty$ at start). The search algorithm, then, is to
\begin{enumerate}
    \item if $R_\text{bulk}=\infty$, return $\vec{d}$; otherwise, if $R_\text{bulk}>p_{\text{grid}}$, update $\vec{p}\,' \to \vec{p}\,' + (R_\text{bulk}-p_\text{grid}) \vec{v}$.
    \item If $|\vec{p}\,'-\vec{p}|$ is greater than a user defined maximum, return $\vec{d}$; otherwise,
    \item search a $w\times w\times w$ grid of voxels centered at the voxel containing $\vec{p}\,'$, calculate the vector distance to the surface of any droplets in this grid and update $\vec{d}$ and $R_0$ appropriately.
    \item If $|\vec{p}\,'-\vec{p}|>R_0+\sigma$, return $\vec{d}$; otherwise,
    \item Find the smallest $n\in\mathbb{N}$ such that $\vec{p}\,' + n\,0.99\,p_\text{\text{grid}}\vec{v}$ lies in a new voxel, define $\vec{\delta}=\delta_x\hat{x} + \delta_y\hat{y} + \delta_z\hat{z}$ to be the difference in voxel positions between the $\vec{p}\,' + n\,0.99\,p_\text{\text{grid}}\vec{v}$ and $\vec{p}\,'$, and update $\vec{p}\,' \to \vec{p}\,' + n\,0.99\,p_\text{\text{grid}}\vec{v}$.
    \item If $R_\text{bulk}>p_{\text{grid}}$, return step $1$; otherwise,
    \item for each non-zero $\delta_\alpha$ for $\alpha\in\{x,y,z\}$, search for droplets the `end caps' $\vec{p}_\text{old}\,'+ (w+1)\delta_\alpha \hat{\alpha} + a\hat{\beta} + b\hat{\gamma}$ for $\hat{\beta},\hat{\gamma}\in\{\hat{x},\hat{y},\hat{z}\}\setminus\hat{\alpha}$ with $\hat{\beta}\neq\hat{\gamma}$ and with $a,b\in[-w,w]\cap\mathbb{Z}$.
    \item If $|\vec{\delta}|>1$, also search for droplets the edges/corners between the end caps.
    \item Update $\vec{d}$ and $R_0$ appropriately. Return to step $4$.
\end{enumerate}

Since the surface distance needs to be calculated to know if a droplet is a candidate (the primary could travel through the droplet's bounding radius but not actually intersect the droplet), the explicit generation of the collection $\mathscr{C}$ and then iteration through $\mathscr{C}$ is unnecessary.

\subsection{Shape Flexibility}
The fastAerosol class accepts any bulk shape in the form of a solid volume object, which is valuable since aerosols often have non-rectangular shapes such as tori\cite{doi:10.3847/0004-637x/818/2/164} and the shapes of atmospheric clouds. Similarly, fastAerosolSolid accepts any droplet shape in form of a solid volume object, with the disclaimer that, while distance calculations to droplets of any shape are accurate, these droplets will be placed as if they were spheres of their bounding radius. Thus fine packing of elongated objects (e.g., sticks) is currently not possible with this package.

By only containing the centers in the fastAerosol helper class, construction of structured droplets is simple: associate multiple fastAerosolSolid objects per fastAerosol helper class. Droplet rotation distributions are allowed.

\subsection{Code Implementation}
The fastAerosol class achieves simpler user implementation by its reliance only on macroscopic properties. While implementation of a parameterised aerosol requires the user to provide droplet positions and rotations at construction, fastAerosol only needs distributions for position and rotation:
\begin{lstlisting}
fastAerosol* aerosol =
    new fastAerosol(``aerosol", bulkShape, boundingR,
                    minSpacing, avgDropNumDens, sigma,
                    positionDistribution);

fastAerosolSolid* solidAerosol =
    new fastAerosolSolid(``aerosolSV", aerosol,
                         dropletShape,
                         rotationDistribution));
\end{lstlisting}
where the first argument to each function is its name and the variables are self-explanatory. The distribution arguments (`positionDistribution' and `rotationDistribution') and `sigma' are optional, taking default values of uniform position/rotation and $\sigma=0$. More detailed implementation may be found in the full source code, available as an advanced example in GEANT4 (as of version 10.7).

\section{Experimental Methods}
To test fastAerosol's accuracy and computational efficiency, a fastAerosol aerosol is generated at the origin with a $150\times150\times5000$ mm bounding box (see \cref{fig:experiment} displaying an ellipsoidal aerosol) with varying droplet size, droplet number density, and bulk/droplet shape (droplet size/shape uniform for all droplets in aerosol). Transport is verified by shooting a variable number of $50$ MeV protons along the $z$-axis from a Gaussian distribution centered at $(0,0,-2612.5)$ mm with standard deviation $\sigma=16.5$ mm in the $x$ and $y$ directions. The global seed is allowed to vary across $s_{global}\in[1,5]\cap\mathbb{N}$ for statistics.

\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{experiment_close.png}
	\caption{Ellipsoidal aerosol geometry with $50$ MeV protons shot through it from $(0,0,-2612.5)$ mm (with $x$ and $y$ spread $\sigma=16.5$ mm each) along the $+z$-direction.}
	\label{fig:experiment}
\end{figure*}

The following experiments were performed on GEANT4 10.04 (patch 02) running on an Amazon Web Services (AWS) `c5.24xlarge' computer. We have provided the source code as an advanced example in GEANT4 (as of version 10.7). The data reflects code compiled in release mode and ran in multi-threaded mode.

\subsection{Droplet Distribution}
\label{sec:distribution}
For the droplet distribution tests, droplets with bounding radius $r=1$ mm and number densities of $n=10^{-6+k/5}$ mm$^{-3}$ for $k\in\mathbb{Z}\cap [0,20]$ were populated into a box-shaped aerosol. Droplet placement in the aerosol was performed by individually populating voxels in order of increasing $o(x_i,y_i,z_i)$, defined as
\begin{equation}
    \label{eq:order}
    o(x_i,y_i,z_i) = z_i + N_z(y_i + N_y x_i)
\end{equation}
for voxel position $(x_i,y_i,z_i)$. This is effectively in order of voxel index (\cref{eq:index}) except with $x$ and $z$ swapped. The distribution of droplets is then saved and analyzed to verify that droplets are roughly uniform in the clouds and thus actually modeling the aerosol of interest.

\subsection{Particle Transport}
Five particle transport experiments were performed to test the following:
\begin{enumerate}
    \item the dependence of particle transport on droplet size,
    \item the dependence of relative computational efficiency between fastAerosol and parameterised aerosols on droplet number density,
    \item the dependence of this relative computational efficiency on primary count,
    \item the dependence on parameterised computational efficiency on its voxelization fineness, and
    \item the dependence of particle transport on bulk/droplet shape.
\end{enumerate}
Four aerosol models are used for these tests: a parameterised model, a fastAerosol model, a fastAerosol model which has droplet positions generated before shooting any primaries (`pre-populated' fastAerosol; as in \cref{sec:distribution}), and a single volume model (`smooth'; density set as average aerosol density). The parameterised aerosol loads droplet positions generated by the pre-populated fastAerosol.

For droplet size investigations, $5000$ protons are shot through parameterised, fastAerosol, pre-populated fastAerosol, and smooth box-shaped aerosols consisting of $n=0.01$ mm$^{-3}$ spherical droplets with radius $0.1\leq r\leq1.50$ mm. The energy deposited into the detector aluminum block is recorded by a $20\times 20\times 1$ scoring grid.

For droplet number density efficiency tests, $5000$ protons are shot through parameterised, fastAerosol, and pre-populated fastAerosol box-shaped aerosols consisting of spherical droplets with radius $r=1$ mm and number density $n=10^{-6+k/5}$ mm$^{-3}$ for $k\in\mathbb{Z}\cap [0,20]$. The total simulation time along with the program's maximum resident set size was measured via `/usr/bin/time'.

For primary count efficiency tests, $1$ to $10000$ protons are shot through parameterised, fastAerosol, and pre-populated fastAerosol box-shaped aerosols consisting of $n=10^{-3}$ mm$^{-3}$ spherical droplets with radius $r=1$ mm. The simulation time draw was measured by `/usr/bin/time'.

For parameterised solid voxelization tests, $5000$ protons are shot through a box-shaped aerosol consisting of $n=10^{-3}$ mm$^{-3}$ spherical droplets with radius $r=1$ mm. The `smartless' parameter, defined as the average number of voxels per geometry object in the parameterised solid, measures the fineness of parameterised solids' voxelization (higher value implies finer voxelization)\cite{bookForApplicationDevelopers}. GEANT4 automatically calculates an `optimal' smartless value and then sets the real smartless value as the minimum of the calculated and user-set values. To test how the voxelization fineness affects simulation time and memory, the user-set smartless parameter is allowed to vary from $2^{-14}$ to $2^5$ and  both simulation time and memory draw were measured by `/usr/bin/time'.

Finally, for droplet shape tests, $5000$ protons are shot through a box-shaped fastAerosol aerosol with $n=0.01$ mm$^{-3}$ droplets shaped like cylinders, boxes, hemispheres, and spheres (all maximally sized to fit in a bounding radius of $r=1$ mm; no rotations applied - aligned along z-axis) and the energy deposited into the aforementioned detector/scoring grid recorded. Similarly, for bulk shape investigations, $5000$ protons are shot through a fastAerosol aerosol containing $n=0.01$ mm$^{-3}$ spherical droplets with radius $r=1$ mm inside a bulk shaped like
\begin{enumerate}
    \item a box ($150\times150\times 5000$ mm),
    \item a cylinder ($r=75$ mm, $\text{length}=5000$ mm; aligned along $z$-axis),
    \item an ellipsoid ($r_\text{minor}=75$ mm, $r_\text{major}=2500$ mm; major radius aligned along $z$-axis), and
    \item a pipe ($r_\text{outer}=75$ mm, $r_\text{inner}=37.5$ mm, $\text{length}=5000$ mm; aligned along $z$-axis).
\end{enumerate}
Again, the energy deposited into the detector/scoring grid was recorded.

\section{Results}
\subsection{Droplet Distribution}
The uniformity of the droplet $x$-position distribution is verified using a one-sample Kolmogorov-Smirnov test with the null hypothesis that the droplet $x$-positions come from a uniform distribution. The distributions from the five seeds were concatenated for better statistics. Under a statistical cutoff of $\alpha=0.05$, the $n=10^{-6}$ mm$^{-3}$ data fails to reject the null hypothesis while the $n=10^{-2}$ mm$^{-3}$ data succeeds in rejecting this hypothesis ($p\approx 44\%$ and $p\approx 2\%$ of drawing a distribution from the null hypothesis more extreme than measured for $n=10^{-6}$ mm$^{-3}$ and $n=10^{-2}$ mm$^{-3}$ respectively). This suggests that, while the $n=10^{-6}$ mm$^{-3}$ data is roughly uniform, the $n=10^{-2}$ mm$^{-3}$ data is not.

This failure of uniformity can be visualized in \cref{fig:distribution_plot}, displaying in a solid black line the observed cumulative distribution function (CDF) of droplet $x$-position (other directions similar) minus the expected uniform distribution. The 95\% lower confidence bound (LCB) and upper confidence bound (UCB) are plotted in dotted blue lines, similarly minus the expected uniform distribution. The CDF for $n=10^{-6}$ mm$^{-3}$ data is displayed on top; the CDF for $n=0.01$ mm$^{-3}$ data is displayed on bottom.

\begin{figure}
	\centering
	\includegraphics[width=0.75\linewidth]{distribution_plot}
	\caption{Cumulative distribution plot of droplet $x$-position minus the expected CDF of a uniform distribution. Top: the $n=10^{-6}$ mm$^{-3}$ $x$-position distribution is generally uniform, with the UCB of the observed distribution consistently higher than the uniform expectation and the LCB lower than the uniform expectation. Bottom: the $n=10^{-2}$ mm$^{-3}$ $x$-position distribution is observably non-uniform, with the UCB and LCB both generally smaller than the uniform expectation and a periodic bunching at frequency $\sim7.4$ mm (matching voxel pitch), suggesting a $\sim0.01\%$ heightened probability of droplet placement at voxel boundaries.}
	\label{fig:distribution_plot}
\end{figure}

For truly uniform data, one would expect the predicted uniform CDF would lay between the observed LCB and UCB, implying that, on \cref{fig:distribution_plot}, the UCB-Uniform curve should be mostly positive while the LCB-Uniform curve mostly negative. This is observed in the $n=10^{-6}$ mm$^{-3}$ data, as expected since our statistical test failed to distinguish this data from a uniform distribution. The $n=10^{-2}$ mm$^{-3}$ data, however, does not show this feature. Namely, the empirical data, LCB, and UCB are all consistently lower than a uniform distribution, except for periodic spikes in the CDF. This indicates that the observed droplet number density is, on average, low, with periodic bunching. The bunching period of $\sim7.4$ mm matches the grid pitch of $7.368$ mm and, suggesting heightened droplet population at voxel boundaries.

Despite this deviation, the overall distribution is still qualitatively uniform: these spikes in probability correspond to an increased likelihood of finding a droplet with $x$-position near a voxel boundary of $\sim 0.01\%$. This level of uniformity is expected to be sufficient for all but the most precise simulations. Regardless, this bunching should be resolved in future work.

\subsection{Droplet Size}
The energy deposited into the detector after a parameterised, fastAerosol, pre-populated fast-Aerosol, and smooth aerosol is shown in \cref{fig:net_plot} (error bars too small for visualization). The energy deposited after the smooth aerosol drops off approaching $r=1$ mm, a feature absent from all of the other simulations. This most likely indicates that, for `large' droplet radii, the granularity is important for the particle transport, as expected since aerosols with large droplets have a heightened probability for multiple interactions in quick succession, something absent from the single volume model.

Above a droplet radius of $r=0.9$ mm, the parameterised simulations show a statistically equal energy deposit to both the dynamically populated fastAerosol and the pre-populated fastAerosol using a p-value of $0.05$ (tested with Welch's t-test). Below $r=0.9$ mm (only two of the five such data points displayed) there is a statistically significant deviation indicating some difference in the simulations. Despite this, the deviation is on the order of $100$s of MeV, which is only around $1\%$ of the total energy deposited.

This discrepancy can likely be explained by the approximations used for determining $\sigma$ (only strictly correct for absolute distance) and for the approximation setting the upper limit in \cref{sec:absoluteClosest} (limit calculated in $2D$ but used for $3D$ geometry).

\begin{figure}
	\centering
	\includegraphics[width=0.75\linewidth]{net_plot}
	\caption{Plot of the energy deposited into the $150\times 150\times 50$ mm aluminum block at $(0,0,2625)$ mm by $5000$ protons ($50$ MeV) after travelling through a box-shaped aerosol ($n=0.01$ mm$^{-3}$). The smooth aerosol displays a drop off at $r=1$ mm which is absent from the other simulations, indicating a difference in underlying simulations.}
	\label{fig:net_plot}
\end{figure}

\iffalse

\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{droplet_histogram_plot}
	\caption{Plots of energy deposited into the $150\times 150\times 50$ mm aluminum block at $(0,0,2625)$ mm. Moving right: the radius increases from $0.1$ mm to $2$ mm; moving downwards: the aerosol build type varies from smooth to parameterised to fastAerosol. The energy deposited after the smooth aerosol is roughly independent on droplet radius (even though the effective density changes) which is drastically different from the parameterised and fastAerosol clouds for which energy deposited quickly drops off with increasing radius. Parameterised and fastAerosol clouds are qualitatively equivalent.}
	\label{fig:droplet_histogram_plot}
\end{figure*}

\fi

\subsection{Number Density Efficiency}
Total simulation time and maximum resident size for the parameterised, fastAerosol, and pre-populated fastAerosol simulations are shown in \cref{fig:efficiency_plot} for varying droplet number density. Note that, for both plots, error bars are typically too small to be visible.

\begin{figure}
    \centering
    \includegraphics[width=0.75\linewidth]{efficiency_plot}
    \caption{Plots of computational efficiency in simulating the transport of $5000$ protons ($50$ MeV) through a box-shaped aerosol with $r=1$ mm spherical droplets. Top: this plot of simulation time displays that, for all number densities, fastAerosol simulations are quicker than (or comparable to) parameterised simulations with significant speedup at large $n$. The kink in fastAerosol's simulation time at $n\approx 10^{-2.6}$ mm$^{-3}$ is likely due to early stopping of protons, reducing the number of generated droplets. Bottom: this plot of memory usage shows that fastAerosol simulations use less memory than (or comparable to) parameterised simulations for all number densities with a reduction of approximately one order of magnitude at large $n$.}
    \label{fig:efficiency_plot}
\end{figure}

\subsubsection{Timing}
The fastAerosol simulations are significantly faster than (or comparable to) the parameterised simulations for all number densities, with a pronounced speedup beyond number densities of $n=10^{-3}$ mm$^{-3}$. Since fastAerosol simulation times include the generation of the droplet positions while the parameterised simulations do not, the pre-populated fastAerosol simulation minus this population time is also plotted (green solid line).

Both the parameterised and pre-populated fastAerosol simulation times increase roughly polynomially with number density beyond $n=10^{-3}$ mm$^{-3}$ with roughly equal powers. The relatively constant $y$-gap of approximately $1$ order of magnitude indicates that parameterised simulations are roughly $10\times$ slower than fastAerosol simulations.

The dependence of fastAerosol's simulation time on droplet number density appears to largely be due to the population of droplets, as is seen in the large decrease in the pre-populated fastAerosol simulation time when the population time is subtracted. Despite this reduction, the simulation time seems to have the same power dependence on droplet number density, indicating that droplet population is not the sole timing load that is dependent on number density.

The drop-off in the dynamically populated fastAerosol's simulation time at $n\approx 10^{-2.6}$ mm$^{-3}$ is likely due to early termination of particles, reducing the effective aerosol volume and thus reducing the effective droplet count. This both provides evidence that droplet population is a large factor in the simulation load and it shows the value of the dynamically populated fastAerosol variant: it can be extremely quick compared to the pre-populated simulation models when the effective aerosol volume is significantly smaller than the bulk volume or when few primaries are shot. Specifically, at the largest number density, the dynamically populated fastAerosol's simulation time is approximately two orders of magnitude less than the parameterised simulation time.

\subsubsection{Memory}
For all number densities, the fastAerosol simulations use significantly less memory than (or comparable to) the parameterised simulations, peaking at a discrepancy of around one order of magnitude reductions at large number densities. With the current data, it is hard to determine the scaling of fastAerosol (pre-populated and dynamically populated) simulations at large number densities and whether these match the parameterised memory scaling

\subsection{Primary Count Efficiency}
Total simulation time for the parameterised, fastAerosol, and pre-populated fastAerosol simulations is shown in \cref{fig:efficiency_numPrimaries_plot} for varying primary count (error bars are too small to be visible).

\begin{figure}
    \centering
    \includegraphics[width=0.75\linewidth]{efficiency_numPrimaries_plot}
    \caption{Plots of simulation time after shooting $50$ MeV protons through a box-shaped aerosol ($n=0.01$ mm$^{-3}$; $r=1$ mm) for primary count ranging from $1$ to $10000$. Regardless of this count, the fastAerosol simulations are around $10\times$ (or more) quicker. Both parameterised and the pre-populated fastAerosol simulations appear relatively independent of the primary count while the dynamically populated fastAerosol seems approximately polynomially dependent on this count.}
    \label{fig:efficiency_numPrimaries_plot}
\end{figure}

The pre-populated fastAerosol maintains a speedup of approximately a factor of $10$ over parameterised simulations for all primary counts, indicating computational gains due to the simpler voxelization optimization utilized in the fastAerosol class (other speedups expected to change with particle count).

While the pre-populated fastAerosol simulation times are roughly independent of primary count (minor upwards slope not visible in plot), the dynamically populated fastAerosol simulation times are roughly polynomially dependent on primary count. This dependence is most likely due to droplet population since, at small primary counts, the dynamically populated fastAerosol can generate significantly fewer droplets than the pre-populated variant. Since these classes are otherwise identical, it is expected that both fastAerosol variants converge to roughly the same simulation time as the number of primaries becomes large enough for the dynamically populated fastAerosol to populate all voxels.

While not visible in the current plot, there is a positive dependence of simulation time on primary count for the pre-populated fastAerosol simulations and no such dependence for parameterised solids. It would be interesting to see how these dependences extend to extremely high primary counts, such as billions of primaries.

\subsection{Smartless Efficiency}
Total simulation time and maximum resident set size for the parameterised simulations is shown in \cref{fig:efficiency_smartless_plot} for varying `smartless' parameter, defined as the average number of voxels per contained geometry object (default is $2$). The fastAerosol (dynamically populated and pre-populated) average simulation time and memory load are shown as horizontal lines (dashed blue and dotted red, respectively).

\begin{figure}
    \centering
    \includegraphics[width=0.75\linewidth]{efficiency_smartless_plot}
    \caption{Plots of computational efficiency after shooting $50$ MeV protons through a box-shaped aerosol ($5000$ primaries; $n=0.01$ mm$^{-3}$; $r=1$ mm) for smartless parameter ranging from $2^{-14}$ to $2^5$. It is observed both that the default smartless value does not give a minimum in simulation time/memory and that, even at the true minimum simulation time, fastAerosol is more efficient.}
    \label{fig:efficiency_smartless_plot}
\end{figure}

Both the memory and simulation time appear dependent on the user-set smartless value, implying that the user-set value is consistently smaller than the calculated value for this aerosol (since the real smartless value is set as the minimum of these two; maximum user-set smartless is $32$). This suggests that the tuning algorithm is overridden by the user-set smartless value for this aerosol and thus an automatically tuned voxelization is not expected, contradicting the user manual's recommendation that no manual tuning is needed\cite{bookForApplicationDevelopers}. This need for additional manual tuning is potentially due to the exotic nature of the simulated geometry: the aerosol contains 12,500,000 droplets, significantly more sub-objects than parameterised typically takes.

The observed dependence of simulation time on smartless is expected: too fine a voxelization and the primary will undergo the overhead of travelling through many empty voxels; too coarse a voxelization and the distance to many droplets that are irrelevant to the primary will be calculated since they will all be lumped into a single voxel. Similarly, the observed memory dependence of parameterised simulations with smartless is as expected: increasing with smartless (since higher smartless implies more voxels) but with a minimum value representing the program's base memory draw.

Even when tuning smartless for minimum simulation time, to a value of approximately $2^{-2}$, both variants of fastAerosol (pre-populated and dynamically populated) maintain a simulation time reduction over the optimized parameterised solids. Similarly, both variants of fastAerosol maintains their more memory reductions over parameterised solids even when tuning smartless for minimum memory usage (even at smartless value $2^{-2}$, approximately where a minimum in simulation time occurs).

A more detailed investigation on tuning the smartless parameter should be done to further verify the computational gains of fastAerosol over parameterised solids for a variety of aerosols, since this test only verified fastAerosol's computational gains over parameterised aerosols for aerosols consisting of $n=0.01$ mm$^{-3}$ spherical droplets with $r=1$ mm. It would be valuable to learn if serious computational gains are allowed for other parameterised geometries by tuning this smartless parameter.

\subsection{Shape Flexibility}
The energy deposited is recorded after particle transport through a box-shaped aerosol with cylindrical, box-shaped, hemispherical, and spherical droplets (top of \cref{fig:shape_histogram_plot}; chosen to maximally fit in bounding sphere; no rotations applied). Additionally, the energy deposited is recorded after particle transport through a ellipsoid, pipe, cylinder, box-shaped aerosol with spherical droplets (bottom of \cref{fig:shape_histogram_plot}). From left to right, the histograms are placed in order of increasing filling ratio in the bounding sphere/box respectively.

\begin{figure*}
	\centering
	\includegraphics[width=\linewidth]{shape_histogram_plot}
	\caption{Plots of energy deposited into the $150\times 150\times 50$ mm aluminum block at $(0,0,2625)$ mm, ordered with increasing droplet/bulk filling ratio in the $r=1$ mm bounding radius and in the $150\times 150\times 5000$ mm bounding box moving right, respectively. The energy deposited decreases with increasing filling ratio for the droplet shape tests and for all bulk shape tests except for the pipe, which displays the largest energy deposited by far. This is expected since the pipe presents no barrier to most of the particles.}
	\label{fig:shape_histogram_plot}
\end{figure*}

With increasing filling ratio of the droplets in the bounding radius, the net energy deposited is observed to decrease as expected. Similarly, with increasing bulk filling ration in bounding box, the net energy deposited is observed to decrease except for the bulk pipe cloud. The exception of the pipe-shaped cloud is expected since the pipe-shaped cloud presents no barrier to the particles near $(x,y)=(0,0)$, where the majority of particles are shot.

\section{Conclusion}
Presented is a GEANT4 solid, fastAerosolSolid, and a helper class, fastAerosol, which enable accurate and efficient simulations of large aerosols. The fastAerosol classes represent aerosols as a collection of droplets, which we demonstrate is necessary when droplet size becomes large. When droplet number densities become large, the presented fastAerosol model is orders of magnitude more efficient in terms of timing and memory compared to the most efficient droplet-level method of building an aerosol currently included in GEANT4, using parameterised solids, as is seen in \cref{fig:efficiency_plot}.

The computational improvements in fastAerosol are achieved by fastAerosol's simpler voxelization optimization and in the allowed dynamical population of droplets. Additionally, the distance calculations are comparable in time to those performed in parameterised solids as is seen in \cref{fig:efficiency_numPrimaries_plot}: the pre-populated fastAerosol simulation times show a roughly $10\times$ speedup over parameterised simulations regardless of primary count and a higher primary count implies a larger number of distance calculations.

The reduction in memory usage is important for allowing simulations of a wide variety of experiments. The peak observed memory load in simulating particle transport through a parameterised cloud, $18.18\pm 0.01$ GB, is not small when compared to current computer hardware. Simulating the same scenario with the (dynamically populated) fastAerosol class reduces the memory load by a factor of $\sim17\times$ to $1.069\pm0.002$ GB. The parameter space of scenarios that can be simulated, then, is drastically larger when building aerosols with fastAerosol as compared to parameterised solids.

As this will be featured as an advanced example in GEANT4 (available as of version 10.7), it is our hope that the simplicity of implementation of fastAerosol along with the large computational gains will make this class attractive whenever an aerosol needs to be simulated.

\subsection{Future work}
There are some areas for future study to be investigated. First, the periodic droplet bunching apparently at voxel boundaries observed in \cref{fig:distribution_plot} should be studied, maybe by modifying the the order (\cref{eq:order}) in which voxels are filled. Second, the study of the approximations (use of $\sigma$ for vectorized distance and the upper limit on spherical voxel radius in \cref{sec:absoluteClosest}) should be done in hopes of removing the discrepancy in energy deposited after particles transport through aerosols containing droplets with radius $r<0.9$ mm. Third, a more detailed investigation of tuning the smartless parameter for parameterised solids should be done, as that appears to have capacity for large computational gains and lessons beyond the scope of simulating aerosols. It appears possible to disable parameterised voxelization optimization completely, so re-doing tests with this optimization disabled would be interesting in parallel to the smartless investigations to determine if this optimization really does hurt the efficiency.

There are features to add for more sophisticated simulations such as randomized droplet sizes, shapes, and materials. This would require alteration of both the droplet placement and search algorithm, but it seems possible.

Another feature extension that does not seem difficult is to allow a user-specified collision detection function. This would allow accurate placement of non-spherical droplets at the cost of slower simulations. In a similar vein, extending fastAerosol to accept user-specified droplet locations could be valuable, but this would be incompatible with dynamic population.

\subsection{Acknowledgments}
This work built on the GEANT4 source code and examples, so we thank those developers and maintainers. In particular, Dr. Makoto Asai, Dr. Gabriele Cosmo, and  Dr. Susanna Guatelli provided invaluable discussions, recommendations, testing, and assistance with integration with the GEANT4 codebase.

\bibliographystyle{elsarticle-num}


\section{Introduction}
The $\alpha$-cluster structure established for different excited states in several
light nuclei like $^{12}$C or $^{16}$O has inspired numerous experimental and
theoretical studies, especially, the direct nuclear reactions measured with $^{12}$C
as projectile and/or target \cite{Freer07}. Given the cluster states above the
$\alpha$-decay threshold of $^{12}$C, some direct reactions with $^{12}$C were shown
to produce both the free $\alpha$ particle and unstable $^{8}$Be in the exit channel
\cite{Munoz10,Soy12,Curtis16}. Consequently, the knowledge of the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\
interaction potentials should be important for the studies of such reactions within
the distorted wave Born approximation (DWBA) or coupled reaction channel (CRC)
formalism.

Given a well established 2$\alpha$-cluster structure of the unbound $^{8}$Be nucleus,
the $^{8}$Be+$^{8}$Be\  interaction potential poses a four-body problem which is a challenge for
the determination of the $^{8}$Be+$^{8}$Be\  optical potential (OP) that cannot be deduced from
a standard optical model (OM) analysis because of the lack of the elastic $^{8}$Be+$^{8}$Be\
scattering data. The knowledge about the $\alpha+^{8}$Be\  and $^{8}$Be-nucleus\  OP's should be also important
for the studies of those direct reaction processes that produce $^8$Be fragments in
the exit channel \cite{Rom09,PSM16,PSM19}. Although the $\alpha$-nucleus\  and nucleus-nucleus\  OP's are proven
to be well described by the double-folding model (DFM) using the accurate ground-state
densities of the colliding nuclei and a realistic density dependent nucleon-nucleon
(NN) interaction (see, e.g., Refs.~\cite{Sat79,Kho97,Kho00,Kho01,Kho16}), the DFM
cannot be used to calculate the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  OP's because of a strongly deformed,
extended two-center density distribution of $^{8}$Be. In general, one could think
of the triple- and quadruple folding models for the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  potentials,
respectively, but these will surely be complicated and involve much more tedious
calculation in comparison with the standard DFM method. Although some phenomenological
OP's are available in the literature for $^7$Li and $^9$Be, two nuclei neighboring
$^8$Be, a strongly (two-center) deformation of the unstable $^8$Be nucleus casts
doubt on the extrapolated use of these potentials for the $^{8}$Be+$^{8}$Be\  and $^{8}$Be-nucleus\  systems.

Given a very loose (unbound) $^{8}$Be nucleus that breaks up promptly into 2 $\alpha$
particles, we determine in the present work the $^{8}$Be+$^{8}$Be\  OP using the Continuum-Discretized
Coupled Channel (CDCC) method which was developed to take into account explicitly
the breakup of the projectile and/or target. A textbook example is a direct reaction
induced by deuteron, which is loosely bound and can be, therefore, easily broken up
into a pair of free proton and neutron. Originally, the deuteron breakup states
were included in terms of a discretized continuum by the CDCC method (see, e.g.,
Refs.~\cite{KYI86,AIK87,YOM12} for reviews).
In the recent version of the CDCC theory, the continuum of deuteron is approximated
by the square-integrable functions corresponding to positive energies. As a result,
this approach can be well extended to study elastic scattering of exotic nuclei
that have rather low breakup energies (typical examples are $^6$He and $^{11}$Be).

The first developments of the CDCC method were done in the framework of a three-body
system where the projectile is seen as a two-body nucleus and target is assumed
to be structureless, being in its ground state. More recently, four-body calculations
were developed, for either a three-body projectile on a structureless target \cite{MHO04},
or a two-body projectile and a two-body target \cite{De18}. The latter approach
is highly time consuming, but was successfully applied to study $^{11}{\rm Be}+d$
scattering in terms of $^{11}{\rm Be}=^{10}$Be$+n$ and $d=p+n$. The goal of the present
study is to determine the $^{8}$Be+$^{8}$Be\  OP based on the elastic scattering matrix predicted
by the 4-body CDCC calculation of four interacting $\alpha$ clusters. While such a
$4\alpha$ model is not appropriate for the spectroscopy of $^{16}$O \cite{De18},
the derived OP for elastic $^{8}$Be+$^{8}$Be\  scattering is expected to be reliable. The only input
for the present $4\alpha$ CDCC calculation is a realistic $\alpha+\alpha $\  interaction potential.

Although $^8$Be is particle-unstable, its half life around $10^{-16}$s is long enough
for the $^8$Be-nucleus OP to contribute significantly to a direct reaction that
produces $^8$Be in the exit channel, like the $\alpha$ transfer
$^{12}$C($\alpha$,$^8$Be) reaction. This particular reaction was shown to be a good
tool for the study of the high-lying or resonance states of $^{16}$O \cite{Che67,Mar72}
and to determine the $\alpha$ cluster configurations of this nucleus
\cite{Che67,Bro76,Cur13}. Because of the unbound structure of $^8$Be, the direct
reaction reactions $A(\alpha$,$^8$Be)$B$ usually have a very low cross section (of a
few tens microbarn), but they are extremely helpful for the study of the $\alpha$-cluster
structure of the target nuclei \cite{Woz76}. In particular, the $\alpha$ spectroscopic
factors of different cluster states were deduced from these measurements at the
$\alpha$ incident energies of 65 to 72.5 MeV.

In the present work, the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  optical potentials deduced from the scattering
wave functions given by the 3- and $4\alpha$ CDCC calculations are used as the core-core
and the exit OP, respectively, in the CRC study of $^{12}$C($\alpha$,$^8$Be) reaction
measured at 65 MeV \cite{Woz76}. The OP of the entrance channel is calculated in the
DFM using the density dependent CDM3Y6 interaction that was well tested in the
mean-field studies of nuclear matter as well as in the OM studies of the elastic $\alpha$-nucleus\
scattering \cite{Kho97,Kho01}, and it accounts well for the elastic $\alpha+^{12}$C\  scattering
data measured at 65 MeV \cite{Yasue83,Gon14}. The $\alpha$ spectroscopic factors
of the $\alpha+^{8}$Be cluster configurations of $^{12}$C are taken from the results
of the complex scaling method (CSM) by Kurokawa and Kato \cite{KaKu07}.

\section{Three- and four-body CDCC methods}
\label{sec2}
We discuss here the CDCC method used to determine the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  optical
potentials, where the $\alpha$ particles are treated as structureless and interacting
with each other through a (real) potential $v_{\alpha\alpha}(r)$. The Hamiltonian
of the $\mbox{$\alpha+\alpha$}$ system is given by
\begin{equation}
 H_{\alpha\alpha}(\bm{r})=T_r+v_{\alpha\alpha}(r),
\label{eq1}
\end{equation}
where $T_r$ is the relative kinetic energy. There are two versions of the $\alpha+\alpha $\  potential
\cite{AB66,BFW77} parametrized in terms of Gaussians amenable for the present CDCC
calculation. In the present work, we have chosen the $\alpha+\alpha $\  potential suggested by Ali
and Bodmer \cite{AB66} (referred to hereafter as AB potential). The AB potential
simulates the Pauli blocking effect by a repulsive core that makes this potential
much shallower than the deep $\alpha+\alpha $\  potential suggested by Buck {\it et al.} \cite{BFW77}.
Both potentials reproduce equally well the $\alpha+\alpha $\  phase shifts, and they were shown
by Baye \cite{Ba87} to be linked by a supersymmetric transformation. The Buck potential,
however, contains some deeply-bound states (two states with $\ell=0$ and one with
$\ell=2$), which do not have physical meaning but simulate the so-called Pauli forbidden
states \cite{Sa69} in the microscopic $\alpha+\alpha$ model. As long as the two-body
$\alpha+\alpha$ system is considered, the choice of either AB or Buck potential is not
crucial. However, when dealing with more than two $\alpha$ clusters, the forbidden
states have to be removed as they produce spurious states in a multi-cluster system
like $\alpha+^{8}$Be\  or $^{8}$Be+$^{8}$Be. There are two methods to remove the forbidden states
in the multi-cluster systems: either to apply the pseudostate  method \cite{Ne71} or
to use the supersymmetric transformation \cite{Ba87}. These two techniques, however,
give rise to a strong angular-momentum dependence of the $\alpha+\alpha $\  potential that cannot be used
in most of the multi-cluster models. Among the 3$\alpha$ models, only the hypersperical
method and Faddeev method are able to use the deep $\alpha+\alpha $\  potential with an exact
removal of the $\alpha+\alpha $\  forbidden states. This is why other studies of the 3$\alpha$
and 4$\alpha$ systems \cite{Og09,De18,Su02,Ti08} have used only the $\ell$-independent
AB potential. In the present work, we perform the CDCC calculation of the $\alpha+^{8}$Be\  and
$^{8}$Be+$^{8}$Be\  optical potentials using the AB potential of the $\alpha+\alpha $\  interaction, so that
the spurious effects arising from the Pauli forbidden states can be avoided.

The present CDCC method is based on the eigenstates $\Phi_{\lambda}^{\ell m}(\bm{r})$
of the Hamiltonian (\ref{eq1}) which can be written as
\begin{align}
\Phi_{\lambda}^{\ell m}({\bm r})=\frac{1}{r}{u_{\lambda}^{\ell}(r)}
 Y_{\ell m}(\hat{\bm r}), \nonumber
\end{align}
where $\ell$ is the relative orbital momentum of the $\alpha+\alpha $\  system. The radial wave
functions $u_{\lambda}^{\ell}(r)$ of the two-$\alpha$ state $\lambda$ are expanded over
a basis of $N$ orthonormal functions $\varphi_i(r)$
\begin{align}
 u_{\lambda}^{\ell}(r)=\sum_{i=1}^N f_{\lambda,i}^{\ell}\varphi_i(r),
\label{eq2}
\end{align}
where $f_{\lambda,i}^{\ell}$ are determined by diagonalizing the eigenvalue problem
\begin{align}
\sum_j f_{\lambda,j}^{\ell}\bigl(\langle\varphi_i|H_{\alpha\alpha}|\varphi_j\rangle
 - E_{\lambda}^{\ell}\delta_{ij}\bigr)=0. \label{eq3}
\end{align}
The eigenvalues with $E_{\lambda}^{\ell}<0$ correspond to the physical bound states,
while those with $E_{\lambda}^{\ell}>0$ are referred to as the pseudostates,
which are used in the present CDCC method to simulate the breakup of $^8$Be.
Note that there is only a small number of physical states in a CDCC calculation
(often one for exotic nuclei). Although $^8$Be is unbound, its energy is very
close to the $\alpha+\alpha $\  threshold, and the lifetime is long enough to use a quasi-bound
approximation for the ground state (g.s.). Equations (\ref{eq2}) and (\ref{eq3})
are general for any choice of the basis functions $\varphi_i(r)$. We use
here a Lagrange-mesh basis \cite{Ba15} derived from the Legendre polynomials, and
the calculation of matrix elements in Eq.~(\ref{eq3}) is fast and accurate. We refer
the reader to Ref.~\cite{Ba15} for more details and the application of the
Lagrange-mesh basis.

The $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems are described by the 3$\alpha$ and 4$\alpha$ Hamiltonians,
respectively, as
\begin{align}
& H_3=H_{\alpha\alpha}(\pmb{r}_1)+T_R+\sum_{i=1}^2 v_{\alpha\alpha}(S_{1i}) \nonumber \\
& H_4=H_{\alpha\alpha}(\pmb{r}_1)+H_{\alpha\alpha}(\pmb{r}_2)+
T_R+\sum_{i,j=1}^2 v_{\alpha\alpha}(\vert \pmb{S}_{1i}-\pmb{S}_{2j}\vert),
\label{eq4}
\end{align}
where $\bm{R}$ is the projectile-target coordinate and $(\bm{r}_1,\bm{r}_2)$ are
the internal coordinates of the $^8$Be nuclei, as illustrated in Fig.~\ref{f1}.
Coordinates $\pmb{S}_{1i}$ and $\pmb{S}_{2j}$ are expressed as a function
of $(\pmb{R},\pmb{r}_1)$ and of $(\pmb{R},\pmb{r}_2)$, respectively.
\begin{figure}[htb]
	\begin{center}
		\epsfig{file=config.eps,width=7.5cm}
		\caption{Alpha configurations and coordinates used for $\alpha+^{8}$Be\  (a) and for $^{8}$Be+$^{8}$Be\  (b).}
		\label{f1}
	\end{center}
\end{figure}

In the CDCC approximation, the $\alpha+^{8}$Be\  wave function for a total angular moment $J$
and parity $\pi$ is written as
\begin{align}
\Psi^{JM\pi}(\bm{R},\bm{r}_1)=\sum_{\lambda\ell L}g^{J\pi}_{\lambda \ell L}(R)
\bigl[\Phi_{\lambda}^{\ell}(\bm{r}_1)\otimes Y_L(\hat{\bm R})\bigr]_{JM}, \label{eq5}
\end{align}
where the summation over the pseudostates $\lambda$ is truncated at a given energy
$\mbox{$E_{\rm max}$}$. The $^{8}$Be+$^{8}$Be\  wave functions involve 4 $\alpha$ clusters and can be expressed as
\begin{align}
&\Psi^{JM\pi}(\bm{R},\bm{r}_1,\bm{r}_2)=\sum_{\lambda_1 \ell_1}\sum_{\lambda_2\ell_2}
 \sum_{IL}g^{J\pi}_{\lambda_1 \ell_1 \lambda_2 \ell_2 I L}(R)\nonumber \\
&\hspace*{1cm} \times \biggl[\bigl[\Phi_{\lambda_1}^{\ell_1}(\bm{r}_1)
\otimes\Phi_{\lambda_2}^{\ell_2}(\bm{r}_2)\bigr]_I
\otimes Y_L(\hat{\bm R}) \biggr]_{JM},
\label{eq6}
\end{align}
where the parity conservation imposes $(-1)^L=(-1)^I$. There are two parameters defining
the CDCC basis: the maximum $^8$Be angular momentum $\mbox{$\ell_{\rm max}$}$ and maximum pseudostate
energy $\mbox{$E_{\rm max}$}$. The physical quantities obtained from the CDCC calculation (scattering
matrix, elastic cross section, and local equivalent OP) must be converged with respect
to these parameters. In practice, the 4$\alpha$ calculations involve many channels and
are, therefore, highly time consuming \cite{De18}.

The relative radial wave functions $\chi^{J\pi}_c(R)$, with the indices
$c=(\lambda\ell L)$ for $\alpha+^{8}$Be\  and $c=(\lambda_1 \ell_1 \lambda_2 \ell_2 I L)$ for $^{8}$Be+$^{8}$Be\ ,
are obtained from the solutions of the following coupled-channel equations
\begin{align}
&\biggl[-\frac{\hbar^2}{2\mu}\biggl(\frac{d^2}{dR^2}-\frac{L(L+1)}{R^2}  \biggr)
+E_{c_1}+E_{c_2}-E_{\rm c.m.} \biggr]\chi^{J\pi}_{c}(R)\nonumber \\
&\hspace*{1cm} +\sum_{c'}V^{J\pi}_{cc'}(R)\chi^{J\pi}_{c'}(R)=0,
\label{eq7}
\end{align}
where $\mu$ is the reduced mass, $E_{\rm c.m.}$ is the center-of-mass energy,
$E_{c_1}$ and $E_{c_2}$ are the excitation energies of the two interacting nuclei,
separated by the distance $R$ as shown in Fig.~\ref{f1}. The coupling potentials
$V^{J\pi}_{cc'}(R)$ are determined by the method explained in Refs.~\cite{RAG08,De18}.
The system of the coupled channel equations (\ref{eq7}) is solved using the $R$-matrix
method which provides explicitly the scattering matrix and the associated wave
functions \cite{DB10,De16a}. Although the AB potential \cite{AB66} is real, it
consistently reproduces the experimental $\alpha+\alpha $\  phase shifts up to about 20 MeV.
In the present CDCC approach, the loss of flux from the elastic scattering channel
is due entirely to the breakup channels, and the local equivalent $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\
optical potentials are therefore complex. Owing to the strong 2$\alpha$ structure
of $^8$Be it is likely that these breakup channels represent the main source
of the absorption.

\section{Results and discussion}
\label{sec3}
\subsection{Local equivalent OP for the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems}
\label{sec3.1}
The main goal of our study is to determine the local ($J$-independent) equivalent
optical potential $U$ for the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems at the considered energies,
based on the scattering wave functions given by the solutions of the CDCC equations
(\ref{eq7}). The main requirement for this procedure is that the solutions
$\chi^{J\pi}$ of the one-channel OM equation with the optical potential $U$
\begin{align}
\bigl[E_{\rm c.m.}-T_R-U(R)\bigr]\chi^{J\pi}(R)=0 \label{eq8}
\end{align}
give the cross section of the elastic $\alpha+^{8}$Be\  or $^{8}$Be+$^{8}$Be\  scattering close to that
given by the 3$\alpha$ or 4$\alpha$ CDCC calculation (\ref{eq7}), especially,
the cross section at forward angles which is sensitive to the surface part
of the 3$\alpha$ or 4$\alpha$ interaction potential. We briefly discuss the
two approaches used in the present work for this purpose.
\begin{figure}[h]\vspace*{-1cm}
	\includegraphics[width=1.0\textwidth]{fig2.eps}\vspace*{-1cm}
\caption{Complex OP for the elastic $^{8}$Be+$^{8}$Be\  scattering at $E_{\rm c.m.}=41.3$ MeV
given by the method (i), and that in the WS form (\ref{eq9}) given by the method
(ii).}  \label{f2}
\end{figure}

(i) The quantum-mechanically consistent method using the matrix inversion was suggested
in Refs.~\cite{TNL89,De18} to derive a local equivalent potential (LEP) that exactly
reproduces the elastic cross section given by the CDCC calculation (\ref{eq7}).
However, this LEP has two major drawbacks that prevent its further use in the direct
nuclear reaction calculation. Namely, the derived LEP strongly depends on the total
angular momentum $J$, and its radial dependence has the singularities caused by the
nodes of the scattering wave functions. These problems can be handled by the method
proposed in Ref.~\cite{TNL89} which averages the obtained LEP over the angular momenta
to obtain a smooth $J$-independent OP without discontinuity that approximately
reproduces the CDCC elastic scattering cross section. The recent 4-body CDCC
calculation \cite{De18} has shown that such an averaging method to a fairly
good approximation determines the local $J$-independent OP. The complex OP
derived using this approach is denoted hereafter as $U_{\rm LEP}$, with its imaginary
part $W_{\rm LEP}$ originating from the breakup channels included in the CDCC
calculation (\ref{eq7}). We have first performed the CDCC calculation (\ref{eq7})
for the elastic $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  scattering at $E_{\rm c.m.}=43.3$ and 41.3 MeV,
respectively, using the AB potential \cite{AB66} for the $\alpha+\alpha $\  interaction.
The maximum angular momentum of the $\alpha+\alpha $\  system is $\mbox{$\ell_{\rm max}$}=2$,
and the maximum pseudostate energy is $\mbox{$E_{\rm max}$}=10$ MeV. These cutoff values were
well tested to ensure the convergence of both the elastic cross section and
$U_{\rm LEP}$. The complex $U_{\rm LEP}$ for the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems were
obtained first in a Lagrange mesh \cite{TNL89}, and then interpolated into the smooth
shapes for use as the external input of the complex OP in the CRC calculation.
\begin{figure}[h]\vspace*{0cm}
 \includegraphics[width=0.75\textwidth]{fig3.eps}\vspace*{-0.5cm}
 \caption{The CDCC prediction for the elastic $^{8}$Be+$^{8}$Be\  (upper panel) and $\alpha+^{8}$Be\  (lower
 panel) scattering at $E_{\rm c.m.}=41.3$ and 43.3 MeV, in comparison  with the
 corresponding OM results given by the optical potentials $U_{\rm LEP}$  and
 $U_{\rm WSD}$ determined by the methods (i) and (ii), respectively.}	\label{f3}
\end{figure}

(ii) The standard OM method can also be used to determine a phenomenological
$J$-independent OP in the conventional Woods-Saxon (WS) form, with its parameters
adjusted to obtain a good OM fit to the elastic cross section given by the 3$\alpha$
or 4$\alpha$ CDCC calculation. We have assumed in the present work the following
(volume+surface) WS form for the complex OP of the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems at the
energies under study, denoted hereafter as $U_{\rm WSD}$
\begin{align}
- & U_{\rm WSD}(R)=V_vf_v(R)-4V_da_d\frac{df_d(R)}{dR}+i\left[W_vf_w(R)
 -4W_da_s\frac{df_s(R)}{dR}\right], \nonumber\\
 & {\rm where}\ f_x(R)=\frac{1}{1+\exp[(R-R_x)/a_x]},\ x=v,d,w,s.
\label{eq9}
\end{align}
The elastic $^{8}$Be+$^{8}$Be\  and $\alpha+^{8}$Be\  cross sections at $E_{\rm c.m.}=41.3$ and 43.3 MeV,
predicted by the 4$\alpha$ and 3$\alpha$ CDCC calculations (\ref{eq7}), respectively,
have been used in the method (ii) as the ``experimental data" with the uniform 10\%
uncertainties for the OM analysis to determine the WS parameters of $U_{\rm WSD}$.
We found quite a shallow WS potential that gives a good agreement of the OM result
with the CDCC elastic cross section (see the OP parameters in Table~\ref{t1}).
\begin{table}[h]
\setlength{\tabcolsep}{0.3em}
\centering
\caption{WS parameters (\ref{eq9}) of $U_{\rm WSD}$ given by the method (ii) based
 on the OM fit to the elastic $^{8}$Be+$^{8}$Be\  and $\alpha+^{8}$Be\  cross sections at $E_{\rm c.m.}=41.3$
and 43.3 MeV, predicted by the 4$\alpha$ and 3$\alpha$ CDCC calculations
(\ref{eq7}), respectively. $J_V$ and $J_W$ are the volume integrals per interacting
nucleon pair of Re~$U_{\rm WSD}$ and Im~$U_{\rm WSD}$, respectively. $J_{V_{\rm LEP}}$
and $J_{W_{\rm LEP}}$ are those of $U_{\rm LEP}$ given by the method (i).}
\label{t1}\vspace{0.5cm}
\begin{tabular}{ccccccccc}
 \hline\hline
 & $V_v(W_v)$ & $R_{v(w)}$ & $a_{v(w)}$ & $V_d(W_d)$ & $R_{d(s)}$ & $a_{d(s)}$ &
 $-J_V(J_W)$ & $-J_{V_{\rm LEP}}(J_{W_{\rm LEP}})$\\
 & (MeV) & (fm) & (fm) & (MeV) & (fm) & (fm) & (MeV~fm$^3$) & (MeV~fm$^3$) \\ \hline
 \multicolumn{9}{c}{$^8$Be+$^8$Be, $E_{\rm c.m.}=41.3$ MeV} \\ \hline
  Real & 10.42 & 3.652 & 0.190 & 5.785  & 4.811  & 0.562 & 95.88 & 100.0 \\
  Imag. & 1.486 & 7.225 & 0.183 & 0.567  & 2.665  & 1.564 & 46.87 & 29.71 \\ \hline
    \multicolumn{9}{c}{$^4$He+$^8$Be, $E_{\rm c.m.}=43.3$ MeV}  \\ \hline
  Real  & 1.535 & 5.182 & 0.123 & 5.397 & 2.584 & 0.997 & 111.0 & 97.62 \\
  Imag. & 10.09 & 1.771 & 0.125 & 8.201 & 3.299 & 0.120 & 24.56 & 33.95 \\
	\hline\hline
	\end{tabular}
\end{table}

The radial shapes of both $U_{\rm LEP}$ and $U_{\rm WSD}$ potentials for the $^{8}$Be+$^{8}$Be\
system at $E_{\rm c.m.}=41.3$ MeV are shown in Fig.~\ref{f2}, the use of the shallow
AB potential of the $\alpha+\alpha $\  interaction is shown to result on quite a shallow potential
$U_{\rm LEP}$. A moderate oscillation of $U_{\rm LEP}$ is seen at small radii that
might originate from the $J$-dependence of the exact LEP discussed above. The best-fit
WS complex OP determined by the method (ii) has the strength of Re~$U_{\rm WSD}$
enhanced slightly at the surface ($R\approx 5$ fm), and a weak and smooth
Im~$U_{\rm WSD}$. One can see in Table~\ref{t1} that the volume integrals of
Re~$U_{\rm WSD}$ and Im~$U_{\rm WSD}$ are close to those of Re~$U_{\rm LEP}$ and
Im~$U_{\rm LEP}$, which indicates that the OP's given by both methods belong to about
the same potential family. The results of the CDCC calculation (\ref{eq7}) for
the elastic $^{8}$Be+$^{8}$Be\  and $\alpha+^{8}$Be\  scattering at $E_{\rm c.m.}=41.3$ and 43.3 MeV, respectively,
are compared in Fig.~\ref{f3} with the results of the OM calculation (\ref{eq8})
using $U_{\rm LEP}$ and $U_{\rm WSD}$. The OM results given by both OP's agree
fairly good with the CDCC prediction at forward angles, while at medium and large angles
the phenomenological $U_{\rm WSD}$ determined by the method (ii) better reproduces
the CDCC cross sections. The agreement with the CDCC results becomes worse at large
angles, and it might be due to the nonlocality effects.

We note further that the method (i) fails to derive a smooth $\ell$-independent
$U_{\rm LEP}$ based on the CDCC results obtained with the deep Buck $\alpha+\alpha $\  potential.
Namely, the obtained $\ell$-independent $U_{\rm LEP}$ turns out to be deeper but
strongly oscillatory, and it gives the elastic cross section substantially different
from that given by the CDCC calculation. Such a failure of the $\ell$-independent
$U_{\rm LEP}$ based on the Buck potential is presumably caused by the Pauli forbidden
states, and this remains an unsolved problem for the present 4$\alpha$ CDCC method.
Therefore, we deem hereafter reliable only the CRC results obtained with $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\
OP's derived based on the CDCC elastic cross section obtained with the AB potential
of the $\alpha+\alpha $\  interaction.

\subsection{CRC study of the $^{12}$C($\alpha,^8$Be) reaction}
\label{sec3.2}
The $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  optical potentials determined by the methods (i) and (ii) have been
further used as the potential inputs for the CRC study of the $\alpha$ transfer
$^{12}$C($\alpha$,$^8$Be) reaction measured at $E_\alpha=65$ MeV \cite{Woz76}.
We briefly recall the multichannel CRC formalism, to illustrate how the OP's of the
$\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  systems enter the CRC calculation of the $\alpha$ transfer cross section.
In general, the CRC equation for the initial channel $\beta$ of the transfer reaction
can be written as \cite{Sat83,ThompsonN}
\begin{eqnarray}
&&\hskip 1cm \bigl[E_{\beta}-T_{\beta}-U_{\beta}({\bm R})\bigr]
 \chi_{\beta}({\bm R})= \nonumber\\
&&\sum_{\beta'\neq\beta}
 \bigl\{\langle\beta|W|\beta'\rangle +\langle\beta|\beta'\rangle
 \bigl[T_{\beta'}+U_{\beta'}({\bm R}')-E_{\beta'}\bigr]\bigr\}
 \chi_{\beta'}({\bm R}'). \label{eq10}
\end{eqnarray}
\begin{figure}[bht]\vspace*{-2cm}
 \hspace*{-1cm}\includegraphics[width=0.8\textwidth]{fig4.eps}\vspace*{-2.0cm}
 \caption{Cluster configurations in the entrance and exit channels of the
 $^{12}$C($\alpha,^8$Be) reaction, and the corresponding coordinates used
 for the inputs of the potentials in the CRC calculation.}
\label{f4}
\end{figure}
Without coupling to the inelatic scattering channels, the indices $\beta$ and $\beta'$
in Eq.~(\ref{eq10}) stand for the initial $\alpha+^{12}$C\  and final $^{8}$Be+$^{8}$Be\  partitions of the $\alpha$
transfer reaction, respectively, as shown in Fig.~\ref{f4}. In the present CRC analysis,
the index $\beta'$ is used to identify both the $^{8}$Be+$^{8}$Be$_{\rm g.s.}$ and
$^{8}$Be+$^{8}$Be$^*_{2^+}$ exit channels of the final partition. The distorted waves
$\chi_{\beta}$ and $\chi_{\beta'}$ are given by the optical potentials $U_{\beta}$
and $U_{\beta'}$ of the $\alpha+^{12}$C\  and $^{8}$Be+$^{8}$Be\  systems, respectively. The $\alpha$ transfer
proceeds via the transfer interaction $W$ which is determined in the post form
\cite{Sat83,ThompsonN} as
\begin{equation}
 W=V_{\alpha - ^{12}\text{C}}(\bm{r})+\bigl[U_{\alpha+^8\text{Be}}(\bm{R}_{\rm cc})
 -U_{^8\text{Be}+^8\text{Be}}(\bm{R}')\bigr], \label{eq11}
\end{equation}
with the radii of the potentials illustrated in Fig.~\ref{f4}. Here
$V_{\alpha-^{12}\text{C}}(r)$ is the potential binding the $\alpha$ cluster
to the $^8$Be core in the g.s. of $^{12}$C. The difference between the core-core OP
and that of the final partition,
$U_{\alpha+^8\text{Be}}(R_{\rm cc})-U_{^8\text{Be}+^8\text{Be}}(R')$,
is the complex remnant term of $W$. The CRC equations (\ref{eq10}) are solved
iteratively using the code FRESCO written by Thompson \cite{Thompson}, with the
complex (nonlocal) remnant term and boson symmetry of the identical $^{8}$Be+$^{8}$Be\
system properly taken into account. One can see that the $^{8}$Be+$^{8}$Be\  OP enters the CRC
calculation of the $\alpha$ transfer $^{12}$C($\alpha$,$^8$Be) reaction as the input
of both the remnant term and the OP of the final partition. Therefore, it can be
tested indirectly based on the CRC description of the $\alpha$ transfer data.

\begin{figure}[bht]\vspace*{-0.5cm}
 \includegraphics[width=\textwidth]{fig5.eps}\vspace*{-0.5cm}
 \caption{CRC description of the elastic $\alpha+^{12}$C  scattering at $E_{\alpha}=65$
 MeV given by the (unrenormalized) real folded OP and imaginary OP chosen
 in the WS form \cite{Kho01}, in comparison with the data taken from
 Refs.~\cite{Yasue83,Gon14}} \label{f5}
\end{figure}
The OP of the initial partition $U_{\alpha+^{12}\text{C}}$ has its real part given
by the double-folding model using the density dependent CDM3Y6 interaction \cite{Kho01},
and imaginary part chosen in the WS shape, with the parameters fine tuned to the
best CRC fit of the elastic $\alpha+^{12}$C\  scattering data measured at $E_\alpha=65$ MeV
\cite{Yasue83,Gon14}. A reasonably good CRC description of the elastic $\alpha+^{12}$C  scattering
data at $E_{\alpha}=65$ MeV (see Fig.~\ref{f5}) is achieved without renormalizing
the strength of the real OP. In principle, we could also think of using the
4$\alpha$ CDCC method to predict the $\alpha+^{12}$C  OP at the considered energy. However,
Suzuki {\it et al.} \cite{Suzuki08} have shown that the use of a \emph{local} $\alpha+\alpha $\  potential
(that properly reproduces the experimental $\alpha+\alpha $\  phase shifts) cannot provide a proper
3$\alpha$ description of both the g.s. and $0^+_2$ excitation (known as Hoyle state)
of $^{12}$C. This problem could only be solved by introducing a microscopically
founded \emph{nonlocal} $\alpha+\alpha $\  force that mimics the interchange of three $\alpha$ clusters
in the phase space allowed by the Pauli principle \cite{Suzuki08}. The use a nonlocal
$\alpha+\alpha $\  interaction remains beyond the scope of the present 4-body CDCC method \cite{De18}.
On the other hand, the elastic $\alpha+^{12}$C  scattering at energies above 10 MeV/nucleon
is proven to be strongly refractive \cite{Kho97,Kho01,Kho07r}, with a far-side
dominant elastic cross section at large angles typical for the nuclear rainbow,
which can be well described by the deep (mean-field type) real OP predicted
by the double-folding model \cite{Kho01,Kho07r}.

For the $\alpha$ transfer reaction, the initial (internal) state of the $\alpha$
cluster bound in $^{12}$C is assumed to be $1s$ state. Then, the relative-motion wave
function $\Phi_{NL}(\bm{r})$ of the $\alpha$+$^{8}$Be configuration in the $^{12}$C
target ($L$-wave state) has the number of radial nodes $N$ determined by the Wildermuth
condition \cite{Sat83}, so that the total number of the oscillator quanta
$\mathcal{N}$ is conserved
\begin{equation}
 \mathcal{N}=2(N-1)+L=\sum_{i=1}^{4}2(n_i-1)+l_i, \label{eq12}
\end{equation}
where $n_i$ and $l_i$ are the principal quantum number and orbital momentum of each
constituent nucleon in the $\alpha$ cluster. $\Phi_{NL}(\bm{r})$ is obtained in the
potential model using $V_{\alpha - ^{12}\text{C}}(r)$ chosen in the WS shape, with its
radius and diffuseness fixed as $R=3.767$ fm and $a=0.65$ fm, and the WS depth
($V=51.6$ MeV) adjusted to reproduce the $\alpha$ separation energy of $^{12}$C.
Because the ground state of $^8$Be is unbound by 92 keV, we have used in the present
CRC calculation a {\it quasi-bound} approximation for $^8$Be similar to that used for
the g.s. of $^8$Be in the CDCC calculation as discussed in Sec.~\ref{sec2}, to
describe the formation of $^8$Be on the exit channel of the $\alpha$ transfer reaction.
For this purpose, the repulsive core of the AB potential was slightly weakened
to give the ($1s$) state $\Phi_\alpha(\bm{r}')$ of the $\alpha$ cluster in $^8$Be
a quasi-binding energy of 0.01 MeV. The cluster wave functions $\Phi_{NL}(\bm{r})$
and $\Phi_\alpha(\bm{r}')$ are used explicitly in the calculation
of the complex nonlocal $\alpha$ transfer form factor
\begin{equation}
\langle\beta'|W|\beta\rangle\sim \langle[\Phi_\alpha(\bm{r}')\otimes
 Y_{L_{\beta'}}(\hat{{\bm R}'})]_{J_{\beta'}}|W|[\Phi_{NL}(\bm{r})\otimes
 Y_{L_{\beta}}(\hat{\bm R})]_{J_{\beta}}\rangle, \label{eq13}
\end{equation}
where $L_{\beta}$ and $L_{\beta'}$ are the relative orbital momenta of the
initial and final partitions. The wave functions of $^8$Be core in the initial
and $^4$He core in the final partitions are omitted in (\ref{eq13}) because they
are spectators and do not contribute to the transfer \cite{ThompsonN}.

The CRC calculation of the $\alpha$ transfer $^{12}$C($\alpha$,$^8$Be) reaction
requires the input of the spectroscopic factor of the $\alpha$ cluster in $^{8}$Be
which is naturally assumed to be unity, and that of the cluster configuration
$\alpha$+$^{8}$Be in $^{12}$C. The latter is determined as $S_\alpha=|A_{NL}|^2$,
where the spectroscopic amplitude $A_{NL}$ is given by the dinuclear overlap
\begin{equation}
\langle{\rm ^{8}Be}|^{12}{\rm C}\rangle=A_{NL}\Phi_{NL}(\bm{r}). \label{eq14}
\end{equation}
Because two different exit channels of the $\alpha$ transfer $^{12}$C($\alpha$,$^8$Be)
reaction were identified, with the emitting $^8$Be being in the g.s. and excited
$2^+$ state \cite{Woz76}, one needs to evaluate (\ref{eq14}) for the two
configurations $\alpha$+$^{8}$Be$_{\rm g.s.}$ and $\alpha$+$^{8}$Be$^*_{2^+}$,
which are associated with $\Phi_{N=3,L=0}(\bm{r})$ ($S$-wave) and
$\Phi_{N=2,L=2}(\bm{r})$ ($D$-wave). In general, one can treat these two $S_\alpha$
values as free parameters to be adjusted by the best DWBA or CRC fit to the $\alpha$
transfer data. Instead of this procedure, we have adopted in the present work the
$S_\alpha$ values predicted for these configurations by Kurokawa and Kato using
the CSM method \cite{KaKu07}. Namely, $S_\alpha ({\rm g.s.})\approx 0.36$ and
$S_\alpha (2^+)\approx 0.38$, which are rather close to the spectroscopic factors
predicted recently by other cluster models \cite{FHO09,Enyo07}. Note that the
$S_\alpha$ values used in our CRC calculation are also close to those extracted
from the DWBA analysis of the $^8$Be transfer reaction
$^{24}$Mg($\alpha$,$^{12}$C)$^{16}$O \cite{BHJ80}. The same $^{8}$Be+$^{8}$Be\  OP has been used
for both  $^8\text{Be}+^8$Be$_{\rm g.s.}$ and $^8\text{Be}+^8$Be$^*_{2^+}$ exit
channels of the final partition (see more discussion below).

\begin{figure}[bht]\vspace*{-1cm}
 \includegraphics[width=\textwidth]{fig6.eps}\vspace*{-0.5cm}
 \caption{CRC description of the $\alpha$ transfer reaction
$^{12}$C($\alpha$,$^{8}$Be)$^{8}$Be$_{\rm g.s.}$ measured at $E_{\alpha}=65$ MeV
\cite{Woz76}, using the $\alpha$ spectroscopic factor $S_\alpha({\rm g.s.})\approx 0.36$
taken from the CSM calculation \cite{KaKu07}. The CRC results obtained with the CDCC-based
optical potentials $U_{\rm WSD}$ and $U_{\rm LEP}$ for the $^{8}$Be+$^{8}$Be\  partition are shown
as the solid and dash-dotted lines, respectively. The dashed line is the CRC result
obtained with the WS potential $U_{\rm Int WS}$, interpolated from the OP's adopted
for the $^{7,8}$Be+$^9$Be systems at the nearby energies \cite{Rom09}.} \label{f6}
\end{figure}

The CRC results for the $\alpha$ transfer $^{12}$C($^4$He,$^8$Be) reaction at
$E_{\alpha}=65$ MeV to the ground state of $^8$Be are compared with the measured data
\cite{Woz76} in Fig.~\ref{f6}. One can see that the CDCC-based optical potentials
of the $^{8}$Be+$^{8}$Be\  partition ($U_{\rm LEP}$ and $U_{\rm WSD}$ determined by the methods
(i) and (ii), respectively) give a good CRC description of the $\alpha$ transfer data
without any adjustment of its strength, using
$S_\alpha({\rm g.s.})=|A_{30}|^2\approx 0.36$ taken from the results of the CSM
calculation \cite{KaKu07}. With a better OM description of the CDCC elastic cross
section given by the $U_{\rm WSD}$ potential (see Fig.~\ref{f3}), the CRC cross section
given by $U_{\rm WSD}$ also agrees slightly better the measured $\alpha$ transfer data.

Because the $^{8}$Be+$^{8}$Be\  OP is unknown so far, a practical assumption is to estimate it
from the phenomenological OP's adopted for the neighboring $^{9,7}$Be isotopes.
For example, the proton transfer reaction $^7$Li($^{10}$B,$^9$Be)$^8$Be was measured
by Romanyshyn {\it et al.} \cite{Rom09}, and a deep WS potential was deduced for
the real OP of the $^8$Be+$^9$Be system at $E_{\rm c.m.}=31.7$ MeV from the DWBA
analysis of transfer data. Interestingly, these authors also found that the
$^8$Be+$^9$Be OP is quite close to that adopted earlier for the $^7$Be+$^9$Be
system (see Fig.~10 of Ref.~\cite{Rom09}). Therefore, one might expect the $^{8}$Be+$^{8}$Be\
OP to be close to the WS optical potentials adopted for the $^{7,8}$Be+$^9$Be systems.
To explore the reliability of this practical approach, we have interpolated
the $^{8}$Be+$^{8}$Be\  OP from those of the $^{7,8}$Be+$^9$Be systems adopted in Ref.~\cite{Rom09}
and denoted it as $U_{\rm Int WS}$, with $V_v=155.0$ MeV, $R_v=3.152$ fm, $a_v=0.768$ fm;
and $W_v=13.5$ MeV, $R_w=5.6$ fm, $a_w=0.768$ fm. The use of $U_{\rm Int WS}$ in the CRC
calculation of the $^{12}$C($\alpha$,$^{8}$Be) reaction completely fails to account
for the data (see dashed lines in Fig.~\ref{f6}). In fact, the spectroscopic factor
$S_\alpha({\rm g.s.})$ taken from Ref.~\cite{KaKu07} must be scaled by a factor
of 25, so that the CRC cross section obtained with $U_{\rm Int WS}$ can be comparable
with the measured $\alpha$ transfer data.
\begin{figure}[bht]\vspace*{-0.5cm}
 \includegraphics[width=\textwidth]{fig7.eps}\vspace*{-0.5cm}
 \caption{The same as Fig.~\ref{f6} but for the $\alpha$ transfer reaction with
 one emitting $^{8}$Be nucleus being in its $2^+$ state ($E_x\approx 2.94$ MeV),
 using the $\alpha$ spectroscopic factor $S_\alpha(2^+)\approx 0.38$ taken from
 the CSM calculation \cite{KaKu07}} \label{f7}
\end{figure}

We have also performed the CRC calculation of the $^{12}$C($\alpha,^8$Be)$^8$Be$^*$
reaction at $E_{\alpha}=65$ MeV with one emitting $^8$Be nucleus being in its $2^+$
state ($E_x\approx 2.94$ MeV). Although the $2^+$ state of $^8$Be is a broad resonance,
its 2$\alpha$-cluster structure remains similar to that of the ground state, and
the $\alpha$ transfer cross section measured for the $2^+$ state is of about the same
strength as that measured for the g.s. as shown in Figs.~\ref{f6} and \ref{f7}.
The  $\alpha$ spectroscopic factors $S_\alpha$ predicted by the CSM calculation
\cite{KaKu07} are also close for both the ground- and $2^+$ states.  It is,
therefore, reasonable to use the same CDCC-based $^{8}$Be+$^{8}$Be\  OP for the partition
$^8$Be+$^8$Be$_{2^+}^*$ in the exit channel. The results of the CRC calculation
are compared with the data \cite{Woz76} in Fig.~\ref{f7}, and one can see that
the (unrenormalized) CDCC-based $^{8}$Be+$^{8}$Be\  OP also delivers a good description
of the $\alpha$ transfer data using $S_\alpha(2^+)=|A_{22}|^2\approx 0.38$ given
by the CSM calculation \cite{KaKu07}. The use of $U_{\rm Int WS}$ for the $^{8}$Be+$^{8}$Be\  OP
in the CRC calculation also strongly underestimates the $\alpha$ transfer data
(see dashed line in Fig.~\ref{f7}).

In conclusion, a good CRC description of the measured $^{12}$C($\alpha,^8$Be) data
has been obtained with the $^{8}$Be+$^{8}$Be\  OP's determined by the methods (i) and (ii) from
the elastic $^{8}$Be+$^{8}$Be\  cross section given by the 4$\alpha$ CDCC calculation, and
with the $\alpha$ spectroscopic factors given by the CSM calculation \cite{KaKu07}.
The fact that no adjustment of the potential strength of $U_{\rm LEP}$ and
$U_{\rm WSD}$ was necessary suggests that the 4$\alpha$ CDCC method is a reliable
approach to study the $^{8}$Be+$^{8}$Be\  system. These results also show that, despite the short
life-time of $^{8}$Be, the $\alpha$ transfer $^{12}$C($^4$He,$^8$Be) reaction is very
sensitive to the OP of the $^{8}$Be+$^{8}$Be\  partition. The measured $\alpha$ transfer data
clearly prefer the shallow OP based on the result of the 4$\alpha$ CDCC calculation
using the AB potential of the $\alpha+\alpha $\  interaction \cite{AB66}, over the deep WS potential
interpolated from those adopted for the $^{7,8}$Be+$^9$Be systems \cite{Rom09}.

\section{Summary}
The 3-body and recently suggested 4-body CDCC methods \cite{De18} have been used
to predict the elastic $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  scattering at $E_{\rm c.m.}=43.3$ and 41.3
MeV, respectively, using the $\alpha+\alpha $\  interaction suggested by Ali and Bodmer \cite{AB66}
that well reproduces the experimental $\alpha+\alpha $\  phase shifts. The elastic cross sections
predicted by the CDCC calculation were used to determine the local OP’s of the $\alpha+^{8}$Be\
and $^{8}$Be+$^{8}$Be\  systems.

The CDCC-based $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  OP's are further used as the inputs of the core-core OP
and that of the final partition, respectively, in the CRC study of the $\alpha$
transfer $^{12}$C($\alpha,^8$Be) reaction measured at $E_\alpha=65$ MeV \cite{Woz76},
with the emitting $^8$Be being in both the g.s. and $2^+$ state. These $\alpha$ transfer
data are well reproduced by the CRC results obtained with the CDCC-based OP’s and
$\alpha$ spectroscopic factors of the $^8$Be+$^8$Be$_{\rm g.s.}$ and
$^8$Be+$^8$Be$_{2^+}^*$ configurations in $^{12}$C taken from the CSM cluster
calculation \cite{KaKu07}.

As alternative to the shallow (surface-type) $^{8}$Be+$^{8}$Be\  OP determined from the
elastic $^{8}$Be+$^{8}$Be\  cross section predicted by the 4$\alpha$ CDCC calculation, a deep
WS potential with parameters interpolated from the OP’s adopted for the
$^{7,8}$Be+$^9$Be systems at the nearby energies \cite{Rom09} has been used
in the CRC calculation, and it grossly underestimates the $\alpha$ transfer data.
This might be due to the fact that $^7$Be and $^9$Be are well bound nuclei, and
the breakup effect is therefore much weaker than that of the unbound $^8$Be.

We conclude that the $\alpha+^{8}$Be\  and $^{8}$Be+$^{8}$Be\  optical potentials can be determined from the
elastic scattering cross section predicted, respectively, by the 3$\alpha$ and
4$\alpha$ CDCC calculations using the realistic $\alpha+\alpha $\  interaction that properly
reproduces the experimental $\alpha+\alpha $\  phase shifts \cite{AB66}. The present CRC study should
motivate further theoretical and experimental studies of the $\alpha$ transfer
$^{12}$C($\alpha,^8$Be) reaction as a probe of the 4$\alpha$ interaction and
the $\alpha$-cluster structure of $^{12}$C.

\section*{Acknowledgement}
The present research has been supported, in part, by the National Foundation for
Scientific and Technological Development of Vietnam (NAFOSTED Project No. 103.04-2016.35).
P.D. is Directeur de Recherches of F.R.S.-FNRS, Belgium, supported by the Fonds
de la Recherche Scientifique - FNRS under Grant Number 4.45.10.08. Computational resources have
been provided by the Consortium des Équipements de Calcul
Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS)
under Grant No. 2.5020.11 and by the Walloon Region. We also thank
Nguyen T.T. Phuc for his helpful discussion.

\bibliographystyle{apsrev4-1}