# Pheno-Robot: An Auto-Digital Modelling System for In-Situ Phenotyping in the
Field

Yaoqiang Pan1,∗, Kewei Hu1,∗, Tianhao Liu2, Chao Chen2, and Hanwen Kang2,# ∗
Equal Contribution1 K.Hu and Y.Pan are with the College of Engineering, South
China Agriculture University, Guangzhou, China2 T.Liu, C.Chen, and H.Kang are
with the Department of Mechanical and Aerospace Engineering, Monash
University, Melbourne, Australia

###### Abstract

Accurate reconstruction of plant models for phenotyping analysis is critical
for optimising sustainable agricultural practices in precision agriculture.
Traditional laboratory-based phenotyping, while valuable, falls short of
understanding how plants grow under uncontrolled conditions. Robotic
technologies offer a promising avenue for large-scale, direct phenotyping in
real-world environments. This study explores the deployment of emerging
robotics and digital technology in plant phenotyping to improve performance
and efficiency. Three critical functional modules: environmental
understanding, robotic motion planning, and in-situ phenotyping, are
introduced to automate the entire process. Experimental results demonstrate
the effectiveness of the system in agricultural environments. The pheno-robot
system autonomously collects high-quality data by navigating around plants. In
addition, the in-situ modelling model reconstructs high-quality plant models
from the data collected by the robot. The developed robotic system shows high
efficiency and robustness, demonstrating its potential to advance plant
science in real-world agricultural environments.

##  I Introduction

Plant phenotyping plays a fundamental role in precision agriculture. It
involves identifying and selecting genetics with advantageous input traits and
complementary output traits [1]. While phenotyping in the laboratory plays a
key role in identifying promising lines for crossbreeding, they are a
surrogate for the primary goal of understanding how a crop will grow in real-
world environments [2]. Uncontrolled ’non-laboratory’ conditions present
significant challenges, particularly in the analysis of the traits responsible
for beneficial responses [3]. Therefore, the mass collection of phenotypic
data in the field is essential for precision agriculture.

At the heart of plant phenotyping is the digital modelling of plant growth and
traits, including both appearance and geometry [4]. This process involves
monitoring the growth, development, and changes in plants, taking into account
factors such as climate, soil properties, pests, and diseases [5]. Image-based
plant modelling is a widely used approach, particularly for analysing
observable traits [4]. However, analysing plants from a single viewpoint poses
challenges, especially when plants overlap [6]. To overcome this limitation,
RGB images acquired from multiple viewpoints are widely used. Recently, 3D
reconstruction technology has become a prominent tool for plant analysis [7].
This new paradigm provides essential information about the geometric aspects
of plant [8], including height [9], volume [10], and even mass [11]. These
advances in modelling technology, coupled with the increasing availability of
information-rich data, contribute significantly to the capture of detailed
plant characteristics.

Figure 1: Pheno-Robot system operates in the greenhouse. Figure 2: System
Overview of the Pheno-Robot system.

Despite significant advances, the current practice heavily relies on manual
operation, which is highly labour-intensive and inefficient, making it
impractical for large-scale farms. Robotics offers the potential for
widespread sampling in the field under authentic agricultural conditions.
However, the realisation of robotic automation for high-quality in-situ plant
phenotyping in the field faces two main challenges. Firstly, there is a lack
of effective methods for perception and autonomously navigating the
intricacies of farm landscapes to perform large-scale and direct phenotyping.
Secondly, there is a crucial gap in the availability of an accurate and
effective digital modelling method capable of generating high-fidelity and
multi-modal plant models. Successfully overcoming these challenges holds the
key to enabling capabilities for repeated and detailed assessments of plants,
potentially rendering a paradigm shift in the development of agri-genetics.

In this research, we present an innovative robotic system designed for
autonomous in-situ phenotyping in the field, synthesising both robotic and
digital technologies. The developed Pheno-Robot system comprises a
comprehensive framework that includes a deep learning model and a motion
planning method for robotic automation. It also incorporates a Neural Radiance
Field (NeRF)-based modelling network, which enhances the system’s ability to
perform detailed and accurate in-situ phenotyping. Our key contributions are:

  * •

develops a novel environmental understanding and the robotic navigation method
tailored for farm environments.

  * •

develops an improved NeRF method to address sparse-view input to achieve good
quality in plant modelling.

  * •

presents a hierarchy mapping method and demonstrates the Pheno-Robot system in
agricultural settings.

The rest of this paper is organised as follows. Section II surveys related
work, followed by the proposed methodologies in Section III. The experiment
results are given in Sections IV. and then the conclusions are given in
Section V.

##  II Related Works

###  II-A Robotic Automation in Agriculture

Robots play a key role in precision agriculture, using smart and intervention
technologies to improve efficiency through sensing and automation. Precision
agriculture addresses spatial and temporal variability in the environment and
plant growth patterns, scaling from traditional farm level [9] to sub-field
precision [12]. Horticulture farms are more complex than crop farms, requiring
robots to navigate semi-structured, challenging terrain in dynamic
environments [13]. Firstly, sensor information is essential for detecting
objects and potential risks in the field to ensure the safe operation of
robotic vehicles. Secondly, machine learning is indispensable for farming
robots to maximize locomotion flexibility (e.g., moving sideways or navigating
narrow spaces between crops) and optimize input utilization with greater
efficacy. The agility of these robots, coupled with the ability to carry
specialized sensors, holds the key to unlocking the full potential of
precision agriculture. Despite notable progress, exemplified by robots
designed for autonomous operations in orchard environments for selective
harvesting [14] or monitoring applications [15], there is still uncharted
territory. Robots capable of precise environmental recognition, robust and
flexible motion, and accurate plant modelling remain an area for exploration.

###  II-B In-situ phenotyping

Phenomics encompasses the study of various phenotypic plant traits, including
growth, yield, plant height, leaf area index and more. Traditional methods
rely on imaging and 3D range sensors to measure various plant traits such as
colour, shape, volume, and spatial structure [8]. For example, the growth rate
of rosette plants such as Arabidopsis [16] and tobacco [6] is often analysed
using images taken from a single viewpoint. On the other hand, 3D range
sensors allow accurate measurement of plant geometry and characteristics such
as height, width and volume [13]. Kang et al. [17] proposed a sensor fusion
system for yield estimation in apple orchards, using deep learning-based
panoptic segmentation algorithms. Wu et al. [2] utilised Multi-View Stereo
(MVS) for the reconstruction of crop geometry in the field. However, 3D range
sensors face limitations in agricultural environments due to strong occlusion,
often resulting in a bleeding effect [18] near discontinuous geometry. In
addition, their point resolution is sparse compared to image-based data. While
MVS-based methods offer better quality with sufficient data, they require
hours to process even simple plants [11], limiting their applicability in
real-time operations. Furthermore, current technologies do not provide a
multi-modal representation of instances, making them inefficient for analysing
plant features when multiple data types (image/geometry) are required.

##  III Methods

The Pheno-Robot system comprises three subsystems: an Environmental Perception
Module (EPM), a Motion Planning Module (MPM), and an In-situ Phenotyping
Module (IPM), as illustrated in Figure 2.

###  III-A Environment Perception

####  III-A1 3D Detection Network

In this study, we utilise a novel neural network, the 3D Object Detection
Network (3D-ODN), as introduced in our previous study [19]. Specifically
designed for processing point cloud maps in Bird’s-Eye View (BEV), the 3D-ODN
consists of three integral network branches: point cloud subdivision, feature
coding, and a detection head. First, the entire point cloud map is uniformly
segmentated. The points within each segment are then projected onto the BEV
perspective, where local semantic features are extracted. These features are
then fed into the main branch, forming a single-stage network architecture
dedicated to processing the point cloud features from the BEV angle. To
improve the learning of multi-scale feature embeddings, the feature pyramid is
used to facilitate the fusion of features of multi-scales. Finally, the
recognition results are projected into 3D space based on the predicted height
values. The identified instances are denoted as
$\mathbf{T}=\\{t_{1},t_{2},\cdots\rvert t_{i}\in\mathbb{R}^{2}\\}$,
representing the instances detected by the 3D-ODN (see Fig.3 (a)).

Figure 3: Illustration of EPM principle.

####  III-A2 Graph Mapping

For each instance, eight nodes are designated to represent the four corners in
two different directions of motion, as shown in Fig.3 (b). This orientation
consideration aims to mitigate turning movements in narrow passages, which
could otherwise lead to a higher failure rate of movement due to trapping into
potential obstacles. Consequently, each instance $o$ is associated with eight
nodes, denoted as $t_{i}=\\{N^{i}_{1},N^{i}_{2},\cdots,N^{i}_{8}\rvert
N^{i}_{j}\in\mathbb{R}^{2}\\}$ (see Fig.3 (b)). To extract the line structure
of farms, a line detection algorithm is used to group instances in the same
row, denoted $\mathbf{g}$. Each row has two common access nodes at both ends,
denoted $N_{l}$ and $N_{r}$ ($N_{l}$ and $N_{r}\in\mathbb{R}^{2}$). These
nodes serve to establish connections between all instances within the group
and to form links with other groups (see Fig.3 (a)). A group is thus
characterised as $g_{i}=\\{t_{1},t_{2},\cdots,N^{i}_{l},N^{i}_{r}\\}$. The map
is then represented as a set of groups denoted by
$\mathbf{G}=\\{g_{1},g_{2},\cdots\\}$.

###  III-B Trajectory Planning

Figure 4: Illustration of MPM principle.

####  III-B1 Global Path Planning

A Greedy-search-based path generation algorithm on the graph map is developed,
as illustrated in Alg.1. The following are the relevant definitions:

  * •

$\mathbf{T}$ represents the goal set that requires phenotyping.

  * •

$\mathbf{V}$ represents the subgroups if their instances in $\mathbf{T}$.

  * •

$\mathbf{N_{start}}$ represents the robot position at the start.

  * •

$\mathbf{\Gamma}$ represents the global path that includes a sequence of
nodes.

Algorithm 1 Global Path Generation

0: $\mathbf{T},\mathbf{G},\mathbf{N_{start}}$

0: $\mathbf{\Gamma}$

1: for $t_{i}$ in $\mathbf{T}$ do

2: $\mathbf{v}_{new}\leftarrow\mathbf{FindParent}(t_{i});$

3: $\mathbf{V}\leftarrow\mathbf{V}\bigcup\mathbf{v}_{new};$

4: end for

5: while $\mathbf{V}\neq\mathbf{\emptyset}$ do

6:
$\mathbf{v}_{j}\leftarrow\mathbf{FindNearestSubgroup}(\mathbf{\Gamma},\mathbf{V});$

7:
$\mathbf{\Gamma}_{new}\leftarrow\mathbf{PlanConnection}(\mathbf{\Gamma},\mathbf{v}_{j});$

8: $\mathbf{\Gamma}\leftarrow\mathbf{\Gamma}\bigcup\mathbf{\Gamma}_{new};$

9: $\mathbf{Delete}\ \mathbf{v}_{j}$ $\mathbf{from}\ \mathbf{V};$

10: end while

11: return $\mathbf{\Gamma}$

Two key sub-functions presented in Alg.1 are described as follows:

  * •

$\mathbf{FindParent}$: this function finds the group that this instance
belongs to.

  * •

$\mathbf{FindNearestSubgroup}$: this function finds the subgroup that has the
instance closest to the current robot location.

The key functions $\mathbf{PlanConnection}$ is described in Alg.2.

Algorithm 2 $\mathbf{PlanConnection}$

0: $\mathbf{\Gamma},\mathbf{v}=\\{t_{1},t_{2},\cdot,t_{n}\rvert
t_{i}\in\mathbf{g}_{i}\\}$

0: $\mathbf{\Gamma}_{new}$

1: while $\mathbf{v}\neq\mathbf{\emptyset}$ do

2: $\mathbf{N}_{new},\
t_{i}\leftarrow\mathbf{FindNearestAndFeasible}(\mathbf{\mathbf{\Gamma},v});$

3:
$\mathbf{\Gamma}_{new}\leftarrow\mathbf{\Gamma}_{new}\bigcup\mathbf{N}_{new};$

4: if $\mathbf{IFFullyCover}(\Gamma_{new},t_{i})$ then

5: $\mathbf{Delete}\ t_{i}$ $\mathbf{from}\ \mathbf{v};$

6: end if

7: end while

8: return $\mathbf{\Gamma}_{new}$

The function $\mathbf{FindNearestAndFeasible}$ finds the feasible and the
nearest node (see Fig.4 (a)) of an instance in the subgroup. It evaluates the
feasibility from two perspectives: traversability and orientation, which will
be detailed in Sec III-B3. Another function $\mathbf{IFFullyCover}$ evaluates
if an instance has been fully covered by the generated path. This is achieved
by checking the number of connected nodes of an instance, as shown in Fig.4
(b).

####  III-B2 Local Trajectory Generation

This step aims to compute the detailed path based on the global path,
including three steps: path generation, optimisation, and interpolation.

Initial path Generation: Given the global path $\Gamma$, the first step is to
generate an initial collision-free path between two nodes. Two scenarios are
considered, including path generation for phenotyping data acquisition (where
two adjacent nodes belong to the same instance) and paths between nodes of
different instances. In the first case, a sequence of collision-free preset
sample positions is established and the RRT algorithm is used to connect these
positions to form the sampling trajectory between nodes. In the second case,
the $\mathbb{A}^{*}$ algorithm is used to determine a collision-free path
between nodes.

Trajectory Optimisation: The initial path $\Phi$ consists of a sequence of
points $\\{Q_{0},Q_{1},\cdots|Q\in\mathbb{R}^{2}\\}$, and we conceptualise a
trajectory $t\in[0,1]\rightarrow\Phi\subset\mathbb{R}^{2}$ as a continuous
function mapping time to robot states. The objective function incorporates
three key aspects of the robot’s motion. It penalises velocities to encourage
smoothness, proximity to the environment to secure trajectories that maintain
a certain distance from obstacles, and the distance between the state and the
preset viewpoint. These terms are represented as $f_{s}$, $f_{c}$, and
$f_{o}$, respectively. The objective function is:

$\Phi^{*}=\mathop{\mathbf{argmin}}_{\Phi}\
\alpha_{s}f_{s}(\Phi)+\alpha_{c}f_{c}(\Phi)+\alpha_{o}f_{o}(\Phi)$ (1)

the initial state $\Phi(0)=Q_{0}$ and final state $\Phi(1)=Q_{1}$ are fixed.
$\alpha_{s}$, $\alpha_{c}$, and $\alpha_{o}$ are weights for each penalty
terms. The detailed formulations of $f_{s}$, $f_{c}$, and $f_{o}$ are as
follows:

$f_{s}(\Phi)=\int_{\Phi}||\frac{d}{dt}\Phi(t)||^{2}$ (2)
$f_{c}(\Phi)=\int_{\Phi}c(\Phi(t))||\frac{d}{dt}\Phi(t)||dt$ (3)
$f_{o}(\Phi)=\int_{\Phi}||\Phi(t)-\Phi_{d}||^{2}dt$ (4)

where $c(\cdot):\mathbb{R}^{2}\to\mathbb{R}$ be the function that penalises
the state near the obstacles, $\Phi_{d}$ is the trajectory that contains the
desired view-position for data acquisition.

We update the trajectory by the functional gradients $\bar{\nabla}f(\Phi_{i})$
using $\Phi_{i+1}=\Phi_{i}-l_{\mathbf{r}}\cdot\bar{\nabla}f(\Phi_{i})$,
following the work [20], where $l_{\mathbf{r}}$ is the learning rate. The
functional gradient of the objective in (2)-(4) is given by

$\bar{\nabla}f_{s}(\Phi)=-\frac{d^{2}}{dt^{2}}\Phi(t)$ (5)
$\bar{\nabla}f_{c}(\Phi)=||\Phi^{\prime}(t)||\cdot[(I-\Phi^{\prime}(t)\Phi^{\prime}(t)^{T})\nabla
c-c\kappa]$ (6) $\bar{\nabla}f_{o}(\Phi)=-(\Phi(t)-\Phi_{d})$ (7)

where $\nabla c$ is the derivative of obstacles to the control points $Q_{j}$
by the $\partial c/\partial Q_{j}$, the definition of $\kappa$ is given as
below.

$\kappa=||\Phi^{\prime}(t)||^{2}((I-\Phi^{\prime}(t)\Phi^{\prime}(t)^{T})\Phi^{\prime\prime})$
(8)

B-spline Interpolation: The optimised trajectory is then parameterised by a
piece-wise B-spline into a uniform curve $\Phi_{B}$. Given the determined
degree $m$, and a knot vector $\\{k_{0},k_{1},\cdot,k_{M}\\}$, where
$M=n_{q}+2m$. The parameterised uniform curve $\zeta(t)$ can be formulated as:

$\zeta(t)=\frac{\sum\limits_{i=1}^{n_{q}}R_{i,m}(t)w_{i}Q_{i}}{\sum\limits_{i=1}^{n_{q}}R_{i,m}(t)w_{i}}$
(9)

where $w_{i}$ is the weight of each control points $Q_{i}$, $R_{i,m}$ is the
basis function for $Q_{i}$ at $m$ degree, $u$ is the control value of the
curve. Then, the TEB-planner [21] is used to find the proper velocity to
follow the trajectory.

####  III-B3 Feasibility checking

The $\mathbf{FindNearestAndFeasible}$ function in Alg.2 evaluates the
feasibility between two nodes from two aspects: orientations and
traversability.

I. Orientation Checking: This term evaluates whether two nodes have similar
directions, as the narrow tunnels between plants do not allow turn-round. The
largest allowed orientation difference between two nodes is $1/3\pi$.

Figure 5: Illustration of terrain analysis, a) the points of robot
surroundings, b) the risk area is converted into polygons.

II. Traversability Analysis: The point cloud between two nodes is projected
onto a 2D grid map. The traversability of each grid is assessed as the
multiple risk factors, including:

(a) Collision risk: A risk factor quantified by the possibility of a grid
belonging to obstacles, denoted $\theta\in[0,1]$. A terrain analysis [22] is
utilised here.

(b) Slope risk: For each point on the terrain $T_{i}$, points on the map are
selected by a cube box with sides of length $l_{s}$, which is represented as
$\Omega_{i}=\\{(p^{j}_{i})_{j=1:N_{i}}\rvert p_{i}^{j}\in\mathbb{R}^{3}\\}$.
The SVD is used to fit a plane $P_{i}$ from the $\Omega_{i}$ and get the
normal vector $n_{z}\in\mathbb{R}^{3}$. The slop angle between terrain $T_{i}$
and vertical direction $n_{z}$ is obtained by
$\mathbf{s}=\mathbf{arccos}\frac{||e_{z}\cdot
n_{z}||}{||e_{z}||\cdot||n_{z}||}$.

(c) Step risk: This term evaluates the height gap between adjacent grid cells.
Negative obstacles can also be detected by checking the lack of measurement
points in a cell. The maximum height gap is denoted as $\lambda$.

The above three risks are combined by a weighted sum to obtain the weighted
traversability value $\Upsilon$, as:

$\Upsilon=\theta+\alpha_{\mathbf{s}}\frac{\mathbf{s}}{\mathbf{s}_{crit}}+\alpha_{\lambda}\frac{\lambda}{\lambda_{crit}}$
(10)

where the $\mathbf{s}_{crit}$ and $\lambda_{crit}$ are the maximum allowed
slope angle and height gap, respectively.

III. Terrain analysis for trajectory optimisation: For robot operation, the
odometry and KD-tree [23] are utilised to build the local terrain map (Fig.5
(a)). The detected obstacle points (Fig.5 (b)) are converted to polygons,
denoted as $\\{\mathcal{O}\subset\mathcal{R}^{2}\\}$. Let the
$\phi(Q_{i},\mathcal{O})$ describe the minimal Euclidean distance between
obstacles and a control point. A minimum separation $\phi_{\mathrm{min}}$
between all obstacles and $Q_{i}$ is found by

$\phi_{\mathrm{min}}=\mathrm{min}[\phi(Q_{i},\mathcal{O}_{1}),\phi(Q_{i},\mathcal{O}_{2}),\cdots]$
(11)

The obstacles $\mathcal{O}$ is updated online for dynamic environments.

###  III-C In-situ Phenotyping Model

####  III-C1 Neural Rendering

NeRF represent the 3D scene as a radiance field which describes volume density
$\sigma$ and view-dependent color $c$ for every point $x$ and every viewing
direction $d$ via a MLP [24]:

$\sigma,c=H_{\Theta}(x,d)$ (12)

The ray tracing-based volume rendering is used to render the parameters of a
3D scene into colours $\hat{C}$ of a ray $r$, expressed as:

$\displaystyle\hat{C}(r)$
$\displaystyle=\sum_{i=1}^{N}T_{i}\left(1-\exp\left(-\sigma_{i}\delta_{i}\right)\right)c_{i},T_{i}$
(13)

Where $T$ is the volume transmittance and $\delta$ is the step size of ray
marching. For every pixel, the squared loss on photometric error is used for
the optimization of MLP. When applied across the image, this loss
$\mathcal{L}_{\mathrm{Rendering}}$ is represented as:

$\mathcal{L}_{\mathrm{Rendering}}=\sum_{r\in R}||\hat{C}(r)-C(r)||_{2}^{2}$
(14)

where $C(r)$ is the ground truth colours.

Figure 6: Few-shot learning from sparse input. a) Ray samples under limited
training views. b) Occlusion regularization for sparse-view rendering.

####  III-C2 Few-shot learning from Sparse Views

The artefact of ”white floaters” caused by rendering distortion is the common
failure mode in NeRF learning. Autonomous data acquisition by robots can
always lead to imperfect density and sparse views with fewer overlapping
regions (Fig. 6 (a)), thus distorting the rendering of these regions. This is
essential because, between these sparse views, NeRF’s training lacks
sufficient information to estimate the correct geometric information of the
scene, leading to significant and dense floats in the region that is close to
the camera. To reduce those artefacts, an optimisation term
$\mathcal{L}_{occ}$, which is designed to regulate the learning behaviour of
NeRF [25], is used to penalize the dense fields near the camera via
“occlusion” regularization, which can be expressed as:

$\displaystyle\mathcal{L}_{occ}$
$\displaystyle=\frac{\boldsymbol{\sigma}_{K}\mathbf{m}_{K}}{K}=\frac{1}{K}\sum_{K}\sigma_{k}\cdot
m_{k}$ (15)

where $\boldsymbol{\sigma_{K}}$ represents the density values of the $K$
points sampled along the ray, ranked in order of proximity to the origin,
$\mathbf{m_{K}}$ is the binary mask vector that determines whether a ray
sector will be penalised or not.

####  III-C3 Geometry Extraction from Field

Given a predefined 3D region of interest, a set of spatial points
$P=\\{p_{1},p_{2},...,p_{n}\\}$ is generated via dense volumetric sampling.
For each point $p_{i}\in P$, it’s evaluated through the NeRF model to obtain
The density values, $\sigma(p_{i})=\mathrm{NeRF}_{\sigma}(p_{i})$, form the
basis for surface extraction. The Marching Cubes algorithm identifies the iso-
surface by threshold of the density values:

$M=\mathbf{MarchingCubes}(P,\sigma_{\text{threshold}})$ (16)

Where $M$ is the resultant mesh and $\sigma_{threshold}$ is an optimal density
value demarcating the object’s boundary.

For vertex $v_{j}$ in mesh $M$, a viewing ray $r_{j}$ is constructed and
queried $\mathbf{c}(v_{j})=\mathrm{NeRF}_{\mathbf{c}}(v_{j},r_{j})$ in NeRF.
Those values derived from the field are mapped onto mesh $M$, assigning colour
to each vertex:

$\operatorname{Color}(v_{j})=\mathbf{c}(v_{j})$ (17)

For a 2D texture representation, the vertex-coloured mesh undergoes UV
unwrapping. To minimize distortion, Least Squares Conformal Mapping (LSCM)
[26] is used, which is to minimize the conformal energy:

$E(u,v)=\int_{\Omega}\left(|\nabla u|^{2}+|\nabla v|^{2}\right)dA$ (18)

where $(u,v)$ are the 2D texture coordinates for each vertex in $M$, $\Omega$
represents the object’s surface, and $dA$ is a differential area element on
the mesh’s surface, indicating that the energy is computed by integrating over
the surface of the mesh.

####  III-C4 NeRF Model

The Instant-NGP [27], which makes use of a multi-resolution hashed encoding
that can represent learned features of the scenes with tiny MLPs, is utilised.
In detail, Instant-NGP operates on the premise that the object to be
reconstructed is enclosed within multi-resolution grids. For any point
$\mathbf{x}\in\mathbb{R}^{3}$ in various resolution grids, it obtains the hash
encoding $h^{i}(\mathbf{x})\in\mathbb{R}^{d}$, where $d$ is the features’
dimension, $i$ is the level of tri-linear interpolation. The hash encodings of
all levels are concatenated to form the multi-resolution feature
$h(\mathbf{x})=\\{h^{i}(\mathbf{x})\\}_{i=1}^{L}\in\mathbb{R}^{L\times d}$.

Training: The primary goal of our study is to produce high-quality rendering
models of instances. We utilise the rendering loss as it quantifies the
discrepancy between the rendered images and the input images, denoted as
${\mathcal{L}_{color}}$. Besides, to improve the MLP learning under sparse
view, we also apply ${\mathcal{L}_{occ}}$ from (12). Therefore, given a set of
posed-images $I_{gt}$ and the predicted renderings from the network
$I_{pred}$, the training loss ${\mathcal{L}_{MLP}}$ is defined as:

${\mathcal{L}_{MLP}}={\mathcal{L}_{color}}+\mathcal{L}_{occ}$ (19)

where
${\mathcal{L}_{color}}=(1/N)\cdot\sum_{i=1}^{N}\|I_{gt}^{(i)}-I_{pred}^{(i)}\|_{2}^{2}$
and $N$ is the total number of images in the datasets.

###  III-D Hierarchy Map Representation

A hierarchy scene representation that combines large-scale mapping with local
high-fidelity rendering is introduced, including structure-level and instance-
level.

I. Structure-level representation: The coarse-level map is designed to extract
high-level information, focusing on the structure of the scene and the plants
in the environment. In this context, we construct a map of the farm by
creating a coloured point cloud map and its corresponding semantic map.

II. Instance-level representation: Each plant in the environment will be
meticulously represented, ensuring that each plant $o_{i}$ is captured in fine
detail through the robot’s visual image stream. In the context of phenotyping,
each detected plant will be associated with a high-fidelity neural-learned
rendering model, accompanied by its corresponding 3D model.

##  IV Experiments

###  IV-A Pheno-Robot Hardware

The Pheno-Robot is equipped with a 32-line LiDAR, a Realsense D435 camera and
a 9-axis IMU (see Fig.7 (a) and (b)) for navigation. Communication between the
different modules is via the common message layer on the ROS server. The
sensors are connected to the computer via ROS-Noetic on Ubuntu 20.04. A 4K
GoPro Hero-11 is mounted on a gimbal stabilizer on the left side, and the
image stream from the GoPro is transmitted to the computer through the GoPro-
ROS-node. Plant modelling by NeRF training is performed remotely, with data
transmission via the 5G wireless network.

Figure 7: (a) and (b): Pheno-Robot systems; (c) and (d): the demonstration of
the EPM results in two test scenarios.

###  IV-B Evaluation on Pheno-Robot system

####  IV-B1 Evaluation on EPM

This section evaluates the performance of the EPM in typical agricultural
environments. Two different settings, shown in Fig. 7 (a) and (b), were
selected for the system evaluation, and their corresponding semantic maps are
shown in Fig. 7(c) and (d), respectively. The results show the precise
instance-level understanding of the EPM of the point cloud map in agricultural
scenarios, with an overall recall and precision for detection of 0.95 and 1.0,
respectively. The average position error and bounding box error are measured
to be 0.06m and 0.02m respectively. Particularly in agricultural landscapes
with row structures, our method demonstrates the ability to detect such
features and generate a graph map for robot autonomy. Furthermore, the method
is versatile and supports both online and offline modes depending on specific
requirements. For this study, we perform the semantic extraction process
offline.

####  IV-B2 Evaluation on MPM

This section evaluates the performance of the MPM for automated phenotyping by
robots. We conducted tests in two different environments, shown in Fig. 8 and
9, corresponding to the environments shown in Fig. 7 (a) and (b),
respectively.

Figure 8: Planning scene and motion planning results of the MPM in orchard-
like environments. Figure 9: Planning scene and motion planning results of the
MPM in greenhouses.

In the experiment, target plants requiring phenotyping were randomly selected.
The results show that the presented MPM achieves a high success rate and meets
the robot’s requirements in both scenarios, as shown in (a) of Fig. 8 and 9.
Compared to conventional $A^{*}$ or Dijkstra planners, which may generate
inappropriate global trajectories leading to turning movements in narrow
channels, causing local motion planning failures for car-like robots, our
graph-based global planner performs better. Beyond the global planner, the
developed local trajectory planner also shows strong robustness in field
environments, as illustrated in (b) and (c) of Figs. 8 and 9, respectively. A
significant factor affecting the performance of the local trajectory planner
is the traversability analysis. Given the complex terrain in agricultural
environments, traversability may be prone to overestimating or underestimating
risk areas during planning, leading to motion planning failures. In the
experiment, the robot was optimised with an additional solid-state lidar at
the front to improve its perception under complex conditions. In addition, a
replanning mechanism is employed for local motion planning, with an update
frequency of 5HZ and a forward planning distance of 10m. Overall, our system
demonstrates robust performance in the field with appropriate parameter
tuning. The maximum speed during sampling is 0.2m/s, while the maximum speed
for other conditions is 1m/s.

Figure 10: Demonstration of hierarchy maps of two scenarios in orchard-like
environments (A) and greenhouse (B). The a) of both A and B are global-level
maps of the scenes, and b) are the detailed models of instances in the
environments.

####  IV-B3 Evaluation on in-situ Phenotype

This section assesses the performance of IPM. We compare the in-situ
phenotyping models of both scenarios by using handheld data acquisition and
robot-automated acquisition. The quality of the results is evaluated by Peak
Signal to Noise Ratio (PSNR), as detailed in Table 1.

TABLE I: Comparison of Modelling quality |  | PSNR(dB) | Time(min)
---|---|---|---
Dataset | Methods | HA | RA | HA | RA
Outdoor | Instant-NGP | 24.8 | 22.5 | 2.4 | 3.2
Ours | 25.2 | 24.2 | 2.3 | 2.7
Greenhouse | Instant-NGP | 23.4 | 21.5 | 3.1 | 5.2
Ours | 23.3 | 22.7 | 2.3 | 4.1
|  |  |  |  |

In the initial experiments, both vanilla Instant-NGP and our model showed fast
convergence, taking less than 2 minutes when using hand-collected samples,
resulting in a PSNR above 24 dB and indicating high-quality results. However,
when using robot-collected samples, Instant-NGP took over 4 minutes to
converge and achieved a PSNR below 23 dB. In contrast, our method achieved a
PSNR above 23.5 dB and converged in 3 minutes, demonstrating the effectiveness
of occlusion regularisation in NeRF training with sparse-view inputs,
effectively mitigating geometric estimation errors in the renderings. Moving
on to the greenhouse experiments, the quality of modelling using both Instant-
NGP and our model showed a decrease, with PSNR values of around 23 dB in
approximately 3 minutes when using handheld collected samples. This reduction
in quality is mainly due to the complex geometry of the canopy leading to
occlusion. When modelling with robotically collected samples, both models
showed a decrease in model quality along with a longer training convergence
time. In comparison, our model performed better under these field conditions.
The hierarchy maps for both scenarios are shown in Figures 10 and
LABEL:fig:enter-label-2, revealing global structure-aware maps of farms and
detailed models for individual plants. These results highlight the ability of
the robotic system to achieve high-quality in-situ phenotyping in complex
agricultural environments, demonstrating resilience despite certain
limitations.

###  IV-C Demonstration of Digital-Modelling for Simulation

Our system can streamline the creation of virtual environments for robot
development using Isaac Sim, an advanced virtual environment known for
realistic physics simulation, rich sensor simulation, and graphical rendering.
The environment model, derived from terrain and trees of hierarchy maps, is
seamlessly integrated (see Fig. 11).

Figure 11: Isaac simulation environments of a) robot navigation and b) robot
harvesting.

The first scenario simulates a car-like robot equipped with LiDAR operating in
orchards (Fig. 11 (a)). The real-world terrain is converted into 2D images and
meshes are constructed based on these images. Tree locations are imported
directly from EPM predictions and tree models are generated by extracting
meshes from the renderings. In the second scenario (Fig. 11 (b) and (c)),
fruit models with fragile joints are added to the trees, creating a
reinforcement learning environment for robotic harvesting tasks. The
harvesting robot, equipped with a UR-5 manipulator, a soft gripper and a
camera on the remote robot base, is trained and evaluated in this virtual
environment.

##  V Conclusion

This study investigates the utilisation of robots in plant phenotyping to
increase efficiency and reduce labour-intensive tasks. The proposed system
consists of three key sub-systems that address environmental information
processing, motion planning for data acquisition, and modelling using data
collected by the robot. Experimental results demonstrate the effectiveness of
the system, particularly in outdoor environments with mild terrains, where the
robot can collect high-quality samples leading to superior plant phenotypic
models. For undulating terrains, typically challenging for plant phenotyping,
our enhancements to the NeRF model also exhibit promise in generating quality
plant phenotypic models. Future research will focus on further refining
phenotyping quality by exploiting the robotic system, incorporating the
robotic arm, and making additional improvements to the NeRF model to minimise
artefacts in challenging environments.

## References

  * [1] L. Fu, F. Gao, J. Wu, R. Li, M. Karkee, and Q. Zhang, “Application of consumer rgb-d cameras for fruit detection and localization in field: A critical review,” _Computers and Electronics in Agriculture_ , vol. 177, p. 105687, 2020.
  * [2] S. Wu, W. Wen, Y. Wang, J. Fan, C. Wang, W. Gou, and X. Guo, “Mvs-pheno: a portable and low-cost phenotyping platform for maize shoots using multiview stereo 3d reconstruction,” _Plant Phenomics_ , 2020.
  * [3] D. M. Deery and H. G. Jones, “Field phenomics: Will it enable crop improvement?” _Plant Phenomics_ , 2021.
  * [4] Y. Zhang and N. Zhang, “Imaging technologies for plant high-throughput phenotyping: a review,” _Frontiers of Agricultural Science and Engineering_ , vol. 5, no. 4, pp. 406–419, 2018.
  * [5] R. Xu and C. Li, “A modular agricultural robotic system (mars) for precision farming: Concept and implementation,” _Journal of Field Robotics_ , vol. 39, no. 4, pp. 387–409, 2022.
  * [6] B. Dellen, H. Scharr, and C. Torras, “Growth signatures of rosette plants from time-lapse video,” _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ , vol. 12, no. 6, pp. 1470–1478, 2015.
  * [7] L. Feng, S. Chen, C. Zhang, Y. Zhang, and Y. He, “A comprehensive review on recent applications of unmanned aerial vehicle remote sensing with various sensors for high-throughput plant phenotyping,” _Computers and electronics in agriculture_ , vol. 182, p. 106033, 2021.
  * [8] S. Paulus, “Measuring crops in 3d: using geometry for plant phenotyping,” _Plant methods_ , vol. 15, no. 1, pp. 1–13, 2019.
  * [9] N. Virlet, K. Sabermanesh, P. Sadeghi-Tehran, and M. J. Hawkesford, “Field scanalyzer: An automated robotic field phenotyping platform for detailed crop monitoring,” _Functional Plant Biology_ , vol. 44, no. 1, pp. 143–153, 2016\.
  * [10] S. Paulus, H. Schumann, H. Kuhlmann, and J. Léon, “High-precision laser scanning system for capturing 3d plant architecture and analysing growth of cereal plants,” _Biosystems Engineering_ , vol. 121, pp. 1–11, 2014.
  * [11] D. Xiong, D. Wang, X. Liu, S. Peng, J. Huang, and Y. Li, “Leaf density explains variation in leaf mass per area in rice between cultivars and nitrogen treatments,” _Annals of Botany_ , vol. 117, no. 6, pp. 963–971, 2016.
  * [12] M. S. A. Mahmud, M. S. Z. Abidin, A. A. Emmanuel, and H. S. Hasan, “Robotics and automation in agriculture: present and future applications,” _Applications of Modelling and Simulation_ , vol. 4, pp. 130–140, 2020.
  * [13] T. Duckett, S. Pearson, S. Blackmore, B. Grieve, W.-H. Chen, G. Cielniak, J. Cleaversmith, J. Dai, S. Davis, C. Fox _et al._ , “Agricultural robotics: the future of robotic agriculture,” _arXiv preprint arXiv:1806.06762_ , 2018.
  * [14] W. Au, H. Zhou, T. Liu, E. Kok, X. Wang, M. Wang, and C. Chen, “The monash apple retrieving system: a review on system intelligence and apple harvesting performance,” _Computers and Electronics in Agriculture_ , vol. 213, p. 108164, 2023.
  * [15] D. Albani, J. IJsselmuiden, R. Haken, and V. Trianni, “Monitoring and mapping with robot swarms for agricultural applications,” in _2017 14th IEEE International Conference on Advanced Video and Signal Based Surveillance (AVSS)_. IEEE, 2017, pp. 1–6.
  * [16] M. Jansen, F. Gilmer, B. Biskup, K. A. Nagel, U. Rascher, A. Fischbach, S. Briem, G. Dreissen, S. Tittmann, S. Braun _et al._ , “Simultaneous phenotyping of leaf growth and chlorophyll fluorescence via growscreen fluoro allows detection of stress tolerance in arabidopsis thaliana and other rosette plants,” _Functional Plant Biology_ , vol. 36, no. 11, pp. 902–914, 2009.
  * [17] H. Kang, H. Zhou, X. Wang, and C. Chen, “Real-time fruit recognition and grasping estimation for robotic apple harvesting,” _Sensors_ , vol. 20, no. 19, p. 5670, 2020.
  * [18] H. Kang, X. Wang, and C. Chen, “Accurate fruit localisation using high resolution lidar-camera fusion and instance segmentation,” _Computers and Electronics in Agriculture_ , vol. 203, p. 107450, 2022.
  * [19] Y. Pan, H. Cao, K. Hu, H. Kang, and X. Wang, “A novel perception and semantic mapping method for robot autonomy in orchards,” _arXiv e-prints, pp. arXiv–2308_ , 2023.
  * [20] M. Zucker, N. Ratliff, A. D. Dragan, M. Pivtoraiko, M. Klingensmith, C. M. Dellin, J. A. Bagnell, and S. S. Srinivasa, “Chomp: Covariant hamiltonian optimization for motion planning,” _The International journal of robotics research_ , vol. 32, no. 9-10, pp. 1164–1193, 2013.
  * [21] C. Rösmann, F. Hoffmann, and T. Bertram, “Kinodynamic trajectory optimization and control for car-like robots,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 5681–5686.
  * [22] W. Zhang, J. Qi, P. Wan, H. Wang, D. Xie, X. Wang, and G. Yan, “An easy-to-use airborne lidar data filtering method based on cloth simulation,” _Remote sensing_ , vol. 8, no. 6, p. 501, 2016.
  * [23] Y. Cai, W. Xu, and F. Zhang, “ikd-tree: An incremental kd tree for robotic applications,” _arXiv preprint arXiv:2102.10808_ , 2021.
  * [24] B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, and R. Ng, “Nerf: Representing scenes as neural radiance fields for view synthesis,” _Communications of the ACM_ , vol. 65, no. 1, pp. 99–106, 2021\.
  * [25] J. Yang, M. Pavone, and Y. Wang, “Freenerf: Improving few-shot neural rendering with free frequency regularization,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 8254–8263.
  * [26] B. Lévy, S. Petitjean, N. Ray, and J. Maillot, “Least squares conformal maps for automatic texture atlas generation,” in _Seminal Graphics Papers: Pushing the Boundaries, Volume 2_ , 2023, pp. 193–202.
  * [27] T. Müller, A. Evans, C. Schied, and A. Keller, “Instant neural graphics primitives with a multiresolution hash encoding,” _ACM Transactions on Graphics (ToG)_ , vol. 41, no. 4, pp. 1–15, 2022.


# The Initial Magnetic Field Distribution in AB Stars

Eoin Farrell Center for Computational Astrophysics, Flatiron Institute, New
York, NY 10010, USA  Adam S. Jermyn Center for Computational Astrophysics,
Flatiron Institute, New York, NY 10010, USA  Matteo Cantiello Center for
Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA
Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA  Daniel Foreman-Mackey Center for Computational Astrophysics,
Flatiron Institute, New York, NY 10010, USA

###### Abstract

Stars are born with magnetic fields, but the distribution of their initial
field strengths remains uncertain. We combine observations with theoretical
models of magnetic field evolution to infer the initial distribution of
magnetic fields for AB stars in the mass range $1.6-3.4M_{\odot}$. We tested a
variety of distributions with different shapes and found that a distribution
with a mean of $\sim$ 800 G and a full-width of $\sim$ 600 G is most
consistent with the observed fraction of strongly magnetized stars as a
function of mass. Our most-favoured distribution is a Gaussian with a mean
$\mu=770\,\mathrm{G}$ and standard deviation $\sigma=146\,\mathrm{G}$.
Independent approaches to measure the typical field strength suggest values
closer to 2 - 3 kG, a discrepancy which could suggest a mass-dependent and bi-
modal initial field distribution, or an alternative theoretical picture for
the origin of these magnetic fields.

Stellar physics (1621); Stellar evolutionary models (2046); Stellar convection
zones (301)

††journal: ApJ

##  1 Introduction

Magnetic fields play a significant role in many aspects of stellar evolution.
They transport angular momentum (Spruit, 2002a; Fuller et al., 2019a). They
affect stellar winds (Weber & Davis, 1967; Ud-Doula et al., 2009), which in
turn can alter planetary evolution (Vidotto et al., 2013). Magnetic fields can
also alter heat transport and produce star spots (Cantiello & Braithwaite,
2011), influence accretion (Bouvier et al., 2007), and enhance (Harrington &
Garaud, 2019) or inhibit (Gough & Tayler, 1966) chemical mixing.

Stars are expected to form with a range of initial magnetic field strengths.
Unfortunately this Initial Field Distribution (IFD) is difficult to directly
observe, as magnetic fields interact with convective regions and can be erased
by subsurface convection layers as a star evolves (e.g. Gough & Tayler, 1966;
Jermyn & Cantiello, 2020). Nonetheless, observations of early-type stars (OBA
stars) can probe the present-day distribution of magnetic field strengths, and
have revealed that it is _bimodal_. Stars are observed either with strong
fossil fields, in excess of $\sim 300\rm G$, or ultra-weak fields, with
amplitudes below a $\sim 3\rm G$, and these two modes appear to be mass-
dependent. Few or no stars are observed with field strengths in between
(Aurière et al., 2007; Grunhut et al., 2017; Fossati et al., 2015).

The bimodal field distribution suggests that there are at least two different
origin stories for early-type stellar magnetism. Strong magnetic fields may be
the relic of fossil field generation during star formation (e.g. Donati &
Landstreet, 2009) or during a stellar merger (Schneider et al., 2016, 2019a),
while weak magnetic fields could be the result of ongoing dynamo processes
(Cantiello & Braithwaite, 2011, 2019; Jermyn & Cantiello, 2021).

Exploring this scenario, Jermyn & Cantiello (2020) suggested that the
difference between the two modes could be explained by the strength of the
initial magnetic field in a star. Strong enough initial magnetic fields can
suppress near-surface convective motions and evolve just under flux
conservation, producing the strong-field mode of the present-day distribution.
Below some critical field strength, however, convection begins and can act to
erase the near-surface magnetic field, producing the observed weak-field mode.
This naturally explains both modes of the distribution as well as the
approximate (mass-dependent) field strengths of the magnetic desert between
the modes.

If this suggestion holds then whether or not a star will be observed in the
strong-field mode can be calculated given its initial magnetic field strength,
mass and evolutionary state. Hence, the fraction of strongly magnetized stars
as a function of mass should depend only on the IFD. Importantly, this
fraction has recently been observed in a volume-limited survey of AB stars by
Sikora et al. (2019a, b).

Our aim is to infer the Initial Field Distribution. We begin in Section 2 by
introducing the observational data (Section 2.1), our forward model of the
magnetic field evolution (Section 2.2) and the Approximate Bayesian
Computation approach (Section 2.3) that we use to infer the IFD from the
observations. We then apply these methods and obtain the IFD in Section 3. We
find that a distribution with a mean of $\sim$ 800 G and a full-width of
$\sim$ 600 G is most consistent with the observations, though we cannot
precisely infer the shape of the distribution. We explore uncertainties in our
analysis in Section 4 and discuss the implications of this IFD in Section 5.

##  2 Methods

###  2.1 Observations

Figure 1: Left: A Hertzsprung–Russell diagram with observed sample (Sikora et
al., 2019a, b) and our MESA stellar evolution models over-plotted with the
mass indicated at the beginning of each track. Non-magnetic AB stars are
plotted in grey and strongly-magnetized stars are shown in red. Right: The
critical magnetic field strength computed using Equation 1.

Sikora et al. (2019a) recently produced a volume-limited survey of all
identified intermediate mass MS stars within a heliocentric distance of 100 pc
as determined using Hipparcos parallaxes. Sikora et al. (2019b) used
spectropolarimetry to measure the magnetic field strengths of the magnetic
chemically peculiar stars (mCP) in the sample. Figure 1 shows their sample in
the Hertzsprung–Russell (HR) diagram.

We then choose to divide the data into five mass bins. Using less than five
bins limits the power of our inference method, while more than five leads to
limited sample sizes of just 1 – 2 stars in the upper mass bins. In each bin
we compute the fraction of strongly-magnetized stars. We identify objects as
strongly magnetic if Sikora et al. (2019b) were able to measure a magnetic
field and as weak otherwise, associating these with the two modes of the
present-day field distribution. Table 1 summarises the number of magnetic
stars nB, non-magnetic stars nnonB and fraction of magnetic stars fB =
nB/nnonB in the five mass bins. These five values of fB, along with the mass
bins, form the data which we use to infer the IFD.

Table 1: Characteristics of stars in the volume limited survey by Sikora et al. (2019a, b) binned by mass as derived from our stellar evolution models. Mass bin ($\mathrm{M}_{\odot}$) |  Number of mCP stars | Total Number | $f_{\rm B}$
---|---|---|---
1.44 - 1.87 | 10 | 1384 | 0.007
1.87 - 2.30 | 19 | 544 | 0.034
2.30 - 2.72 | 13 | 192 | 0.068
2.72 - 3.15 | 7 | 53 | 0.132
3.15 - 3.58 | 6 | 26 | 0.231

###  2.2 Theory

begin by noting that there is a critical vertical magnetic field strength
above which systems are _stable_ to convection. MacDonald & Petit (2019)
compute this as

$B_{\rm crit}^{2}=\frac{4\pi\rho c^{2}_{s}Q(\nabla_{\rm rad}-\nabla_{\rm
ad})}{1-Q(\nabla_{\rm rad}-\nabla_{\rm ad})+d\ln\Gamma_{1}/d\ln p},$ (1)

where $c_{s}$ is the sound speed, $\Gamma_{1}$ is the first adiabatic index,
$\nabla_{\rm ad}$ is the adiabatic temperature gradient, $\nabla_{\rm rad}$ is
the radiative temperature gradient, $\rho$ is the density, and

$Q=\frac{4-3p_{\rm gas}/p}{p_{\rm gas}/p}.$ (2)

Here $p_{\rm gas}$ is the gas pressure and $p$ is the total pressure.

Using this criterion, we model the evolution of the magnetic fields assuming
the following:

  1. 1.

The magnetic field is uniform near the surface of a star.

  2. 2.

So long as $B>B_{\rm crit}$, convection is suppressed and the field evolves
according to flux conservation ($B\propto R^{-2}$, where $R$ is the stellar
radius).

  3. 3.

If the magnetic field is ever $B<B_{\rm crit}$, near-surface convective
motions begin and quickly erase the surface magnetic field.

We calculated stellar evolution tracks for stars ranging from
$1.6-3.4\,M_{\odot}$ using revision 15140 of the Modules for Experiments in
Stellar Astrophysics (MESA, Paxton et al., 2011, 2013, 2015, 2018, 2019)
software instrument. Details on the MESA microphysics inputs are provided in
Appendix A. We swept the mass range in increments of $0.2M_{\odot}$ and used
solar metallicity ($Z=0.02$). The evolutionary tracks of our models are
plotted in the HR diagram in Fig. 1. In these calculations we forced the
temperature gradient to be $\nabla_{\rm rad}$ in the outer envelope where
subsurface convective layers would otherwise form. This allows us to simulate
the inhibition of convection by a strong magnetic field in those regions. This
choice has a very small effect on the surface properties, as shown by Jermyn &
Cantiello (2020). However, it does impact the thermodynamic properties of the
subsurface convection regions, which are relevant for computing $B_{\rm
crit}$.

We evaluate $B_{\rm crit}$ for each model at each point in time using equation
(1) and taking the maximum value inside the subsurface convective layers,
because in our scenario even a small convective region is enough to erase the
surface magnetic field. For the stars observed by (Sikora et al., 2019a) the
relevant subsurface convective regions are driven by hydrogen (H CZ) and
helium ionization (HeI CZ, HeII CZ) (Cantiello & Braithwaite, 2019). We note
that the HeI CZ is frequently stabilized by viscosity and thermal diffusion
(Jermyn et al., 2022), so formally we should exclude its contribution to
$B_{\rm crit}$. However, our models indicate that the HeI CZ never produces
the largest $B_{\rm crit}$ among available convection zones.

Fig. 2 (upper) shows $B_{\rm crit}$ as a function of mass and evolutionary
stage in our models. Similar to Jermyn & Cantiello (2020), we find that the
critical field strength decreases with increasing stellar mass. We also find
that forcing the temperature gradient to the radiative gradient in the
subsurface convective layers typically decreases $B_{\rm crit}$ by about 1 %
relative to what Jermyn & Cantiello (2020) calculated, though this can
increase to 10% for low-mass models at towards the late main-sequence.

Stars expand during the main-sequence (MS). Assuming conservation of magnetic
flux, the surface magnetic field $B$ decrease with radius $R$ as $B\propto
1/R^{2}$. Therefore the initial field strength required to produce a strong
(super-critical) magnetic field for a given mass and MS fractional age is the
cumulative maximum of $B_{\rm crit}$ from the upper panel in Fig. 2 scaled by
the square of the increase of the radius. This _initial_ $B_{\rm crit}$ is
plotted in lower panel of Fig. 2, where the increase as a function of MS age
is due to the expansion of the star. As in the upper panel, the critical field
strength decreases with increasing mass.

Figure 2: Upper panel: Critical magnetic field required to prevent the
formation of subsurface convective layers as a function of normalised main
sequence age for models of different masses. Lower panel: The critical initial
magnetic field required to suppress the formation of subsurface convective
layers and maintain a strong magnetic field until a given age as a function of
normalised main sequence age. The trend in age is due to expansion during the
main sequence.

###  2.3 Approximate Bayesian Computation

We use an Approximate Bayesian Computation (ABC) method (e.g. Alsing et al.,
2018; Sunnåker et al., 2013) to perform likelihood-free inference of the
initial field distribution. The ABC method allows us to reconstruct the IFD as
a probability density function $I(B)=dN/dB$ from the observed fraction of
strongly magnetised stars f${}_{\rm B}(M_{\star})$ using only forward
simulations, free from any likelihood assumptions or approximations. The
advantage of using ABC is that it is relatively easy to to simulate mock data
based on a given IFD, but it is not feasible to analytically propagate an IFD
to an expected distribution of observed fields. We use the elfi python package
(Lintusaari et al., 2018) for ABC to perform our simulations, which take place
as follows.

First, we choose a parameterised, functional form for $I(B)$. We write this as
$I_{\mu,\sigma,...}(B)$, where the subscripts denote the parameterization. For
simplicity, we assume that the fossil field distribution is mass independent
throughout. We discuss the implications of this assumption later on.

Next, we choose an appropriate prior distribution $p(\mu,\sigma,...)$ over
these free parameters. For example, one of our parameterizations is a Gaussian
distribution described by a mean $\mu$ and a standard deviation $\sigma$, and
we choose flat priors over $\mu$ and $\sigma$.

Finally, we enter a loop:

  * •

We sample values for the free parameters of our distribution from the prior.

  * •

We construct a population of stars with initial magnetic field strengths
sampling the resulting IFD ($I_{...}(B)$). These stars have masses and ages
chosen to match those in the observed volume-limited sample described above.

  * •

For each star, we use our stellar evolution models and frozen flux evolution
to determine if the magnetic field is weaker than the critical field at any
point between the pre-main-sequence and the age of the star inferred from
observations. If this occurs, we mark that star as ‘not magnetic’, assuming
that convection will erase the fossil field. Otherwise we mark the star as
‘magnetic’.

  * •

We bin the simulated stars by mass and compute the magnetic fraction $f_{\rm
B}$ in each bin.

  * •

We compare $f_{\rm B}$ to that reported from observations by computing the sum
of the absolute differences between the fractions in each mass bin.

  * •

If this difference is less than some threshold distance, the values of the
parameters chosen at the beginning are saved, otherwise they are rejected.

We perform this loop 105 times for each functional form of $I(B)$. The
accepted samples form our posterior distribution over the free parameters of
each functional form. We verified that the accepted posterior distributions
did not change substantially when we increased the number of iterations in the
loop or the threshold distance between the simulations and observations

##  3 Inferred Initial Field Distribution

Table 2: Details of a range of distributions that we tested, including the free parameters, prior distributions, means and standard deviations ($\sigma$) of the posterior distributions and the ‘distance’ between the threshold 1000th best simulation and the observations (D1000). Distribution | Parameters | Priors | Posterior Mean | Posterior Standard Deviation | D1000
---|---|---|---|---|---
Standard Assumptions
Linear | m (slope) | U(-1, 1) | -0.12 | 0.01 | 0.39
Trapezoidal | m (slope) | U(-1, 1) | -0.086 | 0.48 | 0.07
Bmin | U(300, 850) | 551 | 106
Bmax | U(850, 1350) | 1028 | 53
Gaussian | $\mu$ | U(500, 1000)) | 770 | 44 | 0.06
$\sigma$ | U(10, 230) | 146 | 37
Triangular | Midpoint | U(500, 1000) | 767 | 131 | 0.05
Bmin | U(300, 600) | 445 | 85
Bmax | U(800, 1300) | 1107 | 79
Additional Tests
Gaussian + $f_{0}$ | $\mu$ | U(500, 1000) | 843 | 64 | 0.07
$\sigma$ | U(10, 230) | 115 | 49
f0 | U(0, 1) | 0.34 | 0.20
Gaussian + | $\mu_{0}$ | U(500, 1000) | 731 | 172 | 0.03
mass dependence | $\sigma$ | U(10, 230) | 150 | 47
($\mu=\mu_{0}+\beta*m$) | $\beta$ | U(-100, 100) | 14 | 52

Because magnetohydrodynamics is complicated, we do not have a strong physical
prior for the functional form of the IFD. For simplicity, we begin by assuming
that all stars are born with some non-zero magnetic field and that the
distribution is independent of stellar mass. We discuss the impact of relaxing
these assumptions later on. We tested the following different forms for the
IFD, all of which are summarized in the upper half of Table 2.

  1. 1.

Linear: A distribution over $B\in[300,1500]\mathrm{G}$ of the form
$I(B)=I_{0}+mB$, where $I_{0}$ is determined to normalize the distribution. We
assign the slope $m$ a uniform prior distribution over $[-1,1]$.

  2. 2.

Trapezoidal: A variant of the linear distribution but with variable left and
right endpoints, resulting in two additional free parameters ($B_{\rm min}$,
$B_{\rm max}$), which we assign broad uniform priors.

  3. 3.

Gaussian: A Gaussian distribution characterized by a mean $\mu$ and variance
$\sigma^{2}$. Both $\mu$ and $\sigma$ are given uniform priors.

  4. 4.

Triangular: A distribution with a single peak, dropping to zero linearly on
either side. Both the peak location and the left and right roots are allowed
to vary, with all three having uniform priors.

Figure 3: Comparing the results from the ABC simulations for different
functional forms of the fossil field distribution. Left panel: Plots of the
fossil field distributions with the free parameters equal to the means of
their posterior distributions. Middle panel: Comparison between outputs of
simulations and observations. Right panel: The distribution of the “distance”
between the observations and simulations for each distribution type.

To properly compare these different forms we performed our ABC procedure over
$10^{5}$ distributions of each form, sampling the distribution parameters from
their priors. We chose our acceptance tolerance separately for each form so as
to accept the best $10^{3}$ samples, so that the posterior would be well-
sampled, and used these to construct the posterior distributions for each
parameter. We verified that the posterior distributions we find are not
sensitive to this acceptance threshold. Corner plots of these posteriors are
provided in Appendix B.

The left panel of Fig. 3 shows examples of all four forms, with their
respective parameters set to the means of their respective posterior
distributions. The middle panel then compares the distributions from the 100
best simulations for each functional form with the observations. Finally, the
right panel shows the distribution of the “distance” between the observations
and simulations for each distribution type. The ‘distance’ is the sum of the
distance between each of the observational data points in the middle panel and
the corresponding simulated points.

We find that a linear distribution from 300 G to 1500 G is inconsistent with
the observed data (Fig. 3, middle). The distribution of distances in the lower
panel of Fig. 3 confirms that linear distributions perform poorly compared
with the other distributions. The reason for this is that the linear
distribution produces many more strongly-magnetized $\sim
1.8\mathrm{M}_{\odot}$ stars than are actually observed. By contrast, the
trapezoidal, Gaussian and triangular forms are all able to reproduce the
observations within a much tighter tolerance, and visually reproduce the
observed trend in $f_{\rm B}$ much more closely (Fig. 3, middle).

The trapezoidal form is a useful one for building intuition about our
inference procedure. This form favours a relatively uniform distribution from
500 G to 1100 G. Because the critical magnetic field strength decreases with
increasing stellar mass, the lower bound of the distribution is sensitive to
higher mass stars and the upper bound of the distribution is sensitive to
lower mass stars. The upper bound of the trapezoidal distribution is
relatively strongly constrained due to the very low fraction of strongly
magnetised $\sim\,1.8\mathrm{M}_{\odot}$ stars. The lower bound and slope are
less constrained because the fraction of strongly magnetised stars at higher
masses just tells us how many stars have field strengths below the lowest
$B_{\rm crit}$, not how those stars are distributed.

The Gaussian and triangular forms are very similar in shape, and both rather
different from the trapezoidal form. They both favour a peak around $800\rm G$
and a similar width as that of the trapezoidal distribution (e.g. $\sigma\sim
130\rm G$). This suggests that the initial field distribution could have a
peak at 800 G and a full-width of around 600 G, with very few stars born with
magnetic fields of less than 500 G or greater than 1100 G.

Given that Gaussian, triangular, and trapezoidal all fit the data similarly
well we cannot tell much else about the shape of the distribution. At the same
time, the fact that all three forms obtain similar IFDs is encouraging, and
tells us that features which all three agree on, like the mean and the width
of the IFD, are likely robust.

The trapezoidal and triangular forms, with their sharp kinks and
discontinuities, seem less likely to be found in nature than the Gaussian
Therefore, given the similarity between their performance and that of the
Gaussian, we choose to explore the Gaussian form in more detail below.

Figure 4: The black line indicates the IFD using the means from the posterior
distributions. A selection of samples from the posterior distribution over
Gaussian IFDs are shown in grey.

The posterior distribution for the Gaussian form has mean $\mu=770\pm
44\,\mathrm{G}$ and standard deviation $\sigma=145\pm 77\,\mathrm{G}$. These
are anti-correlated, such that we see that distributions with higher mean tend
to be narrower, and those with higher variance tend to have lower means. Fig.
4 shows samples of distributions drawn from this posterior, as well as the
average over these samples, which broadly tells the same story as above: the
initial field distribution is peaked around 800 G and has a full-width of
around 600 G, with very few stars born below 500 G or above 1100 G.

Figure 5: Posterior distributions for the mean $\mu$ and the standard
deviation $\sigma$ are shown for the Gaussian form.

##  4 Uncertainties

We now consider various uncertainties in our approach.

###  4.1 Bimodal IFDs

We have assumed that all stars were born with some non-zero initial magnetic
field. While this is almost certainly the case, it could well be that the IFD
is itself bimodal, and observations of young stars suggest the possibility of
a second mode at very weak field strengths (Villebrun et al., 2019).

Figure 6: Posterior distribution of the mean $\mu$ and standard deviation
$\sigma$ and fraction of stars born with no strong magnetic field ($f_{0}$)
for a Gaussian initial field distribution.

To explore this possibility we investigated distributions in which some
fraction of stars are born with an initial magnetic field of $B_{\rm init}=0$,
while the rest have initial fields distributed according to a Gaussian form.
We refer to this as the Gaussian+$f_{0}$ form. Fig. 6 shows the posterior
distribution over $f_{0}$, $\mu$, and $\sigma$ for this form. We see that
while $\sigma$ is mostly uncorrelated with $f_{0}$, $\mu$ and $f_{0}$ are
strongly correlated. This makes sense: as the fraction of stars with zero
initial field strength increases, more stars with strong magnetic fields are
needed to match the observed magnetic fraction $f_{\rm B}$. For this reason,
regardless of the precise form we choose, any functional form allowing an
additional mode at weak field strengths will exhibit such a degeneracy. We
cannot exclude such forms, and this degeneracy introduces an additional
uncertainty to our models.

###  4.2 Possible Mass Dependence

We investigated the possibility of a dependence of the distribution of initial
fossil field on mass. To do this, we tested a Gaussian IFD with a mean that
scales linearly with the stellar mass $m$ such that $\mu=\mu_{0}+\beta*m$. We
found that $\mu_{0}$ and $\beta$ are strongly degenerate (Fig. 7).

To understand this, note that our approach is primarily using the fact that
$B_{\rm crit}$ varies with mass, and so differences in the magnetic fraction
between different mass bins tell us about the IFD. As a result if we permit
our parameterized IFD to vary with mass we effectively lose our ability to
constrain the distribution. What little inferential power we retain comes from
the evolution of $B_{\rm crit}$ across the main sequence, but we do not have a
large enough sample of stars to also bin by age and so this is only minimally
constraining. As a result, we emphasize that our results are strongly
contingent on the assumption of a mass-independent IFD over the range from
$1.6-3.4\mathrm{M}_{\odot}$.

Figure 7: Posterior distribution of IFD parameters for a Gaussian form and
allowing a dependence of the mean of the Gaussian with the mass through
$\mu=\mu_{0}+\beta*m$.

###  4.3 Sensitivity to overshooting

To study the sensitivity of our results to the details of our stellar
evolution calculations, we computed a grid of models with a higher value of
core convective overshooting ($f_{\rm ov}=0.020$ compared to $f_{\rm
ov}=0.014$). Convective overshooting affects the size of stellar cores, and
consequently the stellar luminosity and lifetime. We find that the posterior
distributions of the parameters inferred using stellar models with larger
overshooting is consistent with the distribution inferred using the original
models to within 1%.

###  4.4 Magnetic Field Evolution

We have assumed that, so long as the magnetic field is strong enough to shut
off convection, the field strength evolves according to flux conservation from
the beginning of the main sequence phase onward. This neglects a variety of
different possible phenomena. On the main sequence the most important are
magnetic diffusion, differential rotation, and stellar interactions.

We think it is safe to neglect magnetic diffusion, as the diffusion time
across a star is long compared with its main-sequence lifetime (e.g.
Braithwaite & Spruit, 2017). We likewise suspect that differential rotation
can be neglected, as super-critical magnetic fields are likely strong enough
to inhibit any significant shears from developing on the main-sequence.

Stellar interactions, on the other hand, do pose a challenge to our analysis.
Mass accretion can potentially enhance surface stellar magnetic fields
(Grunhut et al., 2013), and very strongly magnetized stars are believed to be
the outcome of stellar mergers (Schneider et al., 2019b). Depending on the
incidence of stellar interactions on the main sequence, these processes could
impact our ability to infer the IFD at the ZAMS. While the expected fraction
of massive stars undergoing interactions during the main sequence is fairly
high ($\sim$ 30%, de Mink et al., 2014), it is much smaller for the mass range
explored here. This is because the fraction of short period binaries is
smaller in these later-type stars (Moe & Di Stefano, 2017). In the future it
would be useful to carefully quantify the impact of stellar interactions on
the magnetization of main-sequence, intermediate mass stars.

##  5 Discussion

We have inferred the initial distribution of magnetic field strengths (IFD)
from a volume-limited sample of AB stars with magnetic field measurements. To
do this we have assumed that the IFD is independent of stellar mass and that
convection erases near-surface subcritical magnetic fields. Our favoured
distribution is a Gaussian with mean $\mu=770\,\mathrm{G}$ and standard
deviation $\sigma=145\,\mathrm{G}$. We now turn to the astrophysical
implications of this distribution.

###  5.1 Magnetic Desert

Jermyn & Cantiello (2020) proposed that the observed distribution of magnetic
fields - and the “magnetic desert” – can be explained in a scenario where some
fraction of stars are born with a smooth distribution of initial fields and
these fields are then modified by convection. The fact that our preferred IFD
reproduces the observed magnetic fraction in AB stars is consistent with this
hypothesis.

It is important to note that this did not need to be the case: we could not
have successfully fit arbitrary data. For instance this story would have
trouble reproducing a non-monotone magnetic fraction as a function of mass,
and would need to invoke more complicated assumptions such as a mass-dependent
IFD.

###  5.2 Connection to Star Formation

It is unclear what physical mechanism sets the IFD. It is possible that the
IFD reflects some frozen-in flux present in the molecular cloud a star formed
from. However, that story runs into trouble explaining the mean of the IFD.
Typical molecular clouds are have strong enough magnetic fields that they
cannot collapse without shedding magnetic energy (Troland & Crutcher, 2008).
If just enough magnetic field is shed to allow the cloud to collapse then we
should expect initial magnetic fields in approximate equipartition with gas
pressure, predicting field strengths which are orders of magnitude stronger
than what is observed.

On the other hand if convection indeed erases any magnetic fields that came
before then we would expect the magnetic field on the ZAMS to be generated by
a convective dynamo as the star descends the Hayashi track. Villebrun et al.
(2019) compiled measurements of magnetic fields on the pre-main-sequence in
their figure 1. For fully convective stars at the base of the Hayashi track
they see a broad distribution of magnetic field strengths ranging from $\sim
100\,\rm G$ to $\sim 800\,\rm G$. These field strengths are in equipartition
with the kinetic energy of convection in the outer layers of a star carrying
luminosity $L\sim L_{\odot}$, and occur for stars with luminosities of that
order, so it is at least plausible that the these fields are dynamo-generated.

From that point on though the story becomes unclear. These relatively strong
fields persist through the development of a radiative interior, though they
mostly vanish by the time early-type stars develop convective cores and
radiative envelopes (Villebrun et al., 2019), so there could well be non-
trivial field evolution happening in the immediate evolution prior to the
ZAMS.

###  5.3 Solar Magnetism

If the IFD we have inferred holds at lower masses, we should expect the Sun to
have been born with a $\sim 800\,\rm G$ magnetic field, and that field could
well have survived below the solar convection zone to this day (Cowling,
1953).

So far as we know, the fossil field in the interior of the Sun has not been
measured. There are, however, various upper bounds which have been obtained.
Reasoning from measurements of the solar oblateness, Friedland & Gruzinov
(2004) found an upper bound of $7\times 10^{6}\,\rm G$, which is certainly
consistent with our prediction. By contrast Boruta (1996) inferred an upper
bound of $\sim 30\,\mathrm{G}$ in the radiative zone based on the lack of bias
in the solar cycle, contingent on some assumptions including that the dynamo
is localized to the tachocline. If their assumptions hold then this tight
limit is evidence against our IFD, and favours more complicated models like
the bimodal Gaussian+$f_{0}$ forms we considered in Section 4.

The results of Stello et al. (2016) also support the idea that the interior of
stars like the Sun is not strongly magnetized. They used asteroseismology to
infer the presence of strong magnetic fields in evolving red giant cores
(Fuller et al., 2019a), and claimed that the incidence of core magnetization
goes to zero for stars below 1.1$M_{\odot}$. These stars are the descendants
of main sequence stars with radiative cores, while our inference used stars
with convective cores on the main sequence. The interaction of convective
cores with a large scale magnetic field could be an important piece of the
puzzle (e.g. Featherstone et al., 2009).

###  5.4 Compact Remnants

Our results predict that some early-type stars with subcritical surface
magnetic fields may have strong fossil fields hiding just beneath their
subsurface convection zones (Jermyn & Cantiello, 2021). Assuming simple flux
conservation , a field of $\sim$ 800 G would give rise to fields of order
$10^{6}\,\rm G$ in red giant cores and white dwarfs (WD) and $10^{15}\,\rm G$
in neutron stars (NS) (e.g. Woltjer, 1964; Duncan & Thompson, 1992).
Asteroseismology has revealed magnetic fields of order $10^{6}\,\rm G$ in
about 20% of red giant cores (Fuller et al., 2015; Stello et al., 2016), and
roughly 20% of white dwarfs have $10^{5-6}\,\rm G$ fields (Landstreet &
Bagnulo, 2019). This does not necessarily mean that all strongly magnetized NS
(magnetars) and magnetic WD inherit their fields directly from the IFD. Indeed
there are other mechanisms for generating or destroying magnetic fields in the
course of stellar evolution (Duncan & Thompson, 1992; Spruit, 2002b; Fuller et
al., 2019b), but the scenario discussed here provides one possible, simple
mechanism for the generation of (strongly) magnetized compact remnants.

The Flatiron Institute is supported by the Simons Foundation. We are grateful
to Gregg Wade for helpful conversations on the IFD.

##  Appendix A MESA

The MESA EOS is a blend of the OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon et
al., 1995), FreeEOS (Irwin, 2004), HELM (Timmes & Swesty, 2000), and PC
(Potekhin & Chabrier, 2010) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime by Buchler & Yueh (1976).
Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al., 2010) plus
additional tabulated weak reaction rates Fuller et al. (1985); Oda et al.
(1994); Langanke & Martínez-Pinedo (2000). Screening is included via the
prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from
Itoh et al. (1996).

The MESA input files are available as a ZENODO repository: Farrell et al.
(2022).

##  Appendix B Posterior distributions

Figs. 8 and 9 show the posterior distributions of the free parameters of the
IFD assuming respectively trapezoidal and triangular functional forms.

Figure 8: Posterior distribution of the slope $m$, the minimum B field value
and maximum B field value for the trapezoidal initial field distribution.
Figure 9: Posterior distribution of the location of the peak of the
distribution, the minimum B field value and maximum B field value for the
triangular initial field distribution.

## References

  * Alsing et al. (2018) Alsing, J., Wandelt, B., & Feeney, S. 2018, MNRAS, 477, 2874, doi: 10.1093/mnras/sty819
  * Aurière et al. (2007) Aurière, M., Wade, G. A., Silvester, J., et al. 2007, A&A, 475, 1053, doi: 10.1051/0004-6361:20078189
  * Boruta (1996) Boruta, N. 1996, ApJ, 458, 832, doi: 10.1086/176861
  * Bouvier et al. (2007) Bouvier, J., Alencar, S. H. P., Boutelier, T., et al. 2007, A&A, 463, 1017, doi: 10.1051/0004-6361:20066021
  * Braithwaite & Spruit (2017) Braithwaite, J., & Spruit, H. C. 2017, Royal Society Open Science, 4, 160271, doi: 10.1098/rsos.160271
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440, doi: 10.1086/154847
  * Cantiello & Braithwaite (2011) Cantiello, M., & Braithwaite, J. 2011, A&A, 534, A140, doi: 10.1051/0004-6361/201117512
  * Cantiello & Braithwaite (2019) —. 2019, ApJ, 883, 106, doi: 10.3847/1538-4357/ab3924
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028, doi: 10.1103/PhysRevD.76.025028
  * Cowling (1953) Cowling, T. G. 1953, in The Sun, ed. G. P. Kuiper, 532
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240, doi: 10.1088/0067-0049/189/1/240
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7, doi: 10.1088/0004-637X/782/1/7
  * Donati & Landstreet (2009) Donati, J. F., & Landstreet, J. D. 2009, ARA&A, 47, 333, doi: 10.1146/annurev-astro-082708-101833
  * Duncan & Thompson (1992) Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9, doi: 10.1086/186413
  * Farrell et al. (2022) Farrell, E., Jermyn, A., Cantiello, M., & Foreman-Mackey, D. 2022, Supporting data for ”The Initial Magnetic Field Distribution in AB Stars”, Zenodo, doi: 10.5281/zenodo.6139141
  * Featherstone et al. (2009) Featherstone, N. A., Browning, M. K., Brun, A. S., & Toomre, J. 2009, ApJ, 705, 1000, doi: 10.1088/0004-637X/705/1/1000
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fossati et al. (2015) Fossati, L., Castro, N., Schöller, M., et al. 2015, A&A, 582, A45, doi: 10.1051/0004-6361/201526725
  * Friedland & Gruzinov (2004) Friedland, A., & Gruzinov, A. 2004, ApJ, 601, 570, doi: 10.1086/380480
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1, doi: 10.1086/163208
  * Fuller et al. (2015) Fuller, J., Cantiello, M., Stello, D., Garcia, R. A., & Bildsten, L. 2015, Science, 350, 423, doi: 10.1126/science.aac6933
  * Fuller et al. (2019a) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019a, MNRAS, 485, 3661, doi: 10.1093/mnras/stz514
  * Fuller et al. (2019b) —. 2019b, MNRAS, 485, 3661, doi: 10.1093/mnras/stz514
  * Gough & Tayler (1966) Gough, D. O., & Tayler, R. J. 1966, MNRAS, 133, 85, doi: 10.1093/mnras/133.1.85
  * Grunhut et al. (2013) Grunhut, J. H., Wade, G. A., Leutenegger, M., et al. 2013, MNRAS, 428, 1686, doi: 10.1093/mnras/sts153
  * Grunhut et al. (2017) Grunhut, J. H., Wade, G. A., Neiner, C., et al. 2017, MNRAS, 465, 2432, doi: 10.1093/mnras/stw2743
  * Harrington & Garaud (2019) Harrington, P. Z., & Garaud, P. 2019, ApJ, 870, L5, doi: 10.3847/2041-8213/aaf812
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Irwin (2004) Irwin, A. W. 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors.  http://freeeos.sourceforge.net/
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411, doi: 10.1086/192264
  * Jermyn et al. (2022) Jermyn, A. S., Anders, E. H., & Cantiello, M. 2022, arXiv e-prints, arXiv:2201.10567.  https://arxiv.org/abs/2201.10567
  * Jermyn & Cantiello (2020) Jermyn, A. S., & Cantiello, M. 2020, ApJ, 900, 113, doi: 10.3847/1538-4357/ab9e70
  * Jermyn & Cantiello (2021) —. 2021, ApJ, 923, 104, doi: 10.3847/1538-4357/ac2d2a
  * Landstreet & Bagnulo (2019) Landstreet, J. D., & Bagnulo, S. 2019, A&A, 628, A1, doi: 10.1051/0004-6361/201936009
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481, doi: 10.1016/S0375-9474(00)00131-7
  * Lintusaari et al. (2018) Lintusaari, J., Vuollekoski, H., Kangasrääsiö, A., et al. 2018, Journal of Machine Learning Research, 19, 1.  http://jmlr.org/papers/v19/17-374.html
  * MacDonald & Petit (2019) MacDonald, J., & Petit, V. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3904, doi: 10.1093/mnras/stz1545
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231, doi: 10.1006/adnd.1994.1007
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Schneider et al. (2019a) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019a, Nature, 574, 211, doi: 10.1038/s41586-019-1621-5
  * Schneider et al. (2019b) —. 2019b, Nature, 574, 211, doi: 10.1038/s41586-019-1621-5
  * Schneider et al. (2016) Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355, doi: 10.1093/mnras/stw148
  * Sikora et al. (2019a) Sikora, J., Wade, G. A., Power, J., & Neiner, C. 2019a, MNRAS, 483, 2300, doi: 10.1093/mnras/sty3105
  * Sikora et al. (2019b) —. 2019b, MNRAS, 483, 3127, doi: 10.1093/mnras/sty2895
  * Spruit (2002a) Spruit, H. C. 2002a, A&A, 381, 923, doi: 10.1051/0004-6361:20011465
  * Spruit (2002b) —. 2002b, A&A, 381, 923, doi: 10.1051/0004-6361:20011465
  * Stello et al. (2016) Stello, D., Cantiello, M., Fuller, J., et al. 2016, Nature, 529, 364, doi: 10.1038/nature16171
  * Sunnåker et al. (2013) Sunnåker, M., Busetto, A. G., Numminen, E., et al. 2013, PLOS Computational Biology, 9, 1, doi: 10.1371/journal.pcbi.1002803
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Troland & Crutcher (2008) Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457, doi: 10.1086/587546
  * Ud-Doula et al. (2009) Ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2009, MNRAS, 392, 1022, doi: 10.1111/j.1365-2966.2008.14134.x
  * Vidotto et al. (2013) Vidotto, A. A., Jardine, M., Morin, J., et al. 2013, A&A, 557, A67, doi: 10.1051/0004-6361/201321504
  * Villebrun et al. (2019) Villebrun, F., Alecian, E., Hussain, G., et al. 2019, A&A, 622, A72, doi: 10.1051/0004-6361/201833545
  * Weber & Davis (1967) Weber, E. J., & Davis, Leverett, J. 1967, ApJ, 148, 217, doi: 10.1086/149138
  * Woltjer (1964) Woltjer, L. 1964, ApJ, 140, 1309, doi: 10.1086/148028


# On Interval Non-Edge-Colorable Eulerian Multigraphs

Petros A. Petrosyan

Department of Informatics and Applied Mathematics,
Yerevan State University, 0025, Armenia
Institute for Informatics and Automation Problems,
National Academy of Sciences, 0014, Armenia
E-mail<EMAIL_ADDRESS>
###### Abstract

An edge-coloring of a multigraph $G$ with colors $1,\ldots,t$ is called an
interval $t$-coloring if all colors are used, and the colors of edges incident
to any vertex of $G$ are distinct and form an interval of integers. In this
note, we show that all Eulerian multigraphs with an odd number of edges have
no interval coloring. We also give some methods for constructing of interval
non-edge-colorable Eulerian multigraphs.

##  1 Introduction

In this note we consider graphs which are finite, undirected, and have no
loops or multiple edges and multigraphs which may contain multiple edges but
no loops. Let $V(G)$ and $E(G)$ denote the sets of vertices and edges of a
multigraph $G$, respectively. The degree of a vertex $v\in V(G)$ is denoted by
$d_{G}(v)$, the maximum degree of $G$ by $\Delta(G)$, and the chromatic index
of $G$ by $\chi^{\prime}\left(G\right)$. A multigraph $G$ is Eulerian if it
has a closed trail containing every edge of $G$. The terms and concepts that
we do not define can be found in [7].

A proper edge-coloring of a multigraph $G$ is a coloring of the edges of $G$
such that no two adjacent edges receive the same color. If $\alpha$ is a
proper edge-coloring of $G$ and $v\in V(G)$, then $S\left(v,\alpha\right)$
denotes the set of colors of edges incident to $v$. A proper edge-coloring of
a multigraph $G$ with colors $1,\ldots,t$ is called an interval $t$-coloring
if all colors are used, and for any vertex $v$ of $G$, the set
$S\left(v,\alpha\right)$ is an interval of integers. A multigraph $G$ is
interval colorable if it has an interval $t$-coloring for some positive
integer $t$. The set of all interval colorable multigraphs is denoted by
$\mathfrak{N}$.

The concept of interval edge-coloring of multigraphs was introduced by
Asratian and Kamalian [1]. In [1, 2], they proved the following result.

Theorem 1. If $G$ is a multigraph and $G\in\mathfrak{N}$, then
$\chi^{\prime}\left(G\right)=\Delta(G)$. Moreover, if $G$ is a regular
multigraph, then $G\in\mathfrak{N}$ if and only if
$\chi^{\prime}\left(G\right)=\Delta(G)$.

Some results on interval edge-colorings of multigraphs were obtained in [5].
In [6], the authors described some methods for constructing of interval non-
edge-colorable bipartite graphs and multigraphs.

In this note we show that all Eulerian multigraphs with an odd number of edges
have no interval coloring. We also give some methods for constructing of
interval non-edge-colorable Eulerian multigraphs.

##  2 Results

Let $G$ be a multigraph. For any $e\in E(G)$, by $G_{e}$ we denote the
multigraph obtained from $G$ by subdividing the edge $e$. For a multigraph
$G$, we define a multigraph $G^{\star}$ as follows:

$V(G^{\star})=V(G)\cup\\{u\\}$, $u\notin V(G)$,

$V(G^{\star})=E(G)\cup\\{uv:v\in V(G)~{}and~{}d_{G}(v)~{}is~{}odd\\}$.

For a graph $G$, by $L(G)$ we denote the line graph of the graph $G$.

We also need a classical result on Eulerian multigraphs.

Euler’s Theorem. ([3]) A connected multigraph $G$ is Eulerian if and only if
every vertex of $G$ has an even degree.

Now we can prove our result.

Theorem 2. If $G$ is an Eulerian multigraph and $|E(G)|$ is odd, then
$G\notin\mathfrak{N}$.

Proof Suppose, to the contrary, that $G$ has an interval $t$-coloring $\alpha$
for some $t$. Since $G$ is an Eulerian multigraph, $G$ is connected and
$d_{G}(v)$ is even for any $v\in V(G)$, by Euler’s Theorem. Since $\alpha$ is
an interval coloring and all degrees of vertices of $G$ are even, we have that
for any $v\in V(G)$, the set $S\left(v,\alpha\right)$ contains exactly
$\frac{d_{G}(v)}{2}$ even colors and $\frac{d_{G}(v)}{2}$ odd colors. Now let
$m_{odd}$ be the number of odd colors in the coloring $\alpha$. By Handshaking
lemma, we obtain $m_{odd}=\frac{1}{2}\sum\limits_{v\in
V(G)}\frac{d_{G}(v)}{2}=\frac{|E(G)|}{2}$. Thus $|E(G)|$ is even, which is a
contradiction. $\square$

Corollary 1. If $G$ is an Eulerian multigraph and $G\in\mathfrak{N}$, then
$|E(G)|$ is even.

Let us note that there are Eulerian graphs with an even number of edges that
have no interval coloring. For example, the complete graph $K_{5}$ has no
interval coloring. On the other hand, there are many Eulerian graphs with an
even number of edges that have an interval coloring. In [4], Jaeger proved the
following result.

Theorem 3. If $G$ is a connected $r$-regular graph ($r\geq 2$),
$\chi^{\prime}\left(G\right)=r$ and $|E(G)|$ is even, then
$\chi^{\prime}\left(L(G)\right)=2r-2$.

Since $G$ is a connected $r$-regular graph ($r\geq 2$) and $|E(G)|$ is even,
we have that $L(G)$ is a connected $(2r-2)$-regular graph with an even number
of edges. Moreover, by Theorems 1 and 3 and Euler’s Theorem, we obtain the
following

Corollary 2. If $G$ is a connected $r$-regular ($r\geq 2$) graph with an even
number of edges and $G\in\mathfrak{N}$, then $L(G)$ is an Eulerian graph with
an even number of edges and $L(G)\in\mathfrak{N}$.

Let us note that Theorem 2 also gives some methods for constructing of
interval non-edge-colorable Eulerian multigraphs from interval colorable
multigraphs.

Corollary 3. If $G$ is an Eulerian multigraph and $G\in\mathfrak{N}$, then for
each $e\in E(G)$, $G_{e}\notin\mathfrak{N}$.

Corollary 4. If $G$ is a connected multigraph with an odd number of edges and
$G\in\mathfrak{N}$, then $G^{\star}\notin\mathfrak{N}$.

## References

  * [1] A.S. Asratian, R.R. Kamalian, Interval colorings of edges of a multigraph, Appl. Math. 5 (1987) 25-34 (in Russian).
  * [2] A.S. Asratian, R.R. Kamalian, Investigation on interval edge-colorings of graphs, J. Combin. Theory Ser. B 62 (1994) 34-43.
  * [3] L. Euler, Solutio problematis ad geometriam situs pertinentis, Commentarii Academiae Sci. I. Petropolitanae 8 (1736) 128-140.
  * [4] F. Jaeger, Sur l’indice chromatique du graphe representatif des aretes d’un graphe regulier, Discrete Math. 9 (1974) 161-172.
  * [5] R.R. Kamalian, Interval edge-colorings of graphs, Doctoral Thesis, Novosibirsk, 1990.
  * [6] P.A. Petrosyan, H.H. Khachatrian, Interval non-edge-colorable bipartite graphs and multigraphs, Journal of Graph Theory, 2013, http://onlinelibrary.wiley.com/doi/10.1002/jgt.21759/pdf
  * [7] D.B. West, Introduction to Graph Theory, Prentice-Hall, New Jersey, 1996.