New insights for the description of magnetic correlations inferred from muSR

Whenever a compound exhibits a spontaneous muSR oscillation, long-range magnetic ordering is usually inferred. Here we show that some caution is required. The coherence length needs not to be large for a spontaneous muon spin precession to be observed. Establishing the incommensurate nature of a magnetic structure, solely based on muSR measurements, may not be reliable. The absence of a spontaneous muon precession at low temperature does not mean that the system under investigation does not display long-range magnetic ordering. The relaxation measured in zero and longitudinal field in the quasi-static limit is usually analyzed in the framework of the strong-collision model, the static polarization function being taken to be the famous Kubo-Toyabe function. This might not be satisfactory if short-range correlation effects are strong. Here we propose a method based on the maximum entropy concept and reverse Monte Carlo technique which gives results consistent with those obtained in 2013 by analytical means for the considered example.


Introduction
In this report we consider problems relevant for the analysis of experimental positive muon spectroscopy (µSR) data recorded for a magnetic compound. First, we discuss the meaning of an observed spontaneous muon precession in relation to the possible long-range nature of the magnetic order. Using published examples, we point out the difficulty of reliably establishing solely with µSR the incommensurate nature of the magnetic structure. In Sec. 4 we present a numerical method for the determination of the field distribution at the muon site for a compound characterized by quasi-static spin dynamics. As an example, we apply the method for the quantum spin-ice system Yb 2 Ti 2 O 7 . A discussion of our results is presented in Sec. 5. An appendix details the reverse Monte Carlo algorithm.

Long-versus short-range magnetic order
A large fraction of the µSR studies are performed with the purpose of determining whether a magnetic material exhibits a magnetic transition. If a spontaneous µSR oscillation is detected, the material is usually proposed to be characterized by a long-range magnetic order. This may not always be justified as it was shown theoretically in Ref. [1]. Let us consider the polarization function P stat Z (t) = exp(−γ 2 µ ∆ 2 Z,m t 2 /2) cos(γ µ B 0 t). It describes a conventional Gaussian damped  Asymmetry: Figure 3. A zero-field µSR spectrum for a powder of NiGa 2 S 4 taken at 2.3 K; adapted from Ref. [2]. The solid line results from a fit.
A magnetic neutron powder diffraction pattern for NiGa 2 S 4 resulting from the difference between the patterns recorded at 1.5 and 15 K. The magnetic reflections are quite broad. They correspond to a correlation length ξ m = 2.5 (3) nm. The solid line results from a fit. Adapted from Ref. [3]. spontaneous oscillation. In Fig. 1 is plotted P stat Z (t) for two values of the ratio ∆ Z,m /B 0 . From this plot we deduce that this ratio should be smaller than ≈ 1 to conclude a spontaneous precession is detected. It is theoretically possible to express the ratio in terms of the correlation length of the magnetic structure ξ m . Assuming a ferromagnetic order in a face-centered-cubic compound with lattice parameter a and the muon in the octahedral site, we compute ∆ Z,m /B 0 versus √ 2ξ m /a. The result is displayed in Fig. 2. A relatively short ξ m of about ten interatomic distances is enough to give rise to a detectable spontaneous field. This result is qualitatively independent of the type of magnetic structure and the muon site. Are ten atomic distances long enough to call the magnetic structure long-range? In Fig. 3 we display a zero-field µSR spectrum recorded for the triangular system NiGa 2 S 4 . A spontaneous oscillation is clearly observed at 2.3 K. Yet, as shown in Fig. 4 the system only exhibits short-range correlations at low temperature. In fact, the lack of long-range magnetic order is the main interest of NiGa 2 S 4 .
To round up our discussion of the meaning of a damped oscillation for a zero-field µSR asymmetry, we point out that we have recently shown that such a signal may in fact reflect shortrange correlations [4]. We recall that the Gaussian component field distributions as assumed by Kubo and Toyabe explicitly neglect them. A simple analytical method to account for them is to extend the component field distribution beyond the Gaussian function. This method successfully describes the spectrum recorded for the spin-ice system Yb 2 Ti 2 O 7 at low temperature [4].

Commensurate versus incommensurate magnetic order
An increasing number of observed spontaneous oscillations are found to be well accounted for by a zeroth-order Bessel function of the first kind, rather than a simple cosine. Usually from this observation it is inferred that the magnetic structure of the compound under study is incommensurate. However, the observation of such oscillations does not guarantee the inference to be correct. To support our argument, we cite two counter examples. A relatively complex yet commensurate structure was found in 1992 to yield a field distribution at a nuclear probe reminiscent of the one expected for an incommensurate modulated structure [5]. Recently, the triangular system Ag 2 NiO 2 previously suggested to be incommensurate using zero-field µSR data was discovered to be commensurate by neutron diffraction [6,7]. While the determination of a magnetic structure solely from µSR is questionable, it is pertinent to test the reality of a neutron-determined structure by zero-field µSR measurements. This is so because of the difficulty sometimes encountered to determine a magnetic structure from neutron diffraction data. In fact, a good practice is to record zero-field µSR and neutron diffraction data for the same sample and to extract the magnetic structure from a combined analysis of the data. We stress that information on ξ m can be derived from neutron diffraction and that the field distribution at the muon site depends on it. This is further discussed in the next section.

Numerical determination of the field distribution at the muon site
Here we further consider the effect of short-range correlations on the longitudinal polarization function P Z (t). There is a priori no reason for limiting the description of the field distribution to the kurtosis as done recently in Ref. [4]. Higher terms could matter. Here we introduce a numerical method which does not have this limitation. At the end of this section it is applied to the analysis of the aforementioned Yb 2 Ti 2 O 7 spectrum.
Taking the vector field distribution at the muon sites, D v (B loc ), to be the product of the three Cartesian component field distributions, where D c (B α loc ) is the component field distribution along the α direction, the zero-field static longitudinal polarization function can be written as follows: with 2 and ω µ = γ µ B (γ µ = 2π × 135.54 MHz/T). It is interesting to examine the P stat Z (t) expression. Let us replace B X loc by −B X loc in Eq. 2. The function P stat Z (t) does not change, but now it is expressed in terms of D c (−B X ) rather than D c (B X ). The same analysis applies for the B Y loc and B Z loc variables. Consequently, . It is easy to extend our scope to the case of a longitudinal-field measurement of the relaxation. By definition, an external field B ext is applied along the Z direction. To account for it we only Thus, we implicitly assume that the system is not modified by the external field. Then it must be noted that, contrary to the zero-field case, it is in principle possible to resolve an asymmetric component field distribution.
In fact, there is always some magnetic dynamics. Assuming the strong-collision model to be valid, P Z (t) can be expressed with an integral equation in terms of the fluctuation rate ν c and P stat Z (t): To proceed further, as for the analytical analysis of P Z (t) given in Ref. [4], we shall assume the three Cartesian component field distributions to be identical. Thus we can write P Z (t) = P Z (t, ν c , {D c (B α loc )}). Experimentally, the measured asymmetry is described as the sum of two components: where a s and a bg are the initial asymmetry for the sample and background, respectively. We stress that experimentally there is always some background, i.e. some implanted muons miss the sample. Our experience shows that it should always be taken into account for a reliable analysis. For the numerical analysis we shall consider N B discrete field values B α loc,i in the finite field range [B min , B max ]. To determine D c (B α loc ) we need to minimize the following chi-square: Here N t is the number of discrete equidistant time channels of the measured asymmetry data, d(t i ) the measured asymmetry at time t i , and σ i the related standard deviation. For a reliable fit χ 2 should be roughly equal to the number of degrees of freedom, i.e. χ 2 ≃ N dof = N t − N p , where N p is the number of free parameters. If we were to solely minimize χ 2 the noise for D c (B α loc ) would be substantial. To proceed we shall use the maximum entropy (ME) concept, and therefore we shall minimize where λ is a Lagrange parameter and E the entropy, i.e.
where δB α loc,i is the field distribution step which is taken uniform here. The basic idea is to minimize χ 2 while maximizing E.
As a first example of the use of the technique described above, we analyzed the Yb 2 Ti 2 O 2 spectrum already considered with the analytical method in Ref. [4]. There it was shown that the small applied longitudinal field did not have any effect on the spectrum. So it will be analyzed as a zero-field spectrum. We take B max = −B min = 20 mT and N B = 23 points to describe the field distribution. Since a symmetrized distribution is considered, it counts only to  Figure 5. Comparison of the results from two methods for the determination of the symmetrized component field distribution at the muon site derived from measurements of the longitudinal asymmetry for a powder sample of Yb 2 Ti 2 O 7 [8] below the temperature at which the specific heat displays a sharp peak. While the result from the maximum entropy concept is obtained in this work, the analytical analysis was previously presented [4]. 12 field points besides the N p = 3 parameters a s , a bg , and ν c . Two methods were tested for the F minimization: the reverse Monte Carlo (RMC) and the conventional Gauss-Newton (GN) methods. We found the first method better in terms of speed of convergence. In addition, it has the ability to converge to the global minimum, i.e. local minima are avoided. The RMC method is described in the appendix to this report. In Fig. 5 we display the symmetrized field distributions extracted from the spectrum shown in Fig. 6. It is rewarding that the previous analytical result [4] is consistent with the ME/RMC data. This provides a justification to limit the description of the field distribution to the kurtosis in this case. Obviously, the analytical analysis is performed much faster than the numerical one.

Discussion
The µSR method is recognized for its ability to easily detect magnetic transitions. However, one should not overinterpret the experimental data, such as proposing a long-range magnetic order when the features of the asymmetry does not allow to do it. Inferring a magnetic structure to be incommensurate with only µSR data at hand does not seem reasonable. To further complicate the state of affairs, the observation of the expected initial asymmetry without a spontaneous oscillation is not a proof of the absence of a magnetic order. The first recognized counter example was published for the ordered spin-ice Tb 2 Sn 2 O 7 in Ref. [9].
The µSR technique is well adapted for the study of quasi-static magnetic dynamics, i.e. with dynamics in the micro-and nanosecond time scale. Such dynamics is seen in exotic magnetic materials such as some heavy fermion and frustrated systems. The longitudinal relaxation is usually found to strongly deviate from the conventional exponential relaxation. The number of analytical formulas available for modeling the data is limited. We have introduced a numerical method to extract from the measured relaxation the field distribution at its origin. It supplements a recently proposed analytical method. The numerical method combines the best of two worlds: the maximum entropy concept first introduced by Rainford and Daniell in the µSR community [10] and the reverse Monte Carlo algorithm of common use for analyzing neutron diffraction data; see Ref. [11] and references therein. The numerical method proposed here can probably be extended to the analysis of spontaneous oscillations, at least in the case of only one muon magnetic site. The computation discussed here involves three-dimensional field integrations; see Eq. 2. If B 0 is sufficiently large relative to ∆ Z.m , the field integrations get restricted to one dimension [1]. Then the number N B of discrete field values B α loc,i can be increased relative to the case considered here such that the asymmetry spectrum in a rotating frame is properly described, while keeping the computing time reasonable.
Acknowledgments AY gratefully acknowledges Dr J. Sugiyama for a useful communication.

Appendix A. The reverse Monte Carlo method
We write F = F a 0 , a bg , ν c , D c (B α loc,i ) ≡ F (r), i.e. the set of parameters for F -a 0 , a bg , ν c , D c (B α loc,i ) -is denoted with the single vector r of dimension N r = N p + N B . Each iteration j of the analysis consists of the following steps (see e.g. [11,12]): (1) A component r i of r is randomly chosen (i = 1, 2, ..., N r ).
(2) We change r i by a factor of (1 + Sǫ), where, in our case, ǫ = 0.03 (i.e. 3% of change) with the sign S = 1 or S = −1 randomly chosen. Thus, we obtain the new vector r new .
(3) F at iteration j, that is F j , is computed using Eq. 6 for r = r new . (4) If ∆F j ≡ F j −F j−1 < 0 the iteration is successful and r new is accepted and stored (r new → r). Otherwise r new is rejected and previous r is kept. The next iteration is repeated from step (1). There is an exception at this stage: for ∆F j positive there is still a possibility to accept the iteration, and therefore r new , with the probability P = exp(−∆F j /p). The parameter p is chosen such that about half of the iterations are accepted [11]. We found p = 0.003 to be a proper choice. The p value depends on ∆F j and consequently on ǫ. This mechanism of accepting the iteration with a positive ∆F j turned out to be important for escaping local minima of F (r). The larger the number of independent parameters N r the higher the probability to converge to a local minimum.
This algorithm depends on several parameters: ǫ, p, and λ. The larger ǫ the faster the convergence at the early stage of iterations. However, one should consider that ǫ determines also the accuracy of the fitted parameters. As a variation of this algorithm one can gradually reduce ǫ with increasing iteration number j. The parameter ǫ is related to ∆F j and consequently to the parameter p, determining the probability to accept an iteration with a positive ∆F j . Obviously the larger ǫ the larger the absolute value of ∆F j . In our implementation of the algorithm we analyze the average number of accepted and rejected iterations within the last 10 iterations and we increase or reduce parameter p to reach the 50% goal for accepted iterations. The value of λ determines the weight of entropy in the function F (r). Entropy determines the probability of a distribution, while χ 2 determines the fit quality. χ 2 is a statistical quantity with an expectation value of N dof and standard deviation √ N dof . Consequently, for a good fit Note that the number of free parameters used for the evaluation of N dof is N p i.e. the parameters entering the entropy calculation are excluded. If χ 2 is too small one should increase λ and conversely for a large χ 2 , λ should be reduced. The result does not change significantly with a reasonable variation of λ for χ 2 ≃ N dof . The particular value of λ must be adapted according to the number of field points N B . To ensure the reliability of the result it is important to check that the condition χ 2 ≃ N dof is satisfied. The value of N B can be arbitrarily chosen, however one should consider that the computation time increases as N 3 B . A minimum number of points in the distribution is obviously required.
In the case of our spectrum, a stable solution was reached after about 500 successful RMC iterations. Its stability was checked by extending the computation up to 2000 successful RMC iterations. The computation was started with an unbiased distribution. i.e. assuming a flat distribution. Note that despite the large number of iterations required to converge, F (r) is evaluated only once per iteration. Thus the convergence rate is quite fast considering the large number of parameters r i .
Obviously, other than RMC methods may be used for a nonlinear minimization such as that of F . Note that the linear optimization method suggested for transverse field µSR analysis in Ref. [10] does not apply in present case. The relation between time and field domains given by Eq. 2 is highly nonlinear and consequently iterative nonlinear minimization algorithms should be used. The obvious merits of RMC are its simplicity, a reasonable computation efficiency, and most importantly, possibility to avoid local χ 2 minima in parameter space due to the mechanism of accepting also the iterations with reduced fit quality. The presented method can be very useful in analysis of µSR data of correlated and magnetically ordered systems in zero or longitudinal field.