# Human-Robot Kinaesthetic Interaction Based on Free Energy Principle

Hiroki Sawada, Wataru Ohata, Jun Tani∗
Cognitive Neurorobotics Research Unit
Okinawa Institute of Science and Technology Graduate University
Okinawa, Japan 904-0302.
{hiroki.sawada1, wataru.ohata<EMAIL_ADDRESS>

###### Abstract

The current study investigated possible human-robot kinaesthetic interaction
using a variational recurrent neural network model, called PV-RNN, which is
based on the free energy principle. Our prior robotic studies using PV-RNN
showed that the nature of interactions between top-down expectation and
bottom-up inference is strongly affected by a parameter, called the meta-
prior, which regulates the complexity term in free energy. The current study
examines how changing the meta-prior $w$ in the interaction phase affects the
counter force generated when an experimenter attempts to induce movement
pattern transitions familiar to the robot through its prior training. The
study also compares the counter force generated when trained transitions are
induced by a human experimenter and when untrained transitions are induced.
Our experimental results indicated that (1) the human experimenter needs
more/less force to induce trained transitions when $w$ is set with
larger/smaller values, (2) the human experimenter needs more force to act on
the robot when he attempts to induce untrained as opposed to trained movement
pattern transitions. Our analysis of time development of essential variables
and values in PV-RNN during bodily interaction clarified the mechanism by
which gaps in actional intentions between the human experimenter and the robot
can be manifested as reaction forces between them.

_K_ eywords Human-robot kinaesthetic interaction $\cdot$ predictive coding
$\cdot$ active inference $\cdot$ free energy principle.

##  1 Introduction

Studies on social human-robot interaction have attracted considerable
attention recently because of their practical applications, especially with
using the linguistic modality [1, 2, 3]. However, investigation of more direct
interaction such as via kinaesthesia should be also indispensable when
considering more embodied aspects or enactivism [4] of human-robot
interactions. Although there have been a reasonable number of practical
studies on human-robot interaction using kinaesthesia, which have contributed
greatly to human-robot joint collaboration, human assistance, and user
interface [5, 6, 7], few studies have attempted to understand essential
mechanisms underlying kinaesthetic interaction in light of embodied cognition,
social cognition, and system-level neuroscience.

In this regard, the current study investigated human-robot kinaesthetic
interaction by applying system neuroscience theory, the free energy principle,
proposed by Friston [8] which is consonant with enactivism [9, 10]. Let us
consider a situation in which a robot and a human dance, holding each other
with both hands, executing memorized dance patterns. If the robot initiates a
particular pattern from memory with strong intention, the human counterpart
might follow it without resisting because of the strong counter force. On the
other hand, if the robot generates a pattern without strong intention, the
human counterpart might be able to shift to a different pattern without
experiencing strong counter force. Let us consider another situation. If the
human counterpart attempts to induce a movement pattern that is familiar to
the robot, such guidance should proceed easily without strong counter force,
since the robot can infer the intended pattern immediately and can move as
anticipated. On the other hand, if the human counterpart attempts to induce a
movement pattern unfamiliar to the robot, such guidance should experience
strong counter force, since the movement is neither inferential nor
predictable for the robot.

We presume that these situations can be well explained by predictive coding
[11] and active inference [12, 13, 14], based on the free energy principle
[8]. Predictive coding provides a theory in which perception is achieved by
inferring hidden causes for sensory observations that minimize the error
between top-down prediction of sensation made by the generative model and the
actual sensation. On the other hand, active inference provides a theory for
action generation in which an optimal action is inferred to minimize the error
between the preferred sensation and the actual sensation resulting from the
action. These two are not separate, but integrated through a sensorimotor loop
in embodied cognition in which a perceptual inference and action generation
can be achieved simultaneously by minimising errors through iterative
interaction between the top-down predictive/generative process and the bottom-
up inference. Also in social embodied cognition, the dynamic interaction
between the top-down pathway for predicting while acting on the others and the
bottom-up pathway for inferring the actional intention of the counterpart
through sensory observation should become a crucial element. Accordingly, we
presume that different interactions in the dance appear, depending on the
balance between the strength of the top-down pathway and that of the bottom-up
pathway.

Recently, the number of studies addressing application of free energy
principle in cognitive robotics has increased [15, 16, 17]. Maselli et al.
[18] showed that the active inference model is able to characterise movements
generated by the agent’s intention to resolve multi-sensory conflict or to
achieve an external goal, such as reaching with its arm to a certain point
with the agent having a VR-vision of its arm that was tilted from the actual
position to confuse the agent. Tschantz et al. [19] presented a working
implementation of active inference to reinforcement learning that demonstrated
efficient exploration and an order of magnitude higher sample efficiency in a
high-dimensional task, such as a mountain-car environment. Also, Pezzato et
al. [20], showed that a robotic task, such as reactive action planning can be
formulated as a free energy minimisation problem by introducing a hybrid
combination of active inference and behaviour trees. However, studies
investigating human-robot interaction using the free energy principle remain
few in number.

The author’s group has conducted neurorobotic studies related to frameworks of
predictive coding and active inference to develop various types of recurrent
neural network (RNN) models [21]. However, an essential problem in these
studies is that RNN models have difficulty in dealing with probabilistic
properties hidden in interactions between robots and environments. To tackle
this problem, our group developed a probabilistic variational RNN model,
called PV-RNN [22], based on the free energy principle.

An indispensable feature of PV-RNN is that free energy is computed as a sum of
the negative accuracy term and the complexity term, weighted by a parameter
$w$, called the meta-prior. Here the complexity term represents the divergence
between the approximated posterior probability distribution and the prior
probability distribution in probabilistic latent variables allocated in PV-
RNN. More intuitively, the complexity term can be understood as the internal
gap between top-down expectation and bottom-up sensory reality. The model can
learn to extract probabilistic structures hidden in data either by embedding
them in nonlinear, deterministic dynamics of chaos by setting a large value of
the meta-prior $w$ or into stochastic processes by setting a small $w$ [22].

By following the above study, Chame et al. [23] conducted a human-robot bodily
interaction experiment using a PV-RNN model and attempted to show possible
effects of $w$ on the interaction dynamics. However, this study is considered
preliminary because only some snapshots of the experiments were shown and no
rigorous analysis of the experimental results was made. Wirkuttis et al. [24]
performed simulation studies on synchronized imitative interaction of dyadic
robots using a PV-RNN model in which robots observed the sensation of
simplified vision and proprioception, but without kinaesthesis. Through
statistical analysis of repeated simulation experiments, these studies show
that a robot with a larger $w$ in both the learning phase and the test phase
tends to lead the counterpart robot set with smaller $w$ by projecting
stronger actional intention to the counterpart.

The current study examined human-robot kinaesthetic interaction by
implementing a PV-RNN model in a real humanoid robot. In particular, the study
attempted to translate top-down intention or prior belief, formulated using
the free energy principle to the force bodily exerted on the counterpart by
conducting statistical analysis on data obtained from repeated experiments. In
this experiment, a humanoid robot controlled by a PV-RNN model is trained to
generate movement patterns with specific transition probability distributions
among them. After the training phase, two types of human-robot kinaesthetic
interaction experiments were conducted. In the first experiment, while the
robot was generating trained patterns, a human experimenter tried to
physically induce a set of trained familiar movement pattern transition in the
robot. We examined how the robot’s rigidity in terms of the counter force in
its response changed depending on $w$, which was set during the interaction
phase. The second experiment compared the counter force during responses in
these two cases, i.e., when the experimenter induced the robot to proceed with
learned movement pattern transitions and when the experimenter initiated
previously unlearned pattern transitions.

Our intensive analysis of time-development of the latent variables, complexity
term, and prediction error observed in these experiments under different
conditions clarifies how bodily interaction between an experimenter and a
robot proceeds and what sorts of neural activities are generated through top-
down and bottom-up interactions during an experimental task.

##  2 Model

###  2.1 Free Energy Principle

As mentioned in the Introduction, the current study used PV-RNN, based on the
free energy principle, as the basic model. The free energy principle assumes
that brains learn as well as inferred for given observations, by minimizing
free energy, defined in Eq. 1.

$\displaystyle\mathcal{F}=\underbrace{D_{\rm
KL}[q_{\phi}({\boldsymbol{z}}|{\boldsymbol{X}})\|p_{\theta}(\boldsymbol{z})]}_{\rm
complexity}-\underbrace{\mathbb{E}_{q_{\phi}({\boldsymbol{z}}|{\boldsymbol{X}})}[\log
p_{\theta}({\boldsymbol{X}}|{\boldsymbol{z}})]}_{\rm accuracy},$ (1)

$p_{\theta}({\boldsymbol{X}}|{\boldsymbol{z}})$ is the likelihood of the
sensory observation $\boldsymbol{X}$, given the probabilistic latent variables
$z$ which is parameterized by $\theta$.
$q_{\phi}({\boldsymbol{z}}|{\boldsymbol{X}})$ is the inference model
parameterised by $\phi$. As shown in Eq. 1, free energy consists of two terms.
The first term indicates the $complexity$, which is the divergence between the
approximate posterior probability distribution and prior probability
distribution for the latent variable $z$, and the second term is the
$accuracy$ in predicting the sensory observation.

###  2.2 PV-RNN implementation

PV-RNN is operated in two distinct phases, a training phase and an interaction
phase. In the training phase, we prepared a dataset (joint angle sequences of
the robot) and trained the PV-RNN model using it. In the interaction phase,
the robot and the human experimenter interact physically, such that the PV-RNN
attempts to drive the robot’s arms by predicting next-time-step target joint
angles while it infers latent variables using actual joint angle readings.
Each operation phase is explained in detail below.

####  2.2.1 Training Phase

First, we start with computation during the training phase. The evidence-free
energy ${\mathcal{F}}^{train}$ as a loss function for training the PV-RNN with
the joint angle observation from time step $1$ to $T$ is shown in Eq.2 (the
exact derivation should refer to [25]).

$\begin{split}{\mathcal{F}}^{train}&=\beta\times\underbrace{D_{\text{KL}}[q_{\phi}(\mathbf{z}_{1}|\mathbf{d}_{0},\mathbf{X}_{1})\|p_{\theta}(\mathbf{z}_{1})]}_{\rm
complexity}\\\
+&\underbrace{w\times\sum_{t=1}^{T-1}\mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{d}_{t},\mathbf{X}_{t:T+1})}\big{[}D_{\text{KL}}[q_{\phi}(\mathbf{z}_{t+1}|\mathbf{d}_{t},\mathbf{X}_{t+1:T+1})\|p_{\theta}(\mathbf{z}_{t+1}|\mathbf{d}_{t})]\big{]}}_{\rm
complexity}\\\
-&\underbrace{\sum_{t=1}^{T}\mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t-1}|\mathbf{d}_{t-1},\mathbf{X}_{t:T})}[\ln
p_{\theta}(\mathbf{X}_{t}|\mathbf{d}_{t})]}_{\rm accuracy},\end{split}$ (2)

where $q(\phi)$ and $p(\theta)$ are the inference model parameterised by
$\phi$ and the generative model parameterised by $\theta$, respectively.

As shown in Eq.2, the meta-prior $w$ weights the complexity term, except for
the initial step. As described previously, the setting of meta-prior in the
training phase strongly affects behavioural characteristics of the trained
network. Prior studies [25, 24] showed that a network trained with a larger
$w$ develops higher precision, i.e., smaller standard deviation, in the prior
distribution, whereas a smaller $w$ develops lower precision in the prior.
More or less precision in the prior means stronger or weaker top-down belief,
respecyively. PV-RNN is composed of two variables $d,z$; in which the former
is a deterministic latent variable, and the latter is a random latent variable
sampled from a Gaussian distribution. During forward computation of PV-RNN,
$d$ in layer $l$ ($l=1$ is the bottom layer, which is closest to the output
layer) at time step $t$ is computed as,

$\displaystyle\mathbf{h}^{l}_{t}=$
$\displaystyle\left(1-\frac{1}{\tau^{l}}\right)\mathbf{h}^{l}_{t-1}+\frac{1}{\tau^{l}}\biggl{(}\mathbf{W}^{ll}_{dd}\mathbf{d}^{l}_{t-1}+\mathbf{W}^{ll}_{zd}\mathbf{z}^{l}_{t}+\mathbf{W}^{ll+1}_{dd}\mathbf{d}^{l+1}_{t}+\mathbf{b}_{h}\biggr{)},$
(3) $\displaystyle\mathbf{d}^{l}_{t}=$
$\displaystyle\tanh(\mathbf{h}^{l}_{t}).$

Here, $\mathbf{h}_{t}^{l}$ denotes the internal state of $\mathbf{d}_{t}^{l}$
before the activation function $\tanh$ is applied, and $\mathbf{b}_{h}$
denotes the bias term of $\mathbf{h}$. $\tau^{l}$ is a time constant unique to
each layer, which encourages the network to process information by following
an intrinsic time scale of the layer [26, 27, 28, 29]. Matrices
$\mathbf{W}_{dd}$ and $\mathbf{W}_{zd}$ express intra- and inter-layer
connectivity weights in the network, respectively. At $t=1$, the input to
$\mathbf{h}$ in the top layer is calculated only from $\mathbf{z}_{1}^{l}$ and
$\mathbf{b}_{h}$.

Each dimension of $\mathbf{z}_{t}^{p}$ is sampled over the prior distribution,
which is a Gaussian distribution of mean $\boldsymbol{\mu}_{t}^{p}$ and
standard deviation $\boldsymbol{\sigma}_{t}^{p}$ individually. At $t=1$,
$\mathbf{z}_{1}^{p}$ is sampled from a standard normal distribution, and
$z_{t}^{p}$ for the following time steps is computed by,

$\displaystyle\boldsymbol{\mu}^{p}_{t}$
$\displaystyle=\tanh(\mathbf{W}^{ll}_{d\mu}\mathbf{d}_{t-1}+\mathbf{b}_{\mu}^{p}),$
(4) $\displaystyle\boldsymbol{\sigma}^{p}_{t}$
$\displaystyle=\exp(\mathbf{W}^{ll}_{d\sigma}\mathbf{d}_{t-1}+\mathbf{b}_{\sigma}^{p}),$
$\mathbf{z}^{p}_{t}={\boldsymbol{\mu}^{p}_{t}}+{\boldsymbol{\sigma}^{p}_{t}}*\boldsymbol{\epsilon}_{t},\quad\textrm{with
}\boldsymbol{\epsilon}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{I}),$ (5)

where $*$ indicates the element-wise product of the two vectors. This equation
follows the idea of the conditional prior [30]. Here,
$\mathbf{W}^{ll}_{d\mu},\mathbf{W}^{ll}_{d\sigma},\mathbf{b}_{\mu}^{p}$ and
$\mathbf{b}_{\sigma}^{p}$ are the weight matrices and the bias terms for
$\boldsymbol{\mu}_{t}^{p}$ and $\boldsymbol{\sigma}_{t}^{p}$, respectively.
$\boldsymbol{\epsilon}$ is a value sampled from the standard normal
distribution. In order to make model parameters differentiable through the
random latent variable, we use the reparametrization trick [31] to sample
$\mathbf{z}_{t}^{p}$.

On the other hand, $\mathbf{z}_{t}^{q}$ in inference model $q_{\phi}$ is
sampled from the approximate posterior probability distribution
$q(\mathbf{z}_{t})$ which is a Gaussian distribution with mean $\mu_{t}^{q}$
and standard deviation $\sigma_{t}^{q}$ which are computed as,

$\displaystyle\boldsymbol{\mu}^{q}_{t}$
$\displaystyle=\tanh(\mathbf{A}^{\mu}_{t}),$ (6)
$\displaystyle\boldsymbol{\sigma}^{q}_{t}$
$\displaystyle=\exp(\mathbf{A}^{\sigma}_{t}),$
$\displaystyle\mathbf{z}^{q}_{t}$
$\displaystyle=\boldsymbol{\mu}^{q}_{t}+\boldsymbol{\sigma}^{q}_{t}*\boldsymbol{\epsilon}_{t},\quad\textrm{with
}\boldsymbol{\epsilon}_{t}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$

$\mathbf{A}^{\mu}_{t},\mathbf{A}^{\sigma}_{t}$ are adaptive variables which
are optimised during training. During the training phase, $\mathbf{z}_{t}$ in
Eq.3 is sampled from the approximate posterior. On the other hand,
$\mathbf{z}_{t}^{p}$ is sampled from the prior distribution to compute the
Kullback-Leibler divergence (KLD) between the approximate posterior and the
prior, which is the complexity term of the loss function Eq.2.

The obtained $\mathbf{d}_{t}^{1}$ is used for computing the output layer, of
which function is,

$\mathbf{o}_{t}=\mathbf{W}_{o}\mathbf{d}_{t}^{1}+\mathbf{b}_{0},$ (7)

where $\mathbf{W}_{o},\mathbf{b}_{o}$ is the weight matrix and the bias of the
output layer, respectively.

Finally, we compute the predicted output $\mathbf{\bar{x}}_{t}$. Each
dimension of the predicted output is represented in a probabilistic
distribution using the Softmax of $N_{soft}$ elements. Therefore, the $j$-th
softmax element of the $i$-th dimension of the predicted output
$\bar{x}_{t}^{i,j}$ is computed by,

$\bar{x}_{t,i,j}=\frac{\exp(o_{t,i,j})}{\sum_{j=1}^{N_{soft}}\exp(o_{t,i,j})}.$
(8)

As shown in Eq.2, the evidence-free energy $\mathcal{F}$, which is the loss
function of this model, is the sum of the KLD between the approximate
posterior $q(\mathbf{z}_{t}|\mathbf{x}_{t})$ and prior
$p(\mathbf{z}_{t}|\mathbf{d}_{t-1})$, and the negative log-likelihood
calculated from the output $\mathbf{\bar{x}}_{t}$ and the perception
$\mathbf{x}_{t}$.

At time step $t$, the KLD between the approximate posterior and prior is
defined as,

$D_{KL}[q(\mathbf{z}_{t}|\mathbf{x}_{t:T})||p(\mathbf{z}_{t}|\mathbf{d}_{t-1})]=\sum_{i}^{N_{z}}q(z_{t}^{i}|\mathbf{x}_{t:T})\ln\frac{p(z_{t}^{i}|\mathbf{d}_{t-1})}{q(z_{t}^{i}|\mathbf{x}_{t:T})},$
(9)

where $z_{t}^{i}$ is the $i$-th dimension of $\mathbf{z}_{t}$, $T$ is the
length of the sequence and $N_{z}$ is the dimension of $\mathbf{z}_{t}$ of
each layer. Here, $q(\mathbf{z}_{t}|\mathbf{x}_{t})$ and
$p(\mathbf{z}_{t}|\mathbf{d}_{t-1})$ are both assumed to follow multivariate
Gaussian distributions with diagonal covariant matrices. Hence, they are
expressed with the mean $\boldsymbol{\mu}_{t}$
($=[\mu_{t,1},\mu_{t,2},...,\mu_{t,N_{z}}]$) and the standard deviation
$\boldsymbol{\sigma}_{t}$
($=[\sigma_{t,1},\sigma_{t,2},...,\sigma_{t,N_{z}]}$) as below,

$q(z_{t}^{i}|\mathbf{x}_{t:T})=\frac{1}{\sqrt{2\pi(\sigma_{t,i}^{q})^{2}}}\exp{\left[-\frac{1}{2}\left(\frac{z_{t}^{i}-\mu^{q}_{t,i}}{\sigma_{t,i}^{q}}\right)^{2}\right]},$
(10)
$p(z_{t}^{i}|\mathbf{d}_{t-1})=\frac{1}{\sqrt{2\pi(\sigma_{t,i}^{p})^{2}}}\exp{\left[-\frac{1}{2}\left(\frac{z_{t}^{i}-\mu_{t,i}^{p}}{\sigma_{t,i}^{p}}\right)^{2}\right]}.$
(11)

Therefore, the KLD loss $r_{t}$ between the approximate posterior and the
prior is computed as,

$r_{t}=\begin{cases}{\sum_{i}D_{KL}[q(z_{t}^{i}|\mathbf{x}_{t})||p(z_{t}^{i}|\mathcal{N}(0,1))]},\quad\textrm{if
}t=1,\\\
{\sum_{i}D_{KL}[q(z_{t}^{i}|\mathbf{x}_{t})||p(z_{t}^{i}|\mathbf{d}_{t-1})]},\quad\textrm{otherwise}.\end{cases}$
(12)

On the other hand, as an equivalent of the negative log-likelihood term of
Eq.2, the KLD between the target and the predicted output is computed. As
mentioned above, the predicted output is represented as a probabilistic
distribution $P(x_{t}^{i}|\mathbf{\theta})$ using Softmax. The target sequence
is also transformed into a probabilistic distribution $Q(x_{t}^{i})$ along the
same line [32]. Therefore, the KLD is computed as,

$e_{t}=\sum_{i}^{N_{x}}D_{KL}[Q(\mathbf{x}_{t,i})||P(\mathbf{x}_{t,i})]=\sum_{i}^{N_{x}}\sum_{j}^{N_{soft}}Q(x_{t,i,j})\ln\frac{Q(x_{t,i,j})}{P(x_{t,i}|\mathbf{\theta})}.$
(13)

The loss function of our model which is equivalent to the evidence-free energy
$\mathcal{F}^{train}$ is computed as the weighted sum of Eq.12 and Eq.13. As
shown in Eq.2, the KLD between the approximate posterior and prior is
regularised with a meta-prior $\beta$ for $t=1$ and $w$ for $t\neq 1$. Here,
$\beta$ is the weighting parameter specifically for $t=1$, which regulates the
initial sensitivity of the PV-RNN model. Therefore, the loss function of our
model during the training phase is computed as,

$\mathcal{F}=\begin{cases}\sum_{t=1}^{T}\left(e_{t}+\frac{N_{x}}{N_{z}}\beta
r_{t}\right),\quad\textrm{if }t=1,\\\
\sum_{t=1}^{T}\left(e_{t}+\frac{N_{x}}{N_{z}}w^{t}r_{t}\right),\quad\textrm{otherwise}\end{cases}$
(14)

where $N_{x}$ and $N_{z}$ are the $\bar{\mathbf{x}}$ dimension and
$\mathbf{z}$ dimension in each layer, respectively. $T$ is the length of the
dataset sequence. The negative log-likelihood, which we denote as _prediction
error_ , $e_{t}$ is proportional to $N_{x}$ and the KLD between the
approximate posterior and prior is proportional to $N_{z}$, under the
assumption that each $\bar{\mathbf{x}}$ dimension and $\mathbf{z}$ dimension
is independent. Therefore, the meta-prior during the training phase $w^{t}$ is
normalised by $N_{x}$ and $N_{z}$ so that it can be compared among PV-RNN
models with different $\bar{\mathbf{x}}$ dimension and $\mathbf{z}$ dimension.
Through training, weight matrices, biases and adaptive variables
$\mathbf{A}_{t}^{\mu},\mathbf{A}_{t}^{\sigma}$ are updated based on back-
propagation through time (BPTT) [33, 34, 35].

####  2.2.2 Interaction phase

We performed a human-robot real-time interaction in the interaction phase
where the pre-trained PV-RNN drives the robot’s movement by predicting joint
angles of the next time step. The forward computation part remains the same as
that in the training phase. However, there are some differences in the loss
function used and in ways of updating variables. Unlike the training phase,
weight matrices and biases are not updated. Only the adaptive variables
$\mathbf{A}_{t}^{\mu}$ and $\mathbf{A}_{t}^{\sigma}$ are updated so that the
approximate posterior can adapt to the ongoing observation of the joint angle
sequence.

PV-RNN predicts the future sensation and infers past latent variables using
the past window spanning from time-step $t_{c}-t_{w}$ to the current time-step
$t_{c}$ with a window size $t_{w}$ as shown in Fig. 1.

Figure 1: Graphical representation of the interaction phase of PV-RNN.
$d_{l,t},z^{p}_{l,t},z^{q}_{l,t},x_{t}$ and $\bar{x}_{t}$ indicates the
deterministic latent variable, prior distribution, approximate posterior,
output and target of layer $l$ and time-step $t$, respectively. Blue and red
arrows indicate forward propagation and backward propagation, respectively.
The future, the past within the window and the past outside the window are
coloured red, blue, and green respectively. The past window size is 2 in this
graphical model. a) shows the network at time-step $t$. b) shows the network
at time-step $t+1$.

The evidence-free energy in the interaction phase $\mathcal{F}^{int}$,
normalized in the same way as Eq.14, is computed inside the past window as:

$\mathcal{F}^{int}=\begin{cases}\sum_{t=t_{c}-t_{w}}^{t_{c}}\left(e_{t}+\frac{N_{x}}{N_{z}}\beta
r_{t}\right),\quad\textrm{if }t=1,\\\
\sum_{t=t_{c}-t_{w}}^{t_{c}}\left(e_{t}+\frac{N_{x}}{N_{z}}w^{i}r_{t}\right),\quad\textrm{otherwise},\end{cases}$
(15)

The adaptive variables $\mathbf{A}_{t}^{\mu}$ and $\mathbf{A}_{t}^{\sigma}$ at
each time-step $t$ from $t_{c}-t_{w}$ to $t_{c}$ are modified so as to
minimize $\mathcal{F}^{int}$ by iterating the forward computation and the
error back-propagation through time for a fixed number of times, called
_epochs_. After modifying the approximate posterior by going through _epochs_
of iteration, next-step joint angles are predicted and fed into the PID
controller of the robot in order to generate its movement. Then, the past
window is shifted one step ahead. $w^{i}$, which is the meta-prior used in the
interaction phase, could be set with different values from $w^{t}$ used in the
training phase. In [36], a simulation experiment using a PV-RNN model shows
that when a PV-RNN trained with a particular $w^{t}$ was reset with a smaller
$w^{i}$ in the later interaction phase, the approximate posterior shifted away
from the prior and the PV-RNN tended to adapt to the observed sensory
sequence. This means that top-down intention in the PV-RNN became weaker. On
the other hand, when reset with a larger $w^{i}$, the PV-RNN tended to ignore
the observed sensory sequence by generating its own intended patterns, which
means that top-down intention becomes stronger. However, the precision
structure in the prior does not change even as $w^{i}$ is changed, since the
prior distribution was developed in the learning phase.

###  2.3 The employed robot controller

The current study uses a humanoid robot, Torobo 111Humanoid-robot Torobo from
Tokyo Robotics: https://robotics.tokyo/products/torobo/ to conduct human-robot
kinesthetic interaction experiments. Torobo is equipped with a built-in force-
feedback controller that enables humans to back-drive joint angles with subtle
force. In the Torobo control system, when a human exerts certain torque on
joints by pushing or pulling the limbs of Torobo, the exerted torque can be
estimated by subtracting the torque inferred as necessary to account for the
current static state as well as a dynamic state of the robot from the actual
torque measured in the joints. By computing the next time-step joint target
positions by adding the current positions with the estimated exerted torque
multiplied by a constant gain and feeding them in the PID controller, Torobo’s
limbs move by following the force exerted on them by the human. This force-
feedback controller is integrated with the PV-RNN, which generates the next
time-step target joint angles (Fig.2).

Figure 2: the human interacting with Torobo, the PV-RNN target joint angle
generator, the inverse model, and the PID joint controller.

The overall control diagram is shown in Fig.2, which includes a human
experimenter interacting with Torobo, the PV-RNN target joint angle generator,
the inverse model, and the PID joint controller. The inverse model and PID
controller were developed by the manufacturer of Torobo. Details are as
follows. Excess torque, $\tau_{t}^{exc}$ which is torque exerted on joints by
the human, can be estimated by extracting $\tau_{t}^{dynm}$ as the torque
inferred for the current position, velocity, and acceleration of the joints by
using the inverse model of Torobo from the current measured torque,
$\tau_{t}$. Next, excess torque, $\tau_{t}^{exc}$ is applied with a threshold
control and $e_{t}$ as excess torque after threshold control is obtained.
Then, $e_{t}$ is time-filtered with the decay parameter $0<\alpha<1$, which
generates $\tilde{e_{t}}$ as the time-filtered excess torque. These operations
are necessary to prevent unnecessary overreaction of the torque estimation
against noise. These lines of processes are described in the following
mathematical form:

$\displaystyle\tau_{t}^{exc}$ $\displaystyle=\tau_{t}-\tau_{t}^{dynm},$ (16)
$\displaystyle e_{t}$
$\displaystyle=\begin{cases}\max(0,\tau_{t}^{exc}-\tau^{th}),\quad\textrm{if
}\tau_{t}^{ext}>0,\\\
\min(0,\tau_{t}^{exc}+\tau^{th}),\quad\textrm{otherwise},\end{cases}$
$\displaystyle\tilde{e_{t}}$
$\displaystyle=\begin{cases}e_{t},\quad\textrm{for }t=1,\\\
\max\Bigl{(}\alpha\times\tilde{e}_{t-1}+e_{t},e_{max}\Bigl{)},\quad\textrm{otherwise},\end{cases}$

where $\tilde{e_{t}}$ has a fixed upper limit $e_{max}$.

The PV-RNN predicts $\bar{\theta_{t}}$ next time-step joint angles while
performing the online inference, which is subtracted by $\theta_{t}$ as the
current joint angle to generate $\Delta\theta_{t}$, the target joint angle
difference shown as:

$\Delta\theta_{t}=\bar{\theta_{t}}-\theta_{t}$ (17)

Finally, $\Delta\theta_{t}$ representing the next move intended by the PV-RNN
and $\tilde{e}_{t}$ representing the estimate of torque exerted by the human
are multiplied by each gain $k_{r}$ and $k_{p}$, respectively and they are
added to the current joint angles to generate the final target joint
positions, $\bar{\theta_{t}}$, which are fed into the PID controller as shown
in the following:

$\accentset{\ast}{\mathbf{\theta}}_{t+1}=\Delta\theta_{t}\times
k_{r}+\tilde{e}_{t}\times k_{p}+\theta_{t},$ (18)

##  3 Experiment

###  3.1 Experiment Setup

####  3.1.1 Train Data-Set Preparation

The proposed PV-RNN model was trained in a supervised manner by preparing
4-dimensional joint angle teaching trajectories. In preparation for teaching
trajectories, we considered four types of cyclic movement patterns (20 time-
steps for each cycle) using Torobo’s shoulder and elbow joint angles in both
arms. Then, it was assumed that cycling movement patterns transit from one to
another following a probabilistic finite state machine (Fig.3). For example,
after a movement pattern _A_ is generated for one cycle with the state at _S1_
, _A_ can be generated for one more cycle with $90\%$ probability staying at
the same state, or _B_ or _C_ can transit to _S2_ or _S3_ with a probability
of $3\%$ and $7\%$, respectively.

Figure 3: Schematic of the probabilistic finite state machine from which
training data was generated.

We prepared 10 sequences that each consisted of 200 cycles of movement
patterns, extending 4000 time-steps.

####  3.1.2 The Network Configuration and Training

To conduct human-robot dyadic interaction experiments, PV-RNN was trained 3
times with identical parameters (Table.1).
$\\#\mathbf{d},\\#\mathbf{z},\tau,\mathbf{w}^{t},\mathbf{w}^{i}$ indicates the
number of d neurons and z neurons, time constant, and meta-prior during the
training phase and interaction phase, respectively. The network was trained
for 50,000 epochs to minimize the evidence-free energy shown in Eq.2, starting
with random weights generated with different seeds for each training.

Table 1: PV-RNN Parameters | $\\#\mathbf{d}$ | $\\#\mathbf{z}$ | $\tau$ | $\mathbf{w}^{t}$ | $\mathbf{w}^{i}$
---|---|---|---|---|---
Layer 1 (Top) | 60 | 6 | $3$ | 0.01 | $[0.01,0.05,0.1]$
Layer 2 (Bottom) | 30 | 3 | $9$ | 0.01 | $[0.01,0.05,0.1]$

Fig.4 shows one of the resultant training processes.

Figure 4:  Resultant time-development from the training phase of one of the
PV-RNN models. Development of KLD with respect to the number of training
epochs is shown for (a) the Top layer (layer 2) and (b) the Bottom layer
(layer 1). Development of the average prediction error over time steps in each
teaching sequence is shown in (c).

It can be seen that prediction error and the KL divergence in both the top and
the bottom layers decreased throughout training. All three training processes
converged in a similar way, achieving prediction errors and KL-divergences
shown in Table.2.

Table 2: Training processes of all three PV-RNN model |  | Epoch
---|---|---
| Model | 2,000 | 5,000 | 50,000
KL divergence | 1 | $2.0\times 10^{-3}$ | $5.4\times 10^{-4}$ | $4.4\times 10^{-5}$
2 | $2.4\times 10^{-3}$ | $9.1\times 10^{-4}$ | $5.6\times 10^{-5}$
3 | $2.5\times 10^{-3}$ | $7.2\times 10^{-4}$ | $6.9\times 10^{-5}$
Pred. Err. | 1 | $1.8\times 10^{-2}$ | $7.2\times 10^{-4}$ | $1.4\times 10^{-4}$
2 | $1.4\times 10^{-2}$ | $6.6\times 10^{-3}$ | $4.6\times 10^{-5}$
3 | $1.3\times 10^{-2}$ | $1.6\times 10^{-3}$ | $1.4\times 10^{-4}$

####  3.1.3 Evaluation of the trained networks

Trained PV-RNNs were evaluated on the basis of how closely the state
transition probability wasreconstructed compared with the target (Fig.3). For
this purpose, we conducted _Prior Generation_ in which the forward computation
by following Eq.3-7 was performed without any external observation for 40,000
time-steps for each trained network. Then, probabilities for all possible
movement pattern transitions were measured during prior generation and the
resulting state transition probability was inferred (Table.3). The state
transition probability computed for all trained networks is quite similar to
the target one. Output and network dynamics of movement pattern transitions
during prior generation are shown in Appendix, Fig.8.

Table 3: Transition probabilities of the _Prior Generation_ and data set Model |  _S1_ $\rightarrow$ _S2_ |  _S1_ $\rightarrow$ _S3_ |  _S2_ $\rightarrow$ _S4_ |  _S3_ $\rightarrow$ _S4_ |  _S4_ $\rightarrow$ _S1_
---|---|---|---|---|---
1 | 2.8 % | 3.5 % | 15.6 % | 12.5 % | 3.1 %
2 | 2.7 % | 9.3 % | 14.1 % | 11.5 % | 7.2 %
3 | 3.1 % | 4.1 % | 9.8 % | 12.3 % | 5.1 %
Target | 3 % | 7 % | 10 % | 15 % | 5 %

####  3.1.4 Human-Robot Interaction

After the training phase, we conducted two types of human-robot interaction
experiments using trained PV-RNNs. In these experiments, while Torobo was
generating movement pattern transitions successively based on the training,
the human experimenter attempted to induce various movement pattern
transitions by grasping both arms of Torobo and exerting force on them. These
movement pattern transitions included trained transitions (_AB, AC, BD, CD,
DA_) and untrained transitions (_AD, DB, DC, BA, CA_) in which _AB_ , for
example, dictates that _A_ pattern is forced to transit to _B_ pattern. Each
transition from one pattern to another requires some guiding force by a human
experimenter, even for trained transitions, since ongoing patterns tend to
repeat another cycle with high probability, more than $85\%$ for all patterns.
This means that there exist some conflicts between movement trajectories
intended by the robot and those by the human experimenter during both trained
and untrained transitions.

Experiment-1 examined the effect of $w^{i}$ settings on the interactions.
While Torobo was generating movement pattern transitions successively for
2,200 time-steps, the human experimenter attempted to induce one of the
trained transition every 200 time-steps starting from $t=200$ which resulted
in 10 trained transitions. The duration of each attempt lasted 100 time-steps
at the most. This experiment was conducted three times for all three trained
PV-RNNs by changing $w^{i}$ with $0.01$, $0.05$, and $0.1$.

Experiment-2 examined the difference between trained and untrained transitions
by setting $w^{i}$ with a fixed value. The experimental procedure was the same
as in Experiment-1, although this time, $w^{i}$ was fixed at $0.01$ for both
layers. This experiment, using all three trained PV-RNNs, resulted in 30
untrained transition attempts. In these experiments, we recorded the time
development of essential values including latent variables, predicted and
observed joint angles, prediction error, and the KL-divergence in both layers
for later analysis of experiment results.

###  3.2 Experiment Results

In this section, we show the results of the aforementioned experiments.

####  3.2.1 Experiment-1

First, we examined time-development of essential values during movement
pattern transitions induced by the experimenter for each case with a different
meta-prior setting. Fig.5 shows an example snapshot of future prediction and
past reflection, which shifted every 7 time-steps of the current time during
the transition _AB_ performed under different settings of meta-prior,
$w^{i}=0.01,0.05,0.1$.

Figure 5: Time-development of excess torque, prediction error, and KL-
divergence between the approximate posterior and the prior in trained movement
transitions in cases with three different meta-prior $w^{i}$ settings with
$0.01$, $0.05$, and $0.1$. The grey area represents the past window where the
head of the window is the current time. The window is shifted 5 times during
the movement transition _AB_.

Each snapshot shows one of the observed joint angles $\theta$ (dotted blue)
and its prediction $\bar{\theta}$ (blue) in the top row, the KL-divergences
between the approximate posterior and the prior in the layer 1 (orange) and in
the layer 2(dark orange) in the second row, the prediction error (green) in
the third row, and the excess torque (black) in the bottom row. The grey area
represents the past window where the approximate posterior in terms of
adaptive variables $\mathbf{A}_{t}^{\mu},\mathbf{A}_{t}^{\sigma}$ is updated.
We provide two supplementary videos for experiment-1 showing the interaction
between the experimenter and Torobo, as well as network dynamics in the case
with the meta-prior $w^{i}$ set to 0.01 (video-link1) and 0.1 (video-link2).

The sequences of time-shifted snapshots in Fig.5, show that excess torque
appears first followed by rises in the prediction error (negative log-
likelihood) and the KL-divergence. Later, the predicted joint angle pattern
shifts from pattern _A_ to _B_ while the gap between the observed joint angle
and the reconstructed joint angle remains in the past window. This is the same
for all three cases with different $w^{i}$ settings.

However, we can see some qualitative differences in the transition process
depending on the $w^{i}$ setting. The prediction error and the excess torque
in the case with a small $w^{i}$ ($w^{i}=0.01$) are smaller than those in the
case with large $w^{i}$ ($w^{i}=0.1$). Also, the error and torque with a small
$w^{i}$ are less persistent than those with a larger $w^{i}$. However, the KL-
divergence in the case with a small $w^{i}$ is larger and persists longer than
that in the case with large $w^{i}$. In order to confirm these observations,
we conducted statistical analysis on the excess torque, KL-divergence, and
prediction error time-averaged for each transition period (100 steps). This
computation was repeated 10 times for each of three different trained networks
set with three different $w^{i}$ values.

The results are shown in Fig.6.

Figure 6: Time-averaged excess torque, prediction error, and the KL-divergence
between the approximate posterior and the prior in the trained movement
transition in cases with three different meta-prior $w^{i}$ settings with
$0.01$, $0.05$, and $0.1$.

Both the time-averaged prediction error and the excess torque measured in the
cases with small $w^{i}$ are significantly smaller than those with large
$w^{i}$. On the other hand, the time-averaged KL-divergence with small $w^{i}$
is significantly larger than that with large $w^{i}$.

By considering these statistical results, movement pattern transitions exerted
by the experimenter require greater force when $w^{i}$ is set larger, since
the approximate posterior distribution strongly follows the prior distribution
representing the current movement intention of PV-RNN, minimizing the KL-
divergence between them while the error $\Delta\theta_{t}$ between the
predicted joint angle $\bar{\theta_{t}}$ and the observed joint angle
$\theta_{t}$ becomes larger. If the experimenter attempts to move the
trajectory of the robot’s joint angles in a direction different from that
predicted by the PV-RNN, this requires a large excess torque $\tilde{e}_{t}$
to counteract the large error $\Delta\theta_{t}$, as derived from Eq. 18.

On the other hand, when $w^{i}$ is set smaller, the approximate posterior
follows the prior only weakly, allowing larger KL-divergence between them
while the prediction error becomes smaller. In this case, only a small amount
of excess torque is necessary to counteract the small error. The top-down
actional intention of the robot became stronger in the case of larger $w^{i}$
settings; therefore, the human experimenter was required to exert more force
on the robot arms to induce a transition, whereas less force was required with
smaller $w^{i}$ settings, since the top-down actional intention of the robot
became weaker.

In the current experiment, the value of $w^{i}$ was set between $0.01$ and
$0.1$. This is because our preliminary experiments showed that the robot’s
behaviour became noisy when $w^{i}$ was set smaller than $0.01$, since the
approximate posterior could easily deviate from the prior by noise sampling.
On the other hand, it became difficult for the experimenter to initiate a
transition when $w^{i}$ was set larger than $0.1$, because of substantially
increased resistance.

####  3.2.2 Experiment-2

Next, we looked at the difference in the excess torque, prediction error, and
KL-divergence between trained and untrained transitions while $w^{i}$ was
fixed at a given value. Both trained and untrained transitions were attempted
10 times for each of the three trained networks with $w^{i}$ set to $0.01$ for
both trained and untrained cases. However, for the untrained case, among 30
attempts, only 24 succeeded in performing a transition. Our preliminary
experiment showed that the untrained transition became more difficult when
$w^{i}$ was set higher than $0.01$ because of the strong resistance when the
robot attempted to lead trained movement transitions. Fig.7 shows the time-
average of the excess torque, prediction error, and KL-divergence for both
trained and untrained transitions attempted by the experimenter.

Figure 7: The amount of excess torque, prediction error exerted and the KL-
divergence between the approximate posterior and the priorc during each
attempt for trained and untrained transitions. The meta-prior $w^{i}$ was set
to $0.01$.

The time-average of the prediction error, the excess torque, and the KL-
divergence are larger in untrained than trained transitions. This means that
untrained transitions require the experimenter to exert more force than for
trained transitions because of the free energy, which is the sum of the
prediction error and the KL-divergence, and which increases more in untrained
transitions. We provide two supplementary videos for experiment-2 showing the
interaction between the experimenter and Torobo as well as the network
dynamics in the case with meta-prior $w^{i}$ set to 0.01, where trained
transitions (video-link1) and untrained transitions (video-link3) are
performed.

##  4 Discussion

The current study investigated human-robot bodily interactions via
kinaesthesia using a PV-RNN model that was developed based on the free energy
principle. Bodily interactions between a human experimenter and a robot were
conducted using Torobo, a humanoid robot equipped with a PV-RNN model that can
sense excess torque exerted by a human counterpart. We especially examined how
the counter force between the robot and the human experimenter changed during
movement pattern transitions physically guided by the human experimenter,
depending on two different condition changes.

In experiment 1, we examined how the setting of a parameter called the meta-
prior $w^{i}$, which regulates the KL-divergence between the approximate
posterior distribution and the prior distribution in the interaction phase,
affects the counter force generated in executed or attempted transitions.
Results of this experiment showed that in the case of a smaller $w^{i}$, while
the KL-divergence between the approximate posterior and the prior distribution
becomes larger, the prediction error (negative log-likelihood) becomes
smaller. Since the prediction error diminishes, the excess torque to
counteract this error also decreases. On the other hand, in the case of a
larger $w^{i}$ setting, while the KL-divergence between the approximate
posterior and the prior distribution becomes smaller, the prediction error
becomes larger, which requires more excess torque for the transition. The
conflict that appeared between the movement intended by the robot and that
executed by the experimenter is distributed to the prediction error and the
KL-divergence between the approximate posterior and the prior in proportions
determined by the meta-prior $w^{i}$. With larger $w^{i}$ the top-down
movement intention of the robot becomes stronger, which results in a stronger
counter force, whereas the top-down intention as well as the counter force
become weaker with smaller $w^{i}$.

The above is consistent with past research from our group [37]. In [37], Tani
conducted simulation studies on synchronized imitative interaction by dyadic,
vision-based robots using a PV-RNN model. That study showed that a robot with
smaller/larger $w^{i}$ tends to follow or lead the other robot set with a
larger or smaller $w^{i}$ with weaker or stronger actional intention in
synchronized imitative interaction. Similarly in the current study, the
experimenter easily led the robot with a smaller $w^{i}$ with a smaller
counter force because of the weaker top-down intention of the robot. On the
other hand, when $w^{i}$ is quite large, such as $>0.1$ for the robot, it was
difficult for the experimenter to lead the robot because of the extremely
strong counter force. In this situation, the experimenter just followed
movement patterns strongly led by the robot while grasping the robot’s hands,
as shown in the preliminary experiment described previously.

In experiment 2, we examined the difference in the counter force required for
the experimenter to execute trained and untrained movement pattern
transitions. These experimental results showed that untrained transitions
require more force since such transitions are accompanied by larger increases
in free energy, the KL-divergence (between the approximate posterior and the
prior), and the prediction error. Trained transitions, on the other hand,
require less force because of smaller increases in free energy, the KL-
divergence, and the prediction error.

As already mentioned, there have been few studies of human-robot interactions
based on the free energy principle. Although [23] showed that the setting of
$w^{t}$ as meta-prior in the training phase could strongly affect
characteristics of human-robot kinaesthetic interactions, the description of
their experiment results did not include rigorous analysis with repeated
experiments.

The major limitation of the present study is that presented human-robot
interactions are not fully interactive since experimenter-induced sequences of
movement pattern transitions were determined a priori. In this regard, Ikegami
and his colleagues [38, 39] investigated underlying mechanisms for turn-taking
that were generated by spontaneous interaction between artificial agents as
well as artificial agents and humans. Recently, Masumori et al. [40] developed
a humanoid robot platform, called Alter3, behaviour of which was controlled by
sub-modules, including a self-simulator, an automatic mimicry unit, and memory
storage, which were perturbed by a specific neurodynamic model for the purpose
of conducting experiments on spontaneous human-robot interactions. They showed
that spontaneous turn-taking between imitator and imitated could be developed
by autonomous switching of information flow between the two sides.

In future studies, we will undertake human-robot kinaesthetic interaction
experiments that assume less a priori. Such experiments should be done not
with experimenters as counterparts of the robots (as in the present study) but
by inviting an adequate number of human participants, since the human side
also needs to be analyzed. Such studies will focus on two research issues. One
is to investigate how spontaneous turn-taking can occur in imitative
interaction based on kinaesthesis by using the active inference framework [12,
14, 41]. Spontaneous turn-taking in this setting means that the role of the
leader to initiate the next shared patterns switches autonomously between the
two sides, such that sometimes the robot may push hard with its own intended
patterns and the human counterpart may do so at other times. This study may
require development of an autonomous $w^{i}$ adaptation scheme, since if
$w^{i}$ on the robot side can shift adaptively by sensing contextual flow in
the interaction, the leader-follower relationship should shift accordingly.

The other focus is to investigate how novel movement patterns can be developed
through repeated kinaesthetic interaction associated with continuous learning
in both robots and human participants, based on the free energy principle. One
assumption is that novel patterns could develop in terms of false memory as
the number of movement patterns memorized distributively in the PV-RNN model
increases. This phenomenon of false memory is due to potential non-linearity
and stochasticity in the PV-RNN model. A study on a deterministic RNN model
demonstrated this property [42]. Novel patterns generated by robots could
enhance improvisation of new pattern generation from human counterparts
through iterative interaction.

## Prior Generation

Here, in Fig.8, we show part of the result from the prior generation,
including two transitions. As shown in Table 1, since the numbers of d neurons
in layers 2 and 1 were 30 and 60, respectively, we performed principle
component analysis and reduced both dimensions to 3.

Figure 8: Neural activity in prior generation while performing a movement
pattern transition. The two columns at the left show neural activities in
layers 2 and 1, respectively. The top row shows the 3 principal components of
d values, and the bottom row shows the mean value $\mu_{t}^{p}$ of the prior.
The right column shows the joint angle at each time step. The transition was
performed at time step $t=40,340$.

## References

  * [1] Sergio Guadarrama, Lorenzo Riano, Dave Golland, Daniel Go, Yangqing Jia, Dan Klein, Pieter Abbeel, Trevor Darrell, et al.  Grounding spatial relations for human-robot interaction.  In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 1640–1647. IEEE, 2013.
  * [2] Takayuki Kanda, Takayuki Hirano, Daniel Eaton, and Hiroshi Ishiguro.  Interactive robots as social partners and peer tutors for children: A field trial.  Human–Computer Interaction, 19(1-2):61–84, 2004.
  * [3] Joshua Wainer, Ben Robins, Farshid Amirabdollahian, and Kerstin Dautenhahn.  Using the humanoid robot kaspar to autonomously play triadic games and facilitate collaborative play among children with autism.  IEEE Transactions on Autonomous Mental Development, 6(3):183–199, 2014.
  * [4] Francisco J Varela, Evan Thompson, and Eleanor Rosch.  The embodied mind, revised edition: Cognitive science and human experience.  MIT press, 2017.
  * [5] Shikhar Bahl, Abhinav Gupta, and Deepak Pathak.  Human-to-robot imitation in the wild.  arXiv preprint arXiv:2207.09450, 2022.
  * [6] Jörg Krüger, Terje K Lien, and Alexander Verl.  Cooperation of human and machines in assembly lines.  CIRP annals, 58(2):628–646, 2009.
  * [7] Luka Peternel, Wansoo Kim, Jan Babič, and Arash Ajoudani.  Towards ergonomic control of human-robot co-manipulation and handover.  In 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids), pages 55–60. IEEE, 2017.
  * [8] Karl Friston.  A theory of cortical responses.  Philosophical transactions of the Royal Society B: Biological sciences, 360(1456):815–836, 2005.
  * [9] Maxwell JD Ramstead, Karl J Friston, and Inês Hipólito.  Is the free-energy principle a formal theory of semantics? from variational density dynamics to neural and phenotypic representations.  Entropy, 22(8):889, 2020.
  * [10] Jelle Bruineberg, Julian Kiverstein, and Erik Rietveld.  The anticipating brain is not a scientist: the free-energy principle from an ecological-enactive perspective.  Synthese, 195(6):2417–2444, 2018.
  * [11] Karl Friston and Stefan Kiebel.  Predictive coding under the free-energy principle.  Philosophical transactions of the Royal Society B: Biological sciences, 364(1521):1211–1221, 2009.
  * [12] Karl J Friston, Jean Daunizeau, James Kilner, and Stefan J Kiebel.  Action and behavior: a free-energy formulation.  Biological cybernetics, 102(3):227–260, 2010.
  * [13] Karl Friston, Thomas FitzGerald, Francesco Rigoli, Philipp Schwartenbeck, Giovanni Pezzulo, et al.  Active inference and learning.  Neuroscience & Biobehavioral Reviews, 68:862–879, 2016.
  * [14] Thomas Parr and Karl J Friston.  Generalised free energy and active inference.  Biological cybernetics, 113(5):495–513, 2019.
  * [15] Alejandra Ciria, Guido Schillaci, Giovanni Pezzulo, Verena V Hafner, and Bruno Lara.  Predictive processing in cognitive robotics: a review.  Neural Computation, 33(5):1402–1432, 2021.
  * [16] Pablo Lanillos, Cristian Meo, Corrado Pezzato, Ajith Anil Meera, Mohamed Baioumy, Wataru Ohata, Alexander Tschantz, Beren Millidge, Martijn Wisse, Christopher L Buckley, et al.  Active inference in robotics and artificial agents: Survey and challenges.  arXiv preprint arXiv:2112.01871, 2021.
  * [17] Tadahiro Taniguchi, Shingo Murata, Masahiro Suzuki, Dimitri Ognibene, Pablo Lanillos, Emre Ugur, Lorenzo Jamone, Tomoaki Nakamura, Alejandra Ciria, Bruno Lara, et al.  World models and predictive coding for cognitive and developmental robotics: Frontiers and challenges.  arXiv preprint arXiv:2301.05832, 2023.
  * [18] Antonella Maselli, Pablo Lanillos, and Giovanni Pezzulo.  Active inference unifies intentional and conflict-resolution imperatives of motor control.  PLoS computational biology, 18(6):e1010095, 2022.
  * [19] Alexander Tschantz, Manuel Baltieri, Anil K Seth, and Christopher L Buckley.  Scaling active inference.  In 2020 international joint conference on neural networks (ijcnn), pages 1–8. IEEE, 2020.
  * [20] Corrado Pezzato, Carlos Hernández Corbato, Stefan Bonhof, and Martijn Wisse.  Active inference and behavior trees for reactive action planning and execution in robotics.  IEEE Transactions on Robotics, 2023.
  * [21] Jun Tani.  Exploring robotic minds: actions, symbols, and consciousness as self-organizing dynamic phenomena.  Oxford University Press, 2016.
  * [22] Ahmadreza Ahmadi and Jun Tani.  A novel predictive-coding-inspired variational rnn model for online prediction and recognition.  Neural computation, 31(11):2025–2074, 2019.
  * [23] Hendry Ferreira Chame and Jun Tani.  Cognitive and motor compliance in intentional human-robot interaction.  arXiv:1911.01753, 2019.  accepted for publication in IEEE ICRA2020.
  * [24] Nadine Wirkuttis and Jun Tani.  Leading or following? dyadic robot imitative interaction using the active inference framework.  IEEE Robotics and Automation Letters, 6(3):6024–6031, 2021.
  * [25] Ahmadreza Ahmadi and Jun Tani.  A novel predictive-coding-inspired variational rnn model for online prediction and recognition.  Neural computation, 31(11):2025–2074, 2019.
  * [26] Yuichi Yamashita and Jun Tani.  Emergence of functional hierarchy in a multiple timescale neural network model: a humanoid robot experiment.  PLoS computational biology, 4(11), 2008.
  * [27] Léo Pio-Lopez, Ange Nizard, Karl Friston, and Giovanni Pezzulo.  Active inference and robot control: a case study.  Journal of the Royal Society Interface, 16, 2016.
  * [28] Guido Schillaci, Alejandra Ciria, and Bruno Lara.  Tracking emotions: Intrinsic motivation grounded on multi-level prediction error dynamics.  10th Joint IEEE ICDL-EPIROB, pages 1–8, 2020.
  * [29] J. Hwang, J. Kim, A. Ahmadi, M. Choi, and J. Tani.  Dealing with large-scale spatio-temporal patterns in imitative interaction between a robot and a human by using the predictive coding framework.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 50(5):1918–1931, 2020.
  * [30] Junyoung Chung, Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron C Courville, and Yoshua Bengio.  A recurrent latent variable model for sequential data.  In Advances in neural information processing systems, pages 2980–2988, 2015.
  * [31] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [32] Ahmadreza Ahmadi and Jun Tani.  How can a recurrent neurodynamic predictive coding model cope with fluctuation in temporal patterns? robotic experiments on imitative interaction.  Neural Networks, 92:3–16, 2017.
  * [33] Paul J Werbos.  Backpropagation through time: what it does and how to do it.  Proceedings of the IEEE, 78(10):1550–1560, 1990.
  * [34] Timothy P Lillicrap and Adam Santoro.  Backpropagation through time and the brain.  Current opinion in neurobiology, 55:82–89, 2019.
  * [35] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  nature, 323(6088):533–536, 1986.
  * [36] Wataru Ohata and Jun Tani.  Investigation of the sense of agency in social cognition, based on frameworks of predictive coding and active inference: A simulation study on multimodal imitative interaction.  Frontiers in Neurorobotics, 14:61, 2020.
  * [37] Nadine Wirkuttis, Wataru Ohata, and Jun Tani.  Turn-taking mechanisms in imitative interaction: Robotic social interaction based on the free energy principle.  Entropy, 25(2):263, 2023.
  * [38] Takashi Ikegami and Hiroyuki Iizuka.  Turn-taking interaction as a cooperative and co-creative process.  Infant Behavior and Development, 30(2):278–288, 2007.
  * [39] Hiroyuki Iizuka and Takashi Ikegami.  Adaptability and diversity in simulated turn-taking behavior.  Artificial Life, 10(4):361–378, 2004.
  * [40] Atsushi Masumori, Norihiro Maruyama, and Takashi Ikegami.  Personogenesis through imitating human behavior in a humanoid robot “alter3”.  Frontiers in Robotics and AI, 7, 2021.
  * [41] Manuel Baltieri and Christopher L Buckley.  Pid control as a process of active inference with linear generative models.  Entropy, 21(3):257, 2019.
  * [42] Jun Tani, Masato Ito, and Yuuya Sugita.  Self-organization of distributedly represented multiple behavior schemata in a mirror system: reviews of robot experiments using rnnpb.  Neural Networks, 17(8-9):1273–1289, 2004.


propositiontheoremproposition lemmatheoremlemma corollarytheoremcorollary
definitiontheoremdefinition exampletheoremexample notetheoremnote
remarktheoremremark

# Linear $q$-difference, difference and differential operators preserving some
$\mathcal{A}$-entire functions

Jiaxing Huang  College of Mathematics and Statistics, Shenzhen University,
Shenzhen, PR China<EMAIL_ADDRESS>and  Tuen-Wai Ng  Department of
Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong
<EMAIL_ADDRESS>
###### Abstract.

We apply Rossi’s half-plane version of Borel’s Theorem to study the zero
distribution of linear combinations of $\mathcal{A}$-entire functions (Theorem
1.1). This provides a unified way to study linear $q$-difference, difference
and differential operators (with entire coefficients) preserving subsets of
$\mathcal{A}$-entire functions, and hence obtain several analogous results for
the Hermite-Poulain Theorem to linear finite ($q$-)difference operators with
polynomial coefficients. The method also produces a result on the existence of
infinitely many non-real zeros of some differential polynomials of functions
in certain sub-classes of $\mathcal{A}$-entire functions.

Keyword Laguerre-Pólya class; q-difference operators; differential polynomial;
real zeros; Nevanlinna Theory

MSC 2010 Primary 30C15; Secondary 30D35; 30D15

##  1\. Introduction and main results

The investigations of linear operators preserving real-rootedness of certain
classes of entire functions of one complex variable has a long history. In the
1870s, the linear operator preserving the class of hyperbolic polynomials
$\mathcal{HP}$ (i.e. polynomials with real coefficients whose zeros are all
real) was initiated by Hermite, and further developed by Laguerre. In 1914,
Pólya and Schur [30] completely described the operators acting diagonally on
the standard monomial basis $1,x,x^{2},\dots,$ of $\mathbb{R}[x]$ and
preserving $\mathcal{HP}$. One may then consider the corresponding
classification problem for some classes of entire functions containing
$\mathcal{HP}$, for example, the classical Laguerre-Pólya class (see [31,
Definition 5.4.11]).

Let $S$ be a subset in the complex plane $\mathbb{C}$. An entire function $f$
is said to be in the _S-Laguerre-Pólya class_ , $f\in\mathcal{LP}(S)$, if

(1.1) $f(z)=h(z)e^{-\alpha z^{2}+\beta z},\quad
h(z)=cz^{n}\prod_{k=1}^{\infty}(1-z/z_{k})e^{tz/z_{k}}$

where $\beta\in\mathbb{R}$, $c\in\mathbb{R}\setminus\\{0\\}$, $\alpha\geq 0$,
$t=\\{0,1\\}$, $n$ is a non-negative integer and $\\{z_{k}\\}$ is a finite or
infinite sequence in $S$ with $\sum_{k}|z_{k}|^{-t-1}<\infty$. By [15, Theorem
1.11] or [22, Theorem 3.8.5], for $M(r,h)=\displaystyle\max_{|z|=r}|h(z)|$, we
have $\log M(r,h)=o(r^{t+1})$ as $r\rightarrow\infty$.

Any function $f$ in _S-Laguerre-Pólya class_ with $t=0$, $\alpha=0$ and
$\beta\geq 0$ is said to be of $S$-type I. Clearly, _S-Laguerre-Pólya class_
and its sub-class $S$-type I are generalisations of the classical Laguerre-
Pólya class (when $S=\mathbb{R}$ and $t=1$) and type I class (when
$S=\mathbb{R}_{-}:=\\{x\in\mathbb{R}|x<0\\}$ or
$S=\mathbb{R}_{+}:=\\{x\in\mathbb{R}|x>0\\}$) respectively. Notice that by a
theorem of Pólya [28] (see [31, Theorem 5.4.12]), $\mathcal{LP}(\mathbb{R})$
is the closure (in the sense of the uniform convergence on compacta) of
polynomials in $\mathcal{HP}$.

To describe our results, we need to introduce some basic definitions in the
Nevanlinna theory (see [15]) and the $\mathcal{A}$-entire functions (which are
called class $\mathcal{A}$ functions in Chapter V of B. Ja. Levin’s book
[23]).

######  Definition .

A sequence $\\{a_{n}\\}$ of $\mathbb{C}$ is called an $A$-sequence if it
satisfies the condition

(1.2) $\sum_{n=1}^{\infty}\left|\mathrm{Im}\frac{1}{a_{n}}\right|<\infty.$

Here, if $a_{n}=0$, then we define $\mathrm{Im}\dfrac{1}{a_{n}}=0$. An entire
function $f$ is in class $\mathcal{A}$ if its zero set $\\{z_{n}\\}$ is an
$A$-sequence and we will also call such $f$ an $\mathcal{A}$-entire function.
If such $f$ is also in $\mathcal{LP}(\\{z_{n}\\})$ or $\\{z_{n}\\}$-type I
class, then we say $f$ is in class $\mathcal{LP}(A)$ or $A$-type I class
respectively. Finally, $\mathcal{LP}(A;2)$ is the subset of $\mathcal{LP}(A)$
which contains $f$ with $\alpha>0$ in (1.1) and $\mathcal{LP}_{t_{0}}(A;2)$
($\mathcal{LP}_{t_{0}}(A)$) is the subset of $\mathcal{LP}(A;2)$
($\mathcal{LP}(A)$) which contains $f$ with $t=t_{0}$ in (1.1).

Clearly, class $\mathcal{A}$ contains entire functions with only real zeros
and hence $\mathcal{HP}\subset\mathcal{LP}(\mathbb{R})\subset\mathcal{A}$.
Also, by definition, we have
$\mathcal{LP}_{t_{0}}(A;2)\subset\mathcal{LP}(A;2)\subset\mathcal{LP}(A)\subset\mathcal{A}$.

For any meromorphic function $f$ on $\mathbb{C}$, let $n(r,f)$ be the number
of poles of $f$ in $|z|\leq r$, and

$N(r,f):=\int_{1}^{r}\frac{n(t,f)-n(0,f)}{t}dt+n(0,f)\log r$

be the counting function of $f$. The proximity function $m(r,f)$ is defined by

$m(r,f):=\frac{1}{2\pi}\int_{0}^{2\pi}\log^{+}|f(re^{i\theta})|d\theta,$

where $\log^{+}a:=\max(\log a,0)$. The Nevanlinna characteristic function
$T(r,f)$ is defined by

$T(r,f):=m(r,f)+N(r,f).$

We also introduce the exponent of convergence of the zeros of $f$,
$\lambda(f)$, and the order $\rho(f)$ of $f$, which are given respectively by

$\lambda(f):=\limsup_{r\rightarrow\infty}\frac{\log N(r,1/f)}{\log
r}=\limsup_{r\rightarrow\infty}\frac{\log n(r,1/f)}{\log r}.$

and

$\rho(f):=\limsup_{r\rightarrow\infty}\frac{\log T(r,f)}{\log r}.$

It is clear that $\lambda(f)\leq\rho(f)$. If $f$ is an entire function, it is
not hard to see that $T(r,f)\leq\log^{+}M(r,f)$, for all $r>0$. Finally, by
$\mathcal{M}_{k}$ we mean the field of meromorphic functions $f$ with
$T(r,f)=o(r^{k})$ as $r\to\infty$ outside a set of finite measure and by
$\mathbb{C}[z]$ we mean the ring of polynomials with complex number
coefficients. The field $\mathcal{M}_{k}$ appears naturally in the studies of
some hypertranscendental functions (see for example [17]).

Using some ideas from Eremenko and Rubel [14] and Ng and Yang [26], and the
half-plane version of Borel’s lemma by Rossi [32] (see Lemma 2), we obtain the
following

######  Theorem 1.1.

Let $f_{1},\ldots,f_{n}$ be linearly independent entire functions over
$\mathcal{M}_{2}$ satisfying $N(r,1/f_{i})=o(r^{2})$ for all $i$. Let
$a_{1},\ldots,a_{n}$ be entire functions in $\mathcal{M}_{1}$. Suppose that
each $f_{i}$ is in class $\mathcal{A}$. Then

$F=a_{1}f_{1}+\dots+a_{n}f_{n}$

is in $\mathcal{A}\backslash\\{0\\}$ if and only if $a_{i}\not\equiv 0$ for at
most one $i$ and this $a_{i}$ is in class $\mathcal{A}\backslash\\{0\\}$.

###  1.1. Linear operators preserving real-rootedness

We now explain how Theorem 1.1 can be used to classify linear ($q$-difference,
difference, differential) operators of finite order preserving some sub-
classes of $\mathcal{A}$-entire functions.

Let $\mathcal{M}$ be the field of meromorphic functions on $\mathbb{C}$.
Consider the linear difference operator
$\Delta_{M_{1},M_{2},h}:\mathcal{M}\rightarrow\mathcal{M}$ defined by

(1.3) $\Delta_{M_{1},M_{2},h}(f)(z):=M_{1}(z)f(z+ih)+M_{2}(z)f(z-ih),$

where $M_{1}$ and $M_{2}$ are complex-valued functions, and $h$ is a complex
number. In 1926, Pólya [29] established that
$\Delta_{1,1,c}(\mathcal{LP}(\mathbb{R}))\subset\mathcal{LP}(\mathbb{R})$ for
every real number $c$ and de Bruijn [13] noticed that actually
$\Delta_{e^{i\theta},e^{-i\theta},c}(\mathcal{LP}(\mathbb{R}))\subset\mathcal{LP}(\mathbb{R})$
for every real number $c$ and $\theta\in[0,2\pi]$. These studies were
continued by Oberschkov [27], Levin [23], Craven and Csordas [11, 12], Walker
[35] and others. On the other hand, Brändén and Borcea [6] have completely
characterised all linear operators on $\mathbb{C}[z]$ preserving
$\mathcal{HP}$ as well as $\mathcal{LP}(\mathbb{R})$ in the entire function
space. However, it is not easy to apply their characterisation to check when
the two term difference operator (1.3) preserves $\mathcal{LP}(\mathbb{R})$.
In 2017, Katkova et al. [18] thoroughly solved this classification problem and
hence generalised the results of Pólya and de Bruijn. They gave the necessary
and sufficient conditions of the operator (1.3) preserving the class
$\mathcal{LP}(\mathbb{R})$ in terms of the explicit expressions of $M_{1}$ and
$M_{2}$ (see [18, Theorem 1.1]).

It is then natural to consider linear difference operators with more than two
terms. When the coefficients are polynomials, by applying Brändén and Borcea’s
characterisation of linear operators on $\mathbb{C}[z]$ preserving
$\mathcal{HP}$ ([6, Theorem 1]), Brändén et al. [8] obtained the following

######  Theorem 1.2 (Brändén-Krasikov-Shapiro [8]).

Let $T:\mathbb{C}[z]\to\mathbb{C}[z]$ be a linear operator defined by

(1.4) $T(p)(z)=q_{0}(z)p(z)+q_{1}(z)p(z-1)+\dots+q_{k}(z)p(z-k),$

where $q_{0},\dots,q_{k}$ are fixed complex-valued polynomials. Then
$T(\mathcal{HP})\subset\mathcal{HP}$ if and only if $q_{i}\not\equiv 0$ for at
most one $i$, and this $q_{i}$ is in $\mathcal{HP}$.

Recall that the well-known Hermite-Poulain theorem [27, page 4] states that a
finite order linear differential operator
$T:=a_{0}+a_{1}d/dx+\dots+a_{k}d^{k}/dx^{k}$ preserves the class
$\mathcal{HP}$ if and only if its symbol polynomial
$Q_{T}(t)=a_{0}+a_{1}t+\dots+a_{k}t^{k}$ is hyperbolic. Thus, Theorem 1.2 is a
somewhat difference operator analog of the Hermite-Poulain theorem. Inspired
by this result of Brändén-Krasikov-Shapiro, we consider finite order linear
$q$-difference, difference and differential operators with possibly non-
polynomial coefficients, and obtain three corollaries of Theorem 1.1.

######  Corollary .

Let $q$ be a nonzero real number such that $|q|\neq 1$ and $p\in\mathcal{M}$.
Let $T_{1}:\mathcal{M}\to\mathcal{M}$ be the finite $q$-difference operator
defined by

$T_{1}(p(z)):=a_{0}p(z)+a_{1}p(qz)+\dots+a_{k}p(q^{k}z),$

where $a_{i}\in\mathcal{M}_{1}$ are entire, for $i=0,1,\dots,k$. Then, the
following statements hold:

  1. (1)

Suppose there exists some $p\in\mathcal{LP}(A;2)$ such that
$T_{1}(p)\in\mathcal{LP}(A;2)$. Then $a_{i}\not\equiv 0$ for at most one $i$.

  2. (2)

The finite $q$-difference operator $T_{1}$ preserves the class
$\mathcal{LP}_{0}(A;2)$ if and only if $a_{i}\not\equiv 0$ for at most one
$i$, and this $a_{i}$ is in the $A$-type I class.

  3. (3)

Suppose $a_{i}\in\mathbb{C}[z]$ for all $0\leq i\leq k$, then
$T_{1}(\mathcal{LP}(\mathbb{R};2))\subset\mathcal{LP}(\mathbb{R};2)$ if and
only if $a_{i}\not\equiv 0$ for at most one $i$, and this
$a_{i}\in\mathcal{HP}$.

######  Corollary .

Let $T_{2}:\mathcal{M}\to\mathcal{M}$ be the finite linear difference operator
defined by

(1.5) $T_{2}(p(z)):=a_{0}p(z)+a_{1}p(z+c_{1})+\dots+a_{k}p(z+c_{k})$

where $c_{i}\in\mathbb{R}\setminus\\{0\\}$. Suppose that each $a_{i}$ is of
the form

(1.6) $a_{i}(z)=g_{i}(z)e^{-\mu_{i}z^{2}+d_{i}z}$

where $\mu_{i}$’s are mutually distinct complex numbers with
$\mathrm{Re}\,\mu_{i}>0$, $d_{i}\in\mathbb{C}$ and $g_{i}$ is entire in
$\mathcal{M}_{1}$, for $i=0,\dots,k$. Then the following assertions hold:

  1. (1)

If there exists some $p\in\mathcal{LP}(A)$ such that
$T_{2}(p)\in\mathcal{LP}(A)$, then $a_{i}\not\equiv 0$ for at most one $i$,
and this $a_{i}\in\mathcal{LP}(A;2)$.

  2. (2)

Let $t=0$ or $1$, the finite linear difference operator $T_{2}$ preserves the
class $\mathcal{LP}_{t}(A)$ if and only if $a_{i}\not\equiv 0$ for at most one
$i$, and this $a_{i}\in\mathcal{LP}_{t}(A;2)$.

######  Corollary .

Let $T_{3}$ be the finite linear differential operator given by

(1.7) $T_{3}(p(z)):=a_{0}p(z)+a_{1}p^{\prime}(z)+\dots+a_{k}p^{(k)}(z).$

Suppose that each $a_{i}$ is of the form (1.6) where $\mu_{i}$’s are mutually
distinct complex numbers with $\mathrm{Re}(\mu_{i})>0$, $d_{i}\in\mathbb{C}$
and $g_{i}$ is entire in $\mathcal{M}_{1}$, for $i=0,\dots,k$. If
$T_{3}(\mathcal{LP}(\mathbb{R}))\subset\mathcal{LP}(\mathbb{R})$, then
$a_{i}\not\equiv 0$ for at most one $i$ and this
$a_{i}\in\mathcal{LP}(\mathbb{R};2)$.

Corollaries 1.1-1.1 give necessary and/or sufficient conditions of such
operators preserving some sub-classes of $\mathcal{A}$-entire functions. In
fact, Corollary 1.1 could be seen as a generalisation of the Hermite-Poulain
theorem (á la work [8]) in the $q$-difference context. Corollary 1.1 is a
transcendental version of Theorem 1.2 and sort of a complement to the results
of [18] and to the Hermite-Poulain theory, developed in [19], for difference
operators of finite order. Finally, Corollary 1.1 is somewhat an extension of
the Hermite-Poulain theorem to some sub-classes of $\mathcal{A}$-entire
functions.

###  1.2. Zeros of differential polynomials

In 1989, Sheil-Small [33] settled a longstanding conjecture of Wiman (1911).
As a consequence, Sheil-Small obtained that _if $f$ is a real entire function
(mapping the real line to itself) of finite order and $ff^{\prime\prime}$ has
no non-real zeros, then $f\in\mathcal{LP}(\mathbb{R})$._ This result also
solves a problem posed by Hellerstein (see [9, p.552, Problem 4.28]). Later,
Bergweiler et al.[2] completed Sheil-Small’s result to the real entire
function $f$ with infinite order: _for every real entire function $f$ of
infinite order, $ff^{\prime\prime}$ has infinitely many non-real zeros_. Thus,
combining the results of Sheil-Small and Bergweiler et al., we have the
following

######  Theorem 1.3 (Sheil-Small [33], Bergweiler-Eremenko-Langley [2]).

Let $f$ be a real entire function and $ff^{\prime\prime}$ has only real zeros,
then $f\in\mathcal{LP}(\mathbb{R})$.

Applying a result of Bergweiler et al. to the function of the form

$f=\exp\int_{0}^{z}g(t)dt,$

one can obtain the following result which also follows from a result of
Bergweiler and Fuchs [4].

######  Theorem 1.4 (Bergweiler-Fuchs [4]; Bergweiler-Eremenko-Langley [2]).

For every real transcendental entire function $g$, the function
$g^{\prime}+g^{2}$ has infinitely many non-real zeros.

In 2005, Bergweiler et al. [3] extended Theorem 1.4 to the real meromorphic
functions, and considered the zeros of $f^{\prime}+f^{m}$ where $m\geq 3$. It
is natural to ask if $f^{\prime}$ can be replaced by any linear differential
polynomial of $f$ or linear difference polynomial of $f$ (Langley [21, page
108] asked a similar question when $f$ is a real entire function with finitely
many non-real zeros). In general, this is not true. For example, let
$f=e^{-z}$ which is a real entire function, then
$f^{\prime\prime}+f^{\prime}+f^{m}=e^{-mz}$ has no zeros for any integer $m$.
However, applying Theorem 1.1, we do have a positive result if we restrict to
certain sub-classes of real entire functions.

######  Corollary .

Let $P$ be a complex polynomial with degree at least two and $P(0)=0$. If $f$
is in class $\mathcal{A}$ with $2\leq\rho(f)<\infty$ and $N(r,1/f)=o(r)$, and
$L$ is a non-constant linear differential operator with coefficients in
$\mathbb{C}$, then $L(f^{\prime})+P(f)$ is not in class $\mathcal{A}$. Hence
$L(f^{\prime})+P(f)$ has infinitely many non-real zeros. In particular, this
result holds for $f$ in $\mathcal{LP}(A;2)$ with $N(r,1/f)=o(r)$.

######  Remark .

The condition that $2\leq\rho(f)<\infty$ is necessary as can be seen from the
above example for $f=e^{-z}$. It would be interesting to see if the condition
$N(r,1/f)=o(r)$ is also necessary.

######  Remark .

Essentially the same proof also works when the above differential operator
$L(f^{\prime})$ is replaced by the linear difference operator
$a_{n}f(z+c_{n})+\cdots+a_{1}f(z+c_{1})$ where $a_{i}\in\mathbb{C}$ and those
$c_{i}$ are mutually distinct nonzero constants.

######  Remark .

In Corollary 1.2, we assume the degree of $P$ is at least two. When $\deg
P=1$, Langley [21, Theorem 1.4] showed that if $f$ is an infinite order real
entire function with finitely many non-real zeros, then
$f^{\prime\prime}+\omega f$ has infinitely many non-real zeros for any
positive $\omega$.

The rest of the paper is organized as follows. In Sect. 2 we state several
results that will be used in our proofs. Then we prove our main result
(Theorem 1.1) in Sect. 3, Corollaries 1.1–1.1 in Sect. 4 and finally Corollary
1.2 in Sect. 5.

##  2\. Some Lemmata

As we will apply the half-plane version of Borel’s lemma by Rossi [32] (see
Lemma 2) to prove Theorem 1.1, we first introduce Tsuji’s characteristic of a
meromorphic function in the upper (lower) half-plane (see [34]).

Let $\mathfrak{n}_{u}(t,f)$ be the number of poles of $f$ in
$\\{z:|z-it/2|\leq t/2,|z|\geq 1\\}$, where $f$ is meromorphic in the open
upper half plane. Define

$\mathfrak{N}_{u}(r,f):=\int_{1}^{r}\frac{\mathfrak{n}_{u}(t,f)}{t^{2}}dt=\sum_{1\leq
r_{k}\leq
r\sin\theta_{k}}\left(\frac{\sin\theta_{k}}{r_{k}}-\frac{1}{r}\right),$
$\mathfrak{m}_{u}(r,f):=\frac{1}{2\pi}\int_{\arcsin(r^{-1})}^{\pi-\arcsin(r^{-1})}\log^{+}|f(r\sin\theta
e^{i\theta})|\frac{d\theta}{r\sin^{2}\theta},$

and

$\mathfrak{T}_{u}(r,f):=\mathfrak{N}_{u}(r,f)+\mathfrak{m}_{u}(r,f)$

where $r_{k}e^{i\theta_{k}}$ are the poles of $f$ in
$\\{z:\mathrm{Im}\,z>0\\}$. Similarly, one can also define
$\mathfrak{m}_{l}(r,f),\mathfrak{N}_{l}(r,f)$ and $\mathfrak{T}_{l}(r,f)$ for
functions meromorphic in the open lower half plane.

######  Lemma  ([24]).

Let $f$ be meromorphic in the open upper (lower) half plane. Define

$m_{\alpha,\beta}(r,f):=\frac{1}{2\pi}\int_{\alpha}^{\beta}\log^{+}|f(re^{i\theta})|d\theta.$

Then

$\int_{r}^{\infty}\frac{m_{0,\pi}(t,f)}{t^{3}}dt\leq\int_{r}^{\infty}\frac{\mathfrak{m}_{u}(t,f)}{t^{2}}dt$
$\left(\int_{r}^{\infty}\frac{m_{\pi,2\pi}(t,f)}{t^{3}}dt\leq\int_{r}^{\infty}\frac{\mathfrak{m}_{l}(t,f)}{t^{2}}dt\right).$

######  Lemma  ([32]).

Let $n\geq 2$, $G=\\{f_{0},\dots,f_{n}\\}$ be a set of meromorphic functions
in $\mathrm{Im}\,z>0$ such that any proper subset of $G$ is linearly
independent over $\mathbb{C}$. If $G$ is linearly dependent over $\mathbb{C}$,
then for all positive $r$ except possibly a set of finite measure,

$\mathfrak{T}_{u}(r)=O\left(\sum_{k=0}^{n}\left(\mathfrak{N}_{u}(r,f_{k})+\mathfrak{N}_{u}\left(r,\frac{1}{f_{k}}\right)\right)+\log\mathfrak{T}_{u}(r)+\log
r\right),$

where $\mathfrak{T}_{u}(r):=\max\\{\mathfrak{T}_{u}(r,f_{i}/f_{j})|\,0\leq
i,j\leq n\\}$.

We also need the following generalisation [25] of Borel’s lemma [7].

######  Lemma  (Corollary 4.5 of [1]).

Let $a_{i}$ and $g_{i}$, $i=1,\dots,k$ be nonzero meromorphic and entire
functions in $\mathbb{C}$ respectively, satisfying

$a_{1}e^{g_{1}}+a_{2}e^{g_{2}}+\dots+a_{k}e^{g_{k}}\equiv 0.$

If $T(r,a_{j})=o(T(r,e^{g_{m}-g_{l}}))$, for any $m\neq l$ and $1\leq j\leq
k$, then $a_{j}(z)\equiv 0$ for all $1\leq j\leq k$.

##  3\. Proof of Theorem 1.1

We may assume that $a_{i}\neq 0$ for any $i=1,\dots,n$ (otherwise, we can
relabel $a_{i}$ so that we can replace $n$ by a smaller number). Let $g_{0}=F$
and $g_{i}=a_{i}f_{i}$ for $i=1,\dots,n$. By the assumption that
$f_{1},\dots,f_{n}$ are linearly independent over $\mathcal{M}_{2}$ and
$a_{i}\in\mathcal{M}_{1}\subset\mathcal{M}_{2}$, we have $g_{1},\dots,g_{n}$
are linearly independent over $\mathbb{C}$. Consider the set
$G=\\{g_{0},\dots,g_{n}\\}$, then $G$ is linearly dependent over $\mathbb{C}$
and any proper subset of $G$ is linearly independent over $\mathbb{C}$.
Therefore, $G$ satisfies the assumptions of Lemma 2.

Suppose that $F=a_{1}f_{1}+\dots+a_{n}f_{n}\in\mathcal{A}\backslash\\{0\\}$.
Since all $g_{0},f_{1},\dots,f_{n}$ are in $\mathcal{A}$ (hence entire), it is
easy to check from the definitions of class $\mathcal{A}$ and
$\mathfrak{N}_{u}(r,*)$ that

$\mathfrak{N}_{u}(r,1/g_{0})=O(1)\ \mathrm{and}\
\mathfrak{N}_{u}(r,1/f_{i})=O(1).$

We also have for $i=1,\dots,n$, $\mathfrak{N}_{u}(r,g_{i})=0$ because each
$a_{i}$ is entire.

Recall that $T(r,a_{i})=o(r)$. Then for each $i=1,\dots,n$,

$\displaystyle\mathfrak{N}_{u}(r,1/g_{i})+\mathfrak{N}_{u}(r,g_{i})$
$\displaystyle\leq\mathfrak{N}_{u}(r,1/a_{i})+\mathfrak{N}_{u}(r,1/f_{i})+O(1)$
$\displaystyle\leq N(r,1/a_{i})+O(1)\leq T(r,a_{i})=O(r^{\epsilon})$

for some positive $\epsilon<1$. We also have
$\mathfrak{N}_{u}(r,1/g_{0})+\mathfrak{N}_{u}(r,g_{0})=O(1)$. Therefore, we
can deduce from Lemma 2, that

$\mathfrak{T}_{u}(r)=O(r^{\epsilon}),\ \text{and hence}\
\mathfrak{T}_{u}(r,g_{i}/g_{j})=O(r^{\epsilon})$

for $i,j=0,\dots,n$. From the definition of $\mathfrak{T}_{u}(r,g_{i}/g_{j})$,
it follows that

$\mathfrak{m}_{u}(r,g_{i}/g_{j})=O(r^{\epsilon}).$

Similarly,

$\mathfrak{m}_{l}(r,g_{i}/g_{j})=O(r^{\epsilon}).$

Notice that
$m(t,g_{i}/g_{j})=m_{0,\pi}(t,g_{i}/g_{j})+m_{\pi,2\pi}(t,g_{i}/g_{j})$. Then
Lemma 2 implies that for any $r>0$,

$\displaystyle\frac{m(r,g_{i}/g_{j})}{2r^{2}}$
$\displaystyle\leq\int_{r}^{\infty}\frac{m(t,g_{i}/g_{j})}{t^{3}}dt\leq
O(r^{\epsilon-1})$

and hence $m(r,g_{i}/g_{j})=o(r^{2})$.

Using the fact that $a_{i}\in\mathcal{M}_{1},N(r,f_{i})=o(r^{2})$ and
$a_{i},f_{i}$ are entire, we have

$N(r,g_{i}/g_{j})\leq N(r,1/a_{j})+N(r,1/f_{j})\leq
T(r,1/a_{j})+N(r,1/f_{j})=o(r^{2}).$

Therefore, $T(r,g_{i}/g_{j})=o(r^{2})$.

Take $i=0$ and $j=1$, then $T(r,\frac{F}{a_{1}f_{1}})=o(r^{2})$. Let
$b=\frac{F}{a_{1}f_{1}}$ so that $F=ba_{1}f_{1}$ and
$T(r,(1-b)a_{1})=o(r^{2})$. Hence

$(1-b)a_{1}f_{1}+a_{2}f_{2}+\dots+a_{n}f_{n}=0.$

As $f_{1},\dots,f_{n}$ are linearly independent over $\mathcal{M}_{2}$, we
must have $a_{i}=0$ for each $i\neq 1$ and $(1-b)a_{1}=0$. Thus, the only
possibility is $n=1$ and $b=1$ so that $F=a_{1}f_{1}$ and
$a_{1}\in\mathcal{A}$.

The converse assertion is clearly true and this completes the proof.

##  4\. Proof of Corollaries 1.1, 1.1 and 1.1

###  4.1. Proof of Corollary 1.1$(\ref{1})$

Since $p\in\mathcal{LP}(A;2)$, we have $p(z)=h(z)e^{-\alpha z^{2}+\beta z}$
where $\alpha>0$ and $h(z)$ is of the form (1.1) so that
$h\in\mathcal{A}\backslash\\{0\\}$ and $T(r,h)=o(r^{2})$. Let
$b_{i}(z)=h(q^{i}z)$ and $f_{i}(z)=b_{i}(z)e^{-\alpha q^{2i}z^{2}+\beta
q^{i}z}$ for $i=0,\dots,k$ where $|q|\neq 0,1$.

Clearly $f_{i}\in\mathcal{A}\backslash\\{0\\}$ if
$b_{i}\in\mathcal{A}\backslash\\{0\\}$. To see
$b_{i}\in\mathcal{A}\backslash\\{0\\}$, let $\\{\beta_{k}\\}$ and
$\\{z_{k}\\}$ be zeros of $b_{i}$ and $h$ respectively. Since $q$ is real, it
follows that

$\displaystyle\sum|\mathrm{Im}\frac{1}{\beta_{k}}|=|q|^{i}\sum|\mathrm{Im}\frac{1}{z_{k}}|<\infty$

and hence $b_{i}\in\mathcal{A}$. Since $\mathcal{LP}(A;2)\subset\mathcal{A}$,
each $f_{i}$ and $T_{1}(p)=a_{0}f_{0}+\cdots+a_{k}f_{k}$ are in class
$\mathcal{A}\backslash\\{0\\}$. Finally, we notice that
$N(r,1/f_{i})=N(r,1/b_{i})=N(|q|^{i}r,1/h)+O(1)\leq
T(|q|^{i}r,h)+O(1)=o(r^{2})$.

To show that at most one $a_{i}\not\equiv 0$ by applying Theorem 1.1, it
remains to prove that $f_{0},\dots,f_{k}$ are linearly independent over
$\mathcal{M}_{2}$. Suppose there exist $c_{i}\in\mathcal{M}_{2}$, for
$i=0,\dots,k$, such that $c_{0}f_{0}+\cdots+c_{k}f_{k}=0$. Let
$g_{i}=e^{-\alpha q^{2i}z^{2}+\beta q^{i}z}$. As

$\frac{g_{i}}{g_{j}}=e^{-(\alpha q^{2i}-\alpha q^{2j})z^{2}+(\beta q^{i}-\beta
q^{j})z}$

and $\alpha q^{2i}-\alpha q^{2j}\neq 0$ for $i\neq j$, we have whenever $i\neq
j$,

$T(r,g_{i}/g_{j})\geq Cr^{2},\ \text{as}\ r\rightarrow\infty,\ C>0.$

On the other hand, one can check easily that
$T(r,h(q^{j}z))=T(|q|^{j}r,h)+O(1)$ (see [5, Page 249]) so that
$T(r,h(q^{j}z))=o(r^{2})$ for all $0\leq j\leq k$ and hence
$T(r,c_{j}b_{j})=o(r^{2})$ as $c_{j}\in\mathcal{M}_{2}$. Since

$c_{0}b_{0}e^{-\alpha z^{2}+\beta z}+\cdots+c_{k}b_{k}e^{-\alpha
q^{2k}z^{2}+\beta q^{k}z}=c_{0}f_{0}+\cdots+c_{k}f_{k}=0,$

by Lemma 2, $c_{j}b_{j}=0$ for all $j=0,\dots,k$ which implies that $c_{j}=0$
for all $j$ and therefore $f_{0},\dots,f_{k}$ are linearly independent over
$\mathcal{M}_{2}$.

###  4.2. Proof of Corollary 1.1$(\ref{2})$

Since
$T_{1}(\mathcal{LP}_{0}(A;2))\subset\mathcal{LP}_{0}(A;2)\subset\mathcal{LP}(A;2)$,
it follows from part one that $a_{i}\not\equiv 0$ for at most one $i$.
Therefore, it remains to prove that this $a_{i}$ is in class $A$-type I.

Without loss of generality, we may assume that $i=0$ and hence
$T_{1}(p(z))=a_{0}(z)p(z)=a_{0}(z)h(z)e^{-\alpha z^{2}+\beta z}$, where
$h(z)=cz^{n}\displaystyle\prod_{k=1}^{\infty}(1-z/z_{k})\in\mathcal{M}_{1}$,
$c\in\mathbb{R}\backslash\\{0\\}$. Since
$T_{1}(p(z))\in\mathcal{LP}_{0}(A,2)$, we have

(4.1) $a_{0}(z)h(z)e^{-\alpha z^{2}+\beta
z}=h_{1}(z)e^{-\alpha_{1}z^{2}+\beta_{1}z},$

where
$h_{1}(z)=c_{1}z^{m}\displaystyle\prod_{k=1}^{\infty}(1-z/w_{k})\in\mathcal{M}_{1}$,
$c_{1}\in\mathbb{R}\backslash\\{0\\}$. Since $a_{0}h$ and $h_{1}$ are in
$\mathcal{M}_{1}$, we can apply Lemma 2 to conclude that $a_{0}h=h_{1}$ and
hence

$a_{0}=c^{\prime}z^{l}\displaystyle\prod_{w_{k}\neq
z_{k}}(1-z/w_{k})\quad\mathrm{for\ some}\
c^{\prime}\in\mathbb{R}\setminus\\{0\\}\ \mathrm{and}\ l\in\mathbb{N}.$

Since

$\sum_{w_{k}\neq
z_{k}}|w_{k}|^{-1}\leq\sum_{k=1}^{\infty}|w_{k}|^{-1}<\infty\quad\mathrm{and}\quad\sum_{w_{k}\neq
z_{k}}|\mathrm{Im}\frac{1}{w_{k}}|\leq\sum_{k=1}^{\infty}|\mathrm{Im}\frac{1}{w_{k}}|<\infty,$

it follows that $a_{0}$ is in class $A$-type I.

The converse assertion is obvious and we complete the proof.

###  4.3. Proof of Corollary 1.1$(\ref{3})$

Since $a_{i}\in\mathbb{C}[z]\subset\mathcal{M}_{1}$ for all $i$, and
$T_{1}(\mathcal{LP}(\mathbb{R};2))\subset\mathcal{LP}(\mathbb{R};2)$, by
Corollary 1.1(1), we conclude that $a_{i}\not\equiv 0$ for at most one $i$ and
there is no harm to assume that only $a_{0}\not\equiv 0$. The property that
$a_{0}\in\mathcal{HP}$ then follows from the identity (4.1).

###  4.4. Proof of Corollary 1.1$(\ref{a})$

Let $p=h(z)e^{-\alpha z^{2}+\beta z}\in\mathcal{LP}(A)$ be expressed in the
form of (1.1) where $T(r,h)=o(r^{2})$. Let $c_{0}=0$ so that we can write
$T_{2}(p)$ as

$T_{2}(p)=b_{0}f_{0}+b_{1}f_{1}+\dots+b_{k}f_{k}$

where $b_{i}=g_{i}$ and
$f_{i}(z)=h(z+c_{i})e^{-(\alpha(z+c_{i})^{2}+\mu_{i}z^{2})+\beta(z+c_{i})+d_{i}z},$
with $\alpha+\mu_{i}\neq 0$ for all $i$ because $\mathrm{Re}(\mu_{i})>0$ and
$\alpha\geq 0$. Since $N(r,1/h)\leq T(r,h)+O(1)=o(r^{2})$, $\lambda(h)\leq 2$.
By [10, Theorem 2.2], it follows that
$N(r,1/f_{i})=N(r,1/h(z+c_{i}))=N(r,1/h)+O(r^{\lambda(h)-1+\epsilon})+O(\log
r)$ for any positive $\epsilon$. Hence $N(r,1/f_{i}))=o(r^{2})$ for all $i$.

Since each $c_{i}$ is real, one can check that if $|z_{k}|>2|c_{i}|$, then

$|\mathrm{Im}\frac{1}{z_{k}-c_{i}}|\leq 4|\mathrm{Im}\,\frac{1}{z_{k}}|$

and hence $f_{i}\in\mathcal{A}$. Now suppose $T_{2}(p)\in\mathcal{LP}(A)$. In
order to apply Theorem 1.1 to show that at most one $a_{i}\not\equiv 0$, we
only need to show that $f_{0},\dots,f_{k}$ are linearly independent over
$\mathcal{M}_{2}$. This suffices to show that each $h(z+c_{i})$ is in
$\mathcal{M}_{2}$ because by Lemma 2 and the fact that $\mu_{i}$’s are
distinct,
$\\{e^{-(\alpha(z+c_{i})^{2}+\mu_{i}z^{2})+\beta(z+c_{i})+d_{i}z}:i=0,\ldots,k\\}$
is linearly independent over $\mathcal{M}_{2}$.

Since $T(r,h)=o(r^{2})$, its order $\sigma$ is at most two. By [10, Theorem
2.1], it follows that $T(r,h(z+c_{i}))=T(r,h)+O(r^{\sigma-1+\epsilon})+O(\log
r)$ for any positive $\epsilon$. Hence $T(r,h(z+c_{i}))=o(r^{2})$ for all $i$
and we are done.

###  4.5. Proof of Corollary 1.1$(\ref{b})$

The argument is similar to that of Corollary 1.1(2), and we omit the details
of the proof.

###  4.6. Proof of Corollary 1.1

Let $p(z)=h(z)e^{-\alpha z^{2}+\beta z}\in\mathcal{LP}(\mathbb{R})$, where $h$
is given by (1.1) with $T(r,h)=o(r^{2})$. Then $p^{(i)}(z)=h_{i}(z)e^{-\alpha
z^{2}+\beta z}$, where $h_{i}(z)=h_{i-1}^{\prime}(z)+(-2\alpha
z+\beta)h_{i-1}(z)$ and $h_{0}(z)=h(z)$. By [16, Theorem A], it follows that
$p^{(i)}\in\mathcal{LP}(\mathbb{R})\subset\mathcal{A}$. Hence
$f_{i}:=p^{(i)}e^{-\mu_{i}z^{2}+d_{i}z}=h_{i}e^{-(\alpha+\mu_{i})z^{2}+(\beta+d_{i})z}\in\mathcal{A}$.
As $T(r,h_{0})=o(r^{2})$ and $T(r,h_{i})=O(T(r,h_{i-1}))$, we have
$T(r,h_{i})=o(r^{2})$ for all $i$. Therefore, $N(r,1/f_{i})=o(r^{2})$ for all
$i$. The rest can be done as in the proof of Corollary 1.1(1).

##  5\. Proof of Corollary 1.2

If $L(f^{\prime})+P(f)\equiv 0$, then $L(f^{\prime})+P(f)\notin\mathcal{A}$
and we are done. So we may assume that $L(f^{\prime})+P(f)\not\equiv 0$ and
write it in the following form

$\displaystyle L(f^{\prime})+P(f)$
$\displaystyle=\sum_{k=1}^{n}a_{k}f^{(k)}+\sum_{m=1}^{l}b_{m}f^{m}=(\sum_{k=1}^{n}a_{k}\frac{f^{(k)}}{f})f+\sum_{m=1}^{l}b_{m}f^{m}$
$\displaystyle=\sum_{k=1}^{l}c_{k}f^{k}$

where $a_{i}\in\mathbb{C}$,
$c_{1}=b_{1}+\sum_{k=1}^{n}a_{k}\frac{f^{(k)}}{f}$, $c_{m}=b_{m}\in\mathbb{C}$
for $m=2,\dots,l$. Let $f_{j}=f^{j}$ for $j=1,\dots,l$, then

$f_{j}\in\mathcal{A}\quad\mathrm{and}\quad N(r,1/f_{j})=jN(r,1/f)=o(r)$

and hence $N(r,1/f_{j})=o(r^{2})$. As $\rho(f)\geq 2$, by Valiron’s Theorem
([20, Theorem 2.2.5]), $f_{1},\dots,f_{l}$ are linearly independent over
$\mathcal{M}_{2}$. Since $f$ is entire with $\rho(f)<\infty$ and
$N(r,1/f)=o(r)$, it follows from the logarithmic derivative lemma and
$N(r,f^{(k)}/f)\leq N(r,1/f)$ that

$\displaystyle T(r,f^{(k)}/f)$ $\displaystyle=m(r,f^{(k)}/f)+N(r,f^{(k)}/f)$
$\displaystyle=o(r).$

Hence $f^{(k)}/f\in\mathcal{M}_{1}$ and so is $c_{1}$. Now suppose
$L(f^{\prime})+P(f)=c_{1}f+\cdots+c_{k}f^{k}\in\mathcal{A}\backslash\\{0\\}$.
By Theorem 1.1, we have $c_{i}\not\equiv 0$ for at most one $i$. As $\deg
P\geq 2$, $b_{l}\neq 0$, therefore, $c_{1}\equiv 0$, i.e.,
$L(f^{\prime})+b_{1}f=a_{n}f^{(n)}+\dots+a_{1}f^{\prime}+b_{1}f=0$. This
implies that $\rho(f)\leq 1$ which is a contradiction to the order of growth
of $f$ is at least $2$.

## Acknowledgments

The authors would like to thank the anonymous referees for their valuable
suggestions and helpful comments.

## References

  * [1] C. Berenstein, D. C. Chang, and B. Q. Li.  A note on Wronskians and linear dependence of entire functions in $\mathbb{C}^{n}$.  Complex Variables, 24:131–144, 1993.
  * [2] W. Bergweiler, A. Eremenko, and J. K. Langley.  Real entire functions of infinite order and a conjecture of Wiman.  Geom. Funct. Anal., 13:975–991, 2003.
  * [3] W. Bergweiler, A. Eremenko, and J. K. Langley.  Zeros of differential polynomials in real meromorphic functions.  Proc. Edinburgh Math. Soc., 48:279–293, 2005.
  * [4] W. Bergweiler and W. H. J. Fuchs.  On the zeros of the second derivative of real entire functions.  J. Analysis, 1:73–79, 1993.
  * [5] W. Bergweiler, K. Ishizaki, and N. Yanagihara.  Meromorphic solutions of some functional equations.  Methods and Applications of Analysis, 5(3):248–258, 1998.
  * [6] J. Borcea and P. Brändén.  Pólya-Schur master theorems for circular domains and their boundaries.  Ann. of Math., 170(1):465–492, 2009.
  * [7] E. Borel.  Lecons sur les fonctions entieres.  Gauthier-Villars, Paris, 1921.
  * [8] P. Brändén, I. Krasikov, and B. Shapiro.  Elements of Pólya-Schur theory in the finite difference setting.  Proc. Amer. Math. Soc., 144(11):4831–4843, 2016.
  * [9] D. A. Brannan and J. G. Clunie, editors.  Aspects of contemporary complex analysis.  Academic Press, 1980.
  * [10] Y. M. Chiang and S. J. Feng.  On the Nevanlinna characteristic of $f(z+\eta)$ and difference equations in the complex plane.  The Ramanujan J., 16(1):105–129, 2008.
  * [11] T. Craven and G. Csordas.  Integral transforms and the Laguerre-Pólya class.  Complex Variables, 12(1-4):211–230, 1989.
  * [12] T. Craven and G. Csordas.  Problems and theorems in the theory of multiplier sequences.  Serdica Math J, 22(4):515–524, 1996.
  * [13] N. G. de Bruijn.  The roots of trigonometric integrals.  Duke Math. J, 17(3):197–226, 1950.
  * [14] A. Eremenko and L. A. Rubel.  On the zero sets of certain entire functions.  Proc. Amer. Math. Soc., 124(8):2401–2404, 1996.
  * [15] W. K. Hayman.  Meromorphic Functions.  Oxford Clarendon Press, 1964.
  * [16] S. Hellerstein and J. Williamson.  Derivatives of entire functions and a question of Pólya.  Trans. Amer. Math. Soc., 227:227–249, 1977.
  * [17] J. Huang and T. W. Ng.  Hypertranscendency of perturbations of hypertranscendental functions.  J. Math. Anal. Appl., 491(2):124390, 2020.
  * [18] O. Katkova, M. Tyaglov, and A. Vishnyakova.  Linear finite difference operators preserving the Laguerre–Pólya class.  Complex Var. Elliptic Equ., 63(11):1604–1619, 2017.
  * [19] O. Katkova, M. Tyaglov, and A. Vishnyakova.  Hermite-Poulain theorems for linear finite difference operators.  Constructive Approximation, 52:357–393, 2020.
  * [20] I. Laine.  Nevanlinna theory and complex differential equations.  Walter de Gruyter, Berlin, 1993.
  * [21] J. K. Langley.  Non-real zeros of linear differential polynomials.  Journal d’Analyse Mathématique, 107(1):107–140, 2009.
  * [22] J. K. Langley.  Postgraduate notes on complex analysis, 2014.
  * [23] B. Levin.  Distribution of zeros of entire functions, volume 5.  American Mathematical Society, 1964.
  * [24] B. J. Levin and I. V. Ostrovskii.  The dependence of the growth of an entire function on the distribution of zeros of its derivatives (Russian).  Sibirsk. Mat. Zh., 1:427–455, 1960.
  * [25] R. Nevanlinna.  Le théorème de picard-borel et la théorie des fonctions méromorphes. paris, gauthier-villars, 174 pp.  Jbuch, 55:773, 1929.
  * [26] T. W. Ng and C. C. Yang.  On the zeros of $\sum a_{i}\exp(g_{i})$.  Proc. Japan Acad. Ser. A Math. Sci., 73(7):137–139, 1997.
  * [27] N. Obreschkov.  Verteilung und Berechnung der Nullstellen reeller Polynome.  Berlin: VEB Deutscher Verlag der Wissenschaften, 1963.
  * [28] G. Pólya.  Über Annäherung durch polynome mit lauter reellen Wurzeln.  Rend. Circ. Mat. Palermo, 36:279–295, 1913.
  * [29] G. Pólya.  Bemerkung über die integraldarstellung der Riemannsche $\zeta$-funktion.  Acta Math., 48:305–317, 1926.
  * [30] G. Pólya and J. Schur.  Über zwei arten von faktorenfolgen in de theorie de algebraischen gleichungen.  J. Reine Angew. Math., 144:89–113, 1914.
  * [31] Q. I. Rahman and G. Schmeisser.  Analytic theory of polynomials.  Number 26. Oxford University Press, 2002.
  * [32] J. Rossi.  A halfplane version of a theorem of Borel.  In D. Drasin, I. Kra, C. J. Earle, A. Marden, and F. W. Gehring, editors, Holomorphic Functions and Moduli I, pages 111–118, New York, NY, 1988. Springer US.
  * [33] T. Sheil-Small.  On the zeros of the derivatives of real entire functions and Wiman’s conjecture.  Ann. of Math., 129:179–193, 1989.
  * [34] M. Tsuji.  On Borel’s directions of meromorphic functions of finite order, I.  Tohuku Math. J., 2:97–112, 1950.
  * [35] P. Walker.  Separation of zeros of translates of polynomials and entire functions.  J. Math. Ann. Appl., 206:270–279, 1997.