Statistical Analysis of Phase-Inversion Neutron Specular Reflectivity

4132

Langmuir 2009, 25, 4132-4144

Statistical Analysis of Phase-Inversion Neutron Specular Reflectivity† N. F. Berk*,‡ and C. F. Majkrzak§ Department of Materials Science and Engineering, UniVersity of Maryland, and NIST Center For Neutron Research, National Institute of Standards and Technology, 100 Bureau DriVe, Stop 6102, Gaithersburg, Maryland 20899-6102 ReceiVed August 26, 2008. ReVised Manuscript ReceiVed October 31, 2008 The phase-inversion approach to neutron specular reflection is elucidated in a formal setting, in order to emphasize its conceptual coherence and to facilitate study of some of its statistical properties in the context of real data. An operational notion of data degradation is introduced and illustrated with the randomizing effects of shot noise (“counting” noise) and the systematic “bias” induced by data truncation. Some basic statistical effects of phase-inversion are worked out in the new formalism and illustrated by simulated examples. A principal is advanced that phase-inversion sets the limit of available information from specular reflection.

1. Introduction Neutron specular reflection probes the laterally averaged scattering length density (SLD) depth profile of thin films, including those of biological interest. Specular reflection can be cast formally as a one-dimensional scattering problem, and the power of the technique ultimately stems from its resulting mathematical coherence, establishing a one-to-one correspondence between the reflection amplitude spectrum and the SLD responsible for it. Thus not only is the reflection amplitude unique to its SLD (the “direct” problem), but the SLD is unique to its specular reflection spectrum (the “inverse” problem) for perfect data. Moreover, practical “phase-inversion” methods have been developed for determining the reflection amplitude from one or more reflectivity measurements and then inverting it, providing a means for reliable interpretations of thin film reflection data. In this work we discuss some basic aspects of the statistical analysis of neutron specular reflection phase-inversion, both analytically and with the use of examples from numerical simulations with a realistic model lipid bilayer. Much of the effort here will be given to developing a symbolic operator language for describing the relationship between the direct and inverse problems, especially in the context of dealing with degraded (i.e., actual) data. Such a formalized approach certainly is not required for the problem at hand, but we believe it provides a good organizational foundation for discussing a topic that may be less familiar to some readers than to others. Thus, in sections 2 and 3, we begin with mathematical outlines of the direct and inverse problems for neutron specular reflection, emphasizing the one-to-one correspondence between the veridical (i.e., genuine) SLD, F(x), and the (perfect) reflection amplitude spectrum, r(k), it produces. A formalized notion of measured or degraded r(k) f rˆ(k) is introduced in section 4, and section 5 defines the meaning of inverting a degraded spectrum to obtain a measured F f Fˆ . The idea of reverting an inverted Fˆ (x) to a reflection amplitude rˆ(k) is described in section 5.1. Some basic statistical properties of phase-inversion are then discussed in section 6, and we give a phase-inversion “principle” in section 6.1. Along the way, in section 5.2, we recast the Born approximation in the new operator language, which, we believe, † Part of the Neutron Reflectivity special issue. * Corresponding author. E-mail: [email protected]. ‡ University of Maryland and NIST. § NIST.

helps to illuminate some of its essential properties. Several illustrative graphs are shown in section 7, derived from numerical simulations using a model thin film. In particular, we focus on the distinction between the effects of noise (“statistical”) and truncation (“bias”) degradation. Measurement uncertainty is discussed in terms of estimated prediction error in section 7.3, and single-measurement “statistics” are briefly addressed in section 8. We then give a brief summary and offer some concluding remarks in section 9, including a wave of the hand at some broader-ranged thoughts. The essentials of phase-inversion for neutron specular reflectivity may be found in our earlier work1 and references therein. A standard reference on inverse scattering problems is the book by Chadan and Sabatier.2

2. Direct Problem Given a scattering length density (SLD) profile F(x), where we assume, without much loss of generality, that F(x) ) 0 for x < 0, we can exactly calculate the unique reflection amplitude r(k) it produces for all k on the real k-line. The reflectivity associated with r(k) is R(k) ) |r(k)|2. Getting r(k) directly from F(x) amounts to computing

r(k) )

4π 2ik

∫0∞ eikxF(x)ψ(x, k)dx

(1a)

where ψ(x,k) is the exact physical solution of the Schro¨dinger equation (SE) for a plane wave incident from the left,

-

∂2ψ(x, k) + 4πF(x)ψ(x, k) ) k2ψ(x, k) ∂x2

(1b)

i.e., the solution with boundary conditions, ψ(0,k) ) 1 + r(k) and ∂ψ(x,k)/∂x|x)0 ) ik(1 - r(k)). Equation 1a follows from the Lippmann-Schwinger equation representation of the physical solution of (1b), viz.,

ψ(x, k) ) eikx + 4π

∫0∞ g0(x - x′, k)F(x′)ψ(x′, k)dx′

(1c)

where g0(x) is the free-space Green’s function, satisfying ([∂2/ ∂x2] + k2)g0(x,k) ) δ(x), and δ(x) is the Dirac delta function. For the problem at hand, (1) Majkrzak, C. F.; Berk, N. F.; Perez-Salas, U. A. Langmuir 2003, 19, 7796. (2) Chadan, K.; Sabatier, P. C. InVerse Problems in Quantum Scattering Theory; Springer-Verlag: New York, 1989.

10.1021/la802779r CCC: $40.75  2009 American Chemical Society Published on Web 01/07/2009

Achemso Statistical Analysis

g0(x, k) )

1 -ikx (e θ(-x) + eikxθ(x)) 2ik

Langmuir, Vol. 25, No. 7, 2009 4133

(1d)

where θ(x) is the Heaviside (unit step) function. For x e 0, (1d) in (1c) leads at once to ψ(x,k) ) eikx + r(k)e-ikx, with r(k) given by (1a). Alternatively, (1a) can be derived from (1b) using the Wronskian method,3 which can be generalized to include polarized neutrons.4 In practice eq 1a is rarely used in this form for numerical computation, but it establishes the existence of a rigorous transformation, F f r, which can be represented formally as

r ) ropF

(2)

In (2) the operator rop, acting on a given F, means the procedures summarized in (1) or their functional equivalents, which can be called the class of rop. When it is useful to identify a particular procedure in the class, we can say that it represents rop. We refer to (1) as the fundamental representation of rop. In particular, exact methods exist for computing r for any F of interest from two or three R’s from F in contact with known references.1 Such procedures define what we will call induced representations of rop, nominally different from the fundamental representation in that their actions on F are implicit. Determination of reflection amplitudes from directly observable reflectivities is called the phase problem, the idea being that |r|2 is directly measured, so that the phase, arg r, must somehow be inferred indirectly from the data, if it is desired at all. There is no general solution to the phase problem for reflectivity in precisely these terms. With only occasional exceptions, the “lost” phase is not retrievable from diffuse |r|, absent convenient circumstances, such as a F exhibiting high symmetry or some other special shape constraint. The induced representations of rop just mentioned simultaneously determine |r| and arg r pointwise at each k and are not limited in any essential way to particular symmetry or shape assumptions about F. For our purposes, we can treat the phase problem as solved and thus, for now, need not distinguish among representations of the rop class. The ultraconcise (2) can be written more explicitly as

r(k) ) rop[F(x)]

(3)

to clarify that rop takes functions of x into functions of k. However, for brevity, normally we will leave out the arguments of r and F. Notice that rop is not a linear operator; i.e.,

rop(F1 + F2) * ropF1 + ropF2

(4)

This follows from (1), since F appears both explicitly in the integral and implicitly in the wave function. Indeed, rop is linear only in the Born approximation (BA), defined in (1) by setting BA F ) 0 in (1b) but not in (1a). The linearity of rop is useful for making some points, as we will see later, since the BA is asymptotically exact as |k| f ∞, and it can serve as a good approximation (except at small k) for weakly scattering films. The BA is discussed below in section 5.2. Equation 2 expresses what is usually called the direct problem in scattering theory; calculate r from a given F. The complementary inVerse problem specifies that F be exactly retrieved from the r it produced. The direct problem is “easy” in the sense already given: for any (reasonable) F there are both formal and feasible methods of solution. The inverse problem is “hard” because, to begin with, a solution does not always exist, and (3) Majkrzak, C. F.; Ankner, J. F.; Berk N. F.; Gibbs, D. In Magnetic Multilayers; Bennet, L. H., Watson, R., Eds.; Elsevier: New York, 2005; p 299. (4) Majkrzak, C.; Berk, N. F.; O’Donovan, K. V. In Neutron Scattering From Magnetic Materials; Chatterji, T. K., Ed.; Elsevier: New York, 2005; p 397. (5) Sacks, P. E. WaVe Motion 1993, 18, 21.

furthermore, when it does, its implementations tends to be more technically involved than direct problem methods.

3. Inverse Problem Let us assume, to begin with, that the inverse problem has a unique solution; i.e., there exists an exact procedure for converting r(k) into the veridical F(x). In an obvious notation, we can express this as

F ) rop-1r

(5)

F ) rop-1r ) rop-1ropF ≡ F

(6a)

r ) ropF ) roprop-1r ≡ r

(6b)

Thus,

and

as appropriate for an inverse. In fact (unique) inverses exist only within a certain class of F: viz., those F(x) which do not have bound states (in the usual quantum mechanical sense) and which satisfy some relatively weak (in physical terms) regularity conditions. We will call this the class of invertible Fsnotwithstanding that actually it is r that is invertedsand continue to use the unmodified symbol F for it. The class F includes all nonnegative F of finite thickness L, i.e., F satisfying F(x) g 0 for all x, but also includes functions with negative portions that are “too weak” to support bound states. The importance of this class is that it supports the logical relation F T r, i.e., it puts F and r into one-to-one correspondence. For invertible F, a procedural definition of rop-1, analogous to that for rop in (1), exists in the Gel’fand-Levitan-Marchenko (GLM) method.2,5,6 Namely, for given r first obtain its Fourier transform (FT)

G(x) )

1 2π

∫-∞∞ e-ikxr(k) dk

(7a)

with G(x) ) 0 for x e 0; next, with the input G(x), solve the GLM integral equation for the kernel K(x,y),

K(x, y) + G(x + y) +

∫-yx K(x, y)G(z + y) dz ) 0

(7b)

and then retrieve F from the diagonal, K(x,x), according to

F(x) )

1 dK(x, x) 2π dx

(7c)

For real-valued SLDs, r(-k) ) r*(k), so the integration in (7a) does not require explicit measurement at negative k. The small imaginary parts of most neutron SLDs do not cause observable problems, but the described inversion formalism does not apply to strong absorbers, which rules out certain materials in neutron reflection and effectively prevents extending phase-inversion, as defined here, to X-ray reflection. The form of the GLM equation given in (7b) applies only to cases where F(x) ) 0 for x < 0, i.e., where films of interest have their leading edges placed at the origin, but this is general enough for practical use. The integral equation in step (7b) can be replaced by a partial differential equation,5 which is better suited to numerical work. We will refer to this modification as GLMS. Also, the FT in (7a) can be reduced to either of the real integrals (6) Berk, N. F.; Majkrzak, C. F. J. Phys. Soc. Jpn. 1996, 65(Suppl. A), 107.

4134 Langmuir, Vol. 25, No. 7, 2009

∫-∞∞ cos(kx) Re r(k) dk 2 ∞ ) - ∫-∞ sin(kx) Im r(k) dk π

G(x) )

Berk and Majkrzak

2 π

(8)

indicating that either Re r or Im r is sufficient data to invert r.5,6 This property follows from the analytic (holomorphic) behavior of r(k) in the upper half of the complex k-plane, which applies to invertible F. We could, therefore, sharpen the notation by distinguishing between versions of rop-1 that work separately on Re r or Im r, but we do not need to do so for the purposes at hand. (The differences between Re r and Im r can be important in implementations, however. Some experimental setups are better suited to obtaining Re r rather than Im r, while, on the other hand, Im r is more sensitive to certain features of F, especially in some thin films of biological interest.)

4. Degradation Operators We defined r ) ropF as the ideal reflection amplitude “caused” by (i.e., calculated from) the veridical F. We refer to this as the veridical r. The veridical r(k), therefore, can be calculated exactly at all k. Indeed the infinite spectrum, r(k) for all 0 e |k| e ∞, is required in the definition of rop-1, and it is only in this sense that (5) is true. (It is sufficient to restrict k to positive values, since Re r(k) is even in k, and Im r(k) is odd.) The “r” we measure, however, is a an imperfect image of the veridical r, degraded by counting noise (“shot” noise), data truncation (i.e., the restriction of r(k) to a finite domain of measurable k), instrumental resolution, and other possible intrusions. In particular, for solutions of the phase problem, the degradation of r is inherited from the underlying degraded reflectivity measurements, depending on the methods used. We will see a specific case of this in section 7, where we show simulated examples. Now with any giVen circumstances leading to degraded rscall the result rˆswe may associate an operator rôp having the property

rˆ ) rôpF

(9)

for any F. In proposing rôp we assume that the phase determination methods devised for ideal data can be applied directly to real data. Since several methods of phase-sensitive reflectometry have been developed, we could affix indices or other notations to rˆ and rôp in order to distinguish among them. For the sake of simplicity, however, we leave the methodological dependence implicit and further assume that any of them correctly produces r in an ideal experiment. One such method directly determines Re r only, but we saw above that this is sufficient information for inversion. So, for ease of notation, we can let r stand for Re r, Im r, or Re r + iIm r, as appropriate. Alternatively, (9) is

rˆ ) Dopr

(10)

Dop ) rôprop-1

(11)

where we will call

the degradation operator for the given measurement setup. The operator rop-1 appears in (11) because F can be formally represented by rop-1r, as in (5), even though r is not directly observable. Thus it is Dop that permits us to speak of arriving at rˆ by way of r. (We return to this point shortly.) In the absence of degradation we have rôp ≡ rop, and then, Dop ≡ roprop-1 ) 1op. From (11) and (6),

rôp ) Doprop

(12)

Dop can be factored into explicit degradation components. For example, say rôp determines rˆ locally in k; i.e., it determines

rˆ(k) at each k only from measurements at the given k. Then the operation of truncating a theoretically infinite data set, say to 0 e k e K, can be applied at any step of the procedure used K to determine rˆ. In other words, if we let T op represent k-truncation, leaving the actual truncation limit implicit, we are able to write K K Dop ) D′opTop ) Top Dop

where Dop ′ represents the “rest” of Dop. Formally, by the action K Top f(k) )

{

(13) T opK

is defined

f(k) for 0 e k e K 0 otherwise

(14)

∞ for any f(k). Thus if we let T op stand for ”truncation” to the ∞ ) 1op. It is important to (semi) infinite k-domain, we have T op bear in mind that in real experiments, k-truncation is not optional: every feasible measurement is restricted to a finite k-domain. One strives, however, to make K as large as possible while minimizing other sources of degradation, such as background and noise. Typically, K is determined after the data are in hand and usually is smaller than the maximum k of the data collected. If Dop is a non-local methodsi.e., one requiring integration over all ksfactorization still is possible but noncommutative; we must truncate first, because the measurements necessarily have a finite domain of k, and this immediately limits the domain of rˆ. In such cases only the first equality of (13) is valid. The very notion of a Dop (or an rôp) presupposes that we are able to identify all of the degrading components of a measurement. For the purposes of simulating data and performing statistical analysis on the results, rôp can be suitably modeled; the computational penalties of such calculations can be incorporated into Dop. With limited and mostly academic exceptions, data degradation is irreversible. In other words, rôp does not have a useful inverse which would allow us to write F ) rôp-1rˆ. For example, consider K , so that rˆ ) the case of truncation-only degradation, Dop ) T op K r. While the analyticity of r formally admits the possibility T op of analytic continuation of this rˆ to the entire k-axisseffectively K with 1opsthe problem usually is so ill-behaved as replacing T op to be unsolvable. The compounding of truncation with resolution broadening and noise only worsens the situation. It follows that ˆ op-1 ) roprôp-1 also does not exist. D 4.1. Degradation of R Before leaving this topic, it may be helpful to distinguish between degradations of the reflection amplitude r and the reflectivity R ) |r|2. In the latter case, we observe a degraded image Rˆ of the ideal reflectivity R, for which we also may introduce a degradation operator Dop such that

Rˆ ) DopR

(15)

Indeed Dop (or “type-D” degradation) likely is what the reader has in mind when we speak of Dop, since degradations due to noise and instrumental resolution are directly observed in Rˆ. Such instrumental effects are naturally expressed as “doing something” to the ideal R, as formalized by (15). For example, a convolution integral associated with instrument resolution is directly represented this way, as is shot noise. If these are the only degradations of concern and we approximate the counting uncertainties by (presumably small) additive noise, we can write

Rˆ(k) ) Dop[R(k′)] )

∫-∞∞ K(k, k′)R(k′)) dk′ + R(k)

(16)

where K(k,k′) is the resolution function (this use of K not to be confused with the real-space GLM kernel in (7b)) and R(k) is the noise. On the other hand, except for k-truncation, degradations of r are harder to describe succinctly, because on standard


Langmuir, Vol. 25, No. 7, 2009 4135

instruments r is reconstructed from the results of one or more reflectivity measurements. Thus degradations of type D are inherited indirectly by rˆ as they propagate through the reconstruction. As suggested earlier, it may not be obvious that rˆ can be viewed as the result of doing something to r, i.e., that a Dop exists such that rˆ ) Dopr. And, as we saw, this is rigorously true only for the class of invertible F, where F ) rop-1r. Dop, however, is definable without referring to the underlying F. It also can be pointed out that generally DopR * |Dopr|2.

Fˆ (x) ) (rop-1rˆ)(x) ≈ F(x) + P (x)

where P ) and calling the inversion result Fˆ . R(k) is far-removed from white noise but nevertheless has a broad “frequency” spectrum P (x). Frequency here lies in the domain of the variable conjugate to k, namely, along the x-axis; low and high frequencies refer to small and large x, respectively. Specifically, we identity the physically relevant low-frequency band of measurement noise R(k) with the support of F(x), viz., x ∈ (0,L), and introduce the corresponding x-truncation operator.

5. Inverting Degraded Data We briefly summarize the discussion so far. The veridical F induces the veridical (or ideal) r, according to the prescription denoted by r ) ropF. In a feasible experiment we measure rˆ, a degraded image of r, which, over the class of invertible F, can be represented by rˆ ) Dopr, where Dop encodes the deleterious effects of the measurement and reconstruction process. For this class, there is an exact solution of the inverse reflection problem, which we represent by rop-1, such that rop-1r ) F retrieves the veridical F from ideal data. Without degradation (Dop ) 1op), rop-1rˆ ) rop-1r ) F. Now it seems reasonable that if rˆ ≈ r in some sense, then we can expect rop-1rˆ ≈ F in a similar sense. There are potential problems with this line of thought, of course. To begin with, it may not be clear in what “sense” rˆ needs to approximate r in order for rop-1rˆ to approximate F to like degree. Data truncation, i.e., the T Kop componentofDop,isagoodsandarguablythemostimportants case. Take Dop ) T opK , so that K rˆ(k) ) Top r)

{

r(k) for 0 e k e K 0 otherwise

(17)

Speaking loosely, we might say that the rˆ in hand, i.e., up to K, is perfect, but this could be misleading, since rop-1 receives all of rˆ, including the “truncation zeros” beyond K and there are no a priori means for accurately accessing the relevant size of the discrepancy between this rˆ and r. We do know, however, that truncation may significantly degrade the output of rop-1rˆ compared to F. The effects of truncation are easily described in the Born approximation, which we will discuss later, but for now let it suffice to recall that the first step in the GLM representation of rop-1 in (7a) is a Fourier transform of the input; so for truncated data, inversion of rˆ begins with a low-pass filter. From a formal perspective, the derivation of the GLM algorithm itself rests on several “good” properties of the veridical reflection amplitude that do not survive truncationsor shot noise, for that matter. Thus, in the narrowest sense, rop-1 is wrong for degraded data. This is an extreme viewpoint; in practice, the GLM and related inversion methods appear to be well-behaved when applied to measured and simulated rˆ. Before we adapt rop-1 to the inversion of rˆ, however, it is useful to consider a situation complementary to truncation-only degradation, viz., noise-only degradation. Thus let Dop ) Nop, where Nop represents the result of propagating the shot noise in the original reflectivity measurements through the reconstruction leading to rˆ ) Nopr. For the sake of argument, let us assume that rop and rop-1 are linear operators, as they are in the Born approximation, and assume further that the noise is additive. Then we can decompose rˆ into its veridical signal and noise component, r(k) and R(k), respectively,

rˆ(k) ≈ r(k) + R(k) and invert this to get

(18)

(19)

rop-1R

L Top g(x) )

{

g(x) for 0 e x e L 0 otherwise

(20)

for any g(x). In terms of (20), L is the smallest L satisfying L F(x) ) F(x). Then T op

Fˆ (x) ) (rop-1rˆ)(x) ≈ F(x) + PL(x) + PH(x)

(21)

with L PL ) Top P ) rop-1RL

(22a)

L PH ) (1op - Top )P ) r-1 op RH

(22b)

and

L In (22), the second equalities define RL ) ropT op P and RH ) L )P, which identify low- and high-frequency rop(1op - T op components of the k-space noise R(k)si.e., the noise in r(k) that generates the high- and low-frequency (out-of-film and in-film) real space noise. RL(k) thus can be viewed as analogous to bandlimited noise from the x-space f k-space perspective. These notions will be used in section 5.1. Only the in-film PL(x) is relevant to accessing the uncertainty in Fˆ , and PH(x) can be discarded. It could be argued that PH(x) produces uncertainty about the support of Fˆ , but the effective support of PH(x) can be much larger than L, a clear artifact of measurement noise. Rounding effects of k-truncation produce the main uncertainty in locating a sharp film edge. Moreover, L usually is accurately known from analyzing well defined Kiessig fringes in the neutron data or from separate X-ray measurements. For films with soft boundaries, where the support may be significantly in doubt, L can be treated as a parameter and results discussed over a range of values. Therefore, we arrive at

L Fˆ ) Top rop-1rˆ(x) ≈ F(x) + PL(x)

(23)

Applying k-truncation induces smearing (convolution broadening) in the x-domain, requiring us to replace (23) with L K L ˜ Fˆ ) Top rop-1rˆ(x) ≈ Top F˜ (x) + Top P (x)

(24)

where F˜ and P˜ represent smearedsi.e., reduced resolutions images of F and P. Some of these effects are more easily described in the BA, as discussed below in section 5.2. These arguments do not rest strongly on the assumption of linearity, and with this in mind, we define the physically relevant inversion of rˆ by the operator L Fop ) Top rop-1

(25)

Fˆ ) Foprˆ

(26)

such that It is tempting to refer to Fop as rˆ-1, but this would be inappropriate. As we saw earlier, when Dop * 1op, rˆ is not invertible, and because of TopL , neither is Fop. By making constituent factors of Fop and rˆ explicit, we are able to refer back to the veridical F and rewrite (26) as

4136 Langmuir, Vol. 25, No. 7, 2009

Berk and Majkrzak

Fˆ ) SopF

(27)

L Sop ) Top rop-1Doprop

(28)

where

is the F-degradation operator. In the absence of r-degradation, Dop ) 1op and then TopL ) 1op in effect, because TopL F ) F. This gives Sop ) 1op, and so Fˆ ) F as expected. However, since the composition of Sop includes noninvertible factors, there is no Sop-1; therefore, we cannot invert (27). In other wordssand unsurprisinglysin the presence of r-degradation, Fˆ is insufficient information for the exact reconstruction of the veridical F. We will recast some of these results explicitly in the BA in section 5.2. 5.1. Reversion. One way of judging the utility of Fˆ is to compare the reflection it produces with the measured reflection rˆ. We call this reVerting Fˆ to a reflection amplitude rˆˆ, according to K rˆˆ ) Top ropFˆ

ψ(x, k) ≈ ψBA(x, k) ≡ eikx

rBA(k) )

K K L sop ) Top ropFop ) Top ropTop rop-1

(31)

It is natural to compare the reverted and measured spectra via the residual

∆ ) rˆ - rˆˆ ) (1op - sop)rˆ ≡ δoprˆ

(32)

Using the linear approximations we discussed near the end of section 5, we easily find for noise-only degradation that

rˆˆ ≈ r + RL

(33a)

∆ ≈ RH

(33b)

and

Therefore sop acts as a low-pass filter on rˆ, transmitting only the information relevant to the support of F, while δop acts as the complementary high-pass filter, identifying the disposable noise. With k-truncation restored, the symbolic generalization of (33) is cumbersome: for we must then consider reversions of F˜ and TopL P˜ and devise notations for those results. Since our goal remains chiefly pedagogical, we sacrifice some precision for simplicity and let r˜ be the reversion of F˜ and RL stand for the reversion L ˜ of T op P . Then, (33) becomes

rˆˆ ≈ r˜ + RL ) r + RL + BIAS

4π 2ik

∫-∞∞ e2ikxF(x) dx

(36)

The integral in (36) is a Fourier transform, and thus if it exists (i.e., is finite everywhere in k), its inverse also exists (Plancherel’s theorem), so that

F(x) )

(30)

where

(35)

i.e., the exact wave function in the film is replaced with the incident wave. With this, we immediately get

(29)

That is, we treat Fˆ as the veridical F, computing its ideal reflection amplitude and truncate the resulting spectrum to match the k-domain of rˆ. (As in our usage of “inversion” and “inverted”, we will call rˆˆ a “reverted” r.) By chaining operators, reversion also can be viewed as a direct transformation of the data, viz.,

rˆˆ ) soprˆ

Let us briefly return to noise-only degradation. We see from (33a) and (18) that rˆˆ effectively is a smoothed version of rˆ, i.e., rˆ without its “rough component,” RH. There are, of course, many techniques available for smoothing data. Smoothing by reversion has the advantage over most of being nonparametric and essentially automatic. There remains the question of how to use a reverted r quantitatively. We return to this in section 6.1. 5.2. Born Approximation. The Born approximation (BA) (or kinematical theory) is well-known in reflectometrysindeed in scattering, in generalsfor the manner in which it “reduces” the direct and inverse problems to Fourier analysis. It is worthwhile, therefore, to work out the operator language in the BA, which we do now. The BA represents first-order perturbation theory and can be derived from (1a) with the approximation

i π

dk ∫-∞∞ e-2ikxkrBA(k) 2π

(37)

In standard definitions of FTsssuch that as in (7a)sthe integral in (36) represents a mapping from x to 2k ) q as conjugate variable, where q is the wavevector transfer in specular reflection. However, we observe the usual practice in dynamical theory and express rBA as a function of k. The validity of the BA rests on r(k) being sufficiently small to justify treating the film as weakly scattering. Because of the factor k-1 in (36), this justification generally breaks down as k f 0, even when F itself is small. On the other hand, the BA becomes exact as k f ∞ because limkf∞ r(k) ) 0. (This limit does not depend on the prefactor but is a general result of Fourier integrals for admissible F.) It is important to realize that the left-hand side of (37) is the veridical F only if the r under the integral is rBA, as written. In typical applications of (37), the measured rˆ(k) is put into the integral, so that the result is Fˆ BA(x), the Born approximation of Fˆ (x). Equation 37 also can be expressed as ∞ -2ikx BA dk e r (k) ∫x F(x′) dx′ ) -1 ∫ -∞ 2π 2π

(38)

which now involves the inverse FT of rBA(k) rather than krBA(k). These formulas also can be economically represented by operator notation. First, introduce the x-derivative and antiderivative operators

dop(·) ) ∂x(·)

(34a)

(39a)

and

and

∆ ≈ r - r˜ + R - RL ) RH + BIAS ˜

for 0 e k < K, where ˜

BIAS ) r˜ - r + RL - RL

dop-1(·) )

(34b)

(34c)

Therefore, the reversion-filters sop and δop also transmit bias induced by k-truncation and, more generally, by other degradations.

∫x dx(·)

(39b)

which satisfy dopdop-1 ) dop-1dop ) 1op. Next, let the relevant FT and FT-1 pair be represented by

Fop(·) ) and

∫-∞∞ e2ikx(·) dx

(40)


Langmuir, Vol. 25, No. 7, 2009 4137

Fop-1(·) ) 2

dk ∫-∞∞ e-2ikx(·) 2π

(41)

with Fop-1Fop ) FopFop-1 ) 1op, as expected. Then, (38) takes the form

-1 -1 BA Fop r 8π2

dop-1F )

(42)

Subsequent application of dop gives

F)

-1 dopFop-1rBA 8π2

(43)

-1 d F -1 2 op op 8π

(44)

and thus

rop-1 )

in the BA. Inverting (44) then gives

rop ) -8π2Fopdop-1

(45)

which re-expresses the inverse of (38). In words, to compute r from F in the BA, take the FT of the antiderivative of F. The seeming awkwardness of this formulation replaces the “nuisance” of the k- 1 prefactor in (36). In particular, the divergence (infinity) of rBA(k) as k f 0 immediately evident in (36) now stems from the constancy of ∫x-∞ F(x′) dx′ as x f ∞, for compactly supported F. Since Fop is a linear operator, it follows at once from (45) that rop is a linear operator in the BA, unlike its exact behavior shown in (4). We also note that (43) is equivalent to the GLM inversion in (7) if K(x,x) in (7c) is replaced by -G(2x), as achieved by dropping the integral term from (7b). One may loosely view the solution of the GLM equation (7b) as “stripping away” the dynamical scattering from r. The effect of k-truncation is easy to work out in the BA, since only Fourier transforms are involved. Thus consider the case of truncation-only degradation. From (27) and (28) Fˆ ) SopF, with L L K Sop ) Top rop-1Doprop ) Top rop-1Top rop

(46)

which in the BA, using (43) and (44), is L K Sop ) Top dopUop dop-1

(47)

where we have introduced K UKop ) Fop-1Top Fop

(48)

From the definitions (40) and (41), it is straightforward to obtain, for any f(x),

UKop f(x) )

∫-∞∞ UK(x - x′)f(x′) dx′

(49)

∫-∞∞ UK(x - x′)F(x′) dx′

UK(x) )

K 2Kx sinc π π

(50)

sin πx πx

(51)

and where

x)

is the “sine cardinal” function. Moreover, since ∂xUK(x - x′) ) -∂x′UK(x - x′), it follows from integration by parts in (49) that K K K ) Uop dopUop dop, i.e., that dop commutes with Uop . Therefore, from (47), L K Sop ) T LopdopUKopdop-1 ) Top Uop

(52)

(53)

for Dop ) TopK . In particular, from the behavior of the sinc function L in (50), UK(x - x′) f δ(x - x′) as K f ∞, and then Fˆ f Top F ) F. The result in (53) is unsurprising but still worth comment. It is, of course, what would be expected for the convolutionbroadening of any pulse transmitted by a flat, low-pass filter. On the other hand, the derivation here is not trivial because of the fact that in the specular reflection problem it is r and dop-1F that comprise the FT pair in the BA, not r and F directly; the K resolves this technical distinction. commutivity of dop and Uop K Furthermore Sop in (52) actually is a filter bank, consisting of Uop L L and Top in series. The sharp cutoff that Top effects in the x-domain introduces “distortions” of its own near the nominal edges of F(x) when the smearing induced by TopK is significant, but these are consistent with information loss in the x-domain in the presence of truncation in the k-domain. Finally, although the BA analysis of truncation-only degradation does not apply strictly to exact inversion, it serves as an intuitive guide and, equally importantly, works reasonably well for K not too small. The effects of TopK are more difficult to summarize in the presence of noise and other degradations, but to first approximation, additive degradations are convolution-smeared in similar fashion to (53). In particular, further distortions near the film edges may result from truncation-induced leakage of the otherwise disposable PH(x) into the support of F.

6. Statistical Analysis A phase-inversion measurement of the function F by neutron specular reflection yields an estimate of it, viz., Fˆ ) SopF. M repeated measurements of the same F, i.e., repeated measurements of the same physical specimen, thus yield a set of estimates,

Fˆ m ) Sop,mF

(54)

for m ) 1,2,...,M, where Sop,m is the real-space degradation operator of the mth measurement. For a given setup, {Sop,m}M depends on m only by way of its chance elements, i.e., on the noisy components of Dop,m. A measurement result Fˆ m can be likened to the outcome of a trial on a chance setup. Normally with such outcomes, one would expect that the greater the discrepancy between the measured result and the veridical “measurand,” the smaller the chance of it being observed. In the present case, however, we first need to separate the two dominant components of Dop,m, viz., truncation-only and noise-only degradation, before we can assert this. For even in the absence of randomness, the inversion of rˆ ) TopK r can not yield F if K is not sufficiently large. Thus let us define K Fˆ T ) FopTop r

with

so in the BA

L Fˆ (x) ) Top

(55)

for the truncation-only case. Relative to the random component of Dop,m, Fˆ T is a surrogate for the veridical F, absent a means of “undoing” the truncation; therefore, we call Fˆ T the T-Veridical F. The discrepancy between Fˆ T and F is a “systematic error”. In traditional statistical arguments, it is theoretically identifiable as the result of averaging the outcomes of an infinite number of trials on the same chance setup, assuming that the noise averages to zero (usually taken as axiomatic). We take this view of the discrepancy between Fˆ T and F. Then it is sensible to assert that the greater the discrepancy between Fˆ m and Fˆ T, the less chance of it being observed. In conventional uncertainty analysis, it is often assumed that systematic errors, once identified, are removable. In our case, this requires retrieving F from Fˆ T, an ill-conditioned problem, as we discussed in section 5.2. Thus

4138 Langmuir, Vol. 25, No. 7, 2009

Berk and Majkrzak

truncation bias can not be “corrected,” other than by increasing Ksi.e., changing the setup. This viewpoint could be considered unnecessarily pessimistic. After all, if we guessed F, say Fg, and were lucky enough to guess correctly, then outcomes of M simulations {Sop,mFg}M would be statistically indistinguishable from M measurements {Fˆ m}M on the setup. In other words, {Fˆ m}M does “tell us something” about F: namely, that any Fg close enough to F is supported by the data. The obvious problem with guessing F is that we would be wrong almost always, and a great many bad guesses also would be statistically indistinguishable from {Fˆ m}M because of the loss of K . Of course, “guessing” F generally information induced by T op means fitting of some sort, and we will return to this topic later on. In particular, we will need to examine whether for any measured rˆ, a fit can reliably produce a better estimate of F than the inversion Fˆ can. In the subsequent discussion we use the following definitions. For a set {fm}M, we denote the arithmetical (“sample”) average by M

∑

jf ) 1 f M m)1 m

(56)

where fm is fm(k) or fm(x), as appropriate. For functions f(k)s“rtype” functionsswe define

〈f〉 )

1 K

∫0K ω(|k|)(|k|)f(k) dk

(57a)

where K is the maximum value of k measured. The weight function ω(k) compensates for the rapid decrease with increasing k of functions like r(k) to produce a relatively uniform dynamic range over the integration; normally we assume ω(k) ) k2. For functions f(x)s“F-type” functionsswe use

〈f〉 )

1 L

∫0L f(x) dx

(57b)

where L is the thickness of veridical F. K and L are specified by the setup. For r-type or for F-type functions, we measure the separation of fm from fn by the Euclidean norm, or “loss” function,

L[fm, fn] ) 〈|fm - fn|2 〉

(58)

which we also denote by Lm,n when the context is clear. 6.1. Phase-Inversion Principle. We propose the following: Principle. In a set of M measurements on a setup, {rˆm}M, no estimate of the Veridical F accounts for any rˆm significantly better than one of the inVersions Fˆ m′, for some 1 e m′ e M, or the arithmetic aVerage F ˆ. Corollary. In a single measurement on the setup, say rˆ1, no estimate of the Veridical F accounts for rˆ1 significantly better than the inVersion Fˆ 1. The gist of the principle is that the chance setup defines a pool of statistically equivalent {Fˆ 1,Fˆ 2,...}. The pool is large but far from arbitrary: every member derives from the veridical F, is K and is systematically degraded by the same data truncation T op degraded by the same noise mechanism. Such lack of arbitrariness in phase-inversion, essentially the absence of undetermined parameters, underlies the principle, regardless of the representations of rop and rop-1. A particular rˆ in the pool, say rˆ1, may be “best” accounted for by its inversion Fˆ 1 or by another member of the pool, say Fˆ 100. Both Fˆ 1 and Fˆ 100 will suffer from the same degree of truncation degradation and will differ only as a result of chance effects. Thus, which of the corresponding reversions, rˆˆ1 or rˆˆ100, is closer to rˆ1 with respect to a given measure, is

largely a matter of chance. More formally, assume that M is large and write

rˆm(k) ) rˆ(k) + δrˆm(k)

(59)

The residual δrˆm(k) is self-defined by (59), so that rˆm(k) ) 0. Also

〈δrˆm 〉 ≈ 0

(60)

The ≈ sign is necessary because the noise in rˆm(k) is neither strictly additive nor “iid” (independently and identically distributed) over k; however, the approximation generally holds well for large enough K. Let us compare reversion rˆˆm to measurand rˆn using the M × M loss matrix

Lm,n ) 〈|rˆˆm - rˆn|2 〉

(61)

Then using (59),

Lm,n )〈|rˆˆm - rˆ|2 〉 + 〈|δrˆn|2 〉 - 2Re〈(rˆˆm - rˆ)*δrˆn 〉 (62) ≈〈|rˆˆ - rˆ|2 〉 + 〈|δrˆ |2 〉 m

n

The first line is exact, while the second results from asserting 〈(rˆˆm - rˆ)*δrˆn〉 ≈ (rˆˆm - rˆ)*〈δrˆn〉 and using (60). The factorization is justified for large M, where rˆ (k) becomes effectively noiseless and smoothly varying on the k-scale of δrˆm(k), which is “all” noise. (Somewhat more abstractly, rˆˆm(k) - rˆ(k) and δrˆm(k) essentially reside in orthogonal manifolds of the function space they share, since they are mostly characterized by low frequency and high frequency noise, respectively.) For a setup with low to moderate noise levels and for large M, we can expect 〈|δrˆn|2〉 to be relatively independent of nsi.e., roughly speaking, the rˆm are equally noisy. Moreover, as discussed in section 7, any rˆˆm is less noisy than the rˆ from which it is derived, since the high frequency noise in the latter is filtered out of the former. Thus, except for very low noise levels, the second line of (62) is

Lm,n ≈ 〈|rˆˆm - rˆ|2 〉 + 〈|δrˆ1|2 〉

(63)

for all n. Therefore, Lm,n is smallest for the reversions rˆˆm that are closest to rˆm. This means, in turn, that most rˆm will tend to lie closer to some rˆˆm′ ) TopK ropFˆ m′ other than their “own” rˆˆm ) TopK ropFˆ m; moreover, most may lie closest to rˆˆm* ) TopK ropFˆ m*, the m ) m* for which 〈|rˆˆm* - rˆ|2〉 is smallest. Extensive simulations reveal these loss differences to be small, however, and we do not consider them statistically useful. In effect, they show that all Fˆ for the setup are statistically equivalent. We think it is easy to accept that members of a statistical pool {Fˆ m} are equivalent among themselves. The principle, however, makes the stronger assertion that no other estimates of veridical Fsi.e., estimates outside the poolsdo significantly better in terms of accounting for a measurement rˆ. The reason is the unavoidable degradation TopK . Because of data truncation at k ) K, a measurement can only reveal structure compatible with the “trend” scale of spatial resolution π/(2K) (i.e., π/Qmax). “Detail” on significantly finer scales entails data at k > K, which are not available. Precise notions of trend and detail in Fˆ can be defined with wavelet analysis.8 Indeed, the data are as compatible with Fˆ T, the T-veridical F, as with F itself. Thus agreement between rˆ and the rmodel calculated from a Fmodel may as good as the agreement generated by an inversion pool member, but we contend that it is unlikely to be significantly superior. (7) Majkrzak, C. F.; Berk, N. F. Phys. ReV. 1998, B 58, 15416. (8) Berk, N. F.; Majkrzak, C. F. Langmuir 2003, 19, 7881.


We do not attribute metaphysical content to an infinite pool of measurements {rˆ1, rˆ2,...}, or to the corresponding inversions {Fˆ 1, Fˆ 2,...}, emanating from a given phase-inversion setup. In our naive view, any {Fˆ m}M speaks for itself, and we do not need to envision it as a sampling of a predefined (or “premeasured”) universal set. Indeed, since we are dealing here with functional data and resultssrˆ and Fˆ are functions, rˆ(k) and Fˆ (x)sthere is no measurable universal set. No two independently measured countable pools {rˆ1, rˆ2,...} and {rˆ1′, rˆ2′,...} are likely to coincide or even contain pointwise common members. A given {Fˆ m}M can be thought of as a subset of an an infinite pool obtained by endlessly adding measurements to it, in the manner of a von Mises “collective,” but we draw no specific implications from this kind of imagining, such as providing long run frequencies. Nevertheless, we can speak of the “true” F common to all such pools, viz., the veridical F which spawns them in the given setup. It is plausible to assert that any pool of measurands is as good as any other of similar size for statistical purposes, at least for estimating measurement uncertainties. Moreover, given the definition of the phase-inversion setup, we may expect the emergence of stable statistical behavior in {Fˆ m}M as M increases without limit, since such behavior can be rationalized as being a property of the particular setup. In practice, of course, M cannot be especially large and often is as small as M ) 1. However, it is possible to simulate large-M sets of measurements more-orless realistically, for sensible models of F, and from these draw useful inferences about the statistical behavior of the setup. As we will see, using a variation of the bootstrap method, such information can be applied to estimating uncertainties of individual measurements. The single measurement corollary to the phase-inversion principle follows from the general statement with the supplement that we can not use information we do not have. That is, without a Fˆ m*1 in hand, we cannot infer that its reversion rˆˆm+1 lies closer to rˆ1 than rˆˆ1 does, even if we admit this possibility in the event. Nonetheless, if all {Fˆ 1, Fˆ 2,...} for the measurement setup are statistically equivalent, it is unlikely that any estimate of F would be significantly better than Fˆ 1.

7. Examples We now show examples of some of the statistical properties of phase-inversion reflectometry described in the previous sections. These are all based on slightly idealized simulations using a surround Variation7 induced representation of rop. Let us assume the film (i.e., F) under study adheres to a substrate of known composition. Well-characterized silicon or sapphire blocks typically are used for such purposes. It should be said that the film “of interest” often represents one layer of material on a sample composed of several layers, as we will see, so there is always the possible issue of what parts of the sample are accurately “known” prior to the measurement. For surround variation, we assume that only the substrate (in this instance) is known, and that everything on the substrate is “unknown,” even if we have (or believe we have) reliable prior knowledge of some components (a gold-chromium layer, for example.) Let us also say that the neutron beam is incident on the film side of the sample in contact with air. In the language of surround variation, we then speak of the substrate as the backing medium and the incident medium as the fronting medium. The surround variation method has been worked out for arbitrary frontings and backings (e.g., the supporting substrate could be the fronting medium), but for these simulations we use the simplest such configuration. The induced representation of rop for this case, such that r ) ropF, is defined by the following equations, valid at every k:

Langmuir, Vol. 25, No. 7, 2009 4139

Re r )

β-R β+R+2

(64a)

where

R)-

b22Σ1 - b12Σ2 b12 - b22

(64b)

and

β)-

Σ1 - Σ2 b12 - b22

(64c)

and where

1 + Rj 1 - Rj

Σj ) 2bj

(64d)

with

bj(k) )

1-

4πFj k2

(64e)

for j ) 1, 2. In (64d), R1 and R2 are the observable reflectivities for F in contact with backings F1 and F2, respectively. This method directly gives Re r(k), which, as mentioned in section 4 is sufficient information for inversion. (Im r ) (2(Rβ - 1)1/2/(β + R + 2).1 For perfect data, the branch of the square root can be inferred from the known general behavior of Im r near k ) 0, but we will not need to make use of this here.) The surround-variation (SV) rop induced by (64), is exact for perfect data. For imperfect (i.e., actual) data, degradation of rop to Doprop, such that rˆ ) Dopr, is induced by applying the same procedures to measured R1 and R2. This process is straightforward to simulate on a computer. For a model F, we calculate R1 ) |r1|2 and R2 ) |r2|2 by applying the fundamental representation of rop to the setup described above. Since R1 and R2 are produced by samples with different substrates, R1 ) R2 ) 1 for k below the larger of the two critical wavevectors, kc,max; therefore, the SVinduced rop can not be applied at these small k-values. For thin films, however, it is easy to interpolate Re r, with reasonable accuracy, below kc,max to its invariant value at the origin, Re r(0) ) -1se.g., by fitting the available small-k behavior using a simple barrier model for the film. More interestingly, theoretically, we have found that the effect of “filling in” the missing data is well-simulated (for sufficiently thin films) simply by computing G(x) in (7a) with the data at hand and then adding a constant to enforce the G(0) ) 0 requirement.9 For the purpose of these simulations, we simply use the known Re r for k < kc,max. Shot noise in R1(k) and R2(k) is realistically simulated by random numbers P(I0R(k))/I0 from a Poisson distribution P(µ) of mean µ, where I0 is the monitor count for the incident beam, taken to be the same in both measurements. Data truncation is trivial, of course, since we are free to choose any range of k. We assume here that the reflectivity measurements have been corrected for background. (It is rare that reflectivity data are not so corrected.) We also do not consider degraded k-resolution here, since usually it is a minor concern in contemporary neutron specular reflectometry. The simulations shown here were made for the hypothetical F shown in Figure 1, which is a reasonably faithful rendition10 of a hybrid bilayer membrane composed of a (deuterated) lipid ˆ (x) is computed with missing data (9) Berk, N. F. unpublished. That is, if G ˆ (x) - G ˆ (0). below kc take G(x) ) G (10) Majkrzak, C. F.; Berk, N. F.; Krueger, S.; Dura, J. A.; Tarek, M.; Tobias, D.; Silin, V.; Meuse, C. W.; Woodward, J.; Plant, A. L. Biophys. J. 2000, 79, 3330.

4140 Langmuir, Vol. 25, No. 7, 2009

Figure 1. Model for veridical F, inverted Fˆ from simulated “truncationonly” surround variation experiment (Qmax ) 2.0 Å-1), and 〈Fˆ 〉M, the average over M experiments, for M ) 100 and I0 ) 104. The terms Fˆ and 〈Fˆ 〉M are virtually indistinguishable on the scale of the figure. The small biases relative to F, visible in some portions, are attributed to implementation effects. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

monolayer (the two peaks to the right of x ) 90 Å) on an alkanethiol leaf (the dip to slightly negative SLD values centered near x ) 80 Å), which has been deposited on a chromium-gold layer (to the left of x ) 65 Å). This, the supporting substrate, which disappears in the application of rop-1, lies along the negative x-axis and is not shown. The support of F (thickness of the model film) is L ) 120 Å. For computations, F(x) was rendered as as NL ) 128 bins of width ∆x ) L/NL ) 0.938 Å. Simulated measurements were made at I0 ) 500, 103, 2 × 103, 5 × 103, 104, 2 × 104, 40 × 104, and 80 × 104 counts. For each I0, “experiments” were performed to generate data sets {rˆm}M, usually for M ) 100, and these were inverted to produce {Fˆ m}M, etc. In the figures, we plot spectra as a function of wavevector transfer Q ) 2k, in keeping with common practice, and show Q2Re r(Q), which displays more uniformly than Re r(Q) over a large Q-range. We specify data truncation cutoffs by Qmax ) K for the truncation operator 2K but continue using the symbol Top to avoid confusion with earlier formulas. 7.1. Noise-Only Degradation. Figure 1 shows the veridical K ropF for noiseless F (the model) and the inversion Fˆ ) rop-1T op data (equivalent to I0 ) ∞) using the SV representation of rop and the GLMS representation of rop-1, with Qmax ) 2.0 Å-1. Fˆ 2.0 ≈ F to good accuracy, overall, indicating that T op ≈ 1 for this model and that both the SV-induced representation of rop and the GLMS representation of rop-1 can be used with reasonable confidence as computational tools. The figure also shows the average over {Fˆ m}100 for I0 ) 104. Fˆ and 〈Fˆ 〉M are nearly indistinguishable on the scale of the figure over most of the support of F. Small biases relative to F in some portions are partly do to data truncation, since they decrease slowly with increasing Qmax (tested out to Qmax ) 3.0 Å-1), but our numerical implementations (particularly of the computationally involved rop-1) likely also contribute to them. We say more about the behavior of 〈Fˆ 〉M in the discussion of Figure 5. Figure 2 shows instances of R1 ) RSil and R2 ) RSap for F on silicon and sapphire substrates, respectively, for I0 ) 104. Simulated data such as these were used to effect SV-induced Doprop. Figure 3 shows {Re rˆm}M, where rˆ ) Dm,opropF, for a run of K Nop, and M ) 100 SV experiments at I0 ) 104, with Dop ) T op -1 Qmax ) 2.0 Å . The thickened white line is 〈Re rˆ〉100, and the dotted line is the veridical r ) ropF out to Qmax ) 2.0 Å-1. To

Berk and Majkrzak

Figure 2. Simulated instances of R1 ) RSil and R2 ) RSap for F on silicon and sapphire substrates, respectively, for I0 ) 104. The unit of the abscissa is Å-1.

Figure 3. {Re rˆm}100, where rˆ ) Dm,opropF, for I0 ) 104 and Qmax ) 2.0 Å-1: (thickened white line) 〈Re rˆ〉100; (dotted line) veridical r ) ropF for Q e Qmax. The units of the abscissa and ordinate are Å-1 and Å-2, respectively.

Figure 4. {Fˆ m}100 for conditions of Figure 2, also showing (2 × ESTDs. The units of the abscissa and ordinate are Å and and 10-6 Å-2, respectively.

good approximation, 〈Re rˆ〉100 ≈ r, so it seems apparent that 〈Re rˆ〉M f r as M f ∞. Figure 4 shows {Fˆ m}100 for a run of M ) 100 samples with I0 ) 104 and Qmax ) 2.0 Å-1. Display conventions are similar to Figure 2. Different runs reproduce effectively the same result. The central white line is 〈Fˆ 〉100, which for all such runs essentially coincides with F, as in Figure 1. The outer dark lines show 〈Fˆ 〉100 ( 2σˆ 100, where σˆ is the estimated standard deviation (ESTD)


Langmuir, Vol. 25, No. 7, 2009 4141

about the estimated mean, 〈Fˆ 〉100. This band contains almost all of {Fˆ m}100 pointwise (rather than 95%, as for a normal distribution), suggesting that that the random function Fˆ (x) is pointwise more narrowly distributed than normal. Similar behavior is exhibited for other monitor counts. Since σˆ M(x) is relatively insensitive to x for M ) 100, or larger, it is useful to characterize the noise-induced variation of Fˆ by the integrated measure NL

σ j ) L-1

∫0L σˆ M(x) dx ≈ NL-1∑ σˆ M(xn)

(65)

n)1

The quantity σ j 2 is similar to the “object dissimilarity” statistic used in statistical learning.11 Specifically, if we assert Fˆ mn as the mth measurement of “attribute” nsviz., the amount of SLD located at xnsthen the object (here Fˆ ) dissimilarity is defined as

j ) NL-1 D

NL

∑ djn

(66a)

Figure 5. (a) 〈Fˆ 〉M for M ) 1, 2,..., 100: (white line) 〈Fˆ 〉10; (dotted line) 〈Fˆ 〉100. (b) ESTDMM for {Fˆ m(x)}M, averaged over x. See text in section 7.1. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

n)1

where, for Euclidean measure of attribute dissimilarity,

djn ) N-2

N

N

∑ ∑

(Fˆ mn - Fˆ m′n)2 ) 2Var[Fˆ (xn)] (66b)

m)1 m′)1

Therefore,

j ) 2σ D j2

(67)

The elementary idea, of course, is that less statistical variation implies more similarity (less dissimilarity) among possible outcomes. Such notions can be generalized to include nonuniform attribute weighting (shape-dependent effects) and different measures of attribute dissimilarities.11 Figure 5 collectively shows 〈Fˆ 〉M for M ) 1, 2,..., 100. For each M, {Fˆ m}M is a randomly chosen subset, without repetition, of a fixed {Fˆ m}100. The white line is 〈Fˆ 〉10 and the dotted line is 〈Fˆ 〉100. It is noteworthy how rapidly 〈Fˆ 〉M f 〈Fˆ 〉100; even 〈Fˆ 〉10 wellapproximates 〈Fˆ 〉100. The inset shows that the estimated standard deviation of the mean for {Fˆ m(x)}M, averaged over x, which we call ESTDMM, decreases with increasing M substantially faster than M-1/2. It is no surprise, perhaps, that the “random variables” Fˆ m(x) are not pointwise independently distributed, as this behavior indicates, but the implication of this for 〈Fˆ 〉M is remarkable, nonetheless. Figure 6 presents {Fˆ m}100 for several monitor counts I0. Each Fˆ m is drawn as a dotted line. For each I0, the pooled ∆Fˆ m ) Fˆ m - 〈Fˆ 〉M, i.e., {∆Fˆ mn}M×NL (with M ) 100 and NL ) 128), are distributed approximately normally. This is illustrated in Figure 7, which shows a quantile plot (“probability plot”) of data quantiles for I0 ) 5 × 103 versus normal quantiles for STD ) 0.2. The straight line has unit slope. Uniform and triangular distributions do not account for the observed quantiles nearly as well. From such simulations, we obtain

σ j ≈ 15/√I0

(68)

for I0 g 500. This behavior depends only weakly on data truncation for not-too-small Qmax. Since, as we mentioned above, σˆ (x) is relatively independent of x, we can take σˆ (x) ≈ σ j . The numerator of (68) depends, of course, on the veridical F. Since the noise in Fˆ ultimately derives from surround variation shot noise, it is not surprising that σ j scales as I0-1/2. (11) Hastie, T.; Tibshirani, R.; Friedman, J. The Elements of Statistical Learning; Springer-Verlag: New York, 2001.

Figure 6. {Fˆ m}100 for indicated I0. See section 7.1 for display convention. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

Figure 7. Quantile plot of ESTDs for {Fˆ m}100, for I0 ) 5 × 103, pooled over x, versus a normal distribution N(0,0.2). The abscissa and ordinate are dimensionless.

Figure 8 depicts the separation of reciprocal-space shot noise into its direct-space in-film (“low-frequency”) and out-of-film (“high-frequency”) components, PL(x) and PH(x) (shown only for x > L), respectively, as in (21). The example is for Qmax ) 0.88 and thus also includes a small amount of truncation degradation. Figure 9 shows the “smoothing” effect in reciprocal space of throwing away PH(x). The gray lines are {Fˆ m}M for M ) 100,

4142 Langmuir, Vol. 25, No. 7, 2009

Figure 8. In-film (“low-frequency”) and out-of-film (“high-frequency”) shot noise. See section 7.1. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

Figure 9. Smoothing effect in reciprocal space of throwing away PH(x). Gray lines are {Fˆ m}M for M ) 100; dark lines are 10 randomly selected reversions, {rˆˆm}10. The units of the abscissa and ordinate are Å-1 and Å-2, respectively.

Berk and Majkrzak

Figure 11. Reversion of Fˆ (Figure 9) for Qmax ) 0.25 Å-1 and L ) 120 Å. The units of the abscissa and ordinate are Å-1 and Å-2, respectively.

recognizable. The alkane-thiol leaf bound to the gold is discernible, but the deuterated lipid leaf remains a single, albeit distorted, peak. When Qmax ≈ 0.5 Å, the major components of the film (gold, alkane, and lipid) are resolved, but the detailed structure in the lipid remains undefined; and, noticeable truncation “wiggles” distort the hydrogenated alkane. For Qmax ≈ 1.0 Å, the detailed structures in the lipid head-group and tail are revealed, and the alkane SLD is nearly flat. For Qmax ≈ 2.0 Å, Fˆ and F are nearly indistinguishable everywhere. Figure 11 shows the reversion of Fˆ ) Fopr to rˆˆ for Qmax ) 0.25 -1 Å and L ) 120 Å. (See Figure 9 for Fˆ .) Clearly, rˆˆ(Q) accounts very well for r(Q) at Q e K, as one would hope. This same rˆˆ(Q) fails for Q > K, but that also is to be expected, since Fˆ ) L K rop-1T op r “knows” nothing of r(Q) when Q > K. That is to T op say, we know nothing of r(Q) at Q > K, and Fˆ inherits our ignorance of absent data. This state of affairs underlies the formulation of the PIP in section 6.1. Any “improvements” to Fˆ made outside phase-inversion can only lead to “better agreement” with absent data. It is appropriate to describe the non-random effects of truncation-degradation on Fˆ as a bias and introduce the pooled truncation-bias statistic

σK2 )

1 L

∫0L (Fˆ K(x) - F(x))2 dx

(69)

where Fˆ K is the inversion result for truncation-only degradation. (Since we have switched the wavevector variable from k to Q, we now can think of K as standing for wavevector “cutoff.” From the results in Figure 9, we find

σK ≈ e-1.6Qmax

Figure 10. Fˆ (truncation-only degradation) for various Qmax. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

and the dark lines are 10 randomly selected reversions, {rˆˆm}10. The comparison makes clear that much of the noise in rˆ is irrelevant. 7.2. Truncation-Only Degradation. We briefly describe K . Figure 10 shows Fˆ for several effects of data truncation, T op values of Qmax. For Qmax < 0.125 Å, Fˆ is rendered merely as two peaks: one for the CrAu layer and the other for the entire lipid complex. For Qmax ≈ 0.25 Å, the CrAu layer is rounded but

(70)

for 0.25 Å e Qmax e 1.5 Å. 7.3. Estimated Prediction Error. A global statistic for gauging agreement with predicted values is the estimated prediction error (EPE), which here can be defined as11

EPE2 )

1 L

∫0L 〈(Fˆ m(x) - F(x))2〉M dx

(71)

if we view veridical F as “predicting” the M measured quantities {Fˆ m}M at every x. In fact, if F(x) is viewed as adjustable, then EPE(F) is minimized pointwise by F(x) ) 〈Fˆ m(x)〉M; i.e., variance is minimized by the estimated mean “conditioned” on x.11 Here, however, we take F as fixed and the EPE as a measure of its uncertainty. By adding and subtracting


Langmuir, Vol. 25, No. 7, 2009 4143

Fj ≡ 〈Fˆ m 〉M

(72)

from the difference in (71), and using (65) and (70), we get

∫

1 L 1 2 EPE2 ) L 0 〈(Fˆ m(x) - Fj(x)) 〉M dx + L )σ j 2 + σK2

∫

L

0

(Fj(x) - F(x))2 dx

≈σ j S2 + σK2 (73) for large M, where we let S denote “shot-noise-only.” In (73), we have made use of the fact that 〈Fm〉M converges exponentially rapidly to 〈Fm〉∞, as seen in Figure 5. That example assumed noise-only degradation, but the rate of convergence with M essentially is independent of Qmax; therefore, Fj “is” the truncationonly result for Fˆ . In the last line of (73), we set σ j)σ j S, since, as mentioned earlier, σ j is insensitive to truncation over a significant range of Qmax. In effect, noise and truncation bias are independent degradations. Thus, with (68) and (70),

EPE(I0, Qmax) ≈ √a ⁄ I0 + e-bQmax

(74)

where for our example, (dimensionless) a ≈ 225.0 and b ≈ 3.2 Å-1. In general, we expect a and b to depend on the representation of rop-1 and on the shape and support of F. However, we also expect that the functional form of the EPE in (74) is somewhat robust over useful ranges of I0 and Qmax. With these properties in mind, and assuming that the measurements {Fˆ m}M actually are at hand for sufficiently large M, we suggest that a reasonable reporting of the “result” is

Fˆ (x) ) Fj(x) ( kσˆ

(75)

where

σˆ ≡ ESTDMM2 + EPE2 ≈σ j S2 + σK2

(76)

In metrology,12,13 kσˆ is called the combined standard uncertainty with coVerage factor k. Typically k ) 1 in scientific literature, but k ) 2 often is recommended in standards contexts. Equation 76 treats σˆ as independent of x. In eq 76, ESTDMM accounts for the error in the convergence of Fj, which is negligible for large M, while EPE pertains to the overt variation or dissimilarity over {Fˆ m}M and to the bias induced by truncation-only degradation, here treated as “uncorrectable” error.12 The exact meaning of (75) is open to debate. The simplest view along traditional lines is that, with k ) 2, roughly 95% of the results in {Fˆ m(x)}M fall completely within the indicated interval. More importantly, however, we have seen that, except for truncation bias, Fj(x) tends to veridical F(x) with increasing I0 and thus deserves consideration as the best possible estimate of F, under the circumstances. Knowledge of σK in (76) requires knowledge of the veridical F, which is not available from {Fˆ m}M alone. However, meaningful estimates of σK could be obtained by comparing Fˆ (x) with models incorporating some degree of “prior knowledge” about F(x). This is important, since in typical contemporary reflection measurements, the EPE is likely to be dominated by σK. For example, using (74) on our demonstration F, and with Qmax ) 0.3 Å-1, one (12) U.S. Guide to the Expression of Uncertainty in Measurement; ANSI/ NCSL Z540-2-1997, National Conference of Standards Laboratories: Boulder, CO, 1997. (13) Taylor, N., Kuyatt C. E. Guidelines for EValuating and Expressing Uncertainty of NIST Measurement Uncertainties of NIST Measurement Results; NIST Technical Note 1297, National Institute of Standards and Technology: Gaithersburg, MD, 1994.

finds EPE(5 × 103, 0.3) ) 0.65 Å-2, while σK ) EPE(∞, 0.3) ) 0.62 Å-2. Before moving on, we note that in the realm of statistical learning, eg., see ref 11, our {Fˆ m(x)}M would be viewed as training data for predicting F(x). The best prediction minimizes EPE(F), which, for Euclidean metric, leads to Fˆ (x) ) Fj(x).11 In conventional model fitting, the training data consists of a single measurement, say Fˆ 1 (recall that in the context of phase-inversion, we denote the result of inversion as the degraded “estimate” Fˆ of veridical F.) Then the ”best” parametrized model F(x|θ), over a parameter set θ of dimension dθ, could be determined, given F1, by minimizing the logarithm of the likelihood function (LLF),

LLF(θ) )

1 L

∫0L (Fˆ 1(x) - F(x|θ))2 dx

(77)

with respect to θ. Equation 77 looks similar to (71) but lacks averaging over a multimember training set. Since LLF g 0, clearly one cannot do better, “errorwise,” than F(x|θˆ ) ) Fˆ 1(x). Of course, in practice, model fitting usually is done in the wavevector domain. In our terms this means “reverting” F(x|θ) to a “predicted” rˆ(Q|θ) ) ropF(x|θ) and then minimizing

LLF(θ) )

∫0Q

max

(rˆ1(Q) - rˆ(Q|θ))2 dQ

(78)

with respect to θ. This will lead to F(x|θˆ ) * Fˆ 1(x), since Fˆ 1 is “model-independent.” On the other hand, if F(x|θ) is sufficiently complex (i.e., dθ is sufficiently high), there likely exists some F(x|θ) ≈ Fˆ 1(x); but generally, this will be true over a significant range of θ-sets, and unregulated least-squares minimization will then almost certainly “over-fit” noisy rˆ1(Q). In the model-selection problem, such a state of affairs is called the bias-Variance tradeoff:11 roughly said, low-dθ (high modeling bias) f poor fits; high-dθ (low modeling bias) f too many (mostly wrong) fitssunless a regularization scheme is employed. We will not pursue this topic further here, except to note that fitting Fˆ 1(x) directly in the x-domain is useful if one seeks to interpret the inverted Fˆ 1 in terms of a parametrized physical model. Then, in terms of phase-inversion, Fˆ 1(x) may be viewed as (the best) surrogate data.

8. Single-Measurement Statistics The value in having {Fˆ m}M is obvious. First, since 〈Fm〉M converges rapidly to 〈Fm〉∞, much of the noise degradation in {Fˆ m}M can be “averaged out” even for modest M ≈ 10 and I0 not too small. Second, {Fˆ m}M provides the means of quantifying uncertainties using statistical methods of choice. Indeed, simply plotting {Fˆ m}M, as in Figures 4 and 6, qualitatively conveys all the essential information about how well Fˆ is determined from the measurements. A sufficiently large set of measurements enables more sophisticated statistical treatments, as well. For example, {Fˆ m}M can be extended to even larger sets with bootstrapping techniques. And from a Bayesian perspective, {Fˆ m}M provides a means of predicting futuresi.e., as yet unmeasuredsoccurances of Fˆ . Our own view of this is that a posterior probability for Fˆ is useful only if everyone concerned measures exactly the same veridical Fsor, at least, similarly prepared films, as in a round-robin standarization effort. But, at present, most interest in reflection from thin films of biological interest has other goals, we believe. In practice, we expect M ) 1 in the laboratory; so, in the absense of other information, the question is how to embed a single measurement, Fˆ 1, in a statistically useful population. Several answers are possible, but the one we put forward is to treat Fˆ 1 ≡ τ as veridical F and simulate the experiment S times, where

4144 Langmuir, Vol. 25, No. 7, 2009

Figure 12. (white lines) Example {Fˆ m}100; light and dark gray lines: {τs}100 for two “surrogate” veridical F randomly selected from the white set. See section. The units of the abscissa and ordinate are Å and 10-6 Å-2, respectively.

S is as large as we like. That is, generate {τˆ s}S and treat this as a surrogate for {Fˆ m}M)S. We give and example of this in Figure 12 for our model F, which was “measured” 100 times with I0 ) 104 and Qmax ) 0.88 Å-1. In the figure, {Fˆ m}100 is shown in white. Two particular Fˆ m values were selected randomly from this set to act as surrogates for F, and 100 measurements were simulated for each choice, giving two examples of {τs}100. Each of these is shown in Figure 12 in black and gray. The white set was plotted first, then the black, followed by the gray, and the lines were drawn with high transparancy to allow views of the partially occluded black and white sets. This rendering somewhat obscures the differences between the black and gray sets, but that is partly the point to be made; each gives a reasonably good depiction of the variation in the white set. Quantitatively, to two significant j S(black) ) 0.21 Å-2, and figures, σ j S(white) ) 0.15 Å-2, while σ -2 σ j S(gray) ) 0.18 Å . All 100 choices of τ from this {Fˆ m}100 gave (over 100 simulated measurements) 〈σ j S〉 ≈ 0.18 Å-2, with an -2 ESTSD ≈ 0.01 Å and with maximum and minimum values j S}100 σ j S ≈ 0.21 and 0.15 Å-2, respectively. (The distribution of {σ was approximately normal, as determined by a quantile plot.) It j S(F), since τ is a degraded surrogate is not surprizing that σ j S(τ) g σ for F. Thus, in general we can expect this approach to singlemeasurement statistics to result in more cautious assesements of uncertainty. Truncation bias must be estimated separately, using (70)sfor Qmax ) 0.88 Å-1, σK ≈ 0.24sor by modeling.

9. Conclusion In this work we have developed a formal languagesa sort of relational calculussfor describing neutron phase-inversion spectroscopy and have used it to elucidate some basic statistical properties of phase inversion, including notions of “disposable” noise and reversion introduced here. In section 6.1, we have asserted a belief, which we call the phase-inversion principle, that phase-inversion measurements and their parameter-free analysis determine the limits of reliable information about veridical F by specular reflection.

Berk and Majkrzak

Section 7 summarizes simulation studies of a realistic model of a biological thin film. Some details in the results certainly depend on the shape of F(x) and to a limited extent on the surrounding variation method of phase determination. But, the model exhibits features common to many such films, and we feel that the statistical behaviours of the examples are largely representative of phase-inversion in general. The most important results for inverted F are the following: (1) Noise and data truncation degradations propagate approximately independently. (2) Propagated shot (counting) noise is distributed approximately normally and scales, as expected, with monitor count I0 as I0-1/2. (3) Over M measurements, the average SLD, 〈Fˆ (x)〉M converges pointwise rapidlysfaster than O(M-1/2)s to Fˆ K(x), the veridical F(x) degraded only data truncation (Fˆ K f F as Qmax ) 2K f ∞). That is, in phase-inversion one can say, “a little averaging goes a long way.” Indeed, in section 7.3, we argued from a different perspective that over a measured population, {Fˆ m}M, no estimate lies significantly closer to Fˆ K than 〈Fˆ m〉M. We also suggested a format for reporting results from a set of measurements that more-or-less follows conventional guidelines for expressing measurement uncertainty. Finally, in section 8, we showed that useful estimates of noise characteristics of a single measurement, Fˆ 1, can be deduced by “bootstrapping” an ensemble of measurements using Fˆ 1 as a surrogate for veridical F. Roughly similar effects could be achieved by “resampling” the original data sets used in the representation of Fop. While our notion of data degradation introduced in section 4.1 is general, we have here limited our study to shot noise and truncation, viz., Dop ) TopK Nop. A more complete treatment of data degradation in phase-inversion certainly would also include effects of background and limited Q-resolution. In this work, we have sidestepped philosophical issues underlying traditional and contemporary statistical analysis. The assertion of veridical F, measurable in principle, seems to be a frequentist notion, to be sure, and without it, we believe, one could hardly speak of phase-inversion as we have done here. After all, something produces {Fˆ m}M, and Fˆ seq 75sconverges to that “something” as Dop f 1op. This can be contrasted with model-based approaches to data analysis. Even allowing for ”phase-information”sideally, possession of true r(Q), itselfswe can not predict with certainty that a given model F(x|θˆ ) is able to fit r(Q) perfectly. (A straight line cannot fit a parabola perfectly.) Thus, there is no guarantee that F(x) ∈ {F(x|θˆ c)}C for a chosen set of models, c ) 1, 2,..., C, or, even if it is true, that actually finding the right model is feasible. The implications of phaseinversion for model-based regression methodssincluding modern ideas of model selection and model averaging11,14sis an interesting subject for subsequent study. LA802779R (14) Hjort, N. L. J. Am. Stat. Assoc. 2003, 98, 21.

Statistical Analysis of Phase-Inversion Neutron Specular Reflectivity

Recommend Documents