# Dandelion-Picking Legged Robot

Sandilya Sai Garimella<EMAIL_ADDRESS>Department of Mechanical
Engineering, University of Michigan, Ann Arbor  Shai Revzen<EMAIL_ADDRESS>Department of Electrical Engineering and Computer Science, Department of
Ecology and Evolutionary Biology, University of Michigan, Ann Arbor

###### Abstract

Agriculture is currently undergoing a robotics revolution, but robots using
wheeled or treads suffer from known disadvantages: they are unable to move
over rubble and steep or loose ground, and they trample continuous strips of
land thereby reducing the viable crop area. Legged robots offer an
alternative, but existing commercial legged robots are complex, expensive, and
hard to maintain. We propose the use of multi-legged robots using low-degree-
of-freedom (low-DoF) legs and demonstrate our approach with a lawn pest
control task: picking dandelions using our inexpensive and easy to fabricate
BigANT robot. For this task we added an RGB-D camera to the robot. We also
rigidly attached a flower picking appendage to the robot chassis. Thanks to
the versatility of legs, the robot could be programmed to perform a “swooping”
motion that allowed this 0-DoF appendage to pluck the flowers. Our results
suggest that robots with six or more low-DoF legs may hit a sweet-spot for
legged robots designed for agricultural applications by providing enough
mobility, stability, and low complexity.

###### Index Terms:

agricultural robotics, mechanical implement, robotic picking, multi-legged
robots, mobile manipulation, computer vision

##  I introduction

The conventional approach to harvesting agricultural produce, such as fruit,
is to use large harvesting machines to produce mechanical vibration to
separate the fruit from the stalk or to use human labor for manual picking.
Both these methods have downsides: the former method fails to reap a third of
edible produce harvested (33.7% market yield remains on the farms) [1] and
there is no guarantee that the plants or picked fruit remain intact due to
damage from machinery [1]. The latter method requires trained pickers and
involves health hazards such as lacerations [2]. Existing wheeled harvesters
such as the Octinion Rubion and the strawberry-picking system from Dogtooth
Technologies require produce to be cultivated in tailored environments like
raised beds [3] and they cannot be used in a field. Other wheeled robotic
pickers depend on sophisticated mechanical design, for example, the apple
picker from FFRobotics operates between 4 and 12 robotic arms and requires
human supervision [3] and the wheeled raspberry harvesting robot from Field
Robotics relies on 4 robotic arms [3]. Purchasing and maintaining several
agricultural machines and robotic pickers is expensive, requires training, and
they are not easily transportable, thus dissuading farmers with small farms
(defined by the USDA as farms with gross cash income under $250,000 annually)
from using them [4].

Any machine designed to support agriculture must either be the size of the
agricultural fields themselves or be able move through them. Conventionally,
such mobility employs large wheels or treads to push against the ground –
thereby regularly trampling a substantial fraction of the potential growing
area. Additionally, wheels and treads encounter trouble when the ground is
soft or flowing, and cannot move over very sharp inclines or discontinuous
terrain (like stepping stones).

Robotic platforms with legs might be a means to avoiding these disadvantages.
Unfortunately, these typically have 3 or more DoF per leg, making even a
quadruped require a minimum of 12 motors. Here we provide an alternative to
such complex multi-legged robots using BigANT, a 6-legged (hexapedal) robot
with only one DoF per leg. Thus, BigANT has only half as many motors as the
typical quadruped robot, yet is far easier to stabilize even when running at
speed - because it is statically stable. The superior value proposition of
multi-legged robots with low-DoF per leg was recently presented in [5, 6],
where the author studied the design, modeling, and control of such robots.
This class of robots offer the advantages of legged systems without the
complexity and cost associated with having numerous actuators.

As a model problem, we considered the task of removing dandelions from a lawn.
Dandelions are a known lawn pest, and they reproduce extremely quickly.
Because their bright yellow flowers are easy to detect on the green background
of a lawn, problems of machine vision and target acquisition are minimized. We
addressed the dandelion-picking problem by generating a series of actions from
a library of existing behaviors (walking, steering, and turning in place),
some of which are parametric. We chose which behavior to perform based on the
RGB-D data which we first reduced to the dandelion’s azimuth, elevation, and
distance. Below we provide background on the BigANT robot and its behaviors
(section I-A), followed by an investigation of the techniques used to cut the
stem of a dandelion (section I-B).

Figure 1: BigANT with cutting appendage attached to the front of the chassis.

###  I-A Previous work on the BigANT hexapedal robot

This project was done entirely on the BigANT robot (Figure 1) which has six
legs, each with 1-DoF. Its properties are described in detail in [5]. To walk
the robot we drove the six legs in an “alternating tripod gait”, with the left
tripod comprising of [F]ront-[L]eft, [M]iddle-[R]ight, and [H]ind-[L]eft legs
(FL-MR-HL), and the right tripod being FR-ML-HR [5]. The toe trajectory
(Figure 2) occupies a one-dimensional manifold with respect to the body frame
of reference and choosing where the leg pairs (F, M, and H pairs) should be on
this trajectory formed the basis of the overall BigANT chassis ‘swooping’
motion which we developed to pick the dandelions. The idea of exploiting the
existing motors to produce the swooping motion came from work on limiting
parasitic vertical chassis oscillations when designing the tripod gait – where
Zhao exploited the timing of leg motions to reduce the vertical motions. By
doing (mostly) the opposite, we exploited the same individual leg trajectories
to lower and raise the robot while pitching, and produce the ‘swooping’ motion
to pick a dandelion.

Figure 2: Toe trajectory in body frame, with points at equal phase intervals.
Each leg pair’s rate modulation and position in the cycle will be used to
define a function for the chassis trajectory. Adapted from [5].

By modifying the gait parameters (defined in [5]) of the mid left and mid
right legs we obtained a steering gait, e.g. slowing down the mid right leg
and speeding up the mid left leg during the ground contact part of the
trajectory results in steering right. We varied a dimensionless quantity
called the turn value $T$, which lies in the range $-0.3\leq T\leq 0.3$, to
produce steering with a particular turn radius.

###  I-B Methods for cutting a dandelion stem: chopping or slicing?

The ways in which a dandelion stem is cut can be distilled to three main
categories: separation due to tensile fracture, applied normal force due to a
blade (chopping), and applied normal and shear forces (slicing). For tensile
fracture to occur, pulling along the stem axis of the dandelion is sufficient.
During chopping the stem undergoes localized normal compression until failure
occurs. Slicing and chopping are dissimilar because the former involves both
normal and shearing deformations whereas the latter involves normal
deformations only, causing global deformations [7]. A soft solid such as a
dandelion stem offered greater resistance to failure under compressive stress
than it did to tensile stress so the slicing technique was more effective due
to localized deformation. Considering which cutting action to use was crucial
because to produce any of the three cutting actions, the cutting device will
need to have at least 1 DoF, either innately or by using BigANT’s DoFs. We
explored the combination of two separation techniques, i.e. tensile stress
with either chopping or slicing by manually testing cutting devices which we
fabricated using rapid prototyping.

We tested the chopping idea by rapidly prototyping a box-like enclosure with a
utility knife blade attached to the closing boundary of one of the parts. Once
the dandelion stem was cut, it remained engulfed by the box, which was closed.
The chopping action occurred when the blade pressed against the closing
boundary of the opposite piece, with the dandelion stem trapped in between;
the closing piece with the knife blade was actuated by a servo motor.

To test the slicing idea, we fabricated a ‘V’ shaped cutting appendage using
two utility knife blades fixed to a foam core frame. Rather than predicting
the performance theoretically, we examined the design by slicing dandelions,
and made improvements through two design iterations (Figure 3). The testing
method was to move the cutting appendage through a concave-up trajectory by
hand (Figures4, 4, and 4) to simulate the BigANT producing the same motion via
its chassis ‘swooping’ trajectory. Through our trials, we identified that
dandelions were being sliced when the cutting appendage approached them at
azimuthal angles near its center (near dashed line in Figure 3). By widening
the two opening sides of the ‘V’ shape, we increased the azimuth range in
which the dandelion was picked successfully. This increased the error
allowance for the swooping trajectory and dandelion position to intersect
after steering towards the dandelion. With repeated use, we improved the first
design’s flaws in the second iteration, namely the lack of restraint for the
dandelion after it is plucked and the cut dandelion occasionally falling
forwards. The second cutting appendage had foam core barriers to prevent the
cut dandelion from falling sideways out of the platform, and a positive
inclination with respect to the BigANT’s chassis base to make it fall
backwards consistently (Figure 4).

Figure 3: (a) Comparison of the first cutting appendage (top) and second
design. Note the improvements in the second iteration: positive inclination
and side barriers for consistent plucking. In (b), the top view of ‘V’ shape
cutting appendage with a dandelion stem of radius $r$ approaching the vertex.

Although the ‘V’ shape cutting appendage has 0-DoF on its own, it harnesses
the BigANT’s forward motion to produce normal and shear stresses required for
slicing. As the dandelion stem approaches the vertex of the ‘V’, the maximum
stem radius which the blades can accommodate decreases. The relationship
between the maximum radius $L$, distance of the stem center to the vertex $d$,
and half the interior angle of the ‘V’ $\theta$, is given by
$L=d\sin(\theta)$. The condition for slicing is satisfied when the stem center
is at critical distance $d_{cr}$ from the vertex; then
$d_{cr}=r/\sin(\theta)$, where $r$ is the stem radius (shown in green in
Figure 3).

Depending on the height of the dandelion, either the tensile fracture or
slicing occurred. Suppose the dandelion height was greater than that of the
cutting appendage blades above ground, slicing occurred because only the stem
interacted with the blades. If the dandelion was short, its head was
positioned above the blades near the vertex until the final stage of the
BigANT’s ‘swooping’ trajectory when enough tensile stress was applied for
separation to occur. This phenomenon is demonstrated in figures 4 and 4.

Figure 4: Approaching the dandelion in early swoop trajectory in (a). The
dandelion head is caught firmly at the vertex in (b). Due to tensile fracture,
the dandelion is cut in (c).

##  II A chassis fixture for picking dandelions

The chosen cutting appendage would have to be move along an appropriate
trajectory to operate. One option was to attach it to an actuated mechanism
designed to produce the needed motion. Instead, we explored a solution that
attached the cutting appendage to the chassis. The existing DoFs which the
cutting appendage could have due to BigANT’s mobility were forwards and
backwards movement (on BigANT’s sagittal plane), yaw due to steering or
turning in place, roll due to height differences on left and right legs, pitch
through differences in height of front and hind legs, and vertical motion
withing the clearance height of the legs.

By performing identical symmetric motions on right and left contralateral leg
pairs, we reduced the complexity of the trajectory design from 3-D to 2-D.
From here on, in describing the picking maneuver, we will discuss it as if it
is a 2-D problem. The ‘swooping’ motion we designed consists of lowing the
cutting appendage, then moving it forward and up, while intercepting the
dandelion near the nadir of the ‘swooping’ arc. We created motions (e.g. the
cutting appendage trajectory in Figure 5) by using the rapid prototyping
abilities developed in [8]: puppeting the robot legs by hand, recording the
motor positions in a gait table in CSV format, and replaying the that table as
a motion primitive (a Plan in the language of [8]).

Figure 5: Swooping motion of the cutting appendage. The feasible toe
trajectories relative to the current body location are shown in Figure 2. We
combined toe motions to produce the cutting appendage motion.

We divided the concave-up trajectory motion into five stages (Figures 6-6) for
understanding the leg pair motion in each stage. The BigANT began in the
‘slack’ stage where all the legs were near the apex of their trajectories. The
hind leg pair then reached the lowest point of the toe trajectory, putting the
chassis at a negative attitude to the horizontal. The front and hind pair
moved simultaneously in opposite directions to make the cutting appendage
attachment point reach the lowest point of the chassis trajectory. The mid
legs were swiftly lowered while the hind legs pushed slightly forward to
complete the ‘swooping’ motion.

Figure 6: Swooping chassis trajectory stages. (a) shows starting stage;
chassis has negative attitude due to lowering hind legs in (b), reaching the
lowest point of concave-up trajectory in (c) and (d), and raising the cutting
appendage in the final portion of concave-up trajectory in (e).

We determined this picking motion ad-hoc, however, given that the robot has
sufficient actuation to control fore-and-aft motion, height above ground, and
pitch angle simultaneously, it is likely straightforward to design a picking
motion trajectory using some form of inverse kinematics computation.

##  III Computer Vision Based Method for Detecting Dandelions

We used an Intel RealSense L515 RGB-D camera for detecting dandelions to
support motion planning. We mounted the camera at a distance of 0.20m from the
vertex of the cutting appendage on the chassis and we accounted for this in
the ‘distance to dandelion’ parameter in our code, so that BigANT would not
continue walking past the dandelion. We processed the color (Red-Green-Blue)
and depth streams concurrently and we obtained them in the same resolution
(640x480) to avoid frame misalignment issues.

For simplicity, we considered a ‘dandelion’ to be any yellow colored globular
blob in the field of view. This also allowed us to test our code with a small
inflated yellow balloon and focus our attention on BigANT’s robotic behaviors,
while setting us a computer vision framework which future schemes could
extend.

Our code used the OpenCV 4.5.2.54 and NumPy 1.21.0 Python libraries. We
identified the ‘dandelion’ using hue and saturation thresholding in a hue-
saturation-value (HSV) representation of the image (Figure 7). We used moments
of inertia along the vertical and horizontal axes to identify the geometric
center of the ‘dandelion’. To distinguish between yellow spherical objects and
yellow non-spherical objects, we defined a circularity variable as $c=4\pi
A/P^{2}$, where $A$ and $P$ denote area and perimeter of the object
respectively. When circularity $c$ was in the range $0.4<c<1.45$, the yellow
object satisfied the condition for being a ‘dandelion’.

We reduced the location of the ‘dandelion’ to azimuth, elevation, and distance
(Figures 7 and 8) using the centroid pixel location and the fact that the
frame had a field of view of 54° (horizontal) by 40° (vertical). We obtained
the remaining distance to the dandelion from the depth frame by taking the
mean of depth values contained in the bounding box obtained from the RGB
frame.

Figure 7: RGB (top) and depth frames show how the dandelion’s azimuth,
elevation, and depth are tracked; the red dot represents the centroid, the
bounding box is in blue, and the contour is green.

Azimuth, elevation, and distance (referred to as ($a$, $e$, $d$) for
concision) can be used to express the vector from BigANT’s onboard camera
position to the dandelion in spherical coordinates. It can be shown that the
($a$, $e$, $d$) coordinates translate to ($r$, $z$, $\theta$) as
$r=d$cos($\gamma)$, $z=d$sin($\gamma)$, and $\theta=a$, where
$\gamma=\arctan(1/(\cos(a)\tan(e)))$ (Figure 8). For the purpose of navigating
toward the dandelion, we used azimuth $a$ and distance $d$; elevation played a
role in determining whether the dandelion was within picking scope of the
cutting appendage or not.

Figure 8: How azimuth, elevation, and distance ($a$, $e$, $d$) are defined.

##  IV Algorithmic Control Scheme for Motion Planning

The two parameters we used to formulate the path to the dandelion were azimuth
and distance. With steering (based on definition of turning versus steering in
[5]) involved, BigANT’s path resembled an arc on which two points represented
BigANT’s initial position and the target, i.e. the dandelion’s position
(Figure 9). For instructing BigANT to steer for a specified distance, we
analyzed two approaches. The first approach entailed dead reckoning by walking
for a time period calculated from a measured stride length and the gait
frequency.

Figure 9: (a) Given azimuth $a$ and distance $d$, a vector in polar
coordinates can be defined from BigANT’s cutting appendage $B$ and dandelion’s
position $T_{d}$. An arc with center $C$ resembles BigANT’s steering path. (b)
shows an example of successive steering arc paths shown to reach the target
dandelion placed at 18°azimuth and 1.50 m distance. The centers of the arcs
which overshoot the target are $C_{1},C_{2},,C_{3}$.

For example, for a gait frequency of $f=0.16$ Hz, the overall gait cycle time
period was $T=1/f=6.25$ s. We experimentally measured that approximately $3$
cycles were needed to walk $1.00$ m, allowing us to estimate $18.75d$ seconds
of walking to go $d$ meters. The drawback of this method was that the number
of cycles to cover a certain distance was not always the same due to
variability in the walking gait. Nevertheless, this naive approach was useful
for initial prototyping.

We replaced this scheme with a feedback based policy using the parametric
steering gait available for the BigANT [5].

Here too we considered a dead-reckoning solution first. We placed a dandelion
at several azimuth points along an arc which was 1.50 m from BigANT’s cutting
appendage (Figure 9) and explored turn values manually according to whether
the BigANT overshot or undershot the dandelion. Overshooting in this context
refers to BigANT steering past the dandelion’s azimuth such that the steering
direction has to be reversed for it to reduce the difference angle between the
cutting appendage and dandelion azimuth. We produce a calibration curve of
turn value as a function of azimuth such that BigANT slightly overshot the
dandelion while steering. This allowed us to produce successive overshooting
paths, which eventually led to the dandelion (Figure 9) by BigANT autonomously
choosing appropriate turn values.

We considered a second steering control scheme which relied on varying a
steering parameter $s$ in the shaft angle equations (3) and (4) to modify the
phase of the middle legs in an anti-symmetric manner [5]. With azimuth and
distance, we can define a target point like before, and identify an arc
(Figure 8 (right) in [5]) that contains the start and target points. Solving
for the steering input $s$ which corresponds to the arc then gives a turn
value $T$ for steering toward the dandelion. This approach was different
because it identified a turn value $T$ using steering arc results generated by
varying a steering parameter $s$, whereas the first approach chose $T$ values
based solely on the dandelion’s azimuth at a distance of 1.50 m. The strategy
involving steering parameter $s$ may give the ability to steer to distances
other than 1.50 m. Both approaches are limited in that steering results
produced on one surface would vary from those on another surface due to
differences in foot-surface interactions.

$\displaystyle\psi_{FL}$ $\displaystyle=\psi_{HL}:=b(\varphi)$ (1)
$\displaystyle\psi_{FR}$ $\displaystyle=\psi_{HR}:=b(\varphi+1/2)$ (2)
$\displaystyle\psi_{ML}$
$\displaystyle=b(1/2+\varphi+sk_{s}\cos(2\pi\varphi))$ (3)
$\displaystyle\psi_{MR}$ $\displaystyle=b(\varphi-sk_{s}\cos(2\pi\varphi))$
(4)

We used a state machine to control the picking process. This consisted of the
following processes: (1) receive ($a$, $e$, $d$) data to recognize the
dandelion’s position, (2) set the turn parameter for steering based on the
dandelion’s azimuth and walk forwards based on the dandelion’s distance, (3)
inspect dandelion’s position with respect to cutting appendage after
completing step (2), (4.1) perform swooping motion if dandelion is in the
picking zone or (4.2) steer while walking backwards to reapproach missed
dandelion, this time by walking in a straight line towards it. Due to the
nature of the distance covering method, where BigANT continuously adjusts its
turn parameter for steering, processes (1) and (2) occur concurrently.

We continuously monitored the remaining distance to reach the dandelion and
once BigANT was approximately 0.20 m from the dandelion, we stopped walking.
The 0.20 m value was to compensate for the distance between the camera and the
cutting appendage’s vertex.

If the dandelion’s centroid was in a box defined by the azimuthal range -7.0°
to 7.0° and an elevation range of -15.0° to 20.0° of the cutting appendage
(the ‘picking zone’) we performed the swooping motion. We determined that
dandelions can be successfully picked within this zone by repeated trials.

When successive overshooting paths failed to move the cutting appendage such
that the dandelion was within the picking zone, BigANT steered backwards on a
path with the smallest turning radius (turn parameter $|T|=0.3$) just until
the dandelion was visible within the picking zone.We then approached the
dandelion by walking forwards in a straight line path, keeping it in the
picking zone.

##  V Conclusions

To have a net positive economic impact on agriculture, robots can evolve in
one of two ways: become cheap, reliable, and moderately productive, or remain
complex and expensive, but be truly spectacularly productive. Here we focused
on the first of these two approaches.

We have shown how the BigANT hexapod robot’s capabilities could perform an
agricultural pest control task in a lab setting that simulated some features
of dandelion picking in an open field. This illustrates that low-DOF multi-
legged robots can perform meaningful field robotics tasks and future work
should consider how the robot morphologies can be used to reap the benefits of
legs without the cost of leg complexity. The solution we presented did not
involve adding any DoFs to the BigANT or the mechanism (0-DoF cutting
appendage) with which it plucks dandelions, illustrating one of the many
advantages of legged systems – their ability to produce novel motions through
coordinated motion.

Our work here was more an illustration of a principle than a complete
demonstration. One improvement which would likely be necessary for any field
use is the detection of multiple dandelions at the same time, and the ability
to track a target dandelion in a field of other dandelions. These kinds of
computer vision problems are well studied, and we believe existing solutions
could be applied.

More interestingly, future work could explore how to control foot placement
while performing the dandelion picking task. One of the great promises of
legged robots for agricultural uses is the ability of legs to trample far less
of the agricultural fields than wheels do. Demonstrating how this could be
done well with low DOF, low complexity robot morphologies is a promising
follow-up to our work here.

##  VI Acknowledgements

Work on this project was funded by the Wang Chu Chien-Wen Summer Research
Award.

## References

  * [1] G. A. Baker, L. C. Gray, M. J. Harwood, T. J. Osland, and J. B. C. Tooley, “On-farm food loss in northern and central california: Results of field survey measurements,” _Resources, Conservation and Recycling_ , vol. 149, pp. 541–549, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0921344919301296
  * [2] S. A. McCurdy, S. J. Samuels, D. J. Carroll, J. J. Beaumont, and L. A. Morrin, “Agricultural injury in california migrant hispanic farm workers,” _American journal of industrial medicine_ , vol. 44, no. 3, pp. 225–235, 2003\.
  * [3] R. Bogue, “Fruit picking robots: has their time come?” _Industrial Robot: the international journal of robotics research and application_ , 2020.
  * [4] N. Akram, M. W. Akram, and W. Hongshu, “Study on the socioeconomic factors affecting adoption of agricultural machinery,” _Journal of Economics and Sustainable Development, 11 (3)_ , 2020.
  * [5] D. Zhao and S. Revzen, “Multi-legged steering and slipping with low dof hexapod robots,” _Bioinspiration & biomimetics_, vol. 15, no. 4, p. 045001, 2020.
  * [6] D. Zhao, “Locomotion of low-dof multi-legged robots,” Ph.D. dissertation, University of Michigan, 2021.
  * [7] E. Reyssat, T. Tallinen, M. Le Merrer, and L. Mahadevan, “Slicing softly with shear,” _Physical review letters_ , vol. 109, no. 24, p. 244301, 2012.
  * [8] S. Revzen, J. Sastra, N. Eckenstien, and M. Yim, “Ckbot platform for the icra 2010 planetary challenge.”


# Spectroscopy of phase transitions for multiagent systems

Niccolò Zagli<EMAIL_ADDRESS>Department of Mathematics, Imperial
College London, London, UK  Department of Mathematics and Statistics,
University of Reading, Reading, UK  Centre for the Mathematics of Planet
Earth, University of Reading, Reading, UK  Valerio Lucarini  Department of
Mathematics and Statistics, University of Reading, Reading, UK  Centre for the
Mathematics of Planet Earth, University of Reading, Reading, UK  Grigorios A.
Pavliotis  Department of Mathematics, Imperial College London, London, UK

###### Abstract

In this Letter we study phase transitions for weakly interacting multiagent
systems. Using linear response theory we probe the emergence of singular
behaviour in the thermodynamic limit by studying the properties of the real
and imaginary parts of the susceptibility. We find clear evidence of the loss
of analyticity due to a pole crossing the real axis of frequencies. Such
behaviour has a high degree of universality, as it does not depend on either
the applied forcing or on the considered observable. We present results
relevant for both equilibrium and nonequilibrium phase transitions by studying
the Desai-Zwanzig and Bonilla-Casado-Morillo models.

††preprint: APS/123-QED

##  I  Introduction

Agent based models are regularly employed to model various phenomena in the
natural sciences, social sciences and engineering Naldi  _et al._ (2010);
Pareschi and Toscani (2013). Multiagent systems are used to model diverse
phenomena such as cooperation Dawson (1983), synchronisation Acebrón  _et al._
(2005), systemic risk Garnier  _et al._ (2013) and consensus formation Wang
_et al._ (2017); Garnier  _et al._ (2017). They are fundamental in developing
algorithms for sampling and optimization Garbuno-Inigo  _et al._ (2020) and
they have also been used for the management of natural hazard Simmonds  _et
al._ (2019) and climate change impact Geisendorf (2018).

Multiagent systems can often exhibit abrupt changes in their behaviour, often
corresponding to critical transitions that occur when a parameter, e.g.
interaction strength or temperature, passes a certain threshold. Such
transitions are often associated to cataclysmic events such climate change,
market crashes etc Scheffer (2009); Sornette (0006). The importance of
developing tools for predicting critical transitions has long been recognized.
One of the main tools used in order to develop early warning signals for
critical transitions is that of linear response theory.

Following the seminar contribution by Kubo, Kubo (1966), linear response
theory represents a very powerful framework for studying the properties of
statistical mechanical systems by investigating how they respond to external
perturbations. Linear response theory has been successfully applied to classic
problems of solid-state physics and optics Lucarini  _et al._ (2005) as well
as plasma physics and stellar dynamics (Binney and Tremaine, 2008, Ch.5); see
applications of the theory in both equilibrium and out of equilibrium
processes Öttinger (2005); Marconi  _et al._ (2008a); Baiesi and Maes (2013);
Cessac (2019); Sarracino and Vulpiani (2019); Gottwald (2020). Linear response
has also been employed to predict the response of the climate system to
external forcings Leith (1975); North  _et al._ (1993); Gritsun  _et al._
(2008); Lucarini  _et al._ (2017); Bódai  _et al._ (2020); Lembo  _et al._
(2020). Rigorous mathematical foundations for linear response theory have been
provided for the case of Axiom A systems Ruelle (1998, 2009) and for diffusion
processes, both in finite and in infinite dimensions Dembo and Deuschel
(2010); Hairer and Majda (2010); see also the interesting viewpoint proposed
by Wormell and Gottwald (2019).

Critical transitions arise when the the spectral gap of the transfer operator
of the unperturbed system shrinks to zero Liverani and Gouëzel (2006);
Chekroun  _et al._ (2014); Lucarini (2016); Tantet  _et al._ (2018) as the
Ruelle-Pollicott poles Pollicott (1985); Ruelle (1986), touch the real axis.
Near criticality, the negative feedbacks of the system become increasingly
ineffective, resulting in arbitrarily large, usually non-Gaussian,
fluctuations and a divergence of correlation properties of the system Dawson
(1983); Shiino (1987); Delgadino  _et al._ (2020).

In the thermodynamic limit, multiagent system can also undergo a qualitative
change of their properties through a different mathematical mechanism, namely
phase transitions Dawson (1983); Shiino (1987), defined as exchange of
stability of nonunique stationary distributions as the parameters of the
systems vary; see a detailed analysis in Carrillo  _et al._ (2020).

In a previous paper Lucarini  _et al._ (2020) we derived linear response
formulas for a system of weakly interacting diffusions described by an
$N-$particle Fokker-Planck equation and have explicitly identified two
qualitatively different scenarios for the breakdown of the linear response,
associated with the previously mentioned critical transitions and phase
transitions. We focus here on the latter case.

Phase transitions are a genuine thermodynamic phenomenon, where the divergence
of the response stems from the coordination taking place, in suitable
conditions, because of the coupling between the infinite number of agents
composing the total system. The coupling among the subsystems results in a
memory effect that leads to obtaining the macroscopic response function of the
system as a renormalised version of its microscopic counterpart Lucarini  _et
al._ (2020), with formal similarities with the well-known Clausius-Mossotti
relation Lucarini  _et al._ (2005); Jackson (1975); Talebian and Talebian
(2013). The role of memory in determining criticality due to endogenous
processes has been emphasised in Sornette and Helmstetter (2003); Sornette
(0006). The link between phase transitions and slow decay of correlations for
interacting particle systems is well established, see. e.g. Yoshida (2003).

###  I.1 This Letter

We relate the occurrence of phase transitions to the emergence of non-
analytical properties of the susceptibility of the system describing the
frequency-dependent response of an observable to a given perturbation due to
the emergence of a pole in the real axis of the frequency. We use a formalism
that closely mirrors spectroscopic techniques that are commonly used for
investigating the optical properties of matter Lucarini  _et al._ (2005). The
presence of the pole becomes very clear as we increase the number of agents,
thus approaching the thermodynamic limit. The attractiveness of the approach
proposed here is that the loss of analyticity depends neither on the choice of
the observable, nor on the applied perturbation (except degenerate, non-
typical cases), but is instead universal in nature Lucarini  _et al._ (2020).
We study the Desai-Zwanzig (DZ) Desai and Zwanzig (1978) and the Bonilla-
Casado-Morillo (BCM) Bonilla  _et al._ (1987) models. The DZ model exhibits a
paradigmatic example of an equilibrium order-disorder phase transition,
analogous to the Ising ferromagnetic transition Shiino (1987); Dawson (1983),
while the BCM model describes an out-of-equilibrium synchronisation transition
of an infinite collection of coupled nonlinear oscillators. As the transition
point is crossed, the order parameter (magnetization) acquires a non vanishing
constant value for the DZ model and is periodically oscillating for the BCM
model.

##  II The general framework

We investigate a system composed of $N$ exchangeable interacting $M$
dimensional sub-systems whose dynamics is determined by the following Itô
stochastic differential equations (SDEs)

$\mathrm{d}\mathbf{x}^{k}=\mathbf{F}_{\alpha}(\mathbf{x}^{k})\mathrm{d}t-\frac{\theta}{N}\sum_{l=1}^{N}\mathbf{\nabla}\mathcal{U}\left(\mathbf{x}^{k}-\mathbf{x}^{l}\right)\mathrm{d}t+\sigma\mathbf{S}(\\{\mathbf{x}^{k}\\})\mathrm{d}\mathbf{W},$
(1)

where $\mathbf{x}^{k}\in\mathbb{R}^{M}$ and $k=1,\ldots,N$.
$\mathbf{F_{\alpha}}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{M}$ is a smooth
vector field, possibly depending on a parameter $\alpha$, and $\mathbf{W}$
denotes a standard $P-$dimensional Brownian motion;
$\mathbf{S}:\mathbb{R}^{M}\rightarrow\mathbb{R}^{M\times P}$ is the volatility
matrix and the parameter $\sigma>0$ controls the strength of the stochastic
forcing, i.e. plays the role of the temperature. We consider a fully coupled
system given by the quadratic (Curie-Weiss) interaction potential
$\mathcal{U}\left(\mathbf{y}\right)=\frac{|\mathbf{y}|^{2}}{2}$. In this case,
the order parameter is known and it is given by the first
moment/magnetization.

The coefficient $\theta$ modulates the intensity of the coupling, which
attempts at synchronising all systems by attracting them to the center of
mass. In the thermodynamic limit $N\rightarrow\infty$ the one-particle
distribution function converges to the distribution $\rho(\mathbf{x},t)$ that
satisfies a nonlinear and nonlocal Fokker-Planck equation Dawson (1983);
Sznitman (1989); Oelschlager (1984); Dawson and Gärtner (1987)

$\begin{split}\partial_{t}\rho(\mathbf{x},t)=&-\nabla\cdot\left[\rho(\mathbf{x},t)\left(\mathbf{F}_{\alpha}(\mathbf{x})+\theta\left(\langle\mathbf{x}\rangle(t)-\mathbf{x}\right)\right)\right]\\\
&+\frac{\sigma^{2}}{2}\nabla\cdot\nabla\cdot(\mathbf{D}\rho(\mathbf{x},t)),\end{split}$
(2)

where $\mathbf{D}=\mathbf{S}\mathbf{S}^{T}$. This mean field PDE, might
support multiple coexisting stationary measures at low temperatures/large
interaction strengths. In particular, in a conservative system described by a
confining potential $\mathbf{F_{\alpha}}(\mathbf{y})=-\nabla
V_{\alpha}(\mathbf{y})$ with additive noise $\mathbf{S}=\mathbf{I}$, the
identity matrix, stationary solutions of Eqn. 2 correspond to local minima of
$V_{\alpha}(\mathbf{y})$. In this case, the thermodynamic limit (2) can be
written in the standard form $\lim_{N\rightarrow+\infty}-\frac{1}{N}\log
Z_{N}=\inf F(\rho)$, where $Z_{N}$ denotes the partition function of the
$N-$particle system and $F(\rho)$ denotes the free energy of mean field system
Helffer (2002). A stationary state is characterised by the order parameter
$\langle\mathbf{x}\rangle_{0}$ and the associated stationary distribution
$\rho_{0}(\mathbf{x})$.

The Fokker-Planck operator that appears on the right hand side of (2)
$\mathcal{L}_{\langle\mathbf{x}\rangle_{0}}$ (evaluated at the stationary
state $\rho_{0}(\mathbf{x})$) can be interpreted as the adjoint of the
generator of the Koopman operator of the stationary dynamics Frank (2004);
Lucarini  _et al._ (2020). As such, correlation properties of the system are
related to the spectrum of $\mathcal{L}_{\langle\mathbf{x}\rangle_{0}}$
Chekroun  _et al._ (2020, 2014).

We now perturb the stationary state by setting
$\mathbf{F}_{\alpha}(\mathbf{x})\rightarrow\mathbf{F}_{\alpha}(\mathbf{x})+\varepsilon\mathbf{X}(\mathbf{x})T(t)$
and we study the response of the system by expanding the distribution function
as
$\rho(\mathbf{x},t)=\rho_{0}(\mathbf{x})+\varepsilon\rho_{1}(\mathbf{x},t)+O(\varepsilon^{2})$.
The response of the order parameter in the frequency domain is written in
terms of a macroscopic (or renormalised) susceptibility
$\tilde{\Gamma}_{i}(\omega)$ as 111Repeated indices are summed.

$\langle
x_{i}\rangle_{1}(\omega)=P_{ij}^{-1}(\omega)\Gamma_{j}(\omega)T(\omega)\equiv\tilde{\Gamma}_{i}(\omega)T(\omega)$
(3)

where $P_{ij}(\omega)=\delta_{ij}-\theta\mathcal{Y}_{ij}(\omega)$ and the
susceptibilities $\Gamma_{i}(\omega)$and $\mathcal{Y}_{ij}(\omega)$ are
respectively the Fourier Transform of the microscopic response functions that
can be written as correlation functions in the unperturbed state as

$\displaystyle G_{i}(\tau)$
$\displaystyle=-\Theta(\tau)\left\langle\frac{\nabla\cdot\left(\rho_{0}(\mathbf{y})\mathbf{X}(\mathbf{y})\right)}{\rho_{0}(\mathbf{y})}\exp\left(\mathcal{L}^{+}_{\langle\mathbf{x}\rangle_{0}}\tau\right)y_{i}\right\rangle_{0}$
(4) $\displaystyle Y_{ij}(\tau)$
$\displaystyle=-\Theta(\tau)\left\langle\partial_{y_{j}}\log\rho_{0}(\mathbf{y})\exp\left(\mathcal{L}^{+}_{\langle\mathbf{x}\rangle_{0}}\tau\right)y_{i}\right\rangle_{0}$
(5)

where $\mathcal{L}^{+}_{\langle\mathbf{x}\rangle_{0}}$represents the adjoint
of $\mathcal{L}_{\langle\mathbf{x}\rangle_{0}}$ and $\langle\cdot\rangle_{0}$
is the expectation value on the unperturbed state $\rho_{0}(\mathbf{x})$. The
renormalisation of the susceptibility derives from the coupling among the
subsystems; note that $\tilde{\Gamma}_{i}(\omega)$ inherits the poles of both
$\Gamma_{i}(\omega)$ and of the matrix $P^{-1}_{ij}(\omega)$.

Away from criticality, both the microscopic and macroscopic susceptibilities
are analytic in the upper half of the $\omega$ complex plane.

Figure 1: Renormalised response functions as a function of time. Left panel:
response $\tilde{G}(\tau;N)$ for the one dimensional order parameter of the DZ
model. Right panel: response $\tilde{G}_{x}(\tau;N)$ of the first component of
the bi-dimensional order parameter for the BCM model. Black and blue lines
correspond to non critical values of the strength of the noise $\sigma$. Red
lines correspond to response functions at the transition point. For each value
of $\sigma$, there are five lines corresponding to different values of $N$,
namely $N=2^{k}\times 10^{3}$ with $k=1,\dots,5$. The arrows and the colour
gradient indicate the direction of increasing $N$

Figure 2: Macroscopic susceptibilities as a function of the frequency
$\omega$. Left panel: susceptibility $\tilde{\Gamma}(\omega)$ for the one
dimensional order parameter of the DS model. Right panel: susceptibility
$\tilde{\Gamma}_{x}(\omega)$ for the first component of the two dimensional
order parameter for the BCM model. The parameter in the lower extreme of the
integral is for both cases $c=0.05$. Blue and black lines in the left panel
have been multiplied by a factor $5$ for visualisation purposes. Colour code
and plotting conventions as in Figure 1.

It is of interest here the case where poles appear in the real $\omega$ axis
for $\tilde{\Gamma}_{i}(\omega)$ because the matrix $P_{ij}(\omega)$ becomes
singular. As discussed in Lucarini  _et al._ (2020), this corresponds to phase
transitions originated by the coupling and do not show a divergence of
correlation properties, because $\Gamma_{i}(\omega)$ is, instead, analytic in
the upper complex $\omega$ plane. At the transition point, the usual
dispersion relations need to be modified Lucarini  _et al._ (2020, 2005). The
conditions underpinning the breakdown of linear response theory do not depend
on the perturbation field $\mathbf{X}(\mathbf{x})$ nor on the choice of the
observable ($\mathbf{x}$, in our case) and can be related to the spectral
properties of a modified transfer operator Lucarini  _et al._ (2020).

##  III Numerical results

Below, we present results for the Desai-Zwanzig (DZ) Desai and Zwanzig (1978)
and the Bonilla-Casado-Morillo (BCM) Bonilla  _et al._ (1987) models. These
are composed by $N$ interacting agents evolving according to Eq. 1. Each
individual agent of the DZ (BCM) model evolves in $\mathbb{R}$
($\mathbb{R}^{2}$). A detailed description of the two models is presented in
the Supplementary Material. We repeat our experiments for various choices of
$N$, in order to detect the emergence of singularities for the combination of
the parameters corresponding to phase transitions. Here, we keep fixed the
values of the internal parameter $\alpha$ and of the coupling intensity
$\theta$ and we vary, instead, the noise intensity $\sigma$. The phase
transition is obtained for a specific value
$\tilde{\sigma}=\sigma(\theta,\alpha)$ that can be calculated analytically for
both models, see Supplementary Material.

Following Marconi  _et al._ (2008b), we perform $n$ simulations where the
initial conditions are chosen according to the unperturbed invariant measure
$\rho_{0}(\mathbf{x})$ and where at time $\tau=0$ we apply a perturbation
proportional to a Dirac $\delta$ function. The average of the response for the
observable $\mathbf{x}$ over the $n$ simulations gives an estimate
$\tilde{G}_{i}(\tau;N)$ of the renormalised response function . Details on the
numerical simulations are also reported in the Supplementary Material. Figure
1 shows the response functions $\tilde{G}_{i}(\tau;N)$ for an additive
perturbation $\mathbf{X}(\mathbf{x})=1$ for the DZ model (left panel) and
$\mathbf{X}(\mathbf{x})=(1,0)$ for the BCM model (right panel). The two
response functions are qualitatively different because, by and large, the one
for the DZ model describes a monotonic decay, whereby the system relaxes
towards the unperturbed state, while the one for the BCM combines the decay
with an oscillatory behaviour taking place at the natural frequency
$\tilde{\omega}=1$ of the unperturbed model.

In the DZ model, the response functions initially undergo a fast and
substantial decay, both far from and at the phase transition, associated with
a time scale of order $1$. However, at the phase transition, a new, much
longer, timescale appears. This timescale increases monotonically with $N$.
The same is observed in the case of the BCM model if one considers the
envelope of the response function rather than the response function itself: at
the transition the decay of the oscillations becomes slower and slower as $N$
increases.
The origin of the new timescales resides in the appearance of simple pole at
$\omega=\omega_{0}$ in the susceptibility $\tilde{\Gamma}_{i}(\omega;N)$, the
Fourier transform of the response function. The pole is located at
$\omega_{0}=0$ for the DZ model and at $\omega_{0}=\tilde{\omega}=1$ for the
BCM model. When considering finite values of $N$, the susceptibilities
describing the response of (virtually) any observable to (virtually) any
external perturbation have a contribution of the form
$\frac{\kappa}{\omega-\omega_{0}+i\gamma(N)}$, where $\gamma(N)\rightarrow
0^{+}$ as $N\rightarrow+\infty$ and $\kappa$ represents the residue of the
pole, because
$\lim_{N\rightarrow\infty}\frac{\kappa}{\omega-\omega_{0}+i\gamma(N)}=-i\pi\kappa\delta(\omega-\omega_{0})+\kappa\mathcal{P}(1/(\omega-\omega_{0}))$.
The quantity $\kappa\in\mathbb{C}$ depends on the choice of observable and of
the perturbation. We next investigate the phase transitions by looking at the
properties of the susceptibilities, see Fig. 2. When
$\sigma\neq\tilde{\sigma}$, the susceptibilities do not show any singularity
nor any remarkable dependence on $N$, thus indicating that the thermodynamic
limit has been reached to a good approximation. As $N$ increases, for both the
DZ model (left panel) and the BCM model (right panel) the resonance at
$\omega=\omega_{0}$ of the real part of the susceptibility approaches the
limiting Dirac function $\pi k\delta(\omega-\omega_{0})$ with coefficient
$k>0$. This singular behaviour is clear from the plot of the primitive
function of the real part of the susceptibility (bottom inset) that tends to
step function. For both models $\kappa=\mathrm{i}k$ is an imaginary number.
Indeed, the imaginary part of the susceptibility behaves exactly as the Cauchy
principal value distribution and can be used to get easily a quantitative
estimate of $k$. The top insets of Figure 2 shows the function
$(\omega-\omega_{0})\mathbf{Im}\\{\tilde{\Gamma}_{i}(\omega)\\}$. As
$N\rightarrow\infty$, this function converges to $k$ everywhere except for
$\omega=\omega_{0}$. An explicit expression for $k$ is known Lucarini  _et
al._ (2020) in the case of the DZ model 222There is a typo in the formula
given in Lucarini  _et al._ (2020). Furthermore, the convention for the
Fourier transform we have used in this paper has the opposite sign, hence
ditto the residue. $k=\frac{\langle
x^{2}\rangle_{0}}{\theta\int_{0}^{+\infty}\mathrm{d}t\langle
x(t)x(0)\rangle_{0}}.$
Using the statistics of the unperturbed runs we obtain $k\approx 0.89$, which
agrees within 2% with the one obtained from the limiting behaviour of the
susceptibility, thus validating our results. In the case of the BCM model, our
procedure allows one to derive a direct estimate $k\approx 0.44$; in this case
no expression for the residue is available in the literature and, following
Bonilla  _et al._ (1987); Lucarini  _et al._ (2020), its evaluation seems
cumbersome. We here observe that, by the evaluation of the susceptibility for
a finite value of $N$, we are able to predict the residue of the pole
$\omega_{0}$ that appears only in the thermodynamic limit. As discussed
earlier, the singular behaviour of the susceptibility is universal; see the
Supplementary Material for an additional example.

##  IV Conclusions

The study of how a large network of identical agents respond to exogenous
perturbations is of the uttermost importance in different areas of science.
One might be interested not only in the smooth response of the system, where
its properties change ever so slightly, but also in the critical, nonsmooth,
regime, where small perturbations can lead to large and possibly undesired
changes. Multiagent systems modeled as weakly interacting Itô diffusions
represent a rich class of models exhibiting such critical behaviour in their
thermodynamic limit, and for which rigorous analysis and careful numerial
investigations can be carried out. Usually these phenomena are accompanied by
a large spatial (among sub-systems) and temporal restructuring of the system
where correlations get highly magnified. Linear response theory provides
powerful and universal tools for detecting and anticipating phase transitions
due to the presence of sufficiently strong coupling between the agents. Using
methods borrowed from spectroscopy, we have been able to predict the
appearance of poles in the susceptibility as well as to obtain a correct
estimate of critical thermodynamic properties, such as the residue of the
poles, based on the knowledge of the response for the finite particle system
in two paradigmatic models describing equilibrium and nonequilibrium phase
transitions. This is an encouraging starting point for improving our ability
to understand and predict transitions in more complex multiagent systems.

###### Acknowledgements.

VL acknowledges the support received by the European Union’s Horizon 2020
research and innovation program through the project TiPES (Grant Agreement No.
820970). The work of GP was partially funded by the EPSRC, grant number
EP/P031587/1, and by J.P. Morgan Chase & Co. Any views or opinions expressed
herein are solely those of the authors listed, and may differ from the views
and opinions expressed by J.P. Morgan Chase & Co. or its affiliates. This
material is not a product of the Research Department of J.P. Morgan Securities
LLC. This material does not constitute a solicitation or offer in any
jurisdiction. NZ has been supported by an EPSRC studentship as part of the
Centre for Doctoral Training in Mathematics of Planet Earth (grant number
EP/L016613/1).

## References

  * Naldi  _et al._ (2010) G. Naldi, L. Pareschi, and G. Toscani,  _Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences_ (Birkhäuser Basel, 2010).
  * Pareschi and Toscani (2013) L. Pareschi and G. Toscani, Interacting multiagent systems: kinetic equations and monte carlo methods (2013).
  * Dawson (1983) D. A. Dawson, Critical dynamics and fluctuations for a mean-field model of cooperative behavior, Journal of Statistical Physics 31, 29 (1983).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77, 137 (2005).
  * Garnier  _et al._ (2013) J. Garnier, G. Papanicolaou, and T. Yang, Large deviations for a mean field model of systemic risk, SIAM Journal on Financial Mathematics 4, 151 (2013), https://doi.org/10.1137/12087387X .
  * Wang  _et al._ (2017) C. Wang, Q. Li, W. E, and B. Chazelle, Noisy hegselmann-krause systems: Phase transition and the 2r-conjecture, Journal of Statistical Physics 166, 1209 (2017).
  * Garnier  _et al._ (2017) J. Garnier, G. Papanicolaou, and T. Yang, Consensus convergence with stochastic effects, Vietnam Journal of Mathematics 45, 51 (2017).
  * Garbuno-Inigo  _et al._ (2020) A. Garbuno-Inigo, N. Nüsken, and S. Reich, Affine invariant interacting Langevin dynamics for Bayesian inference, SIAM J. Appl. Dyn. Syst. 19, 1633 (2020).
  * Simmonds  _et al._ (2019) J. Simmonds, J. A. Gómez, and A. Ledezma, The role of agent-based modeling and multi-agent systems in flood-based hydrological problems: a brief review, Journal of Water and Climate Change 10.2166/wcc.2019.108 (2019), jwc2019108, https://iwaponline.com/jwcc/article-pdf/doi/10.2166/wcc.2019.108/625177/jwc2019108.pdf .
  * Geisendorf (2018) S. Geisendorf, Evolutionary climate-change modelling: A multi-agent climate-economic model, Computational Economics 52, 921 (2018).
  * Scheffer (2009) M. Scheffer,  _Critical Transitions in Nature and Society_, Princeton Studies in Complexity (Princeton University Press, Princeton, 2009).
  * Sornette (0006) D. Sornette, Endogenous versus exogenous origins of crises, in  _Extreme Events in Nature and Society_ , edited by K. H. Albeverio S., Jentsch V. (, Berlin, Heidelberg, 20006) pp. 95–119.
  * Kubo (1966) R. Kubo, The fluctuation-dissipation theorem, Reports on Progress in Physics 29, 255 (1966).
  * Lucarini  _et al._ (2005) V. Lucarini, J. J. Saarinen, K.-E. Peiponen, and E. M. Vartiainen,  _Kramers-Kronig relations in Optical Materials Research_ (Springer, New York, 2005).
  * Binney and Tremaine (2008) J. Binney and S. Tremaine,  _Galactic Dynamics_ , 2nd ed. (Princeton University Press, Princeton, 2008).
  * Öttinger (2005) H. Öttinger,  _Beyond Equilibrium Thermodynamics_ (Wiley, Hoboken, 2005).
  * Marconi  _et al._ (2008a) U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Fluctuation–dissipation: Response theory in statistical physics, Physics Reports 461, 111 (2008a).
  * Baiesi and Maes (2013) M. Baiesi and C. Maes, An update on the nonequilibrium linear response, New Journal of Physics 15, 013004 (2013).
  * Cessac (2019) B. Cessac, Linear response in neuronal networks: From neurons dynamics to collective response, Chaos 29, 103105 (2019).
  * Sarracino and Vulpiani (2019) A. Sarracino and A. Vulpiani, On the fluctuation-dissipation relation in non-equilibrium and non-hamiltonian systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 083132 (2019), https://doi.org/10.1063/1.5110262 .
  * Gottwald (2020) G. A. Gottwald, Introduction to focus issue: Linear response theory: Potentials and limits, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 020401 (2020), https://doi.org/10.1063/5.0003135 .
  * Leith (1975) C. E. Leith, Climate response and fluctuation dissipation, J. Atmos. Sci. 32, 2022 (1975).
  * North  _et al._ (1993) G. North, R. Bell, and J. Hardin, Fluctuation dissipation in a general circulation model, Clim. Dyn. 8, 259 (1993).
  * Gritsun  _et al._ (2008) A. Gritsun, G. Branstator, and A. J. Majda, Climate response of linear and quadratic functionals using the fluctuation-dissipation theorem, J. Atmos. Sci. 65 (2008).
  * Lucarini  _et al._ (2017) V. Lucarini, F. Ragone, and F. Lunkeit, Predicting climate change using response theory: Global averages and spatial patterns, Journal of Statistical Physics 166, 1036 (2017).
  * Bódai  _et al._ (2020) T. Bódai, V. Lucarini, and F. Lunkeit, Can we use linear response theory to assess geoengineering strategies?, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 023124 (2020), https://doi.org/10.1063/1.5122255 .
  * Lembo  _et al._ (2020) V. Lembo, V. Lucarini, and F. Ragone, Beyond forcing scenarios: Predicting climate change through response operators in a coupled general circulation model, Scientific Reports 10, 8668 (2020).
  * Ruelle (1998) D. Ruelle, Nonequilibrium statistical mechanics near equilibrium: computing higher-order terms, Nonlinearity 11, 5 (1998).
  * Ruelle (2009) D. Ruelle, A review of linear response theory for general differentiable dynamical systems, Nonlinearity 22, 855 (2009).
  * Dembo and Deuschel (2010) A. Dembo and J.-D. Deuschel, Markovian perturbation, response and fluctuation dissipation theorem, Ann. Inst. Henri Poincaré Probab. Stat. 46, 822 (2010).
  * Hairer and Majda (2010) M. Hairer and A. J. Majda, A simple framework to justify linear response theory, Nonlinearity 23, 909 (2010).
  * Wormell and Gottwald (2019) C. L. Wormell and G. A. Gottwald, Linear response for macroscopic observables in high-dimensional systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 113127 (2019).
  * Liverani and Gouëzel (2006) C. Liverani and S. Gouëzel, Banach spaces adapted to Anosov systems, Ergodic Theory and Dynamical Systems 26, 189 (2006).
  * Chekroun  _et al._ (2014) M. D. Chekroun, J. D. Neelin, D. Kondrashov, J. C. McWilliams, and M. Ghil, Rough parameter dependence in climate models and the role of Ruelle-Pollicott resonances, Proceedings of the National Academy of Sciences 111, 1684 (2014).
  * Lucarini (2016) V. Lucarini, Response operators for Markov processes in a finite state space: Radius of convergence and link to the response theory for Axiom A systems, Journal of Statistical Physics 162, 312 (2016).
  * Tantet  _et al._ (2018) A. Tantet, V. Lucarini, and H. A. Dijkstra, Resonances in a Chaotic Attractor Crisis of the Lorenz Flow, Journal of Statistical Physics 170, 584 (2018), arXiv:1705.08178 .
  * Pollicott (1985) M. Pollicott, On the rate of mixing of Axiom A flows, Inventiones Mathematicae 81, 413 (1985).
  * Ruelle (1986) D. Ruelle, Resonances of chaotic dynamical systems, Physical Review Letters 56, 405 (1986).
  * Shiino (1987) M. Shiino, Dynamical behavior of stochastic systems of infinitely many coupled nonlinear oscillators exhibiting phase transitions of mean-field type: H theorem on asymptotic approach to equilibrium and critical slowing down of order-parameter fluctuations, Phys. Rev. A 36, 2393 (1987).
  * Delgadino  _et al._ (2020) M. G. Delgadino, R. S. Gvalani, and G. A. Pavliotis, On the diffusive-mean field limit for weakly interacting diffusions exhibiting phase transitions (2020), arXiv:2001.03920 .
  * Carrillo  _et al._ (2020) J. A. Carrillo, R. S. Gvalani, G. A. Pavliotis, and A. Schlichting, Long-time behaviour and phase transitions for the mckean–vlasov equation on the torus, Archive for Rational Mechanics and Analysis 235, 635 (2020).
  * Lucarini  _et al._ (2020) V. Lucarini, G. A. Pavliotis, and N. Zagli, Response theory and phase transitions for the thermodynamic limit of interacting identical systems, Proc. R. Soc. A. 476 (2020), https://doi.org/10.1098/rspa.2020.0688 .
  * Jackson (1975) J. D. Jackson,  _Classical electrodynamics; 2nd ed._ (Wiley, New York, NY, 1975).
  * Talebian and Talebian (2013) E. Talebian and M. Talebian, A general review on the derivation of clausius-mossotti relation, Optik 124, 2324 (2013).
  * Sornette and Helmstetter (2003) D. Sornette and A. Helmstetter, Endogenous versus exogenous shocks in systems with memory, Physica A: Statistical Mechanics and its Applications 318, 577 (2003).
  * Yoshida (2003) N. Yoshida, Phase transition from the viewpoint of relaxation phenomena, Rev. Math. Phys. 15, 765 (2003).
  * Desai and Zwanzig (1978) R. C. Desai and R. Zwanzig, Statistical mechanics of a nonlinear stochastic model, Journal of Statistical Physics 19, 1 (1978).
  * Bonilla  _et al._ (1987) L. L. Bonilla, J. Casado, and M. Morillo, Self-synchronization of populations of nonlinear oscillators in the thermodynamic limit, Journal of Statistical Physics 48, 571 (1987).
  * Sznitman (1989) A. Sznitman,  _Topics in propagation of chaos._ , Hennequin PL. (eds) Ecole d’Eté de Probabilités de Saint-Flour XIX — 1989. Lecture Notes in Mathematics, Vol. 1464 (Springer, Berlin, Heidelberg, 1989).
  * Oelschlager (1984) K. Oelschlager, A martingale approach to the law of large numbers for weakly interacting stochastic processes, Ann. Probab. 12, 458 (1984).
  * Dawson and Gärtner (1987) D. A. Dawson and J. Gärtner, Large deviations from the mckean-vlasov limit for weakly interacting diffusions, Stochastics 20, 247 (1987), https://doi.org/10.1080/17442508708833446 .
  * Helffer (2002) B. Helffer,  _Semiclassical analysis, Witten Laplacians, and statistical mechanics_, Series in Partial Differential Equations and Applications, Vol. 1 (World Scientific Publishing Co., Inc., River Edge, NJ, 2002) pp. x+179.
  * Frank (2004) T. Frank, Fluctuation–dissipation theorems for nonlinear fokker–planck equations of the desai–zwanzig type and vlasov–fokker–planck equations, Physics Letters A 329, 475 (2004).
  * Chekroun  _et al._ (2020) M. D. Chekroun, A. Tantet, H. A. Dijkstra, and J. D. Neelin, Ruelle–pollicott resonances of stochastic systems in reduced state space. part i: Theory, Journal of Statistical Physics 10.1007/s10955-020-02535-x (2020).
  * Note (1) Repeated indices are summed.
  * Marconi  _et al._ (2008b) U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Fluctuation-dissipation: Response theory in statistical physics, Phys. Rep. 461, 111 (2008b).
  * Note (2) There is a typo in the formula given in Lucarini  _et al._ (2020). Furthermore, the convention for the Fourier transform we have used in this paper has the opposite sign, hence ditto the residue.