# Issues with $J$-dependence in the LSDA$+U$ method for non-collinear magnets

Eric Bousquet1,2 Nicola Spaldin3 1Materials Department, University of
California, Santa Barbara, CA 93106, USA  2Physique Théorique des Matériaux,
Université de Liège, B-4000 Sart Tilman, Belgium  3 Department of Materials,
ETH Zurich, Wolfgang-Pauli-Strasse 10 CH-8093 Zurich, Switzerland

###### Abstract

We re-examine the commonly used density functional theory plus Hubbard U
(DFT$+U$) method for the case of non-collinear magnets. While many studies
neglect to explicitly include the exchange correction parameter J, or consider
its exact value to be unimportant, here we show that in the case of non-
collinear magnetism calculations the J parameter can strongly affect the
magnetic ground state. We illustrate the strong J-dependence of magnetic
canting and magnetocrystalline anisotropy by calculating trends in the
magnetic lithium orthophosphate family LiMPO4 (M = Fe and Ni) and difluorite
family MF2 (M = Mn, Fe, Co and Ni). Our results can be readily understood by
expanding the usual DFT$+U$ equations within the spinor scheme, in which the J
parameter acts directly on the off-diagonal components which determine the
spin canting.

first-principles, LDA+U, non-collinear magnetism, magnetocrystalline
anisotropy

Density functional theory (DFT) within the local density (LDA) and generalized
gradient (GGA) approximations is widely used to describe a large variety of
materials with good accuracy. The LDA and GGA functionals often fail, however,
to correctly reproduce the properties of strongly correlated materials
containing d and f electrons. The LDA$+U$ approach – in which a Hubbard U
repulsion term is added to the LDA functional for selected orbitals – was
introduced in response to this problem, and often improves drastically over
the LDA or GGA. Indeed, it provides a good description of the electronic
properties of a range of exotic magnetic materials, such as the Mott insulator
KCuF3liechtenstein1995  and the metallic oxide LaNiO2 lee2004 .

Two main LDA$+U$ schemes are in widespread use today: The Dudarev dudarev1998
approach in which an isotropic screened on-site Coulomb interaction
$U_{eff}=U-J$ is added, and the Liechtenstein liechtenstein1995  approach in
which the $U$ and exchange ($J$) parameters are treated separately. The
Dudarev approach is equivalent to the Liechtenstein approach with $J=0$
baettig2005 . Both the effect of the choice of LDA+$U$ scheme on the orbital
occupation and subsequent properties ylvisaker2009 , as well as the dependence
of the magnetic properties on the value of $U$ savrasov2005 , have recently
been analyzed. There has been no previous systematic study, however, of the
effect of the $J$ parameter of the Liechtenstein approach in non-collinear
magnetic materials. Here we show that neither the approach of not explicitly
considering the $J$ parameter (as in the Dudarev implementation), nor the
assumption that its importance is borderline – a common approximation is to
use $J\simeq 10\%\ U$ without careful testing – within the Liechtenstein
implementation are justified in the case of non-collinear magnets. We
demonstrate that in the case of non-collinear antiferromagnets, the choice of
$J$ can strongly change the amplitude of the spin canting angle (LiNiPO4) or
even modify the easy axis of the system (LiFePO4 and FeF2), with consequent
drastic effects on the magnetic susceptibilities and magnetoelectric
responses.

First we remind the reader how the $U$ and $J$ parameters appear in the usual
collinear spin LSDA$+U$ formalism. The LSDA$+U$ reformulation of the LSDA
Hamiltonian is usually written as:

$\displaystyle H_{LSDA+U}=H_{LSDA}+H_{U}\quad,$ (1)

whith

$\displaystyle
H_{U}^{\sigma}=\displaystyle\sum_{m_{1},m_{2}}P_{m_{1},m_{2}}V^{\sigma}_{m_{2},m_{1}}\quad,$
(2)

where $P$ is the projection operator, $\sigma$ is the spin index, and (on a
given atomic site):

$\displaystyle V^{\uparrow(\downarrow)}_{m_{2},m_{1}}=$
$\displaystyle\displaystyle\sum_{3,4}\left(V^{ee}_{1,3,2,4}-U\delta_{1,2}-V^{ee}_{1,3,4,2}+J\delta_{1,2}\right)n^{\uparrow(\downarrow)}_{3,4}$
$\displaystyle+\left(V^{ee}_{1,3,2,4}-U\delta_{1,2}\right)n^{\downarrow(\uparrow)}_{3,4}+\frac{1}{2}(U-J)\delta_{1,2}$
(3)

Here $V^{ee}_{1,3,2,4}=\left\langle
m_{1},m_{3}\left|V^{ee}_{m_{1},m_{3},m_{2},m_{4}}\right|m_{2},m_{4}\right\rangle$
are the elements of the screened Coulomb interaction (which can be viewed as
the sum of Hartree (direct) contributions $V^{ee}_{1,3,2,4}$ and Fock
(exchange) contributions $V^{ee}_{1,3,4,2}$ and $n^{\sigma}_{i,j}$ are the
$d$-orbital occupancies.

In the case of non-collinear magnetism, the formalism is extended and the
density is expressed in a two-component spinor formulation:

$\displaystyle\rho=$
$\displaystyle\begin{pmatrix}\rho^{\uparrow\uparrow}&\rho^{\uparrow\downarrow}\\\
\rho^{\downarrow\uparrow}&\rho^{\downarrow\downarrow}\end{pmatrix}=\begin{pmatrix}n+m_{z}&m_{x}-im_{y}\\\
m_{x}+im_{y}&n-m_{z}\end{pmatrix}$ (4)

where $n$ is the charge density and $m_{\alpha}$ the magnetization density
along the $\alpha$ direction ($\alpha=x,y,z$). Using the double-counting
proposed by Bultmark et al.bultmark2009 , the LSDA$+U$ potential is then also
expressed in the two-component spin space as:

$\displaystyle V_{i,j}=$
$\displaystyle\begin{pmatrix}V^{\uparrow\uparrow}_{i,j}&V^{\uparrow\downarrow}_{i,j}\\\
V^{\downarrow\uparrow}_{i,j}&V^{\downarrow\downarrow}_{i,j}\end{pmatrix}$ (5)

where $V^{\uparrow\uparrow}$ and $V^{\downarrow\downarrow}$ are equal to Eqs.3
and

$\displaystyle
V^{\uparrow\downarrow(\downarrow\uparrow)}_{m_{2},m_{1}}=\displaystyle\sum_{3,4}\left(-V^{ee}_{1,3,4,2}+J\delta_{1,2}\right)n^{\uparrow\downarrow(\downarrow\uparrow)}_{3,4}$
(6)

For collinear magnets, only $V^{\uparrow\uparrow}$ and
$V^{\downarrow\downarrow}$ (Eqs. 3) are relevant since
$n^{\uparrow\downarrow}$ and $n^{\downarrow\uparrow}$ are equal to zero, and
$J$ affects the potential mainly through an effective $U-J$. However, in the
case of non-collinear magnetism, the $n^{\uparrow\downarrow}$ and
$n^{\downarrow\uparrow}$ and hence the $V^{\uparrow\downarrow}$ and
$V^{\downarrow\uparrow}$ (Eqs. 6) are non-zero. Then it is clear from Eqs. 6
that $J$ acts explicitly on the off-diagonal potential components.

Next, we show the effect of the choice of $J$ parameter in the family of
lithium orthophosphates, LiMPO4 (M = Ni and Fe) and in the family of
difluorites MF2 (M = Mn, Co, Fe and Ni). The orthophosphates crystallize in
the orthorhombic Pnma space group with $C$-type antiferromagnetic (AFM) order.
The difluorites crystalize in the tetragonal P42/mnm rutile structure with AFM
order. We performed calculations within the Liechtenstein approach of the
DFT$+U$ as implemented in the VASP code vasp1 ; vasp2  111We note that
LSDA$+U$ double-couting term taking into accound the magnetization density as
proposed by Bultmark et al.bultmark2009  is mandatory within non-collinear
magnetism calculations. This is not necessarily done in the present
implementation of other codes. with $U$ and $J$ corrections applied to the 3d
orbitals of the M cations. In all cases we relaxed the atomic positions until
the residual forces on each atom were lower than 10 $\mu$eV/Å at the
experimental volume and cell shape reported in Tab. 1, taking into account the
spin-orbit interaction. We found good convergence of the non-collinear spin
ground state with a cutoff energy of 500 eV on the plane wave expansion and a
k-point grid of $2\times 4\times 4$ for the orthophophates and $4\times
4\times 6$ for the difluorites.

| $a$ | $b$ | $c$ | Ref.
---|---|---|---|---
LiFePO4 | 10.332 | 6.010 | 4.692 | streltsov1993,
LiNiPO4 | 10.032 | 5.854 | 4.677 | abrahams1993,
NiF2 | 4.650 | 4.650 | 3.084 | hutchings1970,
FeF2 | 4.700 | 4.700 | 3.310 | dealmeida1989,
MnF2 | 4.650 | 4.650 | 3.084 | oguchi1958,
CoF2 | 4.695 | 4.695 | 3.179 | otoole2001,
Table 1: Experimental cell parameters (Å) used in the simulations of LiMPO4
phosphates and MF2 difluorites.

Figure 1: (a) Calculated LSDA$+U$ canting angle of LiNiPO4 versus $J$ for
$U=5$ eV. The experimental value of the canting angle is equal to 7.8∘
jensen2009 . (b) Energy versus canting angle in LiNiPO4 for $U=5$ eV and $J=0$
eV (red circles), $U=5$ eV and $J=1$ eV (blue triangles), $U_{eff}=4$ eV
(green crosses) and $U=5$ eV and $J=1$ eV but by fixing $J=$ 0 eV in Eqs.6
(pink squares). The zero energy reference is chosen at zero canting angle. (c)
Magnetocrystaline anisotropy energy (MCAE) between the $a$ and $b$
orientations of the magnetic moments of LiFePO4. The experimental $b$
orientations is taken as energy reference.

First, we focus on LiNiPO4, which is known experimentally to be $C$-type AFM,
with an easy-axis along the $c$ direction and a small A-type AFM canting of
the spins along the $a$ direction ($C_{z}A_{x}$ ground state with mm’m
magnetic point group) jensen2009 . Performing calculations within the LSDA$+U$
method with $J=0$, we find that we correctly reproduce the $C_{z}A_{x}$ ground
state with a rather small $U$ sensitivity of the magnetocrystalline anisotropy
energy (MCAE) and the spin canting; this finding is consistent with a previous
report using the GGA functional yamauchi2010 . However, our calculated canting
angle of 1.6∘ for $U=5$ eV and $J=0$ eV severely underestimates the
experimental value of 7.8∘jensen2009 . In Fig.1 (a) we show the evolution of
the canting angle with $J$ at $U=5$ eV. We find that the canting angle is
extremely sensitive to the value of $J$ – in fact it is $\propto J^{3}$ –
changing from 1.6∘ at $J=0$ eV to 7.8∘ at $J=1.7$ eV. To reproduce the
experimental value of the canting angle we need to use the rather large $J$
value of 1.7 eV. The dependence of the canting angle on $J$ is consistent with
Eqs. 6, as the off-diagonal elements $n^{\uparrow\downarrow}$ and
$n^{\downarrow\uparrow}$ are non-zero when the spins cant away from the easy
axis.

In Fig. 1 (b) we report the energy versus the canting angle in LiNiPO4 for
$U=5$ eV and different values of $J$. We see that as $J$ is increased from
$J=0$ eV to $J=1$ eV (red circles and blue triangles) the minimum of the
energy shifts to larger canting angle, with a stronger gain of energy with
respect to the uncanted reference. When performing the same calculation with
$U_{eff}=$4 eV (green crosses in Fig. 1) we obtain results that are very
similar to the case $U=5$ eV and $J=0$ eV, which is formally equivalent to the
Dudarev approach with $U_{eff}=5$ eV. These comparisons confirm that varying
$U$ has a minimal effect on the canting angle in LiNiPO4 and also that the use
of the Liechtenstein treatment of $J$ is extremly important. To further
confirm the direct relationship between the spin canting and the $J$
parameter, we performed the same calculations with $U=5$ eV and $J=1$ eV but
we artificially fixed $J=0$ eV only in Eqs. 6 (pink squares in Fig.1 (b)). We
clearly see that the energy versus canting angle is strongly affected by this
modification and in fact the canting is almost removed.

Similar $J$ dependence of the canting angle was also reported previously for
Ni2+ in BaNiF4 ederer2006 ; in Ref. ederer2006,  it was found that at $U=5$
eV, the canting varies from 2∘ to 3∘ when $J$ is varied from 0 eV to 1 eV. In
both LiNiPO4 and BaNiF4 the Ni ion is divalent, with a $d^{8}$ configuration,
and octahedrally coordinated. To investigate the generality of this behavior,
we next consider the case of the canted-spin antiferromagnet NiF2, in which
the Ni ion is in the same coordination environnement as in BaNiF4.
Experimentally, NiF2 has the spins aligned preferentially in the plane
perpendicular to the c axis with a slight canting from antiparallel alignment
by an estimated $\sim$0.5∘ at low temperatures hutchings1970 . Performing
LSDA$+U$ calculations at the experimental volume and with $U=5$ eV and $J=0$
eV we indeed obtain the easy axis perpendicular to the c axis and a small
canting of 0.3∘, in excellent agreement with the experiments. In contrast to
the case of LiNiPO4, however, we find that the amplitude of the canting angle
is almost insensitive to the value of $J$ with just a small tendency to be
reduced when $J$ increased. This insensitivity of the canting angle to the
value of $J$ in NiF2 can be understood from the fact that in this compound the
magnetism is almost collinear, and therefore the off-diagonal elements of the
occupation matrix, $n^{\uparrow\downarrow}$ and $n^{\downarrow\uparrow}$, are
close to zero. Inspection of Eqs. 3 then shows that the effect of $J$ is
reduced largely to the diagonal part of the potential where the $U$ parameter
is dominant.

To summarize our findings for the Ni-based compounds, in cases where the
experimental canting is large (2-3∘) we find a strong $J$-dependence of the
canting angle, which increases with increasing $J$; when the canting is weak
experimentally the $J$-dependence is much weaker.

Figure 2: Magnetocrystaline anisotropy energy versus the $J$ parameter of (a)
FeF2 (Experimental value from Ref.lines1967, ), (b) NiF2, (c) MnF2
(Experimental value from Ref.gafver1977, ) and (d) CoF2 (“sc“ are calculations
with Co semi-cores while ”no sc” are calculations without Co semi-cores). The
MCAE reported here is the energy between the $a$ and $c$ orientation of the
spins, the energy of the $c$ orientation is taken as reference.

Next we analyse the effect of $J$ on the behavior on the corresponding
divalent iron compounds. We begin with LiFePO4, which is known experimentally
to be a $C$-type AFM with an easy axis along the $b$ direction and no observed
canting of the spins zimmermann2009 ; liang2008  ($C_{y}$ ground state with
mmm’ magnetic point group). Our calculations within the LSDA$+U$ functional at
the commonly used values of $U=4$ eV and $J=0$ eV for Fe2+ yield the correct
$C$-type AFM order but find the easy axis incorrectly along the $a$ direction.
Now we switch to $J\neq 0$ eV and report in Fig. 1.c the MCAE between the $b$
and $a$ directions, calculated by turning all the spins homogenously from the
$C_{y}$ to the $C_{x}$ direction. We find that the MCAE is approximately
linear with $J$, but with rather dramatic qualitative dependence: while at
$J=0$ eV the easy axis is along the $a$ direction (negative MCAE) the MCAE is
almost reduced to zero around $J=0.5$ eV and the easy axis changes to the $b$
direction for $J\gtrsim 0.5$ eV (positive MAE). To reproduce the experimental
easy axis ($C_{y}$) a value of $J$ greater than 0.58 eV is required. In the
cases where the correct easy axis is reproduced ($C_{y}$) we do not observe
any canting of the spins, in agreement with the experimental magnetic point
group mmm’.

As a second example with Fe2+, we analyse the effect of $J$ on the MCAE of
FeF2. Experimentally FeF2 is known to have its spin magnetization parallel to
the tetragonal $c$ axis with a rather large MCAE of about +4800 $\mu$eV
rudowicz1977 ; ohlmann1961 . In Fig. 2.a we report the LSDA$+U$ MCAE energies
with respect to $J$ at four different values of $U$ (3, 4, 5 and 6 eV). All
the calculations with $J=0$ eV give the wrong easy axis (spins are
perpendicular to $c$) with a huge error in the MCA energy (MCAE from -16000 to
-26000 $\mu$eV for $U$ going from 3 to 6 eV). Increasing the value of $J$ in
the range of 0–0.5 eV has the tendency to strongly reduce this error with a
linear increase of the MCAE with $J$ as we found above for LiFePO4. However
beyond $J\simeq 0.5$ the increase of the MCAE is reduced and the evolution
becomes more complex with the appearance of two maxima before a drastic
decrease beyond $J\simeq 1.3$ eV. The correct easy axis (MCAE$>0$) is only
obtained for a very small range of $U$ and $J$ values, and the amplitude of
the MCAE is correct over an even smaller range. This $J$ dependence of the
MCAE is again consistent with Eqs.3-6. From Eq.4 it is clear that when
changing the orientation of the spins from the $z$ axis to the $x$ or $y$ axis
the off-diagonal parts of Eq.4 become non-zero resulting in a direct effect of
$J$ on the MCAE from Eqs.6.

We also performed the same analysis of the MCAE for NiF2 (Fig.2.b), MnF2
(Fig.2.c) and CoF2 (Fig.2.d). MnF2 and CoF2 have the same easy axis as FeF2
while NiF2 has its easy axis perpendicular to the $c$ direction. The easy axis
is well reproduced for all three compounds at $J=0$ eV. As for FeF2, the
amplitudes of the MCAE depend strongly on $J$ but with a completely different
trend in each compound. For MnF2 and FeF2 the experimental value can be
reproduced by adjusting the values of $U$ and $J$. In the case of CoF2 and
NiF2 no experimental values are available. For CoF2 we also performed
calculations with and without Co semi-cores states (Fig.2.d) and find a strong
difference in the magnitude of the MCAE for the two cases. For FeF2 we also
performed calculations within the GGA functional (black pentagons in Fig.2.a)
and obtained a completely different $J$ dependence than those calculated with
the LDA functional. These comparisons illustrate the difficulty of extracting
a general rule about the $J$ dependence of the MCAE.

Our results reveal a problem with the predictability of the LSDA$+U$ method
for non-collinear magnetic materials: A strong dependence of the MCAE and spin
canting angles on the values of $U$ and particularly $J$ that are used in the
calculation. Since properties such as magnetostriction, piezomagnetic
response, magnetoelectric response and exchange bias coupling are directly
related to MCAEs and spin canting, it is of primary importance to reproduce
these quantities accurately. At the moment, the most reliable, although not
entirely satisfactory, option appears to be a fine tuning of the $U$ and $J$
parameters by adjustment to reproduce experimentally measured anisotropies and
canting angles; there is some evidence to suggest that properties such as
magnetoelectric responses are then in turn well reproduceddelaney2010 . Future
studies might explore methodologies for self-consistent calculation of the $J$
parameter, or the predictions of new descriptions of the exchange and
correlation such as the hybrid functionals heyd2004 . On the flip side, it is
clear that non-collinear magnetic systems provide a challenging case for
testing the correctness of new exchange correlation functionals within the
density functional formalism.

Acknowledgments: This work was supported by the Department of Energy SciDAC
DE-FC02-06ER25794. We made use of computing facilities of TeraGrid at the
National Cen-ter for Supercomputer Applications and of the California
Nanosystems Institute with facilities provided by NSF grant No. CHE-0321368
and Hewlett-Packard. EB also acknowledges FRS-FNRS Belgium and the ULg SEGI
supercomputer facilities.

## References

  * (1) A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (Aug 1995)
  * (2) K.-W. Lee and W. E. Pickett, Phys. Rev. B 70, 165109 (Oct 2004)
  * (3) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (Jan 1998)
  * (4) P. Baettig, C. Ederer, and N. A. Spaldin, Phys. Rev. B 72, 214105 (2005)
  * (5) E. R. Ylvisaker, W. E. Pickett, and K. Koepernik, Phys. Rev. B 79, 035103 (Jan 2009)
  * (6) S. Y. Savrasov, A. Toropova, Kat, K. M. I., L. A. I., A. V., and G. Kotliar, Z. Kristallogr. 220, 473 (2005)
  * (7) F. Bultmark, F. Cricchio, O. Grånäs, and L. Nordström, Phys. Rev. B 80, 035121 (2009)
  * (8) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (Oct 1996)
  * (9) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (Jan 1999)
  * (10) We note that LSDA$+U$ double-couting term taking into accound the magnetization density as proposed by Bultmark et al.bultmark2009  is mandatory within non-collinear magnetism calculations. This is not necessarily done in the present implementation of other codes.
  * (11) V. A. Streltsov, E. L. Belokoneva, V. G. Tsirelson, and N. K. Hansen, Acta Cryst. B 49, 147 (Apr 1993)
  * (12) I. Abrahams and K. S. Easson, Acta Cryst. C 49, 925 (1993)
  * (13) M. T. Hutchings, M. F. Thorpe, R. J. Birgeneau, P. A. Fleury, and H. J. Guggenheim, Phys. Rev. B 2, 1362 (Sep 1970)
  * (14) M. J. M. de Almeida, M. M. R. Costa, and J. A. Paixão, Acta Cryst. B 45, 549 (1989)
  * (15) T. Oguchi, Phys. Rev. 111, 1063 (Aug 1958)
  * (16) N. J. O’Toole and V. A. Streltsov, Acta Cryst. B 57, 128 (2001)
  * (17) T. B. S. Jensen, N. B. Christensen, M. Kenzelmann, H. M. Rønnow, C. Niedermayer, N. H. Andersen, K. Lefmann, J. Schefer, M. v. Zimmermann, J. Li, J. L. Zarestky, and D. Vaknin, Phys. Rev. B 79, 092412 (2009)
  * (18) K. Yamauchi and S. Picozzi, Phys. Rev. B 81, 024110 (2010)
  * (19) C. Ederer and N. A. Spaldin, Phys. Rev. B 74, 020401 (Jul 2006)
  * (20) M. E. Lines, Phys. Rev. 156, 543 (Apr 1967)
  * (21) U. Gäfvert, L. Lundgren, P. Nordblad, B. Westerstrandh, and . Beckman, Sol. State Comm. 23, 9 (1977)
  * (22) Zimmermann, A. S., Van Aken, B. B., Schmid, H., Rivera, J.-P., Li, J., Vaknin, D., and Fiebig, M., Eur. Phys. J. B 71, 355 (2009)
  * (23) G. Liang, K. Park, J. Li, R. E. Benson, D. Vaknin, J. T. Markert, and M. C. Croft, Phys. Rev. B 77, 064414 (2008)
  * (24) C. Rudowicz, J. Phys. Chem. Solids 38, 1243 (1977)
  * (25) R. C. Ohlmann and M. Tinkham, Phys. Rev. 123, 425 (1961)
  * (26) K. Delaney, E. Bousquet, and N. A. Spaldin, arXiv:0912.1335v2(2010)
  * (27) J. Heyd and G. E. Scuseria, The Journal of Chemical Physics 120, 7274 (2004)


because the latter is supposed to be completely screened at that depth. The
situation is right opposite in the MW model: the greater $H_{0}$, the greater
$H(=H_{0}$ in the given geometry) and therefore the greater the induced
currents $J_{i}$ circulating all over the sample volume. Correspondingly,
inhomogeneity of the microscopic fields increases with increasing $H$. If so,
the depolarization rate of muons stopped at the depths where $B=0$ should be
increasing with increasing $H_{0}$. This is exactly what we see in Figs. 13a
and b, and 14.

## SUMMARY AND OUTLOOK

Historical experiments of Meissner and Ochsenfeld and of Rjabinin and
Shubnikov, in which the Meissner effect was established, are reviewed.
Theories of Gorter and Casimir and of F. and H. London addressing,
respectively, thermodynamic and electromagnetic properties of superconductors
in the Meissner state are reviewed as well. It is shown that due to
significant difficulties in Londons’ theory, the standard description of the
electromagnetic properties of superconductors is inadequate and requires
revision.

A novel semi-classical micro-whirls model targeted to address this issue is
developed and presented. The model is based on the concept of paired electrons
obeying the Bohr-Sommerfeld quantization condition. Accordingly, each Cooper
pair represents a microscopic circular current like the current due to
orbiting electrons in atoms.

The model is valid in the whole range of the superconducting state for samples
of any shape and at any orientation of the applied field. A fundamentally new
prediction of the model is that the S phase of superconducting ellipsoidal
bodies in a static magnetic field is filled with ordered microscopic whirls of
field induced currents caused by precessing Cooper pairs. The model is free
from disadvantages of the standard theories, simultaneously reproducing their
achievements. All properties predicting by the model can be verified
experimentally.

The model consistently describes equilibrium properties of superconductors,
which includes the Meissner effect, persistency of induced currents, zero
entropy of “superconducting” electrons (i.e. of the ensemble of Cooper pairs),
the London rigidity principle for Cooper pairs, the flux quantization, the S/N
interface energy of two signs and hence two type of superconductivity, and
others. Some well-known experimental facts, such as the Meissner state in non-
spheroidal samples and paramagnetism of the Abrikosov vortices, are
consecutively described in the new model for the first time. Non-equilibrium
properties of superconductors, e.g., the non-existence of the Hall effect and
other properties associated with the total current are explained by the model
as well.

According to the new model, there are two kinds of non-dissipating currents in
superconductors. (1) The current formed by electrons bound in stationary
Cooper pairs; it arises due to precession of the magnetic moments of the pairs
caused by the orbital motion of the coupled electrons about their center of
mass; these microscopic currents dictate thermodynamic (equilibrium)
properties of the singly connected isolated ellipsoidal samples, i.e.,
properties of the Meissner, intermedium and mixed states. (2) A dissipation-
free flow of Cooper pairs forming the total current leading to the non-
equilibrium properties resembling the properties of hypothetical perfect
conductors. Both equilibrium and non-equilibrium properties follow from the
single Bohr-Sommerfeld quantization condition for the paired electrons (Eq.
(16)).

There was a long lasting discussion about what is the prime property of
superconductors, i.e. of the superconducting fraction of conduction electrons
or the ensemble of Cooper pairs, either it is zero resistivity or zero
induction Landafshitz_II ; DeGennes ; Shoenberg . For some reason, this
discussion often overlooks the third fundamental property postulated by Gorter
and Casimir: zero entropy Gorter_Casimir . In the micro-whirls model these
three key characteristics of superconductivity are on an equal footing and
share the same root: quantization of the orbital motion of the paired
electrons.

Contrarily to the standard theories, the model predicts a strong field-
dependent inhomogeneity of the induced microscopic currents in the bulk of the
samples in the Meissner state. The existence of such currents follows from the
requirement of the First law of thermodynamics and were foreseen by Hall Hall
coming from the disappearance of the Hall effect first observed by Onnes and
Hof Onnes-1914 , and by I. Kikoin Isaak-2  based on his and Goobar
measurements of the gyromagnetic ratio Isaak-1 . Available LE-$\mu$SR data
support this prediction. However these data were obtained in the experiment
aimed to another task and therefore require a targeted verification. It can be
also interesting to perform the SANS experiment with samples in the Meissner
state. Potentially, such an experiment can make it possible to measure the
diameter of micro-whirls.

One more way to verify the model can be as follows. According to the model,
the sample in the Meissner state is filled with the ordered and parallel to H
whirls of the field-induced microscopic circular currents with a frequency
$\omega=eH/2mc$, where H is the field intensity inside the sample (Eq. (30)).
Therefore, if it is so, a weak ac transverse field (with respect to H) should
experience a resonance absorption at the linear frequency
$\nu_{s}=\omega/2\pi=1.8H$ MHz, where $H$ is in Oe (numerically the same as
G). The frequency $\nu_{s}$ equals half the standard frequency of the electron
paramagnetic resonance (EPR). Important to stress that the sample should not
be small. For example, it can be a piece of straight superconducting wire with
a diameter of 0.5 mm and a length of 5-10 mm or longer. In such case
$\textbf{H}=\textbf{H}_{0}$ if the applied field $\textbf{H}_{0}$ is parallel
to the sample longitudinal axis, and $\textbf{H}=2\textbf{H}_{0}$, if
$\textbf{H}_{0}$ is perpendicular to this axis. For a spherical sample (e.g.,
4 mm in diameter like the sample used in the experiment of I. Kikoin and
Goobar) $\textbf{H}=3\textbf{H}_{0}/2$. To ensure unambiguity, the sample
should be pure and free of frozen-in flux.

Note that a possible existence of such a diamagnetic resonance can be traced
from measurements of the gyromagnetic ratio in superconducting samples by I.
Kikoin and Goobar Isaak-1 ; Isaak-2  similar as the existence of EPR is traced
from measurements of the gyromagnetic ratio in ferromagnetic samples by
Einstein and de Haas636363In both works the gyromagnetic ratio was measured
using a resonance technique, in which the sample oscillates in resonance with
the oscillating applied field. Einstein .

Regarding the model itself, its story is not complete yet, of course. One of
remaining questions is the penetration depth at the surface of the samples in
the MS and at the S/N interphase boundaries within the samples in the
intermediate and the mixed states.

Perhaps the main feature of the new model is its striking simplicity. As cited
above, a model of such kind was anticipated by Fritz London based on his
discovery of the rigidity principle. On the other hand, being semi-classical,
the model can not and does not address specifics of the lattice polarization
responsible for Cooper pairing.

Finishing, we cite one more excerption from the book of F. London London50 :
”Thus we may conclude that it ought to be sufficient [for construction of the
model] to derive the rigidity of the average momentum for simply connected
isolated superconductor”. It seems we have strong reason to say that the
presented model copes with this task.

## APPENDIX

Thermodynamic argument of H. London.

In 1947 Heinz London suggested a thermodynamic argument according to which the
London penetration depth $\lambda_{L}$ should depend on the applied field
$H_{0}$ HLondon47 ; Pippard50 ; Shoenberg . This argument has played an
important role in the history of superconductivity, but is rarely mentioned in
textbooks. Let us recall and discuss it.

Consider a long superconducting (for definiteness type-I) rod of radius
$\Re\gg\lambda_{L}$ in a parallel field $\textbf{H}_{0}$, i.e., a cylindrical
sample in the MS. A strict thermodynamic relationship reads for this
sample646464In general case, i.e., for any ellipsoidal sample in an arbitrary
oriented field this relationship is
$(\nabla_{\textbf{H}}S)_{T}=(\partial\textbf{M}/\partial T)_{H}$.  (see, e.g.,
VK )

$\left(\frac{\partial S}{\partial H_{0}}\right)_{T}=\left(\frac{\partial
M}{\partial T}\right)_{H_{0}},$ (59)

where $S$ and $M$ are the entropy and magnetic moment of the sample,
respectively, and $T$ is its temperature.

In the London theory the magnetic moment of the chosen sample is (see footnote
(20) above)

$M=-\frac{H_{0}}{4\pi}(V-V_{p})=-\frac{H_{0}}{4\pi}(V-\lambda_{L}A_{sm}),$
(60)

where $V_{p}$ and $A_{sm}$ stand for the volume of penetration layer and the
sample surface area, respectively.

Note that $V_{p}$ represents a kind of “excluded” volume possessing pretty
peculiar properties: it does not contribute into the sample magnetic moment,
but the current in this very volume sets this moment.

Next, neglecting thermal expansion and taking into account that $\lambda_{L}$
depends on $T$, from Eq. (59) it follows

$\left(\frac{\partial S}{\partial
H_{0}}\right)_{T}=\frac{A_{sm}H_{0}}{4\pi}\left(\frac{\partial\lambda_{L}}{\partial
T}\right)_{H_{0}}.$ (61)

Therefore, since $(\partial\lambda_{L}/\partial T)_{H_{0}}\neq 0$ Shoenberg ,
the sample entropy $S$ depends on $H_{0}$ at constant temperature. Important
to point out that this entropy is contained entirely within the penetration
layer656565This is another strange property of $V_{p}$, which is excluded from
the production of magnetic moment. (if $V_{p}=0$ then there is no the entropy
change with the field) and its density (calculated as $(S-S_{0})/V_{p}$, where
$S_{0}$ is the entropy at zero field) is  not small Pippard50 ; London50 ;
Shoenberg .

According to the two-fluid model, the entropy of a superconductor in the MS is
caused by the change of $n_{s}$666666Actually, as we have seen, this is not
quite a lawful extension of the two-fluid model, where $n_{s}$ changes only
with temperature.. Hence, the field dependence of $S$ implies that $n_{s}$
depends on the field at constant temperature. Therefore, $\lambda_{L}(\sim
1/\surd n_{s})$ depends on the field. This is the essence of H. London’s
argument. Recall that the field-independence of $\lambda_{L}$ is one of the
base assumptions of the London theory. So, H. London’s argument poses one more
dilemma in this theory.

This argument was discussed by Pippard Pippard50  who, as mentioned above,
found experimentally that $\lambda_{L}$ is essentially field-independent.
Pippard’s discussion is reproduced in London50  and Shoenberg , therefore here
we note only two points. (i) Pippard dropped from consideration 70% of his
data676767Those are the data at 1.5 K $\lesssim T\lesssim$ 3 K. At these
temperatures the reported variation of the effective penetration depth
$\delta\lambda_{L}=\Delta\lambda_{L}/\lambda_{L}$ at $H_{0}$ changing from
zero to $H_{c}$ decreases from about 0.023 at 1.5 K to 0.002 at 3 K. Above 3 K
$\delta\lambda_{L}$ increases up to 0.027 near $T_{c}$. From Eq. (61) (after
two more steps) it follows that $\delta\lambda_{L}$ increases with increasing
temperature starting from zero at 0 K. On that ground Pippard concentrated on
the data at $T>$ 3 K. thus admitting, in fact, that error bars he provides are
not fully adequate. Indeed, as it was demonstrated by Sridhar and Mercereau
Sridhar , Pippard’s data can be subject to non-equilibrium effects caused by
high-frequency radiation, which have not been accounted. (ii) Basing on the
remaining data and (although not always consequent) interpretation of Eq. (61)
Pippard concluded that F. London’s idea of a long-range ordering of the
superconducting phase is consistent with his observations. Coming from that
Pippard arrived at the concept of nonlocality Pippard53 , which paved the way
for the BCS theory.

Now, let us consider what one can say about $\lambda_{L}$ based solely on
thermodynamics.

By definitions of the condensation energy $E_{c}(T)$ and the thermodynamic
critical field $H_{c}(T)$, the difference between the total free energies of
the sample in the superconducting ($\widetilde{F}_{s}$) and the normal
($\widetilde{F}_{n}$) states in zero field is686868In zero field
$\widetilde{F}(T,\textbf{H}_{0})=\hat{F}(T,\textbf{H})=F(T,\textbf{B})$, where
the last two functions are the Gibbs and Helmholtz free energies, which are
appropriate for samples of the cylindrical and transverse geometries,
respectively. The total free energy $\widetilde{F}$ is suitable for all
geometries VK .

$\widetilde{F}_{n0}-\widetilde{F}_{s0}=\widetilde{F}_{n}-\widetilde{F}_{s0}=E_{c}(T)=V\frac{H_{c}^{2}(T)}{8\pi},$
(62)

where subscript $0$ designates zero field; for the normal state
$\widetilde{F}_{n0}=\widetilde{F}_{n}$ since the magnetic permeability in the
normal state is unity by definition.

Next, by definition of the total free energy VK , $\widetilde{F}_{s}$ of a
cylindrical sample in the field $\textbf{H}_{0}$ at constant $T$ is

$\widetilde{F}_{s}(T,H_{0})=\widetilde{F}_{s0}-\int_{0}^{H_{0}}\textbf{M}d\textbf{H}_{0}=\widetilde{F}_{n}-V\frac{H_{c}^{2}}{8\pi}+V\frac{H_{0}^{2}}{8\pi}.$
(63)

Here we also used the Meissner state definition according to which M of the
sample in the MS does not depend on temperature at constant field, as it takes
place in all other diamagnetics VK . As known (see, e.g., Fig. 3), this is
consistent with experiment.

Now, from Eq. (63) we get

$S_{s}\equiv-\left(\frac{\partial\widetilde{F}_{s}}{\partial
T}\right)_{H_{0}}=S_{n}+\frac{VH_{c}}{4\pi}\frac{dH_{c}}{dT}.$ (64)

Figure 15: A cross section of the field-induced currents in Cooper pairs
$J_{cp}$ inside a thin traverse (with respect to H) slice of a sample in the
MS as follows from the MW model. The sample is a film with thickness $d$, like
the film in Fig. 3. The slice contains currents induced in one Cooper pair of
each micro-whirl of the sample; $\mu_{cp}$ is the magnetic moment induced in
one pair. The whirls represent identical solenoids parallel to the field
$\textbf{H}(=\textbf{H}_{0}$ in this geometry); their rms radius
$r_{i}=2\lambda_{L}$. The graph schematically shows the distribution of the
induction $B$ along the $z$ axis perpendicular to the sample surface;
$\lambda(\neq\lambda_{L})$ is the width of the surface layer in which $B$
decays from $H_{0}$ to zero.

Since none of the terms on the right hand side depend on the field, the
entropy of the superconducting sample in the MS does not depend on the field
(like in all other diamagnetics). This implies that $n_{s}$ and therefore
$\lambda_{L}$ are independent of the field. The same conclusion can be drawn
immediately from Eq. (59) basing on the MS definition. Following to Pippard
Pippard50 , the field independence of $\lambda_{L}$ means that the long-range
order of the superconducting phase extends over the entire volume of the
sample in the MS. Needless to say that this is exactly what follows from the
MW model.

Thus, thermodynamics applied to the London theory requires the field
dependence of $\lambda_{L}$, the agreement of which with experiment is, at
least, questionable. On the other hand, in accordance with experiment, the
pure thermodynamic approach excludes such a dependence. However, it may look
like that the latter approach ignores the presence of the penetration layer
where $B\neq 0$, which is confirmed in countless experiments. Let us now turn
to the MW model.

Fig. 15 schematically shows a cross section of a thin ($\approx$ 6 fm thick)
transverse slice of a sample in the MS, as it follows from the MW model. The
slice contains the induced currents in one Cooper pair $J_{cp}$ of each micro-
whirl of the sample.

As one can see, contrarily to the London theory, there is no excluded volume.
The currents induced in a unit volume just next to the sample boundary (at $z$
near $0$ and $d$) make the same contribution to the sample magnetic moment as
the induced currents in the unit volume in any other part of the sample.
Therefore, the magnetic moment per unit volume $\chi$ is the same everywhere
inside the sample. Then, the sample magnetic moment is

$M=\chi HV=-\frac{1}{4\pi}HV=-\frac{1}{4\pi}\frac{H_{0}}{(1-\eta)}V,$

where $H$ is the field intensity, $\eta$ is the demagnetizing factor, which
equals zero for the samples of cylindrical geometry, and $V$ is the entire
sample volume as it was in the pure thermodynamic approach.

Hence, taking into account that $r_{i}$ is the root mean square value,
implying that the currents in Fig. 15 overlap, we arrive at the distribution
of the averaged microscopic field, i.e. $B$, as schematically shown in the
graph in Fig. 15. Thus, the MW model explains both the validity of the solely
thermodynamic approach as well as the reality of the penetration layer.

Finishing our discussion of H. London’s argument or the second dilemma of the
London theory, we see that its erroneousness, as in the case of the first
dilemma, stems from the assumption of the screening current, which, in
particular, leads to the appearance of an apparently non-physical excluded
volume in this theory.

## ACKNOWLEDGMENTS

It is my pleasant duty to thank Dr. A.-M. Valente-Feliciano, Dr. A. Suter, Dr.
T. Prokscha and Prof. C. Van Haesendonck for their crucial contributions into
the LE-$\mu$SR project temporarily interrupted by the pandemic. I am deeply
grateful to Prof. Chris Van Haesendonck for  his hospitality and many years of
support, and Prof. Joseph O. Indekeu for driving me into studies of
superconductivity. I am pleased to thank Prof. Tony Leggett for an encouraging
response on a short version of the manuscript. Special thanks to Prof.
Vladimir Z. Kresin for important recommendations after reading an original
version of the manuscript and for the invitation to write this review.

  1. (1) W. H. Keesom and J. N. van den Ende, KNAW Proceedings 35, 143 (1932); W. H. Keesom and J. A. Kok, KNAW Proceedings 35, 743 (1932); W. H. Keesom and J. A. Kok, Commun. Phys. Lab. Univ. Leiden, No 230c (1932).
  2. (2) H. Kamerlingh Onnes, KNAW Proceedings: 13, 1274 (1911); 14, 113 (1911); 14, 818 (1912); 16, 673 (1914).
  3. (3) D. van Delft, Freezing physics. Heike Kamerlingh Onnes and the quest for cold, (KNAW, Amsterdam, 2007).
  4. (4) H. Kamerlingh Onnes, Phys. Lab. Univ. Leiden, Suppl. No 50a (1924).
  5. (5) J. Mehra, The Solvay Conferences on Physics (D. Reidel Publishing Co, Dordrecht-Holland, 1975).
  6. (6) F. London, Superfluids v. I (N.Y., Wiley, 1950; N.Y. Dover, 1960).
  7. (7) D. Shoenberg, Superconductivity, 2nd. ed., (Cambridge, University Press, 1962).
  8. (8) V. Kozhevnikov, Thermodynamics of Magnetizing Materials and Superconductors (CRC Press, Boca Raton, 2019).
  9. (9) W. Meissner and R. Ochsenfeld, Naturwissenschaften 21, 787 (1933); for English translation see A. M. Forrest, Eur. J. Phys. 4, 117 (1983).
  10. (10) G. N. Rjabinin and L. W. Shubnikow, Nature 134, 286 (1934).
  11. (11) H. G. Smith and J. O. Wilhelm, Rev. Mod. Phys. 7, 237 (1935).
  12. (12) V. Kozhevnikov, A. Suter, T. Prokscha and C. Van Haesendonck, J. Supercond. Nov. Magnetism 33, 3361 (2020).
  13. (13) K. Mendelssohn and J. D. Babbitt, Nature 133, 459 (1934).
  14. (14) F. G. Tarr and J. O. Wilhelm, Can. J. Research 12, 265 (1935); Trans. Roy. Soc. Canada 28, 61 (1934).
  15. (15) K. Mendelssohn and J. D. Babbitt, Proc. R. Soc. A 151, 316 (1935).
  16. (16) J. C. Maxwell, A Treatise on Electricity and Magnetism, v.II, 2nd ed. (Oxford, Clarendon Press, 1881).
  17. (17) L. D. Landau, E.M. Lifshitz and L. P. Pitaevskii, Electrodynamics of Continuous Media, 2nd ed. (Elsevier, 1984).
  18. (18) L. D. Landau and E. M. Lifshitz, Quantum Mechanics, 3d ed, (Elsevier, Amsterdam, 2003).
  19. (19) C. J. Gorter and H. B. G. Casimir, Phys. Z. 35, 963 (1934).
  20. (20) R de Bruyn Ouboter, D. van Delft and P. H. Kes, in 100 Years of Superconductivity, Eds. H. Rogulla and P. H. Kes, p. 1 (CRC Press, Roca Raton, 2012).
  21. (21) J. A. Kok, Physica 1, 1103 (1934).
  22. (22) J. Bardeen, L. N. Cooper and Schrieffer, Phys. Rev. 108, 1175 (1957).
  23. (23) J. R. Schrieffer, Theory of Superconductivity, (Taylor & Francis Group, Boca Raton, 1999).
  24. (24) F. and H. London, Proc. Roy. Soc. A 149, 71 (1935).
  25. (25) V. L. Ginzburg and L. D. Landau, Zh.E.T.F. 20, 1064 (1950).
  26. (26) H. A. Lorentz, Commun. Phys. Lab. Univ. Leiden, Suppl. No 50b to Nos 157-168, 37 (1924).
  27. (27) E. H. Hall, Proc. Nat. Acad. Sci. Wash. 19, 619 (1933).
  28. (28) A. B. Pippard, Biogr. Mems Fell. R. Soc. 51, 379 (2005).
  29. (29) D. Shoenberg, Biogr. Mems Fell. R. Soc. 31, 326 (1985).
  30. (30) A. B. Pippard, Proc. Roy. Soc. London A 216, 547 (1953).
  31. (31) A. Suter, E. Morenzoni, N. Garifianov, R. Khasanov, E. Kirk, H. Luetkens, T. Prokscha, and M. Horisberger, Phys. Rev. B 72, 024506 (2005).
  32. (32) V. F. Kozhevnikov, C. V. Giuraniuc, M. J. Van Bael, K. Temst, C. Van Haesendonck, T. M. Mishonov, T. Charlton, R. M. Dalgliesh, Yu. N. Khaidukov, Yu. V. Nikitenko, V. L. Aksenov, V. N. Gladilin, V. M. Fomin, J. T. Devreese, and J. O. Indekeu, Phys. Rev. B 78, 012502 (2008).
  33. (33) V. Kozhevnikov, A. Suter, H. Fritzsche, V. Gladilin, A. Volodin, T. Moorkens, M. Trekels, J. Cuppens, B. M. Wojek, T. Prokscha, E. Morenzoni, G. J. Nieuwenhuys, M. J. Van Bael, K. Temst, C. Van Haesendonck, J. O. Indekeu, Phys. Rev. B 87, 104508 (2013).
  34. (34) R. Feynman, The Character of Physical Law (Cox and Wyman LTD, London, 1965).
  35. (35) B. S. Deaver, Jr., and W. M. Fairbank, Phys. Rev. Letters 7, 43 (1961).
  36. (36) L. Onsager, Phys. Rev. Letters 7, 50 (1961).
  37. (37) C. J. Gorter, Rev. Mod. Phys. 36, 3 (1964).
  38. (38) M. Desirant, D. Shoenberg, Proc. Roy. Soc. A 194, 63 (1948).
  39. (39) V. Kozhevnikov, A.-M. Valente-Feliciano, P. J. Curran, G. Richter, A. Volodin, A. Suter, S. J. Bending, C. Van Haesendonck, J. Supercond. Nov. Magnetism 31, 3433 (2018).
  40. (40) I. E. Tamm, Fundamentals of the Theory of Electricity (Mir, Moscow, 1979).
  41. (41) J. D. Jackson, Classical Electrodynamics, 3d ed. (John Wiley & Sons, Inc., Hoboken NJ, 1999).
  42. (42) I. Kikoin and I. Fakidov, Z. Physik 71, 393 (1931); reprinted in I. K. Kikoin - Physics and Fate, Ed. S. S. Yakimov, p. 44 (Nauka, Moscow, 2008).
  43. (43) J. H. Van Vleck, The Theory of Electric and Magnetic Susceptibilities (Clarendon Press, Oxford, 1932).
  44. (44) R. Feynman, R. Leighton, M. Sands, The Feynman Lectures on Physics, v. II (Basic Books, N.Y., 1964).
  45. (45) L. D. Landau and E. M. Lifshitz, The Classical Theory of Field, 4th Edition (Elseiver Ltd, Oxford, 1975).
  46. (46) L. N. Cooper, Phys. Rev. 104, 1189 (1956).
  47. (47) M. Tinkham, Introduction to Superconductivity (McGraw-Hill, 1996).
  48. (48) L. D. Landau, Z. Phys. 64, 629, 1930.
  49. (49) L. D. Landau and E. M. Lifshitz, Statistical Physics, part I, 3d ed., (Elsevier, Amsterdam, 1980).
  50. (50) C. Kittel, Introduction to Solid State Physics, 8th Ed. (John Wiley & Sons, Hoboken, 2005).
  51. (51) H. London, Proc. Roy. Soc. A 155, 102 (1936).
  52. (52) E. M. Purcell, Electricity and Magnetism.  Berkeley Physics Course-v.2 (McGraw-Hill, Boston, 1985).
  53. (53) D. J. Griffiths Introduction to Electrodynamics, 4th Ed. (Cambridge University Press, Cambridge 2017).
  54. (54) I.K. Kikoin and S. V. Goobar, C. R. Acad. Sci. USSR 19, 249 (1938); J. Phys. USSR 3, 333 (1940).
  55. (55) I. K. Kikoin, Zh. Tekh. Fiz. (Technical Physics) 166, 129 (1946); reprinted in I. K. Kikoin - Physics and Fate, Ed. S. S. Yakimov, p. 145 (Nauka, Moscow, 2008), in Russian.
  56. (56) G. S. Landsberg, Optics, 3d Ed. (Nauka, Moscow, 1976).
  57. (57) H. Kroemer, Quantum Mechanics (Prentice-Hall, Inc., New Jersey, 1994).
  58. (58) M. Peshkin, A. Tonomura, The Aharonov–Bohm effect (Springer-Verlag, N.Y., 1989).
  59. (59) I. S. Grigoriev, E. Z. Meilikhov, A. A. Radzig, Eds., Handbook of Physical Quantities (CRC Press, Boca Raton, 1997).
  60. (60) R. P. Huebener, Magnetic Flux Structures in Superconductors, 2nd Ed. (Springler-Verlag, N.Y., 2010).
  61. (61) V. Kozhevnikov, R. J. Wijngaarden, J. de Wit, and C. Van Haesendonck. PRB 89, 100503(R) (2014).
  62. (62) U. Essmann and H. Träuble, Phys. Letters 24A, 526 (1967).
  63. (63) R. P. Feynman, Progress in Low Temperature Physics, v. 1, p. 17 (1955); Rev. Mod. Phys. 29, 205 (1957).
  64. (64) G. M. Androes and W. D. Knight, Phys. Rev. Letters 2, 386 (1959).
  65. (65) F. Reif, Phys. Rev. 106, 208 (1957).
  66. (66) K. Ishida, M. Manago, K. Kinjo and Y. Maeno, J. Phys. Soc. Japan 89, 034712 (2020).
  67. (67) A.B. Pippard, Proc. Roy. Soc. London A 203, 210 (1950).
  68. (68) Report on experiment ID 20180450, Paul Scherrer Institut (2018); V. Kozhevnikov, A.-M. Valente-Feliciano, A. Suter, T. Prokscha and C. Van Haesendonck, to be published.
  69. (69) H. Kamerlingh Onnes and K. Hof, Commun. Phys. Lab. Univ. Leiden, no. 142b, (1914); Proc. Roy. Soc. Netherlands 17, 520 (1914).
  70. (70) L. V. Shubnikov, V. I. Khotkevich, Yu. D. Shepelev, Yu. N. Ryabinin, Zh.E.T.F. 7, 221 (1937).
  71. (71) E. Zeldov, in 100 Years of Superconductivity, Eds. H. Rogulla and P. H. Kes, p. 222 (CRC Press, Roca Raton, 2012).
  72. (72) A. B. Pippard, Math. Proc. Cambridge Phil. Soc. 47, 617 (1951).
  73. (73) H. London, Proc. Roy. Soc. London A 152, 650 (1935).
  74. (74) L. D. Landau, Zh.E.T.F. 7, 371 (1937).
  75. (75) E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics Part 2 ( Butterworth-Heinemann, 1980).
  76. (76) V. Kozhevnikov, A.-M. Valente-Feliciano, P. J. Curran, A. Suter, A. H. Liu, G. Richter, E. Morenzoni, S. J. Bending, and C. Van Haesendonck, Phys. Rev. B 95, 174509 (2017).
  77. (77) L. P. Gorkov, in 100 Years of Superconductivity, p. 72, Eds. H. Rogalla and P. H. Kes (CRC Press, Roca Raton, 2012).
  78. (78) A. B. Pippard, IEEE Trans. Magnetics 23, 371 (1987).
  79. (79) A. Shalnikov, Nature 142, 74 (1938).
  80. (80) L. B. Okun, Physics-Uspekhi 178, 653 (2008); arXiv:0809.2379 (2008).
  81. (81) S. Sridhar and J. E. Mercereau, Phys. Rev. B 34, 203 (1986).
  82. (82) P. G. De Gennes, Superconductivity of Metals and Alloys (Perseus Book Publishing, L.L.C., 1966).
  83. (83) A. Einstein and W.J. de Haas, KNAW Proceedings (Amsterdam) 18, 696 (1915).
  84. (84) H. London, Physical Society Report on a Low Temperature Conference, London, p. 51 (1947), cited from London50 .


# Robust parameter estimation for partially observed second-order diffusion
processes

Jan Albrecht and Sebastian Reich  Institut für Mathematik, Universität
Potsdam, Karl-Liebknecht-Str. 24/25, 14476 Potsdam

###### Abstract.

Estimating parameters of a diffusion process given continuous-time
observations of the process via maximum likelihood approaches or, online, via
stochastic gradient descent or Kalman filter formulations constitutes a well-
established research area. It has also been established previously that these
techniques are, in general, not robust to perturbations in the data in the
form of temporal correlations. While the subject is relatively well understood
and appropriate modifications have been suggested in the context of multi-
scale diffusion processes and their reduced model equations, we consider here
an alternative setting where a second-order diffusion process in positions and
velocities is only observed via its positions. In this note, we propose a
simple modification to standard stochastic gradient descent and Kalman filter
formulations, which eliminates the arising systematic estimation biases. The
modification can be extended to standard maximum likelihood approaches and
avoids computation of previously proposed correction terms.

##  1\. Introduction

We consider the problem of online learning parameters
$\theta\in\mathbb{R}^{d_{\theta}}$ of second-order stochastic differential
equations (SDEs)

(1a) $\displaystyle{\rm d}U_{t}$ $\displaystyle=f(X_{t},U_{t}){\rm
d}t+F(X_{t},U_{t})\theta{\rm d}t+\sqrt{\sigma}{\rm d}W_{t},$ (1b)
$\displaystyle{\rm d}X_{t}$ $\displaystyle=U_{t}{\rm d}t$

from incoming observations $X_{t}^{\dagger}\in\mathbb{R}^{d}$, $t\geq 0$, of
an underlying reference process. Here
$f:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ and
$F:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d_{\theta}}$ are
given functions, $\sigma>0$ is the diffusion constant, and $W_{t}$ denotes
$d$-dimensional Brownian motion.

The data $X_{t}^{\dagger}$ is being generated from a second-order SDE of the
same type, i.e.,

(2a) $\displaystyle{\rm d}U_{t}^{\dagger}$
$\displaystyle=g(X_{t}^{\dagger},U_{t}^{\dagger}){\rm d}t+\sqrt{\sigma}{\rm
d}W_{t}^{\dagger},$ (2b) $\displaystyle{\rm d}X_{t}^{\dagger}$
$\displaystyle=U_{t}^{\dagger}{\rm d}t$

with the drift $g(x,u)$ being unknown. We assume that there exists a unique
$\theta_{\ast}$ such that

(3) $g(x,u)=f(x,u)+F(x,u)\theta_{\ast}.$

While standard online estimation techniques, such as stochastic gradient
descent (SGD) [17] and Kalman filtering [11, 16], can be applied in case both
$(X_{t}^{\dagger},U_{t}^{\dagger})$ from (2) are being observed, the
estimation problem becomes more challenging in case only positions
$X_{t}^{\dagger}$ are being observed and velocities need to be estimated using
finite differencing [15, 5], i.e.,

(4)
$\tilde{U}_{n+1/2}:=\frac{X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger}}{\tau}$

with $\tau>0$ denoting the sampling interval, $t_{n}=\tau n$, for all $n\geq
0$.

As demonstrated previously [9], naive application of standard maximum
likelihood-based (MLE) estimators [10, 13] leads to a systematic bias. In
order to correct for the arising estimation bias, appropriate corrections
terms to the MLE estimator have been proposed in [9, 5, 3]. Related phenomena
have been investigated in the context of multi-scale diffusion [2, 12, 4],
which again can be corrected via an appropriately modified MLE estimator [8].

While these correction terms can be applied to online learning, we follow an
alternative approach in this note, which has been motivated by the related
work [1, 14]. Our approach can easily be implemented both for SGD [17] as well
as Kalman filter [11, 16] formulations.

As demonstrated in the next section by means of a simple example, a
discretised SGD of the form

(5a) $\displaystyle\theta_{n+1}$
$\displaystyle:=\theta_{n}+\alpha_{t}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\Delta I_{n}(\theta_{n}),$ (5b) $\displaystyle\Delta I_{n}(\theta)$
$\displaystyle:=(\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2})-\left(f(X_{t_{n}}^{\dagger},\tilde{U}_{n})+F(X_{t_{n}}^{\dagger},\tilde{U}_{n})\theta\right)\tau,$

with $U_{t_{n+1/2}}$ being replaced by its estimated $\tilde{U}_{n+1/2}$ etc.
and the symmetric approximation

(6)
$\tilde{U}_{n}:=\frac{1}{2}\left(\tilde{U}_{n+1/2}+\tilde{U}_{n-1/2}\right)=\frac{X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}}{2\tau}$

will fail even in the limit $\tau\to 0$. Here

(7) $\alpha_{t}=c_{1}/(c_{2}+t)$

denotes the learning rate for appropriate constants $c_{i}>0$, $i=1,2$.

The discussion from the following section will lead to an appropriate
modification of (5), which removes the undesirable bias in the SGD estimator
(5). See Section 3.

We extend the proposed methodology to Kalman filter based estimation in
Section 4. The unbiasedness of the corrected SGD approach is demonstrated for
a nonlinear second-order diffusion process in Section 5. The note closes with
comments on further directions for research.

##  2\. Motivating example

In order to demonstrate the difficulties the standard SGD implementation (5)
encounters, we consider velocity data $U_{t}^{\dagger}$ being generated by the
linear OU process

(8) ${\rm d}U_{t}^{\dagger}=-\theta_{\ast}U_{t}^{\dagger}{\rm
d}t+\sqrt{\sigma}{\rm d}W_{t}^{\dagger}$

for given diffusion constant $\sigma>0$ and unknown drift parameter
$\theta_{\ast}\geq 0$.

We will set $\theta_{\ast}=0$ for some of the analysis to be conducted later
in this section, which implies that solutions of (8) satisfy

(9) $U_{t}^{\dagger}=U_{s}^{\dagger}+\sqrt{\sigma}W_{t-s}^{\dagger}$

for any $t\geq s$ and $s\geq 0$ and, consequently,

(10)
$X_{t}^{\dagger}=X_{s}^{\dagger}+(t-s)U_{s}^{\dagger}+\sqrt{\sigma}\int_{s}^{t}W_{r-s}{\rm
d}r.$

Here we have introduced the shorthand
$W_{t-s}^{\dagger}=W_{t}^{\dagger}-W_{s}^{\dagger}$.

We first recall some of the properties of SGD when a solution of (8) has in
fact been fully observed. An application of SGD [17] to this problem yields

(11a) $\displaystyle{\rm d}\theta_{t}$
$\displaystyle=-\alpha_{t}U_{t}^{\dagger}{\rm
d}I_{t}=-\alpha_{t}U_{t}^{\dagger}\left({\rm
d}U_{t}^{\dagger}+\theta_{t}U_{t}^{\dagger}{\rm d}t\right)$ (11b)
$\displaystyle=-\alpha_{t}U_{t}^{\dagger}\left((\theta_{t}-\theta_{\ast})U_{t}^{\dagger}{\rm
d}t+\sqrt{\sigma}{\rm d}W_{t}^{\dagger}\right).$

We call

(12) $I_{t}=\int_{0}^{t}\left({\rm
d}U_{s}^{\dagger}+\theta_{s}U_{s}^{\dagger}{\rm
d}s\right)=\int_{0}^{t}\left((\theta_{s}-\theta_{\ast})U_{s}^{\dagger}{\rm
d}t+\sqrt{\sigma}{\rm d}W_{s}^{\dagger}\right)$

the innovation and $\alpha_{t}U_{t}^{\dagger}$ the gain. It is obvious that
$I_{t}$ becomes a martingale at the optimal parameter choice
$\theta_{s}=\theta_{\ast}$, $s\in[0,t]$.

It has been shown that $\theta_{t}\to\theta_{\ast}$ as $t\to\infty$ [17]. In
our context, a key observation is that the integral

(13) $\int_{0}^{t}\alpha_{s}U_{s}^{\dagger}\left({\rm
d}U_{s}^{\dagger}+\theta_{t}U_{s}^{\dagger}{\rm d}s\right)$

also becomes a martingale for any $t>0$ if and only if
$\theta_{s}=\theta_{\ast}$ for all $s\in[0,T]$.

In terms of numerical implementations of the SGD formulation 11), we consider
the discrete time SGD method

(14)
$\theta_{n+1}:=\theta_{n}-\alpha_{t_{n}}U_{t_{n}}^{\dagger}\left(U_{t_{n+1}}^{\dagger}-U_{t_{n}}^{\dagger}+\frac{\theta_{n}}{2}\left(U_{t_{n+1}}^{\dagger}+U_{t_{n}}^{\dagger}\right)\tau\right).$

This discrete-time SGD method still implies that $\theta_{n}$ converges to
$\theta_{\ast}$ as $n\to\infty$ and $\tau\to 0$. More precisely, any finite
sum

(15)
$\sum_{n=0}^{N-1}U_{t_{n}}^{\dagger}\left(U_{t_{n+1}}^{\dagger}-U_{t_{n}}^{\dagger}+\frac{\theta_{n}}{2}\left(U_{t_{n+1}}^{\dagger}+U_{t_{n}}^{\dagger}\right)\tau\right)$

becomes a martingale as $\tau\to 0$ for $T=N\tau$ fixed if and only if
$\theta_{n}=\theta_{\ast}$.

However, the focus of this note is on the case when only positions
$X_{t}^{\dagger}$, $t\geq 0$, can be observed and we recall the velocity
approximations (4), where the sampling interval $\tau>0$ can be made
arbitrarily small. We note that (4) can also be viewed as a local time
average, i.e.,

(16)
$\tilde{U}_{n+1/2}=\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}U_{t}^{\dagger}{\rm
d}t,$

which provides a link to the related work [9, 7].

Using (6) instead of the true $U_{t_{n}}$ etc. lead to the previously stated
(5), which reduces here to

(17a) $\displaystyle\theta_{n+1}$
$\displaystyle:=\theta_{n}-\alpha_{t_{n}}\tilde{U}_{n}\left(\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2}+\theta_{n}\tilde{U}_{n}\tau\right)$
(17b)
$\displaystyle=\theta_{n}-\frac{\alpha_{t_{n}}}{2}\left(X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}\right)\left(\frac{X_{t_{n+1}}^{\dagger}-2X_{t_{n}}^{\dagger}+X_{t_{n-1}}^{\dagger}}{\tau^{2}}+\frac{\theta_{n}}{2}\frac{X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}}{\tau}\right).$

We note that, while the discrete innovations

(18)
$I_{N}=\sum_{n=1}^{N-1}\left(\frac{X_{t_{n+1}}^{\dagger}-2X_{t_{n}}^{\dagger}+X_{t_{n-1}}^{\dagger}}{\tau^{2}}+\frac{\theta_{n}}{2}\frac{X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}}{\tau}\right)\tau$

still converge to a martingale as $\tau\to 0$ for $T=N\tau$ fixed provided
$\theta_{n}=\theta_{\ast}$ for all $n\geq 0$, this is no longer true for

(19)
$\frac{1}{2}\sum_{n=1}^{N-1}\left(X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}\right)\left(\frac{X_{t_{n+1}}^{\dagger}-2X_{t_{n}}^{\dagger}+X_{t_{n-1}}^{\dagger}}{\tau^{2}}+\frac{\theta_{n}}{2}\frac{X_{t_{n+1}}^{\dagger}-X_{t_{n-1}}^{\dagger}}{\tau}\right).$

This violation of the martingale property is perhaps not surprising since (19)
resembles a discretised Stratonovitch integral [13]. We note that such an
interpretation as a discretised stochastic integral concerns only the second-
order differencing term in (19) and we may therefore set
$\theta=\theta_{\ast}=0$ in the subsequent derivations. Recall that

(20)
$\frac{X_{t_{n+1}}-X_{t_{n}}}{\tau}=\tilde{U}_{n+1/2}=U_{t_{n}}^{\dagger}+\frac{\sqrt{\sigma}}{\tau}\int_{t_{n}}^{t_{n+1}}W_{t-t_{n}}^{\dagger}{\rm
d}t$

and we therefore obtain

(21a)
$\displaystyle\frac{1}{2}\mathbb{E}\left[\left(\tilde{U}_{n+1/2}+\tilde{U}_{n-1/2}\right)(\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2})\right]$
$\displaystyle=\frac{1}{2}\mathbb{E}\left[\tilde{U}_{n+1/2}^{2}-\tilde{U}_{n-1/2}^{2}\right]$
(21b)
$\displaystyle=\frac{1}{2}\mathbb{E}\left[(U_{t_{n}}^{\dagger})^{2}-(U_{t_{n-1}}^{\dagger})^{2}\right]=\frac{\tau\sigma}{2},$

which implies an estimation bias induced by the discrete sum (19) even in the
limit $\tau\to 0$.

In order to remove the resulting bias in the SGD formulation (17), we propose
to simply shift the innovation term forward in time such that it does no
longer overlap with the gain; thus regaining an Itô-type stochastic integral
approximation. A related idea has been proposed and investigated in [14] in
the context of multi-scale diffusion. This modification leads to the following
implementation:

(22)
$\theta_{n+1}:=\theta_{n}+\alpha_{t_{n}}\tilde{U}_{n}\left(\tilde{U}_{n+5/2}-\tilde{U}_{n+3/2}-\theta_{n}\tilde{U}_{n+2}\tau\right).$

Let us introduce the random variables

(23)
$\tilde{\Xi}_{n}:=\frac{1}{\sqrt{\sigma\tau}}\left(\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2}\right)$

for $n\geq 1$, which is a rescaling of the innovation $\Delta
I_{n}(\theta_{\ast})$, as defined in (5b), with $\theta_{\ast}$ again set to
$\theta_{\ast}=0$. It is crucial to observe that $\tilde{\Xi}_{n+2}$ is
correlated with both $\tilde{U}_{n+5/2}$ and $\tilde{U}_{n+3/2}$, while being
independent of $\tilde{U}_{n}$. Hence, provided $\tilde{\Xi}_{n}$ is a
martingale for all $n\geq 1$, the SGD scheme (22) will become unbiased in the
limit $\tau\to 0$.

Hence let us analyse the statistical properties of $\tilde{\Xi}_{n}$. We find
that

(24a) $\displaystyle\tilde{\Xi}_{n}$
$\displaystyle=\frac{1}{\sqrt{\sigma\tau}}\left(\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2}\right)$
(24b)
$\displaystyle=\frac{1}{\sqrt{\sigma\tau}}\left(U_{t_{n}}^{\dagger}-U_{t_{n-1}}^{\dagger}\right)+\frac{1}{\sqrt{\tau}}\left(\frac{\int_{t_{n}}^{t_{n+1}}W_{t-t_{n}}^{\dagger}{\rm
d}t-\int_{t_{n-1}}^{t_{n}}W_{t-t_{n-1}}^{\dagger}{\rm d}t}{\tau}\right)$ (24c)
$\displaystyle=\frac{W_{t_{n}-t_{n-1}}^{\dagger}}{\sqrt{\tau}}+\frac{1}{\sqrt{\tau}}\left(\frac{\int_{t_{n}}^{t_{n+1}}W_{t-t_{n}}^{\dagger}{\rm
d}t-\int_{t_{n-1}}^{t_{n}}W_{t-t_{n-1}}^{\dagger}{\rm d}t}{\tau}\right).$

It follows that $\tilde{\Xi}_{n}$ is Gaussian with mean zero and that
$\tilde{\Xi}_{n+2}$ is indeed independent of $\tilde{U}_{n}$.

We recall the following two properties of integrated Brownian motion:

(25) $\mbox{var}\,\left(\int_{t_{n}}^{t_{n+1}}W_{t-t_{n}}^{\dagger}{\rm
d}t\right)=\frac{1}{3}\tau^{3}.$

and

(26)
$\mathbb{E}\left[W_{t_{n}-t_{n-1}}^{\dagger}\left(\int_{t_{n-1}}^{t_{n}}W_{t-t_{n-1}}^{\dagger}{\rm
d}t\right)\right]=\int_{t_{n-1}}^{t_{n}}\mathbb{E}\left[(W_{t-t_{n-1}}^{\dagger})^{2}\right]{\rm
d}t=\frac{1}{2}\tau^{2}.$

Upon applying these results to (24c) when computing the variance of
$\tilde{\Xi}_{n}$, it is found that

(27) $\mbox{var}\,(\tilde{\Xi}_{n})=\frac{2}{3}.$

In other words, we obtain that $\tilde{\Xi}_{n}\sim{\rm N}(0,2/3)$ in
agreement with results from [9, 6]. This result carries over to
$\theta_{\ast}\not=0$ since the arising additional terms vanish as $\tau\to
0$.

In line with the previous result [9], our calculations imply the following
estimator for the diffusion constant $\sigma$:

(28) $\sigma_{N}:=\frac{3}{2\tau
N}\sum_{n=1}^{N}\|\tilde{U}_{n+1/2}-\tilde{U}_{n-1/2}\|^{2}.$

We finally formally investigate the time-continuous limit of the modified
discrete-time SGD formulation (22). We first note that

(29)
$\mathbb{E}[\tilde{\Xi}_{n}\tilde{\Xi}_{n+k}]=\left\\{\begin{array}[]{cc}\frac{1}{6}&|k|=1\\\
0&|k|>1\end{array}\right.,$

which implies that $\sqrt{\tau}\tilde{\Xi}_{n}$ converges to standard Brownian
motion in the limit $\tau\to 0$. When using this result in the SGD formulation
(22) for general $\theta_{\ast}$, we obtain

(30a) $\displaystyle\theta_{n+1}$
$\displaystyle\approx\theta_{n}-\alpha_{t_{n}}\tilde{U}_{n}\left((\theta_{n}-\theta_{\ast})\tilde{U}_{n+2}\tau+\sqrt{\sigma\tau}\tilde{\Xi}_{n+2}\right)$
(30b)
$\displaystyle\approx\theta_{n}-\alpha_{t_{n}}\tilde{U}_{n}\left((\theta_{n}-\theta_{\ast})\tilde{U}_{n}\tau+\sqrt{\sigma\tau}\tilde{\Xi}_{n+2}\right)$

and the limit $\tau\to 0$ leads formally to the time-continuous formulation

(31)
$\dot{\theta}_{t}=-\alpha_{t}U_{t}^{\dagger}\left((\theta-\theta_{\ast})U_{t}^{\dagger}+\sqrt{\sigma}\dot{W}_{t}^{\dagger}\right),$

which coincides with (11b).

##  3\. Bias-corrected SGD

The analysis from Section 2 leads us to the following modification of the
discrete-time SGD implementation (5):

(32)
$\theta_{n+1}=\theta_{n}+\alpha_{t_{n}}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\Delta I_{n+2}(\theta_{n}),$

for $n\geq 1$ with the innovations $\Delta I_{n}(\theta)$ defined by (5b) and
the learning rate $\alpha_{t_{n}}$ by (7) with $t=t_{n}=n\tau$. We recall that
the key point of shifting the innovation forward in time is to retain an Itô-
type approximation in the gain times innovation stochastic integral
approximation. This is the property that maintains the desired unbiasedness of
the resulting sequential learning scheme in the limit $\tau\to 0$.

We note that unbiased formulations are not unique. For example, (32) could be
replaced by

(33)
$\theta_{n+1}=\theta_{n}+\alpha_{t_{n}}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\left(\tilde{U}_{n+5/2}-\tilde{U}_{n+3/2}-\left(f(X_{t_{n}}^{\dagger},\tilde{U}_{n})+F(X_{t_{n}}^{\dagger},\tilde{U}_{n})\theta_{n}\right)\tau\right).$

This freedom is also utilised when formulating the modified MLE estimator (37)
further below in Section 6.

##  4\. Bias-corrected Kalman filter

The SGD method (32) provides an online point estimator for the parameter
$\theta$. However, it is also possible to implement an online Bayesian
estimator using the Kalman filter [11, 16]. The key idea is to treat the
unknown parameter as a Gaussian random variable $\Theta_{n}$ with mean $m_{n}$
and covariance $\Sigma_{n}$, $n\geq 0$. The evolution equation for the mean
and the covariance follow the standard Kalman update formulas:

(34a) $\displaystyle m_{n+1}$ $\displaystyle:=m_{n}+K_{n}\frac{\Delta
I_{n+2}(m_{n})}{\tau},$ (34b) $\displaystyle\Sigma_{n+1}$
$\displaystyle:=\Sigma_{n}-K_{n}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})\Sigma_{n}$

with Kalman gain

(35) $K_{n}:=\Sigma_{n}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\left(F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\Sigma_{n}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})+\frac{\sigma}{\tau}I\right)^{-1}$

for $n\geq 1$. The iteration is initialised from the Gaussian prior
$\Theta_{1}\sim{\rm N}(m_{\rm prior},\Sigma_{\rm prior})$.

In contrast to the SGD formulation (32), the learning rate $\alpha_{t_{n}}$ no
longer appears and has instead been submerged implicitly into the Kalman gain
$K_{n}$. Indeed, the singular values of the covariance matrix $\Sigma_{n}$
decay with rate $1/t_{n}$ under suitable ergodicity assumptions on (2). If the
required diffusion constant $\sigma$ is unknown, it can either be estimated
via (28) or, alternatively, be treated as part of the prior covariance matrix
$\Sigma_{\rm prior}$ and then set to $\sigma=1$ in the Kalman gain (35).

##  5\. Numerical example

We consider the nonlinear second-order SDEs

(36a) $\displaystyle{\rm d}U_{t}$ $\displaystyle=-X_{t}{\rm d}t+U_{t}{\rm
d}t-\theta U_{t}^{3}{\rm d}t+\sqrt{2}{\rm d}W_{t},$ (36b) $\displaystyle{\rm
d}X_{t}$ $\displaystyle=U_{t}{\rm d}t$

with unknown parameter $\theta\in\mathbb{R}$. The data $X_{t}^{\dagger}$ is
generated with $\theta_{\ast}=1$, while the SGD iteration starts from
$\theta_{0}=2$. The sample interval is set to $\tau=0.025$. The data is
generated with a much smaller step-size of $h=\tau/100$. A total of $N=10^{5}$
positions $X_{t_{n}}$, $n=1,\ldots,N$ have been generated. The learning rate
is given by $\alpha_{t_{n}}=0.2/n$. See Figure 1 for the time evolved
estimates $\theta_{n}$, which clearly converge to $\theta_{\ast}$. We also
display results from the standard SGD in Figure 1. While the standard SGD
converges slightly faster, it converges to a wrong parameter value.

Figure 1. Left panel: Estimated parameter value from standard SGD as a
function of the iteration index. The iteration starts at $\theta_{0}=2$ and
the reference value is $\theta_{\ast}=1$. A systematic bias can clearly be
seen. Right panel: Same data is used in the modified SGD. The modified SGD
algorithm converges to the correct reference value.

##  6\. Conclusions

In this note, we have proposed a simple-to-implement modification of standard
SGD and Kalman filter formulations for online parameter estimation from
position data only. The essential idea is to establish an Itô-like
modification that re-establishes the martingale property in the presence of
temporal correlations. The approach has been inspired by related work [1, 14]
on multi-scale diffusion and parameter estimation for the arising reduced
equations. We finally note that a corrected MLE estimator can be based on the
negative log-likelihood function

(37a) $\displaystyle\mathcal{L}(\theta,\sigma)$
$\displaystyle=\frac{1}{2\sigma}\sum_{n=1}^{N-3}\theta^{\rm
T}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})\theta\tau$ (37b)
$\displaystyle\qquad-\,\frac{1}{\sigma}\sum_{n=1}^{N-3}\theta^{\rm
T}F(X_{t_{n}}^{\dagger},\tilde{U}_{n})^{\rm
T}\left(\tilde{U}_{n+5/2}-\tilde{U}_{n+3/2}-f(X_{t_{n}}^{\dagger},\tilde{U}_{n})\tau\right)$
(37c)
$\displaystyle\qquad+\,\frac{3}{4\tau\sigma}\sum_{n=1}^{N-3}\|\tilde{U}_{n+5/2}-\tilde{U}_{n+3/2}\|^{2}+\frac{N-3}{2}\log\sigma.$

Future directions include the investigation of the proposed estimators in the
context multi-scale diffusion when (1) arises from homogenisation of an
underlying multi-scale diffusion process and the inclusion of observation
errors in the positions $X_{t}$. In the former case, a combination of the
approaches taken here and in [14] emerges as a feasible starting point.

### Acknowledgements.

This work has been partially funded by Deutsche Forschungsgemeinschaft (DFG) -
Project-ID 318763901 - SFB1294.

## References

  * [1] Assyr Abdulle, Giacomo Garegnani, Grigorios Pavliotis, Andrew Stuart, and Andrea Zanoni.  Drift estimation of multiscale diffusions based on filtered data.  Foundations of Computational Mathematics, 23:33–84, 2023.  doi:10.1007/s10208-021-09541-9.
  * [2] Yacine Ait-Sahalia, Per A Mykland, and Lan Zhang.  How often to sample a continuous-time process in the presence of market microstructure noise.  The Review of Financial Studies, 18:351–416, 2005.
  * [3] Jan Albrecht, Carsten Beta, Robert Großmann, and Manfred Opper.  Heterogeneity inference for second order Langevin models, 2024.  in preparation.
  * [4] S. Bo and A. Celani.  White-noise limit of nonwhite nonequilibrium processes.  Physical Review E, 88:062150, 2013.  doi:10.1103/PhysRevE.88.062150.
  * [5] D.B. Brückner, P. Ronceray, and C.P. Broedersz.  Inferring the dynamics of underdamped stochastic systems.  Phy. Rev. Lett., 125:28005, 2010.
  * [6] I. Cialenco, H.-J. Kim, and G. Pasemann.  Statistical analysis of discretely sampled semilinear SPDEs: A power variation approach.  Stochastics and Partial Differential Equations: Analysis and Computations, 12:326–351, 2024.
  * [7] Igor Cialenco and Y. Huang.  A note on parameter estimation for discretely sampled SPDEs.  Stoch. Dyn., 20:2050016, 2020.
  * [8] J. Diehl, P. Friz, and H. Mai.  Pathwise stability of likelihood estimators for diffusion via rough paths.  The Annals of Applied Probability, 26:2169–2192, 2016.  doi:10.1214/15-AAP1143.
  * [9] Arnaud Gloter.  Parameter estimation for a discretely observed integrated diffusion process.  Scandinavian Journal of Statistics, 33:83–104, 2006.
  * [10] Yury A Kutoyants.  Statistical inference for ergodic diffusion processes.  Springer Science & Business Media, 2013.
  * [11] Nikolas Nüsken, Sebastian Reich, and Paul J Rozdeba.  State and parameter estimation from observed signal increments.  Entropy, 21(5):505, 2019.  doi:10.3390/e21050505.
  * [12] A. Papavasiliou, G. Pavliotis, and A.M. Stuart.  Maximum likelihood estimation for multiscale diffusions.  Stochastic Processes and their Applications, 19:3173–3210, 2009\.
  * [13] G.A. Pavliotis.  Stochastic Processes and Applications.  Springer Verlag, New York, 2016.
  * [14] Grigorios A Pavliotis, Sebastian Reich, and Andrea Zanoni.  Filtered data based estimators for stochastic processes driven by colored noise, 2023.  arXiv:2312.15975.
  * [15] Jonas N. Pedersen, Liang Li, Cristian Gradinaru, Robert H. Austin, Edward C. Cox, , and Henrik Flyvbjerg.  How to connect time-lapse recorded trajectories of motile microorganisms with dynamical models in continuous time.  Physical Review E, 94:062401, 2016.
  * [16] Sebastian Reich.  Frequentist perspective on robust parameter estimation using the ensemble Kalman filter.  In B. Chapron, D. Crisan, D. Holm, E. Mémin, and A. Radomska, editors, Stochastic Transport in Upper Ocean Dynamics. STUOD 2021. Mathematics of Planet Earth, volume 10, pages 237–258, Cham., 2023. Springer.  doi:10.1007/978-3-031-18988-3_15.
  * [17] J. Sirignano and K. Spiliopoulos.  Stochastic gradient descent in continuous time.  SIAM J. Financial Math., 8:933–961, 2017.  doi:10.1137/17M1126825.