S Zerguine

P Van Isacker

Spin-aligned neutron-proton pairs in N = Z nuclei

Keywords: numbers: 21, 60, Cs, 21, 60, Ev

A study is carried out of the role of the aligned neutron-proton pair with angular momentum J = 9 and isospin T = 0 in the low-energy spectroscopy of the N = Z nuclei 96 Cd, 94 Ag, and 92 Pd. Shell-model wave functions resulting from realistic interactions are analyzed in terms of a variety of two-nucleon pairs corresponding to different choices of their coupled angular momentum J and isospin T . The analysis is performed exactly for four holes ( 96 Cd) and carried further for six and eight holes ( 94 Ag and 92 Pd) by means of a mapping to an appropriate version of the interacting boson model. The study allows the identification of the strengths and deficiencies of the aligned-pair approximation.

I. INTRODUCTION

The study of nuclei with equal numbers of neutrons and protons (N = Z) is one of the declared goals of radioactive-ion-beam facilities, either in operation or under construction. Such nuclei are well known when N and Z are small but, as the atomic mass number A = N +Z increases, they lie increasingly closer to the proton drip line and therefore become more difficult to study experimentally. Nevertheless, several phenomena of interest, such as the breaking of isospin symmetry or the emergence of new collective modes of excitation, are predicted to become more pronounced with increasing A, and this constitutes the main argument for undertaking the difficult studies of ever more heavier N = Z nuclei.

Arguably, the goal of most interest in this quest is the uncovering of effects due to isoscalar (T = 0) neutron-proton (n-p) pairing. In contrast to the usual isovector (T = 1) pairing, where the orbital angular momenta and the spins of two nucleons are both antiparallel (i.e., L = 0 and S = 0), isoscalar pairing requires the spins of the nucleons to be parallel (S = 1), resulting in a total angular momentum J = 1. Collective correlation effects are predicted to occur as a result of isoscalar n-p pairing [1] but have resisted so far experimental confirmation because (i) the jj coupling scheme which is applicable in all but the lightest nuclei disfavors the formation of isoscalar n-p pairs with L = 0 [2], and (ii) the states associated with this collective mode are of low angular momentum J and often hidden among high-J isomeric states which hinders their experimental detection.

Currently, N = Z experiments are approaching 100 Sn, involving studies of nuclei such as 92 Pd [3] where nucleons are dominantly confined to the 1g 9/2 orbit. In the context of these experiments, Blomqvist recently proposed [4] that a realistic description of shell-model wave functions can be obtained in terms of isoscalar n-p pairs which are completely aligned in angular momentum, that is, pairs with J = 9. This proposal is attractive since it is well adapted to the jj coupling scheme, valid in this mass region, and because it should encompass the description of high-J isomeric states. Blomqvist's idea is related to the socalled stretch scheme which was advocated a long time ago by Danos and Gillet [5]. In the stretch scheme, shell-model states are constructed from aligned n-p pairs, treated in a quasi-boson approximation which neglects antisymmetry between the nucleons in different pairs. The latter approximation is absent from Blomqvist's approach.

In this paper we examine Blomqvist's proposal for the N = Z nuclei 96 Cd, 94 Ag, and 92 Pd. We consider several realistic two-body interactions for the 1g 9/2 orbit and analyze the shell-model wave functions, obtained with these interactions, of the four-nucleon-hole system ( 96 Cd), in terms of a variety two-pair states. For the six-and eight-hole nuclei ( 94 Ag and 92 Pd) a direct shell-model analysis in terms of pair states is more difficult, and we prefer therefore to carry out an indirect check by means of a mapping to a corresponding boson model. In these cases our approach is intermediate between that of Blomqvist [4] and of Danos and Gillet [5]. The boson mapping takes care of antisymmetry effects in an exact manner on the level of four nucleons but becomes approximate for more. This paper is organized as follows. First, some necessary concepts and techniques are introduced: the formulas needed to carry out a shell-model calculation in a pair basis are given in Sect. II and two mapping techniques from an interacting fermion to an interacting boson model are reviewed in Sect. III. The results of our analysis of N = Z nuclei are presented and discussed in Sect IV. The conclusions and outlook of this work are summarized in Sect. V.

II. FOUR-PARTICLE MATRIX ELEMENTS

This section summarizes the necessary ingredients to carry out calculations in an isospin formalism. The formulas given are valid for fermions as well as for bosons. Four-particle states are described by grouping the particles in two pairs. These two-pair states can be used, for example, as a basis in a shell-model calculation, facilitating the subsequent analysis of the pair structure of the eigenstates. Furthermore, the two-pair representation of fourparticle states is the natural basis to map the shell model onto a corresponding model in terms of bosons. Once this mapping is carried out, the original interacting fermion problem is reduced to one of interacting bosons, which can also be solved with the results summarized in this section.

In the pair representation of four particles with angular momentum j and isospin t (both integer for bosons and half-odd-integer for fermions), a state can be written as

|(jt) 2 (J 1 T 1 )(jt) 2 (J 2 T 2 )
; JT where particles 1 and 2 are coupled to angular momentum and isospin J 1 T 1 , particles 3 and 4 to J 2 T 2 , and the intermediate quantum numbers J 1 T 1 and J 2 T 2 to total JT . It is convenient to introduce a short-hand notation that takes care simultaneously of angular momentum and isospin quantum numbers and we denote JT as Γ; indices will always be carried over consistently, i.e., Γ i refers to J i T i , γ to jt, etc.

The state |(jt) 2 (J 1 T 1 )(jt) 2 (J 2 T 2 ); JT is then denoted as |γ 2 (Γ 1 )γ 2 (Γ 2 ); Γ . This state is not (anti)symmetric in all four bosons (fermions) but can be made so by applying the (anti)symmetrization operator P ,

|γ 4 [Γ 1 Γ 2 ]Γ ≡ P |γ 2 (Γ 1 )γ 2 (Γ 2 ); Γ = ΓaΓ b [γ 2 (Γ a )γ 2 (Γ b ); Γ|}γ 4 [Γ 1 Γ 2 ]Γ]|γ 2 (Γ a )γ 2 (Γ b ); Γ , (1)
where

[γ 2 (Γ a )γ 2 (Γ b ); Γ|}γ 4 [Γ 1 Γ 2 ]Γ] is a four-to-two coefficient of fractional parentage (CFP).
The notation in square brackets [Γ 1 Γ 2 ] implies that ( 1) is constructed from a parent state [START_REF] Talmi | Simple Models of Complex Nuclei. The Shell Model and Interacting Boson Model[END_REF] with intermediate angular momenta and isospins J 1 T 1 and J 2 T 2 . It is implicitly assumed that all pair quantum numbers are allowed, i.e. that J i is even for T i = 1 and odd for T i = 0.

The four-to-two CFP is given by

[γ 2 (Γ a )γ 2 (Γ b ); Γ|}γ 4 [Γ 1 Γ 2 ]Γ] = 1 N Γ 1 Γ 2 γΓ       δ Γ 1 Γa δ Γ 2 Γ b + (-) J+T δ Γ 1 Γ b δ Γ 2 Γa + 4σ       γ γ Γ a γ γ Γ b Γ 1 Γ 2 Γ             . (2)
where σ = +1 for bosons and σ = -1 for fermions. Furthermore,

N Γ 1 Γ 2 γΓ is a normalization constant, δ Γ 1 Γa ≡ δ J 1 Ja δ T 1 Ta , δ Γ 2 Γ b ≡ δ J 2 J b δ T 2 T b , and       γ γ Γ a γ γ Γ b Γ 1 Γ 2 Γ       = [Γ 1 ][Γ 2 ][Γ a ][Γ b ]            j j J a j j J b J 1 J 2 J                       t t T a t t T b T 1 T 2 T            , (3)
where the symbol in curly brackets is a nine-j symbol and [Γ i ] = (2J i + 1)(2T i + 1). The normalization constant is known in closed form as

N Γ 1 Γ 2 γΓ = 6       1 + (-) J+T δ Γ 1 Γ 2 + 4σ       γ γ Γ 1 γ γ Γ 2 Γ 1 Γ 2 Γ             . ( 4
)
The states (1) do not form an orthonormal basis. One needs therefore to determine the overlap matrix elements which can be written in terms of the CFPs as

γ 4 [Γ 1 Γ 2 ]Γ|γ 4 [Γ ′ 1 Γ ′ 2 ]Γ = ΓaΓ b [γ 2 (Γ a )γ 2 (Γ b ); Γ|}γ 4 [Γ 1 Γ 2 ]Γ][γ 2 (Γ a )γ 2 (Γ b ); Γ|}γ 4 [Γ ′ 1 Γ ′ 2 ]Γ] = 6 N Γ 1 Γ 2 γΓ [γ 2 (Γ 1 )γ 2 (Γ 2 ); Γ|}γ 4 [Γ ′ 1 Γ ′ 2 ]Γ]. (5)
Furthermore, the matrix elements of the two-body part of the Hamiltonian can be expressed as

γ 4 [Γ 1 Γ 2 ]Γ| Ĥ2 |γ 4 [Γ ′ 1 Γ ′ 2 ]Γ = 6 ΓaΛ [γ 2 (Γ a )γ 2 (Λ); Γ|}γ 4 [Γ 1 Γ 2 ]Γ][γ 2 (Γ a )γ 2 (Λ); Γ|}γ 4 [Γ ′ 1 Γ ′ 2 ]Γ]ν Λ , (6)
where ν Λ are two-body matrix elements between normalized two-particle states, ν Λ ≡ γ 2 ; Λ| Ĥ2 |γ 2 ; Λ . The label Λ is a short-hand notation for the two-particles' angular momentum λ and isospin T . For example, in the nuclear shell model the particles are nucleons with half-odd-integer angular momentum j and isospin t = 1/2. We recall that in this case the antisymmetric two-nucleon states are uniquely determined by the total angular momentum λ and that the total isospin is a redundant quantum number which is T = 0 for odd λ and T = 1 for even λ. The two-body matrix elements therefore depend on λ only, ν λ = j 2 ; λ| Ĥ2 |j 2 ; λ . In the interacting boson model, in contrast, the particles are bosons with integer angular momentum j and integer isospin t = 0 or t = 1. In the latter case, the two-particles' isospin T (obtained by coupling t = 1 with t = 1) is not redundant but is needed to fully characterize the two-particle state.

III. METHODS OF MAPPING

The basic and common idea of different boson mapping methods is to truncate the full shell-model space to a subspace which is written in terms of fermion pairs and to establish subsequently a correspondence between the fermion-pair space and an analogous space which is written in terms of bosons. This section recalls briefly two methods, namely the OAI (Otsuka-Arima-Iachello) [START_REF] Otsuka | [END_REF] and the democratic [8] mappings that will be used in the applications.

A. The OAI mapping

We limit ourselves here to the problem of n nucleons with isospin t = 1/2 in a single orbit j which is the case of interest in this paper. A detailed description of the method in its full generality can be found in Ref. [START_REF] Otsuka | [END_REF].

Let us first introduce the pair creation operators

A †(JT )

M J M T = 1 √ 2 [a † × a † ] (JT ) M J M T , (7)
where a † is the creation operator of a t = 1/2 nucleon in orbit j and the square brackets denote coupling to a tensor with angular momentum J and isospin T , and projections M J and M T , respectively. For nucleons in a single orbit j, the pair ( 7) is totally determined by its angular momentum J and its isospin follows from it (i.e., T = 0 for odd J and T = 1 for even J). The starting point of the OAI mapping is a given shell-model Hamiltonian ĤF which is a scalar in J and T , and the selection of a number of fermion pairs J 1 , J 2 , . . . , J p , to each of which will be associated a boson with the same angular momentum J i . (For the bosons this angular momentum is by convention denoted as ℓ i , ℓ i = J i .) For a single-j orbit one may ignore the single-particle term of the Hamiltonian which we therefore take to be of pure two-body character and denote it as ĤF 2 . The one-body term of the mapped boson Hamiltonian ĤB is directly obtained from the matrix elements of the shell-model

Hamiltonian ĤF

2 between the pair states,

ǫ ℓ i ≡ ℓ i | ĤB |ℓ i = J i | ĤF 2 |J i . (8)
These matrix elements are interpreted as single-boson energies.

The determination of higher-body terms of the mapped boson Hamiltonian is more complicated, as we illustrate here for the two-body part. For a given angular momentum J

and isospin T , one first enumerates all possible two-pair states |F i ≡ |J a J b ; JT , where the index i = 1, . . . , d corresponds to pairs (a, b) with 1 ≤ a ≤ b ≤ p. These will be mapped onto the two-boson states

|B i ≡ |ℓ a ℓ b ; JT , |F i -→ |B i , i = 1, . . . , d. (9)
This correspondence cannot be established directly since the boson states are orthogonal while the fermion states are not:

B i |B j = δ ij , o ij ≡ F i |F j = δ ij , (10)
where o ij is the overlap matrix in fermion space. To deal with this problem, the strategy of the OAI mapping is to define for each JT a hierarchy of pair states, to apply a Gram-Schmidt orthogonalization of this ordered sequence, and to associate the boson states with the orthonormalized fermion states. This leads to the definition of the correspondence

| Fi =⇒ |B i , i = 1, . . . , d, (11)
where the states | Fi form an orthonormalized basis. One has the series

| F1 = 1 √ o 11 |F 1 , | F2 = 1 o 22 -(õ 21 ) 2 (|F 2 -õ21 | F1 ), . . . | Fk = 1 √ N k |F k - k-1 i=1 õki | Fi , (12)
until k = p, with

N k = o kk - k-1 i=1 (õ ki ) 2 , õki ≡ F k | Fi . ( 13
)
An efficient algorithm to carry out the Gram-Schmidt procedure involves expanding the orthonormal states | Fk in the non-orthogonal basis |F i as

| Fk = k i=1 a ki |F i , (14)
and calculating the coefficients a ki (needed for k ≥ i) recursively from

a kk = 1 √ N k , a ki = - 1 √ N k k-1 i ′ =i õki ′ a i ′ i , õki = i i ′ =1 a ii ′ o ki ′ , (15)
which as only input requires the knowledge of the overlap matrix elements o ij . The matrix elements of the fermion Hamiltonian in the orthonormal basis | Fk are then given by

Fk | ĤF 2 | Fl = ij a ki a lj F i | ĤF 2 |F j , (16)
and the mapped boson Hamiltonian follows from

ℓ a ℓ b ; JT | ĤB |ℓ a ′ ℓ b ′ ; JT = J a J b ; JT | ĤF 2 | J a ′ J b ′ ; JT . ( 17
)
Note that this equation defines the entire boson Hamiltonian up to and including two-body interactions. To isolate its two-body part ĤB 2 , one should subtract the previously determined one-body terms according to

ℓ a ℓ b ; JT | ĤB 2 |ℓ a ′ ℓ b ′ ; JT = ℓ a ℓ b ; JT | ĤB |ℓ a ′ ℓ b ′ ; JT -(ǫ ℓa + ǫ ℓ b )δ aa ′ δ bb ′ , (18)
always assuming that a ≤ b and a ′ ≤ b ′ . This is the version of the OAI mapping as it will be applied in this paper. The boson Hamiltonian is obtained from the two-and four-particle systems and is kept constant for systems with higher numbers of particles. This is akin to the democratic mapping but the latter has the advantage that no hierarchy of two-pair states is required, as is discussed in the next subsection.

Nevertheless, the hierarchy imposed in the OAI mapping is often based on arguments of seniority which enables an extension to n-particle systems. To illustrate this point, we consider the case of a truncation onto a subspace constructed out of pairs of angular momenta J = 0 (S pair) and J = 2 (D pair). (In this example we assume for simplicity identical nucleons so that isospin can be omitted.) The creation operators S † and D † are defined as

S † = A †(0) 0 , D † M = A †(2) M . (19)
The truncated shell-model space is the SD subspace spanned by the states

|j n vξ; JM J = 1 N F P(S † ) (n-v)/2 [(D † ) v/2 ] (J) M J |o , (20)
where |o refers to the closed core, ξ denotes additional quantum numbers related to the intermediate angular-momentum couplings of the D pairs, and N F is a normalization constant. The operator P is needed to ensure orthogonality of the basis. Its effect is a projection onto a subspace with seniority v, the number of particles not in pairs coupled to J = 0.

As a result the states (20) are normalized and can be mapped onto sd states. This can be expressed by the general correspondence

|j n vξ; JM J =⇒ |n s n d ξ; JM J = 1 N B (s † ) ns [(d † ) n d ] (J) M J |o), (21)
where |o) is the boson vacuum, n s and n d are the sand d-boson numbers, respectively, with

n s = (n -v)/2 and n d = v/2,
n s n d ξ; JM J | ĤB |n ′ s n ′ d ξ ′ ; JM J = j n vξ; JM J | ĤF 2 |j n v ′ ξ ′ ; JM J . (22)
The difference with Eq. ( 17) is that the latter equation is specific to n = 4 particles while the mapping (22) applies to any n. The calculation of an n-particle fermion matrix element is possible because the seniority formalism allows it to be reduced to n = max(v, v ′ ) [START_REF] Talmi | Simple Models of Complex Nuclei. The Shell Model and Interacting Boson Model[END_REF].

Hence, as long as v, v ′ ≤ 4 the fermion matrix element can be computed and fixes the corresponding boson matrix element. The advantage of this procedure is that it determines a boson Hamiltonian which varies with particle number yielding a more adequate description of Pauli effects. Its disadvantage is that a generalization towards any choice of pairs and/or to neutrons and protons is difficult, if not impossible. Therefore, we opt in this paper for the simpler OAI mapping that defines a constant boson Hamiltonian from the two-and four-particle systems.

B. The democratic mapping

Again the starting point is the selection of a number of fermion pairs J 1 , J 2 , . . . , J p and a corresponding series of bosons with energies ǫ ℓ i determined from Eq. ( 8). To establish the correspondence ( 9), complicated by the non-orthogonality of the fermion two-pair states, the democratic mapping relies on the diagonalization of the overlap matrix o ij defined in Eq. ( 10). This diagonalization provides, besides the eigenvalues o k , an orthogonal basis |X k ,

|X k = i c ki |F i , k = 1, . . . , d, (23)
with

X k |X l = ij c ki c lj o ij = o k δ kl , k, l = 1, . . . , d. (24)
We follow the convention of labeling non-orthogonal basis states by i or j, and orthogonal basis states by k or l. The coefficients

c ki ≡ F i |X k ≡ X k |F i are
The states |X k form an orthogonal but non-normalized basis since they satisfy [START_REF] Devi | [END_REF]. To normalize these states, we define

| Fk = 1 √ o k |X k , k = 1, . . . , d, (26)
which indeed satisfy Fk | Fl = δ kl , k, l = 1, . . . , d. (27)
The basis | Fk has the same role as the basis | Fk in the OAI mapping but it is important to realize that both bases are not identical and lead to different boson Hamiltonians. The democratic mapping relies on the definition of a transformation in boson space which is analogous to the one in fermion space

| Bk = i c ki |B i , k = 1, . . . , d. (28)
From the orthogonality of the basis |B i and the properties of the coefficients c ki it can be shown that these states form an orthonormal set,

Bk | Bl = δ kl , k, l = 1, . . . , d. (29)
We have now arrived at fermion and boson bases that are both orthonormal, and we can therefore establish the mapping

| Fk =⇒ | Bk , k = 1, . . . , d, (30)
and determine the matrix elements of the boson Hamiltonian in this basis,

Bk | ĤB | Bl = Fk | ĤF 2 | Fl , k, l = 1, . . . , d. (31)
With use of the inverse of the relation ( 28), of the equality (31), and of Eqs. ( 23) and ( 26), the matrix elements of the boson Hamiltonian in the original basis can be written in terms of those of the fermion Hamiltonian (also in the original basis) as

B i | ĤB |B j = kl i ′ j ′ 1 √ o k o l c ki c ki ′ c lj c lj ′ H F i ′ j ′ , (32)
where

H F i ′ j ′ ≡ F i ′ | ĤF 2 |F j ′ .
Again, as with the OAI mapping, this defines the entire boson Hamiltonian from which the two-body part can be isolated by applying Eq. (18).

In some cases, m fermion vectors |F i are linearly dependent on the dm others, leading to m vanishing eigenvalues of the overlap matrix,

o i = 0, i = 1, . . . , m. (33)
This problem can be solved by excluding from the fermion space m states |F i and by calculating the matrix elements of ĤB from the remaining dm states. All other matrix elements of ĤB are defined such that

B i | ĤB |B i = ∞, i = d -m + 1, . . . , d, B i | ĤB |B j = B j | ĤB |B i = 0, i = j, (34)
where i and j in the last equation take the values i = 1, . . . , d and j = dm + 1, . . . , d. One has still to decide which m vectors to remove from the fermion space. To avoid arbitrariness, one chooses the m fermion states with the smallest overlap with the shell-model states.

IV. APPLICATION TO N = Z NUCLEI

Of particular interest in this work are N = Z nuclei with neutrons and protons in the 1g 9/2 orbit. Recently, Blomqvist [4] has conjectured that a valid interpretation of yrast states in these nuclei can be obtained in terms of n-p T = 0 pairs which are coupled to maximum angular momentum J = 9 and which therefore can be termed aligned n-p pairs. This is contrary to the usual interpretation of such states which involves low-spin pairs with isospin T = 1 and possibly also with T = 0. In this section we examine Blomqvist's proposal with specific reference to the nuclei 96 Cd, 94 Ag, and 92 Pd, corresponding to four, six, and eight holes with respect to the 100 Sn core, respectively. Our study consists of two separate parts: the analysis of shell-model wave functions of 96 Cd in terms of a variety of two-fermion pairs and the mapping of shell-model onto corresponding boson states for 96 Cd and 92 Pd, following the formalism developed in the previous sections. First, an appropriate shell-model interaction should be determined.

A. Shell-model interaction

For the purpose of checking the stability of our results and the reliability of our conclusions, we have carried out the analysis for three different shell-model interactions. The SLGT0 interaction is taken from Serduke et al. [9]. It was used in the more recent analysis of Herndl and Brown [10] where it was found to give satisfactory results for the neutrondeficient nuclei in the mass region A = 86 to 100, which are of interest in the present study. A second shell-model interaction, which shall be named GF, is taken from Gross and Frenkel [11]. Both the SLGT0 and GF interactions are defined in the 2p 1/2 + 1g 9/2 shell-model space and, to carry out an analysis in terms of aligned pairs, it is necessary to renormalize them to the 1g 9/2 orbit. The resulting two-body matrix elements are shown in columns 2 and 3 of Table I. The renormalization only affects the (λ, T ) = (0, 1) and (1, 0) matrix elements since these are the only ones that also occur in the 2p 1/2 orbit.

To avoid the renormalization procedure, we define a third interaction directly for the 1g 9/2 orbit. This is done in the following way. The spectrum of 90 Nb is well known [12] and enables the determination of the particle-hole interaction matrix elements The isospin T is determined from λ, that is, T = 4 for all states except for λ = 0 which

ν ph λ ≡ 1g 9/2 (1g 9/2 ) -1 ; λT | ĤF 2 |1g 9/2 (1g 9/2 ) -1 ; λT . ( 35
)
where E stands for the nuclear binding energy, taken from the 2003 atomic mass evaluation [13] and corrected for the electrons' binding energy [14]. (The minus sign is needed to convert from binding energy to interaction energy.) All levels λ = 0, . . . , 9 of the particlehole multiplet (35) are known in 90 Nb and fix the differences ν ph λν ph 8 . A problem with this procedure concerns the choice of the relevant J π = 1 + state in 90 Nb. Several of them are observed at low excitation energy and only one should be taken as a member of the particle-hole multiplet. We have taken the J π = 1 + level at 2.126 MeV because this state is strongly populated in the β-decay of 90 Mo (log f t = 4.9) [15]. Furthermore, it appears to be the only 1 + state observed in the triton spectrum obtained in the charge-exchange reaction 90 Zr( 3 He, t) 90 Nb [16].

Once the particle-hole matrix elements ν ph λ are derived in this way, the particle-particle (or hole-hole) matrix elements ν pp λ are obtained from an inverted Pandya tranformation [17]. Since absolute (as opposed to relative) matrix elements have been extracted from binding energies, care should be taken to use the appropriate Pandya transformation. Equation (18.63) of Ref. [START_REF] Talmi | Simple Models of Complex Nuclei. The Shell Model and Interacting Boson Model[END_REF] gives the following relation between absolute matrix elements:

ν ph λ = E 0 - λ ′ (2λ ′ + 1) j j λ j j λ ′ ν pp λ ′ , (37)
where E 0 is a constant (i.e., λ-independent) interaction energy given by

E 0 = 1 2j + 1 all λ ′ (2λ ′ + 1)ν pp λ ′ + 2j -1 2j + 1 even λ ′ (2λ ′ + 1)ν pp λ ′ . (38)
To express the particle-particle in terms of the particle-hole matrix elements, as is needed here, the relation (37) can now be inverted in the usual manner, leading to

ν pp λ = E 0 - λ ′ (2λ ′ + 1) j j λ j j λ ′ ν ph λ ′ . (39)
The problem with this relation is that it is of no use as long as the constant E 0 is expressed in terms of particle-particle matrix elements, as in Eq. (38). We need an equation for E 0 in terms of particle-hole matrix elements. This can be obtained by inserting the expression (39) for ν pp λ in Eq. (38) and solving for E 0 , leading to

E 0 = 1 2j(2j + 1) all λ ′ (2λ ′ + 1)ν ph λ ′ - 2j -1 2j(2j + 1) ν ph 0 . (40)
The resulting particle-particle matrix elements are shown in column 4 of Table I under 'Nb90', where the 'pp' index is omitted, ν λ ≡ ν pp λ , as will be done from now on. Note that the matrix elements thus obtained differ by a constant from those of Sorlin and Porquet (see Fig. 6 of Ref. [18]), since we have taken care here of the constant interaction energy E 0 .

It is seen from Table I that the matrix elements of the SLGT0 and GF interactions are similar. The biggest difference concerns the (λ, T ) = (9, 0) matrix element which is more attractive by ∼ 250 keV for GF. This conceivably might influence the approximation in terms of aligned λ = 9 n-p pairs. The interaction derived from 90 Nb is less attractive (or more repulsive), in particular the pairing matrix element with (λ, T ) = (0, 1). All three Before proceeding with the wave-function analysis, we first check to what extent the N = Z nuclei in this region can be described by confining nucleons to the 1g 9/2 orbit. This approximation should be reasonable for nuclei close to (south-west of) 100 Sn but will become increasingly poor as one approaches 80 Zr. The latter nucleus is known to be deformed [19] and hence a single (spherical) orbit will not suffice for a reliable description of nuclei in its neighborhood. To quantify the limitation to the 1g 9/2 orbit, we show in Table II With the formalism developed in Sect. II it is now possible to perform a shell-model calculation for four holes in the 1g 9/2 orbit and analyze the resulting wave functions in terms of two-pair states. We concentrate on the N = Z case which corresponds to two neutrons and two protons. For a given angular momentum J, the overlap matrix between all possible two-pair states |(jt) 4 [(J 1 T 1 )(J 2 T 2 )]JT is constructed using Eq. ( 5). The number d of linearly independent states JT is given by the number of non-zero eigenvalues of the overlap matrix. We then select d linearly independent but otherwise arbitrary two-pair states that shall be denoted in short as

| Fk = i cki |F i , k = 1, . . . , d, (41)
with cki = c ki / √ o k in the notation of Sect. III. This basis can now be used to construct the

energy matrix, Fk | ĤF 2 | Fl = ij cki clj F i | ĤF 2 |F j , (42)
to be computed with use of Eq. ( 6). The diagonalization of this matrix leads to the energy eigenvalues ǫ r , r = 1, . . . , d, with corresponding eigenvectors given by

|E r = k e rk | Fk = ki e rk cki |F i . (43)
Note that the energies ǫ r are independent of the choice of two-pair states |F i as long as the latter span the entire JT space. Although the above procedure might seem rather cumbersome for performing a four-hole shell-model calculation, it has the advantage of providing us directly with the pair structure of the shell-model eigenstates from the overlaps

F i |E r = kj e rk ckj o ij . ( 44
)
We now analyze the yrast eigenstates of the various interactions defined in Table I, that is, the quantities F i |E r 2 for r = 1 and a variety of two-pair states F i . Pairs with angular momentum J are generically denoted as P J and explicitly, following standard spectroscopic

TABLE III: Coefficients a λ [Γ 1 Γ 2 ]
Γ in the expansion (45) of the diagonal energies of the (9/2) 4 system for isospin T = 0. notation, as S, D, G, I, and K for J = 0, 2, 4, 6, and 8. Since this paper deals in particular with the aligned n-p pair with J = 9, we reserve for it the non-standard notation 'B' (for Blomqvist). The central result is shown in Fig. 1 which displays the quantity J 1 |B 2 ; J 2

J = 0 J = 2 J = 4 J = 6 J = 8 J = 10 λ [B 2 ] [S 2 ] [B 2 ] [SD] [B 2 ] [SG] [B 2 ] [SI] [B 2 ] [SK] [B 2 ] [DK] 0 1.
where |J 1 is the yrast eigenstate with angular momentum J and isospin T = 0 of the interactions SLGT0, GF, and Nb90. Most yrast states have a large overlap with B 2 , as was shown by Blomqvist, but this is conspicuously not the case for J ≈ 8. It seems as if the two aligned n-p pairs do not like to couple to a total angular momentum which equals their individual spins.

To acquire some insight in this finding, we note that simple expressions in terms of the two-body matrix elements are available for the diagonal energies of two-pair states,

E(γ 4 [Γ 1 Γ 2 ]Γ) ≡ γ 4 [Γ 1 Γ 2 ]Γ| ĤF 2 |γ 4 [Γ 1 Γ 2 ]Γ = λ a λ [Γ 1 Γ 2 ]Γ ν λ , (45)
where the coefficients a λ [Γ 1 Γ 2 ]Γ of relevance to the present discussion are shown in Table III. According to a recent paper by Talmi [20], they are non-negative rational numbers; since they involve ratios of rather large integers in this case, Table III gives numerical approximations. Of all the coefficients given in this table, the important ones have λ = 0, 1, and 9, because these are the multipolarities of the most attractive interaction matrix elements.

It is seen that the contributions of the aligned matrix element (λ = 9) to the energies Figure 1 shows that the overlaps J 1 |B 2 ; J 2 are very similar for the three interactions.

This finding is at the basis of the fact that the subsequent analysis gives consistent results

for the three interactions. While there can be significant differences in the shell-model results with the different interactions, the approximation in terms of aligned pairs is similar for the three interactions. In other words, if a particular shell-model state is well approximated in terms of aligned pairs for one interaction, it is so for the other two as well; if the approximation is less good, it is so for all three interactions. Although we have carried out the complete analysis for the three interactions SLGT0, GF, and Nb90, we will show in the following only the results of the former since it has a proven track record of satisfactorily reproducing the data in the mass region of interest [10].

In Table IV are shown the amplitudes in percentages for yrast eigenstates of the SLGT0 interaction with even J and T = 0, that is, the quantities 100 × F i |E r 2 for r = 1 and a variety of pair states |F i . The numbers illustrate that, at least at low and at high angular momentum J, the overlaps of the physical eigenstates with |B 2 ; J are more important than those with other pair combinations. The percentages shown in Table IV also illustrate the non-orthogonality of the two-pair basis. For the example, the J = 0 ground state has a 91% overlap with B 2 but also a 80% overlap with S 2 ; this can only be if the overlap

B 2 |S 2
itself is rather large.

C. Boson mappings

Ideally, one would like to perform a similar analysis of shell-model eigenstates for more than four nucleons. That is a challenging problem, however, since it requires the formulation of a nucleon-pair shell model [21,22] in an isospin-invariant formalism. In this paper we choose to extend our analysis toward higher hole number through the boson mapping techniques explained in Sect. III. It is important to stress that this approximation goes beyond the original proposal of Blomqvist since it involves an additional assumption of the boson character of the fermion pairs. The results presented in this subsection therefore do not directly address Blomqvist's conjecture.

Once the mapping is carried out for the two-and four-hole systems according to one of the two procedures described in Sect. III, a Hamiltonian is obtained in terms of the selected bosons ℓ 1 , ℓ 2 , . . . , ℓ p , which can then be applied to systems with two or more bosons. Such a description shall be referred to as ℓ 1 ℓ 2 . . . ℓ p -IBM, where IBM stands for interacting boson model [START_REF] Iachello | The Interacting Boson Model[END_REF]. Note that all versions of IBM thus obtained are isospin invariant; for example, the sd-IBM is in fact the IBM-3 of Elliott and White [25].

To compare the merits of different selections of fermion pairs, the following mappings are considered:

1. A single fermion pair B with J = 9, T = 0, leading to the b-IBM.

2. Two fermion pairs S and B with J = 0, T = 1 and J = 9, T = 0, leading to the sb-IBM.

3. Two fermion pairs S and D with J = 0 and 2, both with T = 1, leading to the sd-IBM.

4. Three fermion pairs S, D, and G with J = 0, 2, and 4, all with T = 1, leading to the sdg-IBM.

The first two cases are inspired by Blomqvist's conjecture, involving the aligned n-p pair, while the next two are the standard choice of the IBM [START_REF] Iachello | The Interacting Boson Model[END_REF] and its most frequently used extension which includes g bosons. (For a review on the latter, see Ref. [START_REF] Devi | [END_REF].)

The results obtained with the various boson Hamiltonians are compared with T = 0 eigenstates of the SLGT0 interaction for four, six, and eight nucleons in Tables V, VI, and VII, respectively. The numerical calculations have been performed with the codes example of the complication mentioned at the end of Sect. III and the solution given there should be applied. In this case it implies that the matrix element B 2 ; 18| ĤB 2 |B 2 ; 18 should be taken infinitely repulsive and it is only by adhering to this procedure that reasonable results can be obtained.

Not much is known experimentally about 94 Ag except for the presence of two isomers, with tentative spin-parity assignments 7 + (presumably the lowest T = 0 state) and 21 + at about 5.78(3) MeV above the 7 + [START_REF] Mukha | [END_REF]. The different shell-model interactions SLGT0, GF, and Nb90 all predict a 7 + as the T = 0 ground state, and a 21 + level at 6.464, 5.948, and 4.632 MeV, respectively. The b-IBM reproduces the shell-model result for the binding energies of these isomers to about 100 keV for the 7 + and less than that for the 21 + (see Table VI for the results of the SLGT0 interaction). This result is valid for the different shellmodel interactions: although the binding energies calculated with the three shell-model interactions vary by several MeV, in each case they are matched by the (appropriately mapped) b-IBM to within about 100 keV, indicating that the B pair incorporates most of the correlations for the 7 + and 21 + states.

The b-IBM results should be contrasted to those obtained with sd-IBM which fails completely to reproduce the spectroscopy of 94 Ag. This is not surprising since it is known from the work of Elliott and Evans [29] that IBM-3 cannot give a satisfactory description of oddodd nuclei which require the addition of isoscalar s and d bosons leading to the so-called IBM-4. While the latter is a realistic model when low-j shell-model orbits are involved (e.g., for sd-shell nuclei [30,31]), the present results seem to indicate that the mapping from a shell-model space with high-j orbits calls for the inclusion of an aligned isoscalar n-p pair with J = 2j. While b-IBM and sb-IBM are largely equivalent for the odd-odd nucleus 94 Ag, this is not the case for the even-even nucleus 92 Pd. As can be seen from Table VII, the s boson provides crucial additional correlation energy which brings the boson result close to its shellmodel equivalent. This nucleus was studied recently in a fusion-evaporation experiment [3].

The excitation energies of the yrast levels, calculated with the SLGT0 interaction, are in reasonable agreement with the observed values of 0.874, 1.786, and 2.535 MeV, respectively.

T F µ (E2) = e ν i∈ν r 2 i Y 2µ (θ i , φ i ) + e π i∈π r 2 i Y 2µ (θ i , φ i ), (46)
where the sums are over neutrons and protons, and each sum is multiplied with the appropriate effective charge. This operator can be written alternatively as a sum of an isoscalar operator, multiplied by (e ν + e π ), and an isovector operator, multiplied by (e νe π ). For the E2 transitions between T = 0 levels of interest here, only the former part contributes.

In second quantization, which is a convenient formalism for carrying out the mapping, the fermion E2 operator can be written as

T F µ (E2) = - 55 3π l 2 ho e ν (a † ν × ãν ) (2) µ + e π (a † π × ãπ ) (2) µ , (47)
where a † ρ creates a neutron (ρ = ν) or a proton (ρ = π) in the 1g 9/2 orbit, and ãjm = (-) j+m a j-m . Furthermore, the factor in front comes from the radial integral over harmonicoscillator wave functions (with length parameter l ho ) involving the 1g 9/2 orbit. The lowest-order bosonic image of the fermion E2 operator is defined by the diagonal (reduced) matrix element in the 9 + state of the 1n-1p system which is given by (1g 9/2 ) 2 ; 9

+ || T F (E2)||(1g 9/2 ) 2 ; 9 + = - 55 3π l 2 ho × 1330 187 (e ν + e π ). (48)
The E2 operator of the b-IBM is of the form

T B µ (E2) = e b (b † × b) (2) µ , (49)
and is necessarily of scalar character in isospin. Since the mapping implies the equality (1g 9/2 ) 2 ; 9 + || T F (E2)||(1g 9/2 ) 2 ; 9

+ = b|| T B (E2)||b , (50)
In the following, the factor 55/3π(e ν +e π )l 2 ho is divided out of all matrix elements, fermionic as well bosonic.

A first test is from E2 transitions between four-nucleon-hole states with T = 0. The shellmodel results obtained with the SLGT0 interaction, shown in Fig. 2, display a characteristic approximation for two ranges of angular momenta, namely J = 6 to 13 and J = 17 to 21. This conclusion agrees, at least qualitatively, with the one drawn on the basis of energies (see Table VII). Since three b bosons cannot couple to total angular momentum J = 2, this state is absent from b-IBM while it is present in the shell model (see Table VII). As a consequence, no 4 → 2 or 3 → 2 transitions occur in b-IBM (see Figs. 3 and4). These transitions exist in the shell model but it is rather striking that the calculated matrix elements are very small indeed.

Finally, in Fig. 5 are shown the E2 transitions between eight-nucleon-hole states with T = 0. A small depletion of the E2 strength calculated in the shell model is perceptible for J ≈ 8 and is absent in b-IBM. Apart from this deviation both calculations agree, indicating that the shell-model wave functions can be adequately represented in terms of a single b boson. We emphasize once more that, the number of states that can be written in terms of B pairs is small (of the order of 10), no effective charges are needed to arrive at the agreement in Fig. 5.

Let us now formulate an overall evaluation of these results. While the b-IBM gives in all cases an reasonable description of the ground-state binding energy, the addition of the standard s boson (with J = 0, T = 1) further improves the agreement. In fact, the energies obtained in sb-IBM (both with the OAI and democratic mappings) agree well with those of the shell-model levels except for (i) the J = 8 level in the four-hole system, (ii) low-J levels of the six-hole system, and (iii) levels with J = 14 to 16 in the six-hole system. The first discrepancy is obviously related to the small overlap of the J = 8 shell-model state with |B 2 ; J = 8 , noted in Table IV and explained in Sect. IV B. The second difference is also understandable since a correct description of the low-J states in the odd-odd nucleus 94 Ag requires the consideration of low-J pairs with T = 0 which have been omitted from the present mapping. The third deviation could be related to an unfavorable coupling of three B pairs to the angular momenta J = 14, 15, and 16, akin to the coupling of two B pairs to J = 8. The results as regards E2 transitions are consistent with what is concluded from the analysis of spectra.

V. CONCLUSIONS

We have shown in this paper that part of the low-energy spectroscopy of N = Z nuclei which have their valence nucleons confined to a single high-j orbit, can be represented in terms of an aligned isoscalar n-p B pair with J = 2j and is further improved by the inclusion of an isovector S pair coupled to J = 0. This was proven explicitly for a four-hole system and indirectly, through a mapping onto a corresponding boson system, for six and eight holes. Some deficiencies were found in this approach. A first concerns states of the four-hole system with angular momentum J ≈ 2j which turn out to be poorly approximated with just S and B pairs. A second deficiency is more generic and pertains to the low-J states in odd-odd N = Z nuclei, the description of which calls for the inclusion of isoscalar n-p pairs with low angular momentum. Nevertheless, it should be noted that the two isomers that have been observed so far in 94 Ag, (7 + ) and (21 + ), are adequately described in terms of b bosons.

These results were obtained for the 1g 9/2 orbit and for three different choices of twobody interaction. To what extent are they valid in general and can they be considered as representative of a system of neutrons and protons confined to a high-j orbit? In essence, two ingredients, geometry and dynamics, determine the outcome of the present pair analysis.

The geometry is defined by the CFPs and, provided j is not too small, is expected to evolve only slowly with j. (It would in fact be an exercise of some interest to perform the pair analysis in the limit of large j.) The dynamics is determined by the two-body interaction which in our study was varied significantly but within reasonable bounds. The matrix elements shown in Table I are typical of what is obtained for a residual interaction with a short-range, attractive character [32] and we may thus expect similar results when we move to orbits other than 1g 9/2 . This work calls for further studies. The pair analysis of the shell-model wave functions should be extended to higher particle numbers which can be achieved through an isospininvariant formulation of the nucleon-pair shell model. Consequently, the present results require further confirmation at higher particle number but one is tempted to conclude at this point that a realistic model can be formulated in terms of s and b bosons. Due to its simplicity, such a model could be of use to elucidate the main structural features of N ∼ Z nuclei in this mass region. These topics are currently under study.

We wish to thank Stefan Heinze for his help with the numerical calculations with

  transformation coefficients from the non-orthonormal basis |F i to the basis |X k and they satisfy i c ki c li = δ kl , k, l = 1, . . . , d.

has T = 5 .ph 8 =

 58 The absolute value of the matrix element (35) for the ground state (λ = 8) is obtained from the binding energies of the surrounding nuclei as ν -E( 90 Zr) + E( 90 Nb) -E( 91 Nb) -E( 89 Zr) = -783.794 -776.895 + 788.942 + 771.825 = 0.078 MeV,

2 ; J 2 FIG. 1 :

 221 FIG.1: Overlaps of the (1g 9/2 ) 4 yrast eigenstates of the interactions SLGT0, GF, and Nb90 with angular momentum J and isospin T = 0, with the two-pair state |B 2 ; J . matrix o ij ≡ F i |F j allows the definition of an orthonormal basis

  of the states |(9/2) 4 [B 2 ]J and |(9/2)4 [SP J ]J remain more or less constant, independent of J. This is not the case for the pairing matrix element (λ = 0) whose contribution to the energy of |(9/2) 4 [B 2 ]J disappears as J increases while it remains important for |(9/2) 4 [SP J ]J . The combined effect of these contributions is that the state |(9/2) 4 [SP J ]J dips below |(9/2) 4 [B 2 ]J around J ≈ 8 and as a result picks up the largest component of the yrast eigenstate. Hence, the feature of the disappearing B dominance around J ≈ 8 is explained by a combination of geometry-the CFPs in the j = 9/2 orbit, and dynamics-the dependence of the interaction matrix elements on λ.

FIG. 2 :

 2 FIG.2: Absolute values of the E2 reduced matrix elements for the transitions J → J -2 between four-nucleon-hole states with T = 0 in the 1g 9/2 orbit ( 96 Cd), calculated with the SLGT0 shellmodel interaction and compared with the mapped b-IBM. Matrix elements are expressed in units 55/3π(e ν + e π )l 2 ho (see text).

  and the since the boson matrix element on the right-hand side equals √ 5e b , we find the following expression of the boson effective charge e b in terms of the shell-model neutron and proton effective charges: + e π ).

FIG. 3 :FIG. 4 :FIG. 5 :

 345 FIG.3: Absolute values of the E2 reduced matrix elements for the transitions J → J -2 between six-nucleon-hole states with T = 0 in the 1g 9/2 orbit ( 94 Ag), calculated with the SLGT0 shellmodel interaction and compared with the mapped b-IBM. Matrix elements are expressed in units 55/3π(e ν + e π )l 2 ho (see text).

ArbModel

  and Aurore Dijon for her help with the shell-model calculations. This work has been carried out in the framework of CNRS/DEF project N 19848. S.Z. thanks the Algerian Ministry of High Education and Scientific Research for financial support. This work was also partially supported by the Agence Nationale de Recherche, France, under contract nr ANR-07-BLAN-0256-03.

TABLE I :

 I Two-body matrix elements in the 1g 9/2 orbit, in units of MeV, derived from the interactions SLGT0 and GF, and from the experimental spectrum of 90 Nb. ν λ ≡ (1g 9/2 ) 2 ; λT | ĤF 2 |(1g 9/2 ) 2 ; λT

	(λ, T )	SLGT0	GF	Nb90
	(0, 1)	-2.392	-2.321	-1.758
	(1, 0)	-1.546	-1.524	-1.225
	(2, 1)	-0.906	-0.937	-0.573
	(3, 0)	-0.747	-0.700	-0.521
	(4, 1)	-0.106	-0.160	0.064
	(5, 0)	-0.423	-0.447	-0.332
	(6, 1)	0.190	0.140	0.266
	(7, 0)	-0.648	-0.640	-0.481
	(8, 1)	0.321	0.241	0.334
	(9, 0)	-1.504	-1.752	-1.376

TABLE II :

 II Fractions of 2p 1/2 + 1g 9/2 shell-model eigenstates of the GF interaction that lie within the 1g 9/2 subspace, expressed in percentages.

	J π i	0 + 1	0 + 2	2 + 1	2 + 2	4 + 1	4 + 2	6 + 1	8 + 1	10 + 1
	96 Cd	96	86	98	97	98	98	98	97	99
	92 Pd	90	95	92	95	94	93	95	94	96
	interactions, while having reasonable characteristics, are sufficiently different to test the
	robustness of our analysis.

TABLE IV :

 IV Overlaps of the (1g 9/2 ) 4 yrast eigenstates of the SLGT0 interaction with angular momentum J and isospin T = 0 with various two-pair states, expressed in percentages.

	J	B 2	SP J	D 2	DG	DI	DK	G 2	I 2	K 2
	0	91	80	35	-	-	-	18	7.4	1.9
	2	97	85	17	22	-	-	1.5	0.0	0.4
	4	89	64	42	11	11	-	0.2	0.2	0.0
	6	55	70	-	43	0.2	4.3	0.0	0.2	0.0
	8	5.3	83	-	-	7.4	24	1.8	0.2	0.1
	10	42	-	-	-	-	58	-	6.1	0.5
	12	88	-	-	-	-	-	-	57	1.5
	14	96	-	-	-	-	-	-	-	31.4
	16	100	-	-	-	-	-	-	-	100

TABLE V :

 V Energies (in MeV) of T = 0 levels for four nucleons in the 1g 9/2 orbit ( 96 Cd) calculated with the shell-model interaction SLGT0 and compared with various versions of IBM obtained by democratic or OAI mapping. E 0 is the binding energy of the ground state.

		SLGT0	b-IBM	sb-IBM	sd-IBM	sdg-IBM
	E 0	9.050	8.643	9.041	8.932	9.050
	0 + 1	0	0	0	0	0
	2 + 1	0.963	0.678	1.077	1.199	1.002
	4 + 1	2.100	1.941	2.339	3.754	2.204
	6 + 1	3.079	3.302	3.700	-	4.034
	8 + 1	3.449	4.425	4.824	-	5.688
	10 + 1	5.227	5.179	5.578	-	-
	12 + 1	5.904	5.572	5.971	-	-
	14 + 1	6.056	5.692	6.091	-	-
	16 + 1	5.904	5.496	5.895	-	-
	18 + 1	-	∞	∞	-	-
	0 + 2	4.594	-	4.613	4.491	4.594
	2 + 2	4.491	-	-	4.730	4.554
	4 + 2	4.390	-	-	-	4.538

TABLE VII :

 VII Energies (in MeV) of T = 0 levels for eight nucleons in the 1g 9/2 orbit ( 92 Pd) calculated with the shell-model interaction SLGT0 and compared with various versions of the IBM obtained with two methods of mapping, democratic (Dem) and OAI. E 0 is the binding energy of

	the ground state.
		SLGT0 b-IBM	sb-IBM	sd-IBM
				Dem	OAI	Dem	OAI
	E 0	18.937	18.135 18.771 18.646 18.624 19.999
	0 + 1	0	0	0	0	0	0
	2 + 1	0.927	0.637	1.170 0.917	0.728 0.762
	4 + 1	1.728	1.104	1.740 1.608	1.561 2.054
	6 + 1	2.512	1.965	2.628 2.441	3.155 4.267
	8 + 1	3.198	2.836	3.501 3.320	5.486 6.861
	10 + 1	4.233	3.683	4.325 4.185	-	-
	12 + 1	5.123	4.414	5.050 4.924	-	-
	D. Electric quadrupole properties
	A further test of the aligned-n-p-pair hypothesis can be obtained from electric quadrupole
	(E2) transition properties. The E2 operator in the shell model is given by

ArbModel [START_REF] Heinze | program ArbModel[END_REF] and IBM-3 [START_REF] Van Isacker | [END_REF]. The former is a general purpose program that can handle systems of fermions and/or bosons with arbitrary spins and can thus be used for the shellmodel as well as the IBM calculations; the latter code is specifically written for the isospininvariant sd-IBM. Alternatively, for three and four identical bosons (i.e., for b-IBM) the calculations can be performed with the expressions given in Sect. II and equivalent ones for the three-hole case.

A few remarks are in order. All results concern absolute energies. In the first line of each table are given the binding energies E 0 of the T = 0 ground states, as obtained in the various mappings, which should be compared with the corresponding quantity in the shell model. In subsequent lines are given the energies of a selected number of states, relative to this ground state. This might lead to some seemingly counterintuitive results. For example, it is seen from Table IV that the four-hole 2 + 1 shell-model state overlaps 97 % with a B 2 configuration. Why, then, should its excitation energy come out rather poorly in b-IBM (0.678 MeV compared with 0.963 MeV in the shell model, see Table V)? The reason is that the absolute energy of the 2 + 1 state is rather well reproduced (it misses only 0.122 MeV of the shell-model correlation energy) while the absolute energy of the ground state is rather worse underbound (by 0.407 MeV) in b-IBM.

The two-boson spectra obtained from mapping the four-hole shell-model results depend on the kind of pairs included in the basis but otherwise they are identical in the OAI and democratic mappings. Therefore, for each of the different IBM versions, there is a unique spectrum shown in Table V which is identical to that of the shell-model hamiltonian when diagonalized in the corresponding (possibly truncated) two-pair basis. While the OAI and democratic mappings yield the same energy spectrum for four holes, they lead to different boson-boson interaction matrix elements. Hence, the OAI and democratic results are different for the six-and eight-hole spectra shown in Tables VI andVII.

If the number of two-pair states equals the number of independent four-hole shell-model states, then the two-boson calculation reproduces the four-hole results exactly. According to Table V this happens, for example, for the three J π = 0 + states with T = 0 which can be exactly described as combinations of |S 2 ; 0 , |D 2 ; 0 , and |G 2 ; 0 . Consequently, the three shell-model 0 + states are exactly reproduced in sdg-IBM.

Because of the Pauli exclusion principle, no four-hole shell-model state exists with J = 18 while this angular momentum is allowed in the coupling of two bosons with J = 9. This is an


Hellen Altendorf
email: hellen.altendorf@mines-paristech.fr

Dominique Jeulin
email: dominique.jeulin@mines-paristech.fr

Random Walk Based Stochastic Modeling of 3D Fiber Systems

For the simulation of fiber systems, there exist several stochastic models: systems of straight nonoverlapping fibers, systems of overlapping bending fibers or fiber systems created by sedimentation. However, there is a lack of models providing dense, non-overlapping fiber systems with a given random orientation distribution and a controllable level of bending. We introduce a new stochastic model in this paper, that generalizes the force-biased packing approach to fibers represented as chains of balls. The starting configuration is modeled using random walks, where two parameters in the multivariate von Mises-Fisher orientation distribution control the bending. The points of the random walk are associated with a radius and the current orientation. The resulting chains of balls are interpreted as fibers. The final fiber configuration is obtained as an equilibrium between repulsion forces avoiding crossing fibers and recover forces ensuring the fiber structure. This approach provides high volume fractions up to 72.0075 %.

I. INTRODUCTION

The increasing interest in fibrous materials expands to a large variety of use cases. The most common fiber reinforced polymers are nowadays contained in the enclosure of aircrafts, boats and cars, but also wound disinfection tissues and thermal insulations make use of fibrous media. The macroscopic properties of these materials are highly influenced by the geometry of the fiber component, in particular by the direction distribution. Physical properties can be optimized by adapting the structural parameters respectively. To this end, we need a realistic stochastic model for the geometry including the main parameters for the fiber structure and taking into account the natural variability of the system.

In the area of stochastic geometry models, we distinguish between softcore and hardcore systems. A hardcore system forbids overlap between the objects, which is strongly required, as such penetrations are not realistic in the case of solid glass or carbon fibers. Concerning random fiber systems, there exist several approaches providing either hardcore fiber systems with a low volume fraction, softcore systems and sedimentation algorithms, providing hardcore systems with restrictions on the orientation distribution.

We recall in details the state-of-the-art in fiber modeling: The classical approach called dilated Poisson line process was introduced in [START_REF] Matheron | Random sets and integral geometry[END_REF] and creates infinite straight cylinders in a softcore network. Since then, the fiber modeling has evolved into more flexible approaches. In order to achieve hardcore systems, the random sequential adsorption (RSA) model was created in [START_REF] Widom | Random sequential addition of hard spheres to a volume[END_REF], which iteratively generates objects and tries to place them in the system, such that they do not overlap with the already existing objects. This approach was applied for cylinders in [START_REF] Feder | Random sequential adsorption[END_REF]. Two other random fiber packing methods have been developed: one for ellipsoids presented in [START_REF] Bezrukov | Simulation and statistical analysis of random packings of ellipsoids[END_REF] and one for spherocylinders in [START_REF] Williams | Random packings of spheres and spherocylinders simulated by mechanical contraction[END_REF]. All three random packing approaches produce straight fibers in a hardcore system and achieve only low volume fractions of about 10% -15% for isotropic orientation distribution and a fiber aspect ratio of 10. For long fibers, as it is the case for fiber reinforced composites (aspect ratio of about 200), the producible volume fraction goes down to 5%, while 15% -55% is requested. Furthermore, long fibers cannot be realized in a periodic window, as fibers tend to overlap themselves. Still, the periodic boundary condition is often required for simulations of physical properties.

There are also more physically motivated approaches of fiber sedimentation, which achieve in general high volume fractions, but are limited in realizing a given orientation distribution. The fiber deposition model in [START_REF] Provatas | Fiber deposition models in two and three spatial dimensions[END_REF] generates fibers oriented in the plane and deposes them with a certain bending parameter on the existing system. In this approach, high volume fractions can be achieved, but the most important parameter, the orientation distribution, is limited to the plane. The sequential deposition algorithm in [START_REF] Coelho | Geometrical and transport properties of random packings of spheres and aspherical particles[END_REF] realizes a fall of random particles until they reach a local minimum of their potential energy. During the fall, rotation and displacement are not restricted. Therefore, the particle can change its orientation arbitrarily and it cannot be assured that the desired orientation distributions is achieved. This approach allows up to 59% volume fraction for cylinders with aspect ratio 10.

Another hardcore bending fiber model was proposed in [START_REF] Pan | Numerical generation of a random chopped fiber composite rve and its elastic properties[END_REF]. The representative volume element is divided in sublayers. Fibers are randomly oriented in the xy-plane and in case of overlap, the newly added fiber changes the sublayer to avoid the existing fiber. The resulting chains of polyhedra with ellipsoid cross section achieve a volume fraction of 35% -40%. Still, the buckled polyhedra appear not very realistic and again the orientation distribution is restricted to the xy-plane.

Moreover, softcore bending fiber systems can be achieved by random walks introduced in [START_REF] Wirjadi | Models and Algorithms for Image-Based Analysis of Microstructures[END_REF] with the von Mises-Fisher distribution, controlling the smoothness of the bending. In this approach it is possible that long fibers bend to a circle, thus the final orientation distribution cannot be controlled. The approach [START_REF] Karkkainen | A stochastic shape model for fibres with an application to carbon nanotubes[END_REF] provided 2D random walks with the multivariate von Mises-Fisher distribution, which controls in addition the deviation from the main fiber orientation, assuring a certain loyalty to a given orientation distribution. Furthermore, in [START_REF] Faessel | 3d modelling of random cellulosic fibrous networks based on x-ray tomography and image analysis[END_REF] a bending fiber model is proposed based on dilated cores, made of curvature points connected by the spline interpolation. The last three mentioned models are softcore systems.

The idea to represent fibers as chains of spheres was already considered in [START_REF] Ch | Monte carlo scheme for generation and relaxation of dense and nearly jammed random structures of freely jointed hard-sphere chains[END_REF]. An irregular assembly of chains of tangent hard spheres (like in a pearl necklace) were studied to build dense packings with Monte Carlo simulation schemes. Hardcore configurations with low density are packed with a "box shrinkage" algorithm and relaxed in Monte Carlo algorithms by localized moves as flip, rotation, reptation and intermolecular reptation. The aim of the approach in [START_REF] Ch | Monte carlo scheme for generation and relaxation of dense and nearly jammed random structures of freely jointed hard-sphere chains[END_REF] is a highly dense packing, which was successfully achieved with a volume fraction of 63.9%, which corresponds to the densest packing of hard spheres. To this end, the restrictions on the fiber structure are very low, and thus, the approach is not suitable for reconstructing fiber reinforced composites.

Material properties are highly influenced by the fiber structure, mainly its orientation distribution. This correlation is studied in [START_REF] Schulgasser | Fibre orientation in machine-made paper[END_REF], [START_REF] Chin | Effects of fiber length and orientation distribution on the elastic modulus of short fiber reinforced thermoplastics[END_REF], [START_REF] Jain | Effect of fiber orientation on the fracture toughness of brittle matrix composites[END_REF], [START_REF] Hine | Measuring the fibre orientation and modelling the elastic properties of injection-moulded long-glass-fibre-reinforced nylon[END_REF], [START_REF] Shao | Effects of fiber length and fiber orientation distributions on the tensile strength of shortfiber-reinforced polymers[END_REF] and [START_REF] Favier | Simulation and modeling of three-dimensional percolating structures: Case of a latex matrix reinforced by a network of cellulose fibers[END_REF]. Furthermore, Berhan ([3] and [START_REF] Berhan | Modeling percolation in high-aspect-ratio fiber systems. ii. the effect of waviness on the percolation onset[END_REF]) studied the influence of model characteristics on percolation simulation for different high-aspect-fiber systems. In [START_REF] Berhan | Modeling percolation in high-aspect-ratio fiber systems. i. soft-core versus hard-core models[END_REF], the influence of softcore and hardcore systems is examined, while [START_REF] Berhan | Modeling percolation in high-aspect-ratio fiber systems. ii. the effect of waviness on the percolation onset[END_REF] cares about the waviness of the fibers. Berhan states that the frequently used straight and overlapping fiber systems are inappropriate to model fiber reinforced materials and he suggests to make use of hardcore bending fiber systems. Those studies enhance the need for stochastic models, creating realistic fiber systems dedicated to fiber composite materials. In [START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF] and [START_REF] Altendorf | Fiber separation from local orientation and probability maps[END_REF], we proposed quantification methods to analyze the most important characteristics of fiber systems. Now, a stochastic model is requested to realize the measured properties. In contrast to physical simulations dedicated to systems of liquid crystals (as for example [START_REF] Cinacchi | Phase behavior of wormlike rods[END_REF]), we are only interested in the final configuration of the stochastic model and its properties, instead of the behavior of the fibers during the simulating process.

The aim of this paper is to provide a random hardcore fiber model with a controllable bending and high volume fractions. For this purpose, random walks are used to create a realistic system of bending fibers. The level of bending is controlled by two parameters in the multivariate von Mises-Fisher distribution. The created points are provided with a radius and the current orientation. The balls, defined by the points and the assigned radii, are connected in chains according to the paths of the random walks. This approach results in a softcore system of bending fibers, modeled as chains of balls. In a second step, we apply a force-biased approach, to achieve a hardcore configuration of the fiber system. Force-biased algorithms on spheres were introduced in [START_REF] Mościński | The force-biased algorithm for the irregular close packing of equal hard spheres[END_REF] and statistically analyzed in [START_REF] Bezrukov | Statistical analysis of simulated random packings of spheres[END_REF]. The forces in our approach were inspired by the energy reducing models known from molecular dynamics [START_REF] Mościński | The force-biased algorithm for the irregular close packing of equal hard spheres[END_REF]. The mentioned forces describe the necessary displacement of the balls to relax the system. They do not act like mechanical forces. Two kinds of forces are applied to the ball centers: repulsion and recover forces. The repulsion force arises in case of a fiber overlap and displaces the balls to a non-overlapping position. The recover force maintains the fiber structure between the balls. It keeps the distance and the angles between a ball and its neighbors, allowing only small deviation.

The approach works in a closed cubic window W ⊂ R 3 with periodic boundary conditions. Overlap checks are effectively computed in subdivisions of the window W , where only neighboring cells need to be compared. Therefore, the algorithm has low computing time. We assured the quality of realization by evaluating the parameters of the final configuration.

Although ball chains are considered as the second type of ideal amorphous solids, they are purely used in modeling of fibrous material. As stated in [START_REF] Laso | Random packing of model polymers: local structure, topological hindrance and universal scaling[END_REF], in the case of non-overlapping neighbor spheres, the maximal possible volume fraction rises up to the one of individual spheres: ≈ 0.64. With decreasing distance between neighbor spheres and increasing alignment of fibers, we approach to a system of parallel cylinders. In this case, the problem reduces to a random packing of 2D discs, for which the maximal volume fraction can reach up to ≈ 0.78 [START_REF] Einar | Random packing of disks in two dimensions[END_REF]. In our approach, the distance between neighbor spheres is ≈ r/2, where r is the sphere radius. This distance is one quarter of the necessary distance for non-overlapping neighbor spheres. Thus, the maximal possible volume fraction may lay between 0.64 and 0.78. As we include higher restrictions to the fiber structure as a random packing of tangent sphere chains [START_REF] Ch | Monte carlo scheme for generation and relaxation of dense and nearly jammed random structures of freely jointed hard-sphere chains[END_REF], it is probable that we cannot reach such high volume fractions for every set of parameters. Still, volume fractions around 50% were achieved for several input parameters. The maximal volume fraction of our experiments was 72.0075 %, achieved for a z-axis preferred orientation distribution and an aspect ratio of 9.

II. METHOD

In this section, we describe the two main steps of the approach: the random walk and the force-biased fiber packing. We explain how the structural parameters of the fiber system are incorporated into the random walk and how we control the bending of the fibers. For the force-biased fiber packing, we introduce the repulsion and recover forces and explain, how they are applied to the system. Furthermore, the stop criterion and the endstep problem is discussed, and implementation details are given.

A. Random Walk

Each fiber of the initial fiber system is modeled as a random walk. A random walk is a Markov process on R 3 producing a chain of points. We assign to every point the current direction and a radius. The result of a random walk is a sequence of points in R 3 × S 2 × R + , describing a fiber:

P = {p 0 , . . . , p l } with p i = (x i , µ i , r i ) ∈ R 3 × S 2 × R + .
The starting point p 0 = (x 0 , µ 0 , r 0 ) and the path length l could be generated from four main distributions for the starting coordinate, the orientation, the radius and the length. In general, the initial coordinate x 0 is uniformly distributed in a cubic window with periodic edge treatment x 0 ∼ U (W ⊂ R 3 ). The other three distributions describe the main parameters of the fiber system. For the orientation distribution, we have applied the β orientation (see [START_REF] Schladitz | Design of acoustic trim based on geometric modeling and flow simulation for non-woven[END_REF] or [START_REF] Ohser | 3d Images of Materials Structures -Processing and Analysis[END_REF]) with a global parameter β ∈ R + \{0}. For β = 1 it results in the uniform distribution on the sphere; for β → 0 the distribution concentrates on the z-axis and for β → ∞ the orientations are distributed isotropically in the xy-plane. The probability density function of the β orientation distribution is

p(θ, φ|β) = β sin θ 4π(1 + (β 2 -1) cos 2 θ) 3/2 ,
where (θ, φ) are the polar coordinates of the orientation µ 0 ∈ S 2 .

The radius and length can be simulated from any distribution on R + . The radius can either be fixed for the whole system, chosen once for each fiber or for every point during the random walk with some strategy for smooth variation. The amount of fibers for the system is chosen dependent on the required volume fraction. The fourth main characteristic of the fiber system, the bending of the fibers is indirectly controlled by two parameters κ 1 and κ 2 in the multivariate von Mises-Fisher distribution.

During the random walk, the generation of a new point p i+1 requires three steps: First, the new orientation µ i+1 is generated according to the multivariate von Mises-Fisher distribution, with the last orientation µ i and the main fiber orientation µ 0 as parameters. Second, the new radius r i+1 is generated. Finally, the new coordinate is calculated as

x i+1 = x i + ri+1
2 µ i+1 . The distance between the points is a trade-off between a representative fiber structure and a treatable amount of points.

In [START_REF] Karkkainen | A stochastic shape model for fibres with an application to carbon nanotubes[END_REF] the multivariate von Mises-Fisher distribution is introduced for the 2D case. In 3D the corresponding probability density function is defined as:

f (x|x 1 , κ 1 , x 2 , κ 2 ) = c(x 1 , κ 1 , x 2 , κ 2 ) e κ1x T 1 x+κ2x T 2 x ,
where x 1 , x 2 ∈ S d-1 are two preferred directions and κ 1 , κ 2 ∈ R + the reliability parameters towards the preferred directions. The factor c(x 1 , κ 1 , x 2 , κ 2 ) serves for the normalization, such that the integral over S 2 is equal to 1.

A useful observation is that every multivariate von Mises-Fisher distribution can be written as a classical von Mises-Fisher distribution with the parameters

κ = |κ 1 x 1 + κ 2 x 2 | and x 0 = κ 1 x 1 + κ 2 x 2 κ ,
which can be easily checked by inserting those parameters in the classical von Mises-Fisher density function. This simplifies the generation of the pseudo random variables to the standard case, which is well described in [START_REF] Fisher | Statistical analysis of spherical data[END_REF]. Furthermore, this defines the normalization factor

c(κ) = κ 2π(e κ -e -κ ) ⇒ c(x 1 , κ 1 , x 2 , κ 2 ) = |κ 1 x 1 + κ 2 x 2 | 2π(e |κ1x1+κ2x2| -e -|κ1x1+κ2x2|
) .

In our approach, the first preferred direction is the main fiber orientation µ 0 and the second preferred direction is the orientation of the previous point µ i-1 . Therefore, κ 1 describes the reliability to the main fiber orientation, where κ 2 describes the reliability to the last orientation and hence specifies the smoothness of the bending. The probability density function for the i th step in the random walk (with i ≥ 1) is

f (µ i |µ 0 , κ 1 , µ i-1 , κ 2 ) = c(µ 0 , κ 1 , µ i-1 , κ 2 ) e κ1µ T 0 µi+κ2µ T i-1 µi .
Realizations of the softcore fiber model created by the random walks with varying choices of the reliability parameters κ 1 and κ 2 are shown in FIG. 1. The mean fiber orientation of a created ball chain P = {p 0 , . . . , p l } is defined as μ(P ) = p l -p0 |p l -p0| ∈ S 2 . The mean orientations of the created fibers are supposed to be realizations of a global orientation distribution for the fiber system. The starting orientations µ 0 are chosen with respect to the global orientation distribution, but the mean orientation differs from µ 0 with a deviation according to the reliability parameters κ 1 and κ 2 . To ensure the global orientation distribution, we rotate the chosen fiber, so that the mean fiber orientation equals the originally chosen main fiber orientation: μ(P ) = µ 0 . The fiber structure (length, radii, bending) remains constant during the rotation. Every ball

p i = (x i , µ i , r i ) is adjusted to p i = (x ′ i , µ ′ i , r i ) with x ′ i = x 0 + Rot(x i -x 0 , n, α) µ ′ i = Rot(µ i , n, α) n = µ 0 × μ(P ) α = ∡(µ 0 , μ(P )) Rot(µ, n, α) = (n • µ)n + cos α((n × µ) × n) + sin α(n × µ),
where "•" notates the scalar product (or inner vector product) and "×" the cross product (or vector product). The random walks create a stochastic fiber model made of overlapping ball chains. The system is defined by the graph (P, E), where

P = {p 1,1 , p 1,2 , . . . , p 1,l1 , p 2,1 , . . . , p n,ln } ⊂ R 3 × S 2 × R +
yields the created points and E = {(p i,j , p i,j+1 ) | i ∈ {1, . . . , n}, j ∈ {1, . . . , l i }} the connections between the points according to the n paths of the random walks. l i is the amount of points for the i th fiber. Furthermore, we define

C = {(q, p, q ′ ) ∈ P 3 | (q, p) ∈ E ∧ (p, q ′ ) ∈ E}
as the joints between two edges. For further steps, we need to remember the original distance between two neighboring points and the original angle in every joint. Thus, we define an original distance and an original angle function

D orig : E → R + (p, q) → |x p -x q | A orig : C → [0, π] (q, p, q ′ ) → ∡(x q , x p , x ′
q ) and fix the mappings before the fiber packing process starts.

B. Force-Biased Fiber Packing

To achieve a hardcore fiber system, the fiber packing applies two kinds of forces to the ball centers given from the softcore fiber model: the repulsion force to separate overlapping balls and the recover force to keep the fiber structure between the balls. The reconfiguration of the ball chains are realized in several steps. In each step, the forces are calculated for the current configuration of the system and displacements are performed respectively.

Repulsion Force

The repulsion force for a pair of overlapping balls describes the necessary displacement to make them nonpenetrating. The overlap between two balls is permitted, if the balls are close neighbors in the same fiber. In this case, we define the two points p, q ∈ P as related: p ∼ q ⇔ ∃ path between p and q with length ≤ 5.

(

The value of 5 is an approximation for the minimal distance of two balls in the same chain, that do not naturally overlap. The distance of the ball centers of two direct neighbor balls is chosen r/2. In this case, balls in the same chain connected with a path with length less than 4, do overlap naturally. According to small deviations in the radius or additional curvature, we increase the minimal path length to 5.

Let

p 1 = (x 1 , d 1 , r 1 ), p 2 = (x 2 , d 2 , r
2 ) ∈ P be two arbitrary points in the fiber system. The overlap between those two points is

Overlap(p 1 , p 2 ) = max(0, r 1 + r 2 -d(x 1 , x 2 ))
with the distance function d(x 1 , x 2 ) respecting the periodic edge treatment. The repulsion force is

F rp (p 1 , p 2 ) = 1 p1≁p2 Overlap(p 1 , p 2 ) 2 x 1 -x 2 |x 1 -x 2 | ,
and is applied, if the points are not related as defined in equation [START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF]. The force on p 1 works in the opposite direction to p 2 with a strength linearly dependent on the overlap. The total repulsion force for the point p 1 is cumulated over all points in P :

F rp (p 1 ) = q∈P F rp (p 1 , q).

Recover Force

The recover force keeps the distance between neighbor points by modeling spring-like forces between them, and it keeps the angle in the joints by modeling open springs between the connections. Open springs describe the fact, that there is a force for decreasing, but none for increasing angles. Thus, a fiber can be straightened, but not folded into a clew.

The recover force is induced by a displacement of the original ball coordinate as effect of the applied repulsion force. By applying the recover forces to one ball it induces recover forces on its neighbors. Even if the strength of the indicated forces is decreasing, a kind of domino effect has started and the process becomes infinite. In order to stabilize the movement, we introduce a smoothing factor, which acts like friction. We define a minimal change x s , which is necessary to induce a force and a change x e , when the force regains its full strength. The friction factor is defined as

f xs,xe (x) =      0 , x < x s 1 2 -1 2 cos( |x|-xs xe-xs π) , x s ≤ x ≤ x e 1 , x e < x .
The curve f xs,xe (x) is given in FIG. 2.

FIG. 2. Friction factor for recover forces.

a. Spring Force: In order to keep the distances between the balls in a fiber, we define a force between the centers of the balls, with a strength linearly dependent on the change of distance. This kind of force can be interpreted as a spring between the ball centers. Let p = (x p , d p , r p ) and q = (x q , d q , r q ) be two neighboring points. The change of distance is ∆D(q, p) = D orig (q, p) -|x q -x p | and the unit vector for the direction of the force is

v(p, q) = x p -x q |x p -x q | .
Thus, the recovery spring force F rs on a point p is defined as

F rs (p) = q∈P (1 E (p, q) + 1 E (q, p))f ds,de |∆D(q, p)| D orig (q, p) (∆D(q, p))v(p, q) .
The spring force is applied to p for every direct neighbor (at most two). It is linear to the distance change and multiplied by the friction factor on the ratio of distance change. For the realizations created in this paper, the friction parameters were chosen as d s = 5 % and d e = 10 %. b. Angle Force: In order to keep the angle between joints, we define a recover force on each point having two direct neighbors. Let (q 1 , p, q 2 ) ∈ C be a joint in the point p = (x, d, r) with the neighbors q 1 = (x 1 , d 1 , r 1 ) and q 2 = (x 2 , d 2 , r 2 ). We define:

α 0 = A orig (q 1 , p, q 2 ) m = x 1 x 2 ∩ Plane (x, ⊥(x 1 -x 2 )) h 1 = |m -x 1 |, h 2 = |m -x 2 |, z = |m -p| α 1 = tan -1 (h 1 /z), α 2 = tan -1 (h 2 /z) tan α = tan(α 1 + α 2 ) = z(h 1 + h 2 ) z 2 -h 1 h 2 z 0 = h 1 + h 2 + (h 1 + h 2 ) 2 + 4h 1 h 2 tan 2 α 2 tan α .
Now, the recovery angle force F ry on a point p is defined as

F ra (p) = 1 ∃q1,q2∈P,(q1,p,q2)∈C f αs,αe (α 0 -α) (z -z 0 )v(m, p).
The angle force is linearly dependent on the necessary displacement to regain the original angle and it is multiplied by the friction factor, dependent on the angle change. For the realizations created in this paper, the friction parameters were chosen as α s = 0.1 • and α e = 0.2 • . The geometrical construction is visualized in FIG. 3.

Application of Force

The total force on an arbitrary point p ∈ P is composed of the sum of all mentioned forces:

F total (p) = q∈P F rp (p, q) + ρF rs (p) + ρF ra (p),
where ρ is a factor in [0, 1] leading to a smooth development of the recover forces, favoring the repulsion forces in the first steps. For the realizations created in this paper, we chose ρ = 0.2. In order to limit the possible displacement of a point, we define an upper bound F max for the final force strength:

F f inal (p) = F total (p) |F total (p)| min(|F total (p)|, F max ).
In each step, the ball center of p is displaced by the limited final force:

x ′ p = x p + F f inal (p).
After the displacement in a step, the forces are recalculated for the new configuration and applied in a following step.

C. Stop Criterion and End Step

During the packing process the forces decrease very fast in the beginning and converge slowly to zero in the end. The process should be stopped by a criterion dependent on the force strength. In our implementation, the algorithm terminates with a solution, if the total force strength p∈P F total (p) falls below a certain limit and the displacements become negligible. For the realizations created in this paper, we chose the limit as 0.002 s n, where n is the amount of fibers and s the side length of the window. That means the mean sum of necessary displacements in a fiber is smaller than 0.2% of a unit size. When the volume fraction is chosen too high, it is possible that there is no configuration with low force strength. In this case, the process will not converge and we need to add a stop criterion. We have chosen an upper bound for the amount of steps of 10 000. If the process reaches this limit, it will terminate unsuccessfully. In the case of a successfully finished process, the total force strength is negligible. Still, one should be aware that the total force on a point is a sum of different forces. Thus, even if the sum of forces is zero, it is not sure that every specific force is negligible in its strength. Regarding the repulsion force as specific force, it may happen that the recover force acts in the opposite direction with exactly the same strength (as shown in Fig. 4(a). As the repulsion force is linear to the overlap, we conclude that the fiber packing has not necessarily a zero overlap, that means the system may not be totally hardcore.

This effect is similar to the maximal random jammed state, known from random packings (see [START_REF] Laso | Random packing of model polymers: local structure, topological hindrance and universal scaling[END_REF]). In our approach, a ball is never jammed as its displacement is defined by reducing forces and overlap in a following configuration is not strictly forbidden. Still, we have a similar situation if the requested volume fraction is too high for the fiber structure.

To assure a non-overlapping system, we accept an outcome only if the maximal overlap falls below a certain limit. Here, we chose 0.1 r min , where r min is the minimal radius in the system. Then, a final step is added, where the radius of each point is reduced according to the maximal overlap with other points. This approach is similar to the idea in the force biased sphere packing. Let p = (x, d, r) ∈ P be an arbitrary point of the fiber system, the radius r is reduced to

r ′ = r - 1 2 max q∈P,q≁p
Overlap(p, q).

If the accuracy of the radius distribution is very important, the end step can be replaced by a second round of the packing algorithm, where only the repulsion force is active. In this case, the theoretical recover force strength should be checked afterwards, to make sure that the fibers are not drawn apart.

D. Implementation Details

The periodic distance function for x = (x x , x y , x z ) T and y = (y x , y y , y z ) T and a window W = [0,

w x ]×[0, w y ]×[0, w z ] is defined as d(x, y) =   min(|x x -y x |, w x -|x x -y x |) min(|x y -y y |, w y -|x y -y y |) min(|x z -y z |, w z -|x z -y z |)   .
The repulsion force on an arbitrary point p ∈ P is defined as sum over all points q ∈ P . Still the repulsion force is equal to zero, if the balls p and q are not overlapping, which is definitely the case if q / ∈ B rp+rmax (x p ), where r max = max q∈P r q . As proposed in [START_REF] Mościński | The force-biased algorithm for the irregular close packing of equal hard spheres[END_REF] we divide the window W in smaller cubic subwindows with a sidelength s i ≥ r max with i ∈ {x, y, z}. For a window W = [0, w x ] × [0, w y ] × [0, w z ] the sidelengths of the subwindow are calculated as

s i = w i n i with n i = ⌈ w i r max ⌉.
Thus, B rp+rmax (x p ) is included in the union of the subwindow including p and its neighbor subwindows. The overlap check is limited to this union.

III. RESULTS

FIG. 5 shows realizations of the presented model with varying input parameters for the fiber aspect ratio χ, the amount of objects n and the main orientation distribution. The parameters and the achieved volume fraction V V is given in the corresponding caption.

The aspect ratio is generally defined as length divided by diameter. In the case of a ball chain p 1 , . . . , p n with fixed radius r the aspect ratio χ can be calculated as

χ = 2r + n-1 i=1 |p i -p i+1 | 2r .
To evaluate, if the final configuration fulfills the required distributions, we observed the changes in the characteristics of the fiber system during the force biased packing for 3 realizations with different choices of β ∈ {0.1, 1, 10}. The remaining parameters are chosen fix: window size = 100 3 , length = 50, radius = 5, amount of fibers = 100, requested density = 44.483%. The evolution of the characteristics is shown in Fig. 6. The estimation of length and radius is obvious. The measured density is the volume fraction of the discretized image. The beta estimate is the numerical approximation of the maximal likelihood estimator, which cannot be resolved theoretically. The estimation of the bending parameters κ 1 and κ 2 is not that simple, but experiments showed that they can be approximated for κ 2 /κ 1 > 2 with the following estimates:

κ1 = 1/Var 1 κ2 = 2/Var 2 Var 1 = 1 n n j=1   1 l j lj i=0 µ j,i -μ(p j,1 , . . . , p j,lj ) 2   Var 2 = 1 n n j=1   1 l j -1 lj -1 i=0 (µ j,i -µ j,i+1 ) 2   .
Even if the estimators for κ 1 and κ 2 are not highly accurate, they are still measures for the global and local bending.

As long as the measures of the real dataset and the measures of the final model configuration are close, the aim of realistically reproducing the fiber structure is fulfilled.

In FIG. 6, the first steps are mostly driven by the repulsion forces, which evokes changes in the local structure and therefore in the parameter measure. When the repulsion force decreases, the recover forces reconstruct the fiber structure, which results in a convergence of the measures to the requested values. Most measures start with the requested values, except for the density and the global bending parameter κ 1 . The low density in the beginning is evident, as the overlapping fiber volume is counted only once. The deviation of the global bending parameter κ 1 or reliability to the mean fiber orientation is caused by the different measure condition. During the realization of the fiber structure, the κ 1 represents the reliability to the chosen fiber orientation µ 0 , whereas the measure of κ 1 reflects the reliability to the actual mean fiber orientation μ(P ), which is slightly too high. The convergence to the requested value of κ 1 is surprisingly good, and probably caused by the trade-off between the recover force to the original fiber structure and the tolerance range of local deviation.

We stress once more that the important configuration is the final one, the geometric characteristics during the process are studied to understand better, what happens during the packing process. The evaluation of the quality of our models is only based on the geometric characteristics of the final configuration, as shown in TABLE I.

In both visualizations, we can see that the parameters of length, radius and density fit perfectly to the requested values. The global and local bending κ 1 and κ 2 are approaching to the requested values, as expected. As the local orientations have some tolerance during the packing process, it is understandable that the final bending values are lower than the initial ones. If one is not satisfied with the final values for bending, there are two possibilities to improve this: firstly one can choose the initial values higher than requested and secondly reduce the tolerance α s and α e in the packing process. The parameter β for the orientation distribution is acceptably approximated for the choices β = 0.1 and β = 1, whereas the girdle orientation with β = 10 could not be maintained that strictly during the packing process. The results would even become worse for higher choices of β. With the described model, we can realize parameters in the range of 0.05 to 5. If more restricted orientation distribution is requested, recover forces for the orientation of a fiber have to be included. This could be obtained by forces on the extremities of a fiber. As we included no restrictions in changing the mean fiber orientation during the packing process, the system will relax in a less restrictive configuration.

IV. DISCUSSION

In this section, we discuss the volume fractions achieved by numerical simulations and the evolution of the repulsion and recover forces during the fiber packing. Furthermore, we analyze the computing time.

A. Volume Fraction

The maximal volume fraction depends on the aspect ratio and on the orientation distribution of the system. FIG. 7 shows the maximal volume fractions, achieved in realizations for the β orientation distribution with β ∈ {0.1, 1, 10} and varying aspect ratio χ ∈ [START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF]33]. We achieved the highest volume fraction 72.0075 % for an aspect ratio of 9 and a z-axis preferred orientation distribution with β = 0.1. In general, the figure supports the intuitive expectation, that fibers can be packed in a particularly dense system, if they are oriented parallel and that short fibers can be packed more easily.

In FIG. 8, the ratio of achieved volume fraction and initial volume fraction is shown. The initial volume fraction is the sum of all fiber volumes after the random walk divided by the window size. The change of the distances between the balls is related to a change of the volume fraction. As we allow small deviations in the distances, we allow also small deviation in the volume fraction. Here, the deviation of the distances is limited to 10 %. As a consequence, there exists also a limit for the deviation of the volume fraction caused by the distance change. It ranges from 90 % to 110 % of the initial volume fraction. Another deviation in the volume fraction is caused by the end step, reducing the radii. The maximal overlap is 0.1 r min , and a radius r is maximally reduced by 1 2 overlap ≤ 0.05 r min ≤ 0.05 r. Therefore, the volume V of a fiber could maximally be reduced to

V ′ ≈ lπ(r ′ ) 2 > (0.95) 2 lπr 2 ≈ 0.9025 V.
Thus, a lower limit for the volume fraction would be 0.9 0.9025 = 81.225% of the initial volume fraction. The upper limit remains at 110%, since the radius is not enlarged in the end step. Note, that this is a rough approximation of the limits. Especially for dense systems, it may happen, that bending or force-canceling effects influence the structural arrangement, so that these limits are exceeded. This happened in FIG. 8 for the point with the highest requested volume fraction of 90%, which is impossible to reach for a fiber system.

The reduction of the volume fraction is expected for dense systems of fibers with a low aspect ratio, where the packing consists of moving and shortening of the objects, which causes a decrease of the volume fraction. The aspect ratio 1 is a special case, as the "fibers" consist only of one ball and shortening is not possible. For a high aspect ratio, the packing is essentially bending the fibers around each other, which causes elongation of the fibers and rising of the volume fraction.

B. Evolution of Forces

FIG. 9 shows the evolution of the repulsion and recover forces during the fiber packing. Initially, the recover force is zero as the fibers still have the initial structure, whereas the repulsion force is high at the beginning, because the fibers can overlap without restriction. The repulsion force is decreasing very fast, while the recover force rises slightly in the first steps, caused by the displacements according to the repulsion forces, and decreases quite smoothly with respect to the parameter ρ.

C. Computing Time

In practical applications the algorithm showed good performance. All time measurements are taken with a desktop PC with an Intel(R) Core(TM)2 processor with CPU X6800, 2.93 GHz and 3.8 GB memory. The creation time for the random walk progresses approximately linearly in the amount of steps, which corresponds to the total amount of ball centers in the system. The mean creation time per ball center is approximately 0.05 ms: t RW ≈ 0.05 ms #(ball center).

The computing time for the fiber packing increases exponentially with the volume fraction V V for a fixed aspect ratio χ. This is caused by the enormous rise of interaction by increasing volume fraction. Still, the realization with the longest computing time took about 41 min for an aspect ratio χ = 25, an isotropic orientation distribution and an achieved volume fraction of V V = 51.43%, which is at the limit of the possible volume fraction for the given parameters. The computing time is acceptable for those circumstances and a standard desktop computer. The computing time for the random walk and the fiber packing are shown in FIG. 10.

V. CONCLUSION AND PERSPECTIVES

We have presented an algorithm generating bending hardcore fibers, with given orientation, radius and length distributions. We showed how to evaluate the quality of the final configuration of the fiber system. For the performed realization we could achieve high convergence to the requested parameters and for the global orientation distribution, we studied the realizable range for the parameter β. In perspective, we will include further recover forces to be able to realize more restricted orientation distributions.

In practical tests, we achieved the highest volume fraction of 72.0075 % for the β orientation distribution (β = 0.1) and an aspect ratio of 9. Practical applications for isotropic orientation distribution have shown a computing time linearly dependent on the amount of points for the random walk and a computing time exponentially increasing with the initial volume fraction for the fiber packing.

A softshell ratio can be easily included, to allow partial overlap or to keep distances between the objects. In the overlap calculation, each radius would be multiplied by the softshell ratio. Further studies will include the observation of additional morphological properties during the packing process, and the application and adaption of the model parameters to data sets on real materials.

  FIG. 1. (Color online) Realizations of the fiber model created by random walks. Common parameters: amount of fibers n = 10, window side length s = 100, fix fiber length l = 300 and fix fiber radius r = 3, isotropic orientation distribution (β = 1).

FIG. 3 .

 3 FIG.3. Geometry for angle force.

  (a) Trade-off between repulsion and recover forces. (b) Solution 1: Reducing the radius. (c) Solution 2: Endstep regarding only the repulsion forces.

FIG. 4 .

 4 FIG. 4. (Color online) Trade-off problem in final configuration with two possible solutions: reducing the radius or applying an endstep regarding only the repulsion forces. Note: Force vectors are elongated for visualization reasons, they do not correspond to the necessary displacements.

  FIG. 5. (Color online) Realizations for packed fiber systems. Common parameters: window side length s = 100, bending parameters κ1 = 10 and κ2 = 100.

10 (distribution with β = 10 FIG. 6 .

 10106 FIG. 6. (Color online) Evolution of estimated model parameters during packing process, for different choices of β for the orientation distribution. Fixed parameters: window size = 100 3 , length = 50, radius = 5, amount of fibers = 100, starting density = 44.483%.

FIG. 7 .

 7 FIG. 7. (Color online) Maximal volume fractions for realizations with varying orientation distribution and aspect ratio.

1 FIG. 8 .

 18 FIG. 8. (Color online) Ratio of the achieved volume fraction and the initial volume fraction for isotropic orientation distribution.

FIG. 9 .

 9 FIG. 9. (Color online) Evolution of force during fiber packing process. Aspect ratio χ = 6. Isotropic orientation distribution.

  )= 0.04955 ms * p (a) Computing time of random walk vs. amount of ball centers. Computing time of fiber packing vs. initial volume fraction.

FIG. 10 .

 10 FIG. 10. (Color online) Computing time of the random walk and the fiber packing. For isotropic orientation distribution and constant fiber volume of 0.005 s 3 , with s = 100 window side length.

TABLE I .

 I Parameter measures for final configuration of experiments from Fig.6.

		Density Length Radius κ1	κ2	β
	Requested 44.48%	50	5	10 100 various
	β = 0.1	43.94% 50.66 4.99 9.91 96.32 0.134
	β = 1	43.90% 50.75 4.98 8.21 80.95 0.87
	β = 10	44.10% 50.90 4.98 8.03 82.39 5.6

ACKNOWLEDGMENTS

We thank Prof. Joachim Ohser, who came up with the idea of applying the force-biased approach to fibers. We also thank Dr. Oliver Wirjadi for the implementation of the VTK Surface Label Rendering, that provided us with the shown visualizations. This work is financially supported by the Institute Carnot M.I.N.E.S. and Fraunhofer ITWM, Kaiserslautern. The authors are grateful to the reviewers for helping to improve the first version of the paper.