# Negative Temperature in Spin Dynamics Simulations

Pui-Wai Ma<EMAIL_ADDRESS>United Kingdom Atomic Energy Authority,
Culham Science Centre, Abingdon, Oxfordshire, OX14 3DB, United Kingdom  Sergei
L. Dudarev<EMAIL_ADDRESS><EMAIL_ADDRESS>United
Kingdom Atomic Energy Authority, Culham Science Centre, Abingdon, Oxfordshire,
OX14 3DB, United Kingdom

###### Abstract

A simple and computationally efficient algorithm enables implementing negative
temperature values in a spin dynamics simulation. The algorithm uses a
Langevin spin dynamics thermostat with a negative damping parameter, enabling
the thermalization of an arbitrary interacting spin system to the Gibbs energy
distribution with a given negative temperature value. Canonical spin dynamics
simulations at a negative temperature are as robust as conventional positive
spin temperature simulations, providing a tool for quantitative dynamic
studies of the physics of highly excited magnetic states. Two simulation case
studies describing spin systems with antiferromagnetic and ferromagnetic
ground states are explored. The phase transitions occurring in the negative
temperature range do not necessarily exhibit similarities with their positive
temperature counterparts. The transition temperatures and the character of
spin alignment vary depending on the spatial range and strength of spin-spin
interactions.

Negative temperature is a well established concept in statistical mechanics
[1, 2, 3, 4, 5]. The occurrence of the negative temperature phenomenon was
first demonstrated experimentally by Purcell and Pound [6] for the nuclear
spin sybsystem of a LiF crystal. Similar studies were later performed using
other materials [7, 8, 9, 10], achieving negative temperatures on the nano-
Kelvin scale.

Negative temperature is a general notion not restricted to a spin ensemble.
Negative temperatures are encountered in the context of plasma states [11] and
optical lattices [12, 13]. A negative temperature value is associated with the
formation of an ordered low entropy high energy state in an ensemble where the
spectrum of energy eigenstates is bounded from above [1, 2, 3, 4, 5]. From the
thermodynamic definition of temperature, $1/T=\partial S(E)/\partial E$, where
$S(E)$ is the entropy of a macroscopic system at energy $E$, negative
temperature corresponds to a range of statistical configurations where entropy
$S$ is a decreasing function of energy $E$. At the point where the derivative
of $S$ changes sign and function $S(E)$ changes from monotonically increasing
to monotonically decreasing, we find that $1/T=0$, or in other words, $T$
changes from $T=\infty^{+}$ to $T=\infty^{-}$. A system with negative
temperature has a higher energy than a system at $T=\infty^{+}$, hence
necessarily requiring that the spectrum of energy eigenstates is bounded from
above. The highest attainable total energy state corresponds to $T=0^{-}$.
Negative temperatures are not encountered in atomistic molecular dynamics
simulations, since for particles moving in real space with kinetic energy
$m{\bf v}^{2}/2$, velocities are not bounded and a higher average energy $E$
implies higher positive temperature.

If a system is at equilibrium, the probability $W_{n}$ of finding it in an
eigenstate $\\{n\\}$ with energy $E_{n}$ is given by the Gibbs distribution

$W_{n}={1\over{\cal Z}}\exp\left(-\frac{E_{n}}{k_{B}T}\right),$ (1)

where ${\cal Z}$ is the partition function ${\cal
Z}=\sum_{n}\exp\left(-{E_{n}}/{k_{B}T}\right)$. Equation (1) imposes no
constraint on the sign of $T$, which can be positive or negative, as long as
the partition function can be computed. Monte Carlo simulations provide a
convenient means for studying statistical mechanical properties of systems
accessible to simulations at negative temperature, but such studies are
limited [11].

In a particular example of nuclear spins [6, 7, 8, 9, 10], since the strength
of spin-spin interactions is low, a negative temperature state can be created
by applying a strong external magnetic field $\mathbf{H}_{ext}$. In this
strong magnetic field limit, the nuclear spin Hamiltonian is dominated by the
Zeeman term ${\cal H}_{Z}=-\mathbf{I}\cdot\mathbf{H}_{ext}$, where
$\mathbf{I}$ is the nuclear spin vector. In accord with Eq. (1), once the
direction of external field is instantaneously reversed, a metastable short-
lived population of inverted spin orientations is formed, characterized by a
negative temperature value.

The external field reversal technique is what has so far been used for
generating negative temperature states, both experimentally and
computationally [6, 7, 8, 9, 10, 12]. The problem with this technique is that
temperature remains undefined during the equilibration process, when the
system dynamically relaxes from the initial highly excited state into a state
of thermodynamic equilibrium at a positive temperature. Besides, the technique
only applies to systems that can be manipulated by macroscopic spatially
homogeneous external fields. For example, a spin system with a
ferromagnetically ordered ground state cannot be transformed into a negative
temperature configuration using an applied external magnetic field technique
since the external field would simply change the overall orientation of the
spins in response to any attempt to change the spin order.

In this study, we propose an algorithm that enables performing a controlled
computationally efficient thermalization of a dynamic spin ensemble to a
negative temperature. The algorithm, involving the use of a Langevin spin
thermostat developed earlier [14, 15, 16], enables performing a detailed
quantitative exploration of equilibrium physics at negative temperatures. The
approach is general and applies to any interacting or non-interacting dynamic
statistical mechanics spin ensemble. Dynamic spin systems exhibiting anti-
ferromagnetic and ferromagnetic ground states are explored, providing
computational examples where the Curie-Weiss phase transition law is observed
in the negative temperature range.

Consider a system of interacting spins described by a general type Hamiltonian
${\cal H}(\mathbf{S}_{1},\mathbf{S}_{2},\ldots,\mathbf{S}_{N})$. The
corresponding Langevin equations of motion for individual spins
$\mathbf{S}_{i}$ are [17, 18, 19, 14, 16]

$\frac{d\mathbf{S}_{i}}{dt}=\frac{1}{\hbar}\left[\mathbf{S}_{i}\times\left(\mathbf{H}_{i}+\mathbf{h}_{i}\right)-\gamma\mathbf{S}_{i}\times\left(\mathbf{S}_{i}\times\mathbf{H}_{i}\right)\right],$
(2)

where $\mathbf{H}_{i}=-\partial{\cal H}/\partial\mathbf{S}_{i}$ is the
effective field acting on $i$-th spin, $\gamma$ is a damping parameter and
$\mathbf{h}_{i}$ is a delta-correlated fluctuating field, satisfying
conditions $\left<\mathbf{h}_{i}(t)\right>=0$ and
$\left<h_{i\alpha}(t)h_{j\beta}(t^{\prime})\right>=\mu\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$
[14, 15, 16]. Subscripts $\alpha$ and $\beta$ denote the Cartesian components
of a vector.

According to the fluctuation-dissipation theorem [20, 21, 22], one can prove
that in a thermal equilibrium, by identifying the energy distribution with the
Gibbs distribution (Eq. 1), parameters $\mu$ and $\gamma$ are related through
temperature $T$ as [23, 24, 18, 25, 15]:

$\mu=2\gamma\hbar k_{B}T.$ (3)

If $T$ is negative and $\gamma$ is positive, a seemingly straightforward idea
is to assign a negative value to $\mu$. However, $\mathbf{h}_{i}$ then becomes
a complex number vector. If one plugs such a complex fluctuating field into
Eq. (2), this makes $\mathbf{S}_{i}$ complex-valued, a clearly unnecessary
complication in the context of a classical spin dynamics considered here.

An alternative, simple yet computationally efficient and robust, solution is
proposed below. We assign a negative value to the parameter $\gamma$ instead,
with $\mu$ remaining positive. If we now examine the corresponding Fokker-
Planck equation [25] corresponding to Langevin equation (2), we see that this
imposes no constraint on the sign of $\gamma$ and still enables attaining a
Gibbs distribution with an arbitrary chosen $T$.

In terms of microscopic fluctuating spin dynamics, we note that if $\gamma$ is
negative, the dissipation term in Eq. (2) pushes the spin vector
$\mathbf{S}_{i}$ away from it being aligned with the effective field
$\mathbf{H}_{i}$ acting on it. This necessarily injects extra energy into the
spin system, in full agreement with the notion that a negative temperature
corresponds to a high-energy excited thermodynamic state of a spin ensemble.
We note that if the method were to be applied to a system where the spectrum
of energy states is not bounded from above, the Langevin term would continue
pumping energy into the system ad infinitum.

A remarkable feature of the Langevin algorithm involving a negative friction
parameter $\gamma$ and negative temperature $T$ is that it is now the friction
term that delivers energy to the statistical system. This implies that the
fluctuation term acts as effective damping, reversing the conventional meaning
of the two (fluctuation and dissipation) terms involved in a Langevin dynamics
simulation [26, 22, 27].

In what follows, we consider two case studies illustrating applications of the
algorithm. Both cases involve an ensemble of interacting spins described by a
general non-collinear Heisenberg Hamiltonian

${\cal H}=-\frac{1}{2}\sum_{i,j}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ (4)

where $\mathbf{S}_{i}=\mathbf{M}_{i}/g\mu_{B}$ is an atomic spin vector,
$\mathbf{M}_{i}$ is the corresponding atomic magnetic moment, $g=2.0023$ is
the electronic g-factor, $\mu_{B}$ is the Bohr magneton and $J_{ij}$ is the
exchange coupling parameter [28] between spins $i$ and $j$. Although there are
no literature data describing negative temperature states involving atomic
spin configurations, it does not appear impossible that an application of
advanced experimental techniques of spin flipping using high energy laser
pulses [29] could generate highy excited non-collinear negative temperature
spin states in an atomic ensemble. We note that the Ruderman-Kittel
interaction for nuclear spins [7] adopts the same functional form for the
Hamiltonian. Therefore, the following discussion is equally applicable to
nuclear as well as to atomic spins. Also, since a similar Hamiltonian
formalism applies also to alloys [30, 31], we note the general applicability
of the negative temperature notion to alloy configurations.

In the first case study, we assume that $J_{ij}\neq 0$ only for the first
nearest neighbours, choosing the Heisenberg parameter and magnetic moment
values of $J_{1}=-69.8$ meV and $|\mathbf{M}_{i}|=0.89$ $\mu_{B}$. In the
second case study, we assume $J_{1}=22.52$ meV, $J_{2}=17.99$ meV and
$|\mathbf{M}_{i}|=2.2$ $\mu_{B}$. Both types of interaction are realized on a
BCC lattice. The two sets of parameters selected above refer to pure chromium
and pure iron that have collinear antiferromagnetic and collinear
ferromagnetic ground states, respectively [32]. We denote them as model 1 and
model 2.

The simulation cells explored in the study involved 16,000 dynamic spins.
Initial configurations were assumed to be the respective ground states,
corresponding to $T=0^{+}$ K. Spin dynamics simulations were performed using
the coupled dynamic equations of motion, Eq. 2, see Refs. [14, 15, 16]. The
temperature of the system during a simulation was monitored using an explicit
expression derived in Ref. [25] for a dynamic spin ensemble. As was noted in
the derivation [25], the resulting value of the dynamic spin temperature can
be positive or negative.

$T=\frac{\sum_{i}|\mathbf{S}_{i}\times\mathbf{H}_{i}|^{2}}{2k_{B}\sum_{i}\mathbf{S}_{i}\cdot\mathbf{H}_{i}}.$
(5)

Since the numerator in the fraction is positive definite, the value of $T$
computed using this formula is negative if the denominator is negative, i.e.
if the direction of a spin $\mathbf{S}_{i}$ is on average anti-parallel to the
direction of the effective field $\mathbf{H}_{i}$ acting on it.

Figure 1: Process of thermalization of dynamic spin ensembles with 16,000
spins from their respective ground states to $-100$ K using a Langevin
thermostat with $\gamma=-0.05$. For model 1, $J_{1}=-69.8$ meV and
$|\mathbf{M}_{i}|=0.89$ $\mu_{B}$. For model 2, $J_{1}=22.52$ meV,
$J_{2}=17.99$ meV and $|\mathbf{M}_{i}|=2.2$ $\mu_{B}$. These parameters
approximately describe chromium and iron with collinear antiferromagnetic and
ferromagnetic ground states, respectively [32, 31]. (top) Evolution
illustrated using the conventional Kelvin temperature scale. (bottom)
Evolution illustrated using a plot involving an alternative temperature scale
[1, 2, 3, 4, 5], where the vertical axis shows values of $-1/T$ in K-1 units.

In Fig. 1a, we illustrate the dynamics of thermalization process for models 1
and 2 as a function of time from $0^{+}K$ to $-100K$, assuming $\gamma=-0.05$.
Temperatures are computed on the fly using Eq. (5), and are shown to approach
the prescribed values at the end of a simulation. Energy is pumped into the
spin system by a heat bath at $-100K$. Simulations show that temperatures
first climb to $\infty^{+}$ and then come back from $\infty^{-}$.

Indeed, a spin system at a negative temperature is hotter than the same system
at any positive temperature. One can adopt an alternative temperature scale
[1, 2, 3, 4, 5] using a parameter $-1/T$, which offers a better description of
how hot a statistical system is. In Fig. 1b, we plotted temperatures during
the thermalization process again, but now on a new temperature scale. Now the
curves are continuous and the functional dependence is monotonic. The average
energies of the spin systems are the same at $T=\infty^{+}$ and $\infty^{-}$.

Figure 2: (top) The equilibrium energy per spin as a function of temperature.
(bottom) The specific heat $C_{S}$ is obtained by performing numerical
differentiation. Peaks can be clearly observed in the negative temperature
region. They correspond to negative Curie and Néel temperature for model 1 and
2 respectively. The asymmetric peaks for model 2 are due to the competition of
$J_{1}$ and $J_{2}$ in an interacting spin system.

Figure 3: (top) Schematic picture for antiferromagnetic type-I. (bottom)
Schematic picture for antiferromagnetic type-II. Figure 4: Atomic spins of
model 2 (Fe) at $T=-1$ K.

By thermalizing the two model systems at various temperatures to equilibrium,
we find the average energy per spin as a function of temperature, plotted in
Fig. 2a. In Fig. 2b we also plot the specific heat $C_{S}$ calculated by
performing the numerical differentiation of data in Fig. 2a. Peaks of specific
heat are observed in Fig. 2b in both positive and negative temperature
regions. These peaks correspond to magnetic phase transitions and exhibit the
characteristic shape commonly referred to as the Curie-Weiss law [33].

In the positive temperature region, the occurrence of magnetic phase
transitions is well known. The peaks of the specific heat correspond to the
Néel temperature for model 1 and to the Curie temperature for model 2,
respectively. At a point where the specific heat diverges, the spin systems
undergoes a phase transition from a magnetically ordered antiferromagnetic or
magnetically ordered ferromagnetic state to a paramagnetic state with no long-
range spin orientational order.

In the negative temperature region, model 1 also exhibits a peak of the
specific heat, illustrated in Fig. 2b, that is symmetric with respect to the
one in the positive temperature region. The occurrence of this peak can be
interpreted as follows. Since the exchange coupling $J_{ij}$ is non-zero only
for the 1st nearest neighbours, if we transform $J_{1}\rightarrow-J_{1}$, the
ground state of model 1 changes from an antiferromagnetic type-I (Fig. 3) to
the ferromagnetic configuration. However, the effective field for any
particular spin does not change, and as a result the absolute value of the
transition temperature stays the same. Alternately, if we make the
transformation $T\rightarrow-T$, while not changing $J_{1}$, this produces the
same energy distribution (Eq. 1) as the reversal of the sign of $J_{1}$. Model
1 then exhibits a phase transition from a paramagnetic to a ferromagnetic
state at a negative Curie temperature.

For model 2, the peaks are asymmetric with respect to the origin of the
temperature axis. Since the exchange coupling constants $J_{ij}$ in model 2
extend to the 2nd nearest neighbour distance, by following the same logic as
above for model 1, we discover that $J_{1}$ and $J_{2}$ are now in competition
with each other at a negative temperature. There is no longer an obvious
preferred relative arrangement of spin orientations that any pair of spins
adopts at a negative temperature.

Dynamic simulations suggest that model 2 becomes an antiferromagnetic type-II
(Fig. 3) system when the temperature is higher than the negative Néel
temperature. Fig. 4 illustrates the configuration of atomic spins at $T=-1$ K.
The direction of alignment of spins shown in the Figure is arbitrary, and does
not correspond to the orientation of any Cartesian axis since the Hamiltonian
is invariant with respect to the rotation of a spin configuration as a whole.

The occurrence of type-II antiferromagnetic ordering can be rationalized from
the energetic perspective. If model 2 adopts an antiferromagnetic type-I
ordering, the energy of a magnetic moment at $T=0^{-}$ is
$8J_{1}-6J_{2}=72.22$ meV. If it adopts type-II ordering, the energy per
moment is $6J_{2}=107.94$ meV. Therefore, type-II ordering corresponds to
higher energy, which is preferable at a negative temperature. The difference
in the strength of the effective field in the positive and negative
temperature regions is the reason for the lack of mirror symmetry in the peak
positions in the plot for $C_{S}$. This result follows naturally from spin
dynamics simulations but cannot be readily obtained if we consider a
conventional mean-field approximation for the ferromagnetic and type-I
antiferromagnetic cases.

This highlights the fact that phase transitions in the negative temperature
region do not necessarily mirror their counterparts at positive temperatures.
The nature of phase transitions and the transition temperatures depend on the
range and strength of spin-spin interactions.

In conclusion, in this study we develop and describe a simple and robust
algorithm for the thermalization of a dynamic spin system to a negative
temperature using a Langevin spin dynamic thermostat. The algorithm enables
thermalizing a spin ensemble to any exact desired value of temperature,
irrespectively of whether it is positive or negative. There is no need to
introduce an external field to flip the directions of selected spins, and the
energy distribution is not perturbed by the Zeeman term. The case studies
explored above suggest that the systems with antiferromagnetic and
ferromagnetic ground states exhibit phase transitions both in the positive and
negative temperature regions. Negative Curie and Néel temperatures were found
in dynamic spin simulations. Finally, we note that although the examples
considered in this study only refer to dynamic spin ensembles, the same
algorithm should enable studying any dynamic statistical ensemble obeying the
fluctuation-dissipation theorem.

###### Acknowledgements.

This work has been carried out within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No 101052200 — EUROfusion) and was partially
supported by the Broader Approach Phase II agreement under the PA of
IFERC2-T2PA02. This work was also funded by the EPSRC Energy Programme (grant
number EP/W006839/1). To obtain further information on the data and models
underlying the paper please contact<EMAIL_ADDRESS>Views and
opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Commission.
Neither the European Union nor the European Commission can be held responsible
for them.

## References

  * Goldman [2008] M. Goldman, Overview of spin temperature, thermal mixing and dynamic nuclear polarization, Applied Magnetic Resonance 34, 219 (2008).
  * Abragam and Goldman [1978] A. Abragam and M. Goldman, Principles of dynamic nuclear polarisation, Reports on Progress in Physics 41, 395 (1978).
  * Ramsey [1956] N. F. Ramsey, Thermodynamics and statistical mechanics at negative absolute temperatures, Phys. Rev. 103, 20 (1956).
  * Klein [1956] M. J. Klein, Negative absolute temperatures, Phys. Rev. 104, 589 (1956).
  * Landsberg [1959] P. T. Landsberg, Negative temperatures, Phys. Rev. 115, 518 (1959).
  * Purcell and Pound [1951] E. M. Purcell and R. V. Pound, A nuclear spin system at negative temperature, Phys. Rev. 81, 279 (1951).
  * Oja and Lounasmaa [1997] A. S. Oja and O. V. Lounasmaa, Nuclear magnetic ordering in simple metals at positive and negative nanokelvin temperatures, Rev. Mod. Phys. 69, 1 (1997).
  * Hakonen and Yin [1990] P. Hakonen and S. Yin, Studies of nuclear magnetism in silver at positive and negative nanokelvin temperatures, Physica B: Condensed Matter 165-166, 785 (1990), proceedings of the 19th International Conference on Low Temperature Physics.
  * Hakonen  _et al._ [1992] P. J. Hakonen, K. K. Nummila, and R. T. Vuorinen, Spin dynamics in highly polarized silver at negative absolute temperatures, Phys. Rev. B 45, 2196 (1992).
  * Hakonen  _et al._ [1994] P. Hakonen, R. Vuorinen, and J. Martikainen, Anomalous spin-lattice relaxation in dilute RhFe at positive and negative nanokelvin temperatures, Physica B: Condensed Matter 194-196, 341 (1994).
  * Krommes and Rath [2003] J. A. Krommes and S. Rath, Monte Carlo sampling of negative-temperature plasma states, Phys. Rev. E 67, 066402 (2003).
  * Sørensen  _et al._ [2010] A. S. Sørensen, E. Altman, M. Gullans, J. V. Porto, M. D. Lukin, and E. Demler, Adiabatic preparation of many-body states in optical lattices, Phys. Rev. A 81, 061603 (2010).
  * Rapp  _et al._ [2010] A. Rapp, S. Mandt, and A. Rosch, Equilibration rates and negative absolute temperatures for ultracold atoms in optical lattices, Phys. Rev. Lett. 105, 220405 (2010).
  * Ma and Dudarev [2011] P.-W. Ma and S. L. Dudarev, Langevin spin dynamics, Phys. Rev. B 83, 134418 (2011).
  * Ma and Dudarev [2012] P.-W. Ma and S. L. Dudarev, Longitudinal magnetic fluctuations in Langevin spin dynamics, Phys. Rev. B 86, 054416 (2012).
  * Ma and Dudarev [2014] P.-W. Ma and S. L. Dudarev, Dynamic magnetocaloric effect in bcc iron and hcp gadolinium, Phys. Rev. B 90, 024425 (2014).
  * Antropov  _et al._ [1997] V. P. Antropov, S. V. Tretyakov, and B. N. Harmon, Spin dynamics in magnets: Quantum effects and numerical simulations (invited), Journal of Applied Physics 81, 3961 (1997).
  * Ma  _et al._ [2008] P.-W. Ma, C. H. Woo, and S. L. Dudarev, Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron, Phys. Rev. B 78, 024434 (2008).
  * Landau and Binder [2009] D. P. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_, 3rd ed. (Cambridge University Press, 2009).
  * Chandrasekhar [1943] S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod. Phys. 15, 1 (1943).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on Progress in Physics 29, 255 (1966).
  * Coffey  _et al._ [2004] W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron,  _The Langevin Equation_, 2nd ed. (World Scientific, 2004).
  * Brown, Jr. [1963] W. F. Brown, Jr., Thermal fluctuations of a single-domain particle, Phys. Rev. 130, 1677 (1963).
  * García-Palacios and Lázaro [1998] J. L. García-Palacios and F. J. Lázaro, Langevin-dynamics study of the dynamical properties of small magnetic particles, Phys. Rev. B 58, 14937 (1998).
  * Ma  _et al._ [2010] P.-W. Ma, S. L. Dudarev, A. A. Semenov, and C. H. Woo, Temperature for a dynamic spin ensemble, Phys. Rev. E 82, 031111 (2010).
  * Mazo [2002] R. M. Mazo,  _Brownian Motion_ (Clarendon Press, 2002).
  * Dudarev [2008] S. L. Dudarev, The non-Arrhenius migration of interstitial defects in bcc transition metals, Comptes Rendus Physique 9, 409 (2008).
  * Turek  _et al._ [2006] I. Turek, J. Kudrnovský, V. Drchal, and P. Bruno, Exchange interactions, spin waves, and transition temperatures in itinerant magnets, Philosophical Magazine 86, 1713 (2006).
  * Koopmans  _et al._ [2010] B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Explaining the paradoxical diversity of ultrafast laser-induced demagnetization, Nature Materials 9, 259 (2010).
  * Nguyen-Manh  _et al._ [2008] D. Nguyen-Manh, M. Y. Lavrentiev, and S. L. Dudarev, The Fe–Cr system: atomistic modelling of thermodynamics and kinetics of phase transformations, Comptes Rendus Physique 9, 379 (2008).
  * Lavrentiev  _et al._ [2016] M. Y. Lavrentiev, J. S. Wróbel, D. Nguyen-Manh, S. L. Dudarev, and M. G. Ganchenkova, Magnetic cluster expansion model for random and ordered magnetic face-centered cubic Fe-Ni-Cr alloy, Journal of Applied Physics 120, 043902 (2016).
  * Lavrentiev  _et al._ [2011] M. Y. Lavrentiev, R. Soulairol, C.-C. Fu, D. Nguyen-Manh, and S. L. Dudarev, Noncollinear magnetism at interfaces in iron-chromium alloys: The ground states and finite-temperature configurations, Physical Review B 84, 144203 (2011).
  * Desai [1986] P. D. Desai, Thermodynamic properties of iron and silicon, Journal of Physical and Chemical Reference Data 15, 967 (1986).


# Key Point-based Orientation Estimation of Strawberries for Robotic Fruit
Picking

Justin Le Louëdec1 and Grzegorz Cielniak1 1Lincoln Centre for Autonomous
Systems, University of Lincoln, Brayford Way, Brayford Pool, Lincoln LN6 7TS,
United Kingdom<EMAIL_ADDRESS>
###### Abstract

Selective robotic harvesting is a promising technological solution to address
labour shortages which are affecting modern agriculture in many parts of the
world. For an accurate and efficient picking process, a robotic harvester
requires the precise location and orientation of the fruit to effectively plan
the trajectory of the end effector. The current methods for estimating fruit
orientation employ either complete 3D information which typically requires
registration from multiple views or rely on fully-supervised learning
techniques, which require difficult-to-obtain manual annotation of the
reference orientation.

In this paper, we introduce a novel key-point-based fruit orientation
estimation method allowing for the prediction of 3D orientation from 2D images
directly. The proposed technique can work without full 3D orientation
annotations but can also exploit such information for improved accuracy. We
evaluate our work on two separate datasets of strawberry images obtained from
real-world data collection scenarios. Our proposed method achieves state-of-
the-art performance with an average error as low as $8^{\circ}$, improving
predictions by $\sim 30\%$ compared to previous work presented in [1].
Furthermore, our method is suited for real-time robotic applications with fast
inference times of $\sim 30$ms.

Automation through robotisation of the agricultural sector is seen as a
promising solution to the socio-economic challenges faced by this industry. A
key application that would benefit most from automation, which still relies on
manual human labour, is selective crop harvesting. Cultivated strawberries
(Fragaria x ananassa) are a perfect example of a crop with a recent
significant increase in demand but affected by labour shortages. Like many
varieties of fruit, they have several characteristics, rendering their
harvesting challenging for robots. To perform precise manipulation and
grasping of the harvestable crop, robotic systems require precise information
about the location and pose of the crop (see Fig. 1). Whilst detecting and
localising strawberry crop has been well-studied in prior work (e.g., [2, 3]),
inferring their precise pose is still an ongoing challenge. Typical approaches
for estimating the location and pose of the fruit rely on a combination of 2D
images together with 3D information (e.g., [4]) and require additional
projection steps and direct operations on point clouds which add significant
computation overhead.

Figure 1: An example strawberry harvesting robot with an example strawberry
indicating the Roll, Pitch, Yaw (RPY) orientation representation assumed in
this work.

In this paper, we present a method that enables precise fruit pose estimation
directly from 2D images allowing for bypassing computationally expensive
transformation and estimation steps. A similar approach has been originally
proposed by [1] which estimates the fruit pose by regressing the orientation
direction vector directly from an image in a fully supervised manner using a
learnt feature extractor and regression layer. The method is characterised by
fast inference times but also reduced accuracy caused by occasional failures
which are difficult to attribute to any particular part of the architecture.
In addition, the supervised approach relies on difficult-to-obtain high-
quality orientation annotations making the method impractical for wider use.
In contrast, the proposed method employs a two-stage approach first predicting
the localisation of two characteristic key points of the strawberry fruit
which are then used for pose estimation. This two-stage approach, inspired by
prior work in human pose estimation (i.e., [5]), leads to improved accuracy
and a simplified training regime thanks to the straightforward annotation of
key points. In particular, the contributions of this paper include:

  * •

a new method for efficient estimation of fruit orientation directly from
colour images based on key-point detectors;

  * •

an improved roll angle estimation method based on learnt image features and
the estimated key point information;

  * •

evaluation of the proposed fruit pose estimation system on two different
datasets of strawberry images consisting of high-resolution reconstructed
models of fruit and data collected directly in the field demonstrating the
improvements to the state of the art.

##  I Related work

Obtaining accurate pose estimation of objects, and especially the orientation
component is essential for enabling precise manipulation and grasping tasks
such as robotic picking [6, 7, 8]. Whilst early approaches utilised
handcrafted image descriptors [9, 10] for that task, most recent research
focuses on learning approaches exploiting the existing extensive datasets. The
examples include work making use of multiple-views for precise orientation and
category estimation such as in [11], through iterative refinement of object’s
pose using Gated Recurrent Unit (GRU) operators as in [12] or by employing
transformer architectures as in [13]. The majority of the state of that art
methods, however, rely on precise and difficult-to-obtain annotations and
therefore there is a big interest in methods which can reduce that requirement
for example through the use of simulation for training as in Sim2Real [14].
When accurate 3D point clouds are available, the use of 3D geometric insights
can be used to improve orientation prediction across the whole categories of
objects [15]. All these techniques, however, do not focus on a specific
application and rely on benchmark datasets of non-organic objects exclusively
(e.g., YCB-V dataset [16]), making their usefulness for agri-robotics
applications more difficult to assess.

Object pose estimation for robotic manipulation and grasping has been
researched thoroughly over the years. For example, the approach in [17] is
proposing the use of simulation and self-supervised training through
manipulating objects to identify their pose. There are also methods relying on
point cloud data and improved descriptors which were applied to obtain better
accuracy of object orientation for robotic bin picking [18]. Crop pose
estimation for robotic harvesting has also gained recent interest with various
applications in different crops. For example, the approach presented in [19]
estimates the pose of guava fruit from point cloud data obtained from an RGBD
sensor by segmenting plant components (i.e. fruit and branches) and combining
their relative pose. Other methods, such as [20], propose the refinement of
the fruit pose by registering a 3D reconstruction of the captured point cloud
to offline templates of the identified fruit [20]. The method presented in
[21] predicts the crop detection bounding boxes, maturity, pose and precise
stem orientation to identify the optimal cutting point for tomatoes, obtaining
more accurate and thorough information for harvesting.

For strawberries, recent work introduced a learning-based regression of the
orientation vectors of the fruits from a single colour and, optionally, depth
image [1]. While achieving state-of-the-art results, the method’s accuracy is
affected by occasional failures which due to the dimensional difference
between the images and the output 3D orientation vector are difficult to
analyse. In addition, the fully-supervised nature of the method requires
accurate annotation of the fruit orientation in images which are complex to
obtain without additional geometrical information about the fruit size and
shape.

In contrast to the state of the art, our method employs key-point detectors
for estimating 3D fruit orientation directly from 2D images. The proposed
technique can work without full 3D orientation annotations but can also
exploit such information for improved accuracy.

##  II The Approach

The proposed method consists of a set of key components which include a learnt
key point detector, orientation calculation and an optional, learnt component
for improved estimation of the orientation. In this work, we assume the “roll-
pitch-yaw” representation as per [1]. In the proposed application (see Fig.
1), the yaw angle is irrelevant due to the symmetry and non-uniformity of the
strawberries. Thus, we simplify the definition of the orientation to two
angles: pitch ($\phi\in\left[-180,180\right]$) and roll
($\theta\in\left[0,90\right]$). Both angles can be derived numerically from
the detected crop key points in 2D, but in addition, we demonstrate that the
calculation of $\theta$ can be regressed from the image and the key points,
leading to improved results when compared to the direct numerical formula. The
main advantage of key points is the ease of their annotation in images which
involves the marking of a single-pixel location only.

In our case, the two “top” and “tip” key points represent projections of the
stem attachment point and the extreme point of the fruit, respectively onto
the image plane. Due to the fruit growing conditions and typical harvesting
robot camera configuration, the tip is always located in space between being
parallel to the image plane or pointing towards the camera. At the same time,
the top might become obstructed by the crop and invisible in the image. In
such a case, its projection on the image plane relative to the centre of mass
is a good indication of the fruit’s size. An example demonstrating the
appearance of key point location under different fruit rotations is presented
in Fig. 3.

###  II-A Key point detection

Figure 2: An hourglass network composed of a succession of residual modules
consisting of convolution layers with skip connection from input to output as
in [22]. The input is the colour RGB image, and the output is two heat maps
indicating the location of each key point.

We propose using a method inspired by human pose estimation literature for
predicting the position of the key points. For each strawberry image, we
estimate two heat maps corresponding to the “top” and “tip” key points. The
top and tip key points correspond to the relative position of the stem
attachment point and outer part of the berry respectively projected onto the
image plane (see Fig. 3). We use the method presented in [5] composed of a
stack of $S$ small hourglass networks (see Fig. 2). For our scenario, we
choose experimentally $S=8$ with one final convolution layer followed by a
sigmoid function for each hourglass network. This last activation function
forces the network outputs within the range $[0,1]$ and can be interpreted as
the likelihood of the location of the key points. The output of an
intermediary hourglass module is combined with its feature map and input for
the next module in the stack.

Figure 3: The location of key points and relative distances used for
orientation calculation under changing values of $\theta$.

During training, we use the average binary cross-entropy as the loss function
$L_{K}=-ylog(f(x))-(1-y)log(1-f(x))$, where$f(x)$ is the output of each
stacked hourglass and $y$ the ground truth heat map for the key point. During
inference, the key point location corresponds to the image coordinates at the
maximum value selected from the $S=8$ predicted heat maps for each key point.
The ground truth heatmaps consist of 2D Gaussian kernels ($\sigma=2.0$)
applied to images containing each annotated key point location.

###  II-B Key point-based orientation

We define the coordinates of the key points as $A$ and $B$ which correspond to
the top and tip of the strawberries respectively. The roll $\phi$ can be
calculated in a straightforward manner: $\phi=\arctan(\frac{y_{BA}}{x_{BA}})$,
where $\vec{BA}=(x_{BA},y_{BA})$ is a vector between the both key points.

We propose to derive pitch $\theta$ from the relative position of the key
points to the centre of mass and also to the geometric outline of the
strawberry. As the fruit rotates, the two key points get closer/further to
each other and it is therefore possible to correlate their relative distance
to the orientation of the berry. To this end, we define the following four
distances with respect to the centre of mass:

  * •

$d_{top}$ and $d_{tip}$ representing the distance to the top and tip key
points, respectively;

  * •

$d_{topside}$ and $d_{tipside}$ is the distance to the fruit’s contour
following the straight line intersecting the top and tip key points,
respectively.

We then use these measures to derive normalised distances to the key points as
ratios $\hat{d}_{top}=\frac{d_{top}}{d_{topside}}$ and
$\hat{d}_{tip}=\frac{d_{tip}}{d_{tipside}}$. The key point locations and
corresponding distances relative to the centre of mass for different values of
$\theta$ are illustrated in Fig. 3.

Due to the specific shape and the skewed centre of mass of strawberries, the
3D rotation of the fruit results in an elliptical trajectory of the key points
in the image space which we model using a simple square root model (see Eq.
1). With the changing pitch angle $\theta$, the $\hat{d}_{top}$ decreases fast
until a specific value ($\sim\theta=50^{\circ}$ in our case), whilst
$\hat{d}_{tip}$ is not affected as much (see Fig. 3). The distance
$\hat{d}_{tip}$ is, however, a good indicator for $\theta$ values above that
threshold. Finally, with the image resolution normalised, we can use the
length of the berry $T$ as a threshold between the phase of importance for
$\hat{d}_{top}$ and $\hat{d}_{tip}$. With these considerations, we define
$\theta$ as:

$\theta=\begin{cases}\sqrt{\hat{d}_{top}}\alpha&\quad\text{if}\quad
d_{tt}>T,\\\
\sqrt{\hat{d}_{tip}}\omega+\sigma&\quad\text{otherwise,}\end{cases}$ (1)

where $d_{tt}$ is the distance between the top and tip key points. The
parameters $\alpha$, $\omega$ and $\sigma$ define the non-linear relationship
between the relative key point distances and $\theta$. The values of these
parameters are correlated with a typical shape of the fruit and in our case,
for strawberries, these were tuned experimentally to $T=170.0$, $\alpha=54.0$,
$\omega=50.0$ and $\sigma=40.0$.

###  II-C Improved estimation of the pitch angle

In general, obtaining the orientation ground truth for images is complex and
not always possible. When this ground truth information is available (eg. with
simulated data where object orientation is known), however, we propose a
supervised method to predict the pitch angle $\theta$ from a strawberry image,
which results in improved estimation results. Similarly to the previous work
presented in [1], we use a pre-trained VGG16 [23] architecture to extract
relevant features from a strawberry image. However, the size of the regressor
used in our architecture is larger to improve the predictive capabilities of
the network. We then combine the extracted feature map with the detected key
point locations before using a two-stage classifier to regress $\theta$. The
architecture used to predict the $\theta$ is presented in Fig. 4.

Figure 4: A network architecture for predicting $\theta$ from the two key
points and an input image. The feature extractor (orange) consists of
convolutions (light shade) and max-pooling (dark shade) layers whilst the
$\theta$ regressor (purple) includes fully-connected layers and the sigmoid
activation function. The red cube represents the key point location input. The
network outputs a single value $\theta\in\left[0,\frac{\pi}{2}\right]$.

During training the loss is expressed as the mean squared error
$L_{\theta}=(f(x)-y)^{2}$, where $f(x)$ is the output of the network and $y$
is the ground truth value for $\theta$.

##  III Evaluation Setup

###  III-A Data collection

To train and evaluate our method for orientation estimation of the strawberry
fruit, we introduce two distinct datasets of strawberry images. $Straw_{2D}$
consists of strawberry images collected in-field and annotated with simple key
point annotations whilst $Straw_{3D}$ includes images generated from high-
quality 3D models of strawberries, which in addition to key point annotations,
include also the full orientation ground truth.

####  III-A1 $Straw_{2D}$

is based on the dataset provided in [24], which includes high-resolution
images of strawberries from 20 different plantations in Spain (see Fig. 5).
From these images, 500 individual strawberries with minimal occlusion and of
different shape, orientation and maturity stage were extracted and
preprocessed so that background and calyx were masked out. The cropped square
area around each berry was resized to the resolution of
$$256\text{\,}\mathrm{px}$\times$256\text{\,}\mathrm{px}$$. The top and tip
key points were then manually annotated by marking their projected location on
the cropped image plane.

|
---|---
Figure 5: Example original images from [24] (top row) and the annotated
instances of individual strawberries from $Straw_{2D}$ together with the
annotation of the top (pink) and tip (blue) key points (bottom row).

####  III-A2 $Straw_{3D}$

is a dataset of strawberry images obtained from high-resolution 3D models. In
addition to annotated key point locations, full orientation annotation is also
available since the models can be rendered at the arbitrary pose. For creating
these models, we have used multi-view stereo-photogrammetry which allows for
reconstructing 3D point clouds from a set of images with unknown poses. The
object is photographed from multiple views covering all sides of the object’s
geometry. For each image, features are computed using standard algorithms such
as SIFT [25] and are used to infer the relative pose between all the views.
Pair-wise stereo depth prediction is then applied to produce a 3D object from
every viewpoint. Photogrammetry creates precise 3D shapes and details without
the need for distance measurements (e.g. typically provided by expensive lidar
sensors) but requires significant computational resources which significantly
increase with the number of images and the required precision. To create 3D
shape ground truth of strawberries using photogrammetry, we use the Agisoft
Metashape software [26].

Our setup for capturing the 3D models consists of a high-resolution camera
(“Olympus E-50” with a focal length of 45 mm) mounted on a tripod and a
manually rotated table carrying a strawberry fixed with two picks preventing
the slippage with a strawberry (see Fig. 6(a)). The images are taken at 30 cm
distances from several viewpoints and at regular rotation intervals. We
capture high-resolution images of
$$3264\text{\,}\mathrm{px}$\times$2448\text{\,}\mathrm{px}$$ allowing for
high-quality textures for realistic rendering. On average, we capture $\sim
80$ images per strawberry. To match the appearance of the models to the field
conditions, we capture the dataset outdoors using sunlight rather than
artificial lights. We also use a white background to easily mask the fruit and
improve the quality of the reconstructed shape.

Our new dataset is significantly larger and offers a greater variability in
the data distribution than previous work, such as [27], with increased texture
and mesh quality. Of the 127 strawberries captured, 36 were from an unknown
species bought at the farmer’s market, and the rest were a mix between Zara
and Katrina. Furthermore, we captured half of the strawberries at the unripe
maturity stage and the other half fully ripe for greater variance in shape and
appearance. The realism of the resulting rendered models is demonstrated in
Fig. 6(a).

(a) The photogrammetry setup used for capturing the 3D models of strawberries.
Figure 6: Comparison between the images of real berries (top row) and their
realistic 3D renders (bottom row).

|  |  |
---|---|---|---

To create the $Straw_{3D}$ dataset, each 3D berry is rendered in 84 different
orientations resulting in 10668 individual views. The images were post-
processed in a similar way as the $Straw_{2D}$ dataset resulting in squared
and masked images with the resolution of
$$256\text{\,}\mathrm{px}$\times$256\text{\,}\mathrm{px}$$. The annotation
consists of the orientation ground truth (known due to the data being
simulated) and the projected 3D key points (top and tip) on their image
locations. The resulting image examples together with annotations are
presented in Fig. 7.

Figure 7: Selected examples from the $Straw_{3D}$ dataset, consisting of post-
processed images and annotations including top (pink) and tip (blue) key
points.

Overall, we have 500 images representing 500 strawberry instances for
$Straw_{2D}$ and 10668 image instances of 127 different strawberries for
$Straw_{3D}$. Both datasets are annotated with ground truth key points, but
only $Straw_{3D}$ has an associated orientation ground truth.

###  III-B Evaluation metrics and training details

To evaluate the performance of our model for predicting the key point
locations, we use the Euclidean distance
$e=\sqrt{(x-x_{gt})^{2}+(y-y_{gt})^{2}}$ between the predicted $(x,y)$ and
ground truth $(x_{gt},y_{gt})$ location expressed in pixels. For evaluating
the predicted 3D pose for the strawberries from $Straw_{3D}$, we compute the
angular distance between the predicted and ground truth pose. The direction
vector $V=RV_{neutral}^{\intercal}$ where $V_{neutral}=\left[0,-1,0\right]$
and $R$ the rotation matrix obtained with $\phi$ and $\theta$. The angular
error $\epsilon$ between the predicted $V_{pred}$ and ground truth $V_{gt}$
direction is then computed as $\epsilon=acos\left(V_{pred}V_{gt}\right)$. We
use a recent method presented in [1] as the baseline for comparisons on the
$Straw_{3D}$ dataset. We use a 10-fold cross-validation evaluation scheme for
all datasets, with a dataset split of 90% for training and 10% for testing.
For key points prediction, we train with a starting learning rate of $110-4$
for ${\sim}100$ epochs for $Straw_{2D}$ and ${\sim}6$ epochs for $Straw_{3D}$.
We train our baseline [1] and the $\theta$ regressor with a starting learning
rate of $110-5$ until convergence (${\sim}6$ epochs).

##  IV Results

We first analyse the quality of key point detection, perform comparisons of
the direct orientation computation and the supervised approach to the baseline
method and present a qualitative analysis of the predictions for the estimated
angles.

dataset |  $e_{tip}$ [px] |  $e_{top}$ [px] | $e_{phi}$
---|---|---|---
$Straw_{2D}$ | 9.25$\pm$4.20 | 14.91$\pm$ 9.45 |  $14.53^{\circ}\pm$ 39.85
$Straw_{3D}$ | 7.61$\pm$2.61 | 8.60$\pm$3.21 |  $12.18^{\circ}\pm$ 27.05
TABLE I: Median estimation errors for the top and tip key point locations in
$Straw_{2D}$ and $Straw_{3}D$ datasets. As well as $\phi$ error, when computed
from these key points.

Table I presents the median prediction errors of our key point detector for
both datasets. For the in-field dataset $Straw_{2D}$, the results are around
$9\text{\,}\mathrm{px}$ (${\sim}4\%$ image size) for the tip and
$19\text{\,}\mathrm{px}$ (${\sim}7\%$ image size) for the top but with a high
variance, which is still relatively good for imperfect fruit images collected
from the field. This higher variance can be attributed to the quality of
manual annotation, which was affected by the difficulties in identifying key
points and strawberry shape for this dataset. Indeed with atypical shapes, and
missing flesh $Straw_{2D}$, presents more challenges for annotation and
prediction. With more accurate annotations present in the $Straw_{3D}$
dataset, however, the algorithm identifies the location of the key points more
accurately which indicates that this is a critical consideration in training
the key point detectors for real applications. Furthermore, the error in
$\phi$ prediction ($e_{phi}$), shows a clear correlation between key-points
accuracy and angle prediction, with the high-variance coming from difficult
cases and ambiguously shaped strawberries.

We further show the correlation between key point localisation error and fruit
orientation in Fig. 8(b) on $Straw_{3D}$. We can see that $e_{top}$ is
consistent across different orientations but with a higher variance. $e_{tip}$
spikes mainly when the tip of the fruit points toward the camera and is harder
to distinguish precisely (${\sim}70^{\circ}$). Indeed, the top key point
aligns with the centre of mass when $\theta$ gets close to $90^{\circ}$, while
the tip key point is often displaced randomly due to the strawberry shape and
growth. It is worth noting that more precise tip point estimation does not
improve orientation results at lower values of $\theta$ as shown in Fig. 3.

Furthermore, we show in Fig. 8(a) that the key points are predicted accurately
in most of the cases. In $Straw_{2D}$, the errors are primarily due to
inaccurate top prediction due to arguable and difficult to annotate precisely.
We see in the last example also key points predicted on the centre of mass
indicating a value of $\theta=90^{\circ}$, probably due to the round shape and
missing flesh information. For $Straw_{3D}$, an higher imprecision with the
tip prediction for high values of $\theta$ can be observed.

(a)

(b)

Figure 8: (a) Examples of accurate (top row) and less accurate (bottom row)
key point predictions from the $Straw_{2D}$ dataset. (b) Key point prediction
error relative to the $\theta$ angle for $Straw_{3D}$.

(a) Computed $\phi$ values for selected strawberry instances.

(b) A comparison of estimation accuracy for $\theta$ and key point
localisation between real strawberries (black background) and their 3D renders
(white background). The tip key point is indicated in blue whilst the
intersection with the contour is indicated in orange.

Figure 9: Qualitative results obtained using the numerical method applied to
the $Straw_{2D}$ dataset.

The example output from the key point detectors applied to the $Straw_{2D}$
dataset together with numerically calculated angles $\phi$ and $\theta$ is
presented in Fig. 9(a) and Fig. 9(b) respectively. We indicate the tip key
point on the rendered strawberries for $\theta$ to compare its localisation
with respect to the reference strawberry. With accurate key points, in Fig.
9(a) the numerical prediction of $\phi$ gives a clear indication of the fruit
orientation within the image plane. The numerical computation of $\theta$ in
Fig. 9(b) shows that using the tip key point and contour of the fruit leads to
accurate orientation estimates, with the template strawberry meshes precisely
aligning with the target images.

(a) Distribution of angular errors and median and mean values.

(b) Median estimation errors with respect to different values of $\theta$.

Figure 10: Quantitative evaluation of the proposed key point-based orientation
estimation, supervised method and [1] on $Straw_{3D}$.

The comparison of the predicted poses from our proposed methods including
direct numerical estimation and the learnt $\theta$ estimator to the baseline
from [1] is presented in Fig. 10. The evaluation was performed on the
$Straw_{3D}$ dataset which includes the full pose annotation. The supervised
orientation method performs significantly better with lower error ($\sim
8^{\circ}$), and less variance than the baseline ($\sim 11^{\circ}$).
Furthermore, as displayed in Fig. 10(b), the numerical method shows viable
results in general with worse performance for the most ambiguous poses
$\theta=0^{\circ}$ and $\theta=49^{\circ}$ only. The latter is due to the
ambiguous orientation lying in between the two parts stages of the formula
used. This is a limitation of the numerical formula based on 2D images which
is difficult to overcome without employing additional 3D shape information.
The orientation prediction performs more accurately for all methods as
$\theta$ gets closer to $90^{\circ}$. This is coherent with the lack of
variation in appearance for values of $\theta$ bellow $30^{\circ}$ as
presented in Fig. 3.

Figure 11: Image projected orientation vector $V_{pred}$ (orange, green,
white) for our method and baseline compared to the ground truth direction
$V_{gt}$ (blue) on $Straw_{3D}$. The projected vector direction corresponds to
$\phi$ and its length to $\theta$. The predicted (triangle) and ground truth
(star) tip (blue) and top (pink) key points are also indicated.

Finally, we visually compare the predicted orientation vectors for all three
methods in Fig. 11. The white direction vectors ([1]) confirm the results from
Fig. 10(a), with the baseline always slightly off the ground truth and
difficult to interpret in some cases. On the other hand, the key points give a
better understanding when our method shows imprecision.

Our method is computationally very efficient and suitable for real-time
robotic applications with inference times of $30.0$ ms for the key point
detector and $1.4$ ms for the supervised regression of $\theta$. The
performance is measured on an NVIDIA GeForce GTX 1880 Ti, with Intel(R)
Core(TM) i7-7700K CPU and 16 GB memory.

##  V Conclusions and future work

In this paper, we presented a novel approach for predicting the orientation of
strawberries from single-view images. Our method exploits the key points and
understanding of the fruit’s shape, with two different techniques for
predicting the rotation angles of the fruit. The experiments indicate that our
approach achieves state-of-the-art results with average errors as low as
$8^{\circ}$. In addition, our key point-based method leads to a better
understanding of failure cases when compared to the baseline, clearly
indicating sources of errors which are directly linked to mislocalised key
points. Our method is suited for robotic strawberry harvesting where the
fruit’s orientation is important for effective end-effector motion planning.

Future work will include further optimisation of the key point prediction
network and supervised regression of $\theta$. For example, sharing the weight
from the feature extraction layers should reduce the training time and
complexity of the models with the number of weights needed. The developed
numerical $\theta$ computation was developed without considering the possible
berry shapes which, if taken into account, should improve the results by
improving the correlation between the key point location and orientation.
Furthermore, this work can easily be extended to other crops by modifying the
numerical formula. Future work would also include considering the external and
self-occlusion of the key points, and adding uncertainty measures for hidden
key points.

## References

  * [1] N. Wagner, R. Kirk, M. Hanheide, and G. Cielniak, “Efficient and robust orientation estimation of strawberries for fruit picking applications,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 13 857–13 863.
  * [2] R. Kirk, G. Cielniak, and M. Mangan, “L* a* b* fruits: A rapid and robust outdoor fruit detection system combining bio-inspired features with one-stage deep learning networks,” _Sensors_ , vol. 20, no. 1, p. 275, 2020.
  * [3] Y. Ge, Y. Xiong, and P. J. From, “Instance segmentation and localization of strawberries in farm conditions for automatic fruit harvesting,” _IFAC-PapersOnLine_ , vol. 52, no. 30, pp. 294–299, 2019.
  * [4] Y. Xiong, Y. Ge, and P. J. From, “An improved obstacle separation method using deep learning for object detection and tracking in a hybrid visual control loop for fruit picking in clusters,” _Computers and Electronics in Agriculture_ , vol. 191, p. 106508, 2021.
  * [5] A. Newell, K. Yang, and J. Deng, “Stacked hourglass networks for human pose estimation,” in _European conference on computer vision_. Springer, 2016, pp. 483–499.
  * [6] R. Harrell, P. D. Adsit, R. Munilla, and D. Slaughter, “Robotic picking of citrus,” _Robotica_ , vol. 8, no. 4, pp. 269–278, 1990.
  * [7] J. Baeten, K. Donné, S. Boedrij, W. Beckers, and E. Claesen, “Autonomous fruit picking machine: A robotic apple harvester,” in _Field and service robotics_. Springer, 2008, pp. 531–539.
  * [8] W. Liu, W. Wang, Y. You, T. Xue, Z. Pan, J. Qi, and J. Hu, “Robotic picking in dense clutter via domain invariant learning from synthetic dense cluttered rendering,” _Robotics and Autonomous Systems_ , vol. 147, p. 103901, 2022\.
  * [9] E. Muñoz, Y. Konishi, V. Murino, and A. Del Bue, “Fast 6d pose estimation for texture-less objects from a single rgb image,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 5623–5630.
  * [10] K. Park, J. Prankl, and M. Vincze, “Mutual hypothesis verification for 6d pose estimation of natural objects,” in _Proceedings of the IEEE International Conference on Computer Vision Workshops_ , 2017, pp. 2192–2199.
  * [11] A. Kanezaki, Y. Matsushita, and Y. Nishida, “Rotationnet: Joint object categorization and pose estimation using multiviews from unsupervised viewpoints,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 5010–5019.
  * [12] L. Lipson, Z. Teed, A. Goyal, and J. Deng, “Coupled iterative refinement for 6d multi-object pose estimation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 6728–6737.
  * [13] T. Jantos, M. Hamdad, W. Granig, S. Weiss, and J. Steinbrener, “PoET: Pose Estimation Transformer for Single-View, Multi-Object 6D Pose Estimation,” in _6th Annual Conference on Robot Learning (CoRL 2022)_.
  * [14] C. Zhong, C. Yang, F. Sun, J. Qi, X. Mu, H. Liu, and W. Huang, “Sim2real object-centric keypoint detection and description,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 36, no. 5, 2022, pp. 5440–5449.
  * [15] Y. Di, R. Zhang, Z. Lou, F. Manhardt, X. Ji, N. Navab, and F. Tombari, “Gpv-pose: Category-level object pose estimation via geometry-guided point-wise voting,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 6781–6791.
  * [16] Y. Xiang, T. Schmidt, V. Narayanan, and D. Fox, “Posecnn: A convolutional neural network for 6d object pose estimation in cluttered scenes,” 2018.
  * [17] X. Deng, Y. Xiang, A. Mousavian, C. Eppner, T. Bretl, and D. Fox, “Self-supervised 6d object pose estimation for robot manipulation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 3665–3671.
  * [18] X. Cui, M. Yu, L. Wu, and S. Wu, “A 6d pose estimation for robotic bin-picking using point-pair features with curvature (cur-ppf),” _Sensors_ , vol. 22, no. 5, p. 1805, 2022.
  * [19] G. Lin, Y. Tang, X. Zou, J. Xiong, and J. Li, “Guava detection and pose estimation using a low-cost rgb-d sensor in the field,” _Sensors_ , vol. 19, no. 2, p. 428, 2019.
  * [20] N. Guo, B. Zhang, J. Zhou, K. Zhan, and S. Lai, “Pose estimation and adaptable grasp configuration with point cloud registration and geometry understanding for fruit grasp planning,” _Computers and Electronics in Agriculture_ , vol. 179, p. 105818, 2020.
  * [21] J. Kim, H. Pyo, I. Jang, J. Kang, B. Ju, and K. Ko, “Tomato harvesting robotic system based on deep-tomatos: Deep learning network using transformation loss for 6d pose estimation of maturity classified tomatoes with side-stem,” _Computers and Electronics in Agriculture_ , vol. 201, p. 107300, 2022.
  * [22] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learningfor image recognition,” _ComputerScience_ , 2015.
  * [23] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” in _International Conference on Learning Representations_ , 2015.
  * [24] I. Pérez-Borrero, D. Marín-Santos, M. E. Gegúndez-Arias, and E. Cortés-Ancos, “A fast and accurate deep learning method for strawberry instance segmentation,” _Computers and Electronics in Agriculture_ , vol. 178, p. 105736, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0168169920300624
  * [25] D. G. Lowe, “Object recognition from local scale-invariant features,” in _Proceedings of the seventh IEEE international conference on computer vision_ , vol. 2. Ieee, 1999, pp. 1150–1157.
  * [26] A. Software, “Agisoft photoscan professional,” _Retrieved from http://www.agisoft.com/downloads/installer/_ , pp. 1534–4320, 2016.
  * [27] J. Q. He, R. J. Harrison, and B. Li, “A novel 3d imaging system for strawberry phenotyping,” in _Plant Methods_ , 2017.


$\displaystyle=-\frac{\gamma}{\lambda}\frac{\hat{i}}{\hat{i}_{E}}+1,$ (6.51)
$\displaystyle\dfrac{d\hat{i}_{E}}{dt}$
$\displaystyle=-1-\frac{\gamma}{\lambda}+\mu\hat{s}-\frac{\hat{i}_{E}}{\hat{s}}+2\frac{\hat{w}}{\hat{s}}-\frac{\omega}{\lambda},$
(6.52) $\displaystyle\dfrac{d\hat{w}}{dt}$
$\displaystyle=\frac{\omega\alpha}{\lambda}\hat{s}-2\frac{\hat{w}}{\hat{s}},$
(6.53)

with initial condition
$(\hat{s}(0),\hat{i}(0),\hat{i}_{E}(0),\hat{w}(0))=(1,0,0,0)$. Note that the
only difference between (6.40)-(6.43) and (6.50)-(6.53) is that the final term
in (6.42) is replaced by $\frac{\omega}{\lambda}$ since now every rewiring
necessarily leads to a drop of an infective edge.

For $t\geq 0$, let
$\hat{\boldsymbol{X}}^{(n)}(t)=n^{-1}(\hat{S}^{(n)}(t),\hat{I}^{(n)}_{E}(t),\hat{W}^{(n)}(t))$
and $\hat{x}(t)=(\hat{s}(t),\hat{i}_{E}(t),\hat{w}(t))$. Let
$\hat{\zeta}=\inf\\{t>0:\hat{i}_{E}(t)=0\\}$. Observe that (6.50), (6.52) and
(6.53) form an autonomous system that is Lipschitz provided that $\hat{s}$ is
bounded away from $0$. Also, in an obvious notation, (6.13) becomes

$\hat{\alpha}^{(n)}_{R}(\xi)=n^{-2}\left(\frac{\gamma}{\lambda}\right)\frac{n^{E}_{2}}{n^{E}_{2}},$

leading to the bound (cf. (6.24))

$\tilde{\alpha}^{(n)}_{R}(\hat{X}^{(n)}(t))\leq
2\frac{\gamma}{\lambda}n^{-1}\log n.$

Consequently, it is straightforward to modify the proof of Theorem 2.2 to show
that, for any $t_{0}\in(0,\hat{\zeta})$,

$\sup_{0\leq t\leq
t_{0}}\left|\hat{\boldsymbol{X}}^{(n)}(t)-\hat{x}(t)\right|\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0\quad\mbox{as }n\to\infty.$ (6.54)

Now $\hat{i}_{E}^{\prime}(0)=\mu-1-\frac{\gamma+\omega}{\lambda}$, so
$\hat{i}_{E}^{\prime}(0)>0$ if and only if $\lambda>\lambda_{c}$ and
$\hat{i}_{E}^{\prime}(0)=0$ if $\lambda=\lambda_{c}$. Further, using (6.50),

$\dfrac{d^{2}\hat{i}_{E}}{dt^{2}}=-\mu-\frac{1}{\hat{s}}\dfrac{d\hat{i}_{E}}{dt}-\frac{1}{\hat{s}^{2}}\hat{i}_{E}+\frac{2}{\hat{s}}\dfrac{d\hat{w}}{dt}+\frac{2\hat{w}}{\hat{s}^{2}},$

so, using (6.53), when $\lambda=\lambda_{c}$,

$\hat{i}_{E}^{\prime\prime}(0)=-\mu+\frac{2\omega\alpha}{\lambda_{c}}=\frac{\mu(\omega(2\alpha-1)-\gamma)-2\omega\alpha}{\omega+\gamma}.$

Note that $\hat{i}_{E}^{\prime\prime}(0)>0$ if and only if the condition in
Theorem 2.4(a) is satisfied.

Let $\hat{\zeta}^{(n)}=\inf\\{t\geq 0:\hat{I}^{(n)}_{E}(t)=0\\}$ and
$\hat{T}^{(n)}=n-\hat{S}^{(n)}(\hat{\zeta}^{(n)})$, so $\hat{T}^{(n)}$ is the
size of the modified epidemic. Suppose that the condition in Theorem 2.4(a) is
satisfied. Then
$\lim_{\lambda\downarrow\lambda_{c}}(1-\hat{x}(\hat{\zeta}))>0$ and it follows
using (6.54) that there exists $\tau_{0}>0$ such that

$\lim_{\lambda\downarrow\lambda_{c}}\lim_{n\to\infty}{\rm
P}(n^{-1}\hat{T}^{(n)}>\tau_{0})=1.$

Theorem 2.4(a) now follows as $\bar{T}^{(n)}$, conditioned on the occurrence
of a major outbreak, is stochastically larger than $n^{-1}\hat{T}^{(n)}$.

Turning to part (b), suppose now that the model $\mathcal{E}^{(n)}$ is
modified so that if a susceptible rewires an edge from an infective to a
recovered individual then the edge to the infective is also retained. In the
construction this means that when an infective sends a warning down an
infective edge and that edge is rewired to a recovered individual then the
infective does not lose that infective edge. This leads to a model whose final
size is stochastically larger than that of $\mathcal{E}^{(n)}$. As before, we
run the model $\mathcal{E}^{(n)}$ until time $t_{n}$ defined above and then
the time-transformed modified model. The corresponding time-transformed
deterministic model is given by (6.50)-(6.53), with (6.52) replaced by

$\dfrac{d\hat{i}_{E}}{dt}=-1-\frac{\gamma}{\lambda}+\mu\hat{s}-\frac{\hat{i}_{E}}{\hat{s}}+2\frac{\hat{w}}{\hat{s}}-\frac{\omega}{\lambda}(1-\alpha+\alpha\hat{s}).$
(6.55)

Note that (6.50), (6.55) and (6.53) form an autonomous system that is
Lipschitz provided that $\hat{s}$ is bounded away from $0$. The proof
parallels that of part (a) in the obvious fashion, the only difference being
that now, when $\lambda=\lambda_{c}$,

$\hat{i}_{E}^{\prime\prime}(0)=-\mu+\frac{3\omega\alpha}{\lambda_{c}}=\frac{\mu(\omega(3\alpha-1)-\gamma)-3\omega\alpha}{\omega+\gamma}.$

Thus $\hat{i}_{E}^{\prime\prime}(0)<0$ if either $\omega(3\alpha-1)<\gamma$ or
$\mu<\frac{3\omega\alpha}{\omega(3\alpha-1)-\gamma}$. Moreover, a further
calculation shows that if $\lambda=\lambda_{c}$ and
$\hat{i}_{E}^{\prime\prime}(0)=0$ then
$\hat{i}_{E}^{\prime\prime\prime}(0)=-\frac{2\omega\alpha}{\lambda_{c}}<0$.
Part (b) of the theorem now follows using a similar argument to that used at
the end of the proof of part (a).

The bounding process used in the proof of Theorem 2.4(b) suggests an approach
to proving the conditions in part(a) are also necessary for the final size to
be discontinuous at $\lambda=\lambda_{c}$. Run the epidemic
$\mathcal{E}^{(n)}$ until time $t_{n}$ and then make the random time-scale
transformation. For the latter, let
$\tilde{\boldsymbol{X}}^{(n)}=\\{\tilde{\boldsymbol{X}}^{(n)}(t):t\geq 0\\}$,
where
$\tilde{\boldsymbol{X}}^{(n)}(t)=\tilde{\boldsymbol{x}}(\tilde{X}^{(n)}(t))$
$(t\geq 0)$ and $\tilde{\boldsymbol{x}}:E^{(n)}\to\mathbb{R}^{4}$ is defined
by (6.2). Let $\tilde{I}^{(n)}(t)$ and $\tilde{I}^{(n)}_{E}(t)$ denote
respectively the numbers of infectives and infective edges in
$\tilde{\boldsymbol{X}}^{(n)}$ at time $t$. Recall from Section 6.5.2 that for
the epidemic $\mathcal{E}^{(n)}$, conditional upon a major outbreak,
$I^{(n)}(t_{n})/I^{(n)}_{E}(t_{n})\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}L$ as $n\to\infty$. This suggests that there exists
$L^{\prime}>L$ and $t_{0}>0$ such that, conditional upon a major outbreak,
$\tilde{I}^{(n)}(t)\leq L^{\prime}\tilde{I}^{(n)}_{E}(t)$ for all
$t\in[0,t_{0}]$ with high probability.

The term $\frac{\omega}{\lambda}(1-\alpha+\alpha(1-\tilde{i}))$ in (6.42)
arises from rewirings to susceptibles or recovered. In particular, the
probability that a rewiring is to a recovered, and the infective edge is
dropped, is $\frac{n-X^{(n)}(t)-I^{(n)}(t)-2}{n-2}$. Modifying the model so
that this probability is replaced by
$\frac{n-X^{(n)}(t)-L^{\prime}I^{(n)}_{E}(t)-2}{n-2}$ leads to an upper
bounding process for the true model and hence also for the time-transformed
model. The latter has deterministic limit which satisfies (6.40)-(6.43) over
$[0,t_{0}]$ but with (6.42) replaced by

$\dfrac{d\tilde{i}_{E}}{dt}=-1-\frac{\gamma}{\lambda}+\mu\tilde{s}-\frac{\tilde{i}_{E}}{\tilde{s}}+2\frac{\tilde{w}}{\tilde{s}}-\frac{\omega}{\lambda}(1-\alpha+\alpha(1-L^{\prime}\tilde{i}_{E})).$
(6.56)

Thus (6.40), (6.56) and (6.43) form an autonomous system, that is Lipschitz
provided that $\hat{s}$ is bounded away from $0$ and hence susceptible to
analysis. Omitting the details, when $\lambda=\lambda_{c}$ we obtain
$\tilde{i}_{E}^{\prime\prime}(0)=-\mu+\frac{2\omega\alpha}{\lambda_{c}}$ and
$\tilde{i}_{E}^{\prime\prime\prime}(0)=-\frac{2\omega\alpha}{\lambda_{c}}<0$
if $\tilde{i}_{E}^{\prime\prime}(0)=0$, which imply that the conditions in
Theorem 2.4(b) are also necessary.

It may be possible to make the argument rigorous using branching processes
that bound $\mathcal{E}^{(n)}$ but that does not seem straightforward since
(i) $t_{0}$ depends on $\lambda$ and $\downarrow 0$ as
$\lambda\downarrow\lambda_{c}$ and (ii) the result is conditional upon a major
epidemic, the probability of which also $\downarrow 0$ as
$\lambda\downarrow\lambda_{c}$.

###  6.6 Proof of Theorems 2.7 and 2.8

The proof of Theorem 2.7 is omitted as it parallels that of Theorem 2.2 in
Section 6.2, with only very minor modification. Turning to Theorem 2.8, we
make the same random time-scale transformation as in Section 6.3, in which the
clock is slowed down by a factor given by the total force of infection.
Recalling (2.12)-(2.15), this leads to the deterministic model given by

$\displaystyle\dfrac{d\tilde{s}}{dt}$ $\displaystyle=-1,$ (6.57)
$\displaystyle\dfrac{d\tilde{i}}{dt}$
$\displaystyle=-\frac{\gamma}{\lambda}\frac{\tilde{i}}{\tilde{i}_{E}}+1$
(6.58) $\displaystyle\dfrac{d\tilde{i}_{E}}{dt}$
$\displaystyle=-1-\frac{\gamma}{\lambda}+\mu\tilde{s}-\frac{\tilde{i}_{E}}{\tilde{s}}+2\frac{\tilde{w}}{\tilde{s}}-\frac{\omega}{\lambda},$
(6.59) $\displaystyle\dfrac{d\tilde{w}}{dt}$
$\displaystyle=\frac{\omega\alpha}{\lambda}-2\frac{\tilde{w}}{\tilde{s}},$
(6.60)

with initial condition
$(\tilde{s}(0),\tilde{i}(0),\tilde{i}_{E}(0),\tilde{w}(0))=(1,0,0,0)$. The
equations for $(\tilde{s},\tilde{i}_{E},\tilde{w})$ form a closed system
having solution

$\displaystyle\tilde{s}(t)$ $\displaystyle=1-t,$ (6.61)
$\displaystyle\tilde{i}_{E}(t)$
$\displaystyle=\tilde{s}(t)\tilde{g}(\tilde{s}(t)),$ (6.62)
$\displaystyle\tilde{w}(t)$
$\displaystyle=\frac{\omega\alpha}{\lambda}\tilde{s}(t)(1-\tilde{s}(t)),$
(6.63)

where

$\tilde{g}(s)=\left(1+\frac{\gamma+\omega(1-2\alpha)}{\lambda}\right)\log\tilde{s}+\left(\mu-2\frac{\omega\alpha}{\lambda}\right)(1-\tilde{s}(t)).$

The solution is derived as follows. Equation (6.61) follows immediately from
(6.57) and the initial condition. Equation (6.63) is obtained by dividing
(6.60) by (6.57) and solving the resulting differential equation using the
integrating factor $\tilde{s}^{-2}$, together with the initial condition.
Finally, (6.62) is obtained by dividing (6.59) by (6.57), substituting for
$\tilde{w}$ from (6.63) and then solving the resulting differential equation
using the integrating factor $\tilde{s}^{-1}$, together with the initial
condition.

Let $\tilde{\zeta}=\inf\\{t>0:\tilde{i}_{E}(t)=0\\}$ and
$\hat{s}=\tilde{s}(\tilde{\zeta})$. Then either $\hat{s}=0$ or
$\tilde{g}(\hat{s})=0$. We thus investigate the solutions of $\tilde{g}(s)=0$
in $[0,1]$. Note that
$\tilde{i}_{E}^{\prime}(0)=\mu-1-\frac{\omega+\gamma}{\lambda}>0$, as
$\lambda>\lambda_{c}$, and $\tilde{s}(1)=0$, so $\tilde{\zeta}\in(0,1]$, so we
investigate the solutions of $\tilde{g}(s)=0$ in $[0,1)$. Now

$\tilde{g}^{\prime}(s)=\left(1+\frac{\gamma+\omega(1-2\alpha)}{\lambda}\right)s^{-1}-\mu+2\frac{\omega\alpha}{\lambda},$

so $\tilde{g}$ has at most one stationary point in $(0,1]$. Further,
$\tilde{g}(1)=0$ and
$\tilde{g}^{\prime}(1)=1-\mu+\frac{\omega+\gamma}{\lambda}<0$, as
$\lambda>\lambda_{c}$, so by considering $\tilde{g}(0+)$ it is seen easily
that $\tilde{g}$ has no root in $(0,1)$ if
$\lambda+\gamma+\omega(1-2\alpha)\leq 0$ and precisely one root in $(0,1)$ if
$\lambda+\gamma+\omega(1-2\alpha)>0$. Note that since
$\lambda_{c}=\frac{\gamma+\omega}{\mu-1}$, a necessary (but not sufficient)
condition for $\tilde{g}$ not to have a root in $(0,1)$ is that
$r(\mu,\gamma,\omega,\alpha)<0$, where $r(\mu,\gamma,\omega,\alpha)$ is
defined at (2.16).

Suppose that a major epidemic occurs. Let $t_{n}=\inf\\{t>0:S^{(n)}(t)\leq
n-\log n\\}<\infty$. Define the time-transformed process
$\tilde{\boldsymbol{X}}^{(n)}$ analogously to the SI model, using the
coordinate functions (6.26), with
$\tilde{\boldsymbol{X}}^{(n)}(0)=(n^{-1}(S^{(n)}(t_{n}),I^{(n)}_{E}(t_{n}),W^{(n)}(t_{n}))$.

Suppose that $r(\mu,\gamma,\omega,\alpha)<0$ and
$\lambda\in(\lambda_{c},\omega(2\alpha-1)-\gamma]$. Then, given any
$\epsilon\in(0,1)$, let $t_{1}=1-\epsilon/2$, so
$\tilde{s}(t_{1})=\epsilon/2$. Let
$\tilde{U}=[\frac{\epsilon}{2},1]\times[-\frac{\epsilon}{2},2\mu]\times[0,2\mu]$.
The vector field $\tilde{b}:\tilde{U}\to\mathbb{R}^{3}$ given by the right-
hand sides of (6.57), (6.59) and (6.60) is Lipschitz and the proof of Theorem
2.2 is easily modified to yield

$\lim_{n\to\infty}{\rm P}\left(\sup_{0\leq t\leq
t_{1}}\left|\tilde{\boldsymbol{X}}^{(n)}(t)-\tilde{x}(t)\right|>\frac{\epsilon}{4}\right)=0,$
(6.64)

where $\tilde{x}(t)=(\tilde{s}(t),\tilde{i}_{E}(t),\tilde{w}(t))$. By a
similar argument to that used at the end of the proof of Theorem 2.5(b) in
Section 6.3, there exists $t_{2}\in(0,t_{1})$ such that
$\tilde{I}^{(n)}_{E}(t)>0$ for all $t\in[0,t_{2}]$ with high probability. It
follows that

$\displaystyle\lim_{n\to\infty}{\rm P}(\bar{T}^{(n)}>1-\epsilon)$
$\displaystyle\geq\lim_{n\to\infty}{\rm P}(\tilde{S}^{(n)}(t_{1})<\epsilon)$
$\displaystyle=1,$

using (6.64) and the definition of $t_{1}$. Theorem 2.8(b) for the case when
$\lambda\in(\lambda_{c},\omega(2\alpha-1)-\gamma]$ now follows as
$\epsilon\in(0,1)$ is arbitrary.

The proof of Theorem 2.8 for the other cases, i.e. when the limit
$\tilde{\tau}<1$, follows a similar argument to that of Theorem 2.5(b) and
hence is omitted.

##  7 Concluding comments

We have presented a construction of an SIR epidemic with preventive rewiring
on an Erdős-Rényi random graph, together with a rigorously justified
deterministic approximation. For the special case of the SI model, these
yielded a detailed analysis of the final outcome of the epidemic, and in
particular a proof of a necessary and sufficient condition for the final size
to be discontinuous at the phase transition $\lambda=\lambda_{c}$. Similar
results were also obtained for the SIR model when rewiring is necessarily to a
susceptible individual. Moreover in that case the behaviour at the phase
transition is very striking in that the final fraction infected by the
epidemic jumps from $0$ to $1$. For the original SIR model only partial
results were proved, in that there is a non-negligible gap between the
necessary and sufficient conditions for a discontinuity, although stronger
results were conjectured together with supporting evidence. Proving those
conjectures is a worthwhile future work.

It seems likely that similar results hold for epidemics on networks
constructed using the configuration model. Indeed [Yao and Durrett(2020)]
provides such results although, as described in Section 4, we believe that
their results are for a model which is different from the original SIR model
and does not provide an exact asymptotic approximation to the original model
in the limit as $n\to\infty$. However it seems hard to extend our construction
to epidemics on configuration-model networks since the construction exploits
symmetries and independence properties of a Erdős-Rényi random graph that are
not present in configuration-model networks. Although not considered here, it
seems likely that the construction can be extended to models that are more
akin to the Erdős-Rényi random graph, such as the stochastic block model.

##  Appendix A Time-transformed deterministic SI model.

###  A.1 Solution of time-transformed deterministic model (6.27)-(6.29)

We solve (6.27)-(6.29) with the initial condition
$(\tilde{s}(0),\tilde{i}_{E}(0),\tilde{w}(0))=(1-\epsilon,\mu\epsilon(1-\epsilon),0)$.
The solution (6.30) for $\tilde{s}(t)$ follows immediately from (6.27) as
$\tilde{x}(0)=1-\epsilon$. Dividing (6.29) by (6.29) yields

$\dfrac{d\tilde{w}}{d\tilde{s}}=-\frac{\omega\alpha}{\lambda}\tilde{s}+2\frac{\tilde{w}}{\tilde{s}},$

which with the above initial condition yields, using the integrating factor
$\tilde{s}^{-2}$,

$\tilde{w}(\tilde{s})=\frac{\omega\alpha}{\lambda}\tilde{s}^{2}\log\left(\frac{\tilde{s}}{1-\epsilon}\right),$
(A.1)

and (6.32) follows. Dividing (6.28) by(6.29) and substituting from (A.1) gives

$\dfrac{d\tilde{i}_{E}}{d\tilde{s}}=1-\mu\tilde{s}+\frac{\tilde{i}_{E}}{\tilde{s}}+2\frac{\omega\alpha}{\lambda}\tilde{s}\frac{\tilde{w}}{\tilde{s}}\log\left(\frac{\tilde{s}}{1-\epsilon}\right)+\frac{\omega}{\lambda}(1-\alpha+\alpha\tilde{s}),$

which, with the above initial condition can be solved using the integrating
factor $\tilde{s}^{-1}$ to yield (6.31).

###  A.2 Final size and $\tilde{i}_{E}^{\prime}(\tilde{\zeta})$ when
$\epsilon>0$

Define $f_{\epsilon}:[\epsilon,1)\to\mathbb{R}$ by

$f_{\epsilon}(x)=\log\left(\frac{1-x}{1-\epsilon}\right)+\frac{(\lambda\mu+\omega\alpha)x-\omega\alpha\epsilon}{\lambda+\omega(1-\alpha)+2\omega\alpha(1-x)},$

so $\tau$ satisfies $F_{\epsilon}(\tau)=0$ if and only if
$f_{\epsilon}(\tau)=0$. From (6.31), we have
$\tilde{i}_{E}(t)=H(\tilde{s}(t))$, where

$H(\tilde{s})=\frac{\tilde{s}}{\lambda}G(\tilde{s}),\text{ with
}G(\tilde{s})=(\lambda+\omega(1-\alpha)+2\omega\alpha\tilde{s})f_{\epsilon}(1-\tilde{s}).$

Now $f_{\epsilon}(\epsilon)>0$ and $f_{\epsilon}(1-)=-\infty$, so
$f_{\epsilon}$ has at least one root in $[\epsilon,1)$, and since
$\tilde{s}(t)=1-\epsilon-t$ $(t\geq 0)$,
$\tilde{\zeta}=\inf\\{t>0:\tilde{i}_{E}(t)=0\\}$ is finite and
$\tilde{s}(\tilde{\zeta})=1-\tau$, where $\tau$ is the smallest solution of
$F_{\epsilon}(x)=0$ in $(\epsilon,1)$. Further, since
$\tilde{s}^{\prime}(\tilde{\zeta})=-1$,

$\displaystyle\tilde{i}_{E}^{\prime}(\tilde{\zeta})$
$\displaystyle=-H^{\prime}(\tilde{s}(\tilde{\zeta}))$
$\displaystyle=-\frac{1}{\lambda}\left[\tilde{s}(\tilde{\zeta})G^{\prime}(\tilde{s}(\tilde{\zeta}))+G(\tilde{s}(\tilde{\zeta}))\right]$
$\displaystyle=-\frac{\tilde{s}(\tilde{\zeta})}{\lambda}G^{\prime}(\tilde{s}(\tilde{\zeta})),$

as $G(\tilde{s}(\tilde{\zeta}))=0$. Now,

$G^{\prime}(\tilde{s}(\tilde{\zeta}))=-(\lambda+\omega(1-\alpha)+2\omega\alpha\tilde{s}(\tilde{\zeta}))f_{\epsilon}^{\prime}(1-\tilde{s}(\tilde{\zeta})),$

since $f_{\epsilon}(1-\tilde{s}(\tilde{\zeta}))=0$, so
$\tilde{i}_{E}^{\prime}(\tilde{\zeta})<0$ if and only if
$f_{\epsilon}^{\prime}(\tau)<0$. It is shown easily that
$f_{\epsilon}^{\prime}(\tau)<0$ if and only if
$F_{\epsilon}^{\prime}(\tau)<0$.

###  A.3 Final size and $\tilde{i}_{E}^{\prime}(\tilde{\zeta})$ when
$\epsilon=0$

Setting $\epsilon=0$ in (6.31) shows that $\tau=1-\tilde{s}(\tilde{\zeta})$ is
given by the smallest solution in $(0,1)$ of $F_{0}(x)=0$, provided there is
at least one solution, otherwise $\tau=1$. We show now that when $R_{0}>1$
there is a unique solution in $(0,1)$. Let $f$ be the function defined at
(6.35), so $\tau$ satisfies $F_{0}(\tau)=0$ if and only if $f(\tau)=0$. Now

$f^{\prime}(x)=-\frac{1}{1-x}+\frac{(\lambda\mu+\omega\alpha)[\lambda+\omega(1+\alpha)]}{[\lambda+\omega(1-\alpha)+2\omega\alpha(1-x)]^{2}},$

so, for $x\in(0,1)$, $f^{\prime}(x)=0$ if and only if $g_{1}(1-x)=g_{2}(1-x)$,
where $g_{1}(x)=(\lambda+\omega(1-\alpha)+2\omega\alpha x)^{2}$ and
$g_{2}(x)=(\lambda\mu+\omega\alpha)[\lambda+\omega(1+\alpha)]x$. Now
$g_{1}(0)>g_{2}(0)$, $\lim_{x\to\infty}g_{1}(x)-g_{2}(x)=\infty$ and a simple
calculation shows that $g_{1}(1)<g_{2}(1)$ if and only if
$\lambda>\frac{\omega}{\mu-1}=\lambda_{c}$. Suppose $R_{0}>1$, so
$\lambda>\lambda_{c}$. Then $g_{1}-g_{2}$ has a root in $(0,1)$ and a root in
$(1,\infty)$. Moreover, $g_{1}-g_{2}$ has precisely one root in each of these
intervals, as it convex on $\mathbb{R}$. Hence $f^{\prime}$ has a unique root
in $(0,1)$. Further, $f(0)=0$,
$f^{\prime}(0)=-1+(\lambda\mu+\omega\alpha)/[\lambda+\omega(1+\alpha)]>0$, as
$\lambda>\lambda_{c}$, and $f(1-)=-\infty$, so the unique stationary point of
$f^{\prime}$ in $(0,1)$ must be a maximum. It follows that $f(x)=0$, and hence
also $F_{0}(x)=0$ has a unique solution $\tau$ in $(0,1)$ and
$f^{\prime}(\tau)<0$. Letting $\epsilon=0$ in the argument in Appendix (A.2)
shows that $\tilde{i}_{E}^{\prime}(\tilde{\zeta})<0$.

## Acknowledgements

Tom Britton is grateful to the Swedish Research Council (grant 2015-05015) for
financial support.

## References

  * [Altmann (1998)] Altmann, M. (1998) The deterministic limit of of infectious disease models with dynamic partners. Math. Biosci. 150, 153–175.
  * [Andersson and Britton(2000)] Andersson, H. and Britton, T. (2000) Stochastic Epidemic Models and Their Statistical Analysis. Springer, New York.
  * [Ball and Britton(2007)] Ball, F. and Britton, T. (2007) An epidemic model with infector-dependent severity. Adv. Appl. Prob. 39, 949–972.
  * [Ball and Neal(2008)] Ball, F. and Neal, P. (2008) Network epidemic models with two levels of mixing. Math. Biosci. 212, 69–87.
  * [Barbour and Reinert(2013)] Barbour, A. and Reinert, G. (2013) Approximating the epidemic curve. Electron. J. Probab. 18(54), 1–30.
  * [Bollobás(1980)] Bollobás, B. (1980) A probabilistic proof of an asymptotic formula for the number of labelled regular graphs. European J. Combin. 1, 311-316.
  * [Britton and O’Neill(2002)] Britton, T. and O’Neill, P.D. (2002) Bayesian inference for stochasric epidemics in populations with random social structure. Scand. J. Stat. 29, 375–390.
  * [Britton and Trapman(2012)] Britton, T. and Trapman, P. (2012) Maximizing the size of the giant. J. Appl. Prob. 49, 1156–1165.
  * [Chinazzi et al.(2020)] Chinazzi, M., Davis, J.T., Ajelli, M., Gioannini, C., Litvinova, M., Merler, S., Pastore y Piontti, A., Mu, K., Rossi, L., Sun, K., Viboud, C., Xiong, X., Yu, H., Halloran, M.E., Longini, I.M. and Vespignani, A. (2020) The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak. Science 368, 395–400.
  * [Darling and Norris(2008)] Darling, R.W.R. and Norris, J.R. (2008) Differential equation approximations for Markov chains. Probability Surveys 5, 37–79.
  * [Durrett(2007)] Durrett, R. (2007) Random Graph Dynamics. Cambridge University Press, Cambridge.
  * [Ethier and Kurtz(1986)] Ethier, S.N. and Kurtz, T.G. (1986) Markov Processes: Characterization and Convergence. Wiley, New York.
  * [Ferguson et al.(2006)] Ferguson, N.M., Cummings, D.A.T., Fraser, C., Cajka, J.C., PC Cooley P.C. and Burke, D.S. (2006) Strategies for mitigating an influenza pandemic. Nature 442, 448-452.
  * [Funk et al.(2010)] Funk, S., Salathe, M., and Jansen, V.A.A. (2010) Modelling the influence of human behaviour onthe spread of infectious diseases: a review. J. R. Soc. Interface, 7:1247–1256, doi:10.1098/rsif.2010.0142.
  * [Gross et al.(2006)] Gross, T., Dommar D’Lima, C. J. and Blasius, B. (2006) Epidemic dynamics on an adaptive network. Phys. Rev. Letters 96, 208701.
  * [van der Hofstad(2016)] van der Hofstad, R. (2016) Random Graphs and Complex Networks Volume 1. Cambridge University Press, Cambridge.
  * [Jacobsen et al.(2018)] Jacobsen, K.A., Burch, M. G., Tien, J. H. and Rempala, G. A. (2018) The large graph limit of a stochastic epidemic model on a dynamic multilayer network. sl J. Biol. Dyn. 12, 746–788.
  * [Janson(2009)] Janson, S. (2009) The probability that a random multigraph is simple. Combin. Prob. Comput. 18, 205–225.
  * [Janson et al.(2014)] Janson, S., Luczak, M. and Windridge, P. (2014) Law of large numbers for the SIR epidemic on a random graph with given degrees. Random Structures Algorithms 45, 726–763.
  * [Jiang et al.(2019)] Jiang, Y., Kassem, R., York, G., Junge, M. and Durrett, D. (2019) SIR epidemics on evolving graphs. arXiv:1901.06568v1.
  * [Kimber(1983)] Kimber, A.C. (1983) A note on Poisson maxima. Probability Theory and Related Fields 63, 551–-552.
  * [Leung et al.(2018)] Leung, K., Ball, F., Sirl, D. and Britton, T. (2018). Individual preventive social distancing during an epidemic may have negative population-level outcomes. J. R. Soc. Interface 15:20180296.
  * [Martin-Löf(1986)] Martin-Löf, A. (1986) Symmetric sampling procedures, general epidemic processes and their threshold limit theorems. J. Appl. Prob. 23, 265–282.
  * [Molloy and Reed(1995)] Molloy, M. and Reed, B. (1995) A critical point for random graphs with a given degree sequence. Random Structures Algorithms 6, 161–179.
  * [Neal(2003)] Neal, P. (2003) SIR epidemics on a Bernoulli random graph. J. Appl. Prob. 23, 265–282.
  * [Nerman(1981)] Nerman, O. (1981) On the convergence of supercritical general (C-M-J) branching processes. Z. Wahrscheinlichkeitsth 57, 365–395.
  * [Newman et al.(2001)] Newman, M.E.J., Strogratz, S.H. and Watts, D.J. (2001) Random graphs with arbitrary degree distributions and their applications. Phys. Rev. E 64, 026118.
  * [Watson(1980)] Watson, R. (1980) A useful random time-scale transformation for the standard epidemic model. J. Appl. Prob. 17, 324–332.
  * [Yao and Durrett(2020)] Yao, D. and Durrett, R. (2020) Epidemics on evolving graphs. arXiv:2003.08534v1.