23-140

# Interstellar Object Uncertainty Evolution and Effect on Fast Flyby Delivery
and Required Delta-V

Declan M. Mages Mission Design and Navigation, Jet Propulsion Laboratory,
California Institute of Technology. Davide Farnocchia11footnotemark: 1
Benjamin Donitz Project Systems Engineering and Formulation, Jet Propulsion
Laboratory, California Institute of Technology.

###### Abstract

Interstellar objects (ISOs) are small bodies that can travel through our solar
system from other star systems. When present in our solar system, they
represent an opportunity to study the properties and origins of these objects,
as well as the potential for cross-pollination of material between star
systems. With current propulsion technology, rendezvous with these objects is
likely infeasible, and thus the maximum science return results from a rapid
response flyby and impactor. However, while trajectories to ISOs may be
feasible, their potentially high ephemeris uncertainties and high-speed
hyperbolic orbits present significant challenges to navigation. In this paper
we assess these challenges by modeling the uncertainties of reachable
synthetic ISOs as a function of time, as derived by measurements from ground
observatories and an approaching spacecraft. From these uncertainties we
derive the final delivery accuracy of fast flyby spacecraft to the ISO and
required statistical delta-v for navigation. We find that these two challenges
can lead to hundreds of meters-per-second or even kilometers-per-second of
required statistical delta-v for navigation, reduce delivery accuracy to
hundreds of kilometers, and make autonomous navigation a requirement.

##  1 INTRODUCTION

With the discovery of 1I/‘Oumuamua [1] and 2I/Borisov [2], interstellar
objects (ISOs) have been thrust into the forefront of planetary science and
the robotic exploration of ISOs is now of significant interest. While there is
significant scientific understanding to be gained from observing these objects
remotely, ground based observations are far, and in situ characterization
remains the most powerful means to studying these objects However, ISOs
present unique challenges to mission design and navigation. They by nature
pass through our solar system on hyperbolic trajectories as depicted in Figure
1. They travel at extreme speeds and spend only a few years inside Saturn’s
orbit before exiting the solar system forever. Even when the bodies are active
and posses comas, they are still dim and thus difficult to detect at distances
greater than 10 AU from the Sun, which means there is limited time to respond
with an intercept mission. Given these extreme velocity trajectories and
limited warning times, rendezvous opportunities with ISOs are likely
prohibitive with current and near term propulsion technology [3]. Thus, the
maximum science return would result from a rapid-response fast flyby with
possible impactor spacecraft, akin to the Deep Impact mission [4] and Comet
Interceptor mission [5].

Figure 1: Trajectories of ’Oumuamua and Borisov.

However, this rapid-response mission architecture presents its own set of
challenges. The hyperbolic orbits drive high relative encounter velocities and
high solar phase angles (the angle between vector from target body to Sun and
vector from target body to spacecraft), while short observation arc lengths
that are intrinsic to a “rapid-response” mission produce large ISO position
uncertainties that can change dramatically during cruise. The motivation for
this work was the concern that high uncertainties during cruise and approach
would require prohibitive amounts of delta-v to correct for and limit science
return during flyby. To assess this we first establish a population of
synthetic ISOs, and then determine which are reachable and the associated
spacecraft trajectories. Then we model the uncertainties of all ISOs and
spacecraft pairs over the course of each mission starting at launch, and
determine the distribution of delta-vs required for navigation and the
distribution of final flyby delivery accuracies.

##  2 POPULATION MODELING

###  2.1 ISO Population

In the generation of the ISO population we leveraged previous work by T.
Engelhardt et al. [6]. From an initial population of 1.7 billion synthetic
ISOs, modeled as both active and inert, they derived the ISOs that would have
appeared in any Pan-STARRS1 [7] and Catalina Sky Survey [8] telescope fields
taken between between 2005-January-1 and 2015-January-1. T. Engelhardt et al.
defined as “detected” any ISO that resulted in a digest score [9] greater than
or equal to 90. We instead assume that any active ISO is “detected” in the
first field in which it has an apparent magnitude brighter than 21, i.e.,
around the current limiting magnitude of the Catalina Sky Survey. In fact, the
detection of cometary activity would ensure that the newly discovered object
would get posted on the Possible Comet Confirmation
Page111https://minorplanetcenter.net/iau/NEO/pccp_tabular.html and in turn
followed up by other observers. On the other hand, for inert ISOs we
additionally require a digest score greater than or equal to 65, which is the
threshold that triggers posting on the NEO Confirmation
Page222https://minorplanetcenter.net/iau/NEO/toconfirm_tabular.html and, in
turn, follow-up.

Of the 1.7 billion synthetic ISOs, 165,378 of them pass within 50 AU of the
Sun, and of those 7,813 ISOs pass within 10 AU. Each ISO is studied for
detectability as if it were active and as if it were inert. The active ISOs
are assumed to become active at 10 AU, becoming significantly brighter and in
all cases reaching an apparent magnitude brighter than 21, and thus becoming
potentially detectable. When farther than 10 AU from the Sun, with no modeled
activity, only a few ISOs reached apparent magnitude 21. All 7813 active ISOs
that passed within 10 AU were detected, versus only 371 inert ISOs detected.
Figure 2 shows the distribution of distances at which the ISOs are first
detected and of the ISOs’ perihelion distance. Active ISOs are significantly
brighter than inert ISOs and so are detected roughly 2-3 times farther out
from the Sun. However, while all the modeled active ISOs maintain detectable
magnitudes within 10 AU, the referenced surveys do not cover the entire sky,
and so there are trajectory geometries such that a bright ISO does not appear
in a field till much later. 2I/Borisov is a good example as it was discovered
by an amateur astronomer specifically looking in a region of sky closer to the
Sun and thus not covered by current surveys [2].

Figure 2: Histograms of detected ISOs’ detection distances and perihelion
distances.

This methodology produces an ISO population with detection times given current
survey telescopes. Future surveys, particularly the Vera C. Rubin
Observatory[10] and NEO Surveyor[11], will significantly improve our detection
capabilities. Vera Rubin will be performing daily full sky surveys down to the
24th magnitude while NEO Surveyor hunts for dim near-Earth objects in the
infrared. Future analysis will also have to consider ISO activity and
increased brightness outside 10 AU in order to better assess potential
detection times at dimmer magnitudes.

###  2.2 Spacecraft Trajectory Design

With the 8,184 total detected ISOs we applied small satellite (SmallSat)
system design limitations to constrain a patched-conic search for trajectories
that depart Earth to intercept each ISO. The constraints are listed in Table 1
and outlined in more detail along with the trajectory design process in
previous work, “Interstellar Object Encounter Trade Space Exploration”[12]. Of
the 8,184 detected ISOs, 205 active and 5 inert ISOs can be reached for a fast
flyby/impact. While this seems like a negligible fraction, only 931 active and
231 inert ISOs passed within 2.5 AU of the Sun and were thus even within the
trajectory design constraints. This gives an intercept chance of 22% and 2%
respectively for active and inert objects that pass within constraints. Three
example ISO and spacecraft trajectories are shown in Figure 3. Inclination
changes are extremely costly in terms of delta-c, and so the feasible
intercept trajectories are largely constrained to the ecliptic. One additional
note is that if the launch timing constraint is completely removed, i.e. we no
longer care when the ISO is detected, a total of 549 ISOs can be reached for
intercept. This is a representative upper bound of possible performance that
could be gained from future more powerful observatories.

Figure 3: Three example ISO and parent feasible spacecraft trajectories. ISO
is plotted starting 1 year out from intercept.

The minimum launch characteristic energy (C3) trajectories produce the
encounter space depicted in Figure 4. Because of their hyperbolic
trajectories, ISOs pass through the inner solar system at extreme heliocentric
speeds and at all angles relative to the ecliptic. This means the high speed
flyby encounters can be both prograde and retrograde, resulting in the wide
range of relative velocities from 20 km/s to nearly 100 km/s.

With their high speeds, ISOs are traveling significantly faster compared to
the flyby spacecraft, given our C3 constraints. As such, the encounter should
not be thought of as the spacecraft going to and flying by the ISO. Instead,
the spacecraft just seeks to get in the path of the ISO, usually where it
pierces or nears the ecliptic. This in turn means that if the ISO is
intercepted pre-perihelion, during the inbound portion of its trajectory, the
spacecraft will be looking out and away from the Sun, observing a well
illuminated target at low solar phase angles. However, if the ISO is
intercepted post-perihelion, while outbound, it will typically approach the
spacecraft from the direction of the Sun and at high solar phase angles, as
highlighted by the right plot in Figure 4. Time is a precious resource for
rapid response missions, and since post-perihelion intercepts inherently give
more time from detection to intercept, they are more likely, representing 61%
of cases, with 50% of all cases having solar phase angles greater than 100∘.
These high solar phase angles will require a highly capable camera and can
make it difficult for the spacecraft to image and detect the ISO.

Table 1: Constraints applied to patched-conic trajectory search Min Time from Detection to Launch | 10 days
---|---
Max Time from Launch to encounter | 3 years
Max orbital revolutions | 4 before encounter
Max $C_{3}$ | 60 km2/s2
Max Phase Angle | 160∘
Max Relative Velocity | 100 km/s
Max Solar Distance at Encounter | 2.5 AU
Flyby Distance | 1000 km
Figure 4: Left, flyby/impact solar phase angle versus relative velocity.
Right, flyby/impact solar phase angle versus intercept date relative to ISO
perihelion. Figure 5: Histogram of spacecraft wait times and response times.

##  3 NAVIGATION ANALYSIS

###  3.1 Small Body Deep Space Navigation

When navigating to a small body in our solar system, the goal is to deliver
the spacecraft to a specific flyby location relative to the target. To
accomplish this, two key pieces of information need to be estimated: the
spacecraft’s trajectory and the small body’s ephemeris. The spacecraft’s orbit
determination (OD) is primarily done through tracking data in the form of two-
way Doppler, range, and delta-differential one-way ranging with the Deep Space
Network (DSN). The target’s ephemeris is then computed using primarily optical
astrometric measurements from ground-based observatories and improves as the
observation arc length is extended. As the spacecraft approaches the target,
optical navigation (opnav), the technique of imaging a target with a camera
on-board the spacecraft[13], is also used. Opnav data can directly measure the
relative state and improve the target’s ephemeris and is powerful due to the
fact it is an angular measurement whose metric resolution increases as the
spacecraft gets closer to the target. During mission development it is an
important step to perform navigation analysis and model the combination of
these measurements to predict the uncertainties of the spacecraft and ISO
trajectories over time, and make adjustments if needed. The uncertainty over
time results are critical to understanding the statistical delta-v required
for correction maneuvers and understanding the final delivery accuracies that
determine science collection.

The results from ISO population formulation and spacecraft trajectory design
ultimately gave 210 unique reachable ISOs and timelines, and for each unique
ISO a set of dozens of possible spacecraft trajectories. Typically, navigation
analysis involves evaluating a single target body ephemeris with an assumed
unchanging ground-based uncertainty, and a single spacecraft trajectory. The
challenge this ISO analysis presents is that there are over 200 target
ephemerides, each with uncertainties that change dramatically over time – and
also multiple possible spacecraft trajectories. To perform the broad
navigation analysis, we first needed to constrain the problem, and so we
nominally selected the lowest C3 spacecraft trajectory associated with each
ISO. Some alternate trajectories selection methodologies are also discussed in
the results. With each spacecraft and ISO pair, the next step was to model
their uncertainties over time. The process and the assumptions used is
outlined in the following sections.

The last step in our process was to determine the 99th percentile required
delta-v (DV99) and final delivery accuracy. This is accomplished by feeding
the estimated uncertainty profiles into LAMBIC (Linear Analysis of Maneuvers
with Bounds and Inequality Constraints), which produces these statistics by
simulating the execution of a sequence of maneuvers through a Monte Carlo
process[14]. The maneuvers are designed utilizing critical plane targeting,
which targets flyby location but not when exactly the flyby will occur. This
approach reduces the maneuver design constraints and therefore the required
delta-v, especially with high uncertainty targets where relative distance and
thus flyby time is poorly known. The coordinate frame depicting delivery
results is the b-plane [15]; where the $Z$-axis is along the spacecraft-ISO
relative velocity vector (here referred to as the downtrack direction), $Y$
axis perpendicular to the relative velocity vector and laying in the ecliptic,
and $X$ completing the right-handed frame ($X$ and $Y$ referred to here as
crosstrack directions).

###  3.2 ISO Uncertainty Modeling

The ephemeris uncertainty for a discovered ISO affects the feasibility of a
mission, especially a rapid response one. We derived the ISO ephemeris
uncertainty by simulating optical astrometric observations assuming the
cadence described in Table 2. The astrometric observations were assumed to
have an uncertainty of 0.5 arcsec, which is in line with the typical
performance of current observing programs [16]. To mimic a realistic evolution
of the ephemeris uncertainty and the real-time uncertainty information that
would be available for navigation, we started at the time of discovery and
progressively advanced the data cutoff with a time step of one week. For each
data cutoff, we computed a best fitting orbit solution and mapped its
covariance to compute the future position and velocity uncertainties.

In terms of covariance analysis, one important consideration is the effect of
nongravitational perturbations (nongravs), which are typically modeled as
radial, transverse, and normal accelerations $A_{1}g(r)$, $A_{2}g(r)$, and
$A_{3}g(r)$, respectively, where the $g(r)$ captures the dependency on the
heliocentric distance [17]. The motion of the first two interstellar objects
was significantly affected by nongravitational perturbations, i.e., for
‘Oumuamua [18] $A_{1}$ was $3\times 10^{-7}$ AU/d2, while for Borisov [19]
$A_{1}$ was $7\times 10^{-8}$ AU/d2. To analyze the dependency of the ISO
uncertainty on nongravitational perturbations, we considered three scenarios,
similarly to what was done for comet C/2013 A1 (Siding Spring) and the
analysis of its close approach to Mars in 2014 [20]. The first scenario is the
optimistic, gravity-only case, i.e., without nongravitational perturbations.
As the baseline, best-guess scenario (referred to as “reference”), we estimate
$A_{1}$, $A_{2}$, and $A_{3}$ as part of the orbital fit by assuming a priori
uncertainties of $10^{-7}$ AU/d2, $10^{-8}$ AU/d2, and $10^{-8}$ AU/d2,
respectively. Finally, as a worst-case scenario (referred to as “wide”), the
$A_{1}$, $A_{2}$, and $A_{3}$ as a priori uncertainties are $10^{-6}$ AU/d2,
$10^{-7}$ AU/d2, and $10^{-7}$ AU/d2, respectively. This setup produces three
distinct uncertainty evolution profiles for each ISO trajectory. We only
consider the gravity-only case for the 5 inert ISOs while the 205 active ISOs
are modeled with both reference and wide nongravitational perturbations
(nongravs).

Figure 6 shows an example of the resulting ISO uncertainty profiles in the
b-plane frame at the spacecraft encounter time versus data-cutoff. This
example is indicative of many ISO uncertainty profiles. They start with high
initial uncertainties that slowly improve as the observation arc length
increases and have unobservable periods than can halt ephemeris improvement. A
more statistical representation of the population of ISO uncertainties is then
depicted in Figure 7. This shows how for both reference and wide nongrav
cases, the 90th percentile uncertainty is roughly 108 km after one month of
observation and roughly 106 km after three months of observation. After that
the reference case uncertainties continue to steadily decline – versus the
wide nongrav uncertainties, whose higher magnitude accelerations limit the
ability to accurately predict the ISOs location at intercept.

Table 2: Inputs to ISO covariance evolution analysis. ISO Covariance Setup
---
Time of Discovery |  First detection where ISO V magnitude $<$ 21 and solar elongation $>$ 90∘
Observation Schedule |  Daily observations if V $<$ 24 and elongation $>$ 90∘, biweekly observations if $24<V<26$ and elongation $>$ 60∘, weekly observations otherwise (as long as $V<26$ and elongation greater than $60^{\circ}$)
Measurement Uncertainty |  0.5 arcsec
Nongrav a priori Uncertainty |  Ref Case: $A_{1}=10^{-7}$ AU/d2, $A_{2,3}=10^{-8}$ AU/d2
Wide Case: $A_{1}=10^{-6}$ AU/d2, $A_{2,3}=10^{-7}$ AU/d2
Spacecraft Covariance Setup
Initial State Uncertainty |  100 m, 1 mm/s - position, velocity 1$\sigma$
Tracking Schedule |  Four 8 hour Doppler and range passes per week
OpNav Schedule |  Every 12 hours between E-50 days and E-1 day
Measurement Weights |  Two-way Doppler 0.0112 Hz, Two-way SRA range 35 ru, OpNav 1 arcsecond
Maneuver Schedule |  1-2 cruise maneuvers timed to minimize $\Delta V$, 3 approach maneuvers at E-40, E-20, E-1 days
Maneuver Data-cuttof |  24 hours

Figure 6: Example ISO 1-sigma ephemeris uncertainty over time when mapped to
b-plane at encounter. Left is reference nongravs, right models wide nongravs.
Figure 7: 1-sigma uncertainties of where each ISO will be at targeted
encounter time versus ISO observation arc-length.

###  3.3 Spacecraft Uncertainty Modeling

Each spacecraft’s uncertainty is modeled by simulated Doppler, range, and
opnav measurements [21]. Typically, with a fixed target and trajectory each
mission would have its own detailed custom DSN tracking schedule, while the
unique geometry and target brightness would determine the opnav imagery
schedule. Then any number of trajectory correction maneuvers would be placed
throughout the trajectory to maximize performance. However, in order to
analyze the 210 spacecraft and ISO pairs, we apply to all a generic data
collection and maneuver timing strategy as outlined in Table 2 and depicted in
Figure 8. The only changes to this schedule were made to investigate the
effect of late opnav detection on navigation given most ISO encounters involve
high solar phase angles. The emphasis of this analysis is also specifically on
how the unique nature of a high speed rapid response mission inflates
uncertainties. We therefore chose to only model a simple and benign spacecraft
– estimating only the dynamic state of the spacecraft and ISO while modeling
only maneuver execution error and bias solar radiation pressure parameters and
no injection error from a launch vehicle.

For the maneuver timeline, we assumed two trajectory correction maneuvers
(TCMs) during cruise – where cruise is defined as the time from Earth
departure to start of opnav. These two maneuvers are initially spaced evenly
apart. A check is then done to see if either maneuver is located during a
period when the ISO has a low elongation angle and can not be observed. If a
maneuver is initially during that period, it is instead placed one week after
the next ISO ground astrometry measurement so to be able to take advantage of
the improved ephemeris knowledge from the new measurements. This is fairly
simple “optimization” and possible improvements to maneuver placement is
discussed later as future work. During approach there are three maneuvers
placed at encounter minus 40 days, 20 days, and 1.5 days, where approach is
defined as the time from start of opnav (nominally 50 days from encounter) to
the last ground in the loop involvement.

Figure 8: Nominal spacecraft timeline.

##  4 ACTIVE ISO RESULTS

###  4.1 Delivery Accuracy and Statistical Delta-V

The resulting delivery accuracies from the final maneuver at 1.5 days from
encounter are depicted in Figure 9 and reveal crosstrack deliveries that are
significantly worse compared to typical flyby missions. This is a result of
degraded opnav measurement power because the high relative velocities lead to
high relative distances. The high crosstrack delivery accuracies are nearly
1:1 correlated with the high relative velocities and this relationship is
emphasized by the fact that the larger ISO ephemeris uncertainties from wide
nongravs have no effect on crosstrack delivery. The opposite is true regarding
the downtrack delivery as opnav cannot measure the distance to a target. The
downtrack delivery is therefore largely un-observable and thus defined by the
ground-based ISO ephemeris uncertainty. In the reference nongrav case this
results in thousands of kilometers of uncertainty, while in the wide nongrav
case the uncertainties can be on the order of tens of thousands of kilometers.

Figure 9: Crosstrack and downtrack delivery uncertainty.

These results show that given the assumed optical navigation camera accuracy,
the crosstrack delivery error can be a significant percent of the targeted
1000 km flyby distance – at the extreme cases more than 50%. The downtrack
error results then show that imaging the ISOs during flyby will likely require
autonomous tracking. At closest approach the downtrack uncertainty is fully
mapped into the imager or other instrument’s FOV – at 1000 km away this means
a downtrack uncertainty of 2000 km represents a roughly 90∘ long strip of sky.
Only with very late autonomous opnav data, when the spacecraft is starting to
pass along-side the ISO and the downtrack uncertainty starts to map into the
image FOV, can the spacecraft then autonomously estimate when close approach
will occur.

Figure 10: Active ISO DV99s.

A histogram of the required DV99 to navigate to active ISOs is depicted in
Figure 10. These results show that navigating to these highly uncertain
targets can require very significant statistical delta-v, particularly in
cases with high nongravitational accelerations. With reference nongravs, the
majority-75% of the cases require below 100 m/s of DV99. However, roughly 15%
of cases are shown to require in excess of 500 m/s of DV99. With wide nongravs
just 25% of cases required less than 100 m/s of DV99. The 75th-percentile DV99
case with wide nongravs is shown to be a little over 500 m/s.

###  4.2 High Delta-V Cases

While much of the cases’ delta-vs fit within a reasonable distribution that
tapers off away from the mean, there are a non-negligible $>$10% number of
cases on the edge of the histogram charts that represent multiple kilometer-
per-second DV99 cases. These cases require likely infeasible amounts of
statistical delta-v. Upon investigation, we determined that the driver of the
large DV99 was early launches with a minimum C3 trajectory that departs from
Earth only weeks after detection. At this time the ISO’s ephemeris is very
poorly known as depicted previously in Figure 7. Therefore after launch the
spacecraft must execute large maneuvers during cruise in order to effectively
re-target the new ISO trajectory as more ground data refines its predicted
path. Figure 11 depicts how requiring longer wait times will remove some of
the extreme delta-v cases but will also remove some otherwise viably reachable
targets.

Figure 11: Active ISO DV99s depending for different wait time constraints.

By imposing a 30 day wait constraint, 25 spacecraft trajectories no longer
minimize C3, and 10 ISOs are no longer reachable. Imposing a 180 day wait
results in 37 spacecraft trajectories that no longer minimize C3, and 36% or
73 ISOs that are no longer reachable. The histogram plots of these new
trajectories’ DV99s clearly shows almost all of the $>$1 km/s cases have been
eliminated and skewed lower. This is completely true with the reference
nongrav ISOs, with all 30 day delayed launches below 200 m/s and all 180 day
delayed launches below 100 m/s. However, for wide nongrav ISOs the enforced
wait times do not remove all of the multi-kilometer-per-second cases, and many
still result in DV99s in the hundreds of meters-per-second range.

An investigation of those high delta-v, long-wait cases highlighted a
secondary reason for large ephemeris uncertainties: ISO low elongation angle
observation constraints. Some ISOs, despite being detected far in advance,
become un-observable for long periods, when the ephemeris cannot be improved.
This happens for a variety of reasons. For example, ISOs approach the inner
solar system from a fairly fixed region of sky and so as the Earth’s continues
on its orbit, the Sun can pass near the ISO’s approach asymptote, thus
preventing Earth-based observations. Poor timing of a low elongation period
can result in high uncertainties at the start of opnav, thus leaving the
spacecraft with a significant ephemeris update to make and only 40 days left
to make a correction. An example ISO that highlights this behavior is shown in
Figure 12.

Figure 12: ISO that enters long period of low elongation angle before
encounter.

###  4.3 Late Detection

Another important navigation challenge to highlight is the possibility of late
opnav detection. Late detection means there is less time to correct the
ground-based ephemeris, which increase the delta-v costs. For all previous
analysis opnav was assumed to start at 50 days from encounter, which requires
that the target ISO be bright enough to detect with the onboard imager at that
time. As 61% of encounters involve high solar phase angles this detection may
be difficult. For active ISOs, their comas scatter light forward[22] and so
should enable early opnav detection even at high phase – however for low-
activity ISOs this may not be the case. We therefore reran the navigation
analysis assuming detection 5 days from encounter along with new maneuver
assumptions, listed in Table 3, to assess how late detection would affect
delta-v.

The 40 day delay in opnav measurement collection results in a modest increase
in DV99 for reference nongrav ISOs and an enormous increase in DV99 for wide
nongrav ISOs as depicted in Figure 13. The reference nongrav increase is
mostly concentrated around 100 m/s, with some large 1 km/s increases coming
from highly uncertain ISOs with substantial gaps in ground observability. For
the wide nongrav ISOs, the DV99 increases are almost all above 1 km/s. This
highlights how nongravitational forces can severely limit the ISO ephemeris
accuracy – and if the spacecraft is left to correct for that with a maneuver
only 4 days from encounter, the delta-v will be massive. Luckily these wide
nongrav ISOs are synonymous with highly active ISOs, and so should have
detectable coma.

Table 3: Late detection new inputs. Spacecraft Covariance Setup
---
OpNav Schedule |  Every 12 hours between E-5 days and E-1 day
Maneuver Schedule |  2 approach maneuvers at E-4 and E-1 days
Figure 13: Increase in DV99 due to opnav detection at 5 days versus 50 days.

##  5 INERT ISO RESULTS

Even with just five ISOs, each case presents unique mission parameters – and
with the small sample size each of their results can be outlined in more
detail. For each inert ISO, with its unique ID, the minimum C3 trajectories
with 10, 30, and 180 day wait constraints were analyzed, with results depicted
in Table 14. The different wait time constraints had no effect except for ISO
464888 – for which a global minimum C3 trajectory had a wait time of 15 days
and still had feasible trajectories after waiting 30 days, unlike for ISO
275123.

ISO 23184 is shown to have a high DV99 even with a wait time of over 30 days.
Investigation into the ISO’s uncertainty profile showed, as expected, a long
gap in ground observation data soon after detection that prevents ephemeris
improvement. ISO 130419 is a reasonable target, with a launch wait time over
30 days and good ISO observability and moderate relative velocity – resulting
in a small DV99. ISOs 275123 is the most challenging target with a 13 day wait
resulting in a DV99 of 588 m/s. It involves a very high solar phase angle
approach which could make opnav imaging extremely difficult. ISO 464888
highlights how waiting 15 days to 30 days can reduce the DV99 from 197 m/s to
41 m/s. And finally ISO 650829 is a nearly idyllic case with a 130 day wait,
only 29 km/s relative velocity, and modest solar phase angle.

Figure 14: Table of results on 5 reachable inert ISOs.

##  6 CONCLUSION

Our results show that navigating a rapid response mission to an ISO can be
extremely challenging. The required delta-v just for statistical trajectory
correction maneuvers is comparable to some missions’ entire deterministic
delta-v. Navigation is rarely considered in early mission development, however
for rapid response missions to ISOs these navigation challenges must be
modeled from the beginning. If not, the statistical corrections have the
potential to break a mission’s delta-v budget. There ultimately may be ISOs
that, while accessible from a trajectory design perspective, may not be
feasible targets due to extreme navigation delta-v.

The high relative velocities are also shown to limit delivery accuracies to
hundreds of kilometers. To significantly improve this, autonomous OD and
maneuver execution is likely needed [23]. High ISO ephemeris uncertainties are
also shown to massively inflate the downtrack uncertainties during flyby.
Meaningful science collection will therefore likely require autonomous
tracking capabilities.

Other conclusions from this analysis include the requirement for a dynamic
navigation team that is closely integrated with the ground-based ISO OD team.
Nongravitational accelerations are also shown to be possibly the most critical
parameter in determining navigability. A highly capable opnav imager is also
needed for long range detection even at high solar phase angles. And finally,
a counter-intuitive result is that for a rapid response mission, rather than
launching to the target as soon as possible, it can instead be more optimal to
launch as late as possible.

##  7 ACKNOWLEDGEMENTS

The research described in this paper was carried out at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration. Copyright (c) 2022 California
Institute of Technology. U.S Government sponsorship acknowledged.

## References

  * [1] K. J. Meech, R. Weryk, M. Micheli, J. T. Kleyna, O. R. Hainaut, R. Jedicke, R. J. Wainscoat, K. C. Chambers, J. V. Keane, A. Petric, L. Denneau, E. Magnier, T. Berger, M. E. Huber, H. Flewelling, C. Waters, E. Schunova-Lilly, and S. Chastel, “A brief visit from a red and extremely elongated interstellar asteroid,” Nature, Vol. 552, Dec. 2017, pp. 378–381, 10.1038/nature25020.
  * [2] P. Guzik, M. Drahus, K. Rusek, W. Waniak, G. Cannizzaro, and I. Pastor-Marazuela, “Initial characterization of interstellar comet 2I/Borisov,” Nature Astronomy, Vol. 4, No. 1, 2020, pp. 53–57.
  * [3] D. Miller, F. Duvigneaud, W. Menken, D. Landau, and R. Linares, “High-performance solar sails for interstellar object rendezvous,” Acta Astronautica, Vol. 200, 2022, pp. 242–252.
  * [4] M. F. A’Hearn and L. N. Johnson, “Deep Impact and Related Missions,” Handbook of Cosmic Hazards and Planetary Defense, pp. 513–534, 2015, 10.1007/978-3-319-03952-7_43.
  * [5] C. Snodgrass and G. H. Jones, “The European Space Agency’s comet interceptor lies in wait,” Nature communications, Vol. 10, No. 1, 2019, pp. 1–4.
  * [6] T. Engelhardt, R. Jedicke, P. Vereš, A. Fitzsimmons, L. Denneau, E. Beshore, and B. Meinke, “An observational upper limit on the interstellar number density of asteroids and comets,” The Astronomical Journal, Vol. 153, No. 3, 2017, p. 133.
  * [7] R. Wainscoat, K. Chambers, E. Lilly, R. Weryk, S. Chastel, L. Denneau, and M. Micheli, “The Pan-STARRS search for Near Earth Objects,” Asteroids: New Observations, New Models (S. R. Chesley, A. Morbidelli, R. Jedicke, and D. Farnocchia, eds.), Vol. 318, Jan. 2016, pp. 293–298, 10.1017/S1743921315009187.
  * [8] E. Christensen, B. Africano, G. Farneth, C. Fuls, A. Gibbs, A. Grauer, H. Groeller, R. Kowalski, S. Larson, G. Leonard, T. Prune, D. Rankin, R. Seaman, and F. Shelly, “The Catalina Sky Survey: 20 Years of NEO Observations,” EPSC-DPS Joint Meeting 2019, Vol. 2019, Sept. 2019, pp. EPSC–DPS2019–1912.
  * [9] S. Keys, P. Vereš, M. J. Payne, M. J. Holman, R. Jedicke, G. V. Williams, T. Spahr, D. J. Asher, and C. Hergenrother, “The digest2 NEO Classification Code,” Publications of the Astronomical Society of the Pacific, Vol. 131, June 2019, p. 064501, 10.1088/1538-3873/ab1157.
  * [10] Ž. Ivezić, S. M. Kahn, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, D. Alonso, Y. AlSayyad, S. F. Anderson, J. Andrew, et al., “LSST: from science drivers to reference design and anticipated data products,” The Astrophysical Journal, Vol. 873, No. 2, 2019, p. 111.
  * [11] T. Hoffman, Y. Liu, M. Lysek, A. Murray, E. Nilsen, P. Peddada, M. Rokey, M. Vaquero, A. Mainzer, A. Wong, et al., “Near-Earth Object Surveyor Overview,” 2022 IEEE Aerospace Conference (AERO), IEEE, 2022, pp. 1–16.
  * [12] B. P. S. Donitz, D. Mages, D. Landau, J. V. Maydan, J. Castillo-Rogez, and D. Farnocchia, “Interstellar Object Encounter Trade Space Exploration,” 2022 IEEE Aerospace Conference (AERO), 2022, pp. 1–15, 10.1109/AERO53065.2022.9843567.
  * [13] W. Owen, T. Duxbury, C. Acton, S. Synnott, J. Riedel, and S. Bhaskaran, “A brief history of optical navigation at JPL,” 31st Annual AAS Guidance and Control Conference. Breckenridge, CO, USA, 2008.
  * [14] E. H. Maize, “Linear statistical analysis of maneuver optimization strategies,” Proceedings of the American Astronautical Society and American Insititute of Aeronautics and Astronautics Astrodynamics Specialist Conference, Kalispell, MT, 1987.
  * [15] W. Kizner, “A method of describing miss distances for lunar and interplanetary trajectories,” Planetary and Space Science, Vol. 7, 1961, pp. 125–131.
  * [16] P. Vereš, D. Farnocchia, S. R. Chesley, and A. B. Chamberlin, “Statistical analysis of astrometric errors for the most productive asteroid surveys,” Icarus, Vol. 296, Nov. 2017, pp. 139–149, 10.1016/j.icarus.2017.05.021.
  * [17] B. G. Marsden, Z. Sekanina, and D. K. Yeomans, “Comets and nongravitational forces. V,” Astronomical Journal, Vol. 78, Mar. 1973, p. 211, 10.1086/111402.
  * [18] M. Micheli, D. Farnocchia, K. J. Meech, M. W. Buie, O. R. Hainaut, D. Prialnik, N. Schörghofer, H. A. Weaver, P. W. Chodas, J. T. Kleyna, R. Weryk, R. J. Wainscoat, H. Ebeling, J. V. Keane, K. C. Chambers, D. Koschny, and A. E. Petropoulos, “Non-gravitational acceleration in the trajectory of 1I/2017 U1 (’Oumuamua),” Nature, Vol. 559, June 2018, pp. 223–226, 10.1038/s41586-018-0254-4.
  * [19] M.-T. Hui, Q.-Z. Ye, D. Föhring, D. Hung, and D. J. Tholen, “Physical Characterization of Interstellar Comet 2I/2019 Q4 (Borisov),” Astronomical Journal, Vol. 160, Aug. 2020, p. 92, 10.3847/1538-3881/ab9df8.
  * [20] D. Farnocchia, S. R. Chesley, P. W. Chodas, P. Tricarico, M. S. P. Kelley, and T. L. Farnham, “Trajectory Analysis for the Nucleus and Dust of Comet C/2013 A1 (Siding Spring),” The Astrophysical Journal, Vol. 790, Aug. 2014, p. 114, 10.1088/0004-637X/790/2/114.
  * [21] B. Schutz, B. Tapley, and G. H. Born, Statistical orbit determination.  Elsevier, 2004.
  * [22] I. Bertini, F. La Forgia, C. Tubiana, C. Güttler, M. Fulle, F. Moreno, E. Frattin, G. Kovacs, M. Pajola, H. Sierks, et al., “The scattering phase function of comet 67P/Churyumov–Gerasimenko coma as seen from the Rosetta/OSIRIS instrument,” Monthly Notices of the Royal Astronomical Society, Vol. 469, No. Suppl_2, 2017, pp. S404–S415.
  * [23] D. Mages, D. Landau, B. Donitz, and S. Bhaskaran, “Navigation evaluation for fast interstellar object flybys,” Acta Astronautica, Vol. 191, 2022, pp. 359–373.


# Separating the effects of experimental noise from inherent system
variability in voltammetry: the $[$Fe(CN)${}_{6}]^{3-/4-}$ process

Martin Robinson<EMAIL_ADDRESS>Department of Computer Science,
University of Oxford, Wolfson Building, Parks Road, Oxford, OX1 3QD, United
Kingdom.  Alexandr N Simonov  School of Chemistry, Monash University, Clayton,
Vic. 3800, Australia.  Jie Zhang  School of Chemistry, Monash University,
Clayton, Vic. 3800, Australia.  Alan Bond<EMAIL_ADDRESS>School of
Chemistry, Monash University, Clayton, Vic. 3800, Australia.  David Gavaghan
<EMAIL_ADDRESS>Department of Computer Science, University of
Oxford, Wolfson Building, Parks Road, Oxford, OX1 3QD, United Kingdom.

###### Abstract

Recently, we have introduced the use of techniques drawn from Bayesian
statistics to recover kinetic and thermodynamic parameters from voltammetric
data, and were able to show that the technique of large amplitude ac
voltammetry yielded significantly more accurate parameter values than the
equivalent dc approach. In this paper we build on this work to show that this
approach allows us, for the first time, to separate the effects of random
experimental noise and inherent system variability in voltammetric
experiments. We analyse ten repeated experimental data sets for the
$[$Fe(CN)${}_{6}]^{3-/4-}$ process, again using large-amplitude ac cyclic
voltammetry. In each of the ten cases we are able to obtain an extremely good
fit to the experimental data and obtain very narrow distributions of the
recovered parameters governing both the faradaic (the reversible formal
faradaic potential, $E_{0}$, the standard heterogeneous charge transfer rate
constant $k_{0}$, and the charge transfer coefficient $\alpha$) and non-
faradaic terms (uncompensated resistance, $R_{u}$, and double layer
capacitance, $C_{dl}$). We then employ hierarchical Bayesian methods to
recover the underlying “hyperdistribution” of the faradaic and non-faradaic
parameters, showing that in general the variation between the experimental
data sets is significantly greater than suggested by individual experiments,
except for $\alpha$ where the inter-experiment variation was relatively minor.
Correlations between pairs of parameters are provided, and for example, reveal
a weak link between $k_{0}$ and $C_{dl}$ (surface activity of a glassy carbon
electrode surface). Finally, we discuss the implications of our findings for
voltammetric experiments more generally.

##  1 Introduction

In a previous paper 1, we described the use of Bayesian inference for
quantitative comparison of voltammetric methods for investigating electrode
kinetics. We illustrated the utility of the approach by comparing the
information content in both dc and ac voltammetry at a planar electrode for
the case of a quasi-reversible one electron reaction mechanism. Using both
synthetic and experimental data, we were able to demonstrate that realistic
levels of purely random experimental (Gaussian) noise have a relatively minor
affect on the inverse problem of recovering both the faradaic (the reversible
formal potential, $E_{0}$, the standard heterogeneous charge transfer rate
constant $k_{0}$, and the charge transfer coefficient $\alpha$) and non-
faradaic (uncompensated resistance, $R_{u}$, and double layer capacitance,
$C_{dl}$) parameters that govern this reaction mechanism. We also demonstrated
the clear advantages in terms of accuracy of parameter recovery of the large
amplitude ac approach. With this ability of being able to recover parameter
values reliably from a single experimental data set in place, we are now in a
position to go on to investigate _system level_ variability i.e. if we repeat
the same experiment multiple times and implement our Bayesian parameter
recovery procedure for each data set independently, how reproducibly do we
recover the governing parameters?

To approach this problem, we again return to the “pathological”
$[$Fe(CN)${}_{6}]^{3-/4-}$ process. It is well known that for this system the
kinetic parameters reported are highly variable even when using apparently
identical electrodes and conditions (see 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14 and references cited therein). Using our new approach we are able to
show that this difficulty is not due to the impact of random experimental
noise, since for each individual data set we are able to fit the mathematical
model to the experimental data extremely accurately. However, we are able to
demonstrate that the recovered values from each individual data set vary
significantly from one another, i.e. the system itself varies between
experimental runs, and so the recovered parameter values are extremely
sensitive to the precise experimental conditions pertaining on that particular
run. Since we have already demonstrated its advantages over the dc approach,
we restrict ourselves to the ac case in this paper.

##  2 Methods

###  2.1 Experimental methods

Details of the experimental data sets used in this paper have been given
previously in 2, 15. In summary, large amplitude ac voltammetry was performed
in a standard three-electrode cell, using a glassy carbon macrodisk (diameter
3 mm) working electrode. All potentials are reported versus an Ag/AgCl/KCl(3
M) reference electrode (hereinafter Ag/AgCl). The frequencies, amplitude and
scan rate were 9.02 Hz, 0.080 V and 0.894 $\text{Vs}^{-1}$ respectively and
data were collected over the potential range of 0.5 to -0.1 V versus Ag/AgCl.
The surface area of the electrode was estimated as $0.070\;\mathrm{cm}^{2}$.
The value of the diffusion coefficient, $D$, of
$\left[\mathrm{Fe(CN)}_{6}\right]^{3-}$ was found to be $7.2\times
10^{-6}\;\mathrm{cm}^{2}\;\mathrm{s}^{-1}$, as described in 15. In this paper
we make use of the ten repeated data sets for ac voltammetry taken from 15 for
the reduction of aqueous 1.0 mM $[$Fe(CN)${}_{6}]^{3-}$ in 3 M KCl aqueous
electrolyte solution. In order to reduce the amount of computation required
for the parameter inference, we use a moving average window of length 21 to
reduce the number of experimental data points within each data set to about
25,000 data points.

###  2.2 Mathematical Modelling

The details of the mathematical and computational modelling approach that we
have taken in this paper were given previously in 15, 1. In summary, our
chosen experimental system is modelled as a quasi-reversible reaction

$\displaystyle A+e^{-}\cee{&<=>[E_{0},k_{0},\alpha]}B$

where species $A$ and $B$ are in solution, and $E_{0}$, $k_{0}$, and $\alpha$
are the reversible formal potential, standard heterogeneous charge transfer
rate constant at $E_{0}$ and the charge transfer coefficient, respectively. We
assume that the Butler-Volmer formalism applies to the electron transfer
process 16, 17, 18. We also assume that both convection and migration can be
neglected since we are using a macrodisk stationary electrode and an excess of
supporting electrolyte, respectively. Since we assume equal diffusion
coefficients for each species $A$ and $B$ ($D_{A}=D_{B}=D$) we need to solve
for the concentration of only one of the species (i.e. the concentrations
$c_{A}$, $c_{B}$ of species $A$ and $B$ satisfy $c_{A}=c_{\infty}-c_{B}$,
where $c_{\infty}$ is the bulk concentration of species $A$) and we choose to
solve for species $A$. We can then use Fick’s second law to model the
variation with time of species $A$ via

$\displaystyle\frac{\partial c_{A}}{\partial t}$
$\displaystyle=D\frac{\partial^{2}c_{A}}{\partial x^{2}},$ (1)

where $x$ is distance from the electrode surface and $t$ is time. The initial
and boundary conditions are

$\displaystyle c_{A}(x,0)$ $\displaystyle=c_{\infty}$ $\displaystyle c_{A}$
$\displaystyle\rightarrow c_{\infty},\quad\text{as}\quad
x\rightarrow\infty,\quad t>0.$ (2)

At the electrode surface, $x=0$, for $t>0$, we have the conservation and flux
conditions

$\displaystyle D\frac{\partial c_{A}}{\partial x}=\frac{I_{f}}{FS},$ (3)

along with the Butler-Volmer condition

$\displaystyle D\frac{\partial c_{A}}{\partial x}=\text{ }$ $\displaystyle
k_{0}\left[(c_{\infty}-c_{A})\exp\left((1-\alpha)\frac{F}{RT}(E_{\mbox{\tiny
eff}}(t)-E_{0})\right)\right.\left.-c_{A}\exp\left(-\alpha\frac{F}{RT}(E_{\mbox{\tiny
eff}}(t)-E_{0})\right)\right].$ (4)

Here, $I_{f}$ is the faradaic current, $S$ is the electrode area, and
$E_{\mbox{\tiny eff}}(t)$ is the effective applied potential (defined below).

We complete the model by defining $E_{\mbox{\tiny{app}}}(t)$ to be the applied
potential, then for the case of an ac voltammetry ramp we have

$\displaystyle
E_{\mbox{\tiny{app}}}(t)=E_{\mbox{\tiny{start}}}\left\\{\begin{array}[]{ll}+vt+\Delta
E\sin{(\omega t)},\qquad 0\leq t\leq t_{\mbox{\tiny reverse}},\\\
-vt+2vt_{\mbox{\tiny reverse}}+\Delta E\sin{(\omega t)},\qquad t_{\mbox{\tiny
reverse}}\leq t\leq 2t_{\mbox{\tiny reverse}}\end{array}\right.$ (7)

where $v$ is the sweep rate, $E_{\mbox{\tiny{start}}}$ is the initial
potential, $t_{\mbox{\tiny reverse}}$ is the time of switching from the
forward to the reverse sweep in cyclic voltammetry, $\omega$ is the radial
frequency and $\Delta E$ is the amplitude of the sine wave. The effective
applied potential can now be defined as

$\displaystyle
E_{\mbox{\tiny{eff}}}(t)=E_{\mbox{\tiny{app}}}(t)-E_{\mbox{\tiny{drop}}}=E_{\mbox{\tiny{app}}}(t)-I_{\mbox{\tiny{tot}}}R_{u}$

where $E_{\mbox{\tiny{drop}}}$ models the effect of uncompensated resistance,
$R_{u}$. $I_{\mbox{\tiny{tot}}}$ is the total (measured) current, and combines
the faradaic current and the background capacitive current, $I_{c}$, which can
be modelled as

$\displaystyle I_{c}$ $\displaystyle=C_{dl}\frac{dE_{\mbox{\tiny{eff}}}}{dt},$
(8)

where $C_{dl}$ is the double layer capacitance (assumed constant in this
work), and then

$\displaystyle I_{\mbox{\tiny{tot}}}=I_{f}+I_{c}.$ (9)

Equations 1–9 are non-dimensionalised as described previously 15, 1. The
resulting non-dimensional system of equations is solved using an implicit
finite difference method with an exponentially expanding grid (again, as
described previously 19). We can now see mathematically that the reaction
mechanism in Eq. 2.2 is governed by five parameters
$(E_{0},k_{0},\alpha,C_{dl},R_{u})$, and we will collectively denote these
parameters by the vector $\boldsymbol{\theta}$. The inverse problem that we
wish to solve can be defined as finding the best possible approximation to
$\boldsymbol{\theta}$ given our measured experimental output trace of the
current $I_{\mbox{\small{\tiny{tot}}}}^{\mbox{\tiny{data}}}$ versus potential.

##  3 Parameter Recovery Methods

In a recent paper we described in detail how methods of Bayesian inference can
be used to solve the inverse problem of parameter recovery 1 from voltammetric
data. We illustrated how these methods yield not only a point estimate for
each parameter of interest, but also a measure of our confidence in that
estimate. Full details of the approach that we adopted, including the
algorithms used can be found in 1, so that here we simply give a brief outline
of the methods.

###  3.1 Bayesian Inference

We denote by $\mathbf{y}=(y_{1},\dots,y_{T})$ the experimental data trace,
that is, the total measured current $I_{\mbox{\tiny{tot,
t}}}^{\mbox{\tiny{data}}}$ at each time point $t$. In our previous work (see
Figure 4 and Table S2 of 15) we demonstrated that, to a very good
approximation, experimental measurements can be assumed to be subject to
normally distributed random noise, typically with zero mean and some standard
deviation which we will denote by $\sigma$ (we showed that typical values of
the standard deviation of the experimental noise are in the range 1 to 2% of
the peak current). Our mathematical model of the system, Eqs 1 to 9 above,
assumes that the observed experimental data $\mathbf{y}$ is a function of the
parameters of interest $\boldsymbol{\theta}$,

$\boldsymbol{\theta}=(E_{0},k_{0},\alpha,C_{dl},R_{u}).$

We now assume further that the parameters $\boldsymbol{\theta}$ are themselves
drawn from a probability distribution. We can then frame our inverse problem
as trying to find this probability distribution for $\boldsymbol{\theta}$
given the observed values of the data $\mathbf{y}$, and denote this
probability distribution as $P(\boldsymbol{\theta}|\mathbf{y})$ (the vertical
line indicates that the values of $\mathbf{y}$ are given). In Bayesian
inference, $P(\boldsymbol{\theta}|\mathbf{y})$ is termed the posterior
probability density or posterior distribution; this is the distribution that
we want to approximate.

We now make use of Bayes’ rule which states

$P(\boldsymbol{\theta}|\mathbf{y})=\frac{P(\mathbf{y}|\boldsymbol{\theta})P(\boldsymbol{\theta})}{P(\mathbf{y})},$
(10)

where $P(\boldsymbol{\theta})$ is called the prior distribution of
$\boldsymbol{\theta}$ and is chosen to capture any prior knowledge we have of
$\boldsymbol{\theta}$ before any experimental observation. The distribution
$P(\mathbf{y}|\boldsymbol{\theta})$ is the probability density of the
experimental data $\mathbf{y}$ given a model parameterised with parameters
$\boldsymbol{\theta}$, and is termed the likelihood of the data; assuming a
known distribution of the error in the data this likelihood can be calculated.
$P(\mathbf{y})$ is a normalising term (which is the integral of all possible
densities
$P(\mathbf{y},\boldsymbol{\theta})=P(\mathbf{y}|\boldsymbol{\theta})P(\boldsymbol{\theta})$
over all values of $\boldsymbol{\theta}$), and ensures that the posterior
density $P(\boldsymbol{\theta}|\mathbf{y})$ integrates to $1$. In practice,
the calculation of this normalising term (which can be very computationally
expensive) is avoided by considering ratios of the likelihood (see below).

###  3.2 Calculating the likelihood

Writing the likelihood as

$L(\boldsymbol{\theta}|\mathbf{y})=P(\mathbf{y}|\boldsymbol{\theta}),$ (11)

we can re-arrange Bayes’ rule in Eq. 10 to give

$P(\boldsymbol{\theta}|\mathbf{y})\propto
P(\boldsymbol{\theta})L(\boldsymbol{\theta}|\mathbf{y}).$ (12)

Since we assume that the errors are independent at each time point, the
conditional probability density of observing the experimental trace from
$t=1,\ldots,T$ given $\boldsymbol{\theta}$ is simply the product of the
probability density functions at each time point, that is, the likelihood is
given by

$\displaystyle L(\boldsymbol{\theta}|\mathbf{y})$
$\displaystyle=\prod_{t=1}^{T}P(y_{t}|\boldsymbol{\theta})=\prod_{t=1}^{T}\mathcal{N}(y_{t}|f_{t}(\boldsymbol{\theta}),\sigma^{2})$
(13)
$\displaystyle=\prod_{t=0}^{T}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{\left(y_{t}-f_{t}(\boldsymbol{\theta})\right)^{2}}{2\sigma^{2}}\right),$
(14)

using the assumption that the experimental noise is normally distributed with
a mean zero and variance of $\sigma^{2}$. For notational simplicity we have
set
$f_{t}(\boldsymbol{\theta})=I_{\mbox{\tiny{tot}},t}^{\mbox{\tiny{model}}}$.

Since, in most experimental situations in electrochemistry, a priori we will
have only a rough idea of what the values of the parameters are likely to be,
we assume an “uninformative” prior and use a uniform distribution for each
parameter across a suitably wide range, that is

$P(\boldsymbol{\theta})=\left\\{\begin{array}[]{l}c,\>\\{\boldsymbol{\theta}\\}\text{\>in\>some\>suitably\>chosen\>5-dimensional\>
hypercube,}\\\ 0,\text{\>otherwise},\end{array}\right.$ (17)

where $c$ is a non-zero finite normalizing constant. The bounds for this
hypercube were set to

$\displaystyle E_{reverse}+0.1\delta E$ $\displaystyle\leq E_{0}\leq
E_{start}-0.1\delta E,$ $\displaystyle 0$ $\displaystyle\leq k_{0}\leq 1\text{
cm s}^{-1},$ $\displaystyle 0.4$ $\displaystyle\leq\alpha\leq 0.6,$
$\displaystyle 0$ $\displaystyle\leq C_{dl}\leq 200\text{ }\mu\text{F
cm}^{-2},$ $\displaystyle 0$ $\displaystyle\leq R_{u}\leq 80\text{ }\Omega,$

where $\delta E=E_{start}-E_{reverse}$.

Note that this prior is only used for the single level Markov Chain Monte
Carlo (MCMC) algorithm; in the hierarchical MCMC algorithm, this is replaced
by a multivariate normal, as described below. However, these bounds are still
used in the hierarchical MCMC algorithm to prevent the lower level samplers
from accepting any samples that lie outside these bounds.

####  3.2.1 Markov Chain Monte Carlo parameter inference

To obtain a sample from the posterior distribution
$P(\boldsymbol{\theta}|\mathbf{y})$ we make use of the Markov Chain Monte
Carlo method. In outline, this involves finding an approximation to the
posterior distribution $P(\boldsymbol{\theta}|\mathbf{y})$ by drawing a finite
(but sufficiently large to be accurate) number of samples from this
distribution. To do this we simulate a Markov Chain whose limiting
distribution is the required posterior distribution using an efficient
implementation of the Metropolis-Hastings algorithm 20, within which candidate
parameter sets are proposed from a proposal distribution
$q(\boldsymbol{\theta_{cand}}|\boldsymbol{\theta_{i}})$ which depends only on
the previously accepted parameter set $\boldsymbol{\theta_{i}}$; we take
$q(\boldsymbol{\theta_{cand}}|\boldsymbol{\theta_{i}})$ to be a multivariate
normal distribution. If the proposed parameter set contains any parameters
outside the range of the prior, then the parameter set is assigned an
acceptance probability of $0$, i.e. it is rejected, and the previously
accepted parameter set is added to the Markov chain — that is,
$\boldsymbol{\theta_{i+1}}=\boldsymbol{\theta_{i}}$. Otherwise, we compare
$\boldsymbol{\theta_{cand}}$ to the current parameter set
$\boldsymbol{\theta_{i}}$ by calculating the ratio of the posteriors of the
two parameter sets. If the candidate parameter set has a greater posterior
density value than the existing parameter set then it will be added to the
Markov chain, that is $\boldsymbol{\theta_{i+1}}=\boldsymbol{\theta_{cand}}$.
Otherwise, (making use of Eq. 12) the proposed parameter is accepted with
probability, $r$, given by

$r=\min\Bigg{\\{}\frac{P(\boldsymbol{\theta_{cand}})L(\boldsymbol{\theta_{cand}}|\mathbf{y})}{P(\boldsymbol{\theta_{i}})L(\boldsymbol{\theta_{i}}|\mathbf{y})},1\Bigg{\\}}.$
(18)

If the proposed parameter set is rejected (with probability $1-r$), then the
previously accepted parameter set is again added to the Markov chain — that
is, $\boldsymbol{\theta_{i+1}}=\boldsymbol{\theta_{i}}$.

Since the likelihoods for large samples are extremely small, in practice we
work with the natural log of the likelihood which reduces to

$l(\boldsymbol{\theta}|\mathbf{y})=-T\log(\sigma)-\frac{1}{2\sigma^{2}}\sum_{t=1}^{T}(y_{t}-f_{t}(\boldsymbol{\theta}))^{2},$
(19)

where terms which are constant in $\boldsymbol{\theta}$ have been removed
(since these will cancel on taking the difference of log-likelihoods in the
Metropolis-Hastings algorithm). More comprehensive descriptions of the theory
of MCMC can be found in the statistics literature (see for example 21).

###  3.3 Practical implementation of the Metropolis-Hastings Algorithm

In practice, we make use of an adaptive covariance matrix version of the
Metropolis-Hastings algorithm which helps identify the directions in parameter
space which have the highest likelihood values 22. At each iteration of the
algorithm, the covariance matrix of the multivariate normal distribution is
updated and a scalar value is also updated to define the width of the
distribution. In the results presented in this paper, we run our MCMC chains
for 10,000 samples and discard the first 5,000 samples as ‘burn in’ (see 21).
To ensure efficiency in our Monte-Carlo sampling, we first find the location
of the optimum, that is, the maximum likelihood estimate of
$\boldsymbol{\theta}$ using a standard global minimisation algorithm (we use
the cma-es algorithm 23), which we use as a seed point for the MCMC algorithm
as described in 1. In the results section, these samples are shown as
histograms which illustrate the nature of the posterior distribution.

###  3.4 Hierarchical Bayesian Inference

In the results section, the MCMC algorithm described above is used to estimate
the parameter values and their posterior distributions from each of the ten
repeat runs of the ac voltammetry experiment for the reduction of aqueous 1 mM
$[$Fe(CN)${}_{6}]^{3-}$, as described in the Experimental methods section.
This will allow us to show the degree of variability in the recovered values
of the parameters across these ten data sets. This will in turn allow us to
postulate that on each repeat of the experiment the parameters
$\boldsymbol{\theta}=\left\\{E_{0},k_{0},\alpha,C_{dl},R_{u}\right\\}$ are
themselves drawn from a multivariate normal distribution with mean hyper-
parameters
$\boldsymbol{\mu}=(\hat{E_{0}},\hat{k_{0}},\hat{\alpha},\hat{C_{dl}},\hat{R_{u}})$
and a covariance hyper-parameter matrix $\boldsymbol{\Sigma}$. Our aim now is
to sample from the distributions of $\boldsymbol{\mu}$ and
$\boldsymbol{\Sigma}$ to enable us to quantify how $\boldsymbol{\theta}$
varies between different experiments.

We first write the posterior distribution of all of our parameters given the
data, taking into account our new hierarchical model structure

$\displaystyle
P(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{n},|\mathbf{y}_{1},...,\mathbf{y}_{n})$
$\displaystyle\propto
P(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{n})\prod_{i=0}^{n}P(\mathbf{y}_{i}|\boldsymbol{\theta}_{i}),$
(20)
$\displaystyle=P(\boldsymbol{\mu},\boldsymbol{\Sigma})\prod_{i=0}^{n}P(\boldsymbol{\theta}_{i}|\boldsymbol{\mu},\boldsymbol{\Sigma})\prod_{i=0}^{n}P(\mathbf{y}_{i}|\boldsymbol{\theta}_{i}),$
(21)

where $\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{n}$ are all the bottom
level parameters for the $n$ different experiments, and
$\mathbf{y}_{1},...,\mathbf{y}_{n}$ are the corresponding measurements.

We choose a multivariate normal distribution for the bottom level parameters
$\boldsymbol{\theta}$

$P(\boldsymbol{\theta}_{i}|\boldsymbol{\mu},\boldsymbol{\Sigma})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}).$
(22)

It now remains to choose a suitable hyper-prior
$P(\boldsymbol{\mu},\boldsymbol{\Sigma})$, which for ease of computation is
generally taken to be a normal-inverse-Wishart distribution (see for example
24)

$P(\boldsymbol{\mu},\boldsymbol{\Sigma})=\mathcal{NIW}(\boldsymbol{\mu}_{0},\kappa_{0},\nu_{0},\boldsymbol{\Psi}).$
(23)

The distributions in Eqs. 22 and 23 are chosen to be conjugate, so that the
conditional distribution of the hyper-parameters can be analytically derived
as 25

$P(\boldsymbol{\mu},\boldsymbol{\Sigma}|\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{n})=\mathcal{NIW}(\frac{\kappa_{0}\boldsymbol{\mu}_{0}+n\hat{\boldsymbol{\theta}}}{\kappa_{0}+n},\kappa_{0}+n,\nu_{0}+n,\boldsymbol{\Psi}+\mathbf{C}+\frac{\kappa_{0}n}{\kappa_{0}+n}(\hat{\boldsymbol{\theta}}-\boldsymbol{\mu}_{0})(\hat{\boldsymbol{\theta}}-\boldsymbol{\mu}_{0})^{T}),$
(24)

where $\hat{\boldsymbol{\theta}}$ and $\mathbf{C}$ are the sample mean and
covariance of the bottom level parameters

$\displaystyle\hat{\boldsymbol{\theta}}$
$\displaystyle=\frac{1}{n}\sum_{i=0}^{n}\boldsymbol{\theta}_{i},$ (25)
$\displaystyle\mathbf{C}$
$\displaystyle=\sum_{i=0}^{n}(\boldsymbol{\theta}_{i}-\hat{\boldsymbol{\theta}})(\boldsymbol{\theta}_{i}-\hat{\boldsymbol{\theta}})^{T}.$
(26)

We can now use our original adaptive MCMC algorithm to sample from the bottom
level parameters, combined with classical Gibbs sampling and Eq. 24 to sample
from the hyper-parameters (see Algorithm 1 for details). For the hyperprior
parameters we use $\kappa_{0}=0$ and $\nu_{0}=1$. We set
$\boldsymbol{\mu}_{0}$ to be the centre point of the 5-dimensional hypercube
in Eq. 17, and $\boldsymbol{\Psi}$ as a diagonal matrix with the standard
deviation of each parameter set to half the width of this same hypercube. To
prevent the bottom level MCMC chains from exploring unphysical parameter
regimes, we automatically reject any proposed point which lies outside of the
hypercube.

Algorithm 1 Metropolis within Gibbs

$s=0$

$\boldsymbol{\mu}_{s},\boldsymbol{\Sigma}_{s}=\mbox{SampleFrom}\mathcal{NIW}(\boldsymbol{\mu}_{0},\kappa_{0},\nu_{0},\boldsymbol{\Psi})$
{using Eq. 23}

repeat

for $i=1$ to $n$ do

$\boldsymbol{\theta}_{i}=\mbox{AdaptiveMCMCStep}(\boldsymbol{\mu}_{s},\boldsymbol{\Sigma}_{s})$
{see Eq. 18 and enclosing section}

end for

$s=s+1$

$\boldsymbol{\mu}_{s},\boldsymbol{\Sigma}_{s}=\mbox{SampleFrom}\mathcal{NIW}(\boldsymbol{\mu}_{0},\kappa_{0},\nu_{0},\boldsymbol{\Psi},\boldsymbol{\theta}_{1},...,\boldsymbol{\theta}_{n})$
{using Eq. 24 & 25}

until finished sampling

###  3.5 Generating synthetic data as a test case

To test our inference procedure and algorithms we make use of “synthetic” test
data i.e. we generate ten synthetic data sets by solving Equations 1 to 9 for
a randomly chosen set of values of
$\boldsymbol{\theta}=(E_{0},k_{0},\alpha,C_{dl},R_{u})$, drawn from a
multivariate normal distribution with mean $\boldsymbol{\theta}_{\mbox{\tiny
true}}=(7.27,2.01,0.53,3.70\times 10^{-3},1.06\times 10^{-2})$ (in non-
dimensional units), and standard deviation $\boldsymbol{\sigma}_{\mbox{\tiny
true}}=(0.06,0.7,0.005,0.7\times 10^{-3},0.3\times 10^{-2})$. We then add
randomly generated Gaussian noise at each time point with zero mean and
standard deviation of 0.3% of the maximum current (chosen to match the
experimental noise level).

##  4 Results

###  4.1 Synthetic data

To verify that we have implemented our hierarchical Bayes algorithms correctly
we briefly describe the results of using synthetic data to test our inference
procedure. Figure 1 shows the sampled distributions of the hyper-parameters.
The dashed lines indicate the expected peak of each distribution, as
calculated from the sample mean and variance of the ten true values of
$\boldsymbol{\theta}_{\mbox{\tiny true}}$. As can be seen the hierarchical
MCMC algorithm obtains the correct peak for all the five different mean and
variance hyper-parameters.

Figure 1: Histograms of sampled hyper-parameter distributions
$\boldsymbol{\mu}$ (left column) and the variances of $\boldsymbol{\Sigma}$
(right column), generated by using Algorithm 1 on 10 sets of synthesised ac
voltammetric data. The dashed lines indicate the sample mean and sample
standard deviation of the true parameters used to generate the synthetic
datasets, and since this aligns correctly with the maximum likelihood point of
the histograms we can be confident that our implementation of Algorithm 1 is
sampling the correct posterior distribution for $\boldsymbol{\mu}$ and
$\boldsymbol{\Sigma}$.

###  4.2 Experimental data

Figure 2 shows an example of the initial fitting process for experimental data
set 1. Once the best fit for each data set was obtained using maximum
likelihood estimation and the cma-es optimisation algorithm, each of the
bottom level MCMC algorithms was initialised at these points. Then Algorithm 1
was used to generate 5,000 samples (10,000 total samples, of which the first
5,000 was discarded as burn-in) of the parameters $\boldsymbol{\theta}_{i}$
and hyper-parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$.

Figure 2: Comparison between the simulated current trace and the experimental
data for the reduction of aqueous 1 mM $[$Fe(CN)${}_{6}]^{3-}$ (dataset 1 from
reference 15). Simulation parameters were obtained by maximum likelihood
estimate, and their values in dimensional units were $E_{0}=0.214$ V,
$k_{0}=0.010$ $\text{ cm s}^{-1}$, $\alpha=0.528$, $C_{dl}=16.9$ $\mu\text{F
}\text{cm}^{-2}$, $R_{u}=0.00$ $\Omega$.

Figure 3 (left column) shows histograms of the samples of
$\boldsymbol{\theta}_{i}$ taken from the bottom level samplers in Algorithm 1.
On the same axis is drawn the posterior predictive distribution for each
variable, calculated by summing the individual Gaussian distributions
described by each sample of the hyper-parameters $\boldsymbol{\mu}$ and
$\boldsymbol{\Sigma}$. The posterior predictive distributions describes the
distribution of each parameter that would be expected from another repeat of
the experiment, given the results of the 10 already observed experiments. As
can be seen, this distribution covers the width of all 10 bottom level samples
of $\boldsymbol{\theta}_{i}$, and illustrates the significantly greater
variation in each parameter that is expected between subsequent experiments.

It is interesting to compare the samples obtained from the hierarchical model
to those taken using the original non-hierarchical model (i.e. just running
ten independent MCMC chains on the ten different data sets), and this is shown
in Figure 3 (right column). As can be seen, the result in this case are almost
identical to the hierarchical model, giving us confidence that we are
capturing the distributions of $\boldsymbol{\theta}_{i}$ correctly in each
case.

Figure 3: Analysis of ten independent ac voltammetric experiments for the
reduction of aqueous 1 mM $[$Fe(CN)${}_{6}]^{3-}$ using the hierarchical MCMC
algorithm 1. The samples obtained from the 10 lower level adaptive MCMC
samplers (i.e. $\boldsymbol{\theta}_{i}$) are shown as histograms with the
axis label $P_{i}(\cdot)$, and the different chains from $i=1,..,10$ are shown
with different colours. The left column plots show the samples obtained using
the hierarchical model, and these also show the posterior predictive
distribution $P(\cdot)$ for each variable, calculated by summing the Gaussian
distributions described by the samples of the hyper-parameters
$\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ (see Figure 4 for histograms of
these hyper-parameter samples). For comparison, the right column plots show
the samples of $\boldsymbol{\theta}_{i}$ with no hierarchical model (and thus
no hyper-parameter samples). This is identical to the analysis done in our
previous paper 1.

The chief benefit of the hierarchical model is that it allows us to quantify
(with the hyper-parameters) the variability of the parameters _between_
different experiments. Figure 4 shows the histograms for the hyper-parameters
samples of $\boldsymbol{\mu}$ (left) and $\boldsymbol{\Sigma}$ (right). Also
shown as vertical dashed lines are the the sample mean and variance of the
concatenated ten $\boldsymbol{\theta}_{i}$ chains. The most obvious feature of
these plots is that these hyper-parameter distributions clearly show a much
wider confidence interval as compared with the individual
$\boldsymbol{\theta}_{i}$ distributions. While the value of each parameter for
each _individual_ experiment is known with high accuracy, once the variability
between experiments is taken into account we see that the possible parameter
range is significantly broadened.

Figure 4: Histograms of hyper-parameter samples of $\boldsymbol{\mu}$ (left)
and $\boldsymbol{\Sigma}$ (right) obtained by applying the hierarchical MCMC
algorithm 1 to the ten experimental ac voltammetric datasets for the reduction
of aqueous 1 mM $[$Fe(CN)${}_{6}]^{3-}$. Note that only the diagonal elements
of $\boldsymbol{\Sigma}$ (i.e. the variances $\sigma^{2}$) are shown. For
comparison, the dashed lines show the sample mean and variance taken across
all ten lower level MCMC samples of $\boldsymbol{\theta}_{i}$.

##  5 Discussion

The Bayesian data analysis approach introduced in this paper provides access
to fundamentally new knowledge that assists in elucidating nuances that have
contributed to the highly non-reproducible electrode kinetic data published
for the “pathologically variable” $[$Fe(CN)${}_{6}]^{3-/4-}$ process.
Traditional heuristic and data analysis optimisation methods produce only
single point values for a limited set of parameters and do not quantify the
system variability, which is crucial information. Using a statistically based
Bayesian inference approach, we are now able to show that the difficulty in
achieving reproducibility in the voltammetry is not associated with the impact
of random noise, since for each data set we are able to fit the experimental
data extremely accurately using a model derived from use of Butler-Volmer
electron transfer kinetics, mass transport by planar diffusion, uncompensated
resistance and double layer capacitance. Thus, while substantial variation in
$k_{0}$ from about 0.002 to 0.018 cm $\text{s}^{-1}$ is evident in 10
individual experiments at a nominally identical electrode surface (Figure 3),
conformance to the quasi-reversible model is exceptionally good for each
individual experiment.

The variation in performance of the now very widely used glassy carbon
electrode was identified as a point of concern soon after the material was
introduced into electroanalytical chemistry 26. In essence, the exact nature
of the glassy carbon, and indeed other carbon based surfaces, has their origin
in the method (e.g. temperature) of manufacture, nature of pre-treatment and
its history as an electrode (see for example 13, 14, 26, 27, 28, 29, 30, 31,
32 and references cited therein). Accordingly, it seems likely that the widely
variable electrode kinetic data associated with published studies on the
$[$Fe(CN)${}_{6}]^{3-/4-}$ process (summarised in Eq. 27)

$\text{Fe(CN)}_{6}^{3-}+\text{e}^{-}\cee{<=>}\text{Fe(CN)}_{6}^{4-}\qquad(E_{0},k_{0},\alpha,C_{dl},R_{u}),$
(27)

mimics the variability of the surface state used in the different
publications.

In practice, carbon electrodes are highly heterogeneous with surface defects
and organic functional groups in abundance when used in aqueous electrolyte
media 28, 29, 30. This seems to translate into an electron transfer process
occurring at a surface consisting of microscopically small and distinctly
different regions that must be sufficiently close so that complete overlap of
diffusion layers occurs on the measurement timescale 33, presumably allowing
the entire surface to be successfully modelled by approximating the mass
transfer by planar diffusion. Thus, even though highly heterogeneous with
variable $k_{0}$ values at the microscopic level, the surface behaves as
though it is fully homogeneous, hence giving rise to measurement of an
apparently single ”averaged” $k_{0}$ value.

Electron transfer reaction mechanisms are often classified into inner and
outer sphere categories in both homogeneous chemical redox reactions that
occur in the solution phase and in heterogeneous reactions in electrochemistry
that take place across an electrode-solution interface33, 16. In a chemical
redox reaction, the homogeneous outer sphere class mechanism is characterised
by weak interactions of the reactant and product, as in Eq. 28,

$^{*}\text{Fe(CN)}_{6}^{3-}\text{(aq)}+\text{Fe(CN)}_{6}^{4-}\text{(aq)}\cee{<=>}^{*}\text{Fe(CN)}_{6}^{4-}\text{(aq)}+\text{Fe(CN)}_{6}^{3-}\text{(aq)},$
(28)

which is the homogeneous analogue of the electrochemical one of interest in
this study, while in the inner class, ligands involving bridging to a common
metal centre may be involved.

In an outer sphere electrochemical process, the plane of closest approach to
the Outer Helmholtz Plane does not allow penetration of the layer of non-
specifically adsorbed or coordinated solvent adhered to the electrode surface
by reactants. In contrast, inner sphere electron transfer processes involve
specifically adsorbed reactants 33, 16, and therefore are anticipated to
exhibit electrode kinetics that are strongly dependent on the chemical nature
of the electrode surface. Thus, the $[$Fe(CN)${}_{6}]^{3-/4-}$ process can be
designated as outer sphere under homogeneous chemical redox reaction
conditions but inner sphere under electrochemical conditions at a glassy
carbon electrode. Presumably, electron transfer at such electrodes is
accompanied by interaction with surface functional groups, facilitated by the
high negative charges associated with the $[$Fe(CN)${}_{6}]^{3-}$ reactant and
$[$Fe(CN)${}_{6}]^{4-}$ product which allow specific binding or electrostatic
attractive and repulsive interactions to accompany the electron transfer
process. On this basis, the $[$Fe(CN)${}_{6}]^{3-/4-}$ electrode kinetics are
strongly dependent on the treatment and origin of the glassy carbon.

Before attempting to identify the factors that may be most relevant to the
variation encountered in the electrode kinetics, it is worthwhile reviewing
the electrode pre-treatment regime. A nominally 3 mm diameter GC working
electrode was purchased from BAS. Prior to each set of measurements the
surface of the GC electrode was thoroughly polished with 0.3 $\mu\text{m}$
alumina powder in the form of an aqueous slurry on a wet polishing cloth
(BAS). After polishing, the electrode was repeatedly washed with high purity
water and subjected to sonication for 10 to 20 s. To provide effective removal
of any residual alumina powder, after the initial sonication, the electrode
was carefully wiped with a clean wet polishing cloth, again washed with water
and sonicated for 10 to 20 s in high purity water. Finally, the working
electrode was dried under a nitrogen stream. As far as possible, each
electrode preparation was undertaken under identical conditions; nevertheless,
the electrode kinetics differ substantially from experiment to experiment.
However, of course the history of the electrode could be important since the
electrode used for experiment ten will have a more extensively polished
electrode than that used in experiment one.

The narrow distribution of each $P_{i}(\cdot)$ in Figure 3 enables us to
conclude that the contribution to parameter variability from noise associated
with each particular experiment is minimal for all five parameters estimated.
From perusal of Figures 3 and 4, we can conclude that while the degree of
experiment-to-experiment variation changes with respect to each parameter, it
is in general significantly greater than the parameter variability from noise
alone (i.e. the width of each $P_{i}(\cdot)$). From the hyper-parameter
samples we can calculate the distribution $P(\cdot)$ of each parameter that
would be expected from another repeat experiment (shown in Figure 3, left
column). Quantitatively, the parameter mean values (with one standard
deviation values in parenthesis) as deduced from this analysis are as follows:
$E_{0}=0.215$ V vs Ag/AgCl (0.002), $k_{0}=0.010$ cm $\text{s}^{-1}$ (0.005)
$\alpha=0.521$ (0.006), $C_{dl}=20.1$ $\mu F$ $\text{cm}^{-2}$ (4.5),
$R_{u}=15.2$ $\Omega$ (21.0). Intriguingly, experiment-to-experiment variation
in $\alpha$ is very small. That is, its value lies in a very narrow range of
about 0.515 to 0.525 with a standard deviation of just 0.006. It is also
notable that while the value of $R_{u}$ covers a range from zero to about 35
$\Omega$, it is always small. This means that the more important ohmic drop
term that can distort voltammograms also is small and hence not highly
significant. The large standard deviation of 21.0 $\Omega$ with a mean value
of 15.2 $\Omega$ allows us to conclude the $R_{u}$ parameter only has a minor
impact on the ac voltammetry (highly conducting 3 M KCl aqueous electrolyte)
and hence on $E_{0}$, $k_{0}$, or $\alpha$ variability. Although the
experiment-to-experiment variability of $E_{0}$ is much greater than predicted
by $P_{i}(E_{0})$, it is still fitted to a tight regime (about 2 mV) when
compared with the total voltage scan range. Indeed, in the absence of
adsorption, $E_{0}$ is theoretically predicted to be completely independent of
electrode material or state, suggesting that the inter-experiment variability
of $E_{0}$ is due to compensating factors from other parameters. $C_{dl}$
displays some variability with values lying within the range of about 16 to 23
$\mu F$ $\text{cm}^{-2}$ which confirms that the electrode is not in an
identical surface state, and hence does have a variable level of activity for
each experiment. $k_{0}$ as noted above encompasses a wide range of about
0.002 to 0.018 cm $\text{s}^{-1}$ (Figure 3) with a mean value of 0.010 cm
$\text{s}^{-1}$ and a standard deviation of 0.005 cm $\text{s}^{-1}$. The
question that arises is whether there is a correlation between $C_{dl}$ or
surface activity and $k_{0}$.

Figure 5 shows quantitatively the correlation between pairs of parameters
using scatterplots of the mean hyper-parameter $\boldsymbol{\mu}$ samples.
Clearly, $R_{u}$ and $\alpha$ are not correlated at all with each other or any
other parameter (i.e. they give “shot gun” correlation plots). However, there
appears to be a weak correlation of $k_{0}$ and $C_{dl}$ with larger $k_{0}$
values being associated with data sets having higher $C_{dl}$ values. There
also appears to be a weak correlation between $k_{0}$ and $E_{0}$, with the
more positive $E_{0}$ values coinciding with the larger $C_{dl}$ values. One
needs to be careful not to over interpret the significance of weak
correlations, as other non-quantified parameters also can be operative.
However, a larger capacitance current caused by a larger activation of the
glassy carbon surface (more surface functionality) may be expected to increase
the rate of an inner sphere process, like the $[$Fe(CN)${}_{6}]^{3-/4-}$ one
probed in this work. The origin of the weak correlation of $k_{0}$ and $E_{0}$
is more problematical. The presence of extremely weak adsorption not
accommodated in the model is one possible explanation, but a small drift in
reference electrode potential cannot be ruled out.

Figure 5: Pairwise correlation plots for the mean hyper-parameter samples
$\boldsymbol{\mu}=(\hat{E_{0}},\hat{k_{0}},\hat{\alpha},\hat{C_{dl}},\hat{R_{u}})$
(see Algorithm 1) obtained for 10 independent ac voltammetric experiments for
the reduction of aqueous 1 mM $[$Fe(CN)${}_{6}]^{3-}$. Each dot in the scatter
plots shows a sample drawn from the posterior distribution
$P(\boldsymbol{\mu},\boldsymbol{\Sigma})$ (Eq. 23), showing the correlation
between parameter pairs. The diagonal plots show histograms of each individual
mean hyper-parameter.

##  6 Conclusion

In summary, the Bayesian inference-inspired strategy introduced in this paper
for data evaluation represents a significant advance in understanding the
contribution of different parameters to voltammetric data and their
significance in experiment-to-experiment variability at a heterogeneous
electrode in a manner that has not been possible in earlier studies. Perhaps
remarkably, each data set in the pathologically variable
$[$Fe(CN)${}_{6}]^{3-/4-}$ process conforms exceptionally well with simulated
data derived from the Butler-Volmer model of electron transfer and mass
transport by planar diffusion even though the variation in $k_{0}$ is quite
substantial. Intriguingly, $\alpha$, unlike $k_{0}$, does not reveal
significant experiment-to-experiment variation in this data analysis exercise.
This may be expected if conformance to the Butler-Volmer model is strong.

In early electrode kinetic studies, the ideal and very homogenous mercury
electrode was used. It is now evident that electrode design is becoming very
sophisticated, particularly when advances in materials science are used to
generate highly heterogeneous electrode materials 34. Thus, there is a
tendency in electrochemistry nowadays to use far more complex electrodes than
the historically important liquid mercury and pure solid metal surfaces. The
new breed of highly heterogeneous electrodes imply that data analysis
strategies in the future will also need to be more sophisticated. That is, use
of models significantly more complex than used in this study with Butler-
Volmer theory for electron transfer and mass transport by planar diffusion, as
well as Bayesian forms of data analysis, will become increasingly essential if
informative experimental versus simulation comparisons are to be reported.

## Acknowledgements

AMB, DJG and JZ would like to acknowledge the support Australian Research
Council through the award of a Discovery Grant DP170101535. JZ and ANS also
acknowledge the financial support through the ARC Centre of Excellence for
Electromaterials Science (ACES). MR gratefully acknowledges research support
from the EPSRC Cross-Disciplinary Interface Programme (EP/I017909/1). AMB also
wishes to thank the Vallee Foundation for travel support that enabled him to
spend time at the University of Oxford.

## References

  * Gavaghan et al. 2018 Gavaghan, D. J.; Cooper, J.; Daly, A. C.; Gill, C.; Gillow, K.; Robinson, M.; Simonov, A. N.; Zhang, J.; Bond, A. M. Use of Bayesian inference for parameter recovery in DC and AC Voltammetry. _ChemElectroChem_ 2018, _5_ , 917–935
  * Bond et al. 2005 Bond, A. M.; Duffy, N. W.; Guo, S.-X.; Zhang, J.; Elton, D. Changing the look of voltammetry. _Anal. Chem._ 2005, _77_ , 186 A–195 A
  * McCreery 2008 McCreery, R. L. Advanced carbon electrode materials for molecular electrochemistry. _Chemical reviews_ 2008, _108_ , 2646–2687
  * McCreery and McDermott 2012 McCreery, R. L.; McDermott, M. T. Comment on electrochemical kinetics at ordered graphite electrodes. _Analytical chemistry_ 2012, _84_ , 2602–2605
  * Patel et al. 2012 Patel, A. N.; Collignon, M. G.; O’Connell, M. A.; Hung, W. O.; McKelvey, K.; Macpherson, J. V.; Unwin, P. R. A new view of electrochemistry at highly oriented pyrolytic graphite. _Journal of the American Chemical Society_ 2012, _134_ , 20117–20130
  * Ji et al. 2006 Ji, X.; Banks, C. E.; Crossley, A.; Compton, R. G. Oxygenated Edge Plane Sites Slow the Electron Transfer of the Ferro-/Ferricyanide Redox Couple at Graphite Electrodes. _ChemPhysChem_ 2006, _7_ , 1337–1344
  * Curulli et al. 2005 Curulli, A.; Valentini, F.; Orlanducci, S.; Terranova, M. Kinetic evaluation of the standard rate constant from cyclic voltammetric data at SWCNT-modified tungsten microelectrodes. _Indian journal of chemistry. Sect. A: Inorganic, physical, theoretical & analytical_ 2005, _44_ , 956–967
  * Mundinamani and Rabinal 2014 Mundinamani, S.; Rabinal, M. Cyclic voltammetric studies on the role of electrode, electrode surface modification and electrolyte solution of an electrochemical cell. _J. Appl. Chem_ 2014, _7_ , 45–52
  * Dekanski et al. 2001 Dekanski, A.; Stevanović, J.; Stevanović, R.; Nikolić, B. Ž.; Jovanović, V. M. Glassy carbon electrodes: I. Characterization and electrochemical activation. _Carbon_ 2001, _39_ , 1195–1205
  * Gonçalves et al. 2006 Gonçalves, E. S.; Rezende, M. C.; Takahashi, M. F. K.; Ferreira, N. G. Electrochemical reversibility of reticulated vitreous carbon electrodes heat treated at different carbonization temperatures. _Materials Research_ 2006, _9_ , 147–152
  * Wightman et al. 1984 Wightman, R. M.; Deakin, M. R.; Kovach, P. M.; Kuhr, W. G.; Stutts, K. J. Methods to improve electrochemical reversibility at carbon electrodes. _Journal of the Electrochemical Society_ 1984, _131_ , 1578–1583
  * Patel et al. 2012 Patel, A. N.; Collignon, M. G.; O’Connell, M. A.; Hung, W. O.; McKelvey, K.; Macpherson, J. V.; Unwin, P. R. A new view of electrochemistry at highly oriented pyrolytic graphite. _Journal of the American Chemical Society_ 2012, _134_ , 20117–20130
  * 13 Kuwana, T. Analytical Electrochemistry: A Laboratory Manual. http://www.asdlib.org/onlineArticles/elabware/kuwanaEC_lab/ec_labmanual1.htm, Accessed: 2018-07-13
  * 14 eDAQ, Cyclic Voltammetry: Ferro/Ferricyanide, fact and fiction. https://www.edaq.com/wiki/Cyclic_Voltammetry:_Ferro/Ferricyanide,_fact_and_fiction, Accessed: 2018-07-13
  * Morris et al. 2013 Morris, G.; Simonov, A.; Mashkina, E.; Bordas, R.; Gillow, K.; Baker, R.; Gavaghan, D.; Bond, A. A comparison of fully automated methods of data analysis and computer assisted heuristic methods in an electrode kinetic study of the pathologically variable $[\mathrm{Fe(CN)}_{6}]^{(3-/4-)}$ process by AC voltammetry. _Analytical chemistry_ 2013, _85_ , 11780–7
  * Bard and Faulkner 2001 Bard, A.; Faulkner, L. _Electrochemical Methods: Fundamentals and Applications_ ; John Wiley & Sons, Inc, 2001
  * Brett and Brett 1993 Brett, C.; Brett, A., Oliveira Principles, Methods, and Applications of Electrochemistry. 1993
  * Pletcher et al. 2001 Pletcher, D.; Greff, R.; Peat, R.; Peter, L.; Robinson, J. _Instrumental methods in electrochemistry_ ; Elsevier, 2001
  * Sher et al. 2004 Sher, A.; Bond, A.; Gavaghan, D.; Harriman, K.; Feldberg, S.; Duffy, N.; Guo, S.-X.; Zhang, J. Resistance, capacitance, and electrode kinetic effects in Fourier-transformed large-amplitude sinusoidal voltammetry: emergence of powerful and intuitively obvious tools for recognition of patterns of behavior. _Analytical chemistry_ 2004, _76_ , 6214–28
  * Hastings 1970 Hastings, W. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. _Biometrika_ 1970, _57_ , 97–109
  * Gilks et al. 1996 Gilks, W.; Richardson, S.; Spiegelhalter, D. _Markov Chain Monte Carlo in Practice_ ; London: Chapman and Hall, 1996
  * Haario et al. 2001 Haario, H.; Saksman, E.; Tamminen, J. An Adaptive Metropolis Algorithm. _Bernoulli_ 2001, 223–242
  * Auger and Hansen 2005 Auger, A.; Hansen, N. A restart CMA evolution strategy with increasing population size. Evolutionary Computation, 2005. The 2005 IEEE Congress on. 2005; pp 1769–1776
  * Huang et al. 2010 Huang, Y.; Wu, H.; Acosta, E. P. Hierarchical Bayesian inference for HIV dynamic differential equation models incorporating multiple treatment factors. _Biometrical Journal_ 2010, _52_ , 470–486
  * Murphy 2007 Murphy, K. P. _Conjugate Bayesian analysis of the Gaussian distribution_ ; 2007
  * Van der Linden and Dieker 1980 Van der Linden, W.; Dieker, J. W. Glassy carbon as electrode material in electro-analytical chemistry. _Analytica Chimica Acta_ 1980, _119_ , 1–24
  * DeClements et al. 1996 DeClements, R.; Swain, G. M.; Dallas, T.; Holtz, M. W.; Herrick, R. D.; Stickney, J. L. Electrochemical and surface structural characterization of hydrogen plasma treated glassy carbon electrodes. _Langmuir_ 1996, _12_ , 6578–6586
  * Dekanski et al. 2001 Dekanski, A.; Stevanović, J.; Stevanović, R.; Nikolić, B. Ž.; Jovanović, V. M. Glassy carbon electrodes: I. Characterization and electrochemical activation. _Carbon_ 2001, _39_ , 1195–1205
  * Ilangovan and Chandrasekara Pillai 1997 Ilangovan, G.; Chandrasekara Pillai, K. Electrochemical and XPS characterization of glassy carbon electrode surface effects on the preparation of a monomeric molybdate (VI)-modified electrode. _Langmuir_ 1997, _13_ , 566–575
  * DeClements et al. 1996 DeClements, R.; Swain, G. M.; Dallas, T.; Holtz, M. W.; Herrick, R. D.; Stickney, J. L. Electrochemical and surface structural characterization of hydrogen plasma treated glassy carbon electrodes. _Langmuir_ 1996, _12_ , 6578–6586
  * Chaisiwamongkhol et al. 2017 Chaisiwamongkhol, K.; Batchelor-McAuley, C.; Compton, R. G. Amperometric micro pH measurements in oxygenated saliva. _Analyst_ 2017, _142_ , 2828–2835
  * Chaisiwamongkhol et al. 2018 Chaisiwamongkhol, K.; Batchelor-McAuley, C.; Palgrave, R. G.; Compton, R. G. Singlet oxygen and the origin of oxygen functionalities on the surface of carbon electrodes. _Angewandte Chemie International Edition_ 2018, _57_ , 6270–6273
  * Compton and Banks 2011 Compton, R.; Banks, C. _Understanding voltammetry_ ; World Scientific, 2011
  * Zhang et al. 2018 Zhang, Y.; Simonov, A. N.; Zhang, J.; Bond, A. M. Fourier transformed alternating current voltammetry in electromaterials research: direct visualisation of important underlying electron transfer processes. _Current Opinion in Electrochemistry_ 2018,

Distributions for all ten lower level adaptive MCMC samples of the reversible
formal potential $E_{0}$. The coloured histograms show the lower level
posterior for each of the 10 repeat experiments $P_{i}(E_{0})$. The solid line
shows the predictive posterior for $E_{0}$, which describes the distribution
of each parameter that would be expected from another repeat of the
experiment.