text:
\section{Introduction}
Radio pulsars are generally highly polarised sources, and the
propagation of their radiation through the magneto-ionised
interstellar medium (ISM) is affected by Faraday rotation. Measurement of
the rotation measure (RM), combined with the dispersion measure (DM)
for each pulsar, provides a good estimate of the average magnetic-field strength in the direction of the line of sight to many of the
known pulsars. With the high sensitivity and broad bandwidths of
current observing instruments, it is possible to accurately measure
RMs within a single band, thus avoiding calibration issues of
observing with different receivers in different bands. Furthermore,
observing pulsars with high temporal resolution allows RM measurements
at different pulse longitudes, provided the signal to noise ratio of
the linearly polarised component is sufficient.

Ramachandran et al. (2004) presented a phase-resolved RM analysis for
PSR B2016+28. In an attempt to explain the variation in RM by
$\pm15$~rad~m$^{-2}$ across the pulse profile, they conducted a single-pulse analysis of the polarization. They found that there are clearly
two, almost orthogonal polarization modes across the pulse profile,
with slightly different spectra (as also postulated in Karastergiou et
al. 2005 and demonstrated in Smits et al. 2006). Superposition of
these modes results in a net change in the polarization position angle
with frequency, leading to the measurement of a variable RM across the
pulse.  Although the effect measured in their study rotates the position
angle (PA) with frequency, they concluded that it is not Faraday rotation
intrinsic to the pulsar magnetosphere. The authors clearly demonstrate that for PSR B2016+28 the position angle does not have a quadratic dependence on wavelength: i.e.~$\Delta{\rm PA}={\rm RM}\cdot\lambda^2$.

A very recent work by \cite{kar09} showed that even small amounts
of interstellar scattering can significantly deform the steep features of PA profiles (e.g.~orthogonal mode jumps), and can provide a simple explanation for the RM variations observed by Ramachandran et al. The latter is caused by the frequency dependence of the scattering constant, which causes a frequency-dependent amount of polarised intensity to be transfered from a given phase bin to later phases. Regions where the PA profile has a steep gradient are worst affected, since it is at those phases where the scattered portion of the polarised intensity has the most impact on the shape of the PA profile: flat profiles can only stay flat, but the kinks and wiggles in PA profiles will be distorted by the redistribution of the polarised intensity.
Hence, Karastergiou noted that the polarization of core components, where PA profiles are steepest, is affected the most. The author concluded that scattering can smear orthogonal polarization modes so that they appear slightly less ``orthogonal''.

In this paper we concentrate on pulsars for which we have recently
published RM values (e.g. Noutsos et al. 2008). Our aim is not only to
fit for pulse phase--resolved RMs, as was done in Ramachandran et al., but to
closely investigate whether the polarization position angles obey
$\lambda^2$ fits. We investigate the possible causes for this
particular PA frequency dependence and discuss implications on RM
measurements performed using different techniques.

\section{Observations and data analysis}
The data used in our analysis were collected with the Parkes 64-m
radio telescope during two observing periods: in 2004, August 31 to
Septemeber 1, and in 2006, August 24 to 27. We observed at the
frequency band of the 20cm H-OH receiver, installed at the focus cabin
of the Parkes telescope. The receiver is equipped with a pair of
orthogonal, linear feeds that are sensitive to polarised emission with
an equivalent system flux density of 43 Jy at 20cm (\nocite{joh02}Johnston 2002). The 2004 data were
observed with a central frequency of 1.375 GHz and 256 MHz bandwidth,
while in 2006 we observed with a central frequency of 1.369 GHz and
256 MHz bandwidth.

The wide-band correlator (WBC) back end that was used in both
observations is capable of splitting the available bandwidth
into 1,024 frequency channels, on each of four polarizations, $I$, $Q$,
$U$ and $V$, and producing 10-s de-dispersed pulse subintegrations
with up to 2,048 phase bins across the pulse profile.

In Noutsos et al. (2008), we outline the method used to derive a
single RM value for each pulsar we observed. Here, we performed the
same analysis of the polarised emission in each longitude bin, across
the total pulse width. Our aim was to produce RM profiles as a
function of longitude for each observed pulsar. Prior to our analysis,
the raw data from the WBC were subjected to the
necessary process of RFI mitigation and polarization calibration. We
processed the calibrated data using the PSRCHIVE software package
\cite{hvm04}.

It is expected that the plane of linearly polarised emission will be
rotated by the magnetised ISM it traverses on its way to the Earth:
this is the well-known interstellar Faraday effect. The amount of
rotation, $\Delta {\rm PA}$, for a given pulsar should be constant for
a specific frequency but scales as $1/f^2$ across the observing
band. By definition, the rotation between frequency $f$ and infinite
frequency gives the PA at frequency $f$:

\begin{equation}
\label{eq:paeq2}
{\rm PA}=c^2{\rm RM}\frac{1}{f^2}
\end{equation}

Since the RM is a physical property of the ISM, it is expected to be
constant for a given ISM configuration between the pulsar and the
Earth. Its value depends on the magnitude and direction of the
magnetic field, $\boldsymbol{B}$, along the propagation path of
the polarised emission, $C$; the free-electron density, $n_e$, along $C$,
and the length of $C$; it is given by the path integral

\begin{equation}
\label{eq:rmeq1}
{\rm RM} \propto \int_{C}n_e\boldsymbol{B} \cdot d\boldsymbol{l}
\end{equation}
where $d\boldsymbol{l}$ is directed from the pulsar to the Earth.

Apart from the pulsar distance, none of the remaining physical
quantities, on which RM depends, are directly measurable. However,
given a finite observing bandwidth, one can derive RM by fitting the
PA across the band with the quadratic function of equation~\ref{eq:paeq2}. The
PA of each bin $i$ in each frequency channel $j$ was calculated as
${\rm PA}_{i,j}=0.5\arctan{(U_{i,j}/Q_{i,j})}$.

Having previously
developed the software capable of such quadratic fits (see
\nocite{njk+08} Noutsos et al. 2008 for details), we performed
bin-wise fitting across the longitude range inside the pulse and
derived a value of RM for each individual longitude bin. In order to
increase the s/n per frequency channel, $j$, before fitting across the
whole available bandwidth we averaged the frequency channels by a
factor 16. For a single averaging factor, as is the case in this
analysis, only the statistical error on each RM value can be
calculated, and the fitting routine gives the
best integer RM for each bin. Finally, we plotted pulse profiles of the total
flux density, the linear- and circular-polarization components, and
the PA and RM profiles across the pulse; we only plotted significant PA
values with ${\rm s/n} > 3.5$ and RM values with $\sigma_{\rm RM}<5$
rad m$^{-2}$.

\section{Results}
 \label{subsec:20cmdata}
The resulting 19 RM profiles are shown in Fig.~\ref{fig:largeRMvar1}
and Fig.~\ref{fig:smallRMvar1}. 9 pulsars show large RM variation
across the pulse and 10 pulsars show a rather constant RM profile. In
these plots, the bottom panel shows the polarization profile, where
solid lines denote total power; dashed lines, linear polarization; and
dotted lines, circular polarization. The middle panels show the
integrated PA profile and the top panels, the RM profile. In the interest of
characterisation, we have included a number of properties of the pulsar emission and the ISM in the direction of the 19 pulsars, in Table~\ref{tab:pulprop}.

\subsection{Pulsars with large RM variation}

{\bf PSR J0835$-$4510 (2004 \& 2006 data).} The Vela pulsar has a
highly polarised, yet relatively simple pulse profile in both the 2004
and 2006 data. Both profiles appear scattered.  The PA shows the
characteristic swing which has been used as a strong argument for the
rotating vector model (\nocite{rc69a}Radhakrishnan \& Cooke 1969). The most striking profile of Vela is its RM
profile, which shows a clear swing across the pulse that resembles a
sideways `S' (peak-to-peak variation $\approx 6$ rad m$^{-2}$). The
shape of the RM profile persists through both 2004 and 2006 data sets,
albeit the earlier of the two profiles displays a somewhat higher peak-to-peak variation of the RM ($\approx 13$ rad m$^{-2}$).

{\bf PSR J0942$-$5552 (2004 data).} The pulse profile of this pulsar
has 3 distinct components, with the peaks of linear polarization being
slightly offset from the corresponding ones in the total power. The PA
profile shows an orthogonal jump at the bridge emission between the
leading and central components and at the exact phase where $L$
vanishes. The RM profile shows a steep decline of high-s/n points
coincident with the leading end of the central component and an
equally steep increase, coincident with the leading component; this
behaviour occurs around the orthogonal PA jump. In both cases, the RM
changes by at least 20 rad m$^{-2}$. At the central peak's maximum the
RM fluctuates within a few rad m$^{-2}$.

{\bf PSR J1056$-$6258 (2004 data).}  The total-power profile of this
pulsar appears very slightly stretched toward later pulse phases. It
has nearly zero circular polarization but $\approx 50$\% linear. The
PA profile sweeps smoothly across the pulse, without showing any
breaks. The RM profile is roughly consistent with the published RM
around the the flux maximum, but it drops suddenly near the pulse edges,
giving it an overall shape of a swing across the entire pulse. \cite{kj06} performed a comparison between the 20cm and 10cm profiles of this pulsar and concluded that the linear polarization, and hence the PA profile, must undergo significant frequency evolution. They supported their argument by deriving a pulse-averaged RM value for PSR B1056$-$6258 that was much closer to the published one (6 rad m$^{-2}$), whereas if no evolution was assumed that value was $-1$ rad m$^{-2}$. However, this pulsar shows a peak-to-peak RM variation across the pulse of almost 100 rad m$^{-2}$ --- the highest in our sample --- which dwarfs any possible RM error due to profile evolution.

{\bf PSR J1243$-$6423 (2004 data).} The central pulse of this pulsar's
profile consists of two Gaussian components, whereas a third, much
weaker component is present at the far leading edge of the
profile. The circular polarization profile resembles an unsusual
double swing, with its minimum value being roughly coincident with the
maximum $L$. There is an orthogonal PA jump after the leading
component. The PA profile corresponding to the central pulse
components displays a full `S'-shaped swing, characterised by wiggles
often seen in PA profiles. The RM profile of this pulsar is extremely
complex: it shows a clear structure that departs from flatness, with a
leading swing followed by a sharp peak and a steep decline of the
RM. The peak-to-peak variation of the profile is $\approx 60$ rad
m$^{-2}$. The maximum RM ($\approx 200$ rad
m$^{-2}$) appears coincident with the inflexion point of the PA swing.

{\bf PSR J1326$-$5859 (2004 \& 2006 data).} This pulsar's profile is the sum of at least two components: a weak leading one and
a strong central one that leads into a trailing tail, perhaps
indicative of a third, very weak component. The linear polarization of
the central component is moderatly high, with more than 30\% of its
total flux being polarised, whereas the leading component is just over
20\% polarised. Even more interesting is the evident swing of the
circular polarization under the central component. Like in the case of
PSR J1157$-$6224, this may well be hinting at emission from near the
magnetic axis. In both epochs, the PA profile is broken by an
orthogonal jump near the peak of the leading component, at the exact
point where $L$ vanishes. At later phases, the PAs form an unusual
swing interrupted mid-way by a bulge that is coincident with the
maximum of the linear polarization. The RM profile clearly shows a
smooth functional dependence on phase that deviates from flatness;
like in the case of Vela, this dependence is exactly reproduced in both
2004 and 2006 data, which argues against this being an instrumental
effect. The form of the RM function resembles a double swing with 2
local minima separated by a local maximum. The trailing of the two
minima is exactly coincident with the peak in total power.

{\bf PSR J1453$-$6413 (2004 data).} The pulse profile of this
pulsar is a composition of a relatively weak leading component and a
much stronger central one. However, in the linear-polarization profile
these two components are of comparable intensity. The PA profile
swings along the entire central component but is disrupted by an
orthogonal jump that coincides with the merging point of the two
components' linear polarization, where $L=0$. The RMs across the pulse
show a moderate amount of scattering that prevents us from being
certain about the existence of a possible structure. On one hand,
there is a clear break in the RMs due to the jump in the PAs. On the
other hand, even away from the jump there seems to be a dip in the RMs
that is coincident with the peak of the central component's linear
polarization.

{\bf PSR J1644$-$4559 (2004 data).} The total power of this pulsar can
be described by a single Gaussian with an extended, trailing scattering
tail (see also analysis by \nocite{joh04}Johnston 2004). A hint at a second, very weak component can be seen near the
leading edge of the main pulse. The maximum of the linear polarization
occurs at slightly earlier phase than that of the total-power
maximum. In addition, this pulsar exhibits a swing in the circularly
polarised profile. The PA profile shows a more complex picture, in
comparison: it features two smooth bumps across the main pulse, which
are only interrupted by an approximately orthogonal jump, coincident
with the bridge emission between the main pulse and the weak leading
component; it also shows the obvious signs of scattering on
polarization shown in Li \& Han (2001)\nocite{lh03}. The RM profile
features a clear structure that comprises at least three local maxima
and an equal number of local minima; none of these extremes is exactly
coincident with either the total flux or the linearly polarised flux
extremes. The fluctuation of the RM across the pulse appears mildly
correlated to the PA profile.

{\bf PSR J1745$-$3040 (2004 data).} There are three, distinct, highly
polarised components in this pulsar's profile. \cite{jkm+08} performed multifrequency observations of this pulsar and noted the significant frequency evolution of its profile as well as the apparent increase of the linear polarization fraction with frequency. The PA profile
resembles an upside-down $\upsilon$ across the trailing and central
components, which is interrupted by an orthogonal jump occuring
between the trailing components and the leading one. The RM profile
shows a seemingly constant downwards slope for the most part of the
two trailing components but kinks downwards at the leading edge of the
central component. Above the leading component, the RMs are scattered,
although an upwards slope can perhaps be inferred from the 6 data points
coincident with this component.

{\bf PSR J2048$-$1616 (2006 data).} Three Gaussian components can be
discerned in the total profile of this pulsar. The linear polarization
is strong in the leading and trailing of the three components but
relatively weak in the central one. The PA profile resembles an
`S'-shaped swing that rotates the PAs smoothly across the profile (see also \nocite{jkk+07}Johnston et al. 2007; \nocite{jkm+08}Johnston et al. 2008). The
RM is rather constant across the leading component but dips by
$\approx 8$ rad m$^{-2}$ at the central, least-polarised component. A
sudden rise of the RMs is observed on the trailing edge of the
trailing component, although the s/n of the linear polarization is
already low in that phase region.

\subsection{Pulsars with low RM variation}

{\bf PSR J0134$-$2937 (2006 data).} The PA profile of this pulsar is
flat everywhere, apart from the phase region around the pulse maximum,
where the profile is broken by an orthogonal jump. The jump is also
coincident with the narrowest of the two distinct components present in
the linearly polarised flux, the other being broad and coincident with
the long leading tail of the total flux. The RM profile of this pulsar
shares the same flatness, which seems to be disturbed by the presence
of the orthogonal jump. However, the RM fluctuation occurs only over a
few phase bins, and its magnitude is of the same order as at least
three other instances away from the orthogonal jump.

{\bf PSR J0738$-$4042 (2004 \& 2006 data).}  The polarization of this
pulsar was shown at two frequencies (1.4 and 3~GHz) in Karastergiou \&
Johnston (2006). The position angle shows two orthogonal-mode
jumps. The RM profile, although generally flat, is below the average
RM value in the leading part of the profile; it then turns above the
average RM value for the rest of the profile, after the first
orthogonal jump. In the dual frequency data of Karastergiou \&
Johnston, the frequency evolution of the position-angle profile is in
agreement with what is seen here: the first segment of the PA profile
has a different frequency dependence to the rest of the profile.

{\bf PSR J0742$-$2822 (2004 \& 2006 data).} This pulsar comprises a
complex, 7-component profile (see e.g.~\nocite{kra94}Kramer
1994; \nocite{kjm05}Karastergiou et al. 2005). Despite the complexity of the pulse itself, the PA and RM
profiles display a relatively simple dependence on phase, with the
former being a simple sweep across the profile and the latter being
rather flat. It is perhaps noteworthy that a comparison
between the profiles of 2004 and 2006 shows that the average RM in the earlier observations appears slightly higher than the published value, whereas the opposite is observed in 2006. The $\langle {\rm RM}\rangle$ difference between the two data sets is $\approx 5$ rad m$^{-2}$.

{\bf PSR J0907$-$5157 (2006 data).} The total-power and linear-polarization profiles of this pulsar have two components. The PA profile
is a smooth, continuous sweep, without any breaks, and one can notice
that it slopes slightly upwards, toward the trailing edge of the
pulse profile. The RM profile is flat and consistent within the errors
with the published RM.

{\bf PSR J1057$-$5226 (2006 data).} This young pulsar ($5.4\times 10^5$ y) is an interesting and a peculiar one, and one of a few
that are visible at all wavelengths (see e.g.~\nocite{vl72}Vaughan \& Large 1972; \nocite{ch83}Cheng \& Helfand 1983;
\nocite{mcb97}Migniani et al.~1997; \nocite{fbc+94}Fichtel et al.~1994). Even more interesting is its strong radio interpulse, which
is separated from the main pulse by roughly 150$^\circ$ and has the same energy as the main pulse over a wide range of radio
frequencies (\nocite{mhak76}McCulloch et al.~1976; Weltevrede \& Wright 2008 {\em MNRAS submitted}); it has been suggested that we observe emission from both magnetic poles
(\nocite{lm88}Lyne \& Manchester 1988; \nocite{wqxl06}Wang et al.~2006). Furthermore, the steepness of the PA swing for the main
pulse and the interpulse differs significantly, which has given rise to the suggestion that the corresponding emission is generated
at different heights in the magnetosphere (\nocite{wqxl06}Wang et al.~2006). Due to the higher degree of linear polarization ($\sim 90\%$) and the modest s/n of the pulse in the 2006 data, we chose to study only the main pulse for RM variations. The PA profile of the main pulse is flat, with no breaks, and shows rotation of no more than 60$^{\circ}$ across the entire pulse. The RM profile is also consistent with a constant RM.

{\bf PSR J1157$-$6224 (2004 data).} This pulsar shows a simple total-power profile, with substantial amounts of linear ($>25$\%) and circular
polarization.
The PA profile exhibits a smooth but
complex swing, broken only by an orthogonal jump half-way along the
trailing edge of the pulse.  The RM profile also swings, although the
edges are too noisy and the magnitude of the effect too small to draw a
certain conclusion.

{\bf PSR J1253$-$5820 (2006 data).} This pulsar's pulse profile can be
fitted with a single Gaussian component, and the same pulse shape is
shared by the linearly polarised flux, which accounts for a large
fraction of the total observable emission. The PA shows a smooth
variation throughout. The RM profile shows little deviation from
flatness, mainly at the leading edge.

{\bf PSR J1359$-$6038 (2004 data).} This pulsar's profile is $>70$\%
linearly polarised and $\approx 15$\% circularly polarised. The PA profile
is rather flat across the pulse. The RM profile is consistent with a
constant value for the most part of the profile, diverging slightly at
the tail end, similarly to PSR J1056$-$6258. However, at that point the s/n is greatly reduced so solid conclusions cannot be extracted.

{\bf PSR J1705$-$1906 (2006 data).} This is another highly polarised
pulsar with a main pulse composed of multiple components: despite having been classified as a three-component pulse, at least five components can be discerned in both the total power and the linear polarization
profiles in our data (\nocite{ran90}Rankin 1990; \nocite{gl98}Gould \& Lyne 1998). In addition, a much simpler, single-peaked interpulse is also present in this pulsar's profile, which is separated by roughly $180^\circ$ from the main pulse. At 20cm, the interpulse is almost 100\% linearly polarised, which has also been noted in previous work (see e.g.~\nocite{kj04}Karastergiou \& Johnston 2004). In this work, we chose to examine only the main pulse for RM variations, due to its larger width and higher s/n. The PA profile is a gradual sweep across the entire main pulse. The RM profile appears flat, on average, although a coherent jitter over a-few-bin scale is evident. The only possible departure from flatness appears to coincide with the leading outrider.

{\bf PSR J1807$-$0847 (2004 data).} Another case of a multi-component,
total-power profile, with low linear polarization and with a central
component surrounded by at least three outriders. The PA profile is
broken twice by orthogonal jumps, which coincide with the minima of the
linear polarization; this has also been noted in a separate analysis by \cite{kj06}. The RM profile is quite noisy, owing to the low
linear polarization, but flat, on average, across the entire pulse.

\subsection{Fit Accuracy Across the Band}
In order to check for the presence of any frequency-dependent, systematic effect that could cause a change of the PA across the frequency band that is different to the $\lambda^2$ law of Faraday rotation, we examined in more detail the fits that produced the calculated RMs for Vela, J1243$-$6423 and J1326$-$5859.
Fig.~\ref{fig:0835rmfits}, \ref{fig:1243rmfits} and
\ref{fig:1326rmfits} show the RM fits for the minimum- and maximum-RM
bins (denoted by arrows) that reside near the central part of the
profile of the three aforementioned pulsars. It is important to note
that the associated errors on the fitted RMs, being only statistical,
would be underestimated if there was an additional, systematic,
non-Faraday component that affects the PAs. In order to reveal the
presence and magnitude of a possible systematic discrepancy between
the fitted function and the data, for each fit we also plotted the fit
residuals and their distribution (shown below the fit plots, in
Fig.~\ref{fig:0835rmfits}, \ref{fig:1243rmfits} and
\ref{fig:1326rmfits}). The quality of the data and the goodness of the
fit are evident in these figures. Upon closer inspection, however, the
residual PAs of all three pulsars, but mostly J1243$-$6423 and
J1326$-$5859, show evidence of some weak, sinusoidal fluctuations around the mean. It is difficult to attribute some significance to this, especially given
the Gaussian-like distributions of the residuals. The conclusion
therefore is that, if there is indeed a physical mechanism that affects the PA rotation, apart from the interstellar Faraday effect, we are unable to distinguish it from the dominant $1/f^2$ law.

\subsection{DM errors}
We also investigated the possibility that the RM variations observed are
caused by erroneous estimation of the DM, with, say, a $\Delta {\rm DM}$
offset from the true value. Such an offset is translated into a
dispersive delay between the pulse arrival time at a reference
frequency $f_0$ (usually the centre frequency of the band) and that at
frequency $f=f_0+\Delta f$, where $\mid\Delta f\mid \ll f_0$. In that
case, the delay in ms is given by

\begin{equation}
\label{eq:ddelay}
\Delta t = 8.3\times 10^6 {\rm ms} \ \ \frac{\Delta f}{f_0^3}\Delta {\rm DM}
\end{equation}
which also corresponds to a phase offset in the pulse window, $\Delta
\phi=\Delta t/P$, with $P$ the pulsar period.

If this phase offset is not accounted for, the calculation of the PA
in each frequency channel will correspond to a pulse longitude that is
shifted with respect to that at the reference frequency; this may
introduce a frequency-dependent shift in the PA of any given phase
bin.
In order to eliminate this as the cause of the observed RM
variations, we de-dispersed our data using the published DM value, as well as two other values that were $5\sigma$ away, either side of the published one, and generated again the RM profile in each case. The published DMs for our pulsars are shown in Table~\ref{tab:pulprop}.
By checking the resulting RM profiles for all pulsars, we concluded that there is no way of recovering a flat RM curve for any reasonable departure from the published DM. We therefore rule out DM errors as a possible cause of RM variation.

\section{Discussion}
Our analysis reveals that the RM profiles of pulsars with large RM
variations are diverse: they can be as simple as the `S'-shaped
profiles of Vela and PSR J1056$-$6258 or as complex as those of PSR
J1243$-$6423 and PSR J1644$-$4559. The peak-to-peak variation of those
pulsars' RMs across their profiles was found to be as little as a few
rad m$^{-2}$ (e.g. Vela in 2006) but also as large as 100 rad m$^{-2}$,
for PSR J1056$-$6258. Despite the hints of correlation between
individual features in the flux and PA profiles and those in the RM profiles (see e.g.~PSR J1326$-$5859), in most cases these correlations
seem coincidental. However, a general observation is that the greatest RM fluctuations across the pulse profile seem coincident with the steepest gradients of the PA profile, whereas pulsars with flat PA profiles show little or none at all RM variation. Good examples that support this argument are PSR J1243$-$6423 and PSR J2048$-$1616, whose RM departure from the mean is maximised around the inflexion point of the PA swing; in contrast, PSR J0907$-$5157 and PSR J1057$-$5226, having flat PA profiles, show insignificant RM variations.

The subset of pulsars with low RM variations adheres to the
theoretical predictions. More than half of those pulsars have
phase-resolved RMs that are scattered around their respective
published RM value, thus showing convincingly longitude-independent RM
values across the pulse. However, there are cases like those of PSR
J1157$-$6224, PSR J1253$-$5820, PSR J1359$-$6038 and PSR J1705$-$1906,
where the RM scatter across the pulse is clearly not
Gaussian: these profiles show a clear structure but with a small RMS
spread. An interesting case is that of PSR J0742$-$2822, where
although no significant RM variations exist across the profile, the
longitude-averaged RM across the pulse is higher in 2004 data than in
2006. Nevertheless, one has to be careful when drawing the above
conclusions; the main statement is that any RM variations with pulse
phase are smaller than the statistical error of the measurement.

We have demonstrated that the RM variations observed are not due to measurement errors, and that the PA rotation in all cases follows the
$\lambda^2$ law of Faraday Rotation. In the following we consider three possible explanations: (a) Faraday Rotation within the pulsar magnetosphere, (b) the superposition and frequency dependence of quasi-orthogonal
polarization modes, and (c) interstellar scattering and its effects on the polarization.

\subsection{Magnetoshperic Faraday Rotation}

Faraday Rotation through the ionised ISM has the effect of rotating
the PA with frequency. Our data, however, comprise of full sets of the
four Stokes parameters, $I$, $Q$, $U$ and $V$, so we can observe the frequency
dependence of the entire polarization vector
($\hat{\boldsymbol{P}}=(Q\hat{\boldsymbol{i}}+U\hat{\boldsymbol{j}}+V\hat{\boldsymbol{k}})/\sqrt{Q^2+U^2+V^2}$)
as a function of frequency, for specific high-s/n phase bins. An
example for one phase bin of PSR
J1326$-$5859 is shown in Fig.~\ref{fig:poincare}, where the tip of the Poincar\'e vector is plotted across
32 frequency channels. This plot demonstrates that, together with the
pronounced effect of PA rotation (around the vertical axis), a small
and rather systematic change in the fraction of circular polarization
can be seen. We obtained similar plots for four additional pulsars in
our sample that displayed large RM variations: namely Vela, PSR
J1243$-$6423, PSR J1644$-$4559 and PSR J2048$-$1616. In nearly all
cases, we found that the trace of $\hat{\boldsymbol{P}}$ on the
Poincar\'e sphere does not strictly follow paths of constant latitude,
as in Fig.~\ref{fig:poincare}.

The wave modes in the cold and ionised ISM are purely circularly
polarised. Therefore, the phase offset introduced by propagation
through this medium results in the well-documented changes of the
PA. On the contrary, in a medium of highly relativistic charged
particles, the natural modes are close to linearly polarised (\nocite{saz69}Sazonov 1969; \nocite{mel97}Melrose 1997a). As a
consequence, different phase offsets between the modes will result in
a conversion between linear and circular polarization. In a pulsar's
magnetosphere, the observed degrees of circular polarization indicate
that the natural wave modes are elliptically polarised, which suggests
that propagation through this medium will generate some frequency-dependent conversion between linear and circular polarization. An
extensive discussion on this can be found in \cite{km98}.

If such a generalised Faraday Rotation is taking place within the
pulsar magnetosphere and giving rise to RM variations as well as
changing the degree of circular polarization in the observed
frequency band, we expect the greatest change in $V$ to coincide with
the greatest variations in RM.
In order to test this, we first calculated the fraction of circularly
polarised intensity, $V/I$, for each frequency channel across the
band. The division by $I$ removes any fluctuations of $V$ caused by
scintillation through the interstellar medium, whereas any change in $V$
that originates in the magnetoshpere should remain
unaffected. Subsequently, we fitted a linear function, $V/I \propto
f$, to the data and calculated the change in $V/I$ across the band,
based on that function. This was done to get a good estimate of the
total effect while minimizing the influence of instrumental noise. We
found linear fits to describe the data sufficiently well.
The resulting profiles of $\Delta(V/I)$ as a function of longitude,
for Vela, PSR J1243$-$6423, PSR J1326$-$5859, PSR J1644$-$4559 and PSR
J2048$-$1616 are shown in Fig.~\ref{fig:VDVprofiles}.  A side-by-side
comparison of the RM and $\Delta(V/I)$ profiles, shown in
Fig.~\ref{fig:VDVprofiles}, makes clear that none of the pairs of profiles seems correlated.

There is certainly little evidence that the maximum RM deviation from the mean and the maximum change in circular polarization across the band occur at the same pulse phase. Therefore, to first order, the assumption that the observed RM variations are caused by magnetospheric Faraday rotation is deemed unlikely.

\subsection{Quasi-orthogonal Emission Modes}
Another factor that could potentially affect the measured RM is the
{\em incoherent}\/ superposition of quasi-orthogonal modes of
polarization (see e.g.~\nocite{rbr+04}Ramachandran et al. 2004). The superposition of strictly
orthogonal modes results in an observed PA of the strongest of the two
modes, regardless of the ratio of the intensity of the modes
(e.g.~\nocite{brc76}Backer et al 1976; \nocite{crb78}Cordes et al. 1978). The reason is that orthogonal modes are represented by
anti-parallel vectors on the Poincare sphere; regardless of their
relative lengths, the addition is always pointing in the direction of
either of the two modes. This is not the case for non-orthogonal
modes: in this case, the sum of the mode vectors points at an angle
that depends on the degree of non-orthogonality and the relative
intensity of the modes. As a consequence, if the ratio of the mode
intensities changes as a function of frequency, the observed PA will
also rotate. The question is under what conditions this could
resemble Faraday Rotation.

In the pulsars where significant RM variations are seen, the
variations often exceed 20 rad m$^{-2}$ across the pulse profile. At
the frequency and bandwidth of these observations, this corresponds to
a variation in PA by $\sim 18.5^\circ$. Therefore, in order for
quasi-orthogonal polarization modes to present a possible explanation
for the observed RM variations, they should be capable of inducing
such a PA rotation.

From polarization studies of 18 pulsars in the European Pulsar Network
(EPN) archive, \cite{sse+06} conclude that the ratio of the average
intensity of the orthogonal modes for those pulsars follows a power
law with spectral index that does not typically exceed
$\gamma \sim -0.5$. Assuming linear modes ($L_1$ and $L_2$; the
degree of circular polarization is generally very low), then at a
reference frequency $f_0$ the mode intensity ratio is
$R_0=L_1/L_2$, and the form of this power law is
\begin{equation}
\label{eq:opmpl}
R=R_0\left(\frac{f}{f_0}\right)^{-0.5}
\end{equation}

For given values of $R_0$, we can compute $R$ across the frequency band
we are interested in: in this case, between 1300 and 1556 MHz. We have
numerically computed the total rotation of the PA across the band (hereafter
$\Delta$PA), for different values of $R_0$ and different angles of
non-orthogonality, $\phi$. The latter angle is defined on the Poincar\`e
sphere as the supplementary to the angle between the two modes. Figure
\ref{fig:dpafunc} shows $\Delta$PA versus $\phi$, for various values of $R_0$,
which equate to certain degrees of linear polarization. The plot on
the left shows examples where $R_0\geq 1$ as opposed to $R_0\leq 1$
for the plot on the right. Both plots show a set of lines
corresponding to the same degrees of polarization , $L$\%, which is
approximately $\mid L_1-L_2\mid/I$.  Two things are immediately
evident: first, the greatest rotation of PA happens for small values
of $\phi$; second, and most important, the only possibility that
$\Delta$PA$>15^\circ$ occurs when the fractional polarization remains under
10\% and only in the specific case where $R_0>1$, i.e. the case where
the mode with the steepest spectral index is brightest at the lowest
frequency $f_0$.  In our data, there are clear examples of high $\Delta$RM
and high $L$\%, the most prominent being Vela and J1326$-$5759.  It is
therefore unlikely that non-orthogonal polarization modes provide a
general explanation for the RM variations in the data.

\subsection{Scattering}
As was mentioned earlier in this paper, some of the pulsars we studied
show signs of interstellar scattering in their flux profiles: vivid
examples are PSR J1056$-$6258, justified by its high DM of 320.3 pc
cm$^{-3}$, and PSR J1644$-$4559, with a DM of 478.8 pc
cm$^{-3}$ --- the highest in
our sample. We checked for a possible correlation between the amount
of interstellar scattering and the magnitude of the observed RM
variations. The effect of interstellar scattering is usually
quantified by the scattering timescale, $\tau_{\rm sc}$. Assuming
a Kolmogorov spectrum for the density fluctuations along the
sightline to the pulsar, $\tau_{\rm sc}$ has a roughly quadratic
relation to DM for low-DM pulsars (i.e.~$\tau_{\rm sc}\propto {\rm DM}^{2}$)
but can be much steeper for high-DM pulsars: i.e.~$\tau_{\rm sc}\propto {\rm DM}^{4}$ (e.g.~\nocite{sut71}Sutton 1971; \nocite{ric77}Ricket 1977). More recently, based on a global fit to 98 scattering times of high-DM pulsars, \cite{bcc+04} showed that a parabolic curve of the form

\begin{equation}
\label{eq:tausc}
\log\tau_{\rm sc}=-6.46+0.154\log{\rm DM}+1.07(\log{\rm DM})^2-3.86\log{f_0}
\end{equation}
where $f_0$ is the observing frequency, is a reasonable description of the trend of $\tau_{\rm sc}$ as a function of DM. Hence, as long as we are refering to the same observing frequency, and since there are only a few available measurements of $\tau_{\rm sc}$ for our pulsar sample (see Table~\ref{tab:pulprop}), the DM is a convenient alternative to the amount of scattering a pulsar's emission is subject to.

We have quantified the magnitude of the RM variations for those
pulsars whose RM profiles show signs of a longitude-dependent
structure, in terms of the peak-to-peak distance of their RM profile,
$\Delta {\rm RM}_{\rm pp}$. Fig.~\ref{fig:dRMmaxDM} shows a scatter
plot of $\Delta {\rm RM}_{\rm pp}$ against DM for 9 pulsars. Although
there is perhaps no clear correlation, it is striking that the pulsars
that typify high RM variations are also high-DM pulsars. This includes
the most prominent examples of RM variation with pulse phase, such as
PSR J1243$-$6423, PSR J1326$-$5859 and PSR J1644$-$4559. In a recent
paper, \cite{kar09} showed that small amounts of scattering and
orthogonal polarization modes of emission in the pulsar magnetosphere
can result in a frequency dependent PA that can be falsely attributed
to Faraday Rotation. This explanation is extremely simple and may well
be the key to understanding the high-DM examples. It is also
strengthened by the observation that most of the low-DM pulsars show
relatively smooth PA profiles, devoid of the bumps and wiggles seen in
the higher DM profiles. The overall correlation in
Fig.~\ref{fig:dRMmaxDM} is perhaps reduced by the fact that RM
measurements can only be made for pulse phases with significant linear
polarization.

\section{Conclusions}

We have shown phase-resolved RM profiles for 19 pulsars, 9 of which show
substantial RM variations across the pulse profile. These variations appear in
most cases to show coherent rather than noisy behaviour. The high s/n of our
data and our robust fitting technique have shown that whatever the cause of the
PA rotation with frequency, it is inseparable from the $1/f^2$ law of
interstellar Faraday Rotation within our observed band.

We have considered a
number of possibilities that may result in the observed phenomena. For example,
we have ruled out the possibility that de-dispersion with a slightly erroneous
DM may be responsible for the data here. We have discussed the possibility that
the RM variations are of magnetospheric origin. Two cases were examined, that of
generalised Faraday Rotation due to the relativistic plasma in the pulsar
magnetospheres, and the existence and superposition of quasi-orthogonal modes of
polarization.

In the case of generalised Faraday Rotation, we examined whether
pulse phases with the highest rotation of the PA also showed the highest change
of the fractional circular polarization. Our first-order investigation turned
out negative results, although it is true that the theory for generalised
Faraday Rotation is not complete and does not provide us with good handles on
the rate at which polarization rotations occur with frequency.

As concerns
quasi-orthogonal modes, we find that it is impossible to obtain a large rotation
of the PA within our observing frequency band and observe large degrees of
linear polarization at the same time. Our analysis indicates that there are
specific cases of the data that can be explained by this process, but it
certainly does not form a general framework for explaining RM variations across
the pulse.

Finally, this work supports the ideas presented in \cite{kar09}, that interstellar scattering can distort the polarization in a frequency-dependent way. This may resemble Faraday Rotation and will generally be
pulse-phase dependent. It is important to emphasise that the possibility of
magnetospheric Faraday Rotation is not totally eliminated by this work; we
merely conclude that it is an unlikely explanation for RM variations of the
observed magnitude.

There is a possibility that both the magnetospheric mechanisms and
interstellar scattering are responsible for the observed RM
variations, along with other possibe mechanisms which we have not
considered. Whichever the explanation, certain care must be taken
during RM-determination experiments from pulsars to correctly
identify the amount of PA rotation that is due to the iono-magnetised
interstellar medium. Especially for the studies of the planar component
of the Galactic magnetic field, where pulsars with high DMs and RMs are used,
the identification of the ISM-induced Faraday rotation and the correct
calculation of the associated RM plays an important role in our efforts to model the interstellar magnetic fields. Hence, in light of the RM profiles presented herein, pulsar-RM determination from specific, high-s/n pulse bins can lead to considerable errors.

\vspace{0.1cm}

\section*{Acknowledgements}
\noindent The Australia Telescope is funded by the Commonwealth of
Australia for operation as a National Facility managed by the CSIRO.
We would like to thank Andrew Lyne, Don Melrose and Roy Smits for useful discussions and comments on this manuscript.


text:
\section{Introduction}
\label{sec:intro}

Active inference (AIF) is prominent in neuroscientific literature as a general mathematical framework of the brain at the computational level~\cite{friston2010free}. According to this theory, the brain learns by interaction a generative model of the world/body that is used to perform state estimation (perception) as well as to execute actions, optimizing one single objective: Bayesian model evidence. This approach, which grounds on hierarchical variational inference and dynamical systems estimation \cite{friston2008dem}, has strong connections with Bayesian filtering~\cite{sarkka2013bayesian} and control as inference~\cite{millidge2020relationship}, as it both estimates the system state and computes the control commands as a result of the inference process.

In the last years, some proof-of-concept studies in robotics have shown that AIF may be a powerful framework to address key challenges in robotics, such as adaptation, robustness and abstraction in goal-directed tasks \cite{lanillos2018adaptive, oliver2021empirical, sancaktar2020end, pezzato2020novel, pio2016active, minju2019goal, baioumy2020active, meera2020free, lanillos2020robot}. For instance, a state estimation algorithm and an AIF-based reaching controller for humanoid robots were proposed in \cite{lanillos2018adaptive} and \cite{oliver2021empirical} respectively, showing robust sensory fusion (visual, proprioceptive and tactile) and adaptability to unexpected sensory changes. However, they could only handle low-dimensional inputs. Meanwhile, AIF-based joint torque control was investigated using proprioceptive inputs \cite{pezzato2020novel}, showing better performance than a state-of-the-art model reference adaptive controller. Latterly, we presented a pixel-based deep AIF controller~\cite{sancaktar2020end} that incorporated generative model learning using convolutional neural networks but limited to one sensor modality (visual) and velocity control. Thus, this work investigates how to design a brain-inspired controller, with the adaptability and robustness properties of AIF, that: 1) scales to high-dimensional inputs, 2) allows multisensory integration and 3) commands in torque.

\subsection{Contribution}
We propose a Multimodal Vartiational Autoencoder Active Inference (MVAE-AIF) torque controller that combines free energy optimization \cite{friston2010action} with generative model learning (Fig. \ref{fig:architecture}).

\begin{figure}[hbtp!]
	\centering
	\subfloat[Architecture]{
		\centering
		\includegraphics[width=0.8\columnwidth, height=150px]{Images/AIFscheme}
		\label{fig:architecture:def}
	}\\
	\subfloat[Camera input $\v{I}$]{
	        \includegraphics[width=0.38\columnwidth]{robot_arm_p}
	  \label{fig:architecture:visual}
    }
    \subfloat[Predicted input $\g_v(\latent)$]{
        \includegraphics[width=0.38\columnwidth]{robot_arm_rp}
         \label{fig:architecture:visual_out}
    }
	\caption{MVAE-AIF architecture. (a) Abstract description of the algorithm and environmental setup. State representation learning is given by the multimodal VAE, and estimation and control are provided by the AIF framework. The Panda robot had access to (noisy) proprioceptive input $\joints$ and a camera placed in front $\v{I}$. (b) Camera input. (c) Predicted visual input from the generative model.}
	\label{fig:architecture}
\end{figure}

We extended our previous AIF formulation \cite{oliver2021empirical, sancaktar2020end} to work with high-dimensional multimodal input at the torque level. Inspired by \cite{Yuge2019Variational}, \cite{sancaktar2020end} and \cite{rood2020deep}, we designed a VAE that provides light coupled multimodal state representation learning \cite{lesort2018state}. One of the main advantages of our approach is that the robot only has to learn the kinematic forward mapping and then is the inference process that produces the right torques for achieving the goal, even in the presence of unmodeled situations or external forces.

We studied, on a simulated\footnote{Planned onsite experiments with the real platform were not possible to finish due to university access restrictions.} 7DOF Franka Emika Panda industrial arm,
\begin{enumerate}
    \item How multimodal encoding improves representation power and thus, state estimation and control accuracy;
    \item How the variational inference nature of our approach provides strong adaptation and robustness against unmodeled dynamics, environment parameter variations and sensory noise.
\end{enumerate}
Results show its improved performance by qualitative and statistical comparison with two baselines: proprioceptive AIF \cite{pezzato2020novel} and the Panda industrial arm built-in controller. For reproducibility, the code is publicly available at~\url{https://github.com/Cmeo97/MAIF}.

\section{AIF general formulation and notation}
\label{sec:formulation}
Here we introduce the standard equations and concepts from the AIF literature \cite{friston2010free}, and the notation used in this paper, framed for estimation and control of robotic systems~\cite{oliver2021empirical}. The aim of the robot is to infer its state (unobserved variable) by means of noisy sensory inputs (observed). For that purpose, it can refine its state using the measurements or perform actions to fit the observed world to its internal model. This is dually computed by optimizing the variational free energy, a bound on the Bayesian model evidence.

\begin{description}
\item[System variables.] State, observations, actions and their n-order time derivatives (generalized coordinates).
\begin{align}
    \Tilde{\obs} &= [\obs,\obs',\obs'',...,\obs^n] &\text{, sensors}\nonumber\\
    \Tilde{\latent} &= [\latent,\latent',\latent'',...,\latent^n] &\text{, multimodal system state}\nonumber\\
    \Tilde{\muvec} &= [\muvec,\muvec',\muvec'',...,\muvec^n] &\text{, proprioceptive state}\nonumber\\
    \Tilde{\xnoise} &= [\xnoise,\xnoise',\xnoise'',...,\xnoise^n] &\text{, sensory noise}\nonumber\\
    \Tilde{\znoise} &= [\znoise,\znoise',\znoise'',...,\znoise^n] &\text{, state fluctuations}\nonumber\\
    \action &= \{a_1,\ldots, a_m\}&\text{, action (m actuators)}\nonumber
\end{align}
Where $\obs' = d\obs/dt$. Depending on the formulation the action $\action$ can be force, torque, acceleration or velocity. In this work action refers to torque. We further define the time-derivative of the state vector $D\Tilde{\latent}$ as:
\begin{align}
    D \Tilde{\latent} = \frac{d}{dt}(\latent,\latent',\ldots,\latent^n) = [\latent',\latent'',\ldots,\latent^{n+1}]\nonumber
\end{align}

\item[Generative models.]
Two generative models govern the robot; the mapping function between the robot's state and the sensory input $\g(\Tilde{\latent})$ (e.g., forward kinematics) and the dynamics of the internal state $f(\Tilde{\latent})$.
\begin{align}
    \Tilde{\obs} &= \g(\Tilde{\latent}) + \Tilde{\xnoise} \\
    D\Tilde{\latent} &= \f(\Tilde{\latent}) + \Tilde{\znoise}
    \label{gen}
\end{align}
where $\xnoise \sim \mathcal{N}(\boldsymbol{0}, \Sigma_x )$ and $\znoise \sim \mathcal{N}(\boldsymbol{0}, \Sigma_z )$ are the sensory and process noise respectively.  $\Sigma_x$ and $\Sigma_z$ are the covariance matrices that represent the controller's confidence about each sensory input and about its dynamics respectively.
\vspace{0.2 cm}
\item[Variational Free Energy (VFE).]
The VFE is the optimization objective for both estimation and control. We use the definition of the $\mathcal{F}$ based on \cite{friston2010action}, where the action is implicit within the observation model $\obs(a)$. Using the KL-divergence the VFE is:
\begin{align}
\mathcal{F} &= \text{KL}\left[ q(\latent) || p(\latent|\obs) \right] - \log p(\obs) \\
&= \text{KL}\left[ q(\latent) || p(\latent,\obs) \right] = -\text{ELBO} \nonumber
\end{align}

\noindent \textit{State estimation} using gradient optimization:
\begin{equation}
   \boldsymbol{\dot{\Tilde{\latent}}} =  D\Tilde{\latent} - k_{z}\nabla_{\latent}\mathcal{F}({\Tilde{\obs}}, {\Tilde{\latent}})
    \label{eq:perception_general}
\end{equation}
\noindent \textit{Control} using gradient optimization:
\begin{equation}
    \dot{\action} = -k_a \sum_x \frac{d \Tilde{\obs}}{da} \cdot \nabla_{\Tilde{\obs}}\mathcal{F}({\Tilde{\obs}}, {\Tilde{\latent}})
    \label{eq:action_general}
\end{equation}
The VFE has a closed form under the Laplace and Mean-field approximations \cite{friston2007variational,oliver2021empirical,buckley2017free} and it is defined as:
\small
\begin{align}
    \mathcal{F}(\Tilde{\latent},\Tilde{\obs}) \simeq& -\ln p(\Tilde{\latent},\Tilde{\obs}) = -\ln p(\Tilde{\obs}|\Tilde{\latent}) - \ln p(\Tilde{\latent})\\
    \simeq& \;(\obs-g(\obs))^T \Sigma_x^{-1} (\obs-g(\obs)) \nonumber\\
    &+ (D\latent-f(\latent))^T \Sigma_z^{-1}(D\latent-f(\latent)) \nonumber\\
    &+ \frac{1}{2} \ln|\Sigma_x| + \frac{1}{2} \ln|\Sigma_z| \label{eq:flaplace}
\end{align}
\normalsize
The first term of Eq. (\ref{eq:flaplace}) is the sensor prediction error and the second term is the dynamics prediction error.

\end{description}

\section{MVAE-AIF controller}
\label{sec:model}
The main novel addition presented in this work is the introduction of high-dimensional representation learning within the AIF optimization framework to provide at the same time multimodal state estimation and torque control. We first describe the architecture, the generative models and the training. Second, we detail the estimation and control equations and finally, we summarize the algorithm.

\subsection{Architecture and design}
Our proposed multimodal variational autoencoder active inference (MVAE-AIF) architecture is depicted in Fig.~\ref{fig:architecture:def}. Figures \ref{fig:architecture:visual} and \ref{fig:architecture:visual_out} shows the camera input and the predicted sensory input from the learnt generative model.

\subsubsection{Sensors}
We assume that the robot has (noisy) proprioceptive sensors that provide the joint angles $\joints$, velocities $\dot{\joints}$ and accelerations $\ddot{\joints}$. Furthermore, there is a camera that provides images $\v{I}$ of size $(w\times h)$ of the robot.

\subsubsection{Robot generative models}

The generative models of the agent are approximations of the real generative models of the world and are composed of a sensory generative function $g(\latent)$, which maps the internal state to the visual and proprioceptive spaces and a dynamics function $f(\latent)$, which describes the evolution the system.
\begin{align}
\obs_v &=   \g_{v} (\latent) + \xnoise_v &\text{\quad visual}\\
\obs_q &= \g_{q} (\latent) + \xnoise_q  &\text{\quad joints}\\
D\latent &= \f(\latent, \rhovec) + \znoise_{z}  &\text{\quad state dynamics}
\end{align}
where $\obs_v$ and $\obs_q$ are the visual and proprioceptive predictions respectively. The processes noise $\xnoise_q$, $\xnoise_v$ and $\znoise_z$ are assumed to be drawn from a multivariate normal distribution with zero mean and covariance matrices $\Sigma_{q}$,  $\Sigma_v$ and $\Sigma_z$ respectively. Note that we included $\rhovec$ into the dynamics function. This variable will describe the desired goal and acts as a prior that drives the system.

\subsubsection{Generative models and state representation}

The defined sensory generative models are approximated by a multimodal variational autoencoder (MVAE) which maps the internal state representation to both proprioceptive and visual sensory spaces into a common latent space. We designed the MVAE using two couples of parallel encoders and decoders. The MVAE latent space describes the state of the system and the output of the decoders are the predicted image $\obs_v$ and joint values $\obs_q$. The state of the system $z$ is defined as:
\begin{equation}
    \latent = encoder_{\joints}(\joints) + encoder_{v}(\v{I})
    \label{z_eq}
\end{equation}
As a result, modalities have a low correlation, allowing sensory integration without reciprocal degradation. In other words, if one of the two modalities is not well reconstructed, the other one can still be inferred. For instance, in the case of visual occlusion, proprioceptive reconstruction is still possible.
The robot can infer the sensory input given the state using the decoders, which act as the generative learnt models:
\begin{align}
    \obs_v = g_v(\latent) = decoder_{v}(\latent) \\
    \obs_q = g_p(\latent) = decoder_{q}(\latent)
\end{align}

As we only need the sensory generative models, the encoders are just employed for the training. Then, in the AIF control loop, only the decoders are used. The detailed description of the network layers can be found in Table \ref{tab:network_specs}.

\noindent For algorithmic purposes we also define the proprioceptive state (joint angles belief) as $\muvec= \{g_q(\latent), \muvec', \muvec'', \muvec'''\}$, where the first component is the predicted joint angles and the others are the higher-order derivatives.

\subsubsection{Generative models learning}
The MVAE was trained with a dataset containing 500000 images with size $(128 \times 128)$ and the associated joints angle values which we created performing motor babbling on the robot arm. To accelerate the training, we included a precision mask $\Pi_{\obs_v}$, computed by the variance of all images and highlighting the pixels with more information. The augmented reconstruction loss employed was:
\begin{align}
    \mathcal{L} = \text{MSE}(\obs_v \!+\! \Pi_{\obs_v}\obs_v,\; \v{I}) + \text{MSE}(\obs_q, \joints)
\end{align}

\subsubsection{Goal definition}
We define the goal of the system following a similar approach as \cite{oliver2021empirical}. The desired goal, described in the sensory space $\rhovec=\obs_d$, is introduced in the internal dynamics as follows:
\begin{align}
f(\latent,\rhovec\!=\!\obs_d) = T(\latent) (\obs_d -g(\latent)) = \frac{\partial \g(\latent)}{\partial \latent} (\obs_d - g(\latent))
\end{align}
Where $T(\latent)$ is a function that maps the error in the sensory space to the latent space. In this case, the partial derivative of $g(\obs)$ with respect to $\latent$. The desired goal in this work is defined by the final image and joints angles $\obs_d=\{\v{I}_d, \joints_d\}$.

\subsection{Estimation and control}
The online estimation and control problem is solved by optimizing the VFE through gradient optimization, computing Equations (\ref{eq:perception_general}) and (\ref{eq:action_general}).

\vspace{2mm}

State \textbf{estimation} uses the Laplace approximation of $\mathcal{F}$, described in Eq. (\ref{eq:flaplace}). We update the multimodal state $\latent$ by means of the partial derivative of the VFE with respect to $\latent$ ($\dot{\latent} = -k\partial_\latent\mathcal{F}$):
\small
\begin{align}
    \dot{\latent} =&  k_{v}\frac{\partial \g_{I}}{\partial z} \Sigma_v^{-1} (\obs_v - \g_{I}(\latent)) +  k_{q}\frac{\partial \g_{q}}{\partial z} \Sigma_q^{-1} (\muvec - \g_{q}(\latent)) \nonumber\\
    &+ k_{v}\frac{\partial \f_{I}}{\partial z} \Sigma_v^{-1} (\boldsymbol{I}_d - \f_{I}(\latent,\rhovec)) +  k_{q}\frac{\partial \f_{q}}{\partial z} \Sigma_q^{-1} (\joints_d - \f_{q}(\latent, \rhovec)) \nonumber\\
    \label{eq:up_z}
\end{align}
\normalsize
As we do not have access to the high order generalized coordinates of the latent space $\latent',\latent''$, we track both the multimodal shared latent space $\latent$ and the higher orders of the proprioceptive (joints) state $\muvec',\muvec''$. This is very useful as we can use the angular velocity and acceleration within the controller. Thus, we update the proprioceptive state velocity and acceleration by computing the partial derivative of $\mathcal{F}$ with respect to $\muvec',\muvec''$, as it was a unimodal proprioceptive controller~\cite{pezzato2020novel}, but the joint angles are predicted by the MVAE.
\small
\begin{align}
     \muvec =& g_q(\latent) \label{eq:update:mu}\\
     \dot{\muvec}' =& \muvec'' + k_{\mu} \left(\Sigma_{\dot{q}}^{-1} ( \dot{\joints} - \muvec') \nonumber\right.\\
     &\left.- \Sigma_{\mu}^{-1} ( \muvec' + \muvec - \joints_{d}) - \Sigma_{\dot{\mu}}^{-1} (\muvec'' + \muvec') \right) \label{eq:update:dotmu}\\
     \dot{\muvec}'' =&  - k_{\mu} \Sigma_{\dot{\mu}}^{-1} (\muvec'' + \muvec') \label{eq:update:ddotmu}
\end{align}
\normalsize
where $\Sigma_{\mu}^{-1},\Sigma_{\dot{\mu}}^{-1}$ are the precision (inverse variance) matrices related to joint angles and velocities beliefs.

\vspace{2mm}

The \textbf{action} (torque) is computed by optimizing the VFE using Eq. (\ref{eq:action_general}). Here we only consider the proprioceptive errors. In the AIF literature this models the reflex arc, where top-down proprioceptive predictions coming from higher cortical levels produce the motoric reflexes~\cite{otto2020deep}. Thus, the torque commands are updated with the following differential equation:
\small
\begin{align}
    \dot{\action} &= - k_a \left(\partial_\action \muvec \Sigma_q^{-1} (\joints - \muvec) + \partial_\action \muvec' \Sigma_{\dot{\mu}}^{-1} (\dot{\joints} - \muvec') \right) \\
    &= - k_a (\Sigma_{q}^{-1}(\joints - \muvec) + \Sigma_{\dot{\mu}}^{-1} (\dot{\joints} - \muvec'))
\end{align}
\normalsize
Although we can compute the action inverse models $\partial_\action \muvec, \partial_\action \muvec'$ through online learning using regressors~\cite{lanillos2018active}, we let the adaptive controller absorb the non-linearities of the relation between the torque and the joint velocity. Thus, as described by \cite{pezzato2020novel} we just consider the sign of the derivatives.

\subsection{Algorithm}
Algorithm \ref{alg:mvae-aic} describes the proposed controller. The first stage solves the free energy optimization using the prediction errors between the generative model and the sensor measurements. We compute $\dot{z}_{q}, \dot{z}_{v}, \dot{z}_{q_d}$ and $\dot{z}_{v_d}$ as the four terms in Eq. \ref{eq:up_z}. Since these terms are equivalent to backpropagation operations over the respective errors, as we proposed in \cite{sancaktar2020end}, they are computed by defining a node with $\dot{z}$ on the network virtual tree and backpropagating errors through the decoders. Once we have the contributions from the sensory and dynamics prediction errors, we update the robot multimodal state and the proprioceptive higher-order derivatives. Finally, we compute the action to minimize the inverse variance weighted proprioceptive prediction errors.

\begin{algorithm}[hbtp!]
\caption{MVAE-AIF}
\label{alg:mvae-aic}
\begin{algorithmic}
\Require $\obs_{desired}=\{\v{I}_d,\joints_d\}$
\While {$\neg goal$}
\vspace{0.1 cm}
\State $\joints,\v{I} \gets$ Sensors (joint values, camera)
\vspace{0.2 cm}
\State $Compute \hspace{0.1 cm} prediction \hspace{0.1 cm} errors \hspace{0.1 cm} in \hspace{0.1 cm} \latent \hspace{0.1 cm} space, \hspace{0.1 cm}Eq. \hspace{0.1 cm}\ref{eq:up_z}$
\State $\obs_{v} = \g_{v}(\latent)   \quad\quad Decoder \hspace{0.1 cm} visual $
\State $\obs_{q} = \g_{q}(\latent)  \quad\quad Decoder \hspace{0.1 cm} joints$
\State $\dot{\latent}_{q} = \g_{q}.backward(\Sigma^{-1}_{q}(\joints- \obs_q))$
\State $\dot{\latent}_{v} = \g_{v}.backward(\Sigma^{-1}_{v}(\v{I} - \obs_v))$
\State $\dot{\latent}_{q_d} = \g_{q}.backward(\Sigma^{-1}_{q}(\joints_{d} - \obs_q))$
\State $\dot{\latent}_{v_d} = \g_{v}.backward(\Sigma^{-1}_{v}(\v{I}_{d} - \obs_v))$
\vspace{0.2cm}
\State $Update \hspace{0.1 cm} state $
\State $\dot{\latent} = z_q+z_v+z_{q_d}+z_{v_d}$
\State $\dot{\muvec}', \dot{\muvec}''  \gets  proprioceptive \hspace{0.1 cm} state \hspace{0.1 cm} Eq. \hspace{0.1 cm} \ref{eq:update:dotmu}, \hspace{0.1 cm}  \ref{eq:update:ddotmu}$
\vspace{0.2 cm}
\State $Active \hspace{0.1 cm} inference$
\State $\dot{\action} = - k_a \Sigma_q^{-1}(\joints - \obs_q) + \Sigma_{\dot{q}}^{-1} (\dot{\joints} - \muvec') $
\vspace{0.2 cm}
\State $Euler \hspace{0.1 cm} integration$
\State $\left[\latent, \muvec', \muvec''\right] \mathrel{+}= \delta_t \left[\dot{\latent},\dot{\muvec}', \dot{\muvec}''\right]$
\EndWhile
\end{algorithmic}
\end{algorithm}

\section{Results}
\label{sec:results}
\subsection{Experiments and evaluation measures}
We performed three different experimental analyses and compared our approach (MAIF) against two state-of-the-art controllers: a proprioceptive torque AIF controller (PAIF)~\cite{pezzato2020novel} and the built-in Panda controller (BPC) with the factory parameters. Video footage of the experiments can be found in the supplementary video: \url{https://github.com/Cmeo97/MAIF}
\begin{enumerate}
    \item \textbf{Qualitative analysis} in sequential reaching (Sec. \ref{sec:results:qualitative}). We evaluated the algorithm accuracy and behavior using several metrics, as well as studied the effect of including the visual modality. Precisely, we analyzed how the multimodal representation power affects accuracy. Besides, we show how our algorithm can work in an imaginary regime by mentally simulating the behavior.
    \item \textbf{Statistical evaluation} in 1000 random goals (Sec. \ref{sec:results:stat}). We compared our proposed controller's average accuracy and robustness to sensory noise against the baselines PAIF and BPC.
    \item \textbf{Adaptation study} (Sec. \ref{sec:results:adapt}). We evaluated the response of the system to unmodeled dynamics and environment variations by altering: the world gravity (Jupiter experiment), the compliance of the robot, and the sensory noise. Here we also compared our approach with the PAIF and the BPC.
\end{enumerate}

In the experiments, we used the following evaluation metrics.
\begin{itemize}
    \item  \textit{Joints perception error}. Error between the inferred (belief) and the observed joint angle. The more accurate are the predictions, the lower will be the perception error.
    \item  \textit{Joints goal error}. Error between the current joint angles and the desired ones (goal).
    \item \textit{Joints belief-goal error}. Error between the inferred joint angle and the desired joints goal. It shows how the perceived joints approach the final goal.
    \item \textit{Image reconstruction error}. Error between the predicted visual input and the observed image. It describes the goodness of the visual generative model.
    \item \textit{End-effector error}. Error between the current end-effector location and the desired one (goal). It provides a single metric to evaluate the precision of the controller.
\end{itemize}

Summarizing, joints perception and image reconstruction errors measure how well the state is estimated, while joints goal, joints belief-goal and end-effector errors give a measure of how well the control task is executed.

\subsection{Experimental setup and parameters}
\label{sec:results:experimental_setup}

 Experiments were performed on the Gazebo simulator and the 7DOF Franka Panda robot arm~\cite{franka} using ROS~\cite{koubaa2019ROS} as interface and Pytorch~\cite{stevens2020pytorch} for the MVAE. A camera model was used to acquire visual grey scaled images with size 1x128x128 pixels. The camera was placed in front of the robot arm with a distance of $1.2$ m. The average frequency of the AIF algorithm was $110$ Hz. Proprioceptive and Visual sensors have publishing frequencies of $1000$ Hz and $10$ Hz respectively. Due to this discrepancy, they were synchronized updating them once per AIF loop iteration.

 The tuning parameters for the AIF controller are:
\begin{itemize}
    \item $\Sigma_{q}, \Sigma_{v}$: Variances representing the confidence about sensory data were set as the variances of the training dataset (joints and images).

    \item $\Sigma_{\mu}\!\!=\!\!2.5, \Sigma_{\mu'}\!\!=\!\!1.0$: Variances representing the confidence of internal belief about the states;
    \item $k_{\mu}\!\!=\!\!11.67,k_q\!\!=\!\!0.6,k_v\!\!=\!\!0.1,k_a\!\!=\!\!900$:  The learning rates for state update and control actions respectively were manually tuned in the ideal settings experiment.
\end{itemize}

Table \ref{tab:network_specs} details the specification of the encoders and decoders used in the MVAE. For training the MVAE we set the number of training epochs to 100, and we used ADAM~\cite{kingma2017adam} as the optimization algorithm. All experiments are simulated on a computer with CPU: AMD ryzen 9 3900X, GPU: Nvidia GeForce RTX 3080 msi. \newpage

\hspace{0.1 cm}
\vspace{0.05 cm}
\hspace{0.1 cm}

\begin{table}[hbtp!]
\scalebox{0.93}{\begin{tabular}{|l|l|l|l|}
\hline
\textbf{Visual Enc.} & \textbf{Prop Enc.} & \textbf{Visual Dec.} & \textbf{Prop Dec.} \\ \hline
\textit{Input} - $\v{I}$   & \textit{Input} - $\joints$    & \textit{Input} - $\latent$  & \textit{Input} - $\latent$ \\
Conv(1, 128, 3)   & FC(7, 64)        & upConv(1, 16, 4)  & FC(256, 128) \\
Conv(128, 64, 3)  & FC(64, 512)      & upConv(16, 32, 4) & FC(128, 64)  \\
Conv(64, 32, 3)   & FC(512, 4096)    & upConv(32, 16, 4) & FC(64, 7)    \\
upConv(32, 32, 3) & Conv(16, 1, 2)   & upConv(16, 1, 4)  & \textit{Output} - $\obs_p$\\
Conv(32 ,16, 2)   & \textit{Output} - $\latent$ & \textit{Output} - $\obs_v$  &              \\
upConv(16, 16, 4) &                &                   &              \\
Conv(16, 1, 2)    &                &                   &              \\
\textit{Output} - $\latent$    &                &                   &              \\ \hline
\end{tabular}}
\caption{Network layers specifications. Type: Convolution (Conv), Trasposed Convolution (upConv), Fully connected (FC). All layers use ReLu activation function, except for the visual decoder output that uses $\tanh(\text{ReLu})$. MaxPool and AvgPool are used to increase network generalization.}
\label{tab:network_specs}
\end{table}

\subsection{Qualitative analysis in a sequential reaching task}
\label{sec:results:qualitative}
To analyse the controller behavior, we designed a sequential reaching task with five desired goals defined by the final joint angles $\{\joints_{d_1}, \joints_{d_2}, \joints_{d_3}, \joints_{d_4}, \joints_{d_5}\}$:
\begin{itemize}
    \item $\boldsymbol{q_{d_1}} = \left[\hspace{0.27 cm} 1.0,\hspace{0.1 cm} 0.5,\hspace{0.1 cm} 0.0,\hspace{0.05 cm} -2.0,\hspace{0.1 cm} 0.0,\hspace{0.1 cm} 2.5,\hspace{0.1 cm} 0.0 \hspace{0.1 cm}  \right]$ rad
    \item $\boldsymbol{q_{d_2}} = \left[\hspace{0.27 cm} 0.0,\hspace{0.1 cm} 0.2,\hspace{0.1 cm} 0.0,\hspace{0.05 cm} -1.0,\hspace{0.1 cm} 0.0,\hspace{0.1 cm} 1.2,\hspace{0.1 cm} 0.9 \hspace{0.1 cm} \right]$ rad
    \item $\boldsymbol{q_{d_3}} = \left[-1.0,\hspace{0.1 cm} 0.5,\hspace{0.1 cm} 0.0,\hspace{0.05 cm} -1.2,\hspace{0.1 cm} 0.0,\hspace{0.1 cm} 1.6,\hspace{0.1 cm} 0.0 \hspace{0.1 cm} \right]$ rad
    \item $\boldsymbol{q_{d_3}} = \boldsymbol{q_{d_2}}$
    \item $\boldsymbol{q_{d_5}} = \boldsymbol{q_{d_1}}$
\end{itemize}
and the desired visual input $\{\v{I}_{d_1},\v{I}_{d_2},\v{I}_{d_3},\v{I}_{d_4}\!\!=\!\!\v{I}_{d_2}, \v{I}_{d_5}\!\!=\!\!\v{I}_{d_1}\}$, depicted in Fig.~\ref{fig:goal_pose}. In all experiments the robot starts in home position ($\joints_{home} \!\!=\!\! \v{0}$ rad), corresponding to the arm straight up.

\begin{figure}[hbtp!]
    \centering
	\subfloat[Home]{
		\centering
		\includegraphics[width=0.21\columnwidth]{Images/starting_pose}
		\label{fig:goal_pose:1}
	}
	\subfloat[$\v{I}_{d_1}$]{
		\centering
		\includegraphics[width=0.21\columnwidth]{Images/goal_pose_1}
		\label{fig:goal_pose:2}
	}
	\subfloat[$\v{I}_{d_2}$]{
		\centering
		\includegraphics[width=0.21\columnwidth]{Images/goal_pose_2}
		\label{fig:goal_pose:3}
	}
	\subfloat[$\v{I}_{d_3}$]{
		\centering
		\includegraphics[width=0.21\columnwidth]{Images/goal_pose_3}
		\label{fig:goal_pose:4}
	}
	\caption{Robot home position and visual goals.}
	\label{fig:goal_pose}
\end{figure}

\subsubsection{Sequential reaching}
Fig. \ref{fig:results} describes the evolution of the evaluation metrics. The joint errors (Fig. \ref{fig:results:1}) show synchronous convergence without overshooting or steady-state errors. Fig. \ref{fig:results:3} and \ref{fig:results:4} show perception and joints belief error respectively. The delay between the convergence of joints beliefs and the real ones is due to the AIF framework. The robot updates its internal belief by approximating the conditional density, maximizing the likelihood of the observed sensations and then it generates an action that results in a new sensory state, which is consistent with the current internal representation. Fig. \ref{fig:results:2} shows the visual generative model image reconstruction errors through the simulation. From one goal to the next one the error drops down. However, some goals can be better reconstructed than others, resulting in different steady errors. The reason is that many $\latent$ solutions lead to similar images. Fig.  \ref{fig:results:5} shows some visual predictions in the 20 seconds simulation. Finally, the generated actions (torques) using MAIF, for the 7 joints during the experiment, are showed in Fig. \ref{fig:torques}. Switching to a new goal generated torques that minimize prediction errors.

\begin{figure}[hbtp!]
	\subfloat[Joints goal error]{
		\centering
		\includegraphics[scale=0.26]{Images/joints_err.pdf}
		\label{fig:results:1}
	}
	\hspace{0.01 cm}
	\subfloat[Image reconstruction error]{
		\centering
		\includegraphics[scale=0.25]{Images/err_Im_attr.pdf}
		\label{fig:results:2}
	}\\
	\subfloat[Joints perception error]{
		\centering
		\includegraphics[scale=0.26]{Images/err_perception.pdf}
		\label{fig:results:3}
	}
	\hspace{0.01 cm}
	\subfloat[Joints belief-goal errors]{
		\centering
		\includegraphics[scale=0.24]{Images/err_attr.pdf}
		\label{fig:results:4}
	}\\
	\subfloat[Sequence of the MVAE predicted visual input $\g_v(\latent)$.]{
		\centering
		\includegraphics[width=0.48\textwidth]{Images/reconstruction}
		\label{fig:results:5}
	}
	\caption{Qualitative analysis of the error measures in the sequential reaching of five goals. All errors present peaks when the new goal is set. (a) Each line represents the error between the i-th joint angle and the desired one. (b) Image reconstruction error. (c) Each line represents the error between the i-th joint belief and the ground truth. (d) Each line represents the error between the i-th joint belief and the desired joint angle. (e) Sequence of the predicted images by the generative model along the trajectory.}
	\label{fig:results}
\end{figure}

\begin{figure}[hbtp!]
\centering
    \includegraphics[width=0.7\columnwidth, height=180px]{Images/torques.pdf}
    \caption{Torque commands for the 7-DOF robot arm using our MAIF in sequential reaching of five goals.}
    \label{fig:torques}
\end{figure}

\subsubsection{Visual modality study}
We evaluated the effect of adding the visual modality. To this end, we tested the algorithm in the presence of visual noise and in the absence of visual input (camera occlusion or broken) by clamping the image to zeros. We compared our algorithm with the proprioceptive-AIF (PAIF) baseline. Visual noise was implemented as additive noise sampled from a Normal distribution $\xnoise_{I} \sim \mathcal{N}(\boldsymbol{0}, \Sigma_I = 0.25 )$.

\begin{figure}[hbtp!]
\centering
    \includegraphics[width=0.8\columnwidth]{Images/err_end_effector_.pdf}
    \caption{Visual influence in sequential reaching in ideal conditions. Lines describe the end-effector location RMSE to goal. Peaks are present when the new goal is set.}
    \label{fig:err_end_effector}
\end{figure}

Fig. \ref{fig:err_end_effector} shows the end-effector location errors using four different algorithm settings: MAIF, MAIF with visual noise, MAIF with visual occlusion and PAIF. MAIF had the best response while PAIF scored the lowest. In the case of visual occlusion, the MAIF approached the PAIF response, and in the case of noisy visual input, accuracy fell in the middle of MAIF and MAIF with visual occlusion.

\subsubsection{Mental simulation}
Unlike most of the controllers present in literature, a great advantage of using our approach is the possibility to perform simulations on the robot's brain exploiting the generative nature of the MVAE. In other words, given $\obs_{desired}$, the entire experiment can be mentally simulated. Since sensory data are not available, the state update law becomes:
\begin{equation}
    \dot{\latent} = \dot{z}_{q_d}+\dot{z}_{v_d}
\end{equation}
As a result, decoding the new internal state, the updated $\{\v{I},\joints\}$ can be computed and the new errors can be back-propagated again, creating a loop that allows the system to do imaginary simulations.

\begin{figure}[hbtp!]
	\centering
	\subfloat[Joints errors]{
		\centering
		\includegraphics[width=0.45\columnwidth]{Images/imagined_joints_.pdf}
		\label{fig:im_results:1}
	}
	\hspace{0.01 cm}
	\subfloat[Image reconstruction error]{
		\centering
		\includegraphics[width=0.45\columnwidth]{Images/reconstruction_imagined.pdf}
		\label{fig:im_results:2}
	}
	\caption{Mental simulation of sequential reaching of five goals. The goal is updated on time steps where peaks are present. (a) Joints errors of an imagined simulation. Each line represents the error of the i-th joint. (b) Image reconstruction errors of an imagined simulation.}
	\label{fig:imagined_results}
\end{figure}

Fig. \ref{fig:im_results:1} and \ref{fig:im_results:2} show respectively imagined joints error and images reconstruction error through the entire simulation. These results show that the errors converge faster to zero than in the normal regime (Fig. \ref{fig:results:4}) as it does not need to accommodate the real dynamics of the robot.

\subsection{Statistical analysis}
\label{sec:results:stat}
A statistical comparison between our MAIF, PAIF and BPC was performed to evaluate the controllers' performances on reaching tasks under normal and noisy conditions. To this end, we created a list of $1000$ random desired goals that we executed in each controller. Then, we compared their average joint errors.
\begin{figure}[hbtp!]
	\subfloat[No noise]{
		\centering
		\includegraphics[width=0.44\columnwidth]{Images/comparison.pdf}
		\label{fig:stat_results:1}
	}
	\hspace{0.01 cm}
	\subfloat[Proprioceptive noise]{
		\centering
		\includegraphics[width=0.44\columnwidth]{Images/comparison_noise.pdf}
		\label{fig:stat_results:2}
	}
	\caption{Statistical comparison of our MAIF (green), PAIF (orange) and the BPC (blue) over 1000 executions with different random goals in ideal conditions (a) and with proprioceptive noise (b). Lines represent the Root Mean Square Error~(RMSE) of the end-effector and shadings describe their standard deviation.}
	\label{fig:stat_results}
\end{figure}
In every execution, the robot started from the home position. Fig. \ref{fig:stat_results} shows the root-mean-square error (RMSE) between the end-effector position and the desired pose over all executions for the three controllers.

In perfect conditions (Fig.\ref{fig:stat_results:1}), BPC was the fastest and had the lowest RMSE on average. MAIF had high accuracy and presented a small overshoot when reaching the goal. PAIF was much slower and had lower accuracy. However, when injecting Gaussian noise ($\xnoise_q \sim \mathcal{N}(\boldsymbol{0}, \Sigma_q = 0.5 )$) in the joint sensor measurements (Fig. \ref{fig:stat_results:2}), our proposed MAIF scored the lowest RMSE and standard deviation. Hence, MAIF presented higher robustness to proprioceptive noise with more stable behavior than the BPC and the PAIF. BPC had the highest standard deviation, showing poor robustness to sensory noise, especially we approaching the desired goal. By visual inspection of the robot behavior, we observed strong oscillations in the BPC. The robustness of MAIF is due to the algorithm structure. The free energy optimization scheme and the decoder interpolation make that the sensory noise has a low effect on the state estimation.

\subsection{Adaptation}
\label{sec:results:adapt}
To investigate our approach adaptability to unmodeled dynamics and environment variations we tested the controllers in three experiments. First, we changed the gravity of the world as if the robot was in Jupiter. Second, we altered the motors stiffness parameter of the Panda robot. Third, we reevaluated the controller in the presence of sensory noise focusing on the robot behavior. Again, we compared our algorithm (MAIF) with the built-in Panda controller (BPC) and the proprioceptive controller (PAIF). All controllers parameters were the same as in the previous experiments. Thus, no retuning was done.

\subsubsection{Jupiter experiment}
This experiment aims to show the controller performance in case of gravity changes. The reasoning behind this experiment is to evaluate MAIF performance under unexpected external conditions, such as external forces.
\begin{figure}[hbtp!]
    \centering
    \includegraphics[width=0.8\columnwidth]{Images/jupiter_exp.pdf}
    \caption{Jupiter experiment. The lines describe the end-effector error to the desired goal under Jupiter gravity in a reaching task with five sequential goals. MAIF (green), PAIF (orange) and BPC (blue). All end-effector errors present peaks when the new goal is set.}
    \label{fig:jupiter_exp}
\end{figure}
To do that, we performed the simulation with Jupiter gravity, which is $g = 24.79 \frac{m}{s^2}$ on the planet's surface. Fig. \ref{fig:jupiter_exp} illustrates joint angle absolute errors for each controller. MAIF was the least affected by the change of gravity. Under BPC the robot arm oscillated around the desired pose and had steady-state errors. PAIF response was also affected by the gravity change. MAIF registered a faster response than PAIF due to multisensory integration: images are not affected by gravity.

\subsubsection{Compliant experiment}
The goal of this experiment is the evaluation of MAIF adaptation under changes on the robot's physical parameters. We altered the robot arm's compliance by changing the motors stiffness values from $0.61 \frac{Nm}{rad}$ to $0.01\frac{Nm}{rad}$. Figure \ref{fig:compliant_exp} shows the absolute joints errors for each controller. As in the previous experiment, the BPC makes the robot arm oscillate around some goal poses, while both MAIF and PAIF do not. With low stiffness, MAIF and PAIF slightly overshot when the desired goal was far away. However, since MAIF updates his state from both proprioceptive and visual sensory data, it had much lower overshooting and faster response than PAIF.

\begin{figure}[hbtp!]
    \centering
    \includegraphics[width=0.8\columnwidth]{Images/compliant_exp.pdf}
    \caption{Compliant experiment. Sequential reaching of five goals when motors stiffness is changed from $0.61 \frac{Nm}{rad}$ to $0.01\frac{Nm}{rad}$.  All end-effector location absolute errors present peaks when the new goal is set. }
    \label{fig:compliant_exp}
\end{figure}
\subsubsection{Sensory noise experiment}
We reevaluated the controller behavior in the presence of proprioceptive noise but in a sequential reaching task, focusing on the adaptation capabilities of the three controllers. Figure \ref{fig:noisy_exp} shows that MAIF was the most adaptive, presenting the smoothest behavior. PAIF also adapted to the noise but presented higher oscillations than MAIF. By contrast, BPC had the worst behavior, showing strong oscillations around the desired goals.

\begin{figure}[hbtp!]
    \centering
    \includegraphics[width=0.8\columnwidth]{Images/noisy_exp.pdf}
    \caption{Sensory noise experiment. Sequential reaching of five goals under proprioceptive noise. All end-effector absolute errors present peaks when the new goal is set.}
    \label{fig:noisy_exp}
\end{figure}

\section{Conclusion}
\label{sec:conclusions}

We presented an improved adaptive torque controller for robot manipulators based on active inference. Our approach makes use of the alleged adaptability and robustness of AIF, taking advantage of previous works and overwhelming some related limitations. We solved state estimation by combining representation learning (with a multimodal VAE) with free energy optimization, improving the representational power and adaptability. As a result, we derived a schema for online proprio-visual torque control, which does not require any dynamic or kinematic model of the robot, is less sensitive to unmodeled dynamics, and can handle high-dimensional inputs. Although in this work we used images it can be generalized to any sensor modality.

Results showed that our proposed algorithm outperforms in terms of adaptation the built-in Panda controller and a state-of-the-art torque AIF baseline. Moreover, it was more accurate in the presence of sensory noise. Our AIF controller was highly adaptive and robust to different contexts, such as sensory noise, visual occlusions, changes in the world dynamics (i.e., gravity) and changes in the robot properties (i.e. motor stiffness). Besides, our approach natively allowed mental simulation of the planned trajectory.

Future work will finalize the physical robot experiments that were not possible due to university restrictions and will investigate strongly coupled cross-modal learning and hierarchical schemes for sequential decision making.

\section*{Acknowledgment}
\addcontentsline{toc}{section}{Acknowledgment}
We would like to thank Professor Martijn Wisse for his fruitful comments on the elaboration of this work.

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv