A Bayesian analysis of neutron spin echo data on polymer coated gold nanoparticles in aqueous solutions

TL;DR

This paper introduces a Bayesian method for analyzing neutron spin echo data of PEG-coated gold nanoparticles, enabling model-free, statistically supported insights into their nanosecond dynamics and polymer corona relaxation.

Contribution

It presents a Bayesian approach to determine the optimal number of exponential terms in dynamical analysis, avoiding pre-selected models and overparametrization.

Findings

Polymer corona relaxation follows a pure exponential decay.

The Bayesian method identifies the most statistically supported dynamical model.

Disentanglement of nanoparticle diffusion from polymer dynamics.

Abstract

We present a neutron spin echo study (NSE) of the nanosecond dynamics of polyethylene glycol (PEG) functionalised nanosized gold particles dissolved in D$_2$O at two temperatures and two different PEG molecular weights. The analysis of the NSE data was performed by applying a Bayesian approach to the description of time correlation function decays in terms of exponential terms, recently proved to be theoretically rigorous. This approach, which addresses in a direct way the fundamental issue of model choice in any dynamical analysis, provides here a guide to the most statistically supported way to follow the decay of the Intermediate Scattering Functions I(Q, t) by basing on statistical grounds the choice of the number of terms required for the description of the nanosecond dynamics of the studied systems. Then, the presented analysis avoids from the start resorting to a pre-selected framework and can be considered as model free. By comparing the results of PEG coated nanoparticles with those obtained in PEG2000 solutions, we were able to disentangle the translational diffusion of the nanoparticles from the internal dynamics of the polymer grafted to them, and to show that the polymer corona relaxation follows a pure exponential decay in agreement with the behavior predicted by coarse grained molecular dynamics simulations and theoretical models. This methodology has one further advantage: in the presence of a complex dynamical scenario I(Q,t) is often described in terms of the Kohlrausch-Williams-Watts function that can implicitly represent a distribution of relaxation times. By choosing to describe the I(Q,t) as a sum of exponential functions and with the support of the Bayesian approach, we can explicitly determine when a finer-structure analysis of the dynamical complexity of the system exists according to the available data without the risk of overparametrisation.

Tables

Table 1. Table 1: Instrumental setups used for the various samples

Sample	$C$[w%]	Instrument	$\lambda$ [Å]	NSE time [ns]	$Q$ range [Å-1]
PEG2000 AuNP	5	IN15	10, 16	0.3-206	0.02-0.21
PEG400 AuNP	2	IN15	8, 16	0.18-206	0.02-0.25
PEG2000	10	IN15	8, 16.8	0.18-240	0.02-0.25
PEG2000	10	IN11C	8.5	0.02-5	0.09-0.27
PEG400	5	IN15	8, 16	0.18-206	0.02-0.25

Table 2. Table 2: Comparison between the values of the relaxation times and their weights obtained by fitting only the IN15 data and by fitting the nISF obtained by gathering together the data from IN11 and IN15.

Q [Å-1]	Instrument	${\tau }_1$ [ns]	${\tau }_2$ [ns]	$A_1$
0.12	IN15	6$\pm$1	16$\pm$4	0.59
	IN15 & IN11	6$\pm$1	12$\pm$ 4	0.60
0.14	IN15	3.7$\pm$ 0.8	10 $\pm$ 2	0.46
	IN15 & IN11	4.7 $\pm$ 0.8	10$\pm$ 4	0.60
0.20	IN15	1.8$\pm$ 0.3	5.3 $\pm$ 0.8	0.49
	IN15 & IN11	1.6 $\pm$ 0.2	5.2$\pm$ 0.8	0.62
0.23	IN15	1.6$\pm$ 0.2	5 $\pm$ 1	0.69
	IN15 & IN11	1.3 $\pm$ 0.2	4.0$\pm$ 0.8	0.69

Table 3. Table S 1: Fitting parameters obtained for the NSE measurements taken on the PEG2000 AuNP at 280 K and at λ = 10 𝜆 10 \lambda=10 Å. P( k | y conditional 𝑘 𝑦 k|y ) is the probability for the most frequently chosen model. The constant γ 𝛾 \gamma (that was not reported in this table and in the next ones) in front of Eq. 10 in the manuscript, resulted always in the range 0.98–1 with few exceptions within a 5% difference from 1.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	${\tau }_{pol}\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.04	1	$430\pm 10$		1	0.82	$1.42\pm 0.06$
0.05	$0.93\pm 0.02$	$300\pm 10$	$16\pm 5$	2	1	$1.28\pm 0.07$
0.06	$0.54\pm 0.05$	$200\pm 20$	$24\pm 5$	2	1	$1.3\pm 0.1$
0.07	$0.25\pm 0.04$	$150\pm 10$	$18\pm 2$	2	0.51	$1.3\pm 0.1$
0.08	$0.10\pm 0.03$	$110\pm 10$	$15\pm 0.7$	2	1	$1.3\pm 0.1$
0.09	$0.07\pm 0.02$	$89\pm 6$	$11.7\pm 0.6$	2	0.56	$1.3\pm 0.1$
0.11	$0.05\pm 0.02$	$69\pm 7$	$9.5\pm 0.5$	2	0.69	$1.3\pm 0.1$
0.12	$0.10\pm 0.3$	$56\pm 6$	$7.8\pm 0.5$	2	0.72	$1.3\pm 0.1$
0.13	$0.11\pm 0.3$	$48\pm 5$	$6.3\pm 0.5$	2	0.73	$1.3\pm 0.1$
0.14	$0.10\pm 0.04$	$39\pm 4$	$5.6\pm 0.7$	2	0.72	$1.3\pm 0.1$
0.15	$0.09\pm 0.03$	$34\pm 3$	$4.6\pm 0.5$	2	0.86	$1.3\pm 0.1$
0.16	$0.08\pm 0.02$	$29\pm 2$	$4\pm 0.2$	2	0.83	$1.3\pm 0.1$
0.19	$0.06\pm 0.02$	$21\pm 2$	$2.3\pm 0.2$	2	0.93	$1.3\pm 0.1$
0.20	$0.06\pm 0.02$	$19\pm 2$	$2.0\pm 0.1$	2	0.93	$1.3\pm 0.1$
0.21	$0.07\pm 0.04$	$17\pm 2$	$1.8\pm 0.2$	2	0.94	$1.3\pm 0.1$

Table 4. Table S 2: Fitting parameters obtained for the NSE measurements taken on the PEG2000 AuNP at 280K and at λ = 16 𝜆 16 \lambda=16 Å.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	${\tau }_{pol}\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.019	1	$1890\pm 50$		1	0.98	$1.41\pm 0.04$
0.026	1	$1250\pm 40$		1	0.97	$1.18\pm 0.04$
0.032	1	$860\pm 40$		1	0.87	$1.11\pm 0.05$
0.044	1	$400\pm 10$		1	0.79	$1.29\pm 0.04$
0.050	$0.90\pm 0.02$	$310\pm 7$	$20\pm 7$	2	0.91	$1.30\pm 0.03$
0.057	$0.66\pm 0.04$	$230\pm 10$	$28\pm 9$	2	0.59	$1.31\pm 0.07$
0.076	$0.22\pm 0.06$	$132\pm 7$	$20\pm 4$	2	0.55	$1.32\pm 0.07$
0.082	$0.09\pm 0.04$	$111\pm 6$	$16\pm 2$	2	0.74	$1.32\pm 0.07$
0.089	$0.19\pm 0.05$	$96\pm 5$	$13\pm 2$	2	0.68	$1.32\pm 0.07$
0.110	$0.10\pm 0.04$	$62\pm 3$	$9.7\pm 0.9$	2	0.74	$1.32\pm 0.07$
0.116	$0.13\pm 0.05$	$56\pm 3$	$9\pm 1$	2	0.78	$1.32\pm 0.07$
0.123	$0.22\pm 0.07$	$50\pm 3$	$8\pm 2$	2	0.70	$1.32\pm 0.07$

Table 5. Table S 3: Fitting parameters obtained for the NSE measurements taken on the PEG2000 AuNP at 318 K and at λ 𝜆 \lambda = 10 Å. By way of example in the 5 th and the 6 th column of this table we report also the average estimated value of β 2 subscript 𝛽 2 \beta_{2} and the probability to have β 2 = 1 subscript 𝛽 2 1 \beta_{2}=1 when the two component model is chosen. ∗ For these datasets we set k = 2 𝑘 2 k=2 , for consistency with the results obtained for the other data sets, even if the model with two components is not the one with the highest posterior. Moreover the joint posterior distribution (conditional on k = 2 𝑘 2 k=2 ) is bimodal. The first maximum in the distribution corresponds to β 2 = 1 subscript 𝛽 2 1 \beta_{2}=1 and values of the other parameters as given in the table. The second maximum corresponds to values of the parameters that are meaningless from a physical point of view. For example, at Q 𝑄 Q =0.051 Å -1 , the second mode corresponds to β 2 = 0.71 subscript 𝛽 2 0.71 \beta_{2}=0.71 , A 1 = 0.27 subscript 𝐴 1 0.27 A_{1}=0.27 , τ D subscript 𝜏 𝐷 \tau_{D} = 131.12 ns and τ p ​ o ​ l subscript 𝜏 𝑝 𝑜 𝑙 \tau_{pol} = 86.55 ns, i.e. values of the parameters completely incongruous with the trend observed for the other datasets. This second mode is probably a spurious solution, related to a local maximum and was, therefore, discarded. Thus, we simply set β 2 = 1 subscript 𝛽 2 1 \beta_{2}=1 for these datasets, in analogy with the values estimated for the remaining datasets. In all the other datasets we reported in the table the average value of β 2 subscript 𝛽 2 \beta_{2} when k = 2 𝑘 2 k=2 as estimated by the algorithm but a value of 1 for the stretching coefficient when the conditional probability P ​ ( β 2 = 1 | k = 2 ) > 0.5 𝑃 subscript 𝛽 2 conditional 1 𝑘 2 0.5 P(\beta_{2}=1|k=2)>0.5 was adopted for the final curve fitting.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	${\tau }_{pol}\left[\text{ns}\right]$	${\beta }_2$	P(${\beta }_2=1|k=2$)	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.041	1	$177\pm 2$				1	0.80	$3.42\pm 0.09$
0.051	$0.86\pm 0.02$	$128\pm 9$	$9\pm 2$	1	1*	2	1	$3.0\pm 0.2$
0.062	$0.55\pm 0.03$	$87\pm 9$	$10\pm 1$	1	1*	2	1	$3.1\pm 0.3$
0.072	$0.28\pm 0.02$	$60\pm 6$	$7.5\pm 0.3$	$1.000\pm 0.005$	0.97	2	0.77	$3.2\pm 0.3$
0.083	$0.18\pm 0.02$	$47\pm 5$	$5.6\pm 0.3$	1	1*	2	1	$3.2\pm 0.3$
0.094	$0.11\pm 0.02$	$36\pm 4$	$4.2\pm 0.2$	$1.000\pm 0.006$	0.97	2	0.84	$3.2\pm 0.3$
0.106	$0.04\pm 0.01$	$27\pm 3$	$3.1\pm 0.1$	$0.99\pm 0.02$	0.83	2	0.98	$3.3\pm 0.3$
0.127	$0.06\pm 0.02$	$19\pm 2$	$2.3\pm 0.1$	$0.999\pm 0.008$	0.96	2	0.98	$3.3\pm 0.3$
0.140	$0.08\pm 0.02$	$15\pm 2$	$1.8\pm 0.1$	$0.99\pm 0.03$	0.86	2	0.98	$3.4\pm 0.3$
0.150	$0.10\pm 0.03$	$13\pm 1$	$1.50\pm 0.09$	$0.99\pm 0.03$	0.80	2	0.98	$3.4\pm 0.3$
0.161	$0.07\pm 0.02$	$12\pm 2$	$1.30\pm 0.08$	$0.99\pm 0.02$	0.91	2	0.98	$3.3\pm 0.3$
0.190	$0.03\pm 0.02$	$8.5\pm 0.9$	$0.95\pm 0.05$	$0.99\pm 0.02$	0.89	2	0.99	$3.3\pm 0.3$
0.200	$0.04\pm 0.02$	$7.6\pm 0.8$	$0.84\pm 0.06$	$0.99\pm 0.02$	0.91	2	0.99	$3.3\pm 0.3$
0.211	$0.03\pm 0.02$	$6.9\pm 0.7$	$0.78\pm 0.06$	$0.99\pm 0.02$	0.91	2	0.99	$3.3\pm 0.3$

Table 6. Table S 4: Fitting parameters obtained for the NSE measurements taken on the PEG2000 AuNP at 318 K and at λ = 16 𝜆 16 \lambda=16 Å. *Refer to note in Table 3 .

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	${\tau }_{pol}\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.019	1	$760\pm 10$		1	0.76	$3.5\pm 0.2$
0.026	1	$488\pm 4$		1	0.58	$3.2\pm 0.3$
0.032	1	$312\pm 4$		1	0.75	$3.0\pm 0.1$
0.044	$0.92\pm 0.02$	$175\pm 6$	$16\pm 5$	2	1*	$3.1\pm 0.3$
0.051	$0.80\pm 0.02$	$135\pm 5$	$12\pm 2$	2	1*	$3.1\pm 0.3$
0.057	$0.64\pm 0.04$	$92\pm 7$	$10\pm 2$	2	1*	$3.4\pm 0.3$
0.076	$0.28\pm 0.06$	$52\pm 5$	$8\pm 1$	2	0.92	$3.4\pm 0.3$
0.082	$0.14\pm 0.05$	$46\pm 5$	$7.1\pm 0.6$	2	0.93	$3.3\pm 0.3$
0.089	$0.10\pm 0.04$	$38\pm 4$	$5.3\pm 0.6$	2	0.97	$3.4\pm 0.3$
0.110	$0.06\pm 0.03$	$26\pm 3$	$3.6\pm 0.2$	2	0.94	$3.2\pm 0.2$
0.116	$0.11\pm 0.04$	$23\pm 2$	$3.1\pm 0.3$	2	0.97	$3.3\pm 0.3$
0.123	$0.10\pm 0.04$	$21\pm 2$	$2.7\pm 0.3$	2	0.85	$3.2\pm 0.3$

Table 7. Table S 5: Fitting parameters obtained for the NSE measurements on the PEG400 AuNP taken at 280 K and at λ 𝜆 \lambda = 8 Å. ∗ For this particular dataset the data are scarcely informative about the number of components and the algorithm assigns relevant posterior probabilities to models with either 1, 2 or 3 components. Therefore, in this case, model choice takes advantage from results obtained for similar data sets, for which the posterior distribution of the number of components is more peaked and the model with 1 component is largely supported.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.047	1	$173\pm 7$	1	0.69	$2.6\pm 0.2$
0.057	1	$128\pm 5$	1	0.35∗	$2.4\pm 0.2$
0.067	1	$90\pm 3$	1	0.57	$2.4\pm 0.2$
0.063	1	$93\pm 3$	1	0.83	$2.6\pm 0.3$
0.073	1	$71\pm 2$	1	0.84	$2.6\pm 0.2$
0.083	1	$58\pm 5$	1	0.59	$2.6\pm 0.3$
0.093	1	$42\pm 5$	1	0.79	$2.4\pm 0.3$

Table 8. Table S 6: Fitting parameters obtained for the NSE measurements on the PEG400 AuNP taken at 280 K and at λ 𝜆 \lambda = 16 Å.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.02	1	$800\pm 20$	1	0.89	$3.1\pm 0.1$
0.03	1	$458\pm 9$	1	0.59	$2.5\pm 0.1$
0.05	1	$191\pm 3$	1	0.78	$2.5\pm 0.1$
0.06	1	$124\pm 2$	1	0.71	$2.4\pm 0.3$
0.07	1	$69\pm 3$	1	0.91	$2.7\pm 0.3$
0.08	1	$54\pm 4$	1	0.92	$2.7\pm 0.2$
0.11	1	$34\pm 5$	1	0.69	$2.4\pm 0.4$

Table 9. Table S 7: Fitting parameters obtained for the NSE measurements on the PEG400 AuNP taken at 318 K and at λ 𝜆 \lambda = 8 Å.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.02	1	$800\pm 20$	1	0.89	$3.1\pm 0.1$
0.037	1	$99\pm 1$	1	0.87	$7.2\pm 0.4$
0.048	1	$64.4\pm 0.8$	1	0.89	$6.8\pm 0.4$
0.057	1	$44.7\pm 0.6$	1	0.91	$6.8\pm 0.5$
0.067	1	$31.6\pm 0.5$	1	0.95	$7.1\pm 0.2$
0.063	1	$35.7\pm 0.4$	1	0.94	$7.0\pm 0.1$
0.073	1	$25.0\pm 0.4$	1	0.96	$7.5\pm 0.3$
0.083	1	$19.3\pm 0.5$	1	0.91	$7.6\pm 0.6$
0.093	1	$17.1\pm 0.7$	1	0.90	$6.9\pm 0.4$
0.100	1	$15.0\pm 0.5$	1	0.92	$6.6\pm 0.3$
0.111	1	$13.2\pm 0.5$	1	0.96	$6.2\pm 0.3$
0.121	1	$10.8\pm 0.4$	1	0.97	$6.4\pm 0.3$
0.131	1	$9.5\pm 0.4$	1	0.96	$6.2\pm 0.3$
0.136	1	$8.7\pm 0.4$	1	0.94	$6.2\pm 0.5$
0.146	1	$7.5\pm 0.4$	1	0.96	$6.3\pm 0.3$
0.155	1	$8.3\pm 0.9$	1	0.83	$5.3\pm 0.6$
0.165	1	$6.6\pm 0.9$	1	0.92	$5.6\pm 0.7$
0.186	1	$4.7\pm 0.8$	1	0.93	$6\pm 1$

Table 10. Table S 8: Fitting parameters obtained for the NSE measurements on the PEG400 AuNP taken at 318 K and at λ 𝜆 \lambda = 16 Å.

Q [Å-1]	$A_1$	${\tau }_D\left[\text{ns}\right]$	k	P($k|y$)	$D_T\left[{\text{Å}}^2{\text{ns}}^{-1}\right]$
0.02	1	$800\pm 20$	1	0.89	$3.1\pm 0.1$
0.02	1	$329\pm 4$	1	0.74	$7.6\pm 0.2$
0.03	1	$173\pm 2$	1	0.87	$6.6\pm 0.4$
0.05	1	$69\pm 1$	1	0.94	$6.8\pm 0.3$
0.06	1	$48\pm 1$	1	0.93	$6.6\pm 0.2$
0.07	1	$28\pm 1$	1	0.94	$6.7\pm 0.3$
0.08	1	$23\pm 2$	1	0.85	$6.5\pm 0.6$
0.11	1	$12\pm 2$	1	0.86	$6.8\pm 0.8$
0.12	1	$10\pm 2$	1	0.87	$6.7\pm 0.9$

Table 11. Table S 9: Fitting parameters obtained for the NSE measurements on the PEG2000 in D 2 O at concentration of 10% in weight taken at 280 K and at λ 𝜆 \lambda = 8 Å.

Q [Å-1]	$A_1$	${\tau }_1\left[\text{ns}\right]$	${\tau }_2$ [ns]	k	P($k|y$)
0.06	1	24.3 $\pm$ 0.4		1	0.75
0.08	1	17.7 $\pm$ 0.2		1	0.65
0.09	$0.5\pm 0.2$	9 $\pm$ 2	23 $\pm$ 7	2	0.64
0.1	$0.6\pm 0.2$	7 $\pm$ 2	20 $\pm$ 5	2	0.64
0.12	$0.6\pm 0.2$	6 $\pm$ 1	16 $\pm$ 4	2	0.65
0.13	$0.6\pm 0.2$	4.9 $\pm$ 0.9	12 $\pm$ 3	2	0.63
0.14	$0.5\pm 0.2$	3.7 $\pm$ 0.8	10 $\pm$ 2	2	0.61
0.15	$0.4\pm 0.1$	2.8 $\pm$ 0.6	7.5 $\pm$ 0.9	2	0.62
0.16	$0.4\pm 0.2$	2.3 $\pm$ 0.6	6 $\pm$ 1	2	0.57
0.19	$0.5\pm 0.2$	1.8 $\pm$ 0.3	5 $\pm$ 1	2	0.63
0.2	$0.5\pm 0.1$	1.8 $\pm$ 0.3	5.3 $\pm$ 0.8	2	0.61
0.21	$0.6\pm 0.1$	1.7 $\pm$ 0.2	6 $\pm$ 1	2	0.61
0.23	$0.69\pm 0.07$	1.6 $\pm$ 0.2	5 $\pm$ 1	2	0.66
0.24	$0.7\pm 0.1$	1.6 $\pm$ 0.2	5 $\pm$ 1	2	0.67
0.25	$0.69\pm 0.05$	1.11 $\pm$ 0.09	6 $\pm$ 2	2	0.73

Table 12. Table S 10: Fitting parameters obtained for the NSE measurements on the PEG2000 in D 2 O at concentration of 10% in weight taken at 318 K and at λ 𝜆 \lambda = 8 Å.

Q [Å-1]	$A_1$	${\tau }_1\left[\text{ns}\right]$	${\tau }_2\left[\text{ns}\right]$	k	P($k|y$)
0.06	1	10.7 $\pm$ 0.1		1	0.55
0.08	$0.5\pm 0.2$	5 $\pm$ 1	11$\pm$ 2	2	0.64
0.09	$0.5\pm 0.2$	3.8 $\pm$ 0.9	7$\pm$ 1	2	0.63
0.1	$0.7\pm 0.1$	3.3 $\pm$ 0.3	8.5$\pm$3	2	0.68
0.12	$0.6\pm 0.2$	2.3 $\pm$ 0.4	4.8$\pm$0.9	2	0.72
0.13	$0.7\pm 0.1$	1.9 $\pm$ 0.2	5 $\pm$ 1	2	0.73
0.14	$0.7\pm 0.2$	1.7 $\pm$ 0.2	3.6$\pm$0.9	2	0.78
0.15	$0.7\pm 0.1$	1.3 $\pm$ 0.2	3.6$\pm$0.9	2	0.81
0.16	$0.8\pm 0.1$	1.2 $\pm$ 0.1	3.5$\pm$0.2	2	0.81
0.19	$0.7\pm 0.2$	0.9 $\pm$ 0.1	2.1$\pm$0.6	2	0.75
0.2	$0.90\pm 0.05$	0.89 $\pm$ 0.02	-	2	0.82
0.21	$0.97\pm 0.03$	0.73 $\pm$ 0.04	-	1	0.53
0.23	$0.93\pm 0.05$	0.64 $\pm$ 0.02	-	2	0.85
0.24	$0.95\pm 0.05$	0.52 $\pm$ 0.02	-	2	0.91
0.25	$0.94\pm 0.05$	0.44 $\pm$ 0.02	-	2	0.93

Equations

$${I}(Q,t)\simeq\exp(-DQ^{2}t)\frac{1}{N_{at}}\sum_{n,n^{\prime}=1}^{N_{at}}\bigl{\langle}b_{n}^{*}b_{n^{\prime}}\bigr{\rangle}\langle\exp[i\vec{Q}\cdotp(\vec{r}_{n}(t)-\vec{r}_{n^{\prime}}(0))]\rangle=\exp(-DQ^{2}t){I}_{\rm int}(Q,t)$$

$$I_{\mathrm{int}}(Q,t)=I_{\mathrm{int}}(Q,\infty )+I_{\mathrm{int}}^{{}^{\prime }}(Q,t).$$

$$\frac{I(Q,t)}{I(Q,0)}=\frac{I^{coh}(Q,t)-\frac{1}{3}{I}^{incoh}(Q,t)}{I^{coh}(Q,0)-\frac{1}{3}{I}^{incoh}(Q,0)}$$

$$I^{coh}(Q,t)=\frac{1}{N_{at}}\sum _{\begin{array}{c}n,n^{\prime }=1\end{array}}^{N_{at}}⟨b_n^{\ast coh}{b}_{n^{\prime }}^{coh}⟩⟨e^{-i\vec{Q}\cdot \overrightarrow{r_n}(0)}{e}^{+i\vec{Q}\cdot \overrightarrow{r_{n^{\prime }}}(t)}⟩$$

$$I^{incoh}(Q,t)=\frac{1}{N_{at}}\sum _{n=n^{\prime }=1}^{N_{at}}⟨(b_n^{incoh}{)}^2⟩⟨e^{-i\vec{Q}\cdot \overrightarrow{r_n}(0)}{e}^{+i\vec{Q}\cdot \overrightarrow{r_n}(t)}⟩$$

$$\frac{{I}_{\rm int}(Q,t)}{{I}_{\rm int}(Q,0)}=\frac{{I}_{\rm int}(Q,\infty)}{{I}_{\rm int}(Q,0)}+(1-\frac{{I}_{\rm int}(Q,\infty)}{{I}_{\rm int}(Q,0)})\sum_{j=2}^{k}I_{j}\exp\Bigg{(}-\frac{t}{\tau_{\rm int,\it j}}\Bigg{)}.$$

$$\frac{{I}(Q,t)}{{I}(Q,0)}\simeq\exp(-DQ^{2}t)\Bigg{[}A_{1}+\sum_{j=2}^{k}A_{j}\exp\Bigg{(}-\frac{t}{\tau_{\rm int,\it j}}\Bigg{)}\Bigg{]},$$

$$\frac{I(Q,t)}{I(Q,0)}\simeq A_1\exp (-\frac{t}{{\tau }_1})+\sum _{j=2}^k{A}_j\exp (-\frac{t}{{\tau }_j}).$$

$$\frac{I(Q,t)}{I(Q,0)}\simeq \sum _{j=1}^k{A}_j\exp (-\frac{t}{{\tau }_j}).$$

$$y_i=\gamma \sum _{j=1}^k{A}_j\exp \left(-{\left(\frac{t_i}{{\tau }_j}\right)}^{{\beta }_j}\right)+{ϵ}_i,\text{mbox}fori=1,\dots ,n$$

$$P(\Theta |y,I)\propto P(y|\Theta ,I)\times P(\Theta |I),$$

$$P(y|\Theta )=\prod _{i=1}^n\varphi \left(y_i;f(t_i),\nu {\sigma }_i^2\right),$$

$${\beta }_j\sim \zeta \mathcal{B}(\kappa ,\psi )+(1-\zeta ){\delta }_{{\beta }_j,1},$$

$$\gamma \sim \xi \mathcal{B}(\phi ,\rho )+(1-\xi ){\delta }_{\gamma ,1},$$

$$\hat{p}(k=\ell |y)=\sum _{m=1}^M{\delta }_{k^{(m)},\ell }/M=M_{\ell }/M$$

$${\hat{A}}_j=\sum _{m:k^{(m)}=\ell }{A}_j^{(m)}/M_{\ell },{\hat{\tau }}_j=\sum _{m:k^{(m)}=\ell }{\tau }_j^{(m)}/M_{\ell },\text{mbox}forj=1,\dots ,\ell$$

$$I(Q,t)/I(Q,0)=A(Q)\cdot F_{NP}(Q,t)+(1-A(Q))\cdot F_{Pol}(Q,t)$$

$$F_{NP}(Q,t)=\exp (-t/{\tau }_D)$$

$$F_{pol}(Q,t)=\exp \left[(-t/{\tau }_{pol}{)}^{\beta }\right]$$

$$exp(-D_R{Q}^{2}t)exp(-\frac{t}{{\tau }_{p=1}^{Rouse}})$$

$${\tau }_{p=1}^{Rouse}(Q)=\frac{3{\pi }^2}{N{b}^2}{\tau }_R\frac{1}{Q^2}=\frac{1}{D_R{Q}^2}$$

$$R=\frac{P(y|{\Theta }^{\prime })}{P(y|\Theta )}\frac{P({\Theta }^{\prime })}{P(\Theta )}\frac{q(\theta |{\theta }^{\prime })}{q({\theta }^{\prime }|\theta )},$$

$$R=\frac{{\prod }_{i=1}^n\varphi \left(y_i;\gamma {\sum }_{j=1}^k{A}_j^{\prime }\exp \left(-{\left(\frac{t_i}{{\tau }_j}\right)}^{{\beta }_j}\right),\nu {\sigma }_i^2\right)}{{\prod }_{i=1}^n\varphi \left(y_i;\gamma {\sum }_{j=1}^k{A}_j\exp \left(-{\left(\frac{t_i}{{\tau }_j}\right)}^{{\beta }_j}\right),\nu {\sigma }_i^2\right)}\times \frac{d(A^{\prime };\lambda )}{d(A;\lambda )}\times \frac{d(A;{ϵ}_A{A}^{\prime })}{d(A^{\prime };{ϵ}_{A}A)},$$

$${\nu }^{-1}\sim \mathcal{G}\left(\iota +\frac{n}{2},\varsigma +\frac{1}{2}\sum _{i=1}^n\frac{{\left(y_i-\gamma {\sum }_{j=1}^k{A}_j\exp \left(-{\left(\frac{t_i}{{\tau }_j}\right)}^{{\beta }_j}\right)\right)}^2}{{\sigma }_i^2}\right).$$

$${\tau }_{j^{\prime }}^{\prime }\sim \mathcal{U}(0,{\tau }_{\max })\text{mbox}and{\beta }_{j^{\prime }}^{\prime }\sim \zeta \mathcal{B}(\kappa ,\psi )+(1-\zeta ){\delta }_{{\beta }_{j^{\prime }}^{\prime },1}$$

$$B$$

Full text

I Introduction

Considerable scientific effort has been devoted in the last years to the study of gold nanoparticles (NP), in view of many important properties peculiar of these systems including quantum size effects and single electron transitions.Alivisatos (1996) Due to their small size (1-20 nm) the electronic bands of the nanoparticles follow quantum-mechanical rules and exhibit a strong dependence of the physical properties on the size and shape of the particles.Brust and Kiely (2002) These properties are further highly dependent upon the presence of an organic ligand coating or capping layer and upon the spatial relationship of one particle to another.

In order to contrast and control their tendency to aggregate, different methods have been developed to synthesise and stabilise such particles. Notable among these is the method developed by Brust et al.Brust et al. (1995) This method has had considerable impact since it has allowed the production of gold NP (AuNP) of controlled size and reduced size dispersity stabilised by alkanethiols. This has enabled a series of functionalisations that allow using the AuNPs in a number of applications including molecular recognitionLabande and Astruc (2000); Lin et al. (2002) and specific binding to biomoleculesReynolds et al. (2000) with significant impact in biological and biomedical studies.Daniel and Astruc (2004); Dreaden et al. (2012)

Functionalisation with thiol-terminated polymers provides great flexibility in modulating the properties of the AuNPs. It was shown that AuNPs protected with 4-(N,N-dimethylamino)pyridine (DMAP) are a convenient precursor of thiol-capped AuNPs. DMAP can be conveniently exchanged with thiol-terminated polymers to obtain AuNP that are soluble in organic solvents or water by choosing the appropriate thiol-terminated polymer including poly(ethylene oxide) (PEG) or polystyrene (PS) over a range of molecular weights.Rucareanu et al. (2008)

The properties of the organic capping layer that protects the AuNPs impact the way the AuNPs interact with other systems such as a surrounding polymer matrix,Genix and Oberdisse (2015) other NPs in complex 3D architectures,Nie et al. (2010) or a biomimetic assembly.Tatur et al. (2013) For practical applications it is crucial to have a very precise knowledge of the molecular details and properties of this protective layer surrounding the NP.

In a previous study, a structural characterisation of 5 nm PEG-functionalised AuNPs was achieved via a combination of experimental methods including neutron and X-ray small angle scattering, transmission electron microscopy, density, and thermogravimetric measurements.Maccarini et al. (2010) The functionalisation was obtained through thiol derivatisation of PEG2000 (45 units). The confinement of the PEG molecule on the highly curved surface of a AuNP was found to considerably affect the conformation and the hydration of the polymer as a function of the distance from the surface of the gold core.

Here, we extend this static study to the dynamical properties of these PEG-functionalised AuNPs. The study of the dynamical processes in nanosized objects like micelles, proteins, and nanoparticle dispersions presents some difficulties due to the presence of multi-scale dynamics occurring at similar time and range of length scales. Methods like Dynamic light scattering (DLS) are useful to measure the translational diffusion of the whole object under study. However, since the wavelength of light is larger than the dimension of these nanosized objects, internal dynamics cannot be explored with this technique. Quasi elastic neutron scattering (QENS) and, in particular, neutron spin echo spectroscopy (NSE) have been shown to be an excellent tool to probe the dynamical processes in micellar dispersions.Matsuoka et al. (2000); Kanaya et al. (2005); Castelletto et al. (2003) Thanks to its wide dynamical range, NSE can probe a variety of dynamical processes, from translational diffusion of the whole micelle to the internal dynamics of the polymer chains.

The polymer coated NPs used in this study have sizes of the same order of those of the micellar aggregates studied in other investigations (50-200 Å).Matsuoka et al. (2000); Kanaya et al. (2005); Castelletto et al. (2003) However, they differ significantly in various aspects. Firstly, unlike the core of a micellar aggregate, the gold core of the NP is a solid. The internal dynamics of gold in the solid state include atomic vibrations and phonons that occur at much faster timescales than those detectable with NSE. Hence, differently from the micellar aggregate comparators, in the systems reported here the dynamical processes present in the polymer corona and in the NPs core occur at significantly different timescales. Moreover, the nonionic surfactants, or block co-polymers of the micelles, aggregate in response to the hydrophobic effect, which favours the confinement of the hydrophobic moieties in the oily bulk. These aggregates are dynamic entities, as the individual surfactant molecules have a finite solubility in the solvent phase and can leave an aggregate to become part of another one, following a kinetic equilibrium. PEG-grafted gold NPs, on the other hand, involve a strong Au-S bond (ca. 160-200 kJ/mol) between the thiol-derivatised PEG to the metal surface. The exchange of material between different NPs, if it in fact occurs, does so at an exceedingly lower rate under ambient condition. NSE was also used to study NPs with various functionalisations dispersed in polymeric materials to elucidate the effect of the NP dispersion on the dynamical properties of the NPs-polymer composite.Senses et al. (2016, 2017); Campanella et al. (2016); Nusser et al. (2011); Mark et al. (2017) These studies differ from that reported here in that the NPs are dispersed in a high molecular weight polymer matrix, whereas the NPs used in our experiments were dispersed in a low molecular weight solvent.

The analysis of the NSE data was undertaken by adopting an approach composed by two building blocks:

we expressed the NSE intermediate scattering function $I(Q,t)$ as a sum of exponentials following the results of general theoretical work that proves that any time correlation function can be expressed as the sum of infinite (generally complex) exponential functionsBarocchi et al. (2012); Barocchi and Bafile (2013); Barocchi et al. (2014). This provided an unbiased model-free expression to interpret the experimental NSE results. 2. 2.

The analysis of the $I(Q,t)$ was supported by a Bayesian approach implemented through the use of a Monte Carlo Markov Chain (MCMC) algorithm integrated with a Reversible Jump (RJ)Green (1995) option and successfully applied recently to Brillouin inelastic neutron scattering De Francesco et al. (2016) and inelastic X-ray scatteringDe Francesco et al. (2018) data.

This methodology addresses one of the most important issues in the analysis of inelastic and quasi-elastic neutron scattering data - the arbitrariness in the choice of the number of components needed to describe the experimental dynamic structure factor. In NSE this issue translates in how many components one chooses to describe the decay of the intermediate scattering function (the time-Fourier transform of the dynamic structure factor). In this approach the number of components is treated as a parameter to be estimated. Since inferential results from this approach are in the form of probabilities, the final choice of the number of components resides on a solid probabilistic basis. The main objective of this work is then to show that the combination of points 1 and 2 above allows a precise and theory complying representation of macromolecular polymeric systems in solutions and yields an unbiased and model-free description of the dynamics of these systems.

Furthermore, it is common practice in the analysis of NSE data to use a stretched exponential function to model the decay of the intermediate scattering function. The use of this function is in some circumstances a heuristic way to account for an ensemble of dynamical processes occurring in the sample with a given distribution of relaxation times and the stretching parameter might measures, in some fashion, the difficulty in describing such a complex dynamics without however being able to give account for this complexity. It is also the results of model based theoretical description of polymer dynamics such as the Zimm or the Rouse Model. Nevertheless, all possible descriptions of the experimental $I(Q,t)$ of the system under study should eventually satisfy the above general property.Barocchi et al. (2012); Barocchi and Bafile (2013); Barocchi et al. (2014) This principle has been exploited to give excellent descriptions of time correlation functions and spectra relevant to simpler liquids.Bellissima et al. (2015, 2017); Guarini et al. (2017) Here, by combining it with the Bayesian approach, we analyse the data in an alternative unbiased and model-free way by assigning a *fine-structure * of relaxation processes to an otherwise undistinguished representation of the decay of the time correlation function as it is provided by the beta-stretching modelling. Furthermore and most importantly, the Bayesian method used here is able to provide the most statistically reliable description of the system supported by the experimental data without any risk of over-parametrisation.

We will also show that our unbiased method will provide results for the polymer corona in agreement with those predicted by coarse grained MD simulations and traditional theoretical models of polymer science.

II Materials and methods

Polyethylene glycol coated gold nanoparticles (PEG AuNPs) were synthesised with the same procedure described in our previous work. Maccarini et al. (2010) Briefly, tetraoctylammonium bromide-Au and DMAP-AuNPs were synthesised using Gittins and Caruso’s procedure. Gittins and Caruso (2001); Gandubert and Lennox (2005) The Au-NPs were subjected to place-exchange reaction of presynthesised DMAP stabilised NPs with thiol-functionalised PEG2000SH and PEG400SH. The gold core had a radius around 25 Å with a polydispersion of 20% as a result of TEM analysis.Rucareanu et al. (2008) However the indetermination in the Au core radius as a result of small angle neutron scattering was around 2%.Maccarini et al. (2010) The detailed synthesis and characterisation are described elsewhere. Rucareanu et al. (2008) PEG2000SH was purchased from Polymer Source Inc. Dorval, QC, Canada (Mw $\sim$ 2000) and PEG400SH (Mw $\sim$ 400) was received from Polypure Oslo, Norway. PEG2000 and PEG400 homopolymers were purchased from Sigma and used without further purification. The translational diffusion coefficients of the NPs were obtained by DLS. These measurements were carried out using an ALV CGS-3 Compact Goniometer equipped with a HeNe Laser with a wavelength of 632.8 nm, a 22 mW output power, and an ALV LSE-5004 Correlator. Samples were measured at a scattering angle of 90*∘, while the sample temperature was controlled via an external water bath circulator. The small angle neutron scattering profile for the PEG2000 AuNP was measured on the D11 SANS instrument at the Institut Laue-Langevin with a standard configuration. The neutron wavelength was 6 Å and two detector distances (10.5 and 1.5 m) were used. This provided a Q range between $6\cdot 10^{-3}$ and 0.43 Å-1* The PEG2000 AuNPs were dissolved in D2O at a concentration of 1% in weight.

II.1 Neutron Spin Echo Spectroscopy

NSE experiments were performed at the IN15 and IN11 neutron spin echo spectrometers at the Institut Laue-Langevin in Grenoble. NSE measurements were carried out in solutions of PEG2000-SH and PEG400-SH coated nanoparticles dispersed in D2O at concentrations of 5% wt and 2% wt, respectively, and in solutions of PEG2000 homopolymers (10% wt) and of PEG400 homopolymers (5% wt), also in D2O. The NSE signal of D2O alone was also measured in the same experimental conditions to allow for background subtraction. All measurements of $I(Q,t)$ were performed in a wide range of $Q$ and time at two temperatures, 280 and 318 K. Different setups were used depending on the sample and the spectrometer. A complete summary of the measured samples and experimental setups is given in Table 1.

II.2 Theory

We briefly present in this section the expressions on which are based the models adopted for analysing the neutron intermediate scattering functions (nISF) $I(Q,t)$ measured in the neutron spin echo experiments.

A number of customary approximation is generally adopted to obtain a workable expression of $I(Q,t)$. First of all the concentration of the nanoparticles in the solutions is supposed to be low enough to neglect any correlation between atoms of different nanoparticles. Then it is assumed that translational, rotational and vibrational motions are decoupled. This hypothesis is strictly valid for very dilute solutions but it is frequently used for systems at higher concentrationsLiu (2017) and quite plausible for our system as we will show later. Under these hypotheses it can be shown that the nISF, after subtraction of the measured solvent contribution, can be expressed as:

[TABLE]

where $D$ is the translational diffusion coefficient of the macromolecule, $N_{at}$ is the total number of atoms in a nanoparticle, $b_{n}$ is the scattering length of the $n$-th atom, and $r_{n}(t)$ is the displacement of an atom from its equilibrium position due to vibration or localised diffusive motionKneller , so that the last factor, labeled with the subscript ‘int’, describes the internal dynamics of the nanoparticle. As far as rotation are concerned it is useful to recall that the combination of independent translational and rotational motions has been shownSharma et al. (2012) to be satisfactorily described in the frequency domain by a single central Lorentzian line, though with a slightly increased width with respect to Fick’s Lorentzian of pure translational diffusion. This suggests that rotations can be accounted for by an effective, small, increase of the translational diffusion coefficient. However, in the specific case of the present samples, as reported in section III, no such increase has been detected, suggesting that the rotational contribution to Eq. (1) is negligible.

In liquid samples the position of a nucleus as $t\to\infty$ becomes uncorrelated with its initial value. Conversely, if the motion of the nucleus is well located in space, or confined to a particular volume (e.g. reorientations of chemical groups in a molecule), the correlation persists and tends to a nonvanishing asymptotic value at long times. The latter case applies to the internal dynamics of the polymer molecules grafted to the nanoparticles, therefore also in the present case it is possible, at any time, to write ${I}_{\rm int}(Q,t)$ as the sum of its time-independent and time-dependent parts:

[TABLE]

Here, the term ${I}_{\rm int}(Q,\infty)$ represents a quantity proportional to the contribution of purely elastic scattering to the differential scattering cross sectionSquires (1978) and is related to the structure of the nanoparticle.

NSE experiments measure a normalized intermediate scattering function (nISF), which is normally written as:

[TABLE]

where in general:

[TABLE]

[TABLE]

Both these terms are implicitly contained in Eq. (1) and contribute to $I(Q,t)$. In the experiment performed here the coherent term prevails especially at low Q ($\leq$ 0.1 Å*-1*). In the region of higher Q the incoherent terms reached 40-50% of the total contribution as it is shown in Figure S 1.

To further model ${I}^{{}^{\prime}}_{\rm int}(Q,t)$, we refer to the already quoted general theoretical results Barocchi et al. (2012); Barocchi and Bafile (2013); Barocchi et al. (2014), showing that any correlation function that decays to zero at long times is exactly represented by an infinite series of exponentials of generally complex argument, each exponential being weighted, in the series, by a generally complex amplitude. When the amplitude and time constant of an exponential term are complex, then the sum of such term with its complex conjugate (also present in the series) describes an exponentially damped oscillation. Differently, if amplitude and time constant are real, then the exponential term depicts a pure exponential decay.111The finding that discrete dynamical modes of exponential type can be identified in an ISF is also present in the work of Kikuchi et al.Kikuchi et al. (2013) However, this appears as the result of a numerical analysis rather than a property demonstrated in full generality. It is important to note that the possibility of expressing the ${I}^{{}^{\prime}}_{\rm int}(Q,t)$ as a sum of exponential terms is valid irrespectively of the fact that it represents, or is obtained from a coherent, incoherent scattering or a combination of the two.

Since we are dealing with quasielastic spectra, we can rather confidently assume that the normalized time correlation function ${I}^{{}^{\prime}}_{\rm int}(Q,t)$ does not require complex modes in the series representation, and that it can be reliably described by real modes only. Moreover, in Refs Bellissima et al. (2015, 2017); Guarini et al. (2017) it was shown that time correlation functions are most typically described with high accuracy by using only a small number (let’s say $k$) of terms in the series, thus finally we modeled the internal dynamics as

[TABLE]

The reason why the index $j$ starts at 2 will be clear immediately below. In Eq. (6) $\tau_{int,j}$ and $I_{j}$ are, respectively, the time constant and the amplitude of the $j$-th decay channel considered in the series (we omitted their obvious dependence on $Q$), and $\sum_{j=2}^{k}I_{j}=1$. The factor in front of the exponential series is imposed by the normalization condition $[{I}_{\rm int}(Q,t)/{I}_{\rm int}(Q,0)]_{t=0}=1$. Thus, by substituting Eq. (6) in Eq.(1), our system can actually be modeled as

[TABLE]

with $A_{1}=I_{\rm int}(Q,\infty)/I_{\rm int}(Q,0)$, $A_{j}=(1-A_{1})I_{j}$ for $j\geq 2$, and the sum rule $\sum_{j=1}^{k}A_{j}=1$. Defining $1/\tau_{1}=DQ^{2}$ and $1/\tau_{j}=1/\tau_{1}+1/\tau_{int,j}$ for $j\geq 2$, Eq. (7) can be re-written as

[TABLE]

As we will see below in the results section, for the AuNP at the lowest values of $Q$ in the explored range, the decay of the time correlation function is fully explained by the translational diffusive term ($A_{1}=1$) alone, therefore by what we will call from now on the relaxation time $\tau_{D}=\tau_{1}$. In the opposite limit of higher Q (i.e. $QR_{h}>>1$, being $R_{h}$ the hydrodynamic radius) the weight $A_{1}$ progressively vanishes and the decay is more and more explained by the internal polymer relaxation dynamics. In the intermediate Q-range the relaxation times $\tau_{j}$ in the second term of Eq. (8) are somehow affected by the translational diffusive contribution of the particle.

Since, as mentioned above, the term $I_{\rm int}(Q,\infty)$ is related to the structure of the nanoparticle, the term $A_{1}$ accounts implicitly for the form factor of the PEG AuNP.

In applying our RJ-MCMC approach (see next section), we will use for the AuNP the equation model (8) where the relaxation times $\tau_{j}$ for $j\geq 2$ will be explicitly referred to as $\tau_{pol,j}$. For the case of the homopolymer solutions of PEG2000 and PEG400, we used a more general expression of the nISF that reflects the fact that we did not make any assumption on the possible character of the first relaxation and that was instead inferred only after the data analysis:

[TABLE]

It is worth noticing that the formulation of a model describing the relaxation dynamics in terms of simple exponential functions, as in Eqs. (8) and (9), is not the standard choice in polymer physics, where stretched exponentials are often assumed. In order to demonstrate how the intermediate scattering function for the considered systems can be fully described as a linear combination of simple exponentials, we generalised the model equations (8) and (9) to sums of Kohlrausch-Williams-Watts exponential functions with stretched coefficients $\beta_{j}$. By means of the Bayesian approach described in Section II.3 applied to this model, we will provide a quantitative estimate of the statistical significance of stretching parameters and we will show that the choice of non-stretched ($\beta_{j}=1$) exponentials is fully adequate to an accurate description of the studied system dynamics and indeed the most valid on a probabilistic basis.

II.3 Bayesian analysis

Equation (9) describes the measured nISF relaxation as a linear combination of pure exponential components. As explained at the end of Section II.2, we generalised this model by replacing the simple with stretched exponential functions. Let $t=(t_{1},\ldots,t_{n})$ denote the vector of discrete times at which data were collected, and $y=(y_{1},\ldots,y_{n})$ the corresponding vector of the nISF measurements. Generalising Eq. (9), the data can be assumed to be sampled from the following expression:

[TABLE]

where $\gamma$ is a proportionality constant to account for a possible adjustement of the data normalisation and $\epsilon_{i}$, for $i=1,\ldots,n$, are independent random noises, accounting for statistical errors in the measurements. We assume these random errors to be independent and identically distributed according to a normal distribution $\mathcal{N}(0,\nu\sigma^{2}_{i})$, where $\sigma_{i}$ is the measurement error corresponding to the $i$-th observation and $\nu$ is a proportionality constant.

In order to choose an adequate number of components and estimate their parameters, we adopted a Bayesian approach. The main advantage of this approach in the present context is that the number of components, i.e. the relaxations in the time correlation function decay, can be treated as an unknown parameter to be estimated, along with the other model parameters. In addition, inferential results are obtained in the form of parameter probability density functions, so that the final choice on the number of relaxations is justified on a probabilistic basis. Finally, the choice of the Bayesian approach intrinsically and automatically applies the Ockam’s razor principleDe Francesco et al. (2016), thus preventing over-parametrisation.

Bayesian inference is based on Bayes theorem (Bayes and Price, 1753), which allows to combine evidence from data and possible prior knowledge $I$, through the equation

[TABLE]

where the symbol $|$ stands for “conditional on”, and $\Theta$ is the whole vector of unknown parameters, i.e. $\Theta=(k,A,\tau,\beta,\gamma,\nu)$, with $A=(A_{1},\ldots,A_{k})$, $\tau=(\tau_{1},\ldots,\tau_{k})$ and $\beta=(\beta_{1},\ldots,\beta_{k})$. We stress, again, that the number $k$ of elementary components in the model is also treated as a parameter.

The first term on the right-hand side of Eq. (11) is the likelihood function, which is the probability of the observed data $y$, conditional on a certain set $\Theta$ of parameter values, and on prior knowledge. In the present case the likelihood of the data is:

[TABLE]

where $\displaystyle f(t_{i})=\gamma\sum_{j=1}^{k}A_{j}\exp\left(-\left(\frac{t_{i}}{\tau_{j}}\right)^{\beta_{j}}\right)$ and \phi\left(\mbox{\boldmath\cdot};f(t_{i}),\nu\sigma_{i}^{2}\right) is the density of the $\mathcal{N}\left(f(t_{i}),\nu\sigma_{i}^{2}\right)$ distribution, evaluated in a specific point “$\cdot$”. The conditioning on $I$ here and in the rest of the paper has been dropped to simplify the notation.

The second term on the right-hand side of Eq. (11) is the prior distribution for the parameters of the model function, before data collection. It incorporates all prior knowledge, if any, on the problem at hand. In a Bayesian perspective, such prior distributions need to be specified for each parameter of the model, including the unknown number of mixture components, $k$.

In order to allow for complete a priori ignorance on $k$, we assume $p(k)=1/k_{\max}$, for $k=1,\ldots,k_{\max}$, where $k_{\max}$ is reasonably chosen to be the maximum possible number of components. The vector of weights, $A$, is assumed to follow a priori a Dirichlet distribution with hyperparameters $\lambda=(\lambda_{1},\ldots,\lambda_{k})$, i.e $A\sim\mathcal{D}(\lambda)$. Complete ignorance can be accommodated by letting $\lambda_{1}=\ldots=\lambda_{k}=1$ so that each component of the mixture has a marginal weight which is uniformly distributed in the interval $(0;1)$. Different prior beliefs on $k$ and $A$ can be easily accounted forDe Francesco et al. (2016). The stretching parameters $\beta_{j}$ are assumed a priori to independently follow a mixed distribution made up of a continuous beta density and a probability mass in 1, that is

[TABLE]

with $\kappa$ and $\psi$ being the hyperparameters of the beta density, and $\zeta$ denoting our prior support in favour of stretched (rather than unstretched) components. Notice that, letting $\zeta=0$ implies assuming a sum of pure exponential components, as the one in Eq. (9), whereas letting $\zeta=1$ determines a sum of stretched components. For any other $0<\zeta<1$, the mixture can include both stretched and unstretched components. In particular, we assumed $\zeta=0.5$, so that the probability of having a stretched or unstretched $j-$th component is, a priori, the same, for all $j=1,\ldots,k$. In addition, we set $\kappa=1$ and $\psi=1$ so that the stretching parameters are assumed uniformly distributed on the interval $(0;1)$. This corresponds to an uninformative prior giving the probability of 0.5 to an unstretched component and the remaining probability of 0.5 to a stretched component, and, for a stretched component, assigning the same density to any value of $\beta_{j}$ in $(0;1)$. Therefore our choice is that of being completely general: the same prior probability is given to any model made up of a sum of just stretched components or a sum of just unstretched components or a sum of stretched and unstretched components. Obviously models with just one stretched or unstretched component are also included when $k=1$. Our choice of prior satisfies both the requirements of the theory of correlation functions Barocchi et al. (2012); Barocchi and Bafile (2013); Barocchi et al. (2014) and the possibility to describe the results in terms of theoretical models used in polymer physicsRichter et al. (2005) without pushing the data in either one or the other direction. Obviously, different priors can be considered. For example, if one wants to give a higher prior to the Zimm model, he could choose $\kappa$ and $\psi$ such that $\kappa/(\kappa+\psi)\approx 0.85$ (so that the a priori expected value of $\beta_{j}$ would be 0.85), with $\psi$ accounting for the degree of prior belief in the Zimm model: the smaller $\psi$, the higher the prior probability assigned to this model. If, instead, one wants to consider both Zimm and Rouse model as a priori more probable than other solutions, a mixture of two beta densities could be considered for $\beta_{j}$, with two peaks centered on 0.85 and 0.5, respectively, and mixture weights proportional to the prior probabilities assigned to the two models.

The relaxation times $\tau_{j}$, for $j=1,\ldots,k$ are assumed to be a priori independent and identically distributed according to uniform densities in the interval $[0;\tau_{\max}]$, with $\tau_{\max}$ reasonably chosen. For the parameter $\gamma$, we also assume a mixed distribution made up of a continuous beta density and a probability mass in 1, that is

[TABLE]

with $\varphi$ and $\rho$ being the hyperparameters of the beta density, and $\xi$ denoting our prior support in favour of a refinement of the data normalisation process. As for $\zeta$, also in this case we let $\xi=0.5$. Finally, for conjugacy reasons, $\nu$ is assumed to have an inverse Gamma prior, $\nu^{-1}\sim\mathcal{G}(\iota,\varsigma)$, parametrised so that the mean and the variance are $\iota/\varsigma$ and $\iota/\varsigma^{2}$.

In the case of PEG AuNP data, where the first component of the model is meant to represent explicitly the centre-of-mass diffusion described as in Eq. (7), we specify slightly different priors. In fact, the stretching parameter is fixed at the value 1 by using a degenerate prior distribution $P(\beta_{1}=1)=1)$. Moreover, for the translational diffusion parameter, we assumed $D\sim\mathcal{N}(\mu_{D},\sigma^{2}_{D})$, with the hyperparameters chosen in such a way that $\mu_{D}$ is equal to the value $D_{LS}$ obtained by independent DLS measurements (see Section III.1) and $\sigma_{D}=0.20\mu_{D}$, so that $D$ can vary by up to $40\%$ from $D_{LS}$, with a probability of $95\%$. This was done to account for possible concentration effects on the translational diffusion, since the DLS experiment was done in a very dilute regime, whereas the NSE was done at a concentration of 5% in weight. Notice that, the relaxation time of the first exponential component, conditional on $D$, is not random anymore and is simply obtained as $\displaystyle\tau_{D}=\frac{1}{DQ^{2}}$.

Bayesian inference is based on the joint posterior distribution of all the random parameters, i.e. the left-hand side in Eq. (11). Due to the complexity of the model at hand, such posterior distribution is known up to a normalisation constant, which would require solving a high-dimensional integral to be evaluatedAndrieu et al. (2001); Razul et al. (2003) (i.e. the integral of the right-hand side of Eq. (11) over the parameters $\Theta$, which is required in the denominator of the right-hand side to have the joint posterior density integrate to one and the proportionality symbol replaced by equality). For this reason, we resort to MCMC algorithms to simulate the joint posterior distribution of interest.

MCMC algorithms allow to draw samples from a target distribution (the joint posterior in Bayesian applications), known up to a normalising constant. In this regard, these algorithms make it feasible to simulate the joint posterior distribution for the parameters of arbitrarily complicated models. MCMC algorithms are based on the construction of an irreducible and aperiodic Markov chain that has the target distribution as invariant distributionGilks et al. (1996). However, in our model the dimension of the parameter space is random, since it depends on the random value of $k$. Therefore, the MCMC algorithm needs to be incremented with a Reversible Jump Green (1995) (RJ) step, which allows to jump between models having different dimensions.

The RJ-MCMC algorithm we propose performs $M$ sweeps and, at each sweep, all the parameters are updated in turn, sampling the new value of a certain parameter conditionally on the data and all the remaining parameters and is able to explore models with any number of components as well as to move between models with stretched or unstretched components or any mixing of the two, thus including also Zimm and Rouse models. At each sweep, we perform a changing dimension move to update the value of $k$ (through RJ sampling), as well as fixed dimension moves to update all the other parameters (either through Gibbs or Metropolis-Hastings sampling (Tierney, 1994)). All the details about the algorithm and the updating moves are exhaustively explained in a dedicated paragraph of the Supplemental Material.Sup

After an initial burn-in period, during which the RJ-MCMC draws are discarded until the convergence of the chain to the target distribution is reached, the algorithm produces at each sweep $m$ a new realization of all the parameters of the model, sampled from their joint posterior distribution. Let $(k^{(m)},A^{(m)},\tau^{(m)},\beta^{(m)},\gamma^{(m)},\nu^{(m)})$, for $m=1,\ldots,M$, be the sample obtained after $M$ sweeps of the algorithm. This sample provides the simulated joint posterior distribution of all parameters and can be used to estimate all the quantities of interest.

The posterior distribution of the number of exponential components can be estimated as the proportion of times each model was visited by the algorithm, i.e.:

[TABLE]

where $M_{\ell}$ is the number of times the model with $\ell$ components was visited. In practice, the model chosen will be the one with the highest posterior probability.

Once a particular model has been chosen, we can estimate the parameters of each exponential component under that model. The estimates of these parameters can be computed as the mean of their simulated marginal posterior distribution, conditional on the model with $\ell$ components. Point estimates for the weights and relaxation times of the $\ell$ exponential components are given respectively by:

[TABLE]

where the sum over $m:k^{(m)}=\ell$ means over all values of $m$ in which the model with $k=\ell$ has been visited.

Note that in case the simulated marginal posterior distribution of a certain parameter is not symmetric, other estimators than the arithmetic mean can be considered (e.g. the simulated posterior mode).

For the stretching parameters, if the simulated posterior probability of $\beta_{j}=1$ is larger than its prior, i.e. $\displaystyle\sum_{m:k^{(m)}=\ell}\delta_{\beta_{j}^{(m)},1}/M_{\ell}\geq 1-\zeta$, then we simply let $\hat{\beta}_{j}=1$, otherwise $\displaystyle\hat{\beta}_{j}=\sum_{m\in S}\beta_{j}^{(m)}/|\mathcal{S}|,$ where $\mathcal{S}$ is the set of sweeps for which $k^{(m)}=\ell$ and $\beta_{j}^{(m)}<1$, and $|\mathcal{S}|$ is the dimension of $\mathcal{S}$.

All the other model parameters do not depend on the number of exponential components and can be estimated from their simulated marginal posterior distributions, averaging over all the sweeps of the algorithm.

III Experimental Results and Analysis

III.1 PEG2000AuNP and PEG400AuNP

A representative set of NSE data obtained for the PEG2000 AuNP is shown in Fig. 2(a). The full set is reported in the Supplemental Material (Fig. S 2).Sup As expected the nISF decays as a function of the NSE time and the decay becomes faster as Q increases.

In the investigated Q range we probed length scales ranging from 30 Å to about 300 Å, which include distances on the same scale of the intra- and inter-polymer chains interactions, as well as distances on the scale typical of translational diffusion of the whole nanoparticle. The translational diffusion coefficient was estimated by DLS measurements carried out on very dilute solution of the NPs in D2O ($\sim$ 0.02 weight %). The values obtained for this parameter were $D_{LS}$= 1.35 and 3.5 Å2/ns at temperatures 280 and 318 K, respectively. From the $D_{LS}$ we can extimate the hydrodynamic radius of the PEG2000 AuNP, $R_{h}\sim 85$ Å  through the relation $D={kT}/{6\pi\eta R_{h}}$ where $k_{B}$ is the Boltzmann constant, T the temperature and $\eta$ the viscosity of D2O  at the temperature T.

At low Q, where the condition $Q\cdot R_{h}\sim 1$ holds, the dynamics should be dominated by the translational diffusion of the AuNP, whereas at higher Q where $Q\cdot R_{h}>>1$ the dynamics should be dominated by the polymer relaxations. At intermediate Q values, a $Q$ dependent combination of these two contribution is expected. With this in mind and in analogy with Arriaga et al. Arriaga et al. (2009) we initially described the decay of ${I}(Q,t)/{I}(Q,0)$ as:

[TABLE]

where

[TABLE]

is the translational diffusion of the NP, and

[TABLE]

is the term approximately associated to the polymer dynamics (see above) often described with a Kohlrausch-Williams-Watts function. In the previous equations $A$ is the $Q$-dependent relative portion of translational diffusion with respect to the polymer dynamics contribution, $\tau_{D}=1/DQ^{2}$ is the characteristic time of the translational diffusion of the nanoparticle in D2O, $\tau_{pol}$ is a characteristic time associated to the polymer dynamics (although it contains also the translational diffusive dynamics (c.f. Eq. (7)), and $\beta$ the stretching exponent, which often accounts for the possibility of having a distribution of relaxation processes. However, we wanted to leave open the possibility to describe the polymer relaxations with a more general model. Consequently, we chose to describe the polymer relaxations with sums of exponential (stretched or simple) functions in order to find the combination that described the data with the greatest statistical significance.

In Eq. (8) we imposed the first component to be a simple exponential (which corresponds, as discussed in Sec.II.3, to set $P(\beta_{1}=1)=1$ as degenerate prior distribution for $\beta_{1}$), since it is associated to the translational diffusion of the NP (see for example Eq. (7)). As already discussed in Sec.II.3, we assumed $D\sim\mathcal{N}(\mu_{D},\sigma^{2}_{D})$. We allowed the remaining components to either be simple or stretched exponentials.

The results of the fitting procedures are listed in Tables S I, S II, S III and S IV in the Supplemental Material.Sup The column labelled $P(k|y)$ lists the probability associated to the most statistically supported model conditional to our data. In Fig. 2 we show a summary of the main results concerning the fitting analysis for the PEG2000 AuNP. For Q $<$ 0.04 Å*-1* the RJ-MCMC algorithm privileged a solution with only the translational diffusive component at both measured temperatures (Fig. 2(b) and 2(c)).

Above this value of Q, the solution with two relaxations was always chosen with a probability $P(k=2|y)$ usually between 70 and 90% especially at high Q and at the higher temperature. The mean $\beta_{2}$ stretching coefficient was estimated by the algorithm to be close to 1 and in any case the probability to have $\beta_{2}=1$ was most of the times greater than 90%.

In Fig. 2(b) the trend in Q of the parameter $A(Q)$ is displayed at the temperature of 280 K. At Q $<$ 0.04 Å*-1* the main contribution to the dynamics came from the translational diffusion of the NP. At $Q$ $>$ 0.1 Å*-1* the dominant contribution was the polymer dynamics. In the intermediate Q regime 0.04 Å*-1* $<$ Q $<$ 0.1 Å*-1*, we observed a transition between these two regimes.

The value of the translational diffusion coefficient estimated by the algorithm was in perfect agreement with the one measured by DLS (see Fig. 2(c), solid red line) which confirms that the contribution of the rotational term can actually be neglected. From Fig. 2(c) it is also clear that the relaxation times estimated at the two wavelengths are consistent with each other. As expected, the relaxation times decreased both with increasing the temperature and Q.

In Fig. 3 we report as example a typical set of results coming from the fitting algorithm obtained at $\lambda=10$ Å and at $Q=0.72$ Å*-1*. Figure 3(a) shows an example of the posterior distribution of the translational diffusion coefficient $D$ at $Q=0.72$ Å*-1*. In this specific case the simulated posterior distribution was a Gaussian, as shown by the fit (red curve). The spread of this distribution is related to the polydispersivity of the NP and to the errors bars of the experimental $I(Q,t)$. Figure 3(b) displays the posterior distribution of the number of components in Eq. (10) probed by the algorithm. In Fig. 3(c) and 3(d) the Trace Plots for the two relaxation times $\tau_{D}$ and $\tau_{pol}$ and for the A parameter are shown, always for the same $\lambda$ and $Q$. In Fig. 4 we report as an example the joint posterior distribution $P(\tau_{D},\tau_{pol}|k=2,y)$ for the two relaxation times marginalized to all the other model parameters as a function of $Q$ for the PEG2000AuNP sample at 318 K. It can be noticed that at small $Q$ values, the posterior distribution of $\tau_{D}$ is more dispersed than that of $\tau_{pol}$, implying a minor accuracy in the estimation of the first parameter, compared to the second one. As $Q$ increases, the difference between the dispersions of the two parameters vanishes and, in addition, the joint distribution becomes more concentrated, implying that more precise estimates of both parameters can be obtained at large $Q$ values. Also notice the absence of correlation between the two parameters.

A representative set of the NSE data measured on the PEG400 AuNPs dispersed in D2O are shown in Fig. 5(a). At both the investigated temperatures the most probable model predicted by the algorithm to describe the decay of ${I}(Q,t)/{I}(Q,0)$ was a single exponential relaxation. At 280 K, even if the quality of data at $Q\geq 0.8$ Å*-1* degraded (see Fig. S4 in the Supplemental MaterialSup ), the probability assigned by the algorithm to a one-component model was, except for a few datasets, significantly greater than 60%. The estimated relaxation time at each $Q$ value is very close to the one expected for the NP translational diffusion as measured by DLS ($D_{LS}=2.2$ Å2/ns; cf. Fig. 5(b), red solid lines). At 318 K the situation was clearer in view of an improved data quality. This allowed an unambiguous fitting of the relaxation curves given that the algorithm assigned to a one-component relaxation model a posterior probability $P(k=1|y)$ almost always above 90%. The translational diffusion coefficient from the fits is in agreement within the errors with that obtained with the DLS ($D_{LS}=6.2$ Å2/ns) and corresponds to a hydrodynamic radius of $\sim$ 48 Å. The complete set of parameters obtained by the fits are listed in Tables S V, S VI, S VII and S VIII in the Supplemental Material.Sup

III.2 Polymer solutions of PEG2000 and PEG400

III.2.1 PEG2000

In order to better characterise the dynamics of the polymer grafted to the NP, we compared the results obtained with the PEG coated NPs and those obtained with a reference system identified as a solution of PEG homopolymers. We used PEG2000 and PEG400 that possess a molecular weight very similar to those tethered to the NPs. As previously shown Maccarini et al. (2010), the concentration of polymer as a function of the distance from the gold surface is not uniform, but varies being rather high in the vicinity of the surface of the NP and lower away from it. Ideally, we should compare the results obtained on the PEG coated NP with several PEG solutions at different concentrations. The PEG2000 was studied in D2O  at weight concentrations of 10% 20% and 45% (the latter is the solubility limit for PEG2000 in D2O) as a reference. However, due to an unfavorable combination of coherent and incoherent components that contribute to the neutron polarization with opposite sign, the measurements at the higher concentrations (20 and 45% in weight) suffered of a too weak spin echo signal and were not taken into account. The choice of the reference system was based on a compromise between the counting statistics and the available measuring time. Ultimately, the PEG2000 and PEG400 were studied at concentrations of 10% and 5% in weight, respectively. A representative set of data obtained for the PEG2000 at 10% is displayed in Fig. 6(a).

The protocol used to analyse the NSE data of the polymers was slightly different from the one used for the NPs, as we did not use any informative prior for the first relaxation. In fact, we left the algorithm free to search for the most statistically probable solution. The prior distribution of the relaxation times was chosen as to be uniform in an interval $0\div\tau_{\max}$.

In a first attempt, we used a combination of simple exponential functions to describe ${I}(Q,t)/{I}(Q,0)$. The results for PEG2000 homopolymers at a concentration of 10% and at $\lambda=8$ Å are shown in Fig. 6(b). The parameters obtained from the fits are listed in Tables S IX and S X in the supporting information.

At T=280 K, the algorithm favors a single exponential for $Q\leq 0.08$ Å*-1* and two exponentials for $Q>0.08$ Å*-1* (see Fig. 6(b) and Fig. 7 top panels). At T=318 K we have a similar behavior up to $Q=0.19$ Å*-1*.

Above this value, the algorithm would still favor the solution with two components. However, while the parameters of the first component are nicely estimated to continue the trend of their respective values at lower $Q$s, the parameter $\tau_{2}$ results scarcely identified. In fact, the simulated posterior distribution of $\tau_{2}$ turns out to be substantially uniform in a large range of time (Fig. S 5(b)). As a consequence, any value of $\tau_{2}$ sufficiently large would define a second exponential function able to fit the decay of the time correlation function at large Spin Echo Time. We, therefore, decided to retain both components in the model (a single component would simply misconsider the last few points; see Fig. S 5(a) in the Supplemental MaterialSup ) but we do not trust estimates for the slower dynamics, since meaningless. The algorithm still provides the best probabilistically supported solution given the experimental data, but at the highest Q values the number of data points available for a precise determination of $\tau_{2}$ is insufficient. The application of our method particularly for datasets as the ones at $Q>0.19$ illustrates the efficacy of the Bayesian analysis since we can infer the reliability of the parameter estimates from the simulated marginal posterior distributions of the different parameters. This represents effectively a diagnosis tool to exploit at maximum the collected experimental data.

We stress that the posterior probability for the chosen values of $k$ is always larger or much larger than 60% (see an example in Fig. 8). Also, the estimate of the relaxation times is accurate for the majority of data sets and providing for the posterior parameter distribution function a Gaussian distribution or at least a well-shaped unimodal distribution although not necessarily symmetric (see Fig. 8 panel b) and c)).

Further measurements on the PEG2000 homopolymers dissolved in D2O  at 10% in weight were carried out on the spin echo spectrometers IN11 at ILL with an incident wavelength of $8.5$ Å. These measurements expanded the accessible time and momentum transfer, covering a time window between 0.02 and 5 ns and a Q range between 0.09 and 0.42 Å*-1*. In order to further validate the scenario so far achieved we joined the data collected with the two spin echo instruments IN11 and IN15 extending the bayesian MCMC data analysis to this enlarged time window. The available data from IN11 were only at 280 K. In Fig. 10 we show the relaxation curves at four selected momentum transfer wavevectors ($Q$=0.12, 0.14, 0.20 and 0.23 Å*-1*) at which we collected the data from both spectrometers.

At $Q$=0.12, 0.20 and 0.23 Å*-1* we re-obtain substantially the same results we found from the analysis of the IN15 data only, both for the characteristic relaxation times and the weights $A_{1}$ and (1-$A_{1}$). At Q=0.14 Å*-1* there is a slight difference for the faster relaxation, but the results are fully consistent within the uncertainties (Tab. 2).

III.2.2 PEG400

The study of the PEG400 homopolymer solution has been undertaken at a concentration of 5% wt. The description once again is rather simple. As in the case of the PEG400AuNP’s, one relaxation is sufficient to give account of the time correlation decays at both the investigated temperatures. The posterior distribution function for $k$ predicts a probability $P(k=1\mid y)$ for a single relaxation model always in the range $80\div 98\%$ for all the spectra analysed, with few exceptions where anyway such probability was larger than $60\%$. The relaxation times obtained at the two temperatures are summarised in Fig. 5(b). These relaxation times are coherent with the values of translational diffusion of PEG400 at the same temperature and concentration (cf. Fig. 5(b), dot red lines). Waggoner et al. (1995) A selection of the experimental NSE curves at the two investigated temperatures and at three selected $Q$ values and fits are displayed in Fig. S 6 of the Supplemental Material.Sup In Fig. S 7 we report two examples of the simulated posterior distribution function for the relaxation time $\tau_{1}$ for two different values of $Q$ at T=318 K. In these cases the posterior distribution functions look unimodal and almost symmetric. The relaxation times obtained for the PEG400 indicate a faster dynamics compared to that found in the PEG400 AuNP (Fig. 5b).

IV Discussion

The fitting analysis of the PEG2000 AuNP data details two relaxation processes that occur in the time window probed by the NSE experiment. The slower relaxation process with characteristic time $\tau_{D}$ corresponds to the translational diffusion of the PEG2000 AuNP. As it appears from Fig. 2(c), it perfectly overlaps with the relaxation times obtained from the expression $\tau_{D}=1/DQ^{2}$ where $D$ is the translational diffusion as derived by DLS measurement of PEG2000 AuNP at infinite dilution and at the respective temperatures. To further check if the interaction of the NP due to the finite concentration used in the experiments had a significant effect on the dynamics, we obtained the structure factor $S(Q)$ by performing SANS experiments at increasing concentrations (see Fig. S 3 in the Supplemental MaterialSup ). Following the de Gennes narrowing, a significant interaction would produce a modulation of the translational diffusion according to the inverse of the structure factor $1/S(Q)$. As it is shown in Fig. S 3 this modulation is barely visible within the error bars. This makes the assumption of negligible interactions between the nanoparticles reasonably applicable in this case.

While at small $Q$ the signal is dominated by the coherent contribution due to the translational diffusion of the PEG AuNP, at high $Q$ the signal is dominated by the chain relaxation ($A_{1}\leq 0.2$) and it is composed by similar proportion of coherent and incoherent (Fig. S1). The observation of only one single exponential decay at those high $Q$-values implies that the coherent and incoherent scattering have a similar relaxation behavior within the limit of the sensitivity of the technique. A possible explanation for this feature is that we are looking to a system of polymers closely tethered to a small surface from the perspective of a relatively large length-scales ($>25$ Å). In these condition it is reasonable to think that the polymers cannot move independently from each other and that the self and coordinated motion of the polymers relax on similar timescales.

The second relaxation process is related to the dynamics of the polymer chains tethered to the NP, as shown by the fact the $\tau_{pol}$ resembled the relaxation times found in solution of PEG2000 in D2O (Fig. 9). The relative contribution of the two relaxations processes is expressed by the $Q$-dependent parameter $A_{1}$, that represents the portion of translational diffusion, with the sum rule $A_{1}+A_{2}=1$. As shown in Fig. 2(b), the translational diffusion of the NP dominates at low $Q$ and the polymer dynamics dominates at high $Q$. Interestingly, the trend of $A_{1}$ as a function of Q suggests a maximum at the same Q where the $I(Q)$ obtained by SANS measurements presents a shoulder (see the shaded region in Fig. 2(b)). In fact the form of the function $A_{1}(Q)$ reflects a modulation of the structure factor. At the $Q$ values where the maximum occurs, the incoherent part of the scattering is around 30%. This incoherent part also contribute to the $A_{1}(Q)$ with a term given by $I_{inc}(Q,\infty)/I_{inc}(Q,0)$. This term is the usually called EISF (Elastic Incoherent Structure Factor), which relates to the volume where the hydrogen atoms are confined in, i.e. the nanoparticle corona.

It is useful to plot the inverse of the relaxation times in terms of $\Gamma=1/\tau$. As detailed in Eq. (7) the \Gamma_{int}=\Bigg{(}DQ^{2}+\frac{1}{\tau_{\rm int,\it j}}\Bigg{)} contains also the term $DQ^{2}$ related to the translational diffusion of the whole NP. In Fig. 11 (a) and (b) we report as a function of $Q^{2}$ the $\Gamma_{D}=1/\tau_{D}$ and the $\Gamma_{pol,cor}$, where the latter corresponds to the internal dynamics corrected from the translational diffusion of the whole NP. Both $\Gamma_{D}$ and $\Gamma_{pol,cor}$ have a clear linear trend that intercepts the origin. A dynamic whose relaxation can be described by a simple exponential decay in time, with a relaxation rate $\Gamma\sim Q^{2}$ has a diffusion-like behaviour. The coefficient diffusion associated to the relaxation time of the tethered polymer chains are around $9\pm 2$ Å2ns*-1* and $25\pm 3$ Å2ns*-1* at temperatures of 280 and 318 K, respectively. Such a behaviour is also associated to a linear trend with $t$ of the mean square displacement. This is in qualitative agreement with what observed in coarse grained molecular dynamics simulations Lo Verso et al. (2013) that described the dynamics of brush polymers on spherical nanoparticles. This is an important result since the fact that the dynamics of a corona is represented by a single exponential in this time and length-scale is not due to intrinsic limitation of the technique or of the analysis that are not able to distinguish a higher number of relaxations but it corroborates computational studies. The linear trend for the mean square displacement found in reference Lo Verso et al. (2013) was qualitatively related to a Rouse description of a polymer brush under the boundary condition that constraints the tethered part of the polymer to the surface.Johner (1993) Under this assumption all the modes but one are quickly damped. The only remaining mode has a relaxation behaviour that can effectively be described with a simple exponential. The diffusion constant associated to this dynamics was shown to be twice the diffusion constant of a free chain in solution.Johner (1993) This in fair agreement for example with the center of mass diffusion found for the PEG2000 solution of our study (see below).

If we tried to frame these dynamical results in a Zimm scenario the relaxation times of the polymers would be related to the viscosity, $\eta_{s}$, of the liquid medium by the following relation $\eta_{s,280K}/\eta_{s,318K}=(280\tau_{1}(Q)^{280K}/318\tau_{1}(Q)^{318K})=R_{z}$.Doi and Edwards (1988) The average ratio between the viscosity of the D2O between the two temperatures 280 and 318 K is $\sim$ 2.6 whereas the $R_{z}\sim 2.3$. Even if the trend of the relaxation times with temperature is roughly compatible with what the Zimm model predicts, the absolute value of the viscosity that it is associated with these relaxation times ($\eta_{solvent}=K_{B}T\tau_{pol}Q^{3}/6\pi$)Kanaya et al. (2005) as predicted by the Zimm model is around 0.004 Pa$\cdot$s at 280 K and 0.0016 Pa$\cdot$s at 318K, i.e. roughly two times bigger than expected. A scaling of the relaxation time based only on the solvent viscosity neglects the important contribution to the dynamics related to the fact that the polymers are tethered to the NP surface. This has a significant impact on the resulting dynamics. This fact suggests that our results cannot be fully represented in the context of this common used model.

This $\Gamma\sim Q^{2}$ behaviour was already reported for amphiphilic diblock copolymer micelles in aqueous solution like in the work of Matsuoka et al. Matsuoka et al. (2000) and Castelletto et al..Castelletto et al. (2003) In the latter, however, the $\Gamma$ corresponding to the internal dynamics of the micelle did not intercept the origin. Our results are similar also to those found for polystirene-butadiene diblock copolymer micelles in that, the dynamical processes are ascribable to two relaxation times, one for the translational diffusion of the whole micelle and one accounting for the dynamics of the polymer corona. However, it differs in the fact that a Zimm behaviour was used to describe the latter. In agreement with all the above cited work the dynamics found for the polymer corona in this case is not compatible with the so called breathing mode, a collective mode of the polymer chains driven by a competition between the osmotic compressibility and the entropic force of the polymer chains predicted by de Gennes. de Gennes (1986) Such a dynamics is characterised by a distribution of relaxation times, at variance with our experimental observation that gives evidence of only one relaxation time associated with the polymer chains.

The fitting analysis of the NSE data for the PEG2000 solution dissolved in D2O at 10% in weight illustrates also the presence of two dynamical processes with relaxation times $\tau_{1}$ and $\tau_{2}$ . By a close inspection of Fig. 9(a) and 9(b) we see that the faster relaxation times $\tau_{1}$ measured for the homopolymer were very close to the relaxation times $\tau_{pol}$ obtained from the analysis of the NSE data of the PEG2000 AuNP at each $Q$. This is particularly evident at the highest temperature (Fig. 9(b)). The slower dynamical process $\tau_{2}$ is compatible to the relaxation time associated to the translational diffusion of a PEG2000 polymer in D2O as measured by NMR experiments at similar concentration.Waggoner et al. (1995) A direct comparison is however not obvious since the data in literature are performed up to 8% concentration and the extrapolation in these range of concentrations and temperatures are not trivial for the PEG2000. Therefore we tentatively ascribed this relaxation time to the translational dynamic of the polymer in the solution. We will come back later to the consistency of this assumption. Previous attempts to fit the data corresponding to the PEG2000 homopolymers with a standard fitting procedure showed that these data could be fitted also using a single stretched exponential in the Q range $0.06\div 0.19$ Å*-1*. By letting the $\beta$-parameter free, the fit returned a value close to 0.9. A $\beta$ parameter smaller then one is a sign of a more complex dynamic scenario involving a distribution of relaxation times. Johnston (2006) Indeed, a distribution of relaxation times is very likely to occur in a macromolecular system such as a polymer in solution. As we have seen previously in this section, instead of considering a continuous distribution of relaxation times, one can represent this complexity in a simplified way by describing the dynamical properties of the system with a finite number of components. But which is then the best way to describe our data?

With the aim of finding the model that support our data with the highest statistical significance, we then run the fitting analysis by letting the RJ-MCMC algorithm use a model consisting of a linear combination of $\beta$-stretching exponentials (Eq. (10)) and leaving the algorithm free to estimate all the parameters. In looking for the best likelihood between the model and the data, the algorithm explored all the different possibilities such as sums of simple and/or stretched exponential functions, and clearly also the single stretched exponential function that, as we have seen, can still describe the data.

The prior distribution function for the $\beta_{j}$ (with $j=1\ldots,n$) described in Sec.II.3 was used. Quite interestingly the algorithm converged to the same combination of simple exponentials described in section III.2.1. The results found by the Bayesian algorithm showed quite clearly that the combination of one (only at very small $Q$) or two simple exponentials was indeed the most probable description of ${I}(Q,t)/{I}(Q,0)$. At each Q, the algorithm reproduced within the errors the results obtained for a combination of simple exponentials. The estimated average values for $\beta_{1}$ and $\beta_{2}$ (when present) were both very close to one (typically $0.98/0.99\pm 0.02$) and the posterior probability $P(\beta_{1}=1)$ and $P(\beta_{2}=1)$ were close to 100% and between 85 and 96%, respectively.

More interestingly, the solution with two single exponential functions that picks out a finite number of relaxation processes was preferred over the one with one single stretched exponential that is often used to describe the presence of multiple relaxation processes. We can conclude that the proposed method has allowed a sort of “fine-structure” analysis of the dynamic scenario for the polymer solutions providing also a probabilistic statistical support to the model adopted for the decription of the time correlation function relaxations. We stress that by choosing arbitrarily to fit the data with just one Kohlrausch-Williams-Watts (KWW) function in aprioristic way, part of the information contained in the data might be lost. Moreover in this particular case the choice of an apparent simpler model (one single stretched exponential rather than two exponentials) would not be even justified by some parameters economy criterion since in both cases the number of fitting parameters would be the same.

This analytical compact form is often used to represent NSE data in the contest of the Zimm and the Rouse models,Richter et al. (2005) where the Langevin equation is solved through an expansion in normal coordinates. A relaxation time is associated to each of these normal modes. The KWW results from an average over the normal modes used to solve the Langevin equation. Although a stretched exponential function is ruled out by the algorithm as highly improbable, this by no means implies that dynamic processes found in our analysis are at variance with what it is predicted by polymer theories.

In fact, an interesting experimental result of our analysis of homopolymers solutions, is that the decay time ratio of the two exponential modes $\tau_{1}$ and $\tau_{2}$ found by the application of Eq. (9) has a constant value of about 2 (to within experimental uncertainties), as displayed in Fig. 12.

In the Rouse model Doi and Edwards (1988); Richter et al. (2005) the $p$=0 normal mode corresponds to the center of mass diffusion with a diffusion coefficient $D_{R}$. All other modes ($p\geq 1$) account for the chain internal motions. Since the Bayesian analysis detected only two modes, including the overall polymer translational diffusion that we tentatively ascribed to the $\tau_{2}$ relaxation, we worked out the $I(Q,t)$Richter et al. (2005) in the context of the Rouse model, taking into account only the first internal mode ($p=1$), that has the longest relaxation time. In particular, it can be shown that $I^{\prime}_{int}(Q,t)$, i.e. the time-dependent part of the internal dynamics in Eq.(7), which now reduces to a single exponential, contributes to $I(Q,t)/I(Q,0)$ with a term:

[TABLE]

with

[TABLE]

where $\tau_{R}=N^{2}b{{}^{2}}\zeta/3\pi^{2}k_{B}T$ is the Rouse timeDoi and Edwards (1988), $N=45$ is the the number of monomer units of PEG2000, $b=3.5$ ÅHansen et al. (2003) is the segment length in the Rouse chain and $\zeta$ is the microscopic friction coefficient. Therefore we finally find that $I(Q,t)/I(Q,0)$ is given by the sum of two exponential terms with decay times expressed as $1/D_{R}Q^{2}$ and $1/2D_{R}Q^{2}$, respectively. Hence, in the context of Rouse Model and under the approximation used the ratio between the relaxation times $\tau_{2}$ and $\tau_{1}$ as expressed in Eq (9) is 2, in fair agreement with the experimental finding. Furthermore, this results provide us an internal consistency proof that $\tau_{2}$ can be confidently attributed to the polymer center of mass diffusion and allow us to calculate from their values an average values of $D_{R}\sim 12\text{\AA}^{2}$/ns and $5\text{\AA}^{2}$/ns at 318 and 280K, respectively, and from the relation $\zeta=k_{B}T/ND_{R}$ a ratio $\zeta_{280K}/\zeta_{318K}\sim 2.1$. If again a Rouse context is assumed to model the internal dynamics of the polymer chains in solution, the same ratio between the friction coefficients at the two investigated temperature can be calculated using the estimated relaxation times $\tau_{int,1}$ (see Eq. 7) and they follow the temperature scaling relationship $280\cdot\tau_{int,1}(Q)^{280K}/318\cdot\tau_{int,1}(Q)^{318K}=\zeta_{280K}/\zeta_{318K}\sim 2.1$. This shows a further internal consistency of the rescaling in temperatures of relaxation times according to the Rouse model.

Unlike for the PEG2000 AuNP, for the PEG 400 AuNP, ${I}(Q,t)/{I}(Q,0)$ could be always described only with one term corresponding to the translational dynamics of the NP, without the need of including a term accounting for the polymer dynamics. This could be due to different reasons. The net amount of the polymer, as well as the volume fraction of polymer with respect to that of the whole NP is significantly lower for PEG400 AuNP compared to PEG2000 AuNP. As a matter of fact if we approximate the volume of the NP to $4/3\pi R_{h}^{3}$ and if we consider that the radius of the gold core is $R_{Au}=~{}25$ Å  we have that in the case of the PEG2000 AuNP the hydrated polymer corona occupies $\sim$ 0.90 of the NP volume, whereas only 0.70 in the case of PEG400 AuNP. Moreover, the concentration of the NP in the solution were 5% in weight for the PEG2000 AuNP and only 2% in weight for the PEG400 AuNP. Clearly, the contribution of the polymer corona in the two experiments was significantly different, in particular the scattering intensity due to the regions where polymer was present was much higher in the case of the PEG2000 AuNP than in the case of the PEG400 AuNP. Furthermore, the condition $Q\cdot R_{h}>>1$ that determines the Q range where the polymer dynamics dominates with respect to the translational diffusion is less effectively fulfilled in the case of the PEG400 AuNP due to the significantly smaller value of $R_{h}$. The combination of these two conditions might have obscured the scattering of the polymers in the case of the PEG400 AuNP. The lack of evidence of a polymer dynamics in our NSE experiment could also be due to an effect of confinement. The short PEG400 chains are anchored to the surface of the NP with an area per molecule around 14 Å2.Maccarini et al. (2010) The effect of the anchor on a solid surface combined with the presence of neighbouring polymers highly packed might slow down the fluctuation of the PEG400 to values outside the experimental window of the instrument. Unlike the PEG400 that extends only up to 25 Å away from the surface, the PEG2000 extends up to 60 Å away from the surface allowing a significant part of the molecule to fluctuate in the larger volume available.

The NSE data measured on the 5% solution of PEG400 in D2O  evidenced only one relaxation process in the time window of the experiment. This relaxation corresponds very well with the NMR results concerning the translational diffusion of the polymer at the same concentration and temperatures. Waggoner et al. (1995) A possible explanation for the absence of internal modes might be related to the fact that the relaxation time corresponding to the faster translational diffusion of PEG400 (compared to PEG2000) might be comparable with the relaxation time of the internal dynamics of PEG400. If the two processes occur at a very similar timescale it would be very difficult to distinguish between them.

V Conclusions

We have exploited a Bayesian approach together with the very important property of decomposition of any time correlation function in terms of exponential functions to explore the nanosecond dynamics of polymer coated nanoparticles dissolved in D2O in comparison with the dynamics of equivalent polymers free in solution measured through NSE spectroscopy. After decomposing the intermediate scattering function in elementary components, we associated each component to a particular kind of dynamics (center of mass translational diffusion or internal polymeric dynamics) both for the grafted nanoparticles and for the homopolymer solutions, by comparing the relaxation times obtained from the NSE analysis to those obtained with DLS and NMR data from literature and checking the agreement with polymer science theory. For the PEGAuNP we disentangled the translational dynamics of the whole nanoparticle from that of the grafted polymer chains. The dynamic phenomena observed depended strongly on the molecular weight of the tethered polymer. In the case of longer PEG2000 grafted chains we could clearly see the translational dynamics of the NPs and the internal dynamics of the polymer chains. The polymers did not follow a Zimm dynamics in the range of times and momentum transfer explored, nor the often invoked breathing modes. Instead our analysis indicates that the polymer corona relaxation follows a pure exponential decay, which appears to be an important confirmation of the behavior predicted by coarse grained molecular dynamics simulations Lo Verso et al. (2013) and consistent with a Rouse description of a tethered polymer brush where all the Rouse modes are damped but the longest one related to a movement over the entire polymer chain length. Johner (1993) In the case of shorter PEG400 grafted chains, only the translational dynamics of the NP was visible. This could be due to the significantly lower amount of polymeric material (and the NSE signal associated with it) in this case, or due to a high mechanical constraint of the PEG400 in proximity of the gold surface.

As for polymer solutions, the analysis of the NSE data evidenced two relaxation times in the case of longer PEG2000 polymers associated to the translational dynamics of the polymer and to the $p$=1 Rouse mode, respectively. The proposed analysis estabilishes that beside the center of mass molecule diffusion only one relaxation time is statistically supported and this is then necessarily described as a simple exponential decay. As an important outcome the more common formal description of the time correlation function decay in terms of stretched exponential, as in the averaged form of Zimm and Rouse theory, is excluded by the joint posterior distribution of the adopted model conditional on the observed data. Nevertheless, the internal relaxation time is absolutely consistent with the dominant Rouse time for a polymer free chain and effectively both estimated relaxations times scale with temperature according the Rouse theory.

In the case of the shorter PEG400 only one relaxation time comparable with the translational diffusion of the polymer was evidenced.

Our methodology addresses one of the most important issues in the analysis of inelastic and quasi-elastic neutron scattering data - the arbitrariness in the choice of the number of components needed to describe the studied dynamical processes. In conventional fitting analysis, a particular model is imposed, based on a reasonable physical hypothesis and without any knowledge of a posterior probability that this choice of the model is statistically significant. Hence, the physical hypothesis, although well motivated, affects substantially the result of the fitting analysis. In the specific case of NSE experiments, nISF is often described in terms of stretched exponential that can account for the superposition of more relaxation mechanisms, in terms of an effective relaxation process modulated by a beta parameter. In view of our lack of knowledge about the most reliable model describing the data, no such invasive hypothesis was done to run the RJ-MCMC algorithm used in this work to explain the experimental data. We have shown for example how it is possible to make explicit the richness of the dynamic properties of a system such the PEG2000 solution that would be instead hidden if nISF’s were modeled as a stretching exponential function. At the same time we provided a probabilistic evidence about the model used to describe the experimental data and without sacrificing the simplicity of the model description. In doing so, we obtained that the output of the Bayesian analysis complies fully with the theoretically demonstrated general result on the exponential expansion of any time correlation function,Bellissima et al. (2015, 2017); Guarini et al. (2017) which has been here applied for the first time to the description of the dynamics of a macromolecular system at the nanosecond time scale typically investigated by means of NSE. We believe that the approach presented here will serve as an effective tool that can be used in general to provide an unbiased interpretation of neutron spin echo data.

VI Acknowledgments

The authors acknowledge the Institut Laue Langevin for neutron beam time allocated for the experiments. We acknowledge Forschungszentrum Jülich at the Heinz Maier- Leibnitz Zentrum (MLZ) in Garching for the allocation of beam time used to perform preliminary NSE experiments. In particular we thank Olaf Holderer and Michaela Zamponi for the assistance on the J-NSE instrument. We thank Ralf Schweins for performing the DLS and SANS measurements. We thanks also Simone De Panfilis and Alessandro Paciaroni for very stimulating discussions. We are immensely grateful to Fabrizio Barocchi who provided insight and expertise that greatly strengthened important aspects of this work. This article is dedicated to the memory of Giuseppe Briganti (1951-2015).

VII coherent and incoherent scattering

VIII PEG2000 AuNP

IX PEG400 AuNP

X PEG2000

XI PEG400

XII Complementary information to the Bayesian Analysis

XII.1 *The RJ-MCMC algorithm

The RJ-MCMC sampler we propose performs $M$ sweeps and, at each sweep $m$, all the parameters, including $k$, are updated in turn. This is performed by drawing the new values of a certain parameter conditionally on the data and all the other parameters. The algorithm uses different fixed-dimension moves to update the model parameters, conditionally on a fixed $k$, plus a variable dimension move to update the number of components $k$.

In the fixed dimension moves, all the parameters, but $\nu$, are updated by means of Metropolis-Hastings moves Tierney (1994). For a generic parameter $\theta$ to be updated, the Metropolis-Hastings move proposes, at each sweep $m$, a new candidate value $\theta^{\prime}$, drawn from a proposal distribution q(\mbox{\boldmath\cdot}|\theta), which depends on the current value of $\theta$ as well as from tuning parameters. This candidate value is then accepted with a probability calculated according to an acceptance rule, which ensures that the Markov chain converges to an equilibrium distribution which is the joint posterior of the parameters, and is given by $\min\{1,R\}$, where

[TABLE]

with $\Theta^{\prime}$ denoting the whole parameter vector, in which the parameter $\theta$ has been substituted by $\theta^{\prime}$, $P(y|\Theta)$ being the likelihood and $P(\Theta)$ the prior on the parameters. Notice that the product of the first two ratios is simply the ratio between the posteriors evaluated in $\Theta^{\prime}$ and $\Theta$, respectively, with the normalizing constants cancelling out. The tuning parameters are tuned during some pilot running of the chain in such a way to guarantee an acceptance rate of the Metropolis-Hastings move of approximately $30\%$. To give an example, we describe in detail the updating of the component weights $A$, the Metropolis-Hastings updating of the other parameters being similar. At the $m$-th sweep of the algorithm, in order to update the current parameter $A$, we draw a candidate $A^{\prime}$ from a proposal $q(A^{\prime}|A)$. A reasonable and natural choice of the proposal distribution for $A^{\prime}$ is the Dirichlet distribution with vector of parameters proportional to the current value of $A$, i.e. $\mathcal{D}(\epsilon_{A}A_{1},\ldots,\epsilon_{A}A_{k})$, where $\epsilon_{A}$ is a tuning parameter. Under this choice, the proposed $A^{\prime}$ fulfils the constraints $A_{j}^{\prime}\geq 0,\forall j$, and $\sum_{j=1}^{k}A_{j}^{\prime}=1$, and the distribution of each $A_{j}^{\prime}$ is centred on the current $A_{j}$, with variance controlled by $\epsilon_{A}$. After drawing $A^{\prime}$, the acceptance probability is calculated as $\min\{1,R\}$, where, accordingly to Eq.(S1),

[TABLE]

with the three ratios being, respectively, the likelihood, the prior and the proposal ratios in Eq.(S1), and d(\mbox{\boldmath\cdot};\lambda) representing the density of the $\mathcal{D}(\lambda)$ distribution evaluated in a specific point “$\cdot$”. Finally, to decide whether the move is accepted or not, a random number $u$ is drawn from a Uniform distribution on the interval $[0;1]$ and compared with $\min\{1,R\}$. If $u<\min\{1,R\}$, the move is accepted and $A^{\prime}$ becomes the new parameter vector, otherwise the move is rejected and the parameter vector stays unchanged.

The parameter $\nu$ is instead updated by means of Gibbs sampling. The prior of this parameter is, in fact, conjugated with the likelihood and this results in a closed form for its full conditional distribution (i.e. the posterior distribution of the parameter, given the data and the value of all the other parameters). Thus, a new value of the parameter $\nu$ can be directly draw from its full conditional distribution, i.e.

[TABLE]

Finally, updating the value of $k$ implies a change of dimensionality for the vectors of weights $A$, stretching parameters $\beta$ and relaxation times $\tau$. We accomplish this by introducing a RJ move, consisting in a random choice between the birth of a new component or the death of an existing one. The probabilities of the birth/death alternatives are $b_{k}$ and $d_{k}=1-b_{k}$, respectively, and depend only on the current value of $k$. Of course, $d_{1}=0$ and $b_{k_{\max}}=0$; otherwise we choose $b_{k}=d_{k}=0.5$, for $k=2,\ldots,k_{\max}-1$. The move starts with a random choice between birth and death, using probabilities $b_{k}$ and $d_{k}$. For a birth, we let $k^{\prime}=k^{m}+1$ and we pick a position $j^{\prime}$ at random among $1,\ldots,k^{\prime}$, with probability $1/k^{\prime}$ for the place to be occupied by the new component $j^{\prime}$. Then, a weight for the new component is drawn as $A^{\prime}_{j^{\prime}}\sim\mathcal{B}(1,k)$ and a new vector $A^{\prime}$ is obtained by rescaling the existing weights, using $A_{j}^{\prime}=A_{j}(1-A^{\prime}_{j^{\prime}})$. Also, values of $\tau^{\prime}_{j^{\prime}}$ and $\beta^{\prime}_{j^{\prime}}$ are drawn from their respective priors

[TABLE]

and collocated in the $j^{\prime}-$th place of vectors $\tau$ and $\beta$ to give $\tau^{\prime}$ and $\beta^{\prime}$, respectively. The parameters $\nu$ and $\gamma$ remain unchanged. The acceptance probability for the birth move is calculated as $\min(1,B)$, where, accordingly to the RJ rule

[TABLE]

where the first three factors in the multiplication represent the likelihood ratio (with $f(t_{i})$ and $f^{\prime}(t_{i})$ being the model function evaluated at the current parameter set and the new parameter set, respectively) and the prior ratios, the fourth factor is the proposal ratio (with be(\mbox{\boldmath\cdot};1,k) denoting the beta density evaluated in “$\cdot$”) and the fifth factor is the Jacobian of the transformation from $A$ to $A^{\prime}$. Thus, the acceptance probability for the RJ move is obtained from the Metropolis-Hastings acceptance probability in Eq.(S1) simply multiplying by the Jacobian of the transformation from the old to the new parameters. For a death move, a component is chosen at random among any existing one, the chosen component is deleted and the remaining weights are rescaled to sum to 1. The acceptance probabilities for a death move is $D=\min(1,B^{-1})$, with $B$ given in Eq.(S3).