Arnaud Derode

Éric Larose

Mickael Tanter

Julien De Rosny

Arnaud Tourin

Michel Campillo

Mathias Fink

Eric Larose

Recovering the Green's function from field-field correlations in an open scattering medium

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The possibility of recovering the Green's function from the field-field correlations of coda waves in an open multiple scattering medium is investigated. The argument is based on fundamental symmetries of reciprocity, time-reversal invariance, and the Helmholtz-Kirchhoff theorem. A criterion is defined, indicating how sources should be placed inside an open medium in order to recover the Green's function between two passive receivers. The case of noise sources is also discussed. Numerical experiments of ultrasonic wave propagation in a multiple scattering medium are presented to support the argument.

Wave propagation in a multiple scattering or reverberating environment has been a subject of interest in a wide variety of domains ranging from solid state physics to optics or acoustics. Ultrasound is particularly interesting because it allows a direct measurement of the field fluctuations, both in amplitude and in phase. In connection with this, a remarkable work by Weaver and Lobkis [START_REF] Weaver | On the emergence of the Green's function in the correlations of a diffuse field: pulse echo using thermal phonons[END_REF][START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF][START_REF] Weaver | Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies[END_REF] recently showed that the Green's function between two points could be recovered from the field-field correlation of a diffuse ultrasonic field. This amounts to doing ''ultrasonics without a source'' since they showed that thermal noise could be used instead of a direct pulse/echo measurement between the two points. The experiment was carried out in an aluminum block, and the theoretical analysis was based on discrete modal expansion of the field, with random modal amplitudes. Applications are promising: it would be possible to recover the Green's function of a complex medium just by correlating diffuse fields received on passive sensors ͑application to shallow water ocean acoustics, where the field is not diffuse but propagates in a wave guide, was also evoked 4 ͒.

However, the basic assumption in the theoretical analysis is that the medium is closed and free of absorption. In a real medium, absorption will tend to cut out the longest scattering ͑or reverberating͒ paths, and discrete modes will not be resolved any more. Similar problems are expected if the medium is open rather than closed ͑actually, in an open medium, the fluctuation-dissipation theorem [START_REF] Weaver | Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies[END_REF] establishes the result, as long as the field is diffuse in the thermal sense͒. The aim of this letter is to examine whether the Green's function can still be recovered from the correlations of an ultrasonic wave field in an open scattering medium, when a discrete expansion on orthogonal modes is no longer relevant and the field is not thermally diffuse.

To that end, we present 2-D numerical experiments of acoustic scattering on rigid inclusions randomly located either in a closed cavity or in a open medium. The wave equation is solved by a finite differences simulation ͑centered scheme͒; the boundary conditions is implemented following Collino's work. [START_REF] Collino | Absorbing boundary conditions for the finite element solution of 3D Maxwell's Equation[END_REF] Naturally, a finite difference scheme shows numerical dispersion. However, the essential point is that the fundamental symmetries of reciprocity and time reversal still hold in the numerical experiments.

To begin with, let us consider two receiving points A and B and a source C placed among a random collection of scatterers, as represented in Fig. 1. The scatterers are in water; only lossless acoustic waves are considered here. At the edges of the grid, the boundary conditions may be either perfectly reflecting ͑Dirichlet͒ as in a closed cavity or absorbing ͑open medium͒. The signal transmitted by C is a pulse with a center frequency 1 MHz and a Gaussian envelope (ϭ0.7 s).

We will note h IJ (t) as the impulse response between I and J, i.e., the wave field sensed in I when a Dirac ␦(t)i s sent by J.I fe(t) is the excitation function in C, then the wave field A and B received in A and B will be respec- tively e(t) h AC (t) and e(t) h BC (t), representing convolution. The cross-correlation of the fields received in A and B is then a͒ Electronic mail: arnaud.derode@espci.fr

1

C AB ͑ t ͒ϭ ͵ A ͑ ͒ B ͑ tϩ ͒ dϭh AC ͑ t ͒ h BC ͑ Ϫt ͒ f ͑ t ͒
with f (t)ϭe(t) e(Ϫt). f (t) depends only on the excitation imposed at the source, whereas the information regarding the impulse response between A and B is hidden in h AC (t) h BC (Ϫt). Indeed, the impulse responses of a closed cavity satisfy a remarkable property, as shown by Carsten Draeger in 1999, [START_REF] Draeger | One-channel time-reversal in chaotic cavities: theoretical limits[END_REF] which he termed the ''cavity equation:''

h AC ͑ t ͒ h BC ͑ Ϫt ͒ϭh AB ͑ t ͒ h CC ͑ Ϫt ͒. ͑1͒
For this relation to hold, the cavity must be lossless and its eigenmodes not degenerate. Note that, in practice, the correlations cannot be performed over an infinite time interval ͑the ring time of a cavity is infinite if it is lossless͒; therefore the cavity equation can be compared to experimental results if the integration time ⌬T is sufficiently large compared to 1/⌬, with ⌬ the characteristic distance between modes, so that the modes are resolved. 1/⌬ is sometimes referred to as the Heisenberg time, or break time. Figure 2 illustrates the validity of the cavity equation; here the impulse responses have been recorded during an integration time of 80 ms (2 ϫ10 6 time steps͒, and the Heisenberg time is ϳ5 ms. From Draeger's cavity equation, the correlation between the fields received in A and B is

C AB ͑ t ͒ϭh AB ͑ t ͒ h CC ͑ Ϫt ͒ f ͑ t ͒.
Therefore, similarly to Weaver's results, [START_REF] Weaver | On the emergence of the Green's function in the correlations of a diffuse field: pulse echo using thermal phonons[END_REF][START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF][START_REF] Weaver | Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies[END_REF] the direct Green's function h AB is present in the correlations of the field within a closed cavity and can be recovered from C AB provided that the h CC term can be properly deconvolved, at least in the frequency domain limited by the spectrum of f (t). Is this valid in an open medium? We have conducted the same numerical experiment, with the same distribution of scatterers, but with absorbing instead of reflecting boundary conditions. As a result, the cavity equation is no longer valid and the correlation of the scattered field C AB shows no resemblance whatsoever with the Green's function ͑the correlation coefficient function between wave forms represented on Fig. 3 is less than 0.5%͒.

However, a physical argument indicates that the Green's function can still be recovered from the correlations C AB , even in an open medium, if several judiciously distributed sources are used instead of a single point C. To that end, we propose to analyze the experiment in terms of time-reversal symmetry. Indeed, there is a strong link between correlations of a diffuse field and time reversal. [START_REF] Derode | Random multiple scattering of ultrasound. II. Is time-reversal a self-averaging process?[END_REF] Because the scatterers do not move and there is no flow within the medium, the propagation is reciprocal, i.e., h IJ (t)ϭh JI (t). When we cross-correlate the impulse responses received in A and B, the result h AC (t) h BC (Ϫt)i s equal to h CB (Ϫt) h AC (t). Now, imagine that we do the following time-reversal experiment: B sends a pulse, C records the impulse response h CB (t), time reverses it and sends it back; the resulting wave field observed in A would then be h CB (Ϫt) h AC (t), which, because of reciprocity, is exactly the cross-correlation h AC (t) h BC (Ϫt) of the impulse responses received in A and B when C sends a pulse. We would like the direct Green's function h AB to appear in the cross correlation. But in the most general case, h CB (Ϫt) h AC (t) has no reason to be equal to h AB , as was shown in Fig. 3. Yet we can go beyond: imagine now that we use several points C to perform the time-reversal operation, and that we place them in such a way that they form a perfect time-reversal device, with no loss of information. Following the Helmholtz-Kirchhoff theorem, such would be the case if the sources C were continuously distributed on a surface surrounding the scattering medium. Then the time-reversal operation should be perfect. During the ''forward'' propagation, B sends a pulse that propagates everywhere in the medium

FIG. 2. Comparison between C AB (t) ͑thick continuous line͒ and h AB (t)

h CC (Ϫt) f (t) ͑dotted line͒. The overall correlation coefficient between the two waveforms is 98.7%.

FIG. 3. Comparison between C AB (t) ͑thick continuous line͒ and h AB (t)

h CC (Ϫt) f (t) ͑dotted line͒ in the open scattering medium. The impulse responses were recorded during 800 s, until they became negligible. The overall correlation coefficient between the two waveforms is 0.48%.

͓including in A where the field received is h BA (t)], it may be scattered many times and is eventually recorded on every point of the time-reversal device, with no loss. When the field is time reversed, since nothing of it has been lost, it should exactly go backward in time ͑and refocus on B at time tϭ0) everywhere in the medium, which implies that the field received in A after the time reversal is exactly h BA (Ϫt), the time-reversed version of the direct Green's function. Then, once the wave has refocused on B ͑at time tϭ0), it does not stop since there is no ''acoustic sink'' in B: [START_REF] De Rosny | Overcoming the diffraction limit in wave physics using a time-reversal mirror and a novel acoustic sink[END_REF] the wave diverges again from B and gives rise, at times tϾ0, to h BA (t) in A.

Thus, if we use a collection of sources C arranged in such a way that they form a perfect time-reversal device, we should have

͚ C h AC ͑ Ϫt ͒ h CB ͑ t ͒ϭh BA ͑ t ͒ϩh BA ͑ Ϫt ͒. ͑2͒
A more detailed analysis, taking into account the monopolar or dipolar nature of the source/receivers, is given by Didier Cassereau. [START_REF] Cassereau | Time-reversal of ultrasonic fields-Part III: theory of the closed time-reversal cavity[END_REF] Equation ͑2͒ implies that the impulse response h BA (t) can still be recovered from the correlation of a diffuse field, even in an open medium, provided that the sources C are distributed judiciously, and all the correlation functions are summed over the source positions. Unlike the case of a closed medium, no additional deconvolution by h CC is needed. From this time-reversal analogy, we deduce a condition for the Green's function to emerge from cross correlations in open media: the sources C must be placed so that they form a perfect time-reversal device.

We have checked this in the numerical experiments depicted in Fig. 4. The results are in excellent agreement with Eq. ͑2͒, as is shown in Fig. 5: the degree of correlation between waveforms is 97.4%. Of course, when the sources are not placed as a perfect mirror, as presented in Fig. 4͑b͒, the results are less good ͑the degree of correlation between waveforms is 81.9%͒ because one part of the waves is not recorded by the time-reversal device due to the presence of scatterers outside the sources. Yet the main features of the Green's function can still be recognized, even at late times. If the number of sources is decreased, the reconstruction of the Green's function is less satisfactory, as shown in Fig. 6. With only 50 sources ͑instead of 250 previously͒ regularly spaced every ϳ5/3 all along a circle as in Fig. 4͑a͒, the correlation coefficients between waveforms is 70%. However, if the 50 sources are gathered together in a 72°angular sector ͑pitch /3͒, it drops to 53%. Indeed, since the coherence length of a diffuse wave field is ϳ, it is useless to place the sources closer.

So far, we have considered that the origin of the field measured in A and B was an active and coherent source transmitting a short pulse ͑or a collection of such sources͒. What if there are no such sources in the medium, but a diffuse continuous noise? The physical origin of this noise may be thermal vibrations. [START_REF] Weaver | Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies[END_REF] In seismology, noise in the seismograms comes from a variety of different sources ͑traffic, sea waves, weather, human activity,...͒ continuously and ͑alleg-edly͒ randomly pumping energy into the earth and essentially exciting surface waves. In ocean acoustics, noise may originate from boats, surf, wind, animals, etc. By definition, the noise sources cannot be controlled. In the light of the discussion above, in order to recover the Green's function from the cross-correlation of the noise received in A and B, the most favorable situation would be that in which noise can be con- sidered as coming from virtual point sources C randomly distributed everywhere in the medium and continuously generating uncorrelated white noises n C (t). In that case, the cross correlation between waveforms sensed in A and B would be

͚ C ͚ CЈ h AC ͑ Ϫt ͒ n C ͑ Ϫt ͒ h C Ј B ͑ t ͒ n C Ј ͑ t ͒.
If the observation time ⌬T is long enough compared to the correlation time of the noise, then n C (Ϫt) n C Ј (t) will con- verge to ␦(t)␦ C,C Ј . Moreover, if the virtual noise sources C are distributed everywhere in the medium ͑in other words, if each degree of freedom is excited randomly and indepen-dently͒, then the C-points would necessarily constitute a perfect time-reversal device, so Eq. ͑2͒ should be verified again.

We have carried out a numerical experiment based on this idea. Two hundred forty sources were distributed at random inside the open scattering medium shown in Fig. 1, and 240 uncorrelated white noises, convolved by e(t), were transmitted by these sources during 40 ms. The resulting wave forms are received in A and B. Their cross correlations C AB (t) is compared to the direct Green's function h AB (t) f (Ϫt): the agreement is still very good ͑61% correlation coefficient͒ even at late times.

The emergence of the Green's function in the field-field correlations in a closed cavity with discrete modes is now well established. [START_REF] Weaver | On the emergence of the Green's function in the correlations of a diffuse field: pulse echo using thermal phonons[END_REF][START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF][START_REF] Weaver | Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies[END_REF] In this letter, we have argued that recovering the Green's function was also possible in an open multiple scattering medium and we have proposed a criterion based on reciprocity, time-reversal symmetry, and the Helmholtz-Kirchhoff theorem: if sources are placed as if they were to form a perfect time-reversal device, then the Green's function can be recovered by summing the cross correlations. This has been validated by numerical experiments. The reduction of the number of sources was also discussed and the possibility of using noise sources was illustrated.

There is still much food for thought, particularly regarding the role of the scatterers in the reconstruction of the Green's function. The argument we developed here is valid for any medium ͑homogeneous, high-order multiple scattering, reverberant, ...͒ where reciprocity and invariance under time reversal hold. The field does not need to be thermally diffuse for the Green's function to emerge from the correlations, as long as there are enough well-positioned sources. Another approach is to consider the scatterers as secondary sources which are necessary to truly randomize the wave field emanating from a single original source. Correlating the late part of the coda would then permit us to reconstruct at least the early arrivals of the Green's function. Recent seismologic results support this idea. [START_REF] Campillo | Long range correlations in the diffuse seismic coda[END_REF] The influence of absorption is also to be investigated.

FIG. 1 .

 1 FIG. 1. Two hundred perfectly rigid scatterers ͑diameter 2.1 mm͒ are randomly distributed over a 7.5ϫ7.5 cm 2 area. A point source is placed in C; A and B are receiving points. The boundary conditions may be perfectly reflecting ͑Dirichlet͒ as in a closed cavity or absorbing ͑open medium͒.

FIG. 4 .

 4 FIG. 4. A and B are receiving points. Two hundred fifty source points are placed regularly around a circle with radius 18.7 mm, 100 scatterers being inside the circle. They completely surround the medium ͑a͒, or only partially ͑b͒. The boundary conditions on the edges of the grid are absorbing ͑open medium͒, in both cases.

ACKNOWLEDGMENTS

The authors wish to thank Ludovic Margerin and Anne Paul.


domain: phys chim

text:
Octave Duros

Hebatalla Elnaggar

Amélie Juhin

Christian Brouder

General expressions for Stevens operators and spin operator equivalents

Keywords:

Several definitions of the crystal electric field (CEF) have been used over time and their variety has lead to many misunderstandings, in both theoretical and experimental literature. Two categories of definitions can be mentioned, the first being the operators equivalents introduced by Stevens in 1952 and the second being the crystal-field operators, introduced later by Wybourne and expanded on the Racah spherical tensors. This paper aims at providing some clarification in this field. We first make a review of several expressions introduced in various references to compute crystal-field operators and we make connections between them. Then, we introduce an explicit way to compute crystal-field operators, in terms of j and m as well as in terms of J 2 and Jz. We eventually give some connections between the Stevens operators equivalents and the crystal-field operators.

I. INTRODUCTION

Operator equivalents O m n were introduced in 1952 by Stevens [START_REF] Stevens | Proc. Phys. Soc. A[END_REF] to parametrize the Hamiltonian describing magnetic properties of rare earths. Indeed the 4f shell of most rare earths can be described by a spectroscopic term 2S+1 L J , where J is a good quantum number. As a result, its ground state belongs to the vector space generated by state vectors |JM (where M = -J, . . . , J) and the action of the crystal field can be written as a linear combination of operators O m n that can be expressed as a polynomial in the spin operators J x , J y and J z . These Stevens operators were then widely used to analyze spectroscopic experiments in EPR, ESR, ENDOR, UV-visible spectroscopy, Mössbauer spectroscopy 2 as well as neutron scattering and x-ray absorption spectroscopy. Explicit expressions for O m n were given in several tables and, because of the complexity of these expressions, each new table listed the misprints in the previous ones [START_REF] Rudowicz | [END_REF] .

On the other hand, it was later realized 3 that the transformation properties of O m n under rotation are very complicated. By constrast, spherical tensors 4 T k µ (with µ = -k, . . . , k) enjoy the same transformation properties as spherical harmonics and are linked by the relation [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] (p. 90): J -, T k µ = (k + µ)(k -µ + 1)T k µ-1 whereas O m n of different m have no simple relations. Moreover, any (2j + 1) × (2j + 1) matrix can be written as a linear combination of T k µ for k = 0, . . . , 2j [START_REF] Grenet | [END_REF] . In particular, spherical tensors can be used to express any finite dimensional Hamiltonian. For all these reasons, many authors replaced Stevens operators by spherical tensors or combinations thereof, leading to new tables. Besides, additional confusion arises from the fact that spherical tensors come with various normalizations in the literature.

The operator equivalents were introduced in 1952 by Stevens was given by Bleaney and Stevens 7 . It was later realized 3 that the transformation properties of O m n would be much simpler if spherical tensors 4 were used. This lead to the introduction of a large number of more or less similar operators, that were reviewed by Rudowicz and coll. [START_REF] Rudowicz | [END_REF]8 .

The absence of an analytical expression for these operator equivalents imposed the publication of extensive tables running over many pages and scattered in several articles and books. Moreover, each single entry can be quite complicated (for example the expression for O [START_REF] Stevens | Proc. Phys. Soc. A[END_REF] 12 is 11 lines long in 2 ). To avoid egregious mistakes some authors proposed to derive these results with a computer through recursion methods 9,10 .

To solve these problems, several closed form expressions were proposed 11 . In this paper we make a critical review of these expressions and we present a new method to translate an analytical expression in terms of the angular moment quantum number j into an analytical expression in terms of the angular momentum operator J 2 . Finally, we describe the connection between the parameters of the effective Hamiltonian and the parameters of the crystal-field Hamiltonian.

The plan of the paper is the following. We start from the definition of spherical tensors T k µ and we introduce the operator P k µ which has the advantage of having only diagonal matrix elements ( jm|P k µ |jm = 0 if m = m ) and from which T k µ can be easily computed. Then, we review the different expressions for P k µ proposed in the literature and we make connections between them. We also provide a way to go from an expression of jm|P k µ |jm in terms of j and m to an expression of P k µ in terms of J 2 and J z . Finally, we give some applications to checking tables, connecting with various Stevens operators and linking effective Hamiltonians to crystal field operators.

II. SPHERICAL TENSORS AND THEIR MATRIX ELEMENTS

A. Definition and main properties General (Irreducible) spherical tensors T k µ (with µ = -k, . . . , k) are defined by the relations 5 (p. 157)

T k k = a k J k + , (1)
J -, T k µ = (k + µ)(k -µ + 1)T k µ-1 , (2)
where a k is a real constant, J ± = J x ± iJ y and J k + = (J + ) k (i.e. k is a true exponent and not an upper index as in T k µ ). These operators satisfy the important properties 5 (p. 656)

(T k µ ) † = (-1) µ T k -µ , (3)
and

k µ=-k (T k µ ) † T k µ = a 2 k (k!) 2 (2k)! k s=1 (4J 2 + 1 -s 2 ), (4)
where J 2 = J 2 x + J 2 y + J 2 z . Moreover if, for a given j, we denote by M k µ the (2j +1)×(2j +1) matrix defined by the matrix elements (M k µ ) mm = jm|T k µ |jm , then we have the orthogonality relation [START_REF] Grenet | [END_REF] :

Tr (M k µ ) † M k µ = δ k,k δ µ,µ (k!) 2 a 2 k (2k + 1)! (2j + k + 1)! (2j -k)! , (5)
where the adjoint † can be replaced by a transpose since the matrix elements are real. This property can be used for instance to determine the expansion of a Hamiltonian over spherical tensors. Indeed, if H is a (2j + 1)×(2j + 1) Hamiltonian matrix, then the parameters B µ k of

H = k,µ B µ k T k µ ,
are given by

B µ k = (2k + 1)! (k!) 2 a 2 k (2j -k)! (2j + k + 1)! Tr (M k µ ) † H .

B. Common conventions

We denote by T k µ the spherical tensor corresponding to 5

a k = (-1) k (2k)! k! .
Another frequently used convention [START_REF] Tuszynksi | Spherical Tensor Operators -Tables of Matrix Elements and Symmetries[END_REF] is

a k = (-1) k 2 -k/2 .
Let us finally define Racah's unit tensor operators 4,13 u k q (j), for which a k depends on the value of j on which the tensor operators are applied:

a k = (-1) k k! (2k)!(2j -k)! (2j + k + 1)! 1/2
.

C. Matrix elements

Any Hamiltonian acting within a vector space generated by |jm m=-j,...,j can be written as a sum

H = k,µ B µ k T k µ .
Please note that the position of the angular momentum quantum number k and the magnetic quantum number µ are traditionally at a different position in B µ k or O µ k and T k µ . The matrix elements of the Hamiltonian can be derived from the matrix elements of T k µ , which are given by the Wigner-Eckart theorem:

jm |T k µ |jm = (-1) j-m j k j -m µ m j||T k ||j .
To compute the reduced matrix element j||T k ||j we only need

T k k |jm = a k J k + |jm = a k (j -m)!(j + m + k)! (j -m -k)!(j + m)! 1/2 |jm + k , (6)
and the explicit expression 14 (p. 268)

j k j -m k m = δ m ,m+k (-1) j-m (2k)!(2j -k)!(j + m + k)!(j -m)! (2j + k + 1)!(k!) 2 (j -m -k)!(j + m)! 1/2 , which give us j||T k ||j = (-1) k a k k! (2j + k + 1)! (2k)!(2j -k)! 1/2 . (7)
In particular, we see that Racah's unit operators derive their names from the reduced matrix element j||u k (j)||j = 1. Note that Eq. ( 7) differs from the reduced matrix element given in 5 (p. 157) by a factor of √ 2j + 1 due to the fact that these authors use a Wigner-Eckart theorem involving Clebsch-Gordan coefficients instead of 3j-symbols. Note that the orthogonality relation (5) follows from the Wigner-Eckart theorem and the orthogonality of 3-j symbols.

III. EXPLICIT EXPRESSIONS FOR T k µ

In this section we assume that µ ≥ 0. The first step to obtain a explicit expression is to write T k µ as the product of J µ + times a polynomial in J 2 and J z . This polynomial will only have non-zero matrix elements between vector states with the same m and j. The second step is to obtain explicit expressions for this polynomial.

A. P k µ and its properties

To derive explicit expression for T k µ we first define P k µ for µ ≥ 0 by 5 (p. 658)

T k µ = a k (k + µ)!(k -µ)! (2k)! J µ + P k µ , (8)
with P k k = 1. We first show that P k µ is a polynomial of J 2 and J z . Indeed, the recursive relation (2) for T k µ leads to the recursive expression for

P k µ (k -µ + 1)P k µ-1 = P µ µ-1 P k µ + J + J -, P k µ ,
where

P µ µ-1 = -µ(µ -1) -2µJ z .
To give a recursive proof of the polynomial nature of P µ , we start from the fact that P k k = 1 is a polynomial in J 2 and J z . Assume now that P k µ is a polynomial in J 2 and J z . Then, by using 9 [J -, J n z ] = J -∇J n z , where ∇J n z = J n z -(J z -1) n is the backward difference, we see that J -, P k µ is the product of J -by a polynomial in J 2 and J z . It follows that J + J -, P k µ and P k µ-1 are also polynomials of this form because J + J -= J 2 -J 2 z -J z . As a consequence, the operator P k µ is a polynomial in J 2 and J z for all µ such that k ≥ µ ≥ 0, which implies that its matrix elements are diagonal in j and m: jm|P k µ |j m = δ jj δ mm jm|P k µ |jm . To be more explicit, if P k µ (J 2 , J z ) is the operator P k µ whose dependence on J 2 and J z is made explicit, we denote by P k µ (J 2 , ∇J z ) the operator P k µ (J 2 , J z ) where each monomial J n z is replaced by ∇J n z . Then, the recursive expression becomes

(k -µ + 1)P k µ-1 (J 2 , J z ) = -(µ(µ -1) + 2µJ z )P k µ (J 2 , J z ) + (J 2 -J 2 z -J z )P k µ (J 2 , ∇J z )
. Finally, we give the formula for T k -µ following Ryabov 10

T k -µ = (-1) k a k (k + µ)!(k -µ)! (2k)! J µ -P k µ (J 2 , -J z ).
If we combine this with (T k µ ) † = (-1) µ T k -µ , J † + = J -and the fact that P k µ (J 2 , J z ) is self-adjoint, we obtain the remarkably simple intertwining formula

P k µ (J 2 , J z )J µ -= (-1) k-µ J µ -P k µ (J 2 , -J z )
. We now list the existing closed-form expressions for the matrix elements of P k µ . Note that each expression corresponds to a formula for the corresponding Clebsch-Gordan coefficients.

B. Coala

The first closed-form expression was given by Coala [START_REF] Coala | [END_REF] jm

|T k µ |jm = (-1) k (2j + 1)(k!) 2 (2j -k)! (2j + k + 1)! 1/2 (k + µ)!(k -µ)! k p=µ (-1) p (j -µ -m)!(j + m)! (j -m -p)!(j -k + m + p)!p!(k -p)!(p -µ)!(k + µ -p)! jm |J µ + |jm ,
where we corrected a misprint. By comparing this expression with Eq. (3.170) in [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] (p. 79) for the Clebsch-Gordan coefficients (jmkµ|jm ) we see that

jm |T k µ |jm = (jmkµ|jm ), so that a k = (-1) k k! (2j + 1)(2k)!(2j -k)! (2j + k + 1)! 1/2 .
Therefore,

jm|P k µ |jm = (-1) k (k!) 2 k p=µ (-1) p (j -µ -m)!(j + m)! (j -m -p)!(j -k + m + p)!p!(k -p)!(p -µ)!(k + µ -p)! .
By using the descending Pochhammer symbol [a] p defined for p ≥ 0 by By using the descending Pochhammer symbol [a] p defined for p ≥ 0 by

[a] p = a! (a -p)! , (9)
we can rewrite this as

jm|P k µ |jm = (-1) k k p=µ (-1) p k p k p -µ [j + m] k-p [j -m -µ] p-µ . (10)
C. Biedenharn and Louck Table 6 of 5 gives the following relation for the non-zero matrix elements of P k µ :

jm|P k µ |jm = k-µ p=0 (-1) p k p k µ + p [j + m] p [j -m -µ] k-µ-p . (11)
where [a] p is the descending Pochhammer symbol for p ≥ 0

[a] p = a! (a -p)! . ( 12
)
Note that, in the rest of their book Biedenharn and Louck use the ascending Pochhammer symbol 5 (p. 352). We notice that the expressions given by Coala and by Biedenharn and Louck become identical if we write q = k -p.

D. Grenet and Kibler

Again for µ ≥ 0, the formula given by Grenet and Kibler is 6

jm|P k µ |jm = k! (2j -k)!(k + µ)! k z=µ (-1) µ+z (2j -z)!(k + z)! z!(k -z)!(z -µ)! (j + m)! (j + m + µ -z)! ,
that we rewrite

jm|P k µ |jm = (k!) 2 (k + µ)! k-µ p=0 (-1) p p! k + µ + p k 2j -µ -p 2j -k [j + m] p . (13)
A remarkable aspect of this formula is that m appears only in one term [j + m] p . This allows us to obtain a relatively simple expression for P k µ as a function of J z . Indeed, if we expand the Pochhammer symbol in terms of (signed) Stirling numbers of the first kind s(p, n):

[a] p = p q=0 s(p, q)a q . ( 14
)
Hence

[j + m] p = p q=0 s(p, q)(j + m) q = p n=0 m n p-n q=0 n + q n s(p, n + q)j q (15)
And we obtain an explicit expression of P k µ in terms of j and J z :

P k µ = k-µ n=0 α k µn (j)J n z , (16)
where

α k µn (j) = (k!) 2 (k + µ)! k-µ p=n (-1) p p! k + µ + p k 2j -µ -p 2j -k p-n q=0 n + q n s(p, n + q)j q . ( 17
)

E. Other expressions

Since matrix elements of T k µ are directly related to Clebsch-Gordan coefficients, every Clebsch-Gordan coefficient formula gives an expression for the matrix elements of T k µ . A large number of such formulas are given in [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF] and Shimpuku gives 168 different formulas in 16 . The formulas used in the literature to calculate T k µ correspond to the cases where there are a small number of terms using the variable m, so that the expression in terms of J z is (relatively) simple. From that point of view, the formula given by Grenet and Kibler cannot be improved.

There are also several semi-explicit expressions using iterated forward differences, quasi-powers or hypergeometric functions [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF] . Within this semi-explicit family, a series of authors [17][18][19][20][21] used a connection between some Clebsch-Gordan coefficients and Chebyshev polynomials of a discrete real variable (not to be confused with the Chebyshev polynomials T n (x) and U n (x) of the first and second kind). Then, they used iterated forward differences of these polynomials to obtain the Clebsch-Gordan coefficients required to calculate the matrix elements of T k µ . We do not elaborate because semi-explicit formulas are not the subject of this paper and because a very nice and extensive article was recently published on this topic 21 .

F. From matrix elements to operators

We proved that P k µ is a polynomial of J 2 and J z . We gave an explicit expression for P k µ as a polynomial of J z with coefficients α k µm (j) which are polynomials of j. It is interesting to use expressions for P k µ as a polynomial of J z and J 2 because these expressions are independent of j and some tables are given in terms of J 2 instead of j. To replace the dependence on j by a dependence on J 2 , we need to find polynomials β k µm (J 2 ) such that

β k µm (J 2 )|jm = β k µm j(j + 1) |jm = α k µm (j)|jm .
We treat this problem in a general way and for notational convenience we replace j by x and we solve the following problem: let P and Q be polynomials related by P (x) = Q x(x + 1) , calculate the coefficients of Q from the coefficients of P . Assume that

P (x) = n i=0 a i x i = Q x(x + 1) = m j=0 b j x(x + 1)) j .
We expand

Q x(x + 1) = m j=0 j k=0 b j j k x 2k x j-k = m j=0 j k=0 b j j k x j+k .
In the last expression we define i = j + k, so that k = i -j. We know that j ≥ k ≥ 0 so that i/2 ≤ j ≤ i. But we also know that 0 ≤ j ≤ m. Thus

Q x(x + 1) = 2m i=0 x i min(m,i) j= i/2 b j j i -j ,
and

a i = min(m,i) j= i/2 b j j i -j .
The term of maximum degree is a 2m = b m , so that n = 2m. If we want to express b j in terms of a i , we need only the first m + 1 values a i . Moreover, we have a 0 = b 0 and the interesting cases are

a i = m j=1 M ij b j ,
where the m × m matrix M = (M ij )

M ij = j i -j for i ≥ j ≥ i/2 , M ij = 0 otherwise, is invertible because it is lower diagonal with diagonal M ii = 1.
Since M is a lower triangular matrix with unit diagonal, its inverse is also a lower triangular matrix with unit diagonal. To prove this, write M = 1 -K, where K is lower triangular with zero diagonal so that K m = 0 and note that M -1 = I + m-1 k=1 K k . Moreover, the first k lines of the inverse of a triangular matrix M depend only on the first k lines of M . In this sense, the inverse of M does not depend on the dimension of M .

It turns out that this problem is solved a theorem given by Riordan [START_REF] Riordan | Combinatorial Identities[END_REF] (p. 50) : If

α i = j (-1) j β j p + qj -j i -j , then β i = j (-1) j α j p + qj -j p + qi -j p + qi -j i -j .
In the case α i = a i , β j = (-1) j b j , p = 0 and q = 2 we recover our equation

a i = j b j j i -j ,
and Riordan's theorem gives the inverse relation:

b 0 = a 0 , b i = i j=1
(-1) i-j j 2i -j 2i -j i -j a j for i ≥ 1.

G. Explicit expressions in terms of j and m

As an application of the formulae presented in this paper, we give explicit expressions for P k µ in terms of j and m for k = 0 to 6. P 0 0 = 1.

P 1 0 = -2m, P 1 1 = 1. P 2 0 = -2j 2 -2j + 6m 2 , P 2 1 = -4m -2, P 2 2 = 1. P 3 0 = 12j 2 m + 12jm -20m 3 -4m, P 3 1 = -3j 2 -3j + 15m 2 + 15m + 6, P 3 2 = -6m -6, P 3 3 = 1. P 4 0 = 6j 4 + 12j 3 -60j 2 m 2 -6j 2 -60jm 2 -12j + 70m 4 + 50m 2 , P 4 1 = 24j 2 m + 12j 2 + 24jm + 12j -56m 3 -84m 2 -76m -24, P 4 2 = -4j 2 -4j + 28m 2 + 56m + 36, P 4 3 = -8m -12, P 4 4 = 1.
Then,

j, m + µ|J µ + P k µ |jm = (j -m)!(j + m + µ)! (j -m -µ)!(j + m)! 1/2 (j + m)!(j -m -µ)! (2j -k)! k-µ p=0 (-1) p k p k µ + p 2j -k j + m -p = ((j -m)!(j + m + µ)!(j + m)!(j -m -µ)!) 1/2 (2j -k)! k-µ p=0 (-1) p k p k µ + p 2j -k j + m -p . Hence jm + µ|T k µ |jm = (-1) k (k + µ)!(k -µ)! k! jm + µ|J µ + P k µ |jm = (-1) k ((k + µ)!(k -µ)!(j -m)!(j + m + µ)!(j + m)!(j -m -µ)!) 1/2 (2j -k)!k! k-µ p=0 (-1) p k p k µ + p 2j -k j + m -p = C jkj mµ,m+µ (2j + k + 1)! (2j + 1)(2j -k)! 1/2
. This gives us the Clebsch-Gordan coefficient

C jkj mµ,m+µ = (-1) k (2j + 1)(k + µ)!(k -µ)!(j -m)!(j + m + µ)!(j + m)!(j -m -µ)! (2j + k + 1)!(2j -k)!(k!) 2 1/2 k-µ p=0 (-1) p k p k k -µ -p 2j -k j + m -p .
We need to find a place where this expression for the Clebsch-Gordan coefficient is given. We first use

C jkj mµ,m+µ = (-1) k C kjj µm,m+µ .
The coefficients of Shimpuku's expression becomes

a 1 = 2j -k, a 2 = k, a 3 = k, b 1 = k -µ, b 2 = j -m, b 3 = j + m + µ, c 1 = k + µ, c 2 = j + m, c 3 = j -m -µ.
Then the sum (p. 413)

α (-1) α a 3 α a 2 b 1 -α a 1 c 2 -α = α (-1) α k α k k -µ -α 2j -k j + m -α ,
is exactly what we need. We must check that the overall factors agree. Indeed:

b1+c1 a3 b2+c2 a3 J+1 a3 b1+c1 b1 b2+c2 c2 b3+c3 c3 = 2k k 2j k 2j+k+1 k 2k k-µ 2j j+m 2j j-m-µ = (2j + 1)(k + µ)!(k -µ)!(j -m)!(j + m + µ)!(j + m)!(j -m -µ)! (2j + k + 1)!(2j -k)!(k!) 2 1/2 .
To summarize. The expression for P k µ cannot be obtained from the usual formulas for the Clebsch-Gordan coefficients. Among the 168 expressions for the Clebsch-Gordan coefficients proposed by Shimpuku 16 , only the first one of p. 413 provides the proper result

C kjj µm,m+µ = 2k k 2j k 2j+k+1 k 2k k-µ 2j j+m 2j j-m-µ 1/2 p (-1) p k p k k -µ -p 2j -k j + m -p .
We saw that

j, m + µ|T k µ |jm = C jm+µ jmkµ (2j + k + 1)! (2j + 1)(2j -k)! 1/2 . Therefore, j, m + µ|T k µ |jm = 1 (2j -k)! (j + m)!(j -m)!(k -µ)!(j + m + µ)! (k + µ)!(j -m -µ)! 1/2 z (-1) k+µ+z (2j -z)!(k + z)! z!(j + m + µ -z)!(k -z)!(z -µ)! .
We know that The identity of Eqs. ( 11) and ( 13) was checked with Mathematica for k = 0 to 12. This can be rewritten

jm|P k µ |jm = (k!) 2 (k + µ)! k-µ s=0 (-1) s k + µ + s k 2j -µ -s k -µ -s j + m s .
On the other hand, Eq. ( 11) can be rewritten

  + µ|T k µ |jm = (-1) k (k + µ)!(k -µ)! k! j, m + µ|J µ + P k µ |jm = (-1) k (k + µ)!(k -µ)! k! j m j, m + µ|J µ + |j m j m |P k µ |jmBy using Eq. (6) and the orthonormality of the basis vectors we obtainj, m + µ|T k µ |jm = (-1) k (k + µ)!(k -µ)! k! (j -m)!(j + m + µ)! (j -m -µ)!(j + m)! 1/2 jm|P k µ |jm .Our final result is now, for µ ≥ 0 and k ≥ 2jjm|P k µ |jm = k! (2j -k)!(k + µ)! z (-1) µ+z (2j -z)!(k + z)! z!(k -z)!(z -µ)! (j + m)! (j + m + µ -z)! , k + µ + s s (2j -µ -s)! (2j -k)!(j + m)! (j + m -s)! .

(- 1 )

 1 m)! (j + m -s)! (j -m -µ)! (j -m -k + s)! = -µ k -µ -s (k -µ -s)! = (k!) 2 (k + µ)! k-µ s=0 (-1) s k + µ µ + s j + m s j -m -µ k -µ -s .The two formulas are identical ifk-µ s=0 s k + µ + s k 2j -µ -s k -

P 5 0 = -60j 4 m -120j 3 m + 280j 2 m 3 + 140j 2 m + 280jm 3 + 200jm -252m 5 -420m 3 -48m, P [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] 1 = 10j 4 + 20j 3 -140j 2 m 2 -140j 2 m -70j 2 -140jm 2 -140jm -80j + 210m 4 +420m 3 + 630m 2 + 420m + 120, P [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] 2 = 40j 2 m + 40j 2 + 40jm + 40j -120m 3 -360m 2 -480m -240, P [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] 3 = -5j 2 -5j + 45m 2 + 135m + 120, P [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] 4 = -10m -20, P [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] 5 = 1.

P 6 0 = -20j 6 -60j 5 + 420j 4 m 2 + 100j 4 + 840j 3 m 2 + 300j 3 -1260j 2 m 4 -1680j 2 m 2 -80j 2 -1260jm 4 -2100jm 2 -240j + 924m 6 + 2940m 4 + 1176m 2 , P [START_REF] Grenet | [END_REF] 1 = -120j 4 m -60j 4 -240j 3 m -120j 3 + 720j 2 m 3 + 1080j 2 m 2 + 1200j 2 m + 420j 2 + 720jm 3 +1080jm 2 + 1320jm + 480j -792m 5 -1980m 4 -4320m 3 -4500m 2 -2808m -720, P [START_REF] Grenet | [END_REF] 2 = 15j 4 + 30j 3 -270j 2 m 2 -540j 2 m -375j 2 -270jm 2 -540jm -390j + 495m 4 + 1980m 3 +4095m 2 + 4230m + 1800, P [START_REF] Grenet | [END_REF] 3 = 60j 2 m + 90j 2 + 60jm + 90j -220m 3 -990m 2 -1790m -1200, P 6 4 = -6j 2 -6j + 66m 2 + 264m + 300, P 6 5 = -12m -30, P 6 6 = 1.

H. Explicit expressions in terms of J 2 and Jz

As an application of the formulae presented in this paper, we give explicit expressions for P k µ in terms of J 2 and J z for k = 0 to 6.

4 = -10(J z + 2), P 5 5 = 1. is proportional to J z +ν. As far as we know, this observation is new and we ignore whether it can be extended to all values of k. As a consequence, P 2m 2ν+1 |j, -2ν+1 2 = 0 and P 2m+1 2ν |j, -ν = 0 for all values of j. This restricts the influence of the corresponding crystal field parameters.

IV. FROM CRYSTAL FIELD PARAMETERS TO OPERATOR EQUIVALENTS

When investigating rare earths, neutron scattering generally uses operator equivalents whereas x-ray absorption spectroscopy uses crystal field parameters. More precisely, the crystal field acting on the angular part of the 4f shell is described by

whereas operator equivalents are

A. One-electron case

As a simple case, we consider the case where the shell is occupied by a single electron. We use formula (40) of 14 (p. 481)

Moreover,

where [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF] (p. 251)

The relation between C k µ acting on | m and T k µ is

B. Many-electron case

We consider now the case where shell is occupied by n electrons (0 ≤ n ≤ 2(2 + 1)). For the case of a shell containing n electrons, we use Eq. (11.53) [START_REF] Cowan | The Theory of Atomic Structure and Spectra[END_REF] (p. 317)

The coefficients of fractional parentage were obtained from the table computed by Velkov in his PhD thesis [START_REF] Velkov | Multi-electron coefficients of fractional parentage for the p, d, and f shells[END_REF][START_REF] Judd | [END_REF] .

V. CONCLUSION

In this paper, we reviewed the existing expressions for the spherical tensor form of Stevens or spin operator equivalents. We gave explicit proofs of these forms, which were not available in the literature and we devised a way to transform expressions in terms of j and m into expressions in terms of J 2 and J z . Along the way, we discovered factorization properties of some of the spherical tensors which seem to be new.

We hope that the present paper will enable the physics and chemistry communities to avoid having rescourse to cumbersome (and generally not fully correct) tables.

Among the perspectives of this work, we might mention the question whether or not the factorization properties extend to general values of k. We also want to mention that some analytical expressions were given for non-diagonal tensor-equivalent operators [25][26][27] , i.e. for the case of matrix elements between different values of j. However, the generalization of our results to the non-diagonal case seems non trivial.

VI. ACKNOWLEDGEMENTS

We are very grateful to Genviève Grenet for pointing out to us that the main formula of her article with Kibler is proved in her PhD thesis. We are also grateful to the Online Encyclopedia of Integer Sequences (https://oeis.org) and to Marjorie Bicknell Johnson whose paper with Hoggatt 28 lead us to Riordan's theorem. We thank David Siminovitch for providing information about his paper 21 . Finally, we are very grateful to Dobromir D. Velkov who sent us a copy of his PhD thesis [START_REF] Velkov | Multi-electron coefficients of fractional parentage for the p, d, and f shells[END_REF] and of its 2891-page-long Appendix A containing the one-, two-, three-and four-body coefficients of fractional parentage, that we used to calculte the matrix elements of U (k) .

VII. PROOFS

As far as we know, no proof of the formulas given by Biedenharn and Louck and by Grenet and Kibler were published, although Genevi'eve Grenet provides a proof in her PhD thesis. Since these proofs depend on rather unusual expressions for the Clebsch-Gordan coefficients, we find it useful to make them explicit. In this section we use the spherical tensors T k µ normalized by 5

A. Proof of Coala's formula Equation ( 14) of Ref. [START_REF] Coala | [END_REF] is, in our notation,

We recall that

Therefore,

To check Coala's formula, we notice that the reference he uses seems to be the symmetric form, Eq. (3.170) in [START_REF] Biedenharn | Encyclopedia of Mathematics and its Applications[END_REF] .

Let us check that. According to Eq. (3.170) in 5 (p. 79).

(j

The relevant coefficient is (jmkµ|jm + µ), so that

Proof of the Biedenharn and Louck formula

We want to prove the expression for jm|P k µ |jm . For this we need an expression for the Clebsch-Gordan coefficient which is adapted to what we want. Recall that,

and that

The formula for the matrix element of P k µ given by Eq. ( 11) is

By using C jkj mµ,m+µ = (-1) k C kjj µm,m+µ we obtain

The Wigner-Eckart theorem for T k µ is expressed by

From the relation