# Higher-order Dirac Semimetal in a Photonic Crystal

Zihao Wang  These authors contributed equally to this work Division of Physics
and Applied Physics, School of Physical and Mathematical Sciences, Nanyang
Technological University, 21 Nanyang Link, Singapore 637371, Singapore
Dongjue Liu  These authors contributed equally to this work School of
Electrical and Electronic Engineering, Nanyang Technological University, 50
Nanyang Avenue, Singapore 639798, Singapore  Hau Tian Teo  Division of Physics
and Applied Physics, School of Physical and Mathematical Sciences, Nanyang
Technological University, 21 Nanyang Link, Singapore 637371, Singapore  Qiang
Wang  Division of Physics and Applied Physics, School of Physical and
Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link,
Singapore 637371, Singapore  Haoran Xue<EMAIL_ADDRESS>Division of
Physics and Applied Physics, School of Physical and Mathematical Sciences,
Nanyang Technological University, 21 Nanyang Link, Singapore 637371, Singapore
Baile Zhang<EMAIL_ADDRESS>Division of Physics and Applied Physics,
School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, Singapore 637371, Singapore  Centre for
Disruptive Photonic Technologies, The Photonics Institute, Nanyang
Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore

###### Abstract

The recent discovery of higher-order topology has largely enriched the
classification of topological materials. Theoretical and experimental studies
have unveiled various higher-order topological insulators that exhibit
topologically protected corner or hinge states. More recently, higher-order
topology has been introduced to topological semimetals. Thus far, realistic
models and experimental verifications on higher-order topological semimetals
are still very limited. Here, we design and demonstrate a three-dimensional
photonic crystal that realizes a higher-order Dirac semimetal phase. Numerical
results on the band structure show that the designed three-dimensional
photonic crystal is able to host two four-fold Dirac points, the momentum-
space projections of which at an edge are connected by higher-order hinge
states. The higher-order topology can be characterised with the calculation of
the $\chi^{(6)}$ topological invariant at different values of $k_{z}$. An
experiment at microwave frequencies is also presented to measure the hinge
state dispersion. Our work demonstrates the physical realization of a higher-
order Dirac semimetal phase and paves the way to exploring higher-order
topological semimetals phases in three-dimensional photonic systems.

††preprint: APS/123-QED

Band topology has been a major topic in condensed matter physics since the
discovery of the quantum Hall effect [1], leading to discoveries of various
topological phases of matter [2, 3]. In recent years, the notion of band
topology has also been utilised to manipulate classical waves, such as in the
emerging fields of topological photonics [4, 5, 6] and topological acoustics
[7, 8]. In both electronic and classical topological systems, the bulk-
boundary correspondence dictates the connection between the bulk topology and
the topological boundary states. Conventionally, the dimension of topological
boundary states is one order lower than that of the bulk. For example, a two-
dimensional (2D) Chern insulator supports chiral edge states that propagate
along its one-dimensional (1D) edges.

Recently, it has been found that the dimension difference between the bulk and
the topological boundary states can be larger than one, leading to the
discovery of higher-order topological phases [9, 10, 11, 12, 13, 14]. For
example, a 2D second-order topological insulator (TI) does not host protected
1D edge states. Instead, it supports zero-dimensional (0D) topological corner
states localized at corners. Over the last four years, various higher-order
topological systems have been theoretically proposed and experimentally
demonstrated [11, 12, 13, 15, 16, 17, 14, 18, 19, 20, 21]. So far, most
studies have focused on higher-order TIs, i.e., systems with a bulk bandgap.

Topological semimetals (TSMs) are another class of topological phases [22].
Different from TIs, TSMs host gapless bulk band structures. One paradigmatic
example is the Weyl semimetal, which supports two-fold linear band crossings
called Weyl points in the bulk, and topological Fermi-arc states on the
surfaces connecting two Weyl points with opposite topological charges [23].
Very recently, studies have found that TSMs can also support higher-order
topological states. Theoretical proposals have predicted topological hinge
states in higher-order Weyl semimetals, Dirac semimetals and nodal lines
semimetals [23, 24, 25, 26, 27, 25, 28, 29, 30, 22, 31, 32, 33, 34]. A few
experiments have been performed on the platform of acoustic metamaterials to
characterize higher-order TSM phases [32, 33, 34].

In this work, we introduce the concept of higher-order TSMs into photonic
systems by presenting the design and realization of a higher-order Dirac
semimetal (HODSM) in a three-dimensional (3D) photonic crystal. Our model
consists of coupled layers of deformed photonic honeycomb lattices. We
identify the Dirac points (DPs), topology of bulk bands, and the existence of
hinge states through full-wave simulations and tight-binding calculations. In
addition, a microwave experiment on the dispersion of the hinge states is also
presented.

Figure 1:  Structure and bulk dispersion of the photonic crystal exhibiting a
higher-order Dirac semimetal phase. (a) Schematic of the unit cell of a
deformed honeycomb lattice with lattice constant $a$. The lattice consists of
grey rods sandwiched by metallic plates. The rods have height $h$, diameter
$d_{1}$ and are separated by a distance $c_{1}$. The two metallic plates on
the top and bottom are drilled with holes to couple the neighbouring layers.
The holes have the depth $d$, diameter $d_{2}$ and are separated by a distance
$c_{2}$. (b) The first Brillouin zone of the crystal. (c) Simulated bulk
dispersion of the photonic crystal. The coloured region from $\Gamma$ to A
indicates the topological transition when $k_{z}$ goes through the Dirac point
at $k_{z}$ = $2\pi/5a_{z}$. (d) Schematic of phase transition in the entire
Brillouin zone at different values of $k_{z}$. The upper and lower parts host
the higher-order topological insulator (HOTI) phase, while the middle part
possesses the normal insulator (NI) phase. The two red dots denote the two
Dirac points.

The designing principles of our photonic HODSM are as follows. First, let us
consider a 2D photonic honeycomb lattice with each lattice site occupied by a
dielectric rod and the surrounding space being air. Then, we choose a larger
cell that contains six rods (see Fig. 1(a)) and moves the rods towards/away
from the centre of the cell [35]. Previous studies have shown that both
deformations will open a bandgap at $\Gamma$ point. When the configuration
changes from one deformation to the other, the bandgap closes and reopens,
inducing the phase transition between a normal insulator (NI) phase and a
higher-order topological insulator (HOTI) phase [36, 37, 38, 39, 40, 41, 42].
As we adjust the distance between rods, the coupling strength between them
changes as well. For the original honeycomb lattice, the distance of two
adjacent rods within a unit cell and that between two unit cells are equal,
thus leading to the equal intracell and intercell coupling. In our photonic
lattice, we adopt the shrunken lattice in which the distance of rods in a unit
cell is smaller, so the ratio between the intercell and intracell coupling is
tuned to be less than one. As a result, a complete bandgap is opened at
$\Gamma$ point in the Brillouin zone (BZ) and the 2D lattice is in the NI
phase. The final step is to extend this 2D lattice into a 3D one by properly
engineering the interlayer couplings. The key of our construction is, as can
be seen from Fig. 1(a), a rod not only couples to rods at the same in-plane
position in adjacent layers but also couples to rods at neighbouring in-plane
positions within the same unit cell in adjacent layers. Thus, a
$k_{z}$-dependent term is introduced as a tuning parameter in the Bloch
Hamiltonian as an additional term, changing the intercell and intracell
coupling strengths as a function of $k_{z}$. In this model, the situation with
equal intercell and intracell coupling strengths is a transition point and it
corresponds to a DP in the momentum space. The NI phase and the HOTI phase are
separated by the DP, forming a TSM. The tight-binding model analysis is
performed in the Supplementary Material I [43] to show the details.

Fig. 1(a) shows the designed unit cell, which consists of six dielectric
cylindrical rods. These rods are made of dielectric material that has the
relative permittivity $\epsilon_{r}$ = 8.5 and relative permeability $\mu_{r}$
= 1. Two metallic plates drilled with air holes are placed on the top and at
the bottom, acting as the perfect electric conductor (PEC) and the air holes
enable wave coupling between adjacent layers. The leftover part of the unit
cell is filled with air with $\epsilon_{r}$ = 1 and $\mu_{r}$ = 1. The in-
plane lattice constant is $a$ = 25 mm. In the vertical direction, the lattice
constant is $a_{z}$ = 10.2 mm and the height of the rods is $h$ = 9 mm. The
metallic plates have thickness $d$ = 0.6 mm and thus the air holes between
layers are 1.2 mm thick. The rods have the diametre $d_{1}$ = 6 mm and the air
holes $d_{2}$ = 6.4 mm. The distance between the adjacent rods centres is
$c_{1}$ = 7.88 mm, and the distance between adjacent air holes is $c_{2}$ =
7.94 mm. Numerically obtained bulk dispersion is shown in Fig. 1(c). The
reciprocal lattice is depicted in Fig. 1(b) with the first BZ identified using
blue lines. In Fig. 1(c), the band dispersion of the left region
($\Gamma$-K-M-$\Gamma$) lies on the $k_{z}$ = 0 plane, holding a topologically
trivial phase, while the right region (A-H-L-A) lies on the $k_{z}$ =
$\pi/a_{z}$ plane, holding a topologically nontrivial phase. The middle region
($\Gamma$-A) shows the transition of $k_{z}$ from $k_{z}$ = 0 to $\pi/a_{z}$.
The DP appears at $k_{z}$ = $2\pi/5a_{z}$, serving as the phase transition
point dividing the left and right halves of the bulk dispersion diagram. Fig.
1(d) depicts the distribution of the DPs along $k_{z}$ in the BZ. The upper
and lower parts of the BZ hold the topologically nontrivial HOTI phase, while
the middle part between two DPs has the topologically trivial NI phase. Both
in Fig. 1(c) and (d), the same colours are used to identify the same phases
along the phase transition process.

Figure 2:  The $\chi^{(6)}$ topological invariants and the band structures for
the model in Fig. 1(a) with different $k_{z}$ values. (a) The $\chi^{(6)}$ =
$([\mathbb{M}],[\mathbb{K}])$ topological invariants calculated from the
tight-binding model in the first Brillouin zone showing the phase transition.
(b-d) The bulk band structure for that in Fig. 1(a) with $k_{z}=$0,
$2\pi/5a_{z}$, and $\pi/a_{z}$ separately. The green bands in (b) and (d) show
the bandgaps when $k_{z}$ deviates from the Dirac point.

To verify the topological nature of the sample, when $k_{z}$ roams over the
first BZ, the $\chi^{(6)}$ topological invariants are calculated to show the
phase transition process and the HOTI phase between two DPs. The calculation
is done on a series of 2D planes with different $k_{z}$. As proposed in [18],
the phase transition can be represented by checking a pair of topological
invariants $\chi^{(6)}$ = $([\mathbb{M}],[\mathbb{K}])$. Since the system
maintains rotational invariance at several high-symmetry momentum points, it
allows us to use the difference of rotation eigenvalues between three momentum
points, M, K and $\Gamma$, to characterise the topological nature of that
plane. The $[\mathbb{M}]$ and $[\mathbb{K}]$ are defined as
$[\mathbb{M}]=\\#\mathrm{M_{1}}-\\#\Gamma_{1}^{(2)}$,
$[\mathbb{K}]=\\#\mathrm{K_{1}}-\\#\Gamma_{1}^{(3)}$, and
$[\mathbb{M}],[\mathbb{K}]\in\mathbb{Z}$. Here, $\\#\mathrm{M_{1}}$ is the
number of bands below the gap in the spectrum which have the rotation
eigenvalue 1 at this rotation invariant momenta point M and similar for
$\\#\Gamma_{1}^{(\alpha)}$, $\alpha=2,3$, and $\\#\mathrm{K_{1}}$ (See
Supplementary Material II [43] for details). As shown in Fig. 2(a), when
$k_{z}$ varies, $[\mathbb{K}]$ will stick to 0 all the time, but
$[\mathbb{M}]$ experiences two times of jumps, which exactly corresponds to
the critical point when the intracell and intercell coupling strengths are the
same. In the HOTI phase, $[\mathbb{M}]$ holds the value 2 while in the NI
phase it jumps to 0. Figs. 2(b-d) show the band structure of the model in Fig.
1(a) when $k_{z}$ equals to 0, $2\pi/5a_{z}$, and $\pi/a_{z}$, respectively.
The bulk topology experiences a transition here along with the bandgap at
$\Gamma$ point closes and reopens, transforming from NI to HOTI as discussed
before.

To demonstrate the HODSM and observe the hinge states, we fabricate a 3D
sample as shown in Fig. 3(a). The sample is a stack of 19 layers in total. The
four long thin rods at the corners are for positioning and alignment. The
inset is the zoom-in view of deformed hexagonal clusters. In the experiment,
in order to observe the hinge state localisation, we have adopted PECs to
cover the two side surfaces around the obtuse-angle corner to form the hard
boundary. (see Fig. 3(b) and Supplemental Material III [43]). Fig. 3(b)
illustrates the experimental setup. When probing the hinge state, a source
(shown as the red star) is fixed near the hinge located at the hard boundary
at the half-height of the whole structure. Then we use a probe (shown as the
blue star) to scan over all the unit cells. After collecting the data, the
Fourier transform is performed and the result is shown in Fig. 3(d).

Fig. 3(c) is the simulation result for a structure with 20 $\times$ 20 unit
cells in the $xy$ plane and periodic boundary condition along $z$ direction.
The black, blue and red dots represent the bulk, surface and hinge states,
respectively. Before the phase transition point at $k_{z}$ = $2\pi/5a_{z}$,
all the states belong to the bulk states, indicating a NI phase. After the
phase transition point, the higher-order hinge states emerge along with the
surface states and the hinge states, as predicted by the $\chi^{(6)}$ index.
Fig. 3(d) shows the experimental results of the hinge states. We also redraw
the numerical hinge dispersion by the green curve. As can be seen, the
experiment result coincides with the simulation result well. Fig. 3(e) is the
field distribution of the hinge state measured from the experiment. Here we
show the result of eleven sample planes with an identical height difference of
20.4 mm between any two adjacent planes. The frequency is chosen to be 6.4
GHz, falling in the range of the hinge states dispersion in Fig. 3(c) and (d).
In Fig. 3(e), the colour and the height of the tiny columns on each layer
denote the amplitude of the field at that site point. At each layer, the
amplitude of the strongest site is taken as the basis for the normalisation
and all the other sites are normalised according to it. As we can find in Fig.
3(e), all the eleven layers show the trend that the field localises on the
hinge on the hard boundary. Compared to the hinge localisation, the states in
the bulk away from the corner decreases rapidly and become invisible. The
experimental field distribution of the bulk distribution can be found in
Supplementary Material IV [43] and obvious contrast between the hinge and bulk
mode are shown.

Figure 3:  Experimental measurement on the higher-order hinge states. (a)
Photograph of the experimental sample. The inset shows the details of the
deformed lattice rods placed on each layer with the top plate removed. (b) The
sketch of the experimental setup. The red line denotes the source pin inserted
from the top air hole. The red star denotes the source dipole fixed at the
half-height of the sample next to hard boundaries. The blue line denotes the
position of the probe. The purple lines denote the side surfaces with perfect
electric conductors and the measurement is performed around the obtuse angle
consisting of these two surfaces. (c) Simulation result of a 20 $\times$ 20
lattice. The black, blue, and red dots denote the bulk, surface, and hinge
state. (d) The measured energy dispersion. The green curve represents the
hinge states redrawn from (c). (e) The measured field distribution showing the
hinge state localised on the hard boundary of the sample.

In our sample, it is a parallelogram in the $x$-$y$ plane, so it hosts two
obtuse-angle corners and two acute-angle ones. Both the obtuse-angle corners
are able to hold the hinge state localisation there. The hinge mode localised
on the corner with hard boundaries is discussed above. On the other obtuse
angle, however, since the source is fixed far away from it, it is difficult to
observe the obvious localisation on that hinge. More than that, the EM wave
propagating to that corner also couple to the atmosphere outside the sample.
Both of these two reasons prevent us from capturing obvious hinge localisation
there. So in the text, only the hinge state localised on the hard boundary is
shown.

In summary, we have designed and demonstrated a HODSM using the photonic
crystal. The design is based on coupling 2D deformed honeycomb lattices with
proper interlayer coupling along $z$ direction. The Dirac points emerge at the
phase transition point between a HOTI phase and a NI phase when $k_{z}$
varies. The hinge states that connect two Dirac points projected on a hinge
are directly observed with field scanning measurement at microwave
frequencies. Our work provides a concrete example of HODSM in photonic
systems. Moreover, the design scheme presented here can also be applied to
construct other types of photonic topological semimetals with higher-order
band topology, such as Weyl semimetals and nodal line semimetals.

We thank Y. Long for discussions. This work is supported by National Research
Foundation Singapore Competitive Research Program No. NRF-CRP23-2019-0007.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Lu  _et al._ [2014] L. Lu, J. D. Joannopoulos, and M. Soljačić, Nat. Photonics 8, 821 (2014).
  * Khanikaev and Shvets [2017] A. B. Khanikaev and G. Shvets, Nat. Photonics 11, 763 (2017).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, _et al._ , Rev. Mod. Phys. 91, 015006 (2019).
  * Yang  _et al._ [2015] Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, Phys. Rev. Lett. 114, 114301 (2015).
  * Ma  _et al._ [2019] G. Ma, M. Xiao, and C. T. Chan, Nat. Rev. Phys. 1, 281 (2019).
  * Benalcazar  _et al._ [2017a] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ [2017b] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b).
  * Langbehn  _et al._ [2017] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Song  _et al._ [2017] Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Schindler  _et al._ [2018] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv 4, eaat0346 (2018).
  * Ezawa [2018] M. Ezawa, Phys. Rev. Lett. 120, 026801 (2018).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nat. Phys. 14, 925 (2018).
  * Serra-Garcia  _et al._ [2018] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Peterson  _et al._ [2018] C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * Noh  _et al._ [2018] J. Noh, W. A. Benalcazar, S. Huang, M. J. Collins, K. P. Chen, T. L. Hughes, and M. C. Rechtsman, Nat. Photonics 12, 408 (2018).
  * Xue  _et al._ [2019] H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Nat. Mater. 18, 108 (2019).
  * Ni  _et al._ [2019] X. Ni, M. Weiner, A. Alu, and A. B. Khanikaev, Nat. Mater. 18, 113 (2019).
  * Zhang  _et al._ [2019] X. Zhang, H.-X. Wang, Z.-K. Lin, Y. Tian, B. Xie, M.-H. Lu, Y.-F. Chen, and J.-H. Jiang, Nat. Phys. 15, 582 (2019).
  * Armitage  _et al._ [2018] N. Armitage, E. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
  * Liu  _et al._ [2014] Z. Liu, B. Zhou, Y. Zhang, Z. Wang, H. Weng, D. Prabhakaran, S.-K. Mo, Z. Shen, Z. Fang, X. Dai, _et al._ , Science 343, 864 (2014).
  * Lu  _et al._ [2015] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, and M. Soljačić, Science 349, 622 (2015).
  * Xu  _et al._ [2015] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, _et al._ , Science 349, 613 (2015).
  * Bernevig [2015] B. A. Bernevig, Nat. Phys. 11, 698 (2015).
  * Burkov [2016] A. Burkov, Nat. Mater. 15, 1145 (2016).
  * Fang  _et al._ [2016] C. Fang, L. Lu, J. Liu, and L. Fu, Nat. Phys. 12, 936 (2016).
  * Lin and Hughes [2018] M. Lin and T. L. Hughes, Phys. Rev. B 98, 241103 (2018).
  * Wieder  _et al._ [2020] B. J. Wieder, Z. Wang, J. Cano, X. Dai, L. M. Schoop, B. Bradlyn, and B. A. Bernevig, Nat. Commun. 11, 627 (2020).
  * Luo  _et al._ [2021] L. Luo, H.-X. Wang, Z.-K. Lin, B. Jiang, Y. Wu, F. Li, and J.-H. Jiang, Nat. Mater. 20, 794 (2021).
  * Wei  _et al._ [2021] Q. Wei, X. Zhang, W. Deng, J. Lu, X. Huang, M. Yan, G. Chen, Z. Liu, and S. Jia, Nat. Mater. 20, 812 (2021).
  * Qiu  _et al._ [2020] H. Qiu, M. Xiao, F. Zhang, and C. Qiu, Higher-order dirac sonic crystals (2020), arXiv:2012.00251 [cond-mat.other] .
  * Wu and Hu [2015] L.-H. Wu and X. Hu, Phys. Rev. Lett. 114, 223901 (2015).
  * Xie  _et al._ [2018] B.-Y. Xie, H.-F. Wang, H.-X. Wang, X.-Y. Zhu, J.-H. Jiang, M.-H. Lu, and Y.-F. Chen, Phys. Rev. B 98, 205147 (2018).
  * Zangeneh-Nejad and Fleury [2019] F. Zangeneh-Nejad and R. Fleury, Phys. Rev. Lett. 123, 053902 (2019).
  * Xie  _et al._ [2019] B.-Y. Xie, G.-X. Su, H.-F. Wang, H. Su, X.-P. Shen, P. Zhan, M.-H. Lu, Z.-L. Wang, and Y.-F. Chen, Phys. Rev. Lett. 122, 233903 (2019).
  * Liu  _et al._ [2019] F. Liu, H.-Y. Deng, and K. Wakabayashi, Phys. Rev. Lett. 122, 086804 (2019).
  * Benalcazar  _et al._ [2019] W. A. Benalcazar, T. Li, and T. L. Hughes, Phys. Rev. B 99, 245151 (2019).
  * Li  _et al._ [2020] M. Li, D. Zhirihin, M. Gorlach, X. Ni, D. Filonov, A. Slobozhanyuk, A. Alù, and A. B. Khanikaev, Nat. Photonics 14, 89 (2020).
  * Xie  _et al._ [2020] B. Xie, G. Su, H.-F. Wang, F. Liu, L. Hu, S.-Y. Yu, P. Zhan, M.-H. Lu, Z. Wang, and Y.-F. Chen, Nat. Commun. 11, 3768 (2020).
  * [43] see Supplemental Material at (URL will be inserted by the publisher) for details on the tight-binding model, topological rotation invariant, experiment sample and experiment results for hinge and bulk modes. .


Also at ]Dept. of Mechanical and Materials Engineering, Florida International
University, Miami, FL33174, United States

Also at ]Nanoelectronics Integrated Systems Center (NISC), Nile University,
Cairo, Egypt and ]Dept. of Electrical and Computer Engineering, University of
Calgary, Calgary, Canada

# Information Encoding/Decoding using the Memory Effect in Fractional-order
Capacitive Devices

Anis Allagui<EMAIL_ADDRESS>Dept. of Sustainable and Renewable Energy
Engineering, University of Sharjah, PO Box 27272, Sharjah, United Arab
Emirates  Center for Advanced Materials Research, Research Institute of
Sciences and Engineering, University of Sharjah, PO Box 27272, Sharjah, United
Arab Emirates  [  Ahmed S. Elwakil  Dept. of Electrical and Computer
Engineering, University of Sharjah, PO Box 27272, Sharjah, United Arab
Emirates  [  [

###### Abstract

In this study, we show that the discharge voltage pattern of a fractional-
order supercapacitor from the same initial steady-state voltage into a
constant resistor is dependent on the past charging voltage profile. The
charging voltage was designed to follow a power-law function, i.e.
$v_{c}(t)=V_{cc}\left({t}/{t_{ss}}\right)^{p}\;(0<t\leqslant t_{ss})$, in
which $t_{ss}$ (charging time duration between zero voltage to the terminal
voltage $V_{cc}$) and $p$ ($0<p<1$) act as two variable parameters. We used
this history-dependence of the dynamic behavior of the device to uniquely
retrieve information pre-coded in the charging waveform pattern. Furthermore,
we provide an analytical model based on fractional calculus that explains
phenomenologically the information storage mechanism. The use of this
intrinsic material memory effect may lead to new types of methods for
information storage and retrieval.

Keywords : Fractional-order capacitors; Memory; Fractional calculus; Constant-
phase element

## Introduction

Many systems and processes in nature exhibit fractional-order behavior, such
as vestibulo-ocular reflex Anastasio (1994), neuronal activity Teka, Marinov,
and Santamaria (2014), ion channel gating Goychuk and Hänggi (2004),
viscoelasticity Mainardi and Spada (2011), and disordered semiconductors Scher
and Montroll (1975), which are models known to capture the existing
short/long-term memory effects in these systems Sabatier and Farges (2014);
Teka, Upadhyay, and Mondal (2017); Yuan, Fu, and Liu (2014); Ventosa-
Santaulària, Heres, and Martínez-Hernández (2014). Supercapacitors are well-
established types of electrochemical capacitive energy storage devices, but
are also known to exhibit non-ideal, time-fractional-order electric behavior
when charged by a power supply and discharged into a load Kant and Singh
(2015); Allagui  _et al._ (2016a); Allagui, Elwakil, and Freeborn (2017);
Fouda  _et al._ (2016); Allagui  _et al._ (2018a, 2017). Their electric
impedance can be modeled as a series resistance ($R_{s}$) with a fractional-
order capacitor (also known as constant-phase element, CPE) over a certain
frequency range; the reduced impedance of such a model is
$Z^{*}(u)=1+{1}/{(\mathrm{j}u)^{\alpha}}$ with
$u=\omega(R_{s}C_{\alpha})^{1/\alpha}$, $C_{\alpha}$ is a fractional-order
capacitance in units of F sα-1, $\alpha$ ($0<\alpha<1$) is a dimensionless
fractional exponent, and $\omega$ is the applied angular frequency in units of
s-1 (see Nyquist plot of a commercially-available NEC/TOKIN supercapacitor in
Fig. 1). In the time-domain, the current-voltage relationship of an
$R_{s}$-CPE-equivalent supercapacitor is expressed by the fractional-order
differential equation Allagui  _et al._ (2016b, 2018b):

$R_{s}C_{\alpha}\frac{d^{\alpha}V_{C_{\alpha}}(t)}{dt^{\alpha}}+V_{C_{\alpha}}(t)=v_{c}(t)$
(1)

where $V_{C_{\alpha}}$ is the voltage across the CPE part ($V_{C_{\alpha}}=0$
for $t\leq 0$), $C_{\alpha}{d^{\alpha}V_{C_{\alpha}}(t)}/{dt^{\alpha}}$ is the
current flowing through the CPE, and $v_{c}(t)$ is the applied charging
voltage. The differentiation of non-integer order in equation 1 is defined as
Podlubny (1999):

$\frac{d^{-\alpha}{V}(t)}{dt^{-\alpha}}=\frac{1}{\Gamma(\alpha)}\int\limits_{0}^{t}{V}(\tau)(t-\tau)^{\alpha-1}d\tau,$
(2)

which can be viewed as a convolution with a hyperbolic function of frequency,
and therefore contains a memory that progressively increases as the fractional
order $\alpha$ decreases Westerlund and Ekstam (1994); Westerlund (1991); Du,
Wang, and Hu (2013). The physical interpretation of such fractional-order
electric behavior of supercapacitors is still under debate, nonetheless, it
has been widely attributed to the surface chemistry and morphological
structure of the electrodes, which are usually composed of high-surface area
and porous materials separated by an ionic conductor Córdoba-Torres, Mesquita,
and Nogueira (2015); Alexander, Tribollet, and Orazem (2015, 2016); Alexander
_et al._ (2017).

Figure 1: Nyquist plane representation of open-circuit spectral impedance of a
NEC/TOKIN supercapacitor (part #FGR0H105ZF, 5.5 V, 1 F). Complex nonlinear
least-squares fitting to $Z(s)=R_{s}+1/C_{\alpha}s^{\alpha}$
($s=\mathrm{j}\omega$) shows two straight line regions giving the values
$(R_{s};C_{\alpha};\alpha)=(6.306\;\Omega;\;0.138\;\text{F\,s}^{\alpha-1};\;0.49)$
from 10 Hz to 20 mHz, and
$(16.87\;\Omega;\;0.524\;\text{F\,s}^{\alpha-1};\;0.88)$ from 20 mHz to 5 mHz.
20 mHz is a critical frequency separating the near-ideal capacitive behavior
from the Warburg region

Despite the fact that fractional-order models involve hereditary effects, very
few robust experimental proofs of this memory have been reported and
investigated Westerlund (1991); Westerlund and Ekstam (1994); Uchaikin  _et
al._ (2016); Uchaikin, Sibatov, and Uchaikin (2009); Allagui, Zhang, and
Elwakil (2018); Allagui  _et al._ (2020). This Letter illustrates the memory
effect in supercapacitors by sequentially encoding information into the
charging pattern of the device, and then uniquely retrieving this code from
the discharge pattern. The purpose here is to investigate this effect
experimentally and mathematically using fractional-order models, which means
that we are leaving the physical interpretation and potential applications of
this effect to another study. The results are obtained on state-of-the-art,
commercially-available devices which offer excellent stability, and relatively
high-voltage ratings which facilitate repeatability and precise measurements.

## Methods

We applied different charging voltage waveforms to a NEC/TOKIN supercapacitor
(part #FGR0H105ZF, 5.5 V, 1 F111This device consists of six 0.917 V-aqueous
electrolytic cells stacked in series; each cell is composed of two symmetric
activated carbon + dilute sulfuric acid electrodes separated by a porous
organic film) following the power-law function:

$v_{c}(t)=V_{cc}\left(\frac{t}{t_{ss}}\right)^{p}\quad(0<t\leqslant t_{ss})$
(3)

where $p$ is an exponent taking values between 0 (step voltage) and 1 (linear
ramp), and $t_{ss}$ is the rise time from 0 V to the steady-state value
$V_{cc}$ ($V_{cc}=5.5\;\text{V}$ for this device), after which the charging
voltage is turned off. The rise time $t_{ss}$ is pre-defined so that the
device will operate either in the capacitive tail from 20 mHz to 5 mHz, or in
the Warburg region from 10 Hz to 20 mHz (see the two quasi-linear regions in
Fig. 1). Prior to each applied charging voltage waveform, the supercapacitor
was fully discharged into a constant 100 $\Omega$ resistor ($R_{p}$) until its
voltage was equal to 1 mV. Depending on the values of $p$ and $t_{ss}$, the
charge waveforms used for information storage and the discharge waveforms used
for information retrieval are given by the letter/number codes in Table 1. All
charging/discharging experiments of the supercapacitor were programmed and
executed sequentially on a Bio-logic VSP-300 electrochemical station using the
EC-Lab control software. The time step for collecting data in all measurements
was set to 10 ms.

Table 1: Table of codes illustrating encoding and decoding waveforms symbols in the form of different values of $p$ and $t_{ss}$ in $v_{c}(t)=V_{cc}({t}/{t_{ss}})^{p}$. We will be using the superscripts $"c"$ and $"d"$ for charge and discharge, respectively. For example $A^{c}_{10}$ refers to the code $A_{10}$ during charge while $A^{d}_{10}$ refers to the same code during discharge. $p$ \$t_{ss}$ | 550 s | 275 s | 110 s | 55 s | 27 s
---|---|---|---|---|---
1.0 | $A_{10}$ | $B_{10}$ | $C_{10}$ | $D_{10}$ | $E_{10}$
0.7 | $A_{07}$ | $B_{07}$ | $C_{07}$ | $D_{07}$ | $E_{07}$
0.4 | $A_{04}$ | $B_{04}$ | $C_{04}$ | $D_{04}$ | $E_{04}$
0.2 | $A_{02}$ | $B_{02}$ | $C_{02}$ | $D_{02}$ | $E_{02}$
0.1 | $A_{01}$ | $B_{01}$ | $C_{01}$ | $D_{01}$ | $E_{01}$

## Results

Figure 2: Charging sequence (left) using a power supply with a voltage
$v_{c}(t)=V_{cc}\left({t}/{t_{ss}}\right)^{p}$, and discharging sequence
(right) into a constant 100 $\Omega$ resistor of an NEC/TOKIN supercapacitor.
(a)-(b) depict the charge/discharge voltage patterns for different values of
$p$ (1.0, 0.7, 0.4, 0.2 and 0.1) and $t_{ss}=550$ s (table 1), i.e.
$A^{c}_{10}$/$A^{d}_{10}$-$A^{c}_{07}$/$A^{d}_{07}$-$A^{c}_{04}$/$A^{d}_{04}$-$A^{c}_{02}$/$A^{d}_{02}$-$A^{c}_{01}$/$A^{d}_{01}$.
(c)-(d) and (e)-(f) show the repeatability of the process for four consecutive
cycles (4x, superimposed on top of each other) of charge/discharge with
$p=1.0$ and $p=0.1$; in (c)-(d) waveforms
$A^{c}_{10}$/$A^{d}_{01}$-$A^{c}_{10}$/$A^{d}_{01}$-$A^{c}_{10}$/$A^{d}_{01}$-$A^{c}_{10}$/$A^{d}_{01}$,
and in (e)-(f) waveforms
$E^{c}_{10}$/$E^{d}_{01}$-$E^{c}_{10}$/$E^{d}_{01}$-$E^{c}_{10}$/$E^{d}_{01}$-$E^{c}_{10}$/$E^{d}_{01}$.
In (g)-(h), we applied the sequence of charge/discharge codes
[$B^{c}_{10}$/$B^{d}_{10}$\- $C^{c}_{05}$/$C^{d}_{05}$\-
$D^{c}_{01}$/$D^{d}_{01}$]–[$B^{c}_{10}$/$B^{d}_{10}$\-
$C^{c}_{05}$/$C^{d}_{05}$\-
$D^{c}_{01}$/$D^{d}_{01}$]–[$C^{c}_{05}$/$C^{d}_{05}$\-
$D^{c}_{01}$/$D^{d}_{01}$\-
$B^{c}_{10}$/$B^{d}_{10}$]–[$C^{c}_{05}$/$C^{d}_{05}$\-
$D^{c}_{01}$/$D^{d}_{01}$\- $B^{c}_{10}$/$B^{d}_{10}$] using three values of
$t_{ss}$ and three values of $p$ (Table 1) demonstrating the possibility of
two-dimensional encoding as well as the repeatability of the process (4x)

The experimental results are shown in Fig. 2. In the first column of the
figure, we show several combinations of the charging voltage waveforms applied
to the supercapacitor by varying the values of $p$ and $t_{ss}$ (see Table 1).
The time scale is shifted by $-t_{ss}$ to represent past events. The second
column of the figure shows the first 60 s of the resulting voltage discharge
waveforms into the same 100 $\Omega$ resistor. Here the potentiostat acts as a
constant resistance by controlling the current to maintain the ratio
voltage/current constant. Fig. 2(b) shows five different voltage discharge
profiles ($A_{10}^{d}$, $A_{07}^{d}$, $A_{04}^{d}$, $A_{02}^{d}$,
$A_{01}^{d}$) after the device was charged with five different voltage
waveforms ($A_{10}^{c}$, $A_{07}^{c}$, $A_{04}^{c}$, $A_{02}^{c}$,
$A_{01}^{c}$) (Fig. 2(a)). The value of $t_{ss}$ is made long enough so that
the device will operate within its low frequency capacitive tail (Fig. 1). All
discharging waveforms show first a quick voltage drop of _ca._ 0.5 V from the
initial voltage $V_{cc}=5.5\;\text{V}$ into the internal series resistance of
the device ($R_{s}$), followed by a non-Debye inverse power law profile (i.e.
$v_{d}(t)\propto t^{-\alpha}$ Jonscher (1977); Westerlund and Ekstam (1994);
Westerlund (1991); Jonscher (1999)). In Fig. 2(d), we show the discharge
profiles of the device for four consecutive cycles (4x) alternating between
code $A_{10}^{c}$ and code $A_{01}^{c}$ (Fig. 2(c)). The corresponding
discharge waveform codes $A_{10}^{d}$ and $A_{01}^{d}$ are perfectly
superimposed on each other demonstrating the repeatability of the
charge/discharge process and the good stability of the device between one
sequence to another. The results in Figs 2(a)-(d) show that the supercapacitor
voltage discharge profile into the same constant resistance depends on the
exponent $p$ in the charging voltage profile, and thus on its prehistory. This
would not be the case for an ideal capacitor for which the exponential
decaying voltage $v_{d}(t)=V_{cc}\exp\left(-{t}/{R_{p}C}\right)$ depends on
the initial voltage $V_{cc}$, but not on how this voltage has been reached..
In other words, information can be encoded in the exponent $p$ being for
example $1.0$ or $0.1$ (Figs. 2(c) and 2(d)) which can be used to encode a
binary data sequence, or in multi-level logic using multiple values of $p$
(Figs. 2(a) and 2(b)).

The effect of the parameter $t_{ss}$ is examined in Figs. 2(e) and 2(f), which
depict the voltage profiles that were obtained in a similar way to the results
shown in Figs. 2(c) and 2(d), respectively, but with a faster charging rate
($t_{ss}=27\,\text{s}$). This value of $t_{ss}$ corresponds to about 37 mHz
frequency which belongs to the Warburg region and not the capacitive tail
anymore (see Fig. 1). It is clear from Fig. 2(f) that the superposition of the
code plots ($E_{10}^{d}$ and $E_{01}^{d}$) obtained in response to the
charging voltage codes $E_{10}^{c}$ and $E_{01}^{c}$ is again impeccable.
Additionally, the difference between the discharge waveforms is more
pronounced than when $t_{ss}$ was set to $550\,\text{s}$ (Fig. 2(d)) at which
point the supercapacitor behaved as a near-ideal capacitor. Although the
values of $R_{s}$ and $C_{\alpha}$ are not same for these two cases, the
correlation between the discharging and the charging voltage waveforms is
stronger as $\alpha$ decreases (equations 1 and 2) which makes $t_{ss}$
another possible information coding dimension.

As for the decoding of the discharge pattern corresponding to a specific
charge pattern, this can be carried out by a simple adaptive thresholding
method either for a fixed time or fixed voltage. For example, as shown in Fig.
2(f), the thresholding can be performed at 5 seconds from the beginning of the
discharge giving the two voltage values of 2.441 and 3.238 V for $p=1.0$ and
$p=0.1$ (i.e. a difference of 0.797 V), respectively. Alternatively, it can be
performed at a fixed voltage value of 2.7 V for example, leading to decoded
time intervals of 3.109 and 14.31 s for $p=1.0$ and $p=0.1$, respectively. It
is important to note that decoding from the discharging waveform should be
carried out way before the device is fully discharged, as shown in decimal
algorithm scale in the insets in Figs. 2(d) and 2(f) depicting the quick
convergence of the voltage-time profiles. The reason is that fractional-order
behavior is space- and time-dependent; i.e. only in the transient time does
the order of the state-space (represented by the values of $\alpha$) manifest
itself leading to the memory effect Allagui, Zhang, and Elwakil (2018);
Allagui  _et al._ (2020). In the steady-state time, fractional-order behavior
becomes space-order independent.

Another set of experimental results are shown in Figs. 2(g) and 2(h), in which
charge/discharge of the device followed the sequence
[$B_{10}$-$C_{05}$-$D_{01}$]–[$B_{10}$-$C_{05}$-$D_{01}$]–[$C_{05}$-$D_{01}$-$B_{10}$]–[$C_{05}$-$D_{01}$-$B_{10}$].
This indicates the possibility of two-dimensional memory encoding using $p$
and $t_{ss}$.

## Discussion

Analytically, the memory effect can be explained as follows. The voltage
discharge response of the supercapacitor, modeled as an $R_{s}$-CPE equivalent
circuit, into a parallel resistance $R_{p}$ is given by Freeborn, Maundy, and
Elwakil (2013):

$v_{d}(t)=v_{d}(0)\frac{R_{p}}{R_{p}+R_{s}}\text{E}_{\alpha,1}\left(-\frac{t^{\alpha}}{\tau_{p}^{\alpha}+\tau_{s}^{\alpha}}\right)$
(4)

where $\tau_{p}=\left(R_{p}C_{\alpha}\right)^{1/\alpha}$,
$\tau_{s}=\left(R_{s}C_{\alpha}\right)^{1/\alpha}$,
$\text{E}_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}z^{k}/\Gamma(\alpha\,k+\beta)$
with ($\alpha>0$, $\beta>0$) is the two-parameter Mittag-Leffler function, and
the initial voltage $v_{d}(0)$ is given by
$v_{d}(0)=V_{cc}-R_{s}i_{c}(t_{ss})$ where $i_{c}(t)$ is the charging current
obtained as $dq_{c}(t)/dt$ (see equation 7). For $R_{p}\gg R_{s}$, equation 4
can be rewritten in a dimensionless form as follows:

$\tilde{v}_{d}{(\tilde{t}_{p})}\simeq
R_{(p,m,\alpha)}\text{E}_{\alpha,1}\left(-\tilde{t}_{p}^{\,\alpha}\right)$ (5)

where $\tilde{v}_{d}{(\tilde{t}_{p})}=v_{d}(t)/V_{cc}$,
$\tilde{t}_{p}=(t/\tau_{p})$ and $m=(t_{ss}/\tau_{s})$. The Mittag-Leffler
function in equation 5 depends on the CPE parameters $\alpha$ and $C_{\alpha}$
(and on $R_{p}$), whereas the parametric function $R_{(p,m,\alpha)}$ is
dependent also on how the device has been charged via the selection of one or
both of the applied waveform parameters $p$ and $t_{ss}$. For an ideal
capacitor, i.e. with $\alpha=1$ and $R_{s}=0$, we verify that equation 5
simplifies to $\tilde{v}_{d}{(\tilde{t}_{p})}=\exp(-\tilde{t}_{p})$ which does
not depend on its charge prehistory, as expected.

Figure 3:  (a) Supercapacitor charging using the voltage waveforms
$v_{c}(t)=V_{cc}\left({t}/{t_{ss}}\right)^{p}$ with $t_{ss}=27\;\text{s}$ and
$p=1.0$ and $0.1$ (i.e. $E_{10}^{c}$ and $E_{01}^{c}$) and the resulting time-
charge profiles. (b) Illustrates the corresponding discharging voltage
waveforms (i.e. $E_{10}^{d}$ and $E_{01}^{d}$). Experimental measurements are
plotted in solid lines and simulations are represented by dashed lines. (c)
Flowchart representing the process and equations of the electric variables
during charge (i.e. applied voltage, accumulated charge, and corresponding
current), and voltage during discharge into a constant resistor (the device is
considered. to be a fractional-order supercapacitor of impedance
$Z=R_{s}+1/C_{\alpha}s^{\alpha}$ and being charged with a voltage function
$v_{c}(t)=V_{cc}\left({t}/{t_{ss}}\right)^{p}$)

Now, to show the memory relationship, we derive the expression for the
electrical charge stored in the device subjected to the charging waveform
$v_{c}(t)=V_{cc}\left({t}/{t_{ss}}\right)^{p}$ as follows Allagui, Zhang, and
Elwakil (2018):

$q_{c}(t)=\left(\frac{C_{\alpha}V_{cc}\Gamma(p+1)}{t_{ss}^{p}\tau_{s}^{\alpha}}\right)t^{p+1}\text{E}_{\alpha,p+2}\left[-\left(\frac{t}{\tau_{s}}\right)^{\alpha}\right]\vspace{.1cm}$
(6)

from which the current is found to be:

$i_{c}(t)=\left(\frac{C_{\alpha}V_{cc}\Gamma(p+1)}{t_{ss}^{p}\tau_{s}^{\alpha}}\right)t^{p}\text{E}_{\alpha,p+1}\left[-\left(\frac{t}{\tau_{s}}\right)^{\alpha}\right]\vspace{.1cm}$
(7)

[Note that Equation 6 is obtained by equating $i(t)=dq/dt$ with
$C_{\alpha}{d^{\alpha}V_{C_{\alpha}}(t)}/{dt^{\alpha}}$ where
$V_{C_{\alpha}}(t)$ is the solution of the fractional-order differential
equation
$V_{cc}\left({t}/{t_{ss}}\right)^{p}=V_{C_{\alpha}}+R_{s}C_{\alpha}{d^{\alpha}V_{C_{\alpha}}(t)}/{dt^{\alpha}}$
obtained using the inverse Laplace transform identity
$\mathscr{L}^{-1}\left(\frac{k!s^{\alpha-\beta}}{(s^{\alpha}+\lambda)^{k+1}}\right)=t^{\alpha\,k+\beta-1}E_{\alpha,\beta}^{(k)}[-\lambda
t^{\alpha}]$ with ($k=0$, $\beta=p+1$, $\lambda=\tau_{s}^{-\alpha}$) Podlubny
(1999)]. At $t=t_{ss}$, the steady-state charge $q_{c}(t_{ss})$ is a function
of both parameters of the charge waveform ($p$ and $t_{ss}$) and those of the
supercapacitor. [Equation 6 simplifies to $q_{c}(t)=CV_{cc}t/t_{ss}$ for an
ideal capacitance ($C_{\alpha}=C$ and $\alpha=1$) using the identity
$E_{1,3}(z)=(\exp(z)-z-1)/z^{2}$ which applies when $p=1$. If $p=0$ (step
voltage), the charge $q_{c}(t)=CV_{cc}$, otherwise $q_{c}(t)$ is a function of
$p$]. This is in line with our recent findings in which we highlighted that
charging a supercapacitor with a voltage input results in a device and
waveform-dependent accumulated electric charge Allagui, Elwakil, and Freeborn
(2017); Allagui, Zhang, and Elwakil (2018); Allagui  _et al._ (2020). In a
dimensionless form, equation 6 looks like this:

$\tilde{q}_{c}(\tilde{t}_{s})=S_{(p,m,\alpha)}\tilde{t}_{s}^{\,p+1}\text{E}_{\alpha,p+2}\left(-\tilde{t}_{s}^{\,\alpha}\right)$
(8)

where $\tilde{t}_{s}=t/\tau_{s}$. In Fig. 3(a), we show the measured and
simulated (using equation 6) charge function $q_{c}(t)$ for a fixed value of
$t_{ss}$ equal to 27 s, and steady-value voltage of 5.5 V. Two values of the
parameter $p$ (i.e. 1.0 and 0.1) are selected. Equation 6 is in excellent
agreement with the experiment using the values ($R_{s},C_{\alpha},\alpha$)
equal to ($10\,\Omega,35\,\text{mF\,s}^{\alpha-1},0.48$) for $p=1.0$, and
equal to ($10\,\Omega,36\,\text{mF\,s}^{\alpha-1},0.49$) for $p=0.1$ obtained
using least-squares fitting. Experimentally, the accumulated charge at time
$t_{ss}$ is found to be 0.68 and 0.99 mA h-1 for $p=1.0$ and $0.1$,
respectively. Given that all parameters in equation 6 are fixed apart from
$p$, it is evident that the information encoded in the value of $p$ is stored
as $q_{c}(t)$. In Fig. 3(b), we show the resulting discharge voltage along
with simulation using equation 4. The best fit is found using the values
($R_{s},C_{\alpha},\alpha$) equal to
($10\,\Omega,803\,\text{mF\,s}^{\alpha-1},0.96$) for $p=1.0$, and equal to
($14.24\,\Omega,855\,\text{mF\,s}^{\alpha-1},0.95$) for $p=0.1$, which can be
improved if a sliding mode fitting is adopted instead. Note that the increase
of the values of $\alpha$ and $C_{\alpha}$ towards those of a near-ideal
capacitance is the result of the slow discharge rate.

## Conclusion

We have shown that the voltage discharge of a supercapacitor that exhibits
fractional-order temporal dynamics depends uniquely on the way by which it was
charged. Figure 3(c) recapitulates the process from charging the device by an
external power supply, to the accumulation of charge and the creation of
current, and their relationship with the discharge voltage function which
depends on the voltage charging parameters $p$ and $t_{ss}$. In other words,
the supercapacitor remembers the pattern by which it was charged, and as a
result, discharges accordingly.

The experimental results reported here were obtained on a commercial, high-
capacitance supercapacitor which resulted in slow information storage and
retrieval. Higher read-write rates should be possible in principle by using
low-capacitance, fractional-order devices Agambayev  _et al._ (2017). However,
it is important in this case to ensure stability over a large number of
cycles, as well as to reduce the effect of any parasitic capacitance for
properly assessing the memory effect.

## References

## References

  * Anastasio (1994) T. J. Anastasio, “The fractional-order dynamics of brainstem vestibulo-oculomotor neurons,” Biol. Cybern. 72, 69–79 (1994).
  * Teka, Marinov, and Santamaria (2014) W. Teka, T. M. Marinov, and F. Santamaria, “Neuronal spike timing adaptation described with a fractional leaky integrate-and-fire model,” PLOS Comput. Biol. 10, e1003526 (2014).
  * Goychuk and Hänggi (2004) I. Goychuk and P. Hänggi, “Fractional diffusion modeling of ion channel gating,” Phys. Rev. E 70, 051915 (2004).
  * Mainardi and Spada (2011) F. Mainardi and G. Spada, “Creep, relaxation and viscosity properties for basic fractional models in rheology,” Eur. Phys. J. Spec. Top. 193, 133–160 (2011).
  * Scher and Montroll (1975) H. Scher and E. W. Montroll, “Anomalous transit-time dispersion in amorphous solids,” Phys. Rev. B 12, 2455 (1975).
  * Sabatier and Farges (2014) J. Sabatier and C. Farges, “Long memory models: a first solution to the infinite energy storage ability of linear time-invariant fractional models,” IFAC Proceedings Volumes 47, 2884–2890 (2014).
  * Teka, Upadhyay, and Mondal (2017) W. W. Teka, R. K. Upadhyay, and A. Mondal, “Fractional-order leaky integrate-and-fire model with long-term memory and power law dynamics,” Neural Networks 93, 110–125 (2017).
  * Yuan, Fu, and Liu (2014) N. Yuan, Z. Fu, and S. Liu, “Extracting climate memory using fractional integrated statistical model: A new perspective on climate prediction,” Sci. Rep. 4, 6577 (2014).
  * Ventosa-Santaulària, Heres, and Martínez-Hernández (2014) D. Ventosa-Santaulària, D. R. Heres, and L. C. Martínez-Hernández, “Long-memory and the sea level-temperature relationship: a fractional cointegration approach,” PloS one 9, e113439 (2014).
  * Kant and Singh (2015) R. Kant and M. B. Singh, “Generalization of Randles-Ershler admittance for an arbitrary topography electrode: application to random finite fractal roughness,” Electrochim. Acta 163, 310–322 (2015).
  * Allagui  _et al._ (2016a) A. Allagui, A. S. Elwakil, B. J. Maundy, and T. J. Freeborn, “Spectral capacitance of series and parallel combinations of supercapacitors,” ChemElectroChem 3, 1429–1436 (2016a).
  * Allagui, Elwakil, and Freeborn (2017) A. Allagui, A. S. Elwakil, and T. J. Freeborn, “Supercapacitor reciprocity and response to linear current and voltage ramps,” Electrochim. Acta 258, 1081–1085 (2017).
  * Fouda  _et al._ (2016) M. E. Fouda, A. S. Elwakil, A. G. Radwan, and A. Allagui, “Power and energy analysis of fractional-order electrical energy storage devices,” Energy 111, 785–792 (2016).
  * Allagui  _et al._ (2018a) A. Allagui, A. S. Elwakil, M. Fouda, and A. G. Radwan, “Capacitive behavior and stored energy in supercapacitors at power line frequencies,” J. Power Sources 390, 142–147 (2018a).
  * Allagui  _et al._ (2017) A. Allagui, Z. Said, M. A. Abdelkareem, A. S. Elwakil, M. Yang, and H. Alawadhi, “DC and AC performance of graphite films supercapacitors prepared by contact glow discharge electrolysis,” J. Electrochem. Soc. 164, A2539–A2546 (2017).
  * Allagui  _et al._ (2016b) A. Allagui, T. J. Freeborn, A. S. Elwakil, and B. J. Maundy, “Reevaluation of performance of electric double-layer capacitors from constant-current charge/discharge and cyclic voltammetry,” Sci. Rep. 6, 38568 (2016b).
  * Allagui  _et al._ (2018b) A. Allagui, T. J. Freeborn, A. S. Elwakil, M. E. Fouda, B. J. Maundy, A. G. Radwanh, Z. Said, and M. A. Abdelkareem, “Review of fractional-order electrical characterization of supercapacitors,” J. Power Sources 400, 457–467 (2018b).
  * Podlubny (1999) I. Podlubny,  _Fractional differential equations_ (Academic Press, 1999).
  * Westerlund and Ekstam (1994) S. Westerlund and L. Ekstam, “Capacitor theory,” IEEE Trans. Dielectr. Electr. Insul. 1, 826–839 (1994).
  * Westerlund (1991) S. Westerlund, “Dead matter has memory!” Phys. Scr. 43, 174 (1991).
  * Du, Wang, and Hu (2013) M. Du, Z. Wang, and H. Hu, “Measuring memory with the order of fractional derivative,” Sci. Rep. 3, 3431 (2013).
  * Córdoba-Torres, Mesquita, and Nogueira (2015) P. Córdoba-Torres, T. J. Mesquita, and R. P. Nogueira, “Relationship between the origin of constant-phase element behavior in electrochemical impedance spectroscopy and electrode surface structure,” J. Phys. Chem. C 119, 4136–4147 (2015).
  * Alexander, Tribollet, and Orazem (2015) C. L. Alexander, B. Tribollet, and M. E. Orazem, “Contribution of surface distributions to constant-phase-element (cpe) behavior: 1. influence of roughness,” Electrochim. Acta 173, 416–424 (2015).
  * Alexander, Tribollet, and Orazem (2016) C. L. Alexander, B. Tribollet, and M. E. Orazem, “Contribution of surface distributions to constant-phase-element (cpe) behavior: 2. capacitance,” Electrochim. Acta 188, 566–573 (2016).
  * Alexander  _et al._ (2017) C. L. Alexander, B. Tribollet, V. Vivier, and M. E. Orazem, “Contribution of surface distributions to constant-phase-element (cpe) behavior: 3. adsorbed intermediates,” Electrochim. Acta 251, 99–108 (2017).
  * Uchaikin  _et al._ (2016) V. Uchaikin, A. Ambrozevich, R. Sibatov, S. Ambrozevich, and E. Morozova, “Memory and nonlinear transport effects in charging–discharging of a supercapacitor,” Tech. Phys. 61, 250–259 (2016).
  * Uchaikin, Sibatov, and Uchaikin (2009) V. Uchaikin, R. Sibatov, and D. Uchaikin, “Memory regeneration phenomenon in dielectrics: the fractional derivative approach,” Phys. Scr. 2009, 014002 (2009).
  * Allagui, Zhang, and Elwakil (2018) A. Allagui, D. Zhang, and A. S. Elwakil, “Short-term memory in electric double-layer capacitors,” Appl. Phys. Lett. 113, 253901–5 (2018).
  * Allagui  _et al._ (2020) A. Allagui, D. Zhang, I. Khakpour, A. S. Elwakil, and C. Wang, “Quantification of memory in fractional-order capacitors,” J. Phys. D 53 (2020).
  * Note (1) This device consists of six 0.917+.1667emV-aqueous electrolytic cells stacked in series; each cell is composed of two symmetric activated carbon + dilute sulfuric acid electrodes separated by a porous organic film.
  * Jonscher (1977) A. K. Jonscher, “The ’universal’ dielectric response,” Nature 267, 673–679 (1977).
  * Jonscher (1999) A. K. Jonscher, “Dielectric relaxation in solids,” J. Phys. D: Appl. Phys. 32, R57 (1999).
  * Freeborn, Maundy, and Elwakil (2013) T. J. Freeborn, B. Maundy, and A. S. Elwakil, “Measurement of supercapacitor fractional-order model parameters from voltage-excited step response.” IEEE J. Emerg. Sel. Topics Circuits Syst. 3, 367–376 (2013).
  * Agambayev  _et al._ (2017) A. Agambayev, S. P. Patole, M. Farhat, A. Elwakil, H. Bagci, and K. N. Salama, “Ferroelectric fractional-order capacitors,” ChemElectroChem 4, 2807–2813 (2017).