$\phi\left(T\right)<\mathcal{C}_{133}.$ (141)

However, $\phi\left(T\right)$ tends to $+\infty$ as $T$ tends to $0$, and as
$T=\log\log H\left(X,\,Y,\,Z\right)$, this happens as
$H\left(X,\,Y,\,Z\right)$ gets arbitrarily large. Thus, there is a value $B$
such that if $H\left(X,\,Y,\,Z\right)>B$, then (141) cannot hold. This gives
an upper bound for values of $H\left(X,\,Y,\,Z\right)$ such that the triple
satisfies the assumptions of the theorem. It thus follows by Northcutt’s
Theorem that there are only finitely many primitive triples
$\left(X,\,Y,\,Z\right)$ satisfying $\left(X+Y=Z\right)$ with

$S\left(X,\,Y,\,Z\right)\leq\log\log H\left(X,\,Y,\,Z\right)\frac{\log\log\log
H\left(X,\,Y,\,Z\right)}{\log\log\log\log
H\left(X,\,Y,\,Z\right)\phi\left(\log\log H\left(X,\,Y,\,Z\right)\right)}$

∎

We note that we could also directly prove that there are only finitely many
primitive integer triples $\left(X,\,Y,\,Z\right)$ satisfying $X+Y=Z$ with

$S\left(X,\,Y,\,Z\right)<c\log\log H\left(X,\,Y,\,Z\right)$

for any constant $c\in\mathbb{R},\,c>0$ using the same method of proof as
above, though this result follows from Theorem 5 as stated previously.

## References

  * [1] J. Berstel and M. Mignotte.  Deux propriétés décidables des suites récurrentes linéaires.  Bulletin de la Société Mathématique de France, 104:175–184, 1976.
  * [2] E. Bombieri and W. Gubler.  Heights in Diophantine geometry.  Number 4. Cambridge university press, 2007.
  * [3] J. Browkin.  The abc-conjecture.  In Number theory, pages 75–105. Springer, 2000.
  * [4] J. Browkin.  The abc–conjecture for algebraic numbers.  Acta Mathematica Sinica, 22(1):211–222, 2006.
  * [5] J. Buchmann.  On the computation of units and class numbers by a generalization of lagrange’s algorithm.  Journal of Number Theory, 26(1):8–30, 1987.
  * [6] Y. Bugeaud and K. Győry.  Bounds for the solutions of unit equations.  Acta Arithmetica, 74:67–80, 1996.
  * [7] N. Childress.  Class field theory.  Springer Science & Business Media, 2008.
  * [8] G. Everest, A. J. Van Der Poorten, I. Shparlinski, T. Ward, et al.  Recurrence sequences, volume 104.  American Mathematical Society Providence, RI, 2003.
  * [9] J.E. Evertse and K. Győry.  Unit equations in Diophantine number theory, volume 146.  Cambridge University Press, 2015.
  * [10] K. Győry.  On the $abc$ conjecture in algebraic number fields.  Acta Arithmetica, 133:281–295, 2008.
  * [11] K. Győry.  Bounds for the solutions of $s$-unit equations and decomposable form equations II.  Publ.Math.Debrecen, 94:507–526, 2019.
  * [12] K. Győry.  $S$-unit equations and Masser’s $abc$-conjecture in algebraic number fields.  Submitted, 2021.
  * [13] K. Győry and K. Yu.  Bounds for the solutions of $s$-unit equations and decomposable form equations.  Acta Arithmetica, 123:9–41, 2006.
  * [14] V. Halava, T. Harju, M. Hirvensalo and J. Karhumäki.  Skolem’s problem–on the border between decidability and undecidability.  Technical Report 683, Turku Centre for Computer Science, 2005.
  * [15] A. J. Harper.  Minor arcs, mean values, and restriction theory for exponential sums over smooth numbers.  Compositio Mathematica, 152(6):1121–1158, 2016.
  * [16] J. C. Lagarias and K. Soundararajan.  Smooth solutions to the $abc$ equation: the $xyz$ conjecture.  Journal de théorie des nombres de Bordeaux, 23(1):209–234, 2011\.
  * [17] E. Landau.  Neuer beweis des primzahlsatzes und beweis des primidealsatzes.  Mathematische Annalen, 56(4):645–670, 1903.
  * [18] S. Lang.  Algebraic number theory, volume 110.  Springer Science & Business Media, 2013.
  * [19] S. Le Fourn.  Tubular approaches to baker’s method for curves and varieties.  Algebra & Number Theory, 14(3):763–785, 2020.
  * [20] R. C. Mason.  Diophantine equations over function fields, volume 96.  Cambridge University Press, 1984.
  * [21] D. W. Masser.  On abc and discriminants.  Proceedings of the American Mathematical Society, 130(11):3141–3150, 2002.
  * [22] D. W. Masser.  Open problems.  In Proceedings of the symposium on Analytic Number Theory, London, 1985. Imperial College, 1985.
  * [23] S. Mochizuki.  Inter-universal teichmüller theory I, II, III, IV.  Publications of the Research Institute for Mathematical Sciences, 57(1):3–723, 2021.
  * [24] S. Natarajan and R. Thangadurai.  Pillars of Transcendental Number Theory.  Springer-Verlag, 2020.
  * [25] J. Neukirch.  Algebraic number theory, volume 322.  Springer Science & Business Media, 2013.
  * [26] J. Oesterlé.  Nouvelles approches du “théoreme” de fermat.  Astérisque, 161(162):165–186, 1988.
  * [27] A. Ostafe and I. E. Shparlinski.  On the Skolem problem and some related questions for parametric families of linear recurrence sequences.  Canadian Journal of Mathematics, 1–24, 2020.
  * [28] J. Ouaknine and J. Worrell.  Decision problems for linear recurrence sequences.  In International Workshop on Reachability Problems, pages 21–28. Springer, 2012.
  * [29] B. Rosser.  The n-th prime is greater than nlogn.  Proceedings of the London Mathematical Society, 2(1):21–44, 1939\.
  * [30] B. Rosser.  Explicit bounds for some functions of prime numbers.  American Journal of Mathematics, 63(1):211–232, 1941.
  * [31] P. Scholze and J. Stix.  Why abc is still a conjecture, 2018.
  * [32] M. Sha.  Effective results on the Skolem problem for linear recurrence sequences.  Journal of Number Theory, 197:228–249, 2019.
  * [33] C. L. Stewart and K. Yu.  On the abc conjecture.  Mathematische Annalen, 291:225–230, 1991.
  * [34] C. L. Stewart and K. Yu.  On the abc conjecture, II.  Duke Mathematical Journal, 108(1):169–181, 2001.
  * [35] A. Surroca.  Sur l’effectivite du theoreme de Siegel et la conjecture abc.  Journal of Number Theory, 124(2):267–290, 2007.
  * [36] P. A. Vojta.  Diophantine approximations and value distribution theory, volume 1239.  Springer, 2006.
  * [37] M. Waldschmidt.  Diophantine Approximation on Linear Algebraic Groups: Transcendence Properties of the Exponential Function in Several Variables.  Springer-Verlag, 2000.
  * [38] K. Yu.  Linear forms in p-adic logarithms.  Acta Arithmetica, 53(2):107–186, 1989.
  * [39] K. Yu.  Linear forms in $p$-adic logarithms. II.  Compositio Mathematica, 74(1):15–113, 1990.
  * [40] K. Yu.  Linear forms in $p$-adic logarithms. III.  Compositio Mathematica, 91(3):241–276, 1994.
  * [41] K. Yu.  P-adic logarithmic forms and group varieties I.  1998\.
  * [42] K. Yu.  p-adic logarithmic forms and group varieties II.  Acta Arithmetica, 89(4):337–378, 1999.
  * [43] K. Yu.  P-adic logarithmic forms and group varieties III.  Forum Mathematicum, 19(2):187–280, 2007.


# CoBaIR: A Python Library for Context-Based Intention Recognition in Human-
Robot-Interaction

Adrian Lubitz1, Lisa Gutzeit1 and Frank Kirchner1,2 *This work was supported
through a grant of the German Federal Ministry of Economic Affairs and Climate
Action (BMWi, FKZ 50 RA 2022)1 Robotics Research Group, Department of Computer
Science, University of Bremen, 28359 Bremen, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>German Research
Center for Artificial Intelligence, 28359 Bremen, Germany
<EMAIL_ADDRESS>
###### Abstract

Human-Robot Interaction (HRI) becomes more and more important in a world where
robots integrate fast in all aspects of our lives but HRI applications depend
massively on the utilized robotic system as well as the deployment environment
and cultural differences. Because of these variable dependencies it is often
not feasible to use a data-driven approach to train a model for human intent
recognition. Expert systems have been proven to close this gap very
efficiently. Furthermore, it is important to support understandability in HRI
systems to establish trust in the system. To address the above-mentioned
challenges in HRI we present an adaptable python library in which current
state-of-the-art Models for context recognition can be integrated. For
Context-Based Intention Recognition a two-layer Bayesian Network (BN) is used.
The bayesian approach offers explainability and clarity in the creation of
scenarios and is easily extendable with more modalities. Additionally, it can
be used as an expert system if no data is available but can as well be fine-
tuned when data becomes available.

##  I Introduction

Our day-to-day lives are becoming increasingly involved with robotic devices.
The industry is currently changing from static robotic environments to dynamic
environments where humans collaborate with robots instead of operating robots.
In private homes, robotic applications like robot vacuums and digital
assistants for home automation are becoming regular household items [1].
Although this movement is drastically changing our society, interactions
between humans and robots are very command-driven and unnatural in contrast to
Human-Human Interaction (HHI) [2, 3]. While Large Language Models (LLM) like
_GPT-3_[4], _BERT_[5], _LLaMA_[6], and _LaMDA_[7] show very promising results
in general language understanding, they have problems with biases, alignment,
uncertainty estimation, and most importantly they lack multimodal
understanding of their surroundings which is a key feature for Intention
Recognition needed in modern HRI applications [8]. In general, data-driven
Intention Recognition systems [9, 10] have the drawback of relying on large
and complex (high dimensional) human-robot interaction data. Because in most
scenarios data is not available or impractical to record, data-driven
approaches are not suitable for this problem. Bayesian context-based Intention
Recognition is an approach to overcome those limitations and offer an expert
system that can be fine-tuned with data when it becomes available. Existing
research in this direction [11, 12, 13, 14, 15, 16] is promising but
introduces very complex network structures which induce the need for the
designer of an HRI scenario to have a profound knowledge of Bayesian
probability theory. Furthermore, it makes adaptations and re-configuration of
the network cumbersome. To overcome the aforementioned limitations of current
systems we propose a two-layer BN for context-based Intention Recognition. The
simple structure enables us to make several optimizations that allow the
designer to concentrate on the HRI scenario instead of Bayesian probability
theory.

In [17] we proposed a first concept on how the design for context-based
Bayesian Intention Recognition in HRI scenarios can be described in a more
compact and intuitive way. In this paper, we introduce Context-Based Intention
Recognition (_CoBaIR_), a python software library that comes with the power to
infer intentions from the current context — context describes every observable
aspect in an HRI scenario. Furthermore, _CoBaIR_ pays great attention to the
design process of HRI scenarios. It provides a configuration format that
decreases the number of values that need to be set during the design process
of the Bayesian Network from an exponential to a linear scale. Additionally,
it provides a Graphical User Interface (GUI) which visualizes the two-layer BN
with its weights and offers an intuitive way of configuring it.

This paper starts by pointing out the key challenges in Intention Recognition
for HRI in Section II-A. In Section III we provide a detailed view of the
proposed two-layer BN structure and highlight how this structure can make the
process of designing HRI scenarios easier and faster. In Section IV we
highlight the advantages of the proposed approach. We point out the most
important features of the python implementation in Section V. In Section VI we
show an example of how the described python library was used in a research
project. Finally, in Section VII we conclude with the interpretation of the
results and an outlook on future work.

##  II Related Work

During the extensive literature review conducted for this paper, it was
observed that existing intention recognition implementations were often
tightly coupled with specific modalities [18, 19] or designed exclusively for
particular scenarios [20, 21]. In some cases, these limitations were found to
coexist [22], further hindering the applicability of such implementations.
However, as part of our objective to provide a generic framework for HRI
within the KiMMI Project, we recognized the need for a solution that could
exhibit high flexibility and adaptability across various modalities and
scenarios.

###  II-A Challenges in Intention Recognition

HRI depends in many ways on the scenario at hand. For Intention Recognition we
define the following challenges that need to be addressed in order to
implement natural and meaningful HRI systems:

#### Hardware constraints

Just like humans, robots come in all shapes and colors. More precisely, social
robots for HRI have different sensing modalities to perceive the human and its
surrounding. While some robots are equipped with stereo-vision or RGBD-camera
systems and multiple microphones for echolocation, as well as LiDAR for
navigation, a simple digital assistant may only be equipped with one
microphone. The designer of the HRI scenario needs to know which sensory
modalities are available for the system in question. [23, 24, 25]

#### Application specifics

Furthermore, the designer needs to know about the application scenario which
can vary from space applications over industrial to domestic applications. In
all of these different application scenarios gestures, voice commands, etc.
can mean different things. [23, 25]

#### Cultural differences

Research on cultural differences in social robotics is an often neglected
topic although social robots will be deployed in multi-cultural place e.g.
airports in the future in a more and more globalized world. One behavior may
have different meanings in different cultures. Therefore the designer must be
aware of the cultural differences, and for different cultures, different HRI
scenarios must be designed. [25, 26]

#### Individual differences

When we think about developing robots that interact with humans in a very
intuitive way we need to ask ourselves what is intuitive for us. Intuitive may
be slightly different from person to person even within one cultural group.
Humans interact slightly different on an individual level based on their
knowledge about the interaction partner. If the interaction partner is not
known a default is chosen which allows for adaptation in the future. While
this behavior is very subtle and unconscious in HHI it is an important factor
while designing HRI scenarios with the possibility for inter-personal
adaptation. [24, 25, 26]

#### Trust & Acceptance

Trust in a robotic system is less scenario specific than the aforementioned
challenges and can therefore not as obviously be integrated into the design
process of an HRI scenario. Trust is primarily connected with the human’s
expectation of the behavior of a robot. If the robot behaves accordingly to
the human’s expectations the human can foresee the behavior and build a model
of trust for the robot’s abilities. Secondarily, it is connected with the
explainability of a behavior. If the human is not able to foresee the robot’s
behavior because it is, e.g. not completely deterministic or too complex to
foresee, the human will seek reasons for the observed behavior. If reasons can
be found trust can still be maintained, while if no reason can be found the
trust will decrease immensely. This model of trust can be very complex in a
way that trust exists for specific abilities but not for others. The trust for
specific abilities weighted with their specific importance for the human
determines the acceptance in the system. [24, 23, 25]
In Section III \- VI we illustrate how we address these challenges, which
advantages arise from the proposed approach, how it is implemented, and how it
can be used to model an HRI scenario.

##  III Architecture of the two-layer Bayesian Network for Context-based
Intention Recognition

In Section II-A we highlighted some of the key challenges in Intention
Recognition in HRI. We propose a two-layer BN to address these challenges in a
computational and data efficient as well as an intuitive way. The general
structure of the BN for context-based Intention Recognition is depicted in
Figure 1. In [12] and [15] Bayesian Networks (BNs) with 3 or more layers are
proposed to additionally model actions. We believe the two-layered structure
comes with several advantages over using three or more layers. From a
usability point of view, it allows us to assume that the designer of the HRI
scenario has no prior knowledge about Bayesian probability in general and BNs
in specific.

In this way every contexts can be treated in the same way, as an observable
phenomenon. Using further layers would give actions, as suggested by [12], a
special meaning.

Figure 1: The architecture of the two-layer Bayesian network allows for a high
degree of flexibility

This special meaning however is not always valid and furthermore, it
introduces a bias towards actions in the design of HRI scenarios. In cases
where external context is of more importance than the performed action, the
introduced bias towards actions can distort the design of the scenario.

In cases where external factors like the time have an influence on the
intention the action itself should not play a predominant role. The action of
grasping a mug could lead to the intention make coffee but at night the
intention becomes more unrealistic and a higher probability should be given to
the intention store mug. This is a very simple example of how the action bias
could lead to the intention make coffee at night where the correct intention
should be store mug.

Modeling relevant observable phenomena as context reduces the action bias and
the complexity of the BN and allows the designer of the scenario to
concentrate on the specifics of the scenario without the need for in-depth
knowledge about the underlying probabilistic modeling.

The final step we took to uncouple probabilistic modeling from the intuitive
design of HRI scenarios is to make some assumptions over the given two-layer
structure. These assumptions reduce the probabilistic notion and drastically
reduce the number of values that need to be set by a human expert to describe
a scenario.

### Independent context and intentions

The basic assumption is that all intentions are independent of each other as
well as all contexts are independent of each other. This allows for the strict
two-layer structure which has only connections between context and intentions
as depicted in Figure 1.

### Binary intentions

We consider all intentions as binary — an intention is either present or not.
This allows us to concentrate on the positive (an intention is present) case
while designing the scenario and consider the negative (an intention is not
present) case as its complement while calculating the probabilistic model.
This small constraint already cuts the number of values that need to be set by
the human expert in half. Contexts on the other hand can have as many discrete
instantiations as needed to create a meaningful HRI scenario.

### Single Condition Assumption

We make the Single Condition Assumption (SCA) which implies that every context
has an individual and independent influence on a specific intention. Using
this assumption, we can approximate the conditional probability
$P(I_{m}|C_{1,l},C_{2,l},...,C_{k,l})$ as the average over all single
condition probabilities $P(I_{m}|C_{k,l})$, where $I_{m}$ is the _m_ -th
intention and $C_{k,l}$ is the _l_ -th instantiation of the _k_ -th context.
This allows us to calculate the conditional probability of the _m_ -th
intention given the first instantiation for all contexts in the following way:

$P(I_{m}|C_{1,l},C_{2,l},...,C_{k,l})=\frac{\sum^{k}_{\hat{k}=1}P(I_{m}|C_{\hat{k},l})}{k}$
(1)

### Influence values on a Likert scale

We generalize the single condition probability $P(I_{m}|C_{k,l})$ as an
influence value $v_{k,l,m}$ on a six-point Likert scale [27] for every
context-intention-tuple $(C_{k,l},I_{m})$ to make them more manageable. The
scale is mapped in the following fashion: $0\mapsto 0\%;1\mapsto 5\%;2\mapsto
25\%;3\mapsto 50\%;4\mapsto 75\%;5\mapsto 95\%$

The above-mentioned assumptions reduce the amount of values that need to be
set by a human expert from an exponential growth given through

$V(i,j,c,n)=\sum^{j}_{\hat{j}=1}c_{\hat{j}}+i\times\prod^{i}_{\hat{i}=1}n_{\hat{i}}\times\prod^{j}_{\hat{j}=1}c_{\hat{j}}$
(2)

to a linear growth given through

$V(i,j,c)=(i+1)\times\sum^{j}_{\hat{j}=1}c_{\hat{j}}$ (3)

where $V$ is the number of values to be set, $i$ is the number of intentions,
$j$ is the number of contexts, $c_{j}$ is the number of context instantiations
for the _j_ -th context, $n_{i}$ is the number of intention instantiations for
the _i_ -th intention. The $\sum^{j}c_{j}$ in Equations 2 and 3 describes the
amount of _a priori_ probabilities for all context instantiations. The
remaining term describes the values needed to fill the Conditional Probability
Tables (CPTs) manually for Equation 2 and automatically for Equation 3. The
SCA contributes in a huge way to the ease of designing HRI scenarios but
ignores cases in which joint probabilities are necessary. An example could be
the handling of voice commands through a robot that is able to estimate
directed speech over Visual Voice Activity Detection (VVAD) as shown in [28].
The robot should infer the intention pick up tool with a higher probability if
the speech command for picking up a tool was emitted AND the speech was
directed towards the robot than one of them individually. Only the speech
command should have a high probability to infer pick up tool but there is
still a chance that the robot was picking up noise or the speech was not
directed towards the robot. The context of directed speech on the other hand
does not have a high probability for any intention individually. For those
special cases, we provide the possibility to set (partially) conditioned
influence values containing multiple contexts that provide more information
when combined. [17]

While the two-layer BN was originally not designed to handle temporal
dependencies, it is possible to model the previously inferred intention as
context. Using a recursive pattern like this, it is possible to model a
temporal dependency under the Markov assumption [29].

##  IV Advantages of a two-layer Bayesian Network for Context-Based Intention
Recognition

_CoBaIR_ uses a two-layer BN to represent the dependencies between contexts
and intentions. In this section we want to highlight the key advantages of
this structure:

### Flexibility

The biggest advantage of the architecture is its flexibility. It allows for
the usage of any algorithm for context creation, whether it be probabilistic,
heuristic, data-driven or any other approach. The generated contexts will be
used as the input to the two-layer BN which fuses the context information to
jointly infer a probability distribution over all possible intentions. On the
one hand, using _CoBaIR_ as an expert system the simplifications explained in
Section III allow the human designer to create a scenario in a fast and
intuitive manner. Fast design and adaptation of HRI scenarios helps
researchers to concentrate on the specifics of an experiment and therefore
reach results faster and more reliably without any deeper knowledge about
Bayesian probability and how to configure BNs. On the other hand the two-layer
BN can be trained or fine-tuned with data, which allows to gradually shift
from a system trained by an expert to a data-driven approach. Furthermore, the
simple structure of the BN allows for the easy removal and addition of
contexts. This makes an iterative prototyping approach, where the HRI scenario
is build up over time, possible.

### Uncertainty Quantification and Explainability

Uncertainty Quantification (UQ) is an often neglected topic, especially in the
field vision based tasks [30]. The Bayesian approach for Intention Recognition
offers the implicit advantage that it comes with a inherently good UQ due to
the probabilistic nature of the model. Additionally, the compact structure of
the two-layer BN allows users to easily identify the contexts that played the
predominant role in the decision making. Using this interpretable compact
structure we are able to generate explanations to understand the decisions
made by a robot using _CoBaIR_ as its Intention Recognition system.
Explainability and UQ in a system strongly increases the trust in the system
and therefore the acceptance to use the system in general [31].

### Modularity

Another advantage of the described architecture is its modularity. Using a
two-layer BN to fuse the output of different modules that provide contexts
makes it possible to use existing solutions for context creation, like PAZ
[32] which provides a large variety of models for visual perception, as well
as use case specific models that need to be trained from scratch. Furthermore,
it is possible to switch seamlessly between different models on the fly for
evaluation and optimization of HRI scenarios.

### Handling missing input

While most data-driven approaches have problems handling missing input [33]
and additional data needs to be recorded or artificially generated, BNs
provide the advantage of defining _a priori_ probabilities for the input. The
_a priori_ probabilities for the context can be estimated by an expert or
calculated from a few observations. In this way knowledge about missing input
can be incorporated and during inference time the missing inputs will handled
accordingly.

The above-mentioned advantages of the two-layer BN highlight why we think a
two-layer BN is suitable to enhance the design, inference and explainability
of Intention Recognition in HRI scenarios.

##  V CoBaIR: a Python Library

_CoBaIR_ is a python library for Context-Based Intention Recognition. The
library allows to create complex HRI scenarios in a fast and intuitive manner.
Furthermore, it provides a GUI which visualizes the underlying two-layer BN
and guides through the configuration procedure.

_CoBaIR_ is divided into two parts:

###  V-A Core library

The core library provides all the key features described in Section III to
make the design of HRI scenarios intuitive and fast. It mainly provides the
class BayesNet which handles the creation of the two-layer BN from a given
configuration. The configuration is provided in YAML and its fields and format
is depicted in Listing 1. The format contains the fields _contexts_ ,
_instantiations_ , _intentions_ and _decision_threshold_. _contexts_ gives
names to the observable phenomena, like _weather_. _instantiations_ are the
discrete instantiations of that phenomenon, like _cloudy, rainy, sunny_. The
binary _intentions_ are the inferable intentions in the scenario, like _turn
on sprinkler_. All _instantiations_ have an a priori probability which needs
to be set and furthermore the _instantiations_ have an influence value for
each _intention_. Additionally, there is a field _decision_threshold_ which
can be a value between 0 and 1 and denotes the threshold an intention’s
likelihood needs to surpass during inference to be considered as the inferred
intention. If the likelihood of the most likely intention is below the
threshold, None is returned as the inferred intention.

Listing 1: Configuration format for Intention Recognition with _CoBaIR_

⬇

contexts:

context 1:

instantiation 1 : float

.

instantiation m_1 : float

context n:

instantiation 1 : float

.

instantiation m_n : float

intentions:

intention 1:

context 1:

instantiation 1: int # one out of [5, 4, 3, 2, 1, 0]

.

instantiation m_1: int # one out of [5, 4, 3, 2, 1, 0]

context n:

instantiation 1: int # one out of [5, 4, 3, 2, 1, 0]

.

instantiation m_n: int # one out of [5, 4, 3, 2, 1, 0]

intention p:

context 1:

instantiation 1: int # one out of [5, 4, 3, 2, 1, 0]

.

instantiation m_1: int # one out of [5, 4, 3, 2, 1, 0]

context n:

instantiation 1: int # one out of [5, 4, 3, 2, 1, 0]

.

instantiation m_n: int # one out of [5, 4, 3, 2, 1, 0]

decision_threshold: float

The core library provides the means to validate, load and save the
configuration. With the information from the configuration a fully defined
two-layer BN will be created. bnlearn [34] is utilized as a backend to handle
the two-layer BNs and do inference on them. The Application Programming
Interface (API) is fully documented an publicly available under https://dfki-
ric.github.io/CoBaIR/API.

### Graphical User Interface

On top of the optimizations we described in Section III, we provide
researchers and practitioners with a helpful GUI depicted in Figure 2.

Figure 2: The GUI of _CoBaIR_ supports the designer of an HRI scenario.

The GUI supports the designer to create a configuration and makes sure it is
always valid and provides helpful insights in the case of an invalid
configuration. Furthermore, it visualizes the two-layer BN in a live view
which helps to keep a good overview of the configuration.
The complete Python software package is available on PyPI. Additionally, the
code is open sourced on GitHub and open for contributions.

##  VI Application of CoBaIR

We applied _CoBaIR_ in an interaction scenario in the project KiMMI
SF111https://robotik.dfki-bremen.de/en/research/projects/kimmi-sf/ to infer
the intentions of a human operator. This user controlled a simulated human
during the interaction with the simulated robotic arm from Universal Robot UR5
mounted on a movable base. The simulated environment is depicted in Figure 3.
The human can be navigated through the environment, which contains shelves in
which different tools are stored, a working station with an inactive robotic
system MANTIS222https://robotik.dfki-bremen.de/en/research/robot-
systems/mantis/ that should be repaired, and on the left back side a dark
area. Based on the human behavior the robotic system should infer the human’s
intention in order to react accordingly. For example, if the human wants to
repair the inactive MANTIS robot, the supporting robot should bring the tool
which is stored in the shelf, which is shown on the right side in Figure 3.

Figure 3: Screenshots of the HRI application scenario, in which a simulated
human interacts with the UR5 robot mounted on a movable base. The scenario
includes a shelf with tools stored in it, a dark area which should be lighted
by the UR5 if the human needs support in the dark area, and a work station
highlighted in green area where the human can work on the robotic system
MANTIS and the UR5 should assist by bringing the desired tool.

The BN used to infer the human intentions is shown in Figure 4. The human
behavior is measured based on four contexts: hand opening, human pose,
location of interest, and speech commands, where the first three are
determined in the simulation and the speech commands are captured from the
operator of the simulation with a microphone. Each context is discretized to
different values, e.g., hand opening can be either open or closed. This is
shown on the left side of Figure 4 for each context. In the presented
scenario, based on the given contexts, the following five intentions should be
inferred: 1\. go work station, i.e., the human wants to go to the work
station; 2\. go dark space, i.e., the human want to go to the dark area; 3\.
robot bring tool, i.e., the human wants the robot to bring the tool stored in
the shelf; 4\. robot stop, i.e., the human wants the robot to stop its current
action; 5\. robot store tool, the human wants the robot to store the tool back
to the shelf.

Figure 4: BN of the application scenario. Based on four different contexts
(left side), five possible human intentions (right side) should be inferred.
Arrows indicate the influence of the possible values of the contexts with
small numbers representing the assigned influence value.

Using _CoBaIR_ and the included GUI shown in Figure 2, the BN shown in Figure
4 could easily be designed and optimized for the described scenario. With the
resulting BN all human intentions could be reliably inferred and the correct
reactions of the UR5 could be triggered to realize a successful interaction
between the human and the robotic system.

During the KiMMI SF project we followed an agile development process which led
to multiple incremental as well as complete changes in the design of the HRI
scenario. _CoBaIR_ enabled us to incorporate these changes fast and
effectively in the development process.

##  VII Conclusion and future work

We presented the python library _CoBaIR_ in this paper. We demonstrated that
the concept from [17] using a two-layer BN with the assumptions highlighted in
III can be effectively implemented. Additionally, we provided a GUI to
visualize and guide the design process for HRI scenarios. In Section VI we
showed that _CoBaIR_ was successfully used in the KiMMI SF project. We advise
practitioners to use the state-of-the-art models for context recognition from
the open source library PAZ [32] which is constantly updated with perception
models for autonomous systems. In the future we plan on providing tutorials
and examples on how to use _CoBaIR_ with PAZ. While it is theoretically
possible to train or fine-tune the two-layer BN with data from the scenario,so
far the fine-tuning and data-driven training from scratch was not tested. In
future works we want to investigate the capability of data-driven training and
fine-tuning within _CoBaIR_. Additionally, we want to quantify the effect in
terms of quality and speed which _CoBaIR_ has on the design of HRI scenarios,
conducting user studies comparing the time and cognitive load to create a HRI
scenario using our solution in contrast to creating a HRI scenario using CPTs
for the creation of the BN.. We will be incorporating _CoBaIR_ in future
projects and by making _CoBaIR_ open source we hope to provide a helpful tool
for researchers and practitioners in HRI all over the world. Everyone is
welcome to use, give feedback and contribute to _CoBaIR_ through GitHub.

## ACKNOWLEDGMENT

This work was supported through a grant of the German Federal Ministry of
Economic Affairs and Climate Action (BMWi, FKZ 50 RA 2022).

## References

  * [1] L. Fortunati, A. Esposito, and G. Lugano, “Introduction to the Special Issue “Beyond Industrial Robotics: Social Robots Entering Public and Domestic Spheres”,” _The Information Society_ , vol. 31, no. 3, pp. 229–236, May 2015.
  * [2] H. Yan, M. H. Ang, and A. N. Poo, “A Survey on Perception Methods for Human–Robot Interaction in Social Robots,” _International Journal of Social Robotics_ , vol. 6, no. 1, pp. 85–119, Jan. 2014.
  * [3] D. Y. Y. Sim and C. K. Loo, “Extensive assessment and evaluation methodologies on assistive social robots for modelling human–robot interaction – A review,” _Information Sciences_ , vol. 301, pp. 305–344, Apr. 2015.
  * [4] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, A. Herbert-Voss, G. Krueger, T. Henighan, R. Child, A. Ramesh, D. M. Ziegler, J. Wu, C. Winter, C. Hesse, M. Chen, E. Sigler, M. Litwin, S. Gray, B. Chess, J. Clark, C. Berner, S. McCandlish, A. Radford, I. Sutskever, and D. Amodei, “Language Models are Few-Shot Learners,” July 2020, arXiv:2005.14165 [cs]. [Online]. Available: http://arxiv.org/abs/2005.14165
  * [5] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding,” May 2019, arXiv:1810.04805 [cs]. [Online]. Available: http://arxiv.org/abs/1810.04805
  * [6] H. Touvron, T. Lavril, G. Izacard, X. Martinet, M.-A. Lachaux, T. Lacroix, B. Rozière, N. Goyal, E. Hambro, F. Azhar, A. Rodriguez, A. Joulin, E. Grave, and G. Lample, “LLaMA: Open and Efficient Foundation Language Models,” Feb. 2023, arXiv:2302.13971 [cs]. [Online]. Available: http://arxiv.org/abs/2302.13971
  * [7] R. Thoppilan, D. De Freitas, J. Hall, N. Shazeer, A. Kulshreshtha, H.-T. Cheng, A. Jin, T. Bos, L. Baker, Y. Du, Y. Li, H. Lee, H. S. Zheng, A. Ghafouri, M. Menegali, Y. Huang, M. Krikun, D. Lepikhin, J. Qin, D. Chen, Y. Xu, Z. Chen, A. Roberts, M. Bosma, V. Zhao, Y. Zhou, C.-C. Chang, I. Krivokon, W. Rusch, M. Pickett, P. Srinivasan, L. Man, K. Meier-Hellstern, M. R. Morris, T. Doshi, R. D. Santos, T. Duke, J. Soraker, B. Zevenbergen, V. Prabhakaran, M. Diaz, B. Hutchinson, K. Olson, A. Molina, E. Hoffman-John, J. Lee, L. Aroyo, R. Rajakumar, A. Butryna, M. Lamm, V. Kuzmina, J. Fenton, A. Cohen, R. Bernstein, R. Kurzweil, B. Aguera-Arcas, C. Cui, M. Croak, E. Chi, and Q. Le, “LaMDA: Language Models for Dialog Applications,” Feb. 2022, arXiv:2201.08239 [cs]. [Online]. Available: http://arxiv.org/abs/2201.08239
  * [8] A. Tamkin, M. Brundage, J. Clark, and D. Ganguli, “Understanding the Capabilities, Limitations, and Societal Impact of Large Language Models,” Feb. 2021.
  * [9] Z. Wang, B. Wang, H. Liu, and Z. Kong, “Recurrent convolutional networks based intention recognition for human-robot collaboration tasks,” in _2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , Oct. 2017, pp. 1675–1680.
  * [10] Z. Liu and J. Hao, “Intention Recognition in Physical Human-Robot Interaction Based on Radial Basis Function Neural Network,” _Journal of Robotics_ , vol. 2019, p. e4141269, Apr. 2019.
  * [11] L. Pereira and T. A. Han, _Intention Recognition via Causal Bayes Networks Plus Plan Generation_ , Oct. 2009, vol. 5816.
  * [12] L. M. Pereira and H. T. Anh, “Intention Recognition with Evolution Prospection and Causal Bayes Networks,” in _Computational Intelligence for Engineering Systems: Emergent Applications_ , ser. Intelligent Systems, Control and Automation: Science and Engineering, A. Madureira, J. Ferreira, and Z. Vale, Eds. Dordrecht: Springer Netherlands, 2011, pp. 1–33.
  * [13] D. Panagopoulos, G. Petousakis, R. Stolkin, G. Nikolaou, and M. Chiou, “A Bayesian-Based Approach to Human Operator Intent Recognition in Remote Mobile Robot Navigation,” Sept. 2021.
  * [14] R. Kelley, A. Tavakkoli, C. King, A. Ambardekar, M. Nicolescu, and M. Nicolescu, “Context-Based Bayesian Intent Recognition,” _IEEE Transactions on Autonomous Mental Development_ , vol. 4, no. 3, pp. 215–225, Sept. 2012.
  * [15] K. A. Tahboub, “Intelligent Human-Machine Interaction Based on Dynamic Bayesian Networks Probabilistic Intention Recognition,” _Journal of Intelligent and Robotic Systems_ , vol. 45, no. 1, pp. 31–52, Jan. 2006.
  * [16] S. Jain and B. Argall, “Recursive Bayesian Human Intent Recognition in Shared-Control Robotics,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Oct. 2018, pp. 3905–3912.
  * [17] A. Lubitz, O. Arriaga, T. Hassan, N. Hoyer, and E. A. Kirchner, “A Bayesian Approach to Context-based Recognition of Human Intention for Context-Adaptive Robot Assistance in Space Missions,” 2022.
  * [18] Y.-M. Jang, R. Mallipeddi, S. Lee, H.-W. Kwak, and M. Lee, “Human intention recognition based on eyeball movement pattern and pupil size variation,” _Neurocomputing_ , vol. 128, pp. 421–432, Mar. 2014.
  * [19] D. Zhang, L. Yao, K. Chen, S. Wang, X. Chang, and Y. Liu, “Making Sense of Spatio-Temporal Preserving Representations for EEG-Based Human Intention Recognition,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 7, pp. 3033–3044, July 2020.
  * [20] Q. Wang, W. Jiao, R. Yu, M. T. Johnson, and Y. Zhang, “Virtual Reality Robot-Assisted Welding Based on Human Intention Recognition,” _IEEE Transactions on Automation Science and Engineering_ , vol. 17, no. 2, pp. 799–808, Apr. 2020.
  * [21] Y. Xing, C. Lv, H. Wang, H. Wang, Y. Ai, D. Cao, E. Velenis, and F.-Y. Wang, “Driver Lane Change Intention Inference for Intelligent Vehicles: Framework, Survey, and Challenges,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 5, pp. 4377–4390, May 2019.
  * [22] C. Zhu, W. Sun, and W. Sheng, “Wearable sensors based human intention recognition in smart assisted living systems,” in _2008 International Conference on Information and Automation_ , June 2008, pp. 954–959.
  * [23] T. B. Sheridan, “Human–Robot Interaction: Status and Challenges,” _Human Factors_ , vol. 58, no. 4, pp. 525–532, June 2016\.
  * [24] B. Hayes and B. Scassellati, “Challenges in Shared-Environment Human-Robot Collaboration,” 2007.
  * [25] C. Jost, B. L. Pévédic, T. Belpaeme, C. Bethel, D. Chrysostomou, N. Crook, M. Grandgeorge, and N. Mirnig, _Human-Robot Interaction: Evaluation Methods and Their Standardization_. Springer Nature, May 2020.
  * [26] C. Nehaniv, K. Dautenhahn, J. Kubacki, M. Haegele, C. Parlitz, and R. Alami, “A methodological approach relating the classification of gesture to identification of human intent in the context of human-robot interaction,” in _ROMAN 2005. IEEE International Workshop on Robot and Human Interactive Communication, 2005._ , Aug. 2005, pp. 371–377.
  * [27] R. Likert, “A technique for the measurement of attitudes,” _Archives of Psychology_ , vol. 22 140, pp. 55–55, 1932.
  * [28] A. Lubitz, M. Valdenegro-Toro, and F. Kirchner, “The VVAD-LRS3 Dataset for Visual Voice Activity Detection,” in _7th International Conference on Human Computer Interaction Theory and Applications_ , Mar. 2023, pp. 39–46.
  * [29] A. A. Markov, “The Theory of Algorithms,” _Journal of Symbolic Logic_ , vol. 18, no. 4, pp. 340–341, 1953, publisher: Association for Symbolic Logic.
  * [30] M. Valdenegro-Toro, “I Find Your Lack of Uncertainty in Computer Vision Disturbing,” Apr. 2021.
  * [31] L. Sanneman and J. A. Shah, “Trust Considerations for Explainable Robots: A Human Factors Perspective,” May 2020.
  * [32] O. Arriaga, M. Valdenegro-Toro, M. Muthuraja, S. Devaramani, and F. Kirchner, “Perception for Autonomous Systems (PAZ),” Oct. 2020.
  * [33] N. Ipsen, P.-A. Mattei, and J. Frellsen, “How to deal with missing data in supervised deep learning?” p. 6, 2020.
  * [34] E. Taskesen, “Learning Bayesian Networks with the bnlearn Python Package.” 1 2020. [Online]. Available: https://erdogant.github.io/bnlearn


# The induced permittivity increment of electrorheological fluids in an
applied electric field in association with chain formation: A Brownian
Dynamics simulation study

Dávid Fertig  Dezső Boda<EMAIL_ADDRESS>Center for Natural Sciences,
University of Pannonia, Egyetem u. 10, Veszprém, 8200, Hungary  István Szalai
Research Centre for Engineering Sciences, University of Pannonia, Egyetem u.
10, Veszprém, 8200, Hungary  Institute of Mechatronics Engineering and
Research, University of Pannonia, Gasparich Márk u. 18/A, Zalaegerszeg, 8900,
Hungary

###### Abstract

We report Brownian Dynamics simulation results for the relative permittivity
of electrorheological (ER) fluids in an applied electric field. The relative
permittivity of an ER fluid can be calculated from the Clausius-Mosotti (CM)
equation in the small applied field limit. When a strong field is applied,
however, the ER spheres are organized into chains and assemblies of chains in
which case the ER spheres are polarized not only by the external field but by
each other. This manifests itself in an enhanced dielectric response, e.g., in
an increase in the relative permittivity. The correction to the relative
permittivity and the time dependence of this correction is simulated on the
basis of a model in which the ER particles are represented as polarizable
spheres. In this model, the spheres are also polarized by each other in
addition to the applied field. Our results are qualitatively similar to those
obtained by Horváth and Szalai experimentally (Phys. Rev. E, 86, 061403,
2012). We report characteristic time constants obtained from bi-exponential
fits that can be associated with formation of pairs and short chains as well
as with aggregation of chains. The electric field dependence of the induced
dielectric increment reveals the same qualitative behavior that experiments
did: three regions with different slopes corresponding to different
aggregation processes are identified.

##  I Introduction

In electrorheological (ER) fluids Winslow (1949) fine non-conducting solid
particles are suspended in an electrically insulating liquid with the
particles having larger relative permittivity than the solvent. Then, an
applied electric field induces polarization charges at the arising dielectric
boundaries that can be expressed as a multipole expansion with dipoles being
the dominant terms.

The interactions of these dipoles lead to a structural change in the ER fluid
known as the ER response. This structural change is the aggregation of ER
particles first into shorter, then into longer chains due to the fact that the
head-to-tail position of two dipoles along the direction of the applied field
is a minimum-energy configuration. In the case of strong applied fields, the
chains form larger clusters, for example, columnar structures.

This structural change results in changes in major physical properties of the
ER fluid. The externally controllable, fast and reversible change in
viscosity, for example, makes ER fluids a central component of various
devices, such as brakes, clutches, dampers, and valves Duclos  _et al._
(1992); Havelka and Filisko (1995).

Another physical quantity whose change can be relatively easily tracked by
measuring the change in the capacitance of a measuring cell when the electric
field is switched on is the relative permittivity, $\epsilon$. Several
experimental studies have been reported for the nonlinear dielectric
properties of ER fluids. Adolf and Garino (1995); Tao and Roy (1994); Wen  _et
al._ (1997, 1998); Rzoska and Zhelezny (2006) Horváth and Szalai Horváth and
Szalai (2012) have proposed a new method that made the measurement of
continuous changes in the increments in permittivity possible. They determined
the field dependence of the change in dielectric permittivity

$\Delta\epsilon(E)=\epsilon(E)-\epsilon(E=0),$ (1)

with $t\rightarrow\infty$ and also the time dependence

$\Delta\epsilon(t)=\epsilon(t)-\epsilon(t=0),$ (2)

where the electric field is switched on from 0 to $E$ at $t=0$. The time
dependence and electric field dependence are shown in Fig. 3 and Fig. 4 of
Ref. Horváth and Szalai, 2012, respectively.

They found that the time dependence can be described with a bi-exponential fit

$\Delta\epsilon(t)=A\left(1-e^{-t/\tau_{1}}\right)+B\left(1-e^{-t/\tau_{2}}\right).$
(3)

They hypothesized that the time constant $\tau_{1}$ can be associated with the
process of formation of pairs (and, perhaps, short chains), because $\tau_{1}$
is very similar to the characteristic time of pair formation derived from a
kinetic rate theory. Baxter-Drayton and Brady (1996) The time constant
$\tau_{2}$, on the other hand, was heuristically corresponded to formation of
long chains and their aggregation and proved to be an order of magnitude
larger than $\tau_{1}$, $\tau_{2}\approx 10\tau_{1}$.

In this paper, we investigate how change in the dielectric permittivity is
associated with structural changes (formation of chains of various lengths,
and the aggregation of chains into columnar structures) by computer
simulations. Although several simulation studies have been reported for
cluster formation Klingenberg  _et al._ (1989); See and Doi (1991); Toor
(1993); Hass (1993); Climent  _et al._ (2004); Domínguez-García  _et al._
(2007), order parameters Hass (1993); Tao and Jiang (1994a, b); Baxter-Drayton
and Brady (1996); Enomoto and Oba (2002), diffusion constant Klingenberg  _et
al._ (1989); Whittle (1990); Hass (1993), pair distribution functions Whittle
(1990); Hass (1993), relaxation times Heyes and Melrose (1990); Toor (1993);
Hass (1993); Cao  _et al._ (2006), aggregation kinetics See and Doi (1991),
and stress under shear Heyes and Melrose (1990); Whittle (1990); Bonnecaze and
Brady (1992); Baxter-Drayton and Brady (1996); Cao  _et al._ (2006), we are
not aware of any paper addressing the dielectric properties of ER fluids.

The field dependence revealed three regimes with different slopes of the
$\Delta\epsilon$ vs. $E$ function (Fig. 4 of Ref. Horváth and Szalai, 2012).
Horváth and Szalai hypothesized that at low electric fields (below a threshold
value $E_{1}$) the induced dipoles are not strong enough to generate chain
formation and $\Delta\epsilon$ increases with $E$, because the probability of
the ER spheres to approach each other becomes larger. Above $E_{1}$ chain
formation begins with a large slope, because the parts of a chain can find one
another easier at large electric fields. Above the threshold value
$E_{\mathrm{h}}$ the chains start to form columnar structures. This process is
less accelerated by the large electric field, because chains can also repulse
each other when they are not in the appropriate mutual configuration with
respect to each other. Here, we support this hypothesis with computer
simulation results.

##  II Models and methods

The ER fluid is modeled as monodisperse dielectric spheres of relative
permittivity $\epsilon_{\mathrm{in}}$ inside the sphere immersed in a carrier
liquid of relative permittivity $\epsilon_{\mathrm{out}}$. The radius of the
spheres is $R$, while their diameter is $d=2R$. If an electric field,
$\mathbf{E}$ is applied on the sphere, a polarization charge density is
induced on the surface of the sphere that can be approximated with an ideal
point dipole placed in the center of the sphere computed as Jackson (1999)

$\bm{\mu}=4\pi\epsilon_{0}\left(\frac{\epsilon_{\mathrm{in}}-\epsilon_{\mathrm{out}}}{\epsilon_{\mathrm{in}}+2\epsilon_{\mathrm{out}}}\right)R^{3}\mathbf{E}=\alpha\mathbf{E},$
(4)

where

$\alpha=4\pi\epsilon_{0}\left(\frac{\epsilon_{\mathrm{in}}-\epsilon_{\mathrm{out}}}{\epsilon_{\mathrm{in}}+2\epsilon_{\mathrm{out}}}\right)R^{3}$
(5)

is the particle polarizability, $E=|\mathbf{E}|$, and $\epsilon_{0}$ is the
permittivity of vacuum. If it is further assumed that the characteristic time
of the rearrangement of the surface charge during the movement of the
particles is much smaller than the characteristic time of the rotation of the
ER particles, the $\bm{\mu}$ dipole always points into the direction of
$\mathbf{E}$ even if the sphere rotates.

If we take a system of $N$ particles at positions $\\{\mathbf{r}_{j}\\}$, the
electric field exerted on dipole $i$ by dipole $j$ is

$\mathbf{E}_{j}(\mathbf{r}_{i})=\frac{1}{4\pi\epsilon_{0}}\frac{3\mathbf{n}_{ij}(\mathbf{n}_{ij}\cdot\bm{\mu}_{j})-\bm{\mu}_{j}}{r_{ij}^{3}},$
(6)

where $r_{ij}=|\mathbf{r}_{ij}|$ and $\mathbf{n}_{ij}=\mathbf{r}_{ij}/r_{ij}$
with $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$.

In Eq. 4, the electric field at $\mathbf{r}_{i}$ is a sum of the applied
field, $\mathbf{E}^{\mathrm{appl}}$ (points to the $z$ direction), and the
electric field produced by all the other dipoles,
$\mathbf{E}(\mathbf{r}_{i})=\sum_{j\neq j}\mathbf{E}_{j}(\mathbf{r}_{i})$. The
total dipole moment

$\bm{\mu}^{\mathrm{tot}}_{i}=\alpha\mathbf{E}^{\mathrm{appl}}+\alpha\mathbf{E}(\mathbf{r}_{i})=\bm{\mu}^{\mathrm{appl}}_{i}+\bm{\mu}^{\mathrm{part}}_{i}$
(7)

then is induced by these two components and is split into the terms
$\bm{\mu}^{\mathrm{appl}}_{i}$ and $\bm{\mu}^{\mathrm{part}}_{i}$ accordingly.
The dipole moment $\bm{\mu}_{i}^{\mathrm{appl}}$ induced by the applied field
is constant, while the dipole moment $\bm{\mu}_{i}^{\mathrm{part}}$ induced by
all the other ER particles needs to be calculated by an iterative procedure.
Vesely (1977)

In the rheological literature, it is usual to ignore the polarization of the
particles by each other ($\bm{\mu}^{\mathrm{part}}_{i}=0$). There are,
however, important exceptions. Tao and Sun (1991); Tao (1993); Tao and Jiang
(1994a); Bonnecaze and Brady (1992); w. Wang  _et al._ (1996); Wang  _et al._
(1997); Martin  _et al._ (1998); Hynninen and Dijkstra (2005) In this work, we
present the full self consistent solution of Eqs. 6–7.

If we introduce the force exerted on dipole $\bm{\mu}_{i}$ by dipole
$\bm{\mu}_{j}$ (irrespective whether they are induced by $E^{\mathrm{appl}}$
or by other particles) is

$\displaystyle\mathbf{f}^{\mathrm{dip}}_{ij}(\mathbf{r}_{ij},\bm{\mu}_{i},\bm{\mu}_{j})=-(\bm{\mu}_{i}\cdot\nabla_{i})\mathbf{E}_{j}(\mathbf{r}_{i})=$
$\displaystyle=\frac{1}{4\pi\epsilon_{0}}\frac{1}{r_{ij}^{4}}\left\\{3\left[\bm{\mu}_{i}(\mathbf{n}_{ij}\cdot\bm{\mu}_{j})+\bm{\mu}_{j}(\mathbf{n}_{ij}\cdot\bm{\mu}_{i})+\right.\right.$
$\displaystyle+\left.\left.\mathbf{n}_{ij}(\bm{\mu}_{i}\cdot\bm{\mu}_{j})\right]-15\mathbf{n}_{ij}(\mathbf{n}_{ij}\cdot\bm{\mu}_{i})(\mathbf{n}_{ij}\cdot\bm{\mu}_{j})\right\\},$
(8)

we can express the force exerted on dipole $\bm{\mu}^{\mathrm{appl}}_{i}$ by
dipole $\bm{\mu}^{\mathrm{appl}}_{j}$ as

$\mathbf{f}_{ij}^{\mathrm{appl}}=\mathbf{f}^{\mathrm{dip}}_{ij}(\mathbf{r}_{ij},\bm{\mu}^{\mathrm{appl}}_{i},\bm{\mu}_{j}^{\mathrm{appl}})$
(9)

and the force exerted on dipole $\bm{\mu}^{\mathrm{appl}}_{i}$ by dipole
$\bm{\mu}^{\mathrm{part}}_{j}$ as

$\mathbf{f}_{ij}^{\mathrm{part}}=\mathbf{f}^{\mathrm{dip}}_{ij}(\mathbf{r}_{ij},\bm{\mu}^{\mathrm{appl}}_{i},\bm{\mu}_{j}^{\mathrm{part}}).$
(10)

The finite size of the ER particles is taken into account with a short-range
repulsive core potential for wich the cut & shifted Lennard Jones (LJ)
potential also known as the Weeks-Chandler-Anderson (WCA) potential is used.
The WCA force is defined as

$\mathbf{f}_{ij}^{\mathrm{WCA}}(\mathbf{r}_{ij})=\left\\{\begin{array}[]{ll}\mathbf{f}_{ij}^{\mathrm{LJ}}(\mathbf{r}_{ij})&\mathrm{if}\;\;r_{ij}<r_{\mathrm{c}}\\\
0&\mathrm{if}\;\;r_{ij}>r_{\mathrm{c}}\end{array}\right.,$ (11)

where

$\mathbf{f}_{ij}^{\mathrm{LJ}}(\mathbf{r}_{ij})=24\varepsilon^{\mathrm{LJ}}\left[2\left(\frac{d}{r_{ij}}\right)^{12}-\left(\frac{d}{r_{ij}}\right)^{6}\right]\frac{\mathbf{r}_{ij}}{r_{ij}^{2}}$
(12)

is the LJ force ($r_{\mathrm{c}}=2^{1/6}d$ is the cutoff distance at the
minimum of the LJ potential).

The trajectories of the particles can be computed from Langevin’s equations of
motion Lemons and Gythiel (1997)

$m\frac{d\mathbf{v}_{i}(t)}{dt}=\mathbf{F}_{i}\left(\mathbf{r}_{i}(t)\right)-m\gamma\mathbf{v}_{i}(t)+\mathbf{R}_{i}(t),$
(13)

where

$\mathbf{F}_{i}=\sum_{j}(\mathbf{f}_{ij}^{\mathrm{WCA}}+\mathbf{f}_{ij}^{\mathrm{appl}}+\mathbf{f}_{ij}^{\mathrm{part}})$
(14)

is the systematic force, $-m\gamma\mathbf{v}_{i}(t)$ is the frictional force,
$\mathbf{R}_{i}(t)$ is the random force, $\mathbf{r}_{i}$, $\mathbf{v}_{i}$,
$m$, and $\gamma$ are the position, the velocity, the mass, and the friction
coefficient of particle $i$, respectively. The friction coefficient can be
computed from Stokes’ law as $\gamma=3\pi\eta d/m$.

To solve this stochastic differential equation we use the GJF-2GJ version
Jensen and Grønbech-Jensen (2019) of a collections of algorithms proposed by
Grønbech-Jensen and Farago: Grønbech-Jensen and Farago (2013); Farago (2019);
Jensen and Grønbech-Jensen (2019)

$v^{n+\frac{1}{2}}=av^{n-\frac{1}{2}}+\frac{\sqrt{b}\Delta
t}{m}f^{n}+\frac{\sqrt{b}}{2m}\left(R^{n}-R^{n+1}\right)$ (15)
$r^{n+1}=r^{n}+\sqrt{b}v^{n+\frac{1}{2}}\Delta t,$ (16)

where $r^{n}=r(t^{n})$ is any position coordinate of any particle,
$v^{n}=v(t_{n})$ is any velocity coordinate of any particle, $t^{n}=n\Delta t$
is the time in the $n$th time step, $\Delta t$ is the time step,
$a=(1-\gamma\Delta t/2)/(1+\gamma\Delta t/2)$ , $b=(1)/(1+\gamma\Delta t/2)$,
$t_{n+\frac{1}{2}}=t_{n}+\Delta t/2$, and $t_{n-\frac{1}{2}}=t_{n}-\Delta
t/2$. The discrete time noise, $R^{n}$, is a random Gaussian number with
properties $\langle R^{n}\rangle=0$ and $\langle R^{m}R^{n}\rangle=2kT\gamma
m\Delta t\delta_{mn}$ with $\delta_{mn}$ being the Kronecker-delta.

The Brownian Dynamics simulations have been performed in a cubic simulation
cell using periodic boundary conditions. The ensemble can be considered
canonical because $V$ and $N$ are fixed, while the temperature is also
constant because the system is thermostated by the Langevin integrator via the
fluctuation–dissipation theorem. The dipolar interactions were truncated at
the half of the cell width. System size dependence has been analyzed in our
previous work. Fertig  _et al._ (2021)

##  III Results and Discussion

Table 1: Reduced quantities defined with $T$, $d$, $m$ or with $T$, $d$, $\rho_{\mathrm{in}}$. Quantity | Reduced quantity
---|---
Time | $t^{*}=t\sqrt{kT/md^{2}}=t\sqrt{6kT/\pi\rho_{\mathrm{in}}d^{5}}$
Distance | $r^{*}=r/d$
Density | $\rho^{*}=\rho d^{3}$
Velocity | $v^{*}=v\sqrt{m/kT}=v\sqrt{\pi\rho_{\mathrm{in}}d^{3}/6kT}$
Energy | $u^{*}=u/kT$
Force | $F^{*}=Fd/kT$
Electric field | $E^{*}=E\sqrt{4\pi\epsilon_{0}d^{3}/kT}$
Dipole moment | $\mu^{*}=\mu/\sqrt{4\pi\epsilon_{0}kTd^{3}}$
Polarizability | $\alpha^{*}=\alpha/4\pi\epsilon_{0}d^{3}$
Friction coefficient | $\gamma^{*}=\gamma\sqrt{md^{2}/kT}=\gamma\sqrt{\rho_{\mathrm{in}}d^{5}/6kT}$

In this work, we use reduced units that are collected in Table 1. They express
physical quantities as dimensionless numbers obtained by dividing a quantity
in a physical unit by a unit quantity, $r^{*}=r/d$, for example. In addition
to $T$ and $d$, we can use either $m$ or $\rho_{\mathrm{in}}$ (mass density of
the ER particles) to define the reduced quantities, because the mass depends
on $m$ and $d$ through $m=\rho_{\mathrm{in}}d^{3}\pi/6$. The diffusion
constant in the high coupling limit can be expressed by Einstein’s relation:
$D=kT/m\gamma$. The square of the reduced dipole moment

$(\mu^{*})^{2}=\frac{\mu^{2}/4\pi\epsilon_{0}d^{3}}{kT}$ (17)

is an important parameter because it relates the ordering effect of the
dipolar energy to the disordering effect of thermal motion. It is proportional
to the $\lambda$ parameter used in the literature. If $(\mu^{*})^{2}$ is
large, the dipolar interactions are strong enough to induce chain formation,
while if it is small, thermal motion prevents chain formation.

In this work, our main goal is to study the dielectric response of the ER
fluid and to analyze the relationship of this response to chain formation in
the system. The relative permittivity of an ER fluid can be computed from the
corrected Clausius-Mosotti (CM) equation that can be derived from a
polarization formula: Neumann (1983)

$\dfrac{\epsilon-1}{\epsilon+2}=\dfrac{1}{3\epsilon_{0}}\dfrac{\left\langle
P\right\rangle}{E^{\mathrm{appl}}},$ (18)

where $P$ is the polarization density obtained from the sum of the dipoles
$\mu^{\mathrm{appl}}=\alpha E^{\mathrm{appl}}$ induced directly by the
external field and the average dipole moment
$\left\langle\mu^{\mathrm{part}}\right\rangle$ induced by other particles:

$\left\langle
P\right\rangle=\frac{N\mu^{\mathrm{appl}}+N\left\langle\mu^{\mathrm{part}}\right\rangle}{V}=\rho\mu^{\mathrm{appl}}\left(1+\frac{\left\langle\mu^{\mathrm{part}}\right\rangle}{\mu^{\mathrm{appl}}}\right),$
(19)

here $V$ is the considered volume, and $N$ the number of particles in it. Eqs.
4, 18 and 19 result in the corrected CM equation

$\dfrac{\epsilon-1}{\epsilon+2}=\dfrac{\alpha\rho}{3\epsilon_{0}}\left(1+\dfrac{\left\langle\mu^{\mathrm{part}}\right\rangle}{\mu^{\mathrm{appl}}}\right),$
(20)

where $\rho=N/V$ is the number density, and the correction factor
$\left\langle\mu^{\mathrm{part}}\right\rangle/\mu^{\mathrm{appl}}$ is the
average induced dipole due to other particles normalized by the dipole due to
the external field. The quantity
$\left\langle\mu^{\mathrm{part}}\right\rangle$ is directly provided by the
simulations. We have investigated the correction to the CM equation in the
case of nonpolar fluids (e.g., carbon dioxide) by Monte Carlo simulations.
Valiskó and Boda (2009)

In the low-field-strength limit, the original CM equation is recovered (the
correction factor is zero), because formally an ensemble of ER particles
corresponds to an ensemble of non-polar, but polarizable molecules. The CM
equation is based on the Lorentz formula Lorentz (1909) for the internal field
and ignores the fact that a particle is also polarized by other particles not
only by the external field. Keyes and Kirkwood (1931)

The number of particles is fixed at $N=256$ in our simulations. This number
proved to be sufficient on the basis of the system size analysis provided in
our previous study. Fertig  _et al._ (2021)

The value of the reduced friction coefficient is fixed at $\gamma^{*}=100$.
This value made simulations feasible, because the system developed relatively
fast at the fixed value of the time step $\Delta t^{*}=0.005$. The effect of
larger values of $\gamma^{*}$ that are more characteristic of typical ER
fluids is analyzed in our previous work. Fertig  _et al._ (2021) The reduced
density is fixed at $\rho^{*}=0.05$. The parameters that we change are the
reduced electric field $E^{*}$ and the reduced polarizability $\alpha^{*}$.

To characterize time dependence, we show values of block averages (denoted by
$\langle\dots\rangle_{\mathrm{b}}$ for various physical quantities as
functions of $t^{*}$. In this work, we performed $M_{\mathrm{b}}{=}5000$ time
steps in a block. We perform $M_{0}=50,000$ time steps ($20$ blocks) in the
absence of applied electric field ($E^{\mathrm{appl}}{=}0$), after which the
electric field is instantaneously switched on. Then we performed
$M_{\mathrm{E}}=450,000$ time steps ($180$ blocks) in the presence of a
constant applied field.

Such a cycle was started over and done several times and averaged to smooth
out noise. When we start a cycle over, we restart from a freshly generated
initial configuration in a completely disordered state without chains. This
way, the subsequent periods are independent and can be averaged.

Fig. 1 shows the time dependence of $\Delta\epsilon$ for $\alpha^{*}=0.03$
(the curves for other values of $\alpha^{*}$ are similar). The numbers near
the curves indicate the values of the reduced electric field, $E^{*}$. This
figure corresponds to Fig. 3 of Horváth and Szalai Horváth and Szalai (2012)
whose data are reproduced in the inset of Fig. 1 for qualitative comparison.

Figure 1: Time dependence of $\Delta\epsilon$ for $\alpha^{*}=0.03$. The
symbols with error bars are simulated data. The error bars have been computed
from the variance of the data in the consecutive and independent periods. The
lines are bi-exponential fits (Eq. 3). The numbers near the curves indicate
the values of the reduced electric field, $E^{*}$. The inset shows the
experimental data of Horváth and Szalai Horváth and Szalai (2012) for
comparison. The numbers in the inset indicate electric field strengths in MV/m
unit.

Direct quantitative comparison with the experimental data is problematic (see
the end of this section) due mainly to the fact that the ER fluid studied by
Horváth and Szalai is polydisperse. They considered nanosized ($10-20$ nm)
silica (SiO2) particles dispersed in silicone oils (polydimethylsiloxane) with
different dynamic viscosities ($0.34$ and $0.97$ Pa s). Qualitative
comparison, however, is possible. By fitting a bi-exponential to our simulated
data (Eq. 3), we obtain the characteristic times (in reduced units) shown in
Fig. 2. Panel A shows the data as functions of $E^{*}$ for different values of
$\alpha^{*}$.

This figure implies that a larger electric field is needed to achieve smaller
$\tau^{*}$ values (faster processes) in the case of smaller $\alpha^{*}$
values. This is obvious, because the important parameter from the point of
view of the dipolar interactions is the induced dipole that is
$\alpha^{*}E^{*}$. Therefore, if we plot $\tau_{1}^{*}$ and $\tau_{2}^{*}$ as
functions of $\alpha^{*}E^{*}$, we obtain a scaling behavior: the curves for
different $\alpha^{*}$ values collapse onto a single curve (Fig. 2B). This
scaling behavior also applies for the time dependence of the normalized
$\Delta\epsilon$; if we plot
$\Delta\epsilon(t^{*})/\Delta\epsilon(t^{*}\rightarrow\infty)$ as a function
of $t^{*}$ for a fixed value of $\alpha^{*}E^{*}$ but for different
combinations of $\alpha^{*}$ and $E^{*}$, the curves collapse onto a single
one.

Figure 2: The characteristic times $\tau_{1}$ and $\tau_{2}$ (the latter is
larger with an order of magnitude) as functions of (A) $E^{*}$ and (B)
$\alpha^{*}E^{*}$ for various values of $\alpha^{*}$. The inset shows the
experimental data of Horváth and Szalai Horváth and Szalai (2012) (from their
Table I.) for the two different values of the viscosity they considered (black
and red colors refer to $\eta=0.34$ and $0.97$ Pa s, respectively). The error
bars of $\tau_{1}$ and $\tau_{2}$ estimated from the Levenberg-Marquardt
algorithm Marquardt (1963) are within the size of the sysmbols.

This figure also shows that $\tau_{2}^{*}$ is an order of magnitude larger
than $\tau_{1}^{*}$ in agreement with the predictions of Hass et al. Hass
(1993) as well as with the experiments of Ly et al. Ly  _et al._ (2001) and
Horváth and Szalai. Horváth and Szalai (2012)

The qualitative behavior of the $\tau_{1}$ vs. $E$ function is also similar to
the experimental behavior as seen from comparison to the data in Table I of
Horváth and Szalai Horváth and Szalai (2012) that are reproduced in the inset
of Fig. 2B.

Figure 3: The equilibrium ($t\rightarrow\infty$) values of the one-particle
dipolar energy, $(u^{\mathrm{dip}})^{*}$ (top-left panel), the diffusion
constant, $D^{*}$ (top-right panel), the induced permittivity increment,
$\Delta\epsilon$ (bottom-left panel), and the average chain length,
$s_{\mathrm{av}}$ (bottom-right panel), for various values of $\alpha^{*}$.
The vertical dashed lines indicate the $E_{\mathrm{h}}^{*}$ field strengths at
which the regions with different slopes meet (see Fig. 4 for more
explanation). The error bars are those computed at $t^{*}=4500$ from the
variance of the values in the periods. The error bars for
$(u^{\mathrm{dip}})^{*}$ and $\Delta\epsilon$ are within the size of the
symbols.

The purpose of this study is to look into the black box and to see how the
dielectric behavior observed in experiment and simulation is related to the
structural changes occurring in the ER fluid as the electric field is
increased. These structural changes can be monitored via various physical
quantities such as the diffusion constant, the dipolar energy, and the average
chain length. Therefore, we plot the equilibrium values of these quantities
(along with $\Delta\epsilon$) as functions of the field strength, $E^{*}$, in
Fig. 3. The equilibrium value of a quantity can be obtained either by running
a long simulation and throwing the equilibration period away or from
substituting $t\rightarrow\infty$ into the bi-exponential. We have chosen the
second option. Fig. 1 shows the fitted functions. The $R^{2}$ coefficient of
the fitting is generally above $0.99$ and the residuals (data not shown)
decrease to zero as $t\rightarrow\infty$. These imply that the equilibrium
values obtained this way are trustable.

The equilibrium values of $\Delta\epsilon$ as functions of $E^{*}$ are shown
in the bottom-left panel of Fig. 3. Also, the curve for $\alpha^{*}=0.03$ is
reproduced in Fig. 4 along with the experimental data in the insets. The
qualitative agreement of the simulations and experimental data is apparent.
The behavior of our curve follows the behavior of the experimental curve in
the respect that there is a initial stage with a small slope. When the
electric field is larger than a threshold value ($E_{1}$), pairs and short
chains start to form and the curve in the interval between $E_{1}$ and the
next threshold value $E_{\mathrm{h}}$ has a larger slope. Above
$E_{\mathrm{h}}$, the chains spanning the simulation cell start to aggregate
with a smaller slope of the $\Delta\epsilon$ vs. $E^{*}$ function.

One notable difference between the simulation and experimental results is that
$E_{1}$ is much smaller relative to $E_{\mathrm{h}}$ in the experiment than in
the model. This is valid for all $\alpha^{*}$ values studied. The reason of
this is not clear, but the results may be system-size and density dependent.
It is also possible that the large-particle fraction of the experimental
polydisperse system can form clusters at lower fields than the average-sized
particles, an effect that is absent in our monodisperse model. The
$E_{\mathrm{h}}$ value, however, appears to be a relatively well defined point
separating two characteristic regions with different slopes, so we indicate
the $E^{*}_{\mathrm{h}}$ vales with vertical dashed lines in all the panels of
Fig. 3.

Note that $\Delta\epsilon(E^{*})$ function exhibits a steep nonlinear increase
at large $E^{*}$ values. This is the result of the strong dipolar attraction
overriding the repulsion of the WCA potential. The overlapping spheres lead to
stronger particle-particle polarization. We consider this behavior an artifact
of the model.

Figure 4: The equilibrium ($t\rightarrow\infty$) value of $\Delta\epsilon$ as
a function of $E^{*}$ for $\alpha^{*}=0.03$. The inset shows the experimental
data from Fig. 4 of Horváth and Szalai Horváth and Szalai (2012). The inset of
the inset shows the results for small $E^{*}$. The solid blue lines indicate
the slopes, while dashed blue lines indicate the crosses of the solid lines
defining the threshold field strengths $E_{1}$ and $E_{\mathrm{h}}$.

Next, let us see how the behavior of the other physical quantities correlates
with the behavior of $\Delta\epsilon$. The one-particle dipolar energy,
$(u^{\mathrm{dip}})^{*}$, does not exhibit the behavior of the three (small-
large-small) slopes (top-left panel of Fig. 3). It decreases with increasing
$E^{*}$ at a continuously increasing rate. The explanation is that the dipolar
energy chiefly depends on the interactions inside the chains. If we increase
$E^{*}$, the dipoles become larger, and also their interactions. Clustering of
chains does not seriously influence this dependence.

The diffusion constant is calculated as the slope of the mean square
displacement (MSD) as a function of time:

$D(t_{\mathrm{b}})=\frac{\langle\mathbf{r}^{2}(t)\rangle_{\mathrm{b}}}{6\Delta
t_{\mathrm{b}}},$ (21)

where $t_{\mathrm{b}}$ is the time at the beginning of a block, and $\Delta
t_{\mathrm{b}}{=}M_{\mathrm{b}}\Delta t$ is the length of the block. This way,
$D$ is characteristic of a block and time dependence can be studied. $D^{*}$
decreases as $E^{*}$ increases as shown in Fig. 3 (top-right panel). The ER
particles lose their mobilities as they are organized into chain-like and
columnar structures. The behavior of $D^{*}$ follows the behavior of
$\Delta\epsilon$. Its value starts with the $E^{*}\rightarrow 0$ limit
($0.01$), breaks down around $E_{1}$, and decreases steeply as longer chains
are formed at higher $E^{*}$ values. Around $E_{\mathrm{h}}$, the chains
aggregate, $D^{*}$ decreases at a lower rate, and saturates into a very small,
but non-zero value. At large $E^{*}$, the spheres move together with their
chains that have much smaller mobility than the single spheres.

These results are closely related to the anomalous diffusion behavior of
dipolar chains that has been studied theoretically and experimentally. Furst
and Gast (2000); Toussaint  _et al._ (2004, 2006) Those studies imply that the
diffusion of chains is reduced compared to the case in the absence of chains.
The degree of reduction is related to the average length of the chains.

Figure 5: Equilibrium limits of the chain length distribution, $n_{s}$ (left),
and radial distribution functions, $g(r^{*})$ (right) for $\alpha^{*}=0.03$.
The curves are obtained by averaging over $10$ blocks at the end of $M_{E}$
simulation periods and averaging over periods. Snapshots are shown in the
middle in front and top view. The electric field strength increases from top
to bottom.

The average chain length, $s_{\mathrm{av}}$, is computed by identifying the
number of chains of length $s$, $n_{s}$, for every configuration, taking the
average

$s_{\mathrm{av}}=\frac{\sum_{s}sn_{s}}{\sum_{s}n_{s}},$ (22)

and then averaging over configurations. Chain length, $s$, is measured in the
number of particles in the chain. Two particles are defined to be in the same
chain if their distance is smaller than $1.2d$. The choice of $1.2$ does not
influence our qualitative conclusions for the dynamics of chain formation.
Other definitions of chains were analyzed in our previous study Fertig  _et
al._ (2020).

The average chain length starts to increase only above the first threshold
value, $E_{1}^{*}$. Above the second threshold value, $E^{*}_{\mathrm{h}}$,
the average length of chains reaches the value $s\approx 18$ that corresponds
to a chain spanning the simulation cell whose length is $L\approx 17.23d$.

As the electric field is increased further well above $E^{*}_{\mathrm{h}}$,
the average number of chains also increases and eventually reaches a limiting
value of 64. At large electric fields characteristically $4$ clusters of
chains are formed each containing on average 64 spheres, but this is just an
average. There is thermal motion, so values different from 64 may occur, but
the average seems solid. Thermal motion at strong field strengths mainly means
the translational motion of chains in the lateral ($x,y$) plane and rotation
of chains about the axes of the chains. Also, the $4$ clusters of chains do
not aggregate further. At high electric fields, the repulsion between these
clusters seems to prevent further aggregation. This finding, however, is
density dependent too.

It is of interest of looking into the black box even deeper to see how these
average chain length values come about. It is done in Fig. 5 in which we show
chain length distributions, $n_{s}$, radial distribution functions, $g(r)$,
and snapshots. The electric field strength increases from top to bottom. The
values are chosen to show the phases in Figs. 3 and 4. We show results for a
value below $E^{*}_{1}$, around $E^{*}_{1}$, between $E^{*}_{1}$ and
$E^{*}_{\mathrm{h}}$, around $E^{*}_{\mathrm{h}}$, and above
$E^{*}_{\mathrm{h}}$.

As $E^{*}$ increases, the chain length distributions show the increased
probability of longer chains. The snapshots clearly show these chains that are
also indicated by peaks in the $g(r)$ functions.

For $E^{*}=33.33$ (below $E^{*}_{1}$), the system is practically a homogeneous
isotropic gas-like fluid regarding the ER particles. For $E^{*}=52.7$ (around
$E^{*}_{1}$), pair formation and, to some degree, formation of short chains
are present. For $E^{*}=57.74$ (between $E^{*}_{1}$ and $E^{*}_{\mathrm{h}}$),
even longer chains and, accordingly, more peaks in $g(r)$ appear.

When we reach $E^{*}_{\mathrm{h}}$ ($E^{*}=66.67$), chains spanning the
simulation cell are clearly present indicated by the peak at $s=18$ in the
$n_{s}$ function. This chain is more stabilized by the periodic boundary
conditions in our cubic simulation cell of length $L\approx 17.23d$ compared
to other chains. In our previous study, Fertig  _et al._ (2021) we analyzed
this behavior in detail.

For even larger electric field strength ($E^{*}=78.17$), we find another peak
at $s=36$ that corresponds to two chains stuck together. Chains are
straighter, and the peaks in $g(r)$ are more pronounced. Also, the $n_{s}$
function is more noisy than at lower electric field strengths that indicates
that the system is “more frozen” or “less fluid”. The evolution of the system
is determined by the movements of the much less mobile chains instead of the
movements of individual particles and short chains.

Figure 6: Snapshots for a very large electric field, $E^{*}=149.07$, for
$\alpha^{*}=0.03$.

If we increase the electric field even further ($E^{*}=149.07$), the chains
aggregate into column-like structures. The resulting $n_{s}$ and $g(r)$
profiles are even more noisy and not very meaningful, so we show only
snapshots in Fig. 6.

Finally, we briefly consider the possibility of a quantitative comparison with
the experimental results. If we accept the hypothesis that the meaning of
$E_{\mathrm{h}}$ is the same in the model and in the experiment (it is a
cornerstone), we can estimate the particle diameter, $d$, that brings the
simulation and the experimental data into correspondence. If we take the value
of $E_{\mathrm{h}}=1.18$ MV/m from Fig. 4 of Ref. Horváth and Szalai, 2012 and
relate it to the $E_{\mathrm{h}}^{*}$ values depicted from Fig. 3
($E_{\mathrm{h}}^{*}=52.5$, $72$, $107$ for $\alpha^{*}=0.02$, $0.03$, and
$0.04$, respectively), we obtain the values $d=673$, $517$, and $418$ nm for
$\alpha^{*}=0.02$, $0.03$, and $0.04$, respectively. These values are much
larger than the $10-20$ nm values specified in the paper of Horváth and
Szalai. Horváth and Szalai (2012)

We can explain this in different ways. First, the ER particles provided by the
manufacturer are polydisperse as opposed to our model that includes particles
of the same size. Also, it was observed in the experiments that the particles
tend to stick together forming larger particles resulting in a larger
effective diameter. Water content that may increase polarizability cannot be
excluded. In any case, the values $10-20$ nm are so small that using them in a
simulation (by transforming them to the corresponding reduced units) does not
result in any kind of chain formation. For these reasons, we regard the
diameters calculated and reported here more realistic than the $10-20$ nm
values.

It is also possible to correspond the $\tau_{1}^{*}$ value depicted from Fig.
2 ($\tau_{1}^{*}\approx 68.5$ for $\alpha^{*}=0.03$) to the experimental
$\tau_{1}$ value depicted from the inset of Fig. 2B ($\tau_{1}\approx 0.38$ s
for viscosity $0.97$ Pa s). The resulting diameter is $d\approx 2470$ nm. The
reason of these large values is that our reduced friction coefficient is small
($\gamma^{*}=100$). Because we can extrapolate to larger values of
$\gamma^{*}$ (see Fig. 6 of our previous work Fertig  _et al._ (2021)), we can
provide the estimation that by increasing $\gamma^{*}$ with two orders of
magnitude, the resulting diameter decreases with about one order of magnitude.
Changing $\gamma^{*}$ does not change the equilibrium value of
$s_{\mathrm{av}}$ and $\Delta\epsilon$, it only influences how fast the system
converges to these values.

To relate $t^{*}$ to $t$ and $E^{*}$ to $E$ in Figs. 1, 2 and 4, we collect
the unit values $E_{0}=E/E^{*}$ and $t_{0}=t/t^{*}$ in Table 2 for three
representative values of the particle diameter: $d=517$ and $2467$ nm obtained
from the procedures of relating $E_{\mathrm{h}}$ to $E_{\mathrm{h}}^{*}$ and
$\tau_{1}$ to $\tau_{1}^{*}$, respectively, and a value in between ($1000$
nm).

Table 2: Unit values of the electric field strength and particle diameter, $E_{0}=E/E^{*}$ and $t_{0}=t/t^{*}$, respectively, for different values of $d$ on the basis of Table 1 ($T=298.15$ K and $\rho_{\mathrm{in}}=2650$ kg/m3). $d$ / nm |  $E_{0}$ / MVm-1 |  $t_{0}$ / s
---|---|---
$516$ | $0.016$ | $0.00011$
$1000$ | $0.0061$ | $0.00058$
$2467$ | $0.0016$ | $0.0055$

##  IV Conclusions

We performed Brownian Dynamics simulations for ER fluids by taking the cross-
polarization among particles into account in a self consistent way and
computed the induced dielectric increment, $\Delta\epsilon$, as a function of
time after the applied field is switched on and as a function of field
strength at the equilibrium limit. Particle-particle polarization is essential
for computing $\Delta\epsilon$ that is a very useful quantity because it is
measurable well and also can be obtained from simulations with a small
statistical error.

Our results are in qualitative agreement with the experimental results of
Horváth and Szalai Horváth and Szalai (2012) and relate the computed data to
structural features in terms of energy, diffusion constant, average chain
length, chain length distribution, radial distribution functions, and
snapshots. The hypotheses about the correlations between dielectric properties
and structural features put forward by Horváth and Szalai have been confirmed
by our calculations.

## Author’s contribution

All authors contributed equally to this work.

## Acknowledgments

This research was supported by the European Union, co-financed by the European
Social Fund, via the project “Research of autonomous vehicle systems related
to the autonomous test track in Zalaegerszeg” (EFOP-3.6.2-16-2017-00002). We
also acknowledge the support of the National Research, Development and
Innovation Office (NKFIH), project No. K124353. We acknowledge KIFÜ for
awarding us access to resource based in Hungary at Szeged.

## Data avalability

The data that supports the findings of this study are available within the
article.

## References

  * Winslow (1949) W. M. Winslow, J. Appl. Phys. 20, 1137 (1949).
  * Duclos  _et al._ (1992) T. G. Duclos, J. D. Carlson, M. J. Chrzan, and J. P. Coulter, in  _Solid Mechanics and Its Applications_ (Springer Netherlands, 1992) pp. 213–241.
  * Havelka and Filisko (1995) K. O. Havelka and F. E. Filisko, eds.,  _Progress in Electrorheology_ (Springer US, 1995).
  * Adolf and Garino (1995) D. Adolf and T. Garino, Langmuir 11, 307 (1995).
  * Tao and Roy (1994) R. Tao and G. Roy,  _Electrorheological Fluids: Mechanisms, Properties, Technology, and Applications_ (World Scientific, 1994).
  * Wen  _et al._ (1997) W. Wen, S. Men, and K. Lu, Phys. Rev. E 55, 3015 (1997).
  * Wen  _et al._ (1998) W. Wen, H. Ma, W. Y. Tam, and P. Sheng, Applied Physics Letters 73, 3070 (1998).
  * Rzoska and Zhelezny (2006) S. Rzoska and V. Zhelezny,  _Nonlinear Dielectric Phenomena in Complex Liquids_, Nato Science Series II: (Springer Netherlands, 2006).
  * Horváth and Szalai (2012) B. Horváth and I. Szalai, Phys. Rev. E 86, 061403 (2012).
  * Baxter-Drayton and Brady (1996) Y. Baxter-Drayton and J. F. Brady, J. Rheology 40, 1027 (1996).
  * Klingenberg  _et al._ (1989) D. J. Klingenberg, F. van Swol, and C. F. Zukoski, J. Chem. Phy. 91, 7888 (1989).
  * See and Doi (1991) H. See and M. Doi, J. Phys. Soc. Japan 60, 2778 (1991).
  * Toor (1993) W. R. Toor, J. Colloid Interf. Sci. 156, 335 (1993).
  * Hass (1993) K. C. Hass, Phy. Rev. E 47, 3362 (1993).
  * Climent  _et al._ (2004) E. Climent, M. R. Maxey, and G. E. Karniadakis, Langmuir 20, 507 (2004).
  * Domínguez-García  _et al._ (2007) P. Domínguez-García, S. Melle, J. M. Pastor, and M. A. Rubio, Phys. Rev. E 76, 051403 (2007).
  * Tao and Jiang (1994a) R. Tao and Q. Jiang, Phys. Rev. Lett. 73, 205 (1994a).
  * Tao and Jiang (1994b) R. Tao and Q. Jiang, Int. J. Modern Phys. B 08, 2721 (1994b).
  * Enomoto and Oba (2002) Y. Enomoto and K. Oba, Physica A 309, 15 (2002).
  * Whittle (1990) M. Whittle, J. Non-Newton. Fluid 37, 233 (1990).
  * Heyes and Melrose (1990) D. M. Heyes and J. R. Melrose, Mol. Sim. 5, 293 (1990).
  * Cao  _et al._ (2006) J. G. Cao, J. P. Huang, and L. W. Zhou, J. Phys. Chem. B 110, 11635 (2006).
  * Bonnecaze and Brady (1992) R. T. Bonnecaze and J. F. Brady, J. Chem. Phys. 96, 2183 (1992).
  * Jackson (1999) J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, New York, 1999).
  * Vesely (1977) F. J. Vesely, J. Comp. Phys. 24, 361 (1977).
  * Tao and Sun (1991) R. Tao and J. M. Sun, Phys. Rev. Lett. 67, 398 (1991).
  * Tao (1993) R. Tao, Phys. Rev. E 47, 423 (1993).
  * w. Wang  _et al._ (1996) Z. w. Wang, Z. f. Lin, and R. b. Tao, Int. J. Mod. Phys. B 10, 1153 (1996).
  * Wang  _et al._ (1997) Z.-w. Wang, Z.-f. Lin, and R.-b. Tao, Chin. Phys. Lett. 14, 151 (1997).
  * Martin  _et al._ (1998) J. E. Martin, R. A. Anderson, and C. P. Tigges, J. Chem. Phys. 108, 3765 (1998).
  * Hynninen and Dijkstra (2005) A.-P. Hynninen and M. Dijkstra, Phys. Rev. Lett. 94 (2005), 10.1103/physrevlett.94.138303.
  * Lemons and Gythiel (1997) D. S. Lemons and A. Gythiel, Am. J. Phys. 65, 1079 (1997).
  * Jensen and Grønbech-Jensen (2019) L. F. G. Jensen and N. Grønbech-Jensen, Mol. Phys. 117, 2511 (2019).
  * Grønbech-Jensen and Farago (2013) N. Grønbech-Jensen and O. Farago, Mol. Phys. 111, 983 (2013).
  * Farago (2019) O. Farago, Physica A 534, 122210 (2019).
  * Fertig  _et al._ (2021) D. Fertig, D. Boda, and I. Szalai, AIP Advances 11, 025243 (2021).
  * Neumann (1983) M. Neumann, Mol. Phys. 50, 841 (1983).
  * Valiskó and Boda (2009) M. Valiskó and D. Boda, J. Chem. Phys. 131, 064120 (2009).
  * Lorentz (1909) H. A. Lorentz,  _Theory of Electrons_ (B. G. Teubner, Leipzig, 1909).
  * Keyes and Kirkwood (1931) F. G. Keyes and J. G. Kirkwood, Phys. Rev. 37, 202 (1931).
  * Marquardt (1963) D. W. Marquardt, J. Soc. Ind. Appl. Math. 11, 431 (1963).
  * Ly  _et al._ (2001) H. V. Ly, K. Ito, H. T. Banks, M. R. Jolly, and F. Reittich, Int. J. Mod. Phys. B 15, 894 (2001).
  * Furst and Gast (2000) E. M. Furst and A. P. Gast, Phys. Rev. E 62, 6916 (2000).
  * Toussaint  _et al._ (2004) R. Toussaint, G. Helgesen, and E. G. Flekkøy, Phys. Rev. Lett. 93 (2004), 10.1103/physrevlett.93.108304.
  * Toussaint  _et al._ (2006) R. Toussaint, E. G. Flekkøy, and G. Helgesen, Phys. Rev. E 74 (2006), 10.1103/physreve.74.051405.
  * Fertig  _et al._ (2020) D. Fertig, D. Boda, and I. Szalai, Hung. J. Ind. Chem. 48, 95 (2020).