© F.-X. Désert 2022

1 Introduction

When modelling the interstellar dust emission it is usually assumed that its spectrum is close to that produced by a grey body. This assumption has been used over the last 50 yr. Usually, the model has been sufficient to explain the dust spectral energy distribution (SED) along a given line-of-sight, be it in the diffuse interstellar medium (ISM; e.g. Hensley & Draine 2021), in molecular clouds, in star-forming regions, or in galaxies (Galliano et al. 2018). This SED is determined by a few broadband photometric measurements with rather large uncertainties (statistical but also systematic), so to the first order a single temperature is enough for the fits to be satisfactory. The intensity, I_ν, is thus modelled by $$I_{\nu }={\tau }_0{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }{B}_{\nu }(T_{\text{d}}),$$(1)

where β is the so-called emissivity index, τ₀ is proportional to the dust column density, and B_ν is the Planck function at frequency ν and temperature T: $$B_{\nu }(T)=\frac{2h{\nu }^3}{c^2}\frac{1}{\exp (\frac{h\nu }{kT})-1},$$(2)

where h is the Planck constant, c the light velocity, and k is the Boltzmann constant.

In orderto have a more physically motivated modelling of dust emission, we can assume that the main dust component has different temperatures along the line-of-sight. Studies have shown that the interstellar radiation field (ISRF), as seen by a grain of dust, can befluctuating due to shadowing or screening effects. That grain integrates an absorption of the ISRF over 4π steradians; therefore, the average temperature should be quite stable, though a fluctuation along that average is unavoidable.

For example, a single slab of attenuating dust placed in front of the ISRF will produce grains at decreasing temperature T_d = T₀e^{−τ∕(4+β)} as we enter the slab with opacity τ. The mass of grains at a given temperature is thus $\frac{\text{d}m}{\text{d}{T}_{}}=\frac{4+\beta }{T_{\text{d}}}\frac{\text{d}m}{\text{d}\tau }$, an inverse power law of the temperature between the standard value, T₀, and the value at the centre of the slab (hence a linear law for the temperature logarithm).

Moreover, grains of different sizes will have slightly different temperatures. While moving forward in complexity, we must not forget that the photometric data for one pixel on the sky are scarce (e.g. nine values for observations by the Planck mission in intensity). Therefore, we must find a flexible modelling without too many parameters. It is also apparent that in many cases a single-temperature and emissivity index model leads to false anti-correlations between these parameters (Shetty et al. 2009).

In this paper we investigate how to cope with the variability in the dust temperature along the line-of-sight using Planck function derivatives. Then we show how to take this effect through a photometric mapping instrument. Considerations of the cosmic microwave background (CMB) and the cosmic infrared background (CIB) are made in this study.

2 The model

The simplest way to deal with a small temperature range is to assume the dust temperature is distributed around a central value, T_m. For reasons discussed below, we prefer to deal with the logarithm variable t ≡ ln T (the normalisation is irrelevant in practice as t will only be used in a differential way). We can think of the probability distribution as a normalised Gaussian function around a mean value, t_m = ln T_m, with a dispersion, σ, for a given line-of-sight: $$p(t)\text{d}t=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{{(t-t_{\text{m}})}^2}{2{\sigma }^2}\right)\text{d}t.$$(3)

A single dimensionless parameter, σ, encapsulates the logarithm temperature spread, and we assume σ ≪ 1. The Gaussian assumption maximises the lack of information on the source of the spread, but it will be shown that what matters is σ as the second-order moment of the distribution. The question becomes about the emission spectrum of such a superposition of modified blackbodies. In the following subsections, we develop the computation of such a spectrum using the derivatives of the Planck function.

2.1 Planck derivatives

In order to simplify the computation, we define the dimensionless reduced frequency as $$x\equiv \frac{h\nu }{kT}=\frac{h\nu }{k{e}^t}$$

and the composite constant as $P\equiv \frac{2h}{c^2}$.

As such, one can reformulate Eq. (2) as $$B_{\nu }(T)=P{\nu }^{3}f(x),$$(4)

with $$f(x)=\frac{1}{e^x-1},$$(5)

which is known as the photon occupation number.

One can then compute the first and second derivatives of the Planck function with the temperature logarithm as $$\frac{\text{d}{B}_{\nu }}{\text{d}t}=P{\nu }^{3}g(x),$$(6)

with $$g(x)=-\frac{\text{d}f}{\text{d}\ln x}=\frac{x{e}^x}{{(e^x-1)}^2},$$(7)

and $$\frac{{\text{d}}^2{B}_{\nu }}{\text{d}{t}^2}=P{\nu }^{3}g(x)(h(x)-1),$$(8)

with $$\frac{{\text{d}}^{2}f}{\text{d}\ln x^2}=g(x)(h(x)-1),$$(9)

and $$h(x)=x\coth \frac{x}{2}=x\frac{e^x+1}{e^x-1}.$$(10)

The Planck derivative functions are displayed in Appendix A.

2.2 Emission spectrum from a Gaussian temperature distribution

We now proceed to combine the temperature spread, as written in Eq. (3), with the Planck function to get the dust emission spectrum: $$I_{\nu }={\tau }_{\nu }{\displaystyle \int \text{d}}tp(t)B_{\nu }(T),$$(11)

where τ_ν is the opacity along the line-of-sight. We now make a Taylor expansion of the Planck function around a reference logarithm temperature, t₀ = ln T₀: $$\begin{array}{ll}{I}_{\nu }\hfill & ={\tau }_{\nu }[B_{\nu }(T_0)+{\displaystyle \int \text{d}}tp(t)(t-t_0)\frac{\text{d}{B}_{\nu }}{\text{d}t}(T_0)\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \int \text{d}}tp(t){(t-t_0)}^2\frac{{\text{d}}^2{B}_{\nu }}{\text{d}{t}^2}(T_0)].\hfill \end{array}$$(12)

We distinguish the reference temperature, T₀, from the Gaussian distribution mean temperature as defined by ln T_m ≡< ln T > in order to be as general as possible. The integral over the temperature can readily be separated from the frequency dependence such that we obtain, by defining Δt ≡ t_m − t₀, $$\begin{array}{ll}{I}_{\nu }\hfill & ={\tau }_{\nu }[B_{\nu }(T_0)+\Delta t\frac{\text{d}{B}_{\nu }}{\text{d}t}(T_0)\hfill \\ \hfill & +\frac{1}{2}\left(\Delta t^2+{\sigma }^2\right)\frac{{\text{d}}^2{B}_{\nu }}{\text{d}{t}^2}(T_0)],\hfill \end{array}$$(13)

which, if T₀ = T_m, simplifies into $$I_{\nu }={\tau }_{\nu }\left[B_{\nu }(T_0)+\frac{{\sigma }^2}{2}\frac{{\text{d}}^2{B}_{\nu }}{\text{d}{t}^2}(T_0)\right].$$(14)

This last expression shows how the dispersion in logarithmic temperature translates into the addition of a new spectrum. This corresponds to the second derivative of the Planck spectrum (Eq. (8)) times half the square of the dispersion (and times the emissivity).

The first-order derivative function disappears in the case of any symmetrical temperature distributions like the Gaussian one. The three spectral components of Eq. (13) are displayed in the example of Fig. 1. While the component of index 1 can be of both signs, the component of index 2 is always positive for all frequencies and mostly influences the Wien part of the spectrum. The three components have the same Rayleigh-Jeans spectral behaviour in ∝ ν² T₀τ_ν. The analogy to the Sunyaev-Zel’dovich (SZ) effect, a CMB spectral distortion in clusters of galaxies, is given in Appendix B. The complete formalism of the interaction of a temperature distribution with the Planck function was developed by Stebbins (2007). He shows that logarithm temperature moments are the key quantities to know, irrespective of the exact distribution, here p(t), as can be seen in our Eq. (12). Pitrou & Stebbins (2014) have argued for t ≡ ln T as the proper variable to use in the context of the CMB distortions (including relativity Lorentz frame invariance). We think that there is an additional reason for that use that is linked to the way photometric data are spread in frequency (see below).

2.3 Finding the underlying distribution of dust temperature from photometry

We now see the interest in the decomposition of the dust emission spectrum into the Planck function and its two derivatives (Eq. (13)). With only a few parameters, we can describe quite different situations. Photometricmeasurements are used to deduce dust properties, abundance, and luminosity. They are obtained from broadband detectors: photoconductors, bolometers, and kinetic inductance detectors (KIDs). The width of the band n, of a given instrument and telescope, described by the transmission curve H_n (ν), is so large that one has to assume a reference spectrum, J_ν, to quote a brightness, $I_{\nu }^n$, at the band centralfrequency, ν_n, and for a given underlying spectrum, I_ν: $$I_n=\frac{J_{{\nu }_n}}{{\displaystyle \int \text{d}}\nu H_n(\nu )J_{\nu }}{\displaystyle \int \text{d}}\nu H_n(\nu )I_{\nu }.$$(15)

To simplify further, we define the set of normalised transmission functions for a given instrument as $$G_{\nu }^n=\frac{J_{{\nu }_n}}{{\displaystyle \int \text{d}}\nu H_n(\nu )J_{\nu }}{H}_n(\nu ),$$

and we can rewrite the bandpass Eq. (15) as: $$I_n=\langle G_{\nu }^n.I_{\nu }\rangle ,$$(16)

where the bracket indicates a simple integration over the frequencies.

The reference spectrum, J_ν, is usually taken as the first derivative of the CMB Planck function for millimetre wavelengths (including the 353 GHz Planck mission band), $$J_{\nu }=\frac{P{\nu }^3}{T_{\text{CMB}}}g(\frac{h\nu }{k{T}_{\text{CMB}}}),$$

and the J_ν ∝ ν⁻¹ Infrared Astronomical Satellite (IRAS) convention at submillimetre wavelengths. The ratio I_ν (ν_n)∕I_n is close to one and is usually called the colour correction.

The bandpass H_n has a typical width of Δν_n and is usually characterised by a resolution of $\frac{{\nu }_{\text{n}}}{\Delta {\nu }_{\text{n}}}~3$ irrespective of ν_n. So the effective spectral resolution is coarse but constant in logarithm terms. It is thus appropriate to match this logarithmic scale with the temperature grid that we seek. For that purpose, it is worth remembering that dust cannot be colder than the CMB. Therefore, we propose discretising the temperature as follows: $$t_l=\log T_l=\log T_{\text{CMB}}+l.s,$$(17)

or T_l = T_CMB.e^l.s, where l is an integer starting at 0 for the CMB temperature. The logarithmic temperature step, s, should match the photometric gridding, so values of $\frac{1}{3}$ or $\frac{1}{2}$ seem appropriate. The grid upper limit depends on the chosen experiments. Defining $x_l\equiv \frac{h\nu }{k{T}_l}$, and using Eq. (13), we can now rewrite the dust brightness as a linear decomposition of known functions that are fixed by the grid of Eq. (17): $$\begin{array}{ll}{I}_{\nu }\hfill & =P{\nu }^3{\sum }_{l=0}{\tau }_{\nu }^l[f(x_l)+\Delta t_{l}g(x_l)\hfill \\ \hfill & +\frac{1}{2}\left(\Delta t_l^2+{\sigma }_l^2\right)g(x_l)(h(x_l)-1)],\hfill \end{array}$$(18)

where the l component is described byits opacity law, ${\tau }_{\nu }^l$, and its temperature distribution, which is solely characterised by two parameters: the logarithm mean temperature shift, $\Delta t_l=t_l^m-t_l$, and the dispersion, σ_l. In order to compare the model with observations, we must take into account the passage of light through an instrument, by combining Eqs. (18) and (16). For that purpose, we normalise the opacity ${\tau }_{\nu }^l$ to a reference opacity, τ_l (say, at a given frequency). We can then define three series of coefficients, which describe the coupling between the dust temperature grid and the set of photometric bands: $$\begin{array}{ll}{K^{\prime }}_{0,n,l}\hfill & \equiv \langle G_{\nu }^n.\left[P{\nu }^3\frac{{\tau }_{\nu }^l}{{\tau }_l}f(x_l)\right]\rangle ,\hfill \end{array}$$ $$\begin{array}{ll}{K}_{1,n,l}\hfill & \equiv \langle G_{\nu }^n.\left[P{\nu }^3\frac{{\tau }_{\nu }^l}{{\tau }_l}g(x_l)\right]\rangle ,\hfill \end{array}$$ $$\begin{array}{ll}{K}_{2,n,l}\hfill & \equiv \langle G_{\nu }^n.\left[P{\nu }^3\frac{{\tau }_{\nu }^l}{{\tau }_l}g(x_l)(h(x_l)-1)\right]\rangle ,\hfill \end{array}$$(19)

and we obtain the photometric expectation in band n: $$\begin{array}{ll}{I}_n\hfill & ={\sum }_{l=0}{\tau }_l[{K^{\prime }}_{0,n,l}+\Delta t_l{K}_{1,n,l}+\frac{1}{2}\left(\Delta t_l^2+{\sigma }_l^2\right)K_{2,n,l}].\hfill \end{array}$$(20)

We now need a digression with respect to the first temperature bin at T_CMB, which is going to modify K′.

Fig. 1 Three elementary spectra of dust (see Eq. (13)) in a log-linear scale. For illustration purposes, we take the case of 20 K dust with an emissivity of ${10}^{-4}{(\frac{\nu }{353\text{GHz}})}^{1.5}$. The component of index 0 is the thermal brightness of dust. The component of index 1 comes from the first derivative of the Planck function (Eq. (6)) with $\frac{\Delta T}{T}=\Delta t=0.1$, and the component 2 is obtained from the second derivative (Eq. (8)) using a logarithmic temperature dispersion ofσ = 0.15.

2.4 The radiative transfer of the CMB through dust

We now want to take into account the dominant background, which goes through interstellar dust before reaching the instruments¹. It happens that dust, although a very thin layer with opacities much smaller than 1, has a non-negligible impact on the CMB and vice versa. The radiative transfer equation is solved for a single-temperature dust to get $$I_{\nu }=B_{\nu }(T_{\text{CMB}})+(1-e^{-{\tau }_{\nu }})\left[B_{\nu }(T_{\text{d}})-B_{\nu }(T_{\text{CMB}})\right],$$

which can be rewritten, for the optically thin case, as the real modified blackbody (RMBB) law: $$I_{\nu d}={\tau }_{\nu }\left[B_{\nu }(T_{\text{d}})-B_{\nu }(T_{\text{CMB}})\right],$$(21)

where we have suppressed the CMB monopole term, which is cancelled out in the measurement by the current differential experiments, except the Far Infrared Absolute Spectrophotometer (FIRAS) on the Cosmic Background Explorer (COBE). The RMBB law is important in the millimetre domain. For example, at 100 GHz the correction amounts to a 7% downward effect for 18 K dust, irrespective of the emissivity law. The RMBB can be compared with the occultation of the CMB by planets, which is taken into account in the calibration of cold planets. Equation (21) is equivalent to saying that in Eq. (20) the first term, ${K^{\prime }}_{0,n,l}$, is special and has to be modified into $K_{0,n,l}\equiv {K^{\prime }}_{0,n,l}-{K^{\prime }}_{0,n,0}$ so that very cold dust (at the CMB temperature) becomes invisible (a consequence of Kirchhoff’s law). The influence of the CMB on the submillimetre spectrum of high-redshift galaxies has also been noted by da Cunha et al. (2013). Here we show that, independent of redshift, the CMB has to be explicitly included in the dust SED.

We can now rewrite Eq. (20) as $$\begin{array}{ll}{I}_n\hfill & ={\sum }_{l=0}{\tau }_l[K_{0,n,l}+\Delta t_l{K}_{1,n,l}+\frac{1}{2}\left(\Delta t_l^2+{\sigma }_l^2\right)K_{2,n,l}].\hfill \end{array}$$(22)

Although K_0,n,0 = 0, very cold dust can still manifest itself with the other two CMB terms: K_1,n,0 and K_2,n,0. We note that it is more likely coming from the CIB than from interstellar dust. Elfgren & Désert (2004) have made a case for part of the CIB being in the spectral form encapsulated by the coefficients K_1,n,0.

3 Discussion

We think that the result described in Eq. (22) is general enough to encompass the description of many dust environments. It includes the modified blackbody (MBB) case if there is only one l component, without dispersion (σ_l = 0). The data at hand in photometric bands $I_n^{\text{meas}}$ (with the associated error bar ${\sigma }_n^{\text{meas}}$) can be easily compared with this model. If we assume a dust emissivity law (identical or not for different l components), the coefficients K can be pre-computed once and for all, thus making the fitting process very fast for many lines of sight: $${\chi }^2={\sum }_n{\left(I_n^{\text{meas}}-{\sum }_{l=0}{\sum }_{i=0}^2{c}_{i,l}{K}_{i,n,l}\right)}^2/{({\sigma }_n^{\text{meas}})}^2.$$(23)

If we keep only the first two coefficients, K_0,n,l and K_1,n,l, the fit is strictly linear, with the condition that c_0,l must be positive. The number of components, as described by the number of useful l components, must be tuned in the process in order to keep $\Delta t_l\equiv t_l^m-t_l=\frac{c_{1,l}}{c_{0,l}}$ small. In that sense, the process is parsimonious: finding one l component that is best at minimising the χ² and testing to see if a second one is really needed, then iterating. If we add the temperature distribution with K_2,n,l and σ_l, and compare Eqs. (22) with (23), we can recover a linear fit if we include Lagrange coefficients to satisfy the condition $c_{2,l}\ge \frac{c_{1,l}^2}{2c_{0,l}}$.

The true non-linear part of the fit is in the emissivity law, which is hidden in the K coefficients. What is also hidden is the absolute calibration uncertainty of the instrument. This second set of values is degenerate with K. Another complication arises with the emissivity law, which is, in principle, different for different l. One can assume the same emissivity for a starter and refine that hypothesis if it happens to be really necessary. From Fig. 1 we see that a temperature spread always produces a positive distortion at all frequencies. The effect is that an MBB fit will tend to overestimate the dust temperature and to undervalue the emissivity index. So not only do we have statistically anti-correlated error bars between these two quantities (Shetty et al. 2009), but the (genuine) temperature spread along the line-of-sight can cause an apparent T-β anti-correlation too.

One could find that this model has ‘too many’ parameters. But we would argue that for the case $s=\frac{1}{3}$, this corresponds to only five discrete temperatures between 8 and 30 K, namely 9.1, 12.2, 16.5, 22.3, and 30 K. Secondly, the fitting procedure can cull the components as the parameters in Eq. (23) have to conform to positivity constraints set by Eq. (22). With the case shown in Fig. 1, a single-temperature logarithmic change by $\frac{s}{2}$ yields a maximum residual of 0.6% relative to the peak emission for three components fitted with Eq. (23) (3.5% if only two components are used).

Another constraint will show up when one finds that the temperature spread is larger than the bin size, s. In that case, the parameter σ should saturate to $\frac{s}{2\sqrt{3}}$ (the equivalent of uniform noise in a digital binning process), and hence ${\sigma }_{\text{sat}}^2~0.009{\left(\frac{s}{0.33}\right)}^2$.

The bolometric luminosity of dust can be directly computed from a linear combination of integrals pre-computed from Eq. (18). Polarisation photometry follows the same rules as intensity. We can expect the degree of dust polarisation to be strong when the temperature distribution is not too spread out (i.e. for small σ).

The CIB fluctuations (Planck Collaboration XXX 2014) may be revisited in light of this model because of the linearisation procedure, as it can lend itself to the powerful statistical tools used for CMB studies, including covariances and model testing with likelihood functions, via the statistics of the c_i,l coefficients.

This work is complementary to the approach by Chluba et al. (2017). Both studies put the emphasis on Planck derivatives and temperature moments. Chluba et al. (2017) argues for $\frac{1}{T}$ as the proper expansion variable. Weprefer ln T because it will fare better, in terms of the number of required steps in the temperature grid, for broad temperaturedistributions (including the CIB and the example in Sect. 1), where more than one l component is necessary. In that case, high-frequency convergence problems, signalled by Chluba et al. (2017), are also alleviated. Moreover, for broadband experiments, it is impossible to go beyond the second Planck derivative as the fitting system becomes too under-determined for limited signal-to-noise ratio measurements.

4 Conclusions

We have shown an analytical development around the Planck function. We have emphasised the role of its second derivative with respect to the logarithm of the dust temperature. This is important in order to include the effect of a very likely distribution of temperatures along the line-of-sight around the main dust temperature. We have also devised an economical way of accounting for a wider temperature range by discretising the temperature ladder in constant logarithm steps, which naturally starts at the CMB temperature. By explicitly including the instrumental measurement process, we have devised a model that can be fit in an almost linear way, from photometry to dust temperature distributions, thus accelerating the computation of the inverse problem. Further work will have to be done to implement these findings for the analysis of dust in the Planck and Herschel missions, in particular the optimal gridding of dust temperatures versus the experimental photometric frequency sampling. This work is complementary to investigations of the possible variations in the emissivity index in some galactic regions (Mangilli et al. 2021; Rigby et al. 2018; Bracco et al. 2017; Tang et al. 2021; Planck Collaboration Int. XIV 2014). Applications to the CIB statistics and dust as the major foreground for CMB studies should be sought too.

Acknowledgements

The author thanks discussions with François-Xavier Hamel, Guilaine Lagache and Nicolas Ponthieu.

Appendix A Planck derivative functions

The third Planck derivative is useful for assessing the errors made by using only two derivatives in the Taylor expansion. For completeness, it is given here: $$\frac{d^3{B}_{\nu }}{d{t}^3}=P{\nu }^{3}g(x)\left[{(h(x)-1)}^2-h(x)+2xg(x)\right].$$(A.1)

It is worth showing the graph of Planck derivative functions. As they all diverge at low frequencies as $\frac{1}{x}$, we can show instead x times the functions. In Fig. A.1 we show the Planck function and its three derivatives (with the index going from 0 to 3).

Fig. A.1 Four Planck functions (times x) as a functionof the dimensionless frequency parameter x. The Planck function has the 0 index, and the three indices label the three consecutive Planck function derivatives. Hence, we have plotted x f(x), x g(x), x g(x) (h(x) − 1), and $xg(x)\left[{(h(x)-1)}^2-h(x)+2xg(x)\right]$, where the functions f, g, and h are defined in Eqs. 5, 7, and 10.

Appendix B The Sunyaev-Zel’dovich effect

The second derivative of the Planck function is linked to the SZ effect, a spectral distortion of the CMB through clusters of galaxies. The SZ effect is quantified by the Compton parameter, y, which is proportionalto the integrated thermal electron pressure along the line-of-sight. Indeed, if we specifically perform the substitution Δt = −3y and σ² = 2y in Eq. 13, where T₀ = T_CMB is the present CMB temperature (Stebbins 2007), we recover the SZ non-relativistic thermal distortion: $$I_{\nu }^{tSZ}=yP{\nu }^{3}g(x)(h(x)-4).$$(B.1)

This gives a direct interpretation of the SZ distortion: the Compton effect (while maintaining the number of photons constant) shifts the CMB average temperature along the cluster line-of-sight by a relative factor of − 3y, and the Maxwellian electron velocity distribution interacts with the CMB photons as if they produced a Gaussian temperature distribution with a dispersion σ. This explains why the quadratic equivalent temperature dispersion in logarithm is proportional to the quadratic Maxwellian velocity distribution (i.e. the electron temperature, T_e): ${\sigma }^2=2\frac{k{T}_e}{m_e{c}^2}\tau$, where m_e is the electron mass and τ is the line-of-sight opacity to the Compton effect.

One of the differences between the SZ spectral distortion (in h − 4) and the Planck second derivative (in h − 1) is that the SZ effect has a negative relative distortion in the Rayleigh-Jeans part of the spectrum, the so-called SZ decrement: $\frac{\Delta I_{\nu }}{I_{\nu }}\simeq -2y$.

References

All Figures

Fig. 1 Three elementary spectra of dust (see Eq. (13)) in a log-linear scale. For illustration purposes, we take the case of 20 K dust with an emissivity of ${10}^{-4}{(\frac{\nu }{353\text{GHz}})}^{1.5}$. The component of index 0 is the thermal brightness of dust. The component of index 1 comes from the first derivative of the Planck function (Eq. (6)) with $\frac{\Delta T}{T}=\Delta t=0.1$, and the component 2 is obtained from the second derivative (Eq. (8)) using a logarithmic temperature dispersion ofσ = 0.15.

In the text

	Fig. A.1 Four Planck functions (times x) as a functionof the dimensionless frequency parameter x. The Planck function has the 0 index, and the three indices label the three consecutive Planck function derivatives. Hence, we have plotted x f(x), x g(x), x g(x) (h(x) − 1), and $xg(x)\left[{(h(x)-1)}^2-h(x)+2xg(x)\right]$, where the functions f, g, and h are defined in Eqs. 5, 7, and 10.
In the text