theorem]Definition

Sample image frames of scenes generated by the tool during runs configured
with four different parameters sets.

# Procedural Generation and Rendering of Realistic, Navigable
Forest Environments: An Open-Source Tool

Callum Newlands0000-0003-4798-9306 and Klaus-Peter Zauner
Electronics and Computer Science, University of Southampton, UK {cn2g18,
<EMAIL_ADDRESS>
###### Abstract

Simulation of forest environments has applications from entertainment and art
creation to commercial and scientific modelling. Due to the unique features
and lighting in forests, a forest-specific simulator is desirable, however
many current forest simulators are proprietary or highly tailored to a
particular application. Here we review several areas of procedural generation
and rendering specific to forest generation, and utilise this to create a
generalised, open-source tool for generating and rendering interactive,
realistic forest scenes. The system uses specialised L-systems to generate
trees which are distributed using an ecosystem simulation algorithm. The
resulting scene is rendered using a deferred rendering pipeline, a Blinn-Phong
lighting model with real-time leaf transparency and post-processing lighting
effects. The result is a system that achieves a balance between high natural
realism and visual appeal, suitable for tasks including training computer
vision algorithms for autonomous robots and visual media generation. The
application and supporting configuration files can be found at
https://github.com/callumnewlands/ForestGenerator. <ccs2012> <concept>
<concept_id>10010147.10010371</concept_id> <concept_desc>Computing
methodologies Computer graphics</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10010147.10010371.10010372</concept_id> <concept_desc>Computing
methodologies Rendering</concept_desc>
<concept_significance>100</concept_significance> </concept> <concept>
<concept_id>10010147.10010371.10010396.10010397</concept_id>
<concept_desc>Computing methodologies Mesh models</concept_desc>
<concept_significance>300</concept_significance> </concept> </ccs2012>
[500]Computing methodologies Computer graphics [100]Computing methodologies
Rendering [300]Computing methodologies Mesh models

##  1 Introduction

Procedural generation and efficient rendering of vegetation and natural
environments is a growing research area in both the entertainment and
commercial industries. Whilst the former is often concerned with visual
attractiveness and generation speed [Mor16, Int20, Lan20], the latter often
focuses on real-world accuracy and generation flexibility, for purposes such
as environmental forecast modelling or simulating real-world locations
[HSAKD98, Kom20, Len20, Vir20].

We here present an open-source Forest Generator (available with sample
configurations at [New22]) for procedurally generating and rendering realistic
forests which can be virtually navigated in real-time, or through asynchronous
connections with external software. The system takes a middle-ground
perspective between visual and real-world accuracy, resulting in an approach
that is ideal, for example, for training the computer-vision systems of forest
navigating robots [HRS13, SNK∗20, TGZ20].

The main contributions of this article are two-fold, firstly, a brief review
of the existing approaches and literature in several aspects of forest
generation and rendering (Section 2)—including vegetation generation,
rendering optimisation, ecosystem simulations and tree-specific
illumination—and secondly, the presentation of a system for generating and
rendering realistic forest environments. Section 3 gives an overview of the
system architecture and control flow—which uses L-systems for representing and
generating a wide range of tree types (Section 4), along with a tree-specific
ecosystem simulation for realistic scene distributions (Section 4.3). The
generated scene is rendered using a Blinn-Phong deferred render pipeline,
along with a leaf-specific translucency model, and post-processing shaders for
shadows and lighting effects (Section 5).

##  2 Related Work

Procedural generation and rendering of virtual forest environments comprises
many aspects; key among them are plant model generation, vegetation
distribution, terrain generation, rendering optimisation, and tree-specific
illumination. The following provides a brief overview of the state-of-the-art
in these areas.

###  2.1 Vegetation

Methods for vegetation generation can be categorised as procedural and non-
procedural. Non-procedural methods allow for high levels of control and can
result in highly detailed models, at the cost of requiring heavy user
input—often at a per-plant level. Key non-procedural methods include manual
modelling, and image-based methods such as photogrammetry [BNB13, Pat18].
Procedural models are constructed algorithmically and controlled through
parameter values without requiring high levels of manual input or expertise.
Several hybrid methods have also been developed to reduce the time taken to
create plants under non-procedural methods, by combining algorithmic
generation with user descriptions (component graphs) [LD99] or sketch-guided
modelling [BvMM11].

Early procedural models for trees include the H-model and P-model [HTF82,
PL90] which defined branch geometry parametrically in terms of their parent
branches through angles and lengths, and further parametric models were
subsequently proposed [WP95, DGCF01, Hew17]. Several methods have been
introduced which take into account the fractal-like (recursive) nature of tree
generation, and operate in terms of repeated sub-modules including Iterated
Function Systems [Bar88] and self-similar branch modules [MHS∗19].

The most well-explored plant modelling technique based on the idea of
recursive hierarchies is the concept of L-systems. L-systems are a parallel
string-rewriting system and were introduced by Lindenmayer in order to model
the growth of cellular organisms [Lin68]. There are many types and extensions
of L-systems [Pru86, PL90, PHHM96] including (deterministic) context-free,
context-sensitive, stochastic, parametric, differential [PHM93],
environmentally-sensitive [PJM94] and open [MP96]. However, often context-
sensitive L-systems, stochastic L-systems, and parametric L-systems are
combined, as is the case here.

Prusinkiewicz was the first to propose a graphical interpretation of L-systems
using a ‘turtle’ [Pru86]. After the desired number of iterations of an
L-system, the resulting string can be interpreted by the turtle by reading
left to right and acting on each symbol in a predefined manner. The turtle
interpretation, combined with the expressive power of L-systems, allows for
the generation of many complex graphics from the specification of simple
rules, and was later extended into 3D and combined with parametric L-systems
to allow for rotation and movement by arbitrary amounts [PL90]. A further key
extension in the field of plant representations was the inclusion of brackets
([ and ]) to allow for branching structures such as trees to be represented
[Lin68] along with their turtle interpretation through the use of a stack
[Pru86].

Another category of procedural plant generation methods are particle-based
methods, which simulate the motion of particles to construct the branch
structure of trees. These include Diffuse-Limited Aggregation [WS81], Space
Colonisation [RLP07] and particle flow modelling [NFD07, NM19].

As well as the generation of plant models, the distribution of plants within a
scene is critical to realism and visual appeal. The simplest approach for
distributing vegetation is to sample a random probability distribution, or use
an algorithm such as Poisson disk sampling [Coo86]. This naïve approach
results in unrealistic plant distributions and specialised methods have been
proposed. These can be grouped according to whether plant positions are
generated from global distribution properties (global-to-local) or whether the
overall distribution arises from simulating interactions between individual
plants (local-to-global) [LP02].

Several global-to-local methods have been proposed, including: modelling
ecotypes from terrain data [Ham01], probability-deformation-kernels [LP02],
halftoning algorithms [DHL∗98], and sample pattern extrapolation [XYY15].
Alsweis and Deussen presented a method for using Wang Tiles, along with the
‘field of neighbourhood’ (FON) distribution model [BH00] to model the
influence of plants on the growth of other plants in their surrounding area
[AD06]. This idea was developed further using Poisson disk distributions
(PDDs) and configurable ecosystem properties [WHDS13], and to accommodate
plants of different sizes without size-related bias [OBRvdK17]. Recently,
Nascimento et al.presented a method which takes advantage of GPU parallelism,
and uses the idea of ecosystems from [Ham01] along with layering, to allow for
different sized plants with per-layer parameter maps and cross-layer growth
interactions [dNFP18].

Local-to-global methods typically allow for higher plant-level interactions,
at the cost of higher compute times, or possible inconsistencies at a global
scale. Early models were based on open L-systems and distributed plants by
simulating resource competition or enforcing species-terrain preferences
[DHL∗98]. This was extended to utilise multiset L-systems and simulate
succession and plant-clustering (using deformation kernels) [LP02]. Beneš
proposed another local-to-global algorithm in the form of a survival
simulation on randomly placed plants [Ben02]. The simulation models the growth
and death of plants, seed dispersal and clustering, resource competition and
removal of colliding plants. Mikuličic and Mihajlović compared random
scattering (with collision detection) to a similar survival simulation and
found that the simulation produced a more natural result [MM16].

###  2.2 Terrain

Terrain generation methods are based on one of two representations:
heightmaps—which store the terrain as a 2D image representing terrain
height—and voxel terrain—where voxels are elements of a 3D volume grid
(typically cubes). The use of heightmaps restricts the terrain to be one
continuous surface with no caves or overhangs, however this can be overcome
with layered heightmaps [BF01] or heightmap transformations [GM01]. Voxel
terrain does not have this limitation but generally has a much higher memory
cost and therefore is typically only used in situations where the terrain
needs to be edited interactively. Here, only procedural, heightmap-based
methods will be discussed, however non-procedural methods include manual
modeling, using real-world data, and sketch-based methods [GMS09, HGA∗10].

Terrain height functions map a coordinate pair to a height value to create
heightmaps before or during render time. One of the most common categories are
noise functions—due to the fact that they can be evaluated without relying on
neighbouring points—and include random noise, coherent noise, fractal noise
and ridged noise [Ols04, Häg06, Dun16, FE17].

Aside from terrain-height functions, several iterative generation methods have
been proposed including midpoint displacement, diamond-square subdivision and
faulting algorithms [BA05, Häg06, Kah13]. At a larger terrain-scale, many
physics-based erosion models have also been proposed [MKM89, CMF98, BF01,
BF02, BTHB06, MDH07, KBKv09, JT11, Kur12, Bey15, BEP18], however, due to the
comparatively small scale of forests, these methods are somewhat unsuitable
for the problem discussed here and are included here only for completeness.

Once the terrain has been generated, post-generation techniques can be applied
to increase realism or reduce render time. Three of the most common mesh-
optimisation methods are: ROAM [DWS∗97], the Lindstrom-Coller Simplification
[LKR∗96] and Geometrical MipMapping [dB00]. Terrain smoothing can be performed
by applying common imaging filters to the heightmaps [SKG∗09], or using agents
which perform random walks and average the height of points [DP10].

###  2.3 Rendering Optimisation

Tree models can have in the order of 100,000 polygons [FUM05] (see
supplementary material), therefore rendering even a few tree models in full
detail can quickly become infeasible if the scene is to be navigable with
reasonable frame rates. One of the key concepts in rendering complex scenes is
Level of Detail (LOD) which involves presenting models using different methods
based on their distance from the viewpoint.

Plant model LOD rendering techniques can be broadly classified into three
categories: geometric simplification, point-based, and image-based. Geometric
simplification methods reduce the number of polygons in the model and include
general polygonal simplification methods (reviewed in [HG97]), as well as
plant-specific models such as FSA [RCB∗02], BFSA and TFSA [DZYJ10]. Point-
based methods partially or totally replace models with point-cloud
representations, as points are typically faster to render than polygon
vertices [RB85, WP95, SD01, WFP∗01, DCSD02]. Image-based methods replace 3D
model meshes with 2D representations—which have much smaller polygon
counts—and result in the highest vertex compression rate at the cost of
reduced detail.

The most common image-based method is a billboard, which replaces a model with
a flat texture representation mapped on a rectangular plane. There are several
types of billboard, which reflect the trade-off between complexity reduction
and accuracy to the original model [DDSD03, DN04, BCF∗05, DCDS05, FUM05,
Häg06, DZYJ10, BLZ∗11, Jak00, BZCJ09], with the most notable types being
camera-aligned billboards (one billboard per model), crossed billboards (2–3
billboards per model) and billboard clouds (a set of $n$ billboards). Among
the limitations of 2D billboard representations are the artefacts when viewed
at shallow angles, and the lack of parallax and depth. As a result, several
other techniques have been presented which consider model depth. The most
notable methods, with regards to walk-through forests, are “Sprites with
Depth” [Sch97, SGHS98], and Layered Depth Images (LDIs) [SGHS98, MO95, Max96,
MDK99]; with other methods including volumetric textures [DN04], volumetric
billboards [DN09], “2.5D imposters” [SK03] and Hierarchies of Bidirectional
Texture Functions (HBTs) [MNP01].

There are several techniques that are used alongside LOD in order to optimise
the rendering further or remove visual artefacts (such as “popping” between
LOD levels); including quadtree data structures [BN08, BN12, OBRvdK17, KS18],
instancing [DHL∗98, KGM14, dV20], GPU mesh construction [BN12, MM16,
OBRvdK17], importance reduction [DCSD02], cross-fading (or alpha-blending)
[Jak00, Häg06], alpha-to-coverage [JRSK09, KCS07], and fizzle LOD [Wha05,
KCS07].

###  2.4 Tree-Specific Illumination Models

Forest environments exhibit unique lighting behaviours, in particular, self-
shadowing between leaves, and complex indirect lighting due to leaf
translucency and light scattering. Accordingly, tree-specific illumination
models have been suggested. Qin et al.presented one of the first models,
limited to a fixed position and offline (non-real-time) rendering [QNTN03].
Both Behrendt et al.and Candussi et al.introduced real-time methods for
rendering trees represented as dynamically rendered billboard clouds, with
differing lighting models to support either accurate external lighting of
large forests [BCF∗05] or realistic, scalable, real-time animation [CCH05].

Wang et al.presented a real-time lighting model which accommodated subsurface
scattering and transmittance of leaves resulting in translucency [WWD∗05]. The
model operates on a per-leaf basis and defines spatially-varying bidirectional
reflectance distribution functions (BRDFs) and bidirectional transmittance
distribution functions (BTDFs) as parameter maps obtained from real leaves.
Habel et al.presented a similar model in [HKW07], but using an 8D
bidirectional scattering-surface reflectance distribution function (BSSRDF)
for the back face of the leaves (as opposed to a 4D BTDF). Rather than
considering individual leaves, as the previous two models do, Boulanger et
al.introduced a lighting model which probabilistically models the spacial
distribution of leaves [BBP08]. The model results in lighting close to Monte-
Carlo PBR rendering, however is limited to a uniform distribution of leaf
positions and normals, so has limited accuracy for some tree types. As ray-
traced methods become feasible in real-time on consumer hardware, the focus is
shifting towards them [GS08, BN12].

##  3 System Design

The Forest Generator constructs and renders interactive visualisations of
forest scenes in four stages (Figure 1).

Figure 1: The four-stage system flow for generating a forest scene.

The generation and rendering are fully configurable through text files (see
[New22]). These configuration files—in human-readable YAML format—facilitate a
wide range of different forest characteristics to be generated, including the
creation of new tree models with differing geometry and textures. The seed for
the pseudo-random generator can also be specified to allow for identical runs
of the application.

Based on the configuration, the application generates the object models and
distributes them within the virtual scene. Next, any pre-rendering operations
are carried out—such as shadow map rendering—and finally, the system enters
the render loop for the final navigable scene.

Movement around the scene at render time can be controlled in two ways:
manually and asynchronously. With manual control, the user can pan the view in
real-time using the mouse and move around the scene using the keyboard. In
asynchronous mode, the application blocks after each frame and waits for a new
line of input on the standard input stream, to parse and react to, before
rendering again. As well as real-time output of the scene, the rendered frames
can be saved as image files. This, combined with the asynchronous input mode,
allows for external software (such as computer vision systems, frame-to-video
converters, or image-to-image style transformers) to be easily connected to
the input and output of the system, as well as allowing the system to be run
in a headless environment.

###  3.1 Scene Construction and Storage

The scene generation and pre-rendering stages are performed offline, which
limits the size of the forest to being finite and has the drawback of large
memory costs—as the entire scene needs to be stored in memory—however, due to
the fact that real-world forests are also finite, and the increased realism
and quality facilitated by this method, this is deemed to be a worthwhile
trade-off.

To store the generated scene, a quadtree data structure is used. In this case,
the tree is complete, as it is fully constructed during the generation phase,
with different depth nodes selected to be visible at render time. This enables
various parts of the scene to be rendered at different levels of detail (LOD)
(Section 5.3) without having to calculate the distance to each object, as well
as enabling efficient view-frustum-culling.

The scene is constructed in the following stages:

  1. 1.

The quadtree data structure is constructed recursively to a max depth $D$.

  2. 2.

At each leaf node, a terrain tile is constructed (Section 4.2) and ground
cover for the tile is generated and distributed (Section 4.2).

  3. 3.

Minimum and maximum tree masks are calculated and trees are distributed with
an ecosystem simulation (Section 4.3).

  4. 4.

Models for trees of the required growth stages are constructed (Section 4.1)
and inserted into the relevant leaf nodes.

As all scene models are constructed before render time, this can result in a
large memory cost. To minimise this, tree meshes are stored in an object pool
and instances are replicated throughout the scene with differing
transformations. The number of variants for each growth stage (iteration
count) of a tree species can be controlled through a system parameter, to
reduce repetition as desired. Ground cover objects, due to their high per-tile
density, are rendered using terrain-tile level GPU instancing. This reduces
the overhead of sending geometry to the GPU, whilst still allowing tiles to be
culled when out of view (Section 5.3).

##  4 Generation

Prior to rendering, the system generates the mesh data for each of the
elements in the forest scene: trees, terrain and ground cover (plants, twigs,
rocks, etc.), and then carries out an ecosystem simulation to distribute the
trees within the scene.

###  4.1 Trees

To generate tree models, parametric, stochastic, context-sensitive L-Systems
(Section 2.1) were implemented along with a 3D turtle interpreter. The
implemented turtle holds its position, heading, and up vectors
$\in\mathbb{R}^{3}$; current radius (to allow for branch tapering); and the
vertices of the most recent cross-section.

With each interpreted symbol—see supplementary material for the full list of
symbols and modules (parameterised symbols)—the turtle’s state is transformed
using standard 3D matrix transformations, and lists of mesh vertices and sub-
model references (position, heading, up-vector and sub-model list index) are
added to where necessary by the interpreter. If tropism is enabled—by seeing a
module of the form $T(v_{1},v_{2},v_{3},e)$—each time the cross-section is
transformed, a (small) rotation (with rotation factor $\alpha$) is applied in
the direction $\bm{v}=[v_{1},v_{2},v_{3}]^{T}$ in the same manner as [PL90] to
simulate the growth of plants dependent on the direction of external forces
(such as gravity or sunlight):

$\alpha=\frac{\theta_{e}}{\theta_{v}}=\frac{e\,\lvert\lvert\bm{h}\times\bm{v}\rvert\rvert}{\arccos\left(\frac{\bm{h}\cdot\bm{v}}{\lvert\lvert\bm{h}\rvert\rvert\,\lvert\lvert\bm{v}\rvert\rvert}\right)}$
(1)

where $\theta_{e}$ is the angle of rotation for the tropism, proportional to
the elasticity ($e$) of the branch, and $\theta_{v}$ is the angle between the
turtle heading $\bm{h}$ and tropism direction $\bm{v}$.

To construct the branch geometry, adjacent cross-sections of $N$ edges are
considered pairwise. For each adjacent pair of vertices in the two cross-
sections, a connecting face $ABCD$ is constructed (see Figure 2), where the
normal vector for a vertex is determined using the cross product of its two
adjacent edges; (e.g. $\widehat{\bm{n_{B}}}=\vv{BD}\times\vv{BA}$). Since the
faces are flat, the tangent vector for each vertex is simply the vector along
one of its adjacent edges. After all connecting faces have been created, the
four normal vectors at each vertex are averaged to create smooth normals
across edges, and similarly for the tangents. This prevents flat shading and
allows for a lower value of $N$ to be used with the same visual appearance.
Each quadrilateral face is divided into two primitive triangles $ABD$ and
$ADC$ for rendering. To save memory, the faces of the interior cross-sections
are not generated.

BADC Figure 2: Diagram of the geometry generated by a turtle interpreting the
commands $F(1)!(0.6)F(1)F(1)$, with $N$ (number of edges) $=6$. Vertices
$A,B,C,D$ outline a connecting face between 2 cross-sections.

To inject a sub-model (such as a leaf or pine cone) into the mesh, given a
sub-model reference (position $\bm{p}$, heading $\bm{h}$, up-vector $\bm{u}$
and model index $i$), the geometry is first rotated so that
$\bm{x}=[1,0,0]^{T}$ aligns with $\bm{h}$, and then rotated around
$\bm{y}=[0,1,0]^{T}$ with angle
$\theta=\operatorname{arctan2}\left((\bm{y}\times\bm{u})\cdot\widehat{\bm{x}},\,\bm{y}\cdot\bm{u}\right)$,
so that the $y$-axis of the sub-model aligns with $\bm{u}$. Next, the geometry
is translated so that the origin of the sub-model lies at $\bm{p}$, and the
resulting vertex data is added to the combined mesh for all injections of
model $i$.

The module $F(d,n_{l},m_{i},\alpha_{r},\alpha_{l})$ is used to inject leaves
which grow out of the branches (such as pine needles or poplar leaves) as sub-
models. This causes the turtle to distribute $n_{l}$ copies of sub-model
$m_{i}$ along a spiral around the constructed prism, controlled by the radial
angle $\alpha_{r}$ and lift angle $\alpha_{l}$. To introduce some pseudo-
random variation into the model placement, the radial angle ($\alpha_{r}$) and
distance between consecutive sub-models ($s$) are not strictly constant, but
instead are offset with a value sampled from a uniform distribution with range
$\pm$ 50% of their (constant) values; resulting in
$\alpha_{r}\in\left[-\frac{\alpha_{r}}{2},\frac{\alpha_{r}}{2}\right)$ and
$s\in\left[-\frac{s}{2},\frac{s}{2}\right)$.

The low LOD representations of the branches are constructed using the same
L-system and parameter values, but using a cross-section with 2 edges to
create a flat version of the tree with the same branches. For the leaves,
adjacent pairs of sub-model references are merged with an average position and
angle, and replaced with a model of double width and 1.1 times height in a
fashion similar to the ideas presented in [RCB∗02].

A branching tree—in this work—is one where the main trunk (and subsequent
branches) split into $n$ branches, and the original trunk ceases. As well as
this, there is the possibility for ‘side branches’ to protrude from the trunk
or branch before it splits (Figure 3).

Figure 3: Examples of a branching tree (left) and monopodial tree (right)
generated by the system.

In the L-systems used, a branching tree species is specified with a number of
parameters—such as initial trunk width or branch offsets—and a set of $n$
(possible) branchings $Br$. The $i$th branching ($Br_{i}$) is a branching into
$m_{i}+1$ sub branches with angles
$\theta_{i}^{0},\,\theta_{i}^{1},\,\dots,\,\theta_{i}^{m_{i}}$ and probability
$p_{i}$ such that $\sum_{i}p_{i}=1$. $\theta^{k}$ is the angle around the
trunk between branches $k$ and $k+1$.

$Br=\left\\{\,\left\langle\theta_{0}^{0},\,\dots,\,\theta_{0}^{m_{0}},\,p_{0}\right\rangle,\,\dots,\,\left\langle\theta_{n}^{0},\,\dots,\,\theta_{n}^{m_{n}},\,p_{0}\right\rangle\,\right\\}$
(2)

Pseudo-random variation is achieved by providing upper and lower bounds
($x_{UB}$ and $x_{LB}$ respectively) for the values of each parameter
(L-system constant), and then, at each usage, a uniformly distributed random
number is sampled from the distribution
$x\sim\mathcal{U}\left[x_{LB},x_{UB}\right)$.

A monopodial tree is one with a single central trunk which extends to the top
of the tree and has side branches protruding outwards from it (Figure 3). Two
L-systems are used to accommodate two different types of ‘side branching’:
alternating and pine-style. Alternating branches appear in a spiral around the
main trunk, with each side branch having a number of 3-level side branches.
Pine-style branches do not appear in a spiral, but rather all appear around
the trunk in discrete layers. The side branches themselves are also different
in that, at each iteration, they split from a single point into two higher-
level side branches.

###  4.2 Terrain and Ground Cover

The scene terrain is constructed from a heightmap which is generated using 2D
fractal simplex noise [Per02] and is configurable through the system
parameters. Details of the terrain mesh construction can be found in the
supplementary material for this paper.

The system has four categories of ground cover: fallen leaves, twigs, crossed
billboards (plants) and external models—with varying levels of customisation.
Fallen leaves are generated by positioning leaf sub-models on (just above) the
ground mesh with controllable density and scale, and twigs are similarly
configurable: generated using a simplified form of the alternating branch
monopodial tree L-system (see supplementary material for the specific
L-system). Their geometry is formed in the same manner as the trees, and then
rotated $\frac{\pi}{2}$ radians around the $x$-axis so the main axis of the
twig lies parallel to the floor.

Crossed billboards (Section 2.3) are rendered as $N$ intersecting planes
(boards). By default, the system generates grass and two fern types as crossed
billboard objects. Each of the boards is scaled in height and width, and
rotated by a random angle
$\phi\in\left[\frac{4\pi}{24},\,\frac{10\pi}{24}\right)$ around a random
vector $\bm{r}$ which lies on the $xz$-plane. This orientates the board at a
random angle whilst keeping the base of the board at the origin. The board is
then rotated by $\theta_{i}$ radians around the y axis where

$\theta_{i}=\frac{i\pi}{N},\quad i\in\left\\{\,0,\,\dots,\,N-1\,\right\\}$ (3)

to distribute the boards around a circle with random angles (Figure 4).

Figure 4: Example of the board rotations for high LOD crossed billboard
objects with four boards (left) and their appearance when rendered (right).

The low LOD representations are constructed using two perpendicular boards.

As well as the three built-in ground cover object types, the system also
allows for any external 3D models to be imported and distributed on the
ground. For example more-detailed leaf piles or fallen tree logs could be
imported from files in common 3D file types. The models (and textures) are
loaded using the Open Asset Import (Assimp) Library [LWJ20], and converted
into the local system format at load time.

As mentioned in Section 3.1, ground cover objects are drawn with OpenGL
instanced rendering. Each instance is given a unique scale, position and
rotation in the form of a model matrix. The density parameter of each object
type ($\rho$) controls the total number of instances ($n_{I}$) to create:

$n_{I}=\frac{d_{0}\cdot\rho}{n_{T}}w^{2}$ (4)

where $w$ is the width of the scene, $d_{0}$ is the default density for the
object and $n_{T}$ is the number of ‘species’ of that object type. The
division by $n_{T}$ ensures that the object type remains at a constant density
regardless of the value of $n_{T}$. For fallen leaves and twigs, $n_{T}=1$.

###  4.3 Ecosystem Simulation

To distribute the trees throughout the scene, an ecosystem simulation was
implemented. This creates a more realistic distribution than random placement,
as natural behaviours emerge such as clustering of species and areas of
negative growth around trees. The simulation is similar to the ones presented
in [Ben02] and [MM16], but has been adapted specifically for trees. Rather
than considering plants as circles on a plane, here, the trees are considered
as 3D masks with a trunk cylinder and a canopy cylinder, which allows for
shorter trees to exist underneath the canopy of larger ones, and enables the
distinction between branches intersecting tree trunks and intersecting other
branches.

To initialise the simulation, a number of plants for each specified tree
species (determined by Equation 4) are created with random ages
$a_{i}\in\left[0,\,a_{M}\right)$—where $a_{M}$ is the maximum age for the
species—and are randomly placed in the scene largest to smallest.

After initialisation, the simulation runs for $N$ steps, with $N_{Y}$ steps
constituting a year. For each step, the following occur:

  1. 1.

If it is the end of the year, all trees seed a number of new plants in a ring
around themselves.

  2. 2.

For each pair of colliding trees, the plant with lower viability is removed.

  3. 3.

Plants which are older than their maximum age are considered dead and are
removed.

  4. 4.

Each plant grows (its age is increased by 1).

The viability ($v$) of each plant is determined similarly to [Ben02] with an
age component and a negative feedback component to prevent exponential growth
of a plant species and simulate environmental feedback. However, here, the age
component threshold ($a_{T}$) is not fixed at 50% of the plants age, but is
instead variable with a system parameter. As well as this, a positive feedback
component was added to increase the viability of plants with a larger radius,
as the initial implementation was found to lead to the decline (and often
species eradication) of larger radius trees. This is likely due to the fact
that trees with larger canopies are more likely to collide with another plant,
however, in nature, larger trees are often hardier than smaller ones.

After the simulation has undergone $N$ iterations, a given plant’s mesh is
generated with $n_{i}$ iterations of the L-system:

$n_{i}=\lfloor\widehat{a}\,(i_{\text{max}}+1-i_{\text{min}})\rfloor$ (5)

where $i_{\text{max}}$ and $i_{\text{min}}$ are the maximum and minimum number
of L-system iterations for the plant’s species, and $\widehat{a}$ is the
normalised age for the plant:

$\widehat{a}=\text{min}\left(\frac{a}{a_{M}},\,1\right),$ (6)

where $a$ is the plant’s current agent and $a_{M}$ is the maximum age for the
plant’s species. Next a mesh scale factor $s$ is used to smoothly interpolate
between iteration counts:

$s=\text{lerp}(s_{\text{min}},\;s_{\text{max}},\;\widehat{a}\;(i_{\text{max}}+1-i_{\text{min}})-n_{i})$
(7)

in which $s_{\text{max}}$ and $s_{\text{min}}$ are the upper and lower bounds
on the scaling factor for the given plant’s species and lerp is the standard
linear interpolation function. Finally, a random rotation around the vertical
axis is applied—as in Section 4.2—to vary the orientations of trees and
prevent the appearance of repetition across model instances.

##  5 Rendering

To convert the scene meshes into screen-space pixel values, the system carries
out six render passes: one shadow pre-pass before the main render loop, one
geometry pass, one lighting pass and three post-processing passes. The system
uses a deferred rendering pipeline: framebuffers store the result of each pass
to be used later. Therefore expensive lighting calculations are only performed
for visible pixels. This also enables pixel-level post-processing effects such
as screen-space ambient occlusion (SSAO) and volumetric light scattering
(Section 5.2).

The geometry pass uses a number of vertex and fragment shaders for the
different object types. The vertex shaders convert the local mesh vertex data
into a global vertex position $\bm{p}^{g}$, a TBN (tangent-bitangent-normal)
matrix $\bm{T}$ and a re-normalised surface normal $\widetilde{\bm{n}}$
[dV20]. Normal mapping is then carried out in the fragment shaders by sampling
fragments from the texture map, which are then multiplied by $\bm{T}$ to
convert from tangent space to world space. For fragments which are part of a
billboard, the map normal ($\widehat{\bm{n^{m}}}$) is negated if the ‘back’
face of the board (the face that the board normal ($\widetilde{\bm{n}}$) is
facing into) is being viewed—to create a board with separate front and back
faces:

$\widehat{\bm{n^{m}}}=\begin{cases}-\widehat{\bm{n^{m}}},&\text{if
}\widetilde{\bm{n}}\cdot\widehat{\bm{v}}<0\\\
\widehat{\bm{n^{m}}},&\text{otherwise}\end{cases}$ (8)

in which $\widehat{\bm{v}}$ is the vector from the fragment position to the
camera position. Finally, alpha-channel fragment discarding is carried out to
enable transparency, and the position, normal, colour and occlusion values are
written to framebuffers.

As well as rendering the scene models, the geometry pass also renders an HDRI
skybox with a ‘clip-space to normalised-device-space’ coordinate
transformation to display the skybox at infinite distance [dV20]. Before the
main render loop, the system converts the equirectangular HDRI projection into
six 2D textures for the skybox faces.

###  5.1 Lighting

In the lighting pass, a Blinn-Phong lighting model [Bli77] is used with
percentage-closer-filtering (PCF) shadow mapping [RSC87], and a specific
lighting model with translucency for leaves. The only light source in the
scene is the sun, which is modelled as a directional light with position
$\bm{p}_{S}$ and colour (strength) $\bm{c}_{S}$.

Whilst the sun position ($\bm{p}_{S}$) can be set manually by a system
parameter, to enable easy changing of light conditions, the position can also
be estimated from the HDRi skybox image. The position is estimated by scanning
over each pixel in the HDRi texture and summing the RGB values for each pixel
in a $3\times 3$ kernel around the pixel. The pixel with the largest sum value
is then converted into a Cartesian coordinate and scaled to a suitable
distance from the scene to give the sun’s position vector.

From the sun, the reflected light $L_{o}$ at a point $\bm{x}$ in direction
$\bm{\omega}_{0}$ is given by:

$L_{o}(\bm{x},\,\bm{\omega}_{0})=L_{a}(\bm{x})+\left(1-s\left(\bm{x},\,\bm{\omega}_{0}\right)\right)\left(L_{d}\left(\bm{x}\right)+L_{s}\left(\bm{x},\,\bm{\omega}_{0}\right)\right)$
(9)

where $s$ is the shadow coefficient (1 in full shadow, 0 in full light),
$L_{d}$ is the standard Phong diffuse component [Pho75], $L_{s}$ is the
standard Blinn-Phong specular (glossy) component [Bli77], and $L_{a}$ is an
ambient component used to approximate indirect lighting,

$L_{a}(\bm{x})=\alpha L_{ao}(\bm{x})\bm{c}_{d}$ (10)

in which $\bm{c}_{d}$ is the diffuse colour of the fragment, $\alpha$ is a
constant ambient factor and $L_{ao}$ is the screen space ambient occlusion
term (Section 5.2).

A single shadow map is generated for the whole scene, in a pre-render pass,
where the scene is rendered with an orthogonal projection matrix $\bm{P_{s}}$
which defines a view cuboid with width (and height) $w$, near plane distance
$d_{near}=0.1$ and far plane distance $d_{far}$:

$w=\text{max}\left(w_{t}\sqrt{2},\,50\right)$ (11)
$d_{far}=\max_{\bm{c}\in\bm{C}}\left(\lvert\lvert\bm{p}_{s}-\bm{c}\rvert\rvert\right)+10$
(12)

where $w_{t}$ is the terrain width, $\bm{p}_{s}$ is the sun position and
$\bm{C}$ is the set of corners of the terrain:

$\displaystyle\bm{C}=\biggl{\\{}\,\frac{1}{2}\begin{bmatrix}-w_{t},&0,&-w_{t}\end{bmatrix}^{T}$
$\displaystyle,\quad\frac{1}{2}\begin{bmatrix}-w_{t},&0,&w_{t}\end{bmatrix}^{T},\,$
(13)
$\displaystyle\frac{1}{2}\begin{bmatrix}w_{t},&0,&-w_{t}\end{bmatrix}^{T}$
$\displaystyle,\quad\frac{1}{2}\begin{bmatrix}w_{t},&0,&w_{t}\end{bmatrix}^{T}\,\biggr{\\}}$

These values of $w$ and $d_{far}$ ensure that the entire scene is captured
within the view cuboid of the sun, and an orthogonal projection matrix is used
to simulate parallel rays of light from the sun. Along with $\bm{P_{s}}$, a
view matrix $\bm{V_{s}}$ is used to orientate the scene from the sun’s
position.

The shadow map is then used during the lighting pass to determine if a
fragment is in shadow or not ($s(\bm{x})$), and PCF is carried out by
averaging the shadow coefficients for the 9 pixels of the shadow map in the
$3\times 3$ kernel around the current fragment.

The lighting model for leaves includes the transparency model proposed in
[HKW07] used alongside the Blinn-Phong shading model. Each leaf is modelled as
a two-sided plane (using normal negation as in Equation 8), and has four
textures associated with each side (front and back):

  1. 1.

An albedo (diffuse) colour texture.

  2. 2.

A translucency colour texture which holds the colour of translucent light
through the leaf.

  3. 3.

A normal map.

  4. 4.

A half-life basis map which holds the (pre-computed) translucency factor for
each point in the leaf represented in the HL2 basis [McT04] created with the
method proposed in [HKW07].

The resulting light from a leaf fragment is given by the following equation
(adapted from [HKW07] for a Blinn-Phong pipeline):

$L_{\text{o-leaf}}(\bm{x},\,\bm{\omega}_{0})=L_{o}(\bm{x},\,\bm{\omega}_{0})+\beta\bm{c}_{t}\odot\bm{c}_{S}\sqrt{\frac{3}{2\pi}}\sum_{i=1}^{3}[\bm{h}_{l}]_{i}(\bm{H}_{i}\cdot\bm{\omega}_{S})$
(14)

where $L_{o}$ is the standard lighting model for the system (Equation 9) ,
$\bm{h}_{l}$ is the half-life basis map value, $\bm{c}_{t}$ is the
translucency colour map value, $\beta$ is a global translucency factor scalar,
and $\bm{H}_{i},\,i\in\\{1\dots 3\\}$ are the HL2 basis vectors [HKW07].

###  5.2 Post-processing

The system uses the implementation of screen-space ambient occlusion (SSAO)
presented in [dV20] where the post-processing shader operates on fragment
depth values and samples a number of values in a normal-aligned hemisphere
around each fragment to determine how exposed the point is to indirect
lighting, giving an ambient occlusion term $L_{ao}$. A second shader is then
used to blur the SSAO texture to reduce noise artefacts.

To implement volumetric light scattering (Figure 5),

Figure 5: Examples of light scattering for different environment maps and
positions with several tree species.

as a screen-space post-processing shader, the system uses the method presented
in [Mit07] adapted for use with HDRi environment maps. During the geometry
pass, an occlusion map is created with the skybox fragments in full (HDR)
colour, and all other fragments as black. The occlusion texture is then passed
through the scattering shader which uses Equation 15 to sample a number of
points along the vector from the sun to each fragment (with an exponential
decay factor). As an adaptation from the original method, the sample
brightness (length of the colour vector) is clamped at a fixed value $s_{max}$
to prevent samples from appearing as bright white circles due to the high
dynamic range of the skybox, giving

$L_{scatter}(\bm{x})=E\sum_{i=0}^{N-1}d^{i}\text{min}\left(\text{occ}\left(\bm{x}^{ss}-\left(i+1\right)\bm{\delta}_{s}\right),\,s_{max}\right)$
(15)

where

$\bm{\delta}_{s}=\frac{\bm{x}^{ss}-\bm{p}_{s}^{ss}}{\rho N},$ (16)

$\text{occ}(\bm{v})$ samples the occlusion texture at coordinate $\bm{v}$,
$\bm{x}^{ss}$ is the screen-space position of $\bm{x}$, $\bm{p}_{s}^{ss}$ is
the screen-space position of the sun, $E$ is an exposure (brightness) scalar,
$N$ is the number of samples to take, $d\in[0,1)$ is the decay coefficient,
and $\rho$ is the sample density (a higher density results in samples being
taken from a smaller area which produces shorter, brighter rays).

This method works with any sun position or environment map and the light rays
disperse well through small gaps in the leaves, with the decay parameter
creating a realistic fall-off in light intensity. Furthermore, the rays update
in position and appearance in a realistic manner as the camera moves around
the scene.

To enable High Dynamic Range (HDR) calculations, the system uses floating
point framebuffers and then tone maps the output to return the values to the
range $[0,1]$. The exposure tone mapping operator presented in [dV20] is used
to allow the user to control the tone mapping with the exposure parameter $E$

$c_{LDR}=1-\exp{(-E\cdot c_{HDR})}$ (17)

where $c_{HDR}$ is the high-dynamic range colour value.

The system also uses gamma correction to convert the textures from sRGB colour
space [IEC99] to linear colour space for calculations and then back to sRGB
space for display on the screen. The gamma correction is carried out alongside
tone mapping at the end of the lighting pass fragment shader.

###  5.3 Rendering Optimisation

The level of detail (LOD) representations for each object type were chosen so
they could be pre-generated (Section Section 4). The representations for each
model are stored in a hash map to enable constant-time look-ups at render
time.

The LOD of a model is determined by its distance from the camera, however, to
prevent calculating the distance to each one individually, a LOD is determined
for each visible quadtree node (quad) based on its depth in the tree, and the
models within that quad are then all represented at that LOD. Visible quads
are determined recursively, where recursion stops if condition 18 is
satisfied. This has the effect that closer quads are subdivided more,
resulting in smaller (deeper) quads closer to the camera.

$\left(\lvert q_{x}-p_{x}\rvert>\tau_{d}w\right)\,\land\,\left(\lvert
q_{z}-p_{z}\rvert>\tau_{d}w\right)\,\lor\,d=D$ (18)

where $(q_{x},q_{z})$ is the centre of the quad on the XZ-plane,
$(p_{x},p_{z})$ is the camera position on the XZ-plane, $w$ is the width of
the quad, $d$ is the depth of the quad in the tree, $D$ is the maximum depth
of the quadtree and $\tau_{d}$ is a threshold scalar which controls the radius
of division around the camera.

(a) $n=8$

(b) $n=10$

(c) $n=13$

Figure 6: Three stages of growth ($n$ iterations) for the ‘oak’ species of
branching tree. The width and length of branches, as well as the number of
branches and leaves, can be seen to increase as the tree ‘grows’.

Following [GH01], view-frustum culling is also performed at this point to
check whether either of two spheres—of radius $2w$ centred at
$(q_{x},0,q_{z})$ and $(q_{x},2w,q_{z})$—intersect with the frustum of the
matrix $\bm{PV}$. Where $\bm{P}$ is the projection matrix for the scene and
$\bm{V}$ is the current camera view matrix. This means that quads which are
not inside the view frustum (plus an offset of 1 quad width to allow for trees
which overlap quads) are not rendered.

The current LOD system has two levels (high and low)—but could be easily
expanded to an arbitrary $N$ levels—with the deepest two quadtree levels ($D$
and $D-1$) rendered with high LOD, as they are closest to the camera.

##  6 Results and Discussion

###  6.1 Scene Quality

The use of tree-specific L-systems[PL90] allows for a broad range of tree
species to be represented and rendered (Figure Procedural Generation and
Rendering of Realistic, Navigable Forest Environments: An Open-Source Tool).
The system generates realistic tree models with highly configurable
parameters, and natural-looking leaf and branch distributions due to the
varying angles and rotations. One slight issue with the use of L-systems is
that branches can occasionally self-intersect in particularly dense areas.
This is only noticeable at close-distance, however, and could be prevented by
using an environmentally-sensitive L-System [PJM94], a branch-level-parametric
tree model [MHS∗19], or an intersection removal technique [XYS∗16], albeit at
a higher computational cost.

In the L-systems used, an increase in iterations results in increasingly later
growth stages of the resulting tree (Figure 6)—with the number, length and
width of branches increasing—enabling tree models of differing ages to be
generated from the same L-system and parameter set.

The LOD reduction techniques used, reduce the number of vertices in the tree
models by an average of 63% (see supplementary information for numerical
breakdown)whilst producing visually similar results in the far and medium
distances (Figure 7) with the overall tree shapes and individual branches
being retained in the low LOD representation, along with the leaves exhibiting
the same lighting behaviour and coverage in both representations.

Figure 7: Comparison of low LOD (top) and high LOD (bottom) representations of
trees and ground cover. Whilst differences between the representations are
discernible at close distance, middle and far-ground objects are visually
similar across both representations.

(a) $n=0$

(b) $n=1000$

Figure 8:  The initial (random) state of the ecosystem simulation, and the
distributed plants after $n=1000$ iterations. Clustering of plant species can
be seen to emerge, and these clusters remain over time, just changing
position.

The use of an ecosystem simulation produces tree distributions which exhibit
natural features, such as plant species clustering and resource consumption,
whilst also resulting in a pleasing and natural-looking forest from within.

Figure 9: Renders of a scene with low quality (LQ) settings (left) compared to
high quality (HQ) settings (right). The scene has a terrain width $w=100$, and
quadtree depth $D=4$, with the LQ scene achieving an average frame render time
of 62.70 ms (and maximum 66.52 ms) across three different views, compared to
the HQ which required an average of 214.85 ms.

Figure 8 shows the distribution produced by the simulation after 1000
iterations, where plant clustering can be seen to have emerged—unlike the
initial (random) distribution in the first image. Furthermore, the simulation
prevents intersecting trees and larger trees inhibit the growth of other
plants around them to simulate resource competition. The system also handles
distributing a single plant species well (Figure 10), resulting in a realistic
distribution from both above and within the forest.

Figure 10: Ground view (left) and distribution (right) for a forest with a
single species.

###  6.2 Performance

Results were measured on an Intel Core i7-8750H CPU, 2.20GHz with an NVIDIA
GeForce GTX 1050 Ti Max-Q at a resolution of $1400\times 780$ px.

The average render-time per frame across 6 runs of varying terrain sizes and
quadtree depths ranged from 0.0643 s (15.5 FPS) to 0.308 s (3.24 FPS). The
most time-consuming part of the render cycle is the geometry pass—taking an
average of 94.8% of the total time across the runs. This suggests that the
main limiting factor for the frame rate is the number of vertices in the
scene. To control this, the system parameters allow for a trade-off between
quality and render time, enabling real-time interactive frame rates in the
manual mode, and higher-quality output in the asynchronous input mode. Two
extremes in quality are shown with their associated average frame render times
in Figure 9—the low-quality frames were rendered three times faster than the
high-quality frames.

##  7 Conclusions

We provide a system—released open-source[New22]—for generating and rendering
realistic virtual forests which can be interactively explored. The system is
fully configurable through the use of text files and generates the full scene
before render time to reduce per-frame computation. Trees are generated with
L-systems and distributed throughout the scene using an ecosystem simulation.
The resulting scene is then rendered with a deferred rendering pipeline, using
a Blinn-Phong shading model adapted for leaf transparency; along with high-
dynamic-range tone mapping, shadow mapping and post-processing shaders for
ambient occlusion and volumetric light scattering.

The system is able to generate an infinite number of tree species, terrains
and overall forest scenes due to the extensive configurability of the system
and flexibility of the techniques implemented. The generated tree models and
distribution are naturalistic and believable, and a high-quality render can be
achieved. The render performance shows that small scenes can be navigated in
real-time and that, for larger scenes, the render quality can be reduced to
enable real-time navigation at reasonable frame rates on a high-end laptop or
PC. The system can also run asynchronously—with non-real-time output—which
allows for high-quality images or video frames to be produced, or for external
software (such as machine learning algorithms) to interface with the Forest
Generator. For example, the tool is currently being used for training computer
vision for autonomous rovers.

A future development—to increase the usability of the application—could be to
adapt the system to utilise higher-level parameters [SGdK∗11, CCGS19] (which
are easier to understand and control), or to apply machine learning techniques
to automatically determine optimised parameter values from inputs such as
photographs of real forests or user sketches [OSKL05, RT20].

## References

  * [AD06] Alsweis M., Deussen O.:  Wang-Tiles for the simulation and visualization of plant competition.  In _Advances in Computer Graphics_ (Hangzhou, China, June 2006), Nishita T., Peng Q., Seidel H., (Eds.), vol. 4035 of _Lecture Notes in Computer Science_ , Springer-Verlag, pp. 1–11.  doi:10.1007/11784203_1.
  * [BA05] Belhadj F., Audibert P.:  Modeling landscapes with ridges and rivers: bottom up approach.  In _Proceedings of the 3rd International Conference on Computer Graphics and Interactive Techniques in Australasia and South East Asia_ (Dunedin, New Zealand, Nov. 2005), GRAPHITE ’05, Association for Computing Machinery, pp. 447–450.  doi:10.1145/1101389.1101479.
  * [Bar88] Barnsley M.:  _Fractals everywhere_.  Academic Press Professional, Inc., USA, 1988.
  * [BBP08] Boulanger K., Bouatouch K., Pattanaik S.:  Rendering trees with indirect lighting in real time.  _Computer Graphics Forum 27_ , 4 (2008), 1189–1198.  doi:https://doi.org/10.1111/j.1467-8659.2008.01257.x.
  * [BCF∗05] Behrendt S., Colditz C., Franzke O., Kopf J., Deussen O.:  Realistic real-time rendering of landscapes using billboard clouds.  _Computer Graphics Forum 24_ , 3 (2005), 507–516.  doi:10.1111/j.1467-8659.2005.00876.x.
  * [Ben02] Beneš B.:  A stable modeling of large plant ecosystems.  In _Proceedings of the International Conference on Computer Vision and Graphics_ (Zakopane, Sept. 2002), ICCVG ’02, pp. 94–101.
  * [BEP18] Backes G. C., Engel T. A., Pozzer C. T.:  Real-time massive terrain generation using procedural erosion on the GPU.  In _Proceedings of 17th Brazilian Symposium on Computer Games and Digital Entertainment (SBGames) 2018_ (Brazil, Oct. 2018), SBGames ’18, pp. 675–678.
  * [Bey15] Beyer H. T.:  _Implementation of a method for hydraulic erosion_.  bathesis, Technische Universität München, Munich, Germany, Nov. 2015.
  * [BF01] Beneš B., Forsbach R.:  Layered data representation for visual simulation of terrain erosion.  In _Proceedings of the 17th Spring Conference on Computer Graphics_ (USA, Apr. 2001), SCCG ’01, IEEE Computer Society, p. 80.
  * [BF02] Beneš B., Forsbach R.:  Visual simulation of hydraulic erosion.  _Journal of WSCG 10_ , 2 (Feb. 2002), 79–86.
  * [BH00] Berger U., Hildenbrandt H.:  A new approach to spatially explicit modelling of forest dynamics: spacing, ageing and neighbourhood competition of mangrove trees.  _Ecological Modelling 132_ , 3 (Apr. 2000), 287–302.  doi:10.1016/S0304-3800(00)00298-2.
  * [Bli77] Blinn J. F.:  Models of light reflection for computer synthesized pictures.  In _Proceedings of the 4th Annual Conference on Computer Graphics and Interactive Ttechniques_ (San Jose, California, July 1977), SIGGRAPH ’77, Association for Computing Machinery, pp. 192–198.  doi:10.1145/563858.563893.
  * [BLZ∗11] Bao G., Li H., Zhang X., Che W., Jaeger M.:  Realistic real-time rendering for large-scale forest scenes.  In _2011 IEEE International Symposium on VR Innovation_ (Los Alamitos, CA, USA, Mar. 2011), IEEE Computer Society, pp. 217–223.  doi:10.1109/ISVRI.2011.5759637.
  * [BN08] Bruneton E., Neyret F.:  Real-time rendering and editing of vector-based terrains.  _Computer Graphics Forum 27_ , 2 (2008), 311–320.  doi:https://doi.org/10.1111/j.1467-8659.2008.01128.x.
  * [BN12] Bruneton E., Neyret F.:  Real-time realistic rendering and lighting of forests.  _Computer Graphics Forum 31_ , 2 (2012), 373–382.  doi:10.1111/j.1467-8659.2012.03016.x.
  * [BNB13] Bradley D., Nowrouzezahrai D., Beardsley P.:  Image-based reconstruction and synthesis of dense foliage.  _ACM Transactions on Graphics 32_ , 4 (July 2013), 74.  doi:10.1145/2461912.2461952.
  * [BTHB06] Beneš B., Těšínský V., Hornyš J., Bhatia S. K.:  Hydraulic erosion.  _Computer Animation and Virtual Worlds 17_ , 2 (Feb. 2006), 99–108.  doi:10.1002/cav.77.
  * [BvMM11] Beneš B., Št’ava O., Měch R., Miller G.:  Guided procedural modeling.  _Computer Graphics Forum 30_ , 2 (2011), 325–334.  doi:https://doi.org/10.1111/j.1467-8659.2011.01886.x.
  * [BZCJ09] Bao G., Zhang X., Che W., Jaeger M.:  Billboards for tree simplification and real-time forest rendering.  In _2009 Third International Symposium on Plant Growth Modeling, Simulation, Visualization and Applications_ (Nov. 2009), Li B., Jaeger M., Guo Y., (Eds.), IEEE, pp. 433–440.  doi:10.1109/PMA.2009.38.
  * [CCGS19] Cook M., Colton S., Gow J., Smith G.:  General analytical techniques for parameter-based procedural content generators.  In _2019 IEEE Conference on Games (CoG)_ (2019), pp. 1–8.  doi:10.1109/CIG.2019.8848024.
  * [CCH05] Candussi A., Candussi N., Höllerer T.:  Rendering realistic trees and forests in real time.  In _26th Annual Conference of the European Association for Computer Graphics, Eurographics 2005 - Short Presentations_ (2005), Dingliana J., Ganovelli F., (Eds.), The Eurographics Association, pp. 73–76.  doi:10.2312/egs.20051027.
  * [CMF98] Chiba N., Muraoka K., Fujita K.:  An erosion model based on velocity fields for the visual simulation of mountain scenery.  _The Journal of Visualization and Computer Animation 9_ , 4 (1998), 185–194.
  * [Coo86] Cook R. L.:  Stochastic sampling in computer graphics.  _ACM Transactions on Graphics 5_ , 1 (Jan. 1986), 51–72.  doi:10.1145/7529.8927.
  * [dB00] de Boer W. H.:  _Fast Terrain Rendering Using Geometrical MipMapping_.  Tech. rep., E-mersion Project, October 2000, Oct. 2000.
  * [DCDS05] Dietrich A., Colditz C., Deussen O., Slusallek P.:  Realistic and interactive visualization of high-density plant ecosystems.  In _Eurographics Workshop on Natural Phenomena_ (2005), Poulin P., Galin E., (Eds.), The Eurographics Association.  doi:10.2312/NPH/NPH05/073-081.
  * [DCSD02] Deussen O., Colditz C., Stamminger M., Drettakis G.:  Interactive visualization of complex plant ecosystems.  In _Proceedings of the Conference on Visualization ’02_ (USA, 2002), VIS ’02, IEEE Computer Society, pp. 219–226.  doi:10.1109/VISUAL.2002.1183778.
  * [DDSD03] Décoret X., Durand F., Sillion F. X., Dorsey J.:  Billboard clouds for extreme model simplification.  _ACM Transactions on Graphics 22_ , 3 (July 2003), 689–696.  doi:10.1145/882262.882326.
  * [DGCF01] Di Giacomo T., Capo S., Faure F.:  An interactive forest.  In _Eurographics Workshop on Computer Animation and Simulation (EGCAS)_ (Manchester, United Kingdom, Sept. 2001), Magnenat-Thalmann N., Thalmann D., (Eds.), Eurographics, Springer, pp. 65–74.  doi:10.1007/978-3-7091-6240-8_7.
  * [DHL∗98] Deussen O., Hanrahan P., Lintermann B., Mech R., Pharr M., Prusinkiewicz P.:  Realistic modeling and rendering of plant ecosystems.  In _Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, July 1998), SIGGRAPH, Association for Computing Machinery, pp. 275–286.  doi:10.1145/280814.280898.
  * [DN04] Decaudin P., Neyret F.:  Rendering forest scenes in real-time.  In _EGSR04: 15th Eurographics Symposium on Rendering_ (Norrköping, Sweden, June 2004), Keller A., Jensen H. W., (Eds.), The Eurographics Association, pp. 93–1202.  doi:10.2312/EGWR/EGSR04/093-102.
  * [DN09] Decaudin P., Neyret F.:  Volumetric billboards.  _Computer Graphics Forum 28_ , 8 (Dec. 2009), 2079–2089.  doi:10.1111/j.1467-8659.2009.01354.x.
  * [dNFP18] do Nascimento B. T., Franzin F. P., Pozzer C. T.:  GPU-based real-time procedural distribution of vegetation on large-scale virtual terrains.  In _Proceedings of 17th Brazilian Symposium on Computer Games and Digital Entertainment (SBGames) 2018_ (Brazil, Oct. 2018), SBGames ’18, pp. 157–166.  ISSN: 2159-6662.  doi:10.1109/SBGAMES.2018.00027.
  * [DP10] Doran J., Parberry I.:  Controlled Procedural Terrain Generation Using Software Agents.  _IEEE Transactions on Computational Intelligence and AI in Games 2_ , 2 (June 2010), 111–119.  doi:10.1109/TCIAIG.2010.2049020.
  * [Dun16] Dunn I.:  _Procedural Generation and Rendering of Large-scale Open-world Environments_.  mathesis, California Polytechnic State University, San Luis Obispo, CA, USA, Dec. 2016.  doi:10.15368/theses.2016.163.
  * [dV20] de Vries J.:  _Learn OpenGL: Learn Modern OpenGL Graphics Programming in a Step-by-step Fashion_.  Kendall & Welling, 2020.  URL: https://learnopengl.com/.
  * [DWS∗97] Duchaineau M., Wolinsky M., Sigeti D. E., Miller M. C., Aldrich C., Mineev-Weinstein M. B.:  ROAMing terrain: Real-time optimally adapting meshes.  In _Proceedings of IEEE Conference on Visualization ’97_ (Oct. 1997), IEEE Computer Society, pp. 81–88.  doi:10.1109/VISUAL.1997.663860.
  * [DZYJ10] Deng Q., Zhang X., Yang G., Jaeger M.:  Multiresolution foliage for forest rendering.  _Computer Animation and Virtual Worlds 21_ , 1 (2010), 1–23.  doi:10.1002/cav.283.
  * [FE17] Freiknecht J., Effelsberg W.:  A survey on the procedural generation of virtual worlds.  _Multimodal Technologies and Interaction 1_ , 4 (2017).  doi:10.3390/mti1040027.
  * [FUM05] Fuhrmann A. L., Umlauf E., Mantler S.:  Extreme model simplification for forest rendering.  In _Proceedings of the First Eurographics Conference on Natural Phenomena_ (Goslar, DEU, 2005), NPH’05, The Eurographics Association, pp. 57–67.
  * [GH01] Gribb G., Hartmann K.:  _Fast Extraction of Viewing Frustum Planes from the World-View-Projection Matrix_.  Tech. rep., June 2001.  https://www.gamedevs.org/uploads/fast-extraction-viewing-frustum-planes-from-world-view-projection-matrix.pdf.
  * [GM01] Gamito M. N., Musgrave F. K.:  Procedural landscapes with overhangs.  In _Proceedings of the 10th Portuguese Meeting of Computer Graphics and Interaction_ (Lisbon, Portugal, Oct. 2001), pp. 33–42.  doi:10.1.1.120.6073.
  * [GMS09] Gain J., Marais P., Straßer W.:  Terrain sketching.  In _Proceedings of the 2009 symposium on Interactive 3D graphics and games_ (Boston, Massachusetts, Feb. 2009), I3D ’09, Association for Computing Machinery, pp. 31–38.  doi:10.1145/1507149.1507155.
  * [GS08] Geist R., Steele J.:  A lighting model for fast rendering of forest ecosystems.  In _2008 IEEE Symposium on Interactive Ray Tracing_ (Los Angeles, CA, USA, Aug. 2008), IEEE, pp. 99–106.  doi:10.1109/RT.2008.4634629.
  * [Häg06] Häggström H.:  _Real-time generation and rendering of realistic landscapes_.  Master’s thesis, University of Helsinki, Helsinki, FIN, Aug. 2006.
  * [Ham01] Hammes J.:  Modeling of ecosystems as a data source for real-time terrain rendering.  In _Digital Earth Moving_ (Berlin, Heidelberg, Sept. 2001), Westort C. Y., (Ed.), vol. 2181 of _Lecture Notes in Computer Science_ , Springer, pp. 98–111.  doi:10.1007/3-540-44818-7_14.
  * [Hew17] Hewitt C.:  _Procedural generation of tree models for use in computer graphics_.  mathesis, Trinity Hall, University of Cambridge, Cambridge, England, May 2017.
  * [HG97] Heckbert P. S., Garland M.:  Survey of polygonal surface simplification algorithms.  In _Multiresolution Surface Modeling - ACM SIGGRAPH 1997 Course 25 Notes_ (May 1997), no. 25 in SIGGRAPH Course Notes, Association for Computing Machinery, pp. 1:1–1:29.
  * [HGA∗10] Hnaidi H., Guérin E., Akkouche S., Peytavie A., Galin E.:  Feature based terrain generation using diffusion equation.  _Computer Graphics Forum 29_ , 7 (2010), 2179–2186.  doi:https://doi.org/10.1111/j.1467-8659.2010.01806.x.
  * [HKW07] Habel R., Kusternig A., Wimmer M.:  Physically based real-time translucency for leaves.  In _Rendering Techniques 2007 (Proceedings Eurographics Symposium on Rendering)_ (June 2007), Kautz J., Pattanaik S., (Eds.), Eurographics, Eurographics Association, pp. 253–263.
  * [HRS13] Hempe N., Rossmann J., Sondermann B.:  Generation and rendering of interactive ground vegetation for real-time testing and validation of computer vision algorithms.  _ELCVIA Electronic Letters on Computer Vision and Image Analysis 12_ , 2 (Apr. 2013), 40–53.  doi:10.5565/rev/elcvia.519.
  * [HSAKD98] House D. H., Schmidt G. S., Arvin S. A., Kitagawa-DeLeon M.:  Visualizing a real forest.  _IEEE Computer Graphics and Applications 18_ , 1 (Jan. 1998), 12–15.  doi:10.1109/38.637262.
  * [HTF82] Honda H., Tomlinson P. B., Fisher J. B.:  Two geometrical models of branching of botanical trees.  _Annals of Botany 49_ , 1 (Jan. 1982), 1–11.  doi:10.1093/oxfordjournals.aob.a086218.
  * [IEC99] _IEC 61966-2-1:1999 - Multimedia systems and equipment - Colour measurement and management - Part 2-1: Colour management - Default RGB colour space - sRGB_.  Standard, International Electrotechnical Commission, Oct. 1999.
  * [Int20] Interactive Data Visualization, Inc:  Speedtree, 2020.  URL: https://store.speedtree.com/.
  * [Jak00] Jakulin A.:  Interactive vegetation rendering with slicing and blending.  In _21st Annual Conference of the European Association for Computer Graphics, Eurographics 2000 - Short Presentations_ (Aug. 2000), de Sousa A. A., Torres J. C., (Eds.), Eurographics Association.  doi:10.2312/egs.20001004.
  * [JRSK09] Jens O., Rezk Salama C., Kolb A.:  GPU-based responsive grass.  _Journal of WSCG 17_ , 1-3 (2009), 65–72.
  * [JT11] Jákó B., Tóth B.:  Fast hydraulic and thermal erosion on GPU.  In _Eurographics 2011 - Short Papers_ (2011), Avis N., Lefebvre S., (Eds.), The Eurographics Association.  doi:10.2312/EG2011/short/057-060.
  * [Kah13] Kahoun M.:  _Realtime library for procedural generation and rendering of terrains_.  Master’s thesis, Charles University in Prague, Prague, CZE, Aug. 2013\.
  * [KBKv09] Krištof P., Beneš B., Krivánek J., Št’ava O.:  Hydraulic erosion using smoothed particle hydrodynamics.  _Computer Graphics Forum 28_ , 2 (2009), 219–228.  doi:https://doi.org/10.1111/j.1467-8659.2009.01361.x.
  * [KCS07] Kharlamov A., Cantlay I., Stepanenko Y.:  _GPU Gems 3_.  No. 3 in GPU Gems. Addison-Wesley Professional, 2007, ch. Chapter 4\. Next-Generation SpeedTree Rendering.
  * [KGM14] Kenwood J., Gain J., Marais P.:  Efficient procedural generation of forests.  _Journal of WSCG 22_ , 1 (2014), 31–38.
  * [Kom20] Komatsu Forest AB:  Komatsu simulators, 2020.  URL: https://www.komatsuforestsimulators.com/.
  * [KS18] Kohek v., Strnad D.:  Interactive large-scale procedural forest construction and visualization based on particle flow simulation.  _Computer Graphics Forum 37_ , 1 (2018), 389–402.  doi:10.1111/cgf.13304.
  * [Kur12] Kurowski M.:  Procedural generation of meandering rivers inspired by erosion.  In _Proceedings of The 20th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision_ (Plzen, Czech Republic, June 2012), Skala V., (Ed.), vol. 1 of _WSCG_ , Vaclav Skala – Union Agency, pp. 79–86.
  * [Lan20] Landin P.:  The road to dungeons: Exploring the random art of procedural levels.  _Minecraft.net Deep Dives_ (May 2020).  URL: https://www.minecraft.net/en-us/article/the-road-dungeons.
  * [LD99] Lintermann B., Deussen O.:  Interactive modeling of plants.  _IEEE Computer Graphics and Applications 19_ , 1 (Jan. 1999), 56–65.  doi:10.1109/38.736469.
  * [Len20] Lenné3D:  Selection of our projects, 2020.  URL: https://www.lenne3d.com/de/alle-projekte/.
  * [Lin68] Lindenmayer A.:  Mathematical models for cellular interactions in development, Parts I and ll.  _Journal of Theoretical Biology 18_ , 3 (Mar. 1968), 280–315.  doi:10.1016/0022-5193(68)90079-9.
  * [LKR∗96] Lindstrom P., Koller D., Ribarsky W., Hodges L. F., Faust N., Turner G. A.:  Real-time, continuous level of detail rendering of height fields.  In _Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, Aug. 1996), SIGGRAPH ’96, Association for Computing Machinery, pp. 109–118.  doi:10.1145/237170.237217.
  * [LP02] Lane B., Prusinkiewicz P.:  Generating spatial distributions for multilevel models of plant communities.  In _Proceedings of the Graphics Interface 2002 Conference_ (Calgary, Alberta, May 2002), vol. Calgary, Canadian Human-Computer Communications Society, pp. 69–80.  doi:10.20380/GI2002.09.
  * [LWJ20] LWJGL: Lightweight Java Game Library, 2020.  URL: https://www.lwjgl.org/.
  * [Max96] Max N.:  Hierarchical rendering of trees from precomputed multi-layer z-buffers.  In _Proceedings of the Eurographics Workshop on Rendering Techniques ’96_ (Vienna, Aug. 1996), Pueyo X., Schröder P., (Eds.), Eurographics, Springer-Verlag, pp. 165–174.
  * [McT04] McTaggart G.:  _Half-Life ® 2 / Valve SourceTM Shading_.  Tech. rep., Valve Corporation, 2004.
  * [MDH07] Mei X., Decaudin P., Hu B.:  Fast hydraulic erosion simulation and visualization on GPU.  In _15th Pacific Conference on Computer Graphics and Applications_ (Maui, Hawaii, Oct. 2007), Alexa M., Gortler S., Ju T., (Eds.), Pacific Graphics, pp. 47–56.  doi:10.1109/PG.2007.15.
  * [MDK99] Max N., Deussen O., Keating B.:  Hierarchical image-based rendering using texture mapping hardware.  In _Proceedings of the Eurographics Workshop on Rendering Techniques ’99_ (1999), Lischinski D., Larson G. W., (Eds.), Eurographics, The Eurographics Association, pp. 57–62.  doi:10.2312/EGWR/EGWR99/057-062.
  * [MHS∗19] Makowski M., Hädrich T., Scheffczyk J., Michels D. L., Pirk S., Pałubicki W.:  Synthetic silviculture: Multi-scale modeling of plant ecosystems.  _ACM Trans. Graph. 38_ , 4 (jul 2019).  doi:10.1145/3306346.3323039.
  * [Mit07] Mitchell K.:  _GPU Gems 3_.  No. 3 in GPU Gems. Addison-Wesley Professional, 2007, ch. Chapter 13\. Volumetric Light Scattering as a Post-Process.
  * [MKM89] Musgrave F. K., Kolb C. E., Mace R. S.:  The synthesis and rendering of eroded fractal terrains.  _ACM SIGGRAPH Computer Graphics 23_ , 3 (July 1989), 41–50.  doi:10.1145/74334.74337.
  * [MM16] Mikulicic N., Mihajlovic v.:  Procedural generation of mediterranean environments.  In _39th International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO)_ (Opatija, HRV, May 2016), IEEE Computer Society, pp. 261–266.  doi:10.1109/MIPRO.2016.7522149.
  * [MNP01] Meyer A., Neyret F., Poulin P.:  Interactive rendering of trees with shading and shadows.  In _Eurographics Workshop on Rendering Techniques 2001_ (Vienna, 2001), Gortler S. J., Myszkowski K., (Eds.), EUROGRAPH, Springer, pp. 183–196.  doi:10.1007/978-3-7091-6242-2_17.
  * [MO95] Max N., Ohsaki K.:  Rendering trees from precomputed z-buffer views.  In _Proceedings of the Eurographics Workshop on Rendering Techniques ’95_ (Vienna, June 1995), Hanrahan P. M., Purgathofer W., (Eds.), Eurographics, Springer, pp. 74–81.  doi:10.1007/978-3-7091-9430-0_8.
  * [Mor16] Morin R.:  Inside the artificial universe that creates itself.  _The Atlantic_ (Feb. 2016).  URL: https://www.theatlantic.com/technology/archive/2016/02/artificial-universe-no-mans-sky/463308/.
  * [MP96] Mech R., Prusinkiewicz P.:  Visual models of plants interacting with their environment.  In _Proceedings of the 23rd annual conference on Computer graphics and interactive techniques_ (New York, NY, USA, Aug. 1996), SIGGRAPH ’96, Association for Computing Machinery, pp. 397–410.  doi:10.1145/237170.237279.
  * [New22] Newlands C.:  ForestGenerator GitHub repository, 2022.  URL: https://github.com/callumnewlands/ForestGenerator.
  * [NFD07] Neubert B., Franken T., Deussen O.:  Approximate image-based tree-modeling using particle flows.  _ACM Transactions on Graphics 26_ , 3 (July 2007), 88.  doi:10.1145/1239451.1239539.
  * [NM19] Nuić H., Mihajlović Ž.:  Algorithms for procedural generation and display of trees.  In _2019 42nd International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO)_ (May 2019), pp. 230–235.  doi:10.23919/MIPRO.2019.8757140.
  * [OBRvdK17] Onrust B., Bidarra R., Rooseboom R., van de Koppel J.:  Ecologically sound procedural generation of natural environments.  _International Journal of Computer Games Technology 2017_ (May 2017), 7057141.  doi:10.1155/2017/7057141.
  * [Ols04] Olsen J.:  _Realtime Procedural Terrain Generation_.  Tech. rep., IMADA, University of Southern Denmark, Odense, DNK, Oct. 2004\.
  * [OSKL05] Ong T. J., Saunders R., Keyser J., Leggett J. J.:  Terrain generation using genetic algorithms.  In _Proceedings of the 7th annual conference on Genetic and evolutionary computation_ (Washington DC, USA, June 2005), GECCO ’05, Association for Computing Machinery, pp. 1463–1470.  doi:10.1145/1068009.1068241.
  * [Pat18] Patel A.:  Using photogrammetry and megascans to plant a forest, 2018.  URL: https://www.exp-points.com/photogrammetry-tree-trunk.
  * [Per02] Perlin K.:  Chapter 2: Noise hardware.  In _Real-Time Shading Languages - ACM SIGGRAPH 2002 Course 36 Notes_ (2002), no. 36 in SIGGRAPH Course Notes, Association for Computing Machinery, pp. 2:1–2:24.
  * [PHHM96] Prusinkiewiczy P., Hananz P., Hammely M., Mech R.:  L-systems: From the theory to visual models of plants.  In _Proceedings of 2nd CSIRO Symposium on Computational Challenges in Life Sciences_ (1996), CSIRO Publishing.
  * [PHM93] Prusinkiewicz P., Hammel M. S., Mjolsness E.:  Animation of plant development.  In _Proceedings of the 20th annual conference on Computer graphics and interactive techniques_ (Anaheim, CA, Sept. 1993), SIGGRAPH ’93, Association for Computing Machinery, pp. 351–360.  doi:10.1145/166117.166161.
  * [Pho75] Phong B. T.:  Illumination for computer generated pictures.  _Communications of the ACM 18_ , 6 (June 1975), 311–317.  doi:10.1145/360825.360839.
  * [PJM94] Prusinkiewicz P., James M., Měch R.:  Synthetic topiary.  In _Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, 1994), SIGGRAPH ’94, ACM Press, pp. 351–358.  doi:10.1145/192161.192254.
  * [PL90] Prusinkiewicz P., Lindenmayer A.:  _The Algorithmic Beauty of Plants_.  Springer-Verlag, Berlin, Heidelberg, 1990.
  * [Pru86] Prusinkiewicz P.:  Graphical applications of L-systems.  In _Proceedings on Graphics Interface ’86/Vision Interface ’86_ (Toronto, CAN, 1986), Canadian Information Processing Society, pp. 247–253.
  * [QNTN03] Qin X., Nakamae E., Tadamura K., Nagai Y.:  Fast photo-realistic rendering of trees in daylight.  _Computer Graphics Forum 22_ , 3 (Sept. 2003), 243–252.  doi:https://doi.org/10.1111/1467-8659.00671.
  * [RB85] Reeves W. T., Blau R.:  Approximate and probabilistic algorithms for shading and rendering structured particle systems.  In _Proceedings of the 12th Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, July 1985), vol. 19 of _SIGGRAPH ’85_ , Association for Computing Machinery, pp. 313–322.  doi:10.1145/325334.325250.
  * [RCB∗02] Remolar I., Chover M., Belmonte O., Ribelles J., Rebollo C.:  Geometric simplification of foliage.  In _23rd Annual Conference of the European Association for Computer Graphics, Eurographics 2002 - Short Presentations_ (2002), Navazo I., Slusallek P., (Eds.), Eurographics Association.  doi:10.2312/egs.20021034.
  * [RLP07] Runions A., Lane B., Prusinkiewicz P.:  Modeling trees with a space colonization algorithm.  In _Proceedings of the Third Eurographics conference on Natural Phenomena_ (Prague, Czech Republic, Sept. 2007), NPH’07, Eurographics Association, pp. 63–70.  doi:10.5555/2381384.2381395.
  * [RSC87] Reeves W. T., Salesin D. H., Cook R. L.:  Rendering antialiased shadows with depth maps.  In _Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, July 1987), vol. 21 of _SIGGRAPH ’87_ , Association for Computing Machinery, pp. 283–291.  doi:10.1145/37401.37435.
  * [RT20] Risi S., Togelius J.:  Increasing generality in machine learning through procedural content generation.  _Nature Machine Intelligence 2_ , 8 (Aug. 2020), 428–436.  doi:10.1038/s42256-020-0208-z.
  * [Sch97] Schaufler G.:  Nailboards: A rendering primitive for image caching in dynamic scenes.  In _Proceedings of the Eurographics Workshop on Rendering Techniques ’97_ (Berlin, Heidelberg, June 1997), Springer-Verlag, pp. 151–162.
  * [SD01] Stamminger M., Drettakis G.:  Interactive sampling and rendering for complex and procedural geometry.  In _Proceedings of the 12th Eurographics Conference on Rendering_ (London, UK, June 2001), EGWR’01, Eurographics Association, pp. 151–162.
  * [SGdK∗11] Smelik R., Galka K., de Kraker K. J., Kuijper F., Bidarra R.:  Semantic constraints for procedural generation of virtual worlds.  In _Proceedings of the 2nd International Workshop on Procedural Content Generation in Games_ (New York, NY, USA, 2011), PCGames ’11, Association for Computing Machinery.  doi:10.1145/2000919.2000928.
  * [SGHS98] Shade J., Gortler S., He L.-w., Szeliski R.:  Layered depth images.  In _Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, July 1998), SIGGRAPH ’98, Association for Computing Machinery, pp. 231–242.  doi:10.1145/280814.280882.
  * [SK03] Szijártó G., Koloszár J.:  Hardware accelerated rendering of foliage for real-time applications.  In _Proceedings of the 19th Spring Conference on Computer Graphics_ (Budmerice, Slovakia, Apr. 2003), SCCG ’03, Association for Computing Machinery, pp. 141–148.  doi:10.1145/984952.984976.
  * [SKG∗09] Smelik R. M., Kraker K. J. d., Groenewegen S. A., Tutenel T., Bidarra R.:  A survey of procedural methods for terrain modelling.  In _Proceedings of the CASA Workshop on 3D Advanced Media in Gaming and Simulation (3AMIGAS)_ (Amsterdam, The Netherlands, June 2009), pp. 25–34.  URL: http://graphics.tudelft.nl/Publications-new/2009/SDGTB09a.
  * [SNK∗20] Song Y., Naji S., Kaufmann E., Loquercio A., Scaramuzza D.:  Flightmare: A flexible quadrotor simulator.  In _4th Conference on Robot Learning (CoRL 2020)_ (Cambridge MA, USA, Nov. 2020), vol. 4.
  * [TGZ20] Tarapore D., Groß R., Zauner K.-P.:  Sparse robot swarms: Moving swarms to real-world applications.  _Frontiers in Robotics and AI 7_ (2020).  doi:10.3389/frobt.2020.00083.
  * [Vir20] Virtueller-Wald.de:  The virtual forest, 2020.  URL: http://www.virtueller-wald.de/en/the-virtual-forest/.
  * [WFP∗01] Wand M., Fischer M., Peter I., Meyer auf der Heide F., Straßer W.:  The randomized z-buffer algorithm: interactive rendering of highly complex scenes.  In _Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, Aug. 2001), SIGGRAPH ’01, Association for Computing Machinery, pp. 361–370.  doi:10.1145/383259.383299.
  * [Wha05] Whatley D.:  _GPU Gems 2_.  No. 2 in GPU Gems. Addison-Wesley Professional, 2005, ch. Chapter 4\. Toward Photorealism in Virtual Botany.
  * [WHDS13] Weier M., Hinkenjann A., Demme G., Slusallek P.:  Generating and rendering large scale tiled plant populations.  _Journal of Virtual Reality and Broadcasting 10_ , 1 (2013).  doi:10.20385/1860-2037/10.2013.1.
  * [WP95] Weber J., Penn J.:  Creation and rendering of realistic trees.  In _Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques_ (New York, NY, USA, Sept. 1995), SIGGRAPH ’95, Association for Computing Machinery, pp. 119–128.  doi:10.1145/218380.218427.
  * [WS81] Witten T. A., Sander L. M.:  Diffusion-limited aggregation, a kinetic critical phenomenon.  _Physical Review Letters 47_ , 19 (Nov. 1981), 1400–1403.  doi:10.1103/PhysRevLett.47.1400.
  * [WWD∗05] Wang L., Wang W., Dorsey J., Yang X., Guo B., Shum H.-Y.:  Real-time rendering of plant leaves.  In _ACM SIGGRAPH 2005 Papers_ (Los Angeles, California, July 2005), SIGGRAPH ’05, Association for Computing Machinery, pp. 712–719.  doi:10.1145/1186822.1073252.
  * [XYS∗16] Xie K., Yan F., Sharf A., Deussen O., Huang H., Chen B.:  Tree modeling with real tree-parts examples.  _IEEE Transactions on Visualization and Computer Graphics 22_ , 12 (2016), 2608–2618.  doi:10.1109/TVCG.2015.2513409.
  * [XYY15] Xu C., Yang G., Yang M.:  Sample-based vegetation distribution information synthesis.  _PLoS ONE 10_ , 8 (Aug. 2015), e0134009.  doi:10.1371/journal.pone.0134009.


11institutetext: $\dagger$Department of Physics, Colorado School of Mines
Golden, CO, 80401 USA
11email<EMAIL_ADDRESS>

$\ddagger$Nuclear and Chemical Sciences Division, Lawrence Livermore National
Laboratory
Livermore, CA 94550, USA
11email<EMAIL_ADDRESS>
# The BeEST Experiment: Searching for Beyond Standard Model Neutrinos using
7Be Decay in STJs

K.G. Leach† and S. Friedrich‡
for the BeEST Collaboration

###### Abstract

Precision measurements of nuclear $\beta$ decay are among the most sensitive
methods to probe beyond Standard Model (BSM) physics in the neutrino sector.
In particular, momentum conservation between the emitted decay products in the
final state is sensitive to any massive new physics that couples to the
neutrino mass. One way to observe these momentum recoil effects experimentally
is through high-precision measurements of electron-capture (EC) nuclear decay,
where the final state only contains the neutrino and a recoiling atom. The
BeEST experiment precisely measures the eV-scale radiation that follows the
radioactive decay of 7Be implanted into sensitive superconducting tunnel
junction (STJ) quantum sensors. STJs are ideally suited for measurements of
this type due to their high resolution at the low recoil energies in EC decay,
and their high-rate counting capabilities. We present the motivation for the
BeEST experiment and describe the various phases of the project.

###### Keywords:

Superconducting Tunnel Junction, Neutrino Physics, Nuclear Decay

††journal: Journal of Low Temperature Physics

##  1 Introduction

The lepton sector of the Standard Model (SM) provides perhaps the best window
into Beyond Standard Model (BSM) physics given the confirmed observation of
non-zero neutrino masses 1, 2, the recent confirmation of the muon $g-2$
anomaly 3, and several other outstanding claims of leptonic BSM physics 4. As
a result, extensions to the SM description of leptons are unavoidable and must
at least account for the fact that neutrinos have non-zero mass eigenstates.
So-called “sterile neutrinos” are are present in almost all SM extensions that
include neutrino mass, and have the characteristic property that they are non-
interacting with respect to the SM forces. These sterile neutrinos are
associated with heavy mass states that do not depend on the Higgs mechanism
and can therefore be on nearly any mass scale. keV-scale sterile neutrinos are
among the most highly motivated, since they have the right cosmological
properties to explain the observed dark matter in our universe 5, 6. In fact,
extensions that include neutrinos in this mass range, such as the $\nu$MSM,
are able to provide a framework that can simultaneously explain the
observational data for dark matter and baryon asymmetry of the Universe 7, 8.

Effective experimental searches for these particles should be model-
independent and cover a large area of the allowed parameter space. The most
powerful approach to search for these heavy neutral leptons (HNLs) in the sub-
MeV mass range is through energy and momentum conservation in weak nuclear
decay. The pure electron capture (EC) decaying nucleus 7Be is the ideal case
for neutrino studies using this approach due to its pure two-body final state,
its relatively large decay energy, and simple atomic and nuclear structure.
Since the quantum information of the $\nu_{e}$ and recoiling 7Li daughter atom
are entangled following EC decay, the superposition of mass eigenstates that
make up the $\nu_{e}$ flavor state can be probed through precision
measurements of the atomic recoil momentum. All possible neutrino mass
eigenstates (including any outside of $\nu_{1,2,3}$) continue to exist as a
superposition until the wavefunction is collapsed through detection of the
7Li. When the wavefunction collapses, a single mass state is selected and the
observed kinetic energy of the recoiling atom differs depending on the mass of
any heavy particle that is emitted. The fraction of events that are shifted
with respect to the SM signal determines the new physics mixing fraction with
$\nu_{e}$ (Fig. 1). This method of decay momentum reconstruction is a simple,
model-independent approach to HNL searches since it relies only on the
existence of a heavy neutrino admixture to the active neutrinos and not on the
model-dependent details of their interactions.

Figure 1: A schematic of the nuclear EC decay of beryllium-7 (left)
illustrating the superposition of neutrino mass states ($\nu_{1,2,3}$) that
exist for the electron neutrino ($\nu_{e}$) – and the potential for heavy BSM
mass states ($\nu_{??}$). The measured lithium-7 recoil energies measured in
the experiment differ depending on the mass-state of the emitted neutrino, and
is schematically shown in the bottom right.

##  2 Nuclear Electron Capture Decay of 7Be

7Be decays by EC primarily to the nuclear ground-state of 7Li with a $Q$ value
of 861.89(7) keV 9 and a half-life of $T_{1/2}=53.22(6)$ days 10. A small
branch of 10.44(4)% results in the population of a short-lived excited nuclear
state in 7Li ($T_{1/2}=72.8(20)$ fs) 10 that de-excites via emission of a
477.603(2) keV $\gamma$-ray 11. In the EC process, the electron can be
captured either from the $1s$ shell ($K$-capture) or the $2s$ shell of Be
($L$-capture). For $K$-capture, the binding energy of the $1s$ hole is
subsequently liberated by emission of an Auger electron whose energy adds to
the decay signal and separates it from the $L$-capture signal. Since the
nuclear decay and subsequent atomic relaxation occur on short time scales, a
direct measurement produces a spectrum with four peaks: two for $K$-capture
and two for $L$-capture into the ground state and the excited state of 7Li,
respectively. Due to the relative spatial overlaps of the $1s$ and $2s$
electron orbitals with the nucleus, the EC process is dominated by $K$-shell
capture with a measured $L/K$ capture ratio of 0.040(6) in HgTe 12 and
0.070(7) in Ta 13. The simplicity of the atomic structure in particular is one
of the primary reasons for using 7Be to perform these BSM physics searches.

##  3 Superconducting Tunnel Junctions (STJs)

Figure 2: A schematic of the BeEST experimental concept.

The detection and characterization of $\sim 20-120$ eV recoil energies from
7Be EC decay at the precision required to provide competitive limits on BSM
physics are difficult to achieve. The BeEST confronts these challenges using
high-rate STJ radiation detectors and is the first instance of adapting them
to search for BSM physics in nuclear decay. STJs are high-rate quantum sensors
that were originally developed for high-resolution X-ray spectroscopy 14, 15,
16 and have been used for many years in astronomy 16 and at synchrotron
facilities 17. STJs are a type of Josephson junction that consists of two
superconducting electrodes separated by a thin insulating tunneling barrier.
The absorption of radiation in one of the electrodes breaks the Cooper pairs
of the superconducting ground state and excites free excess charge carriers
above the superconducting energy gap $\Delta$ in proportion to the absorbed
energy. As these excess charges tunnel across the barrier, they generate a
temporary increase in current that can be directly read out with a field-
effect transistor based preamplifier at room temperature. The high energy
resolution in STJs is due to the fact that the energy to create an excess
charge $\epsilon=1.7\Delta$ scales with the energy gap $\Delta$ 14. For
superconductors, $\Delta$ is of order 1 meV and thus roughly three orders of
magnitude smaller than the band-gap in semiconductors, giving rise to $\sim 1$
eV resolution at the signal energies relevant to the BeEST 18, 19. In order to
prevent thermal excitations across these very small gaps, STJs must be
operated well below the superconducting transition temperature $T_{C}$.

The primary advantages of using STJs for this type of experiment are two-fold:

  1. 1.

The maximum STJ count rates are determined by the time that the excess charges
remain excited above $\Delta$ before they recombine into the superconducting
ground state and again form Cooper pairs, which is typically in the
$10~{}\mu$s range 20. This enables each STJ detector pixel to operate at rates
up to $10^{4}$ counts/s 15, which places them among the highest-rate LTD
technologies with high energy resolution, and makes them ideal to search for
heavy neutrinos with very weak couplings.

  2. 2.

The high energy resolution of STJs can be well characterized for signal
energies in the range of interest for the BeEST with high-rate pulsed lasers
19. Their response is well understood, very predictable, and can be measured
in-situ 18, 19, 21.

##  4 Harvesting High-Intensity 7Be+ Beams

For the BeEST experiment, a pure, high-intensity radioactive 7Be source is
directly implanted into the STJs to ensure that all of the decay products
deposit their energy in the devices. This prevents partial energies from being
recorded due to surface effects where the Auger electrons could leave the
device with some kinetic energy, thus broadening the observed recoil signals
at low energies. The required 7Be is harvested at the TRIUMF Isotope Separator
and Accelerator (ISAC) rare-isotope beam facility on the campus of the
University of British Columbia in Vancouver, Canada 22. TRIUMF-ISAC is able to
produce among the purest ($>99.99\%$), high-intensity ($\geq 10^{8}$ s-1) beam
of 7Be in the world via the isotope separation on-line (ISOL) technique 23. To
perform the implantation, the low-energy ISAC implantation station (IIS) is
used, where laser-purified radioactive ion beams are delivered at energies up
to 50 keV for implantation. The ISAC facility is currently undergoing target
development towards increased yields of 7Be using a new state-of-the-art
graphite ISOL target. Since the 7Be is produced mainly by spallation on
carbon, the prospect of using pure graphite foils is attractive since they are
exempt from the same license restrictions of running proton currents on
uranium-carbide. Thus with a pure graphite target and a proton current of
$100~{}\mu$A, 7Be intensities of $\geq 10^{10}$ s-1 are projected, making the
beam-time requirements for Phase-IV of the BeEST less demanding.

##  5 The BeEST Experiment

Figure 3: Images of the detectors for the various phases of the BeEST
experiment. From left to right - a 10-pixel Ta-based chip used for Phases-I
and -II, a 36-pixel Ta array used for Phase-III, and a 128-pixel Al-based STJ
array prototype for Phase-IV.

The BeEST experiment employs a staged approach that employs four phases to
demonstrate the concept, evaluate systematics of the experiment, scaling the
experiment to multi-pixel arrays, and operation of STJ arrays with new
superconducting material. Phases-I and -II have already been successfully
completed 13, 19, 21, and Phases-III and -IV are currently underway. Briefly
described below are the main objectives and milestones for the 4 phases of the
BeEST experiment.

###  5.1 Phase-I: Proof-of-Concept

The primary goal of Phase-I was to demonstrate the concept of the BeEST and
determine if STJ detectors could be used to perform high-rate nuclear recoil
measurements of EC decay. The Ta-based STJs used were five-layer devices
consisting of Ta (165 nm) - Al (50 nm) - Al2O3 (1 nm) - Al (50 nm) - Ta (265
nm) that were fabricated by photolithography at STAR Cryoelectronics 24 (Fig.
3 (left)). These detectors were used previously for low-energy decay
measurements 25 and were exceptionally well characterized prior to 7Be
implantation 18. The first $A=7$ irradiation at TRIUMF took place over a 16
hour period where a total of $2\times 10^{12}$ 7Be+ and 7Li+ ions were
implanted into 6 STJ detectors at an energy of 26 keV. Of these ions,
approximately $8\times 10^{10}$ were implanted into the $(68~{}\mu$m)2 STJs of
interest for Phase-I - requiring each pixel to count at roughly 800 counts/s.
The results of the high-rate Phase-I data demonstrated the concept, and
provided the platform to initiate Phase-II of the BeEST. However, given the
large number of 7Li atoms that were present in the beam, damage effects were
observed in the STJ that required mitigation in Phase-II.

###  5.2 Phase-II: First BSM Physics Search and Improved Detector
Characterization

The primary goal of Phase-II was to provide the first exclusion limits in our
search for BSM physics in the neutrino sector. To achieve this, we required
the ability to i) evaluate possible ion-beam damage effects suffered in
Phase-I, ii) better evaluate the detector response to nuclear recoils, and
iii) provide a high-precision energy calibration of the recoil spectrum. The
implantation time was limited to 30 minutes to reduce damage effects, which
resulted in $2\times 10^{8}$ 7Be+ ions implanted into the (138 $\mu$m)2 STJ
used for the characterization tests, resulting in an initial event rate per
pixel of only $\sim 10$ counts/s. Data were acquired for $\sim$20 hours/day
over a period of four weeks in the LLNL ADR. Since we implanted $100\times$
less beam in the second implantation, the (138 $\mu$m)2 STJs were used, which
gave the best combination of statistics and resolution. For energy calibration
of the Phase-II data, the STJs were simultaneously exposed to 3.49965(15) eV
photons from a pulsed Nd:YVO4 laser triggered at a rate of 100 Hz. We adjusted
the laser intensity such that multi-photon absorption provides a comb of peaks
over the energy range of interest from 20 - 120 eV in the (138 $\mu$m)2 STJ.
The energy calibration of the spectrum was statistics limited, and resulted in
a precision of 0.005 eV - thus providing a successful demonstration of the
calibration method for future phases of the experiment. This high-level of
detector characterization 19 allowed for the extraction of exclusion limits
from the Phase-II data that were up-to an order of magnitude better than from
all previous laboratory-based limits in the 100 - 860 keV mass range 21.

###  5.3 Phase-III: Scaling of the BeEST to 36- and 112-Pixel Arrays

To provide significant increases in the statistical precision of the BeEST
experiment, large arrays of STJs are required. For Phase-III we are currently
in the process of evaluating scaling the experiment to 36- and 112-pixel
detector arrays of $(200~{}\mu$m)2 Ta-based STJs that were developed for
synchrotron measurements (Fig. 3 (middle)). This will provide two orders of
magnitude increase in statistics over the Phase-II limits for the same running
period. More importantly, these measurements will serve as a test ground to
understand any scaling-related limitations before Phase-IV.

###  5.4 Phase-IV: Final Exclusion Limits with High-Resolution Al-Based STJ
Arrays

For Phase-IV, we have started to fabricate and operate arrays of pure Al (200
nm) - Al2O3 (1 nm) - Al (200 nm) STJs. Al-based STJs have the advantage that
the energy gap $\Delta\sim 0.17$ meV is $4\times$ smaller than the energy gap
in Ta, potentially providing an intrinsic energy resolution of $\sim 0.5$ eV.
Three arrays will be used with pixel sizes of $(70\mu$m)2, $(130\mu$m)2, and
$(200\mu$m)2, respectively, initially allowing us to better understand any in-
medium effects as well as compare data across new detectors. In our current
design, the new STJs will be deposited onto $\sim 1~{}\mu$m thin Si3N4
membranes with the area between pixels covered by a layer of Au to thermalize
high-energy substrate phonons to minimize the number of excess charges in the
STJs resulting from this background process. The final limits from the BeEST
experiment will result from a 128-pixel Al-STJ array counting at 1 000
counts/s/pixel for 100 days.

##  6 Conclusions

The BeEST experiment is a sensitive, model-independent search for any BSM
physics that couples to the neutrino mass using the EC decay of 7Be implanted
in STJs. The high-rate ability of STJs allows for the collection of high
statistics in reasonable experimental periods, and takes advantage of the
high-dose of 7Be that is achievable at TRIUMF-ISAC. The BeEST experiment has
completed the first two phases of the project - including setting the first
model-independent limits on the existence of HNLs in the 100 - 860 keV mass
range from EC decay - and is currently running Phase-III tests. The
demonstration of low-energy nuclear recoil measurements in the first two
phases of the BeEST experiment may also initiate new experiments with other
rare isotopes in our continued search for BSM physics from nuclear decay.

###### Acknowledgements.

The BeEST is supported by the U.S. Department of Energy, Office of Science
under grant No. DE-SC0021245, and the LLNL Laboratory Directed Research and
Development program through grants No. 19-FS-027 and No. 20-LW-006. TRIUMF
receives federal funding via a contribution agreement with the National
Research Council of Canada. This work is performed under the auspices of the
U.S. Department of Energy by Lawrence Livermore National Laboratory under
Contract No. DE-AC52-07NA27344. We thank the ongoing efforts of the entire
BeEST collaboration and acknowledge the extensive efforts that everyone has
put into our program over the past several years!

## References

  * Ahmad  _et al._ 2001 Q. R. Ahmad _et al._ (SNO Collaboration), Phys. Rev. Lett. 87, 071301 (2001).
  * Fukuda  _et al._ 1998 Y. Fukuda _et al._ (Super-Kamiokande Collaboration), Phys. Rev. Lett. 81, 1562 (1998).
  * Abi  _et al._ 2021 B. Abi _et al._ (Muon $g-2$ Collaboration), Phys. Rev. Lett. 126, 141801 (2021).
  * Garisto 2020 D. Garisto, Physics 13, 79 (2020).
  * Boyarsky  _et al._ 2019 A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, Progress in Particle and Nuclear Physics 104, 1 (2019).
  * Adhikari  _et al._ 2017 R. Adhikari _et al._ , Journal of Cosmology and Astroparticle Physics 2017, 025 (2017).
  * Asaka  _et al._ 2005 T. Asaka, S. Blanchet, and M. Shaposhnikov, Physics Letters B 631, 151 (2005).
  * Asaka and Shaposhnikov 2005 T. Asaka and M. Shaposhnikov, Physics Letters B 620, 17 (2005).
  * Wang  _et al._ 2017 M. Wang, G. Audi, F. Kondev, W. Huang, S. Naimi, and X. Xu, Chinese Physics C 41, 030003 (2017).
  * Tilley  _et al._ 2002 D. Tilley, C. Cheves, J. Godwin, G. Hale, H. Hofmann, J. Kelley, C. Sheu, and H. Weller, Nuclear Physics A 708, 3 (2002).
  * Helmer and van der Leun 2000 R. Helmer and C. van der Leun, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 450, 35 (2000).
  * Voytas  _et al._ 2001 P. A. Voytas, C. Ternovan, M. Galeazzi, D. McCammon, J. J. Kolata, P. Santi, D. Peterson, V. Guimarães, F. D. Becchetti, M. Y. Lee, T. W. O’Donnell, D. A. Roberts, and S. Shaheen, Phys. Rev. Lett. 88, 012501 (2001).
  * Fretwell  _et al._ 2020 S. Fretwell, K. G. Leach, C. Bray, G. B. Kim, J. Dilling, A. Lennarz, X. Mougeot, F. Ponce, C. Ruiz, J. Stackhouse, and S. Friedrich, Phys. Rev. Lett. 125, 032701 (2020).
  * Kurakado 1982 M. Kurakado, Nuclear Instruments and Methods in Physics Research 196, 275 (1982).
  * Friedrich  _et al._ 2003 S. Friedrich, A. Vailionis, O. Drury, T. Niedermayr, T. Funk, W. N. Kang, Eun-Mi Choi, Hyeong-Jin Kim, Sung-Ik Lee, S. P. Cramer, Changyoung Kim, and S. E. Labov, IEEE Transactions on Applied Superconductivity 13, 1114 (2003).
  * Verhoeve 2008 P. Verhoeve, Journal of Low Temperature Physics 151, 675 (2008).
  * Carpenter  _et al._ 2014 M. H. Carpenter, S. Friedrich, J. A. Hall, J. Harris, and R. Cantor, Journal of Low Temperature Physics 176, 222 (2014).
  * Ponce  _et al._ 2016 F. Ponce, M. H. Carpenter, R. Cantor, and S. Friedrich, Journal of Low Temperature Physics 184, 694 (2016).
  * Friedrich  _et al._ 2020 S. Friedrich, F. Ponce, J. A. Hall, and R. Cantor, Journal of Low Temperature Physics 200, 200 (2020).
  * Friedrich  _et al._ 2009 S. Friedrich, T. Hertrich, O. B. Drury, N. J. Cherepy, and J. Hohne, IEEE Transactions on Nuclear Science 56, 1089 (2009).
  * Friedrich  _et al._ 2021 S. Friedrich, G. B. Kim, C. Bray, R. Cantor, J. Dilling, S. Fretwell, J. A. Hall, A. Lennarz, V. Lordi, P. Machule, D. McKeen, X. Mougeot, F. Ponce, C. Ruiz, A. Samanta, W. K. Warburton, and K. G. Leach, Phys. Rev. Lett. 126, 021803 (2021).
  * Dilling and Krücken 2014 J. Dilling and R. Krücken, Hyperfine Interactions 225, 111 (2014).
  * Blumenfeld  _et al._ 2013 Y. Blumenfeld, T. Nilsson, and P. V. Duppen, Physica Scripta 2013, 014023 (2013).
  * Carpenter  _et al._ 2013 M. H. Carpenter, S. Friedrich, J. A. Hall, J. Harris, W. K. Warburton, and R. Cantor, IEEE Transactions on Applied Superconductivity 23, 2400504 (2013).
  * Ponce  _et al._ 2018 F. Ponce, E. Swanberg, J. Burke, R. Henderson, and S. Friedrich, Phys. Rev. C 97, 054310 (2018).