Markus Pernpointner
email: markus.pernpointner@pci.uni-heidelberg.de

The four-component two-particle propagator for the calculation of double ionization spectra of heavy-element compounds I. Method

In this work we present a new four-component implementation of the two-particle (2p) propagator as a powerful tool for calculating double ionization spectra of systems containing heavy elements.

The 2p-propagator approach provides immediate access to the complete final state distribution of doubly ionized systems together with a detailed configurational information of each final state. This information is vital for the analysis of electronic decay processes as radiationless relaxation mechanisms after an initial core hole ionization and for the simulation of Auger spectra. The capabilities of the four-component realization is first demonstrated at the xenon atom and the bromine molecule for which experimental results are available. It turns out that a relativistic treatment is required for an adequate description of the experimental findings. In this first paper we focus on the electronic structure described by the 2p-propagator and will later on apply this method to the calculation of Auger spectra and electronic decay processes in heavy systems.

I. INTRODUCTION

In the last few years modern experimental techniques such as double charge transfer spectroscopy [1][2][3], time-of-flight mass spectrometry [4], photoion and photoelectron coincidence spectroscopy [5][6][7] and threshold photoelectron coincidence (TPEsCO) [8] allowed for a detailed analysis of double ionization processes and a characterization of the resulting fragments. In these processes atoms, molecules or clusters are irradiated by photons of sufficient energy to produce at least two photoelectrons and dications that undergo fragmentation in most cases. The analysis of the photoelectrons and fragments yields a detailed insight into the ionization processes from an experimental viewpoint and is of prominent interest for the understanding of atmospheric chemistry, for example, where highly energetic radiation is initiating a multitude of chemical reactions. Of similar importance is the understanding of radiation damage in biological tissue also initiated by highly energetic primary ionization processes. In the case of CF 4 , for example, the double ionization threshold was determined as 37.5 eV [8] and by photoion coincidence measurements the corresponding fragmentation products such as CF + 3 + F + (37.6 eV) and CF + 2 + F + (42.4 eV) could be identified. Simultaneous theoretical calculations with nonrelativistic Green's function techniques in combination with a two-hole (2h) population analysis [9] allowed for an accurate interpretation of the obtained experimental CF 4 spectrum. Other prominent examples for successful theoretical prediction and experimental analysis are the double ionization processes in formaldehyde [10], the interpretation of the spectrum in SF 6 [11], N 2 O [12] and O 2 [13].

For all these systems the application of propagator methods turned out to be beneficial due to the accessibility of all dicationic states and their composition after a single matrix diagonalization. For the CF 4 molecule a final state analysis showed a prevailing hole localization on the fluorine atoms with only minor carbon contributions. Coincidence maps further allow for a classification in direct and indirect ionization processes depending on the photon energy as has been shown for SF 6 [11]. Here the indirect mechanism via a primary ionization of the F 2s level followed by an electronic decay could be established. These decay mechanisms are accurately predicted by the current theoretical methods and play a major role in the correct interpretation of the spectra.

Recently, experimentalists started to analyze double ionization spectra of heavier systems such as molecular bromine [14,[START_REF] Yencha | Proceedings of the XXIV ICPEAC Conference[END_REF] and iodine [START_REF] Edvardsson | [END_REF], the interhalogen compound ICl [17] and the iodine-containing systems HI, CH 3 I and CF 3 I [18]. For the doubly ionized final states of ICl state-specific multi-reference configuration interaction (MRCI) methods were employed by Edvardsson et al. [17] where the X 3 Σ - 0,1 , A 1 ∆ 2 and B 1 Σ + 0 states could be clearly identified in the outermost part of the spectrum between 27 eV and 29 eV. Above 30 eV the assignments were done tentatively. The potential curves give precise information about the boundedness of a specific state but it is more difficult to decide about further electronic decay processes and the nature of states at higher energy ranges.

Another prominent example for double ionization processes in systems with heavy elements are the methyl iodide CH 3 I and trifluoro methyl iodide CF 3 I spectra which were measured by the TOF-PEPECO (time-of-flight photoelectron coincidence) method [18].

Whereas in the methyl iodide spectrum four well resolved peaks of symmetry 3 A 2 (spinorbit split in an A 1 and E state), 1 E and 1 A 1 occur, the CF 3 I 2+ shows a more complicated structure with five peaks including spin-orbit splitting and a broad ionisation edge. The authors in [18] attribute the extra peaks to final states where one hole is each located on the iodine and fluorine atom leading to a dication with considerably lowered Coulomb repulsion compared to the purely iodine-centered dicationic state. Doubly ionized clusters often result from an electronic decay process following a primary ionization on one site. Hereby the Intermolecular Coulombic Decay (ICD) [19] or the Electron Transfer Mediated Decay (ETMD) [20] are the most prominent examples of these decay processes currently investigated. These processes are normally initiated by a inner valence ionization where the vacancy is filled by an outer valence electron from the same site (ICD) or from a neighboring site (ETMD). The released energy is hereby sufficient to ionize the system further where in general many decay channels are open. It sould be kept in mind that the accessibility and characterization of the possible final states will also depend on the inclusion of relativistic effects.

In order to treat the abovementioned systems containing heavy elements adequately a four-component implementation of the two-particle propagator was undertaken and is presented in this work. Hereby access to the full spectrum together with a detailed configuration information of the final states is obtained. Relativity will have substantial influence on the energetics of the resulting final states which in turn determine the structure of the spectra and the accessibility of decay pathways. In the following we first outline the relativistic methodology of the 2h propagator together with a focus on the symmetry handling in the double group-based framework. Afterwards, the double ionization spectra of xenon and Br 2 are calculated in order to demonstrate the capabilities of the method. For Br 2 double group configuration interaction (CI) and multi-reference coupled cluster (MRCC) calculations are also available [21] and allow for a comparison with our results. In subsequent work we will address relativistic effects in the double ionization spectra of CH 3 I, CF 3 I, I 2 and ICl and their impact on ICD/ETMD processes in clusters containing heavy elements.

II. METHODOLOGY IN THE RELATIVISTIC CASE

The two-particle propagator or Green's function describing a double ionization or attachment process reads as

G pq,p q (t 1 , t 2 ; t 1 , t 2 ) = -Ψ N 0 | T c p (t 1 )c q (t 2 )c † q (t 2 )c † p (t 1 ) |Ψ N 0 , (1)
with |Ψ N 0 being the exact N -particle ground state wave function of the Dirac-Coulomb

(DC) operator ĤDC = N i=1 c α i • p i + β i m e c 2 + V ext (i)1 4 + N i<j 1 r ij , (2)
where magnetic interaction (Gaunt) corrections to the Coulomb term and retardation effects are neglected. T is the Wick time-ordering operator and the c † q (t) [c q (t)] denote creation (destruction) operators for one-particle states |φ(q) in the Heisenberg representation. In the language of second quantization in which the algorithm is cast no excitations from negative energy states are considered (see [22] for the no-pair formalism in combination with propagators and [23] for an analysis of the general unboundedness problem in case of the DC operator).

For the practical implementation one starts form the energy-dependent representation of (1) obtained by a Fourier transformation: Π pq,p q (ω) = Π + pq,p q (ω) + Π - pq,p q (ω) .

Π + pq,p q (ω) hereby describes a two-particle attachment whereas Π - pq,p q (ω) is used for a twoparticle detachment (ionization) process. The Lehmann representation [24] of Π - pq,p q (ω) which we are focusing on in this work reads as

Π - pq,p q (ω) = - m x (m) pq x (m) * p q ω + E N -2 m -E 0 N -iη (3)
utilizing a complete set of exact N -2 particle states |Ψ N -2 m with corresponding energies

E N -2 m
. η is an infinitesimal positive number required for the convergence of the backtransformation. The spectral amplitudes

x (m) pq = Ψ N -2 m |c p c q |Ψ N 0
provide information about the composition of the mth final state with respect to the oneparticle states p and q. It should be pointed out that the N -particle ground state wave function |Ψ N 0 can be expanded in a N -particle basis of Slater determinants made of fourcomponent one-particle functions. The no-pair confinement puts us into the advantageous situation that the spin-orbital based expressions derived from diagrammatic perturbation theory need not be modified for a corresponding realization with four-component one-particle functions.

Up to this stage we have not yet obtained actual expressions for the propagator. This is established by a combination of the diagrammatic representation for Π -together with an algebraically formulated perturbation expansion of a nondiagonal representation of (3) termed as Algebraic Diagrammatic Construction (ADC) [25][26][27][28]. In a compact matrix notation (3) reads as

Π -(ω) = X † (ω1 + Ω) -1 X
with 1 being the identity matrix and Ω being the diagonal matrix of double ionization potentials. This diagonal representation is now brought into a nondiagonal form by inserting a so-called intermediate state basis [29] leading to

Π -(ω) = F † (ω1 -Γ) -1 F .
A perturbational expansion of the Γ and F matrices then allows for an order-by-order comparison to the graphs obtained by a diagrammatic perturbation expansion from which the explicit expressions for the matrix elements can be derived [25,30]. The F matrices are needed for the calculation of transition moments and are disregarded for the current inves-

V abcd = φ * a (1)φ * b (2) 1 r 12 φ c (1)φ d (2) ≡ ab|cd (4)
and their totally antisymmetric versions used in the implementation

V ab[cd] = V abcd -V abdc ≡ ab||cd . (5)
The V ab[cd] integrals exhibit the relations

V ab[cd] = -V ab[dc] = -V ba[cd] = V ba[dc] and V cd[ab] = V * ab[cd] . (6)
The indices a, b, c, d hereby cover the occupied and/or the virtual orbital space. The totally antisymmetric MO integrals together with the (relativistic) spinor energies a are the only entities entering the expressions for the DC-ADC matrix elements which were originally derived in the nonrelativistic regime by Schirmer and Barth [25] and restated by Tarantelli [28] with an appropriate sign change in order to obtain positive eigenvalues of the ADC matrix. In [28] both the spin-containing and spin-free formulas can be found where in the latter case two different sets of equations for the singlet and triplet final states result.

This classification according to spin multiplicity is no longer applicable for heavy systems with large spin-orbit coupling and the spin-orbital-based expressions are to be employed in the four-component case. It should be kept in mind that the AO/MO transformation is done exclusively in the positive energy space utilizing the corresponding Dirac-Hartree-Fock spinors. As a consequence, the molecular spinors can not be further relaxed with respect to the negative energy solutions in the ADC step which might yield to changes in the ionization potentials for very heavy systems and core-like final states. In the valence region these additional contributions are assumed to be small and the electron relaxation and correlation effects are completely covered by the propagator. In the following we will often refer to the equations given in [28] and we will therefore repeat them here for convenience in eqs. ( 7) -( 9). The three equations comprise the 2h/2h main block C ij,kl (7), the 3h1p/2h coupling block C rklm,ij (8) and the 3h1p/3h1p satellite block C rijk,slmn (9) of the DC-ADC matrix. In the following the indices i, j, k, . . . refer to occupied (hole) orbitals, the indices p, q, r, . . . to virtual (particle) orbitals and ijk... is an abbreviation for i

+ j + k + • • •. C ij,kl = V ij[kl] + δ jl 1 2 mrs V rs[km] V im[rs] ( rs -km )( rs -im ) rs -m - ik 2 (7)
-{i ↔ j} -{k ↔ l} + {i ↔ j, k ↔ l} - 1 2 rs V rs[kl] V ij[rs] ( rs -kl )( rs -ij ) rs - ijkl 2 -δ ik δ jl ij . C rklm,ij = δ ik V lm[jr] + δ il V mk[jr] + δ im V kl[jr] -(i ↔ j) . (8)
C rijk,slmn = δ rs δ il V jk[mn] + δ im V jk[nl] + δ in V jk[lm] (9)
+ {(i, j, k) → (j, k, i)} + {(i, j, k) → (k, i, j)} + δ jm δ kn V is[rl] + δ jl δ km V is[rn] -δ jl δ kn V is[rm] + {(i, j, k) → (k, i, j)} -{i ↔ j} +δ rs δ il δ jm δ kn ( r -ijk ) .
The abbreviation {i ↔ j} in ( 7) -( 9) hereby symbolizes an i, j swapping of the preceding expression cast in curly brackets and (i, j, k) → (k, i, j) stands for the index replacements as indicated. For the numerical realization the 2-and 4-tuples ij, . . . and rijk, . . . are to be understood as multi-indices addressing rows and columns of the DC-ADC matrix in a unique way. The expressions ( 7) -( 9) are complete up to second order in perturbation theory (DC-ADC2) and contain additional first order contributions in the satellite block not present in the strict second order formulation.

A fundamental property of all matrix elements above is their antisymmetry with respect to hole-index interchange within the bra or ket space (not between these spaces, of course).

In equation (8), for example, we end up with the identities C rklm,ij = -C rlkm,ij = -C rklm,ji = . . .. As a consequence, matrix elements vanish identically if two or more of their hole indices are equal. It is important for the implementation to omit all permutationally related or zero terms in the DC-ADC matrix which can be achieved by following a restricted enumeration scheme for the multi-indices. For example, in the (rklm) multi-index the confinements k > l > m ∀ r ∈ {virt} ∧ k, l, m ∈ {occ} have to be obeyed.

One major factor essentially governing our realization is molecular symmetry which is fully exploited up to all Abelian subgroups of the D * 2h double group. Since the molecular spinors transform according to the molecular double group G * and not according to the single-valued counterpart G the use of double groups containing so-called extra or fermionic irreducible representations (irreps) γ i is mandatory for DC-ADC. The single-valued counterpart G of G * contains the bosonic irreps β i ∈ G. All irreps of Abelian subgroups to G * are one-dimensional and exhibit the following useful property:

γ * i ⊗ γ i = Γ 0 , γ i ∈ G * and β * i ⊗ β i = Γ 0 , β i ∈ G (10)
with Γ 0 being the totally symmetric irrep of G * (G). The complex conjugate of an irrep therefore is its own inverse allowing for a fast computation of target irreps being part of a multiple direct product. Furthermore, the direct product of two fermionic (bosonic) irreps yields a bosonic irrep, whereas the mixed product yields a fermionic irrep. The highest available point group symmetry in the employed DIRAC program package [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC08[END_REF] is D * 2h followed by D * 2 and C * 2v . The bosonic irreps of these groups are one-dimensional and the corresponding bosonic subgroups therefore are Abelian. However, for the corresponding fermionic representations nonseparable two-dimensional irreps are obtained not leading to an overall Abelian group structure anymore. This can be seen, for example, in the C * 2v group where E 1/2 ⊗ E 1/2 yields a four-dimensional reducible representation which can be decom-

posed into A 1 ⊕ A 2 ⊕ B 1 ⊕ B 2 .
Therefore a projection of the non-Abelian double groups to their Abelian subgroups is performed according to

D * 2h → C * 2h , D * 2 → C * 2 and C * 2v → C * 2 .
For these subgroups the number of fermionic and bosonic irreps is identical and their sum equals the number of classes in G * allowing for an efficient symmetry handling.

In our implementation the storage of integrals and the tensorial contractions are symmetry-driven. In order to demonstrate this for the V ij[kl] integrals as a representative example one always follows the all-encompassing requirement

γ * i ⊗ γ * j ⊗ γ k ⊗ γ l ! = Γ 0 , (11)
otherwise V ij[kl] vanishes. For defining a label-free storage sequence we simply follow the generic order (γ 1 , . . . , γ N ) starting with the rightmost irrep γ l ∈ G * . The subsequent irreps γ k , γ * j and γ * i are then looped over accordingly always maintaining relation (11) where the irreps γ i and γ j have to be taken as complex conjugates due to their relation to the bra space.

For a two-indexed array i|j we then obtain the following storage scheme with respect to the irrep sequence:

i|j : γ * 1 |γ 1 , γ * 2 |γ 2 , . . . , γ * N |γ N with γ * q ⊗ γ q = Γ 0 .
In the case of four-index quantities ij||kl the situation is a bit more complex since on the bra and on the ket side the two fermionic irreps first combine to bosonic irreps which in turn are again stored in generic order:

ij||kl : β (ij) * 1 |β (kl) 1 , β (ij) * 2 |β (kl) 2 , . . . β (ij) * N |β (kl) N with β (ij) * q ⊗ β (kl) q = Γ 0 .
Hereby the superscripts (ij) and (kl) indicate the fermionic parent irreps and the subscripts enumerate the bosonic irreps. For a given β q we recursively apply this logic again and obtain a stream sequence like

β (ij) q : γ i 1 |γ 1 , γ i 2 |γ 2 , . . . γ i N |γ N with γ i k ⊗ γ j ! = β (ij) q . ( 12
)
Both fermionic irreps producing β

(ij) q
either stem from the bra or from the ket space and relation ( 12) also applies to the bra space after complex conjugation. As soon as the smallest symmetry cell, γ i k |γ j , is established the individual spinors occurring therein are enumerated generically completing the index-free storage scheme. This label-free storage is also applied for the relativistic coupled cluster [START_REF] Visscher | [END_REF] and configuration interaction module in the DIRAC program suite [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC08[END_REF].

In the general two-electron integral V ab One very important property of the total ADC matrix which holds in the one-and fourcomponent case is its block structure with respect to different bosonic final state symmetries (see Fig. 1). Hereby no coupling between blocks of different symmetry occurs which can be seen as a consequence of relation (11). By this grouping into disjoint blocks the numerical effort for the matrix diagonalization is drastically reduced. It should be noted that the symbols I 2 , J 2 , . . . in Fig. 1 stand for the twofold and fourfold multi-indices of Eqs. ( 7)- (9) in exactly the same index ordering. The symmetry requirements for nonvanishing 2h/2h matrix elements C ij,kl of Eq. ( 7) now read as bra :

β (ij) * 1 = γ * i ⊗ γ * j , ket : β (kl) 2 = γ k ⊗ γ l , β (ij) * 1 ⊗ β (kl) 2 ! = Γ 0 . ( 13
)
After multiplying the last equation by β

(ij) 1
and considering (10) we obtain relations for the bra and ket irreps in order to obtain nonvanishing matrix elements as

β 1 := β (ij) 1 ! = β (kl) 2 =: β 2 .
This holds for all available bosonic irreps β q in G * , q = 1, . . . N leading to the block structure.

The present considerations can be applied in a similar way to the coupling and satellite blocks where care must be taken to perform integral contractions exclusively over corresponding bra/ket pairs. The transition from Abrikosov diagrams to their algebraic counterpart exactly determines the position of a specific spinor index in the bra or ket space. The left-hand-side occurrence of the (rklm) multi-index in C rklm,ij (Eq. 8), for example, does not automatically imply a bra character for them. Instead, from Eq. ( 8) the correct bra/ket sequence for this multi-index can be derived according to

C rklm,ij : (rklm, ij) ↔ | | | |, | | (14)
and the symmetry conditions are established as

γ r ⊗ γ * k ⊗ γ * l ⊗ γ * m ⊗ γ i ⊗ γ j ! = Γ 0 . (15)
The pair γ i , γ j hereby defines the final state symmetry

β f := β (ij) q
= γ i ⊗ γ j of the DC-ADC matrix block. In order not to omit contributions we loop through all possible (γ k , γ l , γ m ) triples and fix γ r after complex conjugation of (15) via

γ * r ⊗ γ klm = β f with γ klm = γ k ⊗ γ l ⊗ γ m =⇒ γ r = β * f ⊗ γ klm .
The numerous delta-conditions in the satellite block (9) test for like spinor labels which are a priori not available in a purely symmetry-driven storage. Since the DC-ADC matrix is to be constructed for each Lanczos iteration and to be stored in a packed row/column form on disk it is more efficient to utilize a row/column-based addressing in the satellite block and recover the spinor index quadruple via a lookup table whenever needed.

III. COMPUTATIONAL DETAILS

For all calculations the Dirac-Coulomb Hamiltonian with inclusion of the (SS|SS) integrals was employed and the relevance of spin-orbit effects was demonstrated by additionally performing ADC2 calculations with an underlying spin-free Hamiltonian by Dyall [34,35]. In the latter case only scalar relativistic effects are taken into account and the symmetry is governed by the normal group instead of the double group. For xenon the dual-type (26s21p16d4f 4g2h) primitive exponent set from [36] was used and is shown in table I. The duality condition hereby restricts the exponent sets according to {d} ⊂ {s}, {f } ⊂ {p}, {g} ⊂ {s} and {h} ⊂ {p} by which the number of (SS|SS) integrals to be determined in a four-component calculation can be considerably reduced. Spinors in the energy range from -10.0 a.u. up to +100.0 a.u. were taken into account for the electron correlation treatment which corresponds to an inclusion of all occupied Xe orbitals with n ≥ 4.

For the bromine molecule the Dyall (23s16p10d1f ) triple zeta basis set [37] was taken in its uncontracted form and all spinors in the energy range of -4.0 a.u. to +10.0 a.u.

were included in the active correlation space. It should be kept in mind that a molecular calculation is much more elaborate than an atomic one and the range of active spinors therefore had to be restricted accordingly. The experimental distance of R = 2.28105 Å [START_REF] Huber | Molecular Spectra and Molecular Structure Constants of Diatomic Molecules[END_REF] was taken for the Br 2 calculations.

All electronic structure calculations have been carried out with the program system DIRAC [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC08[END_REF] which is a general-purpose four-component package. The relativistic implementation of the two-particle propagator DC-ADC2 is available in the developers' version at the moment and will soon be part of the official release. In addition to spatial symmetry time-reversal symmetry is exploited throughout the DC-ADC2 code. In the current implementation DIRAC also supports linear symmetry which allows for a Ω-coupled final state classification of the bromine molecule and also for a M J classification in the xenon atom.

IV. RESULTS AND DISCUSSION

A. Xenon

As a first application we provide theoretical DC-ADC2 results for the Xe 2+ final states and corresponding line spacings together with a simulated double ionization spectrum of Xe in the low energy regime comprising the 5p 4 configuration. Experimental levels and transitions were obtained by Persson et al. [START_REF] Persson | [END_REF] who also performed theoretical calculations based on Hartree-Fock theory with inclusion of some configuration interaction effects. These energy levels also served for the determination of a xenon N 4,5 OO Auger spectrum [40]. In

) 2 (5p 3
2

) 2 , (5p 1 2 ) 1 (5p 3
2

) 3 and (5p 1 2 ) 0 (5p 3
2 ) 4 occupations in the jj-coupled case.

J level Exp. [START_REF] Persson | [END_REF] DC-SCF DC-ADC2 ∆(DC-SCF) ∆(DC-ADC2)

2 ) 2 (5p 3 2 ) 2 3 P 2 reference state to the various excited states with J = 0, 1, 2, 0. Uncorrelated (DC-SCF) and correlated (DC-ADC2) results are included. In the last two columns the difference to the experimental value is shown stressing the relevance of electron correlation contributions.

The inadequacy of a pure LS coupling scheme can be seen in table II where the individual J terms belonging to a specific 2S+1 L manifold are significantly separated (the splitting of the 5p 1/2 and 5p 3/2 spinors in neutral xenon amounts to 1.44 eV at the DC-SCF level). We note that the obtained energetic order in a relativistic description is identical to the experimental one (see table II) and the corresponding level spacings at the correlated DC-ADC2 level are in much better agreement to experiment than they are for an uncorrelated description (DC-SCF).

J level Energy cm

(J = 0), Φ 2 = p 2 1 2 p 2 3 2 (J = 2), Φ 3 = p 1 1 2 p 3 3 2 (J = 2), Φ 4 = p 1 1 2 p 3 3 2 (J = 1), Φ 5 = p 0 1 2 p 4 3 2 (J = 0)
Adapting the jj coupling scheme we expect configurations belonging to the same J value to mix as observed at the Dirac-Hartree-Fock level already. In a relativistic average open shell calculation of the xenon [core](5p) 4 electron configuration (e.g. with GRASP [41,42])

the CSFs for individual J levels are mixtures of different two-particle jj-coupled Slater determinants (see table III for the atomic results). This multi-configurational character will also transfer to the propagator results.

In table IV a relativistic final state analysis for the 5p 4 configuration is presented where the spinors carrying the holes are listed. The calculated final states all possess main state character and genuine 2h configurations are strongly prevailing over 3h1p admixtures. Since the calculations were performed applying linear symmetry we have immediate access to the projection of the total angular momentum J onto the z axis. These different M J values occur in separate irreps allowing for a detailed classification of the final states. For example, a DC-SCF calculation of the J = 2 ground state yields a 0.95p 2 1/2 p 2 3/2 + 0.31p 1 1/2 p 3 3/2 state composition. According to the notation in table IV this corresponds to a p(m j 1 )p(m j 2 ) and a p(m j 1 )p(m j 2 ) 2h state. As in the DC-SCF case the p(m j 1 )p(m j 2 ) configuration remains the dominant part of the Xe 2+ ground state for all M J projections (in the case of M J = 0 there are two realizations with respect to the individual m j projections but the sum of the squares is identical to the square of the M J = 0 contributions) and therefore shows resemblance to the DC-SCF ground state with respect to the configuration mixing. However, the decrease of the pole strength from unity to 0.8818 clearly indicates contributions from other 2h and 3h1p states accounting for electron correlation. Due to their large number and smallness they are not listed in table IV Additionally, the proper inclusion of spin-orbit coupling becomes evident when a scalar relativistic propagator calculation is performed setting out from the spin-free Hamiltonian (table V). Now the double ionization potentials (DIPs) exhibit a ninefold, fivefold and onefold degeneracy corresponding to the 3 P , 1 D and 1 S LS-coupled manifold of a 5p 4 electron configuration and only two transitions at 9630.3 cm -1 and 27087.3 cm -1 can be predicted where just the first one is in reasonable agreement with experiment. Obviously, this spin-free description is inadequate even for a qualitatively correct interpretation of the spectrum. suitable linear combinations within the degenerate manifolds are to be formed. Similar to the four-component calculation the spin-free states resemble a Hartree-Fock state but with additional contributions accounting for electron correlation.

As it was mentioned in section I the propagator method gives immediate access to the complete final state manifold where the pole strengths relate to the peak heights in the spectrum. The coefficients of the participating 2h configurations constituting a specific final state can hereby be obtained from the eigenvector. Employing this information, a theoretical spectrum can be immediately generated. In order to compare the DC-ADC2 stick spectrum to the experimental one possessing a natural linewidth we convolve the theoretical spectrum with Lorentzian envelopes. Hereby the reported spectroscopic resolution serves as the full width at half maximum (FWHM) value which was taken as Γ = 50 meV for the Xe and Br 2 spectrum. The simulated peaks which appear at the energies E i then attain the following functional form:

L i (E) = 1 π Γ/2 (E -E i ) 2 + (Γ/2) 2 . ( 16
)
The total spectrum I tot (E) is obtained by summing over all contributions L i (E) multiplied by the pole strength X i obtained from the DC-ADC2 calculation as

I tot (E) = i L i (E)X i . (17)
The simulated experimental Xe 2+ spectrum is shown in Fig. 2 and the lowest DIP was shifted to the experimental onset. The areas under the peaks exhibit conformity with the statistical weights of the populated levels and the theoretical spectrum well resembles the experimental one in [43].

B. Molecular bromine

Recently, Fleig et al. calculated potential energy curves for the lowest g and u states of Br 2+ 2 and proposed a new double ionization mechanism [21]. For the determination of the dicationic final states relativistic CI [44,45] and MRCC approaches [46,47] together with a large ANO-RCC basis were employed. It is therefore of prominent interest to compare our propagator approach to these wave function-based results. It was not our intention to reproduce the Br 2+ 2 potential energy curves for a characterization of special ionization processes proposed in [21] but to demonstrate the usefulness of the four-component propagator approach with respect to final state and spectrum analysis for this molecule. The orbital energy diagram at the self consistent field level for the outer valence p orbitals is shown in Fig. 3 obtained with the Lévy-Leblond (nonrelativistic), the spin-free (purely scalar relativistic) and Dirac-Coulomb Hamiltonian (scalar relativistic and spin-orbit). One observes an average spin-orbit splitting of 6.9 mH for the π g , π u outer valence orbitals that prevails over scalar relativistic energy changes of 1.0 -2.4 mH. Due to open π g,u shells in Br 2+ 2 spin-orbit coupling will play a role in the corresponding spectra. The classification of the molecular orbitals follows the ΛΣ and Ω coupling schemes, respectively. For the onset of the double ionization Fleig et al. obtained a value of 26.94 eV compared to the experimental onset of 28.55 eV (vertical) and 28.25 eV (adiabatic) [START_REF] Yencha | Proceedings of the XXIV ICPEAC Conference[END_REF] where in the ADC2 approach 27.39 eV (DC-ADC2) and 27.54 eV (SF-ADC2) were obtained. The interplay of electron correlation and relaxation effects together with spin-orbit coupling is rather involved and puts high demands on the theoretical methodology. The deviations of the DC-ADC2 DIPs from experiment mainly result from an insufficient inclusion of electron relaxation effects which are larger for a double ionized final state than for excited or singly ionized final states. This deficiency could be cured by increasing the order of the perturbation treatment but will not distort the physical picture obtained from the current calculations. It should be mentioned that a strict second order DC-ADC where the 3h1p/3h1p particle (satellite) block is included only in zeroth order would be too inaccurate and the additional first order contributions (the various V ab,[cd] terms in eq. ( 9)) considerably improve on the final state energies. It is common practice to use this so-called extended ADC2 routinely as can be seen from eq. ( 9) where these first order contributions are included. The satellite block will not change any further when it comes to an extension to third-order. The energy differences of the outer valence final state energies obtained at the DC-ADC2 and wave function-based level of Fleig et al. exhibit an average deviation of 0.084 eV (see table VI). With the exception of the DIP onset a good agreement between both approaches is therefore observed.

LL DC SF

The second aspect to be addressed is the comparison of the individual final state compositions. Table VII lists the DC-ADC2 DIPs and pole strengths together with the detailed configuration information for the 10 lowest final states in Br 2+ 2 . Hereby the explicit Kramers partners occurring in the degenerate ±Ω pairs are reported separately. All outer valence states can be identified as main states which are states of high pole strength and no significant admixture of 3h1p configurations. This is beneficial for a future two-hole population analysis because for main states it suffices to include the 2h space in the population analysis without loss of accuracy [48]. States with Ω = 0 are doubly degenerate. a) Our work, b) Ref. [21].

Ω ∆ a) ∆ b) Ω ∆ a) ∆ b) 0 g 0.
For a direct comparison of a specific final state configuration one has to account for the different notations employed in [21] and in our work (table VII). We report the spinors where the hole was created (indicated by the -1 superscript) whereas Fleig et al. list the remaining occupied spinors after the double ionization. In this respect we find a dominant π-1

3/2g π -1 3/2g , a less pronounced π-1

1/2g π -1 1/2g and a weak π-1 3/2u π -1 3/2u contribution to the Br 2+ 2 ground state which also seems natural from an energetic point of view (see Fig. 3). Fleig [21]. The remaining states in this energy range exhibit a good resemblance to the wave function-based results. Already for the ground state of the doubly ionized system considerable configuration mixing is observed which can not be reproduced by a single reference self consistent field calculation. This strong configuration mixing is quite characteristic for doubly ionized states and is already observed in light systems such as BF 3 [48]. This is in contrast to outer valence final states of singly ionized systems which in most cases possess only one major contribution in the intermediate state representation [29]. Due to the relevance of spin-orbit coupling for a correct interpretation of the spectral structure spin-free results are not discussed here.

In an analogous way we convoluted the Br 2+ 2 stick spectrum with Lorentzian curves in the energy range of 28 and 32 eV (Fig. 4) in order to simulate the experimental spectrum. The experimental peak positions in table VIII together with their nonrelativistic designations were obtained by Yencha et al. [START_REF] Yencha | Proceedings of the XXIV ICPEAC Conference[END_REF]. Hereby the a 1 ∆ g and the 3 ∆ u states will undergo spin-orbit splitting.

Comparing the theoretical spectrum with the experimental threshold photoelectron coincidence spectrum one finds a substantial discrepancy of the intensities for the lowest u states (indicated by the dashed box in Fig. 4). The intensity of a transition is controlled by the magnitude of the transition matrix element including the initial and final state wave function. In the absence of strong vibronic coupling effects one can in a first approximation write the full wave function as a product of an electronic and nuclear part yielding additional overlap matrix elements over nuclear wave functions, the well-known Franck-Condon States with Ω = 0 are doubly degenerate and the configurations are totally antisymmetric with respect to particle interchange. A Kramers partner is indicated by a barred symbol and π denotes the fourfold degenerate parent orbital before inclusion of spin-orbit coupling.

the lowest excited u states of Br 2+ 2 in contrast to the bound 0 g (1), 1 g (1), 2 g (1) and 0 g (2)

states supporting vibrational levels. It has been shown that in case of an excitation to a decaying state complicated nuclear dynamics lead to a substantially altered spectral structure [49]. This obviously happens in the considered spectral range for the u states of Br 2+ 2 . Additionally, energies and band shapes are strongly influenced by vibrational motion and Franck-Condon factors therefore should be taken into account if one seeks a reliable prediction of these spectra [50,51]. The purely electronic DC-ADC2 results do not incorporate ted as can be seen from the lower part of Fig. 4. The main changes hereby comprise a reduced structure, partial peak shifts and alterated intensities especially for the u states.

The DIP onset for the spin-free spectrum is by 0.2 eV higher than the DC-ADC2 onset and the theoretical spectrum was also shifted to the experimental DIP onset. Obviously the spectral features even for the moderately relativistic bromine molecule are not reproducible by ignoring spin-orbit coupling in the propagator calculations, a situation which will become more pronounced in heavier systems.

V. SUMMARY

In this work we have shown the methodology and capability of the four-component twoparticle propagator for spectrum analysis and final state classification. The efficient exploitation of symmetry in the integral handling and ADC matrix construction was demonstrated and the performance tested on two systems, the xenon atom and the bromine molecule.

These calculations showed that without a proper account for relativistic effects neither of the two spectra can be adequately described. Compared to the high-level four-component multi-reference coupled cluster treatments our results are in excellent agreement with respect to line spacings and state analysis, already at the extended second-order stage of the propagator implementation. The deviations from experiment especially in the intensities for the u states can be understood by not accounting for nuclear dynamical effects occurring after excitation into decaying states. Subsequent analysis will be extended by inclusion of vibrational effects and a two-hole population analysis allowing for a detailed insight into electronic decay processes of clusters containing heavy elements. Already in the case of the moderately relativistic bromine molecule inclusion of spin-orbit coupling is mandatory for a reasonable resemblance of the experimental structure. Due to its second-order character the DC-ADC matrix construction scales similarly to the second-order Møller-Plesset (MP2) method and at the moment the AO→MO integral transformation and the band Lanczos diagonalization represent the computational bottlenecks in large DC-ADC2 calculations.

VI. ACKNOWLEDGMENTS

  [cd] the indices a, b, c, d stand for occupied(O) or virtual (V ) orbitals. It is evident from Eqs. (7)-(9) that only four distinct integral classes, namely the OO||OO , V O||OO , V V ||OO and V O||V O types are needed for the construction of the matrix elements. Hereby the V V ||OO and V O||V O classes are essentially distinct due to the different bra and ket origin of the O and V spinors. Every other order of indices such as e.g. V is[rm] with i, m ∈ {occ}, r, s ∈ {virt} can be brought into one of the four generic forms utilizing the equivalence relations(6). The matrix elements of Eqs. (7)-(9) are therefore transformed to a form suitable for direct use of the four generic integral classes.

FIG. 1 :

 1 FIG. 1:Illustration of the symmetry block formation of the total ADC matrix with respect to the final state irreps β i . The symbols I 2 , J 2 (I 4 , L 4 ) denote the twofold (fourfold) multi-indices from Eqs. (7)-(9).

  ) and pole strengths (PS) for the final states of gerade symmetry derived from the 5p 4 configuration of xenon. The most prominent change compared to the DC-ADC2 is the high degree of degeneracy and the concomitant loss of spectral structure. The configurations are to be understood as being totally antisymmetric with respect to particle interchange.

FIG. 2 :

 2 FIG. 2: Simulated experimental Xe 2+ spectrum obtained from the stick spectrum by convolution with a Lorentzian curve and shifted to the experimental onset.

FIG. 3 :

 3 FIG.3: Orbital energy diagram (in a.u.) for the outer valence p space of neutral Br 2 obtained with the Lévy-Leblond (LL), the spin-free (SF) and Dirac-Coulomb (DC) Hamiltonian.

TABLE I :

 I Dual

primitive exponent set for xenon used in the DC-ADC2 calculations. The notation {s10 -s25} hereby indicates the s exponent range used for the d functions asf.

  table II we collect the experimental data by Persson et al. for the gerade final states and our theoretical results obtained by the DC-ADC2 procedure. Uncorrelated self-consistent field results (DC-SCF) are also given and stress the relevance of electron correlation contributions.

	All states listed emerge from a valence doubly ionized 5s 2 5p 4 configuration (nonrelativistic
	notation) leading to possible (5p 1 2

TABLE II :

 II Experimental

and theoretical transition energies (in cm -1 ) from the lowest (5p

TABLE III :

 III Composition of the J-adapted CSFs at the DC-SCF level. Φ 1 = p 2

	1	p 2 3
	2	2

TABLE IV :

 IV . Double ionization potentials (DIPs, in eV) and pole strengths (PS) for the 2h final states of gerade symmetry in jj coupling. The p(p) symbols denote the p 3 2 (p 1 2 ) spinors and the corresponding individual m j projections are given in parentheses. M J = m j1 + m j2 .

	DIP	PS	Configurations	M J
	32.7964 0.8818 -0.1786 p(-1 2 )p( 1 2 ) -0.1786 p(-1 2 )p( 1 2 ) -0.6395 p(-1 2 )p( 1 2 ) -0.6395 p(-3 2 )p( 3 2 ) 0
	32.7964 0.8818	0.2187 p( 1 2 )p( 1 2 ) +0.1263 p( 3 2 )p(-1 2 ) +0.9044 p( 3 2 )p(-1 2 )	+1
	32.7964 0.8818	-0.1263 p(-3 2 )p( 1 2 ) -0.9044 p(-3 2 )p( 1 2 ) -0.2187 p(-1 2 )p(-1 2 )	-1
	32.7964 0.8818	0.2526 p( 3 2 )p( 1 2 ) -0.9044 p( 3 2 )p( 1 2 )	+2
	32.7964 0.8818	0.2526 p(-3 2 )p(-1 2 ) -0.9044 p(-3 2 )p(-1 2 )	-2
	33.8494 0.8840	-0.4676 p(-1 2 )p( 1 2 ) +0.5759 p(-1 2 )p( 1 2 ) -0.5759 p(-3 2 )p( 3 2 )	0
	33.9900 0.8749	0.4677 p( 1 2 )p( 1 2 ) -0.8100 p( 3 2 )p(-1 2 )	+1
	33.9900 0.8749	0.8100 p(-3 2 )p( 1 2 ) -0.4677 p(-1 2 )p(-1 2 )	-1
	33.9900 0.8749	-0.6614 p(-1 2 )p( 1 2 ) +0.6614 p(-1 2 )p( 1 2 )	0
	34.9468 0.8744	-0.4501 p(-3 2 )p( 1 2 ) +0.2532 p(-3 2 )p( 1 2 ) -0.7796 p(-1 2 )p(-1 2 )	-1
	34.9468 0.8744	0.7796 p( 1 2 )p( 1 2 ) +0.4501 p( 3 2 )p(-1 2 ) -0.2532 p( 3 2 )p(-1 2 )	+1
	34.9468 0.8744	-0.9001 p( 3 2 )p( 1 2 ) -0.2532 p( 3 2 )p( 1 2 )	+2
	34.9468 0.8744	-0.9001 p(-3 2 )p(-1 2 ) -0.2532 p(-3 2 )p(-1 2 )	-2
	34.9468 0.8744 0.6365 p(-1 2 )p( 1 2 ) +0.6365 p(-1 2 )p( 1 2 ) -0.1790 p(-1 2 )p( 1 2 ) -0.1790 p(-3 2 )p( 3 2 ) 0
	37.4880 0.8945	-0.8070 p(-1 2 )p( 1 2 ) -0.3383 p(-1 2 )p( 1 2 ) +0.3383 p(-3 2 )p( 3 2 )	0

TABLE V :

 V The spin-free final state classifications in table V occur with respect to Cartesian p x , p y and p z orbitals obtained in D 2h symmetry. The spin projections (M S values) can be read off directly from the final state configurations in table V and for a proper M L assignment .19955 p xβ p xα -0.73886 p yβ p yα +0.53930 p zβ p zα 0 34.6633 0.8767 0.66203 p zβ p yα +0.66203 p yβ p zα 0 34.6633 0.8767 0.66203 p zβ p xα +0.66203 p xβ p zα 0 34.6633 0.8767 -0.66203 p yβ p xα -0.66203 p xβ p yα 0 34.6633 0.8767 0.73795 p xβ p xα -0.19616 p yβ p yα -0.54179 p zβ p zα 0 36.8277 0.9066 -0.54485 p xβ p xα -0.54485 p yβ p yα -0.54485 p zβ p zα 0 Spin-free (LS coupled) DIPs (in eV

	DIP	PS	Configurations	M S
	33.4693 0.8778	0.66249 p zβ p yα -0.66249 p yβ p zα	0
	33.4693 0.8778	0.66249 p zβ p xα -0.66249 p xβ p zα	0
	33.4693 0.8778	-0.66249 p yβ p xα +0.66249 p xβ p yα	0
	33.4693 0.8778	-0.93690 p zα p yα	1
	33.4693 0.8778	0.93690 p zα p xα	1
	33.4693 0.8778	0.93690 p zβ p yβ	-1
	33.4693 0.8778	0.93690 p zβ p xβ	-1
	33.4693 0.8778	-0.93690 p yα p xα	1
	33.4693 0.8778	0.93690 p yβ p xβ	-1
	34.6633 0.8767 0

TABLE VI :

 VI Energy differences ∆ (in eV) of the individual Ω-coupled states to the lowest state of Br 2+ 2 .

	000 0.000 3 u 1.538 1.418
	1 g 0.126 0.149 2 u 1.882 1.774
	2 g 0.639 0.622 1 u 1.949 1.854
	0 g 1.134 1.059 0 u 2.207 2.112
	0 u 1.491 1.365 1 u 2.229 2.133

(

  FC) factors. By extensive calculations Fleig et al. found a purely repulsive character of

	DIP (eV) PS Final state configuration		Ω
	27.393 0.8574 0.40229π -1 1/2g π -1 1/2g	+0.81293π -1 3/2g π -1 3/2g	0g
	-0.18101π -1 3/2u π -1 3/2u
	27.519 0.8539 -0.90975π -1 3/2g	π-1 1/2g -0.15729π -1 3/2u	π-1 1/2u	+1g
	27.519 0.8539 0.90975π -1 3/2g π -1 1/2g +0.15729π -1 3/2u π -1 1/2u	-1g
	28.032 0.8668 -0.88962π -1 3/2g π -1 1/2g +0.27335π -1 3/2u π -1 1/2u	+2g
	28.032 0.8668 -0.88962π -1 3/2g	π-1 1/2g +0.27335π -1 3/2u	π-1 1/2u	-2g
	28.527 0.8712 -0.78746π -1 1/2g π -1 1/2g	+0.36346π -1 3/2g π -1 3/2g	0g
	+0.26375π -1 1/2u π -1 1/2u -0.21103π -1 3/2u π -1 3/2u
	28.884 0.8777 0.65611π -1 3/2u π -1 3/2g -0.65611π -1 3/2u	π-1 3/2g	0u
	28.931 0.8806 0.93841π -1 3/2u π -1 3/2g			+3u
	28.931 0.8806 -0.93841π -1 3/2u	π-1 3/2g			-3u
	29.275 0.8782 -0.65854π -1 3/2u π -1 1/2g +0.66591π -1 1/2u π -1 3/2g	+2u
	29.275 0.8782 -0.65854π -1 3/2u	π-1 1/2g +0.66591π -1 1/2u	π-1 3/2g	-2u
	29.342 0.8783 -0.65105π -1 3/2u	π-1 1/2g -0.67280π -1 1/2u π -1 3/2g	+1u
	29.342 0.8783 0.65105π -1 3/2u π -1 1/2g +0.67280π -1 1/2u	π-1 3/2g	-1u
	29.600 0.8742 0.65371π -1 1/2u π -1 1/2g -0.65371π -1 1/2u	π-1 1/2g	0u
	29.622 0.8757 0.93404π -1 1/2u π -1 1/2g			+1u
	29.622 0.8757 -0.93404π -1 1/2u	π-1 1/2g			-1u

TABLE VII :

 VII DC-ADC2 DIPs (in eV) and pole strengths (PS) for the 10 lowest 2h final states of Br 2+ 2 .

The author gratefully acknowledges financial support by the Deutsche Forschungsgemeinschaft.

these nuclear effects and lead to deviations for the decaying u states. For the bound states of Br 2+ 2 , however, a good agreement with experiment is observed already without the inclusion of nuclear effects. A substantial alteration of the spectrum takes place if spin-orbit contributions are omit-


F Anzà

B Militello

A Messina

Tripartite thermal correlations in an inhomogeneous spin-star system

Keywords: numbers: 03.65.Ud, 03.67.Mn, 75.10.Jm

Tripartite thermal correlations in an inhomogeneous spinstar system

F Anzà, B Militello, A Messina

Introduction

The notion of thermal entanglement relies on the amount of entanglement possessed by a physical system when it had undergone a thermalization process that has brought it in a thermal state [1]. In fact, in the presence of interaction between different parts of a compound system, even after a complete thermalization, the system can exhibit appreciable quantum correlations since the Hamiltonian eigenstates are, in general, entangled states. The establishment of a relation between temperature and entanglement has quickly brought to the idea of using the entanglement as an order parameter in quantum phase transitions [2,3]. Moreover, it has given a stronger stimulus of searching quantum correlations in the macroscopic world [4,5,6].

Thermal entanglement has been studied in connection with many possible applications and in different physical systems: detailed analysis in spin chains described by Heisenberg models has been given [7] as well as in atom-cavity systems [8] and in simple molecular models [9]. The quantumness of correlations has been singled out in thermalized spin chains [10] also in connection with non local effects [11], and application of thermal entanglement in optimal quantum teleportation protocols has been proposed [12].

Spin systems have been studied in depth not only in the neighbor-interaction configuration, leading to the Heisenberg model, but also in star networks. This kind of systems can be realized in many physical contexts, from Josephson Junctions [13] to trapped ions [14] to solid state physics [15]. In 2004, Hutton and Bose have analyzed a physical system made of a central spin interacting with a set of N outer spins at zero temperature, bringing to light interesting properties which are immediately traceable back to the features of the ground state of this spin-star network [16]. They have shown that the evenness or oddness of the number of outer spins determines the law the entanglement amount scales with. The same spin configuration, in a simplified version involving three outer spins only, has been recently considered by Wan-Li et al [17] in the special case where all the interactions between the central spin and the outer ones are identical. Exploiting the high degree of symmetry of the system, they concentrate on appropriate concurrences [18] to extract information about the presence of quantum correlations. However, the disclosure of quantum correlations in tripartite systems could require tools and concepts more adequate than the simple concurrences.

Sabin and Garcia-Alcaine have contributed to solve this still very open problem [START_REF] Amico | For a review on entanglement in multipartite systems[END_REF] introducing the notion of tripartite negativity [20], which is an effective tool to reveal the existence of quantum correlations traceable back to the impossibility of separating any of the three subsystems from the other two.

In this paper we consider a spin-star system where three outer spins are coupled to the central one with different strengths, and analyze thermal entanglement in different situations. In the next section we present the model under scrutiny, which is characterized by anisotropic spin-spin interactions due to the absence of longitudinal z-z couplings. In the third section we consider the system in the homogeneous case, showing the presence of sharp changes in the amount of entanglement at zero temperature due to abrupt variations of the ground state. In the fourth section we analyze the inhomogeneous model, focusing on two types of inhomogeneity. Finally, in the last section some comments and conclusive remarks are given.

Physical scenario

System and Hamiltonian -In this section we present the spin star-system we have analyzed, which is pictured in figure 1 we indicate the coupling constants of the central spin with the outer spin 1, 2 and 3, respectively. Moreover, the whole system is immersed in a constant uniform magnetic field of modulus B 0 and directed along the z-axis. The Hamiltonian model is then given by (with = 1):

H = H 1 + H 2 , (1)
H 1 = - γ s 2 B 0 i σ i z = ω 0 2 i σ i z , (2)
H 2 = i=1,2,3 c i σ C + σ i -+ σ C -σ i + , (3)
where

σ i z = |+ i +| -|-i -|, σ i + = |+ i -|, σ i -= |-i +|, ω 0 is the unperturbed
Bohr frequency, and c i 's are coupling constants.

The term H 1 , describes the interaction of the system with the magnetic field, while the second one, H 2 , arises from the dipole-like interactions between the central spin and the outer ones. The absence of z-z interaction reflects a certain degree of anisotropy.

Assuming that the system is in a thermal state, its density operator takes the form:

ρ = e -H k B T Tr e -H k B T = e -H k B T Z , (4)
so that ρ and H share the same eigenstates. The Hamiltonian model given in Eq. (1) takes into account possible inhomogeneities, nevertheless it is surely of interest studying it in the homogeneous case (c 1 = c 2 = c 3 = c) which turns out to be simpler from a mathematical point of view, paving the way to the study of the more complicated inhomogeneous model. Wan-Li et al. [17] have analyzed the homogeneous model taking advantage of the concurrence to quantify entanglement. In addition, they have assumed a more general anisotropic interaction. Since in real situations the construction of a system with perfectly homogeneous interactions could be quite difficult, in this paper we investigate the features of thermal quantum correlations in the presence of inhomogeneity.

Tripartite negativity -Remarking that the analysis is performed on the outer subsystem, which is made of three parts, we have the need to study quantum correlations in tripartite systems. There are various proposals of tripartite entanglement quantifiers [21,22,23,24] and witnesses [26,27,28], the latter ones based on the key assertions in [25]. Nevertheless, up to now and to the best of our knowledge, there is not a definitive answer to the request of quantifying tripartite entanglement [START_REF] Amico | For a review on entanglement in multipartite systems[END_REF]. For instance, it has been shown recently that three-tangle [29] could be an improper tool for such a purpose [30]. The situation is much more complicated when we have to consider mixed states.

Indeed, all the recipes mentioned in [21,22,23,24,29] refer to pure states and fail when applied to non-pure states. Since we are studying thermal correlations, therefore dealing with highly mixed states, we need a tool to investigate quantum correlations in mixed states of tripartite systems. Again, to the best of our knowledge, the tripartite negativity introduced by Sabin and Garcia-Alcaine [20] is a good quantity that allows to study quantum correlations in tripartite systems, even in the case of mixed states.

Other quantities are very much related to the specific structure of the analyzed system, since observables that assume values in a certain range when the state is entangled are considered. Sabin and Garcia-Alcaine apply bipartite negativity to all the three possible bipartitions that can be singled out isolating one subsystem and considering the other two as a whole; then they consider the geometric mean of these three quantities.

According to the definition of negativity [31,32], the partial negativity related to the bipartition I -(JK) (in which the two subsystems J and K are considered as a whole) is given by

N I-JK = i |σ i (τ T I )| -1 , (5)
where σ i (τ T I ) is the i-th eigenvalue of τ T I , which is the partial transpose related to subsystem I, of the total (IJK) density matrix. The parameter used to study correlations in tripartite systems is then

N 123 = 3 N 1-23 N 2-13 N 3-12 , (6)
which is the above mentioned tripartite negativity. In [20] the following properties have been proven: i) τ separable or simply bi-separable ⇒ N 123 = 0;

ii) Invariance of N 123 under LU operators;

iii) Monotonicity of N 123 under LOCC operators;

Though N 123 is not a sufficient condition to single out tripartite entanglement, it is an effective tool to study quantum correlations in tripartite systems. In fact,

finding N 123 = 0 implies that none of the three subsystems is separable, hence revealing correlations involving all the three subsystems. It is also important to note that negativity N I-JK does involve all the three spins, and establishes the existence of a correlation between I and the entire couple made of J and K, which is a very different operation from tracing over one spin, say K, and evaluating the degree of correlation between the other two, say I and J.

Homogeneous model

Let us consider our model in the special case c 1 = c 2 = c 3 = c that we have already addressed as the homogeneous case. The complete diagonalization of this Hamiltonian model is given in Appendix A. In order to compute tripartite negativity we need first to find explicit form of the outer-spin density operator ρ 123 , tracing over the degrees of freedom of the central spin.

It is quite clear that because of homogeneity, N 1-23 = N 2-13 = N 3-12 , so the geometric mean reduces to one of these three quantities. In figure 2 we plot N 123 as a function of temperature (T ) and coupling constant (c). One can see some interesting features of N 123 , for example its diminishing as T increases, and the presence of abrupt transitions at low temperature due to fast variations of the ground state. In figure 3,

where N 123 (k B T, c = 6ω 0 ) is plotted, it is better shown this behavior which is due to the fact that the more temperature increases the more ρ C123 approaches a maximally mixed states ( 1 16 I), making also ρ 123 maximally mixed. A more interesting trend of N 123 is observed at low temperature where abrupt variations are present, as shown in figure 4. At k B T = 0.01ω 0 fast entanglement variations are well visible at c ≈ 0.6ω 0 and c ≈ 3.7ω 0 . The reason for this occurrences is that in each of these points the system undergoes a level crossing involving its two lowest levels (ground and first excited states).

Let us call A 1 , A 2 , A 3 the three plateaux of figure 4, corresponding to (c < 0.6ω 0 )(0.6ω 0 < c < 3.72ω 0 ) and (c > 3.72ω 0 ), respectively. The entanglement in these areas is the same of the entanglement possessed by the ground state. In A 1 entanglement is evidently zero and, in fact, the ground state of the whole system is g (A 1 ) = |0000 , corresponding to ρ (A 1 ) 123 = |000 000| (with the notation |σ C σ 1 σ 2 σ 3 ). In regions A 2 and A 3 ground states of four spin system are g

(A 2 ) = |0100 + |0010 + |0001 - √ 3 |1000 / √ 6 and g (A 3 ) = ((|0011 + |0101 + |0110 ) -(|1100 + |1010 + |1001 )) / √ 6,

Model with inhomogeneity of type A

In spite of the lack of some symmetry, the inhomogeneous model of type A is still easily solvable and its diagonalization is given in Appendix A. The trend with respect to temperature is decreasing. The dependence of N 123 on the inhomogeneity parameter, at low temperature, exhibits different trends: initially, smooth changes, then, fast falling around x ≈ 5.5 and zero value from there on. is the same as for the homogeneous case, i.e. the high degree of mixedness. At low temperature, in the 0 < x < 1 area, the smaller is x, the closer to zero is the negativity.

Since one can see that in the Area

In fact, when x is quite smaller than 1, c 2 = c x is much smaller than c 1 = c 3 = c, and hence, roughly speaking, spin 2 can be considered almost decoupled, and then separable. Moreover, at x 0 ≈ 0.43 this negligibility of the coupling between spin C and spin 2 is concomitant to a level crossing between the two lowest energy levels,

E - 2 (c = 6ω 0 , x) = -3ω 0 x + √ 8 + x 2 and E - 4 (c = 6ω 0 , x) = -(6 √ 2 + x 2 + 1) ω 0 , that causes a sharp change at x = x 0 , where E - 2 = E - 4 .
Still, at low temperature, for x quite greater than 1, N 123 decreases as x becomes increasingly grater than 1, because the coupling with spin 2 becomes stronger than the other two, making spin 1 and 3 almost separable from the spin 2. Around x ≈ 5.5 a level crossing occurs. Because of the diminishing of N 123 for x 1 and x 1, there must be a maximum in the intermediate region. Unexpectfully, this maximum is reached at

x = x M ≈ 2.46, instead of x = 1. We think this is an important result because, even though there is no big difference in the values of negativity in a wide contour of the maximum (let's say from x = 0.5 to x = 5.5). In fact, it is conceptually important that the maximum of correlations is not reached in the homogeneous case, where the high symmetry of the system could lead to the idea of stronger correlations between its parts.

Evaluating the eigenvalues of Hamiltonian one can find that the ground state in the nearby of this maximum is Ψ - 2 , which becomes Ψ

ρ 2-(x) = 1 2 |w 1 (x) w 1 (x)| + 1 2 |w 2 (x) w 2 (x)| , (7)
with

|w 1 (x) = ℵ |011 + |110 + √ 8 + x 2 -x 2 |101 , ( 8
)
|w 2 (x) = ℵ |100 + |001 + √ 8 + x 2 -x 2 |010 , (9)
where

ℵ = 2 + ( √ 8 + x 2 -x)/2 2 -1/2 . ( 10
)
The relevant negativity can be evaluated and it turns out to be:

N (ρ 2-(x)) = 1 2 5/3 - 1 (8 + x 2 ) 3/2 x + √ 8 + x 2 -2(20 + x(x + √ 8 + x 2 )) × x -3 √ 8 + x 2 + 8 √ 8 + x 2 (|λ 1 (x)| + |λ 2 (x)| + |λ 3 (x)|) 2 1/3 (11)
with λ i (x)'s the roots of the algebraic equation in λ:

-10x + x 3 + (2 + x 2 ) √ 8 + x 2 + -32x -4x 3 -48 √ 8 + x 2 -4x 2 √ 8 + x 2 λ + 128x + 16x 3 -384 √ 8 + x 2 -48x 2 √ 8 + x 2 λ 2 + 1024 √ 8 + x 2 + 128x 2 √ 8 + x 2 λ 3 = 0 , (12)
where the coefficients, and hence the solutions, depend on x. The behavior of N (ρ 2-(x))

is shown in fig. 6, where the negativity of the thermal state for c = 6ω 0 and T ≈ 0 versus x has also been plotted. It is well visible that the behavior of N (ρ 2-(x)) perfectly reproduces the edge of the top visible in fig. 5 in the low temperature region and for x belonging to the region where Ψ - 2 is the ground state of the system.

It could be interesting to compare the behavior of the tripartite negativity with the values of the three concurrences related to the three couples of spins that can be extracted tracing over one of them: negativity (and hence all the three relevant bipartite Negativity functions) is not. This apparent contradiction reflects the very different meanings of these quantities. Indeed, for example, while the concurrence C 12 measures the entanglement between spins 1 and 2, without considering the spin 3 that has been preliminarily traced over, the negativity N 1-23 takes into account the correlations between spin 1 and the whole couple made of the two spins 2 and 3. Therefore, in that region, we can talk about correlations between spin 1 and the couple (2,3) without correlations between 1 and any of the other two spins, after the third has been traced over. We think this is an important result, since it underlines the difference between N and C. To conclude the analysis of the inhomogeneous model of type A, in figure 9 we have plotted the complete dependence of N 123 on coupling and inhomogeneity parameters, fixing the value of temperature at k B T = 0.01ω 0 . The dependence on the coupling parameter is qualitatively the same as for the previous case (see figure 8). Indeed, we observe a step trend caused by a level crossing between the two lowest energy levels.

C 23 = C(tr 1 ρ 2-), C 12 = C(tr 3 ρ 2-), C 13 = C(tr 2 ρ 2-).
As a function of x, N 123 exhibits both smooth decreasing, due to increasing separability of spins 1 and 3, and fast variations traceable back to a level crossing amongst the two lowest energy levels E - 2 and E - 4 . In fact, we verified that the big step observed in figure 9 corresponds to the locus of points of the plane (c, x) in which the two lowest levels have the same energy. Thanks to the simple form of the eigenvalues of our model, one can find a simple analytical expression of this curve:

c(x) = x 2 + 4 + α(x) 4x + 1 2 x 2 + 5 + α(x) 2 , (13)
where α(x) is

α(x) = √ 16 + 10x 2 + x 4 . ( 14
)

Model with inhomogeneity of type B

The

Conclusions

In this paper we have analyzed the thermal correlations in a spin-star system consisting of a central spin interacting with three outer ones, all immersed in a magnetic field.

In the absence of a tripartite entanglement measure, we have decided to exploit the tripartite negativity, which is at least able to put in evidence the lack of separability and simple bi-separability of a mixed quantum state of a tripartite system.

We started from considering the homogeneous model in which the three peripheral spins interact in the same way with the central one, showing the appearance of a significant degree of inseparability revealed by a non vanishing tripartite negativity N 123 .

Typical behavior of thermal correlations reflects the behavior of N 123 , which decreases to zero-value as temperature increases and, instead, exhibits both appreciable values and abrupt changes at very low temperature. These occurrences show in a very clear way the role of thermal entanglement mediator played by the central spin. Further we have considered the presence of inhomogeneity, meant as differences in the coupling constants between the outer spins and the central one, focusing on two special cases. In the first case (inhomogeneity of type A) two spins interact in the same way with the central one and the third has a different coupling constant. In the second case (inhomogeneity of type B) all the three outer spins are characterized by different coupling strengths.

Though some differences are visible, the qualitative behaviors of N 123 are quite similar for the two inhomogeneous models and even for the homogeneous one. This suggests the idea that the thermal quantum correlations mediated by the central spin are not very much damaged by a certain lack of homogeneity that could characterize a more realistic situation, provided the degree of inhomogeneity is not high.

A remarkable point is that we have singled out the presence of maxima of tripartite negativity corresponding to inhomogeneous models, i.e. for the values of the inhomogeneous parameter quite larger than 1. Though the differences between a maximum value of negativity and its values in the nearby (even up to the homogeneous case, i.e. up to x = 1) is not very big, we think that it is conceptually important the fact that a higher degree of symmetry in the system does not guarantee a higher degree of correlations between all its parts.

Another important point is that we have found some regions of the parameter space where it happens that two concurrences, say C 12 and C 23 , are vanishing, meaning that there is no entanglement between 1 and 2 and between 2 and 3, while the relevant negativity N 2-13 is non-vanishing, meaning that there is a correlation between 2 and the couple made of 1 and 3. This fact points out in a very transparent way the different meanings of concurrence and negativity.

Appendix A.

In this appendix we give eigenvalues and eigenvectors of the Hamiltonian of the inhomogeneous model of type A (c 1 = c, c 2 = c x, c 3 = c), as functions of x > 0, c > 0, ω 0 . The special case x = 1 gives the solutions for the homogeneous model.

Energies Eigenstates

E ± 1 = ±c x , Ψ ± 1 = 1 2 [(|0011 ± |1100 ) -(|0110 ± |1001 )] , E ± 2 = ± c 2 x + (8 + x 2 ) 1 2 , Ψ ± 2 = 1 K 1 [(|0011 ± |1100 ) + (|0110 ± |1001 ) + √ 8+x 2 -x 2 (|0101 ± |1010 ) , E ± 3 = ± c 2 x -(8 + x 2 ) 1 2 , Ψ ± 3 = 1 K 1 [(|0011 ± |1100 ) + (|0110 ± |1001 ) - √ 8+x 2 +x 2 (|0101 ± |1010 ) , E ± 4 = ± c (2 + x 2 ) 1 2 + ω 0 , Ψ ± 4 = 1 K 2 √ 2 + x 2 |0111 ± (|1011 + x |1101 + |1110 )] , E ± 5 = ± c (2 + x 2 ) 1 2 -ω 0 , Ψ ± 5 = 1 K 2 [(|0100 + x |0010 + |0001 ) ± (2 + x 2 ) |1000 , E 6 = -ω 0 , Ψ A 6 = 1 K 3 1 x |0001 + |0010 -1 x + x |0100 , Ψ B 6 = 1 √ 1+x 2 (|0010 -x |0001 ) , E 7 = ω 0 , Ψ A 7 = 1 K 3 1 x |1011 + |1101 -1 x + x |1110 , Ψ B 7 = 1 √ 1+x 2 (|1101 -x |1011 ) , E 8 = -2ω 0 , |Ψ 8 = |0000 , E 9 = 2ω 0 , |Ψ 9 = |1111 ,
where K 1 , K 2 and K 3 are:

K 1 2 = 4 + 2 √ 8 + x 2 -x 2 2 , K 2 2 = 2 2 + x 2 , K 3 2 = 2 x 2 + 3 + x 2 .

  . Numbers 1, 2, 3 indicate the three outer spin 1/2, while the central one is indicated with latin capital letter C. With c 1 , c 2 and c 3

Figure 1 .

 1 Figure 1. The physical system. Spins labeled 1,2 and 3 are coupled to the central (C) one. The magnetic field B 0 is ortogonal to the plan where the four spins lie.

Figure 2 .Figure 3 .

 23 Figure 2. Tripartite negativity N 123 for the homogeneous model as a functions of coupling constant (c) and temperature (k B T ), both in units of ω 0 . Both the diminishing for increasing temperature and the low-temperature transitions with respect to c are well visible.

Figure 4 .

 4 Figure 4. Tripartite negativity N 123 of the homogeneous model versus coupling constants (in units of ω 0 ), at k B T = 0.01ω 0 . Abrupt variations at c ≈ 0.6ω 0 and c ≈ 3.7ω 0 are well visible. Three plateaux, A 1 (c < 0.6ω 0 ), A 2 (0.6ω 0 < c < 3.72ω 0 ), A 3 (c > 3.72ω 0 ) are quite evident.

Figure 5 .

 5 Figure 5.Negativity N 123 versus temperature (k B T , in units of ω 0 ) and inhomogeneity parameter (x), at c = 6ω 0 , for the inhomogeneous model of type A.

  A 3 of figure 4 negativity assumes its greatest value, we start the analysis of the inhomogeneous model of type A showing in figure 5 the complete dependence of N 123 on temperature and inhomogeneity parameter for c = 6ω 0 (this specific value of c belongs indeed to the area A 3 of the homogeneous model). It is quite obvious that the physical reason of the decreasing of N 123 with respect to T

  (|1001 + |0110 ) -0.6435 (|1010 + |0101 ) for x = x M . In this region, the values of N 123 and in particular the appearance of a maximum are essentially determined by the dependence of the structure of the state Ψ - 2 on the inhomogeneity parameter. The relevant density operator is a mixture of two Werner-like states, and can be written in the following way:

Figure 7 Figure 6 .

 76 Figure 7 shows that around x = 0.5 and x = 5.5 there are abrupt changes in the values of the concurrences (in the same points where the tripartite negativity exhibits the same behavior), due to rapid changes of the ground state. It is worth noting that in the region 3.5 -5.5 the concurrences C 12 and C 23 are vanishing while the tripartite

Figure 8 Figure 7 .Figure 8 .

 878 Figure8shows the complete dependence of N 123 on temperature and coupling parameter, for x = 3. We observed (by performing other plots, which we are not reporting here) that the qualitative behavior of N 123 is independent on the specific value of x, provided it is neither too small nor too big. It is quite evident that the trend with respect to temperature is the usual one. At low temperature the dependence of N 123 on c exhibits only the step trend due to level crossings. The highest step involves the energy levels E - 4 (c, x = 3) = -c √ 11 -ω 0 and E - 2 (c, x = 3) = -3 + √ 17 c/2 (see

Figure 9 .

 9 Figure 9. Negativity N 123 , at low temperature (k B T = 0.01ω 0 ), versus coupling parameter (c, in units of ω 0 ) and inhomogeneity parameter (x), for the inhomogeneous model of type A. Both smooth decreasing and fast variations occur with respect to x.Dependence on c exhibits a step trend.

  second type of inhomogeneity we decided to study is characterized by the coupling constants c 1 = c, c 2 = c x and c 3 = c x 2 . Also for this model we have considered the dependence on T , c and x. Surprisingly, what we have found is that the behavior of negativity in the presence of inhomogeneity of type B is not different from what we have found in the presence of inhomogeneity of type A. The diminishing for increasing temperature, the abrupt changes at low temperature in connection with level crossings and the presence of a maximum for x = 1 are all present also in this model of inhomogeneity. Since for type B we do not have an explicit analytical diagonalization of the Hamiltonian, we also miss an analytical expression of the state of the system in the nearby of the maximum. In fact, we have only numerical results, according to which, for example, for T = 0.01 ω 0 /k B and c = 6 ω 0 the maximum of negativity occurs at x = x M ≈ 1.91, where the ground state of the system is approximately: Ψ (B) M = 0.292 (|1100 -|0011 ) + 0.471 (|1001 -|0110 ) + 0.439 (|1010 -|0101 ).

Appendix A), and it is quite easy to find that the crossing occurs at c t ≈ 4.08ω 0 where