Bimodal 1H Double Quantum Build-Up Curves by Fourier and Laplace-like Transforms on Aged Cross-Linked Natural Rubber

The 1H DQ Fourier and Laplace-like spectra for a series of cross-linked natural rubber (NR) samples naturally aged during six years are presented and characterized. The DQ build-up curves of these samples present two peaks which cannot be described by classical functions. The DQ Fourier spectra can be obtained after a numeric procedure which introduces a correction time which depends less on the chosen approximation, spin-½ and isolated CH2 and CH3 functional groups. The DQ Fourier spectra are well described by the distributions of the residual dipolar coupling correlated with the distribution of the end-to-end vector of the polymer network, and with the second and fourth van Vleck moments. The deconvolution of DQ Fourier spectra with a sum of four Gaussian variates show that the center and the width of Gaussian functions increase linearly with the increase in the cross-link density. The Laplace-like spectra for the natural aged NR DQ build-up curves are presented. The centers of four Gaussian distributions obtained via both methods are consistent. The differences between the Fourier and Laplace-like spectra consist mainly of the spectral resolution in the favor of Laplace-like spectra. The last one was used to discuss the effect of natural aging for cross-linked NR.


Introduction
Homonuclear and heteronuclear residual dipolar couplings (RDCs) or quadrupolar interactions in soft solids such as elastomers and biological tissues represent an important source of information about the structure and molecular dynamics [1][2][3][4][5]. Proton residual dipolar couplings in elastomers reflect changes in the cross-link density, temperature, the uniaxial and biaxial extension or compression as well as the presence of fillers and penetrant molecules. Structure-dynamics-function relationships using RDCs can be investigated for the broad class of elastomer materials [2][3][4].
The measurements of RDCs can be carried out using one-dimensional (1D) and twodimensional (2D) NMR methods. In the 1D case, methods were used such as the dipolar correlation effect in combination with Hahn and solid echo [6,7], the stimulated echo [8], the magic echo [9] and magnetization exchange [10,11]. Model free access to RDCs is given by the analysis of multiple-quantum (MQ) build-up [12][13][14][15][16][17][18] and decay [19] curves recorded in the initial regime of the excitation/reconversion periods of the experiment, as well as the accordion magic sandwich technique [20]. Chemically site-selective RDCs can be elucidated by 2D NMR spectroscopy using, for instance, 13 C-1 H heteronuclear residual dipolar couplings and encoded spinning sideband patterns [21]. The NOESY under magic angle sample spinning (MAS) [22] and double-quantum (DQ) MAS NMR spectroscopy are also used [23,24].
In general, residual dipolar couplings are characterized by heterogeneous distributions. In many cases, the local variation in cross-link density and the presence of network Table 1. The values for the correction time τ c obtained as a result of correction procedure of the DQ Fourier spectra based on spin-1 /2 and isolated CH 2 and CH 3 approximations for the series of cross-link density from NR1 to NR7. The second M 2 and fourth M 4 residual van Vleck moments and the maximum of the DQ Fourier spectra obtained for the last approximation.

NMR Measurements
The measurements, for natural aged samples, were performed using a BRUKER MINISPEC mq20 spectrometer working at 19.7 MHz. The sample temperature during all measurement was 35 • C. The DQ five-pulse sequence used to record the build-up curves is presented in Figure S1 (see Supplementary Information, where the efficiency in creating DQ coherences is also discussed). The tipping pulse length was 8-8.5 µs, dependent on the cross-link density. The excitation/evolution period is denoted by τ and was increased in equal steps up to 4 ms. The evolution period t 0 and z-filter t z were kept short of the order of 20 and 50 µs, respectively. The recycle delay was 0.5 s. The double-quantum filtered signals were normalized to the integral intensity of the single-quantum signal measured in the same conditions as DQ build-up curves.

Proton DQ Build-Up Curves
The DQ build-up curve was recorded from the unaged NR1 sample after the vulcanization (see Table 1) and is shown in Figure 1a. A single component can be observed with a maximum around 700 µs. The maximum DQ signal is obtained from the combined effect of an increase in the intensity of the DQ coherences and decay of the single-quantum coherences due to transverse relaxation during the excitation and reconversion periods of the pulse sequence.  1 H DQ build-up curves measured after six years after production as function of cross-link density (samples NR1 to NR7 from Table 1).
The second component is clearly observed for sample NR1, characterized by the smallest cross-link density, and is almost inexistent for sample NR7 (Table 1), characterized by the largest cross-link density (see Figure 1b). The effect of natural aging, revealed as multi-component 1 H DQ build-up curves, decreases with the increase in the cross-link density.

The Multi-Spin van Vleck Moments Approximation
In many of the previous works, the DQ build-up curves were analyzed in the initial excitation/reconversion time regime where the measurements of the 1 H residual dipolar couplings encoded the second van Vleck moments [16,19,30]. Since one cannot deduce an exact analytical expression to approximate the full DQ build-up and decay curve for a multi-spin system, the normalized DQ signal can be phenomenologically described in terms of the residual DQ second van Vleck moments and effective relaxation times * , as [30] 1 H DQ build-up curves measured after six years after production as function of cross-link density (samples NR1 to NR7 from Table 1).
The effect of one year of aging in natural conditions can be observed by the apparition of a new component in the DQ build-up curve as a shoulder shifted to larger time values compared with the maximum of unaged sample (Figure 1a). This new maximum can be associated with the apparition of polymer chain segments characterized by smaller residual dipolar couplings (RDCs) due to chain scissions by oxidative processes. The data obtained for small and large excitation/reconversion times are almost overlapping over the curve registered with one year before, but a drop is observed in the region of maximum that reflects the presence of two polymer networks with different RDCs. The additional five years of natural aging produces a large effect on the DQ curve. Two components are clearly observed, and the aging effect can be quantified as a displacement of a second component to longer excitation/reconversion times. This displacement can be associated with the increase in the transverse relaxation time and/or decrease in the residual dipolar interactions.
The second component is clearly observed for sample NR1, characterized by the smallest cross-link density, and is almost inexistent for sample NR7 (Table 1), characterized by the largest cross-link density (see Figure 1b). The effect of natural aging, revealed as multicomponent 1 H DQ build-up curves, decreases with the increase in the cross-link density.

The Multi-Spin van Vleck Moments Approximation
In many of the previous works, the DQ build-up curves were analyzed in the initial excitation/reconversion time regime where the measurements of the 1 H residual dipolar couplings encoded the second van Vleck moments [16,19,30]. Since one cannot deduce an exact analytical expression to approximate the full DQ build-up and decay curve for a multi-spin system, the normalized DQ signal can be phenomenologically described in and effective relaxation times T * 2 , as [30] S DQ 2τ, M DQ 2 , T * 2 = where τ is the excitation/reconversion time periods and where the residual second van Vleck moment M DQ 2 is due to the spin-pair nature of the dipolar coupling Hamiltonian and to the fact that the DQ filter edits these pairs.
In Equation (1), the transverse relaxation process during free evolution is accounted by the exponential decay with an effective transverse relaxation time T * 2 . The Equation (1) is valid in the limits of τ T * 2 and M DQ 2 1 2 τ 1. In order to describe bimodal DQ curves, as obtained for the aged NR1 sample (Figure 1a), the theoretical expression could be a sum of two ad hoc functions as given by Equation (1) S DQ (2τ) = A 1 where the M DQ 2,1 , M DQ 2,2 , T * 2,1 and T * 2,2 are the residual DQ second order van Vleck moments and the effective transverse relaxation times of polymer chain segments characterized by different mobility. The A 1 and A 2 are two normalization constants where the quantity A i /(A 1 + A 2 ), (i = 1,2) indicates the proportion of two dynamic components. Figure 2a shows the attempts to fit the experimental data using Equation (2) for NR1 sample six years aged. The failure is due to the fact that Equation (1) is not able to describe the experimental DQ build-up curves for the full excitation/excitation time regime. The best fit of the DQ build-up curve for the aged NR1 sample can approximate only the initial regime and fails dramatically for excitation/reconversion times larger than 0.3 ms (see Figure 3a). Moreover, the function described by Equation (2) as the superposition of two individual components, presented with dashed and short dashed lines in Figure 2a, is not showing a distinct double peak feature. As one can observe from this figure the smallest (left) peak can be approximated but for the large component (right) the theoretical peak (dashed line) is much broad. In fact, the fit of the first component can be conducted relatively well, as one can see from Figure 2b. For this fit, we consider only the experimental data up to~0.85 ms of the beginning of the second built-up curve. Even in the absence of the small peak, if the behavior is described by Equation (1) where τ is the excitation/reconversion time periods and where the residual second van Vleck moment is due to the spin-pair nature of the dipolar coupling Hamiltonian and to the fact that the DQ filter edits these pairs.
In Equation (1), the transverse relaxation process during free evolution is accounted by the exponential decay with an effective transverse relaxation time * . The Equation (1) is valid in the limits of ≪ * and ≪ 1.
In order to describe bimodal DQ curves, as obtained for the aged NR1 sample (Figure 1a), the theoretical expression could be a sum of two ad hoc functions as given by Equation (1 where the , , , , , * and , * are the residual DQ second order van Vleck moments and the effective transverse relaxation times of polymer chain segments characterized by different mobility. The A1 and A2 are two normalization constants where the quantity Ai/(A1 + A2), (i = 1,2) indicates the proportion of two dynamic components. Figure 2a shows the attempts to fit the experimental data using Equation (2) for NR1 sample six years aged. The failure is due to the fact that Equation (1) is not able to describe the experimental DQ build-up curves for the full excitation/excitation time regime. The best fit of the DQ build-up curve for the aged NR1 sample can approximate only the initial regime and fails dramatically for excitation/reconversion times larger than 0.3 ms (see Figure 3a). Moreover, the function described by Equation (2) as the superposition of two individual components, presented with dashed and short dashed lines in Figure 2a, is not showing a distinct double peak feature. As one can observe from this figure the smallest (left) peak can be approximated but for the large component (right) the theoretical peak (dashed line) is much broad. In fact, the fit of the first component can be conducted relatively well, as one can see from Figure 2b. For this fit, we consider only the experimental data up to ~0.85 ms of the beginning of the second built-up curve. Even in the absence of the small peak, if the behavior is described by Equation (1), the right peak cannot be approximated by any values of and * (Figure 2c). The best fit of the first peak (b) and the second peak (c) using Equation (2). Figure 2. (a) The best fit of DQ build-up curve for sample NR1 aged six years using Equation (4). The best fit of the first peak (b) and the second peak (c) using Equation (2).
Polymers 2021, 13, x FOR PEER REVIEW 6 of 27 Figure 2. (a) The best fit of DQ build-up curve for sample NR1 aged six years using Equation (4). The best fit of the first peak (b) and the second peak (c) using Equation (2).  The star symbol showed the error in the transformed spectra due to the correction procedure.

The Spin-½ Pair Approximation
For a spin-½ pair all the terms of the Hamiltonian described by Equation (S1) in supplementary information) commute with each other, which allows obtaining an exact time evaluation of the spin system response to the DQ five pulse sequence that is given by and τ is the excitation/reconversion time period and * is the effective transverse relaxation time of the single-quantum coherences. For simplicity, we can assume that the relaxation processes of the DQ coherences are characterized by the same value. In this section we will demonstrate that, due to the particular nature of the problem, this assumption will not significantly change the final result. The ( )  represents the statistical average over the end-to-end vector, and angle β between the direction of the pre-averaged end-to-end vector and external static magnetic field 0 B  . Figure 3a presents a simulation of a DQ signal (continuous line) for an ideal polymer characterized by a single value of the residual dipolar constant ̅ = 2 ⁄ = 2 kHz and an envelope of NMR signal (dashed line) characterized by the transverse relaxation time * = 2 ms. We will use this simulated curve to demonstrate the simplest procedure based on the Fourier transform, (hereafter we use the notation FT {…}) to obtain the distribution of the residual dipolar couplings. The negative Fourier transform, −FT {…}, of the DQ filtered signal is presented in Figure 3b. The spectrum can be easily interpreted if we will rewrite the Equation (3) The negative spectra obtained from Equation (4) can be written as a sum of two terms with the same weight:  (6) with a single residual dipolar coupling ω D /2π = 2 kHz and an effective relaxation time T * 2 = 2 ms; (b) the Fourier spectrum of the DQ signal shown in (a) and (c) the corrected DQ Fourier spectra from (b). The star symbol showed the error in the transformed spectra due to the correction procedure.

The Spin-1 /2 Pair Approximation
For a spin-1 /2 pair all the terms of the Hamiltonian described by Equation (S1) in supplementary information) commute with each other, which allows obtaining an exact time evaluation of the spin system response to the DQ five pulse sequence that is given by where ω D = √ 3/2ω D and τ is the excitation/reconversion time period and T * 2 is the effective transverse relaxation time of the single-quantum coherences. For simplicity, we can assume that the relaxation processes of the DQ coherences are characterized by the same value. In this section we will demonstrate that, due to the particular nature of the problem, this assumption will not significantly change the final result. The (. . .) represents the statistical average over the end-to-end vector, → R and angle β between the direction of the pre-averaged end-to-end vector and external static magnetic field → B 0 . Figure 3a presents a simulation of a DQ signal (continuous line) for an ideal polymer characterized by a single value of the residual dipolar constant ν D = ω D /2π = 2 kHz and an envelope of NMR signal (dashed line) characterized by the transverse relaxation time T * 2 = 2 ms. We will use this simulated curve to demonstrate the simplest procedure based on the Fourier transform, (hereafter we use the notation FT { . . . }) to obtain the distribution of the residual dipolar couplings. The negative Fourier transform, −FT { . . . }, of the DQ filtered signal is presented in Figure 3b. The spectrum can be easily interpreted if we will rewrite the Equation (3) The negative spectra obtained from Equation (4) can be written as a sum of two terms with the same weight: i.e., a negative Lorentzian peak centered in origin (described by the first term) and a second positive Lorentzian peak centered at double value of residual dipolar coupling constant given by the last term in Equation (5). One can observe that, regardless of the value 2 ⁄ , the first term is always negative centered in origin while the desired spectral amplitudes are disperse over all 2 ⁄ values. A simple numerical procedure can be designed to cancel the negative peak independent on the residual dipolar interactions. Due to the long wings of Lorentzian function the proper action is to be applied in the time i.e., a negative Lorentzian peak centered in origin (described by the first second positive Lorentzian peak centered at double value of residual dipo constant given by the last term in Equation (5). One can observe that, rega value 2 ⁄ , the first term is always negative centered in origin while the tral amplitudes are disperse over all 2 ⁄ values. A simple numerical proc designed to cancel the negative peak independent on the residual dipolar Due to the long wings of Lorentzian function the proper action is to be applie i.e., a negative Lorentzian peak centered in orig second positive Lorentzian peak centered at dou constant given by the last term in Equation (5). O value 2 ⁄ , the first term is always negative ce tral amplitudes are disperse over all 2 ⁄ value designed to cancel the negative peak independe Due to the long wings of Lorentzian function the p i.e., a negative Lorentzian peak centered in origin (described by the first term) and a second positive Lorentzian peak centered at double value of residual dipolar coupling constant given by the last term in Equation (5). One can observe that, regardless of the value ω D /2π, the first term is always negative centered in origin while the desired spectral amplitudes are disperse over all ω D /2π values. A simple numerical procedure can be designed to cancel the negative peak independent on the residual dipolar interactions. Due to the long wings of Lorentzian function the proper action is to be applied in the time domain instead on the frequency domain. Therefore, the Fourier transform procedure is applied on a new function defined in the time domain where an exponential decay with a correction time and amplitude 1 2 is added to the negative DQ signal i.e., a negative Lorentzian peak centered in origin (described by the first term) and a second positive Lorentzian peak centered at double value of residual dipolar coupling constant given by the last term in Equation (5). One can observe that, regardless of the value 2 ⁄ , the first term is always negative centered in origin while the desired spectral amplitudes are disperse over all 2 ⁄ values. A simple numerical procedure can be designed to cancel the negative peak independent on the residual dipolar interactions. Due to the long wings of Lorentzian function the proper action is to be applied in the time domain instead on the frequency domain. Therefore, the Fourier transform procedure is applied on a new function defined in the time domain where an exponential decay with a correction time and amplitude ½ is added to the negative DQ signal (6) where is the correction time. This parameter has to be obtained from the best fit. The negative sign used before has the role to produce the final result into a form that is easily identifiable, by means of a positive Fourier spectral distribution of the residual dipolar constants.
A dedicated program was written in C++ to implement the steps presented in Equation (6). First, the correction time is chosen by the well-known secant method which systematically splits in half an interval described by two extreme values. For each particular , a single point where ε is a small positive constant which represents the chosen fit error. The final interval quantifies also the error in the determination of the correction time value . The corrected spectrum of a doubled residual dipolar constant is presented in Figure 3c. The efficiency of the algorithm is demonstrated by the RDC distribution spectrum with a small residual contribution of around zero value marked by a star in Figure 3c. The Fourier transform with a correction time procedure was applied on the DQ build-up data normalized at the SQ amplitude for the entire NR sample series (see Figure 1b) and the corresponding normalized distributions of the residual dipolar couplings were obtained.
The correction times c τ are presented in Table 1.

The DQ Fourier Spectra
By applying the same steps presented for the spin-½ pair in Equations (3)-(6), we can extend the Fourier analysis, also as an approach to the cases of isolated CH3 functional groups. In this case, the spin system response on the DQ five pulse sequences has the same mathematical form as those presented in Equation (3), but the residual dipolar constant 3 CH D ω is specific to the isolated CH3 functional group. From here, one can ex- i.e., a negative Lorentzian peak centered in origin (described by th second positive Lorentzian peak centered at double value of residu constant given by the last term in Equation (5). One can observe th value 2 ⁄ , the first term is always negative centered in origin wh tral amplitudes are disperse over all 2 ⁄ values. A simple numeri designed to cancel the negative peak independent on the residual Due to the long wings of Lorentzian function the proper action is to b domain instead on the frequency domain. Therefore, the Fourier tra applied on a new function defined in the time domain where an expo correction time and amplitude ½ is added to the negative DQ signal is the correction time. This parameter has to be obtained f negative sign used before has the role to produce the final result into identifiable, by means of a positive Fourier spectral distribution of constants.
A dedicated program was written in C++ to implement the steps tion (6). First, the correction time is chosen by the well-known secan tematically splits in half an interval described by two extreme v ticular , a single point where ε is a small positive constant which represents the chosen fi terval quantifies also the error in the determination of the correction corrected spectrum of a doubled residual dipolar constant is present efficiency of the algorithm is demonstrated by the RDC distributio small residual contribution of around zero value marked by a star in rier transform with a correction time procedure was applied on th normalized at the SQ amplitude for the entire NR sample series (see corresponding normalized distributions of the residual dipolar coup The correction times c τ are presented in Table 1.

The DQ Fourier Spectra
By applying the same steps presented for the spin-½ pair in E can extend the Fourier analysis, also as an approach to the cases o tional groups. In this case, the spin system response on the DQ five the same mathematical form as those presented in Equation (3), but i.e., a negative Lorentzian peak centered in origin (described by second positive Lorentzian peak centered at double value of resid constant given by the last term in Equation (5). One can observe t value 2 ⁄ , the first term is always negative centered in origin w tral amplitudes are disperse over all 2 ⁄ values. A simple numer designed to cancel the negative peak independent on the residual Due to the long wings of Lorentzian function the proper action is to domain instead on the frequency domain. Therefore, the Fourier tr applied on a new function defined in the time domain where an exp correction time and amplitude ½ is added to the negative DQ signa is the correction time. This parameter has to be obtained negative sign used before has the role to produce the final result int identifiable, by means of a positive Fourier spectral distribution o constants.
A dedicated program was written in C++ to implement the ste tion (6). First, the correction time is chosen by the well-known seca tematically splits in half an interval described by two extreme ticular , a single point where ε is a small positive constant which represents the chosen terval quantifies also the error in the determination of the correctio corrected spectrum of a doubled residual dipolar constant is presen efficiency of the algorithm is demonstrated by the RDC distribut small residual contribution of around zero value marked by a star i rier transform with a correction time procedure was applied on t normalized at the SQ amplitude for the entire NR sample series (s corresponding normalized distributions of the residual dipolar cou The correction times c τ are presented in Table 1.

The DQ Fourier Spectra
By applying the same steps presented for the spin-½ pair in can extend the Fourier analysis, also as an approach to the cases tional groups. In this case, the spin system response on the DQ five the same mathematical form as those presented in Equation (3), bu where τ c is the correction time. This parameter has to be obtained from the best fit. The negative sign used before has the role to produce the final result into a form that is easily identifiable, by means of a positive Fourier spectral distribution of the residual dipolar constants.
A dedicated program was written in C++ to implement the steps presented in Equation (6). First, the correction time is chosen by the well-known secant method which systematically splits in half an interval described by two extreme τ c values. For each particular τ c , a single point i.e., a negative Lorentzian peak centered in origin (described by the f second positive Lorentzian peak centered at double value of residual d constant given by the last term in Equation (5). One can observe that, r value 2 ⁄ , the first term is always negative centered in origin while t tral amplitudes are disperse over all 2 ⁄ values. A simple numerical p designed to cancel the negative peak independent on the residual dipo Due to the long wings of Lorentzian function the proper action is to be ap domain instead on the frequency domain. Therefore, the Fourier transfo applied on a new function defined in the time domain where an exponen correction time and amplitude ½ is added to the negative DQ signal

FT FT FT
where is the correction time. This parameter has to be obtained from negative sign used before has the role to produce the final result into a fo identifiable, by means of a positive Fourier spectral distribution of the constants.
A dedicated program was written in C++ to implement the steps pre tion (6). First, the correction time is chosen by the well-known secant me tematically splits in half an interval described by two extreme value ticular , a single point where ε is a small positive constant which represents the chosen fit err terval quantifies also the error in the determination of the correction tim corrected spectrum of a doubled residual dipolar constant is presented i efficiency of the algorithm is demonstrated by the RDC distribution s small residual contribution of around zero value marked by a star in Fig  rier transform with a correction time procedure was applied on the DQ normalized at the SQ amplitude for the entire NR sample series (see Fig   {} ω D =0 is considered until the best value is found i.e., a negative Lorentzian peak centered in origin (described by the second positive Lorentzian peak centered at double value of residua constant given by the last term in Equation (5). One can observe that value 2 ⁄ , the first term is always negative centered in origin whil tral amplitudes are disperse over all 2 ⁄ values. A simple numerica designed to cancel the negative peak independent on the residual di Due to the long wings of Lorentzian function the proper action is to be domain instead on the frequency domain. Therefore, the Fourier trans applied on a new function defined in the time domain where an expon correction time and amplitude ½ is added to the negative DQ signal

FT FT FT
where is the correction time. This parameter has to be obtained fro negative sign used before has the role to produce the final result into a identifiable, by means of a positive Fourier spectral distribution of th constants.
A dedicated program was written in C++ to implement the steps tion (6). First, the correction time is chosen by the well-known secant m tematically splits in half an interval described by two extreme val ticular , a single point where ε is a small positive constant which represents the chosen fit terval quantifies also the error in the determination of the correction corrected spectrum of a doubled residual dipolar constant is presente efficiency of the algorithm is demonstrated by the RDC distribution small residual contribution of around zero value marked by a star in F where ε is a small positive constant which represents the chosen fit error. The final interval quantifies also the error in the determination of the correction time value τ c . The corrected spectrum of a doubled residual dipolar constant is presented in Figure 3c. The efficiency of the algorithm is demonstrated by the RDC distribution spectrum with a small residual contribution of around zero value marked by a star in Figure 3c. The Fourier transform with a correction time procedure was applied on the DQ build-up data normalized at the SQ amplitude for the entire NR sample series (see Figure 1b) and the corresponding normalized distributions of the residual dipolar couplings were obtained. The correction times τ c are presented in Table 1.

The DQ Fourier Spectra
By applying the same steps presented for the spin-1 2 pair in Equations (3)- (6), we can extend the Fourier analysis, also as an approach to the cases of isolated CH 3 functional groups. In this case, the spin system response on the DQ five pulse sequences has the same mathematical form as those presented in Equation (3), but the residual dipolar constant ω CH 3 D is specific to the isolated CH 3 functional group. From here, one can extend the approximation to the isolated CH 3 and CH 2 functional groups. The manner of treatment is unitary; the only thing which will be different is the prefactor of the exponential correction term from Equation (6) Polymers 2021, 13, x FOR PEER REVIEW 8 of 27 i.e., a negative Lorentzian peak centered in origin (described by the first term) and a second positive Lorentzian peak centered at double value of residual dipolar coupling constant given by the last term in Equation (5). One can observe that, regardless of the value 2 ⁄ , the first term is always negative centered in origin while the desired spectral amplitudes are disperse over all 2 ⁄ values. A simple numerical procedure can be designed to cancel the negative peak independent on the residual dipolar interactions. Due to the long wings of Lorentzian function the proper action is to be applied in the time domain instead on the frequency domain. Therefore, the Fourier transform procedure is applied on a new function defined in the time domain where an exponential decay with a correction time and amplitude ½ is added to the negative DQ signal i.e., a negative Lorentzian peak centered in origin (described by the firs second positive Lorentzian peak centered at double value of residual dip constant given by the last term in Equation (5). One can observe that, reg value 2 ⁄ , the first term is always negative centered in origin while the tral amplitudes are disperse over all 2 ⁄ values. A simple numerical pro designed to cancel the negative peak independent on the residual dipola Due to the long wings of Lorentzian function the proper action is to be appl domain instead on the frequency domain. Therefore, the Fourier transform applied on a new function defined in the time domain where an exponentia correction time and amplitude ½ is added to the negative DQ signal where A i is the desired distribution function, and N and k are two constants specific to each approximation. Then, the corrected Fourier transform is: i.e., a negative Lorentzian peak centered in origin (described by the first term) and a second positive Lorentzian peak centered at double value of residual dipolar coupling constant given by the last term in Equation (5). One can observe that, regardless of the value 2 ⁄ , the first term is always negative centered in origin while the desired spectral amplitudes are disperse over all 2 ⁄ values. A simple numerical procedure can be designed to cancel the negative peak independent on the residual dipolar interactions.
where k from Equation (9) is related to the N and k constants from Equation (8).  Figure 4b presents the dependence of the τ c on the cross-link density obtained using the Equation (9). The correction time τ c characteristic to the NR series, with the exception of NR1, decays linearly with the increase in the cross-link density. Since all Fourier spectra are broadened due to the relaxation processes, the residual dipolar coupling distributions present two unresolved peaks (see Figure 4a). These two peaks are more evident for the NR1 sample and merge totally for the NR7 sample. With the increase in the cross-link density the distributions of ω D /2π becomes broader. At a simpler inspection, the residual dipolar coupling distributions look similar to those obtained in the case of spin-1 2 approximation, and this is the reason why the Fourier spectra are not presented for spin-1 2 approximation. In fact, in the following section, we will demonstrate the differences between the obtained spectra, considering the previous described approximations are so small that they can be neglected, compared with the differences between the DQ Fourier spectra obtained for different samples.
where Ai is the desired distribution function, and N and k are two constants specific to each approximation. Then, the corrected Fourier transform is: (9) where k from Equation (9) is related to the N and k constants from Equation (8). Figure 4a presents the distributions of residual dipolar coupling as a corrected Fourier transform of normalized DQ build-up curves for the entire cross-linked series of NR samples while the Figure 4b presents the dependence of the on the cross-link density obtained using the Equation (9). The correction time characteristic to the NR series, with the exception of NR1, decays linearly with the increase in the cross-link density. Since all Fourier spectra are broadened due to the relaxation processes, the residual dipolar coupling distributions present two unresolved peaks (see Figure 4a). These two peaks are more evident for the NR1 sample and merge totally for the NR7 sample. With the increase in the cross-link density the distributions of 2 ⁄ becomes broader. At a simpler inspection, the residual dipolar coupling distributions look similar to those obtained in the case of spin-½ approximation, and this is the reason why the Fourier spectra are not presented for spin-½ approximation. In fact, in the following section, we will demonstrate the differences between the obtained spectra, considering the previous described approximations are so small that they can be neglected, compared with the differences between the DQ Fourier spectra obtained for different samples. In order to be visually compared, in Figure 5 the distributions for both the spin-1 /2 pair approximation (continuous black line) and for the isolated CH 3 and CH 2 functional groups approximation (dashed grey line) were represented together for the first and last series samples, which are NR1 and NR7. In order to evaluate the deviation of a curve from another curve, the merit value χ 2 can be defined as where y 1,2 are the amplitude values of arbitrary (1) and (2) distribution curves. Then, it will be interesting to evaluate the error obtained in the case of bad choice of the approximation compared with the differences between two distributions obtained for samples with successive value of the cross-link density. For example, the merit value χ 2 for the NR1 Fourier spectra obtained for spin-1 /2 pair approximation and with isolated CH 3 and CH 2 functional groups approximation is χ 2 NR1 /χ 2 NR1,NR2 ∼ = 5.4% from the merit value χ 2 obtained by comparing the NR1 with NR2 Fourier spectra. This percentage of the merit value increases at χ 2 NR6 /χ 2 NR6,NR7 ∼ = 30.1% for NR6 Fourier spectra compared with the differences between NR6 and NR7 spectra. In conclusion, since the differences between any two samples are much largest than the error due to the bad choice of the model, the Fourier spectra of a series of aged cross-linked natural rubber samples can be well characterized by any model, presented earlier as approximations.
six-years aged cross-linked NR samples. (b) The dependence of the correction effective relaxation time functions of cross-link density. The dashed lines represent the linear fit of data for samples NR2 up to NR7.
In order to be visually compared, in Figure 5 the distributions for both the spin-½ pair approximation (continuous black line) and for the isolated CH3 and CH2 functional groups approximation (dashed grey line) were represented together for the first and last series samples, which are NR1 and NR7. In order to evaluate the deviation of a curve from another curve, the merit value χ 2 can be defined as  (1) and (2) distribution curves. Then, it will be interesting to evaluate the error obtained in the case of bad choice of the approximation compared with the differences between two distributions obtained for samples with successive value of the cross-link density. For example, the merit value χ 2 for the NR1 Fourier spectra obtained for spin-½ pair approximation and with isolated CH3 and CH2 functional groups approximation is pared with the differences between NR6 and NR7 spectra. In conclusion, since the differences between any two samples are much largest than the error due to the bad choice of the model, the Fourier spectra of a series of aged cross-linked natural rubber samples can be well characterized by any model, presented earlier as approximations. Figure 5. Comparison between DQ Fourier spectra obtained with spin-1 2 pair approximation and CH 2 and CH 3 isolated group's approximation corresponding to (a) NR1 and (b) NR7-aged samples. In order to be compared, the Fourier spectra were renormalized to have the same integral area.

The Distributions of End-to-End Distance and Residual Dipolar Coupling
The Fourier spectra of DQ curves for the aged natural rubber present the general features observed earlier for the distribution of residual dipolar coupling described in references [15,17,40]. As in the case of a 2 wt.% PDMS, see ref. [40], our DQ Fourier spectra consist of a high-narrow peak at lower residual dipolar coupling values and a tail which decays slowly with the increase in the residual dipolar coupling values.
This resemblance between our Fourier spectra and the distribution of the residual dipolar coupling obtained from FTIKREG based on Tikhonov regularization [15] and then a fit [41][42][43] allow us to assume that the Fourier spectra features are mainly due to the distribution of the residual dipolar couplings and the relaxation times T 2 will just broaden the obtained spectra. Therefore, in the following sections we characterize the DQ Fourier spectra for the cross-linked natural rubber samples in terms of the distribution of the residual dipolar couplings related to the distribution of the end-to-end vector of the polymer chains.
The 3D Gaussian probability distribution of the end-to-end vector, → R for a cross-linked polymer network, is given by where R 2 is the average square of the end-to-end distance and the distribution over the length of the end-to-end vector, R, which satisfies the normalization relation The Gaussian distribution of the length of the end-to-end vector can be written as From Equations (11)- (13), the corresponding distribution of the residual dipolar coupling constant is given by the Г function [15,17]: where ω D is the residual dipolar coupling and ω D is the mean residual dipolar coupling. This distribution fast increases around ω D /2π = 0 and then slowly after a maximum decay with a long tail at large ω D /2π. By its nature, the Г distribution is a broad one; it does not consider the powder average (see simulations from supplementary information) and therefore cannot explain the narrow peak observed at lower ω D values. This narrow peak corresponds to a Gaussian distribution of residual dipolar coupling, which is defined as where ω 0 D is the maximum value of the Gaussian distribution and ∆ω D is the width of Gaussian function. This function leads to a distribution of the end-to-end vector of the form In conclusion, the distribution of residual dipolar coupling may be described by the superposition of Г function and Gaussian function described by Equations (14) and (15)

The Characterization of DQ Fourier Spectra
As discussed earlier, the DQ Fourier spectra of aged cross-linked natural rubber consist of two components. The fits of theses Fourier spectra with a sum of Г and Gaussian functions for each component leads to unsatisfied results. In fact, even by increasing the number of components, the presence of the Г function leads to unsatisfied fits of the DQ Fourier spectra. The best fit (continuous line) was found when the DQ Fourier spectrum (open circles) was deconvoluted with four Gaussian functions (continuous line, dashed line, dash-dot line and small dash line). The deconvolution of 1 H DQ Fourier spectra for the aged NR1 and NR7 are presented in Figure 6. The deconvolution of the DQ Fourier spectra corresponding to the aged NR1 sample approximate the spectra well while some inconsistencies can be observed for NR7. From these deconvolutions, one can observe that there are two Gaussian distribution (continuous and dashed lines in Figure 6) responsible for the fit of the peak located at small residual dipolar coupling values, and two Gaussian distribution (dash-dot line and small dash line) responsible for the fit of the small peak located at larger residual coupling values.

( )
In conclusion, the distribution of residual dipolar coupling may be described by the superposition of Г function and Gaussian function described by Equations (14) and (15), respectively. More details about the mediation over azimuthally angle β and end-to-end vector can be found in the Supplementary Information.

The Characterization of DQ Fourier Spectra
As discussed earlier, the DQ Fourier spectra of aged cross-linked natural rubber consist of two components. The fits of theses Fourier spectra with a sum of Г and Gaussian functions for each component leads to unsatisfied results. In fact, even by increasing the number of components, the presence of the Г function leads to unsatisfied fits of the DQ Fourier spectra. The best fit (continuous line) was found when the DQ Fourier spectrum (open circles) was deconvoluted with four Gaussian functions (continuous line, dashed line, dash-dot line and small dash line). The deconvolution of 1 H DQ Fourier spectra for the aged NR1 and NR7 are presented in Figure 6. The deconvolution of the DQ Fourier spectra corresponding to the aged NR1 sample approximate the spectra well while some inconsistencies can be observed for NR7. From these deconvolutions, one can observe that there are two Gaussian distribution (continuous and dashed lines in Figure  6) responsible for the fit of the peak located at small residual dipolar coupling values, and two Gaussian distribution (dash-dot line and small dash line) responsible for the fit of the small peak located at larger residual coupling values.  The DQ Fourier spectra can be characterized by the residual van Vleck moments. The second and fourth van Vleck moment for the DQ Fourier spectra obtained with the approximation of isolated CH2 and CH3 functional groups for the entire series of aged NR samples are presented in Table 1. The second M2 and fourth M4 residual van Vleck moments present a monotone variation with the increase in the cross-link density, with the exception of NR5. This is due to the fact that the van Vleck moment is calculated from DQ Fourier spectra after the determination of spectral maximum position 2 ⁄ , presented in Table 1 in the last column. The DQ Fourier spectra of aged NR samples consist The DQ Fourier spectra can be characterized by the residual van Vleck moments. The second and fourth van Vleck moment for the DQ Fourier spectra obtained with the approximation of isolated CH 2 and CH 3 functional groups for the entire series of aged NR samples are presented in Table 1 moments present a monotone variation with the increase in the cross-link density, with the exception of NR5. This is due to the fact that the van Vleck moment is calculated from DQ Fourier spectra after the determination of spectral maximum position ω max D /2π, presented in Table 1 in the last column. The DQ Fourier spectra of aged NR samples consist of two compo-nents; therefore, they are asymmetrical. From sample NR1 to sample NR4, the maximum was found for the first peak, which is characterized by small values of the residual dipolar coupling. Starting with NR5, the maximum jumps towards the second peak, which affects the dependence of the of the M 2 and M 4 van Vleck moments of the DQ Fourier spectra function of cross-link density. The second M 2 and fourth M 4 residual van Vleck moments for all Gaussian distributions obtained from deconvolution of the DQ Fourier spectra for the entire series of aged natural rubber samples are listed in Table 2. Fit errors were smallest than 5%.

The Effect of Cross-Link Density
The Gaussian distribution (14) is characterized by the center of distribution ω 0 D , where one can find the maximum probability and by the width of residual dipolar coupling distribution ω D (see Table 3). Linear dependences for all of these parameters which describe the distribution of residual dipolar couplings with the cross-link density were found (see Figure 7). The effect of an increase in the cross-link density was to proportionally increase the mean residual dipolar coupling for Gaussian distributions (see Figure 7a). At the same time, the increase in the number of cross-linking points leads to a larger inhomogeneity of the polymer network, which is observed in Figure 7b from the increase in the distribution width.  In all cases, the lines represent the best fit of the data.

The Laplace-like Analysis of Bimodal 1 H DQ Build-Up Curves
If we will take a look at the final expressions of DQ signals in Equation (4), we reach the conclusion that the differences between the mentioned approximations are small. Therefore, in any of the previous approximations we must analyze a signal of the form and since ω D (R, β) is a function of R and β, this can be rewritten as where we assumed an effective averaged transverse relaxation time. A complete mediation over end-to-end vector → R or/and azimuthally angle β, and the distribution of the residual dipolar coupling constant D res , is presented in the Supplementary Information with many details. The best fit of the DQ build-up curve measured for natural NR1 aged during six years of analyses by Laplace-like inversion using Equation (S12) with the kernel presented in SI.13 (see Supplementary Information) and with the best transverse relaxation times T * 2,1 = 0.6 ms and T * 2,2 = 2.6 ms, respectively, is presented in Figure 8a. Unfortunately, the best fit curve (dashed red line) cannot explain the bimodal character of the measured data for six years of naturally aged NR1. The overlap of the two signals originating from a large residual dipolar coupling constant D res (olive dotted line) characterized by T * 2,1 = 0.6 ms and from a small residual dipolar coupling constant D res (olive dotted line) characterized by T * 2,2 = 2.6 ms present just a small shoulder at the initial time regime. This fit is similarly with that presented in Figure 2a and analyzed in terms of second van Vleck moments M DQ 2 analyzed with Equation (4). The distributions of the small and large values of residual dipolar coupling constant D res are presented in Figure 8b. Slightly asymmetrically, these distributions are similarly with those reported in several papers for the residual dipolar coupling constant D res [41,42]. Then, the complete mediation over → R and β was proved not to be able to describe the measured bimodal DQ build-up curve for natural aged NR samples, but could be used to describe a bimodal DQ build-up curve for a non-aged PDMS1 sample [25] (for more details, see the Supplementary Information).  As we mentioned before, the 1 H DQ Fourier spectra are affected by the effective transverse relaxation time. Another method, which can provide us with a good resolution in the distribution of residual dipolar coupling, can be based on the Laplace inversion [26,27,31,32,39,44]. In fact, we must underline that the true Laplace inversion is characterized by an exponential kernel specific to magnetization relaxation processes, i.e.,  As we mentioned before, the 1 H DQ Fourier spectra are affected by the effective transverse relaxation time. Another method, which can provide us with a good resolution in the distribution of residual dipolar coupling, can be based on the Laplace inversion [26,27,31,32,39,44]. In fact, we must underline that the true Laplace inversion is characterized by an exponential kernel specific to magnetization relaxation processes, i.e., The quantity of interest is the distribution function f (T 2 ). This can be obtained using the fast inversion algorithm which assumes a problem written into a matrix form as [31,32] where the matrix M contain the measured data, K is the kernel, E stores the measurement noise and F is the desired distribution. We have adapted the problem by changing the exponential integrand kernel from Equation (19) into an sin 2 (ω D τ) as in Equation (18), and which can be solved using the FMI (Fast Matrix Inversion) algorithm. In this case, hereafter the inversion problem will be called a Laplace-like problem. Additionally, we must implement a method to obtain the effective transverse relaxation time, T * 2 . For that, we have considered a minimum T * 2,min and a maximum T * 2,max value for effective transverse relaxation time, and for values between these two limits we tested the fit of the experimental data, or more specifically the merit function χ 2 as it is described by Equation (10).
The experimental 1 H DQ build-up curves and fits with Laplace-like inversion procedure are presented in Figure 9 for natural aged samples NR1, NR4 and NR7. The inverse Laplace-like procedure can fit well the 1 H DQ build-up curve recorded for the aged NR1 sample up to the maximum (τ ∼ = 1.4 ms), but there are some deviations after this maximum (see Figure 9a). Nevertheless, this is a much better fit compared with the approximation obtained before with multi-moments method (see Figure 2), described in Section 3, since now both build-up components are fitted. The deviation observed at larger τ values is due to the fact that the special kernel used in this case sin 2 (ω D τ) is periodical and we assumed a single effective relaxation time T * 2 .
Polymers 2021, 13, x FOR PEER REVIEW 17 of 27 (20) where the matrix M contain the measured data, K is the kernel, E stores the measurement noise and F is the desired distribution. We have adapted the problem by changing the exponential integrand kernel from Equation (19) into an as in Equation (18), and which can be solved using the FMI (Fast Matrix Inversion) algorithm. In this case, hereafter the inversion problem will be called a Laplace-like problem. Additionally, we must implement a method to obtain the effective transverse relaxation time, * . For that, we have considered a minimum , * and a maximum , * value for effective transverse relaxation time, and for values between these two limits we tested the fit of the experimental data, or more specifically the merit function 2 χ as it is described by Equation (10).
The experimental 1 H DQ build-up curves and fits with Laplace-like inversion procedure are presented in Figure 9 for natural aged samples NR1, NR4 and NR7. The inverse Laplace-like procedure can fit well the 1 H DQ build-up curve recorded for the aged NR1 sample up to the maximum (τ ≅ 1.4 ms), but there are some deviations after this maximum (see Figure 9a). Nevertheless, this is a much better fit compared with the approximation obtained before with multi-moments method (see Figure 2), described in Section 3, since now both build-up components are fitted. The deviation observed at larger τ values is due to the fact that the special kernel used in this case is periodical and we assumed a single effective relaxation time * . A better approximation of the measured 1 H DQ build-up curves are observed for the naturally aged rubber samples with higher cross-link density (see Figure 9b,c). In these cases, the inverse Laplace-like procedure can fit well the entire experimental 1 H DQ build-up curves. Some extremely small oscillations can be observed due to the periodicity of the used kernel and probably due to the fact that the second component in the distribution of RDCs associated with the aging effect is much smallest in these cases.

Proton DQ Laplace-like Spectra of Aged Natural Rubber
The normalized distributions ( ) D f ω of the residual dipolar couplings, or with other words, the 1 H DQ Laplace-like spectra, are presented in Figure 10a for all naturally aged rubber samples. All 1 H DQ Laplace-like spectra consist of four well-resolved peaks with one exception: the middle peaks for the aged NR1 sample. For small and large values of cross-link density, the main peak is characterized by a reduced value of averaged residual dipolar couplings ( This is not a surprising result, since in some of our previous results [19,39], the NR4 sample was also showing a different behavior for various measured microscopic NMR parameters and elasticity modulus. A better approximation of the measured 1 H DQ build-up curves are observed for the naturally aged rubber samples with higher cross-link density (see Figure 9b,c). In these cases, the inverse Laplace-like procedure can fit well the entire experimental 1 H DQ buildup curves. Some extremely small oscillations can be observed due to the periodicity of the used kernel and probably due to the fact that the second component in the distribution of RDCs associated with the aging effect is much smallest in these cases.

Proton DQ Laplace-like Spectra of Aged Natural Rubber
The normalized distributions f (ω D ) of the residual dipolar couplings, or with other words, the 1 H DQ Laplace-like spectra, are presented in Figure 10a for all naturally aged rubber samples. All 1 H DQ Laplace-like spectra consist of four well-resolved peaks with one exception: the middle peaks for the aged NR1 sample. For small and large values of cross-link density, the main peak is characterized by a reduced value of averaged residual dipolar couplings (ω D /2π = 100 − 250 Hz). From this point of view, the samples NR4 and NR5 are not in the range and present the first peak at ω D /2π = 200 − 250 Hz. This is not a surprising result, since in some of our previous results [19,39], the NR4 sample was also showing a different behavior for various measured microscopic NMR parameters and elasticity modulus. A better approximation of the measured 1 H DQ build-up curves are observed for the naturally aged rubber samples with higher cross-link density (see Figure 9b,c). In these cases, the inverse Laplace-like procedure can fit well the entire experimental 1 H DQ build-up curves. Some extremely small oscillations can be observed due to the periodicity of the used kernel and probably due to the fact that the second component in the distribution of RDCs associated with the aging effect is much smallest in these cases.

Proton DQ Laplace-like Spectra of Aged Natural Rubber
The normalized distributions ( ) D f ω of the residual dipolar couplings, or with other words, the 1 H DQ Laplace-like spectra, are presented in Figure 10a for all naturally aged rubber samples. All 1 H DQ Laplace-like spectra consist of four well-resolved peaks with one exception: the middle peaks for the aged NR1 sample. For small and large values of cross-link density, the main peak is characterized by a reduced value of averaged residual dipolar couplings ( This is not a surprising result, since in some of our previous results [19,39], the NR4 sample was also showing a different behavior for various measured microscopic NMR parameters and elasticity modulus. The best-fitted effective relaxation times * for all aged NR1-NR7 samples are presented in Figure 10b as a function of cross-link density. As in the case of correction time obtained by 1 H DQ Fourier spectra (see Figure 4), a linear dependence of * as a function of cross-link density can be observed, with the exception of the NR1 sample. From our observation, the value of * can have a certain influence on the 1 H DQ Laplace-like spectra (Figure 10a), but in the fitting error limit (see Figure 10b) this can be negligible compared with the differences between DQ Laplace-like spectra. Figure 11 presents a comparison between the normalized 1 H DQ Laplace-like spectra corresponding to aged and unaged NR1 samples. Additionally, the unaged 1 H DQ Laplace-like spectrum is composed of four peaks but with a different distribution. Two large distribution peaks are located at lower residual dipolar coupling values ( 2 ⁄ = 200 − 400 Hz), and two smallest peaks are located at larger residual dipolar coupling values ( 2 ⁄ ≈ 1.075 − 1.3 kHz ). The peaks located at lower residual dipolar coupling values are in agreement with the results obtained by Nie et al. [34] and with the distributions reported in [35,36,45]. After an aging in natural conditions, many of the NR1 polymer chains characterized by a small value of 2π ω D become more mobile and have a reduced residual dipolar coupling (see the left large peak). The polymer chains with a small value of 2 ⁄ together with the polymer chain with 2 ⁄ ≈ 350 Hz become more rigid at 2 ⁄ = ≈ 450 − 650 Hz. After aging, the two peaks located at larger residual dipolar coupling values collapse into a small peak located between them at 2 ⁄ ≈ 1175 Hz. It is evident from Figure 11 that the aging process induces a much broader heterogeneity in RDC and is characterized by larger mobility of polymer segments. The best-fitted effective relaxation times T * 2 for all aged NR1-NR7 samples are presented in Figure 10b as a function of cross-link density. As in the case of correction time obtained by 1 H DQ Fourier spectra (see Figure 4), a linear dependence of T * 2 as a function of cross-link density can be observed, with the exception of the NR1 sample. From our observation, the value of T * 2 can have a certain influence on the 1 H DQ Laplace-like spectra (Figure 10a), but in the fitting error limit (see Figure 10b) this can be negligible compared with the differences between DQ Laplace-like spectra. Figure 11 presents a comparison between the normalized 1 H DQ Laplace-like spectra corresponding to aged and unaged NR1 samples. Additionally, the unaged 1 H DQ Laplace-like spectrum is composed of four peaks but with a different distribution. Two large distribution peaks are located at lower residual dipolar coupling values (ω D /2π = 200 − 400 Hz), and two smallest peaks are located at larger residual dipolar coupling values (ω D /2π ≈ 1.075 − 1.3 kHz ). The peaks located at lower residual dipolar coupling values are in agreement with the results obtained by Nie et al. [34] and with the distributions reported in [35,36,45]. After an aging in natural conditions, many of the NR1 polymer chains characterized by a small value of ω D /2π become more mobile and have a reduced residual dipolar coupling (see the left large peak). The polymer chains with a small value of ω D /2π together with the polymer chain with ω D /2π ≈ 350 Hz become more rigid at ω D /2π =≈ 450 − 650 Hz. After aging, the two peaks located at larger residual dipolar coupling values collapse into a small peak located between them at ω D /2π ≈ 1175 Hz. It is evident from Figure 11 that the aging process induces a much broader heterogeneity in RDC and is characterized by larger mobility of polymer segments.

Proton DQ Fourier and DQ Laplace-Like Spectra of Aged Natural Rubber
The 1 H DQ Fourier and DQ Laplace-like spectra can be directly compared (Figure 12) for the sample NR1, which is most affected by aging in natural condition. For this purpose, both spectra were renormalized, having the maximum amplitude 1 obtained at almost the same residual dipolar coupling value ω D /2π ≈ 100 Hz (see also the left vertical dashed line in Figure 12). Up to the maximum value, one can observe an excellent superposition of both DQ Fourier and Laplace-like spectra, then the DQ Laplace-like peak decays faster. The second group of two joined peaks at median values of ω D /2π in the 1 H DQ Laplacelike spectra can have a correspondent (see the second and third vertical dashed lines in Figure 12) with a shoulder in the 1 H DQ Fourier spectra, but the four and smallest peak from DQ Laplace-like spectra hardly can be associated with some shoulder in the DQ Fourier spectra. It is obvious from Figures 6 and 12 that the deconvolution of 1 H DQ Fourier spectra, with four Gaussian functions, can not be so well matched with the peak distributions from 1 H DQ Laplace-like spectra. Moreover, the peaks maximum from DQ Laplace-like spectra will not have a linear dependence function of cross-link density.

Proton DQ Fourier and DQ Laplace-Like Spectra of Aged Natural Rubber
The 1 H DQ Fourier and DQ Laplace-like spectra can be directly compared ( Figure  12) for the sample NR1, which is most affected by aging in natural condition. For this purpose, both spectra were renormalized, having the maximum amplitude 1 obtained at almost the same residual dipolar coupling value 2 ⁄ ≈ 100 Hz (see also the left vertical dashed line in Figure 12). Up to the maximum value, one can observe an excellent superposition of both DQ Fourier and Laplace-like spectra, then the DQ Laplace-like peak decays faster. The second group of two joined peaks at median values of 2 ⁄ in the 1 H DQ Laplace-like spectra can have a correspondent (see the second and third vertical dashed lines in Figure 12) with a shoulder in the 1 H DQ Fourier spectra, but the four and smallest peak from DQ Laplace-like spectra hardly can be associated with some shoulder in the DQ Fourier spectra. It is obvious from Figures 6a and 12 that the deconvolution of 1 H DQ Fourier spectra, with four Gaussian functions, can not be so well matched with the peak distributions from 1 H DQ Laplace-like spectra. Moreover, the peaks maximum from DQ Laplace-like spectra will not have a linear dependence function of cross-link density.

Proton DQ Fourier and DQ Laplace-Like Spectra of Aged Natural Rubber
The 1 H DQ Fourier and DQ Laplace-like spectra can be directly compared ( Figure  12) for the sample NR1, which is most affected by aging in natural condition. For this purpose, both spectra were renormalized, having the maximum amplitude 1 obtained at almost the same residual dipolar coupling value 2 ⁄ ≈ 100 Hz (see also the left vertical dashed line in Figure 12). Up to the maximum value, one can observe an excellent superposition of both DQ Fourier and Laplace-like spectra, then the DQ Laplace-like peak decays faster. The second group of two joined peaks at median values of 2 ⁄ in the 1 H DQ Laplace-like spectra can have a correspondent (see the second and third vertical dashed lines in Figure 12) with a shoulder in the 1 H DQ Fourier spectra, but the four and smallest peak from DQ Laplace-like spectra hardly can be associated with some shoulder in the DQ Fourier spectra. It is obvious from Figures 6a and 12 that the deconvolution of 1 H DQ Fourier spectra, with four Gaussian functions, can not be so well matched with the peak distributions from 1 H DQ Laplace-like spectra. Moreover, the peaks maximum from DQ Laplace-like spectra will not have a linear dependence function of cross-link density.

Comparison between Fourier and Laplace-like Methods
The 1 H DQ Fourier and Laplace-like spectra were applied to analyze a series of bimodal DQ build-up curves characteristic to aged samples since the classical methods fail to produce acceptable results. Nevertheless, these methods also have some disadvantages: the main criticism of 1 H DQ Fourier analysis is that the obtained spectrum is affected by the relaxation processes, which broadens the individual lines; therefore a supplementary deconvolution has to be applied. For both methods, we have to consider first a model in order to calculate the DQ signal function of excitation/reconversion times and then an automatic correction procedure has to be applied to obtain the 1 H DQ Fourier spectra and another procedure must be applied to calculate the effective relaxation time to obtain the 1 H DQ Laplace-like spectra. The Fourier analysis is well defined, which means that with the exception of a coefficient~1 2 in front of the exponential correction term, which cancels the negative peak, the Fourier spectra are uniquely defined. The procedure for Laplacelike inversion deal with an ill conditioned problem and the results depend on various internal parameters. Moreover, since in our case the particular kernel is periodic, special attention has to be paid to the upper limits of the residual dipolar coupling values. Another disadvantage of the use of 1 H DQ Fourier and Laplace-like analysis is the long measurement times. In order to obtain a DQ Fourier spectrum with no wiggles, due to the truncation of experimental data, we had to measure many points at large excitation/reconversion time τ. Despite all these, we found that the 1 H DQ Fourier and Laplace-like analysis give complementary results and can be used successfully to analyze multi-component DQ build-up curves, such as those recorded for the cross-linked natural rubber aged in natural conditions over six years.

The ad Hoc Abragam-like Function for the Distribution of Residual Dipolar Coupling Constant D res
In recent years, an empirical function, so-called Abragam-like, became preferred to be used to describe the DQ build-up curve on isotropic samples [41,42]. Moreover, K. Saalvächter and co. state that the Abragam-like build-up curve can be used in the fitting of any kind of DQ build-up data for homogeneous or inhomogeneous samples [41]. Then, in this case, the DQ signal can be described as where the distribution function was labeled with P (D A−l res ), where D A−l res can be viewed similar to a second-moment-type quantity and the specific kernel is given by [41,42] K D A−l res , τ = The prefactors of D res and the Weibull coefficient were considered from ref. [42,43]. In general, these coefficient and exponential factors can be optimized by fitting the DQ build-up curve with the Kernel (Equation (22)) multiplied with the exponent that represents the transverse relaxation process during the excitation and reconversion of double-quantum coherences.
The best fits of the experimental DQ build-up curves obtained using Equation (21) with the Abragam-like kernel are presented in Figure 13 for NR1, NR4 and NR7 natural rubber samples aged for six years. In all cases, one can remark a good fit of the experimental data. Among these, as expected, the largest value of the merit function (see Equation (10)) was obtained for NR1. This is due mostly to the fact that the experimental data are not so well approximated in the initial time regime up to τ ∼ = 1 ms. With the increase in cross-link density, the effect of aging is reduced and the DQ build-up curves are better approximated on the entire time scale (see Figure 13b,c).
In Figure 14a, the distributions of D A−l res for the six-year-aged cross-linked natural rubber series (NR1-NR7) obtained using Equation (21) with the Abragam-like kernel (Equation (22)) are presented. Compared with the series of distributions of the residual dipolar couplings ω D (shown in Figure 7a), in this case: (i) the main peak (with the largest integral area) is located at low D A−l res /2π value (~53.2 Hz for NR1 to~86.1 for NR7); (ii) a series of four (an additional one) peaks are observed at large D A−l res /2π value all with a more small amplitude; (iii) the variation of D A−l res -distribution in function of cross-link density is more smooth than in the case of ω D -distribution. Regarding the number of components and the particular domain, the obtained distribution of D A−l res /2π is similar to those reported by Chassé et al. in ref. [42] for a mixture of NR-C2 samples. In the experimental error limit, the variation of the effective transverse relaxation time T * 2 , with cross-link density, decay linearly for the entire series NR1 to NR7 (see Figure 14b). crease in cross-link density, the effect of aging is reduced and the DQ build-up curves are better approximated on the entire time scale (see Figure 13b,c). In Figure 14a, the distributions of for the six-year-aged cross-linked natural rubber series (NR1-NR7) obtained using Equation (21) with the Abragam-like kernel (Equation (22)) are presented. Compared with the series of distributions of the residual NR7); (ii) a series of four (an additional one) peaks are observed at large /2 value all with a more small amplitude; (iii) the variation of -distribution in function of cross-link density is more smooth than in the case of -distribution. Regarding the number of components and the particular domain, the obtained distribution of /2 is similar to those reported by Chassé et al. in ref. [42] for a mixture of NR-C2 samples. In the experimental error limit, the variation of the effective transverse relaxation time * , with cross-link density, decay linearly for the entire series NR1 to NR7 (see Figure 14b). Comparing the best fit of DQ build-up curves measured for aged NR1 and analyzed spin-½ pair approximation (blue continuous line in Figure 15a) where the curves are well approximated in the initial time regime, inclusive of the maximum doublet, in the case of an Abragam-like kernel (red dashed line in Figure 15a), the experimental data are much better approximated for a medium and large excitation/reconversion time, τ, but not so well in the region of the maximum doublet. Due to this fact, one can consider that the Comparing the best fit of DQ build-up curves measured for aged NR1 and analyzed spin-1 /2 pair approximation (blue continuous line in Figure 15a) where the curves are well approximated in the initial time regime, inclusive of the maximum doublet, in the case of an Abragam-like kernel (red dashed line in Figure 15a), the experimental data are much better approximated for a medium and large excitation/reconversion time, τ, but not so well in the region of the maximum doublet. Due to this fact, one can consider that the spin-1 /2 pair approximation (leading to so-called Laplace spectrum) is more precise in order to describe the effects of natural aging of natural rubber on the measured bimodal DQ build-up curve than in the case of the use of the Abragam-kernel. Nevertheless, one can expect similarly results. In order to test this hypothesis, the so-called Laplace-like and so-called Abragam-like distributions are represented together in Figure 15b for the NR1 natural sample aged for six years. The spectra are renormalized so that the maximum amplitude is 1. The Abragam-like distributions are rescaled with a factor of 0.583 present in Equation (22), where the Abragam-like kernel was defined. In the presented range, both distributions can be considered similarly. DQ build-up curve than in the case of the use of the Abragam-kernel. Nevertheless, one can expect similarly results. In order to test this hypothesis, the so-called Laplace-like and so-called Abragam-like distributions are represented together in Figure 15b for the NR1 natural sample aged for six years. The spectra are renormalized so that the maximum amplitude is 1. The Abragam-like distributions are rescaled with a factor of 0.583 present in Equation (22), where the Abragam-like kernel was defined. In the presented range, both distributions can be considered similarly. Figure 15. The comparison of DQ Laplace-like (continuous blue line) and Abragam-like (dashed red line) of (a) DQ build-up curve and (b) spectra measured for the NR1 natural rubber sample, naturally aged during six years. The spectra were renormalized in order to have the maximum value equal to one. Laplace-like spectrum is represented in function of averaged residual dipolar coupling, /2 while the Abragam-like spectrum is represented in function of /2 /0.583.

Conclusions
The DQ build-up curves for a series of cross-linked natural rubber aged in natural conditions for six years were characterized by 1 H DQ Fourier and Laplace-like spectra. For that, a numerical program was written in C++ to perform a correction with an effective relaxation time, which allows us to reveal the spectral distributions of the residual dipolar couplings. The 1 H DQ Fourier spectra was treated in terms of a superposition of four Gaussian distributions of the residual dipolar coupling. The parameters which are obtained as a result of correction with an effective relaxation time and as result of spectral Figure 15. The comparison of DQ Laplace-like (continuous blue line) and Abragam-like (dashed red line) of (a) DQ build-up curve and (b) spectra measured for the NR1 natural rubber sample, naturally aged during six years. The spectra were renormalized in order to have the maximum value equal to one. Laplace-like spectrum is represented in function of averaged residual dipolar coupling, ω D /2π while the Abragam-like spectrum is represented in function of D A−l res /2π/0.583.

Conclusions
The DQ build-up curves for a series of cross-linked natural rubber aged in natural conditions for six years were characterized by 1 H DQ Fourier and Laplace-like spectra. For that, a numerical program was written in C++ to perform a correction with an effective relaxation time, which allows us to reveal the spectral distributions of the residual dipolar couplings. The 1 H DQ Fourier spectra was treated in terms of a superposition of four Gaussian distributions of the residual dipolar coupling. The parameters which are obtained as a result of correction with an effective relaxation time and as result of spectral deconvolution seem to depend linearly on the cross-link density. The same measured 1 H DQ build-up curves were used to obtain the 1 H DQ Laplace-like spectra for the aged natural rubber samples. Four resolved Gaussian-like peaks were obtained in the 1 H DQ Laplace-like spectra. Three methods (based on Fourier analysis, spin-1 2 pair approximation and on the ad hoc Abragam-like kernel) were presented to successfully fit the DQ build-up curves affected by six years of natural aging of a series of cross-linked natural rubber, whereas the classical ones fail. Finally, the aging effects on the dynamics of the NR polymer chain were discussed, and we find out that the aging of NR in natural conditions will increase the mobility of the majority of mobile polymer segments. At the same time, while aging, part of these polymer chains will become more rigid.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/polym13203523/s1, Figure S1: (a) General set-up of a double quantum (DQ) NMR experiment; (b) five pulses DQ experiment and (c) seven pulses DQ experiment with two refocusing pulses in the middle of excitation and reconversion periods; Figure S2: Double quantum (DQ) build-up curves measured for a 6 years natural aged NR1 sample using the 5 pulse DQ sequence (open red square) and 7 pulse DQ sequence (open blue circle); Figure S3: (a) The simulated bar plot distribution of ω D given by the Equation (S10) for constant → q 2 D res considering an isotope angular dependence described by sin function of azimuthal angle β shown in the small inclusion; (b) the comparison between the positive simulated bar plot distribution of ω D from Figure S3 (a) and a specific powder spectrum; Figure S4: The simulated bar plot distribution of ω D given by the Equation (S10) for constant D res P 2 {cos(β)} considering only the distribution of dimensionless squared end-to-end vector; Figure S5: The simulated bar plot distribution of ω D given by the Equation (S10) for constant D res for a combined isotope angular dependence described by sin function of azimuthal angle β and (a) a 3D Gaussian distribution of the end-to-end vector and (b) an Gaussian distribution of dimensionless squared end-to-end vector with the center → q 0 = 2.5; Figure S6: (a) Simulated bar plot distributions of residual dipolar coupling constant P(D res ) for a: (I) mono-Gaussian with D 0 res = 500 rad/s; (II) mono-Gaussian with D 0 res = 2000 rad/s and (III) bi-Gaussian with D 0 res,1 = 500 rad/s and D 0 res,2 = 2000 rad/s. (b) The simulated DQ build-up curves using Equation (S12) and kernel (SI.13) for the simulated distributions presented in a) and T * 2 = 2 ms; Figure S7: (a) Simulated bar plot bi-Gaussian distributions of residual dipolar coupling P(D res ) with D 0 res,1 = 500 rad/s and D 0 res,2 = 5000 rad/s and the resulting residual dipolar coupling constant distribution P(D res ) (continuous line) after Laplace-like analysis based on Equation (S14) with T * 2,1 = 2 ms and T * 2,2 = 0.2 ms. (b) The simulated DQ build-up curve (open circles) using Equation (S14) and the fitting of simulated DQ build-up data (continuous line); Figure S8: (a) The DQ-build-up curve (open circle) similarly with those measured by Bertmer et al. in [24] for a PDMS1 sample, the fitting curve obtained after the Laplace-like inversion using Equation (S13) (dashed line), the DQ NMR signals corresponding to the small D res values (continuous orange line) and to the large D res values (dotted olive line); (b) The distributions of residual dipolar coupling constants resulted from the analysis of data presented in (a) by Laplace-like inversion using Equation (S14) with T * 2,1 = 0.22 ms and T * 2,2 = 3.92 ms.