# Phase sensitivity approaching quantum Cramér-Rao bound in a modified SU(1,1)
interferometer

Jian-Dong Zhang<EMAIL_ADDRESS>School of Mathematics and
Physics, Jiangsu University of Technology, Changzhou 213001, China  Chenglong
You  Quantum Photonics Laboratory, Department of Physics & Astronomy,
Louisiana State University, Baton Rouge, LA 70803, USA  Chuang Li  Research
Center for Quantum Sensing, Intelligent Perception Research Institute,
Zhejiang Lab, Hangzhou 310000, China  Shuai Wang  School of Mathematics and
Physics, Jiangsu University of Technology, Changzhou 213001, China

###### Abstract

SU(1,1) interferometers, based on the usage of nonlinear elements, are
superior to passive interferometers in phase sensitivity. However, the SU(1,1)
interferometer cannot make full use of photons carrying phase information as
the second nonlinear element annihilates some of the photons inside. Here, we
focus on improving phase sensitivity and propose a new protocol based on a
modified SU(1,1) interferometer, where the second nonlinear element is
replaced by a beam splitter. We utilize two coherent states as inputs and
implement balanced homodyne measurement at the output. Our analysis suggests
that the protocol we propose can achieve sub-shot-noise-limited phase
sensitivity and is robust against photon loss and background noise. Our work
is important for practical quantum metrology using SU(1,1) interferometers.

##  I Introduction

Phase estimation is an important means of precision measurement, which
provides a route to estimate many physical quantities that cannot be directly
measured via conventional methods. Related to this, optical interferometers
are an effective tool for phase estimation. In general, a conventional
interferometer utilizes a coherent light source and intensity-based
measurement; as a result, the phase sensitivity is bounded by the shot-noise
limit $\Delta{\theta_{\rm{SNL}}}={1/{\sqrt{N}}}$ with $N$ photons inside the
interferometer on average Caves (1981). To improve the phase sensitivity, two
kinds of approaches are proposed. The first approach is to deploy quantum
resources, including exotic non-classical states Pezzé and Smerzi (2008);
Anisimov  _et al._ (2010); Joo  _et al._ (2011); Pezzé and Smerzi (2013);
Israel  _et al._ (2014) and photon-number-resolving measurements Plick  _et
al._ (2010a); Seshadreesan  _et al._ (2013); You  _et al._ (2020); the second
one is to change the evolution process through the use of nonlinear elements.
A typical interferometer designed in terms of the second approach is usually
known as an SU(1,1) interferometer, introduced by Yurke _et al_. Yurke  _et
al._ (1986), where two beam splitters in the conventional interferometer are
replaced by two optical parametric amplifiers (OPAs).

In recent years, a large number of theoretical and experimental studies have
focused on phase estimation using SU(1,1) interferometers. Plick _et al_.
demonstrated that an SU(1,1) interferometer with two coherent states as inputs
can achieve the sub-shot-noise-limited phase sensitivity by using conventional
intensity measurement Plick  _et al._ (2010b). Li _et al_. proposed an SU(1,1)
interferometer fed by coherent and squeezed vacuum states, the phase
sensitivity can approach the Heisenberg limit with balanced homodyne
measurement Li  _et al._ (2014) and parity measurement Li  _et al._ (2016). An
SU(1,1) interferometer using squeezed coherent state and coherent state was
presented by Hu _et al_., they studied the phase sensitivity with balanced
homodyne measurement and discussed the optimal squeezing fraction Hu  _et al._
(2017). Meanwhile, there has been an interest in modified SU(1,1)
interferometers, such as pumped-up and truncated SU(1,1) interferometers
Szigeti  _et al._ (2017); Anderson  _et al._ (2017a, b). In addition, the
phase sensitivity limit of an SU(1,1) interferometer was studied. Gong _et
al_. analyzed the quantum Cramér-Rao bounds of various Gaussian inputs Gong
_et al._ (2017). More recently, You _et al_. discussed the conclusive
sensitivity limit for a protocol using single-mode inputs. You  _et al._
(2019).

Indeed, the phase sensitivities of most protocols cannot achieve or approach
the corresponding quantum Cramér-Rao bound, although they are capable of
reaching the Heisenberg limit Gong  _et al._ (2017). A major reason for this
result originates from the drawback of SU(1,1) interferometers: many photons
stimulated by the first OPA are absorbed by the second OPA. Therefore, only a
part of photons carrying phase information contribute to the measurement,
which limits the phase sensitivity. In this paper, we present a new protocol
for phase estimation using a modified SU(1,1) interferometer; the second OPA
in the conventional SU(1,1) interferometer is replaced by a beam splitter. Two
coherent states are used as inputs, and balanced homodyne measurement is
implemented at the output. The phase sensitivity of our protocol is compared
with the previous SU(1,1)-type protocol and its corresponding quantum Cramér-
Rao bound. For practical purposes, we also analyze the effects of photon loss
and background noise on phase sensitivity of our protocol.

The paper is organized as follows. Sec. II introduces model of the modified
SU(1,1) interferometer and analyzes the phase sensitivity in an ideal
scenario. In Sec. III, we compare the phase sensitivity of our protocol with
that of the previous SU(1,1)-type protocol, and analyze the optimality of our
protocol by calculating the quantum Cramér-Rao bound. The effects of photon
loss and background noise on phase sensitivity are discussed in Sec. IV.
Finally, we summarize our results in Sec. V.

##  II Model and phase sensitivity of modified SU(1,1) interferometer

###  II.1 Model

The schematic of the modified SU(1,1) interferometer is shown in Fig. 1. Two
coherent states, $\left|\alpha\right\rangle$ and $\left|\beta\right\rangle$,
are injected into a pumped OPA. Then the generated states undergo phase shift
of $\theta$ on one of the paths. After the unknown phase shift, these states
are recombined by a beam splitter with transmittance of $\eta$. Finally,
balanced homodyne measurement is implemented at one of the outputs.

Figure 1: Schematic of phase estimation protocol based on a modified SU(1,1)
interferometer. A beam splitter, BS1, takes the place of the second OPA in a
conventional SU(1,1) interferometer. One of the paths is then measured using
balanced homodyne detection. Here the transmittance of BS1 and BS2 is given by
$\eta$ and 1/2, respectively. The creation and annihilation operators for each
stage are shown and given in the Appendix B. OPA, optical parametric
amplifier; BS, beam splitter; LO, local oscillator; D, detector; S,
subtracter.

We assume that phase and gain parameter of the OPA are 0 and $g$,
respectively; accordingly, the complex amplitudes of two coherent states can
be expressed as $\alpha=\left|\alpha\right|{e^{i\varphi}}$ and
$\beta=\left|\beta\right|{e^{i\delta}}$. Since the OPA changes the mean photon
number of the input state, we use the mean photon number inside our
interferometer as the benchmark. Thus the mean photon number is given by

$\displaystyle{N}={\kern 0.5pt}$
$\displaystyle\langle{\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{b}_{1}^{\dagger}\hat{b}_{1}}\rangle$
$\displaystyle={\kern 1.0pt}{\kern 1.0pt}$
$\displaystyle({{{\left|\alpha\right|}^{2}}+{{\left|\beta\right|}^{2}}}){\cosh}\left({2g}\right)+2{\sinh^{2}}g$
$\displaystyle+2\left|{\alpha\beta}\right|\sinh\left({2g}\right)\cos\left({\varphi+\delta}\right).$
(1)

It is straightforward to show that when the total mean photon number of the
input ${N}_{\text{in}}=\left|\alpha\right|^{2}+\left|\beta\right|^{2}$ is
constant, $\left|\alpha\right|=\left|\beta\right|$ and $\varphi+\delta=0$ are
the optimal conditions to obtain maximal mean photon number ${N}_{\text{max}}$
in our protocol. Throughout this paper, we assume such conditions are always
met. That is, the mean photon number in our protocol is

$\displaystyle{N}=2{\left|\alpha\right|^{2}}{e^{2g}}+2{\sinh^{2}}g.$ (2)

###  II.2 Phase sensitivity

The measurement strategy in our protocol is balanced homodyne measurement,
which can measure the quadrature of output by mixing with a strong coherent
state. The phase therefore can be estimated from the difference between the
two detectors D1 and D2. Here we consider amplitude quadrature $\hat{X}$ as
measured observable, and the measurement operator can be written as

$\hat{X}=\hat{a}_{2}^{\dagger}+\hat{a}_{2}.$ (3)

The creation operator of the output modes is given by (See Appendix B)

$\displaystyle\hat{a}_{2}^{\dagger}={\kern 1.0pt}{\kern 1.0pt}$
$\displaystyle\sqrt{\eta}{e^{-i\theta}}({\hat{a}_{0}^{\dagger}\cosh
g+\hat{b}_{0}\sinh g})$
$\displaystyle-i\sqrt{1-\eta}({\hat{b}_{0}^{\dagger}\cosh g+\hat{a}_{0}\sinh
g}).$ (4)

By using Eqs. (3) and (4), the expectation value of measurement operator is
found to be

$\langle{\hat{X}}\rangle=2{e^{g}}\sqrt{\eta}\left|{\alpha}\right|\cos\theta.$
(5)

Using the definition of visibility Dowling (2008),

${\cal
V}=\frac{{{{\langle\hat{X}\rangle}_{\max}}-{{\langle\hat{X}\rangle}_{\min}}}}{{|{{{\langle\hat{X}\rangle}_{\max}}}|+|{{{\langle\hat{X}\rangle}_{\min}}}|}},$
(6)

we can find that the visibility of $\langle{\hat{X}}\rangle$ is 100% for any
transmittance.

Furthermore, we can calculate that

$\displaystyle\langle{{{\hat{X}}^{2}}}\rangle={\kern 1.0pt}{\kern 1.0pt}$
$\displaystyle
2\eta{\left|\alpha\right|^{2}}{e^{2g}}\left[{\cos\left({2\theta}\right)+1}\right]+\cosh\left({2g}\right)$
$\displaystyle-2\sinh\left({2g}\right)\sqrt{\eta\left({1-\eta}\right)}\sin\theta.$
(7)

The phase sensitivity of our protocol can be calculated via error propagation
formula, defined as Sparaciari  _et al._ (2016)

$\Delta\theta=\frac{{\sqrt{\langle{{{\hat{X}}^{2}}}\rangle-{{\langle{\hat{X}}\rangle^{2}}}}}}{|{{{\partial{\langle{\hat{X}}\rangle}}\mathord{\left/{\vphantom{{\partial
X}{\partial\theta}}}\right.\kern-1.2pt}{\partial\theta}}}|}.$ (8)

Plugging Eqs. (5) and (7) into Eq. (8), we get

$\Delta\theta=\frac{{\sqrt{\cosh\left({2g}\right)-2\sqrt{\eta\left({1-\eta}\right)}\sinh\left({2g}\right)\sin\theta}}}{{2\sqrt{\eta}{e^{g}}\left|{\alpha\sin\theta}\right|}}.$
(9)

Clearly, the phase sensitivity depends on the phase shift $\theta$. Evidently,
we obtain the optimal phase sensitivity at $\theta=\pi/2$.

Additionally, from Eq. (9), the phase sensitivity also depends on the value of
transmittance $\eta$. We can obtain the optimal transmittance
$\eta_{\rm{opt}}$ by taking the derivative of Eq. (9) with respect to $\eta$.
When a fixed gain parameter $g$, the optimal transmittance $\eta_{\rm{opt}}$
is given by

${\eta_{{\rm{opt}}}}=\frac{{{{\cosh}^{2}}\left({2g}\right)}}{{{{\sinh}^{2}}\left({2g}\right)+{{\cosh}^{2}}\left({2g}\right)}}=\frac{1+\text{sech}(4g)}{2}.$
(10)

Figure 2: The optimal transmittance $\eta_{\rm{opt}}$ of the beam splitter as
a function of gain parameter $g$.

This reveals that the optimal transmittance depends solely on the gain
parameter $g$, regardless of the mean photon number of the input. If the gain
parameter is known, then the transmittance $\eta$ can be chosen such that the
phase sensitivity reaches its minimum. We note that this effect is pronounced
at low gain parameters. In Fig. 2, we illustrate the relationship between the
optimal transmittance $\eta_{\rm{opt}}$ and gain parameter $g$. From Eq. (10)
and Fig. 2, we conclude that the optimal transmittance approaches 0.5 when $g$
is large. This suggests that a 50-50 beam splitter is nearly optimal for our
protocol with a high gain parameter. Trivially, when $g=0$, our protocol
reduces to a classical interferometer, which is consistent with the result
reported in Ref. Zhang  _et al._ (2020).

In the following sections, the transmittance is fixed to be $1/2$ for the sake
of simplicity. Related to this, the optimal phase sensitivity turns out to be

$\Delta{\theta_{\min}}=\frac{1}{{\sqrt{2}{e^{2g}}\left|\alpha\right|}}.$ (11)

According to the definition of mean photon number in Eq. (2), the shot-noise
limit can be expressed as

$\displaystyle\Delta{\theta_{\rm{SNL}}}=\frac{1}{{\sqrt{2{{\left|\alpha\right|}^{2}}{e^{2g}}+2{{\sinh}^{2}}g}}}.$
(12)

In general, high-intensity coherent states ($\left|\alpha\right|^{2}\gg 1$)
will be used as inputs, i.e.,
${{{\left|\alpha\right|}^{2}}{e^{2g}}\gg{{\sinh}^{2}}g}$. Then the shot-noise
limit is approximately recast as

$\displaystyle\Delta{\theta_{\rm{SNL}}}\simeq\frac{1}{{\sqrt{2}{e^{g}{{{\left|\alpha\right|}}}}}}.$
(13)

By comparing Eqs. (11) and (13), one can find that the phase sensitivity of
our protocol is superior to the shot-noise limit by a factor of $e^{g}$.

##  III Comparison with the SU(1,1)-type protocol and the quantum Cramér-Rao
bound

In the previous section, we calculated the phase sensitivity of our protocol.
In this section, we analyze the advantage of our protocol over the previous
SU(1,1)-type protocol. Then we compare the phase sensitivity with the quantum
Cramér-Rao bound.

###  III.1 Comparison with the SU(1,1)-type protocol

In Ref. Li  _et al._ (2016), Li _et al_. used different measurement strategies
to analyze the phase sensitivity of an SU(1,1)-type protocol with two equal
coherent states as inputs. Their results indicate that balanced homodyne
measurement gives the optimal phase sensitivity, which has the following form

$\displaystyle\Delta{\theta_{\rm{SU(1,1)}}}=\frac{1}{{2e^{g}\left|\alpha\right|\sinh{g}}}.$
(14)

By comparing Eqs. (11) with (14), one can find that the phase sensitivity of
our protocol outperforms the previous SU(1,1)-type protocol. To quantify the
improvement, we calculate the improvement factor ${\cal K}$, given by

$\displaystyle{\cal K}=\frac{{\sqrt{2}{e^{g}}}}{{{e^{g}}-{e^{-g}}}}.$ (15)

For a high gain parameter, the phase sensitivity is improved by a factor of
$\sqrt{2}$. With the decrease of gain parameter, the improvement factor of
phase sensitivity increases rapidly. That is, our protocol can give at least a
factor of $\sqrt{2}$ improvement in phase sensitivity in comparison with an
SU(1,1)-type protocol.

In addition, noticing that the complexity of the measurement system is
simplified. The modified SU(1,1) interferometer circumvents problems
encountered by the SU(1,1) interferometer, such as pump phase-matching and
gain-matching Li  _et al._ (2014). Hence, our protocol remains superior in
terms of sensitivity and complexity when compared with an SU(1,1)-type
protocol.

###  III.2 Comparison with the quantum Cramér-Rao bound

Now we focus on the analysis of phase sensitivity limit of our protocol. The
ultimate limit of phase sensitivity is given by the quantum Cramér-Rao bound,
$\Delta\theta_{\rm{QCRB}}=1/{\sqrt{\cal F}}$ with $\cal F$ being the QFI. In
Ref. Gong  _et al._ (2017), Gong _et al_. demonstrated that three different
phase configurations would give three different QFIs with two coherent states
as inputs. This is due to the asymmetric gain of the OPA; in other words, we
have
$\langle{\hat{a}_{1}^{\dagger}\hat{a}_{1}\rangle\neq\langle\hat{b}_{1}^{\dagger}\hat{b}_{1}}\rangle$
for two unequal coherent states ($\left|\alpha\right|\neq\left|\beta\right|$).
For every phase configuration, the number of photons experiencing the
estimated phase is different. Therefore, the QFIs calculated from different
phase configurations are different.

In particular, when $\left|\alpha\right|=\left|\beta\right|$, all phase
configurations give the same QFI, which is found to be

$\displaystyle{\cal
F}=2{\left|\alpha\right|^{2}}({{e^{4g}}+1})+{\sinh^{2}}\left({2g}\right).$
(16)

The corresponding quantum Cramér-Rao bound is given by

$\displaystyle\Delta\theta_{\rm{QCRB}}=\frac{1}{{\sqrt{2{{\left|\alpha\right|}^{2}}\left({{e^{4g}}+1}\right)+{{\sinh}^{2}}\left({2g}\right)}}}.$
(17)

When $\left|\alpha\right|^{2}\gg 1$ and $e^{4g}\gg 1$, we have

$\displaystyle\Delta\theta_{\rm{QCRB}}\simeq\frac{1}{{\sqrt{2}{e^{2g}}{{{\left|\alpha\right|}}}}}.$
(18)

One can find that Eqs. (11) and (18) are the same, meaning that the phase
sensitivity of our protocol approaches the quantum Cramér-Rao bound with the
increase of $\left|\alpha\right|^{2}$ and $g$. That is, for high-intensity
coherent states and high gain parameters, balanced homodyne measurement is
approximately optimal strategy for our protocol.

In Fig. 3, we show the ratio between the ultimate phase sensitivity
$\Delta\theta_{\rm{QCRB}}$, given by quantum Cramér-Rao bound, and the optimal
phase sensitivity $\Delta\theta_{\rm{min}}$ as a function of gain parameter
$g$ and mean photon number $\left|\alpha\right|^{2}$ in the coherent state. We
can see that based on the current experimental techniques Manceau  _et al._
(2017), it is feasible to achieve phase sensitivity close to the one governed
by the quantum Cramér-Rao bound.

Figure 3: Ratio $\Delta\theta_{\rm{QCRB}}/\Delta\theta_{\rm{min}}$ as a
function of gain parameter $g$ and mean photon number
$\left|\alpha\right|^{2}$ in the coherent state. The ratio is given such that
a bigger ratio implies a better performance.

##  IV Analysis of robustness regarding photon loss and background noise

In Sec. II, we analyzed the phase sensitivity in an ideal scenario. In
practical applications, photon loss and background noise are inevitable and
they will limit the attainable phase sensitivity. In this section, we discuss
the effects of these two factors on phase sensitivity.

The process of photon loss can be modeled by adding two fictitious beam
splitters with reflectivity of $L_{1}$ and $L_{2}$ in two paths PL . Photons
reflected by the fictitious beam splitters are regarded as photon loss, and
the lossy rate in each path is equal to the reflectivity of corresponding
fictitious beam splitter. Background noise refers to photons in the
surrounding environment, which enter the interferometer through another input
port of the fictitious beam splitter.

In the optical regime, the mean photon number at room temperature is
approximately $10^{-20}$ Gard  _et al._ (2017) and the environmental mode is
treated as vacuum. However, when the measurement system is immersed in the
environment with a small amount of stray light, the environmental mode is
treated as a thermal state with $n_{\rm{th}}$ photons on average. Generally,
the mean photon number of background noise is much weaker than those of the
input states. In this paper, $n_{\rm{th}}=1$ is used to simulated the
intensity of background noise. Related to this, $n_{\rm{th}}=0$ indicates that
there is only photon loss, and $n_{\rm{th}}=1$ indicates that photon loss and
background noise exist simultaneously.

For a lossy scenario, the operator relationship between the input and output
modes is rewritten as

$\displaystyle\hat{a}_{2}^{\dagger}={\kern 0.4pt}$
$\displaystyle\sqrt{\frac{{1-{L_{1}}}}{2}}{e^{-i\theta}}({\hat{a}_{0}^{\dagger}\cosh
g+\hat{b}_{0}\sinh g})-i\sqrt{\frac{{{L_{1}}}}{2}}\hat{c}^{\dagger}$
$\displaystyle-i\sqrt{\frac{{1-{L_{2}}}}{2}}({\hat{b}_{0}^{\dagger}\cosh
g+\hat{a}_{0}\sinh g})-\sqrt{\frac{{{L_{2}}}}{2}}\hat{d}^{\dagger},$ (19)

where $\hat{c}^{\dagger}$ and $\hat{d}^{\dagger}$ represent the creation
operators for environmental modes.

We can recalculate that

$\langle{{{\hat{X}}_{\rm{L}}}}\rangle=\sqrt{2\left({1-{L_{1}}}\right)}{e^{g}}\left|\alpha\right|{\rm{cos}}\theta,$
(20)

and

$\displaystyle\langle{\hat{X}_{\rm L}^{2}}\rangle={\kern 0.5pt}$
$\displaystyle
2{\left|\alpha\right|^{2}}{e^{2g}}\left({1-{L_{1}}}\right)\cos^{2}\theta+\left({{L_{1}}+{L_{2}}}\right){n_{\rm{th}}}+1$
$\displaystyle-\sqrt{\left({1-{L_{1}}}\right)\left({1-{L_{2}}}\right)}\sinh\left({2g}\right)\sin\theta$
$\displaystyle+\left({2-{L_{1}}-{L_{2}}}\right){{\sinh}^{2}}g.$ (21)

In the limit of $L_{1}=L_{2}=1$, the inputs are completely lossy and the
outputs are thermal states. Combining Eqs. (20) and (21), we have
$\langle{\hat{X}_{\rm
L}^{2}}\rangle-\langle{{{\hat{X}}_{\rm{L}}}}\rangle^{2}=2{n_{\rm{th}}}+1$,
which matches the variance of a thermal state without information on the
estimated phase. In addition, Eq. (20) means that the visibility remains 100%
regardless of lossy rates.

Since the lengths of two paths are close in a practical scenario, we consider
$L_{1}=L_{2}=L$ from now on. The phase sensitivity is given by

$\displaystyle\Delta\theta_{\rm
L}=\frac{{\sqrt{\left({1-L}\right)\left[{2{{\sinh}^{2}}g-\sinh\left({2g}\right)\sin\theta}\right]+2L{n_{\rm{th}}}+1}}}{{\sqrt{2\left({1-L}\right)}{e^{g}}\left|{\alpha\sin\theta}\right|}}.$
(22)

When the estimated phase sits at the optimal phase point $\theta=\pi/2$, the
phase sensitivity reaches the minimum, which can be expressed as

$\displaystyle\Delta{\theta_{\rm
L,\min}}=\frac{{\sqrt{1+\gamma}}}{{\sqrt{2}{e^{2g}}\left|\alpha\right|}}$ (23)

with

$\displaystyle\gamma=\frac{L}{{1-L}}\left({2{n_{{\rm{th}}}}+1}\right){e^{2g}}.$
(24)

It can be seen that $\gamma$ is the increment of measurement variance induced
by photon loss and background noise. As has been previously pointed out, phase
sensitivity in an ideal scenario can outperform the shot-noise limit. Here we
analyze the maximum value of lossy rate when the phase sensitivity drops below
the shot-noise limit. By solving the equation with respect to Eqs. (12) and
(23),

$\displaystyle\frac{{{1-L_{\rm{max}}+L_{\rm{max}}\left({2{n_{{\rm{th}}}}+1}\right){e^{2g}}}}}{{{2\left({1-L_{\rm{max}}}\right)}{e^{4g}}\left|\alpha\right|^{2}}}=\frac{1}{{{2{{\left|\alpha\right|}^{2}}{e^{2g}}+2{{\sinh}^{2}}g}}},$
(25)

one can calculate the maximal tolerable lossy rate $L_{\rm{max}}$.

We first consider the case of only photon loss, Fig. 4 gives the relationship
between the maximal tolerable lossy rate and gain parameter with different
mean photon numbers in coherent states. The mean number of photons we use is
easy to engineer under the current experimental techniques; meanwhile, the
condition $|\alpha|^{2}\gg 1$ is satisfied. One can find that the effect of
mean photon number on the maximal tolerable lossy rate is negligible. With the
increase of gain parameter, the maximal tolerable lossy rate increases
rapidly; after about $g=2$, this rate tends to the limit, which is equal to
50%. This result indicates that increasing gain parameter can improve the
robustness of our protocol.

Figure 4: Maximal tolerable lossy rate as a function of gain parameter in the
presence of photon loss, where $n_{\rm{th}}=0$.

In Fig. 5, we show the dependence of maximal tolerable lossy rate on gain
parameter in the presence of both photon loss and background noise. Similarly,
mean photon number has almost no effect on the maximal tolerable lossy rate.
The maximal tolerable lossy rate tends to the limit after about $g=2$, and the
limit is 25%.

Figure 5: Maximal tolerable lossy rate as a function of gain parameter in the
presence of photon loss and background noise, where $n_{\rm{th}}=1$.

By comparing Fig. 4 with Fig. 5, it can be seen that the effect of background
noise is significant. Although $n_{\rm{th}}\ll|\alpha|^{2}$, the maximal
tolerable lossy rate decreases to half of that without background noise. To
determine the limit of the maximal tolerable lossy rate, we can calculate the
combination of Eqs. (13) and (23) as follows

$\displaystyle\frac{{1-L_{\rm{max}}+L_{\rm{max}}\left({2{n_{{\rm{th}}}}+1}\right){e^{2g}}}}{{2\left({1-L_{\rm{max}}}\right){e^{4g}}{{\left|\alpha\right|}^{2}}}}=\frac{1}{{2{{\left|\alpha\right|}^{2}}{e^{2g}}}}.$
(26)

The approximate solution for Eq. (26) is found to be

$\displaystyle\mathop{\lim}\limits_{g\gg 1,\left|\alpha\right|\gg
1}{L_{\max}}=\frac{1}{{2\left({{n_{{\rm{th}}}}+1}\right)}}.$ (27)

Obviously, we have ${L_{\max}}=0.5$ and ${L_{\max}}=0.25$ for $n_{\rm{th}}=0$
and $n_{\rm{th}}=1$, respectively.

##  V Conclusion

In conclusion, we propose a new protocol for phase estimation using a modified
SU(1,1) interferometer. The second OPA in the conventional SU(1,1)
interferometer is replaced by a 50-50 beam splitter. We use two coherent
states as inputs and perform balanced homodyne measurement at the output. In
an ideal scenario, the phase sensitivity is superior to the shot-noise limit
by a factor of $e^{g}$ with coherent inputs. Additionally, our protocol shows
the potential to surpass conventional SU(1,1) interferometer in phase
sensitivity. The improvement factor of the phase sensitivity is $\sqrt{2}$ at
high gain parameters and increases rapidly with the decrease of the gain
parameter. Finally, we analyze the effects of photon loss and background noise
on phase sensitivity. The results indicate that our protocol is robust against
these two factors. For high gain parameters, the phase sensitivity of our
protocol can outperform the shot-noise limit even when the lossy rate is
relatively high. We believe that our protocol is feasible to achieve phase
sensitivity close to the quantum Cramér-Rao bound based on current
experimental techniques.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
(No. 11905160 and No. 11404040) and sponsored by Qing Lan Project of the
Higher Educations of Jiangsu Province of China.

##  Appendix A The exact expression for phase sensitivity calculated from
classical Fisher information

For many measurement strategies, error propagation formula is a simple but
approximate method to calculate phase sensitivity. An exact expression for
phase sensitivity equals the inverse of the square root of the classical
Fisher information, i.e., $\Delta\theta=1/\sqrt{{\cal F}_{\rm c}}$. As for our
protocol–Gaussian inputs and Gaussian measurement–the classical Fisher
information is described as

$\displaystyle{{\cal
F}_{\rm{c}}}=\frac{{({\partial\langle{\hat{X}}\rangle}/{\partial\theta})^{2}}}{\langle{\hat{X}^{2}}\rangle-\langle{\hat{X}}\rangle^{2}}+\frac{{2(\partial\Delta{\hat{X}}/{\partial\theta})^{2}}}{\langle{\hat{X}^{2}}\rangle-\langle{\hat{X}}\rangle^{2}}$
(28)

with
$\Delta{\hat{X}}={\sqrt{{\langle{\hat{X}^{2}}\rangle-\langle{\hat{X}}\rangle^{2}}}}$.

One can find that the first term is related to error propagation formula. The
second term at the optimal phase point, $\theta=\pi/2$, is 0 since the
derivative in the numerator can be decomposed into the product of a prefactor
and a cosine function, $\cos(\theta)$. This indicates that error propagation
formula and the classical Fisher information give the same result for the
optimal phase sensitivity. For other phase points, the phase sensitivity
calculated from error propagation formula is greater than true value as the
second term is non-negative. These results hold true for the lossy scenario.
Overall, it is sufficient to calculate the phase sensitivity via error
propagation formula, for we merely are interested in the optimal phase
sensitivity of our protocol.

##  Appendix B The expectation values of measurement operators with ideal and
lossy scenarios

In an ideal scenario, the input-output field operator transformations of a
beam splitter and those of an optical parametric amplifier are given by

$\displaystyle\hat{a}_{2}$
$\displaystyle=\sqrt{\eta}\hat{a}_{1}e^{i\theta}+i\sqrt{1-\eta}\hat{b}_{1},$
(29) $\displaystyle\hat{b}_{2}$
$\displaystyle=i\sqrt{1-\eta}\hat{a}_{1}e^{i\theta}+\sqrt{\eta}\hat{b}_{1},$
(30) $\displaystyle\hat{a}_{1}$ $\displaystyle=\hat{a}_{0}\cosh
g+\hat{b}_{0}^{\dagger}\sinh g,$ (31) $\displaystyle\hat{b}_{1}$
$\displaystyle=\hat{b}_{0}\cosh g+\hat{a}_{0}^{\dagger}\sinh g.$ (32)

Regarding two equal coherent states ${\rho_{A}}\otimes{\rho_{B}}$
(${\rho_{A}}={\rho_{B}}=\left|\alpha\right\rangle\left\langle\alpha\right|$),
we have

$\displaystyle\langle{\hat{a}_{2}^{\dagger}\hat{a}_{2}^{\dagger}+\hat{a}_{2}\hat{a}_{2}}\rangle=2{\left|\alpha\right|^{2}}{e^{2g}}[{\eta\cos\left({2\theta}\right)-2\sqrt{\eta\left({1-\eta}\right)}\sin\theta-\left({1-\eta}\right)}]-2\sinh\left({2g}\right)\sqrt{\eta\left({1-\eta}\right)}\sin\theta$
(33)

and

$\displaystyle\langle{\hat{a}_{2}^{\dagger}\hat{a}_{2}}\rangle={\left|\alpha\right|^{2}}{e^{2g}}[{1+2\sqrt{\eta\left({1-\eta}\right)}\sin\theta}]+{\sinh^{2}}g.$
(34)

In a lossy scenario, the input-output field operator transformation of the
whole interferometer is given by Eq. (19). Related to this, the inputs are
${\rho_{A}}\otimes{\rho_{B}}\otimes{\rho_{E}}\otimes{\rho_{E}}$, where
${\rho_{E}}$ represents a thermal state with $n_{\rm{th}}$ photons on average.
Further, we have

$\displaystyle\langle{\hat{a}_{2}^{\dagger}\hat{a}_{2}^{\dagger}+\hat{a}_{2}\hat{a}_{2}}\rangle={\kern
1.0pt}$
$\displaystyle{\left|\alpha\right|^{2}}{e^{2g}}[{\left({1-{L_{1}}}\right)\cos\left({2\theta}\right)-2\sqrt{\left({1-{L_{1}}}\right)\left({1-{L_{2}}}\right)}\sin\theta-\left({1-{L_{2}}}\right)}]$
$\displaystyle-\sqrt{\left({1-{L_{1}}}\right)\left({1-{L_{2}}}\right)}\sinh\left({2g}\right)\sin\theta$
(35)

and

$\displaystyle\langle{\hat{a}_{2}^{\dagger}\hat{a}_{2}}\rangle=\left({1-\frac{{{L_{1}}+{L_{2}}}}{2}}\right)({{{\left|\alpha\right|}^{2}}{e^{2g}}+{{\sinh}^{2}}g})+{\left|\alpha\right|^{2}}{e^{2g}}\sqrt{\left({1-{L_{1}}}\right)\left({1-{L_{2}}}\right)}\sin\theta+\frac{{{L_{1}}+{L_{2}}}}{2}{n_{\rm{th}}}.$
(36)

## References

  * Caves (1981) Carlton M. Caves, “Quantum-mechanical noise in an interferometer,” Phys. Rev. D 23, 1693–1708 (1981).
  * Pezzé and Smerzi (2008) Luca Pezzé and Augusto Smerzi, “Mach-Zehnder interferometry at the Heisenberg limit with coherent and squeezed-vacuum light,” Phys. Rev. Lett. 100, 073601 (2008).
  * Anisimov  _et al._ (2010) Petr M. Anisimov, Gretchen M. Raterman, Aravind Chiruvelli, William N. Plick, Sean D. Huver, Hwang Lee, and Jonathan P. Dowling, “Quantum metrology with two-mode squeezed vacuum: Parity detection beats the Heisenberg limit,” Phys. Rev. Lett. 104, 103602 (2010).
  * Joo  _et al._ (2011) Jaewoo Joo, William J. Munro, and Timothy P. Spiller, “Quantum metrology with entangled coherent states,” Phys. Rev. Lett. 107, 083601 (2011).
  * Pezzé and Smerzi (2013) Luca Pezzé and Augusto Smerzi, “Ultrasensitive two-mode interferometry with single-mode number squeezing,” Phys. Rev. Lett. 110, 163604 (2013).
  * Israel  _et al._ (2014) Yonatan Israel, Shamir Rosen, and Yaron Silberberg, “Supersensitive polarization microscopy using NOON states of light,” Phys. Rev. Lett. 112, 103604 (2014).
  * Plick  _et al._ (2010a) William N Plick, Petr M Anisimov, Jonathan P Dowling, Hwang Lee, and Girish S Agarwal, “Parity detection in quantum optical metrology without number-resolving detectors,” New J. Phys. 12, 113025 (2010a).
  * Seshadreesan  _et al._ (2013) Kaushik P. Seshadreesan, Sejong Kim, Jonathan P. Dowling, and Hwang Lee, “Phase estimation at the quantum Cramér-Rao bound via parity detection,” Phys. Rev. A 87, 043833 (2013).
  * You  _et al._ (2020) Chenglong You, Mingyuan Hong, Peter Bierhorst, Adriana E. Lita, Scott Glancy, Steve Kolthammer, Emanuel Knill, Sae Woo Nam, Richard P. Mirin, Omar S. Magana-Loaiza, and Thomas Gerrits, “Multiphoton quantum metrology with neither pre- nor post-selected measurements,” (2020), arXiv:2011.02454 [quant-ph] .
  * Yurke  _et al._ (1986) Bernard Yurke, Samuel L. McCall, and John R. Klauder, “SU(2) and SU(1,1) interferometers,” Phys. Rev. A 33, 4033–4054 (1986).
  * Plick  _et al._ (2010b) William N Plick, Jonathan P Dowling, and Girish S Agarwal, “Coherent-light-boosted, sub-shot noise, quantum interferometry,” New J. Phys. 12, 083014 (2010b).
  * Li  _et al._ (2014) Dong Li, Chun-Hua Yuan, Z Y Ou, and Weiping Zhang, “The phase sensitivity of an SU(1,1) interferometer with coherent and squeezed-vacuum light,” New J. Phys. 16, 073020 (2014).
  * Li  _et al._ (2016) Dong Li, Bryan T. Gard, Yang Gao, Chun-Hua Yuan, Weiping Zhang, Hwang Lee, and Jonathan P. Dowling, “Phase sensitivity at the Heisenberg limit in an SU(1,1) interferometer via parity detection,” Phys. Rev. A 94, 063840 (2016).
  * Hu  _et al._ (2017) Xiao-Yu Hu, Chao-Ping Wei, Ya-Fei Yu, and Zhi-Ming Zhang, “Enhanced phase sensitivity of an SU(1,1) interferometer with displaced squeezed vacuum light,” Front. Phys. 11, 114203 (2017).
  * Szigeti  _et al._ (2017) Stuart S. Szigeti, Robert J. Lewis-Swan, and Simon A. Haine, “Pumped-up SU(1,1) interferometry,” Phys. Rev. Lett. 118, 150401 (2017).
  * Anderson  _et al._ (2017a) Brian E. Anderson, Prasoon Gupta, Bonnie L. Schmittberger, Travis Horrom, Carla Hermann-Avigliano, Kevin M. Jones, and Paul D. Lett, “Phase sensing beyond the standard quantum limit with a variation on the SU(1,1) interferometer,” Optica 4, 752–756 (2017a).
  * Anderson  _et al._ (2017b) Brian E. Anderson, Bonnie L. Schmittberger, Prasoon Gupta, Kevin M. Jones, and Paul D. Lett, “Optimal phase measurements with bright- and vacuum-seeded SU(1,1) interferometers,” Phys. Rev. A 95, 063843 (2017b).
  * Gong  _et al._ (2017) Qian-Kun Gong, Dong Li, Chun-Hua Yuan, Ze-Yu Qu, and Wei-Ping Zhang, “Phase estimation of phase shifts in two arms for an su(1,1) interferometer with coherent and squeezed vacuum states,” Chin. Phys. B 26, 94205 (2017).
  * You  _et al._ (2019) Chenglong You, Sushovit Adhikari, Xiaoping Ma, Masahide Sasaki, Masahiro Takeoka, and Jonathan P. Dowling, “Conclusive precision bounds for SU(1,1) interferometers,” Phys. Rev. A 99, 042122 (2019).
  * Dowling (2008) Jonathan P. Dowling, “Quantum optical metrology — the lowdown on high-N00N states,” Contemp. Phys. 49, 125–143 (2008).
  * Sparaciari  _et al._ (2016) Carlo Sparaciari, Stefano Olivares, and Matteo G. A. Paris, “Gaussian-state interferometry with passive and active elements,” Phys. Rev. A 93, 023810 (2016).
  * Zhang  _et al._ (2020) Jian-Dong Zhang, Zi-Jing Zhang, Long-Zhu Cen, Jun-Yan Hu, and Yuan Zhao, “Classical inputs and measurements enable phase sensitivity beyond the shot-noise limit,” (2020), arXiv:2001.00160 .
  * Manceau  _et al._ (2017) Mathieu Manceau, Gerd Leuchs, Farid Khalili, and Maria Chekhova, “Detection loss tolerant supersensitive phase measurement with an SU(1,1) interferometer,” Phys. Rev. Lett. 119, 223604 (2017).
  * (24) In an SU(1,1) interferometer, internal loss has a more detrimental effect on the phase sensitivity than external loss as the vacuum noise is amplified by the second OPA. Owing to the fact that a beam splitter does not amplify vacuum noise, the effects of internal loss and external loss on the phase sensitivity are the same in the modified SU(1,1) interferometer.
  * Gard  _et al._ (2017) Bryan T Gard, Chenglong You, Devendra K Mishra, Robinjeet Singh, Hwang Lee, Thomas R Corbitt, and Jonathan P Dowling, “Nearly optimal measurement schemes in a noisy Mach-Zehnder interferometer with coherent and squeezed vacuum,” EPJ Quantum Technol. 4, 4 (2017).


# Time-Fractional Approach to the Electrochemical Impedance: The Displacement
Current

G. Barbero1,2, L. R. Evangelista1,3, E. K. Lenzi4 1Dipartimento di Scienza
Applicata del Politecnico di Torino,
Corso Duca degli Abruzzi 24, 10129 Torino, Italy.
2National Research Nuclear University MEPhI (Moscow Engineering Physics
Institute), Kashirskoye shosse 31, 115409 Moscow, Russian Federation.
3 Departamento de Física, Universidade Estadual de Maringá, Avenida Colombo,
5790, 87020-900 Maringá, Paraná, Brazil.
4 Departamento de Física, Universidade Estadual de Ponta Grossa, Avenida
Carlos Cavalcanti, 4748, 87030-900, Ponta Grossa, Paraná, Brazil.  E-Mails:
<EMAIL_ADDRESS>(G.B.)
<EMAIL_ADDRESS>(L.R.E.)
<EMAIL_ADDRESS>(E.K.L.)

###### Abstract

We establish, in general terms, the conditions to be satisfied by a time-
fractional approach formulation of the Poisson-Nernst-Planck model in order to
guarantee that the total current across the sample be solenoidal, as required
by the Maxwell equation. Only in this case the electric impedance of a cell
can be determined as the ratio between the applied difference of potential and
the current across the cell. We show that in the case of anomalous diffusion,
the model predicts for the electric impedance of the cell a constant phase
element behaviour in the low frequency region. In the parametric curve of the
reactance versus the resistance, the slope coincides with the order of the
fractional time derivative.

American Chemical Society, LaTeX

The Poisson-Nernst-Planck (PNP) diffusional model is formulated to consider a
neutral species that can dissociate into positive and negative charges of
arbitrary mobilities [1, 2]. In the presence of the external electric field
$\mathbf{E}({\textbf{r}},t)$, these charges move, giving rise to currents of
neutral as well as of positive and negative charges. If we denote the bulk
density of these particles, respectively with $n_{n}({\textbf{r}},t)$,
$n_{p}({\textbf{r}},t)$, and $n_{m}({\textbf{r}},t)$ and, likewise, the
current density as ${\textbf{j}}_{n}({\textbf{r}},t)$,
${\textbf{j}}_{p}({\textbf{r}},t)$, and ${\textbf{j}}_{m}({\textbf{r}},t)$,
the continuity equations will be written as

$\displaystyle\frac{\partial}{\partial t}n_{p}({\textbf{r}},t)$
$\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{p}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\frac{\partial}{\partial t}n_{m}({\textbf{r}},t)$
$\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{m}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\frac{\partial}{\partial t}n_{n}({\textbf{r}},t)$
$\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{n}({\textbf{r}},t)-S({\textbf{r}},t),$
(1)

where $S({\textbf{r}},t)$ accounts for a source term related to the
dissociation of neutral particles and recombination of the ions. In the case
of a full dissociation the neutral particles desappear, as well as the source
term. To obtain the set of the fundamental equations of the PNP model, we have
to consider also the Poisson equation, written as

$\nabla\cdot\mathbf{E}({\textbf{r}},t)=\frac{q}{\varepsilon}\left[n_{p}({\textbf{r}},t)-n_{m}({\textbf{r}},t)\right],$
(2)

where $\varepsilon$ is the dielectric permittivity of the insulating medium
free of ions. This equation connects the bulk density of ions of positive and
negative charges, of absolute value $q$, to the actual profile of the electric
field across the sample.

The total electric current is formed by the conduction and the displacement
currents:

${\textbf{j}}({\textbf{r}},t)=q\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]+\varepsilon\frac{\partial\mathbf{E}({\textbf{r}},t)}{\partial
t}.$ (3)

If we now combine the first two of Eqs. (Time-Fractional Approach to the
Electrochemical Impedance: The Displacement Current) with Eq. (3), we conclude
that

$\nabla\cdot{\textbf{j}}=-q\left\\{\frac{\partial}{\partial
t}\left[n_{p}({\textbf{r}},t)-n_{m}({\textbf{r}},t)\right]\right\\}+\varepsilon\nabla\cdot\frac{\partial\mathbf{E}({\textbf{r}},t)}{\partial
t}=0,$ (4)

when the equation of Poisson, Eq. (2), is used. Indeed, the total current has
to be solenoidal, i.e,

$\nabla\cdot{\textbf{j}}({\textbf{r}},t)=\nabla\cdot\left\\{q\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]+\varepsilon\frac{\partial\mathbf{E}({\textbf{r}},t)}{\partial
t}\right\\}=0.$ (5)

In this case, for a one-dimensional problem, the current density $\bf j$ is
position independent. Only in this framework, the concept of electrical
impedance, defined as the ratio between the difference of potential applied to
the system and the total electric current flowing across it, is meaningful
[3]. To proceed with the application of the PNP model to deal with
experimental data obtained in electrolytic cells, the set of Eqs. (Time-
Fractional Approach to the Electrochemical Impedance: The Displacement
Current) and (2) has to be solved for well-defined boundary conditions.
Anyway, the total current flowing though the cell is subjected to the general
condition established in Eq. (5), i.e., the total current has to be
solenoidal.

Perhaps the simplest extension of the problem using derivatives operators of
arbitrary order, commonly called fractional derivatives, should read as the
one we have proposed some years ago [4]:

$\displaystyle\tau^{\gamma-1}\frac{\partial^{\gamma}n_{p}({\textbf{r}},t)}{\partial
t^{\gamma}}$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{p}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\tau^{\gamma-1}\frac{\partial^{\gamma}n_{m}({\textbf{r}},t)}{\partial
t^{\gamma}}$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{m}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\tau^{\gamma-1}\frac{\partial^{\gamma}n_{n}({\textbf{r}},t)}{\partial
t^{\gamma}}$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{n}({\textbf{r}},t)-S({\textbf{r}},t),$
(6)

where $\tau$ is an intrinsic time of the problem. In Eqs. (Time-Fractional
Approach to the Electrochemical Impedance: The Displacement Current) we have
replaced the time derivative present in Eqs. (Time-Fractional Approach to the
Electrochemical Impedance: The Displacement Current) with the Riemann-
Liouville fractional derivative, by promoting the change:

$\frac{\partial n_{p,m}({\textbf{r}},t)}{\partial
t}\longrightarrow\frac{\partial^{\gamma}n_{p,m}({\textbf{r}},t)}{\partial
t^{\gamma}}\longrightarrow{{}_{t_{0}}\textrm{D}}_{t}^{\gamma}n_{p,m}({\textbf{r}},t),$
(7)

and used the definition of the Riemann-Liouville operator as

$\\!\\!\\!\\!\\!\\!\;{}_{t_{0}}{{\textrm{D}}}_{t}^{\gamma}n_{p,m}({\textbf{r}},t)=\frac{1}{\Gamma\left(m-\gamma\right)}\frac{d^{m}}{dt^{m}}\int_{t_{0}}^{t}d\overline{t}\frac{n_{p,m}({\textbf{r}},\overline{t}\;)}{\left(t-\overline{t}\right)^{\gamma+1-m}},$
(8)

where $m-1<\gamma<m$, with $m$ an integer, and $t_{0}$ is related to the
conditions initially imposed to the system. For an analysis of the kind we are
proposing here, focused on the impedance spectroscopy, without loss of
generality we can consider $t_{0}=-\infty$ to study the response of the system
to the periodic applied potential, as required by the impedance spectroscopy
technique. If we now combine the new set of fundamental equations of the
extended model, which are Eqs. (Time-Fractional Approach to the
Electrochemical Impedance: The Displacement Current) and (2), following the
procedure used before, it is mandatory to conclude that
$\nabla\cdot{\textbf{j}}({\textbf{r}},t)=0$, i.e.,

$\nabla\cdot\left\\{q\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]+\varepsilon\tau^{\gamma-1}\frac{\partial^{\gamma}\mathbf{E}({\textbf{r}},t)}{\partial
t^{\gamma}}\right\\}=0.$ (9)

This means that a physically sound extension of the PNP model using fractional
derivatives requires an extended expression in which the displacement current
is also defined in terms of a derivative of arbitrary order. This fact has
profound implications and requires discussing also possible extensions of the
Maxwell equations to the field of the differential operators of arbitrary
order. Before analyzing this important aspect of the problem, to keep the
approach as general as possible, let us consider a further step in the
generalization of the PNP model using the time-fractional derivative.

One possible way to generalize the previous equations is to use the so-called
time-fractional derivative operator of distributed order, defined as [5, 6]:

$\int_{0}^{1}d\gamma p(\gamma)\\!\;_{-\infty}{\rm
D}_{t}^{\gamma}\left(\cdots\right),$ (10)

where $p(\gamma)$ is a distribution of $\gamma$, with $\int_{0}^{1}d\gamma
p(\gamma)=1$, instead of using only a single fractional operator
${}_{-\infty}{\textrm{D}}_{t}^{\gamma}\left(\cdots\right)$. After implementing
this generalization, the set of fundamental equations of the extended PNP
model will be now formed by the following continuity equations:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{0}^{1}\\!\\!d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}n_{p}({\textbf{r}},t)$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{p}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{0}^{1}\\!\\!d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}n_{m}({\textbf{r}},t)$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{m}({\textbf{r}},t)+S({\textbf{r}},t),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\int_{0}^{1}\\!\\!d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}n_{n}({\textbf{r}},t)$ $\displaystyle=$
$\displaystyle-\nabla\cdot{\textbf{j}}_{n}({\textbf{r}},t)-S({\textbf{r}},t),$
(11)

and the equation of Poisson, Eq. (2). Combining again, as before, these
equations, it is possible to show that

$\displaystyle\int_{0}^{1}d\gamma p(\gamma)\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\tau^{\gamma}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}\left[n_{p}({\textbf{r}},t)-n_{m}({\textbf{r}},t)\right]$ (12)
$\displaystyle+$
$\displaystyle\nabla\cdot\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]=0,$

and, consequently, that

$\displaystyle\nabla\\!\cdot\\!\left[q\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]\\!+\\!\frac{\varepsilon}{\tau}\int_{0}^{1}\\!\\!\\!\\!d\gamma
p(\gamma)\tau^{\gamma}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}\mathbf{E}({\textbf{r}},t)\right]=0.$

Equation (Time-Fractional Approach to the Electrochemical Impedance: The
Displacement Current) is now more general than Eq. (9) and implies that the
total electric current

$\displaystyle{\textbf{j}}({\textbf{r}},t)$ $\displaystyle=$ $\displaystyle
q\left[{\textbf{j}}_{p}({\textbf{r}},t)-{\textbf{j}}_{m}({\textbf{r}},t)\right]$
(14) $\displaystyle+$ $\displaystyle\varepsilon\int_{0}^{1}d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}\mathbf{E}({\textbf{r}},t)$

is solenoidal, as in the previous cases. If we consider
$p(\gamma)=\delta(\gamma-1)$, then Eqs. (Time-Fractional Approach to the
Electrochemical Impedance: The Displacement Current) reduce to Eqs.(Time-
Fractional Approach to the Electrochemical Impedance: The Displacement
Current), as well as Eq. (Time-Fractional Approach to the Electrochemical
Impedance: The Displacement Current) reduces to Eq. (5). Other forms of
$p(\gamma)$ may express a superposition of different diffusive regimes. For
instance, $p({\gamma})=A\delta({\gamma}-1)+B\delta({{\gamma}}-\alpha)$, where
$A$ and $B$ are connected with characteristic times, represents a
superposition of a normal diffusion equation (${\gamma}=1$) with a diffusion
equation of arbitrary order ${\gamma}=\alpha$, with $0<\alpha\leq 1$.

In the prototypical case of time-fractional derivative of order $\gamma$, as
in Eq. (9), the displacement current assumes the form:

${\textbf{j}}_{D}({\textbf{r}},t)=\varepsilon\tau^{\gamma-1}\frac{\partial^{\gamma}\mathbf{E}({\textbf{r}},t)}{\partial
t^{\gamma}}.$ (15)

For consistency, the definition of a modified displacement current in terms of
a fractional time-derivative operator requires an appropriated formulation of
a time-fractional electrodynamics, which, needless to say, is a hard task.
Recently, some approaches have been presented to deal with these kinds of
formulation. Analogs for Maxwell’s equations using time-fractional derivatives
in the Riemann-Liouville and Caputo sense have been obtained by introducing
time integro-differentiation of arbitrary order into the fundamental equations
of the electrodynamics of a material media [7]. The analysis shows that the
stochastic nature of charged particles motion in a medium, properly described
in terms of fractional operators, influences the dynamics of an
electromagnetic field. More recently, a formulation of time-fractional
electrodynamics was derived based on the Riemann-Silberstein vector [8]. The
use of this vector together with fractional-order derivatives indicates a way
to write Maxwell’s equations in terms of time-fractional derivatives in a
compact form, which allows for modeling of energy dissipation and dynamics of
electromagnetic systems with memory. These generalizing approaches require a
displacement current exactly in the form we are proposing here. They open the
possibility to interpret the displacement currents, arising in the formulation
of the time-fractional PNP model, in the appropriate electrodynamics’s
framework.

The PNP model is a very useful framework to interpret impedance spectroscopy
data regarding a wide variety of systems and electrolytic cells. In the range
of frequency for which it works particularly well in describing the response
of the ions in the sample to the stimulus of an external field, the electric
current used to determine the electric impedance of the cell has to include
the displacement current related to the time variation of the electric field.
If the diffusion is normal, it is given by ${\bf
j}_{D}=\varepsilon\partial{\bf E}/\partial t$. If the diffusion is anomalous
it is given by (15). The important point is that, in the framework of
anomalous diffusion the displacement current cannot be used in its usual,
i.e., non-fractional form. To analyze this aspect of the problem in more
detail, let us consider a typical problem dealing with the impedance of an
electrolytic cell. The proper application of the technique requires an
electric field in the form:

${\textbf{E}}({\textbf{r}},t)={\textbf{F}}({\textbf{r}})e^{i\omega t},$ (16)

where $\omega=2\pi f$ is the circular frequency of the applied field. From
Eqs. (5) and (Time-Fractional Approach to the Electrochemical Impedance: The
Displacement Current), and again using the equation of Poisson, Eq. (2), we
may write

$\displaystyle\nabla\cdot{\textbf{j}}=-\varepsilon\nabla\cdot\left[\tau^{\gamma-1}\frac{\partial^{\gamma}{\textbf{E}}({\textbf{r}},t)}{\partial
t^{\gamma}}-\frac{\partial{\textbf{E}}({\textbf{r}},t)}{\partial t}\right].$
(17)

Now, if we use the definition of the electric field as in Eq. (16), we obtain:

$\nabla\cdot{\textbf{j}}=-\varepsilon\left(\nabla\cdot{\textbf{F}}\right)(i\omega)\left[(i\tau\omega)^{\gamma-1}-1\right],$
(18)

which is solenoidal only if $\gamma=1$. However, if the difference
$(i\tau\omega)^{\gamma-1}-1$ may be considered negligible, the requirement of
solenoidal current flowing through the system is essentially fulfilled, and
the time-fractional approach to the PNP model may be fruitfully employed to
analyze experimental data of impedance spectroscopy, as done in the last
decades [9]. Some special attention has to be devoted to the $\omega\to 0$
limit, since $0\leq\gamma\leq 1$. Consequently the condition on the solenoidal
character of ${\bf j}$ is violated in the dc limit.

In addition, when this condition is satisfied, a connection between the
impedance spectroscopy response of an anomalous Poisson-Nernst-Planck (PNP)
diffusional model and of equivalent circuits containing constant phase
elements (CPEs) may be established in general terms for a typical electrolytic
cell [10]. The analysis, carried out in the limit of low frequency in order to
highlight the surface effects, is shown to yield a natural presence of a CPE-
like behavior in the diffusional model as analogous to what is considered as a
tool used in the approaches based on equivalent electrical circuits.

A simple, but representative, example is in order here to illustrate this
important tool provided by the extensions of the PNP model involving
fractional operators. To explore it, we first define the electrolytic cell,
i.e., the region where the charges are diffusing, in the shap of a slab o
thickness $d$ in a Cartesian frame with the electrodes positioned at $z=\pm
d/2$, where $z$ the coordinate normal to the surfaces, which have area $A$.
This way, the problem becomes one-dimensional. For simplicity, we consider
also the case of full dissociation, and we assume that the ions have the same
diffusion coefficients. In this way, the equations to be solved are obtained
from Eqs. (Time-Fractional Approach to the Electrochemical Impedance: The
Displacement Current), rewritten as

$\displaystyle\int_{0}^{1}d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}n_{p}(z,t)+\frac{\partial}{\partial z}j_{p}(z,t)$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\int_{0}^{1}d\gamma
p(\gamma)\tau^{\gamma-1}\frac{\partial^{\gamma}}{\partial
t^{\gamma}}n_{m}(z,t)+\frac{\partial}{\partial z}j_{m}(z,t)$ $\displaystyle=$
$\displaystyle 0,$ (19)

together with the equation of Poisson, Eq. (2). The cell is subjected to a
time-dependent electric potential $V\left(\pm
d/2,t\right)=\pm\left(V_{0}/2\right)e^{i\omega t}$ and we limit our
considerations to the periodic-state. For perfect blocking electrodes, i.e.,
when the boundary conditions $j_{p}(\pm d/2,t)=j_{m}(\pm d/2,t)=0$ hold, the
electrical impedance spectroscopy may be investigated in the asymptotic limit
of small AC-signal limit. This condition permits one to search for solutions
of the fundamental equations, Eqs. (Time-Fractional Approach to the
Electrochemical Impedance: The Displacement Current) and (2) in the form

$n_{p}(z,t)=N+\eta_{p}(z)e^{i\omega t}\quad{\textrm{and}}\quad
n_{m}(z,t)=N+\eta_{m}(z)e^{i\omega t}$ (20)

with $N>>|\eta_{p}(z)|$ and $N>>|\eta_{m}(z)|$, where $N$ represents the
number of ions per unit volume in equilibrium, and $V(z,t)=\phi(z)e^{i\omega
t}$ (see Ref. [9], for details). We obtain for the electrical impedance the
following analytical expression [10]:

$\displaystyle
Z=\frac{2}{\Phi(i\omega)A\varepsilon\beta^{2}}\left\\{\frac{1}{\lambda^{2}\beta}\tanh\left(\beta
d/2\right)+\frac{d}{2{D}}\Phi(i\omega)\right\\},$ (21)

where $\beta=\left(1/\lambda\right)\sqrt{1+\Phi(i\omega)\lambda^{2}/D}$,
$\Phi(i\omega)=(1/\tau)\int_{0}^{1}d\gamma p(\gamma)(i\omega\tau)^{\gamma}$,
and $\lambda=\sqrt{\varepsilon k_{B}T/\left(2q^{2}N\right)}$ is the Debye
screening legnth. Figures 1, 2, and 3 show some prediction obtained from Eq.
(21) for different values of $\gamma$. We observe that in the low frequency
limit the fractional time derivative plays an important role and is
responsible for a different behavior from the normal one. It is characterized
by $Z\sim 1/(i\omega)^{\gamma}$, which has a direct connection with the
presence of constant phase elements (CPE) when the problem is formulated in
terms of equivalent circuits [10].

Figure 1: Nyquist diagram for different values of $\gamma$, with
$\Phi(i\omega)=(i\omega)^{\gamma}$. The black (solid), green (dotted), and
blue (dashed) lines correspond to the cases $\gamma=1.0$, $\gamma=0.9$, and
$\gamma=0.7$, respectively. We consider, for illustrative purposes,
$D=8.0\times 10^{-8}\,$m2/s, $A=3.1415\times 10^{-4}\,$m2,
$\varepsilon=90\varepsilon_{0}$ ($\varepsilon_{0}=8.85\times 10^{-12}\,$ F/m),
$d=1.33\times 10^{-3}\,$m, $\tau=1\,$s, and $\lambda=1.19\times 10^{-7}\,$m.
Figure 2: The imaginary part of the impedance, $X={\textrm{Im}}(Z)$, for
different values of $\gamma$, with $\Phi(i\omega)=(i\omega)^{\gamma}$. The
black (solid), green (dotted), and blue (dashed) lines correspond to the cases
$\gamma=1.0$, $\gamma=0.9$, and $\gamma=0.7$, respectively. The parameters are
the same as in Fig. 1.  Figure 3: The real part of the impedance,
$R={\textrm{Re}}(Z)$, for different values of $\gamma$, with
$\Phi(i\omega)=(i\omega)^{\gamma}$. The black (solid), green (dotted), and
blue (dashed) lines correspond to the cases $\gamma=1.0$, $\gamma=0.9$, and
$\gamma=0.7$, respectively. The parameters are the same as in Fig. 1.

This can be checked by considering in more detail the asymptotic behavior of
the impedance in the low frequency limit, which, for Eq. (21), is given by

$Z\approx\frac{\lambda^{2}d}{\varepsilon AD}+\frac{2\lambda}{\varepsilon
A}\frac{1}{\Phi(i\omega)}\;.$ (22)

The first addendum in Eq. (22) is directly related to the resistance of the
cell, it is a bulk contribution and it scales with the thickness of the cell,
$d$. The other term depends on frequency, contributing to the real and
imaginary parts of the impedance in the low frequency limit similar to the CPE
contribution to the total impedance of the cell. It takes origin from the
surface layer whose thickness is of the order of $\lambda$. In the case of
normal diffusion, it reduces to a surface capacitance. If, instead, the
diffusion is anomalous, the corresponding term takes into account also
dissipative effects, introducing a surface electrical resistance. In
particular, for the cases illustrated in Figs. 1, 2, and 3, where
$\Phi(i\omega)=(i\omega)^{\gamma}$, the contribution to the impedance of the
cell is
$1/(i\omega)^{\gamma}=\omega^{-\gamma}[\cos(\gamma\pi/2)-i\sin(\gamma\pi/2)]$
as predicted by the constant phase element [11, 12, 13]. The suitable
consideration for the displacement current lead us to a different behavior
from the one obtained by considering the standard form of the displacement
current [4], which implies in the following asymptotic limit

$\displaystyle Z\approx\frac{\lambda^{2}d}{i\omega\varepsilon
AD}\Phi(i\omega)+\frac{2\lambda}{i\omega\varepsilon A}\;,$ (23)

with a different behavior for the impedance. In particular, the behavior
related to the constant phase element present in Eq. (22) is absent in Eq.
(23) for the imaginary part of the impedance in the very low frequency limit.

To sum up, we have analyzed the role of the displacement current in a
generalized time-fractional approach to the PNP model for electrical
impedance. We have shown that, in order to keep the concept of electrical
impedance meaningful, the total electric current arising in this extended
model has to be solenoidal. This requirement implies that also the expression
of the displacement current has to be modified by extending it to the
fractional domain too. The proposed form of the displacement current coincides
with the one used recently in some generalizations of the Maxwell equations
for material media involving fractional order operators. The time-fractional
extension of PNP model naturally exhibits a CPE-like behavior of the
electrical impedance in the low frequency region of the spectra and is shown
to be a powerful tool to face a wide variety of complex behaviour found in
electrolytic cells. From the analysis reported above, we conclude that the CPE
behaviour observed in the frequency dependence of the impedance of
electrolytic systems, can have different origin [13]. However, in the case of
media presenting a porous structure, it could be related to the anomalous
diffusion of the ions in it.

Acknowledgements. This work was supported by the MEPhI Academic Excellence
Project (G. B.) and by the Program of Visiting Professor of Politecnico di
Torino (L. R. E.). E.K.L. thanks partial financial support of the CNPq under
Grant No. $302983/2018-0$.

## References

  * [1] J. R. Macdonald, Theory of ac space–charge polarization effects in photoconductors, Semiconductors, and Electrolytes Phys. Rev. 92, 4–17 (1953).
  * [2] E. M. Trukhan, Dispersion of the dielectric constant of heterogeneous systems, Sov. Phys. Solid State (Engl. Transl.) 4, 2560 (1963).
  * [3] I. Lelidis and G. Barbero, Role of the displacement current on Warburg-type behavior.Phys. Rev. E 95,052604 (2017).
  * [4] E. K. Lenzi, L. R. Evangelista, and G. Barbero, Fractional diffusion equation and impedance spectroscopy of electrolytic cells, J. Phys. Chem. B 113, 33, 11371?11374 (2009).
  * [5] A. V. Chechkin, R. Gorenflo, I. M. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Phys. Rev. E 66, 046129 (2002).
  * [6] F. Mainardi, G. Pagnini, and R. Gorenflo, Some aspects of fractional diffusion equations of single and distributed order, App. Math. Comp. 187, 295–305 (2007).
  * [7] A. N. Bogolyubov, A. A. Potapov, and S. S. Rehviashvili, An approach to introducing fractional integro-differentiation in classical electrodynamics, Mosc. Univ. Phys. Bull. 64, 365?368 (2009).
  * [8] T. P. Stefański and J. Gulgowski, Formulation of time-fractional electrodynamics based on Riemann-Silberstein vector, Entropy 2021, 23, 987.
  * [9] L. R. Evangelista and E. K. Lenzi, Fractional diffusion dquations and anomalous diffusion (Cambridge University Press, Cambridge, 2018).
  * [10] E. K. Lenzi, J. L. de Paula, F. R. G. B. Silva, and L. R. Evangelista, A connection between anomalous Poisson-Nernst-Planck model and equivalent circuits with constant phase elements, J. Phys. Chem. C 117, 117, 23685?23690 (2013).
  * [11] L. Nyikos and T. Pajkossy, Fractal dimension and fractional power frequency-dependent impedance of blocking electrodes, Electrochim. Acta 30, 1533–1540 (1985).
  * [12] H. S. Liu, Fractal model for ac response of a rough interface, Phys. Rev. Lett. 55, 529–532 (1985).
  * [13] J. Jorcin, M. E. Orazem, N. Pébére, and B. Tribollet, CPE analysis by local electrochemical impedance spectroscopy, Electrochim. Acta 51, 1473–1479 (2006).