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

Waiting time is a common source of waste in business processes~\cite{delias2017positive}.
Although it is impractical to completely eliminate waiting times in business processes, there are various approaches to reduce them~\cite{jansen2005business}. Waiting times typically arise during transitions between activities, i.e., when the execution of a case moves from one activity to another. There are different reasons why waiting times occur during activity transitions. For instance, when two consecutive activities in a case are executed by different resources (a.k.a.\ a handoff~\cite{van2007business}), the processing of the case is put on hold until the next resource becomes available to execute it. In this scenario, the cause of the waiting time is resource contention, i.e.\ a resource is not available to execute an activity instance because they are busy executing other activity instances~\cite{Dumas2013fundamentals}. Waiting times may also be caused by data exchanges~\cite{ramakrishnan2008defining}, coordination issues~\cite{rummel2005activity}, or synchronization points~\cite{rummel2005activity}.

Process mining techniques allow us to analyze data generated by business process executions, a.k.a.\ \emph{event logs}, to unveil performance and conformance issues, and associated improvement opportunities~\cite{Aalst2016}.
In particular, process mining techniques support the discovery of sources of waste, including waiting times~\cite{Aalst2016}.
However, while existing process mining techniques enable analysts to visualize activity transitions with high waiting time (i.e., bottlenecks), they provide limited support for identifying the causes of waiting times and how to reduce them.

To tackle this gap, this paper addresses the following research questions (RQ): (1) \textit{``What causes of waiting times between pairs of activity instances can be identified from event logs?''}, (2) \textit{``How can these causes of waiting time be identified from event logs?''}, and (3) \textit{``How can improvement opportunities, expressed as inefficiencies due to waiting times, be identified from event logs?''}

The contribution of the paper is two-fold.
First, we conceptualize the causes of waiting time arising from activity transitions.
Second, we propose a process mining approach to (1) discover waiting times associated with activity transitions; (2) identify their causes; and (3) analyze their impact w.r.t.\ a well known measure of temporal efficiency called Cycle Time Efficiency (CTE): the ratio of processing time to cycle time in a process~\cite{Dumas2013fundamentals}.

The proposed approach has been implemented as a software tool and  empirically evaluated to verify its ability to discover different causes of waiting time using synthetic event logs. In addition, the applicability of the approach has been tested by applying it to a real-life event log.

The rest of the paper is structured as follows.
Section~\ref{sec:related-work} introduces background and related work.
Section~\ref{sec:discovery-analysis} presents the proposed approach.
Section~\ref{sec:evaluation} describes the empirical evaluation, and Section~\ref{sec:conclusion} concludes the paper.

\section{Background and Related Work\label{sec:related-work}}

In this section, we introduce relevant conceptual foundations and notations from the fields of business process management and process mining, and position our contribution w.r.t.\ existing approaches to discover and analyze waiting times.

\subsection{Business Processes and Temporal Performance Measures}

A business process is a collection of events, activities, and decisions that lead from a customer need to an outcome that is of value to this customer~\cite{Dumas2013fundamentals}. Each execution of a business process is called a \emph{case}.
A common measure of temporal performance of a business process is its \emph{cycle time}: The time between the moment a case of the process starts, and the case ends, aggregated to the level of the set of cases of a process observed during a period of time.
The cycle time of a process consists of \emph{processing time} (i.e., the time during which a case is being processed) and \emph{waiting time} (i.e., the time when a case waits to be processed)~\cite{Dumas2013fundamentals}.

Waiting time may be caused by resource contention, i.e., no suitable resource is available to execute an activity instance~\cite{jansen2005business}.
Other causes of waiting time include: synchronization between resources within a process~\cite{rummel2005activity} or across multiple processes~\cite{Dumas2013fundamentals}, coordination between resources executing different activities~\cite{rummel2005activity}, data transfer~\cite{Dumas2013fundamentals}, batching~\cite{lashkevich2022data}, handoffs~\cite{rummel2005activity}, and external inputs~\cite{AndrewsW17}.

As stated earlier, the temporal efficiency of a process can be measured by its CTE: the ratio of processing time to cycle time.
When CTE is close to 1, there is relatively little waiting time. Conversely, if the CTE is close to 0, the waiting times are longer relative to processing times and there are opportunities to improve the CTE by reducing waiting times~\cite{Dumas2013fundamentals}.

\subsection{Event Logs and Activity Instance Logs}

Modern IT systems record and store process execution data in \emph{event logs}, i.e., sets of timestamped events capturing the execution of the activities in a process~\cite{Dumas2013fundamentals}.
An event log contains information about state changes of each activity instance (e.g.\ enablement, start, end, or cancellations of activity instances). In this paper, we use the concept of \emph{activity instance log} to represent the execution of a set of cases in a process.
An activity instance log is an event log in which each entry contains information about the start time and end time of an activity instance~\cite{Martin2021}. Below, we introduce several notations used in the paper, leading to a definition of an activity instance log.

We consider a business process that involves a set of \emph{activities} $A$.
We denote each of these activities with $\alpha$.
An \emph{activity instance} $\varepsilon = (\varphi, \alpha, \tau_e, \tau_s, \tau_c, \rho)$ denotes one execution of activity $\alpha$, where $\varphi$ identifies the case (represented by a \emph{process trace}) to which this execution belongs to, $\tau_e$, $\tau_s$, and $\tau_c$ denote, respectively, the instants in time in which this activity instance was \emph{enabled}, \emph{started}, and \emph{completed}, and $\rho$ identifies the \emph{resource} that processed the activity.
Accordingly, we use $\varphi(\varepsilon_{i})$, $\alpha(\varepsilon_{i})$, $\tau_e(\varepsilon_{i})$, $\tau_s(\varepsilon_{i})$, $\tau_c(\varepsilon_{i})$ and $\rho(\varepsilon_{i})$ to denote, respectively, the process trace, the activity, the enablement time, the start time, the completion time, and the resource associated with the activity instance $\varepsilon_{i}$.
We use $(\tau_{i}, \tau_{j})$ to denote the interval of time starting in $\tau_{i}$ and ending in $\tau_{j}$.
In this way, $\omega(\varepsilon_{i}) = (\tau_e(\varepsilon_{i}), \tau_s(\varepsilon_{i}))$ denotes the \textit{waiting time} of $\varepsilon_{i}$, representing the interval since $\varepsilon_{i}$ became available for processing ($\tau_e(\varepsilon_{i})$), until its recorded start ($\tau_s(\varepsilon_{i})$).
Accordingly, the processing time of $\varepsilon_{i}$ is $pt(\varepsilon_{i}) = (\tau_s(\varepsilon_{i}), \tau_c(\varepsilon_{i}))$.
We use $(\tau_{i},\tau_{j}) \in (\tau_{k},\tau_{l})$ to denote that the interval $(\tau_{i},\tau_{j})$ is contained in $(\tau_{k},\tau_{l})$, i.e.\ $\tau_{i} \geq \tau_{k}$ and $\tau_{j} \leq \tau_{l}$.
With $(\tau_{i},\tau_{j}) \perp (\tau_{k},\tau_{l})$ we denote that both intervals (partially or fully) overlap, i.e., $\exists (\tau_{m},\tau_{n}) \in (\tau_{i},\tau_{j}) \mid (\tau_{m},\tau_{n}) \in (\tau_{k},\tau_{l})$.

Given the above, an \textit{activity instance log} $L$ is a collection of activity instances recording the data of the execution of a set of traces of a business process.

\subsection{Related Work}
Process mining is a research field that encompasses a wide range of techniques that enable process discovery and efficiency, quality, compliance, predictive, and prescriptive analysis \cite{milani2022process}.
Some techniques address the question of how to discover and analyze waiting times in specific application domains. For instance, Uysal et al.~\cite{uysal2020process} present a case study where process mining is used to identify bottlenecks and reduce the cycle time in a production process.
Similarly, Erdogan et al.~\cite{erdogan2022multi} apply process mining techniques to identify waiting times in a hospital emergency process, while Yampaka \& Chongstitvatana~\cite{yampaka2016application} describe an application of process mining combined with a queuing system to analyze and improve temporal performance in a healthcare process. Similarly, Antunes et al.~\cite{antunes2019solution} combine process mining with discrete event simulation to optimize waiting time in an emergency department.
However, none of these domain-specific studies considers the question of how to attribute waiting times to their causes, which is the focus of this paper.

Ferreira \& Vasilyev~\cite{ferreira2015using} present a technique to identify why some cases in a process take longer time to complete.
They identify case characteristics correlating with higher delays, e.g., when a given activity occurs in a case, or when a given resource is involved, the case is likely to have higher waiting time. Likewise,
De Leoni \& van der Aalst~\cite{de2016general} combine some of the existing correlation analysis techniques to identify how different process characteristics correlate with the process performance, e.g., if process deviations cause delays.
Similarly, Hompes et al.~\cite{hompes2017discovering} propose an approach based on time series analysis to detect cause\mbox{-}effect relations between process characteristics and performance indicators, e.g., if the waiting time for the receipt of a payment depends on the time of day.
Toosinezhad et al.~\cite{toosinezhad2020detecting} introduce an approach to detect event patterns that frequently precede, i.e. lead to, dynamic bottlenecks.
While these studies take a correlation-based approach to analyze waiting times, we classify the causes of waiting times and consider their impact on process performance.

Some process mining techniques support the identification of waiting times caused by queuing effects, i.e., when activity instances wait in a queue until a resource becomes available~\cite{senderovich2015queue}.
Similarly, in~\cite{lashkevich2022data}, the authors present an approach to discover waiting times caused by batch processing.
However, these techniques focus on identifying a singular cause of waiting time. In contrast, in the present paper we seek to decompose the observed waiting time into multiple causes.

\section{Waiting Time Discovery and Analysis\label{sec:discovery-analysis}}

In this section, we describe our approach to discover and analyze waiting times in a business process.
The approach takes an event log as input and, as a result, produces a report comprising the causes of waiting times in a business process and their impact on the CTE.
\figurename~\ref{fig:mesh1} depicts an overview of the three main steps of the proposed approach.

\begin{figure*}[b]
    \centering
    \includegraphics[width=0.99\linewidth]{images/approach.png}
    \caption{Overview of the proposed approach.}\par
    \label{fig:mesh1}
\end{figure*}

In the first step, we discover the transitions between activities and their characteristics -- total frequency, case frequency, and total waiting time -- from the event log.
In the second step, we identify the causes of waiting time of each transition.
In the third and final step, we analyze the impact of each cause on the temporal efficiency of the process.

\subsection{Activity Transition Discovery}

The first step of our approach is to discover transitions between activities and their waiting times.
To do this, we define an \textit{activity transition instance} as a pair of activity instances $\langle a_{1}, b_{1} \rangle$ in a single case, such that the completion of $a_{1}$ enables $b_{1}$, i.e., $b_{1}$ cannot be executed before $a_{1}$ is completed.
In a sequential process, each activity instance of a case is enabled by the completion of the preceding activity instance.
However, concurrency is common in real-life processes.
\figurename~\ref{fig:mesh2} shows an example of the execution of a sub-process formed by four activities with concurrency between two of them.
In this case, even though the order of the activity instances is $a_1$, $b_1$, $c_1$, and $d_1$, both $b_1$ and $c_1$ are enabled when the activity instance $a_1$ completes.
In the same way, the activity instance $d_1$ is enabled only when $c_1$ is completed.
Thus, there are only three activity transitions in this example, namely $\langle a_1, b_1 \rangle$, $\langle a_1, c_1 \rangle$, and $\langle c_1, d_1 \rangle$.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.95\linewidth]{images/flow-parallel.drawio.png}
    \caption{Waiting time in a case fragment with concurrent activity instances.}\par
    \label{fig:mesh2}
\end{figure}

To detect the concurrency relations, e.g., when the activities $b$ and $c$ are concurrent (see \figurename~\ref{fig:mesh2}), we use the concurrency oracle from the Heuristics Miner~\cite{DBLP:conf/cidm/WeijtersR11}.
This method computes, for each observed relation between two activities, a degree of confidence whether it is a concurrent or directly-follows relation, based on the percentage of occurrences in each order -- $b$ followed by $c$, or $c$ followed by $b$.
Then, based on a set of defined thresholds, it retrieves the concurrent relations of the process.
In this way, we consider each activity instance being enabled by its closest non-concurrent preceding activity instance.
We assume that once an activity instance is enabled, it remains enabled until its processing starts.

We call the first element of an activity transition instance the \emph{source activity instance}, while the latter is the \emph{target activity instance}.
An \emph{activity transition} is a set of activity transition instances with the same source and target activities, where the \emph{source activity} is the activity executed in all its source activity instances, and the \emph{target activity} is the activity executed in all its target activity instances.
For example, the activity transition $\langle a, b \rangle$ is composed of the set of activity transition instances $\{\langle a_1, b_1 \rangle, \langle a_2, b_2 \rangle, ..., \langle a_n, b_n \rangle\}$ in the event log.
For the sake of simplicity, we will refer to activity transitions as \emph{transitions}, and to activity transition instances as \emph{transition instances}.

Once the transition instances are discovered, we calculate their \emph{duration}, i.e., the waiting times they induce.
The waiting time of a target activity instance in a transition instance is the interval between its enablement and its start time.
For example, in \figurename~\ref{fig:mesh2}, the waiting time in the transition $\langle c_1, d_1 \rangle$ corresponds to the interval $(t_{enabled\_d1}, t_{start\_d1})$.
In this way, we identify the waiting time of each target activity instance per transition.

Finally, for each identified transition (composed of all its transition instances), we compute the following characteristics.
\emph{Case frequency} illustrates the proportion of process cases from total number of cases where this transition is observed.
\emph{Total frequency} indicates the number of occurrences of this transition in the process.
\emph{Total duration} is the sum of the waiting times of all transition instances.
The output of this step is a report depicting all identified transitions and their characteristics, sorted by total duration in descending order.
Based on this information, the analysts can see what transitions cause the highest waiting times and how frequently they are executed.

\subsection{Waiting Time Cause Discovery\label{subsec:wt-type-discovery}}

Once the activity transitions and their characteristics are discovered, we analyze the waiting time of each transition instance and identify their causes.
We target five causes of waiting time (WT): WT due to batching, resource contention, prioritization, resource unavailability, and extraneous factors.
In this section, we define these five causes of waiting times (RQ1) and describe how they can be identified from an event log (RQ2).

Within a single transition instance, waiting time can be explained by one or several causes --e.g.\ the resource was working in another activity half of the waiting time, and out of their working hours the other half.
If there are several waiting time causes in a transition, we propose to classify them in non-overlapping time intervals by identifying them in strict order using the dichotomies' classification illustrated in \figurename~\ref{fig:mesh3}.
Thus, when an interval of the waiting time is classified as a certain cause, it cannot be classified as any other.

We first identify if any intervals of the transition duration are caused by batching.
Then we seek intervals of waiting time caused by resource contention and prioritization, followed by resource unavailability and extraneous factors (see \figurename~\ref{fig:mesh3}).
The order of identification is based on the dominance of a certain cause over the others, i.e., although several phenomena might be observed for the same interval, only the cause with the highest dominance is considered.
For instance, if during the same waiting time interval, both batching and resource unavailability exist, we consider the primary cause for the instance waiting, i.e., batching in this case.
Although the resource might be unavailable, the processing cannot start until the batch is accumulated.
Therefore, when multiple causes are observed, only one is considered as the direct cause of waiting time.
In addition, this approach ensures the addition property of the waiting time for each transition instance -- i.e.\ the sum of the waiting time causes is equal to the total waiting time of the transition instance.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{images/flow-extended-1.png}
    \caption{Overview of the waiting time cause discovery process and their definitions.}\par
    \label{fig:mesh3}
\end{figure}

\vspace{10pt}

\noindent\textbf{Waiting time due to Batching.} The first cause that we identify is \emph{waiting time due to batching}.
Batch processing occurs when a set of instances of the same activity are accumulated to be processed together (either simultaneously or one after the other)~\cite{lashkevich2022data}.
In this context, a batch is a set of activity instances $\mathcal{E}_b \subseteq L$, such that all $\varepsilon_{i} \in \mathcal{E}_b$ record the execution of the same activity, performed by the same resource, and processed as a batch (i.e.\ all of them were enabled before any of them started, and they were processed as a group).
We use the technique proposed in~\cite{lashkevich2022data} to identify batch processing.
In batch processing, when an activity instance is enabled and ready to be processed, it can wait for other instances of the same activity until the batch is accumulated, i.e., until all instances that are a part of a batch are collected.

Accordingly, when the target activity instance of a transition instance is detected as part of a batch, the waiting time interval from its enablement time to the batch accumulation time is classified as \emph{waiting time due to batching} (see \figurename~\ref{fig:mesh4}).
The waiting time due to batching of an activity instance is then defined as follows:

\begin{defn}[Waiting Time Due to Batching\label{def:batching-wt}]

    Given an activity instance log $L$, a batch $\mathcal{E}_b \subseteq L$, and an activity instance $\varepsilon_{i} \in \mathcal{E}_b$. The WT due to batching of $\varepsilon_{i}$ is $\omega_{ba}(\varepsilon_{i}) = (\tau_{e}(\varepsilon_{i}), \tau_{bc}) \mid \tau_{bc} = max(\{\tau_{e}(\varepsilon_{j}) \mid \varepsilon_{j} \in \mathcal{E}_b\})$, i.e.\ the interval of time between the enablement of $\varepsilon_{i}$ and the last enablement of the activities in the batch.

\end{defn}

Following with the example in \figurename~\ref{fig:mesh4}, the activity $B$ is a batch processed activity where the resource accumulates instances and then processes them one by one (sequential batch processing)~\cite{Martin2021}.
The batch is accumulated until the activity instance $b_3$ is enabled (case C3 is ready to be processed ($t_{enabled\_b3}$).
Therefore, if we analyze the transition instance $\langle a_1, b_1 \rangle$ of case \texttt{C1}, its WT due to batching corresponds to the interval of the waiting time between $t_{enabled\_b1}$ and $t_{enabled\_b3}$.

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.5]{images/batching.png}
    \caption{Waiting time due to batching.}\par
    \label{fig:mesh4}
\end{figure}

\vspace{10pt}

\noindent\textbf{Waiting time due to resource contention.}
There are situations where the resources that have to process a certain activity are busy processing another activity instances that were enabled earlier than the waiting one and, thus, it's understood that they start processing them before the current one (following a first-in-first-out order).\footnote{
We assume that resources work only on one activity at a time, i.e., there is no multitasking.
Thus, we foresee that the proposed estimation technique will not be suitable for event logs with a high proportion of multitask activity instances.
}
When this situation occurs, we classify as \emph{waiting time due to resource contention} those intervals in which the resource that performed the activity instance was working in other activity instances enabled before it (see \figurename~\ref{fig:mesh5}).
Therefore, the waiting time due to resource contention of an activity instance is defined as follows:

\begin{defn}[Waiting Time Due to Resource Contention\label{def:contention-wt}]

    Given an activity instance log $L$, and an activity instance $\varepsilon_{i} \in L$. The WT due to resource contention of $\varepsilon_{i}$ is $\Omega_{rc}(\varepsilon_{i}) =
    \{
    (\tau_{i}, \tau_{j}) \mid \
      \tau_{i} = max(\tau_e(\varepsilon_{i}), \tau_s(\varepsilon_{j})) \wedge
      \tau_{j} = min(\tau_s(\varepsilon_{i}), \tau_c(\varepsilon_{j})) \wedge
      \varepsilon_{j} \in L \wedge
      \varepsilon_{j} \neq \varepsilon_{i} \wedge
      \rho(\varepsilon_{j}) = \rho(\varepsilon_{i}) \wedge
      \tau_e(\varepsilon_{j}) \leq \tau_e(\varepsilon_{i}) \wedge
      pt(\varepsilon_{j}) \perp \omega(\varepsilon_{i})
    \}
    $,
    i.e., the set of intervals of processing time of all $\varepsilon_{j}$ of $L$ (executed by the same resource as $\varepsilon_{i}$, and enabled before it) overlapping with the WT of $\varepsilon_{i}$.

\end{defn}

\vspace{5pt}

\noindent\textbf{Waiting time due to prioritization.}
However, the resources might not always follow the FIFO policy.
In some situations, the resources might give priority to certain activity instances over others.
We call this behavior \emph{prioritization}, meaning that an activity instance is processed out of turn w.r.t. a FIFO policy, thus causing other activity instances to wait longer.
When this situation occurs, we classify as \emph{waiting time due to prioritization} those intervals in which the resource that performed the activity instance was working in other activity instances enabled after it (see \figurename~\ref{fig:mesh5}).
Therefore, the waiting time due to prioritization of an activity instance is defined as follows:

\begin{defn}[Waiting Time Due to Prioritization\label{def:prioritization-wt}]

    Given an activity instance log $L$, and an activity instance $\varepsilon_{i} \in L$. The WT due to prioritization of $\varepsilon_{i}$ is $\Omega_{prior}(\varepsilon_{i}) =
    \{
    (\tau_{i}, \tau_{j}) \mid \
      \tau_{i} = \tau_s(\varepsilon_{j}) \wedge
      \tau_{j} = min(\tau_s(\varepsilon_{i}), \tau_c(\varepsilon_{j})) \wedge
      \varepsilon_{j} \in L\ \wedge\
      \varepsilon_{j} \neq \varepsilon_{i} \wedge
      \rho(\varepsilon_{j}) = \rho(\varepsilon_{i})\ \wedge\
      \tau_e(\varepsilon_{j}) > \tau_e(\varepsilon_{i})\ \wedge\
      pt(\varepsilon_{j}) \perp \omega(\varepsilon_{i})
    \}
    $,
    i.e., the set of intervals of processing time of all $\varepsilon_{j}$ of $L$ (executed by the same resource as $\varepsilon_{i}$, and enabled after it) overlapping with the WT of $\varepsilon_{i}$.

\end{defn}

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{images/cont+prioritization.png}
    \caption{Waiting time due to resource contention and due to prioritization.}\par
    \label{fig:mesh5}
\end{figure}

\vspace{5pt}

\noindent\textbf{Waiting time due to resource unavailability.}
The fourth cause of waiting time that we propose to consider is resource unavailability, which corresponds to the intervals in time in which the resource is not available to work due to their working schedules.
To identify this waiting time, we need to first discover the working schedules of the resources.
We propose to use the technique presented in~\cite{DBLP:conf/bpm/Lopez-PintadoD22} to discover calendars over time granules not fully described by the input data.
This technique analyzes the instants in time when each resource interacted with the system (i.e.\ the start and end of each activity instance) to build a weekly working calendar composed of time intervals in which there was enough evidence, based on given support and confidence values, that the resource is working.\footnote{Although we use this resource calendar discovery algorithm, the approach can be applied with other calendar discovery algorithms or with manually defined calendars.}
Given the weekly calendars of each resource, we transform them to absolute time intervals to compare them with the waiting times observed in the log.
Then, we classify as \emph{waiting time due to resource unavailability} those intervals where the resource is not available for working.
Therefore, the waiting time due to resource unavailability of an activity instance is defined as follows:

\begin{defn}[Waiting Time Due to Resource Unavailability\label{def:unavailability-wt}]

    Given an activity instance log $L$, an activity instance $\varepsilon_{i} \in L$, and being $cal_{av}(\rho) = \{(\tau_{avs}, \tau_{avc})\}$ a resource availability calendar with the set of time intervals in which the resource is available to work.
    The WT due to resource unavailability of $\varepsilon_{i}$ is
    $\Omega_{unav}(\varepsilon_{i}) =
    \{
    (\tau_{i}, \tau_{j}) \mid \
      (\tau_{i}, \tau_{j}) \in \omega(\varepsilon_{i}) \wedge
      \nexists (\tau_{k}, \tau_{l}) \in cal_{av}(\rho(\varepsilon_{i})) \mid (\tau_{i},\tau_{j})
      \perp (\tau_{k},\tau_{l})
    \}
    $,
    i.e., the set of intervals of the waiting time of $\varepsilon_{i}$ that do not overlap with the availability calendar of the resource $\rho(\varepsilon_{i})$ that executed $\varepsilon_{i}$.

\end{defn}

\begin{figure}[t]
    \centering
    \includegraphics[scale=0.5]{images/unavailability.png}
    \caption{Waiting time due to resource unavailability.}\par
    \label{fig:mesh7}
\end{figure}

\vspace{5pt}

\noindent\textbf{Waiting time due to extraneous factors.}
The last cause of waiting time that we propose to consider is ``extraneous factors''.
We propose to classify as \emph{waiting time due to extraneous factors} the waiting time intervals caused by the external effects that cannot be identified from the event log -- e.g.\ the resource is working on another process, the activity instance cannot start until some unrecorded event has occurred, fatigue effects, or context switch.
Thus, the waiting time of an activity instance due to extraneous factors is defined as follows.

\begin{defn}[Waiting Time Due to Extraneous Factors\label{def:extraneous-wt}]

    Given an activity instance log $L$ and an activity instance $\varepsilon_{i} \in L$.
    The WT due to extraneous factors of $\varepsilon_{i}$ is
    $\Omega_{extr}(\varepsilon_{i}) =
    \{
        (\tau_{i},\tau_{j}) \mid \
        (\tau_{i},\tau_{j}) \in \omega(\varepsilon_{i}) \wedge
        \nexists (\tau_{k},\tau_{l}) \in (
            \{\omega_{ba}(\varepsilon_{i})\} \bigcup
            \Omega_{rc}(\varepsilon_{i}) \bigcup
            \Omega_{prior}(\varepsilon_{i}) \bigcup
            \Omega_{unav}(\varepsilon_{i})
        ) \mid (\tau_{i},\tau_{j}) \perp (\tau_{k},\tau_{l})
    \}$,
    i.e., the set of intervals within the WT of $\varepsilon_{i}$ that does not overlap with waiting times due to batching, resource contention, prioritization, or resource unavailability.

\end{defn}

Therefore, the second step of the approach results in a report discovering and measuring the total waiting time for each transition according to their causes.
These results reveal the cause of each waiting times and by how much.

\subsection{Waiting Time Analysis}

The final step is to analyze how much each waiting time cause contributes to the temporal performance of the process and their impact to the CTE.
With that purpose, we propose to compute the percentage of time that each cause induces in the CTE of the process.
In this context, we measure the CTE as the processing time (PT) divided by the sum of the processing and the waiting time (PT $+$ WT), where PT is the sum of the processing time of all activity instances in the process, and WT is the sum of their waiting time.
The impact of waiting times per cause is calculated as the difference between the original process CTE and the CTE if a particular waiting time is eliminated.
In this way, we can measure (1) the impact that each waiting time cause has on the process CTE, (2) the impact that each transition has on the process CTE, and (3) the impact that each waiting time cause has in each transition.
This metrics can indicate the potential CTE improvement if a particular cause of waiting time is addressed.

As a result of the proposed analysis, we can analyze the discovered waiting time causes and their impact on CTE.
With this information, process analysts can identify where the inefficiencies due to waiting times are localized, choose which transitions and/or waiting time causes to address, and which redesign alternatives to apply.

\section{Evaluation\label{sec:evaluation}}

In this section, we present the evaluation of our approach.
We evaluate the approach by addressing two evaluation questions (EQ): (1) \textit{To which extent is the technique able to detect the presence or absence of certain waiting time causes?}, (2) \textit{To what extent is the technique able to correctly quantify the amount of waiting time waste per each cause?}
In this experimentation we use synthetic data to validate the ability of the technique to accurately discover transitions, their waiting times, and the waiting time causes known to be present in the event logs.
Then, we demonstrate the approach's applicability on a real-life event log.
The approach implementation, the event logs and experiment results, are all available on GitHub.\footnote{\url{https://github.com/AutomatedProcessImprovement/waiting-time-analysis/tree/caise2023}}

\subsection{Evaluation on Synthetic Data}

To answer EQ1, we used a business process simulation model (BPS model) of a loan application process to simulate a set of event logs with different combinations of waiting time causes.
To simulate WT due to resource contention, we set a low number of available resources in the BPS model, so that in some cases there were no resources available to process an enabled activity instance.
To create WT due to resource unavailability, we set resource working calendars so that some resources worked from Monday to Wednesday, and others from Thursday to Friday.
To simulate WT due to extraneous factors, we added timer events before some of the activities in the BPS model, thus delaying their start.
WT due to batching and prioritization cannot be injected by modifying the simulation parameters, as current BPS engines do not support them.
Therefore, in case of batching, we added a set of new traces delaying the start of some activity instances so that they are processed as a batch.
To simulate prioritization, we added a set of new traces changing the order of execution of some activity instances, so they are processed following a prioritization order.
Combining these modifications, we simulated a set of 32 event logs with all the combinations of causes of waiting time
and measured the performance using precision and recall.
True positives and false positives stand for the discovery of a waiting time cause that, respectively, was and was not injected in the event log.
True and false negatives denote an undiscovered waiting time cause that, respectively, was not and was injected.

\begin{table*}[t]
    \centering
    \scriptsize
    \caption{Results for the simulated event logs with all waiting time causes, depicting the true positives with '\cmark`, the false positives with '\xmark`, and the true negatives with an empty cell (there are no false negatives).

    }
    \label{tab:simulated-results}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{l c c c c c c c c c c c c c c c c}
                            & S01    & S02    & S03    & S04    & S05    & S06    & S07    & S08    & S09    & S10    & S11    & S12    & S13    & S14    & S15    & S16    \\ \toprule
        Batching            &        & \cmark &        &        &        &        & \cmark & \cmark & \cmark & \cmark &        &        &        & \xmark &        &        \\
        Prioritization      &        &        & \cmark &        &        &        & \cmark &        &        &        & \cmark & \cmark & \cmark & \xmark & \xmark &        \\
        Res. Contention     &        &        &        & \cmark &        &        &        & \cmark &        &        & \cmark &        &        & \cmark & \cmark &        \\
        Res. Unavailability &        &        &        &        & \cmark & \xmark &        &        & \cmark & \xmark &        & \cmark & \xmark & \cmark & \xmark & \cmark \\
        Extraneous factors          &        &        &        &        &        & \cmark &        &        &        & \cmark &        &        & \cmark &        & \cmark & \cmark \\ \midrule
                            & S17    & S18    & S19    & S20    & S21    & S22    & S23    & S24    & S25    & S26    & S27    & S28    & S29    & S30    & S31    & S32    \\ \toprule
        Batching            & \cmark & \cmark & \cmark & \cmark & \cmark & \cmark & \xmark &        &        & \xmark & \cmark & \cmark & \cmark & \cmark & \xmark & \cmark \\
        Prioritization      & \cmark & \cmark & \cmark &        & \xmark &        & \cmark & \cmark & \cmark & \xmark & \cmark & \cmark & \cmark & \xmark & \cmark & \cmark \\
        Res. Contention     & \cmark &        &        & \cmark & \cmark &        & \cmark & \cmark &        & \cmark & \cmark & \cmark &        & \cmark & \cmark & \cmark \\
        Res. Unavailability &        & \cmark & \xmark & \cmark & \xmark & \cmark & \cmark & \xmark & \cmark & \cmark & \cmark & \xmark & \cmark & \cmark & \cmark & \cmark \\
        Extraneous factors          &        &        & \cmark &        & \cmark & \cmark &        & \cmark & \cmark & \cmark &        & \cmark & \cmark & \cmark & \cmark & \cmark \\ \bottomrule
    \end{tabular}
    }
\end{table*}

\tablename~\ref{tab:simulated-results} depicts the results for the simulated event logs used to evaluate EQ1.
Our approach discovered the injected waiting time causes in all the event logs (no false negatives) resulting in a recall of 100\%.
However, due to the influence that some waiting time causes have between them, the results contain 17 false positives (precision of 83\%).
False positives in WT due to resource unavailability (S06, S10, S13, S15, S19, S21, S24, and S28) are caused by the presence of extraneous waiting time, combined with limited data for the discovery of resources' working calendars.
To simulate these logs, we set 24/7 working calendar for all resources (high availability).
However, when a resource has low occupation (executes few events), there is not enough data for the calendar discovery to identify a 24/7 calendar and some intervals are interpreted as non-working time.
When these non-working intervals overlap with extraneous waiting time, the tool classifies them as WT due to resource unavailability.
This limitation is inherent to the discovery of the resource calendar, if there is no data showing that a resource was active during a period of time, it cannot be assumed that they were working.

The injection of WT due to extraneous factors also induced false positives of prioritization (S15, S21, S26, S30).
When an activity instance is enabled but waiting due to extraneous factors, the resource might execute other activities enabled after the waiting one, being detected as WT due to prioritization.
These false positives are due to the absence of an explicit indicator of the WT due to extraneous factors (i.e.\ extraneous waiting time is only detected when no other causes are identified).

False positives due to batching (S14, S23, S26, S31) are caused by the appearance of a batch processing behavior that was not intentionally added.
To create WT due to resource contention and unavailability, in some scenarios we assigned working calendars that produced low resource availability.
In such cases, while instances were waiting until the resource became available, they were collected and then, executed by the same resource.
This resulted in an unassigned but correctly discovered batch processing behavior.

EQ2 aimed at assessing if our approach is able to correctly identify waiting time intervals and their causes.
We could not use the set of event logs used for EQ1 to answer EQ2, as the logs were created through stochastic business process simulation -- we know we have introduced a certain cause of waiting time, but we don't know to what extent.
Therefore,
we manually created a set of 5 event logs (with up to 15 traces) with activity transitions having different waiting time causes.
Due to the low number of events per resource, we manually defined the resource working calendars for this experiment.
We ran our technique over these logs and obtained accurate waiting time intervals and their causes.

\subsection{Evaluation on a Real-Life Log}

To evaluate the applicability of the proposed approach in a real-life scenario, we used an event log of a manufacturing production process~\cite{dafnalevy_2014}.
The event log has 225 traces (cases),
recording the execution of 24 activities in a total of 4,503 activity instances,
executed by 46 resources.

First, we discovered activity transition instances and their characteristics: 91 transitions with 3421 transition instances (i.e., executions of each transition).
For each transition, we identified its characteristics (case frequency, total frequency and total duration, i.e., total waiting time).

Second, we discovered the causes of waiting times.
As a result, we obtained a report that captures how much of the waiting time is induced by each cause in every activity transition.
The highest waiting times originated from the self-loop transitions, i.e., when the same activity is executed for the case twice in a row.
For instance, transitions between two activity instances of ``Turning \& Milling Q.C.'' induced a total waiting time of 1101d 5h 34m. This made it the greatest contributor to the total waiting time in the process.
The largest portion of the waiting time in this transition (54.74\%) was caused by resource unavailability and equaled 602d 18h 59m.

Finally, we analyzed the waiting times by calculating how much each cause contributed to the total WT of the process and how they affected the CTE.
Thus, resource unavailability was the major source of WT as it caused 57\% of the total WT.
It is followed by batching (22\%), prioritization (9\%), extraneous (8\%) and resource contention (4\%).
Then, we measured the impact of the waiting times by cause on the process CTE (CTE = 6.81\%).
The highest CTE increase (up to 14.62\%) can be achieved if the WT due to resource unavailability is eliminated.
From the transition perspective, the highest improvement opportunity lies in addressing the ``Turning \& Milling Q.C'' self-loop transition.
If the WT is this transition is addressed, the CTE could increase to 7.69\%.

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

This paper outlines a process mining-based approach for identifying causes of waiting times and their impact.
To address research question RQ1, we propose a method to attribute the waiting time of activity transitions to one of five causes: \textit{batching}, \textit{resource contention}, \textit{prioritization}, \textit{resource unavailability}, and \textit{extraneous factors}.
To address RQ2, we outline how these five causes of waiting time can be discovered, for each activity transition, from an event log.
Finally, we propose to measure the impact of each waiting time cause on the CTE of a process, to help analysts to prioritize opportunities to reduce these waiting times (RQ3).
The empirical evaluation shows that our approach can accurately classify the waiting times in a process into the five causes, in (synthetic) event logs where the causes of waiting time are known.
Finally, we illustrated the applicability of our approach using an event log of a production process.

In its current form, the proposed approach only considers waiting times associated with transitions between activity instances. Yet waiting times may also arise in at least two other settings: (i) between case creation and start of the first activity instance; and (ii) within an activity instance due to interruptions (e.g. the resource interrupts their work and resumes it later). The first of these  waiting times could be analyzed by applying methods that estimate the inter-arrival time of each trace~\cite{DBLP:conf/bpm/MartinDC15,DBLP:conf/icpm/BerkenstadtGSSW20}. The second approach requires new methods for modeling and inferring interruptions, possibly using additional attributes in an event log or other additional data.

Another limitation of the approach is that it does not consider multitasking. This could be addressed by inferring multitasking patterns from the log, and using this information to estimate at what point in time a resource would normally have started an activity instance, given their past multitasking behavior.

In future work, we plan to develop a method for discovering business process simulation models from event logs, which takes into account the causes of waiting times considered in this paper. Such simulation models could be used to support analysts in identifying combinations of redesign options to optimize CTE, while also considering other performance dimensions (e.g.\ cost).

\bibliographystyle{splncs04}


text:
\section{Introduction}
\label{SectionI}

The stability of a dispersion against coagulation has been one of the primary
subjects of research in colloid science. Colloidal stability can be obtained
in polar solvents by means of ionic molecules. In this case an electrical
double layer forms at the interface between the colloidal particles and the
solvent, and the dependence of its structure on the electrolyte properties
has been the key ingredient to explain electrostatic-induced colloidal
stability in the past \cite{Shaw}.
In practice however electrostatic stability is difficult to obtain, due to
extreme sensitivity of the process to the electrolyte
conditions. A viable alternative to electrostatic stabilization is to
absorb or end-graft non-ionic polymers on the surface of colloidal
particles  \cite{Vincent,Milner1}.
The stabilization mechanism obtained in this case is steric in nature.
The theory of stabilization of colloids by adsorbed neutral polymers
\cite{Alexander,deGennes} is one of the first problems where self-consistent
field theory has been applied \cite{Dolan,Dolan2}.
Self-consistent field theory, introduced in the context of neutral
polymer systems by Edwards \cite{Edwards} and Helfand \cite{Helfand},
together with its discrete, lattice formulation counterpart \cite{Scheutjens},
has been applied and extended to a variety of inhomogeneous systems,
including block-copolymer systems \cite{Matsen1}, polymers at interfaces,
as well as non-ionic polymer brushes \cite{Orland,Netz,Matsen2}.
The study of polymer chains terminally attached to an interface has been
branching from the colloidal science area to many other domains of research.
Polymer brushes are a topic of interest in surface science,
in view of their ability to modify surface properties \cite{Taunton,Singh,Carignano}.
In particular chain-end-functionalized polymer brushes, with
uncharged \cite{Huang}, as well as charged \cite{Lahann,Gong,Carignano2} end-groups,
have been used to design and control surface properties in a variety of
different ways. Theoretical studies of charged-end-group brushes \cite{Zhulina3},
and their response to an electrical external field \cite{Andelman}, have been
discussed recently.

Polyelectrolytes are macromolecules containing ionizable groups.
Their tendency to dissociate in solution into charged groups and
low-molecular-weight counter-ions makes polyelectrolytes an interesting
problem, as both steric interactions of the high-molecular-weight structure
and electrostatic interactions between counter-ions and charged groups,
coexist at different length scales.
Due to the intricate interplay of these different type of interactions,
the properties of polyelectrolyte solutions, as well as the behavior of
polyelectrolytes near an interface, have been considerably more difficult
to explain than those of non-ionic polymers.
Extensions of the lattice self-consistent field theory for neutral polymers
\cite{Scheutjens} have been discussed \cite{Lyklema}.
The self-consistent theory of inhomogeneous polyelectrolyte systems has
been formulated \cite{Shi,Borukhov}, and recently reconsidered in the
context of symmetric di-block polyelectrolytes and polyelectrolyte
blends \cite{Wang,Wang2}.
As discussed in the past \cite{Pincus,Marcelja,Misra}, for a system
of end-grafted polyelectrolytes, the charge distribution is localized at
the interface and symmetry considerations, as well as considerations about
the double layer that forms at the surface, simplify their study.
We consider both the case of end-grafted polyelectrolyte chains and end-grafted polymer
chains with charged-end-groups at a flat interface, by means of self-consistent field theory.
Related studies have been presented in the past \cite{Muthukumar2,Orland2,Zhulina3,
Witte} and a detailed comparison of the results obtained within self-consistent
field theory with the predictions of analytic strong-stretching theory \cite{Zhulina2}
are discussed in this work. The conformational properties of the system in the presence of a
uniform external field \cite{Andelman} are then addressed.
The response of a polyelectrolyte brush to an external electric field
represents an important problem with several practical applications.
Its role in the actuation of nano-size cantilevers, to cite one of
the many possible applications, has been recently reported \cite{Huck}.
We study a system of end-confined charged polymer chains and its
response to a uniform external field. We discuss the conformational changes of
the monomer-density and counter-ion distribution profiles, for both
the case of a similarly and oppositely charged grafting surface.
We propose a new mathematical framework, that generalizes
strong-stretching theory for charged grafted polymers in the presence of
a uniform field and we discuss the results of our findings.

\section{Theory}
\label{SectionII}

We introduce in this section the theoretical methods to describe the
properties of a polyelectrolyte system of chains attached to a flat interface.
The system is composed of $n_p$ polymers chains of length $N_p$, grafted to a
flat, charged substrate of area $ { \cal A}$ with grafting density
$ \sigma= n_p / { \cal A}$, in the presence of a second surface, at a
distance $d$ and in the presence of $n_c$ counter-ions of valence $z_c$ in solution.
We consider the case of strongly dissociating polyelectrolytes, so that the
number and position of the charges along the polymer chains is fixed and
determined by the ionization degree parameter $f$, the total charge per
chain being $fN_p$. The charge distribution along the polymer will be
referred as $z_p(s)$, $s$ being the parameter that measures distance along
the polymer backbone. Both the uniformly charged as well as the neutral
brush with charged-end-groups are described by the present formalism.
The ionization degree, polymer length and grafting density determine the
counter-ion concentration, according to the electro-neutrality condition,
$n_cz_c+fN_p n_pz_p=0$.
The presence of mono-valent salt may easily be included in the notation
introduced below, as well as the presence of other 'small' and polymeric
charged species, according to the general formulation \cite{Shi} discussed
in Appendix \ref{APP:A}.
Each polymer occupies a volume $N_p/ \rho_0$ and has a natural end-to-end
length of $aN_p^{1/2}$, $a$ being the Kuhn segment. We assume all charged
species, namely counter- and co-ions, to be point-like particles, neglecting
steric contributions due to non-polymeric species.
A general formulation of the theory, where several species are present and
the size of the counter-ions and co-ions is explicitly taken into account,
is discussed in Appendix \ref{APP:A}.
The system is immersed in a solvent. As in the case of non-ionic systems,
the corresponding degrees of freedom can be integrated out, so that the
effective interaction between the monomers is described by the effective,
short-ranged excluded-volume interaction parameter $v$.
The Bjerrum length $l_B=e^2/ \varepsilon ({ \bf r}) k_BT$, together with the
reference density $ \rho_0$ and the excluded-volume parameter $v$,
characterizes the properties of the system.

We consider in this section a single polyelectrolyte species in the presence of counter-ions, so that $N_p=N$
and $n_p=n$.
The dimensionless counter-ion concentration is defined as
\begin{equation}
\hat{ \phi}_c( { \bf r}) = \frac{1}{ \rho_0}  \sum_{\alpha=1}^{n_c}
\delta \big ( { \bf r} - { \bf r}_c^{ \alpha} \big ),
\label{density1}
\end{equation}
and the monomer concentration is defined as
\begin{equation}
\hat{ \phi}_p( { \bf r}) = \frac{N}{ \rho_0}  \sum_{\alpha=1}^{n}
\int_0^1 \delta \big ( { \bf r} - { \bf r}_p^{ \alpha}(s) \big ) ds,
\label{density2}
\end{equation}
where the configuration of the $\alpha$-th chain, belonging to the polymer
species $p$ is defined by the vector ${ \bf r}^{ \alpha}_p$ and where the
position of the $ \alpha$-th counter-ion is given by ${ \bf r}_c^{ \alpha}$.

The partition function associated with a system of polyelectrolyte chains
in the presence of counter-ions, tethered to a grafting surface with
density $ \sigma$, can be converted in a statistical field theory
\cite{Parisi}, according to the general formulation of
model $F$ and model $K$, as defined in the classical inhomogeneous polymer
system classification \cite{Fredrickson}.
After integrating out the solvent degrees of freedom, and considering the
case of strong polyelectrolytes, so that the polymer charge distribution
is independent on the location of each monomer and smeared along the chain
with a uniform charge ratio $f$, the effective hamiltonian of the system reads

\begin{eqnarray}
&&~~~\beta { \cal H}[w_p,\psi] = \frac{1}{2} \int d { \bf r} \Big [ \frac{1}{v}
w_p({ \bf r})^2+\frac{1}{4 \pi l_B }  | \nabla \psi( { \bf r})|^2 \Big ] \nonumber \\
&&- \sigma \int d { \bf r}_{ \perp} \ln { \cal Q}_p({ \bf r}_{\perp}, iw_p + iz_p
\psi )- n_c \ln { \cal Q}_c ( i \psi), \\ \nonumber
\end{eqnarray}
where $w_p({ \bf r})$ and $ \psi({ \bf r})$ are the two fields related to monomer concentration
and electrostatic potential, where $ { \cal Q}_p$ and $ { \cal Q}_c$ are the
single-chain and counter-ion partition functions, defined and in Appendix ~\ref{APP:A}
and where a uniform distribution for the grafting points is assumed.
The excluded-volume parameter, also introduced in Appendix \ref{APP:A}, represents the
effective interaction between monomers.
As both the electrostatic and monomer density fields are uniform in the $x$
and $y$ directions, parallel to the grafting surface \cite{Muller}, the
single-chain polymer energy contribution, along the direction perpendicular
to the grafting surface reads
\begin{equation}
\frac{ E [ z^{ \alpha};0,s]}{k_BT} = \int_0^s \Big [ \frac{3}{2 a^2N}|
\dot{z}^{ \alpha} (t)|^2+w_p(z^{ \alpha}(t)) \Big ]dt.
\label{energy}
\end{equation}
The grafting density $ \sigma$ can be absorbed
introducing rescaled polymer and counter-ion concentrations,
\begin{equation}
\phi_p(z)= \frac{a \rho_0}{ \sigma N^{1/2}} \Big \langle \hat { \phi}_p(z)
\Big \rangle ~~~
\end{equation}
and
\begin{equation}
\phi_c(z) = \frac{a \rho_0}{ f \sigma N^{1/2}} \Big \langle \hat { \phi}_c(z)
\Big \rangle,
\end{equation}
that are now normalized according to
\begin{equation}
\int_0^d \phi_p(z) dz = a N^{1/2}
\label{normalization1}
\end{equation}
\begin{equation}
\int_0^d \phi_c(z) dz= aN^{1/2}.
\label{normalization2}
\end{equation}
The corresponding field equations are
\begin{eqnarray}
\label{field}
w_p(z) &=& \Lambda \phi_p(z)+ N_c \psi(z) \nonumber \\
w_c(z) &=& -N \psi(z), \\ \nonumber
\end{eqnarray}
where the reduced interaction parameter is defined as \cite{Matsen3},
\begin{equation}
\Lambda \equiv \frac{v \sigma N^{3/2}}{a \rho_0},
\end{equation}
and is related to the brush height $L$ according \cite{Kim1} to
$L/aN^{1/2}=(4 \Lambda/ \pi^2)^{1/3}$, where we introduce the number of
ions per polymer parameter $N_c=fN$, and where all length scales are
expressed in units of $a N^{1/2}$.

\section{ Self-Consistent Field Theory}
\label{SectionIII}
The crux of self-consistent field theory is the one-dimensional polymer
partition function,
\begin{eqnarray}
\label{pf}
&& q(z,z_0,s) = \int { \cal D} z_{ \alpha}
\exp \Big ( - \frac{E [ z^{ \alpha};0,s]}{k_BT} \Big ) \nonumber \\
&&~~~ \times~~ \delta(z_{ \alpha}(s)-z) \delta(z^{ \alpha}(0)-z_0), \\ \nonumber
\end{eqnarray}
for a section of $sN$ segments with its ends at $z^{ \alpha}(0)=z_0$
and $z^{ \alpha}(s)=z$.
The polymer partition function satisfies the modified diffusion
equation \cite{Matsen1}
\begin{equation}
\label{spe1}
 \frac{ \partial}{ \partial s} q(z,z_0,s) =
\Big [ \frac{a^2N}{6} \frac{ \partial^2}{ \partial z^2}
- w_p(z) \Big ]q(z,z_0,s)
\end{equation}
where $q(z,z_0,0)=aN^{1/2}\delta(z-z_0)$ and $q(0,z_0,s)=0$,
and where the field equations (\ref{field}) include a contribution from the
dimensionless electrostatic potential $\psi(z)$, that satisfies
a non-linear Poisson-Boltzmann equation
\begin{equation}
\frac{1}{L_B} \frac{ \partial^2}{ \partial z^2} \psi(z)= \phi_c(z)-\phi_p(z),
\label{spe2}
\end{equation}
with boundary conditions
\begin{eqnarray}
\label{bc2}
\Lambda_{GC} \frac{ \partial}{ \partial z} \psi(z) \Big |_{z=0}&=&-1\nonumber \\
\Lambda_{GC} \frac{ \partial}{ \partial z} \psi(z) \Big |_{z=d}&=&+1,\\\nonumber
\end{eqnarray}
that describe an external surface charge, where
\begin{equation}
L_B \equiv [(ze)^2/\varepsilon a N^{1/2}k_bT]/ \sigma N_c
\end{equation}
is the rescaled Bjerrum length parameter and where the
Gouy-Chapman length is defined as
\begin{equation}
\Lambda_{GC} \equiv zf \sigma N / \Sigma L_B.
\end{equation}
The counter-ion distribution is given by
\begin{equation}
\phi_c(z)=\frac{a N^{1/2}}{{ \cal Q}_c} \exp ( - w_c(z)/N),
\end{equation}
where
\begin{equation}
{ \cal Q}_c= \int_0^d \exp( - w_c(z)/N) dz.
\end{equation}
Once the solution $q(z,z_0,s)$ to the modified diffusion equation (\ref{spe1})
has been found, the concentration profile for a chain with its free-end
at $z=z_0$ can be computed according to
\begin{equation}
\phi_p(z;z_0)= \int_0^1 \frac{ q( \epsilon,z,1-s)q(z,z_0,s)}
{q( \epsilon,z_0,1)}ds,
\label{spe3}
\end{equation}
where $\epsilon$ is a small finite distance from the grafting surface
introduced according \cite{Muller} to the standard formulation of model $K$.
The free-end distribution is defined as
\begin{equation}
g(z_0)= \frac { aN^{1/2}}{ { \cal Q}_P}
\exp \big ( - \frac{f_0(z_0)}{k_BT} \big ),
\end{equation}
where
\begin{equation}
f_0(z_0)=- \ln q(0,z_0,1)
\end{equation}
and where finally
\begin{equation}
{ \cal Q}_p= \int_0^d \exp \big ( - \frac{f_0(z_0)}{k_BT} \big ) dz_0,
\end{equation}
and the monomer density distribution satisfies
\begin{equation}
\phi_p(z)= \frac{1}{ a N^{1/2}} \int_0^d g(z_0) \phi_p( z;z_0) dz_0.
\label{concentration}
\end{equation}
The set of self-consistent field equations (\ref{spe1})-(\ref{spe2})
have been solved numerically,
for different values of the four fundamental parameters, namely the
natural brush height $L/aN^{1/2}$, the number of ions per polymer
$N_c\equiv fN$, the rescaled Bjerrum Length $L_B$ and the relative surface
charge parameter
\begin{equation}
f_e=\Sigma/ zf \sigma N,
\end{equation}
where $\Sigma$ is the total charge on the grafting surface in elementary
charge units.
The free-energy of the system can be computed as
\begin{eqnarray}
\label{fe}
&&~~~\frac{F}{k_BTn_c}= \frac{1}{ aN^{1/2} N_c} \int_0^{d} g(z_0)
\big [ \frac{f_e(z_0)}{k_BT}+\ln g(z_0) \big ] dz_0 \nonumber \\
&& + \frac{1}{aN^{1/2}} \int_0^{ d} dz \phi_c(z) \ln \phi_c(z) dz + \frac{ \Lambda}{ 2 aN^{1/2}N_c} \int \phi_p^2(z)dz \nonumber \\
&& ~~~~~~~~~~+ \frac{1}{2 aN^{1/2}}
\int_0^{d} \psi(z) \big [ \phi_p(z)-\phi_c(z) \big ]dz, \\ \nonumber
\end{eqnarray}
where the first term is the polymer free-energy and is obtained as the sum
of the translational entropy of the chain ends and the average of the
stretching energy
\begin{equation}
\frac{f_e(z_0)}{k_BT}=- \ln q( \epsilon,z_0,1)- \frac{1}{ aN^{1/2}}
\int w(z) \phi(z,z_0) dz,
\label{feb}
\end{equation}
the second term represents the translational entropy of the mobile
counter-ions and the third and fourth terms the interaction energy, due to
excluded-volume interactions between polymer segments and the electrostatic
energy related to the charge distribution in the system respectively.
Details on the derivation of equation (\ref{fe}) can be found in
Appendix \ref{APP:B}.

\section{Strong-stretching theory}
\label{SectionIV}
When the polymer chains are stretched it is possible to exploit the
analogy between polymer field theory and quantum mechanics and approximate
the single-chain partition function in equation (\ref{pf}), defined above within the
path-integral formalism, by considering the dominant contribution of the
classical path and the fluctuations around it \cite{Milner2,Johner}.
The problem in this limit relates to the classical mechanics problem of a
particle moving down an incline \cite{Matsen2}.
In the classical limit, the partition function $q(0,z_0,1)$ is dominated by
the path $z_{\alpha}(s)$ that minimizes the energy in equation (\ref{energy}) and an
estimate for the free-energy $f_0(z_0)$ of a chain extending to $z=z_0$, can
be obtained accordingly:
\begin{equation}
q(0,Z_0,1) \propto \exp \Big ( - E[ Z^{ \alpha};1]/k_BT] \Big ),
\end{equation}
where, from here on, we will use dimensionless units $Z=z/aN^{1/2}$.
Following the general principle of strong-stretching theory \cite{Milner2,Zhulina}
we write,
\begin{equation}
\frac{3}{2}S^2(Z,Z_0) \equiv N_c \big [ \frac{3}{2} v^2(Z_0)+U(Z)-U(Z_0) \big ],
\end{equation}
where the potential $U(Z)$ is to be identified with the polymer field $w_p(z)$ introduced in
equation (\ref{field}). The chain-end tension can be written \cite{Johner} as
\begin{equation}
v(Z_0)= - \frac{ d \ln g(Z_0)}{d Z_0},
\label{classical2}
\end{equation}
and, in dimensionless units, the two constraints
\begin{equation}
\int_0^{ Z_0} \frac{1}{S (Z,Z_0)}dZ=1,
\label{isochrone2}
\end{equation}
and
\begin{equation}
\int_Z^{D} \frac{ g(Z_0)}{ S(Z,Z_0)} dZ_0 = \phi_p(Z),
\label{isochrone}
\end{equation}
apply. This yields an expression for the dimensionless speed
at position $Z$, $S(Z,Z_0)$ and the normalization condition in equation (\ref{isochrone2}), corresponding to the isochronicity constraint for the effective particle
rolling down an incline \cite{Matsen3}.
 In the absence of excluded-volume interactions \cite{Zhulina,Zhulina2}, the normalization condition in equation (\ref{isochrone2}) can be written as an
Abel-Volterra integral equation of the first kind \cite{Tricomi}. Its
solution leads to a parabolic form for the potential $U(Z)$ \cite{Milner2,Semenov}.
Neglecting the term related to the end-monomer tension, i.e. setting $v(Z_0) \approx 0$, one finds the integral equation
\begin{equation}
\int_{ U}^{ U_0} dt f(t) (t-U)^{-1/2}=1,
\end{equation}
and its solution, in terms of the Riemann-Liouville fractional derivative
of order $1/2$, as discussed in Appendix \ref{APP:C}, reads
\begin{equation}
f(U)= \sqrt{ \frac{3}{2}} \frac{ dZ}{dU}=
 - \frac{1}{ \pi} \frac{ d }{ d U} \int_{ U}^{U_0} dt ( t - U)^{-1/2}.
\label{sst}
\end{equation}
The associated potential form is obtained as
\begin{equation}
U(Z)= \frac{ 3 \pi^2}{8} ( H^2-Z^2),
\label{sst2}
\end{equation}
and the electrostatic potential $\psi^a(Z)=N_c^{-1/2}U(Z)$ follows a
parabolic form, according to strong-stretching theory \cite{Zhulina}.
We express distances in dimensionless rescaled units
$Z'=Z/N_c^{1/2}$, where $H$ represents the brush height in units of $Z$.
The relative charge ratio per chain reads $Q_1=1-\frac{H}{ \gamma}$, where
\begin{equation}
\gamma= \frac{ 4}{3 \pi^2}L_B N_c^{3/2}.
\end{equation}
The rescaled counter-ion distribution is given by
\begin{equation}
N_c^{1/2} \phi_c^a(Z)=\frac{3 \pi^2 H^2}{8 \gamma} \exp \{ \psi^a(Z) \},
\end{equation}
and the electrostatic potential outside the brush is written as
\begin{equation}
\psi^b(Z)=-2 \ln \Big (  \gamma^{1/2}(Z-H+1/H)\Big ),
\end{equation}
where the brush height \cite{Zhulina} is given by
\begin{equation}
H + \frac{\sqrt{ \pi}}{2} \sqrt{ \frac{3 \pi^2}{8}}H^2
e^{\frac{3 \pi^2}{8}H^2}  erf (\sqrt{\frac{3 \pi^2}{8}}H)= \gamma,
\end{equation}
the monomer density distribution,
\begin{equation}
N_c^{1/2}\phi_p(Z)=N_c^{1/2}\phi_c^a(Z)+1 /\gamma,
\label{phip}
\end{equation}
and the counter-ion distribution outside the brush is given by the
Gouy-Chapman form,
\begin{equation}
 \phi_c^b(Z)=\frac{1}{\gamma}(Z-H+1/H)^{-2}.
\end{equation}
The size of the brush can be expressed as
\begin{equation}
\frac{\langle Z \rangle }{H} =\frac{H^2}{2 \gamma} e^{  \frac{3 \pi^2}{8} H^2}.
\label{thickness}
\end{equation}
Finally the end-monomer distribution can be obtained, as detailed in the
Appendix \ref{APP:C}, solving \cite{Zhulina} a similar Abel-Volterra
equation of the first kind to equation (\ref{isochrone2}), that corresponds to
the constraint of equation (\ref{isochrone}). One finds,

\begin{eqnarray}
\label{gz}
&& \gamma N_c^{1/2}g(Z)= \sqrt{ \frac{3 \pi^2}{8}} Z (1+\frac{3  \pi^2H^2}{8})
\big ( \psi^a(Z) \big )^{-1/2} \nonumber \\
&& ~~~~+ \sqrt{ \pi} \frac{ 3 \pi^2}{8} H^2
\exp( \psi^a(Z)) { \it erf} \big ( \sqrt { \psi^a(Z)} \big ). \\ \nonumber
\end{eqnarray}

In the following section we will compare the strong-stretching theory
predictions with the results obtained within self-consistent field theory,
where the Poisson-Boltzmann and modified diffusion equations
(\ref{spe1})-(\ref{spe2}) are solved explicitly by finite difference methods.
We consider both the case of excluded-volume and electrostatic interactions
being present. It is important to note at this point that typical synthetic
polymers considered in experimental conditions are relatively insensitive
to solvent conditions and the long range electrostatic interactions will
dominate in some cases \cite{Zhulina2}.

\section{Self-Consistent Field Theory-Results}
\label{SectionV}

In the previous section we discussed the main assumptions and results of
strong-stretching theory and introduced the fundamental equations for the
monomer density and electrostatic potential in a system of uniformly charged
polyelectrolyte chains, tethered to a flat interface, in the absence of
the external field.
The numerical method, used to solve the mean-field saddle point equations
discussed within self-consistent field theory, is introduced in this
section and the results throughout discussed and compared with the analytical
predictions of strong-stretching theory \cite{Zhulina}.
In order to solve the coupled set of saddle-point equations
(\ref{spe1})-(\ref{spe2}), we used a real-space, finite difference second-order
method \cite{Henrici}. The modified diffusion equation (\ref{spe1}) was solved
using a discrete Crank-Nicholson real-space method, as discussed in detail in
the past\cite{Matsen3}.
The full Poisson-Boltzmann equation (\ref{spe2}) was solved using an iterative
Newton-Raphson method, where the partition function for the small species is
calculated at each iteration step.
The algorithm we implemented in order to solve the modified
diffusion equation (\ref{spe1}) includes the Anderson mixing method \cite{Thompson}.
We observe our algorithm to converge efficiently in a broad range of
parameter values.
It is important to mention that we observe the range of convergence, namely
the highest number of ions per chain the algorithm can reach for a given
value of the Bjerrum length and brush thickness parameters,
to depend in a non trivial way on the Anderson mixing method. Avoiding
the Anderson mixing method and decreasing the mixing parameter $ \mu$ to
very small values, typically $ \mu =0.01$, improves the range of convergence
but increases the time of convergence significantly.
A typical threshold value we obtained and report for comparison:
$N^{max}_c \simeq 200$, for values of $L=2$ and $L_B = 0.1$ of the
brush thickness and Bjerrum length parameters respectively.

A first set of numerical results has been obtained for values of the brush
thickness parameter $L=2$ and Bjerrum length values $L_B=0.1$ and $L_B=0.003$.
The results for the monomer density and counter-ion distributions,
for increasing ionic strength $N_c=10, \cdots , 100$ and  $N_c=30, \cdots 180$ for two different values of the Bjerrum length parameter $L_B=0.1$ and
$L_B=0.001$ are shown in Fig~\ref{fig1}. Another set of calculations have been
performed at the values of $L \simeq 1.88$, $L_B \simeq 0.68$ to confirm
the agreement with closely related calculations \cite{Orland2}.
We note that, when exceeding the values of the number of ions per
chain $N_c =30$, considered in that paper as a limiting threshold,
a peak at low values of the distance appears, already discussed \cite{Taniguchi}.
Our findings resolve the inconsistency reported \cite{Orland2} for
the monomer density profiles at small values of the grafting distance
$z/aN^{1/2}~ \simeq 0.05$, where a peak, as can be seen by close inspection
of Fig.~\ref{fig1} appears. The presence of this peak has already been
reported \cite{Taniguchi} when considering high values of the ion-per-chain
parameter $N_c>50$.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig1.eps}
\caption{ The two upper plots represent the monomer density profile for
the value of the brush thickness $L/aN^{1/2}=2$ and for increasing values
of the number of ions per chain $N_c=10 \cdots 90$ and
$N_c=30,60, \cdots 180$ , for the values of the Bjerrum Length parameter
$L_B=0.1$ and $ L_B=0.001$ respectively. The two lower plots represent
the corresponding counter-ion concentration profiles for $N_c=20 \cdots 100$,
$N_c=30,60, \cdots, 180$.
\label{fig1}}
\end{figure}

All calculations of this and the following sections were performed
considering value of the rescaled distance $D=d/aN^{1/2}$ between the grafting surface and
the electrode to be large enough, so to measure a vanishing counter-ion
density at the boundaries.
In order to obtain such condition, at the lowest values of the
Bjerrum length parameter considered, e.g. $L_B=0.0003$, we needed to consider
values of order $D=10^3$.
This was possible as the two grids we consider to solve the modified
diffusion equation (\ref{spe1}) and equation (\ref{spe2}) did not have the same size,
the second grid exceeding the size of the first, depending on the
Bjerrum length values considered.
This can be seen looking at the range of distances $z/aN^{1/2}$ we plot in Fig.~\ref{fig1}
for different values of the Bjerrum length parameter $L_B$.

Another method introduced in the past to avoid the numerical self-consistent field
theory approach is to linearize \cite{Marcelja} the Poisson-Boltzmann equation (\ref{spe2})
and to modify the strong-stretching theory accordingly.
As a third, intermediate option, we considered explicitly the Poisson-Boltzmann equation
(\ref{spe2}) and combined its numerical solution with the analytical form of the
potential, as predicted by strong-stretching theory.
We hence assume a parabolic form for the effective potential
$w(Z)$ and solve numerically the Poisson-Boltzmann
equation (\ref{spe2}), and an estimate for the monomer concentration
profile, together with the electrostatic potential, follows directly.
The procedure is iterated until convergence is achieved.
A comparison with the results obtained by the full solution of
equations (\ref{spe1}) and (\ref{spe2}) is shown in Fig.~\ref{fig3}.
Closely related methods and results have also been discussed \cite{Misra}.
\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig2.eps}
\caption{Solid lines represent the self-consistent field theory results for
the monomer concentration profiles for values of the brush thickness parameter
$L/aN^{1/2}=2$ and Bjerrum length $L_B=0.1$, for increasing values of the ion
number parameter $N_c=10, \cdot \cdot , 90$. The dashed lines are the results
obtained with the numerical strong-stretching approach discussed in the text.
  \label{fig2}}
\end{figure}

We now consider the case where excluded-volume interactions are neglected \cite{Zhulina2}. The validity of this assumption \cite{Zhulina} is tested and comparison between polyelectrolyte scaling and strong-stretching theory and the numerical results obtained within self-consistent field theory, when solving the saddle point
equations (\ref{spe1})-(\ref{spe2}) is presented. We assume here the relative strength of the excluded-volume interactions to be small with respect to the average electrostatic interaction strength, so that the solvent can be considered as a $ \theta$ solvent for uncharged chains.
This assumption, already discussed and justified in detail \cite{Zhulina2},
will be considered in the results of Figs.~\ref{fig3}-~\ref{fig6}, so to compare with scaling and strong-stretching theory.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig3.eps}
\caption{ Brush thickness $h/aN^{1/2}$ as a function of the number of ions per chain $N_c$,
for different values of the Bjerrum length $L_B=0.003,0.01,0.03,0.1$,
neglecting excluded-volume interactions.
The two dotted lines are the scaling predictions corresponding to the osmotic
and Pincus regimes \cite{Pincus} respectively.
\label{fig3}}
\end{figure}

In Fig.~\ref{fig3} we present the results for the brush thickness $h/aN^{1/2}$, as
obtained from the first moment of the monomer concentration profile,
for increasing values of the number of ions per chain parameter $N_c$ and for
vanishing values of $L/aN^{1/2}$.
The classical scaling predictions, discussed in the literature \cite{Pincus},
are clearly confirmed, given in particular the large values of the number
of ions per chain parameter our algorithm was capable to achieve.
The data fit the asymptotic behavior expected at large values of the
number of ions per chain parameter $N_c$, as shown by the upper dotted line of Fig.~\ref{fig3},
where the condition of static equilibration of the osmotic pressure dictates \cite{Pincus},
in the strongly charged regime, the scaling law $h/aN^{1/2} = N_c^{1/2}$.
According to the results presented in Fig.~\ref{fig3}, an
intermediate cross-over regime, the so-called Pincus-regime \cite{Pincus}, is
present.
The curves, corresponding to different values of the Bjerrum length parameter
$L_B=0.003,0.01,0.03,0.1$, cross-over to the intermediate scaling regime,
according \cite{Pincus} to the scaling expression $h/aN^{1/2} = N_c^{2}L_B$,
corresponding to the slope of the second dashed line in Fig.~\ref{fig3}.
The region for this crossover to occur is quite narrow: all curves
in Fig.~\ref{fig3} converge to the neutral brush thickness \cite{Netz}
value $h_o/aN^{1/2} \simeq 0.54281$, thus explaining why evidence of the
Pincus regime in numerical simulations  turns out to be an elusive
task \cite{Netz2}.
In Fig.~\ref{fig4} we compare the numerical results obtained by
the numerical self-consistent field theory with the predictions of
strong-stretching theory, as described by equation (\ref{phip}),
for the four values of the strong-stretching parameter
$L_BN_c^{3/2}=3,10,30,100$.
At large values of $ \gamma$ the agreement between the numerical results and
the theory becomes very good, except at small values of the
rescaled distance $Z/N_c^{1/2} < 0.05$, and at the edge of the brush, as expected.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig4.eps}
\caption{ Monomer concentration profile as a function of rescaled distance
from the grafting surface $Z/aN_c^{1/2}$, for different values of the
strong-stretching parameter $L_BN_c^{3/2}=3,10,30,100$. Solid lines are
the results obtained with self-consistent field theory. The dashed
lines are the predictions of strong-stretching theory.
Agreement between strong-stretching and the self-consistent field theory
becomes very good for large values of the parameter $\gamma$.
\label{fig4}}
\end{figure}

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig5.eps}
\caption{ Brush thickness as obtained from the first moment of the monomer
concentration profile predicted by self-consistent field theory (solid-lines),
for four different values of the Bjerrum length $L_B=0.0003,0.001,0.003,0.1$
and the prediction of strong-stretching theory in equation (\ref{thickness}) as a
function of the strong-stretching parameter $\gamma$.
The agreement becomes very good at large values of $\gamma$.
A small discrepancy though is observed. The reason for this discrepancy
is not clear.
\label{fig5}}
\end{figure}

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig6.eps}
\caption{ The left upper plot represents the end-segment distribution
$g(z_0)$. The solid curve corresponds to the results obtained within
self-consistent field theory, at values of the strong-stretching
parameter $L_BN_c^{3/2}=3,10,30,100$. The dashed line corresponds
to the prediction of equation (\ref{gz}) for $L_BN_c^{3/2}=100$.
The agreement is good for all values of the distance from the grafting
surface $z/aN^{1/2}$, except at large values where the prediction of
equation (\ref{gz}) diverges.
In the lower left plot we show the chain-end tension for the four values
of the strong-stretching parameter above.
The upper right plot shows the segment profiles $\phi(z,z_0)$ for
$L_BN_c^{3/2}=100$ for the four values of the end position
$z_0/aN^{1/2}=1,2,3,4$. The corresponding prediction of strong-stretching
theory is represented by the dashed lines.
The corresponding average polymer trajectories are shown in the lower right
plot together with the prediction of strong-stretching theory, shown by the
dashed line.
  \label{fig6}}
\end{figure}

A few conclusive remarks, while closing this section, are here in order.
As shown in Fig.~\ref{fig6}, the fundamental assumption of strong-stretching
theory becomes more accurate for increasing values of the
strong-stretching parameter $L_B N_c^{3/2}$.
At values of $L_BN_c^{3/2}=100$, the chain-end tension, as measured within
self-consistent field theory, approaches zero, for all values of the
distance $z/aN^{1/2}$ inside the brush. Fig.~\ref{fig6} also
shows the segment profiles $\phi(z,z_0)$, for chains with different
end-positions $z_0$ at values of $L_BN_c^{3/2}=100$ and the corresponding
average polymer trajectories $z_{\alpha}(s)/z_o$ as obtained within
self-consistent field theory, together with the prediction of
strong-stretching theory.
We finally present the results obtained within self-consistent field
theory, when considering a uniform charge on the grafting surface. We
obtained results for both the case of similarly and oppositely charged surface.
The monomer-density and counter-ion distribution profiles are shown in
Fig.~\ref{fig7} and Fig.~\ref{fig8}.
\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig7.eps}
\caption{ Monomer concentration $\phi_p(z)$ and counter-ion distribution
profiles $\phi_c(z)$, for brush thickness values of $L/aN^{1/2}=2$ and
for two values of the Bjerrum length parameter $L_B=0.1$ and $L_B=0.001$,
as a function of the distance from the grafting surface, for positive vales
of the relative surface charge parameter $f_e=0.5$, for increasing values
of the number of ions per chain parameter $N_c=10, \cdots , 40$ and
$N_c=30,60, \cdots,150$ respectively.
\label{fig7}}
\end{figure}

The results were obtained for both positive and negative
values of the relative surface charge parameter $f_e$, for a brush thickness
parameter value $L/aN^{1/2}=2$. The boundary condition in equation (\ref{bc2}) has been
changed accordingly. As expected, for the case
of a similarly charged surface, the polymer brush swells, and the monomer-density
profile obtained at small values of the relative surface charge parameter changes its
shape and deviates from the parabolic form. Similarly, when considering an oppositely
charged grafting surface, as in Fig.~\ref{fig8}, the brush contracts and the counter-ion
distribution changes self-consistently. Related results, for the case of
a similarly charged grafting surface, have been reported in the
past \cite{Zhulina2}.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig8.eps}
\caption{Monomer concentration $\phi_p(z)$ and counter-ion distribution
profiles $\phi_c(z)$, for brush thickness values of $L/aN^{1/2}=2$ and
for two values of the Bjerrum length parameter
$L_B=0.1$ and $L_B=0.001$, as a function of the distance from the
grafting surface, for negative vales of the relative surface charge
parameter $f_e=-0.5$ and for increasing values of the number of ions per
chain parameter $N_c=10, \cdots , 60$ and $N_c=30,60, \cdots,240$ respectively.
\label{fig8}}
\end{figure}

\section{Charged-end-group polymer brushes}
\label{SectionVI}

In the previous section we presented numerical self-consistent field theory
results for uniformly charged, end-grafted polyelectrolyte systems,
both in the presence and in the absence of excluded-volume interactions and
discussed the effects of a uniform electric field. We discussed the results for the chain-end,
counter-ion and monomer-density distributions, for different values of the
relative surface charge parameter $f_e$. The analysis suggested that the
parabolic form of the monomer density profile is drastically affected by the
electric field, both in the similarly and the oppositely charged case.

We consider in this section the case of a polymer brush in solution,
characterized by charged end-groups \cite{Zhulina3,Andelman} in the presence of counter-ions
and excluded-volume interactions between the monomers.
Rather than considering a uniformly charged polyelectrolyte, as discussed in section \ref{SectionV}, we study in this section the case of grafted neutral chains, with a functional charged end-group, represented in our model
by the terminal, non-grafted $\chi N$ monomers along the chain. As discussed above, we consider the case
of strong polyelectrolytic functional groups, so that the charges are assumed to stay bounded to the chain-ends.
As for the case of uniformly charged polyelectrolytes, the $N_c$ ions per chain
are uniformly smeared along the terminal, non-grafted section of the chain, with a rescaled charge density $f/ \chi$
so that, consistently with the notation of Appendix \ref{APP:A}, we consider the charge
distribution along the chains to be given by $z_p(s)=f/ \chi \Theta( s- \chi)$.
We solve the saddle-point equations (\ref{spe1})-(\ref{spe2}) that now
involve both the charged and uncharged monomer densities and the
effective field $w_p(z)$ is also defined accordingly.
In Fig.~\ref{fig9} we show the total (charged and uncharged) monomer density distribution of a brush
of thickness $L/aN^{1/2}=4$, for increasing values of the number of ions per chain-end parameter $N_c=1,10,20$.
For small values of the ions per chain parameter, namely  $N_c=1$, the electrostatic interactions
between the chain-ends are negligible, as can be seen in Fig.~\ref{fig9} and the monomer distribution
is close to the parabolic prediction of strong-stretching theory, the brush thickness parameter $L/aN^{1/2}$
being relatively large, as expected.
Deviations from the parabolic profiles can be observed and are due to the inter-chain-end electrostatic
interactions at higher values of the ion-per-chain-end parameter $N_c$.
The calculation has been performed for the specific value of $ \chi=1/30$, so that the size of the
polymer chain is about thirty times the size of the chain-end group. We performed several test runs and
checked that for smaller sizes of the chain-end group the monomer density profiles did not change significantly.
For increasing values of $N_c$ the brush swells, due to the repulsive electrostatic chain-end interactions
and the counter-ion distribution changes accordingly.
We now discuss the response of the system to an electric field. In our analysis we consider small values of the number of ions per chain-end parameter, namely $N_c=1$, so to minimize the interaction between the chain-ends.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig9.eps}
\caption{Monomer density profiles for a brush with thickness parameter $L/aN^{1/2}=4$, Bjerrum length parameter value
$L_B=0.1$ for increasing values of the number of ion-per-chain-end parameter $N_c=1,10,20$. The dashed line
represents the prediction of strong-stretching theory, in the absence of charged-end-groups.
\label{fig9}}
\end{figure}

As discussed in section \ref{SectionIV}, we consider the brush to be grafted to a flat surface that acts
as an electrode, while the second electrode is kept at a distance $D \gg L$. In order to study the response of the
system to the uniform electric field, we solved the saddle point equations (\ref{spe1}) and (\ref{spe2}) changing
the boundary conditions for the electrostatic potential as well as modifying the effective field accordingly. Fig.~\ref{fig11} shows the monomer density
profiles for the brush of thickness $L/aN^{1/2}=4$, Bjerrum length $L_B=0.1$ and number of ions per
chain-end parameter $N_c=1$ for both positive and negative values of the surface charge $f_eL_B= \pm 2$.
The neutral surface charge case $f_e=0$ is also shown. As expected, we observe deviations from the parabolic
profile for finite values of the electric field strength. Differently, for vanishing values of the surface charge
$f_e=0$, the profile is very close to the strong-stretching parabolic prediction, the value of the brush thickness
$L/aN^{1/2}=4$ being relatively large \cite{Matsen3}. In the following section we discuss how to generalize strong-stretching theory for the case of a uniform electric field, assuming the relative charge at the
chain-ends to be small.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig10.eps}
\caption{ Chain-end distribution $g(z_0)$ for a brush with thickness parameter $L/aN^{1/2}=4$, Bjerrum length
parameter value $L_B=0.1$ for increasing values of the number of ion-per-chain-end parameter $N_c=1,10,20$.
The dashed line represents the prediction of strong-stretching theory, in the absence of charged-end-groups.
\label{fig10}}
\end{figure}

\section{Strong-stretching theory of charge-end-group polymer brushes}
\label{SectionVII}

The mathematical foundations of strong-stretching theory stem from the
fact that the isochronicity constraint in equation
(\ref{isochrone2}) may be regarded as an Abel-Volterra integral equation
of the first kind.
This determines uniquely the form of the potential
\cite{Milner2}. To determine how a uniform tension at the chain-ends,
induced by the electric field, affects the potential form $U(Z)$, we
evaluated the corresponding corrections, solving the isochronicity constraint
condition in equation (\ref{isochrone2}). The presence of the electric field
induces a chain-end uniform tension and converts the associated
integral equation from an Abel-Volterra equation of the first kind to an
Abel-Volterra equation of the second kind, as also discussed in detail in
Appendix \ref{APP:C}. This leads to evaluate the main corrections to the
parabolic potential, due to the presence of the field, as we now derive.
Let us consider directly the constraint in equation (\ref{isochrone2}),
that can be written, following the few steps in Appendix \ref{APP:C},
\begin{equation}
\int_{ \eta}^{U_0} dt (t- \eta)^{-1/2}= \int_{ \eta}^{U_0} dt f(t)
 \Big \{ \pi- \alpha (t-\eta)^{-1/2} \Big \},
\label{SST2a}
\end{equation}
where $f(U)= \sqrt{\frac{3}{2}}\frac{dZ}{dU}$ is related to the first
derivative of the inverse function $Z(U)$.
We then find corrections to the strong-stretching results.
Namely, after a few steps, we obtain the following Abel-Volterra integral
equation of the second kind:
\begin{eqnarray}
\frac{df( \eta)}{d \eta} &-& \frac{ \pi}{\alpha^2} f( \eta)
=  \frac{ d}{ d \eta} F( \eta)  \nonumber \\
F( \eta)= g( \eta)&-& \frac{1}{ \alpha} \int_{ \eta}^{U_0} dt g( \eta)
(t- \eta)^{-1/2}, \\ \nonumber
\end{eqnarray}
where $g( \eta)=-1/ \alpha$, and where $ \alpha = \sqrt{6} V_0$ measures
the uniform end-monomer force induced by the electric field at the chain ends.
The solution to this Abel-Volterra integral equation of the second kind
can be obtained as follows:
\begin{equation}
f( \eta)= \frac{ dZ}{ d \eta}= \frac{1}{ \alpha} e^{ \frac{\pi}{ \alpha^2}(U_0- \eta)}
\frac{1}{ \sqrt{\pi}} \Gamma \Big( \frac{1}{2},\frac{\pi}{ \alpha^2}
(U_0- \eta) \Big),
\end{equation}
$\Gamma( \alpha,x)$ being the incomplete gamma function of order $\alpha$.
The above expression can be integrated and the corresponding equation
for the potential reads

\begin{equation}
\sqrt{\frac{3}{2}} Z= \frac{2}{\pi} (U_0-U)^{1/2}+\alpha- \frac{ \alpha}{\sqrt{\pi}} e^{ \frac{\pi}{ \alpha^2}(U_0- \eta)}
\Gamma \Big( \frac{1}{2},\frac{\pi}{ \alpha^2}
(U_0- \eta) \Big).
\label{SST2p}
\end{equation}
The expression above reduces to the parabolic strong-stretching potential
when the chain-end tension vanishes, as expected.
From the potential in equation (\ref{SST2p}) , we determined the monomer-density and
chain-end distributions by means of the field equation (\ref{field}), in the
limit of weakly charged end-groups. This assumption has been justified and
discussed in Fig.~\ref{fig10}.
A similar analysis, detailed in Appendix \ref{APP:C}, leads to an explicit
expression for the chain-end distribution $g(Z_0)$, for small values of the
electric field strength $ \alpha$. We find, after a few steps detailed in
Appendix \ref{APP:C},
\begin{equation}
g( \eta) \frac{ dZ}{ d \eta} = - \frac{2}{ \pi} \eta^{ 1/2} + \frac{ \alpha}{2}- \frac{ \alpha}{ \pi^{3/2}} e^{ \frac{ \pi}{ \alpha^2} \eta} \Gamma \Big( \frac{1}{2},\frac{\pi}{ \alpha^2} \eta \Big).
\label{SST2g}
\end{equation}
Equation (\ref{SST2g}) reduces to the well known expression for the chain-end distribution \cite{Milner2}, when no tension
is present at the chain-ends.
The value of the chain-end prefactor $ \delta \simeq 0.263$ has been obtained comparing the results of numerical self-consistent field theory, as in the inset of Fig. \ref{fig12}, and the results obtained from equation (\ref{SST2p}) for the brush height $h/aN^{1/2}$. It would be interesting to see if the prefactor can be computed explicitly, and if it depends on the brush thickness parameter $L/aN^{1/2}$.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig11.eps}
\caption{Monomer density profile $ \phi_p(z)$ for a brush with thickness parameter $L/aN^{1/2}=4$, Bjerrum length
parameter value $L_B=0.1$, for values of the number of ion-per-chain-end parameter $N_c=1$, for neutral surface
charge $f_e=0$, as well as similarly and oppositely charged grafting surface,with rescaled surface charge parameter values
$f_eL_B= \pm 2.0$. The dotted line represents the prediction of strong-stretching theory, in the absence of charged-end-groups.
The dashed lines represent the prediction of the modified strong-stretching theory (\ref{SST2p}), for values of the parameter $\alpha =\delta f_eL_B$. The inset shows the monomer-density at the reference distance $z_m=0.378$, as measured within
self-consistent field theory, shown by the solid line, and as predicted by equation (\ref{SST2p}), for different values of the chain-end tension parameter $ \alpha$, shown by the star symbols.
\label{fig11}}
\end{figure}

Equations (\ref{SST2p}) and (\ref{SST2g}) generalize the MWC strong-stretching theory to problem of charge-end functionalized brushes in the presence of a
small, uniform external electric field. A comparison with the results
of numerical self-consistent field theory is shown in Fig.~\ref{fig11}
and Fig.~\ref{fig12}.

\begin{figure}
\includegraphics*[angle=0,scale=0.7]{fig12.eps}
\caption{Chain-end distribution $ g(z_0)$ for a brush with thickness parameter $L/aN^{1/2}=4$, Bjerrum length
parameter value $L_B=0.1$ for values of the number of ion-per-chain-end parameter $N_c=1$, for neutral surface
charge $f_e=0$, as well as similarly and oppositely charged grafting surface,with rescaled surface charge parameter values
$f_eL_B= \pm 2.0$. The dotted line represents the prediction of strong-stretching theory, in the absence of charged-end-groups.
The dashed lines represent the prediction of the modified strong-stretching theory in equation (\ref{SST2g}), for different values of the chain-end tension parameter $\alpha= \delta f_e L_B$.
The inset shows the size of the brush $h/aN^{1/2}$, as predicted from numerical self-consistent field theory, shown by the solid line, and as obtained according to the modified strong-stretching theory in equation (\ref{SST2p}), shown by the dotted line. The dashed line shows the results of the modified strong-stretching theory, shifted by the zero-field value, obtained subtracting the brush size as obtained from self-consistent field theory and the MWC theory prediction.
\label{fig12}}
\end{figure}

We conclude observing that the monomer-density profile is affected by the uniform electric field in a different way that what seen for uniformly charged polyelectrolyte brushes. Note however that, when comparing the results of this section with the monomer-density profiles of Fig.~\ref{fig7} and Fig.~\ref{fig8} one should keep in mind that the results discussed for the polyelectrolyte problem were obtained at relatively large values of the number of ions per chain parameter $N_c$, the electrostatic interactions between monomers being relatively strong.
It is interesting to observe that for a uniformly charged polyelectrolyte system,at low values of the number of ions per chain parameter, e.g. $N_c=1$, the effect of the uniform electric field is very similar to what observed for charged-end-group chains and the monomer-depletion effect close to the surface,
seen in the analysis of Fig.~\ref{fig7} for a similarly charged surface, is not present, as in the case of grafted charge-end-group polymer chains discussed in this section. This suggests that the strong-stretching theory corrections discussed in this work might apply reasonably well for the problem of weakly charged
polyelectrolyte chains in a uniform electric field.

\section{Conclusions}
\label{SectionVIII}
In this work we studied the properties of charged polymers grafted to a flat
interface, in the presence of an external electrical field. Different charge
configurations have been discussed. We firstly discussed the properties of
polyelectrolyte brushes, where the charges are uniformly distributed
along the chains, while counter-ions are in solution and free to diffuse
between the grafting surface and the second, oppositely charged electrode.
The monomer-density and counter-ion distributions
have been obtained by numerical self-consistent field theory and comparison
has been made with the analytical strong-stretching theory predictions.
The effect of a uniform electric field, both in the oppositely
and similarly charged case, has been studied. We considered
charge-end-functionalized, grafted polymer chain systems. As for the case
of polyelectrolyte brushes, we studied the conformational and electrostatic
properties of the brush by means of numerical self-consistent field theory.
In the case of weakly charged, end-group-functionalized polymers we have
been able to generalize the predictions of strong-stretching theory to the
case of a uniform electric field.
We have shown how the corrections to the MWC theory for the monomer-density and chain-end
distributions can be derived in closed form. The presence of an external field
converts the corresponding Abel-Volterra equation from an integral equation of the
first kind to an integral equation of the second kind. Deviations from the parabolic
profile have been obtained. A comparison of the theory with numerical self-consistent
field theory has been performed and agreement has been observed.

\begin{acknowledgments}
This research was founded by EPSRC under grant number EP/F068425/1.
I acknowledge very useful discussions with Mark Matsen who participated to
the early stages of this work and Jaeup Kim for providing the algorithm
that implements the self-consistent field calculation for neutral polymer
brushes.
\end{acknowledgments}

redpajama_set_name: RedPajamaArXiv