2021 Egypt Hurgada

# Searching for new physics through neutrino non-standard interactions

Yong Du1  1CAS Key Laboratory of Theoretical Physics, Institute of Theoretical
Physics, Chinese Academy of Sciences, Beijing 100190, P. R. China

###### Abstract

Due to the absence of any definite signals of new physics at colliders and
from precision measurements, it has gradually become more and more popular in
the community to utilize the effective field theory (EFT) framework in
searching for new physics in a model-independent manner. In this letter,
working in the EFT framework and focusing on neutrino non-standard
interactions (NSIs), we report our most recent results on these NSIs from
considering terrestrial neutrino oscillation experiments Daya Bay, Double
Chooz, RENO, T2K and NO$\nu$A, and precision measurements of $N_{\rm eff}$
from Planck and CMB-S4.

###### keywords:

Non-standard interactionsEffective field theories$N_{\rm eff}$

††journal: BSM

##  1 Introduction

Though very successful and has been very precisely tested from various
experiments, the Standard Model (SM) of particle physics can only be a low-
energy effective theory as, for example, it cannot explain neutrino masses
that is essential for neutrino oscillations, dark matter that makes up about
25% of the energy budget of our Universe, as well as the baryon asymmetry of
our Universe that is essential for understanding the existence of beings. For
this reason, the SM has to be extended to enlarge its particle content to
account for aforementioned phenomena. Tremendous models have been invented and
intensively studied in the past decades such as SUSY, 2HDM, seesaw models
etc., but unfortunately, all the experiments so far have reported a null
result for any new physics, implying that if existing, the new physics could
be too heavy to be directly produced from current experiments. As a
consequence, EFTs could serve as an ideal framework in searching for new
physics in a model-independent manner.

The SM EFT (SMEFT), obtained by integrating out the heavy newly introduced
particles above the weak scale $\Lambda_{W}$, respects the SM gauge group $\rm
SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$. Besides the SM Lagrangian, the
SMEFT also consists of a tower of higher dimensional operators that each
respects the SM gauge group. The SMEFT can be parameterized as follows

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm SM}+\sum\limits_{j\geq
5}\frac{C_{j}}{\Lambda^{j-4}}\mathcal{O}^{(j)},$ (1)

where $C_{j}$’s are the Wilson coefficients, $\Lambda$ is the characteristic
scale of new physics and $j$ the dimension of operators. In the case of $j=5$,
one has the well-known Weinberg operators Weinberg:1979sa  that are
responsible for neutrino masses. Now at a much lower scale below $\Lambda_{W}$
and above 2 GeV, one can further integrate out the top quark, the Higgs boson
and the SU(2) gauge bosons in the SMEFT, resulting in the Low-energy EFT
(LEFT) which now respects the $\rm SU(3)_{c}\times U(1)_{\rm EM}$ gauge group.

These higher dimensional operators in the LEFT introduce NSIs to the processes
predicted in the SM. Among these operators, neutrino NSIs have gain
significant attention in the past decades due to the observation of neutrino
oscillations. In general, these neutrino NSIs can be classified into two
categories: Charge-Current (CC) ones and Neutral-Current (NC) ones. On one
hand, both CC and NC NSIs could modify neutrino oscillation probabilities and
thus change the fitting of the oscillation parameters. These parameters are
currently very precisely measured from neutrino oscillation experiments such
as Daya Bay An:2012bu ; An:2013zwz ; Adey:2018zwh , Double Chooz
Ardellier:2006mn ; DoubleChooz:2019qbj , RENO Ahn:2012nd ; Bak:2018ydk , T2K
Abe:2019vii ; Abe:2011ks  and NO$\nu$A Ayres:2007tu ; Acero:2019ksn  in the
standard three active neutrino oscillation picture. Given these precision
measurements, these NSIs, or equivalently the new physics, would have the
chance to unveil themselves indirectly through causing some deviations from
the fiducial SM prediction that can be measured from these neutrino
oscillation experiments.

On the other hand, in the early Universe where the active particles are
electrons, photons and neutrinos, NC NSIs between neutrinos and
electrons/photons or NC NSIs among neutrinos of different flavors could also
modify neutrino decoupling, which in turn would result in a different
prediction of $N_{\rm eff}$, viz., the effective number of relativistic
species in the early Universe. However, since $N_{\rm eff}$ is currently very
stringently constrained by Planck Aghanim:2018eyx  and LEP ALEPH:2005ab  and
would be measured at the percent level in the future by CMB-S4
Abazajian:2016yjj , one naturally expects that $N_{\rm eff}$ could also serve
as an indirect probe in investigating new physics in the LEFT framework.

In this letter, based on our recent work Du:2020dwr ; Du:2021idh , we will
discuss the sensitivity on new physics from both CC and NC neutrino NSIs. For
the former, we will mainly focus on terrestrial neutrino oscillation
experiments Daya Bay, Double Chooz, RENO, T2K and NO$\nu$A, and for the
latter, we study their corrections to $N_{\rm eff}$ and take Planck and CMB-S4
into account for obtaining the constraints on new physics. We find that, for
CC NSIs, terrestrial neutrino oscillation experiments are currently
approaching new physics scale at $\sim$20 TeV, while precision measurements of
$N_{\rm eff}$ are now sensitive to new physics at the $\sim$200 GeV scale.

##  2 Parameterization of Neutrino non-standard interactions

###  2.1 CC NSIs

Since neutrino oscillation experiments are performed at a very low energy
scale, the LEFT, which is obtained by integrating out the heavy new particles
in the UV theory, the top quark, the Higgs boson and the SU(2) gauge bosons,
can be used for studying terrestrial neutrino oscillation experiments model-
indenpendently. The Lagrangian relevant for neutrino oscillations from CC NSIs
can be parameterized as Falkowski:2019kfn

$\displaystyle\mathcal{L}_{\rm CC}\supset$
$\displaystyle-\frac{2V_{ud}}{v^{2}}\left\\{\left[\mathbf{1}+\epsilon_{L}\right]^{ij}_{\alpha\beta}\left(\bar{u}_{i}\gamma^{\mu}P_{L}d_{j}\right)\left(\bar{\ell}_{\alpha}\gamma_{\mu}P_{L}\nu_{\beta}\right)+\left[\epsilon_{R}\right]^{ij}_{\alpha\beta}\left(\bar{u}_{i}\gamma^{\mu}P_{R}d_{j}\right)\left(\bar{\ell}_{\alpha}\gamma_{\mu}P_{L}\nu_{\beta}\right)\right.$
$\displaystyle+\frac{1}{2}\left[\epsilon_{S}\right]^{ij}_{\alpha\beta}(\bar{u}_{i}d_{j})\left(\bar{\ell}_{\alpha}P_{L}\nu_{\beta}\right)-\frac{1}{2}\left[\epsilon_{P}\right]^{ij}_{\alpha\beta}\left(\bar{u}_{i}\gamma_{5}d_{j}\right)\left(\bar{\ell}_{\alpha}P_{L}\nu_{\beta}\right)$
$\displaystyle+\left.\frac{1}{4}\left[\epsilon_{T}\right]^{ij}_{\alpha\beta}\left(\bar{u}_{i}\sigma^{\mu\nu}P_{L}d_{j}\right)\left(\bar{\ell}_{\alpha}\sigma_{\mu\nu}P_{L}\nu_{\beta}\right)+\mathrm{h.c.}\right\\},$
(2)

where $V_{ud}$ is the “ud” component of the CKM matrix, $i,j=\\{1,2,3\\}$ the
flavor of quarks, $\alpha,\beta=\\{e,\mu,\tau\\}$ that of charged leptons and
neutrinos and $P=P_{L,R}$ are the chiral projection operators. The neutrino
NSIs are parameterized by $\epsilon_{L,R,S,P,T}$, and to be generic, we also
include interactions besides the $V-A$ type ones in the above Lagrangian.

Upper bounds of the NSI parameters associated with the neutrino production and
detection. All constraints are given at 95% CL.  NSI parameter Upper bound
Experiments $\left|\epsilon^{s}_{\mu e}\right|$ 0.004 T2K Abe:2020vdv ;
Abe:2019vii ; Zarnecki:2020yag , NO$\nu$A Acero:2019ksn
$\left|\epsilon^{s}_{\mu\mu}\right|$ 0.021
$\left|\epsilon^{s}_{\mu\tau}\right|$ 0.080 $\left|\epsilon^{d}_{ee}\right|$
0.007 $\left|\epsilon^{d}_{\mu e}\right|$ 0.018 $\left|\epsilon^{d}_{\tau
e}\right|$ 0.021 Daya Bay An:2012bu ; Adey:2018zwh , Double Chooz
Ardellier:2006mn ; DoubleChooz:2019qbj , $\left|\epsilon^{s}_{ee}\right|$
0.007 and RENO Ahn:2012nd ; Bak:2018ydk  $\left|\epsilon^{s}_{e\mu}\right|$
0.018 $\left|\epsilon^{s}_{e\tau}\right|$ 0.021

To connect these CC NSIs to the neutrino oscillation probabilities that are
formulated in the quantum mechanical framework in terms of the production
parameter $\epsilon^{s}$ and the detection parameter $\epsilon^{d}$, one can
follow the procedure developed in Ref. Falkowski:2019kfn  for a consistent
matching between $\epsilon^{s,d}$ in the quantum mechanical formalism and the
NSI parameters $\epsilon_{L,R,S,P,T}$ in the LEFT formalism. For reactor-type
neutrino oscillation experiments like Daya Bay, Double Chooz and RENO and
long-baseline neutrino oscillation experiments like T2K and NO$\nu$A, since
the NSI parameters are very stringently constrained Du:2020dwr  as summarized
in Table 2.1, the linear order matching between these two formalisms would be
enough for our purpose. To be specific, for the connection between
$\epsilon^{s,d}$ and the CC NSI parameters $\epsilon_{L,R,S,P,T}$ for beta,
inverse beta and pion decay, we cite the results in Ref. Falkowski:2019kfn
below and refer the readers to the Appendix of Ref. Du:2020dwr  for a detailed
example for the calculation:

$\displaystyle\epsilon_{e\beta}^{s}=$
$\displaystyle\left(\epsilon_{L}-\epsilon_{R}-\frac{g_{T}}{g_{A}}\frac{m_{e}}{f_{T}\left(E_{\nu}\right)}\epsilon_{T}\right)_{e\beta}^{*},\quad{\text{~{}($\beta$
decay)}}$ (3) $\displaystyle\epsilon_{\beta e}^{d}=$
$\displaystyle\left(\epsilon_{L}+\frac{1-3g_{A}^{2}}{1+3g_{A}^{2}}\epsilon_{R}-\frac{m_{e}}{E_{\nu}-\Delta}\left(\frac{g_{S}}{1+3g_{A}^{2}}\epsilon_{S}-\frac{3g_{A}g_{T}}{1+3g_{A}^{2}}\epsilon_{T}\right)\right)_{e\beta},{\text{~{}(inverse
$\beta$ decay)}}$ (4) $\displaystyle\epsilon_{\mu\beta}^{s}=$
$\displaystyle\left(\epsilon_{L}-\epsilon_{R}-\frac{m_{\pi}^{2}}{m_{\mu}\left(m_{u}+m_{d}\right)}\epsilon_{P}\right)_{\mu\beta}^{*},\quad{\text{~{}(pion
decay)}}$ (5)

with $g_{S}$, $g_{A}$ and $g_{T}$ the scalar, axial-vector and tensor charges
of the nucleon, $\Delta\equiv m_{n}-m_{p}$ the difference between neutron and
proton masses, $E_{\nu}$ the neutrino energy, and $f_{T}(E_{\nu})$ the nucleon
form factor resulting from the tensor-type NSI in eq. (2).

With the matching formulae in eqs. (3-5), one can readily transfer the upper
bounds on $\epsilon^{s,d}$ from neutrino oscillation experiments summarized in
Table 2.1 onto these of $\epsilon_{L,R,S,P,T}$. Recall that
$\epsilon_{L,R,S,P,T}$ are the LEFT Wilson coefficients obtained by
integrating out the top quark, the Higgs boson and the SU(2) gauge bosons from
the SMEFT, therefore, ultimately, $\epsilon_{L,R,S,P,T}$ are functions of the
UV scale $\Lambda$ and the Wilson coefficients of the SMEFT in eq. (1). Thus,
from the upper bounds on $\epsilon_{L,R,S,P,T}$, one can obtain the lower
bounds on the new physics scale $\Lambda$ once the Wilson coefficients are
fixed. In addition, due to the large energy gap between the LEFT and the
SMEFT, one needs to take the running effects into account to correctly gain
any information on $\Lambda$. The running effects can be systematically
included by solving the renormalization group equations. We show in Figure 1
the details of this procedure described above.

Figure 1: A schematic description of how to relate the UV scale $\Lambda$ to
the neutrino NSI parameters. The vertical axis defines the EFT scales. Figure
adopted from Ref. Du:2020dwr .

###  2.2 NC NSIs

The effects of NC NSIs become significant for long-baseline neutrino
oscillation experiments where matter effects start to impact. In literatures,
the dimension-6 LEFT NC NSI operators are usually parameterized as follows:

$\displaystyle\mathcal{L}_{\rm NSI}^{\rm
NC}=-2\sqrt{2}G_{F}\sum_{\alpha,\beta,f,P}\epsilon_{\alpha\beta}^{f,P}\left(\bar{\nu}_{\alpha}\gamma_{\mu}P_{L}\nu_{\beta}\right)\left(\bar{f}\gamma^{\mu}Pf\right),$
(6)

with $f=e,u,d$ the charged fermions, $\alpha,\beta=e,\mu,\tau$ the flavor of
neutrinos and $G_{F}$ the Fermi constant.

For a detailed study on the effects on neutrino oscillations from NC NSIs, see
Ref. Du:2021rdg . In this letter, as discussed in the introduction, we will
study these NC NSIs from a different aspect by studying their effects on
neutrino decoupling in the early Universe. To that end, we consider all
possible LEFT NC NSI operators up to dimension-7. These operators are
enumerated in Table 2.2.

LEFT operators relevant for $N_{\rm eff}$ up to dimension-7 with
$\alpha,\beta,\alpha^{\prime},\beta^{\prime}=e,\mu,\tau$, the neutrino flavor
indices, and $f=e$. The symbol “$c$” along with related operators means charge
conjugation. Our convention for the Wilson coefficients are shown in the last
column.  Dimensions Operators Wilson coefficients dimension-5
$\mathcal{O}_{1}^{(5)}=\frac{e}{8\pi^{2}}\left(\bar{\nu}_{\beta}\sigma^{\mu\nu}P_{L}\nu_{\alpha}\right)F_{\mu\nu}$
$C_{1}^{(5)}$ dimension-6
$\mathcal{O}_{1,f}^{(6)}=\left(\bar{\nu}_{\beta}\gamma_{\mu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma^{\mu}f\right)$
$C_{1,f}^{(6)}$
$\mathcal{O}_{2,f}^{(6)}=\left(\bar{\nu}_{\beta}\gamma_{\mu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma^{\mu}\gamma_{5}f\right)$
$C_{2,f}^{(6)}$
$\mathcal{O}_{3}^{(6)}=\left(\overline{\nu^{c}}_{\beta}P_{L}\nu_{\alpha}\right)\left(\overline{\nu^{c}}_{\beta^{\prime}}P_{L}\nu_{\alpha^{\prime}}\right)$
$C_{3}^{(6)}$
$\mathcal{O}_{4}^{(6)}=\left(\bar{\nu}_{\beta}\gamma_{\mu}P_{L}\nu_{\alpha}\right)\left(\bar{\nu}_{\beta^{\prime}}\gamma_{\mu}P_{L}\nu_{\alpha^{\prime}}\right)$
$C_{4}^{(6)}$
$\mathcal{O}_{5}^{(6)}=\left(\overline{\nu^{c}}_{\beta}\sigma^{\mu\nu}P_{L}\nu_{\alpha}\right)\left(\overline{\nu^{c}}_{\beta^{\prime}}\sigma^{\mu\nu}P_{L}\nu_{\alpha^{\prime}}\right)$
$C_{5}^{(6)}$ dimension-7
$\mathcal{O}_{1}^{(7)}=\frac{\alpha}{12\pi}\left(\overline{\nu^{c}}_{\beta}P_{L}\nu_{\alpha}\right)F^{\mu\nu}F_{\mu\nu}$
$C_{1}^{(7)}$
$\mathcal{O}_{2}^{(7)}=\frac{\alpha}{8\pi}\left(\overline{\nu^{c}}_{\beta}P_{L}\nu_{\alpha}\right)F^{\mu\nu}\widetilde{F}_{\mu\nu}$
$C_{2}^{(7)}$
$\mathcal{O}_{5,f}^{(7)}=m_{f}\left(\overline{\nu^{c}}_{\beta}P_{L}\nu_{\alpha}\right)(\bar{f}f)$
$C_{5,f}^{(7)}$
$\mathcal{O}_{6,f}^{(7)}=m_{f}\left(\overline{\nu^{c}}_{\beta}P_{L}\nu_{\alpha}\right)\left(\bar{f}i\gamma_{5}f\right)$
$C_{6,f}^{(7)}$
$\mathcal{O}_{7,f}^{(7)}=m_{f}\left(\overline{\nu^{c}}_{\beta}\sigma^{\mu\nu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\sigma_{\mu\nu}f\right)$
$C_{7,f}^{(7)}$
$\mathcal{O}_{8,f}^{(7)}=\left(\overline{\nu^{c}}_{\beta}i\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}_{\mu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma^{\mu}f\right)$
$C_{8,f}^{(7)}$
$\mathcal{O}_{9,f}^{(7)}=\left(\overline{\nu^{c}}_{\beta}i\stackrel{{\scriptstyle\leftrightarrow}}{{\partial}}_{\mu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma^{\mu}\gamma_{5}f\right)$
$C_{9,f}^{(7)}$
$\mathcal{O}_{10,f}^{(7)}=\partial_{\mu}\left(\overline{\nu^{c}}_{\beta}\sigma^{\mu\nu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma_{\nu}f\right)$
$C_{10,f}^{(7)}$
$\mathcal{O}_{11,f}^{(7)}=\partial_{\mu}\left(\overline{\nu^{c}}_{\beta}\sigma^{\mu\nu}P_{L}\nu_{\alpha}\right)\left(\bar{f}\gamma_{\nu}\gamma_{5}f\right)$
$C_{11,f}^{(7)}$

Each operator in Table 2.2 would modify the scattering and/or annihilation
rates of neutrinos in the early Universe, thus also affect the energy
injection rate (through neutrino-electromagnetic scattering) into the neutrino
sector from the electromagnetic plasma or energy redistribution (through
neutrino self-interactions) in the neutrino sector in the early Universe. As a
consequence, the predicted effective number of relativistic species, $N_{\rm
eff}$Shvartsman:1969mm ; Steigman:1977kc ; Mangano:2001iu  defined below, will
also change:

$\displaystyle\rho_{R}=\left[1+\frac{7}{8}\left(\frac{4}{11}\right)^{\frac{4}{3}}N_{\rm
eff}\right]\rho_{\gamma}$ (7)

where $\rho_{\gamma}$ is the photon energy density, and $\rho_{R}$ the total
energy density from all relativistic species during the epoch of neutrino
decoupling.

To obtain the exact prediction of $N_{\rm eff}$ with the inclusion of these
operators in Table 2.2, one shall solve $\rho_{R}$ and $\rho_{\gamma}$ from
the coupled Boltzmann equations, for which we follow Refs. Escudero:2018mvt ;
Escudero:2020dfa . It turns out that the most challenging part in solving
these Boltzmann equations comes from the collision term integrals that are
extremely expensive to evaluate both numerically and analytically. For this
reason and for the benefit of the community in the future, a complete, generic
and analytic dictionary is provided in Ref. Du:2021idh , which makes the
Boltzmann equations simple ordinary differential equations that are
numerically trivially solvable.

##  3 Results

###  3.1 Constraints on dimension-6 SMEFT operators from CC NSIs

Following the procedure described in section 2.1 and depicted in Figure 1, we
study constrains on dimension-6 SMEFT operators from the CC NSIs by
investigating their impacts on terrestrial neutrino oscillation experiments.
We show our results on the UV scale $\Lambda$ in Figure 2 by fixing the Wilson
coefficients at unity and considering one non-vanishing SMEFT operator at a
time. The upper panel is for long-baseline type neutrino oscillation
experiments T2K and NO$\nu$A, and the lower one for reactor-type ones Daya
Bay, Double Chooz and RENO.

From the upper panel, it is clear that currently the long-baseline type
neutrino oscillation experiments are already exploring new physics around the
20 TeV scale as indicated by the $\mathcal{O}_{lequ_{1211}}^{(1)}$ and
$\mathcal{O}_{lequ_{2211}}^{(1)}$ operators. In contrast, as seen from the
lower panel of Figure 2, the reactor-type ones are barely approaching the 5
TeV scale. However, it is worth pointing out that these two types of
experiments are sensitive to different subsets of the SMEFT operators. For
illustration, while the reactor-type experiments constrain the
$\mathcal{O}_{lq_{1111}}^{(3)}$ operator to be above about 5 TeV, the long-
baseline neutrino oscillation experiments are insensitive to this operator. We
thus conclude that these two types of neutrino oscillation experiments are
complementary to each other in searching for new physics.

---
Figure 2:  Upper panel: Constraints on dimension-6 SMEFT operators in long-
baseline neutrino oscillation experiments T2K and NO$\nu$A at 95% CL. Lower
panel: Constraints on dimension-6 SMEFT operators in reactor-type neutrino
oscillation experiments Daya Bay, Double Chooz and RENO at 95% CL. All
constraints are obtained by fixing the Wilson coefficients at one. Plots
adopted from Ref. Du:2020dwr .

###  3.2 Constraints on NC NSIs from $N_{\rm eff}$

To obtain constraints on the NC NSIs, we first define the change in $N_{\rm
eff}$ as

$\displaystyle\Delta N_{\rm eff}=N_{\rm eff}^{\rm SM+EFT}-N_{\rm eff}^{\rm
SM},$ (8)

where $N_{\rm eff}^{\rm SM+EFT}$ is the theoretical prediction of $N_{\rm
eff}$ from the LEFT with the NC NSI operators in Table 2.2, and $N_{\rm
eff}^{\rm SM}=3.044$ Akita:2020szl ; Froustey:2020mcq  is that from the SM.
For Planck, we use the current result $N_{\rm eff}=2.99^{+0.34}_{-0.33}$
Aghanim:2018eyx  at the 95% CL, and $\Delta N_{\rm eff}<0.06$ at 95% CL for
CMB-S4 Abazajian:2016yjj ; Abazajian:2019tiv ; Abitbol:2017nao ;
Abazajian:2019eic . As discussed in section 2.2, after solving the Boltzmann
equations with the help of the complete dictionary presented in Ref.
Du:2021idh  one can readily obtain the corrections to $N_{\rm eff}$ from these
LEFT NC NSI operators listed in Table 2.2. We show our results in Figure 3 by
requiring the change in $N_{\rm eff}$ to be within the uncertainties from
current Planck data and the planned precision goal of CMB-S4. The results are
shown in orange for Planck and purple for CMB-S4 in Figure 3.

From Figure 3, one concludes that the dimension-6 $\nu-e$ contact interacting
operators $\mathcal{O}_{1,e}^{(6)}$ and $\mathcal{O}_{2,e}^{(6)}$ are
relatively more constrained to be above $\sim$200 GeV and $\sim$86 GeV
respectively. In addition, the constraint on the vector-type operator
$\mathcal{O}_{1,e}^{(6)}$ is stronger than that on the axial-vector operator
$\mathcal{O}_{2,e}^{(6)}$ since the former interfere with the SM
constructively while the latter destructively: In the constructive case,
neutrinos and electrons are effectively more tightly coupled such that the
neutrinos would decouple earlier from the electromagnetic plasma than that in
the destructive case, thus the neutrinos would be hotter at the time of
decoupling and therefore a larger prediction of $N_{\rm eff}^{\rm SM+EFT}$.
This explains why the lower bound on $\Lambda$ from the
$\mathcal{O}_{1,e}^{(6)}$ operator is larger compared to that from the
$\mathcal{O}_{2,e}^{(6)}$ operator.

The dimension-7 operators in Figure 3 are only constrained to be below
$\sim$10 GeV from the current results of Planck as well as the planned
precision goal of CMB-S4. The reason is that these operators are suppressed by
one more power of $\Lambda$ compared with the dimension-6 ones. Nonetheless,
these constraints on the dimension-7 operators shown in Figure 3 are firstly
obtained in Ref. Du:2021idh  and would serve as the starting point in the
future when much higher precision measurements become avaliable.

Last but not least, the constraints on $\mathcal{O}_{8,e}^{(7)}$ and
$\mathcal{O}_{10,e}^{(7)}$ turn out to be the same as can also be seen from
Figure 3. This can be understood from their analytical expression for these
two operators:

$\displaystyle\rho_{\nu}^{\rm tot.}(\mathcal{O}_{8,e}^{(7)})-\rho_{\nu}^{\rm
tot.}(\mathcal{O}_{10,e}^{(7)})=\frac{960}{\pi^{5}\Lambda^{6}}\left(C_{8,e}^{(7)}-C_{10,e}^{(7)}\right)\left(C_{8,e}^{(7)}+C_{10,e}^{(7)}\right)\times\mathcal{F}(T_{\gamma},T_{\nu_{e}},T_{\nu_{\mu}},\mu_{\nu_{e}},\mu_{\nu_{\mu}}),$
(9)

where we omit the irrelevant tiny neutrino masses here and
$\mathcal{F}(T_{\gamma},T_{\nu_{e}},T_{\nu_{\mu}},\mu_{\nu_{e}},\mu_{\nu_{\mu}})$
is a function that only depends on photon temperature, neutrino temperatures
and neutrino chemical potentials.111The explicit expression of $\mathcal{F}$
is lengthy and does not matter here. For its full expression, see the
supplementary material of Ref. Du:2021idh . It can be shown that these two
operators are related by the equation of motion for the neutrinos as shown in
Ref. Li:2021phq , we thank Xiao-Dong Ma for pointing this out to us. Clearly,
these two operators contributes exactly the same to the total energy density
of neutrinos when their Wilson coefficients are fixed to be the same. We point
out that a similar observation also holds for the $\mathcal{O}_{9,e}^{(7)}$
and the $\mathcal{O}_{11,e}^{(7)}$ operators.

Figure 3: Constraints on the UV scale $\Lambda$ from $\Delta N_{\rm
eff}=N_{\rm eff}^{\rm SM+EFT}-N_{\rm eff}^{\rm SM}$, where $N_{\rm eff}^{\rm
SM}$ is the SM prediction of $N_{\rm eff}$, and $N_{\rm eff}^{\rm SM+EFT}$ is
that from SM and the LEFT operators in Table 2.2. All histograms here are
obtained by fixing the Wilson coefficients at unity and considering only one
non-vanishing operator only at a time. This plot is adopted from Ref.
Du:2021idh .

##  4 Conclusions and outlook

We investigate both CC and NC neutrino NSIs in this letter in the EFT
framework. For CC NSIs, we study the effects of dimension-6 SMEFT operators on
terrestrial neutrino oscillation experiments. By fixing the Wilson
coefficients at unity and assuming only one operator dominate at a time, we
find that long-baseline neutrino oscillation experiments T2K and NO$\nu$A are
currently approaching new physics at the 20 TeV scale while reactor neutrino
oscillation experiments Daya Bay, Double Chooz and RENO are sensitive to new
physics at $\sim$5 TeV scale. Even though these two types of neutrino
oscillation experiments are sensitive to new physics at much different scales,
we stress that they are complementary to each other in the sense that they are
sensitive to different subsets of SMEFT operators as indicated in Figure 2.

For NC NSIs, we study in detail their rule in neutrino decoupling in the early
Universe for LEFT operators up to dimension-7. We find that the dimension-6
operators are more constrained compared with the dimension-7 ones as they are
less suppressed by the UV scale $\Lambda$. Quantitatively, from the
$\mathcal{O}_{1,e}^{(6)}$ operator, $\Lambda$ is constrained to be above about
200 GeV while that from $\mathcal{O}_{2,e}^{(6)}$ is about 86 GeV. The
difference in these two operators comes from the fact the former interfere
with the SM constructively while the latter destructively. Furthermore,
constraints on dimension-7 operators as shown in Figure 3 are one order of
magnitude smaller than those on the dimension-6 operators, but these results
are firstly obtained from our work in Ref. Du:2021idh  and would serve as the
starting point in the future when higher precision measurements become
possible.

We comment on that, in the future when long-baseline experiments like JUNO and
DUNE start operation, the matter effects from these NC NSIs would start to
impact. In this case, to obtain the correct bounds on the UV physics in the
EFT framework, one has to include the matter effects from these NC NSIs. For a
detailed study on this topic, see Ref. Du:2021rdg . We also point out that the
correlation among multiple operators may change the constraints on the UV
scale $\Lambda$ by orders of magnitude Du:2020dwr . In both the CC and NC NSI
cases, it is worth obtaining a global fitting on multiple operators to fully
understand the correlation among multiple operators and their impact on
constraining the UV physics.

## Acknowledgements

The author is indebted to Hao-Lin Li, Jian Tang, Sampsa Vihonen and Jiang-Hao
Yu for the collaboration. The author also wishes to thank the committee of the
“Beyond Standard Model: From Theory to Experiment (BSM-2021)” workshop. This
work is supported by the National Natural Science Foundation of China (NSFC)
under Grants No. 12022514 and No. 11875003 and CAS Project for Young
Scientists in Basic Research YSBR-006, by the National Natural Science
Foundation of China (NSFC) under Grant No. 12047503, by the National Key
Research and Development Program of China under Grant No. 2020YFC2201501, and
also by the internal funds for postdocs from the CAS Key Laboratory of
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences.

## References

  * [1] Steven Weinberg.  Baryon and Lepton Nonconserving Processes.  Phys. Rev. Lett., 43:1566–1570, 1979.
  * [2] F.P. An et al.  Improved Measurement of Electron Antineutrino Disappearance at Daya Bay.  Chin. Phys. C, 37:011001, 2013.
  * [3] F.P. An et al.  Spectral measurement of electron antineutrino oscillation amplitude and frequency at Daya Bay.  Phys. Rev. Lett., 112:061801, 2014.
  * [4] D. Adey et al.  Measurement of the Electron Antineutrino Oscillation with 1958 Days of Operation at Daya Bay.  Phys. Rev. Lett., 121(24):241805, 2018.
  * [5] F. Ardellier et al.  Double Chooz: A Search for the neutrino mixing angle theta(13).  6 2006.
  * [6] H. de Kerret et al.  Double Chooz $\theta_{13}$ measurement via total neutron capture detection.  Nature Phys., 16(5):558–564, 2020.
  * [7] J.K. Ahn et al.  Observation of Reactor Electron Antineutrino Disappearance in the RENO Experiment.  Phys. Rev. Lett., 108:191802, 2012.
  * [8] G. Bak et al.  Measurement of Reactor Antineutrino Oscillation Amplitude and Frequency at RENO.  Phys. Rev. Lett., 121(20):201801, 2018.
  * [9] K. Abe et al.  Constraint on the matter–antimatter symmetry-violating phase in neutrino oscillations.  Nature, 580(7803):339–344, 2020.
  * [10] K. Abe et al.  The T2K Experiment.  Nucl. Instrum. Meth. A, 659:106–135, 2011.
  * [11] D.S. Ayres et al.  The NOvA Technical Design Report.  10 2007.
  * [12] M.A. Acero et al.  First Measurement of Neutrino Oscillation Parameters using Neutrinos and Antineutrinos by NOvA.  Phys. Rev. Lett., 123(15):151803, 2019.
  * [13] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys., 641:A6, 2020.
  * [14] S. Schael et al.  Precision electroweak measurements on the $Z$ resonance.  Phys. Rept., 427:257–454, 2006.
  * [15] Kevork N. Abazajian et al.  CMB-S4 Science Book, First Edition.  10 2016.
  * [16] Yong Du, Hao-Lin Li, Jian Tang, Sampsa Vihonen, and Jiang-Hao Yu.  Non-standard interactions in SMEFT confronted with terrestrial neutrino experiments.  JHEP, 03:019, 2021.
  * [17] Yong Du and Jiang-Hao Yu.  Neutrino non-standard interactions meet precision measurements of Neff.  JHEP, 05:058, 2021.
  * [18] Adam Falkowski, Martín González-Alonso, and Zahra Tabrizi.  Consistent QFT description of non-standard neutrino interactions.  JHEP, 11:048, 2020.
  * [19] K. Abe et al.  T2K measurements of muon neutrino and antineutrino disappearance using $3.13\times 10^{21}$ protons on target.  8 2020.
  * [20] G. Zarnecki.  Most recent T2K results on CPV in neutrino sector.  Nuovo Cim. C, 42(6):246, 2020.
  * [21] Yong Du, Hao-Lin Li, Jian Tang, Sampsa Vihonen, and Jiang-Hao Yu.  Exploring SMEFT Induced Non-Standard Interactions from COHERENT to Neutrino Oscillations.  6 2021.
  * [22] V.F. Shvartsman.  Density of relict particles with zero rest mass in the universe.  Pisma Zh. Eksp. Teor. Fiz., 9:315–317, 1969.
  * [23] G. Steigman, D.N. Schramm, and J.E. Gunn.  Cosmological Limits to the Number of Massive Leptons.  Phys. Lett. B, 66:202–204, 1977.
  * [24] G. Mangano, G. Miele, S. Pastor, and M. Peloso.  A Precision calculation of the effective number of cosmological neutrinos.  Phys. Lett. B, 534:8–16, 2002.
  * [25] Miguel Escudero.  Neutrino decoupling beyond the Standard Model: CMB constraints on the Dark Matter mass with a fast and precise $N_{\rm eff}$ evaluation.  JCAP, 02:007, 2019.
  * [26] Miguel Escudero Abenza.  Precision early universe thermodynamics made simple: $N_{\rm eff}$ and neutrino decoupling in the Standard Model and beyond.  JCAP, 05:048, 2020.
  * [27] Kensuke Akita and Masahide Yamaguchi.  A precision calculation of relic neutrino decoupling.  JCAP, 08:012, 2020.
  * [28] Julien Froustey, Cyril Pitrou, and Maria Cristina Volpe.  Neutrino decoupling including flavour oscillations and primordial nucleosynthesis.  JCAP, 12:015, 2020.
  * [29] Kevork Abazajian et al.  CMB-S4 Decadal Survey APC White Paper.  Bull. Am. Astron. Soc., 51(7):209, 2019.
  * [30] Maximilian H. Abitbol et al.  CMB-S4 Technology Book, First Edition.  6 2017.
  * [31] Kevork Abazajian et al.  CMB-S4 Science Case, Reference Design, and Project Plan.  7 2019.
  * [32] Tong Li, Xiao-Dong Ma, Michael A. Schmidt, and Rui-Jia Zhang.  The implication of $J/\psi\to\gamma+{\rm invisible}$ for the effective field theories of neutrino and dark matter.  4 2021.


Design Rules for DMI-Stabilised Skyrmions

L. Ranno1*, M. Alfonso Moro1

1 Institut Néel, CNRS/Université Grenoble Alpes, Grenoble, France

*<EMAIL_ADDRESS>
## Abstract

Magnetic skyrmions are topological objects, which have recently been observed
in thin films at room temperature. Sub 100-nm sizes and spin polarised current
manipulation make them candidates for high density information storage and
processing. Besides material thickness, skyrmion stability involves many
contributions: ferromagnetic exchange, magnetic anisotropy, Dzyaloshinskii-
Moryia exchange interaction (DMI) and demagnetising energy. Multidimensional
phase diagrams have been reported in the litterature. In this work we propose
a simple two-parameter model, exact for thin films in the negligible
magnetisation limit, e.g. ferrimagnets close to compensation and synthetic
artificial antiferromagnets. It allows to define design rules to stabilize
small diameter skyrmions at room temperature with explicit estimates of their
nucleation and collapse energies. Comparisons to micromagnetic calculations
allow to assess the model quality for non-zero magnetisation and show good
agreement.

###### Contents

  1. 1 Introduction
  2. 2 Model
    1. 2.1 Equations
    2. 2.2 Numerical Solution
  3. 3 Results
  4. 4 Discussion
    1. 4.1 Theoretical DMI to anisotropy energy ratio
    2. 4.2 Local Anisotropy Approximation
    3. 4.3 Neglecting volume charges
    4. 4.4 Design Rules
  5. 5 Conclusion

##  1 Introduction

Following the $\mu$m-size magnetic bubble studies from the 1970s [1, 2, 3],
magnetic skyrmions with a similar topology but smaller sizes have been
experimentally evidenced, at first at low temperature and in bulk samples [4,
5] and then in metallic ultrathin films at room temperature [6, 7]. Magnetic
bubbles are stabilised by their demagnetising energy whereas the skyrmion
stabilisation energy is based on an asymmetric exchange interaction such as
the Dzyaloshinskii Moriya Interaction (DMI) which favors whirling spin
textures. This interaction also selects a whirling chirality and makes
skyrmions homochiral objects [8]. Due to their small size and the metallic
character of the materials, skyrmions are being widely studied as candidates
for high density information storage associated with current driven
manipulation [9].

The main characteristics of $\mu$m-size bubbles (radius, energy) can be
evaluated using a few reduced quantities. For example, their characteristic
lengthscale is the material length $L=\frac{4\sqrt{AK}}{\mu_{0}M_{s}^{2}}$,
where $L$ is the ratio between the Bloch domain wall surface energy and the
magnetostatic characteristic energy. The most stable bubbles [10] have then a
diameter 8$L$, in a film with preferred thickness t=4$L$ i.e., an aspect ratio
2. Diameters in the sub-micron range were achieved [1] but further downscaling
was not expected and semiconducting memories rapidly outperformed bubble
memories.

Several skyrmion-bubble phase diagrammes have recently been reported[11, 12,
13]. However some restrictions are usually present, for example modelling only
Bloch skyrmions[11] or Néel skyrmions with imposed profiles such as double
walls[12] or zero-thickness walls [13], non physical zero-radius limit [12].
When long-range demagnetising effects are present, these diagrammes evidence
the coexistence of magnetostatically-stabilised skyrmionic ”large” bubbles and
DMI-stabilised ”small” skyrmions, with a possible continuous path to transform
one into the other [13]. These results lack simplicity due to the multi
dimensional phase space and require heavy demagnetising field computing. Here
we propose a simple model with good physical insight and apply it to define
design rules to guide the fabrication of the relevant multilayers where
thickness t, exchange constant A, DMI constant D, magnetisation Ms and
uniaxial anisotropy Ku are the parameters to be chosen in order to stabilise
skyrmions with 10 nm diameter, stable at room temperature.

##  2 Model

The modeled material is an ultra-thin magnetic film, with thickness t
(typically 0.5 to 5 nm, along the z-axis). It is ferromagnetic at room
temperature (exchange constant A, typically 10 pJ/m), its spontaneous
magnetisation is $M_{s}$ and its uniaxial perpendicular anisotropy is
Ku(J/m3). For a thin film the demagnetising energy is characterised by
$K_{d}=\frac{\mu_{0}M_{s}^{2}}{2}$. Spontaneous perpendicular magnetisation
implies to keep the quality factor $Q=\frac{K_{u}}{K_{d}}$ larger than 1. DMI
energy is characterised by the DMI constant D(J/m2). D values larger than a
critical value $D_{c}$ favor a non ferromagnetic cycloidal ground state.

Demagnetising energy, when non-zero, is the most complicated energy to
calculate due to its non local character. For ultrathin films, a local
anisotropy assumption [14] can be proposed to rescale the anisotropy
$K=K_{u}-K_{d}$ . This assumption neglects the stabilising long range
character of the demagnetising energy. For a Néel spin structure, it also
neglects the demagnetising contribution from the charged domain wall. This
assumption will be discussed and compared to micromagnetic calculations (see
section 4.2).

Treating demagnetising energy as a local anisotropy is correct for some
recently proposed low-magnetisation systems such as ferrimagnetic films close
to their compensation composition in the chosen working temperature range [15]
or synthetic antiferromagnetic structures [16].

We assume a cylindrical symmetry for magnetisation and no z-dependence of the
magnetisation because of the ultrathin (few atomic planes) thickness compared
to the exchange length $l_{ex}=\sqrt{\frac{A}{K_{d}}}$. Magnetisation is
$\vec{M}(r)=M_{s}\vec{m}$, where $\vec{m}$ is a unit vector and $M_{s}$ is
uniform. Micromagnetic expressions are used to describe the film energy, so a
continuous magnetisation field is assumed and we restrict the study to
lengthscales larger than a few nm in order not to need an atomic model.
However, our results agree with atomic spins on a lattice calculations at the
small radius limit [17]. Due to the absence of a long range dipolar
contribution, a relevant lengthscale for this problem is the Bloch length
$l_{B}=\sqrt{\frac{A}{K}}$.

For large bubbles, $l_{B}$ scales with the Bloch domain wall width
$\delta=\pi\sqrt{\frac{A}{K}}=\pi l_{B}$. The domain wall energy in the planar
limit (large radius compared to width) for a Bloch spin structure is
$\gamma_{Bloch}=4\sqrt{AK}$, which becomes $\gamma_{DMI}=4\sqrt{AK}-\pi D$
when DMI is present assuming a Néel spin structure. It allows to define
D${}_{c}=\frac{4\sqrt{AK}}{\pi}$, which is the maximum value for D keeping a
ferromagnetic ground state. In the absence of volume magnetic poles (local
anisotropy assumption), no Neel to Bloch transition is expected when
decreasing D. A Neel structure is favored as soon as D is non zero and its
chirality is determined by the sign of D.

###  2.1 Equations

Following Kiselev et al.[11], micromagnetism energies are written using
cylindrical coordinates and assuming a cylindrical symmetry.

The ferromagnetic exchange energy density is written using the usual isotropic
term: $\epsilon_{ex}=A((\frac{\partial\theta}{\partial
r})^{2}+\frac{sin^{2}\theta}{r^{2}})$. Introducing the reduced radius
$\rho=\frac{r}{l_{B}}$, the exchange energy is:

$E_{ex}=2\pi
At\int_{0}^{\infty}((\frac{\partial\theta}{\partial\rho})^{2}+\frac{sin^{2}\theta}{\rho^{2}})\rho
d\rho$

This exchange energy is scale-invariant, so Eex does not depend on the object
radius. When it is the only energy, the minimiser is known, it is the Belavin-
Polyakov (BP) skyrmion profile [18]. Its energy is $E_{ex}=8\pi At$. Defining
$r_{0}$ as the skyrmion radius when $m_{z}(r_{0})$=0, the BP profile is
$m_{z}=\frac{r_{0}^{2}-r^{2}}{r_{0}^{2}+r^{2}}$ and
$m_{r}=\frac{2r_{0}r}{r_{0}^{2}+r^{2}}$. The
$(\frac{\partial\theta}{\partial\rho})^{2}$ planar wall term and the
$\frac{1}{\rho^{2}}$ curvature term contribute equally to the BP skyrmion
exchange energy.
The magnetocrystalline anisotropy has an interfacial origin for ultra-thin
films and it is modelled as a uniaxial one, with energy density
$\epsilon_{mc}=K\sin^{2}\theta(r)$, which gives

$E_{mc}=2\pi At\int_{0}^{\infty}\sin^{2}\theta(\rho)\rho d\rho$

In a cylindrically symmetric film, DMI is generated by the non equivalent
interfaces between for example cobalt and the buffer (Pt) and the capping
layer (AlOx). The lack of inversion symmetry makes an assymetric exchange
interaction possible. Using an atomic model, the energy can be written as
$E_{DMI}=\vec{D}_{12}.\vec{S}_{1}\wedge\vec{S}_{2}$, where $\vec{D}_{12}$ lies
in the film plane, perpendicular to the $\vec{r}_{12}$ direction between both
considered spins $\vec{S}_{1}$ and $\vec{S}_{2}$. It can be written as a
micromagnetic energy [8] :

$\epsilon_{DMI}=D(\frac{d\theta}{dr}+\frac{sin\theta cos\theta}{r})$

The first term selects the wall chirality, the second term is due to the wall
curvature and changes sign along the radial direction. It stabilises the
skyrmion in the core part ($r<r_{0}$) but it is a positive energy contribution
outside the core.
Normalising D to its maximum value allows to define
$\xi=\frac{D}{D_{c}}=\frac{\pi D}{4\sqrt{AK}}$. We are interested in $\xi$
spanning the $[0,1]$ range (i.e. 0% to 100%). We have :

$E_{DMI}=2\pi
At\int_{0}^{\infty}\frac{4\xi}{\pi}(\frac{d\theta}{d\rho}+\frac{sin\theta
cos\theta}{\rho})\rho d\rho$

Finally the skyrmion total energy is :

$E=2\pi
At\int_{0}^{\infty}[((\frac{\partial\theta}{\partial\rho})^{2}+\frac{sin^{2}\theta}{\rho^{2}})+sin^{2}(\theta(\rho))+\frac{4\xi}{\pi}(\frac{\partial\theta}{\partial\rho}+\frac{sin2\theta}{2\rho})]\rho
d\rho$ (1)

E is calculated with respect to the perpendicularly saturated ferromagnetic
state. Only demagnetising energy is non-zero and the energy density of the
perpendicularly saturated film is $\epsilon=-K_{d}$. With our assumptions,
ultrathin film (no z-dependence), cylindrical symmetry, no long range
demagnetising term, only one equation must be solved and $\xi$ is the only
parameter to determine the $\theta(\rho)$ profile. A numerical calculation was
used to solve equation (1) in the range $\xi=20\%$ to $95\%$.

Figure 1:  Magnetisation profile $\theta(\rho)$ of skyrmion solutions for
$\xi$=20% to 90%, step 10%.

###  2.2 Numerical Solution

The previous 1D equation was numerically solved using its Euler equation.

Space is discretised over 104 nodes from $\rho=0$ to $\rho_{max}=$10 or 20.
The boundary conditions are $\theta(0)=\pi$ and
$\theta(\rho_{max})=2e^{-\rho_{max}}$, which is the asymptotic tail for a 1D
Bloch wall. Energies were calculated in $2\pi At$ units. The radius $\rho_{0}$
of the skyrmion is still defined when $m_{z}=0$, so
$\theta(\rho_{0})=\frac{\pi}{2}$. We define the skyrmion wall width $\delta$
from $\theta(\rho)$ as $\delta=-\frac{\pi
l_{B}}{\frac{d\theta}{d\rho}(\rho_{0})}$

Figure 2:  Energies (exchange, anisotropy, DMI and total) of the skyrmion
solutions as a function of the skyrmion reduced radius. Equation (1) was
solved for the range $\xi$=[20%-95%], step 2.5%.

##  3 Results

The magnetisation profiles of some skyrmions are plotted in Fig(1). As
expected, a larger $\xi$ i.e., a larger DMI parameter D, stabilises a larger
radius skyrmion. The magnetisation profiles evolves from a Belavin-Polyakov
skyrmion profile, without a core, to a bubble-like skyrmion for larger $\xi$
with a core and a domain wall where $\theta(\rho)$ evidences an inflection
point. A cycloidal structure becomes the ground state when D reaches $D_{c}$
($\xi$=100%).

The detailed energies of the equilibrium structures are plotted in Fig(2) for
$\xi$=[20%-95%]. The total energy decreases with increasing $\xi$ i.e.,
increasing D. Since D is a stabilising energy, increasing its value leads to a
lower skyrmion energy even if the equilibrium size of the skyrmion is larger.
In the calculated range, the equilibrium structure is always metastable
compared to the ferromagnetic ground state where only the demagnetising energy
contributes to the total energy. At large radius, it mirrors the fact that the
effective domain wall energy stays positive whatever $\xi$ and at small
radius, a finite (BP) exchange energy is required. Stable equilibrium
structures with respect to the uniform perpendicularly magnetised state need
long-range demagnetising stabilisation energy (multidomain demagnetised
state).

Figure 3:  Energy of the skyrmion (in $2\pi At$ units) as a function of its
radius $\rho_{0}$ (in lB units). $\xi$ varies from =25$\%$ to 95$\%$. The
continuous line is $E=4(1-0.231\rho_{0})$

.

The exchange energy converges to $8\pi At$ for small radii. This is the size-
invariant BP skyrmion energy. Exchange and anisotropy energy tend to similar
linear asymptotes for large radii. This asymptote is half the Bloch domain
wall surface energy $\gamma_{Bloch}$, so it can be written as $2\sqrt{AK}2\pi
rt=4\pi At\rho_{0}$.

The DMI stabilising energy is quasi-linear with radius with a quadratic
crossover for small radii. For a converged spin structure, the DMI energy to
anisotropy energy ratio is -2, as theoretically expected (see proof, section
4.1). The skyrmion total energy (Fig.3) can be expressed in the range 20% to
50%Dc as :

$E=8\pi At(1-0.231\rho_{0})$

Figure 4: Skyrmion reduced radius $\rho_{0}$ (full circle) and skyrmion wall
width $\delta$ (cross) as a function of $\xi$. The continuous line is
$\rho_{0}=1.35\xi^{2}/\sqrt{(}1-\xi^{2})$.

Fig. 4 shows the reduced skyrmion radius $\rho_{0}$ and the wall width
$\delta$ as a function of $\xi$. The wall width $\delta$ evolves from
$\delta=2\rho_{0}$ for a BP profile (small $\xi$) up to the Bloch domain wall
width $\delta=\pi$ (in $l_{B}$ unit) for larger objects.

The skyrmion reduced radius can be expressed as
$\rho_{0}=1.35\frac{\xi^{2}}{\sqrt{1-\xi^{2}}}$, which is correct within 5% in
the range $\xi=20\%$ to $90\%$.

##  4 Discussion

###  4.1 Theoretical DMI to anisotropy energy ratio

As mentioned previously a simple relation between calculated DMI and
anisotropy energies is found. Here is the proof. Let’s assume that
$\theta(\rho)$ minimises the total energy (Equation 1), which can be written :

$E=E_{ex}+E_{K}+E_{D}$

Let’s homothetically expand radially the spin structure by a factor $\lambda$
: $\theta(\rho^{\prime})=\theta(\lambda\rho)$. The total energy becomes
$E^{\prime}=E_{ex}+\lambda^{2}E_{K}+\lambda E_{D}$. It must be larger than the
original energy which was that of the optimised spin structure.
$E^{\prime}(\lambda)$ is minimum for $\lambda=\frac{-E_{D}}{2E_{K}}$ and it
should correspond to our original spin structure, so $\lambda=1$.

Finally one has that $E_{D}=-2E_{K}$ i.e. the integrated DMI energy is always
twice the opposite to the anisotropy energy. Since the proof is based on
homogeneity, adding higher order anisotropy terms will give the same relation.
This identity can be used to assess numerical models. Note that this result is
correct for the integrated energies but not for the local energy densities.
Bogdanov et al.[19] has a similar result for a thick magnetic vortex,
invariant along z, which is mathematically equivalent.

###  4.2 Local Anisotropy Approximation

In order to take into account demagnetising terms, without calculating them,
the local anisotropy approximation was used to rescale anisotropy. In order to
test the validity of this approach for large magnetisations, a micromagnetic
solver (Mumax3) [20, 21] has been used. Fig. 5 represents the equilibrium
skyrmion radius for t=1nm, A=10pJ/m, Ku=105J/m3, D=0.76mJ/cm2, calculated on a
512x512 mesh with 1nm isotropic mesh cellsize. The blue line is the 1-D model
with an effective anisotropy $K=K_{u}-K_{d}$. When Q is larger than 3, it
superimposes the full Mumax calculation. When $M_{s}$ increases, the effective
anisotropy K decreases which makes $D_{c}$ decrease. Since D is constant,
$\xi$ goes to 100% and the skyrmion radius diverges before reaching the in-
plane to out-of-plane magnetic transition (Q=1). Boundary effects are present
in the MuMax3 calculation (512nmx512nm simulation box). Part of the finite
size effect is corrected by adding as an external field the demagnetising
field due to the saturated magnetisation out of the simulated box. However,
repulsion between the skyrmion and the DMI tilted spins at the box edges are
still present and limit the radius increase. Up to demagnetising anisotropies
$K_{d}$ which are significant compared to $K_{u}$, the analytical model is a
useful quantitative model expanding its use to non-zero-magnetisation systems.

Figure 5:  Skyrmion radius calculated as a function of the film spontaneous
magnetisation Ms keeping all parameters constant. Mumax3 micromagnetic solver
(red triangles) and our 1D model (blue circles) with effective local
anisotropy $K=K_{u}-K_{d}$. The x-axis is $M_{s}^{2}$ normalised to be 1/Q.
The chosen material parameters are A=10 pJ/m, Ku=105 J/m3, D=0.76mJ/m2, i.e.
$0.6D_{c}$ when $M_{s}$=0 and t=1nm. Increasing $M_{s}$ corresponds to
decreasing K and $D_{c}$ and increasing $\xi$ (top scale).

###  4.3 Neglecting volume charges

Representing the demagnetising energy as a local anisotropy term means
neglecting the volume charges that are associated to the Neel-domain wall spin
structure. Since the effective K model agrees with the full demagnetising
calculation (micromagnetic solver), it shows that neglecting the volume
charges does not modify significantly the spin structure. The micromagnetic
solver allows to extract in plane $H_{dr}$ and out-of-plane $H_{dz}$
components of the demagnetising field. By extracting the $H_{dz}$ and $H_{dr}$
components for a 10 nm skyrmion, as expected from the radius/thickness ratio,
there is one order of magnitude between them since $r_{0}>>10t$ and $H_{dr}$
is a minor correction. [not shown].

When the magnetisation is finite, our model neglects the expected Néel to
Bloch transition when D decreases. It happens when $E_{DMI}$ and the volume
charge contribution to $E_{demag}$ compensate. An estimate for $E_{DMI}$ is :
$E_{DMI}=D\frac{\pi}{2r_{0}}V$, where $D\frac{\pi}{2r_{0}}$ is an estimate for
the DMI density of energy and $V$ is the skyrmion volume. An estimate for
$E_{demag}$(Néel) is $K_{d}\frac{t}{r_{0}}V$, where $\frac{t}{r_{0}}$ is a
demagnetising factor for the in-plane demagnetising field. A cross over is
expected for $D_{N-B}=\frac{2}{\pi}K_{d}t\approx 0.6K_{d}t$. Recently a more
quantitative prediction $D_{N-B}=0.92K_{d}t$ using a more elaborate
calculation[22] has been reported. It is worth noticing that our rough
estimate identify the relevant physics and provides the numerical prefactor
within a factor 2. So, when Ms increases, or when D decreases, a Neel to Bloch
spin structure transition is expected. Unlike the planar domain wall case[23],
where the wall keeps its length across the Bloch to Néel transition, such a
transition for a skyrmion could be associated with a collapse of the spin
structure if the stabilising demagnetisaton energy does not compensate the
spin structure energy cost. For small values of D, this Bloch-Néel spin
structure transition is expected and our analytical model should not be used
anymore.

###  4.4 Design Rules

We will now apply the model to define design rules for stable room temperature
10-nm-diameter skyrmions. Since $r_{0}=\rho_{0}l_{B}=5$nm, and $\rho_{0}$ is
typically of order 1, a small Bloch length of 5-20 nm is desirable. The study
of the nucleation and collapse energies will help to adjust the choice of
material parameters.

First, the nucleation energy allows to create the skyrmion. Starting from a
ferromagnetic uniform state, the first skyrmion nucleus requires a finite
exchange energy $8\pi A.t$ even for a vanishingly small nucleus. This energy
does not scale with $r_{0}$ for small radii. In the atomic limit, which should
be considered to describe the initial nucleus, this conclusion still
holds[17]. Taking typical values A=10 pJ/m and t=1nm gives $8\pi
A.t=60kT_{300K}$. Such a nucleation barrier is large enough to prevent
spontaneous nucleation from thermal fluctuations at room temperature. Assisted
nucleation using a heat source (current pulse and associated Joule effect) or
a magnetic field pulse (Zeeman pressure) may be required. The nucleation
barrier could also be decreased using for example an electric gate [24]. A
nucleation pad with reduced thickness (lower t but also lower A due to finite
size effect for nanometer-thick films) could also be used to create the
skyrmion which could then be displaced to the relevant region using a spin-
polarised current induced motion.

The collapse barrier is the difference between the skyrmion metastable energy
and the energy maximum at small radii.

$E_{collapse}=8\pi A.t-8\pi A.t(1-0.23\rho_{0})=8\pi A.t.0.23\rho_{0}$

Choosing a lifetime $\tau$ of 1 second for operational use of the skyrmion and
using an Arrhenius law:

$\tau=\tau_{0}e^{\frac{\Delta E}{kT}}$

The try time $\tau_{0}$ is usually taken in the 0.1-1 ns range[17], which
gives a barrier $\Delta E=20kT$
Which brings $0.23\rho_{0}$ to be around 1/3, i.e. $\rho_{0}\approx 1$. For
10-nm diameter skyrmions, i.e. $r_{0}=5$nm, it implies $l_{B}=5$ nm. Using
typical material properties, A=10 pJ/m and K=4.105J/m3 give
$l_{B}=\sqrt{\frac{A}{K}}=5$ nm. So materials with large perpendicular
anisotropies are preferred. As regards DMI, small $\xi$ are preferred to get
small structures. Given the previous choice of exchange and anisotropy
constants, D should be kept smaller than Dc which is 2.5 mJ/m2 with these A
and K values.

What are the limits to these conclusions ? For 10-nm-diameter skyrmions, the
local character of the demagnetising field is still a good assumption. The
assumption of circular shape is even more valid for smaller objects where
fluctuations would be at the nm scale which represents large exchange costs.
Comparison to calculations performed on an atomic lattice show that the
micromagnetic conclusions do not collapse when sizes are typically 5 atomic
planes in thickness and 50 spins in diameter (typically 10nm). A more
difficult assumption is the uniform character of all magnetic properties.
Property fluctuations are less averaged out on a 10nm scale. Due to the
polycrystalline nature of most metallic layers, 10-nm skyrmions would
experience material fluctuations at a similar scale (grain size) and grain
boundaries could become efficient pinning sites.

Zeeman energy (applied magnetic field) was not considered in this work. When
$M_{s}$ is large enough and depending on the choice of a stabilising or
destabilising applied magnetic field direction, Zeeman energy could be added
to manipulate the skyrmion. It is a good external parameter to modify the
radius and to control the density of the objects.

##  5 Conclusion

The model we developed allows to discuss the size, profile, nucleation and
collapse barriers of DMI-stabilised skyrmions not using a 5-degrees of freedom
phase diagram (A, K, D, Ms, t) but simply 2 parameters, the Bloch length and
the normalised DMI. This gives a simpler picture of the role of the material
parameters, with physical arguments to discuss all trends. As a consequence,
simple rules to design skyrmions with certain size and stability are proposed.
The model agrees quantitatively with existing models and experimental results.

## Acknowledgements

#### Funding information

This project has received funding from the DARPA TEE program. M.A.M.
acknowledges an Erasmus scholarship from Universidad Autonoma de Madrid to
Université Grenoble Alpes.

## References

  * [1] P. Chaudhari, J. J. Cuomo and R. J. Gambino,  _Amorphous Metallic Films for Bubble Domain Applications_ ,  IBM Journal of Research and Development 17(1), 66 (1973),  10.1147/rd.171.0066.
  * [2] A. Bobeck, P. Bonyhard and J. Geusic,  _Magnetic Bubbles - Emerging New Memory Technology_ ,  Proc. IEEE 63(8), 1176 (1975),  10.1109/PROC.1975.9912.
  * [3] A. H. Eschenfelder,  _Magnetic Bubble Technology_ ,  Springer Series in Solid-State Sciences. Springer-Verlag, Berlin Heidelberg, 2 edn.,  ISBN 978-3-540-10790-3 (1981).
  * [4] S. Muhlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii and P. Boni,  _Skyrmion Lattice in a Chiral Magnet_ ,  Science 323(5916), 915 (2009),  10.1126/science.1166767.
  * [5] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa and Y. Tokura,  _Real-space observation of a two-dimensional skyrmion crystal_ ,  Nature 465(7300), 901 (2010),  10.1038/nature09124.
  * [6] W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen, S. G. E. te Velthuis and A. Hoffmann,  _Blowing magnetic skyrmion bubbles_ ,  Science 349(6245), 283 (2015),  10.1126/science.aaa1442.
  * [7] O. Boulle, J. Vogel, H. Yang, S. Pizzini, D. de Souza Chaves, A. Locatelli, T. O. Menteş, A. Sala, L. D. Buda-Prejbeanu, O. Klein, M. Belmeguenai, Y. Roussigné _et al._ ,  _Room-temperature chiral magnetic skyrmions in ultrathin magnetic nanostructures_ ,  Nature Nanotechnology 11(5), 449 (2016),  10.1038/nnano.2015.315.
  * [8] A. Bogdanov and A. Hubert,  _Thermodynamically stable magnetic vortex states in magnetic crystals_ ,  Journal of Magnetism and Magnetic Materials 138(3), 255 (1994),  10.1016/0304-8853(94)90046-9.
  * [9] J. Sampaio, V. Cros, S. Rohart, A. Thiaville and A. Fert,  _Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures_ ,  Nature Nanotechnology 8(11), 839 (2013),  10.1038/nnano.2013.210.
  * [10] W. J. DeBonte,  _Properties of thick‐walled cylindrical magnetic domains in uniaxial platelets_ ,  Journal of Applied Physics 44(4), 1793 (1973),  10.1063/1.1662450.
  * [11] N. S. Kiselev, A. N. Bogdanov, R. Schäfer and U. K. Rößler,  _Chiral skyrmions in thin magnetic films: new objects for magnetic storage technologies?_ ,  Journal of Physics D: Applied Physics 44(39), 392001 (2011),  10.1088/0022-3727/44/39/392001,  ArXiv: 1102.2726.
  * [12] F. Büttner, I. Lemesh and G. S. D. Beach,  _Theory of isolated magnetic skyrmions: From fundamentals to room temperature applications_ ,  Scientific Reports 8(1) (2018),  10.1038/s41598-018-22242-8.
  * [13] A. Bernand-Mantel, L. Camosi, A. Wartelle, N. Rougemaille, M. Darques and L. Ranno,  _The skyrmion-bubble transition in a ferromagnetic thin film_ ,  Scipost Physics 4(5), UNSP 027 (2018),  10.21468/SciPostPhys.4.5.027,  WOS:000437036900002.
  * [14] S. Tarasenko, A. Stankiewicz, V. Tarasenko and J. Ferré,  _Bloch wall dynamics in ultrathin ferromagnetic films_ ,  Journal of Magnetism and Magnetic Materials 189(1), 19 (1998),  10.1016/S0304-8853(98)00230-3.
  * [15] L. Caretta, M. Mann, F. Büttner, K. Ueda, B. Pfau, C. M. Günther, P. Hessing, A. Churikova, C. Klose, M. Schneider, D. Engel, C. Marcus _et al._ ,  _Fast current-driven domain walls and small skyrmions in a compensated ferrimagnet_ ,  Nature Nanotechnology (2018),  10.1038/s41565-018-0255-3.
  * [16] X. Zhang, Y. Zhou and M. Ezawa,  _Magnetic bilayer-skyrmions without skyrmion Hall effect_ ,  Nature Communications 7(1), 10293 (2016),  10.1038/ncomms10293.
  * [17] S. Rohart, J. Miltat and A. Thiaville,  _Path to collapse for an isolated Néel skyrmion_ ,  Physical Review B 93(21) (2016),  10.1103/PhysRevB.93.214412.
  * [18] A. Belavin and A. Polyakov,  _Metastable States Of 2-Dimensional Isotropic Ferromagnets_ ,  JETP Lett. 22(10), 245 (1975).
  * [19] A. Bogdanov and A. Hubert,  _The Properties of Isolated Magnetic Vortices_ ,  physica status solidi (b) 186(2), 527 (1994),  10.1002/pssb.2221860223.
  * [20] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez and B. Van Waeyenberge,  _The design and verification of Mumax3_ ,  AIP Advances 4(10), 107133 (2014),  10.1063/1.4899186.
  * [21] J. Mulkers, B. Van Waeyenberge and M. V. Milošević,  _Effects of spatially-engineered Dzyaloshinskii-Moriya interaction in ferromagnetic films_ ,  Physical Review B 95(14), 144401 (2017),  10.1103/PhysRevB.95.144401.
  * [22] A. Bernand-Mantel, C. B. Muratov and T. M. Simon,  _Unraveling the role of dipolar versus Dzyaloshinskii-Moriya interaction in stabilizing compact magnetic skyrmions_ ,  Physical Review B 101(4), 045416 (2020),  10.1103/PhysRevB.101.045416,  ArXiv: 1906.05389.
  * [23] A. Thiaville, S. Rohart, E. Jué, V. Cros and A. Fert,  _Dynamics of Dzyaloshinskii domain walls in ultrathin magnetic films_ ,  EPL 100(5) (2012),  10.1209/0295-5075/100/57002.
  * [24] T. Srivastava, M. Schott, R. Juge, V. Krizakova, M. Belmeguenai, Y. Roussigne, A. Bernand-Mantel, L. Ranno, S. Pizzini, S.-M. Cherif, A. Stashkevich, S. Auffret _et al._ ,  _Large-Voltage Tuning of Dzyaloshinskii-Moriya Interactions: A Route toward Dynamic Control of Skyrmion Chirality_ ,  Nano Letters 18(8), 4871 (2018),  10.1021/acs.nanolett.8b01502,  WOS:000441478300035.


# Diffusion across a concentration step: Strongly nonmonotonic evolution into
thermodynamic equilibrium

Hans R. Moser
Physik-Institut, University of Zürich, Winterthurerstrasse 190, CH-8057
Zürich, Switzerland
E-mail<EMAIL_ADDRESS>
###### Abstract

Dynamical and statistical behavior of the ionic particles in dissolved salts
have long been known, but their hydration shells still raise unsettled
questions. We engineered a “diffusion tunnel diode” that is structurally
analogous to the well-known Esaki diode, but now concentration gradients serve
as generalized voltages and the current means particle flow. In an
equipartition sense, the hydrated ions enter a cavity as individual particles
and later, upon increase of their concentration therein, they lose water
molecules that henceforth are particles of their own. These temporarily
attached water molecules thus are the tunnel current analogue. Unlike the
original tunnel diode, our negative differential resistance has implications
for the second law of thermodynamics, due to thermal effects of changes in the
hydration shells.

Keywords: Generalization of tunnel diode; Negative differential resistance;
Properties of hydration shells; Equipartition law; Entropy

##  I Introduction

The transfer of electrical conduction mechanisms to situations with other
driving forces has a long-standing tradition. Heat conduction driven by
temperature gradients, downhill motion of water due to gravity, and diffusion
caused by concentration gradiens (our present instance) are common examples. A
somewhat less intuitive case is magnetism, where field strength, flux density,
and permeability correspond to driving force, current density, and
conductivity, respectively.

Usually such issues are studied in a linear regime in the spirit of a
generalized Ohm’s law. However, the tunnel diode esaki  departs dramatically
from linearity, due to a current-voltage characteristic that involves a range
of negative differential resistance $dU/dI$ that is able to undamp an
oscillatory circuit. Clearly, this cannot end up in a resonance catastrophe,
since amplification stops as this range gets exceeded by the current or
voltage amplitude. In any case, there is also a positive static (that is,
ohmic) resistance contribution that dissipates energy from the power supply.

As indicated previously, we intend to translate the essentials of the tunnel
diode to the case where potential levels are established by particle
concentrations, giving rise to diffusive particle currents. At first glance
this is doomed to failure, since the tunnel current has no classical analogue
(and we disregard de Broglie wave tunneling of particles as heavy as hydrated
ions). We pursue quite another course by looking for a purely structural
correspondence that in substance is entirely different. Sometimes such
analogues are said to be isomorphic.

In the description of the experiment we shall see how the additional particle
flow (beyond ordinary diffusion) comes about. We choose NaCl for practical
reasons, among them its high solubility in water. Essentially, we consider
ions that complete their hydration shell outside of a cavity and later they
release the water molecules inside again. In short, we use that there is
maximum degree of hydration at minimum concentration of the ions. The
occurrence of a negative differential resistance then is straightforward,
since the “voltage” between a reservoir with saturated solution and the
interior of said cavity is at its maximum for zero ionic particle
concentration therein (see below that we additionally organize a zone of low
salt concentration, thus enhanced hydration, right outside the cavity). But
then there is not yet any water molecules release, and so the transfer rate of
separate (independent) water molecules into the cavity peaks at some
intermediate voltage. We are aware of other contributions that involve
negative differential resistances as well, sometimes also in quite a
generalized context with respect to the original Esaki diode, we mention the
experiments renk ; kannan . Further, Ref. li  presents a simulation not too
far from a possible experimental implementation, where the driving force is no
longer an electric field strength but a temperature gradient. However, all
these authors largely focused on issues we do not primarily have in mind.

The delicate behavior of hydration shells is already a worthy research subject
in our view, but we have an even more ambitious goal. Due to energy
conservation, hydration is attended with thermal effects, but this concerns
the (generalized) tunnel current only. Whereas, in terms of the equipartition
law, the “ordinary” particle current, that is, the flow of ionic carriers for
the polar water molecules is just made of individual particles, regardless of
their hydration status. This has exciting consequences that we intend to
investigate.

##  II Setup for generalized tunnel diode

A device with the desired properties must necessarily have an initial particle
concentration gradient, a barrier of geometrically restricted (or, conversely,
admitted) moderate diffusion conductance in order to separate regions of
different concentrations, and also a possibility to restore the initial
situation. Then, very importantly, there should be a built-in equipment to
measure small temperature differences and thus, by virtue of the hydration
heat, indirectly also concentration profiles.

Fig. 1 presents the required items and, we think, their meaningful assembly,
as will emerge later. Bottom-up we see a petri dish that serves as the
reservoir for saturated NaCl solution, a toothed structure that stands for a
grid to establish a clearance to the thermoelectric generator (TEG) just
above, and then a porous glass fiber paper for the diffusion barrier already
mentioned. Then there is a plastic frame with a square hole, covered with PTFE
tape (Polytetrafluoroethylene, often referred to as Teflon). Above that we use
a glass plate that carries a funnel containing glass granules soaked with
water by capillary forces, and also the hole or cavity underneath is filled
with plain water. All of this is fairly well isolated from the environment by
simple means, as only temperatures very close to ambient occur (see below),
and so there is little heat exchange. Uppermost there is another funnel for
(optional) water inlet, i.e., to flush or replace the water in the aforesaid
square cavity, as will be explained later. But this new external water just
enters the lower funnel, while flushing takes place with temperature-adapted
water that has been there for a long time. Moreover, the water marginally
tends to leak out between plastic frame and the moderately hydrophilic glass
plate, but the Teflon layer there reliably prevents that.

Figure 1: Left: Quasi-realistic illustration of the entire setup, see text.
Not all of the proportions are completely realistic, as we optimize for
intelligibility rather than practical construction details. In particular, we
slightly tend towards an “exploded view” to recognize the stacked items.
Right: Top view restricted to selected parts, namely glass plate that supports
the lower funnel, top rim of either funnel, as well as the thermoelectric
generator (TEG) connected to a voltmeter.

The commercial TEG is a (waterproof sealed) 13 by 13 mm square with an
efficiency of some 0.06 V/K, and the hole in the plastic frame is around the
same area and 1.5 mm in height. Such lengths matter primarily because of
diffusion velocities, as we shall see in the measurements section. This also
means that only within certain limits we might proportionally enlarge our
setup. Further, the glass fiber paper provides the barrier or obstacle to
diffusion from the NaCl solution reservoir into the cavity inside the plastic
frame. We use a filter paper around 0.45 mm thick with 93% porosity and 1.5
$\mu$m pore size in terms of particle retention, but this turns out to be
surprisingly uncritical. To see that, we tried greatly different substitutes,
among them metal meshes of various thickness and porosity, which altered
things just quantitatively. In this respect at least, our setup is quite
robust.

With these preliminaries we may now anticipate what our particular arrangement
in Fig. 1 is good for, why it is the analogue of a tunnel diode and what we
actually want to measure. We place the actual device (with lower funnel and
cavity in the plastic frame filled with plain water) into the prepared
saltwater reservoir, assemble the isolation from the surroundings and connect
the TEG to the voltmeter. Then, at the very beginning of what follows, the
temperatures at the TEG should be as balanced as possible (that is, little
offset in any sign). Once all the material is provided, this is most easily
accomplished by trial and error, as recording such TEG voltage offsets is not
overly elaborate. We have already briefly outlined what then happens, namely
diffusion out of the reservoir into the cavity, attended with thermally
relevant hydration effects. Outside the cavity, predominantly in the filter
paper out there (that has also been previously soaked with plain water), we
get a “blurred zone” with by far less salt concentration than in the
reservoir. There, the hydration increases significantly with respect to the
reservoir, which causes heat (or enthalpy) of dilution. Such dilution
processes occasionally can even be violent, as is known from the dilution of
acids. The reverse then happens in the cavity; the concentration therein grows
over time, causing the hydration to decrease and hence the temperature to
drop. Just this mechanism will ultimately be used in the subsequent figures.
In our view, it is instructive to recognize that in saturated NaCl solution
there are only some 4.5 times more water molecules than ions, and so these
“few” molecules must account for both hydration as well as liquid between the
ions. Provisionally we only note that most likely every water molecule gets
involved in hydration, at least temporarily.

Further, the geometry in Fig. 1 deserves some remarks in order to understand
the tunnel diode analogy. The hydrated particles diffuse into said cavity from
the outside and then proceed further inwards. Clearly, the generalized
potential difference is largest between the reservoir and the innermost part
of the cavity. Here the intended analogy is certainly less intuitive than in a
geometry where the particles perform a transit, i.e., enter the cavity on one
side and exit on the opposite one. However, this causes significant
engineering problems, and so we stick to the more stable and much simpler
version of Fig. 1. We merely have to accept a minor additional generalization
step with regard to the archetypal electronic tunnel diode.

##  III Measurements and their physical meaning

First, we highlight some key properties of diffusion, since they are supremely
relevant to experimental geometries in the spirit of Fig. 1 due to rates of
change in particle concentrations. At the same time, we also want to look at
the Brownian motion of individual particles, primarily because together they
make up diffusion as a whole, and also regarding the hydration effects that
rather incompletely may be described as a macroscopic overall phenomenon. The
interplay of many particles in statistical mechanics and thermodynamics is a
vast research domain, we refer to kubo ; blundell  for many relevant aspects.
We start out with the Langevin equation for a single particle

$m\ddot{{\bf x}}(t)=-\beta\dot{{\bf x}}(t)+{\bf F}_{c}(t)+{\bf F}_{ext}$ (1)

that, besides friction term and possible external force, comprises the random
driving force ${\bf F}_{c}(t)$ due to collisions. Frequently, external force
means gravitation, but here this is negligible since the hydrated Na+ and Cl-
ions do not sediment to the ground. Also, we do not apply an electric field.
The collision force is thought to be a long-term function of time, i.e., an
irregular sequence of many different collisions.

This situation may be reconsidered as a momentary view of the sum over many
particles that perform a random walk according to Eq. (1), once we accept
ergodicity to apply in this case. Conveniently, isotropy removes the random
term from this sum, and as a final prerequisite we make use of equipartition
that introduces temperature. This way, after some steps of rearrangement we
arrive at our preferred version of Eq. (1), namely an equation that also
covers the few just above stated clear-cut assumptions. This reads

$\frac{d^{2}}{dt^{2}}\langle{\bf
x}(t)^{2}\rangle_{x}+\frac{\beta}{m}\frac{d}{dt}\langle{\bf
x}(t)^{2}\rangle_{x}=\frac{2k_{B}T}{m},$ (2)

where the spatial mean value of squared distances the particles have departed
(in some fixed direction) from their starting point is our new variable. The
stationary solution of Eq. (2) is a fairly known result, namely the linearly
growing mean value

$\langle{\bf x}(t)^{2}\rangle_{x}=2Dt,\qquad D=\frac{k_{B}T}{\beta}.$ (3)

Here $D$ is the diffusion constant as it enters the two major statements for
unconstrained diffusion (in large volume liquids) according to Fick’s laws,
namely

${\bf j}_{n}({\bf x},t)=-D\nabla n({\bf x},t),\qquad\frac{\partial n({\bf
x},t)}{\partial t}=D\Delta n({\bf x},t).$ (4)

The vector ${\bf j}_{n}({\bf x},t)$ denotes the particle current density, and
Eqs. (4) describe Brownian motion as an overall diffusion phenomenon of many
particles with density $n({\bf x},t)$. We recognize that already unrestricted
diffusion constitutes a quasi-equilibrium of Brownian particles with their
host medium, since there exists a macroscopic time-evolution while the density
gradients tend to even out. Admittedly the process is slow, and so the
equipartition law stays almost unaffected. Apart from the time-evolution of
$n({\bf x},t)$ in Eqs. (4), we urgently point to the formal analogy of the
first equation there with Ohm’s law, as this is a major premise of our current
approach.

Now we are ready to present Fig. 2 that encompasses many of the issues raised
so far. Just the upper funnel in Fig. 1 for water inlet is not yet used here.
At $t=0$, the hydrated ions begin to fill up the all-important cavity in Fig.
1 after traversing the blurred zone (see above). Then, as the concentration in
this cavity increases, their degree of hydration descends and gradually
approaches the one of saturated solution. This consumes energy and the TEG
measures the corresponding temperature drop with respect to the reservoir
(note that the warming outside the cavity is not subject to our measurements).
The circles in Fig. 2 represent measured Seebeck voltages normalized to the
saturation concentration $c_{0}$ at large times, which corresponds to around
0.9 mV or $\Delta T\approx 0.015$ K across our TEG. After a characteristic
time for our chosen geometry $\tau\approx 155$ s the concentration in the
cavity amounts to $c_{0}(1-1/e)$, that is, just an $1/e$ fraction of $c_{0}$
is still missing.

Based on the just above stated characteristic time $\tau$, we can also plot
the approximate generalized voltage $\Delta c(t)=c_{0}e^{-t/\tau}$. In view of
the filter paper that acts as a particle inlet barrier, an exponential
function with measured decay constant is much more convenient and trustworthy
than direct calculation of $n({\bf x},t)$ using Eqs. (4). In sum, we record a
spatial average of concentration-dependent temperature changes (drops due to
hydration decrease) over the TEG area, measured as a momentaneous overall
thermovoltage, which then provides both curves in Fig. 2. After all, at the
top right we see that the vertical error bars (voltage measurement plus
inherent process instability for our setup) are fairly well given by the
diameter of the circles, and at small times the precision is even better.

By virtue of the diffusion dynamics, Fig. 2 depends, e.g., on the dimensions
of the TEG. We think it is worthwhile to look at Fig. 3 that renders things
now on the more natural abscissa of the generalized voltages $\Delta c$. Thus,
the exponential in Fig. 2 degenerates to a straight line with slope unity, as
this is now simply a function of itself. The most important issue is the
strongly unequal distribution of the previously equidistant circles.

Figure 2: The circles represent measured Seebeck voltages across the TEG
normalized to saturation concentration $c_{0}$, see text for experimental
procedure, efficiency as well as precision. Higher ordinate values indicate
cooling in the cavity above the TEG through progressing hydration loss of the
dissolved Na+ and Cl- ions. Also plotted is the exponentially modeled
generalized voltage $\Delta c(t)$ with a characteristic time ($1/e$ decay) of
$\approx 155$ s, as obtained from the fitted curve for the measured TEG
voltages converted into concentrations $c(t)$.  Figure 3: For the sake of
transparency, we replot the entire Fig. 2 with the natural abscissa of the
generalized voltages $\Delta c$ that normalizes away much of the geometrical
arbitrariness in the setup. The “voltages” $\Delta c$ of the ordinate then are
simply a function of themselves. Note on the left that both curves do not
start pecisely at zero voltage, as complete concentration equilibration in
Fig. 2 is only achieved asymptotically.

Only now, in Fig. 4, we come to the analogues of the currents in the Esaki
diode, namely the ordinary diffusive particle flow, and finally the tunnel
current that interests us the most. We conveniently set the volume of the
cavity to be unity, so that we do not need to distinguish particle numbers and
concentrations. Admittedly, our ordinary particle current is somewhat ill-
defined in the sense that the hydrated ions enter the cavity and then do not
instantaneously reach its innermost part, and so, in a simple linear model for
the mean flow, we consider half of the particle entries per unit time to be
the actual current. However, we emphasize that this will have no relevance at
all for our subsequent investigations, and therefore we prefer a
straightforward definition of this certainly less important part of the total
current. Then, based on the assumption that just one water molecule per ion
gets released from the hydration shell (see below), the cooling in Figs. 2 and
3 that is normalized to saturation concentration $c_{0}$ measures the number
of particles that have contributed to this thermal effect. The ordinary
current is proportional to our already stated generalized voltage $U(t)=\Delta
c(t)=c_{0}e^{-t/\tau}$, and the integral of this current over all times for
our above plain model amounts to $c_{0}/2$, which is half of the saturation
totality $c_{0}$ of particles that finally have arrived in the cavity (recall
that we equate particle numbers and concentrations, as the volume is thought
to be unity). Altogether, the ordinary (or ohmic in our generalized
perception) current as a function of time reads

$I_{ohmic}(t)=\frac{c_{0}}{2\tau}e^{-t/\tau}=\frac{U(t)}{2\tau},$ (5)

which is normalized as

$\int^{\infty}_{0}I_{ohmic}(t)dt=\frac{c_{0}}{2},$ (6)

and in Fig. 4 we merely need the proportionality factor or actually
conductance $1/(2\tau)$ in $I_{ohmic}(U)$ from the second part of Eq. (5). In
the general case, that is, no arrangement specified, this factor corresponds
to the diffusion constant $D$ in Eqs. (4). This is intuitive, since a lower
characteristic time $\tau$ increases the diffusion rate in Fig. 2 and also
means an enlarged factor $1/(2\tau)$.

Next we proceed to our major concern in this article, namely the generalized
tunnel current. On this score, we should care for the concentration-dependent
degree of hydration of the ions. What really matters is the increase in
hydration as the ion concentration drops and vice versa, while the actual
number of attached water molecules affects Fig. 4 only quantitatively. It is
not that easy to extract the quantitative extent of hydration from the
literature, since the determined hydration status of an ion significantly
depends on experimental as well as computational methods. We present a choice
of references that have particularly attracted our interest, namely (in
chronological order) Refs. soper ; onori ; barnett ; rempe ; carrillo ;
korolev ; sahle . These findings may be summerized very briefly as follows:
most authors obtain a primary shell with more strongly bound water molecules,
namely $4\pm 1$ and $2\pm 1$ of them around Na+ and Cl-, respectively, and
some authors do not distinguish inner and outer shells. In addition, an upper
limit of involved water molecules is somewhat arbitrary, as this depends on
where the too-weak interaction is truncated.

Figure 4: The two full curves are the ordinary or ohmic current (straight
line) and tunnel-only current, respectively, as described in the text. Dashed
curve: the sum of the two said currents, that is, the analogue of the total
forward current in an electronic tunnel diode which, however, has a nonlinear
$I(U)$ characteristic for the ordinary current share. The undulation to the
left is due to minor stability problems visible in Fig. 2 and becomes
accentuated by the time derivative.

As already mentioned, we assume that, averaged over all the ions, only a
single water molecule per ion finally gets released in the cavity. This number
(one molecule per ion) depends on the efficiency of our setup, e.g., how well
the blurred zone outside the cavity (see above) is able to complete the
hydration of ionic particles arriving from the reservoir. More importantly, it
strongly matters where in the hydration shell the molecules come from, and so
we actually should weight them accordingly. Thus it is not really useful to
question this number, as it depends on how it was determined. In this respect,
Fig. 4 might then be subject to quantitative changes.

And so, the generalized tunnel current is just the time derivative of the
momentaneous number of particles accumulated in the cavity (that has unit
volume, see above). We also plot the sum (dashed line in Fig. 4) of ohmic and
tunnel current in order to show the structural analogy to the “true”
electronic tunnel diode. However, this is not really what we are after. Our
main subject of investigation is the tunnel current only, and by virtue of the
measured hydration heat, this may well be separated from the ohmic
contribution that is an indirectly obtained particle flow. This is in fact a
qualitative difference from the original Esaki diode that exhibits an electric
tunnel current without such aforesaid thermal effects, just as in the
conventional current share based on conduction mechanisms of any semiconductor
diode. Note also that more than half of the abscissa in Fig. 4 is attended
with a negative differential resistance $dU/dI$ for the tunnel current alone
curve.

The thermal effects are completely decoupled from the ohmic governed particle
current, since the released and now independent water molecules in the cavity
did not come in as individual particles from the outside, but were an attached
part of a hydrated ion. However, this is of paramount importance, as this
implies a non-dissipative heat flow out of the cavity, giving rise to the
measured cooling inside. Hence, the range of negative differential resistance
may be harnessed to undamp thermal oscillations that do not dissipate energy
from the driving power supply. Clearly, this range should never be exceeded by
the (generalized) voltage or current amplitude in the spirit of Fig. 4, and so
we must permanently harvest the thermal oscillation energy. We think this is a
possible second law of thermodynamics violation, as the said heat flow out of
the cavity spontaneously constitutes a temperature difference, and the input
is only the externally powered particle concentration gradient that does not
involve the thermal energy attended with unbalanced temperatures. This must
then be organized recurrently in order to maintain the oscillations.

We show this relationship between thermal and concentration oscillations in
our Fig. 5, which this time also involves the upper funnel in the Fig. 1
setup. Again, the circles are measured points, and the fit still features
moderate smoothing, as is easiest seen at the tops of the three structures
(maximum cooling in the cavity). The steep rises indicate the flushing events,
that is, we pour 0.2 cm3 of plain water into the upper funnel. The basic
level, thus apart from the peaks, compares to the 0.9 mV of Fig. 2, but the
three peak heights indicate lower temperatures. This is expected, as the
liquid in the cavity is somewhat cooler than the one right above the TEG,
where thermal conduction out of the reservoir has an influence (however,
thermal equilibration in Fig. 2 takes place on a much longer time scale). The
cooler saltwater from the bulk region of the cavity then flows over the TEG
and causes the measured “cold shocks”. A larger amount of flushing water
results in even stronger shocks and, also importantly, improves the salt
depletion in the blurred zone (see above) outside the cavity. But, in a
cautious attitude, we stay with an amount somewhat less than the cavity
volume.

Figure 5: The three steep ascents are flushing processes at well-taken time
lags in order to recurrently restore the initial situation, see text. This
time we directly record the TEG voltages that measure the cooling with respect
to the reservoir, where every flushing event brings along the cooler liquid
from the interior or bulk area of the cavity under consideration. Again, the
circles are measured points and the fit comprises moderate smoothing, thus
slightly understates the peak heights.

The actual peak heights depend, among other things, on the current conditions
in the upper funnel. Further, even their long-term average is of minor
importance, since this also depends on the geometrical details in the overall
setup. Fig. 5 is just intended to demonstrate how the initial situation with a
concentration drop between reservoir and the cavity interior (extended by the
blurred zone) can be organized in a periodically recurrent manner.

Moreover, it may be helpful to mind that the situation in Fig. 5 withstands a
complete reversal of concentrations, that is, the reservoir now contains plain
water, and the cavity as well as lower funnel are initially filled with
saturated NaCl solution. Most importantly, “flushing” is also performed with
saturated solution. Then, the liquid inside the cavity gets diluted upon
concentration equilibration and thus warms up. The flushing processes cause
the warmer liquid from the more central area of the cavity to flow over the
TEG again, and so we experience “heat shocks” of similar magnitude as the
previous coolings.

We have already touched on the possibility of exceptions to the second law,
and therein we disregard the issues of a Maxwell’s demon and of time reversal.
To our knowledge, there are no such exceptions in statistical mechanics based
on point-like particles, that is, in a somewhat hypothetical world where these
points are nevertheless realistically spaced for fluids or solids (the ideal
gas is a common example). Interestingly, the domain of continuum mechanics (no
inner structure or constituents) admits extensions where these exceptions or
violations can occur, see ostoja . In the present article, however, we pursue
quite another course, namely we scrutinize some rather hidden properties of
our working medium, hence of the dissolved hydrated Na+ and Cl- ions.

Regarding particle size only, we summarized our ideas in moser , where we
studied a porous medium with hollow sizes that are reasonably related to the
particle extent. There we experienced small but undoubtedly measurable
entropic consequences. What is most striking in such measurements is the
dependence of the diffusion constant $D$ on the geometrical conditions, as is
shown particularly clearly in lancon . Essentially, we dealt with
“thermodynamics” of small systems, namely the pores or hollows in such media,
and these weakly coupled subsystems then add up to a macroscopic volume. For
an interesting perspective on such small systems, see sekimoto  and very
recently watanabe . We point to linear response theory that provides a
convenient way to theoretically tackle such issues, as the forces of the inner
walls on the particles may be naturally incorporated into a generalized
external (scalar) force $f(t)$. Let $x(t)$ be a (for the sake of simplicity
again scalar) unperturbed system’s signal. We state the progressing time-
average of its response to said generalized force as

$\langle
x_{resp}(t)\rangle_{t}=\int^{\infty}_{-\infty}\Theta(t-\tau)\chi(t-\tau)f(\tau)d\tau.$
(7)

Then, the Fourier transform $\tilde{\chi}(\omega)$ of the linear response
function or susceptibility $\chi(t)$ in Eq. (7) enters

$\tilde{P}(\omega)=\frac{2k_{B}T}{\omega}\rm{Im}\tilde{\chi}(\omega)$ (8)

that is the fluctuation-dissipation theorem, see kubo ; blundell . The
Heaviside function $\Theta(t-\tau)$ in Eq. (7) assures causality with respect
to the effect of $f(t)$, and $\tilde{P}(\omega)$ in (8) means the Fourier
transform of $x(t)^{2}$ or power spectral density. For the Brownian particles
that make up diffusion, theorem (8) states that the dissipation caused by
$f(t)$, thus essentially $\rm{Im}\tilde{\chi}(\omega)$, has the same physical
origin as the Brownian fluctuations (characterized by the $\tilde{P}(\omega)$
distribution) themselves. Further, depletion effects are a fairly known
research domain where particle size may cause entropic phenomena, see, e.g.,
karzar . Clearly, in our present approach there are no such peculiarities of
porous media (apart from the diffusion barrier between reservoir and cavity,
see above), and also depletion effects are minor. However, we point to the
strong overlap with our current work, as nonzero particle size is really
crucial. A theory where the water molecules for a hydration shell just
disappear into point-like particles is, in our view, somewhat strange.

A final remark is devoted to Jarzynski’s equality jarzynski  that bypasses the
inequality in the second law of thermodynamics; in plain terms, entropy tends
to increase (or remains constant in the Hamiltonian, that is, conservative
limit). The statement reads

$\langle e^{-W/k_{B}T}\rangle=e^{-\Delta F/k_{B}T}$ (9)

and relates the Helmholtz free energy difference $\Delta F$ to the work $W$
supplied to a system, regardless whether or not we slowly pursue a reversible
path of quasi-static states starting at equilibrium situation of reservoir
temperature $T$. Here, in $F=U-TS$ with $U$ being the internal energy, the
term due to dissipation $TS$ of entropy $S$ may well grow. The inequality
(entropy remains or rises, or also $\langle W\rangle\geq\Delta F$) is removed,
and this actually supports strict validity of the second law that this way
might relate to other physical statements. Jarzynski’s equality goes on to be
an active research field, and exciting generalizations have also emerged, see,
e.g., Ref. knipschild . However, just as in the particle size discussion
above, this is not our primary goal, since we tackle the scope of the second
law from a different perspective. By virtue of hydration, our particles have
changed size and mass along their way, and this is not covered in Eqs. (8) and
(9), neither are the attended thermal effects that are the major subject of
this article. Nevertheless, the particles are always small enough that
equipartition holds true.

##  IV Conclusions and outlook

Overall, we have investigated an analogue of the tunnel diode, with the
driving force now being a concentration gradient of dissolved ionic particles
(we used NaCl in the present work). Clearly, to reproduce the negative
differential resistance is imperative, and we point out its long and steep
range in Fig. 4. Diffusion processes are notoriously slow, but we envision
thermal applications where this is not overly disadvantageous.

The tunnel current analogue consists of water molecules of a hydration shell
that in a sense are “stowaways”, as they later get released in a cavity (upon
growth of concentration therein). But this process, only after a maximum
transfer rate of such water molecules, gradually comes to an end due to
concentration equilibration, which also explains the negative differential
resistance.

Here it is of great importance that the thermal hydration effects only concern
the generalized tunnel current, and not the usual (in our analogy ohmic)
particle flow attended with dissipation. This is a beneficial additional
feature of our approach compared to the original Esaki diode where the tunnel
current cannot be separated from the conventional forward current. Definitely,
this has entropic consequences, as the resulting heat flow out of the cavity
(see above) is opposite to the inward particle flow driven by a concentration
gradient. To our knowledge, the interdependence between hydration heat and a
negative differential resistance has not yet been investigated, or maybe
rather marginally and indirectly.

There is certainly a lot of potential to improve our setup, above all since we
used the energy in a hydration shell only to a small extent. On the
theoretical side, we think statistical mechanics of such over time varying
hydrated ions has not been studied particularly well either, as hydration
involves both statistics and deterministic dynamics. Molecular dynamics
simulations may be a promising tool to anticipate the theory, but there are
still many unsettled questions.

Acknowledgments

The author thanks P. Robmann and U. Straumann for constructive discussions.

## References

  * (1) L. Esaki, Phys. Rev. 109, 603 (1958).
  * (2) K. F. Renk and B. I. Stahl, Phys. Lett. A 375, 2644 (2011).
  * (3) V. Kannan, H.-S. Kim, and H.-C. Park, Phys. Lett. A 380, 3743 (2016).
  * (4) F. Li, J. Wang, G. Xia, and Z. Li, Nanoscale 11, 13051 (2019).
  * (5) R. Kubo, M. Toda, and N. Hashitsume, Statistical Physics II (Nonequilibrium Statistical Mechanics), 2nd ed. (Springer, Berlin, 2003).
  * (6) S. J. Blundell and K. M. Blundell, Concepts in Thermal Physics, 2nd ed. (Oxford University Press, New York, 2010).
  * (7) A. K. Soper, G. W. Neilson, J. E. Enderby, and R. A. Howe, J. Phys. C: Solid State Phys. 10, 1793 (1977).
  * (8) G. Onori, J. Chem. Phys. 89, 510 (1988).
  * (9) R. N. Barnett and U. Landman, Phys. Rev. Lett. 70, 1775 (1993).
  * (10) S. B. Rempe and L. R. Pratt, Fluid Phase Equilibria 183-184, 121 (2001).
  * (11) M. Carrillo-Tripp, H. Saint-Martin, and I. Ortega-Blake, J. Chem. Phys. 118, 7062 (2003).
  * (12) V. P. Korolev, J. Structural Chem. 52, 737 (2011).
  * (13) C. J. Sahle, E. de Clermont Gallerande, J. Niskanen, A. Longo, M. Elbers, M. A. Schroer, C. Sternemann, and S. Jahn, Phys. Chem. Chem. Phys. 24, 16075 (2022).
  * (14) M. Ostoja-Starzewski, Entropy 19, 78 (2017).
  * (15) H. R. Moser, arXiv:1602.04249v1 (2016).
  * (16) P. Lançon, G. Batrouni, L. Lobry, and N. Ostrowsky, Physica A 304, 65 (2002).
  * (17) K. Sekimoto, F. Takagi, and T. Hondou, Phys. Rev. E 62, 7759 (2000).
  * (18) G. Watanabe and Y. Minami, Phys. Rev. Res. 4, L012008 (2022).
  * (19) M. Karzar-Jeddi, R. Tuinier, T. Taniguchi, and T.-H. Fan, J. Chem. Phys. 140, 214906 (2014).
  * (20) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  * (21) L. Knipschild, A. Engel, and J. Gemmer, Phys. Rev. E 103, 062139 (2021).