Taking advantage of multiplet structure for lineshape analysis in Fourier space

Lineshape analysis is a recurrent and often computationally intensive task in optics, even more so for multiple peaks in the presence of noise. We demonstrate an algorithm which takes advantage of peak multiplicity (N) to retrieve line shape information. The method is exemplified via analysis of Lorentzian and Gaussian contributions to individual lineshapes for a practical spectroscopic measurement and benefits from a linear increase in sensitivity with the number N. The robustness of the method and its benefits in terms of noise reduction and order of magnitude improvement in run-time performance are discussed.


I. INTRODUCTION
Analysis of signal lineshapes is a prominent problem and a theme of importance in physics, chemistry and biomedicine. Ranging from spectroscopy [1][2][3] to scattering techniques [4][5][6][7], the lineshape can reveal underlying physical processes. For example, relaxation dynamics very commonly give rise to exponential decays in time which correspond in spectroscopy to Lorentzian peaks in frequency, while static disorder and instrumental effects typically induce Gaussian distributions of characteristic frequencies. Together, the two phenomena yield signal lineshapes which are convolutions of Lorentzians with Gaussians, objects referred to as Voigt lineshapes. The accurate parametrization of a Voigt lineshape retrieves the Lorentzian contribution which e.g. quantifies correlation lengths in X-ray and neutron scattering [7] and coherence times of quantum systems in frequency-dependent spectroscopy [8][9][10]; and the Gaussian part which is due to extrinsic factors such as grain size distribution in x-ray diffraction [6,11] and to spatial inhomogeneity in optical spectroscopy [2,12]. Approaches to lineshape determination deal with efficient numeric approximations [13,14] or tackle the deconvolution of single peaks in Fourier space [15]. However, the calculation and fitting of a Voigt profile remains computationally intensive. In case of finite-periodic signals, frequently encountered e.g. for electronic multiplets of ions in solids [16] or rotational spectra of molecules [17][18][19], the problem is exacerbated by the need to sum over Voigt profiles. In this article we extend previous work on single Voigt analysis (see [15] for an overview) and exploit the regular spacing for a robust and direct determination of the individual lineshape in Fourier space. We fit an exponential to the envelope of the Fourier transform of the signal, which is parametrized only by the Lorentzian and Gaussian contribution to the individual line width. The method quantifies these contributions without the need for involved multi-Voigt profile analysis. Results are also readily checked by direct visual inspection of the Fourier transform. Furthermore, our method is computationally less expensive and it's sensitivity increases linearly in N , the number of profiles. Although we focus on the common case of Voigt-shaped profiles, we expect our method will prove advantageous for other line-profiles, such as listed in refrence [20], as well. We first derive the mathematical foundations of the method in section II. In section III we apply the procedure to typical experimental data from solid-state spectroscopy, in this case a rare-earth doped crystal. Finally, in section IV we discuss the method's performance when the initial assumption of regular spacing is relaxed.

II. DERIVATION OF THE METHOD
We consider a real, finite-periodic signal S(x) = i S i (x), consisting of N profiles S i (x) spaced by ∆x. We focus on the common cases of individual signal profiles S i (x) of Lorentzian, Gaussian and Voigt shapes. The general case for other classes of lineshapes is addressed subsequently.

A. Dirac delta model
For simplicity we first model S(x) as a series of N Dirac delta functions δ(x), with N even, symmetrically disposed around x = 0. Definitions and a detailed derivation are in section V A and V B of the Appendix. The Fourier transform F : x → k applied to the Dirac model S D (x) is given by where A represents the integrated area of an individual signal profile. Figure 1 shows the amplitude spectrum |F D | , which exhibits a long periodicity T s = 2π/∆x modulated by a short periodicity T f = T s /N . Note the frequency doubling of |F D | w.r.t F D . The occurrence of T s in |F D | stems from the regular spacing of the S i (x), and T f from the length of the multiplet signal, which can be understood as the product of a rectangular window function with a Dirac comb (details in Appendix V C) in analogy to the result of an N −periodic diffraction grating [21]. We denote the i-th maximum of |F D | as M i , i.e. |F D (i × T s )| = M i ∀i ∈ N 0 . The inferior local maxima are labeled by L j . From equation 1 we obtain M i ∝ N for ∆xk/2 mod π = 0 with the rule of Bernoulli-de l'Hôspital and L j,min = 1/ √ 2π does not depend on the multiplet number N . M i ∝ N represents the scaling of the sensitivity of our Fourier lineshape analysis (FLA) method.

B. Lorentzian and Gaussian lineshapes
We continue for Lorentzian-shaped signal lines S i (x) = L(x), as illustrated in figure 2(a) and 2(b) for N = 8. The shape is given by with full-width-at-half-maximum (FWHM) Γ, center peak coordinate x 0 and integrated area A = ∞ −∞ L(x)dx. An N -fold repeated signal with spacing ∆x is represented by and its Fourier transform F L = F[S L (x)] as for x 0 = 0. We see that F L (eq. 4) differs from F D (eq. 1) but by an exponentially decaying prefactor but the amplitudes M i for i > 0 decay with E L (k) = exp(− 1 2 Γk) and the sensitivity M i ∝ N of F D still holds. In the same way we obtain for a Gaussian model S G (x) which again exhibits the same periodicity T s as equation 1 and equation 4. The exponential prefactor defining the envelope E G (k) of F G is now given by E G (k) = exp(− 1 2 σ 2 k 2 ). The amplitude reduction of M i for i > 0 thus depends on σ which is related to the Gaussian FWHM C. Voigt profile model The Voigt lineshape V (x) is defined as the convolution ( * ) of a Lorentzian with a Gaussian function: and equivalently We find for and observe that the sensitivity M i ∝ N still holds. Due to the convolution theorem the envelope E V (k) defining the amplitudes of the M i is given by which is only parametrized by the Lorentzian line width Γ and the Gaussian σ. For illustration, figure 3 shows this for N = 8, in addition to the Lorentzian and Gaussian cases. The Voigt FWHM f V can be approximated with an accuracy of 0.02% [22] by Given a finite-periodic signal S(x) consisting of N Voigt-profile peaks S i (x), spaced by ∆x, the line width contributions Γ, σ of the individual signal peaks S i (x) can thus be determined by fitting E V (k) to the M i , allowing to distinguish the Lorentzian contribution to the line width from that due to the Gaussian. From equation 10 and figure 3 the behavior of F V is dominated by E L (k) for small k and by E G (k) for large k, which allows for a quick qualitative analysis of the Lorentzian contribution by examination of the first few M i -even by eye. In principle this applies also to the strong suppression or the slow decay of the tail by E G (k) and E L (k) respectively, but might not be unambiguously determined due to noise.

D. Sensitivity and generalization of the method
The Fourier transform simplifies the convolution integral to a product and thus, the number of profiles can be understood as the number of terms of a discrete Fourier transform leading to an increase of the sensitivity, linearly in N (see equation 21 in the Appendix). This point of view can be taken in analogy to the N −periodic diffraction grating the M i become more expressed and sharpen as N increases with the limit of becoming the Dirac comb for N → ∞. This allows us to extract values for M i even if L j is below the noise level in the experimental setup. Thus, our result is particularly interesting for experimental applications where N is large, although the case for small N or even N = 1 (c.f. reference [15]) is possible, but not recommended. Thanks to the simplification of the convolution integral and the linear increase in sensitivity with N , our method applied to other classes of line-profiles [20] which are convolutions of two or more functions, should be as beneficial as for the Voigt profile. However, to which accuracy would have to be analyzed for each case separately.

III. EXAMPLE
In this section we demonstrate the application of the proposed method to a multiplet signal, which in this case is a typical absorbance spectrum measured in transmission for a LiY 1−x Ho x F 4 single crystal with x = 0.25% at a temperature T of 3.8 K with a Fourier transform infrared (FTIR) spectrometer. The Bruker IFS125 spectrometer is based at the infrared beamline of the Swiss Light Source at Paul Scherrer Institut in Villigen, Switzerland [23]. The spectrum in figure 4(a) shows the hyperfine-split ground state to second excited state transition at ∼ 23.3 cm −1 (700 GHz) in spectroscopic units [cm −1 ]. Thanks to the ultra-high resolution FTIR in combination with the highly collimated, high-brillance infrared beam we achieved 0.001 cm −1 (30 MHz) resolution. We refer to Matmon et al. [16] for details on the general experimental setup and FTIR technique as well as for more extensive, spectroscopic work on LiY 1−x Ho x F 4 .

A. Procedure
Equation 10 leads to an algorithmic procedure for extracting the shape-and line width of a finite-periodic signal: 1. Calculate the Fourier transform of the finite-periodic signal and take the modulus.
2. Determine the maxima M i defining the envelope E V (k), which occur with periodicity T s .

Fit the envelope function
to the M i .
We replace the analytical prefactor |F D | of the envelope with P as the global scaling factor, combining factors from the normalization of the Fourier transform, the single signal line area A and experimental sensitivities. The offset constant c accounts for nonzero white noise in the experiment. We note that if ∆x is known, the M i might be more precisely determined thanks to T s = 2π/∆x. This simple procedure allows for algorithmic implementation.

B. Application to experimental data
We apply the above procedure to the transmission spectrum in figure 4(a). More details on the numerical discrete Fourier transform (DFT) procedure are found in the Appendix V E. Figure 4(b) shows the modulus of the DFT coefficients c k of the unapodized (blue) and apodized (red) spectrum (details of apodization are in Appendix V E) in figure 4(a) for k < 600. The characteristic pattern of equation 10 is evident in figure 4. The inset displays all Fourier coefficients. We observe reduced noise frequencies in k-space in figure 4(b) for the apodized spectrum and a different global scaling factor P with respect to the unapodized spectrum. However, the decay constants Γ, σ are unchanged. Figure 5 shows the DFT of the apodized spectrum and the M i . We apply a simple algorithm detecting the evenly spaced maxima M i with period T s ∼ 45 Fourier coefficients. We set a conservative threshold at k = 520 to ensure a fiducial selection of the M i . We drop M k,0 as it is strongly affected by the Fourier transforms zero center peak. The result of the envelope E(Γ, σ, P ) fit for a Lorentzian, Gaussian and Voigt lineshape to the M i is displayed in figure 5. We observe that the Voigt profile yields the best results with Γ = 1.1 × 10 −2 ± 0.12 × 10 −3 cm −1 and σ = 3.5 × 10 −3 ± 0.6 × 10 −3 cm −1 . With equation 7 we find f G = 8.5 × 10 −3 ± 1.4 × 10 −3 cm −1 and with equation 12 f V = 15.7 × 10 −3 ± 1.6 × 10 −3 cm −1 . The same values within the fit precision are obtain by a direct Voigt profile fit to the individual signal peaks. The method determines that f L , the homogeneous part of the line width, is the slightly larger contribution.

IV. ROBUSTNESS, NOISE STABILITY AND RUN TIME
We discuss the robustness of the method under deviations to the initial assumption of regular spacing and compare its performance to a direct least-squares fit method (DFM).

A. Robustness
We examine the sensitivity of the fitting method to deviations from our initial assumption of finite-periodicity. In the case of unequal areas A i , the individual area A i of the signal lines appear as an overall scaling factor which is already considered in section III. This allows the application of the method in cases where the individual signal line varies strongly in intensity e.g. for rotational/vibrational spectra [19]. The extracted line width in case of differing individual line widths Γ i , σ i is a weighted average. As regular spacing is the underlying assumption, a variation of the line spacing ∆x directly modifies the M i and thus the observed envelope E(Γ, σ, P ). Because of the line spacing variation the N frequencies exhibit phase differences resulting in destructive interference and less expressed M i (c.f. Appendix V B). This leads to an overestimation of the width and deviation from actual profile shape. We note that a second order hyperfine interaction in the example spectrum (for details see [16]) results in a 2.7 % deviation from equidistancy of the eight hyperfine lines. Nevertheless we obtain equivalent results for line width and shape, showing the robustness of the method against small (on the scale of ∆x/N ) relaxation in finite-periodicity of 3%. Furthermore, if the variations ∆x i are known, then E(Γ, σ, P ) could be adjusted to depend on the ∆x i .

B. Noise-stability
The noise spectrum is explicitly revealed with the FLA method. In case of a few isolated frequencies, as in e.g. the inset in figure 4, the noise does not affect the envelope fit in k-space and the method naturally implements a low-pass frequency filter. A particular strength of the method is that if there are a few frequencies between the M i where there are what appear to be noise peaks, as is the case in figure 4 at e.g. k = 70. Such noise peaks will affect the DFM results and low-pass filters can not be applied. Consideration only of the primary peaks at M i provides an impartial filter of the associated noise Assuming Gaussian RMS noise for artificial spectra with up to 50% of the signal amplitude, we find the FLA method to be comparable to DFM within the statistical parameter errors (c.f. Appendix V F). This holds true for different ratios of Γ/σ ∈ 0.1, 1, 10. We find that the precision for both, FLA and DFM strongly decreases for Γ/σ < 0.1 and Γ/σ > 10. We observe a systematic error for the FLA method in the form of an increase of the Gaussian line width contribution with increasing RMS noise amplitude. We attribute this to the fact that our analysis in Fourier space reflects the noise spectrum convolved with the signal and thus adds to the Gaussian line width contribution. Although the DFM method may occasionally exhibit a more precise performance, the FLA method extracts line width and shape information algorithmically, with less computational power and fewer data points, as detailed in the next subsection.

C. Run time
We note the significantly different computational intensity of both methods. On our system (Intel(R) Xeon(R) E5-2670 0 @2.6 GHz, 2 processors with 192 GB usable RAM) DFM is at least one to two orders of magnitude slower than the FLA method. A test on 1000 data sets with N = 8 peaks, 1000 data points, and varying Gaussian RMS noise yielded 1.6 s run time for the FLA method and 2543 s for DFM. Runtimes are dataset specific and further (algorithm) optimization may reduce the run time discrepancy. For fundamental reasons, however, a run time difference will persist as the FLA method allows to reduce the parameter space to four (Γ, σ, P, c), whereas a direct fit method has to consider at least six (Γ, σ,area A, position x 0 , spacing s, offset c) under the assumption of identical areas A i , which is rarely fulfilled. Any individual consideration of a parameter in the direct fit scales the parameter space with N . This increases the run time for the DFM method significantly. Furthermore, the final envelope fit is performed on a fraction of the initial data points. In case of insufficient precision of the FLA method, it may well serve as a fast pre-analysis for good DFM start parameters.

V. CONCLUSION
We have introduced and demonstrated a method of extracting the line-width and shape of a multiplet signal. The basic idea and mathematics are actually straightforward generalizations of the Debye-Waller effect [24] in scattering techniques, where the intensities of Bragg peaks shrink with increasing order in proportion to an envelope function decaying as a Gaussian whose width in reciprocal space is inversely proportional to the uncertainty in the positions of the atoms in real space. We offer a direct and computationally efficient way of quantifying the Lorentzian and Gaussian parts of multiple Voigt profiles by harnessing the regular spacing of the signal. The method benefits from a linear increase in sensitivity with the number of profiles of the multiplet. Furthermore, we highlighted that the FLA method is comparable to a direct least square fit, but is less computationally intensive and exhibits significant advantages in the presence of low-frequency noise. Finally, the FLA method might serve as a fast and algorithmic preanalysis for starting parameter determination cases where the assumption of periodic repetition of the same lineshape is broken. Potential applications are ample. X-ray scattering may benefit where quasiperiodic structures are sampled within finite real space windows for problems such as integrated circuit microscopy [25,26]. Our method could also be applied to optical comb-based spectroscopies [27] for precision measurements in atomic and solid-state physics.

A. Definitions
The Fourier transform we use throughout the article is defined as We use the following definition of a Dirac comb X: where δ(x) denotes the Dirac delta function. The Fourier series of the Dirac comb is given by which directly yields another interpretation of equation 1: It is the sum of the first N terms of the Fourier series of the Dirac comb X ∆x (x) for odd N . The rectangular window function of height h and width w is defined as where we use the standard definition for the unit rectangular function rect(x) = 1 ∀x ∈ {−1/2, 1/2} and 0 otherwise. The model where the individual signal peaks S i (x) = G(x) are Gaussian (c.f. figure 2(a) and (b) ) is given by with standard deviation σ, center signal line position x 0 , spacing ∆x and integrated area A. The Fourier transform of S G (x) (with x 0 = 0) is given by

B. Detailed derivation
The Dirac model S D (x) is given by assuming symmetry around x = 0. With {A, n, ∆x, k} ∈ R we obtain for F D which, for N even yields and for odd N For even N we further use and we leave the case for the odd N to the reader. From equation 22 we conclude that Equation 25 shows that the number of signal lines N of S(x) determines the amplitude of the M i = N A/ √ 2π ∀i. This is determined with the rule of Bernoulli-de l'Hôspital where numerator [sin(N ∆xk/2)] and denominator [sin(∆xk/2)] vanish which is the case for ∆xk/2 mod π = 0. Furthermore, L j ∼ A/ √ 2π ∀j as L j,min denotes the point where numerator and denominator are equal to one. Note that for even N this condition is never exactly met and L j,min slightly smaller than 1 but still independent of N . These results lead to which is the scaling of the sensitivity of the FLA method with the number N of individual signal profiles S i (x), in analogy to a finite diffraction grating. This result can be alternatively understood as the approximate ratio of the main maximum to the first side lobe of a sinc-function, as shown in the next subsection.

C. Alternative derivation
A sum of N delta functions centered around 0 and spaced by ∆x can be written as or as For F[S D (x)] ≡ F D follows with the application of the convolution theorem: which reveals equation 1 to be an infinite sum of sinc functions.

D. Detailed derivation of Voigt model
We present the detailed derivation for F V .
which for x 0 = 0 simplifies to We leave the case of odd N to the reader.

F. RMS noise test
For the RMS noise test we create artificial data, where we add different amplitudes of RMS noise and apply the FLA as well as the DFM method with a Voigt lineshape model, for the case of N = 8 peaks. The results are shown in figure 6. Gridlines denote the initial values before RMS noise application. For contribution ratios Γ/σ = 0.1 the FLA method tends to overestimate the overall line width. We attribute this to a systematic error, which pushes the ratio in the limit where the methods precision strongly decreases.