# DECal, a Python tool for the efficiency calculation of thermal neutron
detectors based on thin-film converters

Álvaro Carmona Basáñez  Kalliopi Kanaki  Francesco Piscitelli  European
Spallation Source ERIC, P.O. Box 176, SE-221 00 Lund, Sweden

###### Abstract

The Detector Efficiency Calculator (DECal) is a series of Python functions and
tools designed to analytically calculate, visualise and optimise the detection
efficiency of thermal neutron detectors, which are based on thin-film
converters. The implementation presented in this article concerns 10B-based
detectors in particular. The code can be run via a graphical user interface,
as well as via the command line. The source code is openly available to
interested users via a GitHub repository.

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: DECal
Licensing provisions: Free for non-commercial use (as per terms in LICENSE
file)
Programming language: Python Version 2.7 and 3.3
Operating system: OS X (preferred) and Linux
Version: 1.0.0

Keywords:
Detector efficiency calculator, neutron scattering, Boron-10, neutron
detector, efficiency, thin film, converter, Python, PyQt5,

Nature of problem:

Implementation of a Python-based tool to calculate and optimise detection
efficiency for thin-film thermal neutron detectors via a graphical user
interface and a command line interface.

##  1 Introduction

The European Spallation Source (ESS) ERIC [1] is a joint European organisation
committed to the construction and operation of the world’s leading facility
for research using thermal neutrons. The ESS is designed to be the world’s
brightest neutron source and the instantaneous neutron flux on detectors at
ESS full power will be without precedent. In general, neutron scattering
facilities are going toward higher fluxes and this translates into higher
demands on instrument and detector performance: higher counting rate
capability, better timing and finer spatial resolution are requested among
others.

The current thermal neutron detector technology is reaching fundamental limits
and most of the neutron sources in the world, including ESS, are pushing the
development of their detector technologies [2, 3, 4, 5], to tackle the
increased flux available and the scarcity of $\mathrm{{}^{3}He}$, the so-
called “Helium-3 crisis” [6, 7, 8]. $\mathrm{{}^{10}B}$ along with the
$\mathrm{{}^{3}He}$ and $\mathrm{{}^{6}Li}$ isotopes are the main actors in
thermal neutron detection due to their large neutron absorption cross sections
[9]. Due to its favourable properties, $\mathrm{{}^{3}He}$ (a rare isotope of
$\mathrm{He}$) has been dominating thermal neutron detection for decades.

Nowadays, the importance of solid conversion layers is increasing as this
technology appears to be a viable and promising alternative to
$\mathrm{{}^{3}He}$. Recently, high-quality, low-cost production of square
metres of $\mathrm{{}^{10}B_{4}C}$-coated substrates [10, 11, 12] became
possible. The detection efficiency of a single thin layer of
$\mathrm{{}^{10}B_{4}C}$ is limited to a few percent at thermal neutron
energies compared to the very high efficiencies of $\mathrm{{}^{3}He}$-based
detectors. The Multi-Blade [13, 14, 15], the Jalousie detector [16, 17], the
A1CLD [18, 19] and many others [20, 21], are examples of the detector
developments which exploit solid neutron converters operated at a grazing
angle in order to increase the detection efficiency. A different way to
increase the efficiency is to arrange many layers in sequence. Examples of
detector developments arranging several layers are the Multi-Grid [22, 23, 24,
25, 26], CASCADE [27, 28] and many others [29, 30, 31, 32, 33, 34, 35, 36,
37].

The difference between the physical processes of a gas-based (such as
$\mathrm{{}^{3}He}$) and a solid-converter-based detector (such as
$\mathrm{{}^{10}B_{4}C}$) is described in detail in [38, 39]; a detector based
on solid converter has several parameters that can be tuned to increase the
detection efficiency. A simple model of single and multi-layer thermal neutron
detectors is explored mainly analytically to help optimize the design in
different circumstances. Several theorems are deduced that help guide the
design. Using powerful simulation software has the advantage of including many
effects and potentially results with high accuracy [40]. On the other hand it
does not always give the insight an equation can deliver. The equations in
[39, 41] and the tool described in this manuscript focus on
$\mathrm{{}^{10}B_{4}C}$-based thermal neutron detectors but can be extended
to any detector based on thin-film converters.

The theoretical details of the calculations which this software implements in
Python are described in [39, 41] and more specifically in the Ph.D. thesis
[42]. The authors explain a series of equations to calculate the efficiency of
detectors based on thin film neutron converters given the specific geometry of
the detector and the neutron beam characteristics.

##  2 Detector geometry and neutron beam configuration

A detector can be made of a single or multiple layers. A layer consists of a
substrate material, usually Aluminium, on which the 10B4C converter is coated,
either on one or both sides. A double-coated layer is referred to as blade.
The number of converter layers, their thickness and their composition matter.
In addition, a detector can contain layers of the same or varying converter
thickness. This parameter can be optimised accordingly for a single neutron
wavelength or for a distribution of neutron wavelengths, a scenario which is
closer to reality in a neuron scattering instrument. The material and
thickness of the substrate are not considered in the calculations presented
here and will be the topic of a future improvement.

The incoming neutrons can enter this geometrical arrangement at an angle or
perpendicularly, depending on the needs of the application (see Fig. 1). Both
the neutron incident angle and its wavelength (energy) affect the efficiency
and thus enter the respective calculation function.

Fig. 1 summarizes all parameters that impact the detection efficiency. The
neutron beam hits the detector with an incident angle $\Theta$. A numbered
series of layers follow. A layer is called back-scattering when neutrons are
incident from the gas-converter interface and the escaping particles are
emitted backwards into the gas volume; it is called a transmission layer when
neutrons are incident from the substrate-converter interface and the escaping
fragments are emitted in the forward direction in the sensitive volume. In
Fig. 1 the detection efficiency of every blade is the sum of the back-
scattering and transmission layer efficiencies, as the substrate holds two
converter layers, one in back-scattering mode and one in transmission mode.

Figure 1: Depiction of detector geometry and neutron beam arrangement. The
beam hits the parallel blades with an incident angle $\Theta$. The blades
consist of a substrate on which the converter is coated on both sides.

##  3 Software overview

DECal has been developed as a cross-platform software for Mac OS X (version
10.11 and higher) and Linux (tested on CentOS 7 and Ubuntu 17.4). It has been
written in Python 3. The objective of this project has been to develop an open
source tool that allows the user to calculate the detection efficiency of a
10B-based multi-blade neutron detector and optimise the respective geometry
parameters.

The software features allow the user to set the parameters of the neutron beam
and detector with intuitive visualisation features. The package containing the
efficiency calculation functions is distributed via Python Package Index
(pip).

The configurations can be exported in JSON format, and this JSON file can be
imported with graphical interface as well as in the functions of the library
as input.

###  3.1 Software availability

The software is divided in two different repositories. DECal is available in
https://github.com/DetectorEfficiencyCalculator/dg_efficiencyCalculator In the
repository https://github.com/alvcarmona/neutronDetectorEffFunctions the users
can find the library that holds the calculation functions. It is also
available via easy install via the following command:

     $   pip install neutron_detector_eff_functions

###  3.2 Requirements and dependencies

This is a list of hardware and software requirements for the users to access
all the features of the tool:

  * 1.

OS X $\geq$10.10, Linux Ubuntu 17.4 and CentOS 7 are the OS where the software
has been tested

  * 2.

Multi-core CPU

  * 3.

1200$\times$610 screen resolution or higher

  * 4.

Python v2.7 or 3.3 and higher

  * 5.

Qt5 for the DECal application

The Python libraries available via pip needed to execute the DECal software
are included in the requirements.txt file found in the repository. These are
some of them:

  * 1.

SciPy [43]

  * 2.

NumPy [44]

  * 3.

Matplotlib $\geq$ 1.5.0 [45]

  * 4.

PySide $\geq$ 1.2.4

  * 5.

PyQt5 (not available via pip in some systems)

For more information about the environment needed the user can refer to the
requirements.txt file.

##  4 Software Design

The classes used to represent the Detector and its parts are the following:
Detector, Blade, B10. The Detector class is the object that represents a
detector configuration and has the plot functions. To deal with the
construction of the Detector entity with the different parameters there is a
class that is built with the Builder pattern [46]. The builder pattern builds
a complex object using simple objects and using a step by step approach. The
builder function builds a detector entity depending on the parameters given.

###  4.1 File structure

Script functions are in scripts.py. In this folder there are other basic
Python files like the requirements file and the license. The exports folder
contains some examples of exported neutron configurations and wavelengths. The
Detector class is written in Detector.py. Efftools holds the core efficiency
calculation functions and the metadata functions for plotting. The entire tool
is built around efficiency4boron and efficiency2particles functions. These
functions are the ones used to calculate the theoretical efficiency of a
neutron detector using B10 as converter.

    ============
    neutronDetectorEffFunctions library file structure
    ============

    README.md
    LICENSE.txt
    neutron_detector_eff_functions/
        Aluminium.py
        B10.py
        Blade.py
        Converter.py
        Detector.py
        Detector_meta.py
    data/
        Aluminium/
            AlCrossSect_(n,g).py
        B10/
            10B4C220/
                IONIZ_Alpha06.py
                IONIZ_Alpha94.py
                IONIZ_Li06.py
                IONIZ_Li94.py
            10B4C224/
                IONIZ_Linkoping_Alpha06.py
                IONIZ_Linkoping_Alpha94.py
                IONIZ_Linkoping_Li06.py
                IONIZ_Linkoping_Li94.py
            B10CrossSect_(n,a).py

    efftools.py
    scripts.py # examples of use
    requirements.txt
    setup.cfg
    setup.py
    tests/
        Detector_test.py
        test_B10.py

The data subpackage contains the files the application uses to calculate the
cross-section of different converter materials. In the current version of the
software there are two available configurations of B10. These files have been
generated with a tool called SRIM (The Stopping and Range of Ions in Matter)
[47]. This data is the stopping power of the materials. The user can find the
description of the materials in the header lines of the files
efficiencyCalculator/data/B10/. This is an example of this information:

    ============= TARGET MATERIAL ======================================
    Layer 1 : Layer 1
    Layer Width = 50000.E+00 A ;
      Layer # 1- Density = 12.74E22 atoms/cm3 = 2.207 g/cm3
      Layer # 1-  B = 77.6 Atomic Percent = 74.4 Mass Percent
      Layer # 1-  C = 20   Atomic Percent = 23.0 Mass Percent
      Layer # 1-  B = 2.4  Atomic Percent = 2.53 Mass Percent
    ====================================================================
     Total Ions calculated =000999.00
    ===============================================================
      Ionization Energy Units are  >>>>  eV /(Angstrom-Ion)  <<<<
    ===============================================================

Fig. 2 represents the class diagram of the DECal project.

Figure 2: Class diagram of the DECal project.

##  5 High level functions

A set of functions has been developed for simplifying the use of the classes
included in the project. The functions can be called using the  scripts.py
file located in the root folder of the library. The explanations listed in the
following serve documentation purposes.

    calculate_eff_multiblade(nb, converterThickness, substrateThickness,
     wavelength, angle, threshold, single, converter)

The function returns a list of two positions. The first contains a list of the
efficiency for each blade in depth order (the last is the deepest) and the
second value is the total efficiency of the detector. This is an example of
usage:

    calculate_eff_multiblade(10, 1, 0, [[1.8,100]], 90, 100, False,
    ’10B4C 2.24g/cm3’)

The function call returns the efficiency of a neutron detector configuration
of 10 blades with a converter thickness of 1 micrometre, 0 micrometres of
substrate, a neutron wavelength of 1.8 Å (monoenergetic), 90 degrees of
incident neutron angle, an applied threshold of 100 keV on the total energy
deposition from the ionized conversion products, double-coated layers and a
10B4C converter with a density of 2.24 g/cm3 [10].

    calculate_eff_json(path)

The function calculates the efficiency of a configuration loaded from a JSON
file. The user can generate this configuration via the graphical user
interface.

    plot_eff_vs_thick(path)

This function plots the detector efficiency as a function of different
converter thicknesses. The argument is the path to a JSON configuration file.
A typical result is presented in Fig. 3.

Figure 3: Detector efficiency vs. converter thickness for a set of default
beam and detector parameters, e.g. 1.8Å wavelength, 1.3 $\mu$m converter
thickness, read from a typical JSON file in GitHub.

    plot_eff_vs_wave(path)

Providing a path to a JSON file as parameter, this function plots the
efficiency as a function of a monochromatic neutron wavelength. The resulting
plot is displayed in Fig. 4.

Figure 4: Detector efficiency vs. neutron wavelength, produced from a typical
JSON file in GitHub.

    optimize_config_same_thick(originPath, destinyPath)

The function returns a list of converter thicknesses optimised to maximise
detector efficiency, with the condition that all converter thicknesses are
identical. originPath points to a configuration from a JSON file. destinyPath
is the destination where the new configuration is saved in JSON format.

for each blade the same converter thickness is applied on both substrate
sides.

    optimize_config_diff_thick(originPath, destinyPath)

Same as before but with the condition that different blades can have different
converter thicknesses.

##  6 Graphical user interface tool

A graphical user interface (GUI) has been developed for easy access to the
features of the software. The GUI can be executed running the launch.py
script, located in the root folder of the project. This invokes the main
window as shown in Fig. 5.

Figure 5: Main GUI window, where detector configurations are listed. New
configurations can be created from here or existing ones can be edited,
duplicated and imported.

The detector configuration list is empty to start with. To create a new
configuration click on the New configuration button and a new window pops up
(see Fig. 6). Here the user can configure the neutron detector and neutron
beam characteristics. When all the parameters are set, the user can proceed to
calculate the efficiency, get an optimised configuration, export the data from
the plots or export the configuration to a JSON file. The full list of use
cases can be found in the use cases section.

Figure 6: Detector configurator view. This is the window allowing the user to
define the parameters, display plots and detector efficiency values.

##  7 Use cases

The following list describes the different functionalities available in the
software.

Create a configuration: There are different parameters that can be changed
with the user input: name, substrate, converter, and energy threshold. Then
the user can Add Wavelength and Add Blades. The configuration can be saved for
later use in this session by clicking Save, and can be exported by clicking on
Export Configuration.

Add wavelength: The user can add wavelength ($\lambda$) in two ways, either by
manually setting a value and its respective ratio (0–100%) or via an imported
text file. The latter should have a two-column format as  in file
2gaussdistr.txt. The tool will display a list of wavelengths and weights (see
Fig. 7), and a plot (see Fig. 8).

Figure 7: Example of an imported wavelength list with respective weights.
Figure 8: Example of a plot from an imported wavelength list with a double
Gaussian shape. The red line is at the barycenter of the wavelength
distribution.

Add blades: The user can set the converter thickness and the number of
substrate layers. The GUI provides visualisation aid, i.e. a list with the
added blades and plots of the thickness vs. the blade number. The user can set
a single layer detector to separately see the backscattering and transmission
efficiency values. To change the blade configuration the Delete Blades button
has to be pushed or the desired value in the blade list can be modified (see
Fig. 9).

Calculate efficiency: When the configuration is complete, the user can
calculate the efficiency pressing the Calculate Efficiency button or calling
the function from the scripts.py file. When the efficiency is calculated, the
window displays the total efficiency for double-coated layers or both
backscattering and transmission values for single layer configurations. The
window displays the respective efficiency of a blade in a multi-blade
configuration. When a multi-blade detector has all the layers with the same
thickness, a plot of the efficiency vs. the blade number and another plot with
the total efficiency vs. converter thickness are displayed (see Fig. 3).

Optimisation: The list of converter thicknesses of detector can be optimised
for maximum efficiency. The optimisation is done by calling the Calculate
Efficiency function changing the thickness until the maximum is found. In a
multi-blade detector, which has to be optimised for any distribution of
neutron wavelengths or for a single wavelength, each blade has to hold two
converter layers of the same thickness ([39], page 11). Naturally converter
thicknesses on different blades can be distinct.

The most complex operation is the optimisation of a multi-blade detector with
different coating thicknesses and polychromatic wavelength. This operation
could take several minutes for a multi-core CPU but as it has been
demonstrated ([39], page 16), it is a sufficient approximation to optimise
using the barycenter of the wavelength distribution.

Figure 9: Example of a list of blade thicknesses optimised as a function of
depth with the aim of maximizing detector efficiency. Figure 10: Example of
blade thicknesses optimised as a function of depth with the aim of maximizing
detector efficiency.

Exporting data: The active neutron detector configuration can be exported from
the detector dialog tab. The configuration is exported in the JSON format.
This JSON file can be imported in the main window or used as a parameter for a
script function. The plots are exported with the desired format using the Save
feature in the toolbar under the plots. Data from the plot is exported with
the Export Data button under some of the plots. This is an example of a
configuration exported as a JSON file:

⬇

{

"angle": 90,

"blades": [

{

"backscatter": 4.0,

"inclination": 0,

"substrate": 0.0,

"transmission": 0

}

],

"converter": "10B4C 2.24g/cm3",

"name": "Detector 1",

"single": true,

"threshold": 100,

"wavelength": [

{

"%": 100,

"angstrom": 1.8

}

]

}

##  8 Conclusions

Based on prior analysis of the thin-film neutron detector characteristics an
open source software is developed, in order to calculate and optimise detector
efficiency. The focus of this work is on B10-based neutron detectors. The
Python code produced to this end is accessible via a GitHub repository for the
neutron scattering community to access, use and modify according to the users’
needs. To this end, the calculations are made available as a separate Python
library but can be run via a GUI application as well. Various utilities are
offered that make the software quick to start with and intuitive to use.

The software is available at
https://github.com/DetectorEfficiencyCalculator/dg_efficiencyCalculator and
https://github.com/alvcarmona/neutronDetectorEffFunctions. Part of the DECal
functionality is implemented as a web application and can be found in
https://github.com/alvcarmona/efficiencycalculatorweb.

## Acknowledgments

This work was supported by the EU Horizon 2020 framework, BrightnESS project
676548. Álvaro Carmona Basáñez would like to thank the Basque government and
Novia Salcedo for their partial financial support via the Global training
grant.

## References

  * Peggs et al. [2013] S. Peggs, et al., ESS Technical Design Report, ESS 2013-001, http://esss.se/scientific-technical-documentation, 2013\.
  * Cooper et al. [2003] R. Cooper, et al., A program for neutron detector research and development - A White Paper based on a workshop held at Oak Ridge National Laboratory, Proceedings of the workshop held at Oak Ridge National Laboratory (2003).
  * Gebauer [2004] B. Gebauer, Nuclear Instruments and Methods in Physics Research Section A 535 (2004) 65 – 78. doi:10.1016/j.nima.2004.07.266.
  * Zeitelhack [2012] K. Zeitelhack, Neutron News 23 (2012) 10–13. doi:10.1080/10448632.2012.725325.
  * Kirstein and et al. [2014] O. Kirstein, et al., arXiv:1411.6194 Proceedings of the 23rd International Workshop on Vertex Detectors, 15-19 September, Macha Lake, The Czech Republic (2014).
  * Kouzes [2009] R. T. Kouzes, The 3He Supply Problem, Technical Report, Pacific Northwest National Laboratory, Richland, WA, 2009.
  * Kramer [2011] D. Kramer, Physics Today 64 (2011) 20–23.
  * Shea and Morgan [2010] D. A. Shea, D. Morgan, Congressional Research Service (2010).
  * Dianoux and Lander [2013] A.-J. Dianoux, G. Lander, Neutron Data Booklet, Second Edition, 2013.
  * Höglund et al. [2012] C. Höglund, J. Birch, K. Andersen, T. Bigault, J.-C. Buffet, J. Correa, P. van Esch, B. Guerard, R. Hall-Wilton, J. Jensen, A. Khaplanov, F. Piscitelli, C. Vettier, W. Vollenberg, L. Hultman, Journal of Applied Physics 111 (2012). doi:10.1063/1.4718573.
  * Schmidt et al. [2016] S. Schmidt, C. Höglund, J. Jensen, L. Hultman, J. Birch, R. Hall-Wilton, Journal of Materials Science 51 (2016) 10418–10428. doi:10.1007/s10853-016-0262-4.
  * Höglund et al. [2015] C. Höglund, K. Zeitelhack, P. Kudejova, J. Jensen, G. Greczynski, J. Lu, L. Hultman, J. Birch, R. Hall-Wilton, Radiation Physics and Chemistry 113 (2015) 14–19. doi:10.1016/j.radphyschem.2015.04.006.
  * Piscitelli [2015] F. Piscitelli, The European Physical Journal Plus 130 (2015) 1–9. doi:10.1140/epjp/i2015-15027-3.
  * Piscitelli et al. [2014] F. Piscitelli, J. C. Buffet, J. F. Clergeau, S. Cuccaro, B. Guérard, A. Khaplanov, Q. L. Manna, J. M. Rigal, P. V. Esch, Journal of Instrumentation 9 (2014) P03007.
  * Piscitelli et al. [2017] F. Piscitelli, F. Messi, M. Anastasopoulos, T. Bryś, F. Chicken, E. Dian, J. Fuzi, C. Höglund, G. Kiss, J. Orban, P. Pazmandi, L. Robinson, L. Rosta, S. Schmidt, D. Varga, T. Zsiros, R. Hall-Wilton, Journal of Instrumentation 12 (2017) P03013.
  * Henske et al. [2012] M. Henske, M. Klein, M. Köhli, P. Lennert, G. Modzel, C. Schmidt, U. Schmidt, Nuclear Instruments and Methods in Physics Research Section A 686 (2012) 151 – 155. doi:10.1016/j.nima.2012.05.075.
  * Modzel et al. [2014] G. Modzel, M. Henske, A. Houben, M. Klein, M. Köhli, P. Lennert, M. Meven, C. Schmidt, U. Schmidt, W. Schweika, Nuclear Instruments and Methods in Physics Research Section A 743 (2014) 90 – 95. doi:10.1016/j.nima.2014.01.007.
  * Kampmann [2012a] R. Kampmann, IEEE NSS, Anaheim N1-109 (2012a).
  * Kampmann [2012b] R. Kampmann, First Tests of Novel Neutron Detectors with Thin Conversion Layers in Inclined Geometry at REFSANS, Oral presentation at the Second International 10B-BF3 Detectors Workshop, 13-14 March, ILL, Grenoble (2012b). URL: https://www.ill.eu/10bbf3.
  * Cippo et al. [2015] E. P. Cippo, G. Croci, A. Muraro, A. Menelle, G. Albani, M. Cavenago, C. Cazzaniga, G. Claps, G. Grosso, F. Murtas, M. Rebai, M. Tardocchi, G. Gorini, Journal of Instrumentation 10 (2015) P10003.
  * Croci and et al. [2014] G. Croci, et al., EPL 107 (2014) 12001. doi:10.1209/0295-5075/107/12001.
  * Anastasopoulos et al. [2017] M. Anastasopoulos, R. Bebb, K. Berry, J. Birch, T. Bryś, J.-C. Buffet, J.-F. Clergeau, P. Deen, G. Ehlers, P. van Esch, S. Everett, B. Guerard, R. Hall-Wilton, K. Herwig, L. Hultman, C. Höglund, I. Iruretagoiena, F. Issa, J. Jensen, A. Khaplanov, O. Kirstein, I. L. Higuera, F. Piscitelli, L. Robinson, S. Schmidt, I. Stefanescu, Journal of Instrumentation 12 (2017) P04030.
  * Andersen et al. [2013] K. Andersen, T. Bigault, J. Birch, J. Buffet, J. Correa, R. Hall-Wilton, L. Hultman, C. Höglund, B. Guérard, J. Jensen, A. Khaplanov, O. Kirstein, F. Piscitelli, P. V. Esch, C. Vettier, Nuclear Instruments and Methods in Physics Research Section A 720 (2013) 116 – 121. doi:10.1016/j.nima.2012.12.021.
  * Birch et al. [2014] J. Birch, J.-C. Buffet, J.-F. Clergeau, J. Correa, P. van Esch, M. Ferraton, B. Guerard, J. Halbwachs, R. Hall-Wilton, L. Hultman, C. Höglund, A. Khaplanov, M. Koza, F. Piscitelli, M. Zbiri, Journal of Physics: Conference Series 528 (2014) 012040.
  * Birch et al. [2013] J. Birch, J. C. Buffet, J. Correa, P. van Esch, B. Guérard, R. Hall-Wilton, C. Höglund, L. Hultman, A. Khaplanov, F. Piscitelli, IEEE Transactions on Nuclear Science 60 (2013) 871–878. doi:10.1109/TNS.2012.2227798.
  * Guerard and Buffet [2011] B. Guerard, J. Buffet, Ionizing radiation detector - patent no. 20110215251., 2011\. URL: https://www.google.com/patents/US20110215251, US Patent App. 13/038,915.
  * Klein and Schmidt [2011] M. Klein, C. J. Schmidt, Nuclear Instruments and Methods in Physics Research Section A 628 (2011) 9–18. doi:10.1016/j.nima.2010.06.278, {VCI} 2010Proceedings of the 12th International Vienna Conference on Instrumentation.
  * Klein [2000] M. Klein, Experimente zur Quantenmechanik mit ultrakalten Neutronen und Entwicklung eines neuen Detektors zum ortsaufgeloesten Nachweis von thermischen Neutronen auf grossen Flächen, Ph.D. thesis, University of Heidelberg, 2000.
  * Stefanescu et al. [2013a] I. Stefanescu, Y. Abdullahi, J. Birch, I. Defendi, R. Hall-Wilton, C. Höglund, L. Hultman, D. Seiler, K. Zeitelhack, Nuclear Instruments and Methods in Physics Research Section A 727 (2013a) 109 – 125. doi:10.1016/j.nima.2013.06.003.
  * Stefanescu et al. [2013b] I. Stefanescu, Y. Abdullahi, J. Birch, I. Defendi, R. Hall-Wilton, C. Höglund, L. Hultman, M. Zee, K. Zeitelhack, Journal of Instrumentation 8 (2013b) P12003.
  * Pfeiffer et al. [2015] D. Pfeiffer, F. Resnati, J. Birch, R. Hall-Wilton, C. Höglund, L. Hultman, G. Iakovidis, E. Oliveri, E. Oksanen, L. Ropelewski, P. Thuiner, Journal of Instrumentation 10 (2015) P04004.
  * Lacy et al. [2009] J. L. Lacy, L. Sun, C. S. Martin, A. Athanasiades, T. D. Lyons, in: 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC), pp. 1117–1121. doi:10.1109/NSSMIC.2009.5402421.
  * Lacy et al. [2002] J. L. Lacy, A. Athanasiades, N. N. Shehad, R. A. Austin, C. S. Martin, IEEE NSS, Norfolk 1 (2002) 392–396.
  * Lacy [2006] J. Lacy, IEEE Nuclear Science Symposium Conference Record 1 (2006) 20–26.
  * Lacy et al. [2011] J. L. Lacy, A. Athanasiades, L. Sun, C. S. Martin, T. D. Lyons, M. A. Foss, H. B. Haygood, Nuclear Instruments and Methods in Physics Research Section A 652 (2011) 359 – 363. doi:10.1016/j.nima.2010.09.011, symposium on Radiation Measurements and Applications (SORMA) {XII} 2010.
  * Lacy et al. [2012] J. L. Lacy, et al., IEEE NSS, Anaheim He-1-6 (2012).
  * Lacy et al. [2013] J. L. Lacy, A. Athanasiades, C. S. Martin, L. Sun, G. L. Vazquez-Flores, IEEE Trans. Nucl. Sci. 60 (2013) 1140–1146.
  * McGregor et al. [2003] D. McGregor, M. Hammig, Y.-H. Yang, H. Gersch, R. Klann, Nuclear Instruments and Methods in Physics Research Section A 500 (2003) 272 – 308. doi:10.1016/S0168-9002(02)02078-8.
  * Piscitelli and Esch [2013] F. Piscitelli, P. V. Esch, Journal of Instrumentation 8 (2013) P04020.
  * Kanaki et al. [2017] K. Kanaki, T. Kittelmann, X. X. Cai, E. Klinkby, E. B. Knudsen, P. Willendrup, R. Hall-Wilton, Arxiv preprint hep-th (2017). arXiv:1708.02135.
  * Piscitelli and Esch [2013] F. Piscitelli, P. V. Esch, Journal of Instrumentation 8 (2013) P04020.
  * Piscitelli [3133] F. Piscitelli, Boron-10 layers, Neutron Reflectometry and Thermal Neutron Gaseous Detectors, Ph.D. thesis, Institut Laue-Langevin and University of Perugia, 2014 - arXiv:1406.3133.
  * scipy [????] scipy, ????  URL: https://www.scipy.org/.
  * numpy [????] numpy, ????  URL: http://www.numpy.org/.
  * matplotlib [????] matplotlib, ????  URL: https://matplotlib.org/.
  * tutorialspoint [????] tutorialspoint, ????  URL: https://www.tutorialspoint.com/design_pattern/builder_pattern.htm.
  * Ziegler [????] J. F. Ziegler, ????  URL: http://www.srim.org/.


# Contact-GraspNet: Efficient 6-DoF Grasp Generation
in Cluttered Scenes

Martin Sundermeyer${}^{1,2,3}\qquad$ Arsalan Mousavian${}^{1}\qquad$ Rudolph
Triebel${}^{2,3}\qquad$ Dieter Fox1,4 ∗This work is done while the first
author was an intern at NVIDIA. 1NVIDIA<EMAIL_ADDRESS>2German Aerospace Center (DLR<EMAIL_ADDRESS>3Technical University of
Munich (TUM), 4University of Washington

###### Abstract

Grasping unseen objects in unconstrained, cluttered environments is an
essential skill for autonomous robotic manipulation. Despite recent progress
in full 6-DoF grasp learning, existing approaches often consist of complex
sequential pipelines that possess several potential failure points and run-
times unsuitable for closed-loop grasping. Therefore, we propose an end-to-end
network that efficiently generates a distribution of 6-DoF parallel-jaw grasps
directly from a depth recording of a scene. Our novel grasp representation
treats 3D points of the recorded point cloud as potential grasp contacts. By
rooting the full 6-DoF grasp pose and width in the observed point cloud, we
can reduce the dimensionality of our grasp representation to 4-DoF which
greatly facilitates the learning process. Our class-agnostic approach is
trained on 17 million simulated grasps and generalizes well to real world
sensor data. In a robotic grasping study of unseen objects in structured
clutter we achieve over 90% success rate, cutting the failure rate in half
compared to a recent state-of-the-art method. Video of the real world
experiments and code are available at
https://research.nvidia.com/publication/2021-03_Contact-GraspNet%3A--
Efficient.

##  I Introduction

The ability to grasp objects is one of the fundamental capabilities required
in most robot manipulation tasks. Grasping involves reasoning about the 3D
geometry and physics properties of the object such as mass and friction, and
also reasoning about complex contact physics. It is studied in two main
directions: Model-based grasping where the 3D model or category of the object
is known and model-free grasping where there is no prior knowledge about the
object. Model-based grasping circumvents reasoning about the physics of
contact and grasp generation by pre-defining a set of grasps in the object
frame and transform those grasps according to the 6-DoF object pose [1, 2, 3,
4] or detected keypoints of the objects [5, 6]. The downside of model-based
approaches is that they only work on a limited subset of known objects or
categories, and any errors in detecting 6-DoF object pose or object keypoints
degrade the grasping performance.

Figure 1: Contact-GraspNet efficiently predicts diverse and stable grasps in
cluttered scenes while avoiding collisions.

Model-free approaches do not make any strong assumptions about the category or
shape of the object, and they learn a shared representation for all object
shapes and sizes. However, having one shared representation for all objects in
addition to the large SE(3) space for the grasp poses makes the learning
problem quite challenging. As a result, a large body of work in data-driven
grasping constraints the space of possible grasps to planar grasping, where
grasps are represented by oriented rectangles around each pixel that define
the grasp frame [7, 8, 9]. Such a representation needs the camera to view the
scene perpendicularly and thus limits 3D reasoning and applications
significantly. A large number of possible grasps and the full kinematic
capabilities of the robot are also neglected. To address the limitations of
planar grasping, there has been a recent interest in tackling the problem of
6-DoF grasping of unknown objects [10, 11, 12, 13, 14]. In this paper, we
tackle 6-DoF grasping of unknown objects in cluttered space from a partial
point cloud observation of the scene.

Grasping objects from cluttered scenes with structure introduces extra
challenges. The target objects must be grasped successfully, while at the same
time any collision with other objects must be avoided to prevent damages or
transformations into other undesired states. This is particularly important in
home robotics and healthcare applications. Additionally, it is crucial to
generate a diverse set of grasps for the object due to robot kinematic
constraints. Depending on the relative pose between the object and the robot,
a different subset of grasps is kinematically feasible.

Our method is closely related to the work of Murali _et al._ [12], where the
goal is to generate collision-free diverse grasps for a designated target
object from a partial point cloud of the scene, and the objects are segmented
using a pre-trained unknown object instance segmentation model [15, 16].
Murali _et al._ [12] use a multi-stage process that synthesizes grasps for the
target objects from the segmented object point cloud with no context around
it, and then filters out the colliding grasps using another learned model.
This leads to three issues: 1) Sensitivity to instance segmentation errors. 2)
Grasps are generated just from the target object point cloud and do not
leverage geometric cues in the scene such as table points and surrounding
object points. 3) Grasps are predicted in the large, unconstrained 6-DoF pose
space. To address these issues, our method instead directly processes a full
scene point cloud or a local region around a target object. Therefore, the
quality of our generated grasps is not depending on an accurate mask and
collisions can be directly taken into account during generation. Instance
segmentation can then subsequently be used to filter grasps belonging to a
target object. Thus, our main contributions are the following:

  * •

A new end-to-end method for 6-DoF grasping of unknown objects in cluttered
real world scenes where we achieve $90\%$ grasp success rate. This is $10\%$
higher than [12] in equal settings.

  * •

A new grasp pose representation that projects 6-DoF grasps to their contact
points in an observed point cloud. Our representation has only 4-DoF which
facilitates the learning problem significantly.

  * •

Comprehensive ablation studies in a physics simulator to evaluate the effects
of different loss functions and training data.

##  II Related Work

As a fundamental problem in robotics, grasping has been studied for decades
[17, 18, 19, 20]. We review related literature in the context of data-driven
methods.

End-to-end policy learning: One line of work for grasping and manipulation of
objects employs an end-to-end policy that learns to generate actions from raw
input pixel values [21, 22]. This results in a monolithic model that
concurrently reasons about perception, planning, grasping, and controlling the
robot. A large group of these works learn from interactions of the robot with
the environment through reinforcement learning. These approaches have mostly
shown promise in bin picking, in (quasi-) planar grasping and in small,
insensible workspaces that do not require complicated motion planning in the
robot configuration space. Few works [23] have demonstrated iterative 6-DoF
grasping approaches with a monolithic policy by combining imitation learning
and reinforcement learning. A common drawback of these methods is the limited
generalization to novel environments, because the perception and control are
learned indirectly at the same time. In addition, these methods are not easily
steerable towards grasping a specific object as the reward function encourages
grasping any object. In contrast, our method learns to generate diverse 6-DoF
grasps on novel objects and scenes for specifiable target objects while just
using simulated training data. Additionally, it can be integrated with other
perception and motion planning algorithms.

3D reconstruction: A complete 3D reconstruction enables traditional grasp
planning. However, learned single-view reconstructions are often ambiguous,
coarse and require class-conditioning [24, 25, 26]. Multiple views for 3D
scanning are beneficial [27] but not always obtainable, take additional time
and typically assume a static scene. In our approach a full explicit 3D
reconstruction is not required.

Discriminative methods: Discriminative methods for grasping train a classifier
that evaluates the quality of existing grasps [28, 29, 7]. They use different
sampling strategies to generate potential candidates. For planar grasping,
cross entropy is widely used since it can converge to the final grasp location
by iteratively evaluating the quality of grasps in different locations [7].
However, the cross-entropy method does not work well in the higher dimensional
6-DoF grasp space. To overcome the sampling complexity issue, grasp locations
are often sampled using geometric heuristics [10, 29].

Generative methods: Learning-based generative grasp methods aim to overcome
the limitations of geometric heuristics and generate meaningful 6-DoF grasps
often from experience in a physics simulator [11, 12]. The main challenge is
the large, multi-modal search space of 6-DoF grasps. Instead of sampling some
potential candidates using heuristics and ranking them, these models directly
predict a per-point graspability score and approach direction in SO(3) space
[30, 14, 31]. One problem with predicting approach directions is that they
cannot easily capture high curvature areas such as mug rims or handles and
also can not represent grasps encompassing hollow structures. Furthermore,
successful approach directions are quite ambiguous to learn as multiple ones
are possible for a single contact. Instead we propose to predict 6-DoF grasps
densely projected to their much less ambiguous contact points. While grasps
without full surface contact are plausible, e.g. through the handle of a mug,
the knowledge about the object state and therefore the ability to steadily
place the object again is lost. Therefore, in this work we are aiming to
generate stable grasps for unknown objects with full surface contact. Our
novel loss formulation further improves convergence by accounting for the
discontinuities, imbalance and multi modality of the grasp distribution.
Unlike other methods [31], our proposed method is independent of category
labels and has no assumption of grasps being always perpendicular to a
surface. Instead we learn a grasp semantic purely from a wide variety of grasp
annotated training shapes [32].

##  III Method

Figure 2: Training Data Pipeline. We place object meshes with dense grasp
annotations from the ACRONYM dataset [32] at random stable poses in scenes.
Grasp poses that produce gripper model collisions are removed. Resulting
grasps are mapped to their contacts on the mesh surface. During training, we
sample virtual cameras to render point clouds from the scenes. We consider
recorded points (yellow) as positive contacts if there exists a mesh contact
(blue) in a 5mm radius and associate the grasp transformation belonging to the
closest mesh contact to them. These per-point annotations are used to
supervise the Contact Grasp Network. Figure 3: Our grasp representation: $c$
depicts an observed contact point. $\mathbf{a}$ and $\mathbf{b}$ constitute
the 3-DoF rotation, $w$ is the predicted grasp width, $d$ the distance from
baseline to base frame. In pink we show the five gripper points $\mathbf{v}$
that we used in the $l_{add-s}$ loss.

We consider the problem of generating 6-DoF grasps from any viewpoint on
structured clutter consisting of unknown objects. Our approach takes in a raw
depth image, optionally with object masks, and generates 6-DoF grasp proposals
together with corresponding grasp widths. Our goal is to predict grasps that
are robust, diverse and non-colliding from an only partially observable scene.

From a learning perspective, generating the distribution of successful 6-DoF
grasps is quite challenging, because the distribution is multi-modal,
discontinuous, imbalanced and ambiguous due to (self-) occlusions.
Furthermore, direct regression in high dimensional output spaces like $SE(3)$
has been shown to be difficult in grasping [11] and also in related fields
such as object pose estimation [33].

###  III-A Grasp Representation

For these reasons, finding an efficient grasp representation is crucial to
solve this task using learning-based methods. This representation should
generalize well to unseen objects and handle the high-dimensional output space
well.

Contact Grasp Representation:  We observe that for most predictable two-finger
grasps at least one of the two contacts is visible prior to grasping. In
contrast, grasps without any visible contact are often ambiguous or do not
preserve the initial object pose after grasping. Therefore, we map a
distribution of successful 6-DoF ground truth grasps $g\in G$ to their
corresponding contact points $c\in\mathbb{R}^{3}$. Since visible contact
points are bound to lie on surfaces that we can observe with a depth sensor,
we can represent their 3D location by nearby points in a recorded point cloud.

Given that we can predict whether observed points are suitable grasp contacts,
we can thus reduce the 6-DoF grasp learning problem to estimating the 3-DoF
grasp rotation $R_{g}\in\mathbb{R}^{3\times 3}$ and grasp width
$w\in\mathbb{R}$ of a parallel-yaw gripper.

Starting from a contact point $\mathbf{c}\in\mathbb{R}^{3}$, where the gripper
baseline intersects the mesh, we depict a 6-DoF grasp pose $g\in G$ defined by
$(R_{g},t_{g})\in SE(3)$ and grasp width $w\in\mathbb{R}$ as

$\mathbf{t}_{g}=\mathbf{c}+\frac{w}{2}\mathbf{b}+d\mathbf{a}$ (1)
$R_{g}=\begin{bmatrix}|&|&|\\\
\mathbf{b}&\mathbf{a}\times\mathbf{b}&\mathbf{a}\\\ |&|&|\\\ \end{bmatrix},$
(2)

where $\mathbf{a}\in\mathbb{R}^{3},||\mathbf{a}||=1$ is the approach vector,
$\mathbf{b}\in\mathbb{R}^{3},||\mathbf{b}||=1$ is the grasp baseline vector,
and $d\in\mathbb{R}$ is the constant distance from the gripper baseline to the
gripper base. Our grasp representation is depicted in Figure 3.

The reduced dimensionality greatly facilitates the learning process compared
to methods that estimate grasp poses in unconstrained $SE(3)$ space. It also
increases the pose accuracy of predicted grasps as they are bound to the
geometry of the observed scene. In contrast to axis-angle representations, our
rotation representation has neither ambiguities nor discontinuities. Moreover,
at test time we can sample grasp proposals by sampling contact points that
cover the whole observable surface of the scene/object and thus represent the
modes of the 6-DoF grasp distribution well. While a 3D view on the scene is
preferable, even a frontal view on a box produces reasonable grasps due to the
radial mapping.

Point Set Networks such as PointNet++ [34] effectively process point clouds
and hierarchically aggregate points and their feature representations in local
3D neighborhoods. Their predictions can be directly associated to 3D points in
the input point cloud and our proposed grasp representation exploits this
ability.

###  III-B Data Generation

To learn the full distribution of stable 6-DoF grasps, diverse and dense grasp
pose annotations are required. We used the ACRONYM dataset [32], which
consists of 8872 meshes from the Shapenet dataset [35] and $17.7$ million
simulated grasps under varying friction. An overview of our offline and online
training data generation is given in Fig. 2.

During training we render a scene point cloud
$\mathcal{P}=\\{\mathbf{p}_{1},\dots,\mathbf{p}_{n}\\}\subset\mathbb{R}^{3}$
and assign a point-wise grasp success

$\forall i=1,\dots,n\quad
s_{i}=\begin{cases}1&\min_{j}||\mathbf{p}_{i}-\mathbf{c}_{j}||_{2}<r\\\
0&otherwise,\end{cases}$ (3)

where $\mathbf{c}_{j}\in\mathcal{P}$ are the mesh contact points of non-
colliding ground truth grasps $g_{j}\in G$ in camera coordinates and
$r\in\mathbb{R}$ is their maximum propagation radius. Thus, $\mathcal{P}$ can
be split into points $\mathcal{P}^{-}:=\\{\mathbf{p}_{i}|s_{i}=0\\}$, where no
feasible grasp contact is found within a radius of $r=5mm$, and
$\mathcal{P}^{+}:=\\{\mathbf{p}_{i}|s_{i}=1\\}$, containing points suitable
for a contact. To the latter ones $\mathbf{p}^{+}_{i}\in\mathcal{P}^{+}$ we
assign the closest grasp as

$\displaystyle\begin{bmatrix}w_{g,i}\\\ R_{g,i}\\\
\mathbf{t}_{g,i}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}w_{g,j}\\\
R_{g,j}\\\
\mathbf{p}^{+}_{i}+\frac{w_{j}}{2}\mathbf{b}_{j}+d\mathbf{a}_{j}\end{bmatrix}$
(4)

with

$\displaystyle
j=\operatorname*{arg\,min}_{k}||\,\mathbf{p}^{+}_{i}-\mathbf{c}_{k}||_{2}$ (5)

Given sufficient coverage we can thereby project the ground truth distribution
of 6-DoF grasps densely on the recorded point cloud.

###  III-C Network

We employ the set abstraction and feature propagation layers proposed in
PointNet++ [34] to build an asymmetric U-shaped network. The network takes
n=20000 random points $p\in\mathbb{R}^{20000\times 3}$ as input and predicts
grasps for only m=2048 farthest points of the input to make sure the inference
fits in GPU memory and predicted grasps have good coverage over the scene. The
network has four heads with two 1D-Conv layers each and per-point outputs
$s\in\mathbb{R},\mathbf{z}_{1}\in\mathbb{R}^{3},\mathbf{z}_{2}\in\mathbb{R}^{3},\mathbf{o}\in\mathbb{R}^{10}$,
from which we form our grasp representation. The predicted grasp width
$\mathbf{\hat{w}}_{i}\in[0,w_{max}]$ is split into 10 equidistant grasp width
bins $\mathbf{\hat{o}}\in\mathbb{R}^{10}$ to counteract data imbalance. Then,
$\mathbf{\hat{w}}_{i}$ is represented by the center value of the bin(s) with
the highest confidence. The approach direction $\mathbf{a}\in\mathbb{R}^{3}$
and the baseline direction $\mathbf{b}\in\mathbb{R}^{3}$ are orthonormal by
definition. We inject this property into training by coupling the predictions
$\mathbf{\hat{a}},\mathbf{\hat{b}}$ through an in-network Gram Schmidt
orthonormalization

$\displaystyle\mathbf{\hat{b}}=$
$\displaystyle\frac{\mathbf{z}_{1}}{||\mathbf{z}_{1}||}$
$\displaystyle\mathbf{\hat{a}}=\frac{\mathbf{\mathbf{z}_{2}-\langle\hat{b},\mathbf{z}_{2}\rangle\hat{b}}}{||\mathbf{z}_{2}||}$
(6)

Thus, we perform a projection and only predict $\mathbf{\hat{a}}$ as the
component that is orthonormal to $\mathbf{\hat{b}}$. The orthonormalization
further reduces the dimensionality of our predicted grasp representation and
facilitates the regression of 3D rotations [36].

Figure 4: Full Inference Pipeline: We segment unknown objects from an RGB-D
image using [15]. Our Contact-GraspNet processes the full scene point cloud or
a local region of interest around a target object. Predicted 6-DoF grasps are
then associated to object segments by filtering their contact points. On the
right we show the predicted 6-DoF grasp distribution and, in bold, the most
confident grasp per segment.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0$$0.2$$0.4$$0.6$$0.8$$1$CoverageSuccess
RateSimulator Success-CoverageContact GraspNetw/o $l_{add-s}$w/o binned
$l_{width}$

Figure 5: Loss Ablations: Without weighted binning in the grasp width loss
$l_{width}$ both, success rate and coverage decrease. The $l_{add-s}$ loss
leads to increased success rates at high confidence contacts (Coverage
$\in[0,0.1]$) and to slightly decreased success rate in the low-confidence
regime. This confidence calibration is important, since it determines which
grasp is eventually executed.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0$$0.2$$0.4$$0.6$$0.8$$1$CoverageSuccessSimulator
Success-CoverageContact GraspNetw/ noise $\sigma=0.001$w/o local regionsours
trained on [11][12] trained on [11]

Figure 6: Data Ablations: Training with Gaussian noise has similar performance
in simulation but helps generalization to noisy sensor data. Predicting grasps
directly on full scenes without extracting local regions yields a similar
average success rate, but significantly lowers grasp coverage. Training on the
small grasp dataset from [11] with 5 categories is not sufficient to
generalize to arbitrary objects and shows the importance of ACRONYM [32]

###  III-D Target Losses

The contact grasp success predictions $\hat{s}\in\mathbb{R}$ are evaluated at
all output points $\mathbf{p}_{i}\in\mathbb{R}^{3}:\forall i\in[0,m]$ using
binary cross entropy. We only backpropagate the top-k point predictions with
the largest errors $l_{bce,k}$, with k=512, to counteract data imbalance. The
other predictions concerning the geometry of grasps are only evaluated at
positive contact points $\mathbf{p}^{+}_{i}$. Instead of supervising all
network heads in isolation, we propose to combine the predictions to the 6-DoF
grasp pose $\hat{g}\in G$ given in Eq. (1) and (2) already during training. We
define five 3D points $\mathbf{v}\in\mathbb{R}^{5\times 3}$ representing the
6-DoF gripper pose, as shown in Fig. 3, and transform these using all ground
truth and predicted grasp poses defined in Eq. (4)

$\displaystyle\mathbf{v}^{gt}_{i}=$
$\displaystyle\mathbf{v}R_{g,i}^{T}+\mathbf{t}_{g,i}$
$\displaystyle\mathbf{v}^{pred}_{i}=\mathbf{v}\hat{R}_{g,i}^{T}+\mathbf{\hat{t}}_{g,i}$
(7)

We formulate the 6-DoF grasp loss $l_{add-s}$ as a weighted minimum average
distance between gripper points $\mathbf{v}^{gt}$ and $\mathbf{v}^{pred}$
where we take the symmetry of the gripper into account.

$\displaystyle
l_{add-s}=\frac{1}{n^{+}}\sum^{n^{+}}_{i}\hat{s_{i}}\min_{u}||\mathbf{v}^{pred}_{i}-\mathbf{v}^{gt}_{u}||_{2},$
(8)

where $n^{+}$ is the size of $\mathcal{P}^{+}$. We weight each distance to the
closest ground truth grasp points with the predicted contact success
confidence $\hat{s}_{i}$.

Our proposed loss formulation has several advantages: (1) We can learn the
different modes of the ground truth grasp distribution, e.g. different
predicted grasp approach directions $\hat{\mathbf{a}}$ can produce a small
error. (2) The point-wise weighting with $\hat{s}_{i}$ couples the contact
point classification with the grasp pose predictions. Contact confidence can
only increase if the network predicts a 6-DoF grasp pose close to a ground
truth pose. (3) Wrongly predicted grasps in regions far away from any ground
truth grasp, e.g. at artificial edges from occlusions, produce a high loss and
are thus avoided.

On the grasp width bin predictions, we optimize a weighted, multi-label binary
cross entropy loss $l_{width}$. Since small grasp widths are highly over-
represented, we weight the bin losses anti-proportional to bin size. Our total
loss is $l=\alpha l_{bce,k}+\beta l_{add-s}+\gamma l_{width}$ with
$\alpha=1,\beta=10,\gamma=1$.

###  III-E Implementation Details

We use the Adam optimizer with an initial learning rate of 0.001 and a step-
wise decay to 0.0001. Our set abstraction layers have 3 parallel branches with
query ball radii [0.02,0.04,0.08], [0.04,0.08.0.16] and [0.08,0.16,0.32]. For
inference the point cloud is centered at its mean in camera coordinates. For
training we generate 10000 table top scenes by placing 8-12 grasp annotated
ShapeNet models [32] at random stable poses. We use rejection sampling to
avoid collisions. We train with a batch size of 3 for 144.000 iterations which
takes $\sim 40$ hours on a single Nvidia V100 GPU. Convergence is
significantly faster than on previous methods [12, 14, 11] which take up to
one week on a single GPU for training. This also reflects the effectiveness of
our proposed grasp representation.

##  IV Experimental Evaluation

We evaluate our method in a grasping study with a Franka robot where we pick
unknown objects in cluttered scenes. We also compare different variations of
our method and of our data by executing a large number of predicted grasps in
the FleX physics simulator [37].

###  IV-A Inference

Our inference pipeline is shown and described in Fig. 4. The Contact-GraspNet
can also be applied to raw depth images by itself, but most robotic tasks
require some kind of instance detection/segmentation to specify a target.

Local regions of interest can be optionally extracted around the 3D centroid
of point cloud segments in order to maximize the number of potential contact
points. In our experiments, we extract cubes with an edge length set to twice
the largest spanning dimension, but at least $0.3m$ and at most $0.6m$.

Run time: The Contact-GraspNet has a run time of $0.28s$ for a full scene or
$\sim 0.19s$ for a local region around a target object. Compared to other
6-DoF grasp generation methods this is quite fast and enables applications
requiring reactive closed loop grasping.

Grasp Selection: At test time we select grasps by setting a contact confidence
threshold of $0.23$ and then use farthest point sampling on the (filtered)
contact points to ensure broad grasp coverage. If the number of predicted
grasps for an object is too low, we reduce the confidence threshold to $0.19$.
In the end we execute the most confident grasp that is kinematically reachable
and where the robot does not collide with the scene [38].

Figure 7: One advantage of our method is that it does not rely on an accurate
segmentation of unknown objects. Here, successful grasp contacts are still
found on the driller despite severe under-segmentation.

###  IV-B Evaluation Metrics

In our robotic experiments we report the number of successful grasps and the
number of trials. The latter is often disregarded when picking small objects
from a bin. However, grasping in only one or two trials is crucial in
cluttered scenes (e.g. in households) with large, densely packed objects where
collisions should be avoided and stable grasp opportunities can vanish after
objects tip over. We limit ourselves to a maximum of two grasp trials per
object without rearrangements and report the success rate after a single trial
as well.

Our simulator experiments allow us to also evaluate the diversity of grasps
and ablate variations of our method. Here, we evaluate the success rate and
coverage of the generated grasps following [11]. A grasp is considered
successful if (1) the open gripper does not collide with the object/scene and
(2) the object is still in the gripper after grasping and a shaking motion.
This is a conservative measure, as most real world grasps can slightly collide
and do not undergo a shaking motion. Coverage is the percentage of ground
truth grasps (including occluded ones) whose base coordinates are within 2cm
of any of the generated grasps.

###  IV-C Real robot grasp experiments

Setup: Our physical setup consists of a 7-DoF Franka Panda robot with a
parallel-jaw gripper. We closely replicate the 9 cluttered scenes defined in
[12] with a total of 51 unseen objects. The task is to pick the objects from
the cluttered scene and place them into a bin. We manually select target
objects and grasp them in the same random order as in [12]. In our
experiments, we use the Intel Realsense L515 LiDAR camera mounted on a tripod
for both RGB and depth data. Robot motions are generated using [38].

Results: Table I shows our grasp evaluation results on the robot. We observe a
significantly higher grasp success rate of our method compared to [11] and
[12] which themselves outperform other learning-based methods and
analytic/heuristic baselines. Furthermore, our method strongly improves the
grasp success at first trial and thereby reduces the number of re-grasps. We
also addressed the shortcomings of cropping objects from the point cloud using
potentially imprecise segmentation masks. Fig. 7 shows an imprecise
segmentation example where cropping would be catastrophic but where our grasp
filtering method can still extract successful grasps.

###  IV-D Ablations

Optimization Targets: In Fig. 5 we first investigate the effect of our loss
targets. The weighted loss on the grasp width bins $l_{width}$ is crucial to
deal with the imbalanced widths in our grasp dataset. Without weighting the
bins, the predictions mostly collapse into narrow grasp widths. Weighting also
performs better than oversampling in our experiments. The average distance
loss $l_{add-s}$ improves the success rate of high confidence contacts which
is important because most grasps that we execute lie in the first decimal of
coverage. The connection of contact confidence with the grasp pose results in
an overall improved calibration.

Data: In Fig. 6 we examine the effects of different training and test data.
Zooming into local regions allows the network to concentrate potential contact
points on the object and thus increases coverage. We also show the importance
of a large and diverse grasp dataset like ACRONYM [32]. Training on a small
grasp datasets with $110$ objects from $5$ categories [11] is not sufficient
for out-of-category generalization irrespective of the method.

Failure Cases: We observe some failure cases for thick objects that only allow
grasps almost at maximum grasp width. Here, grasp predictions are less
confident presumably because of the discontinuous decision boundary. Injecting
noise during training reduces this effect. Finally, small objects sometimes
have contact points with low confidence possibly because of their small impact
on the total loss.

TABLE I: Cluttered Scene Grasping: We achieve a clear improvement over recent state-of-the-art grasping pipelines | Success | First attempt | #Attempts
---|---|---|---
6-DOF GraspNet[11] | 62.7 | - | -
[11] +CollisionNet[12] | 80.39 | 68.63 | 67
Contact-GraspNet | $\boldsymbol{90.20}$ | $\boldsymbol{84.31}$ | $\boldsymbol{59}$

##  V Conclusions

We considered the fundamental problem of grasping unknown objects in
structured clutter with a parallel jaw gripper. We proposed an efficient,
accurate and simplifying 6-DoF grasp generation method called Contact-
GraspNet. By transforming the hardly tractable 6-DoF grasp estimation problem
into a grasp contact point classification and a grasp rotation estimate, we
greatly limit the predicted pose space and facilitate the learning process.
Through tailored optimization targets that take into account the multi-
modality, imbalance and sparsity of the 6-DoF grasp distribution, our network
learns to generate diverse grasps covering the whole graspable surface in a
recorded scene. Gripper collisions are effectively avoided by considering them
during training and by predicting grasps directly in scenes. Our approach can
incorporate segmentation predictions as well but is not dependent on accurate
masks itself. It is also complementary to grasp ranking methods that use
gripper and/or robot models as input. Grasping successfully with a single
attempt is crucial in sensible environments. Our method showed strong advances
in that regard and is a step towards reaching the required grasp reliability.

## References

  * [1] X. Deng, A. Mousavian, Y. Xiang, F. Xia, T. Bretl, and D. Fox, “Poserbpf: A rao-blackwellized particle filter for 6d object pose tracking,” in _Robotics: Science and Systems (RSS)_ , 2019.
  * [2] M. Sundermeyer, Z.-C. Marton, M. Durner, and R. Triebel, “Augmented autoencoders: Implicit 3d orientation learning for 6d object detection,” _International Journal of Computer Vision_ , vol. 128, no. 3, pp. 714–729, 2020.
  * [3] X. Deng, Y. Xiang, A. Mousavian, C. Eppner, T. Bretl, and D. Fox, “Self-supervised 6d object pose estimation for robot manipulation,” in _International Conference on Robotics and Automation (ICRA)_ , 2020.
  * [4] C. Wang, D. Xu, Y. Zhu, R. Martín-Martín, C. Lu, L. Fei-Fei, and S. Savarese, “Densefusion: 6d object pose estimation by iterative dense fusion,” 2019.
  * [5] L. Manuelli, W. Gao, P. Florence, and R. Tedrake, “kpam: Keypoint affordances for category-level robotic manipulation.” _International Symposium on Robotics Research (ISRR)_ , 2019.
  * [6] K. Fang, Y. Zhu, A. Garg, A. Kuryenkov, V. Mehta, L. Fei-Fei, and S. Savarese, “Learning task-oriented grasping for tool manipulation from simulated self-supervision,” _Robotics: Science and Systems (RSS)_ , 2018.
  * [7] J. Mahler, J. Liang, S. Niyaz, M. Laskey, R. Doan, X. Liu, J. A. Ojea, and K. Goldberg, “Dex-net 2.0: Deep learning to plan robust grasps with synthetic point clouds and analytic grasp metrics,” _RSS_ , 2017.
  * [8] J. Mahler, M. Matl, V. Satish, M. Danielczuk, B. DeRose, S. McKinley, and K. Goldberg, “Learning ambidextrous robot grasping policies,” _Science Robotics_ , vol. 4, no. 26, p. eaau4984, 2019.
  * [9] I. Lenz, H. Lee, and A. Saxena, “Deep learning for detecting robotic grasps,” _IJRR_ , 2015.
  * [10] A. ten Pas, M. Gualtieri, K. Saenko, and R. Platt, “Grasp pose detection in point clouds,” _The International Journal of Robotics Research_ , vol. 36, no. 13-14, pp. 1455–1473, 2017.
  * [11] A. Mousavian, C. Eppner, and D. Fox, “6-dof graspnet: Variational grasp generation for object manipulation,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 2901–2910.
  * [12] A. Murali, A. Mousavian, C. Eppner, C. Paxton, and D. Fox, “6-dof grasping for target-driven object manipulation in clutter,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6232–6238.
  * [13] S. Song, A. Zeng, J. Lee, and T. Funkhouser, “Grasping in the wild: Learning 6dof closed-loop grasping from low-cost demonstrations,” _Robotics and Automation Letters_ , 2020.
  * [14] H.-S. Fang, C. Wang, M. Gou, and C. Lu, “Graspnet-1billion: A large-scale benchmark for general object grasping,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 444–11 453.
  * [15] Y. Xiang, C. Xie, A. Mousavian, and D. Fox, “Learning rgb-d feature embeddings for unseen object instance segmentation,” _Conference on Robotic Learning(CORL)_ , 2020.
  * [16] C. Xie, Y. Xiang, A. Mousavian, and D. Fox, “The best of both modes: Separately leveraging rgb and depth for unseen object instance segmentation,” in _Conference on Robot Learning (CoRL)_ , 2019.
  * [17] D. Prattichizzo and J. J. Trinkle, _Grasping_ , 01 2008, pp. 671–700.
  * [18] A. Bicchi and V. Kumar, “Robotic grasping and contact: a review,” in _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)_ , vol. 1, 2000, pp. 348–353 vol.1.
  * [19] K. Hang, J. A. Stork, and D. Kragic, “Hierarchical fingertip space for multi-fingered precision grasping,” in _2014 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2014, pp. 1641–1648.
  * [20] G. Du, K. Wang, S. Lian, and K. Zhao, “Vision-based robotic grasp detection from object localization, object pose estimation to grasp estimation: A review,” _arXiv preprint arXiv:1905.06658_ , 2019.
  * [21] S. Levine, P. Pastor, A. Krizhevsky, and D. Quillen, “Learning hand-eye coordination for robotic grasping with deep learning and large-scale data collection,” _International Symposium on Experimental Robotics (ISER)_ , 2016\.
  * [22] D. Kalashnikov, A. Irpan, P. Pastor, J. Ibarz, A. Herzog, E. Jang, D. Quillen, E. Holly, M. Kalakrishnan, V. Vanhoucke, and S. Levine, “Qt-opt: Scalable deep reinforcement learning for vision-based robotic manipulation,” _Conference on Robot Learning_ , 2018.
  * [23] A. Zeng, S. Song, S. Welker, J. Lee, A. Rodriguez, and T. Funkhouser, “Learning synergies between pushing and grasping with self-supervised deep reinforcement learning,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 4238–4245.
  * [24] X. Yan, M. Khansari, J. Hsu, Y. Gong, Y. Bai, S. Pirk, and H. Lee, “Data-efficient learning for sim-to-real robotic grasping using deep point cloud prediction networks,” 2019.
  * [25] X. Yan, J. Hsu, M. Khansari, Y. Bai, A. Pathak, A. Gupta, J. Davidson, and H. Lee, “Learning 6-dof grasping interaction via deep geometry-aware 3d representations,” 2018.
  * [26] W. Agnew, C. Xie, A. Walsman, O. Murad, C. Wang, P. Domingos, and S. Srinivasa, “Amodal 3d reconstruction for robotic manipulation via stability and connectivity,” 2020.
  * [27] M. Breyer, J. J. Chung, L. Ott, R. Siegwart, and J. Nieto, “Volumetric grasping network: Real-time 6 dof grasp detection in clutter,” 2021.
  * [28] D. Fischinger, A. Weiss, and M. Vincze, “Learning grasps with topographic features,” _The International Journal of Robotics Research_ , vol. 34, no. 9, pp. 1167–1194, 2015. [Online]. Available: https://doi.org/10.1177/0278364915577105
  * [29] H. Liang, X. Ma, S. Li, M. Görner, S. Tang, B. Fang, F. Sun, and J. Zhang, “PointNetGPD: Detecting grasp configurations from point sets,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2019.
  * [30] Y. Qin, R. Chen, H. Zhu, M. Song, J. Xu, and H. Su, “S4g: Amodal single-view single-shot se (3) grasp detection in cluttered scenes,” in _Conference on robot learning_. PMLR, 2020, pp. 53–65.
  * [31] P. Ni, W. Zhang, X. Zhu, and Q. Cao, “Pointnet++ grasping: Learning an end-to-end spatial grasp generation algorithm from sparse point clouds,” _arXiv preprint arXiv:2003.09644_ , 2020.
  * [32] C. Eppner, A. Mousavian, and F. Dieter, “Acronym: A large-scale grasp dataset based on simulation,” _IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [33] M. Sundermeyer, Z.-C. Marton, M. Durner, M. Brucker, and R. Triebel, “Implicit 3D orientation learning for 6D object detection from RGB images,” in _ECCV_ , 2018.
  * [34] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,” _Neural Information Processing Systems (NeurIPS)_ , 2017.
  * [35] A. X. Chang, T. Funkhouser, L. Guibas, P. Hanrahan, Q. Huang, Z. Li, S. Savarese, M. Savva, S. Song, H. Su, _et al._ , “Shapenet: An information-rich 3d model repository,” _arXiv preprint arXiv:1512.03012_ , 2015.
  * [36] Y. Zhou, C. Barnes, J. Lu, J. Yang, and H. Li, “On the continuity of rotation representations in neural networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 5745–5753.
  * [37] M. Macklin, M. Müller, N. Chentanez, and T.-Y. Kim, “Unified particle physics for real-time applications,” _ACM Transactions on Graphics (TOG)_ , vol. 33, no. 4, pp. 1–12, 2014.
  * [38] M. Danielczuk, A. Mousavian, C. Eppner, and D. Fox, “Object rearrangement using learned implicit collision functions,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.


# GRBoondi: A code for evolving Generalized Proca theories on arbitrary
backgrounds

Shaun David Brocus Fell<EMAIL_ADDRESS>Institute for
Theoretical Physics, Universität Heidelberg , Philosophenweg 16, 69120
Heidelberg, Germany  Lavinia Heisenberg<EMAIL_ADDRESS>Institute for Theoretical Physics, Universität Heidelberg , Philosophenweg 16,
69120 Heidelberg, Germany

###### Abstract

While numerical simulations offer unparalleled precision and robustness in
studying complex physical systems, their execution is often hindered by
complexity, costliness, and time consumption due to the intricate equations
involved. This challenge is already encountered in General Relativity, where
non-flat spacetimes exacerbate the computational burden. This complexity is
further intensified when dealing with additional degrees of freedom. To
address this challenge head-on, we introduce GRBoondi , a groundbreaking
fixed-background numerical relativity code designed to provide a unified
interface for numerically solving Generalized Proca theories. GRBoondi grants
users the ability to make arbitrary modifications to the Proca equations of
motion on any background, providing a robust and versatile tool for exploring
diverse classes of Generalized Proca theories. This letter serves as part of
the submission of GRBoondi to the Journal of Open Source Software. For access
to the code, please visit https://github.com/ShaunFell/GRBoondi.git.

##  I Introduction

Generalized Proca theories [1, 2] provide a rich landscape to search for
solutions to deep, fundamental questions, such as the nature of dark energy
and dark matter. Moreover, strong gravity regimes, like those surrounding
dense, compact astrophysical entities, offer novel avenues to explore
fundamental fields. To effectively probe these domains, precise models are
imperative for sifting through vast quantities of data. Crafting concrete
models for the entirety of Generalized Proca theories represents a formidable
endeavor, typically relying on numerical methods. These numerical methods are
usually written to be case dependent, especially for a particular background.
This makes it difficult to generalize the computational code to account for
higher-order couplings or different types of backgrounds. The difficulty is
especially amplified when considering the full landscape of Generalized Proca
theories and beyond [3].

Generalized Proca theories represent the most general ghost-free theories for
a massive vector field that preserve the second-order equations of motion [1,
2, 4, 5, 6]. Being the most general construction, it naturally contains
numerous highly non-linear derivative interactions of the Proca field. Its
defining Lagrangian can be split into six separate interactions determined by
the order of the derivatives of the Proca field. Already at the order of the
quadratic Lagrangian, many interesting properties of non-linear Proca theories
become manifest [7, 8, 9, 10, 11, 12]. On the astrophysical side, interest in
Generalized Proca theories as dark matter and dark energy candidates has seen
recent growth [13, 14, 15, 16, 17, 18, 19, 20, 21]. In the pursuit of
investigating this class of theories within strong gravity regimes, numerical
simulations emerge as an indispensable tool. Additionally, it becomes
imperative to constrain the landscape of allowed interactions and parameters
associated with these theories.

GRBoondi 111Consistent with the GRChombo naming convention, Boondi is a Koori
word for a multi-purpose tool, used for hunting and digging. provides a
unified interface for computing the evolution of any Generalized Proca model
on an arbitrary background. Given a specific background with known expressions
for the metric variables and initial data for the Proca field, GRBoondi
numerically computes the time evolution of the Proca field, given user-
specified Generalized Proca equations of motion. While the simultaneous
evolution of the metric and matter fields offers the most comprehensive
depiction of their temporal progression, there are instances where the density
of the Proca field is negligible in comparison to the background curvature. In
such cases, fixed background evolution routines serve as an excellent
approximation to the complete mutual evolution. A core feature of GRBoondi is
that it simplifies many of the boiler-plate code required to begin a
simulation from scratch. The only additions a user needs to input are the
initial conditions, the modifications to the equations of motion (EOM), and
the background functions. GRBoondi will automatically compute various
diagnostic quantities and plot files. Moreover, GRBoondi is incredibly modular
and modifications to any of the in-built functions is effortless.

The core challenge in numerical computations of Generalized Proca theories
lies in formulating the EOM’s in a manner conducive to time evolution-based
solvers. This is usually done by performing a (3+1)-decomposition of the
equations of motion and solving for the evolution operator acting on the state
variables. For configuration space variables $\vec{Q}$, this amounts to
solving the EOM’s for the quantities $\mathfrak{L}_{m}\vec{Q}$, where
$\mathfrak{L}_{m}$ is the Lie derivative along the hypersurface evolution
vector, henceforth called the evolution operator. Once this quantity has been
solved for, it can be split into a time evolution piece and a hypersurface-
transverse piece,
$\mathfrak{L}_{m}\vec{Q}=\frac{1}{\alpha}\left(\dot{\vec{Q}}-\mathfrak{L}_{\beta}\vec{Q}\right)$,
where $\alpha$ is the lapse function and $\mathfrak{L}_{\beta}$ is the Lie
derivative along the shift vector. Once the evolution operator has been solved
for, the evolution equation can be easily recast in the form

$\dot{\vec{Q}}=F(\vec{Q},\vec{B})\;,$ (1)

where $F$ is some known function of all the other configuration space
variables and the background variables $\vec{B}$. For known expressions for
the background variables $\vec{B}$ and specifying $\vec{Q}$ to be the Proca
(3+1)-variables, GRBoondi numerically solves Eq.(1) using standard finite-
differencing techniques. For more specifics on the (3+1)-splitting, see for
example [22, 23].

The fixed background approximation provides a huge speedup compared to the
full evolution, since only a small fraction of the variables actually require
a time integration. Not only do simulations now require far less resources,
the speed up can be up to hundreds of times faster. GRBoondi is based on the
publicly available numerical relativity (NR) code GRChombo [24, 25], which
itself is based on the open source adaptive mesh refinement (AMR)-based
differential equation solver Chombo [26]. Some minor pieces of GRBoondi are
also based off the recently released open source software GRDzhadzha [27],
which is also based off of GRChombo.

##  II Statement of Need

In practice, any NR software library can be used to evolve Generalized Proca
theories. Prominent examples of NR codes encompass the extensive Einstein
toolkit [28] and its related Cactus framework [29], Kranc [30], LEAN [31], and
Canuda [32]. Other libraries worth mentioning are the non-public BAM [33],
AMSS-NCKU [34], PAMR [35] and HAD [36]. The non-exhaustive list can be
expanded with SPeC [37], which is a pseudo-spectral code that uses generalized
harmonic coordinates. Similarly, SpECTRE [38, 39, 40] uses the Galerkin
methods. NRPy [41] is a python library that aims for use on non-high
performance computing clusters. There are also cosmological simulation codes
like CosmoGRaPH [42] and GRAMSES [43]. Simflowny [44] is a magneto-
hydrodynamic simulation software. GRAthena++ [45] uses the oct-tree AMR for
maximum scaling. GRDzhadzha [27] is another fixed background code, based off
the GRChombo [24] framework.

This extensive list underscores the abundance of NR libraries at one’s
disposal. However, none of them provide a tailored unified interface for
studying the vast landscape of gravity theories, like Generalized Proca. Using
any of the numerous frameworks requires significant work in order to evolve
even a single Generalized Proca theory. In order to gauge the simplicity of
using GRBoondi , we compare the number of lines-of-code (LOC) needed to setup
a standard Proca simulation in GRChombo, GRDzhadzha, and GRBoondi . GRChombo
requires over 1500 lines of code, spread across 17 different files, including
source and header files. GRDzhadzha requires over 600 LOC across 14 files, and
GRBoondi requires just 300 LOC across 6 files. These simple metrics provide a
decent gauge of the complexity faced by the user for setting up a simple
standard Proca simulation. On top of this, GRBoondi offers catered plotting
routines for viewing data, leveraging the highly parallelizable VisIt [46]
analysis tool. In addition, since the metric variables and their derivatives
are computed exactly at each grid point, the adaptability of the AMR grid can
be focused solely on the matter variables.

While backreaction is neglected in GRBoondi , the error incured by this
approximation is estimated by computing the norm of the energy-momentum
tensor. This estimation takes into account the relativistic nature of the
matter field and gauges the error in the evolution at the level of the
Einstein equations.

Since GRBoondi is very similar to the GRChombo code, simulations performed
with GRBoondi can easily be ported to GRChombo, should full NR simulations be
required. This is particularly useful if the backreaction is found to be
significant at some point in the simulation. The data from a GRBoondi
simulation can be used as initial data for a GRChombo simulation, potentially
yielding much better initial data for the full NR simulation.

##  III Key Features of GRBoondi

  * •

Ease of use: A central pillar of GRBoondi is its ease of use, relative to
other NR software. Many of the basic boiler-plate code and complexities are
kept within the source code, allowing the user to have only a basic
understanding of the code in order to start researching their problem.

  * •

Arbitrary choice of background spacetime: The main parts of GRBoondi are
extremely modular, allowing for any arbitrary background to be plugged in,
even a numerically computed one. This allows for massive versatility in the
source code. Users can swap in and out backgrounds with ease. GRBoondi comes
pre-equipped with four background classes ready for use, along with testing
suites to verify convergence of each class. These include:

    * –

Minkowski space

    * –

Boosted Schwarzschild black hole

    * –

Kerr black hole

    * –

Kerr-de Sitter black hole

  * •

Arbitrary modifications to base Proca: GRBoondi uses specific coding idioms
that allow arbitrary modifications to the hard-coded equations of motion. This
allows for numerical simulations of any Generalized Proca theory. The base
model is the standard electromagnetic model, subject to the usual Proca
constraint. Any Generalized Proca theory can be added on top of this by simply
adding the additional pieces from the Generalized Proca Lagrangian. Examples
of simulating standard Proca and a simple non-linear model [47, 9, 12] are
included in the codes repository.

  * •

Accuracy: The metric values and their derivatives are computed exactly at each
point, while the matter variables are evolved using 4th-order Runge-Kutta time
integration and their spatial derivatives calculated using the same finite
difference stencils as GRChombo (currently, 4th- and 6th-order stencils are
available).

  * •

Various boundary conditions: Since GRBoondi inherits from GRChombo, all the
boundary conditions in GRChombo are also straight-forwardly applicable to
GRBoondi . GRBoondi thus supports:

    * –

Periodic boundary conditions

    * –

Sommerfeld out-going wave conditions

    * –

Static boundary conditions

    * –

Extrapolation at linear and zeroth order

    * –

Reflective symmetry (e.g. simulating only the top-half of the computational
box due to the background symmetry)

  * •

Checkpointing: Since GRBoondi inherits from Chombo, it inherits the
checkpointing feature. Simulations can be restarted from a checkpoint file
saved during the previous run, allowing for arbitrarily long simulations, even
on restricted computing clusters. Moreover, should a simulation fail for any
reason, it can be restarted using a previous checkpoint.

  * •

Diagnostics: GRBoondi comes equipped with many different diagnostic quantities
that can be calculated during the course of a simulation and additional tools
for users to calculate their own quantities.

    * –

_Computation of conserved quantities_. Computation of conserved quantities,
such as energy and angular momentum densities and their associated fluxes
across a surface can be toggled. GRBoondi can also calculate the energy-
momentum tensor trace and its square, for studying the importance of
backreaction. Some of the included diagnostic quantities are the square of the
Proca field, conserved energy density, Eulerian energy density, conserved
angular momentum, conserved energy density flux, conserved angular momentum
flux, trace of the energy-momentum tensor, and square of the energy-momentum
tensor. The conserved energy and angular momentum correspond to approximate
time and rotational killing vectors. These can be disabled if the spacetime
does not possess these.

    * –

_Excision of diagnostics in user-specified regions_. Currently, only spherical
excision is implemented, where the user specifies the minimum and maximum
radius for the diagnostic quantities. The user can specify which quantities
they wish to excise as well, as specified in a parameter file. Generalization
to more complicated diagnostic regions is case dependent and left up to the
user.

    * –

_Integration of quantities across the computational grid_. The user can
specify which quantities they wish to integrate.

    * –

_Integration of quantities across spherical surfaces_. The user can specify
which quantities they wish to integrate over a spherical surface and the radii
of those surfaces as well.

  * •

Tailored post-processing tools: GRBoondi comes equipped with custom-built
post-processing scripts, leveraging VisIt’s python interface. Any of the
diagnostic quantities computed during the evolution can be plotted using these
tools. It also contains a simple integral plotter, for quickly visualizing the
various integrals computed during a simulation, leveraging python’s Matplotlib
package [48].

Figure 1: Some examples of the usage of GRBoondi . The left image is the
energy density of a Proca cloud superradiantly excited around a rapidly
spinning black hole, superimposed with streamlines of the spatial Proca vector
field. The right image is a plot in the xz-plane of the square of the Proca
4-vector in the same background. The spin axis of the black hole is orientated
along the z-axis. Superimposed on both images is a slice of the computational
grid, clearly showing the refinement hierarchy.

##  IV Summary

GRBoondi provides a unified interface for computing the evolution of any
Generalized Proca model on an arbitrary background spacetime. It is a powerful
tool that opens up many avenues for possible research studies, including
cosmological, superradiant, modified gravity, and more. Moreover, having
GRChombo as a foundation allows its solutions to be easily ported as initial
conditions for a full Einstein+Generalized Proca evolution. We intend to
explore this avenue in a future follow-up project. Its modularity allows for
great user freedom, thus allowing GRBoondi to simulate the full landscape of
Generalized Proca theories. Its in-built post-processing routines allows for
quick visualization of simulation data, both during and after computations.
GRBoondi will thus be a vital tool in studying Generalized Proca theories,
opening the doors to probing astrophysically relevant models of dark matter
and dark energy.

##  V Acknowledgements

LH is supported by funding from the European Research Council (ERC) under the
European Unions Horizon 2020 research and innovation programme grant agreement
No 801781. LH further acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES
Excellence Cluster).

SF is indebted to Katy Clough for her stimulating discussions and extensive
input to the code. Katy was an invaluable source of aid for the nuances of NR
simulations and some of the finer details of GRChombo.

The simulations performed as part of this release were carried out on the
Baden-Württemberg high-performance computing cluster. The authors acknowledge
support by the state of Baden-Württemberg through bwHPC.

## References

  * Heisenberg [2019] L. Heisenberg, Physics Reports 796, 1–113 (2019), arXiv: 1807.01725.
  * Heisenberg [2014] L. Heisenberg, Journal of Cosmology and Astroparticle Physics 2014, 015–015 (2014), arXiv: 1402.7026.
  * Heisenberg  _et al._ [2016] L. Heisenberg, R. Kase, and S. Tsujikawa, Phys. Lett. B 760, 617 (2016), arXiv:1605.05565 [hep-th] .
  * Beltran Jimenez and Heisenberg [2016] J. Beltran Jimenez and L. Heisenberg, Phys. Lett. B 757, 405 (2016), arXiv:1602.03410 [hep-th] .
  * Allys  _et al._ [2016] E. Allys, P. Peter, and Y. Rodriguez, JCAP 02, 004 (2016), arXiv:1511.03101 [hep-th] .
  * Tasinato [2014] G. Tasinato, JHEP 04, 067 (2014), arXiv:1402.6450 [hep-th] .
  * Aoki and Minamitsuji [2022] K. Aoki and M. Minamitsuji, arXiv (2022).
  * Barausse  _et al._ [2022] E. Barausse, M. Bezares, M. Crisostomi, _et al._ , Journal of Cosmology and Astroparticle Physics 2022, 050 (2022).
  * Clough  _et al._ [2022] K. Clough, T. Helfer, H. Witek, _et al._ , Physical Review Letters 129 (2022), 10.1103/physrevlett.129.151102.
  * Coates and Ramazanoğlu [2022a] A. Coates and F. M. Ramazanoğlu, Physical Review Letters 129 (2022a), 10.1103/physrevlett.129.151103.
  * Coates and Ramazanoğlu [2023] A. Coates and F. M. Ramazanoğlu, Phys. Rev. D 107, 104036 (2023).
  * Ünlütürk  _et al._ [2023] K. I. Ünlütürk, A. Coates, and F. M. Ramazanoğlu, Physical Review D 108 (2023), 10.1103/physrevd.108.044022.
  * De Felice  _et al._ [2016a] A. De Felice, L. Heisenberg, R. Kase, _et al._ , JCAP 06, 048 (2016a), arXiv:1603.05806 [gr-qc] .
  * De Felice  _et al._ [2016b] A. De Felice, L. Heisenberg, R. Kase, _et al._ , Phys. Rev. D 94, 044024 (2016b), arXiv:1605.05066 [gr-qc] .
  * De Felice  _et al._ [2016c] A. De Felice, L. Heisenberg, R. Kase, _et al._ , Phys. Rev. D 93, 104016 (2016c), arXiv:1602.00371 [gr-qc] .
  * Heisenberg  _et al._ [2017a] L. Heisenberg, R. Kase, M. Minamitsuji, _et al._ , JCAP 08, 024 (2017a), arXiv:1706.05115 [gr-qc] .
  * Heisenberg  _et al._ [2017b] L. Heisenberg, R. Kase, M. Minamitsuji, _et al._ , Phys. Rev. D 96, 084049 (2017b), arXiv:1705.09662 [gr-qc] .
  * de Felice  _et al._ [2017] A. de Felice, L. Heisenberg, and S. Tsujikawa, Phys. Rev. D 95, 123540 (2017), arXiv:1703.09573 [astro-ph.CO] .
  * Geng  _et al._ [2021] C.-Q. Geng, Y.-T. Hsu, J.-R. Lu, _et al._ , Physics of the Dark Universe 32, 100819 (2021).
  * Gómez [2023] G. Gómez, Phys. Rev. D 107, 123535 (2023), arXiv:2202.07027 [gr-qc] .
  * de Rham  _et al._ [2022] C. de Rham, S. Garcia-Saenz, L. Heisenberg, _et al._ , Journal of Cosmology and Astroparticle Physics 2022, 053 (2022).
  * Zilhão  _et al._ [2015] M. Zilhão, H. Witek, and V. Cardoso, Classical and Quantum Gravity 32, 234003 (2015), arXiv:1505.00797 [gr-qc, physics:hep-th].
  * Gourgoulhon [2012] E. Gourgoulhon,  _3+1 Formalism in General Relativity_ , Lecture notes in physics (Springer, New York, NY, 2012).
  * Andrade  _et al._ [2021] T. Andrade, L. A. Salo, J. C. Aurrekoetxea, _et al._ , Journal of Open Source Software 6, 3703 (2021).
  * Clough  _et al._ [2015] K. Clough, P. Figueras, H. Finkel, _et al._ , Class. Quant. Grav. 32, 245011 (2015), arXiv:1503.03436 [gr-qc] .
  * Adams  _et al._ [2015] M. Adams _et al._ , (2015).
  * Aurrekoetxea  _et al._ [2024] J. C. Aurrekoetxea, J. Bamber, S. E. Brady, _et al._ , Journal of Open Source Software 9, 5956 (2024).
  * Cupp  _et al._ [2023] S. Cupp _et al._ , “The einstein toolkit,”  (2023).
  * Goodale  _et al._ [2003] T. Goodale, G. Allen, G. Lanfermann, _et al._ , in  _Vector and Parallel Processing – VECPAR’2002, 5th International Conference, Lecture Notes in Computer Science_ (Springer, Berlin, 2003).
  * [30] Kranc, “Kranc: Kranc assembles numerical code,” .
  * Sperhake [2007] U. Sperhake, Phys. Rev. D 76, 104015 (2007), arXiv:gr-qc/0606079 .
  * Witek  _et al._ [2021] H. Witek, M. Zilhao, G. Bozzola, _et al._ , “Canuda: a public numerical relativity library to probe fundamental physics,”  (2021).
  * Bruegmann  _et al._ [2008] B. Bruegmann, J. A. Gonzalez, M. Hannam, _et al._ , Phys. Rev. D 77, 024027 (2008), arXiv:gr-qc/0610128 .
  * Galaviz  _et al._ [2010] P. Galaviz, B. Bruegmann, and Z. Cao, Phys. Rev. D 82, 024005 (2010), arXiv:1004.1353 [gr-qc] .
  * East  _et al._ [2012] W. E. East, F. Pretorius, and B. C. Stephens, Phys. Rev. D 85, 124010 (2012), arXiv:1112.3094 [gr-qc] .
  * Neilsen  _et al._ [2007] D. Neilsen, E. W. Hirschmann, M. Anderson, _et al._ , in  _11th Marcel Grossmann Meeting on General Relativity_ (2007) pp. 1579–1581, arXiv:gr-qc/0702035 .
  * Pfeiffer  _et al._ [2003] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, _et al._ , Comput. Phys. Commun. 152, 253 (2003), arXiv:gr-qc/0202096 .
  * Deppe  _et al._ [2021] N. Deppe, W. Throwe, L. E. Kidder, _et al._ , “Spectre,”  (2021).
  * Kidder  _et al._ [2017] L. E. Kidder _et al._ , J. Comput. Phys. 335, 84 (2017), arXiv:1609.00098 [astro-ph.HE] .
  * Cao  _et al._ [2018] Z. Cao, P. Fu, L.-W. Ji, _et al._ , Int. J. Mod. Phys. D 28, 1950014 (2018).
  * Ruchlin  _et al._ [2018] I. Ruchlin, Z. B. Etienne, and T. W. Baumgarte, Phys. Rev. D 97, 064036 (2018), arXiv:1712.07658 [gr-qc] .
  * Mertens  _et al._ [2016] J. B. Mertens, J. T. Giblin, and G. D. Starkman, Phys. Rev. D 93, 124059 (2016), arXiv:1511.01106 [gr-qc] .
  * Barrera-Hinojosa and Li [2020] C. Barrera-Hinojosa and B. Li, JCAP 01, 007 (2020), arXiv:1905.08890 [astro-ph.CO] .
  * Palenzuela  _et al._ [2018] C. Palenzuela, B. Miñano, D. Viganò, _et al._ , Class. Quant. Grav. 35, 185007 (2018), arXiv:1806.04182 [physics.comp-ph] .
  * Daszuta  _et al._ [2021] B. Daszuta, F. Zappa, W. Cook, _et al._ , Astrophys. J. Supp. 257, 25 (2021), arXiv:2101.08289 [gr-qc] .
  * Childs  _et al._ [2012] H. Childs, E. Brugger, B. Whitlock, _et al._ , “High Performance Visualization–Enabling Extreme-Scale Scientific Insight,”  (2012).
  * Coates and Ramazanoğlu [2022b] A. Coates and F. M. Ramazanoğlu, “The intrinsic pathology of self-interacting vector fields,” (2022b), arXiv:2205.07784 [gr-qc] .
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).