# Hybridization driving distortions and multiferroicity in rare-earth
nickelates

Luca Binci  Theory and Simulation of Materials (THEOS), and National Centre
for Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Michele
Kotiuga  Theory and Simulation of Materials (THEOS), and National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Iurii
Timrov  Theory and Simulation of Materials (THEOS), and National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Nicola
Marzari  Theory and Simulation of Materials (THEOS), and National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland  Laboratory
for Materials Simulations, Paul Scherrer Institut, 5232 Villigen PSI,
Switzerland

###### Abstract

For decades transition-metal oxides have generated a huge interest due to the
multitude of physical phenomena they exhibit. In this class of materials, the
rare-earth nickelates, $R$NiO3, stand out for their rich phase diagram
stemming from complex couplings between the lattice, electronic and magnetic
degrees of freedom. Here, we present a first-principles study of the low-
temperature phase for two members of the $R$NiO3 series, with $R=$ Pr, Y. We
employ density-functional theory with Hubbard corrections accounting not only
for the on-site localizing interactions among the Ni–$3d$ electrons ($U$), but
also the inter-site hybridization effects between the transition-metals and
the ligands ($V$). All the U and V parameters are calculated from first-
principles using density-functional perturbation theory, resulting in a fully
_ab initio_ methodology. Our simulations show that the inclusion of the inter-
site interaction parameters $V$ is necessary to simultaneously capture the
features well-established by experimental characterizations of the low-
temperature state: insulating character, antiferromagnetism and bond
disproportionation. On the contrary, for some magnetic orderings the inclusion
of on-site interaction parameters $U$ alone completely suppresses the
breathing distortion occurring in the low-temperature phase and produces an
erroneous electronic state with a vanishing band gap. In addition – only when
both the U and V are considered – we predict a polar phase with a
magnetization-dependent electric polarization, supporting very recent
experimental observations that suggest a possible occurrence of type-II
multiferroicity for these materials.

##  I Introduction

Rare-earth nickelates _R_ NiO3 (_R_ = rare-earth element) are a class of
materials exhibiting a particularly rich phase diagram Medarde (1997). Except
LaNiO3, these materials undergo a metal-insulator phase transition (MIT)
coinciding with a symmetry-lowering distortion of the crystal structure Alonso
_et al._ (1999, 2000). This distortion changes the space group from $Pbnm$, in
which all Ni sites are crystallographically equivalent, to $P2_{1}/n$ – which
possesses two inequivalent Ni sites – because of the activation of a breathing
mode, that contracts and expands the NiO6 octahedra in a rock-salt pattern.
The MIT temperature, $T_{\mathrm{MIT}}$, is a monotonically decreasing
function of the size of the rare-earth cation
($\text{Lu}\rightarrow\text{La}$), that coincides with the Nèel temperature
$T_{\text{Nèel}}$ for $R$ = Nd, Pr, and vanishes for $R$ = La. In the last
years rare-earth nickelates have attracted a great deal of attention,
especially for their application within epitaxial heterostructures Middey  _et
al._ (2016), which have been shown to possess crossovers between different
magnetic orders Hepting  _et al._ (2018), the possibility to manipulate
$T_{\mathrm{MIT}}$ through electric fields Scherwitzl  _et al._ (2010), and
exotic types of phases transitions Post  _et al._ (2018); Liu  _et al._
(2013). Recently, superconductivity has been induced via apical oxygen
deintercalation together with Sr doping in $R$NiO3 with $R$ = Nd, Pr Li  _et
al._ (2019); Osada  _et al._ (2020). Finally, recent experiments have
suggested a low-temperature multiferroic phase signaled by additional Raman-
active modes compatible with the loss of the inversion center in the
antiferromagnetic (AFM) phase Ardizzone  _et al._ (2021).

Figure 1: (a) Unit cell of $R$NiO3 in the $P2_{1}/c$ space group doubled along
the direction of the lattice vector $\mathbf{a}$ in order to accommodate the
different AFM orderings studied here. Blue (orange) octahedra contain the
short (long)-bond nickels, NiS (NiL), respectively. The grey and the green
atoms are the rare-earth elements ($R$ = Pr or Y) and the O atoms,
respectively. The other panels show the three AFM orderings investigated: (b)
$T$-type, (c) $B$-type, and (d) $S$-type; here, for clarity, only the Ni atoms
are shown. The arrows represent the directions of the magnetic moments on the
Ni atoms and their lengths are proportional to their modulus. In the
$S/T$-type orderings all of the Ni atoms have a magnetic moment, while in the
$B$-type ordering the nickels with shorter Ni$-$O bond lengths are
nonmagnetic. The green and red planes are guides for the eyes to visualize
planes of spin-up and spin-down moments, respectively. Rendered using VESTA
Momma and Izumi (2008).

Besides the MIT, the nature of the low-temperature insulating phase has been
the subject of intense study. An orbital ordering induced by a Jahn-Teller
distortion was ruled out in favour of charge disproportionation Mazin  _et
al._ (2007); this picture had also reached experimental consensus Alonso  _et
al._ (1999); Medarde  _et al._ (2008). However, such charge disproportionation
has not manifested in first-principles calculations Park  _et al._ (2012),
which show a negligible difference in the $3d$ electronic occupations of the
two crystallographically inequivalent NiS and NiL (defined below), both of
which have approximately 2 electrons in the $e_{g}$ states, and therefore 8
electrons in the $3d$ manifold – an effect clearly explained within the
negative-feedback charge regulation mechanism discussed in Ref. Raebiger  _et
al._ (2008). Hereafter we use the subscripts S/L, to indicate the Ni ion in
the smaller/larger NiO6 octahedrons, which are contracted/expanded by the
breathing-mode distortion. As rare-earth nickelates are considered to be
negative charge-transfer insulators Zaanen  _et al._ (1985), it has been
proposed that, starting from a configuration with one ligand hole
($\underline{L}$) on every nickel, the transition
$(3d^{8}\underline{L})\,(3d^{8}\underline{L})\rightarrow(3d^{8}\underline{L}^{2})_{\mathrm{S}}\,(3d^{8})_{\mathrm{L}}$
takes place at the MIT Johnston  _et al._ (2014). According to this picture,
one NiO6 octahedron contracts and the two $e_{g}$ electrons on NiS couple with
the two oxygen holes, forming a low-spin state. The other NiO6 octahedron
expands and ends up in a high-spin state owing to the Hund’s coupling. Similar
conclusions have been drawn from an analysis based on density-functional
theory (DFT) + dynamical mean field theory (DMFT) Park  _et al._ (2012); Haule
and Pascut (2017) (that named the state a “site-selective” Mott state),
DFT+$U$ studies Varignon  _et al._ (2017); Hampel and Ederer (2017), multi-
band many-body Hamiltonians solved with Hartree-Fock Mizokawa  _et al._
(2000); Johnston  _et al._ (2014) and exact diagonalization methods Green  _et
al._ (2016).

Despite such intense theoretical effort, several open questions still remain.
In particular, the precise magnetic order remains elusive. Indeed, there is
wide experimental agreement only on the magnetic propagation vector
$\mathbf{Q}$, which encodes the periodicity of the magnetic order with respect
to the primitive cell of the space group. However, the precise values of the
magnetic moments on the Ni atoms have not yet been unambiguously determined
García-Muñoz  _et al._ (1994); Alonso  _et al._ (1999). From first-principles
calculations it has been reported that only small values of Hubbard $U\lesssim
2$ eV successfully stabilize an AFM order that is both compatible with the
experimental Q and energetically favourable, with respect to the ferromagnetic
(FM) order Varignon  _et al._ (2017); Hampel and Ederer (2017); Mercy  _et
al._ (2017). Still, the values for the Hubbard parameter were tuned by
scanning different values, and were not calculated rigorously with a well-
determined theoretical protocol, thus preventing predictions to be entirely
nonempirical. More importantly, these calculated magnetic orders have
systematically zero magnetic moments on NiS Badrtdinov  _et al._ (2021). As it
is shown here, this feature is not compatible with a possible multiferroic
nature for $R$NiO3, which was first computationally predicted Giovannetti  _et
al._ (2009) and more recently experimentally highlighted by Raman spectroscopy
measurements, which revealed the emergence of additional phonon modes in the
low-temperature phase, compliant with an inversion symmetry breaking required
for a type II multiferroic phase Ardizzone  _et al._ (2021). In order to
computationally demonstrate such multiferroic phase, it is necessary to
establish a coupling between a nonvanishing electric polarization P and the
magnetic degrees of freedom, i.e. to prove a dependence of P on the simulated
magnetic order. A group theory analysis Perez-Mato  _et al._ (2016) has shown
that a nonzero P can be obtained if the magnetic moments on the NiS atoms are
different from zero, for instance within the $S$\- and/or $T$-types magnetic
orderings Giovannetti  _et al._ (2009) (see Fig. 1). Unfortunately, for
$U\simeq 2$ eV these latter magnetic patterns are stable only in a range of
amplitudes for the breathing mode distortion that are much smaller than the
experimental ones Badrtdinov  _et al._ (2021). Larger values of $U$ may lead
to converged $S/T$-types orderings, but they could also cause a FM instability
Varignon  _et al._ (2017); Hampel and Ederer (2017).

In this paper we use a fully first-principles DFT+$U$+$V$ approach that
includes both on-site ($U$) interactions to capture the electronic
localization and inter-site ($V$) interactions to account for hybridization
effects Campo and Cococcioni (2010). We calculate the Hubbard parameters self-
consistently using density-functional perturbation theory (DFPT) Timrov  _et
al._ (2018, 2021), thus resulting in a parameter-free theory. Our calculated
Hubbard $U$ parameters range from 8 to 9 eV and are numerically different for
the two crystallographically inequivalent
$\mathrm{Ni}_{\mathrm{S}}/\mathrm{Ni}_{\mathrm{L}}$. Using these high values
of $U$, after an iterative approach aimed to determine the self-consistent
crystal structure and Hubbard parameters, we find that within DFT+$U$ Anisimov
_et al._ (1991); Dudarev  _et al._ (1998) – having only on-site electronic
contributions – the investigated $R$NiO3 structures would be found in a FM
ground state, in agreement with other theoretical studies Varignon  _et al._
(2017); Hampel and Ederer (2017) and in contrast with experimental
observations. Moreover, the DFT+$U$ calculated AFM orders of the $S/T$-types
are also found to have a vanishing breathing mode distortion and zero band
gap, in disagreement with experiments. Conversely, by including the inter-site
Hubbard $V$ interactions between the Ni and O sites, the results change
considerably: in particular, we find qualitatively different final states for
the $S/T$ orderings. In addition to the fact that some of them exhibit the
expected AFM instability against FM ordering, with DFT+$U$+$V$ we always
recover an insulating character, consistently with experiments. This feature
allows us to evaluate the electric polarization P, which is found to be
magnetization-dependent, thus providing additional support for the emergence
of multiferroicity Giovannetti  _et al._ (2009); Ardizzone  _et al._ (2021).
We also simulate the magnetic order with zero magnetic moments on NiS (i.e.
the $B$-type ordering, see Fig. 1) and show that this is not able to generate
a polar instability, thus keeping the materials in a non-polar phase with a
vanishing P. Our results for PrNiO3 and YNiO3 – sitting on opposite ends of
the $R$NiO3 phase diagram – are qualitatively similar, suggesting that these
numerical findings could also be extended to other members of the rare-earth
nickelate series.

Table 1:  Structural and electronic properties of PrNiO3 and YNiO3 in the low-temperature phase for different magnetic orderings computed with different DFT+Hubbard schemes and as measured in experiments. Positive $\Delta E$ means that FM is lower in energy than AFM, and vice versa. Three types of the AFM order are considered: $T$, $S$, and $B$ types (see Fig. 1). The $V$ interaction parameters are different within a given Ni$-$O octahedron because of the different Ni$-$O bond lengths yielding three inequivalent sites (O1, O2 and O3). However, the differences between the different values of $V$ within a given octahedron are very small and for conciseness we report their average values; The full results are available in the Materials Cloud Archive Binci  _et al._ . The metallic (M) or insulating (I) character for each case is also indicated. Method | AFM |  $U_{\mathrm{S}}/U_{\mathrm{L}}$ (eV) |  $\langle V_{\mathrm{S}}/V_{\mathrm{L}}\rangle_{\mathrm{avg}}$ (eV) |  $m_{\mathrm{S}}/m_{\mathrm{L}}$ ($\mu_{\mathrm{B}}$) |  $\langle\ell_{\mathrm{S}}/\ell_{\mathrm{L}}\rangle_{\mathrm{avg}}$ (Å) |  $\Delta E$ (meV/f.u.) | Met/Ins
---|---|---|---|---|---|---|---
PrNiO3
DFT+$U$ | $S$ | 9.40/9.40 |  | 1.54/1.53 | 2.00/2.00 | 156 | M
| $T$ | 9.32/9.33 |  | 1.52/1.53 | 1.99/2.00 | 137 | M
| $B$ | 9.19/8.19 |  | 0.00/1.68 | 1.88/2.03 | 200 | I
DFT+$U$+$V$ | $S_{(a)}$ | 8.76/9.61 | 0.76/0.90 | 0.64/1.65 | 1.89/2.03 | $-26$ | I
| $T_{(a)}$ | 8.77/9.60 | 0.76/0.90 | 0.69/1.68 | 1.90/2.03 | $-23$ | I
| $S_{(b)}$ | 9.78/9.33 | 1.02/0.95 | 1.14/1.61 | 1.94/2.02 | 71 | I
| $T_{(b)}$ | 9.78/9.42 | 1.01/0.96 | 1.15/1.61 | 1.94/2.02 | 75 | I
| $B$ | 9.60/8.49 | 0.84/0.71 | 0.00/1.67 | 1.88/2.03 | 90 | I
Expt. Gawryluk  _et al._ (2019) |  |  |  |  | 1.92/1.97 |  | I
YNiO3
DFT+$U$ | $S$ | 9.38/9.38 |  | 1.51/1.51 | 2.00/2.00 | 115 | M
| $T$ | 9.21/9.20 |  | 1.49/1.49 | 2.00/2.00 | 72 | M
| $B$ | 9.17/8.06 |  | 0.00/1.68 | 1.88/2.04 | 143 | I
DFT+$U$+$V$ | $S_{(a)}$ | 8.38/9.49 | 0.65/0.80 | 0.15/1.70 | 1.88/2.05 | $-94$ | I
| $T_{(a)}$ | 8.40/9.50 | 0.66/0.81 | 0.18/1.70 | 1.88/2.05 | $-91$ | I
| $S_{(c)}$ | 9.31/9.25 | 0.82/0.79 | 1.32/1.32 | 1.99/1.99 | 7 | I
| $T_{(c)}$ | 9.35/9.30 | 0.84/0.82 | 1.33/1.33 | 1.99/1.99 | 6 | I
| $B$ | 9.54/8.35 | 0.80/0.68 | 0.00/1.67 | 1.88/2.04 | 61 | I
Expt. Alonso  _et al._ (2001) |  |  |  |  | 1.92/2.01 |  | I

##  II Results

The AFM magnetic orderings studied in this paper are shown in Fig. 1. All of
these require a supercell of the 20-atom unit cell of the monoclinic
$P2_{1}/n$ structure, which in turn is a $\sqrt{2}\times\sqrt{2}\times 2$
supercell of the 5-atom cubic perovskite structure. Previously these orderings
have been investigated using a 80-atom 2$\times$1$\times$2 supercell of the
$P2_{1}/n$ space group, in order to reproduce the experimentally observed AFM
propagation vector
$\mathbf{Q}_{P2_{1}/n}=(\nicefrac{{1}}{{2}},0,\nicefrac{{1}}{{2}})$ Hampel and
Ederer (2017); Varignon  _et al._ (2017). To reduce the computational costs,
we have used the $P2_{1}/c$ setting of the space group 14, which has the same
$\mathbf{a}$ and $\mathbf{b}$ primitive lattice vectors as $P2_{1}/n$, but
$\mathbf{c}_{P2_{1}/c}=\mathbf{a}+\mathbf{c}_{P2_{1}/n}$ Ardizzone  _et al._
(2021). In this setting the magnetic propagation vector becomes
$\mathbf{Q}_{P2_{1}/c}=(\nicefrac{{1}}{{2}},0,0)$, only requiring a doubling
of the 20-atoms $P2_{1}/c$ unit cell (instead of quadrupling it, as for the
20-atoms $P2_{1}/n$).

In Table 1, we report for each of the AFM configurations the results of the
calculated magnetic moments $m$, total energy difference with respect to the
FM configuration $\Delta E$, the average bond lengths for the Ni$-$O octahedra
$\langle\ell\rangle_{\mathrm{avg}}$, and the computed values of the Hubbard
interaction parameters $U$ and $V$. As shown in Table 1, DFT+$U$+$V$ yields
two states for each of the two $S/T$ types of orderings. Both of them are
energetically degenerate in pairs (within numerical accuracy): the lowest-
energy ones, that we denote as $(S/T)_{(a)}$, and the higher-energy
(metastable) states named $(S/T)_{(b)}$ and $(S/T)_{(c)}$. As described in the
Methods section these states are found with an iterative self-consistent loop
alternating calculations of Hubbard parameters and crystal structure
optimizations until convergence of the two is reached simultaneously. In
general we find stable points within this workflow, except for the
$(S/T)_{(a)}$ states. For these latter, at each iteration the protocol
produces a periodic reversal of the breathing mode amplitude, to which
equivalent switches of the Hubbard parameters and magnetic moments are
associated. Since the numerical values of $U$, $V$, Ni–O octahedral volumes
and magnetic moments associated to NiS/NiL are numerically always the same,
together with the fact that the total energy does not exhibit any variation
but a smooth convergence with the iterations, we argue that this does not
create physical ambiguities, as all the main features of the final result are
well-defined (for more information, see Supplementary Information Sup ). This
fact, however, underlines the need of a formulation for Hubbard-corrected DFT
functionals able to take into account the variation of the Hubbard parameters
with respect to atomic positions and strains, inasmuch as this dependence can
be quite significant Kulik and Marzari (2011).

Figure 2: pDOS of the relaxed crystal structure for each magnetic ordering.
(a), (b): B-type and $T$-type, respectively, for PrNiO3. (c), (d):
$T_{(a)}$-type and (e), (f): $T_{(b)}$-type and $T_{(c)}$-type, for both
PrNiO3 and YNiO3. The Hubbard functionals employed to produce (a), (c)–(f)
include both the $U$ and the $V$ interaction parameters, while in (b) only the
on-site $U$ is present. In the plots the Ni–$e_{g}$ states and the full O–$2p$
projection are shown. We averaged over inequivalent O atoms and over the two
Ni($e_{g}$) states for each spin channel. Given the symmetry between spin-
up/spin-down channels for AFM orderings, the notation for majority/minority
(maj/min) spin channels is: for a Ni with positive magnetic moments ($m>0$
$\mu_{\mathrm{B}}$), $\sigma_{\mathrm{maj}}=\,\uparrow$ and
$\sigma_{\mathrm{min}}=\,\downarrow$; vice versa for Ni with $m<0$
$\mu_{\mathrm{B}}$. The zero of energy corresponds to the valence band
maximum.

The total energies reported in the last column of Table 1 are calculated at
the relaxed crystal structure with the corresponding self-consistent, site-
dependent $U$ and $V$ Hubbard parameters, both for the AFM and the FM
configurations. Since the final $U$ and $V$ are not numerically the same, it
follows that the reference energies are different. This could in principle
create ambiguities when comparing total energies calculated with functionals
containing numerically different parameters. However, we regard the Hubbard
term as a correction to the self-interaction error of the DFT base functional,
selectively applied to the Ni–$3d$ (and O–$2p$) manifold of the corresponding
ground state. According to this point of view it is natural that different
electronic ground states will be characterized by different self-interaction
errors (i.e. different curvatures with respect to the $3d$ electron occupancy,
see Ref. Cococcioni and de Gironcoli (2005)). As a consequence we follow Ref.
Cococcioni and Marzari (2019) and argue that, for each of electronic state,
the corresponding $U$ and $V$, aimed to remove such errors, must necessarily
be different.

We highlight that we use orthogonalized atomic states taken from the
pseudopotential as Hubbard projectors (see Methods), while the calculations
with $U\simeq 2$ eV by other works mentioned above in the text were carried
out with orbital projectors derived from the projector augmented-wave (PAW)
formalism Bengone  _et al._ (2000). In general same values for the Hubbard
parameters can have a different impact on the electronic structure when
applied to different Hubbard projectors Wang  _et al._ (2016). Therefore, to
establish the consistency of two computational setups in the same range of $U$
parameters, we have numerically verified for PrNiO3 that with the empirical
$U\simeq 2$ eV we stabilize a $B$-type AFM order – even starting from the
$S/T$-type orderings – that is lower in energy than the FM configuration, in
compliance with previous studies carried out with PAW projectors.

In Fig. 2 we show the projected density of states (pDOS) of the $T_{(a)}$ and
$T_{(b)}$ types of orderings with and without the inter-site $V$, and the
$B$-type ordering computed using DFT+$U$+$V$. Equivalent results for the
$S$-type ordering are almost identical to the $T$-type ordering. Furthermore,
the results for the $B$-type ordering obtained with and without inter-site $V$
present no qualitative differences and, thus for sake of conciseness we report
only the DFT+$U$+$V$ results. In particular we plot the Ni–$e_{g}$ and the
O–$2p$ states. The Ni $t_{2g}$ bands, corresponding to localized, fully-
occupied states, yield narrow peaks located below $-5$ eV with respect to the
highest occupied states and for clarity are not reported in the same figure
(they can be found in SI Sup ).

In Table 2 we show the results of the calculated electric polarization P using
the Berry phase approach King-Smith and Vanderbilt (1993); Umari and
Pasquarello (2002). Because of the polar instability, during the structural
optimization the crystal structure lowers its symmetry, and thus the space
group changes. The newly obtained space groups are reported in the same table.
For the $S/T$ orderings, we mention that we obtained a nonzero P also at the
experimental crystal structure, in agreement with Ref. Giovannetti  _et al._
(2009), as the inversion symmetry-breaking is induced at the level of the
magnetic space group.

##  III Discussion

###  III.1 Electronic, crystal and magnetic structures

#### DFT+$U$

The $B$-type AFM ordering presents no relevant differences when examined with
the two methods, i.e. with or without the inter-site parameter $V$: they show
a band gap of $\simeq 1.4$ eV, very similar magnetic moments of $|m|\simeq
1.7$ $\mu_{\mathrm{B}}$ on NiL and it is always found to be higher in energy
than the FM configuration. As shown in Table 1, the breathing mode of the
$S$\- and $T$-type orderings obtained with on-site interactions only (DFT+$U$)
disappears, as the Ni$-$O bond lengths of the two structurally inequivalent
NiO6 octahedra become equal. In addition, as shown in Fig. 2, the electronic
states associated with these orderings display a semimetallic behaviour, with
the $3d$ states of the two Ni atoms becoming equivalent. This state also
exhibits almost the same values of $U$ interaction parameters and magnetic
moments for the two inequivalent Ni sites. For this reason we conclude that
within this _ab initio_ model, at the self-consistent crystal structure and
Hubbard parameters, purely on-site interactions can not sustain the breathing-
mode-based structural distortion.

#### DFT+$U$+$V$

Upon the inclusion of the inter-site Hubbard interactions between the nickels
and the oxygens, the materials recover the experimentally observed insulating
nature. The lowest-energy, semi-degenerate states $(S/T)_{(a)}$ share similar
features for PrNiO3 and YNiO3: they display significantly different magnetic
moments and Ni$-$O bond lengths for the two inequivalent Ni sites. This
difference in $m_{\mathrm{S}}/m_{\mathrm{L}}$ and
$\langle\ell_{\mathrm{S}}/\ell_{\mathrm{L}}\rangle$ is larger for YNiO3 than
for PrNiO3 mirroring the larger disproportionation observed for $R$NiO3 with a
smaller rare-earth ionic radius. Importantly for these states, the numerically
larger (smaller) Hubbard parameters are applied to the larger (smaller) Ni–O
octahedra. The opposite applies to the $B$ and $(S/T)_{(b)}$–types orderings.
This difference turned out to be essential for stabilize $(S/T)_{(a)}$ states
as lower in energy than the FM state. The metastable $(S/T)_{(b)}$ for PrNiO3
states exhibit a substantially smaller
$\mathrm{Ni}_{\mathrm{S}}/\mathrm{Ni}_{\mathrm{L}}$ inequivalence because of a
significant reduction of the breathing mode. They are different from the
corresponding $(S/T)_{(c)}$ states obtained for YNiO3, as in these latter the
breathing mode amplitude shrinks to zero (see Table 1). In addition, both
$(S/T)_{(c)}$ display different electronic occupations in the $e_{g}$ bands
(see Supplementary Information) and different polar distortions with respect
to $(S/T)_{(a)}$ states, exemplified as opposite signs of some components of P
between $S_{(c)}$ and $S_{(a)}$ (see Table 2).

The pDOS in Fig. 2 shows that the materials are negative charge-transfer
insulators, since – according to definition – the $2p$ electronic states of
the oxygens stand between the upper/lower Hubbard bands of the $3d$ states of
the two nickels. The negative character of charge-transfer is due to the
presence of residual O$-{2p}$ states in the lowest conduction bands (the so-
called hole-doped ligands) Bisogni  _et al._ (2016). In all the panels in Fig.
2, the $e_{g}$ states of NiL [$e_{g}(\mathrm{Ni}_{\mathrm{L}})$] present
concentrated spectral weights 5 eV above (minority-spin channel,
$\sigma_{\mathrm{min}}$) and 7 eV below (majority-spin channel,
$\sigma_{\mathrm{maj}}$) the valence-band maximum (VBM). The situation is
similar for $e_{g}(\mathrm{Ni}_{\mathrm{S}})$ in DFT+$U$, as the two nickel
sites become structurally equivalent [see Fig. 2 (b)]. A reduction in the
magnetization on NiS produces a gradual downshift of the
$e_{g}(\mathrm{Ni}_{\mathrm{S}})$ $\sigma_{\mathrm{min}}$ states, until they
finally merge with the $e_{g}(\mathrm{Ni}_{\mathrm{S}})$
$\sigma_{\mathrm{maj}}$ in a unique group of electronic bands centered 2 eV
above the VBM in the limit of vanishing magnetic moments ($B$-type ordering).
The $T_{(a)}$ states of PrNiO3 and YNiO3 exhibit different conduction states:
for PrNiO3 two separate peaks in the pDOS can be distinguished at 1 and 3 eV
above VBM, while in YNiO3 they partially overlap at 2 eV above the same level.
This effect is sharply reflected in the magnetization of NiS, that decreases
substantially from YNiO3 to PrNiO3. We attribute this change in the electronic
structure to the larger structural disproportion for the two Ni that is
present in the YNiO3 compared to PrNiO3. The $T_{(b)}$ and $T_{(c)}$ states
display a strongly reduced band gap. Interestingly in $T_{(c)}$ for YNiO3, the
system still remains insulating despite the vanishing breathing mode, in
contrast with DFT+$U$ where instead the gap closes.

###  III.2 Magnetism-induced electric polarization

Table 2: Electric polarization $\mathbf{P}=\mathbf{P}_{\mathrm{e}}+\mathbf{P}_{\mathrm{I}}$ (electronic + ionic) in cartesian axes (the $x$ direction is aligned along $\hat{\textbf{a}}$) in units of $\mu$C/cm2 within the different AFM phases calculated using the DFT+$U$+$V$ method for PrNiO3 and YNiO3, together with the space group of the relaxed structures. The polarization quanta are $(38.3,19.4,27.2)$ $\mu$C/cm2 for PrNiO3 and $(38.8,21.0,27.6)$ $\mu$C/cm2 for YNiO3. Nonzero components are highlighted in bold. State | $P_{\mathrm{e},x}$ | $P_{\mathrm{e},y}$ | $P_{\mathrm{e},z}$ | Space group (#)
---|---|---|---|---
| $P_{\mathrm{I},x}$ | $P_{\mathrm{I},y}$ | $P_{\mathrm{I},z}$ |
PrNiO3
$B$ | $0.00$ | $0.00$ | $0.00$ | $P2_{1}/c\,(14)$
| $0.00$ | $0.00$ | $0.00$ |
$S_{(a)}$ | 0.00 | 0.21 | $0.00$ | $P2_{1}\,(4)$
| 0.00 | 0.90 | $0.00$ |
$T_{(a)}$ | 4.25 | 0.00 | 2.74 | $Pc\,(7)$
| $-\textbf{0.86}$ | 0.00 | $-\textbf{2.68}$ |
$S_{(b)}$ | 0.00 | 1.72 | 0.00 | $P2_{1}\,(4)$
| 0.00 | 1.58 | 0.00 |
$T_{(b)}$ | 11.7 | 0.00 | 6.14 | $Pc\,(7)$
| $-\textbf{1.15}$ | 0.00 | $-\textbf{4.40}$ |
YNiO3
$B$ | $0.00$ | $0.00$ | $0.00$ | $P2_{1}/c\,(14)$
| $0.00$ | $0.00$ | $0.00$ |
$S_{(a)}$ | 0.00 | 0.07 | 0.00 | $P2_{1}\,(4)$
| 0.00 | 0.77 | 0.00 |
$T_{(a)}$ | 1.27 | 0.00 | 0.68 | $Pc\,(7)$
| 0.13 | 0.00 | 0.99 |
$S_{(c)}$ | 0.00 | $-\textbf{5.57}$ | 0.00 | $P2_{1}\,(4)$
| 0.00 | 3.18 | 0.00 |
$T_{(c)}$ | $-\textbf{38.0}$ | 0.00 | 1.02 | $Pc\,(7)$
| $-\textbf{0.17}$ | 0.00 | 1.70 |

The $P2_{1}/c$ crystal structure of rare-earth nickelates is centrosymmetric.
However, a group theory analysis shows that the magnetic space group of the
$S/T$-type magnetic orderings breaks inversion symmetry, thus allowing
symmetry-lowering _polar_ crystal distortions to occur, resulting in space
groups $Pc$ for $T$-type and $P2_{1}$ for the $S$-type Giovannetti  _et al._
(2009); Ardizzone  _et al._ (2021); Perez-Mato  _et al._ (2016). As described
in detail in Ref. Ardizzone  _et al._ (2021), the loss of centrosymmetry in
the $T$-type ordering is induced by breaking the symmetry of the two-fold
screw-axis parallel to $\mathbf{b}$, while for the $S$-type ordering it is due
to the loss of a glide plane consisting of a reflection in the plane
perpendicular to $\mathbf{b}$ plus a translation along $\mathbf{c}$ of
$|\mathbf{c}|/2$. It has been predicted that the $T$-type ordering would
induce an electric polarization $\mathbf{P}$ with nonzero component along
$\mathbf{a}$ and $\mathbf{c}$, while the $S$-type ordering an electric
polarization predominantly along $\mathbf{b}$ Giovannetti  _et al._ (2009).
The proposed microscopic mechanism consists in a partial charge
disproportionation of the oxygens Oi ($i=1,2,3$) forming the octahedra (thus
creating two inequivalent Oi and $\mathrm{O}_{i}^{\prime}$), which produces a
shift from the Ni-site-centered charge disproportionation to a Ni-O-bond-
centered charge disproportionation, and therefore allowing a nonzero P to
develop. This charge disproportionation of the oxygens is caused by the
presence of neighboring nickels with different magnetic moments: either with
parallel $\mathrm{Ni}_{\sigma}$–$\mathrm{O}_{i}$–$\mathrm{Ni}_{\sigma}$ or
antiparallel
$\mathrm{Ni}_{\sigma}$–$\mathrm{O}_{i}^{\prime}$–$\mathrm{Ni}_{-\sigma}$ spins
$\sigma=\,\uparrow,\downarrow$. Applying the same argument to the $B$-type
ordering, the charge disproportionation of the oxygens can not occur because
every O is surrounded by a magnetic and a nonmagnetic Ni atom and this results
in a zero electric polarization.

We analyze the space groups of the relaxed crystal structures and calculate
from first-principles the electric polarization of the magnetic orders
investigated using DFT+$U$+$V$ method (the $+U$ functional, having vanishing
band gaps for $S$ and $T$ types, would yield an ill-defined polarization). As
shown in Table 2, the relaxed crystal structure of the $B$-type ordering
presents a centrosymmetric $P2_{1}/c$ space group and its electronic ground
state does not exhibit any electric polarization, in agreement with our
previous considerations. On the contrary, the $S/T$-type AFM orderings display
a nonvanishing P, whose magnitude increases with
$m(\mathrm{Ni}_{\mathrm{S}})$. This trend is correlated to the proximity of a
(semi)metallic state in the limit of $m(\mathrm{Ni}_{\mathrm{S}})\rightarrow
m(\mathrm{Ni}_{\mathrm{L}})$. Indeed, an increase of the magnetic moment $m$
is caused by an increase of the difference of the $e_{g}(\mathrm{Ni})$
occupations for the two spin channels, and thus an enhancement of the energy
separation of the associated $e_{g}(\sigma_{\mathrm{min}})$ and
$e_{g}(\sigma_{\mathrm{maj}})$ peaks in the pDOS. Therefore, upon enlarging
$m_{\mathrm{S}}$, the $e_{g}(\mathrm{Ni}_{\mathrm{S}})\,\sigma_{\mathrm{min}}$
states downshift toward the VBM reducing the band gap.

###  III.3 Summary

We studied the structural, electronic and magnetic properties of two members
of the rare-earth nickelates series: PrNiO3 and YNiO3. We used DFT functionals
augmented with on-site $U$ and inter-site $V$ Hubbard corrections, and we
determined the crystal and electronic structures via an iterative self-
consistent procedure alternating structural optimizations and evaluations of
the Hubbard parameters, until convergence between the two is reached.
Calculations were performed without any adjustable parameters, as both Hubbard
$U$ and $V$ were calculated from first-principles using density-functional
perturbation theory; and we showed that the results change qualitatively upon
the inclusion of inter-site Hubbard interactions. We analysed three different
antiferromagnetic orderings: B, S and T-types. When employing purely on-site
electronic interactions (DFT+$U$), the investigated antiferromagnetic patterns
are found to be metastable with respect a ferromagnetic configuration, and
some of these (S-type and T-type orders) yield a semimetallic electronic
structure and lose the crystallographic inequivalence (bond
disproportionation) of the two Ni sites. Within DFT+$U$+$V$, the results for
B-type order are very similar to the ones obtained with DFT+U. On the other
hand, for the $S/T$ antiferromagnetic orderings two different self-consistent
states for each magnetic pattern are reached: two metastable states with
respect the ferromagnetic solution [$(S/T)_{(b)}$ and $(S/T)_{(c)}$ for PrNiO3
and YNiO3, respectively] with very similar magnetic moments on the
inequivalent Ni sites, reduced breathing mode distortion and small band gaps;
and two lowest-energy states with respect to all the magnetic patterns
analysed [$(S/T)_{(a)}$ ] with substantially different Ni–O bond lengths and
magnetic moments for the crystallographic inequivalent nickel atoms, and
larger band gaps. In the iterative loop over the calculation of Hubbard
parameters and structural optimization, these latter states oscillate between
two identical crystal structures with inverted breathing mode amplitudes; we
believe that this behaviour might be due to the neglect of the derivative of
the Hubbard parameters with respect to the atomic positions and strains Kulik
and Marzari (2011), and in general for the lack of a functional formulation of
the method involving the stationarity with respect to the Hubbard parameters
themselves. Finally, without resorting to any empirical parameter, within the
DFT+U+V approach we naturally predict the occurrence of a magnetically-induced
electric polarization, thus supporting recent experimental observations
pointing to the emergence of a multiferroic phase in this class of materials.

##  IV methods

In our first-principles calculations we employ both the rotationally-invariant
DFT+$U$ formulation of Dudarev _et al._ Dudarev  _et al._ (1998), and the
extended DFT+$U$+$V$ one Campo and Cococcioni (2010). The parameters $U$ and
$V$ are calculated from first principles using DFPT Baroni  _et al._ (2001);
Timrov  _et al._ (2018, 2021) as implemented in the HP code Timrov  _et al._
(2022). The localized atomic-like set of functions defining the occupation
matrices (defined in Ref. Cococcioni and de Gironcoli (2005)) employed in our
calculations consists of atomic orbitals as found in the pseudopotentials
orthogonalized through the Löwdin algorithm Löwdin (1950). For the calculation
of the Hubbard contribution to forces and stresses we employ a recently-
developed method that uses the solution of the Lyapunov equation to evaluate
the derivative of the inverse square root of the overlap matrix Timrov  _et
al._ (2020).

All the first-principles calculations are carried out using the Quantum
ESPRESSO distribution Giannozzi  _et al._ (2009, 2017, 2020) . We use a spin-
polarized generalized-gradient approximation (GGA) with the PBEsol
prescription for the exchange-correlation functional Perdew  _et al._ (2008),
as well as a PAW pseudopotential (PP) Blöchl (1994) for O from the Pslibrary
v0.3.1 Dal Corso (2014) and ultrasoft (US) PPs Vanderbilt (1990) from the GBRV
(v1.2) and GBRV (v1.4) libraries Garrity  _et al._ (2014) for Y and Ni,
respectively, as recommended from the SSSP (efficiency) library Prandini  _et
al._ (2018). For the Pr we used the US PP from the Pslibrary1.0.0 (Pr.pbesol-
spdn-rrkjus_psl.1.0.0.UPF) having the $4f$ states frozen in the core. In all
our calculations we use 50 Ry as the plane-wave kinetic-energy cutoff for the
wavefunctions, 400 Ry kinetic-energy cutoff for the charge density, Gaussian
smearing of 0.005 Ry to converge narrow-gap magnetic states and a uniform
$\Gamma$-centered 5$\times$8$\times$6 grid of $\mathbf{k}$ points to sample
the Brillouin zone. Calculations of the electric polarization are performed
using the Berry-phase approach King-Smith and Vanderbilt (1993); Umari and
Pasquarello (2002) and using a finer 7$\times$12$\times$8 k mesh. The
calculations of the Hubbard parameters are performed using 1$\times$2$\times$2
q points, which corresponds to a supercell of an equivalent size Timrov  _et
al._ (2018). We adopt the iterative procedure detailed in Cococcioni and
Marzari (2019); Timrov  _et al._ (2021) to incorporate in the parameters $U$
and $V$ the dependence of both the electronic structure and the crystal
environment: this consists in a series of alternating self-consistent
structural optimizations and linear-response evaluations of the Hubbard
parameters until both the crystal structure and the $U$ and $V$ parameters are
converged self-consistently.

##  V data availability statement

The data used to produce the results of this paper are available on the
Materials Cloud Archive Binci  _et al._ .

##  VI acknowledgements

We gratefully acknowledge Marisa Medarde for fruitful discussions. This
research was supported by the NCCR MARVEL, a National Centre of Competence in
Research, funded by the Swiss National Science Foundation (grant number
182892). Computer time was provided by CSCS (Piz Daint) through project No.
s1073.

##  VII competing interests

The authors declare no competing interests.

##  VIII author contributions

L.B. performed and analyzed the first-principles calculations, M.K., I.T. and
N.M. helped to interpret and understand the data, while N.M. conceived
project. L.B. wrote the initial manuscript and all the authors participated to
the final version of it.

## References

  * Medarde (1997) M. L. Medarde, Journal of Physics: Condensed Matter 9, 1679 (1997).
  * Alonso  _et al._ (1999) J. A. Alonso, J. L. García-Muñoz, M. T. Fernández-Díaz, M. A. G. Aranda, M. J. Martínez-Lope, and M. T. Casais, Phys. Rev. Lett. 82, 3871 (1999).
  * Alonso  _et al._ (2000) J. A. Alonso, M. J. Martínez-Lope, M. T. Casais, J. L. García-Muñoz, and M. T. Fernández-Díaz, Phys. Rev. B 61, 1756 (2000).
  * Middey  _et al._ (2016) S. Middey, J. Chakhalian, P. Mahadevan, J. Freeland, A. Millis, and D. Sarma, Annual Review of Materials Research 46, 305 (2016).
  * Hepting  _et al._ (2018) M. Hepting, R. J. Green, Z. Zhong, M. Bluschke, Y. E. Suyolcu, S. Macke, A. Frano, S. Catalano, M. Gibert, R. Sutarto, F. He, G. Cristiani, G. Logvenov, Y. Wang, P. A. van Aken, P. Hansmann, M. L. Tacon, J.-M. Triscone, G. A. Sawatzky, B. Keimer, and E. Benckiser, Nature Physics 14, 1097 (2018).
  * Scherwitzl  _et al._ (2010) R. Scherwitzl, P. Zubko, I. G. Lezama, S. Ono, A. F. Morpurgo, G. Catalan, and J.-M. Triscone, Advanced Materials 22, 5517 (2010).
  * Post  _et al._ (2018) K. W. Post, A. S. McLeod, M. Hepting, M. Bluschke, Y. Wang, G. Cristiani, G. Logvenov, A. Charnukha, G. X. Ni, P. Radhakrishnan, M. Minola, A. Pasupathy, A. V. Boris, E. Benckiser, K. A. Dahmen, E. W. Carlson, B. Keimer, and D. N. Basov, Nature Physics 14, 1056 (2018).
  * Liu  _et al._ (2013) J. Liu, M. Kargarian, M. Kareev, B. Gray, P. J. Ryan, A. Cruz, N. Tahir, Y.-D. Chuang, J. Guo, J. M. Rondinelli, J. W. Freeland, G. A. Fiete, and J. Chakhalian, Nature Communications 4 (2013), 10.1038/ncomms3714.
  * Li  _et al._ (2019) D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature 572, 624 (2019).
  * Osada  _et al._ (2020) M. Osada, B. Y. Wang, B. H. Goodge, K. Lee, H. Yoon, K. Sakuma, D. Li, M. Miura, L. F. Kourkoutis, and H. Y. Hwang, Nano Letters 20, 5735 (2020).
  * Ardizzone  _et al._ (2021) I. Ardizzone, J. Teyssier, I. Crassee, A. B. Kuzmenko, D. G. Mazzone, D. J. Gawryluk, M. Medarde, and D. van der Marel, Phys. Rev. Research 3, 033007 (2021).
  * Momma and Izumi (2008) K. Momma and F. Izumi, J. Appl. Crystallogr. 41, 653 (2008).
  * Mazin  _et al._ (2007) I. I. Mazin, D. I. Khomskii, R. Lengsdorf, J. A. Alonso, W. G. Marshall, R. M. Ibberson, A. Podlesnyak, M. J. Martínez-Lope, and M. M. Abd-Elmeguid, Phys. Rev. Lett. 98, 176406 (2007).
  * Medarde  _et al._ (2008) M. Medarde, M. T. Fernández-Díaz, and P. Lacorre, Phys. Rev. B 78, 212101 (2008).
  * Park  _et al._ (2012) H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. Lett. 109, 156402 (2012).
  * Raebiger  _et al._ (2008) H. Raebiger, S. Lany, and A. Zunger, Nature 453, 763 (2008).
  * Zaanen  _et al._ (1985) J. Zaanen, G. A. Sawatzky, and J. W. Allen, 55, 418 (1985).
  * Johnston  _et al._ (2014) S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 112, 106404 (2014).
  * Haule and Pascut (2017) K. Haule and G. L. Pascut, Scientific Reports 7 (2017), 10.1038/s41598-017-10374-2.
  * Varignon  _et al._ (2017) J. Varignon, M. N. Grisolia, J. Íñiguez, A. Barthélémy, and M. Bibes, npj Quantum Materials 2 (2017), 10.1038/s41535-017-0024-9.
  * Hampel and Ederer (2017) A. Hampel and C. Ederer, Phys. Rev. B 96, 165130 (2017).
  * Mizokawa  _et al._ (2000) T. Mizokawa, D. I. Khomskii, and G. A. Sawatzky, Phys. Rev. B 61, 11263 (2000).
  * Green  _et al._ (2016) R. J. Green, M. W. Haverkort, and G. A. Sawatzky, Phys. Rev. B 94, 195127 (2016).
  * García-Muñoz  _et al._ (1994) J. L. García-Muñoz, J. Rodríguez-Carvajal, and P. Lacorre, Phys. Rev. B 50, 978 (1994).
  * Mercy  _et al._ (2017) A. Mercy, J. Bieder, J. Íñiguez, and P. Ghosez, Nature Communications 8 (2017), 10.1038/s41467-017-01811-x.
  * Badrtdinov  _et al._ (2021) D. I. Badrtdinov, A. Hampel, and C. E. Dreyer, Phys. Rev. B 104, 054403 (2021).
  * Giovannetti  _et al._ (2009) G. Giovannetti, S. Kumar, D. Khomskii, S. Picozzi, and J. van den Brink, Phys. Rev. Lett. 103, 156401 (2009).
  * Perez-Mato  _et al._ (2016) J. M. Perez-Mato, S. V. Gallego, L. Elcoro, E. Tasci, and M. I. Aroyo,  28, 286001 (2016).
  * Campo and Cococcioni (2010) V. L. Campo and M. Cococcioni, Journal of Physics: Condensed Matter 22, 055602 (2010).
  * Timrov  _et al._ (2018) I. Timrov, N. Marzari, and M. Cococcioni, Phys. Rev. B 98, 085127 (2018).
  * Timrov  _et al._ (2021) I. Timrov, N. Marzari, and M. Cococcioni, Phys. Rev. B 103, 045141 (2021).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * Dudarev  _et al._ (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * (34) L. Binci, M. Kotiuga, I. Timrov, and N. Marzari, Materials Cloud Archive (2022), to be published, available upon request.
  * Gawryluk  _et al._ (2019) D. J. Gawryluk, Y. M. Klein, T. Shang, D. Sheptyakov, L. Keller, N. Casati, P. Lacorre, M. T. Fernández-Díaz, J. Rodríguez-Carvajal, and M. Medarde, Phys. Rev. B 100, 205137 (2019).
  * Alonso  _et al._ (2001) J. A. Alonso, M. J. Martínez-Lope, M. T. Casais, J. L. García-Muñoz, M. T. Fernández-Díaz, and M. A. G. Aranda, Phys. Rev. B 64, 094102 (2001).
  * (37) Supplementary information describing the spin-resolved projected density of states, convergence of Hubbard parameters and total energy differences, and band gaps.
  * Kulik and Marzari (2011) H. J. Kulik and N. Marzari, The Journal of Chemical Physics 135, 194105 (2011).
  * Cococcioni and de Gironcoli (2005) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
  * Cococcioni and Marzari (2019) M. Cococcioni and N. Marzari, Phys. Rev. Materials 3, 033801 (2019).
  * Bengone  _et al._ (2000) O. Bengone, M. Alouani, P. Blöchl, and J. Hugel, Phys. Rev. B 62, 16392 (2000).
  * Wang  _et al._ (2016) Y.-C. Wang, Z.-H. Chen, and H. Jiang, J. Chem. Phys. 144, 144106 (2016).
  * King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  * Umari and Pasquarello (2002) P. Umari and A. Pasquarello, Phys. Rev. Lett. 89, 157602 (2002).
  * Bisogni  _et al._ (2016) V. Bisogni, S. Catalano, R. J. Green, M. Gibert, R. Scherwitzl, Y. Huang, V. N. Strocov, P. Zubko, S. Balandeh, J.-M. Triscone, G. Sawatzky, and T. Schmitt, Nature Communications 7 (2016), 10.1038/ncomms13017.
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Timrov  _et al._ (2022) I. Timrov, N. Marzari, and M. Cococcioni, Comput. Phys. Commun. 279, 108455 (2022).
  * Löwdin (1950) P.-O. Löwdin, J. Chem. Phys. 18, 365 (1950).
  * Timrov  _et al._ (2020) I. Timrov, F. Aquilante, L. Binci, M. Cococcioni, and N. Marzari, Phys. Rev. B 102, 235159 (2020).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ (2020) P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, J. Chem. Phys. 152, 154105 (2020).
  * Perdew  _et al._ (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Dal Corso (2014) A. Dal Corso, Computational Materials Science 95, 337 (2014).
  * Vanderbilt (1990) D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
  * Garrity  _et al._ (2014) K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Computational Materials Science 81, 446 (2014).
  * Prandini  _et al._ (2018) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, npj Computational Materials 4 (2018).


# Deep Universal Blind Image Denoising

Jae Woong Soh  Department of ECE, INMC
Seoul National University
Seoul, Korea
Email<EMAIL_ADDRESS>Nam Ik Cho  Department of ECE, INMC
Seoul National University
Seoul, Korea
Email<EMAIL_ADDRESS>
###### Abstract

Image denoising is an essential part of many image processing and computer
vision tasks due to inevitable noise corruption during image acquisition.
Traditionally, many researchers have investigated image priors for the
denoising, within the Bayesian perspective based on image properties and
statistics. Recently, deep convolutional neural networks (CNNs) have shown
great success in image denoising by incorporating large-scale synthetic
datasets. However, they both have pros and cons. While the deep CNNs are
powerful for removing the noise with known statistics, they tend to lack
flexibility and practicality for the blind and real-world noise. Moreover,
they cannot easily employ explicit priors. On the other hand, traditional non-
learning methods can involve explicit image priors, but they require
considerable computation time and cannot exploit large-scale external
datasets. In this paper, we present a CNN-based method that leverages the
advantages of both methods based on the Bayesian perspective. Concretely, we
divide the blind image denoising problem into sub-problems and conquer each
inference problem separately. As the CNN is a powerful tool for inference, our
method is rooted in CNNs and propose a novel design of network for efficient
inference. With our proposed method, we can successfully remove blind and
real-world noise, with a moderate number of parameters of universal CNN.

##  I Introduction

Image denoising aims to recover the latent clean image from an observed noisy
image. It has been a longstanding fundamental problem because noise
intervention is inevitable during the image acquisition process, which
degrades the visual quality of the acquired image and lowers the performance
of computer vision tasks. The overall noise is the accumulation of multiple
different noise sources, such as capturing sensors, in-camera pipeline, and
data transmission media. Such a noise generation process is too complicated to
address, therefore with the belief of the central limit theorem, the noise is
usually assumed as additive white Gaussian noise (AWGN).

Specifically, an observed noisy image $\mathbf{y}$ has been assumed to be the
sum of the latent clean image $\mathbf{x}$ and an AWGN $\mathbf{n}$ as

${\mathbf{y}}={\mathbf{x}}+{\mathbf{n}}.$ (1)

Traditionally, image denoising has been addressed as a statistical inference
problem, which focused on constructing a maximum a posteriori (MAP) inference
model. In the Bayesian perspective, the MAP inference can be divided into the
data fidelity term (log-likelihood) and regularization term (log-prior) as

$\displaystyle\mathbf{\hat{x}}$ $\displaystyle=\arg\max_{{\mathbf{x}}}\log
p({\mathbf{x}}|{\mathbf{y}}),$ (2) $\displaystyle=\arg\max_{{\mathbf{x}}}\log
p({\mathbf{y}}|{\mathbf{x}})+\log p({\mathbf{x}}).$ (3)

To obtain a plausible model in this form and its solution, many classical
researches focused on the design of image prior terms based on experience and
knowledge of the data. These approaches include total variation regularization
[1], sparsity [2], low-rank constraints [3, 4], non-local self-similarity [5,
6, 7, 8], and so on. These methods _explicitly_ defined the data fidelity term
with the assumption of AWGN and appropriate image prior term. They are shown
to be superior in terms of interpretability and flexibility, but they are
limited when the noise deviates from spatially uniform i.i.d. Gaussian noise.
Also, their optimization implementation mostly depends on iterative solvers,
which takes a long computation time.

Recently, the upsurge of powerful deep learning in computer vision introduced
deep learning-based denoisers, which _implicitly_ learns the MAP estimate
based on the supervised learning framework [9, 10, 11, 12, 13, 14, 15, 16,
17]. These methods surpass the conventional non-learning methods by large
margins, but they lack in flexibility, especially for blind denoising. For
example, if a network is trained to denoise a specific level of noise
variance, the network would not work well for an unseen noise-level. Hence, a
universal network is needed, which is trained with the whole expected range of
noise-level for blind denoising unless the noise-level is known. In this case,
it is still inferior to the specific noise-level denoiser, and discriminative
learning methods suffer from the averaging problem, which means that the
network targets better to the mid-level noise [9, 10]. Moreover, knowledge
from human experience cannot be easily injected into the deep learning
framework. Although there are some works to merge non-local self-similarity,
one of the strong priors, with deep networks [12, 15, 16], introducing other
priors is not easy.

In this paper, we propose a convolutional neural network (CNN)-based universal
blind denoiser, dubbed Deep Universal Blind Denoiser (DUBD), which leverages
the advantages of MAP inference and the power of deep-learning. The design of
DUBD is motivated by the advantages and disadvantages of the above-stated
approaches, how to insert human knowledge priors to the deep-learning
framework without harming the power of the network. In particular, for
practicality, the proposed DUBD can handle a wide range of noise-level and
spatially varying noises without knowing the noise-level. Also, our method can
address spectrally variant noise, as evidenced in [4]. Moreover, our DUBD
outperforms other methods, including blind and non-blind denoisers, while
requiring a comparable number of parameters. Our contributions can be
summarized as:

  * •

We propose a CNN-based universal blind denoiser that can handle a wide range
of noise-level including spatially and spectrally varying noise.

  * •

Our DUBD can explicitly incorporate prior knowledge, which further lifts the
performance of the network.

  * •

Our DUBD outperforms other denoisers with a comparable number of parameters,
which eventually brings better practicality.

  * •

Our DUBD can be applied to real-world noisy images and also shows outstanding
performance compared to the other methods.

##  II Probabilistic View

For blind image denoising, a naïve approach is to train a universal CNN with
the whole expected range of noise-level. However, it is proven not to be a
good choice, lacking performance and flexibility [9, 10]. Hence, we approach
the blind image denoising problem with the divide-and-conquer scheme.
Precisely, we reformulate the log-posterior by introducing a new random
variable ${\mathbf{c}}$ that contains prior based on human knowledge as

$\displaystyle\log p({\mathbf{x}}|{\mathbf{y}})$
$\displaystyle=\log\int_{{\mathbf{c}}}p({\mathbf{x}}|{\mathbf{y}},{\mathbf{c}})p({\mathbf{c}}|{\mathbf{y}})d{\mathbf{c}},$
(4) $\displaystyle\approx\log
p({\mathbf{x}}|{\mathbf{y}},\mathbf{\hat{c}})p(\mathbf{\hat{c}}|{\mathbf{y}}),$
(5) $\displaystyle=\log p({\mathbf{x}}|{\mathbf{y}},\mathbf{\hat{c}})+\log
p(\mathbf{\hat{c}|y}).$ (6)

Note that Eq. 4 is a marginal expression by introducing a random variable
${\mathbf{c}}$, but the integration is mostly not tractable. Hence, for
approximating the integration, we use the point estimate for ${\mathbf{c}}$,
which is
$\mathbf{\hat{c}}=\arg\max_{{\mathbf{c}}}p({\mathbf{c}}|{\mathbf{y}})$. When
$p({\mathbf{c}}|{\mathbf{y}})$ has a unimodal distribution with a sharp peak,
the approximation is quite fair. Then, we can solve the MAP estimate with the
given point estimate $\mathbf{\hat{c}}$. Formally, the MAP inference problem
can be reformulated into solving two sub-problems as

$\displaystyle\mathbf{\hat{c}}$ $\displaystyle=\arg\max_{{\mathbf{c}}}\log
p({\mathbf{c}}|{\mathbf{y}}),$ (7) $\displaystyle\mathbf{\hat{x}}$
$\displaystyle=\arg\max_{{\mathbf{x}}}\log
p({\mathbf{x}}|{\mathbf{y}},\mathbf{\hat{c}}).$ (8)

As neural networks are experts at inference, we employ CNNs for our two
inference problems, described as

$\displaystyle\mathbf{\hat{c}}$ $\displaystyle\approx
g_{\theta}({\mathbf{y}}),$ (9) $\displaystyle\mathbf{\hat{x}}$
$\displaystyle\approx f_{\phi}({\mathbf{y}};\mathbf{\hat{c}}),$ (10)

where $g_{\theta}(\cdot)$ and $f_{\phi}(\cdot)$ are CNNs with parameters
$\theta$ and $\phi$, respectively. By dividing the original problem into sub-
problems, we decrease the modalities of multi-modal distribution of the log-
posterior. Thus the network conquers separate conditionals and eventually
facilitates solving the whole blind denoising inference. Concretely, for
solving the second inference problem, we design a tunable CNN according to
${\mathbf{c}}$. Importantly, we must introduce an appropriate ${\mathbf{c}}$
based on our prior knowledge.

##  III Conditional Estimation Network

Figure 1: The network architecture of conditional estimation network (CENet).
We adopt $3\times 3$ convolution layers for all layers, and the number of
channels is $64$ except for the last one. We decrease the spatial size of
feature maps with the strided convolution layers. At the last stage, we
bilinearly interpolate the map to recover its spatial size back to the same
size of the input image. Figure 2: The network architecture of tunable
denoiser. We adopt $3\times 3$ convolution layers for the main stem of the
network, including conditional affine transform blocks (CATBlocks) and
residual blocks (ResBlocks). Also, the number of channels is $64$ for the main
stem, and $1\times 1$ convolution layers are employed for the condition
encoder. In condition encoder, the number of channels for the first layer is
$128$, and the others are $64$.

In this section, we introduce the choice of ${\mathbf{c}}$, and how we design
a network for the point estimate $\mathbf{\hat{c}}$. Notably, we consider the
cases that $p({\mathbf{c}}|{\mathbf{y}})$ is unimodal with a sharp peak, for
tractable solution, as ${\mathbf{c}}$ may belong to any category unless it
satisfies such properties. Fortunately, a noise-level of a given noisy image
is deterministic, which means that its distribution is unimodal with a
relatively sharp peak. In other words, only a single noise-level corresponds
to the given noisy image. Therefore, we choose the noise-level as
${\mathbf{c}}$. In addition to the blind scenario, we can also handle the non-
blind case by solving the sub-problem of Eq. 8 with known information
${\mathbf{c}}$. Also, we can manually control ${\mathbf{c}}$ to manipulate the
network to produce better outputs.

To estimate the noise-level of a given image with an unknown and spatially
varying noise variance, we exploit a CNN architecture shown in Figure 1. In
addition to the choice of ${\mathbf{c}}$, we add another prior knowledge that
the variation of noise-level is spatially smooth in designing the network. We
decrease the size of feature maps using the strided convolutions and apply
bilinear interpolation at the last stage, rather than using a smoothness term
such as total variation regularization. By doing so, we have two advantages
compared to the additional loss term. First, we can decrease the computational
complexity by decreasing the feature map size, and second, we do not need to
tune the hyper-parameter that balances the original loss and the
regularization term. In summary, the network is trained with the loss:

$\mathcal{L}_{c}(\theta)=\mathbb{E}[||\sigma-
g_{\theta}({\mathbf{y}})||^{2}_{2}].$ (11)

The output of $g_{\theta}({\mathbf{y}})$ is the estimated noise-level map
$\hat{\sigma}\in\mathbb{R}^{H\times W\times C}$, where $H$, $W$, and $C$
denote height, width and number of channels of the noisy image respectively.

##  IV Tunable Denoising Network

In this section, we present a novel network architecture where the input is a
noisy image with the intervention of an additional tunable parameter
${\mathbf{c}}$ to emulate $\arg\max
p({\mathbf{x}}|{\mathbf{y}},\mathbf{\hat{c}})$. The overall architecture of
our tunable denoiser is shown in Figure 2 where its design directly reflects
$f_{\phi}(y;\hat{c})$. The denoising network mainly consists of two parts: the
main stem and condition encoder.

###  IV-A Main Stem

The main stem of our tunable denoising network is based on the cascade of $D$
conditional affine transform blocks (CATBlocks). To incorporate conditional
information efficiently, we present a CATBlock which applies the affine
transformation to feature maps as

$F_{o}=\gamma\otimes F_{i}+\beta,$ (12)

where $F_{i}$ and $F_{o}$ denote feature maps before and after the
transformation, respectively, and $\otimes$ denotes element-wise
multiplication. In the CATBlock, the output feature maps of the last
convolution layer are adjusted by affine transformation with respect to the
condition-encoded affine parameters $\gamma$ and $\beta$, and the encoded
parameters are shared along with other CATBlocks. Similar approaches have been
applied in various tasks [18, 19, 20, 21].

The CATBlock includes a cascade of $N$ residual blocks (ResBlocks), which is
shown in the upper-right corner of Figure 2. Additionally, we adopt the
residual learning scheme [9], which is to learn noise rather than the clean
image itself. Notably, our CATBlocks and ResBlocks adopt residual skip-
connection to bypass short- to long-term information.

###  IV-B Condition Encoder

The condition encoder network is a simple two-layers of $1\times 1$
convolution. It takes the conditional variable and outputs the condition-
encoded parameters $\gamma$ and $\beta$. These parameters adjust the feature
values according to the condition. It is notable that the number of parameters
of the condition encoder is negligible compared to the whole network.

###  IV-C Implementation Details

Now we specify the implementation details of our DUBD. For the number of
blocks, we set $D=5$ and $N=5$. The DIV2K [22] dataset is used for the
training, which is a recently released high-resolution dataset. It consists of
high-quality images of $800$ training, $100$ validation, and $100$ test
images. The noise-levels of Gaussian noise are randomly sampled from the range
of $[5,70]$ uniformly. We extract $96\times 96$ size of patches for training.
We adopt the ADAM optimizer [23], and the initial learning rate is set to
$2\times 10^{-4}$ and halved once during training. Mean squared error (MSE)
loss function, which is expressed as

$\mathcal{L}_{dn}(\phi)=\mathbb{E}[||{\mathbf{x}}-f_{\phi}({\mathbf{y}};{\mathbf{c}})||^{2}_{2}]$
(13)

is adopted for denoising loss.

TABLE I: The average PSNR results on test sets. The best results are denoted
in red and the second best in blue. ‘*’ next to the name of the method means
blind denoising.

Noise-level | Dataset | CBM3D [8] | TNRD [24] | RED [11] | MemNet* [25] | CDnCNN* [9] | FFDNet [10] | UNLNet* [13] | ATDNet* [14] | VDN* [17] | DUBD-NB (Ours) | DUBD-B* (Ours)
---|---|---|---|---|---|---|---|---|---|---|---|---
$\sigma=10$ | CBSD68 | 35.91 | - | 33.89 | 28.52 | 36.13 | 36.14 | 36.20 | 36.29 | 36.29 | 36.35 | 36.33
| Kodak24 | 36.43 | - | 34.73 | 29.70 | 36.46 | 36.69 | - | 36.98 | 36.85 | 37.03 | 37.02
| Urban100 | 36.00 | - | 34.42 | 29.44 | 34.61 | 35.78 | - | 36.31 | 35.97 | 36.32 | 36.23
$\sigma=30$ | CBSD68 | 29.73 | - | 28.45 | 28.39 | 30.34 | 30.32 | 30.21 | 30.61 | 30.64 | 30.65 | 30.62
| Kodak24 | 30.75 | - | 29.53 | 29.55 | 31.17 | 31.27 | 31.18 | 31.72 | 31.67 | 31.75 | 31.75
| Urban100 | 30.36 | - | 28.84 | 28.93 | 30.00 | 30.53 | 30.41 | 31.48 | 31.14 | 31.46 | 31.43
$\sigma=50$ | CBSD68 | 27.38 | 25.96 | 26.34 | 26.33 | 27.95 | 27.97 | 27.85 | 28.33 | 28.33 | 28.35 | 28.31
| Kodak24 | 28.46 | 27.04 | 27.42 | 27.51 | 28.83 | 28.98 | 28.86 | 29.48 | 29.44 | 29.51 | 29.50
| Urban100 | 27.94 | 25.52 | 26.25 | 26.53 | 27.59 | 28.05 | 27.95 | 29.20 | 28.86 | 29.16 | 29.14
$\sigma=70$ | CBSD68 | 26.00 | - | 25.09 | 25.09 | 25.66 | 26.55 | - | - | 26.93 | 26.96 | 26.89
| Kodak24 | 27.09 | - | 26.16 | 26.24 | 26.36 | 27.56 | - | - | 28.05 | 28.12 | 28.11
| Urban100 | 26.31 | - | 24.58 | 24.93 | 25.24 | 26.40 | - | - | 27.31 | 27.59 | 27.58

TABLE II: The number of parameters and PSNR for several CNN-based methods. The
PSNR is measured on Urban100 [26] set with $\sigma=50$. For RNAN [16], we
brought the number from the paper.

Methods | RED [11] | CDnCNN* [9] | FFDNet [10] | ATDNet* [14] | RNAN [16] | VDN* [17] | DUBD-NB (Ours) | DUBD-B* (Ours)
---|---|---|---|---|---|---|---|---
Parameters | 4,135 K | 668 K | 825 K | 9,453 K | 7,409 K | 7,817 K | 2,088 K | 2,239 K
PSNR (dB) | 26.25 | 27.59 | 28.05 | 29.20 | 29.08 | 28.86 | 29.16 | 29.14

Ground Truth

Noisy

CBM3D [8]

TNRD [24]

RED [11]

MemNet [25]

CDnCNN [9]

FFDNet [10]

UNLNet [13]

ATDNet [14]

DUBD-NB (Ours)

DUBD-B (Ours)

Figure 3: Visualized examples of denoising for $\sigma=50$.

Ground Truth

Noisy

CBM3D [8]

TNRD [24]

RED [11]

MemNet [25]

CDnCNN [9]

FFDNet [10]

UNLNet [13]

ATDNet [14]

DUBD-NB (Ours)

DUBD-B (Ours)

Figure 4: Visualized examples of denoising for $\sigma=50$.

##  V Experimental Results

We test denoising performances on the images corrupted with noise-levels
$\sigma=10,~{}30,~{}50,~{}70$, on three famous denoising datasets for color
images: CBSD68 [27], Kodak24, and Urban100 [26]. Importantly, considering the
practical usage, we mainly aim color image denoising where gray image
denoising can also be done with the change of input channels.

We compare our method with several denoising algorithms from conventional non-
learning methods to recent CNN-based state-of-the-arts: CBM3D [8], TNRD [24],
RED [11], MemNet [25], DnCNN [9], FFDNet [10], UNLNet [13], and ATDNet [14].
We present two results of our DUBD, namely DUBD-NB and DUBD-B. The DUBD-NB is
a non-blind model where we assume known noise-levels, and thus the output
corresponds to $f_{\phi}({\mathbf{y}};{\mathbf{c}})$. On the other hand,
DUBD-B is a blind model where the denoiser takes the estimated ${\mathbf{c}}$
from CENet, and hence the output is $f_{\phi}({\mathbf{y}};\mathbf{\hat{c}})$.
It is notable that we train only _one_ universal network for our DUBD that can
be used both as blind and non-blind, whereas other methods include blind and
non-blind models separately for each purpose. We evaluate the results by PSNR
(peak signal-to-noise ratio) on color images.

The overall results are listed in TABLE I, which shows that our methods are on
par with or exceeding previous state-of-the-art methods. Notably, our DUBD-B
shows much better performances than other non-blind ones. Interestingly,
DUBD-B and DUBD-NB show negligible performance gaps, even though the
estimation is not quite perfect. In other words, our method is robust or not
very sensitive to the conditional estimations. ATDNet [14] blind model shows
comparable results to ours or sometimes shows the best results among the
others. But the applicable range of noise-level of ATDNet is narrower than
ours. Also, an interesting fact is that CBM3D [8], which is a conventional
approach exploiting non-local self-similarity with 3D transformed Wiener
filter, shows quite great results better than most of the CNN-based methods in
Urban100 [26] set that have many recurrent patterns.

For visual comparisons, we show some denoising results in Figure 3 and 4. As
shown, ours show a better separation of content and noise, _i.e_ , noise
reduction. For other methods, as shown in Figure 4, they often show cloudy
patterns or entangled noise, and sometimes show color artifacts. On the other
hand, our DUBDs achieve plausible detail preservation and naturalness.

##  VI Analysis

We further analyze our DUBD with other recent state-of-the-art methods,
focusing on the complexity and robustness to spatially or spectrally varying
noises, which are important factors for practicability.

Figure 5: The number of parameter vs. PSNR on CBSD68 with $\sigma=50$.

###  VI-A The Number of Parameters

Since the computational resource is limited in a real-world situation, the
number of parameters is an important factor for the practicability of CNNs.
TABLE II and Figure 5 show the summary of network complexity vs. performance
for the methods compared earlier. Compared to CDnCNN [9] and FFDNet [10], our
methods show PSNR gain more than $1.0$ to $1.5$ dB by doubling or tripling the
number of parameters. Our DUBD outperforms RED [11] with about half the number
of parameters. Also, compared to ATDNet [14] and RNAN [16], our methods need
much fewer parameters, while showing comparable PSNRs. Importantly, RNAN [16]
is a non-blind method which requires each separate model for a specific noise-
level. Additionally, in our methods, changing from non-blind method to blind
method requires only 151 K parameters for CENet. We also visualize the
comparisons of several methods in terms of parameters versus performance in
Figure 5. It is noteworthy that although we both include blind and non-blind
methods in the comparisons, our DUBD shows outstanding performance with a
comparable number of parameters. In summary, our DUBD is practical and
competitive compared to the state-of-the-art methods.

###  VI-B Dealing with Spectrally-Spatially Variant Noise

In a real-world situation, the noise in an image is rarely spatially uniform.
Sometimes, the noise-level varies in different color channels, as pointed out
in [10, 4]. Thus, in this section, we show the results of spectrally-spatially
varying noise with our DUBD-B.

####  VI-B1 Spectrally Variant Noise

To test on spectrally variant noise, we feed an input as shown in Figure 6(d),
where the extent of noise corruption differs in each color channel. As shown
in Figure 6(a-c), the noise corruption is severe in B-channel while it is mild
in R-channel. As shown in Figure 6(e), addressing this type of noise with only
a single averaged noise-level does not show plausible results. It can be seen
that the noise in B-channel (appears in the bluish region of the image)
remains after the denoising, while there is loss-of-details of the red bird.
On the other hand, since our method can handle each channel separately, it
shows visually promising results as in Figure 6(f).

(a) Noise in R

(b) Noise in G

(c) Noise in B

(d) Input noisy image

(e) Output with average noise-level

(f) Our result

Figure 6: Input and results on spectrally varying noise.

####  VI-B2 Spatially Variant Noise

We feed an input image with spatially variant noise as shown in Figure 7(b).
The noise-level of each spatial position is shown in Figure 7(a) and the
results are shown in Figure 7(c) and (d). As shown in the results, denoising
with a single noise-level shows over-smoothed regions, as in the red bird, and
there are some regions that the noise is not removed yet. On the other hand,
our method handles spatially variant noise well within the blind denoising
scheme, showing visually pleasing results.

(a) Noise-level in spatial dimension

(b) Input noisy image

(c) Output with average noise-level

(d) Our result

Figure 7: Input and results on spatially varying noise.

###  VI-C Traversing Conditional Variable

(a) PSNR results on different input ${\mathbf{c}}$.

(b) Visualized results of denoised images with different input ${\mathbf{c}}$.

Figure 8: Results by traversing the conditional variable ${\mathbf{c}}$.

In this section, we show the results by traversing the conditional variable
${\mathbf{c}}$ for two reasons. First, we show that our tunable denoiser can
also be manually controlled in accordance with user preferences. By
controlling the ${\mathbf{c}}$, the denoising strength can be tuned. Second,
we attempt to investigate the effect of conditional variable and indirectly
trace the conditional posterior $p({\mathbf{c}}|{\mathbf{y}})$.

The results are shown in Figure 8, where we feed an image corrupted by i.i.d.
Gaussian noise with $\sigma=30$ by changing the input ${\mathbf{c}}$ of our
tunable denoiser. As shown in the Figure 8(a), the PSNR graph achieves its
peak at ${\mathbf{c}}=31$, near the true $\sigma$. Also, our blind result with
conditional estimation shows quite a competent result. Notably, the graph
shape is asymmetric unimodal, where higher ${\mathbf{c}}$ than the true value
shows much less sensitivity then feeding the lower ${\mathbf{c}}$ than the
true value. Its coherent tendency can also be found in the corresponding
visualizations on Figure 8(b). Starting from small ${\mathbf{c}}$, the noise
remains in the image, but by tuning ${\mathbf{c}}$ to higher values, the noise
gradually diminishes. Then, increasing ${\mathbf{c}}$ higher than the true
noise-level, the image gets blurry, losing details and textures.

##  VII Results on Real Noisy Images

As the noise from the real-world largely deviates from the AWGN, it is
important to validate whether the network will also work well for the real-
world noises, with possible modification or retraining. Real-noise originates
from the shot-noise and read-noise from the sensors and then spreads
throughout the camera pipeline. Hence, the noise is not i.i.d. any more, and
it is signal-dependent and correlated with neighboring pixels due to
demosaicking. Hence, we also train our network to deal with real-noise images.
To train our network for real image denoising, we use Smartphone Image
Denoising Dataset (SIDD) [28], which includes pairs of clean and noisy images
from smartphone cameras. For the choice of ${\mathbf{c}}$, since we do not
know the noise-level, we use the average pixel value of the noisy image of
$4\times 4$ region as our ${\mathbf{c}}$. It can be expressed as

${\mathbf{c}}=Avgpool_{4\times 4}({\mathbf{y}}).$ (14)

Then, the CENet becomes just a simple average pooling operation, and our
tunable denoiser targets different color valued regions separately. That is,
it can be operated within the divide-and-conquer scheme. We choose a widely-
used benchmark for the evaluation: Darmstadt Noise Dataset (DND) [29], which
consists of 50 images with real-noise from 50 scenes and the scenes are
further cropped by the provider which results in 1,000 small patches. The
ground-truth image is not available, and the evaluation can only be done
through online submission.

###  VII-A Analysis

We first analyze the effect of the choice of ${\mathbf{c}}$. As our ablation
investigation, we train a baseline DUBD model within an end-to-end scheme
without any constraints on the conditional variable ${\mathbf{c}}$. Then, we
trained our DUBD-R (R stands for the model for real noise) by using the
spatial average color value as our conditional variable. TABLE III shows the
summary. It is shown that by setting ${\mathbf{c}}$ as Eq. 14, the PSNR has
been improved compared to the baseline model. Also, we note that DUBD-R has a
smaller number of parameters compared to the baseline model by replacing the
CENet by a simple average pooling method.

TABLE III: Ablation investigation for the real-noise denoising.

Dataset | SIDD [28] | DND [29]
---|---|---
Baseline | 39.13 | 39.13
DUBD-R | 39.27 | 39.38

###  VII-B Results

TABLE IV: Average PSNR and SSIM results on DND [29] benchmark. The best results are in red and the second bests are in blue. Method | Parameters | PSNR | SSIM
---|---|---|---
DnCNN+ [9] | 668 K | 37.90 | 0.9430
FFDNet+ [10] | 825 K | 37.61 | 0.9415
CBDNet [30] | 4,365 K | 38.06 | 0.9421
ATDNet [14] | 9,453 K | 39.19 | 0.9526
RIDNet [31] | 1,499 K | 39.26 | 0.9528
VDN [17] | 7,817 K | 39.38 | 0.9518
DUBD-R (Ours) | 2,088 K | 39.38 | 0.9526
DUBD-R+ (Ours) | 2,088 K | 39.44 | 0.9530

In TABLE IV, we demonstrate several comparisons on two real-noise image
benchmark datasets. The comparisons are DnCNN+ [9], FFDNet+ [10], CBDNet [30],
ATDNet [14], RIDNet [31], and VDN [17]. To achieve further performance gain,
we adopt a self-ensemble approach to our method and denote it as “DUBD-R+.” As
shown in the table, our method achieves competent results. Compared to a
recently published blind denoiser VDN [17], our method shows competent results
while requiring a much smaller number of parameters. We also visualize
comparisons on Figure 9 for qualitative evaluation. We expect further
investigation on the choice of ${\mathbf{c}}$ can boost the performance for
real-world noisy images.

Noisy
23.55/0.5185

DnCNN+ [9]
34.51/0.9457

FFDNet+ [10]
34.47/0.9510

CBDNet [30]
35.43/0.9469

ATDNet [14]
36.03/0.9506

RIDNet [31]
37.17/0.9596

VDN [17]
37.34/0.9619

DUBD-R+
37.61/0.9637

Figure 9: Visualization of results on a real-noise image from DND [29]. The
PSNR/SSIM results are written below.

##  VIII Conclusion

In this paper, we have proposed a universal blind denoiser, which can reduce
noise from various environments, including the ones encountered in a real
situation. Based on the idea of divide-and-conquer, we split the original
denoising MAP problem into two inference sub-problems, and we designed new CNN
architectures as sub-problem solvers. Concretely, we introduced an auxiliary
random variable to divide and approximate the original problem. Moreover, on
the choice of the auxiliary random variable, we explicitly reflected our prior
knowledge to augment implicit priors learned from a large-scale dataset. With
the experiments, we have shown that our method is an efficient one in terms of
performance vs. complexity, and shows promising results in various situations
of noise corruption, such as spectrally and spatially varying noises without
variance information. As a result, it also works robustly to the real-world
noise from cameras. The codes are publicly available at
https://www.github.com/JWSoh/DUBD.

## Acknowledgment

This research was supported in part by Samsung Electronics Co., Ltd..

## References

  * [1] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” _Physica D: nonlinear phenomena_ , vol. 60, no. 1-4, pp. 259–268, 1992.
  * [2] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” _IEEE Transactions on Image processing_ , vol. 15, no. 12, pp. 3736–3745, 2006.
  * [3] S. Gu, L. Zhang, W. Zuo, and X. Feng, “Weighted nuclear norm minimization with application to image denoising,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 2862–2869.
  * [4] J. Xu, L. Zhang, and D. Zhang, “A trilateral weighted sparse coding scheme for real-world image denoising,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 20–36.
  * [5] A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image denoising,” in _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)_ , vol. 2. IEEE, 2005, pp. 60–65.
  * [6] W. Dong, L. Zhang, G. Shi, and X. Li, “Nonlocally centralized sparse representation for image restoration,” _IEEE transactions on Image Processing_ , vol. 22, no. 4, pp. 1620–1630, 2012.
  * [7] J. Mairal, F. R. Bach, J. Ponce, G. Sapiro, and A. Zisserman, “Non-local sparse models for image restoration.” in _ICCV_ , vol. 29. Citeseer, 2009, pp. 54–62.
  * [8] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-d transform-domain collaborative filtering,” _IEEE Transactions on image processing_ , vol. 16, no. 8, pp. 2080–2095, 2007.
  * [9] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, “Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising,” _IEEE Transactions on Image Processing_ , vol. 26, no. 7, pp. 3142–3155, 2017.
  * [10] K. Zhang, W. Zuo, and L. Zhang, “Ffdnet: Toward a fast and flexible solution for cnn-based image denoising,” _IEEE Transactions on Image Processing_ , vol. 27, no. 9, pp. 4608–4622, 2018.
  * [11] X. Mao, C. Shen, and Y.-B. Yang, “Image restoration using very deep convolutional encoder-decoder networks with symmetric skip connections,” in _Advances in neural information processing systems_ , 2016, pp. 2802–2810.
  * [12] S. Lefkimmiatis, “Non-local color image denoising with convolutional neural networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 3587–3596.
  * [13] ——, “Universal denoising networks: a novel cnn architecture for image denoising,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 3204–3213.
  * [14] Y. Kim, J. W. Soh, and N. I. Cho, “Adaptively tuning a convolutional neural network by gate process for image denoising,” _IEEE Access_ , vol. 7, pp. 63 447–63 456, 2019.
  * [15] D. Liu, B. Wen, Y. Fan, C. C. Loy, and T. S. Huang, “Non-local recurrent network for image restoration,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 1673–1682.
  * [16] Y. Zhang, K. Li, K. Li, B. Zhong, and Y. Fu, “Residual non-local attention networks for image restoration,” in _ICLR_ , 2019.
  * [17] Z. Yue, H. Yong, Q. Zhao, D. Meng, and L. Zhang, “Variational denoising network: Toward blind noise modeling and removal,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 1688–1699.
  * [18] X. Huang and S. Belongie, “Arbitrary style transfer in real-time with adaptive instance normalization,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 1501–1510.
  * [19] E. Perez, F. Strub, H. De Vries, V. Dumoulin, and A. Courville, “Film: Visual reasoning with a general conditioning layer,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [20] X. Wang, K. Yu, C. Dong, and C. Change Loy, “Recovering realistic texture in image super-resolution by deep spatial feature transform,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 606–615.
  * [21] T. Park, M.-Y. Liu, T.-C. Wang, and J.-Y. Zhu, “Semantic image synthesis with spatially-adaptive normalization,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 2337–2346.
  * [22] E. Agustsson and R. Timofte, “Ntire 2017 challenge on single image super-resolution: Dataset and study,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2017, pp. 126–135.
  * [23] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [24] Y. Chen and T. Pock, “Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 39, no. 6, pp. 1256–1272, 2016.
  * [25] Y. Tai, J. Yang, X. Liu, and C. Xu, “Memnet: A persistent memory network for image restoration,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 4539–4547.
  * [26] J.-B. Huang, A. Singh, and N. Ahuja, “Single image super-resolution from transformed self-exemplars,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 5197–5206.
  * [27] D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” in _Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001_ , vol. 2. IEEE, 2001, pp. 416–423.
  * [28] A. Abdelhamed, S. Lin, and M. S. Brown, “A high-quality denoising dataset for smartphone cameras,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 1692–1700.
  * [29] T. Plotz and S. Roth, “Benchmarking denoising algorithms with real photographs,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , July 2017.
  * [30] S. Guo, Z. Yan, K. Zhang, W. Zuo, and L. Zhang, “Toward convolutional blind denoising of real photographs,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1712–1722.
  * [31] S. Anwar and N. Barnes, “Real image denoising with feature attention,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 3155–3164.