Julien Toulouse
email: julien.toulouse@upmc.fr

Wuming Zhu
email: wuming@lct.jussieu.fr

János G Ángyán
email: janos.angyan@crm2.uhp-nancy.fr

Andreas Savin
email: savin@lct.jussieu.fr

Range-separated density-functional theory with random phase approximation: Detailed formalism and illustrative applications

Using Green-function many-body theory, we present the details of a formally exact adiabatic-connection fluctuation-dissipation density-functional theory based on range separation, which was sketched in Toulouse, Gerber, Jansen, Savin and Ángyán, Phys. Rev. Lett. 102, 096404 (2009). Range-separated density-functional theory approaches combining short-range density functional approximations with long-range random phase approximations (RPA) are then obtained as well-identified approximations on the long-range Green-function selfenergy. Range-separated RPA-type schemes with or without long-range Hartree-Fock exchange response kernel are assessed on rare-gas and alkaline-earth dimers, and compared to range-separated second-order perturbation theory and range-separated coupled-cluster theory.

I. INTRODUCTION

Range-separated density-functional theory has emerged as a powerful approach for improving the accuracy of standard Kohn-Sham (KS) density-functional theory [1,2] applied with usual local or semi-local density-functional approximations, in particular for electronic systems with strong (static) or weak (van der Waals) correlation effects. Based on a separation of the electron-electron interaction into long-range and short-range components, it permits a rigorous combination of a long-range explicit many-body approximation with a shortrange density-functional approximation (see, e.g., Ref. 3 and references therein). Several many-body approximations have been considered for the long-range part: configuration interaction [4,5], multi-configuration self-consistent-field theory [6][7][8], second-order perturbation theory [9][10][11][12][13], coupledcluster theory [14][15][16][17][18], multi-reference second-order perturbation theory [19], and several variants of the random phase approximation (RPA) [20][21][22][23][24].

In the context of the recent revived interest in RPA-type approaches to the electron correlation problem in atomic, molecular and solid-state systems , several range-separated approaches using long-range RPA-type approximations have indeed been proposed and show promising results, in particular for describing weak intermolecular interactions. Toulouse et al. [20] have presented a range-separated RPA-type theory including the long-range Hartree-Fock exchange response kernel. Janesko et al. [21][22][23] have proposed a simpler rangeseparated RPA scheme with no exchange kernel and in which the RPA correlation energy has been rescaled by an empirical coefficient. Paier et al. [24] have added the so-called second-order screened exchange to the latter scheme, which appears to correct the self-interaction error. In all these cases, range separation tends to improve the corresponding fullrange RPA-type approach, avoiding the inaccurate description and slow basis-set convergence of short-range correlations in RPA.

In Ref. 20, only the main lines of range-separated densityfunctional theory with long-range RPA were presented. In this work, we give now all the missing details of the theory. Using Green-function many-body theory, we construct a formally exact adiabatic-connection fluctuation-dissipation density-functional theory based on range separation, without the need of maintaining the one-particle density constant. Range-separated RPA-type schemes are then obtained as wellidentified approximations on the long-range Green-function self-energy. The range-separated RPA-type methods with or without long-range Hartree-Fock exchange response kernel are assessed on rare-gas and alkaline-earth dimers, and compared to range-separated second-order perturbation theory and range-separated coupled-cluster theory. The most tedious details of the theory are given in the appendices.

II. THEORY A. Range-separated density-functional theory

In range-separated density-functional theory (see, e.g., Ref. 3), the exact ground-state energy of an N-electron system is expressed as the following minimization over multideterminant wave functions Ψ E = min

Ψ Ψ| T + Vne + Ŵlr ee |Ψ + E sr Hxc [n Ψ ] , (1)
where T is the kinetic energy operator, Vne is the nuclei-electron interaction operator, Ŵlr ee = (1/2) dr 1 dr 2 w lr ee (r 12 )n 2 (r 1 , r 2 ) is a long-range (lr) electron-electron interaction written with w lr ee (r) = erf(µr)/r and the pair-density operator n2 (r 1 , r 2 ), and E sr Hxc [n] is the corresponding µ-dependent short-range (sr) Hartreeexchange-correlation (Hxc) density functional that Eq. ( 1) defines. The parameter µ in the error function controls the range of the separation. The minimizing wave function, denoted by Ψ lr , yields the exact density. Several approximations [3,7,14,18,[49][50][51] have been proposed for the short-range exchange-correlation (xc) functional E sr xc [n], and an approximate scheme must be used for the long-range wave function part of the calculation. In a first step, the minimization in Eq. ( 1) is restricted to single-determinant wave functions Φ, leading to the rangeseparated hybrid (RSH) approximation [9]

E RSH = min Φ Φ| T + Vne + Ŵlr ee |Φ + E sr Hxc [n Φ ] , (2)
which does not include long-range correlation. The minimizing determinant Φ 0 is given by the self-consistent Euler-Lagrange equation

Ĥ0 |Φ 0 = E 0 |Φ 0 , (3)
where E 0 is the Lagrange multiplier for the normalization constraint and Ĥ0 is the RSH reference Hamiltonian

Ĥ0 = T + Vne + Vlr Hx,HF [Φ 0 ] + Vsr Hxc [n Φ 0 ], (4)
x,HF = dx 1 dx 2 v lr x (x 1 , x 2 )n 1 (x 2 , x 1 ) written with v lr x (x 1 , x 2 ) = -w lr ee (r 12 ) Φ 0 |n 1 (x 1 , x 2 )
|Φ 0 and the oneparticle density-matrix operator n1 (x 1 , x 2 ) expressed with space-spin coordinates x 1 = (r 1 , s 1 ) and x 2 = (r 2 , s 2 ).

The RSH scheme does not yield the exact energy and density, even with the exact short-range functional E sr Hxc [n]. Nevertheless, the RSH approximation can be used as a reference to express the exact energy as

E = E RSH + E lr c , (5)
defining the long-range correlation energy E lr c , for which we will now give an adiabatic connection formula. We introduce the following energy expression with a formal coupling constant λ

E λ = min Ψ Ψ| T + Vne + Vlr Hx,HF [Φ 0 ] + λ Ŵlr |Ψ +E sr Hxc [n Ψ ] , (6)
where the minimization is done over multideterminant wave functions Ψ, Ŵlr is the long-range Møller-Plesset-type fluctuation perturbation operator

Ŵlr = Ŵlr ee -Vlr Hx,HF [Φ 0 ], (7)
and E sr Hxc is the previously-defined λ-independent short-range Hxc functional. The minimizing wave function, denoted by Ψ lr λ , is given by the self-consistent Euler-Lagrange equation

Ĥlr λ |Ψ lr λ = E lr λ |Ψ lr λ , (8)
where E lr λ is the Lagrange multiplier for the normalization constraint and Ĥlr λ is the long-range interacting effective Hamiltonian along the adiabatic connection

Ĥlr λ = T + Vne + Vlr Hx,HF [Φ 0 ] + Vsr Hxc [n Ψ lr λ ] + λ Ŵlr , = Ĥ0 + λ Ŵlr + Vsr Hxc [n Ψ lr λ ] -Vsr Hxc [n Φ 0 ] . (9)
For λ = 1, Eq. ( 6) reduces to Eq. ( 1), and so the physical energy E = E λ=1 and density are recovered. For λ = 0, the minimizing wave function is the RSH determinant Ψ lr λ=0 = Φ 0 and the Hamiltonian of Eq. ( 9) reduces to the RSH reference Hamiltonian, Ĥlr λ=0 = Ĥ0 . Note that, because the density at λ = 0 is not exact, the density necessarily varies along this adiabatic connection. Taking the derivative of E λ with respect to λ, noting that E λ is stationary with respect to Ψ lr λ , and reintegrating between λ = 0 and λ = 1 gives

E = E λ=0 + 1 0 dλ Ψ lr λ | Ŵlr |Ψ lr λ , (10)
with

E λ=0 = Φ 0 | T + Vne + Vlr Hx,HF [Φ 0 ]|Φ 0 +E sr Hxc [n Φ 0 ] = E RSH - Φ 0 | Ŵlr |Φ 0 .
Thus, the long-range correlation energy is

E lr c = 1 0 dλ Ψ lr λ | Ŵlr |Ψ lr λ -Φ 0 | Ŵlr |Φ 0 , (11)
or, equivalently,

E lr c = 1 2 1 0 dλ dx 1 dx 2 dx ′ 1 dx ′ 2 w lr (x 1 , x 2 ; x ′ 1 , x ′ 2 ) ×P lr c,λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ), (12)
where w lr (x 1 , x 2 ;

x ′ 1 , x ′ 2 ) = w lr ee (r 12 )δ(x 1 -x ′ 1 )δ(x 2 -x ′ 2 ) - 1/(N -1) v lr H (r 1 )δ(x 1 -x ′ 1 ) + v lr x (x 1 , x ′ 1 ) δ(x 2 -x ′ 2 )
is the potential corresponding to the perturbation operator Ŵlr and P lr c,λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ) is the correlation part of the two-particle density matrix along the adiabatic connection.

B. Long-range many-body perturbation theory

We now derive a formally exact many-body perturbation theory to calculate the long-range correlation two-particle density matrix P lr c,λ . Details are given in Appendix A. The one-particle Green function G lr λ (1, 2) along the adiabatic connection of Eq. ( 9) in terms of space-spin-time coordinates 1 = (x 1 , t 1 ) and 2 = (x 2 , t 2 ) satisfies the following Dyson equation

G lr λ -1 (1, 2) = G -1 0 (1, 2) -Σ lr λ (1, 2) -∆Σ sr λ (1, 2), ( 13
)
where G 0 (1, 2) is the reference Green function corresponding to the RSH Hamiltonian Ĥ0 , Σ lr λ (1, 2) is the self-energy corresponding to the long-range perturbation operator λ Ŵlr and ∆Σ sr λ (1, 2) is the self-energy correction associated with the short-range potential variation term Vsr Hxc [n

Ψ lr λ ] -Vsr Hxc [n Φ 0 ]
due to the variation of the density [START_REF]The short-range self-energy correction ∆Σ sr λ is wrongly missing in Eq[END_REF]. The long-range selfenergy corresponding to the perturbation operator λ( Ŵlr ee -Vlr Hx,HF [Φ 0 ]) is decomposed into Hartree, exchange and correlation contributions as

Σ lr λ (1, 2) = Σ lr Hxc,λ [G lr λ ](1, 2) -Σ lr Hx,λ [G 0 ](1, 2) = λ Σ lr Hx [G lr λ ](1, 2) -Σ lr Hx [G 0 ](1, 2) +Σ lr c,λ [G lr λ ](1, 2), (14)
where Σ lr Hx [G](1, 2) is the sum of a long-range Hartree selfenergy

Σ lr H [G](1, 2) = -i d3 d4 w lr ee (1, 3)δ(1, 2)δ(3, 4)G(4, 3 + ) = -iδ(1, 2) d3 w lr ee (1, 3)G(3, 3 + ) = δ(1, 2) dr 3 w lr ee (r 13 )n(r 3 ) = δ(1, 2)v lr H [n](r 1 ), (15)
with the instantaneous electron-electron interaction w lr ee (1, 3) = δ(t 1t 3 )w lr ee (r 13 ) and the density extracted from the Green function n(r 3 ) = -i s 3 G(3, 3 + ) (where 3 + stands for x 3 t + 3 with t + 3 = t 3 + η and η is an infinitesimal positive shift), and a long-range exchange self-energy

Σ lr x [G](1, 2) = i d3 d4 w lr ee (1, 3)δ(1, 4)δ(2, 3)G(4, 3 + ) = iw lr ee (1, 2)G(1, 2 + ) = -δ(t 1 -t 2 )w lr ee (r 12 )n 1 (x 1 , x 2 ) = δ(t 1 -t 2 )v lr x [n 1 ](x 1 , x 2 ), (16)
with the one-particle density matrix extracted from the Green function n 1 (x 1 , x 2 ) = -iG(x 1 t 1 , x 2 t + 1 ). The short-range selfenergy correction corresponding to the operator Vsr Hxc [n Ψ lr λ ] -Vsr Hxc [n Φ 0 ] is written as

∆Σ sr λ (1, 2) = Σ sr Hxc [G lr λ ](1, 2) -Σ sr Hxc [G 0 ](1, 2), ( 17
)
where

Σ sr Hxc [G](1, 2) = δ(1, 2)v sr Hxc [n](r 1 )
is the local shortrange Hxc self-energy.

The long-range four-point polarization propagator χ lr λ (1, 2; 1 ′ , 2 ′ ) along the adiabatic connection is given by the solution of the following Bethe-Salpeter-type equation which can be derived from the Dyson equation ( 13) by considering variations with respect to G lr λ [see Appendix A, Eq. (A13)]

χ lr λ -1 (1, 2; 1 ′ , 2 ′ ) = χ lr IP,λ -1 (1, 2; 1 ′ , 2 ′ ) -λ f lr Hx (1, 2; 1 ′ , 2 ′ ) -f lr c,λ (1, 2; 1 ′ , 2 ′ ), (18)
where χ lr IP,λ (1, 2;

1 ′ , 2 ′ ) = -iG lr λ (1, 2 ′ )G lr λ (2, 1 ′
) is an independent-particle (IP) polarization propagator, and

λ f lr Hx (1, 2; 1 ′ , 2 ′ ) = iλ δΣ lr Hx [G lr λ ](1, 1 ′ )/δG lr λ (2 ′ , 2) and f lr c,λ (1, 2; 1 ′ , 2 ′ ) = i δΣ lr c,λ [G lr λ ](1, 1 ′ )/δG lr λ (2 ′ , 2
) are longrange Hartree-exchange and correlation kernels. Note that these kernels only stem from the self-energy term Σ lr Hxc,λ [G lr λ ] in Eq. ( 13) that corresponds to the two-electron interaction λ Ŵlr ee , the other self-energy contributions which come from the one-electron terms are absorbed in the definition of χ lr λ (1, 2; 1 ′ , 2 ′ ). The Hartree kernel is obtained from Eq. ( 15)

f lr H (1, 2; 1 ′ , 2 ′ ) = w lr ee (1, 2)δ(1, 1 ′ )δ(2, 2 ′ ) = w lr ee (r 12 )δ(t 1 -t 2 )δ(1, 1 ′ )δ(2, 2 ′ ), (19)
while the HF-like exchange kernel is obtained from Eq. ( 16)

f lr x (1, 2; 1 ′ , 2 ′ ) = -w lr ee (1, 2)δ(1, 2 ′ )δ(1 ′ , 2) = -w lr ee (r 12 )δ(t 1 -t 2 )δ(1, 2 ′ )δ(1 ′ , 2). ( 20
)
The fluctuation-dissipation theorem is then used to express P lr c,λ as [see Appendix A, Eq. (A24)]

P lr c,λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ) = - ∞ -∞ dω 2πi e iω0 + χ lr λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ; ω) -χ 0 (x 1 , x 2 ; x ′ 1 , x ′ 2 ; ω) +∆ lr λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ), (21)
where χ lr λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ; ω) is the frequency-dependent Fourier transform of the one-time-interval polarization propagator

χ lr λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ; τ = t 1 -t 2 ) = χ lr λ (x 1 t 1 , x 2 t 2 ; x ′ 1 t + 1 , x ′ 2 t + 2 ), χ 0 (x 1 , x 2 ; x ′ 1 , x ′ 2 ; ω)
is the equivalent quantity for the RSH reference Hamiltonian (at λ = 0), and

∆ lr λ (x 1 , x 2 ; x ′ 1 , x ′ 2 )
is the contribution coming from the variation of the one-particle density matrix along the adiabatic connection. The expression of ∆ lr λ in terms of the Green functions G lr λ and G 0 is straightforward but it is sufficient to write it as

∆ lr λ = Γ[G lr λ ]-Γ[G 0 ] where Γ is a known functional given in Appendix A [Eq. (A22)].
So far, the theory is in principle exact. In the following we consider two possible approximations. The RPA approximation

Σ lr xc,λ = 0, (22)
corresponds to neglecting long-range exchange-correlation in all one-electron properties. Indeed, with this approximation, one can check that G lr λ = G 0 is a solution of the Dyson equation (13), i.e. the Green function remains unchanged along the adiabatic connection. It follows that ∆ lr λ = 0, f lr xc,λ = 0 and

χ lr IP,λ (1, 2; 1 ′ , 2 ′ ) = -iG 0 (1, 2 ′ )G 0 (2, 1 ′ ) = χ 0 (1, 2; 1 ′ , 2 ′ ). Sim- ilarly, the RPAx approximation Σ lr c,λ = 0, (23)
corresponds to neglecting long-range correlation only in all one-electron properties. Again, this approximation implies that the Green function remains unchanged along the adiabatic connection, i.e. G lr λ = G 0 and it follows that ∆ lr λ = 0, f lr c,λ = 0 and χ lr IP,λ = χ 0 . As different terminologies are used in the quantum chemistry and condensed-matter physics literature, let us stress that what we call RPA here corresponds to a response equation (18) with no exchange-correlation kernel (and it is also sometimes called linear response time-dependent Hartree theory or direct RPA), and what we call RPAx corresponds to a response equation with an additional HF-like exchange kernel (and it is also sometimes called linear response timedependent Hartree-Fock theory or full RPA).

C. Expressions in an orbital basis

The RPA or RPAx equations in an orbital basis are derived in details in Appendix B. In the basis of RSH spin orbitals, the long-range RPA or RPAx correlation energy writes

E lr c = 1 2 1 0 dλ ia, jb ib| ŵlr ee |a j (P lr c,λ ) ia, jb , (24)
where i and j refer to occupied spin orbitals, and a and b to virtual spin orbitals, ib| ŵlr ee |a j are the two-electron integrals with long-range interaction, and (P lr c,λ ) ia, jb are the matrix elements of the correlation two-particle density matrix. The one-electron terms v lr H and v lr x in the perturbation operator in Eq. ( 12) do not contribute to E lr c because of the occupiedvirtual/occupied-virtual structure of the two-particle density matrix in RPA or RPAx. Following the technique proposed by Furche [26], P lr c,λ can be obtained as

P lr c,λ = A lr λ -B lr λ 1/2 M lr λ -1/2 A lr λ -B lr λ 1/2 -1. ( 25
)
with

M lr λ = A lr λ -B lr λ 1/2 A lr λ + B lr λ A lr λ -B lr λ 1/2
, and the orbital rotation Hessians

A lr λ ia, jb = (ǫ a -ǫ i )δ i j δ ab +λ ib| ŵlr ee |a j -ξ ib| ŵlr ee | ja , (26a)
B lr λ ia, jb = λ ab| ŵlr ee |i j -ξ ab| ŵlr ee | ji . (26b)
where ǫ i are the RSH orbital eigenvalues, and ξ = 0 or ξ = 1 for RPA and RPAx, respectively. For spin-restricted closedshell calculations, the correlation energy writes in terms of spatial orbitals

E lr c = 1 2 1 0 dλ ia, jb ib| ŵlr ee |a j ( 1 P lr c,λ ) ia, jb , (27)
where i and j now refer to occupied spatial orbitals, and a and b to virtual spatial orbitals, and 1 P lr c,λ is the spin-singletadapted correlation two-particle density matrix obtained as

1 P lr c,λ = 2 1 A lr λ -1 B lr λ 1/2 1 M lr λ -1/2 1 A lr λ -1 B lr λ 1/2 -1 ,( 28
)
with 1 M lr λ = 1 A lr λ -1 B lr λ 1/2 1 A lr λ + 1 B lr λ 1 A lr λ -1 B lr λ 1/2
, and the singlet orbital rotation Hessians

1 A lr λ ia, jb = (ǫ a -ǫ i )δ i j δ ab +λ 2 ib| ŵlr ee |a j -ξ ib| ŵlr ee | ja , (29a) 1 B lr λ ia, jb = λ 2 ab| ŵlr ee |i j -ξ ab| ŵlr ee | ji . ( 29b
)
Only singlet excitations contribute to Eq. ( 27), since the twoelectron integrals involved vanish for triplet excitations.

In Eq. ( 25), it is assumed that A lr λ + B lr λ and A lr λ -B lr λ are positive definite. In RPA, this is always the case. On the contrary, in RPAx, this is not always the case, i.e. instabilities can be encountered, and Eq. ( 25) can fail. In spin-restricted closedshell formalism, one may encounter singlet instabilities in the RPAx theory defined here, for example when dissociating a bond, but not triplet instabilities since triplet excitations do not contribute at all. In practice, singlet instabilities are usually not encountered for weakly-interacting closed-shell systems. Note that other variants of RPA-type correlation energy expressions using a HF exchange response kernel, such as the plasmon formula [38,[START_REF] Mclachlan | [END_REF]54] or the equivalent ring coupledcluster-doubles theory [38], require contributions from both singlet and triplet excitations, and are thus subject to triplet instabilities (e.g. in a system such as Be 2 ).

Similarly to the notation used in Ref. 20, the rangeseparated method obtained by adding to the RSH energy the long-range RPAx correlation energy [ξ = 1 in Eqs. (26) or Eqs. (29] will be referred to as RSH+lrRPAx. For consistency, the range-separated method obtained by adding to the RSH energy the long-range RPA correlation energy [ξ = 0 in Eqs. (26) or Eqs. (29)] will be referred to as RSH+lrRPA, although it is equivalent to the method called "LC-ωLDA+dRPA" in Refs. 21-24 in the special case of the short-range LDA functional. At second order in the electronelectron interaction, the RSH+lrRPAx method reduces to the range-separated method of Ref. 9 based on long-range secondorder Møller-Plesset perturbation theory, to which we will refer as RSH+lrMP2. Since RPA approaches can be seen as simple approximations to coupled-cluster theory [38], the RSH+lrRPA and RSH+lrRPAx methods bear some resemblance to the range-separated method of Ref. 14 where the long-range correlation energy is evaluated by coupled-cluster theory (with single, double and perturbative triple excitations), to which we will refer as RSH+lrCCSD(T).

We note that one can develop long-rang many-body perturbation theories starting from other references than the RSH reference. For example, starting from the usual (approximate) Kohn-Sham reference could be appropriate for solid-state systems. For the finite systems considered here, RSH is a good reference, as confirmed by other authors [23]. [START_REF] Werner | Molpro, version 2008.2, a package of ab initio programs[END_REF], implementing equations ( 27)- (29). We first perform a self-consistent RSH calculation with the short-range PBE xc functional of Ref. 14 (this RSH calculation could also be referred to as "lrHF+srPBE", a notation closer to the one used by other authors [14]) and add the long-range MP2, RPA, RPAx or CCSD(T) correlation energy calculated with RSH orbitals. For RPA or RPAx, the λ-integration in Eq. ( 27) is done by a 7-point Gauss-Legendre quadrature [26]. The range separation parameter is taken at µ = 0.5 bohr -1 , in agreement with previous studies [START_REF] Gerber | [END_REF], without trying to adjust it for each system. To show the dependence on the orbitals, the full-range RPA calculations have been done with PBE [57] and HF orbitals, which will be denoted by PBE+RPA and HF+RPA, respectively [58]. The full-range MP2, RPAx and CCSD(T) calculations have been done with HF orbitals, and thus, for notation consistency, will be denoted by HF+MP2, HF+RPAx and HF+CCSD(T), respectively. We use large Dunning basis sets [59][60][61][62][63][64][65]. Core electrons are kept frozen in all the full-range and rangeseparated MP2, RPA, RPAx and CCSD(T) calculations (i.e. only excitations of valence electrons are considered). The basis set superposition error (BSSE) is removed by the counterpoise method. For the alkaline-earth dimers, it has been checked than adding diffuse basis functions or core excitations do not change significantly the results. Extrapolations to the complete basis set (CBS) limit have also been considered for some systems. For the full-range methods, the standard three-point exponential formula for the HF (or KS) reference E HF (n) = E HF (CBS) + Ae -Bn with the cardinal number n = 3, 4, 5, and two-point formula for the correlation energy E c (n) = E c (CBS)+C/n 3 with n = 4, 5 have been used. For the range-separated methods, we have also used these two formulas for the RSH reference and the long-range correlation energy, even though in this case the dependence on the cardinal number would deserve a detailed study.

III. COMPUTATIONAL DETAILS

All calculations have been performed with a development version of MOLPRO 2008

For each dimer interaction energy curve, we choose 16 to 20 intermolecular distances, with denser sampling around the equilibrium distance. A third-order polynomial is used for interpolation. The hard core radius is taken as the distance where the interaction energy is 0, and the equilibrium distance and binding energy are from the minimum of the interpolated interaction energy curve. The harmonic vibrational frequency is obtained from the second-order derivative of the energy curve at the equilibrium distance. For C 6 dispersion coefficients, the interaction energy E int is calculated at seven extra distances R i from 30 to 60 bohr, and the coefficient is estimated by averaging with the following formula

C 6 = exp         1 7 7 i=1 (ln |E int (R i )| + 6 ln(R i ))         , (30)
similarly to what has been done in Ref. 22.

IV. APPLICATIONS A. Basis set dependence

The convergence of the equilibrium binding energy of Ar 2 with respect to the basis set size up to the CBS limit for the full-range methods HF+MP2, PBE+RPA, HF+RPA, HF+CCSD(T) and for the range-separated methods RSH+lrMP2, RSH+lrRPA, RSH+lrRPAx, RSH+lrCCSD(T) is represented in Fig. 1. Full-range RPA with PBE orbitals has a very strong dependence on the basis size, as already noted (e.g. Refs. 20,26). Full-range RPA with HF orbitals has a bit weaker basis dependence, similar to full-range HF+MP2, HF+RPAx and HF+CCSD(T). All the range-separated methods have essentially identical, very favorable basis set convergence. Since the slow convergence of full-range methods is FIG. 1: (Color online) Basis set dependence of the equilibrium binding energy of Ar 2 for different full-range and range-separated methods, presented as the percentage of the binding energy recovered with respect to the CBS limit (aVTZ, aVQZ and aV5Z stand for aug-cc-pVTZ, aug-cc-pVQZ and aug-cc-pV5Z, respectively). related to the explicit description of short-range correlation, it is not surprising that range-separated methods have a faster convergence because they leave the description of short-range correlation to the short-range density functional. These results are consistent with other studies, e.g. Refs. 22, 24. Note that, with the aug-cc-pV5Z basis set, all the range-separated methods are essentially converged (98-99% of the CBS binding energy), therefore we will not use CBS extrapolations in the following. However, one should keep in mind that with this basis set the full-range methods are not yet fully converged, with about 90% of the CBS binding energy.

B. Rare-gas dimers

In Fig. 2, the interaction energy curves of He 2 , Ne 2 , Ar 2 and Kr 2 , obtained with the full-range and range-separated methods are compared. As already known, full-range HF+MP2 underestimates the interaction energy for the smallest systems He 2 and Ne 2 , and overestimates it for the largest systems Ar 2 and Kr 2 . Full-range PBE+RPA gives an almost dissociative curve for He 2 , and largely underestimates the interaction energy for Ne 2 , Ar 2 and Kr 2 . Using HF orbitals in full-range RPA drastically improves the interaction energy curve for He 2 , and to a least extend for Ne 2 , but gives less binding for Ar 2 and Kr 2 . Full-range HF+RPAx significantly improves over full-range HF+RPA, but still gives underestimated interaction energies. It can be noted that full-range HF+RPAx yields interaction energy curves almost identical to the full-range HF+MP2 curves for He 2 and Ne 2 , and almost identical to the full-range PBE+RPA curves for Ar 2 and Kr 2 . Full-range HF+CCSD(T) gives systematically quite accurate interaction energies. Quite similarly to full-range HF+MP2, the range-separated RSH+lrMP2 underestimates the interaction energy for He 2 and Ne 2 , and overestimates it for Ar 2 and Kr 2 . RSH+lrRPA tends to improve over both full-range PBE+RPA and HF+RPA but still leads to significantly underestimated interaction energies. RSH+lrRPAx improves over TABLE I: Hard-core radii σ (bohr), equilibrium distances R e (bohr), equilibrium binding energies D e (mhartree), harmonic vibrational frequencies ω e (cm -1 ) and dispersion coefficients C 6 for ten homonuclear and heteronuclear rare-gas dimers from different full-range and rangeseparated methods with aug-cc-pV5Z basis. Mean absolute percentage errors (MA%E) are also given. both RSH+lrRPA and full-range HF+RPAx; it still systematically underestimates the interaction energy at equilibrium, but appears quite accurate at medium and large distances. On the contrary, RSH+lrCCSD(T) systematically overestimates the interaction energy at medium and large distances.

The hard-core radii, equilibrium distances, equilibrium binding energies, harmonic vibrational frequencies and dispersion coefficients C 6 for ten homonuclear and heteronu- clear rare-gas dimers calculated with the full-range and rangeseparated methods are given in Table I. The trends seen in Fig. 2 are confirmed. Full-range RPA (with PBE or HF orbitals) yields very inaccurate equilibrium properties. Fullrange HF+RPAx improves over full-range HF+RPA (with the exception of C 6 coefficients which turn out to be quite good in PBE+RPA for these systems) but the errors remain large. Range separation largely improves RPA and RPAx. RSH+lrRPAx gives much better equilibrium properties than RSH+lrRPA, with mean absolute percentage errors smaller by more than a factor of two, while these two methods give similar accuracy for C 6 coefficients. Full-range HF+MP2 is reasonably accurate and range separation has a much smaller impact on it. For these systems, RSH+lrMP2 gives an over-all similar accuracy than RSH+RPAx, although the C 6 coefficients tend to be globally more accurate in RSH+lrRPAx. Full-range HF+CCSD(T) gives the best results. Surprisingly, range separation tends to deteriorate the accuracy of CCSD(T), especially for C 6 coefficients. Nevertheless, among the range-separated methods, RSH+lrCCSD(T) still gives the best equilibrium properties.

C. Alkaline-earth dimers

In Fig. 3, the interaction energy curves of Be 2 , Mg 2 and Ca 2 , obtained with the full-range and range-separated methods are compared. These systems have static corre- lation effects, especially Be 2 , and are thus more challenging for the single-reference methods tested here. Full-range PBE+RPA gives unphysical interaction energy curves, with a large bump for Be 2 , and with essentially no bond for Mg 2 and Ca 2 . Full-range HF+RPA yields an almost dissociative curve for Be 2 with no bump (which is consistent with Ref. 43), and physically reasonable curves for Mg 2 and Ca 2 . Full-range HF+RPAx moderately improves over fullrange HF+RPA. Among the full-range methods, HF+MP2 and HF+CCSD(T) clearly give the best interaction energy curves. As for rare-gas dimers, RSH+lrRPA always largely underestimates the interaction energy. RSH+lrMP2 and RSH+lrRPAx give much less underestimated interaction energies, with RSH+lrMP2 being a bit more accurate for Mg 2 and Ca 2 . While RSH+lrCCSD(T) largely overestimates the interaction energy for Be 2 , it is remarkably accurate for Mg 2 and Ca 2 . We note that RSH+lrCCSD(T) could be made more accurate for Be 2 by choosing a larger range-separation parameter µ [71].

The hard-core radii, equilibrium distances, equilibrium binding energies, harmonic vibrational frequencies and dispersion coefficients C 6 for Be 2 , Mg 2 and Ca 2 are given in Table II. It is confirmed that range separation largely improves the equilibrium properties of RPA and RPAx. Again, RSH+lrRPAx is much more accurate than RSH+lrRPA, with mean absolute percentage errors smaller by about a factor of two. Range separation also overall brings a significant improvement in MP2. Among the range-separated methods, RSH+lrCCSD(T) gives the best equilibrium properties.

V. CONCLUSIONS

We have expounded the details of a formally exact adiabatic-connection fluctuation-dissipation densityfunctional theory based on range separation. Range-separated density-functional theory with random phase approximations including or not the long-range Hartree-Fock exchange response kernel (referred to as RSH+lrRPA and RSH+lrRPAx, respectively) are then obtained as well-identified approximations on the long-range Green-function self-energy [Eqs. (22) and (23)]. The long-range Green function does not vary along the adiabatic connection at the RSH+lrRPA and RSH+lrRPAx levels, which makes these schemes relatively simple compared to the exact theory. In practice, RSH+lrRPA and RSH+lrRPAx have been applied in a spin-restricted closed-shell formalism, in which both schemes only include spin-singlet orbital excitations, and thus are not subject to triplet instabilities. These range-separated RPA-type schemes have been tested on rare-gas and alkaline-earth dimers, featuring challenging weak (van der Waals) interactions. Both range separation and inclusion of the exact Hartree-Fock response kernel largely improve the accuracy of RPA. The RSH+lrRPAx method appears as a reasonably accurate method for weak interactions, but globally less accurate for equilibrium properties than the more intensive range-separated coupled-cluster method. Although, for the small systems considered here, rangeseparated second-order perturbation theory (RSH+lrMP2) turns out to yield results similarly as accurate as those from RSH+lrRPAx (and in fact more accurate for Mg 2 and Ca 2 ), a recent investigation [72] shows that RSH+lrRPAx corrects the overestimation of the binding energy in RSH+lrMP2 for larger weakly-interacting stacked complexes, such as the benzene dimer.

Adiabatic connection

We consider the following adiabatic connection defined by the λ-dependent energy

E λ = min Ψ Ψ| K0 + λ Ŵ|Ψ + F[n Ψ ] , (A1)
where K0 is an arbitrary one-particle Hamiltonian, Ŵ is a perturbation operator (generally, the sum of a two-particle operator Ŵee and an one-particle operator) and F[n] is a λindependent density functional. The minimizing multideterminant wave function Ψ λ satisfies the Euler-Lagrange equation

Ĥλ |Ψ λ = E λ |Ψ λ , (A2)
where E λ is the Lagrange multiplier for the normalization constraint, and Ĥλ is the effective Hamiltonian along the adiabatic connection

Ĥλ = K0 + λ Ŵ + Vλ , (A3)
where Vλ = dr δF[n Ψ λ ]/δn(r) n(r) is a self-consistent oneparticle potential operator. Note that Ĥλ=1 is not necessarily the physical Hamiltonian. This adiabatic connection links the energy of interest E λ=1 to the reference energy E λ=0 = Φ 0 | K0 |Φ 0 + F[n Φ 0 ] calculated with the single-determinant wave function Φ 0 = Ψ λ=0 of the reference Hamiltonian Ĥ0 = K0 + V0 . The one-particle density is not kept constant with respect to λ.

An adiabatic connection formula for E λ=1 is found by taking the derivative of E λ with respect to λ, noting that E λ is stationary with respect to Ψ λ , and reintegrating between λ = 0 and λ = 1

E λ=1 = E λ=0 + 1 0 dλ Ψ λ | Ŵ|Ψ λ . (A4)
The correlation energy, defined as E c = E λ=1 -E λ=0 -(dE λ /dλ) λ=0 where (dE λ /dλ) λ=0 = Φ 0 | Ŵ|Φ 0 is the firstorder energy correction, is thus given by

E c = 1 0 dλ Ψ λ | Ŵ|Ψ λ -Φ 0 | Ŵ|Φ 0 , (A5)
or, equivalently, in the representation of space-spin coordinates x = (r, s)

E c = 1 2 1 0 dλ dx 1 dx 2 dx ′ 1 dx ′ 2 w(x 1 , x 2 ; x ′ 1 , x ′ 2 ) ×P c,λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ), (A6)
where w(x 1 , x

2 ; x ′ 1 , x ′
2 ) is the interaction potential corresponding to the operator Ŵ and P c,λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ) is the correlation part of the two-particle density matrix along the adiabatic connection.

This exposition encompasses both standard full-range many-body theory and range-separated density-functional theory. Indeed, if K0 is the Hartree-Fock Hamiltonian (i.e., K0 = T + Vne + VHx,HF ), Ŵ is the standard Møller-Plesset fluctuation perturbation operator (i.e., Ŵ = Ŵee -VHx,HF ) and F[n] = 0 then Eq. (A6) yields the full-range many-body correlation energy, defined with respect to the Hartree-Fock energy. Similarly, with the corresponding long-range operators K0 = T + Vne + Vlr Hx,HF and Ŵ = Ŵlr ee -Vlr Hx,HF and the shortrange density functional F[n] = E sr Hxc [n], Eq. (A6) yields now the long-range correlation energy, defined with respect to the RSH energy [Eq. ( 5)].

One-particle Green function

The one-particle Green function along the adiabatic connection is defined as

G λ (1, 2) = -i Ψ λ |T [ ψλ (1) ψ † λ (2)]|Ψ λ , (A7)
where 1 = (x 1 , t 1 ) and 2 = (x 2 , t 2 ) refer to space-spin and time coordinates, ψλ (1) = e i Ĥλ t 1 ψ(x 1 )e -i Ĥλ t 1 and ψ † λ (2) = e i Ĥλ t 2 ψ † (x 2 )e -i Ĥλ t 2 are the annihilation and creation operators in the Heisenberg picture, and T is the Wick time-ordering operator.

A Dyson-type equation connects the inverse of G λ to the inverse of the Green function associated with the one-electron Hamiltonian K0 + Vλ , denoted by G V,λ ,

G -1 λ (1, 2) = G -1 V,λ (1, 2) -Σ λ (1, 2), ( A8
)
which can be considered as the definition of the self-energy Σ λ . In turn, the inverse of G V,λ can be expressed from the inverse of the Green function G 0 of the reference Hamiltonian

Π λ (ω) -1 = Π 0 (ω) -1 -F λ = - A λ B λ B * λ A * λ -ω 1 0 0 -1 , (B7)
where A λ and B λ are the so-called orbital rotation Hessians (A λ ) ia, jb = (ǫ aǫ i )δ i j δ ab +λ ib| ŵee |a jξ ib| ŵee | ja , (B8a) (B λ ) ia, jb = λ ab| ŵee |i jξ ab| ŵee | ji .

We need to consider the linear response non-Hermitian eigenvalue equation

A λ B λ B * λ A * λ X n,λ Y n,λ = ω n,λ 1 0 0 -1 X n,λ Y n,λ , (B9)
whose solutions come in pairs: positive excitation energies ω n,λ with eigenvectors X n,λ , Y n,λ , and opposite (de-)excitation energies -ω n,λ with eigenvectors Y * n,λ , X * n,λ . Choosing the normalization of the eigenvectors so that X † n,λ X m,λ -Y † n,λ Y m,λ = δ nm , the supermatrix Π λ (ω) can be expressed as the following spectral representation (where the sum is over eigenvectors with positive excitation energies)

Π λ (ω) = n 1 ω -ω n,λ + i0 + X n,λ Y n,λ X † n,λ Y † n,λ (B10)
- 1 ω + ω n,λ -i0 + Y * n,λ X * n,λ Y * † n,λ X * † n,λ . (B11)
The fluctuation-dissipation theorem [Eq. (A24)] leads to the supermatrix representation of the correlation part of the twoparticle density matrix P c,λ (using contour integration in the upper half of the complex plane) where out of the integrals ps| ŵ|qr associated with the general perturbation operator only the integrals of the type ib| ŵee |a j associated with the two-electron contribution of the perturbation operator survive because of the occupiedvirtual/occupied-virtual structure of the two-particle density matrix. Using now real spin orbitals, the correlation energy can be simplified to where remains only the contribution from the spin-singletadapted correlation two-particle density matrix ( 1 P c,λ ) ia, jb = σ 1 =↑,↓ σ 2 =↑,↓ (P c,λ ) iσ 1 aσ 1 , jσ 2 bσ 2 , which can be calculated similarly as before

P c,λ = - ∞ -∞ dω 2πi e iω0 + [Π λ (ω) -Π 0 (ω)] = n       Y * n,λ Y * † n,λ Y * n,λ X * † n,λ X * n,λ Y * † n,λ X * n,λ X * † n,λ       - 0 0 0 1 , (
E c = 1
1 P c,λ = 2        n 1 x n,λ + 1 y n,λ 1 x n,λ + 1 y n,λ T -1        = 2 1 A λ -1 B λ 1/2 1 M -1/2 λ 1 A λ -1 B λ 1/2 -1 , (B26)
where 1

M λ = 1 A λ -1 B λ 1/2 1 A λ + 1 B λ 1 A λ -1 B λ 1/2 .

2 FIG. 2 :

 22 FIG. 2: (Color online) Interaction energy curves of He 2 , Ne 2 , Ar 2 and Kr 2 calculated by different full-range (left) and range-separated (right) methods. The basis is aug-cc-pV5Z. The accurate curves are from Ref. 66.

2 FIG. 3 :

 23 FIG. 3: (Color online) Interaction energy curves of Be 2 , Mg 2 and Ca 2 calculated by full-range (left) and range-separated (right) methods. The basis is cc-pV5Z. The accurate curves are from Refs. 67, 68 and 69.

  B12) the simple contribution coming from Π 0 (ω) resulting from its diagonal form [Eqs. (B3)], and the correlation energy [Eq. (A6)] has the following expression in spin-orbital basis ŵee |a j (Y n,λ ) * ia (Y n,λ ) jb+ i j| ŵee |ab (Y n,λ ) * ia (X n,λ ) jb + ab| ŵee |i j (X n,λ ) * ia (Y n,λ ) jb + a j| ŵee |ib (X n,λ ) * ia (X n,λ ) jbδ i j δ ab ,(B13)

  ib| ŵee |a j (P c,λ ) ia, jb , (B14)Performing the sums over spins in the correlation energy expression of Eq. (B14), one gets, for real spatial orbitals, ib| ŵee |a j ( 1 P c,λ ) ia, jb , (B25)

TABLE II :

 II Hard-core radii σ (bohr), equilibrium distances R e (bohr), equilibrium binding energies D e (mhartree), harmonic vibrational frequencies ω e (cm -1 ) and dispersion coefficients C 6 for Be 2 , Mg 2 and Ca 2 from different full-range and range-separated methods with cc-pV5Z basis. Mean absolute percentage errors (MA%E) are also given.

		HF+MP2 PBE+RPA HF+RPA HF+RPAx HF+CCSD(T) RSH+lrMP2 RSH+lrRPA RSH+lrRPAx RSH+lrCCSD(T) Estimated exact
	Be 2 σ	4.44	4.34	5.59	5.30	4.16	4.25	4.50	4.27	3.87	4.01 a
	R e D e	5.15 1.92	4.60 0.58	7.48 0.39	7.17 0.56	4.71 2.70	4.92 2.95	5.08 1.24	4.92 2.81	4.54 6.92	4.63 a 4.31 a
	ω e	139	297	34	37	242	199	152	198	315	267 a
	C 6	256	164	138	180	195	232	149	213	274	214 d
	Mg 2
	σ	6.44	8.30	7.02	6.83	6.29	6.40	6.98	6.49	6.13	6.10 b
	R e	7.66	10.72	8.28	8.11	7.48	7.59	8.23	7.68	7.31	7.35 b
	D e	1.62	0.09	0.70	0.96	1.67	1.43	0.65	1.24	1.92	1.93 b
	ω e	47	7.9	31	35	48	45	30	42	52	51.1 b
	C 6	686	405	364	485	616	571	349	494	671	627 d
	Ca 2
	σ	7.29	-	7.57	7.49	7.07	7.04	7.33	7.11	6.85	6.88 c
	R e	8.57	-	8.76	8.72	8.30	8.25	8.47	8.30	8.05	8.09 c
	D e ω e	3.85 56	--	2.37 44	2.78 47	4.71 64	4.03 60	2.48 50	3.55 57	5.10 68	5.02 c 63.7 c
	C 6	2574	1335	1301	1710	2311	2090	1173	1617	2224	2221 d
	MA%E (%)
	σ	7.4	-	22	18	3.2	4.4	11	5.4	1.5	0.0
	R e	7.1	-	28	24	2.0	3.8	8.7	4.5	1.0	0.0
	D e	32	-	69	61	19	26	63	33	21	0.0
	ω e	23	-	53	48	5.3	14	35	18	9.1	0.0
	C 6	15	33	40	21	5.0	7.7	41	16	12	0.0
	a From Ref. 67
	b From Ref. 68
	c From Ref. 69
	d From Ref. 70

Acknowledgments

We thank J. F. Dobson and T. Gould for numerous discussions on RPA during the FAST (French-Australian Science and Technology) workshop at Griffith University, Australia. We also thank G. Jansen for discussions. This work was partly supported by ANR (French national research agency) via contract number 07-BLAN-0272 (Wademecom).

Appendix A: Adiabatic-connection fluctuation-dissipation density-functional theory

In this appendix, we outline a general, formally exact adiabatic-connection fluctuation-dissipation densityfunctional theory, using Green-function many-body theory. For further details on standard Green's function theory, see e.g. Refs. [73][START_REF] Gross | Many-particle theory[END_REF][START_REF] Onida | [END_REF][76].

, where v λ and v 0 are the one-electron potentials associated with Vλ and V0 , respectively.

For time-independent Hamiltonians, the Green function only depends on the time difference τ = t 1t 2 , so one defines G λ (x 1 , x 2 ; τ) = G λ (x 1 t 1 , x 2 t 2 ), which has a discontinuity at τ = 0. The one-particle density matrix n 1,λ (x 1 , x 2 ) = Ψ λ |n 1 (x 1 , x 2 )|Ψ λ , with n1 (x 1 , x 2 ) = ψ † (x 2 ) ψ(x 1 ), can be obtained from the limit τ → 0 - n 1,λ (x 1 , x 2 ) = -iG λ (x 1 , x 2 ; τ = 0 -).

(A9)

Four-point polarization propagator

The four-point polarization propagator along the adiabatic connection is defined as

) where G 2,λ is the two-particle Green function

A11) Alternatively, using the Schwinger derivative technique, χ λ can be expressed as the functional derivative of the oneparticle Green function with respect to the two-point potential v λ (see, e.g., Refs. 73, 76)

The four-point polarization propagator satisfies a so-called Bethe-Salpeter equation that directly stems from the Dyson equation of Eq. (A8). Considering variations with respect to iG λ (achieved through variations of v λ ) yields

The term on the left-hand side of Eq. (A13) gives straightforwardly

where

) is a so-called independent-particle (IP) polarization propagator [START_REF]The inverse of a 4-point function χ(1, 2; 1 ′ , 2 ′ )[END_REF]. The first term on the right-hand side of Eq. (A13) gives the inverse of the four-point polarization propagator, according to Eq. (A12),

and the second term is the so-called Bethe-Salpeter four-point kernel

and finally, using Eqs. (A14)-(A16) in Eq. (A13), the Bethe-Salpeter equation for χ λ writes

Fluctuation-dissipation theorem

Similarly to the expression of the one-particle density matrix in terms of the one-particle Green function [Eq. (A9)], the two-particle density matrix can be extracted from the polarization propagator. Defining χ λ (x 1 , x

2 ), i.e. the polarization propagator with times t ′ 1 → t + 1 and t ′ 2 → t + 2 which depends only on the time difference τ = t 1t 2 , it is easy to check that in the limit τ → 0 -, after applying the time-ordering operator in Eq. (A11) and using Eq. (A9), one has the following relation

The two-particle density matrix n 2,λ (x 1 , x 2 ;

The correlation part of the two-particle density matrix P c,λ = n 2,λn 2,λ=0 is thus

where χ 0 is the polarization propagator of the non-interacting reference system for λ = 0, and ∆ λ is a term coming from the variation of the one-particle density matrix along the adiabatic connection

. (A21) Using Eq. (A9), one can also express this term with the Green function as

where we define the functional Γ as

Finally, introducing the Fourier transform of χ λ (x 1 , x 2 ; x ′ 1 , x ′ 2 ; τ) in terms of the frequency ω,

we arrive at the form of the fluctuation-dissipation that we use

Appendix B: Random phase approximation in an orbital basis

In this appendix, we give the working equations in an orbital basis resulting from the many-body theory outlined in Appendix A, in the special case of a random phase approximation (RPA)-type simplification. For further details, see e.g. Refs. 26, 53, 78, 79.

Expressions in a spin-orbital basis

In the RPA and RPAx approximations, the Green function does not vary along the adiabatic connection, i.e. G λ = G 0 , which implies that the independent-particle polarization propagator [Eq. (A14)] is just the non-interacting reference polarization propagator, χ IP,λ (1, 2; 1 ′ , 2 ′ ) = -iG 0 (1, 2 ′ )G 0 (2, 1 ′ ) = χ 0 (1, 2; 1 ′ , 2 ′ ), and in the fluctuation-dissipation theorem of Eq. (A24) the term coming from the variation of the oneparticle density matrix vanishes, ∆ λ = 0.

The frequency-dependent non-interacting polarization propagator has the following well-known Lehmann representation

where φ p (x) and ǫ p are the spin orbitals and corresponding eigenvalues of the reference system, and i and a run over occupied and virtual spin orbitals, respectively. Hence, χ 0 can be completely represented in the basis of spin-orbital products, φ * p (x ′ 1 )φ q (x 1 ), where p refer to an occupied orbital and q to a virtual orbital, and vice versa, with matrix elements

Assuming orthonormality of the spin orbitals, the matrix elements are easily calculated

where both i and j refer to occupied orbitals and both a and b to virtual orbitals. The matrix is thus diagonal, and the inverse of χ 0 has the following 2 × 2 supermatrix representation

where ∆ǫ ia, jb = (ǫ aǫ i )δ i j δ ab , each block matrices being reindexed with the composite indices ia and jb.

In the RPA and RPAx approximations, the Bethe-Salpeter kernel of Eq. (A16) is approximated as the frequencyindependent Hartree(-Fock) form [Eqs. (19) and (20)]

where w ee (r 12 ) is a two-particle interaction, and ξ = 0 or ξ = 1 for RPA and RPAx, respectively. This kernel has the following supermatrix elements

where qr| ŵee |ps are the two-electron integrals. The supermatrix representation of the interacting polarization propagator χ λ is then found from the Bethe-Salpeter equation [Eq. (A17)] written in the spin-orbital basis where

or, in matrix form,

Using the well-known fact that, if A λ + B λ and A λ -B λ are positive definite, the non-Hermitian eigenvalue equation (B9) with real spin orbitals can be transformed into the following half-size symmetric eigenvalue equation

where

, and using the spectral decomposition M -1/2 λ = n ω -1 n,λ Z n,λ Z T n,λ , the correlation two-particle density matrix P c,λ can be expressed as

(B18)

2. Expressions for spin-restricted closed-shell calculations

For spin-restricted closed-shell calculations, the eigenvectors (X n,λ , Y n,λ ) can be transformed into spin-singlet excitation/diexcitation vectors

and spin-triplet excitation/diexcitation vectors

( 3,-1 y n,λ ) ia = Y n,λ i↓a↑ , (B20d)


Mohamed Ali Hajjaji
email: hajjaji@etu.u-bourgogne.fr

El-Bay Bourennane

Abdessalem Ben Abdelali

Abdellatif Mtibaa

Combining Haar Wavelet and Karhunen Loeve Transforms for Medical Images Watermarking

is paper presents a novel watermarking method, applied to the medical imaging domain, used to embed the patient's data into the corresponding image or set of images used for the diagnosis. e main objective behind the proposed technique is to perform the watermarking of the medical images in such a way that the three main attributes of the hidden information (i.e., imperceptibility, robustness, and integration rate) can be jointly ameliorated as much as possible. ese attributes determine the e ectiveness of the watermark, resistance to external attacks, and increase the integration rate. In order to improve the robustness, a combination of the characteristics of Discrete Wavelet and Karhunen Loeve Transforms is proposed. e Karhunen Loeve Transform is applied on the subblocks (sized 8 × 8) of the di erent wavelet coe cients (in the HL , LH , and HH subbands). In this manner, the watermark will be adapted according to the energy values of each of the Karhunen Loeve components, with the aim of ensuring a better watermark extraction under various types of attacks. For the correct identi cation of inserted data, the use of an Errors Correcting Code (ECC) mechanism is required for the check and, if possible, the correction of errors introduced into the inserted data. Concerning the enhancement of the imperceptibility factor, the main goal is to determine the optimal value of the visibility factor, which depends on several parameters of the DWT and the KLT transforms. As a rst step, a Fuzzy Inference System (FIS) has been set up and then applied to determine an initial visibility factor value. Several features extracted from the Cooccurrence matrix are used as an input to the FIS and used to determine an initial visibility factor for each block; these values are subsequently reweighted in function of the eigenvalues extracted from each subblock. Regarding the integration rate, the previous works insert one bit per coe cient. In our proposal, the integration of the data to be hidden is bits per coe cient so that we increase the integration rate by a factor of magnitude .

Introduction

Medical imaging is an important tool and essential in the diagnosis and decisions made by health professionals. In this context, several techniques and imaging models have been proposed by many authors. Among these techniques, the most common are MRI, Echographic images, Radiographic images, and Mammographic images. ese techniques have proven to be very successful in diagnosing, and so forth.

However, these diagnoses are o en insu cient or inconclusive due to the complexity of the diseases or the limitation of the imaging techniques themselves. erefore, developing new tools enable physicians, usually located in di erent regions or countries, to collaborate remotely, in order to get a better diagnosis which has become increasingly widespread and necessary. is trend is known as telemedicine. e main problem of these endeavors arises when managing the integrity and con dentiality of data on the internet against pirates. Several solutions based on the use of access control techniques exist, but they remain elusive and hence the appearance of watermarking techniques in order to ameliorate the security control of the network in which those medical images [ -] are shared [ ].

In the general case, digital watermarking is a technique that consists in hiding information in a digital document (indelible or invisible depending in the nature of the document and hence the name) to ensure security in the intended service (copyright, integrity, and nonrepudiation). A special feature of watermarking compared to other techniques is that the watermark is intimately linked to the associated document and resistant to attacks. erefore, the watermarking is theoretically independent of the le format and it can be detected or extracted even if the document has changed or is incomplete.

is paper is organized as follows. Section presents the problem and criteria for a digital watermarking system. Related work is provided in Section . Section deals with the proposed watermarking technique. Section presents a methodology for adjusting of the visibility factor. e insertion and extraction of the watermark are, respectively, represented in the Sections and . e simulation, validation, and the evaluation of the proposed algorithm are represented in the Section . Finally, Section concludes the paper.

Problem and Criteria for a Digital Watermarking System

Watermarking techniques require various features according to their areas of application and their underlying objectives. e hidden watermark in an image must satisfy three basic conditions: robustness, imperceptibility, and capacity. ey are not independent of each other; rather they are closely related [ ].

. . Imperceptibility.

e watermark should not a ect the quality of the original image a er any watermarking operation. Cox et al. [ ] de ne the imperceptibility as visual similarity between the original and watermarked images.

e watermark must be inserted in a way that is completely invisible to the Human Visual System (HVS) [ ].

e insertion process must not damage the host image; that is, the watermarked image has to be visually equivalent to the original image. Not only the image should not be distorted, but also the watermark must be hidden and invisible, otherwise it could be easily removed. In general, the more the imperceptibility is low, the more the robustness and capacity are strong.

. . Capacity.

e ability of a watermarking system refers to the ratio of the "Amount of data" to hide on the "size of the host document" [ ]. Sometimes the size of the watermark is limited just to bit. In general, the more the capacity is low, the more the robustness and the imperceptibility are strong.

. . Robustness. Robustness [ ] is the resistance of the watermark system against intentional transformations of the watermarked image.

ese transformations can be of geometric type (rotation, cropping); they can also include all types of degradation of the image frequencies (lossy compression, high pass lter, or low pass lter). In general, the more the robustness is low, the more the capacity and the imperceptibility are strong.

Many ongoing e orts on watermarking techniques are being carried out with the aim of optimizing these three aforementioned criteria. Moreover, these three parameters are closely related to a pattern of watermarking images so any modi cation of one of these factors in uences directly the others (Figure ).

Related Works

In this section we present an in-depth review of digital image watermarking techniques. It describes the previous works which had been done on digital watermarking by using DWT technique.

In [ ] Anuradha and Singh proposed a watermarking system aimed at protecting the copyright and the control of the integrity related to the digital products. e multiresolution space, obtained by Haar Wavelet Transform, is used for hiding the watermark in the rd level of decomposition. e horizontal, vertical, and diagonal subbands are processed for inserting the totality of the watermark.

In [ ] Kashyap proposed a blind watermarking system in the rd level Discrete Wavelet Transform. e idea consists in hiding a watermark in the low frequency subband. For insertion, Kashyap de ned the visibility factor in function of scaling factor of the subbands of the original image and the watermark. e insertion step follows the next formula:

WMI = × LL 3 + × WM 3 , ( )
where WMI is low frequency component of watermarked image, LL 3 is the low frequency component of the original image obtained by -level DWT, and WM 3 is the low frequency component of Watermark image. and are the scaling factors of the original image and the watermark, respectively.

In [ ] Peter Cika describes a new watermarking method based on two-dimensional wavelet transform and the Singular-Value Decomposition. e diagonal matrix (a er decomposition on SVD) is used for inserting the watermark a er decomposition on wavelet transform.

In [ ] Kannammal and Subha Rani proposed a watermarking system for medical images authentication. e proposed space, for hiding information, is the low frequency subband a er the rst level of the wavelet transform decomposition. e components chosen to hide the watermark are selected by the random numbers generated, which have the integers value from to . e watermark is the hash value of the MSBs of di erent coe cients selected to hide the watermark. e output of the hash function is embedded into the selected coe cients, and it is combined with the MSBs to get the watermarked coe cients.

In [ ] Rawat and Raman proposed a dual watermarking scheme based on Discrete Wavelet Transform (DWT), Wavelet Packet Transform (WPT) with best tree-and Singular-Value Decomposition (SVD). Using subsampling of the cover image, four subimages are obtained and they calculate their SVD values. ey chose and two subimages which contain the highest sum of SVD for embedding two di erent watermarks. e watermarks W and W are embedded with two di erent methods. For embedding the rst watermark, they apply the WPT in l level; a erwards the best tree is chosen. e insertion step consists in combining the two diagonals matrices of the SVD transforms of the watermark and the best tree. For the second watermark, the DWT in l level is applied. en they calculate the SVD of di erent subbands. For the insertion step, they applied the same method for embedding W .

In [ ] Yang et al. presented a watermarking algorithm based on the Integer Wavelet Transform (IWT). e lowhigh (LH) and high-low (HL) subbands are used for hiding the watermark. e simulations show good results in terms of the watermarked image quality. e proposed schema is robust against several attacks such as brightness, inversion, and compression attack.

In [ ] Makbol and Khoo presented a watermarking algorithm based on the Redundant Discrete Wavelet Transform (RDWT) and the Singular-Value Decomposition. e watermark to be inserted is a gray scale image. e insertion step is applied directly in the SVD components of the RDWT subbands of the host image.

In [ ] Latif proposed a new watermarking method based on the Parametric Slant-Hadamard Transform. To ameliorate the imperceptibility factor, the authors have exploited the masking characteristics of the Human visual System using Fuzzy Inference System. e experimental results show that the proposed technique has a high imperceptibility as well as a high robustness against several attacks.

For controlling the authentication, Al-Saif et al. [ ] propose a new watermarking method applied on the gray scale image. e proposed space is the Karhunen Loeve Transform.

e eigenvalues are used for hiding the watermark. Most methods found in the literature fail to optimize the compromise existing in the digital watermarking (robustness, imperceptibility, and capacity). Moreover, all previous works use a logo to be hidden in the original image which is not the case concerning medical imaging.

In this paper, we propose a new simple methodology applied on medical imaging.

e basic idea consists in preserving the visual quality of the original image to increase the rate related to data to be inserted and to enhance the robustness of our method against many attacks.

Proposed Watermarking Technique

In this section, we explain the proposed algorithm for embedding the totality of the patient's data in a medical image. Before the insertion, the patient's data undergoes many steps with the aim at increasing the level of integration and better ensuring its extraction a er applying di erent attacks. e insertion procedure is done by adding to the image frequency components, proper to the original image, data related to the patient.

e following expression explains the integration procedure of the watermark on the image frequency values:

ὔ ( ) = ( ) + × ( ) , ( ) with (i
) ( ) being the th old coe cient to support the watermark;

(ii) ὔ ( ) being the th new watermarked coe cient;

(iii) ( ) being the th bit to hide;

(iv) being named o en the visibility factor. e visibility factor is an important factor in the watermarking system. If is big, we win in terms of the robustness but we lose in terms of imperceptibility and vice versa.

In order to achieve an optimal imperceptibility, it is necessary to go through a stage of preparation of the insertion space. e rationale behind our contribution is based on the use of the subbands obtained a er applying the second level wavelet decomposition then followed by the subdivision stage in which they are further subdivided into smaller subblocks sized 8 × 8. A er that, we calculate the Karhunen Loeve Transform on each subblock and embed the watermark. Finally, the inverse transforms are performed to obtain watermarked image.

To evaluate the performance of our algorithm; an extraction phase is applied on the watermarked image. In addition, the Normalized Cross-Correlation (NC) is carried out to evaluate the similarity between the original and the extracted data. e proposed process is summarized in Figure .   . . Determining the Regions to Be Watermarked. In our method, we propose to use a combination of two transformations on the original image for embedding the watermark.

First, the second level Haar Wavelet Transform [ ] is carried out. In this stage, each subband is formed of subblocks of 8 × 8 coe cients. Second, we chose the frequency subbands, such as HH , HL , and LH . e principal characteristic of the high frequency subbands is that the edges and textures usually are found here more prominently.

ey are used with the aim to accord the visibility factor with psychovisual characteristics of these subbands. ird, we subdivide HH , HL , and LH into subblocks sized 8 × 8.

Finally, the Karhunen Loeve Transform is applied. Figure presents the steps for the preparation of the regions where the watermark will be hidden.

. . Preparation of the Watermark. Most of the algorithms, found in the literature for embedding watermarks, propose to integrate directly binary data in frequency components chosen in the original image, a er multiplication by a visibility factor. However, it is necessary in the medical imaging domain to take into account the integration rate given that the amount of information to be hidden in the host image (name, age, sex, diagnosis, and signature) can be signi cant. e idea consists in using the data from the patient (the watermark) as an index for insertion, a er its transformation into an octal sequence. is index services as an address to access a row among of each sub-block. Algorithm illustrates the typical example of the data from a patient and which can be inserted in the image.

e computation steps for the index table are performed as follows.

(i) For controlling the authenticity of the patient's record, a signature owned by the hospital center is generated. For this purpose, the SHA-hash function [ ] is used.

(ii) Convert the rest of data le into a binary sequence and concatenate it with the binary signature. us, we obtain a message which will be coded in the following steps.

(iii) For improving the reliability and the detection of the hidden data, we apply the concept of channel coding. So, the message is coded by the serial Turbocode [ ].

(iv) e whole of coded binary sequence is assumed to be equal to "3 × , " where is less or equal to 3 × . 3 × is the total number of subblocks in subbands (LH , HL , and HH ). e binary sequence is then converted into the octal representation format. We then generate three addressing tables of identical sizes equal to . We start by lling the rst table to contain addresses for indexing LH subband by the rst values of the octal presentation sequence. By the remaining octal indexing values, we ll the second table related to HL and the third table related to HH , respectively.

(v) For the insertion process, we start by the subblocks of 8 × 8 values each; extracted from HL , the rst values of the rst indexing table are selected for their uses during the insertion in this subband. A erwards, the same process is applied on the second values of the indexing table related to LH . e rest of indexing table data is then used for the insertion in the third subband HH .

Figure summarizes the di erent steps to follow to prepare the index table related to the original data to be inserted.

ese tables are subsequently used for indexing a pseudorandom sequence generated with a secret key for embedding and extracting of the watermark.

Adjustment of the Visibility Factor

e principal problem in the watermarking domain is the compromise between imperceptibility, robustness, and the integration rate. However, the problem of the integration rate can be solved by collecting the di erent data to be inserted, with the goal of minimizing the necessary modi cations on the image. Concerning the robustness and the imperceptibility, they are directly related to the visibility factor known as " . " erefore, it is very important to take into account this compromise. A good compromised can be achieved by weighting the visibility factor according to the Human Visual System HVS [ ]. One of the main contributions of this work consists in proposing a method in which data is inserted in the LH , HL , and HH subbands a er their decompositions into Karhunen Loeve components [ ]. e idea is to weight the value of alpha according to each subblock.

A er the multiresolution transform, the totalities of the subbands (HL , LH , and HH ) are divided into subblocks sized 8 × 8.

For the adjusting of the visibilities factors values, a twostep algorithm is applied. In the rst step, the visibility factors are de ned in function of contrast sensibility, entropy sensibility, and homogeneity factor, which are extracted from the subblocks of the di erent subbands (LH , HL , and HH ). Proportionally to the di erent subbands, these visibility factors are named, respectively, LH 2 , HL2 , and HH2 . In the next step, the sum of LH 2 , HL2 , and HH2 is used to determine the nal visibility factor named LH 2 Final , HL2 Final , and HH2 Final .

eir weights are based on the properties of the Karhunen Loeve space. Totality of the watermark with maximum size equal to 3 × P (P less or equal to 3 × N)

011011110• • • 100001001• • • 011010110• • • For indexing in the LH2 3 3 6• • • 4 1 1• • • 3 2 6• • • For indexing in the HL2
For indexing in the HH2

3 3 6• • • 4 1 1• • • 3 2 6• • • F : Preparation of the indexing tables.
e remainder of this section details how the visibilities factors are obtained, depending on LH , HL , and HH and on the eigenvalues extracted using the Karhunen Loeve Transform.

. . Fuzzy Inference System and Frequency Subbands for Determination of the Visibility Factor. e entropy sensitivity measures the degree of disorder on the image. It attains high values for a random texture; thus, the more the entropy is high, the more the reinforcement of the visibility factor is possible.

e contrast sensitivity is o en important when moving from a low gray level (or high) to a high gray level (or low).

is is the case in the transition regions in an image. In this case, the more the contrast sensitivity is high, the more the reinforcement of is possible.

e homogeneity sensitivity re ects the fact that we o en encounter or not the same pair of pixels separated by the translation " " (a er application of the cooccurrence matrix). Its value is high if the image is a gray uniform. erefore, when the homogeneity sensitivity is small, it is possible to reinforce the visibility factor . In the proposed method, we integrate an intelligent process which is able to automate the value of the visibility factor in terms of the psychovisual quality of the corresponding insertion space.

is process is used in the so-called Fuzzy Inference Systems: in such a system, the contrast sensitivity , homogeneity sensitivity , and the entropy sensitivity are taken as inputs for estimating the adaptive weight visibility factors .

Using the Fuzzy Inference System enables us to increase the visibility factor in the less sensitivity areas (High texture, important contrast, and small homogeneity), while at the same time decreasing the value of the visibility factor in more sensitive areas (important homogeneity, small contrast, and small entropy). For this purpose, the fuzzy logic [ ] is used. F : e steps for determining the primary visibility factor for each subband. In general, it is based on the idea of the human experts, by their subjective and qualitative descriptions of behavior of watermarking method with natural language. e principle of the fuzzy logic is similar to the human behavior. It is based on the linguistic variables related to the human language. Moreover these variables are determined by some empirically experiments.

As shown in Figure , a Fuzzy Inference System is composed of three principal blocks:

(i) a Fuzzi cation step converts the numeric values to degrees of membership of di erent fuzzy set;

(ii) an Inference engine step that contains di erent rules;

(iii) Defuzzi cation step to generate a net worth for which represents the output of the fuzzy system.

When applying fuzzy logic to image watermarking, it is very important to determine the essential elements to nd the optimal value for adjusting . ose elements are the fuzzy variables, inference rules, and the membership functions.

e Fuzzy variables, also called linguistic variables, do not only take binary values but also have an in nite number of possible values between the "logic true" and the "logic false. "

e Fuzzy variables are involved in the description of certain situations, phenomenon, or process generally containing fuzzy quali ers. For example, (as shown in Figure ) for the entropy sensitivity we use the following fuzzy variables: Low texture, Medium textured, and High texture.

Inference rules and fuzzy rules are used for linking the di erent variables of the fuzzy system with its input variables and fuzzy outputs. ese rules come in the following form: If (condition ) and/or condition ( ) then (action on the outputs).

We summarize that these rules make the experience of the expert and they are usually not uniquely de nable as each individual creates his own rules. To do this, we de ne two notions:

(i) the membership functions that de ne the degree of truth of fuzzy variable depending on the input;

(ii) the fuzzy intervals which determine the number of fuzzy variables.

e input and output membership functions exploited are shown in Figure . It must be noted that this approach enables us to adjust the entropy (or homogeneity and the contrast) membership function in such a manner that best ts to the properties of the image.

In consequence, the approximations of the inferred values are optimized and used to generate an adaptive value strength for the inserted watermark. e membership functions used in our algorithm are the triangular and trapezoidal functions.

Concerning the evaluation of the output system, in fuzzy logic, the Defuzzi cation phase is used for translating the fuzzy values into numerical values.

is step is done by using the membership functions. In our approach, the inference results are subsequently computed by means of the minimum-maximum Defuzzi cation method. In this manner, we determine the initial visibility factors noted LH 2 , HL2 , and HH2 .

. . Impact of the Karhunen Loeve Transform on the Visibility Factor. Given an image (in our case sized 8 × 8) formed by columns of rows each. Let us call those columns for = {1, 2, . . . , 8}; we calculate the covariance matrix of the image [ ]:

= - - ,
indicates the matrix transpose.

( ) Since is real and symmetric, it is always possible to nd a set of orthonormal eigenvectors. Let V and , = {1, 2, . . . , 8}, be the eigenvectors and corresponding eigenvalues of , arranged in decreasing order so that ≥ +1 for = {1, 2, . . . , 8}. Let be a matrix whose rows are formed of the eigenvectors of , arranged so that the rst row of is the eigenvector corresponding to the largest eigenvalue, and the last row is the eigenvector corresponding to the smallest eigenvalue. e Karhunen Loeve Transform (also known as the Principal Components Transform) and its inverse (IKLT) may then be de ned as KLT = × , IKLT = × KLT, with : eigenvector matrix extracted from : eigenvector transpose matrix.

( ) is operation has a few key features. First, the Karhunen Loeve Transform decorrelates the signal components of KLT, suggesting that we could reconstruct each Karhunen Loeve component separately in the Karhunen Loeve domain as a sequence of independent reconstructions. Second, the Karhunen Loeve Transform tends to compact the original block content with the eigenvectors stemming of the eigenvalues. is advantage will be used for the weighting again of the visibility factors (already de ned in terms of multiresolution space).

At the Karhunen Loeve Transform, di erent eigenvectors , for = {1, 2, . . . , 8}, present the main directions that carry energy in the image.

However, the eigenvectors are closely related to the eigenvalues our idea is to weight the value of alpha in function to di erent eigenvalues .

Whether ≥ +1 for = {1, 2, . . . , 8}, in this case the value of is inversely proportional to values; that is, ≤ +1 ( = {1, 2, . . . , 8}). In this case if, for example, the insertion will take place on 5 , the value of is equal to 5 .

In conclusion, a er determining by the psychovisual characteristics, it will be reweighted (depending on the column just before the calculation of the Karhunen Loeve transform matrix) for forming the nal value of the visibility factor . .

Proposed Insertion Algorithm

= V 1 V 9 V 17 V 25 V 33 V 41 V 49 V 57 V 2 V 10 V 18 V 26 V 34 V 42 V 50 V 58 V 3 V 11 V 19 V 27 V 35 V 43 V 51 V 59 V 4 V 12 V 20 V 28 V 36 V 44 V 52 V 60 V 5 V 13 V 21 V 29 V 37 V 45 V 53 V 61 V 6 V 14 V 22 V 30 V 38 V 46 V 54 V 62 V 7 V 15 V 23 V 31 V 39 V 47 V 55 V 63 V 8 V 16 V 24 V 32 V 40 V 48 V 56 V

( )

A er performing the matrix multiplication for × , the components KLT 1 , KLT 2 , . . . , KLT 63 , and KLT 64 are given by

KLT 1 = V 1 × 1 + V 9 × 9 + V 17 × 17 + V 25 × 25 + V 33 × 33 + V 41 × 41 + V 49 × 49 + V 57 × 57 , KLT 2 = V 2 × 2 + V 10 × 10 + V 18 × 18 + V 26 × 25 + V 34 × 34 + V 42 × 42 + V 50 × 50 + V 58 × 58 , . . . KLT 63 = V 7 × 7 + V 15 × 15 + V 23 × 23 + V 31 × 31 + V 39 × 39 + V 47 × 47 + V 55 × 55 + V 63 × 63 , KLT 64 = V 8 × 8 + V 16 × 16 + V 24 × 24 + V 32 × 32 + V 40 × 40 + V 48 × 48 + V 56 × 56 + V 64 × 64 .
( )

e
BioMed Research International = V 1 V 9 V 17 V 25 V 33 V 41 V 49 V 57 V 2 V 10 V 18 V 26 V 34 V 42 V 50 V 58 V 3 V 11 V 19 V 27 V 35 V 43 V 51 V 59 V 4 V 12 V 20 V 28 V 36 V 44 V 52 V 60 V 5 V 13 V 21 V 29 V 37 V 45 V 53 V 61 V 6 V 14 V 22 V 30 V 38 V 46 V 54 V 62 V 7 V 15 V 23 V 31 V 39 V 47 V 55 V 63 V 8 V 16 V 24 V 32 V 40 V 48 V 56 V 64 × KLT 1 KLT 9 KLT
IKLT 1 = V 1 × KLT 1 + V 2 × KLT 2 + V 3 × KLT 3 + V 4 × KLT 4 + V 5 × KLT 5 + V 6 × KLT 6 + V 7 × KLT 7 + V 8 × KLT 8 , IKLT 2 = V 9 × KLT 1 + V 10 × KLT 2 + V 11 × KLT 3 + V 12 × KLT 4 + V 13 × KLT 5 + V 14 × KLT 6 + V 15 × KLT 7 + V 16 × KLT 8 , . . . IKLT 63 = V 49 × KLT 57 + V 50 × KLT 58 + V 51 × KLT 59 + V 52 × KLT 60 + V 53 × KLT 61 + V 54 × KLT 62 + V 55 × KLT 63 + V 56 × KLT 64 , IKLT 64 = V 57 × KLT 57 + V 58 × KLT 58 + V 59 × KLT 59 + V 60 × KLT 60 + V 61 × KLT 61 + V 62 × KLT 62 + V 63 × KLT 63 + V 64 × KLT 64 .
( )

If we modify the eight coe cients in ({V \ 1 ≤ ≤ 8}), the distortions will be spread on the totality of the KLT ({KLT \ 1 ≤ ≤ 64}) matrix. is means that the underlying image will be modi ed. When the IKLT matrix is computed, signi cant modi cations in the components of the rst row are produced. is feature enables us to easily identify the exact location of the watermark. erefore, it can be concluded that a change in ({ \ 1 ≤ ≤ 8}) produces a noticeable (and thus detectable) change in the rst column of the IKLT matrix.

. . Embedding Steps. Our approach is based on the combination of Haar Wavelet and Karhunen Loeve Transforms for hiding the watermark. e embedding process is described as follows.

(i) By using the Haar Wavelet Transform, we decompose the original image into second level subband. A er that, LH , HL , and HH are extracted and decomposed into subblocks sized 8 × 8.

(ii) e Karhunen Loeve Transform is applied on the di erent subblocks.

(iii) Using a secret key, we generate a pseudorandom sequence Key (key \ 1 ≤ ≤ 8) in which each number can take a value of either or . en this sequence is multiplied by the visibility factor, related to the subblock to be watermarked.

(iv) With the th block to be watermarked, the th value of the index table, and the binary sequence generated by the key (a er multiplication by the nal value of ), the embedding process can be initiated. Depending on the value " " ({ \ (0 to 7}) of the index table, the values of the ( th + ) column of eigenvector ({ \ (1 to 8}) are combined with Key multiplied by . is step will be stopped a er we nish the values of di erent indexing tables.

Proposed Extraction Algorithm

Generally, the extraction phase follows the reverse steps with respect to the insertion. We have the original image " " and the watermarked image " . " (i) e second Haar Wavelet Transform is applied on the original and watermarked images.

(ii) Extract the subbands of the second decomposition LH , HL , and HH . (iii) Decompose the di erent subbands into subblocks sized 8 × 8. (iv) Compute the di erence, all the blocks (two by two) between subblocks own to I and those related to . (v) In the proposed method, the inserted watermark is the binary sequence generated by a secret key a er their multiplying by the visibility factor. For extraction, the resulting di erence 8 × 8 matrices ( ) are compared, term by term, with threshold noted ( is determined empirically):

BioMed Research International

if , ≥ ⇒ , = 1 else , ≺ ⇒ , = 0,
, \ 1 to 8 .

( ) (vi) Calculate the correlation between the secret key and all row of . e number of the column containing the maximum correlation value indicates the extracted index value ({index value \ 0 ≤ ≤ 7}).

Figure

illustrates the di erent steps to follow for extracting the index table related to LH (the same step is applied on the watermarked subbands HL and HH ).

(vii) Concatenate all indexing tables extracted from di erent subbands, and convert them to binary data.

(viii) Apply the decoder algorithm (serial Turbocode) [ ].

Finally, the patient data is extracted, veri ed, and eventually corrected.

Validation of the Proposed Method

For medical images, the watermark must be imperceptible.

e watermarked image should be widely similar to the original image so as not leading to a misdiagnosis.

e validity of any watermarking algorithm can become more important than testing it against various attacks types.

For this, we subject the watermarked medical image to a series of attacks and test the sensitivity of the watermark and its ability to detect any change in the image. A er application of each attack, the entirety of embedded watermark is extracted and compared through similarity analysis with the original marks ( original , extracted ) to ensure that these marks are not damaged by the attacks applied on the image.

. . Watermark Detection Tools.

e measure of "degree of reliability" of the detected watermark is accomplished by the "calculation of distances" between the inserted and detected watermark. is measure is carried out using the Normalized Cross-Correlation [ ]. e Normalized Cross-Correlation (NC) of two signals consists in computing their dependence.

e NC is de ned as

NC , = ∑ ∑ - - ∑ ∑ - 2 ∑ ∑ - 2 , ( )
where , ̄ indicate, respectively, the original and the extracted watermark and and ̄ correspond, respectively, to the mean of the original and extracted watermark. In the literature, a NC value which is equal or above . denotes an acceptable extracted watermark [ ].

. . Peak Signal to Noise Ratio (PSNR). Among the most important distorting measures in image processing is the Peak Signal to Noise Ratio, PSNR [ ]. It is an assessment of the decibel di erence between the original image and one that is processed. In fact, a PSNR below dB image can be considered useless. e PSNR is de ned by the following formula:

(PSNR) dB = 10 log 10 × max 2 , ∑ , , -,

where and are the number of rows and columns of the image which contains × pixels, is the host image, and is the watermarked image. . . Experimental Results. Regardless of the domain or the method for hiding of the watermark, it is very important to have a good PSNR and WPSNR values; this to di erent Radiographic images sized of 512 × 512 pixels of resolution of bits/pixel. We start analysing our experimental results by a preliminary study, that is to say, a study without the application of attacks on the watermarked images. Figure shows the original medical images and watermarked ones. We notice that the Human Visual System does not distinguish the di erence caused by the marking.

e rst test of robustness for an attack is the application of the JPEG compression. It must be noted that the image compression algorithms are particularly aggressive for watermarked images. We have chosen to apply di erent image compression ratios to the watermarked images as shown in Table which presents the results of simulations showing the NC values between the original and extracted watermark a er image compression attacks.

A er applying the JPEG image compression attacks, we remark that even when varying the compression rate factor between % and %, the NC factor remains equal to . We conclude that the proposed approach makes the inserted watermarks resistant to this type of attacks.

e second kind of attacks tests is the application of several types of digital lters. In our experiments, we applied median, wiener, and low pass digital ltering attacks. We notice that our proposed method is very e ective against these types of attacks (NC is equal to regardless of the size of the lter).

e third attack is the application of two types of noise: Gaussian noise and Salt & Pepper noise. Tables and illustrate the PSNR, WPSNR, and the NC before watermarking and a er applying these attacks.

In general, it is necessary to test our schema against noise. e applied attacks are Gaussian noise (with di erent variance factors) and Salt & Pepper noise attack (varying the density factors). e di erent tests show that the proposed method attains good results, with NC = 1 for most of the cases.

e fourth type of attack applied is the geometric transform such as the rotation attack and cropping. Among the most dangerous attacks applied on a several watermarking algorithms are the geometric transforms. e proposed tests are the cropping attack (with di erent windows) and the rotation attack. e obtained results give us an NC equal or close to .

. . Evaluation of the Proposed Algorithm. In this section, the proposed watermarking method is investigated by comparing our results to those cited in the subsection of the related work.

Comparing the psychovisual quality of the original image and the watermarked image, the proposed algorithm yields very good results. In the absence of attacks, the PSNR equal to .

and WPSNR equal to . are obtained, yielding results that are approximately equal or o en better than those algorithms previously cited in the related section of the works.

A er applying many attacks, it is necessary to evaluate the Normalized Cross-Correlation (NC). Among the most serious attacks is the attack by image compression such as JPEG . In present paper, our applied algorithm is very e ective against this kind of attacks. e NC value stills equal to when the rate of image compression goes from % to %. is obtained result is more accurate than all results quoted precisely where the NC decreases with the increase ). e obtained results of the NC are equal to .

To evaluate the proposed method against the noise, we attacked watermarked images by two types of noises such as the Salt & Pepper and Gaussian noise. We obtained very promising results; the NC is always equal to .

We also tested our method against geometric transforms attacks (cropping and rotation).

e obtained NC values are very close to . Compared to the previous works, our proposed method gives results near to those found in the literature.

Conclusion

e present work is a new robust watermarking algorithm combining the Haar Wavelet and the Karhunen Loeve Transforms. e main contribution of this paper consists in improving the three principal factors existing in all watermarking systems (robustness, imperceptibility, and integration rate). To do so, we came across many steps.

In order to improve the factor of imperceptibility, we used the high frequency (second subband of the Haar Wavelet Transform) to hide the watermark. e Fuzzy Inference System is used to determine the visibility factor according to the proper characteristics of the insertion plan. Each subband (LH , HL , and HH ) is subdivided into subblocks sized 8 × 8. en, the Karhunen Loeve Transform is applied in order to decorrelate the di erent wavelets coe cients. Indexing tables are then used to choose the location of components supporting the watermark. At this stage, the visibility factors determined by the FIS ( LH 2 , HL2 , and HH2 ) are adapted according to the weights of the eigenvalues to determine the nal visibility factors ( LH 2 Final , HL2 Final , and HH2 Final ).

Our contribution concerning the robustness is the use of the ECC by means of the serial Turbocode. We obtained good results in terms of the extracted watermark which is similar to the original. Our principal idea to increase the integration rate by a factor consists in inserting the octal representation of the watermark.

To evaluate the performance of our method, the proposed system is applied on medical images. Several tests are performed, such as digital ltering, JPEG compression, adding noise, and geometric transformation. e results show that our method is very robust against these attacks. It supports image compression attacks such as JPEG up to % compression ratio.

Our experiments have also shown that our method is resistant to digital ltering attacks. We noticed that the extracted watermark is similar to the original watermark.

To evaluate the resistance of our proposal against the geometric transformation attacks, image rotation techniques and cropping are applied on watermarked images. Here again, the watermark extraction was faithful.

For

  Eigenvalues: i : i = 1 : 8 Eigenvectors: V i : i = 1 : 8

  Figure describes the general steps to follow for setting the visibilities factors values.

  determination of the visibility factor. the procedure for determining the primary values of the optimal visibility factors related to the visual characteristics of the di erent subbands.In the HVS, among the most important factors, three properties are employed in the watermarking algorithm: contrast sensitivity, entropy sensitivity, and homogeneity sensitivity. ese factors are extracted by Cooccurrence matrix (CM) [ ].

  (v) Perform the Inverse Karhunen Loeve Transform. (vi) Perform the Inverse Haar Wavelet Transform to obtain the watermarked image.(vii) Display watermarked image.

Figure

  Figure shows the insertion algorithm applied on the second subband LH (the same steps are applied on the HL and HH ).Figure shows the insertion step in the matrix 8 × 8 on the Karhunen Loeve component.

  Figure shows the insertion algorithm applied on the second subband LH (the same steps are applied on the HL and HH ).Figure shows the insertion step in the matrix 8 × 8 on the Karhunen Loeve component.

  Functions and mapping of their input/output variables to fuzzy sets.

F:

  Insertion step in the matrix 8×8 applied on the Karhunen Loeve components.

  . . Weighted Peak Signal to Noise Ratio (WPSNR).e Peak Signal to Noise Ratio PSNR is based on comparing pixel to pixel the original image and the received watermarking image. e WPSNR proposed by Voloshy Noviskiand and al.[ ] is de ned by the following formula: , ) representing the local variance of pixel at location ( , ), ( , ) is the intensity value of the pixel ( , ) in the original image, and ( , ) the intensity value of the pixel in the image under test. and are, respectively, the height and width of the image.

  Choice of Components to Support the Watermark. Suppose that the input subblock, sized 8 × 8, is represented by a matrix . e Karhunen Loeve Transform can be represented by KLT = × , KLT 1 KLT 9 KLT 17 KLT 25 KLT 33 KLT 41 KLT 49 KLT 57 KLT 2 KLT 10 KLT 18 KLT 26 KLT 34 KLT 42 KLT 50 KLT 58 KLT 3 KLT 11 KLT 19 KLT 27 KLT 35 KLT 43 KLT 51 KLT 59 KLT 4 KLT 12 KLT 20 KLT 28 KLT 36 KLT 44 KLT 52 KLT 60 KLT 5 KLT 13 KLT 21 KLT 29 KLT 37 KLT 45 KLT 53 KLT 61 KLT 6 KLT 14 KLT 22 KLT 30 KLT 38 KLT 46 KLT 54 KLT 62 KLT 7 KLT 15 KLT 23 KLT 31 KLT 39 KLT 47 KLT 55 KLT 63 KLT 8 KLT 16 KLT 24 KLT 32 KLT 40 KLT 48 KLT 56 KLT 64

. .

  IKLT 9 IKLT 17 IKLT 25 IKLT 33 IKLT 41 IKLT 49 IKLT 57 IKLT 2 IKLT 10 IKLT 18 IKLT 26 IKLT 34 IKLT 42 IKLT 50 IKLT 58 IKLT 3 IKLT 11 IKLT 19 IKLT 27 IKLT 35 IKLT 43 IKLT 51 IKLT 59 IKLT 4 IKLT 12 IKLT 20 IKLT 28 IKLT 36 IKLT 44 IKLT 52 IKLT 60 IKLT 5 IKLT 13 IKLT 21 IKLT 29 IKLT 37 IKLT 45 IKLT 53 IKLT 61 IKLT 6 IKLT 14 IKLT 22 IKLT 30 IKLT 38 IKLT 46 IKLT 54 IKLT 62 IKLT 7 IKLT 15 IKLT 23 IKLT 31 IKLT 39 IKLT 47 IKLT 55 IKLT 63 IKLT 8 IKLT 16 IKLT 24 IKLT 32 IKLT 40 IKLT 48 IKLT 56 IKLT 64

		Inverse Karhunen Loeve Transform can be repre-
		sented as
	IKLT =	× KLT,
	IKLT 1

  17 KLT 25 KLT 33 KLT 41 KLT 49 KLT 57 KLT 2 KLT 10 KLT 18 KLT 26 KLT 34 KLT 42 KLT 50 KLT 58 KLT 3 KLT 11 KLT 19 KLT 27 KLT 35 KLT 43 KLT 51 KLT 59 KLT 4 KLT 12 KLT 20 KLT 28 KLT 36 KLT 44 KLT 52 KLT 60 KLT 5 KLT 13 KLT 21 KLT 29 KLT 37 KLT 45 KLT 53 KLT 61 KLT 6 KLT 14 KLT 22 KLT 30 KLT 38 KLT 46 KLT 54 KLT 62 KLT 7 KLT 15 KLT 23 KLT 31 KLT 39 KLT 47 KLT 55 KLT 63 KLT 8 KLT 16 KLT 24 KLT 32 KLT 40 KLT 48 KLT 56 KLT 64 , IKLT 2 , . . . , IKLT 63 , and IKLT 64 are given by

		.
		( )
	A er performing the matrix multiplication for KLT = × × , the pixels IKLT 1	×

  is especially true in the medical imaging domain. is work has been applied

		Extracted	LH2	Original	LH2
					Key () is 8 bits
					in length
		Redivision into subblock sized 8 × 8
	Fragmentation into	images sized 8 * 8			00011010
			Extraction algorithm
			1
					Extracted
				1	indexing table
		Conversion the binary	own to HL2
		data to octal
				100001001
			F	: Extraction Step.

  Table shows PSNR, WPSNR, and the Normalized Cross-Correlation (NC) factor of Radiographic watermarked image robustness of our watermarking schema against these attacks with di erent window sizes of the lters.

	T	: Evaluation of the PSNR, WPSNR, and the Normalized
	Cross-Correlation values of the watermarked and attacked images
	by a JPEG	image compression Algorithm.
	Rate (%)	PSNR (dB)	WPSNR (dB)	NC
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.

  Table shows the PSNR, WPSNR, and the NC values a er applying the rotation attack according to the rotation angle. Table shows the PSNR, WPSNR, and the NC values a er applying the cropping attack with various sized window.

  PSNR, WPSNR, and Normalized Cross-Correlation for watermarked and attacked by Gaussian noise.

	Variance	PSNR (dB)	WPSNR (dB)	NC
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.
	.	.	.	.

T

:

  : PSNR, WPSNR, and Normalized Cross-Correlation for watermarked and attacked by image rotation transform.

	Rotation angle ( ∘ )	PSNR (dB)	WPSNR (dB)	NC
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
			.	.	.
	T	: PSNR, WPSNR, and Normalized Cross-Correlation for
	watermarked and attacked by "cropping. "
	Window size [ × ] [ × ] [ × ] [ × ] [ × ]	PSNR (dB) . . . . .	WPSNR (dB) . . . . .	NC . . . . .

T

Acknowledgment

is work has been supported by the Pharmimage Project by Regional Council of Burgundy, France.

LH2

KLT inverse for J = 0 to 7 V(j, k) = V(j, k) + kLH2 Final × key(j)

Conflict of Interests

e authors declare that there is no con ict of interests regarding the publication of this paper.