# Maximum Entropy Principle underlying the dynamics of automobile sales

D. Villuendas  Data Science Management and Market Intelligence, SEAT, S.A.,
Martorell, Spain  A. Hernando<EMAIL_ADDRESS>Social Thermodynamics
Applied Research (SThAR SA), EPFL Innovation Park Bât. C, Lausanne,
Switzerland  Departament FQA, Facultat de Física, Universitat de Barcelona,
Barcelona, Spain  M. Sulc  Social Thermodynamics Applied Research (SThAR SA),
EPFL Innovation Park Bât. C, Lausanne, Switzerland  R. Hernando  Social
Thermodynamics Applied Research (SThAR SA), EPFL Innovation Park Bât. C,
Lausanne, Switzerland  R. Seoane  Centres Científics i Tecnològics (CCiT),
Universitat de Barcelona, Barcelona, Spain  A. Plastino  Instituto de Física
La Plata-CCT-CONICET, Universidad Nacional de La Plata, La Plata, Argentina
Physics Department and IFISC-CSIC, University of the Balearic Islands, Palma
de Mallorca, Spain

###### Abstract

We analyze an exhaustive data-set of new-cars monthly sales. The set refers to
10 years of Spanish sales of more than 6500 different car model configurations
and a total of 10M sold cars, from January 2007 to January 2017. We find that
for those model configurations with a monthly market-share higher than 0.1%
the sales become scalable obeying Gibrat’s law of proportional growth under
logistic dynamics. Remarkably, the distribution of total sales follows the
predictions of the Maximum Entropy Principle for systems subject to
proportional growth in dynamical equilibrium. We also encounter that the
associated dynamics are non-Markovian, i.e., the system has a decaying memory
or inertia of about 5 years. Thus, car sales are predictable within a certain
time-period. We show that the main characteristics of the dynamics can be
described via a construct based upon the Langevin equation. This construct
encompasses the fundamental principles that any predictive model on car sales
should obey.

The automobile industry is experiencing a deep economic and technological
change. It is moving from (i) a fossil-fueled, private-ownership based,
manually-driven to (ii) an electric, sharing based, driver-less onepwc ;
mckinsey . Understanding the supply-and-demand dynamics of automobile sales
could help to achieve a smooth economic transition between (i) and (ii).
Indeed, 88.1 million of both cars and light commercial vehicles were sold
worldwide in 2016macquarieresearch , involving billions of dollars per year.
Unplanned circumstances or forecast failures could generate immense loses in
this market and its ties, as previously seen in past crisis and economic
transitions as that after 2007wiki .

Thanks to the rise of both Big Data and digital tools we have today access to
exhaustive databases that allow us to analyze socioeconomic systems using
ideas of statistical physics sociophysics1 ; sociophysics2 ; sociophysics3 .
Examples run from allometric laws on city populationscities1 ; cities2 ;
cities3 , evolution of firm sizesfirms , transportation networkstransport1 ;
transport3  and human mobilitymobility , to popularity of digital
productslogistic , the structure of the Internetnet  or even diffusion of
memesmemes . In such an environment we analyze here the sales of new
commercial vehicles in Spain from January 2007 to January 2017. The
corresponding exhaustive database is published by the Spanish Directorate-
General of Traffic (DGT) and contains registrations of more than 6 500
different configurations (model + body shape + engine) and more than 10
million sold carsdata . To analyze the evolution of car sales, we apply, in
particular, the procedure used in Refs. firms, ; entropy1, ; entropy2, ;
entropy3,  on the aggregated data of monthly registrations (or sales) per each
car configuration (model) as provided by Ref. data2 . The procedure is
summarized as follows:

  1. 1.

We first identify the main dynamical variables. In our case, the variable is
the total number of cars sold for the $i$-th automobile model at time $t$,
$x_{i}(t)$.

  2. 2.

We compare $x_{i}(t)$ versus its time derivative $\dot{x}_{i}(t)$ so as to
find indications of any possible scaling rule in the underlying equation of
motion, of the form

$\dot{x}_{i}(t)=v_{i}(t)[x_{i}(t)]^{q},$ (1)

where $v_{i}(t)$ is the growth rate at time $t$ and $q$ is the exponent that
parameterizes the dynamics. Due to the stochastic nature of the growth rates,
it is convenient to compare $x_{i}$ with the variance of $\dot{x}_{i}$,
obtained from all those $i$’s endowed with a similar value of $x$. Thus, we
fit the variance to an expression of the form
$\text{Var}[\dot{x}]=T_{q}x^{2q}$, where $T_{q}$ is used to measure the size
of the fluctuations, that might be called a “temperature”thermo2 ; thermo3 .
If one finds contributions of several independent components with the form of
Eq. (1) with different exponents —as
$\dot{x}_{i}(t)=\sum_{q}v_{q,i}(t)[x_{i}(t)]^{q}$— we will define a
temperature associated to each term from the fit
$\text{Var}[\dot{x}]=\sum_{q}T_{q}x^{2q}$.

  3. 3.

Next, we independently analyze the density distribution of total car sales per
model $p(x)dx$. The principle of Maximum Entropy (MaxEnt) with dynamical
information states that the natural variable for measuring the entropy is the
one that linearizes the equation of motion, transforming the underlying
symmetry into a translational oneMaxEntDyn . For an equation of motion as Eq.
(1), this is achieved via the Tsallis logarithm and its inverse function the
Tsallis exponentialtsallis  defining the new variable $u=\log_{q}(x/x_{0})$
—or $x=x_{0}\exp_{q}(u)$— where $x_{0}$ is any value of reference (needed to
keep the arguments dimensionless). After this transformation, the equation of
motion becomes $\dot{u_{i}}(t)=v_{i}(t)$, which is indeed translationally
invariant. If the system is in equilibrium, the most probable distribution is
the one that maximizes the entropy under the system’s observable constraints.
If the form of the empirical distribution fits the form predicted by MaxEnt,
we can perform an independent quantitative measure of $q$ directly from the
equilibrium distribution.

  4. 4.

If the value of $q$ one gets from the equation of motion, and the value of $q$
measured from the distribution-fit agree, we assume that our procedure is
correct, proving that the system is in dynamical equilibrium and obeys the
Maximum Entropy Principle.

  5. 5.

As an additional test for our procedure, we also perform numerical experiments
by defining microscopic equations of motion equivalent to those observed in
the empirical system. We then analyze the evolution of the simulated system
following the same steps as we traversed for the empirical one. If we obtain
equivalent results, we would reconfirm that our procedure is indeed correct.

Figure 1:  Top: Evolution of car sales for a sample of car-models of the
database. Saturation (as that described by the logistic equation) is observed.
Middle: Evolution of the number of different models available per month during
the ten years recorded in the database. Bottom: Monthly new car sales in
Spain, with a total of $10\,062\,193$ sales during the last ten years.

As stated in step 1, our relevant variable is the total of car sales $x_{i}$
for the $i$th automobile model. We display in Fig. 1 a sample of some
trajectories for the ten years covered by our dataset, where a logistic
profile is clearly visible (i.e., growth with saturation). Saturation occurs
when the popularity of a model drops because of the availability of a newer
version, or if a competitor advances on the same market niche. We also show in
Fig. 1 the evolution of the number of different models sold per month during
the ten years included in the dataset.

Figure 2:  Top: Empirical variance of the sales’ growth versus total car sales
(gray points). We compare the moving average (red line) with the analytical
fit to Eq. (1) (black line). We find three regimes with different exponents.
Bottom: same as top, but with growth corrected by the logistic equation (see
text), that fits to the form of Eq. (3).

When we compare $x$ vs. $\text{Var}[\dot{x}]$ so as to fit Eq. (1), we find
what looks like a three-component equation of motion, with an exponent i) for
low sales of $q=0.75\pm 0.01$, ii) for mid-sales of $q=0.63\pm 0.01$, and iii)
for high sales of $q=0.33\pm 0.01$. R2 adopts a value of $0.9997$ (see Fig.
2). A naive interpretation might be that Eq. (1) is indeed the underlying
equation of motion. The problem would be to explain the origin of the values
for three exponents mentioned above. Even worse, an independent assessment of
$q$ derived from the equilibrium distribution will not match these 3 values,
as we will show later. However, paying attention to the high sales in Fig. 2
(arrow), we appreciate the existence of points with low variance. Indeed, here
growth rates drop to such an extent that one approaches saturation and this
would generate an underestimation of the exponents. To get deeper insights
into the actual underlying dynamics, we need to correct for each car model the
effects of saturation. To this aim, we need to reconsider the form of the
logistic equationlogistic

$\dot{x_{i}}(t)=v_{i}(t)x_{i}(t)\left[1-x_{i}(t)/X_{i}\right],$ (2)

where $X_{i}$ is the final number of total car sales. If we reconsider the
trajectories by first defining
$\dot{x_{i}}^{\prime}(t)=\dot{x_{i}}(t)/[1-x_{i}(t)/X_{i}]$, we should recover
an expression of the form of Eq. (1) as $\dot{x_{i}}^{\prime}=v_{i}x_{i}$ that
corrects growth rates for the effects of saturation. Comparing in this way $x$
vs. $\text{Var}[\dot{x}^{\prime}]$, we obtain the curve displayed in Fig. 2,
which nicely describes a three-terms function: one has
$\text{Var}[\dot{x}^{\prime}]=T_{1}x^{2q}+T_{1/2}x+T_{0}$, with $q=1.016\pm
0.002$, and an R2 coefficient equal to $0.9998$ [$\log(T_{1})=-3.41\pm 0.05$,
$\log(T_{1/2})=1.60\pm 0.06$, and $\log(T_{0})=5.73\pm 0.06$]. This new result
vouches for an equation of motion of the form

$\dot{x}_{i}^{\prime}(t)=v_{1,i}(t)x_{i}(t)+v_{1/2,i}(t)\sqrt{x_{i}(t)}+v_{0,i}(t).$
(3)

These three components have been already observed for firm sizesfirms , city
populationsentropy2  and are well-understood. The first term with $q=1$ is
Gibrat’s law of proportional growthgibrat , which is expected for
multiplicative processes. Such behavior indicates that the popularity of an
automobile model grows as more cars are sold, following a rich-get-richer
mechanism. The second term, with exponent $q=1/2$, emerges from the
proportional growth as a side effect of the central limit theorem (or finite-
size effects), as shown in Ref. entropy3, . It has been shown that this
component is characterized by uncorrelated noise and no extra physics is
expected. Finally, the last term in Eq. (3) correspond to linear forces
($q=0$), independent of the value of $x$. This last term becomes relevant only
for small number of sales.

Figure 3:  Scalability (defined as the contribution of the proportional term
to the total growth in Eq. (3)) as a function of the monthly market share:
with 0.005% of the market share, the proportional term starts to be dominant,
reaching around 90% of the total contribution with 0.03% of monthly sales and
a 95% with 0.1% of the market (dashed blue lines and shaded areas). The
average maximum market share for a single model is around 2% (dashed gray
line).

In our procedure we will focus only on the proportional regime, where the most
interesting physics takes place. However, one may ask: which is the critical
number of sold cars so as to consider that the pertinent trajectory would be
accommodated by the proportional regime? A first estimation can be obtained
via the overlap between the first and second terms of Eq. (3), obtaining
$x_{0}=T_{1/2}/T_{1}=129$ cars. At this critical value we have equal
contributions from each of the two terms. It is then natural to ask when we
can consider that the dynamics is fully governed by proportional or scalable
growth. We can combine our two previous fits to compare the market share or
relative number of sold cars per month $\chi=\dot{x}/\sum_{j}{\dot{x_{j}}}$
with the relative contribution of the proportional term to the total growth,
that we call scalability:
$\kappa(x)=\sqrt{T_{1}x^{2}/\text{Var}[\dot{x_{i}}^{\prime}]}$. As seen in
Fig. 3, an automobile model is $\kappa=95\%$ into the scalable growth regime
when it reaches a $\chi=0.1\%$ of monthly relative sales or market-share.

Figure 4:  Left: Time correlation of the growth rates as a function of the
time interval. A clear decay is found that becomes slightly negative after 5
years. Thus, predictability is possible, but for no longer than 5 years.
Right: Empirical cumulative distribution of total sold cars per model (gray
dots) compared with the MaxEnt prediction [see Eq. (8)] (black line).

Another relevant question regarding growth rates becomes legitimate at this
point. It is crucial for predictability: Are car sales a non-Markovian
process? In other words, are the sales of today correlated with those of
yesterday, and thus predictable? Buying a car is an individual decision,
undertaken according to a plethora of various considerations. We are looking
at a complex decision-process that takes its time and relies on decisions
taken previously by other individuals. As we did for cities’ populations in
Refs. cities3, ; memory,  which unraveled the existence of a “memory” in
cities, we have attempted to define here the putative system’s memory as the
time correlation $c(\tau)=\text{Cor}[\dot{x}(t),\dot{x}(t+\tau)]$, where
$\tau$ is the correlation time interval. We find an exponential-like decay,
with slightly negative values after $\tau=5$ years as shown in Fig. 4. This in
turn establishes the limits of predictability: no accurate prediction on sales
can be done for longer than that time-period.

We continue now with the step 3 outlined on page 3.. For maximizing the
entropy $S$, it was shown in Ref. logistic,  that the only observable and
objective constraints for systems with logistic growth are, due to the
limitation in units of resources, the total available units (in our case,
total number of sold cars $X$) and the number of elements sharing these units
(for us, the number of available automobile models $N$). Considering only
those car-models the sales of which exceeded $x_{0}=130$ units in the ten
years (the limit of proportional growth as derived before) we have
$X=9\,957\,537$ and $N=3\,084$. Following Refs. firms, ; entropy1, ; entropy2,
; entropy3, , we write our thermodynamic potential as:

$\Omega=-T\,S-\mu N+\Lambda X$ (4)

where $T$, $\mu$, and $\Lambda$ are the concomitant Lagrange multipliers for
the general variational problem of the thermodynamic potential
$\Omega(S,N,X)$. For a general equation of motion with arbitrary $q$, we
assume that an underlying probability density $p(u)$ governs our process,
where $u=\log_{q}(x/x_{0})$, $N=\int_{0}^{\infty}dup(u)$ and
$X=\int_{0}^{\infty}dup(u)\exp_{q}(u)$. We explicitly cast the MaxEnt problem
as a $p-$variational one:

$0=\delta_{p}\Omega=\delta_{p}\\!\\!\int_{0}^{\infty}\\!\\!\\!dup(u)\left\\{T\log[p(u)/N]-\mu-\Lambda\exp_{q}(u)\right\\}$
(5)

where $\delta_{p}$ represent variations with respect to $p(u)$. We find the
solution $p(u)=zN\exp[-\Lambda^{*}\exp_{q}(u)]$ which in terms of the variable
$x$ is a power law with exponential cut-off:

$p(x)=zN\frac{\exp(-\Lambda^{*}x/x_{0})}{x^{q}},$ (6)

where $\Lambda^{*}=\Lambda/T$ and $z=\exp(\mu/T)$. The first constraint
normalizes $p(u)$ to the number of car models $N$, and the second to the total
car sales, yielding the equations of state:

$\displaystyle z$
$\displaystyle=\left(x_{0}\Lambda^{*}\right)^{1-q}/\Gamma(1-q,\Lambda^{*}),$
(7a) $\displaystyle x_{0}N/X$
$\displaystyle=\Lambda^{*}\Gamma(1-q,\Lambda^{*})/\Gamma(2-q,\Lambda^{*}).$
(7b)

The cumulative distribution can be written as

$P(x)/N=1-\frac{\Gamma(1-q,\Lambda^{*}x/x_{0})}{\Gamma(1-q,\Lambda^{*})}.$ (8)

Passing now to step 4 of page 1, we show in Fig. 4 the comparison of the
MaxEnt prediction with the empirical cumulative distribution of total car
sales. We have considered only sizes larger than $x_{0}=130$ and fit Eq. (8)
via $\log(x_{0})$, $q$, and $\log(\Lambda^{*})$ (logarithms are used for
numerical stability). We find a remarkable fit with $\log(x_{0})=4.884\pm
0.001$, $q=1.048\pm 0.001$, and $\log(\Lambda^{*})=-4.560\pm 0.003$. Our R2
value is $0.99993$ —the relative difference with the empirical values
$\log(130)$, $1$, and $-4.8712$ (taken from the values of $q$, $X$ and $N$)
are $0.3\%$, $5\%$, and $8.8\%$, respectively. Thus, we can safely claim that
the dynamics of car sales are very close to equilibrium and the distribution
of total sales obeys the Maximum Entropy Principle for scale-free systems.

We finally proceed to step 5 on page 5. by defining the equations needed for
performing microscopic simulations of the empiric system. For writing these
equations, we need the following considerations:

  1. 1.

Geometric Brownian walkers are known to obey Gibrat’s law and reproduce states
of dynamical equilibrium following entropic laws. However, Brownian walkers
are Markovian. We will use instead a Langeving-like equation with a viscous
termthermo2 ; thermo3 , which is known to reproduce time correlations with
exponential decay.

  2. 2.

The constraints in the thermodynamic potential of Eq. (4) are introduced as
forces in the Langevin equation, in analogy to what is done in molecular
simulations.

  3. 3.

In addition to the proportional growth $q=1$, the model should include linear
forces ($q=0$) and finite size fluctuations ($q=1/2$). We include them as
independent thermal baths with no inertia.

  4. 4.

The number of sold cars can grow, but does not decrease (un-selling cars is
not allowed). Thus, a Maxwell demonthermo3  is introduced to select only
positive values from the fluctuations;

Figure 5:  Left: $\Lambda^{*}$ as measured from the equilibrium distribution
versus the value of $\Lambda/T_{1}$ used in the simulations. Each dot is an
independent realization using the range of values for $\Lambda$ and $T_{\eta}$
described in the text. Colors represent realizations with the same value of
$\Delta t$.

After these considerations, we write the full set of microscopic equations as:

$\displaystyle\dot{v}_{1,i}(t)$ $\displaystyle=\eta_{i}(t)-\Lambda
x_{i}(t)-\gamma v_{1,i}(t),$ (9) $\displaystyle\dot{x}_{i}(t)$
$\displaystyle=v_{1,i}(t)x_{i}(t)+v_{1/2,i}(t)\sqrt{x_{i}(t)}+v_{0,i}(t)\vee
0,$

where $\gamma$ is the dumping or viscosity which controls the inertia, and the
thermal bath obeys
$\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2T_{\eta}\gamma~{}\delta_{ij}\delta(t-t^{\prime})$
with $T_{\eta}$ defining its temperature. The other two sources of noise are
described as White Noise with $\langle
v_{1/2,i}(t)v_{1/2,j}(t^{\prime})\rangle=2T_{1/2}/\Delta
t~{}\delta_{ij}\delta(t-t^{\prime})$, and $\langle
v_{0,i}(t)v_{0,j}(t^{\prime})\rangle=2T_{0}/\Delta
t~{}\delta_{ij}\delta(t-t^{\prime})$ with $T_{1/2}$, and $T_{0}$ their
respective temperatures. Here, $\delta$ is the Dirac’s Delta and $\Delta t$ is
the discretized interval of time used to numerically solve our equations.
Finally, the symbol $\vee$ indicates choosing the maximum between the quantity
on the left and on the right, preventing negative values of the growth
$\dot{x_{i}}$.

Since our aim is to prove that the entropic procedure properly predicts the
equilibrium distribution in the proportional regime, we have reduced Eqs. (9)
for the sake of simplicity by removing the finite size term $q=1/2$ (as done
in Ref. firms, ) and considering over-dumping dynamics by setting
$\gamma=1/\Delta t$. (For the interested reader, a full exhaustive exercise
solving and exploring Eqs. (9) will be published elsewhere.)

We have solved the simplified version of Eqs. (9) for the range of values
$0.1\leq T_{\eta}\leq 500$ and $10^{-4}\leq\Lambda\leq 0.5$ with $N=3000$,
$T_{0}=100$, and intervals of times in a range of $10^{-3}\leq\Delta t\leq
10^{-2}$, starting at an initial configuration where every walker has
$x(0)=1$. Eventually, the growth rate of the walker becomes zero due to the
force induced by the $\Lambda$ term, reaching saturation. So as to keep the
number of active walkers constant, we add a new one at $x=1$ every time an
older walker saturates. After some time-intervals have elapsed, the system
reaches a thermodynamical equilibrium.

We study the scalable growth by selecting only walkers with values
$x_{i}(t)>100$ and measuring the effective temperature as the variance of the
growth rates $T_{1}=\text{Var}[\dot{x}/x]$ at equilibrium. We find that the
system’s temperature and that of the thermal bath obey $\log(T_{1})=(0.52\pm
0.06)+(0.65\pm 0.02)\log(T_{\eta})$. The equilibrium distribution $p(x)dx$
fits the form of the MaxEnt prediction in Eq. (8): when fitting via $q$ and
$\Lambda^{*}$, we find $q=1.01\pm 0.01$ and $\log(\Lambda^{*})=(3.0\pm
0.2)+\log(\Lambda/T_{1}$) —independently of the other parameters $\Delta t$
and $T_{0}$, see Fig. 5— confirming the correctness of our entropic procedure.
We show in Fig. 5 the cumulative distribution for a simulation using the
empirical temperature $T_{1}=\exp(-3.41)$ ($T_{\eta}=0.0024$) and
$\Lambda=1.35\times 10^{-5}$, obtaining $\Lambda^{*}=0.01$.

Summing up, we have shown that the dynamics underlying car sales corresponds
to a non-Markovian logistic growth with scale symmetry, and the dynamical
equilibrium obeys the maximum entropy principle. We have also found that the
sales become scalable up to a 95% when a model reaches a 0.1% of the monthly
market share. Numerical experiments reproduce indeed the macroscopic features
of the empirical system and confirm the theoretical procedure proposed here.
In view of such results, our macroscopic theoretical framework should be
recommendable in any simulation of car sales that aims to be predictive.
Today, the industry uses algorithms based on, for example, deep learning on
time series that are very useful to fit internal correlations between every
involved variable. Our approach can be used to constrain the degrees of
freedom in such methods, which should further improve their efficiency. We
hope that accurate forecasts will help to describe controlled evolution of the
market and to reduce the economic risk inherent to profound technological
transitions, such as the one that the automobile market is nowadays
experiencing.

## References

  * (1) R. Hanna, F. Kuhnert, and H. Kiuchi, Re-inventing the wheel: Scenarios for the transformation of the automotive industry, PwC, 2015. (www.pwc.com/auto)
  * (2) H.W. Kaas, et al., Automotive revolution – perspective towards 2030, Advanced Industries, McKinsey & Company, January 2016.
  * (3) M. Turner, C. Hamilton, J. Lennon, V. Lloyd, L. Zhao, and I. Roper, Commodities Comment, Macquarie Research, January 2017.
  * (4) Automotive industry crisis of 2008-10. In Wikipedia. Retrieved May 17, 2017, from en.wikipedia.org/wiki/Automotive_industry_crisis
  * (5) A. Pentland, Social Physics: How Good Ideas Spread-the Lessons from a New Science, Penguin Press, 2014.
  * (6) J.L. McCauley. Dynamics of Markets: Econophysics and Finance. Cambridge: Cambridge University Press. 2004.
  * (7) J. Kemeny and J. L. Snell, Mathematical Models in the Social Sciences. Cambridge: MIT Press, 1978.
  * (8) E. Arcaute, E. Hatna, P. Ferguson, H. Youn, A. Johansson, M. Batty, J. Roy. Soc. Interface 12, 102 (2014).
  * (9) L. Bettencourt, G. West. Nature 467, 912–913 (2010).
  * (10) A. Hernando, R. Hernando, A. Plastino, E. Zambrano, J. Roy. Soc. Interface 12, 102 (2015).
  * (11) E. Zambrano, A. Hernando, A. Fernández-Bariviera, R. Hernando, A. Plastino, J. Roy. Soc. Interface 12, 112 (2015).
  * (12) J. R. Banavar, A. Maritan, and A. Rinaldo. Nature 399, 130-132 (1999).
  * (13) M. Barthelemy and A. Flammini, Phys. Rev. Lett. 100, 138702 (2008).
  * (14) M.C. González, C.A. Hidalgo, and A.L. Barabási, Nature 453, 779–782 (2008).
  * (15) A. Hernando and A. Plastino, Phys. Lett. A 377, 176 (2013).
  * (16) A.-L. Barabasi, R. Albert, Rev. Mod. Phys. 74, 47 (2002).
  * (17) J.P. Gleeson, J.A. Ward, K.P. O’Sullivan, and William T. Lee , Phys. Rev. Lett. 112, 048701 (2014).
  * (18) For detailed microdata on registrations: https://sedeapl.dgt.gob.es/WEB_IEST_CONSULTA (Directorate-General of Traffic, Ministry of the Interior, Spain). For comprehensive analysis: www.anfac.com (Spanish Organization of Car and Truck Builders ANFAC).
  * (19) Aggregated data can be delivered upon request from IHS Markit (https://ihsmarkit.com).
  * (20) M. Gell-Mann and C. Tsallis, Eds. Nonextensive Entropy: Interdisciplinary applications, Oxford University Press, Oxford, 2004; C. Tsallis, Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World, Springer, New York, 2009.
  * (21) H. Rozenfeld, et al., _Proc. Nat. Acad. Sci._ 105, 18702 (2008).
  * (22) A. Hernando, A. Plastino, A.R. Plastino, Eur. Phys. J. B 85, 147 (2012)
  * (23) A. Hernando, A. Plastino, Eur. Phys. J. B 85, 293 (2012).
  * (24) A. Hernando, R. Hernando, A. Plastino, A.R. Plastino, J. R. Soc. Interface 10, 20120758 (2013).
  * (25) A. Hernando and A. Plastino, Phys. Rev. E 86, 066105 (2012).
  * (26) A. Hernando, R. Hernando, A. Plastino, J. R. Soc. Interface 11, 20130930 (2014).
  * (27) F. Reif, Fundamentals of statistical and thermal physics. Waveland Press. 2008.
  * (28) R. Balian, From microphysics to macrophysics, vols I and II. Springer. 2006.


11institutetext: Center for Data Science, Peking University, Beijing, China
22institutetext: Canon Medical Systems, Beijing, China 33institutetext: School
of Mathematical Sciences, Peking University, Beijing, China
33email: {xuju, mcmong<EMAIL_ADDRESS>
# Automatic Data Augmentation for 3D Medical Image Segmentation

Ju Xu 1 Equal contributions.1 Equal contributions. Mengzhang Li 112
11footnotemark: 12 11footnotemark: 1 Zhanxing Zhu 3311

###### Abstract

Data augmentation is an effective and universal technique for improving
generalization performance of deep neural networks. It could enrich diversity
of training samples that is essential in medical image segmentation tasks
because 1) the scale of medical image dataset is typically smaller, which may
increase the risk of overfitting; 2) the shape and modality of different
objects such as organs or tumors are unique, thus requiring customized data
augmentation policy. However, most data augmentation implementations are hand-
crafted and suboptimal in medical image processing. To fully exploit the
potential of data augmentation, we propose an efficient algorithm to
automatically search for the optimal augmentation strategies. We formulate the
coupled optimization w.r.t. network weights and augmentation parameters into a
differentiable form by means of stochastic relaxation. This formulation allows
us to apply alternative gradient-based methods to solve it, i.e. stochastic
natural gradient method with adaptive step-size. To the best of our knowledge,
it is the first time that differentiable automatic data augmentation is
employed in medical image segmentation tasks. Our numerical experiments
demonstrate that the proposed approach significantly outperforms existing
build-in data augmentation of state-of-the-art models.

###### Keywords:

Medical Image Segmentation Data Augmentation AutoML

##  1 Introduction

In the past few years, deep neural network has achieved incredible progress in
medical image segmentation tasks and promoted booming development of computer
assisted intervention. This has benefitted research and clinical treatment of
disease diagnosis, treatment design and prognosis evaluation [13, 14]. Given
the training data, researchers proposed various 2D/3D medical image
segmentation models for supervised or semi-supervised tasks [8, 6]. However,
the performance of deep learning models heavily relies on large scale well-
labeled data. Currently, data augmentation is a widely used and effective
technique to increase the amount and diversity of available data, and thus
improving models’ generalization performance. In the domain of natural image
processing, typical data augmentation strategies include _manually_ cropping,
rotating or adding random noise to the original images. Besides thess ad-hoc
approaches, generative models [7] and unsupervised learning models [15] are
also employed for generating extra data. Unfortunately, those augmentation
techniques might not be optimal for a specific task, and thus the customized
data augmentation strategy is required. Recently, researchers proposed to
search the augmentation policy by reinforcement learning [5] or density
matching [10], inspired by previously works of automatic machine learning
(AutoML) on neural architecture search (NAS [11, 12, 18]).

For medical image segmentation tasks, data augmentation techniques are also
used in UNet and its variants nnUNet [8], R2U-Net [2], etc. However, these
methods are simple and hand-made, and the improvement of segmentation accuracy
is limited. In [16], the authors proposed to utilize reinforcement learning to
search for augmentation strategies. However, it costs 768 GPU hours and it
only searchs the probability of each augmentation strategy in [16]. Moreover,
the difference between natural and medical images such as spatial contextual
correlation, smaller scale of dataset and unique pattern of specified organs
or tumor makes the augmentation strategies adopted in natural images difficult
be transferred to medical domains.

In this paper, we propose an _automatic data augmentation_ framework (ASNG)
through searching the optimal augmentation policy, particularly for 3D medical
image segmentation tasks. It’s the first automatic data augmentation work in
whole semantic segmentation filed. The contributions of our paper are as
follows:

  * •

It’s the first time that we formulate the auto-augmentation problem into a bi-
level optimization problem and apply an approximate algorithm to solve it

  * •

The designed search space in medical image field is novel. Different from
previous methods which searched for a fixed magnitude of operations, we search
for an interval of magnitude

  * •

Different from previous method which searched for a fixed augmentation
strategy, the searched augmentation strategy of our method is dynamically
changing during the training. Besides, we don’t need to retrain the target
network after the searching process

  * •

Experiments demonstrate that our ASNG can indeed achieve the SOTA of the
performance

##  2 Method

In our method, we formulate the problem of finding the optimal augmentation
policy as a discrete search problem. Our method consists of two components:
the designed of search space and search algorithm. The search algorithm
samples a data augmentation policy $S$ from the search space consisting of
proposed operations, and then decides the magnitude of the operation and the
probability of applying this operation. The framework of our method can be
seen in Fig.1. We will elaborate the two components in the following.

Figure 1: The framework of our proposed method. $D_{train}$, $D_{val}$
represent training dataset, validation dataset, respectively. $p_{\theta}$ is
the distribution of $c$.

###  2.1 Search Space of Data Augmentation for 3D Medical Images

Since it is the first work for applying AutoAugment strategies in medical
image area, we have to design the search space for our ASNG algorithm. In our
search space, a policy consists of seven image operations to be applied in a
sequential manner. Each image operation is associated with two
hyperparameters: 1) the probability of applying this operation, and 2) the
interval of magnitude for the image operation.

The seven image operations we used in our experiments are from
batchgenerators, a pupular Python image library111https://github.com/MIC-
DKFZ/batchgenerators, including Scale, RoateX/Y/Z, Alpha (magnitude of the
elastic deformation), Sigma (scale of the elastic deformation), Gamma (same as
gamma correction in photos or computer monitors). In order to increase the
diversity of augmentation policies, we do not fix a specific magnitude for an
operation like previous works [5], but set an interval of magnitude for an
operation. Therefore we should decide the left boundary of interval (LB) and
the right one (RB). To decrease the search complexity, we discretize the range
of magnitude into 11 values with uniform spacing so that we can use a discrete
search algorithm to find them. Besides the magnitude of transformation
operation, we also search for the probability of conducting these
transformations, i.e. the probability of applying scale transformation,
rotation, gamma transformation, and elastic deformation, denoted as
$p_{scale}$, $p_{rot}$, $p_{gamma}$, $p_{eldef}$, respectively. Similarly, we
also discretize the probability of applying that operation into 11 values with
uniform spacing. Table 1 summarizes the range of magnitudes and possibilities
for the seven operations. Fig. 2 shows one example of augmented image and
label based on our method, in which the image transformations are from the
defined search space.

We can easily observe that naively searching one augmentation strategy becomes
a search problem with $11^{11}$ possibilities. The search space is so huge
that an efficient algorithm is required, as proposed in the following.

Table 1: The range of parameters in strategies we will search. Operation |  LB |  RB |  Probability | Range
---|---|---|---|---
Scale |  [0.5, 1.0] |  [1.0, 1.5] |  $p_{scale}$ | [0, 1]
RotationX |  [$\frac{-\pi}{6}$, 0] |  [0, $\frac{-\pi}{6}$] |  $p_{rot}$ | [0, 1]
RotationY |  [$\frac{-\pi}{6}$, 0] |  [0, $\frac{-\pi}{6}$] |  $p_{rot}$ | [0, 1]
RotationZ |  [$\frac{-\pi}{6}$, 0] |  [0, $\frac{-\pi}{6}$] |  $p_{rot}$ | [0, 1]
Alpha |  [0, 450] |  [450, 900] |  $p_{eldef}$ | [0, 1]
Sigma |  [0, 7] |  [7, 14] |  $p_{eldef}$ | [0, 1]
Gamma |  [0.5, 1] |  [1, 1.5] |  $p_{gamma}$ | [0, 1]

###  2.2 Stochastic Relaxation Optimization of Policy Sampling

We denote $f(w,c)$ as the objective function, where $w\in\mathcal{W}$ are
network parameters and $c\in\mathcal{C}$ are data augmentation strategies.
$f_{train}$ and $f_{val}$ are the training and the validation loss,
respectively. Both losses are determined not only by the augmentation policy
$c$, but also the weight $w$. The goal for augmentation strategy search is to
find $c^{*}$ that minimizes the validation loss $f_{val}(w^{*},c^{*})$, where
the weights are obtained by minimizing the training loss
$w^{*}=\text{argmin}_{w}f_{train}(w,c^{*})$. Thus augmentation strategy search
is a bi-level optimization problem, we can write the problem as follows:

$\displaystyle\min_{c}f_{val}(w^{*}(c),c)$ (1) $\displaystyle s.t.\quad
w^{*}(c)=\mathop{\text{argmin}}\limits_{w}f_{train}(w,c)$ (2)

Solving the above problem is not easy, since we cannot obtain the gradient
w.r.t. $c$, thus it is hard to optimize $c$ via gradient descent. Though
simple grid search or reinforcement learning proposed in [5] can be utilized
to search for $c$, the computational cost is extremely high if we evaluate the
performance of every $c$. To this end, we propose to solve this optimization
problem efficiently first by stochastic relaxation and then applying natural
gradient descent [3], as described in the following.

Figure 2: Visualization of the proposed data augmentation on certain 2D
section of Task01 BrainTumour dataset. Top: original image and label, Bottom:
the image and label generated by the searched data augmentation policy, in
which the operations are from Table 1.

####  2.2.1 Stochastic Relaxation

We turn the original optimization problem into an optimization of
differentiable objective $J$ by stochastic relaxation [1]. The basic idea of
stochastic relaxation is: instead of directly optimizing w.r.t $c$, we
consider a distribution $p_{\theta}(c)$ over $c$ parametrized by $\theta$, and
minimize the expected value of the validation loss $f_{val}$ w.r.t $\theta$,
i.e.,

$\displaystyle\min\limits_{\theta}J(\theta)=\int_{c\in\mathcal{C}}f_{val}(w^{*}(c),c)p_{\theta}(c)dc$
(3) $\displaystyle s.t.\quad
w^{*}(c)=\mathop{\text{argmin}}\limits_{w}f_{train}(w,c)$ (4)

The stochastic relaxation makes $J$ differentiable w.r.t both $w$ and
$\theta$. Therefore we can update $w$ and $\theta$ by gradient descent.
However, the gradient $\nabla_{w}J(w,\theta)$ is not tractable because we can
not evaluate the mean performance of $c\in\mathcal{C}$ in a closed-form way.
Here we estimate the gradient by Monte-Carlo (MC) using
$\nabla_{w}J(w^{t},c_{i})$ with i.i.d. samples $c_{i}\sim
p_{\theta^{t}}(c),\,i=1,\dots,N_{w}$, namely:

$\displaystyle
G_{w}(w^{t},\theta^{t})=\frac{1}{N_{w}}\sum\limits_{i=1}^{N_{w}}\nabla_{w}f_{train}(w^{t},c_{i})$
(5)

Now we can approximate $\nabla_{w}J(w,\theta)$ with the stochastic gradient
$G_{w}(w^{t},\theta^{t})$, $w^{t}$ can be updated as

$\displaystyle w^{t+1}=w^{t}-\epsilon_{w}G_{w}(w^{t},\theta^{t}),$ (6)

where $\epsilon_{w}$ is the learning rate for network parameters. Due to that
the distance between two probability distribution is not Euclidean, updating
$\theta$ directly by gradient descent like $w$ is not appropriate. We then
resort to natural gradient (NG [3]) designed for parametric probability
distributions,

$\displaystyle\theta^{t+1}=\theta^{t}-\epsilon_{\theta}F(\theta_{t})^{-1}\nabla_{\theta}J(w,\theta),$
(7)

where $F(\theta_{t})$ is the Fisher matrix, $\epsilon_{\theta}$ is the
learning rate. The probability distribution we utilized for $c\in\mathcal{C}$
is multinomial distribution. How to calculate the Fisher matrix can be seen in
[1]. Similiar with [1], we utilize adaptive step-size $\epsilon_{\theta}$ to
make the learning process faster. Monte-Carlo is also adopted to approximate
$\nabla_{\theta}J(w,\theta)$, and then

$\displaystyle\theta^{t+1}=\theta^{t}-\epsilon_{\theta}F(\theta_{t})^{-1}\frac{1}{N_{\theta}}\sum\limits_{j=1}^{N_{\theta}}\nabla_{\theta}f_{val}(w^{t+1},c_{j})\ln
p_{\theta}(c_{j})$ (8)

We summarize the procedure of our proposed approach in Algorithm 1.

Algorithm 1 ASNG

1: Input:  $w^{0}$, $\theta^{0}$, $\epsilon_{w}$. $\epsilon_{\theta}$,
$N_{w}$, $N_{\theta}$

2: Input:  Training dataset $D_{train}$, validation dataset $D_{val}$, test
dataset $D_{test}$.

3: for i=1 to epoch do

4: for t=1 to T do

5: Generate $N_{w}$ policys according to $p_{\theta_{t}}$

6: Augment training data from $D_{train}$ with $N_{w}$ policys, respectively;

7: Obtain the loss $f_{train}(w_{t},c_{i})$ ($i=1,\dots,N_{w}$) on
$D_{train}$;

8: Update $w_{t}$ according to Equation 6, then obtain $w_{t+1}$;

9: Generate $N_{\theta}$ policys according to $p_{\theta_{t}}$;

10: for j=1 to $N_{\theta}$ do

11: Augment training data according to policy $c_{j}$;

12: Update $w_{t}$ to obtain $\hat{w_{t}}$;

13: Obtain the validation loss $f_{val}(\hat{w_{t}})^{j}$ on $D_{val}$;

14: Restore the network parameters, $\hat{w_{t}}=w_{t}$;

15: end for

16: Utilize validation loss $f_{val}(\hat{w_{t}})^{j}$, policys $c_{j}$
($j=1,\dots,N_{\theta}$) to update $\theta_{t}$ according to equation 8;

17: end for

18: end for

19: Test the network on $D_{test}$;

20: return final networks.

##  3 Implementation and Experiments

###  3.1 Datasets and Implementation Details

Datasets: We conduct the proposed method on three 3D segmentation tasks used
in the medical segmentation decathlon challenge
(MSD222http://medicaldecathlon.com/): (1) Task01 Brain Tumour (484 labeled
images, 3 classes), (2) Task02 Heart (20 labeled images, 1 class) and (3)
Task05 Prostate (32 labeled images, 2 classes). Each dataset is collected for
a specified task, their various input sizes, voxel spacings and foreground
targets are suitable for demonstrating our algorithm’s generalization. We
evaluate the performance with 5-fold cross validation as [4, 9] since the
ground truth labels for test dataset are not publicly available.
Compared Methods include 3D U-ResNet [17], SCNAS [9], nnUNet without data
augmentation (nnUNet NoDA) and 3D nnUNet [8] with default data augmentation
strategy333https://github.com/MIC-DKFZ/nnUNet/. 3D U-ResNet utilizes residual
blocks and attention gates; and SCNAS is the latest neural architecture search
model for 3D medical image segmentation, which applies a scalable gradient-
based optimization to find the optimal model architecture. The method proposed
in SCNAS can’t utilized for differentiable autoaugmentation strategies search.
In SCNAS, a mixed operation is created by adding all these operations in
search space based on the importance of each operation. However, we can’t add
the transformation results of each augmentation strategy. There we don’t apply
the proposed method of SCNAS to our augmentation strategies search tasks. Note
that the code of nnUnet has already implemented random augmentation. In
nnUnet, the magnitude of operation is sampled from a predefined interval in
every training epoch. AutoAugment [5] costs 5000 GPU hours to search for a
policy. FastAutoAugment [10] needs to spilt the training dataset as K folds.
However, dataset in medical image area is quite small. Training model in a
small dataset will overfit. Therefore, we don’t compare our method with
AutoAugment and FastAutoAugment. The prediction result is inferenced using a
sliding window with half the patch size ensuring 50% overlapping, i.e., each
voxel is inferenced at least two times at test.
Implementation Details We preprocess the data with same pipeline used in 3D
nnUNet. We unify the identical voxel spacing values by proper interpolation
due to different spacing values of each case, i.e. resampling them to 0.7 mm
$\times$ 0.7 mm $\times$ 0.7 mm firstly. We apply Z-score normalization of
voxel value for each input channel separately; and grip the input patch whose
size is set as 128 $\times$ 128 $\times$ 128, and its foreground ratio is set
larger than 1/3 ensuring UNet variants could learn features of foreground.

Following the default data augmentation policy in nnUNet, we utilize scaling,
rotation, elastic and gamma transformation both in training and test. The
parameters of the operation and probability of conducting that operation are
both within search space of our ASNG algorithm.

The code is implemented using PyTorch 1.0.0. The ADAM optimizer is utilized
for training where the learning rate and weight decay are initialized as
$3\times 10^{-4}$ and $3\times 10^{-5}$, respectively, where it is reduced by
80% if the training loss is not reduced over 30 epochs. Besides ASNG, the
training process of other benchmarks would last for 500 epochs if the learning
rate is larger than $10^{-7}$. Following [9] and [8], the loss function for 3D
U-ResNet and SCNAS is Jaccard distance, for nnUNet and ASNG is sum of minus
Dice similarity and Cross Entropy. Considering the training time, ASNG is
trained for 50, 200 and 200 epochs on Brain Tumour, Heart and Prostate,
respectively, with batch size of 2. It takes about 10 days on one NVIDIA TITAN
RTX GPU, compared with that one integrated nnUNet training procedure takes
about 3 days. The sampling times $T$ of ASNG is 2 because of the limited
memory of GPU, though larger $T$ could produce better numerical results. Our
codes can be found here 444https://github.com/MengzhangLI/ASNG.

###  3.2 Experimental Results

Our result is shown in Table 2. Because of unavailable labels of test set and
restricted online submission times, those Auto ML models on 3D Medical Image
Segmentation tasks are all evaluated on validation set. In this paper we still
follows this metric for fair comparison. ASNG outperforms other architectures
especially 3D nnUNet, which is the best medical image segmentation framework
with default data augmentation. It should be noted that since Heart and
Prostate only have 20 and 32 labeled images, in [9] the obtained architecture
of SCNAS based on the first fold of 484 labeled Brain Tumour images is
transferred to Heart and Prostate tasks to avoid overfitting. Remarkably, our
method, applied only on the basic network architecture, could still achieve
best prediction accuracy. This clearly demonstrates the necessity and
effectiveness of data augmentation policy search in 3D medical image
segmentation.

Figure 3 shows the example of segmentations results and validation loss w.r.t.
number of epochs in the Prostate task. We can observe that our method ASNG can
produce better prediction and achieve more stable improvement during training
than other compared methods.

In this paper [16], the proposed method utilizes reinforcement learning to
search for augmentation strategies, which costs 768 GPU hours while our method
costs less than 100. And their result (Dice 0.92) on task 02 is worse than
ours (0.933). Besides, the first paper only searched the probability of each
augmentation strategy while our method not only search the probability but
also the magnitude.

Figure 3: Results on Task05 Prostate of selected architectures. Left: Example
of inference, green mask represents peripheral zone and red mask represents
transitional zone. Right: The trend of loss on validation set.

Label | Brain Tumour | Heart | Prostate
---|---|---|---
Edema |  | Non-
---
Enhancing
Enhancing | Average | Left atrium | Peripheral | Transitional | Average
U-ResNet | 79.10 | 58.38 | 77.37 | 71.61 | 91.48 | 48.37 | 79.17 | 63.77
nnUNet NoDA | 81.27 | 60.92 | 77.90 | 73.36 | 92.85 | 58.61 | 83.61 | 71.11
nnUNet | 81.68 | 61.29 | 77.97 | 73.65 | 93.21 | 63.14 | 86.53 | 74.84
SCNAS | 80.41 | 59.85 | 78.50 | 72.92 | 91.91 | 53.81 | 82.02 | 67.92
ASNG | 81.94 | 61.85 | 79.35 | 74.38 | 93.27 | 67.40 | 87.05 | 77.22

Table 2: Average Dice similarity coefficients (%) for Brain tumor, Heart, and
Prostate 3D segmentation tasks of MSD.

##  4 Conclusion

We have proposed an automatic data augmentation strategy to accommodate 3D
medical image segmentation tasks. By configuring proper search space followed
by gradient-based optimization, the customized data augmentation strategy for
each task could be obtained. The numerical results for different segmentation
tasks show it could outperform the state-of-the-art models that are widely
used in this area. Furthermore, the proposed approach shows that, compared
with searching network architectures, searching for optimal data augmentation
policy is also important. As for future work, designing better search space
and accelerating the search process can be considered.

## References

  * [1] Akimoto, Y., Shirakawa, S., Yoshinari, N., Uchida, K., Saito, S., Nishida, K.: Adaptive stochastic natural gradient method for one-shot neural architecture search. arXiv preprint arXiv:1905.08537 (2019)
  * [2] Alom, M.Z., Hasan, M., Yakopcic, C., Taha, T.M., Asari, V.K.: Recurrent residual convolutional neural network based on u-net (r2u-net) for medical image segmentation. arXiv preprint arXiv:1802.06955 (2018)
  * [3] Amari, S.I.: Natural gradient works efficiently in learning. Neural computation 10(2), 251–276 (1998)
  * [4] Bae, W., Lee, S., Lee, Y., Park, B., Chung, M., Jung, K.H.: Resource optimized neural architecture search for 3d medical image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 228–236. Springer (2019)
  * [5] Cubuk, E.D., Zoph, B., Mane, D., Vasudevan, V., Le, Q.V.: Autoaugment: Learning augmentation strategies from data. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 113–123 (2019)
  * [6] Ganaye, P.A., Sdika, M., Triggs, B., Benoit-Cattin, H.: Removing segmentation inconsistencies with semi-supervised non-adjacency constraint. Medical image analysis 58, 101551 (2019)
  * [7] Huang, S.W., Lin, C.T., Chen, S.P., Wu, Y.Y., Hsu, P.H., Lai, S.H.: Auggan: Cross domain adaptation with gan-based data augmentation. In: Proceedings of the European Conference on Computer Vision (ECCV). pp. 718–731 (2018)
  * [8] Isensee, F., Petersen, J., Kohl, S.A., Jäger, P.F., Maier-Hein, K.H.: nnu-net: Breaking the spell on successful medical image segmentation. arXiv preprint arXiv:1904.08128 (2019)
  * [9] Kim, S., Kim, I., Lim, S., Baek, W., Kim, C., Cho, H., Yoon, B., Kim, T.: Scalable neural architecture search for 3d medical image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 220–228. Springer (2019)
  * [10] Lim, S., Kim, I., Kim, T., Kim, C., Kim, S.: Fast autoaugment. In: Advances in Neural Information Processing Systems. pp. 6662–6672 (2019)
  * [11] Liu, H., Simonyan, K., Yang, Y.: Darts: Differentiable architecture search. arXiv preprint arXiv:1806.09055 (2018)
  * [12] Pham, H., Guan, M.Y., Zoph, B., Le, Q.V., Dean, J.: Efficient neural architecture search via parameter sharing. arXiv preprint arXiv:1802.03268 (2018)
  * [13] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [14] Tajbakhsh, N., Shin, J.Y., Gurudu, S.R., Hurst, R.T., Kendall, C.B., Gotway, M.B., Liang, J.: Convolutional neural networks for medical image analysis: Full training or fine tuning? IEEE transactions on medical imaging 35(5), 1299–1312 (2016)
  * [15] Xie, Q., Dai, Z., Hovy, E., Luong, M.T., Le, Q.V.: Unsupervised data augmentation. arXiv preprint arXiv:1904.12848 (2019)
  * [16] Yang, D., Roth, H., Xu, Z., Milletari, F., Zhang, L., Xu, D.: Searching learning strategy with reinforcement learning for 3d medical image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 3–11. Springer (2019)
  * [17] Yu, L., Yang, X., Chen, H., Qin, J., Heng, P.A.: Volumetric convnets with mixed residual connections for automated prostate segmentation from 3d mr images. In: Thirty-first AAAI conference on artificial intelligence (2017)
  * [18] Zoph, B., Vasudevan, V., Shlens, J., Le, Q.V.: Learning transferable architectures for scalable image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 8697–8710 (2018)


# Recent scaling properties of Bitcoin price returns

T Takaishi  Hiroshima University of Economics, Hiroshima 731-0192, JAPAN tt-
<EMAIL_ADDRESS>
###### Abstract

While relevant stylized facts are observed for Bitcoin markets, we find a
distinct property for the scaling behavior of the cumulative return
distribution. For various assets, the tail index $\mu$ of the cumulative
return distribution exhibits $\mu\approx 3$, which is referred to as ”the
inverse cubic law.” On the other hand, that of the Bitcoin return is claimed
to be $\mu\approx 2$, which is known as ”the inverse square law.” We
investigate the scaling properties using recent Bitcoin data and find that the
tail index changes to $\mu\approx 3$, which is consistent with the inverse
cubic law. This suggests that some properties of the Bitcoin market could vary
over time. We also investigate the autocorrelation of absolute returns and
find that it is described by a power-law with two scaling exponents. By
analyzing the absolute returns standardized by the realized volatility, we
verify that the Bitcoin return time series is consistent with normal random
variables with time-varying volatility.

##  1 Introduction

In 2008, Nakamoto[1] proposed Bitcoin as the first practical cryptocurrency
based on the blockchain technology. His proposal was quickly accepted, and in
2009, the Bitcoin network was launched as a peer-to-peer payment network.
Since then, a number of digitized cryptocurrencies have been proposed and
created111 e.g., see https://coinmarketcap.com/ for the current market
capitalizations of cryptocurrencies., and growing cryptocurrency markets have
a strong impact on contemporary financial markets. Bitcoin has recently
attracted the interest of researchers and various properties and aspects of
Bitcoin have been investigated. For example, volatility analysis [2, 3], price
clustering[4, 5], adaptive market hypothesis [6], transaction activity [7],
multifractality[8], liquidity and efficiency[9], Taylor effect[10], structural
breaks [11], long memory effects[12], rough volatility[13], power-law cross-
correlation[14], asymmetric multifractal analysis[15], and so forth.

Universal properties across various financial time series are classified as
stylized facts[16], which include (i) fat-tailed distributions, (ii)
volatility clustering, and (iii) slow decay of autocorrelation in absolute
returns, among others. While stylized facts are also observed in Bitcoin, for
example, [17, 8], a distinct property has been observed for the tail index of
return distributions. For other assets such as stocks, it is known that the
tail of the cumulative return distribution is described by a power-law
function and its tail index $\mu$ is obtained as $\mu\approx 3$, known as ”the
inverse cubic law”[18, 19, 20, 21]. For the Bitcoin market before 2014, the
tail index is estimated to be $\mu\approx 2$, which differs from the known
stylized fact for other assets and is referred to as ”the inverse square
law”[22]. Similar tail indices have also been reported in Ref.[23]. Since this
is a newly emerging market, market properties of Bitcoin could change as
trading activity increases over time. In fact, the market efficiency measured
by the Hurst exponent of the return time series varies over time and at the
early stage of the Bitcoin market, the Hurst exponent is observed to be less
than 1/2, which indicates that the time series is anti-persistent[24]. This
anti-persistent behavior could be related to the illiquidity of the Bitcoin
market[9, 25]. As the liquidity of the Bitcoin market improves, the Hurst
exponent moves to 1/2[24, 25].

We investigate the possible change in scaling properties of the return
distribution for the recent Bitcoin market. We also investigate the cumulative
distributions and autocorrelation of the absolute returns and the return
standardized by the realized volatility.

##  2 Data set

The data are Bitcoin tick data (in dollars) traded on Bitstamp from September
11, 2011 to June 21, 2020222Due to a hacking incident, no data are available
from January 4, 2015 to January 9, 2015. For these missing data, we treat them
as the price is unchanged. and downloaded from
Bitcoincharts333http://api.bitcoincharts.com/v1/csv/. We construct 1-min price
data from the Bitcoin tick data and then calculate 1-min returns
$R_{t},t=1,2,...,N$ as

$R_{t}=\log p_{t}-\log p_{t-1},$ (1)

where $p_{t}$ is the 1-min price at $t$ (min). We further standardize the
returns data using

${\bar{R_{t}}}=\frac{R_{t}-\nu}{\sigma}.$ (2)

where $\nu$ and $\sigma$ represent the average and standard deviation of
returns $R_{t}$, respectively. We select the following two data periods from
the data, which separate low and high liquidity periods: I) September 11,
2011-December 31, 2013 and II) January 1, 2015-June 21, 2020. Dataset I (II)
contains low (high) liquidity data[25].

##  3 Results

First, we show the cumulative distributions of returns in Figure 1(a):
positive tail and (b): negative tail, in a log-log plot, and recognize that
the tail behavior is described by a power law. While there is no significant
difference between positive and negative tails, a considerable difference is
seen between data periods, that is, a heavier tail is observed for the
cumulative distribution of dataset (II), which contains high liquidity data.
We obtain tail indices of the cumulative distributions by a regression fit to
a power-law function of $\kappa\bar{R}^{-\mu}$, where $\kappa$ and $\mu$ are
fitting parameters. The results of the tail indices are reported in Table 1.
The tail index for dataset I is $\mu\approx 2$, which is consistent with the
previous results obtained for earlier periods[22]. On the other hand, the tail
index for II is $\mu\approx 3$, which is consistent with the tail index
obtained for other assets such as stocks. For the entire data set, we find a
tail index between I and II. Our findings suggest that the scaling properties
of the return distributions of the Bitcoin market could change over time, For
the recent liquidity of the Bitcoin market, the tail index is in agreement
with the well-known stylized facts for other assets.

Table 1: Results of tail index. | (I) 2011-2013 | (II) 2015-2020 | 2011-2020
---|---|---|---
Fitting region | [2,20] | [2,20] | [1,10]
positive tail index $\mu$ | -2.06 | -3.33 | -2.35
negative tail index $\mu$ | -2.08 | -3.23 | -2.38

Figure 1:  (a) Cumulative distribution for positive tail. (b) Cumulative
distribution for negative tail.

Next, we show the autocorrelation function (ACF) of the absolute returns. The
ACF of a time series $x_{i}:i=1,...,N$ is defined by

$ACF(t)=\frac{\langle(x_{i}-\rho)(x_{i+j}-\rho)\rangle}{\sigma_{x}^{2}}.$ (3)

where $\rho$ and $\sigma_{x}^{2}$ represent the average and variance of
$x_{i}$, respectively, and $\langle O_{i}\rangle$ means taking an average over
$O_{i}$. Figure 2 displays the ACF of the absolute returns in the log-log plot
and we see that the ACF is also well described by a power law. A remarkable
characteristic of the ACF is that the power-law exponent seems to change at
$t\sim 1000$. The power-law exponent for $t>1000$ is found to be larger than
that for $t<1000$. There is no clear explanation for this change. The results
of the power-law exponents are reported in Table 2. Furthermore, we find that
there is no substantial difference for the power-law exponents between the
three data sets. Thus, the long-memory property of Bitcoin remains unchanged
over time.

Table 2:  Results of the power-law exponent. | (I) 2011-2013 | (II) 2015-2020 | 2011-2020
---|---|---|---
fitting region | [3,500] | [3,500] | [3,500]
exponent $\mu$ | -0.121 | -0.116 | -0.120
fitting region | [15000,10000] | [1500,10000] | [1500,10000]
exponent $\mu$ | -0.197 | -0.235 | -0.206
Figure 2:  ACF of absolute returns. The symbols at t=100 and 1000 show typical
one-sigma errors calculated using the jackknife method.

The power law observed in the absolute return means that the time series of
the absolute return has long memory, which sharply contrasts to returns for
which the ACF of the return time series is short ranged. The long memory in
the absolute returns is considered to be related to the time correlation of
volatility. Let us assume that the return time series is described by
$R_{t}\equiv\sigma_{t}\epsilon_{t}$, where $\sigma_{t}$ is the volatility at
$t$ and $\epsilon_{t}$ is the standard normal random variable. Namely, the
return $R_{t}$ is assumed to behave as a normal random variable with time-
varying volatility. This assumption of returns automatically satisfies no
autocorrelation in returns. On the other hand, the autocorrelation of absolute
returns depends on the time structure of $\sigma_{t}$. If we standardize the
(observed) returns $R_{t}$ by $\sigma_{t}$, that is,
$R_{t}/\sigma_{t}\equiv\tilde{R_{t}}$, we expect to obtain the standard normal
random time series $\epsilon_{t}$. To test whether $\tilde{R_{t}}$ satisfies
this expectation we need to estimate the volatility $\sigma_{t}$ of returns in
the proper manner because volatility is not observable in the financial
markets.

In empirical finance, the standard technique to estimate volatility is to use
volatility models such as GARCH-type models[26, 27, 28, 29, 30, 31, 32, 33].
The drawback of using such models is that volatility estimates are model-
dependent, and which estimate is better is not well defined. Recent
availability of high-frequency financial data enable us to obtain a model-free
estimate known as ”realized volatility (RV)”[34, 35] constructed by the sum of
the square of intraday returns. The returns standardized by the RV are
examined for the exchange rate[36, 37] and return[38], and the normality of
the standardized return $\tilde{R_{t}}$ is established. However, detailed
studies report the existence of deviation from the normality caused by the low
sampling frequency in the RV[39, 40, 41, 42, 43]. Furthermore because
$\tilde{R_{t}}$ is expected to be a random variable, we should observe no
autocorrelation not only for $\tilde{R_{t}}$ but also for the absolute
$\tilde{R_{t}}$, that is, $|\tilde{R_{t}}|$. The autocorrelation of the
absolute $\tilde{R_{t}}$ is examined for Japanese stocks[44], and it is
verified that no autocorrelation is observed for the absolute
$\tilde{R_{t}}$444This behavior is also confirmed for returns simulated by the
artificial spin model and standardized by GARCH volatility[45].. We examine
the nonexistence of autocorrelation in the absolute $\tilde{R_{t}}$ of
Bitcoin. Figure 3 shows the autocorrelation of the absolute $\tilde{R_{t}}$
standardized by the RV at a 5-min sampling frequency555In the empirical
analysis, the 5-min sampling frequency for the RV calculations is optimal[46,
47]., and we confirm that there is no significant autocorrelation in the
absolute $\tilde{R_{t}}$.

Figure 3:  ACF of the absolute returns standardized by RV. For the calculation
of the ACF, we take the data from January 10, 2015 to June 21, 2020.

##  4 Conclusions

We investigate the scaling properties of the Bitcoin market and find that the
tail index $\mu$ of the cumulative return distribution changes from
$\mu\approx 2$ to $\mu\approx 3$, which indicates that the recent Bitcoin
market exhibits the well-established scaling law, that is, the inverse cubic
law. Our findings suggest that some properties of the Bitcoin market could
change over time.

We find that the autocorrelation of absolute returns shows a power law with
two scaling exponents separated at approximately $t=1000$min. The ACF of the
absolute returns standardized by RV shows no autocorrelation, which indicates
that the return time series is consistent with a normal random variable with
time-varying volatility.

## Acknowledgment

Numerical calculations for this work were carried out at the Yukawa Institute
Computer Facility and at the Institute of Statistical Mathematics. This work
was supported by JSPS KAKENHI Grant Number JP18K01556.

## References

## References

  * [1] Nakamoto S 2008 Bitcoin: A peer-to-peer electronic cash system https://bitcoin.org/bitcoin.pdf
  * [2] Bouoiyour J and Selmi R 2016 Economics Bulletin 36 1430–1440
  * [3] Dyhrberg A H 2016 Finance Research Letters 16 85–92
  * [4] Urquhart A 2017 Economics Letters 159 145–148
  * [5] Li X, Li S and Xu C 2018 Finance Research Letters
  * [6] Khuntia S and Pattanayak J 2018 Economics Letters 167 26–28
  * [7] Koutmos D 2018 Economics Letters 167 81–85
  * [8] Takaishi T 2018 Physica A 506 507–519
  * [9] Wei W C 2018 Economics Letters 168 21–24
  * [10] Takaishi T and Adachi T 2018 Economics Letters 172 5–7
  * [11] Thies S and Molnár P 2018 Finance Research Letters 27 223–227
  * [12] Phillip A, Chan J and Peiris S 2019 Finance Research Letters 28 95–100
  * [13] Takaishi T 2020 Finance Research Letters 32 101379
  * [14] Takaishi T 2020 EPL (Europhysics Letters) 129 28001
  * [15] Kristjanpoller W, Bouri E and Takaishi T 2020 Physica A 545 123711
  * [16] Cont R 2001 Quantitative Finance 1 223–236
  * [17] Chu J, Chan S, Nadarajah S and Osterrieder J 2017 Journal of Risk and Financial Management 10 17
  * [18] Gopikrishnan P, Meyer M, Amaral L N and Stanley H E 1998 The European Physical Journal B-Condensed Matter and Complex Systems 3 139–140
  * [19] Gopikrishnan P, Plerou V, Amaral L A N, Meyer M and Stanley H E 1999 Physical Review E 60 5305
  * [20] Plerou V, Gopikrishnan P, Amaral L A N, Meyer M and Stanley H E 1999 Physical review E 60 6519
  * [21] Pan R K and Sinha S 2007 EPL (Europhysics Letters) 77 58004
  * [22] Easwaran S, Dixit M and Sinha S 2015 Econophysics and data driven modelling of market dynamics (Springer) pp 121–128
  * [23] Begušić S, Kostanjčar Z, Stanley H E and Podobnik B 2018 Physica A 510 400–406
  * [24] Urquhart A 2016 Economics Letters 148 80–82
  * [25] Takaishi T and Adachi T 2020 Asia-Pacific Financial Markets 27 145–154
  * [26] Engle R F 1982 Econometrica: Journal of the Econometric Society 987–1007
  * [27] Bollerslev T 1986 Journal of Econometrics 31 307–327
  * [28] Glosten L, Jaganathan R and Runkle D 1993 Journal of Finance 48 1779–1801
  * [29] Nelson D 1991 Econometrica 59 347–370
  * [30] Sentana E 1995 Review of Economic Studies 62 639–661
  * [31] Takaishi T 2017 Physica A 473 451–460
  * [32] Takaishi T 2018 Quantitative Finance and Economics 2 127–36
  * [33] Bollerslev T, Chou R Y and Kroner K F 1992 Journal of econometrics 52 5–59
  * [34] Andersen T G and Bollerslev T 1998 International economic review 885–905
  * [35] McAleer M and Medeiros M C 2008 Econometric Reviews 27 10–45
  * [36] Andersen T G, Bollerslev T, Diebold F X and Labys P 2000 Multinational Finance Journal 4 159–179
  * [37] Andersen T G, Bollerslev T, Diebold F X and Labys P 2001 Journal of the American statistical association 96 42–55
  * [38] Andersen T G, Bollerslev T, Diebold F X and Ebens H 2001 Journal of financial economics 61 43–76
  * [39] Andersen T G, Bollerslev T and Dobrev D 2007 Journal of Econometrics 138 125–180
  * [40] Takaishi T 2012 Procedia-Social and Behavioral Sciences 65 968–973
  * [41] Takaishi T 2014 JPS Conf. Proc. 1 019007
  * [42] Takaishi T and Watanabe T 2016 Journal of Physics: Conference Series 710 012010
  * [43] Takaishi T 2018 Physica A 500 139–154
  * [44] Takaishi T, Chen T T and Zheng Z 2012 Progress of Theoretical Physics Supplement 194 43–54
  * [45] Takaishi T 2013 Journal of Physics: Conference Series 454 012041
  * [46] Bandi F M and Russell J R 2006 Journal of Financial Economics 79 655–692
  * [47] Liu L Y, Patton A J and Sheppard K 2015 Journal of Econometrics 187 293–311


# Macroscopic quantum tunneling in spin filter ferromagnetic Josephson
junctions

D. Massarotti<EMAIL_ADDRESS>Universitá degli Studi di Napoli
”Federico II”, Dipartimento di Fisica, via Cinthia, 80126 Napoli (Na) Italy
CNR-SPIN UOS Napoli, Complesso Universitario Monte Sant’Angelo via Cinthia,
80126 Napoli (Na) Italy  A. Pal  Department of Materials Science and
Metallurgy, University of Cambridge, 27 Charles Babbage Road, Cambridge, CB3
0FS, United Kingdom  G. Rotoli  Seconda Universitá degli Studi di Napoli,
Dipartimento di Ingegneria Industriale e dell’Informazione, via Roma 29, 81031
Aversa (Ce) Italy  L. Longobardi  Seconda Universitá degli Studi di Napoli,
Dipartimento di Ingegneria Industriale e dell’Informazione, via Roma 29, 81031
Aversa (Ce) Italy  American Physical Society, 1 Research Road, Ridge, New York
11961, USA  M. G. Blamire  Department of Materials Science and Metallurgy,
University of Cambridge, 27 Charles Babbage Road, Cambridge, CB3 0FS, United
Kingdom  F. Tafuri<EMAIL_ADDRESS>CNR-SPIN UOS Napoli, Complesso
Universitario Monte Sant’Angelo via Cinthia, 80126 Napoli (Na) Italy  Seconda
Universitá degli Studi di Napoli, Dipartimento di Ingegneria Industriale e
dell’Informazione, via Roma 29, 81031 Aversa (Ce) Italy

###### Abstract

The interfacial coupling of two materials with different ordered phases, such
as a superconductor (S) and a ferromagnet (F) is driving new fundamental
physics and innovative applications. For example, the creation of spin-filter
Josephson junctions and the demonstration of triplet supercurrents have
suggested the potential of a dissipationless version of spintronics based on
unconventional superconductivity. Here we demonstrate evidence for active
quantum applications of S-F-S junctions, through the observation of
macroscopic quantum tunneling in Josephson junctions with GdN ferromagnetic
insulator barriers. We prove a clear transition from thermal to quantum regime
at a crossover temperature of about 100 mK at zero magnetic field in junctions
which demonstrate a clear signature of unconventional superconductivity.
Following previous demonstration of passive S-F-S phase shifters in a phase
qubit, our result paves the way to the active use of spin filter Josephson
systems in quantum hybrid circuits.

Introduction

Superconductor (S)/Ferromagnet (F) interfaces are a special class of hybrid
systems where different ordered phases meet and generate new nanoscale
phenomena and new forms of global orderbuzdin_rew ; efetov_rew . A
ferromagnetic barrier in a Josephson junction (JJ) generates novel
physicsryazanov ; robinson , and represents a technological key for advances
in weak superconductivity, spintronics and quantum computationgolubov ; Buzdin
; kawabata . Recent interest in ultra-low-power, high-density cryogenic
memories has spurred new efforts to simultaneously exploit superconducting and
magnetic properties so as to create novel switching elements having these two
competing ordersbaek . S-F-S junctions are expected to shed light on several
aspects of unconventional superconductivity, including transport through spin
aligned triplet Cooper pairsbergeret_prl ; fominov ; robinson ; robinson3 ;
birge ; sprungmann ; aarts .

The NbN-GdN-NbN junctions investigated in these experiments are spin filter
devicesblamire  with an unconventional predominant second harmonic current-
phase relation (CPR)pal1 . By changing the thickness of the GdN ferromagnetic
insulator (FI) barrier it is possible to change its magnetic properties and
hence the spin filter efficiency (SFE)pal1 ; pal2 . In this paper we report
measurements on NbN-GdN-NbN junctions for very low and very high SFE values,
from almost a few percent up to 90$\%$, with an intermediate value of 60$\%$.
By increasing the thickness of GdN we also obtain junctions with lower values
of the IcRn product, where Ic is the critical current and Rn the normal state
resistance. Following the previous demonstration by Feofanov et al.ustinov  of
passive S-F-S phase shifters in a phase qubit, the experiment reported here
demonstrates the active quantum potential of S-F-S JJs via the occurrence of
macroscopic quantum tunneling (MQT) in spin filter devices.

Results

Transport properties of spin filter Josephson junctions.

The hallmark of the spin-filter effect is the decrease in resistance in the R
vs T curves below the ferromagnetic transition, as shown in Fig. 1a
(T${}_{Curie}\simeq$ 33 K for GdN). SFE is defined as the percentage
difference in the tunnelling probability for up/down spin electrons owing to
the difference in barrier heights of the up/down spin channel in the FI caused
by exchange splitting, so that 100% SFE corresponds to pure tunnelling of one
spin sign. See Methods for details concerning the calculation of SFE. Table 1
collects the parameters of the measured junctions. All the NbN-GdN-NbN JJs
present hysteresis larger than 90$\%$ in the current-voltage (IV)
characteristics, in a wide range of Josephson coupling energies E${}_{J}=$
Ic0$\phi$${}_{0}/2$$\pi$ (where Ic0 is the critical current in absence of
thermal fluctuations and $\phi$${}_{0}=$ h/2e is the quantum flux). For all
the measured devices the junction area is about 7$\mu$m $\times$ 7$\mu$m
$\simeq$ 50 $\mu$m2.

In Fig. 1d, IV curves, measured as a function of temperature T for the
junction with the highest SFE, are reported. The dependence of Ic as a
function of H at T = 4 K is reported in Fig. 1c. The blue curve is the first
measurement of the magnetic pattern, after nominal zero-field cooling. Ic(H)
has then been measured both in the downward direction of the magnetic field
sweep and in the upward direction (black and red curves, respectively). The
black and red curves show a distinctive shift of the absolute maximum of Ic
from $-1$ mT to $+1$ mT respectively arising from the hysteretic reversal of
the FI barrierblamire_sust . The period of Ic(H) in non-spin filter junctions
for the same geometry is twice as large (3 mT), pointing to a largely
predominant second harmonic in the CPR in spin filter JJs, as discussed in
Ref.pal1 .

Measurements of switching current distributions.

We have studied the escape rate of the superconducting phase $\varphi$ as a
function of temperature and for different values of the magnetic field,
through standard switching current distribution (SCD) measurementsdevoret1985
; martinis1987 ; bauch_2005 ; luigi2 . SCDs have been performed for all the
samples reported in Table 1. According to the Resistively and Capacitively
Shunted Junction (RCSJ) modelbarone ; likharev , for a JJ with a conventional
CPR the dynamics is equivalent to that of a particle of mass mφ moving in a
washboard potential $U(\varphi)=-E_{J}(\cos\\!\varphi+\gamma\varphi)$. The
particle mass is given by $m_{\varphi}=C(\phi_{0}/2\pi)^{2}$ with C the
capacitance of the junction. The normalized bias current $\gamma$$=$I/Ic0
determines the tilt of the potential.

The second harmonic component in the CPR, $I=I_{1}\sin\varphi+I_{2}\sin
2\varphi$, leads to a modified washboard potential
$U(\varphi)=-E_{1}(\cos\\!\varphi+\frac{g}{2}\cos\\!2\varphi+\frac{I_{c0}}{I_{1}}\gamma\varphi)$,
E${}_{1}=$$\hbar$I1/2e, which may assume the form of a double well for values
of g$=$I${}_{2}/$I1 larger than 0.5tza03 ; goldobin1 ; goldobin2 ; goldobin3
(see Fig. 2a). The presence of two wells in the washboard potential may result
in the observation of two critical currents in the IV characteristics, since
when tilting back the washboard potential the phase particle may be retrapped
in one of the potential wells with finite probabilitygoldobin1 ; the case of
”$\varphi$ JJs” with g $<-1/2$ has been recently studiedgoldobin2 ; goldobin3
. Measurements of two-well distinguished critical currents constitute a very
direct criterion to estimate the g factorgoldobin1 ; goldobin2 . Although spin
filter JJs have a strong second harmonic component in the CPRpal1  we have not
found evidence of two critical currents, and hence the case of $\varphi$
junction with negative values of g is not considered in this work.

The $\sin\\!2\varphi$ term in the CPR on average lowers the barrier height of
the washboard potential without significantly altering the asymptotic
expression of the potential barrier for $\gamma$ close to one. The height of
the potential barrier is given by: $\Delta
U(\gamma,g)=\frac{4}{3}E_{J}a(g)\left(1-\gamma\right)^{3/2}$, where
$a(g)^{2}=\frac{\sqrt{32g^{2}+1}+3}{2\sqrt{32g^{2}+1}}$, and the plasma
frequency is given by:
$\omega_{p}(\gamma,g)=\sqrt[4]{2/a(g)^{2}}\sqrt{2eI_{c0}/\hbar
C}(1-\gamma^{2})^{1/4}$. In addition, the second harmonic component does not
modify the power law of $\Delta U$. As a consequence in the thermal regime the
standard deviation $\sigma$ of the SCDs is expected to scale as T2/3, as in
the case of the standard CPRdevoret1985 ; martinis1987 .

The motion of the particle is subject to damping given by Q-1, where
$Q=\omega_{p}RC$ is the quality factor and R is the resistance of the
junction. When the bias current is ramped from $\gamma$ $=0$ to $\gamma$ $<1$,
the junction is in the zero voltage state in absence of thermal and quantum
fluctuations, corresponding to the particle being localized in a potential
well. At finite temperature the junction may switch into the finite voltage
state for a bias current $\gamma$ $<1$. This corresponds to the particle
escaping from the well either by thermally activated processes or by tunneling
through the barrier potential (see Fig. 2b). In the thermal activation (TA)
regime, the escape rate for weak to moderate damping (Q $>1$) is determined
bykramers  $\Gamma_{t}=a_{t}\frac{\omega_{p}}{2\pi}\exp\left(-\frac{\Delta
U}{k_{B}T}\right)$, where the thermal prefactor is
$a_{t}=4\left[\left(1+Qk_{B}T/1.8\Delta U\right)^{1/2}+1\right]^{-2}$ but83 .
The escape rate will be dominated by MQT at low enough temperaturescal81 . For
Q $>1$ and $\gamma$ close to 1 the escape rate in the quantum regime is:
$\Gamma_{q}=a_{q}\frac{\omega_{p}}{2\pi}\exp\left(-\frac{\Delta
U}{\hbar\omega_{p}}\right)$, where $a_{q}=\left(864\pi\Delta
U/\hbar\omega_{p}\right)^{1/2}$. The crossover temperature between the thermal
and quantum regimes is given bygrabert

$T_{cross}=\frac{\hbar\omega_{p}}{2\pi
k_{B}}\left\\{\left[1+\left(\frac{1}{2Q}\right)^{2}\right]^{1/2}-\frac{1}{2Q}\right\\}\>.$
(1)

The experimental probability density of switching is related to the escape
rate through the following equationfulton :

$P(I)=\frac{\Gamma(I)}{\Delta I/\Delta
t}\exp\left[-\int_{0}^{I}{\frac{\Gamma(I^{\prime})}{\Delta I^{\prime}/\Delta
t}dI^{\prime}}\right]$ (2)

where $\Delta$I/$\Delta$t is the current ramp rate.

The measurements have been performed in a dilution refrigerator with reaches a
base temperature of 20 mK. A full description of the apparatus is given in
detail in the Methods section. The bias current of the junction is ramped at a
constant sweep rate $\Delta$I/$\Delta$t $=2$ mA s-1 and at least 104 switching
events have been recorded using a standard techniqueprbLuigi .

Fig. 3a shows a set of SCDs as a function of temperature for the high SFE JJ
reported in Fig. 1. The thermal behavior of the SCDs is typical of underdamped
JJs and the standard deviation $\sigma$, which is proportional to the width of
the switching histograms, increasing with temperature as expected. Fig. 3b
shows the SCDs measured below 1 K (black circles). The dependence of the
standard deviation $\sigma$ on temperature is reported in Fig. 4a (right
axis), along with the thermal behavior of the mean value of the SCDs, Imean,
below 0.5 K (left axis). When decreasing the temperature, Imean increases
while $\sigma$ decreases and both saturate at a crossover temperature of about
100 mK. Below this crossover the histograms overlap and the escape process is
no longer regulated by thermal fluctuations, indicating the transition to the
MQT regimedevoret1985 ; martinis1987 .

A further confirmation of the observation of MQT comes from measurements of
SCDs in magnetic field. The behavior of $\sigma$(T) at H $=1.1$ mT is shown in
Fig. 4b. At this value of the magnetic field, which reduces Ic to half of the
value measured at zero field (see the blue squares in Fig. 1c), lower values
of $\sigma$ have been measured and Tcross is reduced by a factor $\sqrt{2}$,
down to about 70 mK, in agreement with MQT theorygrabert . In both cases of 0
mT and 1.1 mT, Tcross has been determined by the intersection of the T2/3
curve in the TA regime (dashed green lines in Fig. 4b) and the mean value of
$\sigma$ in the MQT regime (black full lines in Fig. 4b). The measurements in
presence of magnetic field prove that the flattening of $\sigma$ at H $=0$ mT
is a quantum effect and is not due to noise or heating in the measurement
setupdevoret1985 ; martinis1987 ; bauch_2005 ; luigi2 .

Discussion

In the literature there are no measurements of SCDs on junctions with a
dominant second harmonic component in the CPR. Numerical simulations of the
phase dynamics as a function of the damping parameter Q, the g factor and the
temperature T give the conditions for which a double well potential
effectively behaves in the escape process as a single well (for details see
Supplementary Fig. 1, Supplementary Fig. 2 and Supplementary Note 1). Namely,
for values of Q $\simeq 10$ SCDs with a single peak have been obtained for g
$\leq 2$ or in the limiting case of pure second harmonic CPR (I${}_{1}=0$). In
this case the washboard potential changes its periodicity but assumes the form
of a single well potential. Instead, for g $\geq 2$, two critical currents
should be observed. In fact the heights of the two barriers approach each
other when increasing the g factor, and the phase may be retrapped in both the
potential wells with a finite probability, resulting in a bimodal switching
distribution when counting many escape eventsgoldobin1 ; goldobin2 ; goldobin3
(see Supplementary Fig. 3 and Supplementary Discussion). Since this is not
observed, a pure second harmonic is the only possible explanation consistent
with both measurements of magnetic field patternpal1  (see Fig. 1c) and SCDs.

Ic0 can be obtained by fitting the probability density of switching P(I) (red
lines in Fig. 3b) in the thermal and quantum regime, I${}_{c0}=30.41\pm 0.05$
$\mu$A. A quite accurate value of the capacitance C can be obtained from the
crossover temperature Tcross, which depends on Q, C and Ic0, see equation (1).
By inserting the values of Ic0 and Q in the expression for Tcross we get
C$=4.5\pm 0.9$ pF. These values lead to $\omega$${}_{p}\approx 14$ GHz.
Nevertheless, as shown in Fig. 3c, the function a(g) in the equations for
$\Delta U$ and $\omega$p is a slowly varying function for g $>>1$, thus the
junction parameters weakly depend on the g factor for high values of g.

We expect in future additional insights coming from a comparative analysis
with samples with lower SFE. For junction cross sections of about 50 $\mu$m2
the Ic values for such junctions, are too high to be in the conditions to
observe unambiguously the transition from the thermal to the quantum regime as
commonly occurring also in standard S-Insulator-S junctionsdevoret1985 ;
martinis1987 ; webb3 ; han ; martinis2002 ; Berkley . Only further advances in
fabrication able to insert S-FI-S junctions in cavities and qubit
architectures will give more refined feedback on the modes of the dissipative
domains of the junction, and on the triplet component. However we can infer
that the Q values of spin filter junctions are relatively higher than one
would naively expect on the basis of the properties of the low SFE samples,
which are characterized by higher values of Ic. Concerning the Q values of
spin filter JJs, the reduction in Ic is compensated by the increase of Rn, as
reported in Table 1.

In conclusion, we have demonstrated the occurrence of MQT in NbN-GdN-NbN spin
filter JJs. Spin filtering drives the S-F-S junction in the underdamped regime
and in the appropriate window of junction parameters to observe MQT. The
switching current distributions, together with the period of Ic(H) modulation
provides direct evidence for a pure second harmonic current-phase relation in
the junction where MQT was observed. This is clear evidence of unconventional
superconductivity, and it is possible that transport occurs by means of a pair
of spin aligned triplet Cooper pairspal1  which may suppress magnetic sources
of decoherencehouzet ; trifunovic . Demonstration of macroscopic quantum
phenomena in spin filter devices gives promise for their application in
quantum hybrid circuitsustinov , and also possibly as quiet memories.

Methods

Determination of spin filter efficiency.

Spin filter efficiency (SFE) at a particular temperature is calculated from
the R vs T curve. SFE at any temperature below the Curie temperature (TCurie)
of GdN and above the superconducting transition temperature (Tc ) of NbN is
calculated by defining SFE
$=\frac{\sigma\uparrow-\sigma\downarrow}{\sigma\uparrow+\sigma\downarrow}$
where $\sigma$$\uparrow$, $\sigma$$\downarrow$ are current densities of the up
and down spin channels respectively. By invoking WKB approximation, we can
derive SFE $=\tanh\left(\cosh^{-1}\left(R^{*}/R\right)\right)$ where R∗ is the
measured value of resistance, and R is the value of resistance in the absence
of spin filtering. R is estimated by fitting an exponential to the R vs T
curve to temperatures above 50 K, and extrapolating the exponential to
temperatures below the TCurie of GdN. For calculating SFE below Tc of NbN, an
exponential is fitted to the temperature dependence of SFE between TCurie and
Tc and the fitted exponential is extrapolated to lower temperatures to obtain
an estimate of SFE.

Setup for switching current distributions.

The SCDs have been measured by thermally anchoring the samples to the mixing
chamber of He3/He4 Oxford dilution refrigerator. The bias current is ramped at
a constant sweep rate $\Delta$I/$\Delta$t of about 2 mA s-1, the voltage is
measured using a low noise differential amplifier and is fed into a threshold
detector which is set to generate a pulse signal when the junction switches
from the superconducting state to the finite voltage state. This signal is
used to trigger a fast volt meter to record the value of the switching
current. This procedure is repeated at least $10^{4}$ times at each
temperature, which allows to compile a histogram of the switching currents.
Filtering is guaranteed by a room temperature electromagnetic interference
filter stage followed by low pass RC filters with a cut-off frequency of 1.6
MHz anchored at 1.5 K, and by a combination of copper powder and twisted pair
filters thermally anchored at the mixing chamber of the dilution refrigerator.

Acknowledgments

We acknowledge financial support from COST Action MP1201 [NanoSC COST], by
Progetto FIRB HybridNanoDev RBFR1236VV001 and by Regione Campania through POR
Campania FSE 2007/2013, progetto MASTRI CUP B25B09000010007.

Author Contributions

D.M., A.P., L.L., M.G.B. and F.T. conceived the experiments, A.P. and M.G.B.
designed and realized the junctions; D.M. carried out the measurements; D.M.,
G.R. and F.T. worked on the theoretical modelling and data analysis; D.M.,
M.G.B. and F.T. co-wrote the paper. All authors discussed the results and
commented on the manuscript.

Competing financial interests

The authors declare no competing financial interests.

## References

  * (1) Buzdin, A. I. Proximity effects in superconductor-ferromagnet heterostructures. _Rev. Mod. Phys._ 77, 935-976 (2005).
  * (2) Bergeret, F. S., Volkov, A. F. & Efetov, K. B. Odd triplet superconductivity and related phenomena in superconductor-ferromagnet structures. _Rev. Mod. Phys._ 77, 1321-1373 (2005).
  * (3) Ryazanov, V. V., Oboznov, V. A., Rusanov, A. Yu., Veretennikov, A. V., Golubov, A. A. & Aarts, J. Coupling of two superconductors through a ferromagnet: Evidence for a $\pi$ junction. _Phys. Rev. Lett_ 86, 2427-2430 (2001).
  * (4) Robinson, J. W. A., Witt, J. D. S. & Blamire, M. G. Controlled injection of spin-triplet supercurrents into a strong ferromagnet. _Science_ 329, 59-61 (2010).
  * (5) Golubov, A. A., Kupriyanov, M. & Il’Ichev, E. The current-phase relation in Josephson junctions. _Rev. Mod. Phys._ 76, 411-469 (2004).
  * (6) Buzdin, A. I. Peculiar properties of the Josephson junction at the transition from 0 to $\pi$ state. _Phys. Rev. B_ 72, 100501 (2005).
  * (7) Kawabata, S., Kashiwaya, S., Asano, Y., Tanaka, Y. & Golubov, A. A. Macroscopic quantum dynamics of $\pi$ junctions with ferromagnetic insulators. _Phys. Rev. B_ 74, 180502 (2006).
  * (8) Baek, B., Rippard, W. H., Benz, S. P., Russek, S. E. & Dresselhaus, P. D. Hybrid superconducting-magnetic memory device using competing order parameters. _Nat. Comm._ 5, 3888 (2014).
  * (9) Bergeret, F. S., Volkov, A. F. & Efetov, K. B. Long-range proximity effects in superconductor-ferromagnet structures. _Phys. Rev. Lett._ 86, 4096-4099 (2001).
  * (10) Fominov, Y. V., Golubov, A. A. & Kupriyanov, M. Y. Triplet proximity effects in FSF trilayers. _JETP Lett._ 77, 510 (2003).
  * (11) Halász, G. B., Blamire, M. G. & Robinson, J. W. A. Magnetic-coupling dependent spin-triplet supercurrents in helimagnet/ferromagnet Josephson junctions. _Phys. Rev. B_ 84, 024517 (2011).
  * (12) Khaire, T. S., Khasawneh, M. A., Pratt, W. P. Jr. & Birge, N. O. Observation of spin-triplet superconductivity in Co-based Josephson junctions. _Phys. Rev. Lett._ 104, 137002 (2010).
  * (13) Sprungmann, D., Westerholt, K., Zabel, H., Weides, M. & Kohlstedt, H. Evidence for triplet superconductivity in Josephson junctions with barriers of the ferromagnetic Heusler alloy Cu2MnAl. _Phys. Rev. B_ 82, 060505 (2010).
  * (14) Anwar, M. S., Veldhorst, M., Brinkman, A. & Aarts, J. Long range supercurrents in ferromagnetic CrO2 using a multilayer contact structure. _Appl. Phys. Lett._ 100, 052602 (2012).
  * (15) Senapati, K., Blamire, M. G. & Barber, Z. H. Spin-filter Josephson junctions. _Nat. Mater._ 10, 849-852 (2011).
  * (16) Pal, A., Barber, Z. H., Robinson, J. W. A. & Blamire, M. G. Pure second harmonic current-phase relation in spin-filter Josephson junctions. _Nat. Commun._ 5, 3340 (2014).
  * (17) Pal, A., Senapati, K., Barber, Z. H. & Blamire, M. G. Electric field dependent spin polarization in GdN spin filter tunnel junctions. _Adv. Mater._ 25, 5581-5585 (2013).
  * (18) Feofanov, A. K. _et al._ Implementation of superconductor/ferromagnet/ superconductor $\pi$-shifters in superconducting digital and quantum circuits. _Nature Phys._ 6, 593-597 (2010).
  * (19) Blamire, M. G., Smiet, C. B., Banerjee, N. & Robinson, J. W. A. Field modulation of the critical current in magnetic Josephson junctions. _Supercond. Sci. Technol._ 26, 055017 (2013).
  * (20) Devoret, M. H., Martinis, J. M. & Clarke, J. Measurements of macroscopic quantum tunneling out of the zero-voltage state of a current-biased Josephson junction. _Phys. Rev. Lett._ 55, 1908-1911 (1985).
  * (21) Martinis, J. M., Devoret, M. H. & Clarke, J. Experimental tests for the quantum behavior of a macroscopic degree of freedom: the phase difference across a Josephson junction. _Phys. Rev. B_ 35, 4682-4698 (1987).
  * (22) Bauch, T. _et al._ Macroscopic quantum tunneling in d-wave YBaCuO Josephson junctions. _Phys. Rev. Lett._ 94, 087003 (2005).
  * (23) Longobardi, L. _et al._ Direct transition from quantum escape to a phase diffusion regime in YBaCuO biepitaxial Josephson junctions. _Phys. Rev. Lett._ 109, 050601 (2012).
  * (24) Barone, A. & Paternò, G. _Physics and Applications of the Josephson Effect_ (John Wiley and Sons, 1982).
  * (25) Likharev, K. K. Superconducting weak links. _Rev. Mod. Phys._ 51, 101-159 (1979).
  * (26) Tzalenchuk, A. Y., Lindström, T., Charlebois, S. A, Stepantsov, E. A., Ivanov, Z. & Zagoskin, A. M. Mesoscopic Josephson junctions of high-Tc superconductors. _Phys. Rev. B_ 68, 100501 (2003).
  * (27) Goldobin, E., Koelle, D., Kleiner, R. & Buzdin, A. I. Josephson junctions with second harmonic in the current-phase relation: Properties of $\varphi$ junctions. _Phys. Rev. B_ 76, 224523 (2007).
  * (28) Sickinger, H. _et al._ Experimental evidence of a $\varphi$ Josephson junction. _Phys. Rev. Lett._ 109, 107002 (2012).
  * (29) Goldobin, E., Kleiner, R., Koelle, D. & Mints, R. G. Phase retrapping in a pointlike $\varphi$ Josephson junction: the butterfly effect. _Phys. Rev. Lett._ 111, 057004 (2013).
  * (30) Kramers, H. A. Brownian motion in a field of force and the diffusion model of chemical reactions. _Physica_ (Utrecht) 7, 284-304 (1940).
  * (31) Büttiker, M., Harris, E. P. & Landauer, R. Thermal activation in extremely underdamped Josephson-junction circuits. _Phys. Rev. B._ 28, 1268-1275 (1983).
  * (32) Caldeira, A. O. & Leggett, A. J. Influence of dissipation on quantum tunneling in macroscopic systems. _Phys. Rev. Lett._ 46, 211-214 (1981).
  * (33) Grabert, H. & Weiss, U. Crossover from Thermal Hopping to Quantum Tunneling. _Phys. Rev. Lett._ 53, 1787-1790 (1984).
  * (34) Fulton, T. A. & Dunkleberger, L. N. Lifetime of the zero-voltage state in Josephson tunnel junctions. _Phys. Rev. B_ 9, 4760-4768 (1974).
  * (35) Longobardi, L. _et al._ Thermal hopping and retrapping of a Brownian particle in the tilted periodic potential of a NbN/MgO/NbN Josephson junction. _Phys. Rev. B_ 84, 184504 (2011).
  * (36) Washburn, S., Webb, R. A., Voss, R. F. & Farris, S. M. Effects of dissipation and temperature on macroscopic quantum tunneling. _Phys. Rev. Lett._ 54, 2712-2715 (1985).
  * (37) Li, S. X., Yu, Y., Zhang, Y., Qiu, W., Han S. & Wang, Z. Quantitative study of macroscopic quantum tunneling in a dc SQUID: a system with two degrees of freedom. _Phys. Rev. Lett._ 89, 098301 (2002).
  * (38) Martinis, J. M., Nam, S., Aumentado, J. & Urbina, C. Rabi oscillations in a large Josephson-junction qubit. _Phys. Rev. Lett._ 89, 117901 (2002).
  * (39) Berkley, A. J. _et al._ Entangled macroscopic quantum states in two superconducting qubits. _Science_ 300, 1548-1550 (2003).
  * (40) Richard, C., Houzet, M. & Meyer, J. S. Superharmonic long-range triplet current in a diffusive Josephson junction. _Phys. Rev. Lett._ 110, 217004 (2013).
  * (41) Trifunovic, L. Long-range superharmonic Josephson current. _Phys. Rev. Lett._ 107, 047001 (2011).

Table 1:  Device parameters. Spin filter efficiency has been determined at 4.2 K (see Methods for details). Ic0 has been estimated from the fits of the switching current distributions and Rn from the IV curves. Barrier thickness (nm) | Spin filter efficiency |  Ic0 ($\mu$A) |  EJ (meV) |  IcRn (mV) |
---|---|---|---|---|---
1.5 |  $<$30% | 820 | 1700 | 1.0 |
1.7 |  $<$30% | 280 | 570 | 0.9 |
1.8 | 60% | 120 | 250 | 0.7 |
3.0 | 90% | 30 | 60 | 0.1 |
Figure 1: Spin filter Josephson junctions. (a) Thermal behavior of the
resistance is a hallmark of spin filter NbN-GdN-NbN Josephson junctions (JJs).
Semiconducting behavior has been observed at high temperatures, while below
the Curie temperature of GdN ($\simeq$ 33 K) the resistance decreases when
decreasing the temperature. In non spin filter NbN-GdN-NbN JJs a
semiconducting behavior has been observed down to the transition temperature
of NbN. In panel b a picture of the device area is shown. The black arrow is
the scale bar, 10 $\mu$m. (c) Magnetic field pattern of a spin filter JJ in
the range [-5 mT, 5 mT], measured at 4 K. The blue curve is measured after
nominal zero field cooling, the black and red curves are the magnetic pattern
in the downward and upward direction of the magnetic field sweep,
respectively, as indicated by the coloured dashed arrows. Blue squares
indicate the values of the magnetic field, for which we present the
measurements of switching current distributions. (d) Current-voltage
characteristics of the same junction are shown as a function of the
temperature, along with the voltage threshold (gray line) used for the
switching current measurements. Figure 2: Phase dynamics with strong second
harmonic component. Washboard potential in the case of a strong second
harmonic component in the current-phase relation, for different values of the
bias current $\gamma$. Both panels a and b refer to the case g $=2$, for which
two potential wells are present. A and B indicate the lower and the higher
potential well, respectively. The blue rectangle in a is magnified in panel b,
in which thermal activation (black arrow) and quantum tunneling processes (red
arrow) from the well B at $\gamma$ $=0.9$ are shown. $\omega$p ($\gamma$, g)
and $\Delta U$ ($\gamma$, g) are the plasma frequency and the height of the
potential barrier respectively. Figure 3: Switching current distributions. (a)
Measurements of switching current distributions (SCDs) from 4 K down to 0.3 K.
When increasing the temperature, the histograms move to lower values of
switching currents Isw and the standard deviation $\sigma$ increases. Each
histogram collects 104 switching events. (b) Thermal behavior of the SCDs from
0.9 K down to 20 mK. The red lines are fits of the probability density of
switching, according to equation (2) in the case of a dominant second harmonic
component in the current-phase relation (CPR). Below 100 mK the histograms
overlap, indicating the transition to the quantum regime, definitely confirmed
by the their dependence in magnetic field. In panel c the function a(g) is
plotted (black curve) and the red line is the value of the function in the
limiting case of a pure second harmonic (CPR). Figure 4: MQT in spin filter
S-F-S Josephson junctions. (a) Thermal behavior of the mean switching current
Imean (green up triangles, left axis) and of the standard deviation $\sigma$
(violet circles, right axis) in absence of magnetic field. The red dashed line
indicates the crossover temperature between macroscopic quantum tunneling
(MQT) regime (dark yellow background) and thermal activation (pink
background). In the MQT regime both Imean and $\sigma$ saturate. (b)
Comparison between thermal dependences of $\sigma$ at 0 mT (violet circles)
and at 1.1 mT (orange right triangles). In presence of magnetic field, lower
values of $\sigma$ have been measured and Tcross is reduced, according to
equation (1). Tcross has been determined by the intersection of the T2/3 curve
in the thermal activation regime (dashed green lines) and the mean value of
$\sigma$ in the MQT regime (black full lines). In (a) and in (b) the error
bars have been determined by using the theory for propagation of statistical
errors, from the equations defining the mean and the standard deviation of the
switching histograms.