# Evolutionary Carrier Selection for Shared Truck Delivery Services

Rakpong Kaewpuang1, Suttinee Sawadsitang2∗, Dusit Niyato1, and Han Yu1
∗Corresponding author (e-mail: suttinee.s@cmu.ac.th). 1School of Computer
Science and Engineering, Nanyang Technological University
2College of Art, Media and Technology, Chiang Mai University

###### Abstract

With multiple carriers in a logistics market, customers can choose the best
carrier to deliver their products and packages. In this paper, we present a
novel approach of using the stochastic evolutionary game to analyze the
decision-making of the customers using the less-than-truckload (LTL) delivery
service. We propose inter-related optimization and game models that allow us
to analyze the vehicle routing optimization for the LTL carriers and carrier
selection for the customers, respectively. The stochastic evolutionary game
model incorporates a small perturbation of customers’ decision-making which
exists due to irrationality. The solution of the stochastic evolutionary game
in terms of stochastically stable states is characterized by using the Markov
chain model. The numerical results show the impact of carriers’ and customers’
parameters on the stable states.

###### Index Terms:

Evolutionary game, evolutionary equilibrium, integer programming, vehicle
routing problem.

##  I Introduction

Product and package delivery management is a classical but important problem
in the logistics industry. Especially, with the tremendous popularity of
online shopping, package delivery demand increases exponentially. Together
with fierce competition among logistics businesses to attract customers,
various challenging issues emerge from a classical vehicle routing problem
(VRP). Specifically, instead of just focusing on finding optimal vehicle
routes, e.g., minimizing traveling costs, competition and cooperation must be
considered when considering logistics resource/vehicle management. Game theory
becomes a suitable mathematical tool to analyze various logistics scenarios in
which each logistics company can adapt its strategy to reach the best payoff,
for example, defined in terms of utility or profit, given the strategies of
other companies. Some applications of game theory in logistics exist, e.g.,
[1] [2]. In the logistics industry, multiple carriers own vehicles/trucks for
package delivery. Customers can freely choose the carrier to meet their
demands and requirements. While the major factor in deciding on the carrier to
use is a price or service fee, recently performance and reliability of the
delivery service become increasingly concerning. For example, customers may
not always choose the cheapest carrier if their package delivery is time-
sensitive. Therefore, customers can rationally select the carrier based on
their preference, defined as a utility. The aforementioned scenario can be
generally called the carrier selection problem. Again, the game theory appears
to be a promising approach for analyzing and obtaining the solution of this
problem. However, due to a large number of customers, the standard
noncooperative game, which is tractable only for a few players, may not be
suitable, and we have to resort to the other game model for solving the
carrier selection problems of the customers.

In this paper, we focus on the carrier selection problem in the less-than-
truckload (LTL) environment where a small freight or freight that does not
require a full space of a vehicle is delivered [3]. The customers, i.e.,
players, can adaptively select LTL carriers111In the rest of the paper, we use
“LTL carrier” and “carrier” interchangeably. to achieve the highest payoff
given the decision of the other customers. The payoff depends on the delivery
performance, i.e., delay, which is again related to the delivery route.
Considering this fact, we introduce the stochastic evolutionary game to
analyze the carrier selection problem. The evolutionary game is applied as it
has been shown to be tractably applicable for a large number of players.

Additionally, a small perturbation due to the irrational decision-making of
the players can be also captured in the model. The game model can show the
impact of the number of customers choosing a certain carrier on the delivery
performance, i.e., average delay. Moreover, the impact of the physical
parameter, i.e., the location of a carrier, on the stochastically stable state
of the game is presented.

The major contributions of this paper can be summarized as follows:

  * •

We propose inter-related optimization and game models to analyze the vehicle
routing optimization for the LTL carriers and carrier selection for the
customers, respectively.

  * •

In the proposed models, we formulate and solve the integer programming (IP)
model to obtain the minimum total traveling cost of carriers. We apply the
stochastic evolutionary game model to achieve stable states that customers
obtain the minimum delivery delay and service fee.

  * •

We show the superiority of the proposed models by evaluating the performance
with a real trace data and comparing the solution of the proposed models with
the baseline model.

##  II Related Work

Optimizing the performance of package delivery from a depot or warehouse to
customers and vice versa, e.g., to minimize the cost and delay, is a classical
problem known as a vehicle routing problem with private fleet and common
carrier (VRPPC) [3]. For example, the author in [3] proposed a truckload and
less-than-truckload heuristic algorithm to solve the VRPPC. The objective of
the algorithm is to minimize the total transportation cost by selecting
private and outsourcing trucks for delivery. To tackle the problems, various
heuristic algorithms were proposed. However, the carriers have limited
resources in terms of vehicles and drivers. Game theory has been applied to
address various issues in logistics. For example, in [1], the authors
presented the routing game in a heterogeneous environment, i.e., the players
as drivers and vehicles. The players choose the number of flows to reach the
Nash equilibrium. A similar flow game was analyzed in [2] in which the impact
of GPS information available to the drivers is taken into account.

However, the carrier selection problem by customers to deliver their packages
with a minimum delay has not been studied before. Especially, when the
customers can adjust their decisions adaptively considering the decisions of
other customers in the population. The evolutionary game is a promising tool
for analyzing this problem, and thus it is the main focus of this paper.

##  III System Model and Assumptions

(a) LTL carrier selection.

(b) Reference locations.

Figure 1: (a) LTL carrier selection and (b) Reference locations of two
carriers and three terminals.

We consider the carrier selection scenario that the customers choose an LTL
carrier to deliver packages to the terminals, which is illustrated in Fig.
1(a). In particular, the carrier collects packages from customers, e.g., from
a collection point. Then the carrier loads the packages into an LTL truck and
delivers the packages to the terminals. The carrier can optimize the delivery
of the packages, which can be considered as the vehicle routing problem (VRP)
with the fixed number of trucks. The objective is to minimize the delivery
cost for the carrier. Alternatively, there are multiple carriers, and the
customers select the carrier based on their (dis)utility which depends on
service fee and delivery delay. Specifically, the customers prefer the carrier
with the smallest fee and delay.

For the carrier selection, we define sets and notations as follows. Let
$\mathcal{T}$ denote the set of terminals that LTL carriers deliver packages
to, i.e., $\mathcal{T}=\\{1,2,3,\dots,T\\}$ where $T$ is the maximum number of
terminals. Let $\mathcal{N}$ denote the set of LTL carriers, i.e.,
$\mathcal{N}=\\{1,2,3,\dots,N\\}$ where $N$ is the maximum number of LTL
carriers. Let $\mathcal{S}$ denote the set of customers, i.e.,
$\mathcal{S}=\\{1,2,3,\dots,S\\}$ where $S$ is the maximum number of
customers. In the carrier selection scenario, we are interested in analyzing
how the customers choose the carrier. Therefore, in the following sections, we
introduce the vehicle routing optimization model for the carriers and the
stochastic evolutionary game model for the customers. Given the number of
customers choosing a certain carrier, the vehicle routing optimization model
is used to obtain the optimal route and the corresponding delivery delay of
the carrier. Given the delivery delay, the stochastic evolutionary game model
is applied to obtain the decision of customers selecting any carrier.

##  IV Vehicle Routing Optimization

We propose the IP model for vehicle routing optimization for a certain LTL
carrier. The set of vehicles of the carrier is denoted by $\mathcal{V}$ and
the set of terminals to be visited by the vehicles is denoted by
$\mathcal{L}$. The IP model is presented in (1)-(13). The objective function
in (1) is to minimize the total traveling cost of the carrier.

$\displaystyle\min_{x_{v},y_{i,j,v},z_{i,v}}\sum_{v\in\mathcal{V}}x_{v}C^{\mathrm{inc}}_{v}$
(1)
$\displaystyle+\sum_{v\in\mathcal{V}}\sum_{i\in\mathcal{L}}\sum_{j\in\mathcal{L},i\neq
j}y_{i,j,v}C^{\mathrm{tvc}}_{i,j,v}+\sum_{v\in\mathcal{V}}\sum_{i\in\mathcal{U}}z_{i,v}C^{\mathrm{mcn}},$
subject to $\displaystyle\sum_{j\in\mathcal{L}}y_{i,j,v}=z_{i,v},\forall
i\in\mathcal{U},\forall v\in\mathcal{V},$ (2)
$\displaystyle\sum_{j\in\mathcal{L}}y_{j,i,v}=z_{i,v},\forall
i\in\mathcal{U},\forall v\in\mathcal{V},$ (3)
$\displaystyle\sum_{k\in\mathcal{V}}z_{0,v}\leq|\mathcal{V}|,$ (4)
$\displaystyle\sum_{v\in\mathcal{V}}\sum_{i\in\mathcal{L},i\neq
j}y_{i,j,v}=1,\forall j\in\mathcal{U},$ (5)
$\displaystyle\sum_{v\in\mathcal{V}}\sum_{j\in\mathcal{L},j\neq
i}y_{i,j,v}=1,\forall i\in\mathcal{U},$ (6) $\displaystyle
D^{\mathrm{dem}}_{i}\leq q_{i}\leq C^{\mathrm{cap}}_{v},\forall
i\in\mathcal{U},\forall v\in\mathcal{V},$ (7) $\displaystyle q_{i}\leq
C^{\mathrm{cap}}_{v}+y_{0,i,v}(D^{\mathrm{dem}}_{i}-C^{\mathrm{cap}}_{v}),\forall
v\in\mathcal{V},\forall i\in\mathcal{U},$ (8) $\displaystyle q_{j}\geq
q_{i}+D^{\mathrm{dem}}_{j}-C^{\mathrm{cap}}_{v}+(y_{i,j,v}C^{\mathrm{cap}}_{v})$
(9)
$\displaystyle+y_{j,i,v}(C^{\mathrm{cap}}_{v}-D^{\mathrm{dem}}_{j}-D^{\mathrm{dem}}_{i}),\forall
v\in\mathcal{V},\forall i,j\in\mathcal{U},i\neq j,$ $\displaystyle
x_{v}=z_{0,v},\forall v\in\mathcal{V},$ (10) $\displaystyle q_{i}\geq
0,\forall i\in\mathcal{U},$ (11) $\displaystyle z_{i,v}\in\\{0,1\\},\forall
i\in\mathcal{U},\forall v\in\mathcal{V},$ (12) $\displaystyle
x_{v},y_{i,j,v}\in\\{0,1\\},\forall v\in\mathcal{V},\forall
i,j\in\mathcal{L},i\neq j.$ (13)

The objective function in (1) consists of three terms. The first term
$\sum_{v\in\mathcal{V}}x_{v}C^{\mathrm{inc}}_{v}$ is the total initial cost.
$C^{\mathrm{inc}}_{v}$ is the vehicle- and driver-based cost [10]. $x_{v}$ is
the binary decision variable in which $x_{v}=1$ if vehicle $v\in{\mathcal{V}}$
of the carrier is used to deliver packages, and $x_{v}=0$ otherwise. The
second term
$\sum_{v\in\mathcal{V}}\sum_{i\in\mathcal{L}}\sum_{j\in\mathcal{L},i\neq
j}y_{i,j,v}C^{\mathrm{tvc}}_{i,j,v}$ is the total traveling cost.
$C^{\mathrm{tvc}}_{i,j,v}$ is the traveling cost of vehicle $v$ moving from
terminal $i$ to terminal $j$. This traveling cost is calculated from
$C^{\mathrm{tvc}}_{i,j,v}=P^{\mathrm{ppm}}D^{\mathrm{dst}}_{i,j}$, where
$P^{\mathrm{ppm}}$ is price per mile [10] and $D^{\mathrm{dst}}_{i,j}$ is a
distance from terminal $i$ to terminal $j$. $y_{i,j,v}$ is the binary decision
variable in which $y_{i,j,v}=1$ if vehicle $v$ delivers packages to terminal
$j$ after terminal $i$, and $y_{i,j,v}=0$ otherwise. The third term
$\sum_{v\in\mathcal{V}}\sum_{i\in\mathcal{U}}z_{i,v}C^{\mathrm{mcn}}$ is the
total fixed cost. $C^{\mathrm{mcn}}$ is the miscellaneous cost of visiting a
terminal by a vehicle, e.g., package handling cost. $z_{i,v}$ is the binary
decision variable in which $z_{i,v}=1$ if vehicle $v$ delivers packages to
terminal $i$, and $z_{i,v}=0$ otherwise.

The constraints are described as follows: (2) and (3) ensure that terminal $i$
is visited only once by vehicle $v$. For these constraints, the vehicle $v$
has to deliver all packages required by customers to terminal $i$ [4],[5] .
(4) guarantees that the number of vehicles used for package delivery does not
exceed the number of available vehicles of a carrier. (5) and (6) guarantee
that only one vehicle travels from terminal $i$ to terminal $j$. (7) indicates
that the delivery to terminal $i$ (i.e., $q_{i}$) has to meet the demand
[4],[5]. Additionally, the delivery must not exceed the capacity of vehicle
$v$. The demand of terminal $i$ and the capacity of vehicle $v$ are denoted by
$D^{\mathrm{dem}}_{i}$ and $C^{\mathrm{cap}}_{v}$, respectively. $q_{i}$ is
the total number of packages to be delivered on the route that terminal $i$ is
in. Note that, in practice, the customers’ demands are packed into packages
before delivery. In addition, the vehicle’s capacity can be examined by the
driver of the vehicle. If the vehicle’s capacity is not enough for packages,
an additional vehicle will handle packages. (8) is to indicate whether or not
terminal $i$ is the first terminal, the package of which is to be delivered on
the route. (9) applies, if terminal $i$ is not the first terminal of the
route, to ensure that the number of delivered packages at terminal $j$ (i.e.,
$q_{j}$) must be more than or equal to the sum of the number of delivered
packages on the route from the depot to terminal $i$ and the number of
packages from terminal $j$ (i.e., $D^{\mathrm{dem}}_{j}$). (10), (12), and
(13) represent the binary integer constraints of the corresponding variables.
(11) is the non-negative integer constraint of the corresponding variables.

The complexity of the IP model presented in (1)-(13) is
$\mathcal{O}(|\mathcal{L}|^{2})$ [6] where $\mathcal{O}(\cdot)$ and
$|\mathcal{L}|$ are big $\mathcal{O}$ notation and the maximum number of
terminals to be visited by the vehicles, respectively. This complexity means
that the computational time of solving the IP model to achieve the solution
increases exponentially when the number of terminals increases.

The optimal solutions of (1)-(13) are $x_{v}^{\ast},y_{i,j,v}^{\ast}$, and
$z_{i,v}^{\ast}$ which yield the lowest delivery cost of the carrier. The
traveling delay is then calculated as follows:

$\Phi(\mathcal{L})=\sum_{i,j\in\mathcal{L}}\frac{D^{\mathrm{dst}}_{i,j}}{V_{i,j,v}^{\mathrm{spd}}}y_{i,j,v}^{\ast},$
(14)

where $V_{i,j,v}^{\mathrm{spd}}$ is the average speed of vehicle $v$ traveling
from terminal $i$ to terminal $j$, and $D^{\mathrm{dst}}_{i,j}$ is the
distance between terminals $i$ and $j$. This traveling delay of the carrier’s
vehicle is defined as the function of the set of terminal $\mathcal{L}$. The
traveling delay is one of the factors that the customers are concerned about
in selecting the carrier with which the stochastic evolutionary game to
analyze the decision will be presented.

##  V Stochastic Evolutionary Game Formulation

In this section, we formulate the stochastic evolutionary game to model the
carrier selection of customers. The game has four components as follows:
Players are the customers. Population is the set of customers. Strategy of
customers is to select one of the LTL carriers. The set of carriers, i.e.,
strategy space, is denoted by ${\mathcal{C}}$. Negative payoff is defined by
its disutility which is a function of the service fee and average package
delivery delay. The disutility of strategy $c\in{\mathcal{C}}$, i.e., a
customer selects carrier $c$, is defined as follows:

$U_{c}(n_{c})=F_{c}+\pi(n_{c}),$ (15)

where $n_{c}$ is the number of customers selecting carrier $c$, $F_{c}$ is the
service fee, and $\pi(n_{c})$ is the average delivery delay of carrier $c$. By
economic implication, this disutility implies that customers will obtain
benefits (e.g., high reliability due to trustfulness and preference from
terminals) by selecting appropriate carriers. In addition, carriers will
achieve the minimum traveling cost due to balanced workloads.

The average delivery delay given $n$ customers choosing the carrier is
expressed as follows:

$\pi(n)=nU+\overbrace{\sum_{i_{1}=1}^{T}\sum_{i_{2}=1}^{T}\cdots\sum_{i_{n}=1}^{T}}^{n}\left(\prod_{j=1}^{T}\left(\mathbb{P}_{j}\right)^{M_{j}}\Phi(\mathcal{L}_{c}(\\{i_{1},\ldots,i_{n}\\})\right)$
(16)

where

$\displaystyle M_{j}$ $\displaystyle=$
$\displaystyle\sum_{s=1,i_{s}=j}^{n}1,\quad\quad j\in\mathcal{T},$ (17)
$\displaystyle\mathcal{L}_{c}(\\{i_{1},\ldots,i_{n}\\})$ $\displaystyle=$
$\displaystyle\Big{\\{}\\{j\\}\Big{|}j=1,\ldots,T;M_{j}>0\Big{\\}}.$ (18)

$U$ is the package unloading time per customer, and hence $nU$ is the total
unloading time for $n$ customers. In this paper, we assume that the package
unloading time is known in advance from evaluating statistical information
[8]. $i_{s}$ is the terminal to which the package of customer $s$ is to be
delivered. $M_{j}$ is the number of customers, the packages of which are to be
delivered to terminal $j$. $\Phi(\mathcal{L}_{c})$ is the delivery delay of
carrier $c$ given the set of terminals to visit $\mathcal{L}_{c}$ as it is
obtained from (14) by solving vehicle routing optimization defined in
(1)-(13). $\mathbb{P}_{j}$ is the probability that the customer wants to
deliver a package to terminal $j$.

In the game, each customer compares its disutility with the average disutility
of all carriers. The customer switches to the new carrier if its current
disutility is higher than that of the new carrier. With this strategic action
of the customer population, we can model the strategic evolution of the
customers using a Markov chain. The state space of the Markov chain is defined
as follows:

$\Omega=\Big{\\{}(n_{1},\ldots,n_{N})\Big{|}\sum_{c=1}^{N}n_{c}=S\Big{\\}},$
(19)

where again $n_{c}$ is the number of customers selecting carrier $c$. To deal
with the large state space of the Markov chain, the most probable state
generation algorithm and techniques of decomposition, lumping, and truncation
can be applied [7]. The transition rate from state $(n_{1},\ldots,n_{N})$ to
state $(n^{\prime}_{1},\ldots,n^{\prime}_{N})$ is expressed as follows:

$Q_{(n_{1},\ldots,n_{N}),(n^{\prime}_{1},\ldots,n^{\prime}_{N})}=\left\\{\begin{array}[]{ll}Kn_{c}(U_{c}-U_{c^{\prime}}),&\mbox{if
}U_{c}>U_{c^{\prime}},\\\ \epsilon,&\mbox{otherwise},\end{array}\right.$ (20)

where $K$ is a certain constant related to the decision-making rate of
customers, and $\epsilon$ is a small perturbation rate related to the
irrational decision-making of the players. In particular, the customer can
make a small mistake to choose the carrier which yields higher disutility.
From (20), the transition rate depends on the difference between the
disutilities of two carriers, and the number of customers currently choosing
an inferior carrier, i.e., with larger disutility. Note that the self-
transition rate of each state [4] is expressed as follows:

$\displaystyle Q_{(n_{1},\ldots,n_{N}),(n_{1},\ldots,n_{N})}=$
$\displaystyle\quad\;\;-\sum_{(n^{\prime}_{1},\ldots,n^{\prime}_{N})\neq(n_{1},\ldots,n_{N})}Q_{(n_{1},\ldots,n_{N}),(n^{\prime}_{1},\ldots,n^{\prime}_{N})}.$
(21)

Let ${\mathbf{Q}}$ denote the transition matrix of the Markov chain, the
element of which is
$Q_{(n_{1},\ldots,n_{N}),(n^{\prime}_{1},\ldots,n^{\prime}_{N})}$, the
stationary probability of state $(n_{1},\ldots,n_{N})$ of the Markov chain can
be obtained by solving
$\vec{\boldsymbol{\beta}}^{\top}{\mathbf{Q}}=\vec{\mathbf{0}}^{\top}$ and
$\vec{\boldsymbol{\beta}}^{\top}\vec{\mathbf{1}}=1$, where
$\vec{\boldsymbol{\beta}}$ is the vector of stationary probability,
$\vec{\mathbf{0}}$ is the vector of zeros, and $\vec{\mathbf{1}}$ is the
vector of ones. For $\epsilon\rightarrow 0$, if the stationary probability of
state $(n^{\dagger}_{1},\ldots,n^{\dagger}_{N})$ does not approach zero, then
the state is stochastically stable. In other words, we can consider the state
$(n^{\dagger}_{1},\ldots,n^{\dagger}_{N})$ to be a stable point of the game,
i.e., an equilibrium.

##  VI Performance Evaluation

###  VI-A Parameter Setting

To simplify the performance of our proposed approach, we consider the carrier
selection of 16 customers. There are two LTL carriers, and three terminals
(e.g., an airport, train station, and bus terminal). Note that our proposed
approach can handle many customers (e.g., 100 customers) and terminals (e.g.,
50 terminals). The locations of the two carriers and three terminals are shown
in Fig. 1(b) in which the unit is in kilometer. We adopt the data provided by
I.logistics Pte Ltd company in Singapore [8]. Carriers 1 and 2 have one LTL
vehicle. The vehicle capacities of carriers 1 and 2 are 20 and 30 packages,
respectively. The package unloading time is 5 minutes per customer. The
average speed of the vehicles is 40 km/h [8]. We implement and solve the
vehicle routing optimization model of each LTL carrier using the GAMS/CPLEX
solver [9]. The traveling cost of the LTL vehicle is 0.982$ per kilometer
[10]. We assume that all the carriers charge the same fee due to market
competition. Thus, the customers focus only on the average delivery delay. For
the Markov chain, we set $K=1$ and $\epsilon=10^{-3}$.

TABLE I: Traveling delay of carriers 1 and 2 to different sets of terminals.

| Traveling delay (minutes)
---|---
Terminal | Carrier 1 | Carrier 2
1 | 2.400 | 18.228
2 | 8.625 | 17.376
3 | 9.450 | 17.340
$1\rightarrow 2$ | 13.425 | 29.253
$1\rightarrow 3$ | 14.250 | 30.078
$2\rightarrow 3$ | 26.700 | 35.415
$1\rightarrow 2\rightarrow 3$ | 31.500 | 47.328

Table I shows the traveling delay of each carrier visiting different sets of
terminals. Clearly, the traveling delay of carrier 1 is smaller than those of
carrier 2 due to carrier 1’s location which is closer to terminals.

(a) Average delay without rational decision.

(b) Average delay under different customers.

(c) Stochastically stable states of the evolutionary game.

(d) Average number of customers choosing different carriers.

(e) Total traveling cost comparison of carrier 1

Figure 2: (a) Average delay of carriers 1 and 2 without rational decision-
making of customers, (b) Average delay of carriers 1 and 2 under the different
number of customers choosing carrier 1, (c) Stochastically stable states of
the evolutionary game, (d) The average number of customers choosing different
carriers when the distance from carrier 2 to carrier 1 is varied, and (e) The
total traveling cost comparison between the proposed model and the baseline
model.

Fig. 2(a) shows the average delay of carriers under different carrier
locations when all the customers irrationally select either carrier 1 or
carrier 2 to deliver packages to terminals. In other words, all the customers
only select carrier 1 or carrier 2 since they do not consider the performance
of carriers. In Fig. 2(a), the average delay of carriers for package delivery
dramatically increases when the number of customers increases. This is due to
the fact that all the customers make the irrational decision to only select
carrier 1 or carrier 2, and consequently the average delay of the carrier to
deliver customers’ packages to the terminals increases. Additionally, the
total package unloading time of the carrier increases.

From the result shown in Fig. 2(a), this is a reason why our proposed approach
based on the evolutionary game is beneficial to analyze the carrier selection
problem to minimize the average delay of carriers for package delivery.

We next consider the population of customers. Fig. 2(b) shows the average
delivery delay from both carriers when the number of customers choosing the
carriers is varied. For example, if the number of customers choosing carrier 1
is 2, then the number of customers choosing carrier 2 will be 14. From Fig.
2(b), when the number of customers choosing each carrier increases, the
average delivery delay increases. This is due to the fact that, with more
customers, the chance that the carrier has to visit more terminals increases.
Additionally, the total package unloading time increases. Therefore, the
decision of the customer to choose one carrier to minimize its disutility also
depends on the decisions of other customers. It is important to analyze the
decision-making of the customers and obtain stable states of the carrier
selection.

From the evolutionary game for the carrier selected by customers, the
stochastically stable states, i.e., the number of customers choosing carriers,
are shown in Fig. 2(c). We make the following observations. Firstly,
stochastically stable states have the non-zero stationary probability as
$\epsilon\rightarrow 0$. Secondly, the stochastically stable states associate
with the number of customers choosing the carriers that result in a similar
average delivery delay. Thirdly, there are multiple stochastically stable
states. This is due to the fact that these states result in slightly higher or
lower disutility for the customers. Therefore, the customers can switch among
the states to achieve identical average disutility. In stochastically stable
states shown in Fig. 2(c), the first stable state means that eight customers
select carrier 1 and eight customers select carrier 2. The second stable state
means that seven customers select carrier 1 and nine customers select carrier
2.

Next, we evaluate the impact of the location of a carrier on the decisions of
customers. We assume that initially carrier 2 is at the same location as
carrier 1 (see Fig. 1(b)). Then, the location of carrier 2 is shifted to the
east (the right side) of carrier 1 horizontally, which implies that carrier 2
becomes farther away from the terminals. Fig. 2(d) shows the average number of
customers choosing the carriers when the distance from carrier 1 to carrier 2
is varied. Clearly, as the distance becomes larger, carrier 2 takes more time
to reach all three terminals, and hence the delivery delay increases.
Consequently, due to poorer performance, customers deviate to choose carrier
1. This result shows evidently the usefulness of the proposed stochastic
evolutionary game model that can capture the physical parameters of the
vehicle routing problems.

Finally, we compare the traveling costs of carrier 1 from the proposed model
and the baseline model. Figure 2(e) shows the performance comparison of the
proposed model and the baseline model by varying the number of customers. For
the baseline model, we consider that all customers make the irrational
decision to select carrier 1 only. Clearly, the proposed model can achieve a
significantly lower traveling cost of carrier 1 than the baseline model.

##  VII Conclusion

Due to limited resources, i.e., vehicles of carriers, package delivery
performance can depend on the number of customers served by the carrier for
package delivery. In this paper, we have considered the less-than-truckload
services so that the customers can select the carrier with the best
performance. As the customers can adapt their decision dynamically, we have
presented the stochastic evolutionary game model to analyze the decision-
making of the population of customers. We have first formulated vehicle
routing optimization for obtaining optimal delivery routes to minimize the
cost of a carrier. We have then formulated the evolutionary game model to
obtain the stochastically stable states in terms of the number of customers
choosing different carriers. The disutility of the customers in the
evolutionary game is a function of delivery delay obtained from the vehicle
routing optimization problem. The numerical results have verified the
existence of the stochastically stable states for which all customers receive
identical average disutility. Interestingly, in one of the stochastically
stable states, seven customers select carrier 1 while nine customers select
carrier 2 due to the identical average disutility.

For future work, we will consider multiple populations of customers, the
(dis)utility of which are defined differently. The service fee can be
optimized to maximize the profits of the carriers. In addition, we will
consider the effect of the average speed of the vehicle on the equilibrium
solution.

## Acknowledgment

This research is supported by the National Research Foundation Singapore and
DSO National Laboratories under the AI Singapore Programme (AISG Award No:
AISG2-RP-2020-019); the National Research Foundation (NRF), Singapore and
Infocomm Media Development Authority under the Future Communications Research
Development Programme (FCP); Energy Research Test-Bed and Industry Partnership
Funding Initiative, part of the Energy Grid (EG) 2.0 programme; DesCartes and
the Campus for Research Excellence and Technological Enterprise (CREATE)
programme; Alibaba Group through Alibaba Innovative Research (AIR) Program and
Alibaba-NTU Singapore Joint Research Institute (JRI); and Nanyang
Technological University, Nanyang Assistant Professorship.

## References

  * [1] F. Farokhi, W. Krichene, A. M. Bayen and K. H. Johansson, “A heterogeneous routing game,” Annual Allerton Conference on Communication, Control, and Computing, Monticello, IL, pp. 448-455, 2013.
  * [2] A. Aswani and C. Tomlin, “Game-theoretic routing of GPS-assisted vehicles for energy efficiency,” Proceedings of the 2011 American Control Conference, San Francisco, CA, pp. 3375-3380, 2011.
  * [3] C.-W. Chu, “A heuristic algorithm for the truckload and less-than-truckload problem”, JOR, vol. 165, issue 3, pp. 657-667, 2005.
  * [4] R. Kaewpuang, D. Niyato, P. -S. Tan, and P. Wang, “Cooperative Management in Full-Truckload and Less-Than-Truckload Vehicle System,” in IEEE TVT, vol. 66, no. 7, pp. 5707-5722, 2017.
  * [5] X. Li, 2015, Capacitated Vehicle Routing Problem with Time Windows: A Case Study on Pickup of Dietary Products in Nonprofit Organization, Arizona State University, Arizona, USA.
  * [6] J. K. Lenstra and A. H. G. Rinnooy Kan, “Complexity of vehicle routing and scheduling problems”, Networks, vol. 11, no. 2, pp. 221–227, 1981.
  * [7] E. de Souza e Silva and P. M. Ochoa, “State space exploration in Markov models”, ACM SIGMETRICS joint international conference on Measurement and modeling of computer systems, pp. 152–166, 1992.
  * [8] I.Logistics, http://www.ilogpteltd.com.sg , 2022.
  * [9] GAMS Solvers, https://www.gams.com , 2022.
  * [10] An analysis of the operational costs of trucking: 2015 update, http://atri-online.org/wp-content/uploads/2015/09/ATRI-Operational-Costs-of-Trucking-2015-FINAL-09-2015.pdf, 2015.


# Persistent electrical doping of Bi2Sr2CaCu2O8+x mesa structures

Holger Motzkau  Thorsten Jacobs  Sven-Olof Katterwe  Andreas Rydh  Vladimir M.
Krasnov<EMAIL_ADDRESS>Department of Physics, Stockholm
University, AlbaNova University Center, SE – 106 91 Stockholm, Sweden

###### Abstract

Application of a significantly large bias voltage to small Bi2Sr2CaCu2O8+x
mesa structures leads to persistent doping of the mesas. Here we employ this
effect for analysis of the doping dependence of the electronic spectra of
Bi-2212 single crystals by means of intrinsic tunneling spectroscopy. We are
able to controllably and reversibly change the doping state of the same single
crystal from underdoped to overdoped state, without changing its chemical
composition. It is observed that such physical doping is affecting
superconductivity in Bi-2212 similar to chemical doping by oxygen impurities:
with overdoping the critical temperature and the superconducting gap decrease,
with underdoping the $c$-axis critical current rapidly decreases due to
progressively more incoherent interlayer tunneling and the pseudogap rapidly
increases, indicative for the presence of the critical doping point. We
distinguish two main mechanisms of persistent electric doping: (i) even in
voltage contribution, attributed to a charge transfer effect, and (ii) odd in
voltage contribution, attributed to reordering of oxygen impurities.

###### pacs:

74.72.Gh, 74.62.Dh, 74.55.+v, 74.72.Kf

##  I Introduction

High temperature superconductivity (HTSC) in cuprates occurs as a result of
doping of a parent antiferromagnetic Mott insulator and properties of cuprates
change significantly with doping. Superconductivity in overdoped cuprates is
fairly well described by the conventional BCS-type second-order phase
transition TallonPhC ; SecondOrder ; MR . But properties of underdoped
cuprates are abnormal due to the persistence of the normal state pseudogap,
strong superconducting fluctuations, or possibly preformed pairing, Kaminski
and magnetism Kerr_Shen  at $T>T_{\text{c}}$. There are indications that the
transition from the normal to the abnormal behavior occurs abruptly at a
critical doping point TallonPhC ; Doping ; Balakirev ; Kartsovnik ;
LeBoeuf2011 . This may be a consequence of the quantum phase transition - a
phase transition, which occurs at $T=0$, in frustrated systems as a result of
a competition of coexisting order parameters. The coexistence of
superconductivity at $T<T_{\text{c}}$ with the pseudogap KrasnovPRL2000 ;
Lee2007 ; Bernhard , charge and spin density order Kartsovnik ; LeBoeuf2011 ;
Nernst2010  was indeed reported by several techniques. Clearly, detailed
doping-dependent studies are needed both for understanding the puzzling nature
of HTSC in cuprates and for the development of novel HTSC materials.

Usually, the mobile carrier concentration is controlled by chemical doping via
chemical substitution, or in case of cuprates also by variation of the oxygen
content via appropriate annealing and subsequent quenching to room
temperature. This allows an accurate control of the chemical composition, but
less so of the local arrangement of impurities and disorder, which is equally
important for cuprates DisorderBi2201 . For example, it is well established
that properties of the $\mathrm{YBa_{2}Cu_{3}O_{6+x}}$ compound strongly
depend not only on the average oxygen concentration but also crucially on the
order/disorder of oxygen atoms in Cu-O-Cu chains Chains ; Raman . Therefore,
analysis of the doping phase diagram of cuprates requires accurate control of
both, the concentration and the microscopic structure of impurities.

The carrier concentration can be also varied via two physical doping
processes, well established for semiconductors: photo-doping Kudinov ;
Schuller ; PhotoScience2011  and through the electric-field effect Ahn ;
Moeckly ; Sorkin ; Tulina ; KovalCurrentInj ; KovalMemres . In case of
cuprates, physical doping may be persistent at low temperatures in a sense
that it is relaxing very slowly after removing the light Kudinov ; Schuller
or field Moeckly ; Tulina ; KovalCurrentInj ; KovalMemres . Recently, a
persistent electric doping via strong current injection was discovered
KovalCurrentInj ; KovalMemres . It is resembling a resistive switching
phenomenon in memristor devices Memristor1 ; Memristor2  and is related to
previous similar observations in point contact experiments on Bi2Sr2CaCu2O8+x
(Bi-2212) Tulina . Such an electric doping is reversible, reproducible and
easily controllable. It opens a possibility to analyze the doping dependence
of HTSC on one sample without changing its chemical composition
PhotoScience2011 . Despite that, there were very few direct spectroscopic
studies of cuprates employing physical doping techniques.

Intrinsic Tunneling Spectroscopy (ITS) provides an unique opportunity to probe
bulk electronic properties of HTSC SecondOrder . This technique utilizes weak
interlayer ($c$-axis) coupling in quasi two dimensional HTSC compounds, in
which mobile charge carriers are confined in CuO2 planes separated by some
blocking layer (e.g. SrO-2BiO-SrO in case of Bi-2212). This leads to a
formation of atomic scale intrinsic tunnel junctions, and to an appearance of
the intrinsic Josephson effect at $T<T_{\mathrm{c}}$ Kleiner94 ;
KrasnovPRL2000 ; Suzuki ; Katterwe2009 ; Superluminal .

In this work we employ the persistent electric doping for analysis of the
doping dependence of electronic spectra of Bi-2212 single crystals by means of
intrinsic tunneling spectroscopy SecondOrder ; Doping ; MR ; KrasnovPRL2000 ;
KatterwePRL2008 . Controllable and reversible persistent physical doping is
achieved by applying a $c$-axis voltage of a few volts to small Bi-2212 mesa
structures. Thus we are able to change the doping state of Bi-2212 single
crystals without changing its chemical composition. A wide doping range from a
moderately underdoped to strongly overdoped state could be reached. It is
observed that the physical doping is affecting the intrinsic tunneling
characteristics of Bi-2212 similar to chemical doping Doping . With overdoping
the critical temperature and the superconducting gap decrease. With
underdoping the pseudogap rapidly increases, indicative for the presence of
the quantum critical doping point in the phase diagram, and the $c$-axis
critical current density rapidly decreases, indicating a progressively more
incoherent interlayer tunneling. We distinguish two main mechanisms of
persistent electric doping: (i) an even in voltage contribution, attributed to
a charge transfer effect, and (ii) an odd in voltage contribution, attributed
to reordering of oxygen impurities.

The paper is organized as follows. In Sec. II we make a brief overview of
physical doping mechanisms of cuprates. Sec. III provides experimental
details. In Sec. IV we present the main experimental results and in Sec. V we
discuss possible mechanisms of persistent electric doping, followed by
conclusions.

##  II Physical doping of cuprates

###  II.1 Photo-doping

Photo-doping allows a wide-range variation of doping in the same sample
Kudinov ; Schuller ; PhotoScience2011 . The ordinary non-equilibrium photo-
doping is quickly relaxing because of a very short life time ($\sim$ps) of
photoinduced charge carriers PhotoScience2011 . However, in underdoped
$\mathrm{YBa_{2}Cu_{3}O_{6+x}}$ and some other cuprates a different type of
persistent photo-doping takes place Kudinov ; Schuller . It involves
significant energies $\sim$$1\,\mathrm{eV}$, which makes it metastable at low
temperatures. Several mechanisms are contributing to the persistent photo-
doping PhDopTheory , such as charge transfer, which changes the redox state of
the impurity atom Kudinov , and ordering of oxygen impurities in the lattice
Schuller . Photo-doping always leads to an increase of the doping level with
respect to the initial state.

###  II.2 Electric field effects

Electric fields may both increase or decrease the number of mobile charge
carriers, depending on the direction of the applied field Ahn . The ordinary
electric field effect is not persistent and exists only during the time an
electric field is applied. Since the electric field penetrates only to the
Thomas-Fermi charge screening length, $\lambda_{\mathrm{TF}}\lesssim
1\,\mathrm{nm}$, just a thin surface layer can be modified Shapiro .

A persistent electrostatic field-effect due to net electric polarization or
trapped charges can be realized at the interface between a superconductor and
a ferroelectric SuperFerro  or polar insulator ElStatic . This is also a
surface phenomenon, but in case of layered cuprates, which represent stacks of
metallic CuO planes sandwiched between polar-insulating layers Polariton ,
electrostatic charging of insulating layers may in principle lead to the bulk
persistent electrostatic field-effect.

Another type of a persistent and bulk electric-field effect has been observed
at large current densities Moeckly . Similar to photo-doping, it was
attributed to a charge transfer Salluzzo2008  and reordering of oxygen
impurities Sorkin . Significant oxygen mobility in intense electric fields
also leads to a resistive switching phenomenon Tulina .

###  II.3 Resistive switching in complex oxides

The resistive switching phenomenon occurs in many complex oxides and is the
basis for the development of resistive memory devices. Several mechanisms may
be involved in the resistive switching phenomenon Memristor1 ; Memristor2 ,
including a change of the redox state of some of the elements, oxygen
migration, and filament formation. Resistive switching has been observed on
depleted surfaces of Bi-2212 cuprates Tulina  and attributed to oxygen
migration. Recently, it was demonstrated that the resistive switching-like
phenomenon can be used for controllable and reversible doping of small Bi-2212
micro-structures over a wide doping range KovalCurrentInj ; KovalMemres .

##  III Experimental

Mesas were fabricated on top of freshly cleaved Bi-2212 single crystals by
means of optical lithography, Ar ion milling and focused ion beam trimming.
Four batches of crystals were used: pure near optimally doped (OP) Bi-2212,
pure strongly underdoped (UD) Bi-2212, lead-substituted
$\mathrm{Bi_{1.75}Pb_{0.25}Sr_{2}CaCu_{2}O}_{8+\delta}$ [Bi(Pb)-2212], and
yttrium-doped $\mathrm{Bi_{2}Sr_{2}Ca_{1-x}Y_{x}Cu_{2}O}_{8+\delta}$
[Bi(Y)-2212]. Mesas of different sizes from $5\times 5$ to $1\times
0.5~{}\muup\mathrm{m}^{2}$ and with a different number of junctions $N=8-56$
were studied. Details of the sample fabrication can be found in Ref. Submicron
. All studied mesas exhibited a persistent doping effect upon application of a
sufficiently large bias voltage.

The samples were placed in a flowing gas cryostat and measured in a three-
probe configuration with a common top gold contact. The ground contacts for
current and voltage were provided through other mesas on the same crystal. A
Keithley K6221 current source and a FPGA-based arbitrary waveform generator
and lock-in amplifier were used to bias and measure the samples. Biasing was
done at pseudo-constant voltage, with an optional small superimposed ac
voltage to simultaneously measure the high-bias differential resistance in
addition to the dc resistance. Positive bias is defined as current going into
the mesa through the common top contact, as sketched in Fig. 1 (a).

Figure 1: (Color online). Dynamics of electric doping (a-c) for a Bi(Pb)-2212
mesa at $T=135\,\mathrm{K}$ and (d-f) for an optimally doped Bi-2212 mesa at
$T=2\,\mathrm{K}$. Panel (a) shows time evolution of the dc resistance at a
bias of $\approx 1.9\,\mathrm{V}$. The dashed line represents a stretched
exponential decay. The inset schematically shows a sketch of a mesa structure
and the direction of electric field at positive bias. (b) Demonstration of
odd-in-voltage doping: the mesa resistance increases/decreases at
positive/negative bias. Panel (c) shows the predominantly odd voltage
dependence of the logarithmic rate of the dc resistance change $\mathrm{d}\ln
R_{\mathrm{dc}}/\mathrm{d}t$. Different symbols represent different runs. The
inset in (c) demonstrates the sign change of the doping direction at higher
bias. Panel (d) demonstrates a gradual change from a negative to positive
doping rate with increasing bias voltage and time. It also demonstrates a
history dependence of the doping rate, i.e., a different sign of resistance
change at the same bias voltage, depending on the former doping treatment. The
history dependence upon several sweeps of the bias voltage are shown in (e).
Panel (f) demonstrates the predominantly even in voltage doping for OP
Bi-2212. However an asymmetry indicates the presence of a subdominant odd in
voltage doping. Note that despite a significant difference between the two
mesas, the threshold doping voltage is similar $\sim\pm 1.7\,\mathrm{V}$ (c).

##  IV Results

From chemical (oxygen) doping studies it is known that doping/undoping of
Bi-2212 is accompanied by a proportional decrease/increase of the $c$-axis
resistivity Doping . Therefore, we can control the doping state by tracing the
mesa resistance.

###  IV.1 Dynamics of electric doping

The basic features of the dynamics of the persistent electric doping are shown
in Fig. 1. Panel (a) shows the time-evolution of the Bi(Pb)-2212 mesa
resistance at a bias of $\approx$$1.9\,\mathrm{V}$ at $T=135\,\mathrm{K}$. It
is seen that the mesa resistance is decreasing with time, indicating a gradual
doping of the mesa. The doping rate decreases with time, following a stretched
exponential decay, $R=R_{0}+R_{\mathrm{d}}\exp[-(t/\tau)^{\beta}]$, shown by
the dashed line in in Fig. 1 (a), which is also typical for persistent photo-
doping Kudinov ; Schuller . The doping can be equally well performed at any
temperature from 4 to 300 K, but the rate is increasing with $T$. In most
cases we perform doping at low $T$ in order to be able to immediately probe
the superconducting characteristics. Upon reduction of the bias below the
threshold voltage the state of the mesa remains stable even at room
temperature on the time scale of several days.

The resistive change is reversed upon voltage reversal, as illustrated in
Figs. 1 (b) and (c) for the same Bi(Pb)-2212 mesa. The resistance decreases
for positive bias (electric field into the crystal) and increases for negative
bias (electric field towards the top contact), which indicates that we can
controllably and reversibly dope and undope the mesa. The doping rate and
direction depend both on the sign and the absolute value of bias voltage.

Figure 2: (Color online) Intrinsic tunneling characteristics of the same OP
Bi-2212 (a-d) and Bi(Pb)-2212 (e-g) mesas as in Fig. 1 at different doping
states. (a) $T$ dependence of the ac resistance $R_{\text{ac}}(T)$ in the
initial high-resistive state (1) and subsequent low resistance doping states
(4-6). The inset illustrates the significant enhancement of $T_{\text{c}}$ of
the depleted surface junction. Panels (b-d) represent different parts of $I-V$
characteristics in the doping states 1, 2, 3 and 6 at $T=2$ K. A variation of
the tunnel resistance, sum-gap kink and the critical current is clearly seen.
Panels (e) and (f) show $R_{\text{ac}}(T)$ for Bi(Pb)-2212 in the initial
high-resistive state (1) and subsequent low resistance doping states (4) and
(6). This mesa with a large number of junctions $N=56$ exhibits a significant
spread in the $T_{\text{c}}$ of individual junctions, seen as small resistance
drops in (f). Panel (g) shows $I-V$s of the Bi(Pb)-2212 mesa in the doping
states 2, 4 and 5. A progressive increase of the critical current and decrease
of the sum-gap kink voltage is seen.

Figure 1 (c) summarizes the bias dependence of the doping rate for the
Bi(Pb)-2212 mesa. Below the threshold voltages,
$\left|V_{\mathrm{dc}}\right|\lesssim 1.7\,\mathrm{V}$, the mesa resistance is
stable. Upon increasing the bias voltage, the resistance of the mesa starts to
gradually change at a rate that increases drastically up to $|V|\sim 2.2\,V$
as shown in the inset of Fig. 1 (c). A further voltage increase reduces the
rate and then reverses the resistance alteration rate (see the inset at
$V=-2.32\,\mathrm{V}$). The sign change of the alteration rate at high bias is
in agreement with the observations by Koval _et al._ KovalCurrentInj . The
behavior in this regime is, however, history dependent, as may be seen from
Fig. 1 (c), and the final state depends on how long time the mesas was biased
at every bias voltage. The doping process for the Bi(Pb)-2212 mesa, Fig. 1
(c), is predominantly odd in bias voltage, i.e., the direction of doping is
changed when the sign of the bias voltage is changed.

Figures 1 (d) and (e) show a detailed view of the time and voltage dependence
of doping for a near optimally doped pure Bi-2212 mesa. The top panel of (d)
shows the time evolution of the resistance of an OP Bi-2212 mesa for different
bias voltages, shown in the bottom panel. It is seen that the resistance is
constant at $V=1.7\,\mathrm{V}$, and starts to increase slowly at 1.9 V.
However, at 2 V the resistance initially increases but then starts to decrease
after a few minutes. This clearly shows that there are two counteracting
processes: a positive rate mechanism that saturates quickly and a mechanism
with negative rates that is dominating at longer times and at higher voltages.
The second process also saturates with time, which is clear from Fig. 1 (a)
and history dependent rates of Fig. 1 (e). At larger voltage, the resistance
steadily decreases at a rate which is strongly bias dependent, as shown in
Fig. 1 (f).

Figures 1 (e,f) show the bias dependence of doping and a history dependence
upon sequential voltage sweeps (e). It is seen that for the OP Bi-2212 mesa
the electric doping is predominantly even in voltage, i.e., the direction of
doping does not depend on the sign of the bias voltage. However, certain
asymmetry of the doping rate versus bias voltage characteristics in Figs. 2
(c) and (f) indicates the presence of sub-dominant even- and odd-in-voltage
contributions for Bi(Pb)-2212 and OP Bi-2212 mesas, respectively.

Figure 3: (Color online). Doping of a strongly underdoped Bi-2212 mesa. (a)
$I$-$V$ characteristics in the initial state 1 (black) and two successive
doping states 2 (red) and 3 (green). The increase of the critical current by a
factor five is seen. (b) The corresponding $\mathrm{d}I$-$\mathrm{d}V$
characteristics. The superconducting sum-gap peak moves to slightly lower
voltages and becomes sharper with doping. The $c$-axis pseudogap hump voltage
rapidly decreases with doping. (c) Zero bias ac resistivity
$R_{\text{ac}}(T)$. The $T_{\text{c}}$ has increased by about 15 K in the
state 3.

###  IV.2 Doping dependence of ITS characteristics

Using the described method, the superconducting properties of the mesas have
been altered to different intermediate doping states denoted by a successive
number. The electric doping changes all mesa characteristics: the $c$-axis
resistivity, the critical temperature $T_{\mathrm{c}}$, the $c$-axis critical
current density $J_{\mathrm{c}}$, the superconducting energy gap, $\Delta$,
the $c$-axis pseudogap and the $c$-axis resistivity in a manner very similar
to chemical (oxygen) doping Doping .

Figure 2 shows temperature dependencies of low bias ac-resistances for the
initial, high resistance state (HRS), and doped, lower resistance states
(LRS), for (a) an OP Bi-2212 and (e,f) a Bi(Pb)-2212 mesa. The $I$-$V$
characteristics of those mesas at different doping states are presented in
panels (b-d) and (g), respectively.

From Fig. 2 (a) it is seen that the initial state was characterized by the
strong thermal-activation-type increase of resistance with decreasing $T$,
typical for underdoped Bi-2212 Katterwe2009 . The general shape of the
resistive transition was described in Ref. SecondOrder . At
$T_{\mathrm{c}}\sim 82\,\mathrm{K}$ the resistance dropped to the top-contact
resistance, which originates from the first deteriorated junction between the
top CuO plane, shortly exposed to atmosphere after cleavage, and the second,
un-deteriorated CuO plane SecondOrder . Initially this junction had a very low
$T_{\mathrm{c}}^{\prime}$ and a very small critical current, $I_{\mathrm{c}}$,
as can be seen from the corresponding $I$-$V$ in panel (d). After electric
doping, the resistance in the normal state dropped almost three times and
became less semiconducting. The main $T_{\mathrm{c}}$ of the mesa changed only
slightly, indicating that the doping was changing around the optimal doping
level with the flat $T_{\mathrm{c}}$ vs. doping dependence. However, the
properties of the top junction changed drastically: the
$T_{\mathrm{c}}^{\prime}$ increased to $\sim 50\,\mathrm{K}$, and the critical
current increased ten-fold as shown in panel (d), even though it still remains
$\sim 20$ times smaller than for the rest of the junctions in the mesa, as can
be seen from panel (c). This indicates that the surface CuO plane was
initially strongly underdoped and, therefore, responded much stronger to
variation of doping, due to the steep $T_{\mathrm{c}}$ vs. doping dependence
at the underdoped side of the doping phase diagram of cuprates Doping . A
similar trend was also observed in photo-doping Kudinov ; Schuller .

Figure 4: (Color online). Intrinsic tunneling characteristics of a small
Bi(Y)-2212 mesas at low $T$ before (slightly underdoped, HRS) and after
(optimally doped, LRS) doping by a short voltage pulse. Panels (a) and (b)
show $I$-$V$ curves at a large scale and at the quasiparticle branches,
respectively. A reduction of the high bias tunnel resistance and a
simultaneous increase of the critical current in the doped state is seen.
Panel (c) shows $\mathrm{d}I/\mathrm{d}V$ ITS characteristics. A change in the
shape of the curves is seen: the peak becomes sharper and the dip-hump less
pronounced in the optimally doped state, compared to the initial underdoped
state.

Another possible reason for the stronger response of the underdoped top
junction is the larger $c$-axis resistivity of underdoped intrinsic Josephson
junctions Doping . Because of that the electric field is not uniformly
distributed along the mesa but is larger in the high-resistive top junction.
This together with the strong voltage dependence of the electric doping leads
to a faster doping of the top junction and the doping may even go in the
opposite direction with respect to the rest of the mesa.

The effect of non-uniform doping along the height of the mesa becomes more
pronounced in higher mesas with a larger number of intrinsic Josephson
junctions. This is seen from $R(T)$ for the Bi(Pb)-2212 mesa from Fig. 2 (e),
which contained a fairly large number of junctions $N\approx 56$. It is seen
that after doping some junctions retained the initial $T_{\mathrm{c}}\sim
90\,\mathrm{K}$, but some were very strongly overdoped to $T_{\mathrm{c}}\sim
30\,\mathrm{K}$. Fig. 2 (f) shows $R(T)$ at the intermediate doping state 4.
Small drops represent critical temperatures of individual junctions in the
mesa. Apparently there is a gradual distribution of $T_{\mathrm{c}}$ along the
height of the mesa.

Doping of Bi-2212 leads to a rapid increase of the $c$-axis critical current
density Doping . This is clearly seen from $I$-$V$ curves for the Bi-2212 (OP)
mesa shown in Figs. 2 (b-d). The increase of $I_{\mathrm{c}}$ of the surface
junction is shown in panel (d). One-by-one switching of the rest of the
junctions into the resistive state at $I>I_{\mathrm{c}}$ leads to the
appearance of multiple-quasiparticle (QP) branch structures in the $I$-$V$
curves. The corresponding critical current at the first QP branch is shown in
panel (c). It also strongly increases with doping. The same effect doubles
$I_{\mathrm{c}}$ in the Bi(Pb)-2212 mesa, in Fig. 2 (g).

From Fig. 2 (b) it is seen that the $I$-$V$ curves exhibit a kink at large
bias, followed by an ohmic tunnel resistance. The kink represents the sum-gap
singularity in superconducting tunnel junctions at $V=2\Delta/e$ per junction
SecondOrder ; MR . This is the basis of the ITS technique, which allows
analysis of the superconducting energy gap $\Delta$ in the bulk of the Bi-2212
single crystal.

Accurate analysis of the electronic spectra with the ITS technique requires
mesas with a small area and a small number of identical junctions. This is
needed for avoiding possible artifacts, associated with self-heating, in-plane
non-equipotentiality, and spread in junction parameters Heating_PRL2005 ;
SecondOrder . This is particularly important for the analysis of the
genuineshape of tunneling characteristics, which remains a controversial issue
Comment2 . Even though the $I$-$V$ characteristics in Figs. 2 (b) and (g) are
distorted by self-heating, as evident from a back-bending at large bias, the
general trend for a variation of the sum-gap kink with doping is clearly seen:
The superconducting gap decreases and the sum-gap kink becomes sharper with
(over-)doping. This qualitative conclusion is not affected by self-heating
because the dissipation power at the kink decreases with decreasing resistance
and becomes smaller with subsequent doping. Thus, with over-doping the
superconducting sum-gap singularity becomes sharper and moves to lower
voltages despite the progressive reduction of self-heating. This clearly
reveals the doping-variation of the genuine $c$-axis tunneling characteristics
Comment2 . A similar tendency was observed by other techniques, including the
angular resolved photoemission spectroscopy ARPES_doping , scanning tunneling
spectroscopy STS_inhomo , and tunneling spectroscopy on point contacts
Zasad_doping , as well as in previous ITS studies involving chemical (oxygen)
doping Doping ; KatterwePRL2008 ; SecondOrder .

Figure 3 shows the electric doping of an initially strongly underdoped Bi-2212
mesa. A five-fold increase in critical current in the $I$-$V$ characteristics
(a), and an increase of $T_{\mathrm{c}}$ of about 15 K (c) can be seen. The
tunneling conductance $\mathrm{d}I/\mathrm{d}V(V)$ curves are shown in panel
(b). It is seen that the sum-gap peak shifts to slightly lower voltages and
becomes sharper with doping. The hump voltage, attributed to the $c$-axis
pseudogap SecondOrder , rapidly decreases with increasing conductance,
pointing towards an abrupt opening of the pseudogap at the critical doping
point, consistent with chemical doping studies Doping ; Balakirev ; Bernhard .

Note that in all studied cases the electric doping has lead to a significant
increase of the critical current, while the superconducting gap was
decreasing. The anticorrelation between $I_{\text{c}}R_{\text{n}}$ and
$\Delta$ in underdoped Bi-2212 has been reported before Doping  and was
attributed to progressively more incoherent $c$-axis transport in combination
with the $d$-wave symmetry of the superconducting order parameter.

###  IV.3 Short-pulse doping

Figure 5: (Color online) Short pulse resistive switching in a strongly
underdoped Bi-2212 mesa. (a) $I$-$V$ characteristics in the initial high-
resistive state (HRS) and the doped low-resistive state (LRS) at $T=2$ K.
Panels (b) and (c) demonstrate the resistive switching sequence
$R_{\text{ac}}$ versus time at $T=100$ K. The switching was made by a train of
short positive and negative pulses, shown in the bottom of panel (b). It is
seen that the HRS is stable, but the LRS is initially relaxing and then
saturates at a resistance below the HRS. Several hundred reproducible
resistive switching events can be achieved without visible degradation (c).

So far we were discussing a gradual electric doping of Bi-2212 mesas at the
time scale of an hour, as shown in Fig. 1. Such a long time doping allows a
very strong variation of the doping state, but often leads to inhomogeneous
doping within the mesa height, as shown in Fig. 2. Koval et.al. KovalMemres
demonstrated that a short-pulse doping strategy leads to highly reversible and
reproducible doping, similar to resistive switching in point contacts Tulina .
This is probably related to the lack of significant electromigration during
the short pulse, which may eventually lead to an irreversible destruction of
the crystal structure Moeckly .

Figure 4 represents the ITS characteristics at $T\sim 30\,\mathrm{K}$ for a
small $\sim 2\times 2\,\mathrm{\muup m^{2}}$ Bi(Y)-2212 mesa with a small
amount of junctions and small self-heating Heating_PRL2005 . In the initial
high-resistive state the mesa is slightly underdoped with $T_{\mathrm{c}}\sim
91\,\mathrm{K}$. The mesa was switched to a low-resistive state by a short
voltage pulse $V\gtrsim 2\,\mathrm{V}$ of about a milliseconds width. The
periodicity of QP branches in Fig. 4 (b) demonstrates that after switching
into the LRS the mesa remains highly uniform. From comparison of $I$-$V$
curves in panels (a) and (b) it is seen that the decrease of resistance by
$\sim 1/3$ is accompanied by an almost four-fold increase of $I_{\mathrm{c}}$.
The $T_{\mathrm{c}}$ increased to $\sim 93\,\mathrm{K}$ being an indication
that the mesa became near optimally doped.

Fig. 4 (c) represents the tunneling conductance $\mathrm{d}I/\mathrm{d}V$ in
the high-resistive (underdoped) and the low-resistive (optimally doped)
states. The following main changes in ITS spectra are seen: the
superconducting sum-gap peak voltage decreased in the doped low resistive
state and the shapes of spectra are changed. The relative sum-gap peak height
$\mathrm{d}I/\mathrm{d}V(V_{\mathrm{p}})R_{\mathrm{n}}$ is increased by about
50% in the low-resistive state. The high-resistive state exhibits a peak-dip-
hump structure, which is less obvious in the low-resistive state. All this is
similar to the slow-doping case, Figs. 2 and 3, and consistent with the change
of doping from the slightly underdoped to near optimally-doped state Doping ,
in accordance with other spectroscopic studies involving chemical (oxygen)
doping ARPES_doping ; STS_inhomo ; Zasad_doping .

As already discussed above, the increase of sharpness of the sum-gap kink in
the $I$-$V$ curve and the peak in the $\mathrm{d}I/\mathrm{d}V$ curve in the
low-resistive state, reported in Figs. 2 and 4, and the change of the shape of
the peak-dip-hump feature in Fig. 4 (c) can not be attributed to self-heating,
because the actual power dissipation at the peak is decreasing in the low-
resistive state. For example, the corresponding powers at the peaks in Fig. 4
(c) are $P=0.21$ and 0.17 mW, respectively. This observation supports the
conclusion of Ref. Comment2  that the appearance and the shape of the peak-
dip-hump structure in the ITS characteristics of small mesas is determined
primarily by the doping level.

A detailed analysis of short pulsed resistive switching has been performed on
strongly underdoped Bi-2212 mesas. For negative pulses with
$50\,\mathrm{\muup{}s}$ length and compliance voltages up to
$-2.5\,\mathrm{V}$ there was no resistive switching. At $-3.0\,\mathrm{V}$, a
reduction of the quasiparticle resistance to 99.3% is observed (LRS). Nine
subsequent pulses reduce the resistance further to 98.4%. A single positive
pulse switches the resistance back to the initial HRS, while more positive
pulses do not result in additional changes. A single negative pulse with a
higher compliance voltage of $-3.5\,\mathrm{V}$ instead reduces the resistance
of the LRS to about 93.4% of the HRS, and the corresponding positive pulses
switches it back. Pulses with lower compliance voltage lead to a partial
switching but additional pulses with the same compliance do not lead to a
significant change. Fig. 5 (a) shows $I$-$V$ curves of another UD Bi-2212 mesa
at $T=2\,\mathrm{K}$ in the HRS and the LRS obtained with a $3.5\,\mathrm{V}$
pulse of $50\,\mathrm{\muup{}s}$ width. The general difference between
$I$-$V$s after the pulsed doping in Figs. 4 (a) and 5 (a) is the same as for
the slow doping in Figs. 2 (a) and 3 (a).

Fig. 5 (b) demonstrates a reproducible resistive switching between HRS and LRS
at elevated $T=100\,\mathrm{K}$ for another UD Bi-2212 mesa. The switching was
performed using a similar positive and negative pulse sequence with $\pm
3.5\,\mathrm{V}$, a pulse width of $100\,\mathrm{\muup{}s}$ and an interval
between pulses of 3 s, shown in the bottom panel of Fig. 5 (b). Panels (b) and
(c) show the corresponding time sequence of the measured zero-bias resistance
of the mesa. It is seen that negative voltage pulses lead to switching into
the LRS while a subsequent positive pulse switches the mesa back into the HRS.
The corresponding resistance change rates are of the order of $\mathrm{d}\ln
R/\mathrm{d}t\approx\mp 1000\,\mathrm{s^{-1}}$, which are somewhat higher than
for the slow doping shown in Fig. 1 (d), but not inconsistent with that data,
taking into account that the compliance voltage is also significantly higher.
It is seen that the HRS is stable and shows no visible relaxation at
$T=100\,\mathrm{K}$, while the LRS is initially relaxing with the
characteristic time $\tau_{\mathrm{LRS}}=(0.77\pm 0.30)\,\mathrm{s}$ and then
saturates before reaching the HRS, as shown in panel (b). At
$T=280\,\mathrm{K}$, the behavior is similar with a time constant of
$\tau_{\mathrm{LRS}}=(0.75\pm 0.34)\,\mathrm{s}$.

##  V Discussion: mechanisms of persistent electric doping

In Sec. II we briefly reviewed known mechanisms of persistent physical doping
of cuprates. It is likely that some of them are playing a role in persistent
electric doping, studied here. Indeed, the phenomenon is clearly related to
the persistent electric-field effect Moeckly , observed in
$\mathrm{YBa_{2}Cu_{3}O_{6+x}}$, which in turn is clearly related to
persistent photo-doping observed for various cuprates Kudinov ; Schuller .

To identify possible mechanisms, we first summarize characteristic features of
the persistent electric doping:

i) Observed different voltage dependencies (odd and even) indicate that
several distinct mechanisms are involved.

ii) The doping rate shows a threshold-like behavior as a function of bias
voltage (see Fig. 1). The threshold voltage depends on $T$ in a thermal-
activation manner, i.e. decreases with increasing $T$. Remarkably, the
threshold voltage is weakly dependent on the number of junctions in the mesa,
consistent with previous reports KovalCurrentInj ; KovalMemres , and is
comparable to that for a single point contact Tulina . This suggests that the
phenomenon is connected to some characteristic energy $\sim 1\,\mathrm{eV}$,
rather than directly to the electric field. Indeed, for a given voltage, the
latter should scale inversely proportional to the number of junctions in the
mesa, i.e., would not be universal for different mesas (that is why we
hesitate to refer to the phenomenon as a persistent electric field effect, and
rather call it persistent electric doping).

iii) However, the role of the electric field should not be underestimated.
Indeed, the displacement field is given by $D=V\varepsilon_{\mathrm{r}}/Nt$,
where $\varepsilon_{\mathrm{r}}$ is the dielectric constant and $t$ is the
thickness of the insulating barrier between CuO bi-layers. The ratio
$t/\varepsilon_{\mathrm{r}}\simeq 0.1\,\mathrm{nm}$ was estimated from an
analysis of Fiske (geometrical resonance) step voltages Superluminal .
Therefore, the displacement field in intrinsic junctions is

$D\simeq\frac{V}{N}\times 10^{8}\,\mathrm{cm^{-1}}.$ (1)

For $V\simeq 2\,\mathrm{V}$ and $N\simeq 10$, which corresponds to the case of
Fig. 4, one would get a very large value $D\simeq 2\times
10^{7}\,\mathrm{V/cm}$, which is certainly capable of seriously polarizing and
displacing ions in complex oxides Moeckly ; Tulina ; Memristor1 .

iv) The phenomenon is not associated with a net change of the oxygen content,
which may only decrease within the cryostat. To the contrary, mesas can be
repeatedly and reversibly doped and undoped, as shown in Figs. 1 (b) and 5
(b).

###  V.1 Charge transfer and electrostatic charge trapping

If an injected electron has a high enough energy, it may join one of the ions,
leading to a change of the redox state and the effective doping. It was
suggested that such a “charge transfer” mechanism is involved both in photo-
doping of cuprates Kudinov  and in the resistive switching phenomenon in other
complex oxides Memristor1 .

Alternatively, the electron may be trapped (localized) in dielectric parts,
leading to electrostatic charging of the sample. The Bi-2212 compound has a
layered structure with metallic CuO planes sandwiched between polar insulator
BiO layers Polariton . In this case the electrostatic charging will take place
in BiO layers, which may affect the doping state of the neighboring CuO planes
via the electrostatic field-effect. The electrostatic charging of means by
current injection takes place uniformly within the whole structure. Therefore,
unlike the conventional electric field effect Ahn  and the electrostatic field
effect at the interface between a superconductor and a ferroelectric material
SuperFerro ; ElStatic , the current injection may lead to a persistent bulk
electrostatic field-effect doping of Bi-2212. Koval et. al. KovalCurrentInj
emphasized the similarity of the phenomenon with the floating-gate effect
utilized in Flash memory devices.

Both types of charging effects have common similarities:

i) The charge transfer requires a certain energy ($\sim$eV for
$\mathrm{YBa_{2}Cu_{3}O_{6+x}}$), rather than electric field.

ii) The sign of the current and the direction of the electric field does not
matter. Therefore, such doping should be even with respect to the voltage
sign, consistent with the observations in Ref. KovalCurrentInj .

Therefore, we attribute even in voltage persistent electric doping to charge
transfer and/or charge trapping mechanisms.

###  V.2 Oxygen reorientation and reordering

It is well established that the doping state of cuprates depends not only on
the amount of off-stoichiometric oxygen, but also on the relative orientation
of the oxygen bonds Chains . Therefore, the doping state can be changed by
oxygen reordering. Since the required energy is large $\sim$$\mathrm{eV}$,
compared to thermal energies, oxygen reordering is a slow process and does not
take place spontaneously at low enough $T$. Oxygen reordering is considered as
one of the main mechanisms of the persistent photo Schuller  and electric
field Moeckly  doping.

In the case of persistent electric-field doping, the oxygen reordering is
steered by the polarization. Therefore, the direction of doping should depend
on the direction of the electric field, i.e., should be odd with respect to
the bias voltage. We clearly see such a contribution in our experiment, see
Fig. 1 (b). Note that the odd-in-voltage contribution was reported in the
point-contact case Tulina , but not reported in previous related works made on
zig-zag type Bi-2212 microstructures KovalCurrentInj ; KovalMemres . The
geometry of the latter samples is symmetric with respect to the electric field
direction (changing the sign of the electric field is equivalent to flipping
their sample up-side down). This is not the case in point contacts and mesa
structures, studied here, for which the fields down (into the crystal) and up
(into the top electrode) are not equivalent. Therefore, the difference may
partly be due to the difference in sample geometry, or to the observed sample-
dependence of the relative strength of odd and even in voltage doping
contributions, as shown in Figs. 1 (d) and (e).

Thus, we attribute the odd in voltage persistent electric doping mechanism to
field-induced oxygen reorientation/reordering.

###  V.3 Irreversible processes: electromigration, filament and arc formation

An increase of the bias voltage above $\sim 2.5-3.0\,\mathrm{V}$ leads to a
gradual increase of the current and irreversible change of the mesa
properties. A similar phenomenon was observed in
$\mathrm{YBa_{2}Cu_{3}O_{6+x}}$ thin films and attributed to electromigration
and field-induced diffusion of oxygen, which is even in bias voltage. The
increase of conductance is probably due to a dielectric breakthrough in the
insulating BiO layers, which leads to a pin-hole and filament formation. Thus,
we attribute the slow and irreversible drift of the mesa characteristics at
large bias voltages to electromigration in the mesas. This destructive process
is, however, distinctly different from the reversible and reproducible
electric doping effect, reported above.

After deterioration by electromigration, the mesa characteristics become
similar to resistive switching characteristics for point contacts on top of
oxygen-depleted, Bi-2212 surfaces Tulina . At even higher bias the resistance
becomes very high (infinite). But an inspection in a microscope shows that the
mesa itself remains intact. There is no physical evaporation of material or a
crater at the place of the mesa, as in the case of a violent electric
discharge. Instead there are clear indications of an arc formation at one of
the sharp corners of the mesa, which probably leads to delamination of the
structure and mechanical disattachment of the mesa from the base crystal.

###  V.4 Mechanisms of energy accumulation

The most puzzling property of the reported persistent electric doping is that
the required bias voltage is weakly dependent on the number of junctions in
the mesa KovalCurrentInj ; KovalMemres . This is clearly seen from the
presented data, for which the threshold voltage is always $\sim 2\,\mathrm{V}$
for $N=9$ in Fig. 4, $N=11$ in Fig. 1 (f) and $N\sim 56$ in Fig. 1 (d), which
is also similar to that for a single point-contact $\sim 1.5$ V Tulina . It
is, therefore, clear that electric doping requires a certain electron energy,
rather than electric field. However, for stacked tunnel junctions, the energy
acquired by the injected electron in every tunneling event is proportional to
the voltage drop across the junction, $\delta E\simeq eV/N$, and is
significantly smaller than the threshold energy. The main question is,
therefore, how the electrons accumulate a sufficiently large energy, required
for doping.

In Ref. KovalMemres  it was suggested that an electron can accumulate energy
upon sequential tunneling through several junctions without relaxation.
However, in this case only electrons in the last junction will have enough
energy. This would result in a strongly inhomogeneous doping in different
junctions. Moreover, the probability of sequential tunneling without
relaxation is small, because of the small ratio of relaxation time
($\tau\sim\mathrm{ps}$) Cascade  to the tunnel time, $t_{\mathrm{tun}}$,
$t_{\mathrm{tun}}/\tau\gg 1$. The probability of sequential tunneling through
$N$ junctions is decreasing rapidly $\propto(\tau/t_{\mathrm{tun}})^{N}$ with
increasing $N$. This should lead to a dramatic increase of the doping time
with increasing $N$. Indeed, suppose that it takes $t_{N}\sim
10\,\mathrm{min}$ for $N$ junctions at $I=1\,\mathrm{mA}$ to dope the mesa.
This will involve $N_{e}=It_{N}/e$ tunneling events in each junction ($e$ is
the electron charge). Since the probability of sequential tunneling through
$2N$ junctions is decreasing quadratically, it would require $N_{e}^{2}$
tunneling events per junction, which will take $t_{2N}=N_{e}t_{N}=3.75\times
10^{19}\,\mathrm{min}$. However, such a dramatic increase of the doping time
with increasing mesa height is inconsistent with experiments.

For the sequential tunneling scenario to be relevant, the ratio
$t_{\mathrm{tun}}/\tau$ should rapidly drop with increasing electron energy
and become of the order of unity at $E\sim 1\,\mathrm{eV}$. This may be caused
by resonant tunneling, which increases the tunneling rate of quasiparticles
with a certain energy, and/or by a drastic slowing down of the high energy
quasiparticle relaxation, which may be caused by a rapid decrease of the
Eliashberg’s electron-boson spectral function and a gap in the corresponding
bosonic density of states at high energies Cascade . Such a scenario is
interesting to investigate because it may give an information about the
bosonic spectrum, involved in Cooper pairing, and thus provide a clue about
the electron-boson coupling mechanism, responsible for high-$T_{\mathrm{c}}$
superconductivity in cuprates Cascade .

We also want to propose an alternative mechanism for the energy accumulation
of electrons: the formation of electric-field domains in the natural atomic
superlattice formed by the mesa. Electric field domains are well studied in
semiconducting superlattices Superlattice . They appear in weakly coupled
superlattices close to the resonant tunneling condition. The corresponding
non-linearity leads to an instability and a multiple-valued current-voltage
characteristics. As a result, the electric field distribution in the
superlattice becomes nonuniform and is concentrated in one or several
junctions.

A possibility of a formation of electric field domains in Bi-2212 mesas is not
just a hypothesis. In fact, the multiple-branch $I$-$V$ of Bi-2212 mesas due
to one-by-one switching of intrinsic junctions from the superconducting to the
resistive state, shown in Fig. 4 (b), is due to a formation of electric field
domains in individual tunnel junctions. The formation of electric field
domains in Bi-2212 mesas at high bias would explain many of the features of
the studied persistent electric doping. In this case electrons in the domain
may gain an energy close to eV without sequential tunneling through the whole
mesa. Furthermore, since domains are typically dynamic and propagate through
the whole superlattice Superlattice , this would also explain the uniformity
of doping in the whole mesa, and not just in the outermost junction.

##  VI Conclusions

We have studied the effect of persistent electric doping on intrinsic
tunneling characteristics of small Bi-2212 mesa structures. It was shown that
the application of a sufficiently large voltage to the mesas leads to a
controllable and reversible physical doping of the mesas, without a
modification of their chemical composition. This allows the analysis of bulk
electronic spectra in Bi-2212 in a wide doping range on one and the same mesa.
This physical doping has the same effect as chemical (oxygen) doping on the
intrinsic tunneling characteristics of Bi-2212: the $c$-axis resistivity
decreases, the critical current increases and the energy gap is decreasing
together with $T_{\text{c}}$ with over-doping. The anticorrelation between
$I_{\text{c}}R_{\text{n}}$ and $\Delta$ indicates that the $c$-axis transport
becomes progressively more incoherent at the underdoped side of the phase
diagram. An analysis of the doping variation of the intrinsic tunneling
characteristics of the same mesa provides a clue about its genuine shape: with
subsequent doping, the sum-gap peak in the tunneling conductance becomes
sharper and the pseudogap hump rapidly decreases with doping, suggesting the
presence of a critical doping point, in agreement with previous chemical
doping studies Doping .

By analyzing the bias- and time dependence we could identify different
mechanisms involved in the persistent electric doping: i) The even-in-voltage
process via charge transfer and/or charge trapping. ii) The odd-in-voltage
process via oxygen reordering. Those are distinct from the irreversible
electromigration and oxygen electrodiffusion, observed at higher bias.

We confirm the previous report KovalMemres  that the threshold voltage for the
electric doping is weakly dependent of the number of junctions in the mesas
and is similar to that for a single surface point contact Tulina . This
indicates that it is the energy of injected electrons, rather than electric
field, that determines the phenomenon. We suggest that the required energy
accumulation by tunnel electrons may be due to a formation of electric field
domains in the natural atomic superlattice formed in the Bi-2212 single
crystal.

Acknowledgments We are grateful to the Swedish Research Council and the SU-
Core Facility in Nanotechnology for financial and technical support.

## References

  * (1) J. L. Tallon, and J. W. Loram, Physica C 349, 53 (2001).
  * (2) V. M. Krasnov, _Phys. Rev. B_ 79, 214510 (2009).
  * (3) V. M. Krasnov, H. Motzkau, T. Golod, A. Rydh, S. O. Katterwe, and A. B. Kulakov, _Phys. Rev. B_ 84, 054516 (2011).
  * (4) T. Kondo, Y. Hamaya, A. D. Palczewski, T. Takeuchi, J. S. Wen, Z. J. Xu, G. Gu, J. Schmalian, and A. Kaminski, _Nature Phys._ 7, 21 (2011).
  * (5) R.-H. He, M. Hashimoto, H. Karapetyan, J. D. Koralek, J. P. Hinton, J. P. Testaud, V. Nathan, Y. Yoshida, H. Yao, K. Tanaka, W. Meevasana, R. G. Moore, D. H. Lu, S.-K. Mo, M. Ishikado, H. Eisaki, Z. Hussain, T. P. Devereaux, S. A. Kivelson, J. Orenstein, A. Kapitulnik, and Z. X. Shen, _Science_ 331, 1579 (2011).
  * (6) V. M. Krasnov, _Phys. Rev. B_ 65, 140504(R) (2002).
  * (7) F. F. Balakirev, J. B. Betts, A. Migliori, I. Tsukada, Y. Ando, and G. S. Boebinger, _Phys. Rev. Lett._ 102, 017004 (2009).
  * (8) T. Helm, M. V. Kartsovnik, M. Bartkowiak, N. Bittner, M. Lambacher, A. Erb, J. Wosnitza and R. Gross Phys. Rev. Lett. 103, 157002 (2009).
  * (9) D. LeBoeuf, N. Doiron-Leyraud, B. Vignolle, M. Sutherland, B. J. Ramshaw, J. Levallois, R. Daou, F. Laliberte, O. Cyr-Choiniere, J. Chang, Y. J. Jo, L. Balicas, R. Liang, D. A. Bonn, W. N. Hardy, C. Proust, and L. Taillefer Phys. Rev. B 83, 054506 (2011).
  * (10) V. M. Krasnov, A. Yurgens, D. Winkler, P. Delsing and T.Claeson, Phys. Rev. Lett. 84, 5860 (2000).
  * (11) W. S. Lee, I. M. Vishik, K. Tanaka, D. H. Lu, T. Sasagawa, N. Nagaosa, T. P. Devereaux, Z. Hussain and Z. X. Shen, Nature 450, 81 (2007).
  * (12) L. Yu, D. Munzar, A. V. Boris, P. Yordanov, J. Chaloupka, Th. Wolf, C. T. Lin, B. Keimer and Ch. Bernhard, Phys. Rev. Lett. 100, 177004 (2008).
  * (13) R. Daou, J. Chang, D. LeBoeuf, O. Cyr-Choiniere, F. Laliberte, N. Doiron-Leyraud, B. J. Ramshaw, R. Liang, D. A. Bonn, W. N. Hardy, and L. Taillefer, Nature 463, 519 (2010).
  * (14) K. Fujita, T. Noda, K. M. Kojima, H. Eisaki, and S. Uchida, Phys. Rev. Lett. 95, 097006 (2005).
  * (15) N. Chandrasekhar, O. T. Valls, and A. M. Goldman, Phys. Rev. B 49, 6220 (1994).
  * (16) M. Käll, M. Osada, M. Kakihana, L. Börjesson, T. Frello, J. Madsen, N. H. Andersen, R. Liang, P. Dosanjh, and W. N. Hardy, Phys. Rev. B 57, 14072(R) (1998).
  * (17) V. I. Kudinov, I. L. Chaplygin, A. I. Kirilyuk, N. M. Kreines, R. Laiho, E. Lähderanta, and C. Ayache, Phys. Rev. B 47, 9017 (1993).
  * (18) A. Gilabert and A. Hoffmann and A. G. Medici and I. K. Schuller, _J. Supercond. Novel Magn._ 13, 1 (2000).
  * (19) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, _Science_ 331, 189 (2011).
  * (20) C. H. Ahn, A. Bhattacharya, M. Di Ventra, J. N. Eckstein, C. D. Frisbie, M. E. Gershenson, A. M. Goldman, I. H. Inoue, J. Mannhart, A. J. Millis, A. F. Morpurgo, D. Natelson, and J.-M. Triscone, _Rev. Mod. Phys._ 78, 1185 (2006).
  * (21) B. H. Moeckly, D. K. Lathrop, and R. A. Buhrman, Phys. Rev. B 47, 400 (1993); S. H. Huerth, H. D. Hallen and B. Moeckly, Phys. Rev. B 67, 180506(R) (2003).
  * (22) P. Konsin, and B. Sorkin, Phys. Rev. B 58, 5795 (1998).
  * (23) N. A. Tulina, G. A. Emelchenco, and A. B. Kulakov, Physics Lett. A 204, 74 (1995); N. A. Tulina, A. M. Ionov, and A. N. Chaika, Physica C 366, 23 (2001); A. Plecenik, M. Grajcar, P. Seidel, S. Takacs, A. Matthes, M. Zuzcak, and S. Benacka, Physica C 301, 234 (1998).
  * (24) Y. Koval, X. Jin, C. Bergmann, Y. Simsek, L. Özyüzer, P. Müller, H. Wang, G. Behr, and B. Büchner, _Appl. Phys. Lett._ 96, 082507 (2010).
  * (25) Y. Koval, F. Chowdhury, X. Jin, Y. Simsek, F. Lichtenberg, R. Pentcheva, P. Müller, _Phys. Status Solidi A_ 208, 284 (2011).
  * (26) R. Waser, and M. Aono, _Nature Materials_ 6, 833 (2007); R. Waser, R. Dittmann, M. Salinga, and M. Wuttig, _Sol. St. Electr._ 54, 830 (2010).
  * (27) J. J. Yang, M. D. Pickett, X. Li, D. A.A. Ohlberg, D. R. Stewart, and R. S. Williams, _Nature Nanotechn._ 3, 429 (2008).
  * (28) R. Kleiner, P. Müller, _Phys. Rev. B_ 49, 1327 (1994).
  * (29) K. Anagawa, Y. Yamada, T. Watanabe, and M. Suzuki, Phys. Rev. B 67, 214513 (2003).
  * (30) S. O. Katterwe and V. M. Krasnov, _Phys. Rev. B_ 80, 020502(R) (2009).
  * (31) S. O. Katterwe, A. Rydh, H. Motzkau, A. B. Kulakov, and V. M. Krasnov, _Phys. Rev. B_ 82, 024517 (2010).
  * (32) S. O. Katterwe, A Rydh, and V. M Krasnov _Phys. Rev. Lett._ 101, 087003 (2008).
  * (33) N. Kristoffel, and P. Rubin, _Physica C_ 418, 49 (2005).
  * (34) B. Y. Shapiro, _Sol. St. Commun._ 56, 149 (1985).
  * (35) R. Aidam, D. Fuchs, and R. Schneider, _Physica C_ 328, 21 (1999).
  * (36) A. Walkenhorst, C. Doughty, X. X. Xi, Q. Li, C. J. Lobb, S. N. Mao, and T. Venkatesan, _Phys. Rev. Lett._ 69, 2709 (1992); X. X. Xi, _J. Supercond._ 7, 137 (1994).
  * (37) S. O. Katterwe, H. Motzkau, A. Rydh, and V. M. Krasnov, Phys. Rev. B 83, 100510(R) (2011).
  * (38) M. Salluzzo, G. Ghiringhelli, J. C. Cezar, N. B. Brookes, G. M. DeLuca, F. Fracassi, and R. Vaglio, Phys. Rev. Lett.  100, 056810 (2008).
  * (39) V. M. Krasnov, T. Bauch, and P. Delsing, _Phys. Rev. B_ 72, 012512 (2005).
  * (40) V. M. Krasnov, M. Sandberg, and I. Zogaj, Phys. Rev. Lett. 94, 077003 (2005).
  * (41) V. M. Krasnov, _Phys. Rev. B_ 84, 136501 (2011).
  * (42) H. Ding, J. R. Engelbrecht, Z. Wang, J. C. Campuzano, S.-C. Wang, H.-B. Yang, R. Rogan, T. Takahashi, K. Kadowaki, and D. G. Hinks, _Phys. Rev. Lett._ 87, 227001 (2001).
  * (43) K. M. Lang, V. Madhavan, J. E. Hoffman, E. W. Hudson, H. Eisaki, S. Uchida, and J. C. Davis, _Nature_ 415, 412 (2002).
  * (44) L. Ozyuzer, J. F. Zasadzinski, K. E. Gray, C. Kendziora, and N. Miyakawa, _Europhys. Lett._ 58, 589 (2002).
  * (45) V. M. Krasnov, _Phys. Rev. Lett._ 97, 257003 (2006); _ibid._ , 103, 227002 (2009).
  * (46) L. L. Bonilla, and H. T. Grahn, Rep. Prog. Phys. 68, 577 (2005).