text:
# Two-photon absorption in potassium niobate

## I Introduction

For cw laser second harmonic generation in a build-up cavity, thermal issues in the nonlinear crystal inevitably come into play. This is particularly true at high laser powers. In typical applications, a Fabry-Perot build-up cavity increases the fundamental laser intensity by one or two orders of magnitude compared to the single-pass intensity. The nonlinear crystal absorbs a fraction of the fundamental power, heats up, and changes the cavity properties. Thermal lensing changes the cavity’s spatial mode. Nonlinear absorption changes the cavity’s finesse. Thermal expansion and changes in the crystal’s index of refraction change the cavity’s optical path length.

These effects are well known, and have been observed by several groups working in this field (for example, Arie92 ; Polzik91 ). Thermally-induced changes in the build-up cavity sometimes limit the second-harmonic generation efficiency. In many cases, they are nuisances that can be worked around. However, these changes can be used as a sensitive diagnostic tool for accurately characterizing nonlinear crystal properties.

For gases and atomic vapors, measuring optical properties of materials inside Fabry-Perot cavities has a long history Siegman86 , particularly in relation to optical bistability, saturable absorption, and index of refraction measurements. Thermal and collisional properties of gases have also been measured in this way Ye96 , as well as optical properties of cavity mirrors An97 . However, measurements of optical properties of crystals in Fabry-Perot cavities is relatively new Douillet99 .

In some ways, Fabry-Perot measurements are similar to laser calorimetry Busse94 . In those calorimetry measurements, a laser passes through a thermally isolated absorbing crystal and heats it slightly. From the temperature rise, laser power, thermal mass, and thermal conductivity it is possible to determine an optical absorption coefficient for the crystal. Comparing calorimetry to our Fabry-Perot measurements, both are photothermal methods. Both can have high sensitivities to small absorptions. In contrast to calorimetry, our “thermometer” is the power-induced shift in the cavity resonance frequency. Fabry-Perot measurements can typically reach much higher intensities, making nonlinear effects more apparent. Furthermore, frequency-based techniques, such as the one presented in this paper, usually offer a higher level of precision.

This paper presents a detailed study of the thermal response of KNbO3 inside a Fabry-Perot build-up cavity. We use this response to determine a new value for the optical absorption coefficients in KNbO3 at 846 nm and 423 nm, and to measure the effects of blue-light-induced infrared absorption (BLIIRA). We also measure the temperature-dependence of the index of refraction for light polarized along the $`b`$ and $`c`$ crystal axes. These measurements are compared to values from the literature.

## II Experimental Setup

Our experimental setup is typical for cw second-harmonic generation in an external build-up cavity (see Figure 1). This laser system is designed for a calcium laser cooling experiment at 423 nm. It produces typically 75 mW at 423 nm, measured outside the cavity.

The laser system consists of a single frequency extended cavity diode laser at 846 nm amplified in a single-pass through a tapered laser diode NewFocus . The extended cavity diode laser is frequency-stabilized relative to a separate optical cavity using the Pound-Drever method Pound83  to provide both short and long-term stability. After optical isolation, polarization correction, and mode matching, we inject 200 mW into the Gaussian mode of a four-mirror folded (bowtie) Fabry-Perot frequency doubling cavity. Higher-order modes in the cavity are less than 5% of the Gaussian mode. The input coupler reflectivity is 96.6%, and the other three mirror reflectivities are $`>99.5`$% in the infrared. The radius of curvature for the two curved mirrors is 100 mm, and they are separated by a distance of 117 mm. The angle of incidence on the curved mirrors is $`6.5^\text{o}`$, minimizing the ellipticity of the beam waist inside the crystal (41.1 and 40.7 $`\mu `$m in the tangential and saggital planes, as determined from the cavity geometry). The round-trip cavity length is 618 mm. The KNbO3 crystal VLOC  is $`a`$-cut and antireflection-coated ($`R<`$0.25%) at 846 nm. The non-critical phase-matching temperature is $`11.5^\text{o}`$C, for light polarized along the $`b`$-axis. The crystal is mounted on a Peltier-cooled copper block. An AD590 temperature sensor and a thermister embedded in the copper mount monitor the crystal temperature. The entire bowtie cavity is enclosed in a sealed aluminum box, filled with dry oxygen at atmospheric pressure. For low laser powers, and when the crystal temperature is far from the phase matching temperature, the cavity finesse with the crystal in place is 124. Without the crystal, the cavity finesse is 140.

## III Transmission Lineshape Without Blue Light

When laser light incident on the bowtie cavity is polarized along the $`c`$-axis, it is impossible to meet the phase-matching condition for non-critical second harmonic generation at $`\lambda =846`$ nm, and no blue light is generated. Previous measurements in KNbO3 have shown that even small amounts of blue light dramatically change the infrared absorption properties Goldberg93 ; Mabuchi94 ; Shiv95 . We address this blue-light-induced-infrared-absorption (BLIIRA) later in the paper. However, in this section we discuss infrared absorption in the absence of blue light, i.e.: when the incident laser beam is polarized along the $`c`$-axis.

One of the flat mirrors in the bowtie cavity is mounted on a piezoelectric crystal, and we use this crystal to change the optical length of the cavity. A small amount of infrared light is transmitted through the cavity mirrors, and we measure this light while sweeping the cavity length. When the piezo-mounted mirror passes through the position for cavity resonance (i.e.: when the roundtrip pathlength in the cavity is an integer multiple of the laser wavelength), the power circulating in the cavity, and transmitted through the cavity end mirror, increases. The transmission lineshape is Lorenzian in the absence of disturbing effects:

$$P(\nu )=\frac{P_0}{1+F\left(\nu \nu _0\right)^2},$$
(1)
where $`\nu _0`$ is the resonant frequency of the cavity, $`2/\sqrt{F}`$ is the full width at half-maximum (FWHM), and $`P_0`$ is the infra-red power at the peak of the Lorenzian.

The KNbO3 crystal inside the bowtie cavity absorbs a small amount of the laser power and heats up. This changes the crystal’s length and index of refraction, which in turn changes the optical path length of the bowtie cavity. We model the overall effect by replacing $`\nu _0`$ with $`\nu _0+\eta P/P_0.`$ Accordingly, Equation 1 becomes

$$y=\frac{1}{1+F\left(\nu \eta y\right)^2}.$$
(2)
where $`y=P/P_0`$ is the normalized transmission, and the frequency $`\nu `$ is measured relative to the unshifted frequency $`\nu _0`$. Transmission lineshapes of this kind are typical of optically bistable devices. A similar treatment of this effect is given in the references Douillet99 .

The peak in the transmission curve occurs when $`y=1`$ and $`\nu =\eta `$. Figure 2 shows the transmission lineshape for both low and high laser powers. On the same figure we plot a fit of the data to Equation 2. The asymmetry in the transmission lineshape at high powers is due to thermal changes in the optical pathlength. When sweeping the cavity from low to high frequencies as in Figure 2, the power-dependent cavity resonance is effectively “dragged along” from $`\nu =\nu _0`$ to $`\nu =nu_0+\eta `$ as the cavity length changes. This is known as “thermal self-locking” Ye96 . The value of $`\eta `$ is called the “self-locking range,” and it is related to the “thermo-optic constant” of reference Douillet99 .

For a given power incident on the cavity, we determine the self-locking range $`\eta `$ using a nonlinear least-squares fit of the measured transmission lineshape to Equation 2. For the highest laser intensities, it appears that a better overall fit is obtained by replacing $`\eta y`$ in Equation 2 by a higher order polynomial. However, the fit in these high intensity data is hampered because the fit function is not single-valued, as can be seen in Figure 2. For these measurements, we determine the values of $`\eta `$ visually by finding the frequency at the peak of the transmission curves.

Our values for $`\eta `$ as a function of infrared power in the bowtie cavity are shown in Figure 3a. The error bars indicate the 1 $`\sigma `$ statistical spread in 5 measurements at each power setting. For linear absorption, the cavity shift, $`\eta `$, should be a linear function of the circulating power. However our measurements indicate a quadratic dependence of the cavity shift on the circulating power, suggesting that the dominant absorption mechanism is two-photon absorption in the infrared.

We convert $`\eta `$ into a temperature rise inside the crystal by measuring the temperature change necessary to shift the cavity resonance frequency by one free spectral range. In these calibration measurements, we attenuate the incident laser power to 10 $`\mu `$W to avoid laser-heating in the crystal. The cavity mirror positions remain fixed, and the laser frequency is constant. Using the Peltier device, we change the temperature of the crystal by several degrees while monitoring the cavity transmission. At $`T=11.5^\text{o}`$C, the resonant frequency of the bowtie cavity shifts at a rate of 231 MHz $`^\text{o}`$C-1 for light polarized along the $`b`$-axis and 331 MHz $`^\text{o}`$C-1 for light polarized along the $`c`$-axis.

We calibrate our photodiode signal using a power meter, and measure the efficiency of the optical path to convert this IR signal to a power circulating inside the cavity. The conversion from cavity shift to temperature rise can also be derived from measurements of the change in the index of refraction with temperature. From the published values of $`dn/dT`$ Ghosh94  and the geometry of our bowtie cavity, we calculated what the temperature-dependent shift ought to be. The calculation agrees with our data for the $`b`$-axis. However, we find a small difference between the calculation and our data for the $`c`$-axis, as discussed later in this paper.

Using our measurements, we convert $`\eta `$ into $`\mathrm{\Delta }T`$, the difference in temperature between the crystal axis where the laser propagates and the crystal wall which is kept constant by our temperature control circuit. Figure 3b shows a plot of the temperature rise as a function of the infrared power circulating inside the cavity. It is relatively simple to convert this temperature rise into a measure of the optical absorption coefficient of KNbO3 by solving the heat transfer equation, which we do below.

## IV Solving the Heat Transfer Equation

We solve the heat transfer equation in cylindrical coordinates. The laser beam is in the TEM00 Gaussian mode, with a 41 $`\mu `$m waist. The confocal parameter of the beam is 6.2 mm, and at the face of the crystal, the beam waist is only 7% larger than in the center of the crystal. So we approximate the laser beam as axially symmetric and Gaussian with a 41 $`\mu `$m waist along the entire 10 mm length of the crystal. Our crystal has a square cross section, 3 mm on a side. However, because it is so much larger than the laser beam waist, we assume that the crystal is cylindrical also, with a 1.5 mm radius. With these assumptions, the equilibrium temperature inside the crystal is determined by solving the heat equation in the form:

$$\frac{1}{r}\frac{d}{dr}\left(r\frac{du}{dr}\right)=\frac{G}{K},$$
(3)
where $`u(r)`$ is the temperature as a function of radius, $`G`$ is the power density, and $`K`$ is the thermal conductivity. Following Sutherland Sutherland96 , the absorbed power density inside the crystal can be written as

$$G=\alpha I+\beta I^2$$
(4)
where $`I=(2P/\pi a^2)\mathrm{exp}(2r^2/a^2)`$ is the laser intensity, $`P`$ is the laser power incident on the crystal, $`a`$ is the Gaussian beam waist, $`\alpha `$ is the linear optical absorption coefficient, and $`\beta `$ is the two-photon absorption (TPA) coefficient.

We can simplify Equation 3 by letting $`r=ay/\sqrt{2}`$, $`A=(P\alpha /2\pi K)`$, and $`B=(P^2\beta /2\pi ^2a^2K)`$. With these substitutions, Equation 3 integrates once to

$$\frac{du}{dy}=\frac{A}{y}\left[\mathrm{exp}(y^2)1\right]+\frac{B}{y}\left[\mathrm{exp}(2y^2)1\right],$$
(5)
where the constant of integration is chosen to avoid a singularity as $`y0`$. The solution to this differential equation is

$$u(y)=\frac{A}{2}E_1(y^2)\frac{A}{2}\mathrm{ln}(y^2)\frac{B}{2}E_1(2y^2)\frac{B}{2}\mathrm{ln}(y^2)+C$$
(6)
where $`E_n(x)`$ is the exponential integral Abramowitz64 , and $`C`$ is the constant of integration. This constant is determined by requiring the temperature at the wall of the crystal to be $`T_c`$ (i.e.: $`u(r=b)=T_c`$), which we control in the experiment. With this constraint, we can evaluate the temperature rise on axis inside the crystal. The series expansion for the exponential integral has a logarithmic term that exactly cancels the $`\mathrm{ln}(y^2)`$ terms in Equation 6.

$`\mathrm{\Delta }T`$ $`=`$ $`u(0)T_c={\displaystyle \frac{A}{2}}\left[\mathrm{ln}\left({\displaystyle \frac{2b^2}{a^2}}\right)+\gamma \right]+{\displaystyle \frac{B}{2}}\left[\mathrm{ln}\left({\displaystyle \frac{2b^2}{a^2}}\right)+\gamma +\mathrm{ln}(2)\right]`$ (7)
$`=`$ $`{\displaystyle \frac{\alpha }{4\pi K}}\left[\mathrm{ln}\left({\displaystyle \frac{2b^2}{a^2}}\right)+\gamma \right]P+{\displaystyle \frac{\beta }{4\pi ^2a^2K}}\left[\mathrm{ln}\left({\displaystyle \frac{4b^2}{a^2}}\right)+\gamma \right]P^2,`$
where $`\gamma =.57721\mathrm{}`$ is Euler’s constant.

We can analyze the data in Figure 3 to extract the linear and two-photon absorption coefficients. A weighted fit to the function $`T(P)=c_1P+c_2P^2`$, gives $`c_1=(5.7\pm 3.6)\times 10^4^\text{o}`$C W-1 and $`c_2=(1.09\pm 0.12)\times 10^3^\text{o}`$C W-2. The uncertainties in these numbers are the statistical uncertainty from the weighted fit to the data. Taking the thermal conductivity from Busse94  $`K=4.0`$ W m$`{}_{}{}^{1}^\text{o}`$C-1, $`a=41\times 10^6`$ m, $`b=1.5\times 10^3`$ m, and $`\gamma =0.57721`$. We determine the one-photon (linear) absorption coefficient to be $`\alpha =(0.0034\pm 0.0021)`$ m-1 and the two-photon absorption coefficient to be $`\beta =(3.16\pm 0.35)\times 10^{11}`$ m W-1. The uncertainties in these absorption coefficients are only statistical. We have a systematic uncertainty in the power measurements of 10%. Our assumption of cylindrical symmetry in the crystal for solving the heat equation also introduces an error. The error is known, but related to this is the fact that our laser does not propagate exactly down the center of the crystal. These errors show up in the $`\mathrm{ln}(b^2)`$ term in Equation 7, and probably add $`8`$% uncertainty to the measurements. Our best numbers for the absorption coefficients are therefore $`\alpha =(0.0034\pm 0.0022)\text{m}^1`$ and $`\beta =(3.2\pm 0.5)\times 10^{11}`$ m/W. These numbers do not reflect the uncertainty in the thermal conductivity, which we do not know. To our knowledge, this is the first determination of the two-photon absorption coefficient at any wavelength in KNbO3.

Earlier work on the linear absorption coefficient for KNbO3 Busse94  found $`\alpha =0.001\text{cm}^1=0.1\text{m}^1`$ at 860 nm. This value is much higher than this work. It may be that some of the absorption measured in previous work was two-photon absorption. That work used laser calorimetry to measure the absorption coefficient at specific laser lines from 457 nm to 1064 nm. At 860 nm, they focused a 900 mW laser beam into a 5 mm crystal, and it is likely that the intensity was high enough for multiphoton absorption to be important.

Our method for determining the one and two-photon absorption coefficients is quite general. It is independent of measurements of light scattering inside the crystal. It is also independent of reflection measurements at the crystal faces. This method is particularly well suited to measurements of optical absorption where the absorption coefficients are comparable to or smaller than scattering and reflection coefficients.

## V Transmission Lineshape With Blue Light

When infrared light incident on the bowtie cavity is polarized along the $`b`$-axis of the KNbO3 crystal, it is possible to meet the phase-matching conditions required for generating the second harmonic at 423 nm. The presence of blue light in the crystal dramatically alters the transmission lineshape. Not only is the blue light itself absorbed by the crystal, it also significantly increases the absorption of the infrared light. This is called “blue-light-induced-infrared-absorption” (BLIIRA), and it has been studied at length in the literature Busse94 ; Goldberg93 ; Mabuchi94 ; Shiv95 . Previous work has found BLIIRA to be significant at blue light intensities down to $`7\times 10^4`$ W/cm2. BLIIRA is minimized at longer wavelengths and at higher crystal temperatures Shiv95 . However, for our calcium work, we require high powers at 423 nm.

As before, we fix the laser frequency and scan the cavity length through the resonance condition. Both blue light at 423 nm and a small amount of infrared light at 846 nm exits the cavity. We separate these wavelengths using dichroic mirrors, which transmit 95% in the blue and reflect 99.5% in the infrared. To make sure that no blue light reaches the infrared laser beam detector, we use four of these mirrors in series in the infrared beam path after the doubling cavity (see Figure 1).

In these experiments, we do not independently control the infrared and blue beam intensities. Rather, we optimize the crystal temperature to maximize the blue light production for each infrared power setting in the steady state. Losses inside the doubling cavity limit the maximum infrared power circulating inside and therefore the maximum blue light produced.

The transmission lineshape for the cavity when blue light is present is shown in Figure 4a for our highest IR powers. The temperature rise inside the crystal is significant. At maximum power, as the cavity approaches the resonance condition, only a small amount of blue light is initially generated because the crystal temperature is too low to meet the phase matching condition for second harmonic generation exactly (see Figure 4b). Closer to the cavity resonance condition, the circulating power increases, the crystal heats up, and more blue light is generated, which heats the crystal even more. This positive feedback continues, and the cavity demonstrates thermal self-locking for up to 90 MHz, as shown in Figure 4a.

The blue light signal in Figure 4b has a shoulder on the low-frequency side of the maximum. This is a Maker fringe. As stated previously, we initially optimize the crystal temperature for optimum blue light production in the steady state. At the highest intensities, the temperature shift ($`0.4^\text{o}`$C) due to thermal self-locking of the cavity is a few times the temperature phase matching bandwidth.

Because of these complicated thermal conditions, and especially because the blue light production is not constant across the transmission lineshape, we are reluctant to fit the lineshape to any model for the highest laser intensities. Instead, we find the self-locking range graphically. It is the falling step at the far right hand side of the infrared transmission peak shown in Figure 4. This choice is valid as long as there is no significant thermally-induced change in the cavity coupling efficiency. We can monitor these changes by measuring the light reflected from the cavity. For the relatively modest intensities in this study, the minimum reflected light (which corresponds to the maximum light inside the cavity) changes by only a few percent, and only at the highest intensity measurements. The statistical scatter in the data is larger than this, and we neglect this small systematic error.

At lower laser intensities, where BLIIRA and other absorption processes are less severe, we use a least-squares method similar to Section 3 to find the self-locking range $`\eta `$. The self-locking range is plotted as a function of infrared laser power in Figure 5a. Note that except for the two highest power measurements, the frequency doubling process is not saturated (see Figure 4c).

With blue light present in the crystal, the temperature rise inside the crystal has three sources: infrared absorption (one and two photon), blue light absorption, and BLIIRA. We can easily subtract off the infrared contribution to the self-locking range using the data in Figure 3 as a look-up-table. The remaining contribution to the self-locking range is due to temperature rise from blue light absorption and BLIIRA.

Because the blue light intensity and the temperature are not constant along the length of the axis of the crystal, we must be careful about how we define the temperature rise and the blue light intensity. We will define the absorption coefficient at 423 nm, $`\delta `$, by the equation

$$\mathrm{\Delta }P=_VI_b\delta 𝑑V=\delta P_bL$$
(8)
where $`I_b`$ is the blue light intensity inside the crystal, $`P_b`$ is the measured blue power, and $`L`$ is the length of the crystal. In the non-depleted pump plane-wave approximation, the blue light intensity grows quadratically with distance $`z`$ in the crystal. With the definition from Equation 8, the blue light intensity can be written as $`I_b(r,z)=(12P_bz^2/\pi a^2L^2)\mathrm{exp}(4r^2/a^2)`$, where $`a`$ is the Gaussian waist for the infrared beam in the crystal and we explicitly assume that the Gaussian waist for the blue beam is $`a/\sqrt{2}`$.

The heat equation for blue light absorption and BLIIRA is written as

$`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}\left(r{\displaystyle \frac{u}{r}}\right)+{\displaystyle \frac{^2u}{z^2}}`$ (9)
$`=`$ $`z^2\left(Ce^{4r^2/a^2}+De^{6r^2/a^2}\right).`$
The first term on the right hand side of the equation is the blue light absorption term where $`C=12P_b\delta /(\pi a^2L^2K)`$. The second term on the right hand side of the equation is the BLIIRA term, where $`D=24P_bP_{ir}\xi /(\pi ^2a^4L^2K)`$, $`P_{ir}`$ is the measured infrared power, assumed to be constant along the crystal length, and $`\xi `$ is the BLIIRA absorption coefficient. As before, $`u(r,z)`$ is the temperature inside the crystal and $`K`$ is the thermal conductivity.

The solution of this equation for $`r=0`$ is

$`\mathrm{\Delta }T(z)`$ $`=`$ $`{\displaystyle \frac{3\delta }{4\pi K}}\left[\mathrm{ln}\left({\displaystyle \frac{4b^2}{a^2}}\right)+\gamma \right]{\displaystyle \frac{z^2}{L^2}}P_b`$ (10)
$`+`$ $`{\displaystyle \frac{\xi }{a^2\pi ^2K}}\left[\mathrm{ln}\left({\displaystyle \frac{6b^2}{a^2}}\right)+\gamma \right]{\displaystyle \frac{z^2}{L^2}}P_bP_{ir},`$
but we have to be careful about how we define the temperature rise. In Equation 9, the radial derivatives are much larger than the longitudinal derivatives, and the $`^2u/z^2`$ term can be treated perturbatively. This approximation produces a temperature function that increases quadratically with $`z`$ in the crystal.

We can write the temperature rise on axis as $`\mathrm{\Delta }T(z)=T^mz^2/L^2`$, where $`T^m`$ is the maximum temperature rise in the crystal. We can calculate the change in the optical pathlength inside the crystal with this temperature profile. It is equal to the change in pathlength due to an equivalent temperature $`\mathrm{\Delta }T^{}=T^m/3`$. More expressly, we can write $`T(z)=T^mz^2/L^2=3T^{}z^2/L^2`$, and rewrite Equation 10 as

$`T^{}`$ $`=`$ $`{\displaystyle \frac{\delta }{4\pi K}}\left[\mathrm{ln}\left({\displaystyle \frac{4b^2}{a^2}}\right)+\gamma \right]P_b`$ (11)
$`+`$ $`{\displaystyle \frac{\xi }{3a^2\pi ^2K}}\left[\mathrm{ln}\left({\displaystyle \frac{6b^2}{a^2}}\right)+\gamma \right]P_bP_{ir},`$
Knowing this, we can convert the measured self-locking range of the cavity in the presence of blue light into a temperature rise. The data is plotted in Figure 5. For this data, the “temperature rise” is the average temperature rise in the crystal, i.e.: $`T^{}`$. The “blue power” is the measured blue power, i.e.: the power at the end of the crystal, $`P_b`$.

Again, we fit a function of the form

$$\mathrm{\Delta }T=a_1P_b+a_2P_bP_{ir}$$
(12)
to the data in Figure 5b using a multivariable least-squares method and find $`a_1=2.41^\text{o}`$C W-1 and $`a_2=1.04^\text{o}`$C W-2. Using equation 11, we determine the absorption coefficients $`\delta =13.3\text{m}^1`$ and $`\xi =2.2\times 10^8\text{m/W}`$. Apparently, there are no measurements in the literature of the KNbO3 absorption coefficient at 423 nm. However, there are measurements close to this wavelength Busse94 , and a value of 13.3% per cm is reasonable. For this determination of the absorption coefficient, we have the same kinds of uncertainties as before, namely in the beam position in the crystal and in the power measurement. We estimate the error in this measurement to be 15 %. Our final number for the absorption coefficient at 423 nm is $`\delta =(13\pm 2)`$ m-1.

This uncertainty estimate is perhaps conservative. We have a second KNbO3 crystal. In a simple transmission measurement at 423 nm, using the blue light after the frequency doubling cavity, we find the linear absorption coefficient to be 13% per cm, in agreement with the work above.

We verified the validity of our perturbative approximation by solving Equation 9 numerically. In the numerical solution, the temperature profile is quadratic in $`z`$ except at the very end of the crystal, where is flattens out slightly to meet the boundary condition of no heat flowing out the crystal face. With our values of $`\delta `$, $`\xi `$, $`a`$, $`P_b`$, and $`P_{ir}`$, the calculated temperature $`T^{}`$ agrees with our data to about one percent.

It is interesting that the above treatment of BLIIRA does not distinguish between blue-light induced infrared absorption and infrared light induced blue-light absorption. This model of BLIIRA only postulates that there is a two-color two-photon absorption cross-section. This treatment of two-color two-photon absorption in KNbO3 is perhaps simplistic. However, it serves the purpose in this paper of demonstrating that thermal self-locking of a Fabry-Perot cavity can be used to sensitively characterize the linear and nonlinear optical properties of an absorbing crystal.

## VI New Measurements of $`dn/dT`$

From the calibration measurements that convert $`\eta `$ in Figures 3 and 5 to a temperature rise, we can extract the temperature-dependent change in the index of refraction for KNbO3. When the crystal temperature changes, the optical pathlength of the bowtie cavity also changes because both the crystal length and the index of refraction depend on temperature:

$$\frac{d}{dT}=L_c\left[\frac{dn_x}{dT}+\alpha _T\left(n_xn_0\right)\right],$$
(13)
where $``$ is the cavity optical pathlength, $`L_c`$ is the crystal length (10 mm), $`\alpha _T`$ is the coefficient of linear expansion for KNbO3, and $`n_x`$ and $`n_0`$ are the indices of refraction for the crystal and air, respectively.

In a relatively simple measurement, we can determine $`dn/dT`$ for the crystal across a wide range of temperatures. We fix the laser frequency and the cavity mirror positions. We attenuate the laser intensity to 10$`\mu `$W to avoid thermal self-locking. The crystal temperature can be independently controlled because it is mounted on a Peltier device. We measure the cavity transmission while slowly sweeping the crystal temperature over several degrees. At certain temperatures, the cavity optical pathlength meets the resonance condition, and we measure a peak in the cavity transmission. The temperature between transmission peaks in these measurements corresponds to changing the optical pathlength by one optical wavelength, 846 nm. We repeat these measurements over a wide temperature range to determine $`d/dT`$ from -12 $`^\text{o}`$C to about 50 $`^\text{o}`$C.

Using equation 13, we convert $`d/dT`$ to $`dn/dT`$. With $`n_c=2.13`$, $`n_b=2.28`$, $`\alpha _T=5\times 10^6{}_{}{}^{\text{o}}\text{C}_{}^{1}`$, and $`d/dT=\lambda /\mathrm{\Delta }T`$, we determine $`dn/dT`$ for several temperatures. Our values are shown in Figure 6 for light polarized along the $`c`$-axis and $`b`$-axis, as labeled. Also shown in the figure is the data from Ghosh Ghosh94 . The line from Ghosh is a fit to his measurements of the index of refraction between 0 and 140 $`^\text{o}`$C. We extrapolate that fit to lower temperatures to compare with our data.

For light polarized along the $`c`$-axis, our measurements agree with those of Ghosh over our entire temperature range. For light polarized along the $`b`$-axis, in the range 5 to 50 $`^\text{o}`$C, which is roughly the range for which our data overlaps Ghosh’s data, our measurements also agree. However, below $`5^\text{o}`$C, our data departs slightly from Ghosh’s.

To check our data, we made an additional measurement. When the polarization of the light incident on the cavity is rotated slightly relative to the $`c`$-axis, the transmission lineshape has two features, one for each polarization. Because the index of refraction is different for the two polarizations ($`b`$ and $`c`$), the optical pathlength is also different, and each polarization is transmitted at a different cavity length. When we change the temperature of the crystal, the two peaks in the transmission lineshape shift by different amounts in different directions, as shown in Figure 7. We measured these shifts for a set of small temperature changes near -11.5 $`^\text{o}`$C and determined the relative changes in the indices of refraction, $`\left(dn_c/dT\right)/\left(dn_b/dT\right)`$. These measurements confirm our low temperature data in Figure 6. Following the treatment of Ghosh, we fit our $`dn/dT`$ data for the $`b`$-axis to a second-order polynomial, and find

$`{\displaystyle \frac{dn_b}{dT}}`$ $`=`$ $`1.55809\times 10^9T^2`$ (14)
$`+`$ $`4.06912\times 10^8T33.76305\times 10^6`$
where the units of $`dn/dT`$ are $`^\text{o}`$C-1, and $`T`$ is measured in Celsius.

## VII Conclusion

We have demonstrated how thermal self-locking of a Fabry-Perot cavity can be used to determine optical properties of nonlinear crystals. We have determined new values for the one- and two-photon absorption coefficients in KNbO3 at 846 nm. We have also determined new values for the one photon (linear) absorption coefficient at 423 nm, and a BLIIRA absorption coefficient under specific circumstances. Finally, we present new measurements of the temperature-dependent change in the indices of refraction for the $`b`$ and $`c`$ axes.

These measurements all derive from the thermal response of an absorbing Fabry-Perot cavity to light circulating in the cavity. The methods are general and can be widely applied to other systems, particularly those in which the absorption coefficients are small compared to scattering and reflection coefficients.

## VIII Acknowledgements

We express our appreciation to Ross Spencer for his computational assistance. This work is supported in part by a grant from the Research Corporation and from the National Science Foundation under Grant No. PHY-9985027.


format: text

text:
# Radiative Bulk Viscosity

## 1 Introduction

Both the thermal and mechanical effects of radiation on matter play a significant role in the fluid dynamics of hot objects, from the early universe to early stars. The thermal effects may be reasonably well accounted for by the use of the Eddington approximation \[Unno & Spiegel 1966\], though special circumstances may call for improved descriptions. In particular, when travel times of photons are not negligible, retardation effects may be significant \[Delache & Froeschlé 1972\], as indeed they are when we go to treat the dynamical effects of the radiation field.

An early description of the dynamical effects of radiation on matter was provided by L. H. Thomas \[Thomas 1930\] who expanded the specific intensity in terms of the mean free path of photons, as Hilbert had previously done for the one-particle distribution in kinetic theory \[Kogan 1969\]. Each of them assumed that the relevant mean free path is much smaller than the characteristic macroscopic length scales of the material medium. Though Thomas’s treatment included only emission and absorption, Masaki later extended it to include the effects of Thomson scattering \[Masaki 1971, Hsieh & Spiegel 1976\] and Compton scattering \[Masaki 1981\].

The effectiveness of the dissipation in those studies was characterized by the radiative transport coefficients — shear viscosity and conductivity — about whose evaluation there has been general agreement since Thomas’s derivation of the stress tensor for a matter-radiation mixture with short photon mean free paths \[Weinberg 1972, Mihalas & Mihalas 1984\]. But in recent years, particular attention has been given to the bulk viscous effects, which are of interest for applications to cosmology and, for that case, different expressions have been proposed.

A simple relativistic gas has bulk viscosity \[Stewart 1972\] because the ultrarelativistic particles and particles with small energies respond to volume changes at different rates and so are driven out of equilibrium with each other under expansion. Analogously, in radiative fluid dynamics, the temperatures of the matter and the radiation may not be equal in nonequilibrium processes. This too can be a cause of a bulk viscosity \[Weinberg 1971, Mihalas & Mihalas 1984\], even though photons are ultrarelativistic. Weinberg gave an expression for the radiative bulk viscosity in the Thomas case with only emission and absorption \[Weinberg 1971\]. Subsequent work \[Anderson & Kox 1977, Stewart 1972\], also for media with only emission and absorption but based on other methods, gave different expressions for the forms of the transport coefficients. However, those results have been shown to be equivalent to Weinberg’s in the first order in photon mean free path.

Additional studies have included scattering but, in those more complicated calculations, the bulk viscosity has either not been derived or not extracted explicitly \[Masaki 1971, Hsieh & Spiegel 1976, Straumann 1976, Masaki 1981, Thorne 1981\]. Hence our aim here is to determine the bulk viscosity coefficient from a calculation of the entropy generation in the comoving frame of the matter, once we have obtained suitable approximations for the radiative stress tensor. Since the previous studies have been limited to cases where the mean free paths of photons are much smaller than the prevalent characteristic scales, the problem of treating situations where transparent regions form in the medium has been left open. This issue must be confronted at the edges stars and disks, in intense turbulence and in photon bubbles. A case of particular interest arises in models of the early universe when the mean flight times of photons are comparable to the age of the universe.

The difficulties posed by long mean free paths are less severe for thermal than for dynamical problems since, as we mentioned, the radiative smoothing of temperature fluctuations is reasonably well described by the Eddington approximation, both in optically thick and thin regions. A reason for this is that the Eddington approximation represents a summation of contributions from terms of all orders in a mean free path expansion \[Unno & Spiegel 1966\]. However, since the Eddington approximation leads to a diagonal stress tensor, it fails completely to describe shear viscosity \[Anderson & Spiegel 1972\]. Moreover, it does not give a good representation of the bulk viscosity that opposes changes in volume.

Methods based on truncation of the moment hierarchy at higher moments than the pressure tensor do lead to shear viscosity \[Thorne 1981, Struchtrup 1997\], but they give complicated equations and are difficult to use. Another approach is to use a more general resummation procedure than that leading to the Eddington approximation, one that provides a nonlocal viscous shear tensor. Though this has been done for the case of pure absorption and emission \[Chen & Spiegel 2000\], it has not been carried out with scattering included. Nor is this issue to be considered here. Our purpose in this paper is rather to confront the problem of bulk viscosity for a medium where the photon mean free path is not short and where we also have to deal with scattering processes. For this, we shall adopt a procedure that has recently been introduced in classical kinetic theory \[Chen et al. 2000\] to generalize the usual Navier-Stokes equations. The generalized equations exhibit resistance to volume changes when the particle mean free paths are long, even for a gas of hard spheres.

To test the results obtained with the new method, we have calculated the thicknesses of shock waves and the propagation speeds of ultrasonic waves for classical gases. The results found are in good agreement with experiment, whereas those obtained from the Navier-Stokes equations do not fare well. Though the limits of validity of the new approach are still being probed, the results seem good enough to make the application to the problems of radiative fluid dynamics seem worthwhile. Since the methods that have been used heretofore to compute the radiative dissipation terms parallel the Chapman-Enskog procedure \[Kogan 1969, Cercignani 1988\], we may expect to encounter differences from previous results in this case as well.

In what follows, we describe an expansion in photon mean free path for solutions of the transfer equation, including the effects of Thomson scattering, and develop a stress tensor for use in media whose photon mean free paths need not be small. From this, we compute the rate of comoving entropy generation and provide a formula for bulk viscosity. In a later paper, we shall use this result to estimate the entropy generation in the standard model of cosmology.

## 2 Transfer Theory

### 2.1 The radiation field

We consider radiative transfer in a fluid medium with density $`\rho `$ and velocity $`u^\mu `$, each depending on location in space-time, $`x^\mu `$, where $`\mu =0,1,2,3`$. There are two reference frames of basic interest in this work, an inertial frame, or system frame, and the comoving frame of the matter. In the system frame, we have

$$u^\mu =\gamma (1,𝐯),\gamma =(1𝐯)^{1/2}$$
(1)
while in the comoving frame $`u^\mu =(1,\mathrm{𝟎})`$. Here, as in the following, we assume units in which the speed of light in vacuo and Planck’s constant are each unity. We also adopt the signature in which the Minkowski metric is diag$`(1,1,1,1)`$.

The four-momentum, $`p^\mu `$, of a photon satisfies $`p_\mu p^\mu =0`$, so that we may define a basic null vector,

$$n^\mu =p^\mu /\nu =(1,𝐧),$$
(2)
that characterizes the direction of the photon’s motion in spacetime. With $`\nu `$ as the frequency in the inertial frame, the rest frequency (as seen comoving locally with the medium) is given by the Doppler formula

$$\stackrel{~}{\nu }=u_\mu p^\mu =\gamma \nu (1𝐯𝐧).$$
(3)
The coordinates of the phase space of photons are $`x^\mu `$ and $`p^\mu `$, which we shall sometimes denote as $`x`$ and $`p`$ for brevity. We shall describe the density of photons in phase space by the one-particle distribution function, $`f(p,x)`$. This quantity is a scalar and its integral over the invariant volume in momentum space gives the photon number density in space. To perform this integration, we introduce the invariant volume in phase space, $`dP=\nu d\nu d\mathrm{\Omega }`$, where $`d\mathrm{\Omega }`$ is the element of solid angle \[Landau & Lifshitz 1984\].

A principal quantity of interest in this work is the stress-energy tensor of the radiation field,

$$T^{\mu \nu }=p^\mu p^\nu f𝑑P=n^\mu n^\nu 𝑑\mathrm{\Omega }𝑑\nu ,$$
(4)
where the last term on the right introduces the specific intensity, defined as

$$=\nu ^3f.$$
(5)
A revealing way of writing the stress tensor in terms of basic moments of the radiation field is obtained by defining the directional vector with respect to the moving matter,

$$l^\mu =n^\mu u^\mu ,$$
(6)
so that

$$ł_\mu l^\mu =1\mathrm{and}ł_\mu u^\mu =0.$$
(7)
Hence $`l^\mu =(0,𝐥)`$, where $`𝐥`$ is a unit three-vector in the spatial direction of photon motion in the comoving frame. The angular moments of the radiation field are then

$$=𝑑\mathrm{\Omega },^\mu =l^\mu 𝑑\mathrm{\Omega },𝒫^{\mu \nu }=l^\mu l^\nu 𝑑\mathrm{\Omega },$$
(8)
and their frequency-integrated forms are

$$E=_0^{\mathrm{}}𝑑\nu F^\mu =_0^{\mathrm{}}^\mu 𝑑\nu P^{\mu \nu }=_0^{\mathrm{}}𝒫^{\mu \nu }𝑑\nu .$$
(9)
These moments together make up the energy-momentum stress tensor which may now be written as

$$T^{\mu \nu }=Eu^\mu u^\nu +F^\mu u^\nu +F^\nu u^\mu +P^{\mu \nu }.$$
(10)
### 2.2 The Transfer equation

When both absorption and scattering occur, the transfer equation takes the general form \[Simon 1963\]

$$p^\mu _\mu f=\rho (\alpha \beta f)+\rho \mathrm{}(p,p^{})f(p^{})𝑑P^{}\rho \mathrm{}(p^{},p)f(p)𝑑P^{}$$
(11)
where $`\alpha `$ and $`\beta `$ are scalars characterizing absorption and emission and the kernel $`\mathrm{}(p^{},p)`$ is the differential cross-section for scattering a photon from four-momentum $`p^{}`$ into $`p`$. For brevity, we have not indicated the dependence of $`f`$ on $`x`$ in the transfer equation and have omitted the effects of gravity. We here consider only the case of Thomson scattering for which we have \[Masaki 1971, Tucker 1975\]

$$\mathrm{}(p,p^{})=\frac{3}{4}[1+(𝐥𝐥^{})^2]\delta (\stackrel{~}{\nu }\stackrel{~}{\nu }^{})\frac{\sigma }{4\pi }$$
(12)
where $`\sigma `$ is the Thomson cross-section.

We may factor $`f`$ out of the second integral of (11) and note that with $`i=1,2,3`$ we have

$$l^il^j𝑑\mathrm{\Omega }=\frac{4\pi }{3}\delta ^{ij}$$
(13)
where $`\delta ^{ij}`$ is the Kronecker symbol. Then the last term in (11) reduces to $`\rho \sigma f`$. Since $`f`$ is a scalar, we may use (5) to find the transformation law for the specific intensity. Also, we may compare the covariant form of the transfer equation with the usual one for the comoving frame, as Thomas \[Thomas 1930\] did, to find $`\beta =\stackrel{~}{\nu }\kappa (\stackrel{~}{\nu })`$, where $`\kappa `$ is the absorption coefficient, and $`\alpha =\stackrel{~}{\nu }^2j(\stackrel{~}{\nu })`$, where $`j`$ is the emissivity. Since $`\alpha `$, $`\beta `$ and $`\stackrel{~}{\nu }`$ are scalars, we then know the transformation rules for $`j`$ and $`\kappa `$ once we introduce the (relativistic) Doppler formula (3).

The equation of transfer may now be written as \[Hsieh & Spiegel 1976\]

$$\epsilon n^\rho _{,\rho }=\widehat{\kappa }𝒮+\frac{3\widehat{\sigma }}{16\pi }\left(+l_\rho l_\sigma 𝒫^{\rho \sigma }\right),$$
(14)
where

$$𝒮=\nu ^3\alpha /\beta ,$$
(15)
$$\epsilon =[\rho (\kappa +\sigma )]^1,$$
(16)
$$\widehat{\kappa }=\frac{\kappa }{\kappa +\sigma },\widehat{\sigma }=\frac{\sigma }{\kappa +\sigma }$$
(17)
and a comma with subscript $`\mu `$ indicates differentiation with respect to $`x^\mu `$. We have replaced $`\delta _{ij}𝒫^{ij}`$ by $`\delta _{\rho \sigma }𝒫^{\rho \sigma }`$, to which it is equal, for cosmetic reasons.

Expressions for the moments that appear in the transfer equation may be computed from the equation itself. In doing this, we shall assume that the medium is grey, although we could save appearances by introducing suitable mean absorption coefficients. We shall not follow that practice here since it leads to serious complications in moving media, such as tensorial absorption coefficients, and has not repaid the effort involved.

Next we derive expressions for the monochromatic moments defined in (8) by taking suitable moments of (14). For this purpose, we define

$$h^{\rho \sigma }=\frac{3}{4\pi }l^\rho l^\sigma 𝑑\mathrm{\Omega },$$
(18)
which, in terms of the basic tensors of the problem, may be expressed as

$$h^{\rho \sigma }=\eta ^{\rho \sigma }u^\rho u^\sigma .$$
(19)
Similarly, we have

$$l^\mu l^\nu l_\rho l_\sigma 𝑑\mathrm{\Omega }=\frac{4\pi }{15}(h^{\mu \nu }h_{\rho \sigma }+h_\rho ^\mu h_\sigma ^\nu +h_\sigma ^\mu h_\rho ^\nu ).$$
(20)
Then, when we integrate (14) over all solid angle and note that $`h_{\mu \nu }𝒫^{\mu \nu }=`$, we find that

$$=4\pi 𝒮\frac{\epsilon }{\widehat{\kappa }}n^\rho _{,\rho }𝑑\mathrm{\Omega }.$$
(21)
Further, on applying the second definition in (8) to (14), we get

$$^\mu =\epsilon l^\mu n^\rho _{,\rho }𝑑\mathrm{\Omega }.$$
(22)
Finally, on multiplying (14) by $`l^\mu l^\nu `$ and integrating over angle we find that

$$𝒫^{\mu \nu }=\frac{40\pi \widehat{\kappa }}{3(10\widehat{\sigma })}𝒮h^{\mu \nu }\frac{3\widehat{\sigma }}{10\widehat{\sigma }}h^{\mu \nu }\frac{10\epsilon }{10\widehat{\sigma }}l^\mu l^\nu n^\rho _{,\rho }𝑑\mathrm{\Omega }.$$
(23)
With the moments expressed in this way, we obtain the expression

$$+l_\rho l_\sigma 𝒫^{\rho \sigma }=\frac{16\pi }{3}𝒮\frac{2\epsilon (5+\widehat{\sigma })}{\widehat{\kappa }(10\widehat{\sigma })}n^\rho _{,\rho }𝑑\mathrm{\Omega }\frac{10\epsilon }{10\widehat{\sigma }}l_\mu l_\nu l^\mu l^\nu n^\rho _{,\rho }𝑑\mathrm{\Omega }.$$
(24)
Hence we may write the transfer equation in the compact form

$$=𝒮\epsilon $$
(25)
where we have introduced the linear operator

$$=n^\rho _\rho +c_1𝑑\mathrm{\Omega }n^\rho _\rho +c_2l_\mu l_\nu 𝑑\mathrm{\Omega }l^\mu l^\nu n^\rho _\rho $$
(26)
with $`_\rho =/x^\rho `$,

$$c_1=\frac{3\widehat{\sigma }(5+\widehat{\sigma })}{8\pi \widehat{\kappa }(10\widehat{\sigma })}$$
(27)
and

$$c_2=\frac{15\widehat{\sigma }}{8\pi (10\widehat{\sigma })}.$$
(28)
When there is no scattering, $`\widehat{\sigma }=0`$ and the foregoing equations reduce to those for the case of pure absorption.

## 3 The Radiative Pressure Tensor

### 3.1 The Expansion Procedure

We expand the intensity $``$ as

$$=\underset{m=0}{\overset{\mathrm{}}{}}_{(m)}\epsilon ^m.$$
(29)
From $`_{(m)}`$, we can calculate the coefficients of the expansions of the monochromatic moments and their frequency-integrated counterparts from the definitions

$$_{(m)}=_{(m)}𝑑\mathrm{\Omega },_{(m)}^\mu =_{(m)}l^\mu 𝑑\mathrm{\Omega },𝒫_{(m)}^{\mu \nu }=_{(m)}l^\mu l^\nu 𝑑\mathrm{\Omega }$$
(30)
and

$$E_{(m)}=_0^{\mathrm{}}_{(m)}𝑑\nu F_{(m)}^\mu =_0^{\mathrm{}}_{(m)}^\mu 𝑑\nu P_{(m)}^{\mu \nu }=_0^{\mathrm{}}𝒫_{(m)}^{\mu \nu }𝑑\nu .$$
(31)
When we introduce the expansion (29) into the transfer equation (25) and demand that the expanded equation is satisfied term by term, we get

$$_{(0)}=𝒮,$$
(32)
as in comoving local thermodynamic equilibrium. For the higher orders with $`m1`$, we obtain

$$_{(m)}=_{(m1)}$$
(33)
where $``$ is defined in (26).

The problem studied is really one of a mixture of radiation and matter and a more complete description would begin with coupled transport equations. In the present version, we simply proceed as if the material properties are given, as did Thomas \[Thomas 1930\]. In doing this, we consider that the properties of the medium are expressed, not as functions of $`x^\mu `$, but in terms of the basic fields, $`T`$, $`\stackrel{~}{\nu }`$, $`n^\rho `$. Therefore, in the simplest case, where we operate on a function of $`T`$, $`\stackrel{~}{\nu }`$ and $`n^\rho `$ only, we may write

$$_\mu =T_{,\mu }_T+\stackrel{~}{\nu }_{,\mu }_{\stackrel{~}{\nu }}+n_{,\mu }^\sigma _{n^\sigma }.$$
(34)
On recalling that $`\stackrel{~}{\nu }=u_\mu p^\mu `$ and that $`p^\mu `$ does not depend on coordinate, we obtain

$$\stackrel{~}{\nu }_{,\mu }=p_\nu u_{,\mu }^\nu .$$
(35)
Similarly, on differentiating $`p^\mu =n^\mu \stackrel{~}{\nu }`$ and making use of (35) and the identity $`u^\mu u_{\mu ,\rho }=0`$, we get

$$n_{,\nu }^\mu =n^\mu n^\rho u_{\rho ,\nu }.$$
(36)
So (34) takes the form

$$_\mu =T_{,\mu }_T+\stackrel{~}{\nu }n^\rho u_{\rho ,\mu }_{\stackrel{~}{\nu }}n^\sigma n^\rho u_{\rho ,\mu }_{n^\sigma }.$$
(37)
If we further assume that, like the Planck function, $`𝒮`$ depends on $`\stackrel{~}{\nu }`$ and $`T`$ only through the combination $`\stackrel{~}{\nu }/T`$, we deduce from (37) that

$$𝒮_{,\mu }=\left(T_{,\mu }Tn^\rho u_{\rho ,\mu }\right)_T𝒮.$$
(38)
We note that in the higher orders some of the derived quantities will depend explicitly on the derivatives of $`T`$ and $`u^\rho `$ as well as on the fields themselves. Thus, more generally, when we go to higher order developments, as we did in computing the shear viscosity \[Chen & Spiegel 2000\], we must allow for a dependence on those derivatives. However, here we go only to first order in photon mean free path and that complication does not arise; (37) is therefore sufficient for our purposes.

Returning to the calculation of $``$, we see from (21) that $`_{(0)}=4\pi 𝒮`$ and that, for $`m1`$,

$$_{(m)}=\frac{1}{\widehat{\kappa }}n^\rho _\rho _{(m1)}d\mathrm{\Omega }.$$
(39)
Also, from (22) we have $`_{(0)}^\mu =0`$ and, for $`m1`$,

$$_{(m)}^\mu =l^\mu n^\rho _\rho _{(m1)}d\mathrm{\Omega }.$$
(40)
Then from (23) we find $`𝒫_{(0)}^{\mu \nu }=\frac{4\pi }{3}𝒮h^{\mu \nu }`$ and, for $`m1`$,

$$𝒫_{(m)}^{\mu \nu }=\frac{3\widehat{\sigma }}{10\widehat{\sigma }}_{(m)}h^{\mu \nu }\frac{10}{10\widehat{\sigma }}l^\mu l^\nu n^\rho _\rho _{(m1)}d\mathrm{\Omega }.$$
(41)
### 3.2 The First-Order Stress Tensor

From the recursion relation (33) we have

$$_{(1)}=_{(0)}=𝒮.$$
(42)
Then from (38) and (39) we get

$$_{(1)}=\frac{4\pi }{\widehat{\kappa }}(\frac{1}{3}T\theta +\dot{T})𝒮_{,T}$$
(43)
where

$$\theta =u_{,\mu }^\mu $$
(44)
and $`\dot{T}=u^\mu T_{,\mu }`$. On integrating over frequency, we find

$$E_{(1)}=\frac{4\pi }{\widehat{\kappa }}(\frac{4}{3}S\theta +\dot{S}),$$
(45)
where

$$S=_0^{\mathrm{}}𝒮𝑑\nu =aT^4$$
(46)
and $`\dot{S}=u^\mu S_{,\mu }`$.

In a similar way, we find for the frequency-integrated first-order radiative flux,

$$F_{(1)}^\mu =\frac{4\pi }{3}[S_{,\rho }4S\dot{u}_\rho ]h^{\mu \rho },$$
(47)
while for the frequency-integrated pressure tensor in first order, we obtain

$$P_{(1)}^{\mu \nu }=\frac{40\pi }{10\widehat{\sigma }}\left[\frac{3\widehat{\sigma }}{40\pi }h^{\mu \nu }E_{(1)}+\frac{h^{\mu \nu }}{3}\dot{S}+\frac{4S}{15}(\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }+\frac{5}{3}h^{\mu \nu }\theta )\right]$$
(48)
where

$$\tau ^{\mu \nu \rho \sigma }=h^{\mu \rho }h^{\nu \sigma }+h^{\mu \sigma }h^{\nu \rho }\frac{2}{3}h^{\mu \nu }h^{\rho \sigma }.$$
(49)
We may then introduce (45) and so write the pressure tensor to this order as

$$P^{\mu \nu }=\frac{4\pi }{3}h^{\mu \nu }\left[S\frac{\epsilon }{\widehat{\kappa }}(\frac{4}{3}S\theta +\dot{S})\right]+\mathrm{\Xi }^{\mu \nu }$$
(50)
where

$$\mathrm{\Xi }^{\mu \nu }=\mu \tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }$$
(51)
is the shear tensor and the shear viscosity coefficient is that found by Thomas and Masaki,

$$\mu =\frac{32\pi S}{3\rho (10\kappa +9\sigma )}.$$
(52)
As mentioned, we have elsewhere computed a nonlocal expression for the shear tensor for the case of pure absorption \[Chen & Spiegel 2000\] and expect that a similar generalization may apply in the case with scattering, though we have not yet carried out the necessary calculations to confirm this.

From (47), we may identify the thermal conductivity as the coefficient of $`T_{,\rho }`$:

$$\chi =\frac{16\pi S}{3T\rho (\kappa +\sigma )}.$$
(53)
The second term in the flux, with coefficient $`\xi =T\chi `$, has been considered to be a purely relativistic effect \[Mihalas 1983\], but there is a classical counterpart when the particles have long mean free paths \[Chen et al. 2000\]. We shall see in the next section how these transport coefficients, as well as the bulk viscosity, can be deduced from the entropy generation rate.

## 4 Entropy Generation and Bulk Viscosity

To clarify the roles of the various terms in the radiative stress tensor, we compute their contributions to the entropy generation. In leading order, this is related to the deviation of the stress tensor from its equilibrium value, which we approximate to first order in $`\epsilon `$ as

$$\mathrm{\Delta }T^{\mu \nu }=T^{\mu \nu }T_{(0)}^{\mu \nu }=\epsilon T_{(1)}^{\mu \nu },$$
(54)
where

$$T_{(1)}^{\mu \nu }=E_{(1)}u^\mu u^\nu +F_{(1)}^\mu u^\nu +F_{(1)}^\nu u^\mu +P_{(1)}^{\mu \nu }.$$
(55)
The rate of entropy production is the four-divergence of the entropy flux, $`\mathrm{\Sigma }^\mu `$, whose expression as given by Weinberg \[Weinberg 1972\] and recalled in Appendix A is

$$T\mathrm{\Sigma }_{,\mu }^\mu =\mathrm{\Delta }T^{\mu \nu }\left[u_{\mu ,\nu }u_\mu (\mathrm{log}T)_{,\nu }\right].$$
(56)
For the components of $`T_{(1)}^{\mu \nu }`$ we introduce (45), (47) and (48) and obtain, after some simple rearrangements,

$$\rho (\kappa +\sigma )T\mathrm{\Sigma }_{,\mu }^\mu =E_{(1)}\left[u^\nu (\mathrm{log}T)_{,\nu }+\frac{1}{3}\theta \right]+\mathrm{\Xi }^{\mu \nu }u_{\mu ,\nu }+F_{(1)}^\mu \left[u^\nu u_{\mu ,\nu }(\mathrm{log}T)_{,\mu }\right].$$
(57)
Since $`S=aT^4`$, we can rewrite (45) as

$$E_{(1)}=\frac{16\pi S}{\widehat{\kappa }}\left[u^\rho (\mathrm{log}T)_{,\rho }+\frac{1}{3}\theta \right],$$
(58)
which we use in the first term of (57). To simplify the second term in (57), we note that, on account of the form of $`\mathrm{\Xi }^{\mu \nu }`$, we encounter the expression $`u_{\mu ,\nu }\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }`$. We may use the identity $`h_{\mu \nu }\tau ^{\mu \nu \rho \sigma }=0`$ to show that

$$u_{\mu ,\nu }\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }=\frac{1}{2}(u_{\mu ,\nu }+u_{\nu ,\mu }\frac{2}{3}h_{\mu \nu }\theta )\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }.$$
(59)
Then we find that

$$u_{\mu ,\nu }\tau ^{\mu \nu \rho \sigma }u_{\rho ,\sigma }=\frac{1}{4}(u_{\mu ,\nu }+u_{\nu ,\mu }\frac{2}{3}h_{\mu \nu }\theta )\tau ^{\mu \nu \rho \sigma }(u_{\rho ,\sigma }+u_{\sigma ,\rho }\frac{2}{3}h_{\rho \sigma }\theta ).$$
(60)
All this adds up to

$`\mathrm{\Sigma }_{,\mu }^\mu `$ $`=`$ $`\eta h^{\mu \nu }\left({\displaystyle \frac{T_{,\mu }}{T}}u^\rho u_{\mu ,\rho }\right)\left({\displaystyle \frac{T_{,\nu }}{T}}u^\rho u_{\nu ,\rho }\right)+\zeta \theta ^2`$ (61)
$`+`$ $`{\displaystyle \frac{2\mu }{T}}h^{\mu \rho }h^{\nu \sigma }(u_{\mu ,\nu }+u_{\nu ,\mu }{\displaystyle \frac{2}{3}}h_{\mu \nu }\theta )(u_{\rho ,\sigma }+u_{\sigma ,\rho }{\displaystyle \frac{2}{3}}h_{\rho \sigma }\theta )`$
where

$$\chi =\frac{16\pi S}{3T\rho (\kappa +\sigma )},$$
(62)
$$\mu =\frac{32\pi S}{3(10\kappa +9\sigma )}$$
(63)
and

$$\zeta =\frac{16\pi S}{\rho \kappa T}\left[\frac{1}{3}+\frac{u^\alpha T_{,\alpha }}{T\theta }\right]^2.$$
(64)
To see the meaning of these results, we note that for any four-vector $`v^\mu `$, whose components in the comoving frame are $`(v^0,v)`$, we have

$$h^{\mu \nu }v_\mu v_\nu =\delta ^{ij}v_iv_j=\stackrel{}{v}^20.$$
(65)
Similarly, for any rank-2 tensor $`w^{\mu \nu }`$, we can show that $`h^{\mu \rho }h^{\nu \sigma }w_{\mu \nu }w_{\rho \sigma }0`$. So we find that the rate of entropy generation is never negative as long as all these coefficients are positive.

The first term in (61) is the entropy generation produced by the drag of the radiative flux, which is an effect of radiative thermal conductivity, with due allowance for the acceleration of the matter. The second term represents the entropy generation produced by volume changes and the coefficient in this term, $`\zeta `$, is the radiative bulk viscosity. Finally, the third term is a rate of entropy increase caused by shear, with the coefficient of radiative shear viscosity, $`\mu `$.

We see from (64) that the bulk viscosity depends on $`u^\rho T_{,\rho }=\dot{T}`$. In depending on a rate, our bulk viscosity differs qualitatively from that derived by Weinberg \[Weinberg 1972\] who found a bulk viscosity of the form

$$\zeta _W=\frac{16\pi S}{\rho \kappa T}\left[\frac{1}{3}\left(\frac{P}{e}\right)_n\right]^2.$$
(66)
For a radiation-dominated situation, in which $`\left(\frac{P}{e}\right)_n\frac{1}{3}`$, we get $`\zeta _W0`$.

The principal reason for the difference between the results (64) and (66) is that Weinberg followed the standard practice in transport theory of using a lower order condition (here $`T_{(0),\nu }^{\mu \nu }=0`$) to eliminate the time derivative from the pressure tensor \[Weinberg 1972, Bernstein 1988\]. Thus he obtained $`\dot{S}=4S\left(\frac{p}{e}\right)_n\theta `$ (see Appendix B) and was led to formula (66) for the bulk viscosity. While this does simplify the derivation of (66), it also shrinks its domain of validity. In avoiding the use of this condition in our derivation, we obtain a result that is capable of describing relatively rapid processes and narrow structures. Moreover, we find that the bulk viscous pressure can be large even for the radiation-dominated case.

It is interesting that with both approaches, as the absorption coefficient gets smaller, the bulk viscosity grows larger. This is not so surprising since, for large absorption coefficient, the temperature difference between radiation and matter is small. Still, the tendency toward large bulk viscosity at small absorption coefficient is noteworthy and is reminiscent of the sharply different behaviors of fluids with zero viscosity and those with viscosity tending to zero. Of course, since the effect of Compton scattering will be like that of absorption, the limit of zero absorption strictly does not apply to the general case. Nevertheless, we may expect to find effective heating of expanding media in some situations.

## 5 Conclusion

From its very beginnings, the subject of radiative viscosity has been vexatious, as several of the founders of modern radiative transfer theory — Milne, Eddington, Jeans — found in their attempts to compute the shear viscosity \[Mihalas 1983\]. It was not until Thomas \[Thomas 1930\] did the problem relativistically that agreement was reached on the value of the shear viscosity, to leading order. What is interesting is that the agreement seems to have been reached on purely theoretical grounds since there do not appear to be direct measurements confirming Thomas’s results.

In the matter of bulk viscosity, the situation is even more troubling in that the calculations are subtler and there are by now many contributors to the growing literature on this issue. It is conceivable that disagreements over bulk viscosity may also be settled by common consent, but we feel this is unlikely. As we have noted, in previous work, the simplification of introducing equilibrium approximations into the lower order theory is used, in the spirit of the Chapman-Enskog procedure \[Kogan 1969, Cercignani 1988\]. For gases of classical particles, this approximation leads to the conclusion that a simple classical gas does not have a bulk viscosity. However, when this restrictive approximation is not made, we have found \[Chen et al. 2000\], on using the relaxation (or BGKW) model \[Kogan 1969, Cercignani 1988\], that a classical rarefied gas may have an effective bulk viscosity and that the gas obeys macroscopic equations that differ from the Navier-Stokes equations when the particle mean free paths become long. In the classical case, one has recourse to experiments to test the conclusions and these have lent to support to our conclusion. In the radiative problem, we may hope for a rough empirical check in the photon density of the present universe. This is a calculable quantity, though it is model dependent, and we shall report elsewhere on the results obtained for this case (or see \[Chen 2000\]).

But this is not the whole story. There has been some previous disagreement on the correct form of transport coefficients \[Anderson & Kox 1977\], especially, the bulk viscosity. The calculations are usually done in the frame defined by Eckart and the coefficients are identified by comparison with the general form of the stress tensor. But there are usually two temperatures in the expression of the stress tensor (for matter and radiation) and the difference between them is of the same order as the bulk viscosity. Hence the derived bulk viscosity is dependent on which temperature is actually adopted in the stress tensor. But with our present method, we needed to go only to first order to get good results, so the temperature difference does not figure significantly in the outcome. A calculation of the entropy generation rate then leads to an unambiguous determination of the transport coefficients.

Another way to attack this problem is to use the moment method, which was first formally developed in kinetic theory by Grad \[Grad 1963\]; in transfer theory, it may be traced back to Krook’s \[Krook 1955\] formulation of Eddington’s methods. In modern times, the moment method has been elaborated for use in radiative fluid dynamics by several authors in attempting to deal with viscous effects, especially nonlocal ones \[Anderson & Spiegel 1972, Thorne 1981, Kato et al. 1993, Uday & Israel 1982, Struchtrup 1997\]. However, in the context of classical kinetic theory, the introduction of higher moments beyond those used by Grad does not lead to rapid improvement of the macroscopic description of rarefied gases \[Chen et al. 2000\]. The reason is that the microscopic theory admits only a few slow variables, corresponding to the lowest moments of the distribution function.

Higher moments are fast variables and, as we know from the experience of dynamical systems theory, their introduction does not lead quickly to improved representation of the dynamics. Rather, it is best to attempt to improve the description in terms of the slow variables. Therefore, we have pursued here an approach in which we compute only approximations for the energy density, flux and pressure tensor that may be derived by an expansion in the style of Thomas. The next step in this approach would be to seek a closure relation for the pressure tensor in terms of the radiative variables rather than in terms of a material quantity like the source function $`S`$. We shall take this up in a later paper where we shall see that, to improve our description in terms of the slow variables, we should include their derivatives in the closure relation, since the derivatives of slow variables are also slow variables.

## Appendix A Derivation of entropy generation rate

Consider a medium made up of $`N`$ interacting species in a state of near thermal equilibrium. Let us represent the one-particle distribution of the $`i\mathrm{th}`$ species as

$${}_{}{}^{i}f=^if_0+^if_1$$
(67)
where

$${}_{}{}^{i}f_{0}^{}=exp[(\beta _\mu p^\mu ^i\alpha )^iϵ]$$
(68)
is the equilibrium distribution, with $`{}_{}{}^{i}ϵ=1,1`$ for a Bose-Einstein or a Fermi-Dirac gas, respectively. Also,

$$\beta _\mu =\frac{u_\mu }{T}$$
(69)
and $`{}_{}{}^{i}\alpha `$ is the chemical potential for the $`i\mathrm{th}`$ species, where $`i=1,2,\mathrm{},N`$. The corresponding 4-current and stress tensor are defined as

$${}_{}{}^{i}n_{}^{\mu }^ig\frac{p^\mu }{p^0}^ifd^3p,^iT^{\mu \nu }^ig\frac{p^\mu p^\nu }{p^0}^ifd^3p$$
(70)
while the corrections are, similarly,

$${}_{}{}^{i}n_{1}^{\mu }^ig\frac{p^\mu }{p^0}^if_1d^3p,^iT_1^{\mu \nu }^ig\frac{p^\mu p^\nu }{p^0}^if_1d^3p$$
(71)
where $`{}_{}{}^{i}g`$ is a normalization constant depending on the nature of the statistics of the particles.

The entropy four-vector is

$${}_{}{}^{i}\mathrm{\Sigma }_{}^{\mu }=^ig\frac{d^3p}{p^0}p^\mu [^if\mathrm{ln}^if^iϵ(1+^iϵ)^if\mathrm{ln}(1+^iϵ^if)].$$
(72)
If we put equation (67) into this definition and calculate the entropy generation rate, accurate to second order in $`{}_{}{}^{i}f_{1}^{},`$ we have

$${}_{}{}^{i}\mathrm{\Sigma }_{,\mu }^{\mu }=^i\alpha ^in_{,\mu }^\mu \beta _\mu ^iT_{,\nu }^{\mu \nu }+^iϵ^i\alpha _{,\mu }^in_1^\mu +^iϵ\beta _{\nu ,\mu }^iT_1^{\nu \mu }.$$
(73)
So the total entropy generation rate of the mixture is

$$\mathrm{\Sigma }_{,\mu }^\mu =\underset{i}{\overset{i}{}}\alpha ^in_{,\mu }^\mu \underset{i}{}\beta _\mu ^iT_{,\nu }^{\mu \nu }+\underset{i}{\overset{i}{}}ϵ[^i\alpha _{,\mu }^in_1^\mu +\beta _{\nu ,\mu }^iT_1^{\nu \mu }].$$
(74)
Chemical equilibrium leads to $`_i^i\alpha ^in_{,\mu }^\mu =0`$ and total energy conservation results in $`_i\beta _\mu ^iT_{,\nu }^{\mu \nu }=0`$. With these two constraints, equation(74) reduces to:

$$\mathrm{\Sigma }_{,\mu }^\mu =\underset{i}{\overset{i}{}}ϵ[^i\alpha _{,\mu }^in_1^\mu +\beta _{\nu ,\mu }^iT_1^{\nu \mu }].$$
(75)
With formula (75) we can evaluate the entropy generation rate for radiating media by considering only two species in our mixture. The material medium, is one of these and is assumed to be in thermal equilibrium in the absence of photons. For this component, $`f_1=0`$. The other constituent is the radiation field, whose chemical potential, $`\alpha `$, is zero. So we obtain for the entropy generation formula of the mixture

$$\mathrm{\Sigma }_{,\mu }^\mu =\beta _{\nu ,\mu }T_1^{\nu \mu },$$
(76)
where $`T_1^{\nu \mu }`$ is the correction to the radiative stress tensor. With the definition for $`\beta _\nu `$, equation (76) readily reduces to the formula reported in the text except that, in this appendix, we have not introduced a fiducial $`\epsilon .`$

## Appendix B The relation between $`\dot{S}`$ and $`\theta `$

In order to evaluate the bulk viscosity, we need to express $`\dot{S}`$ in terms of $`\theta `$. This can be done by analyzing the zeroth order entropy of the entire system — medium plus radiation — and using the conservation laws. We have the thermodynamic relation

$$Td\sigma =d(ev)+pdv,$$
(77)
where $`p`$ is the total pressure, $`\sigma `$ is the entropy per particle, $`v`$ is the volume per particle and $`e`$ is the total energy density.

If $`T`$ and $`v`$ are chosen to be the independent variables and $`F`$ is the Helmholtz energy per particle, we can write

$$dF=\sigma dTpdv.$$
(78)
Since $`dF`$ is a perfect differential, we have the Maxwell relation

$$\left(\frac{\sigma }{v}\right)_T=\left(\frac{p}{T}\right)_v.$$
(79)
On writing $`F=evT\sigma `$, we find $`\sigma =\frac{1}{T}(evF)`$. If we put this into the left of the foregoing equation and notice that from (78), $`\left(\frac{F}{v}\right)_T=p`$, we get

$$v\left(\frac{e}{v}\right)_T+e+p=T\left(\frac{p}{T}\right)_v.$$
(80)
We also have the continuity equation

$$\dot{n}u^\mu n_{,\mu }=nu_{,\mu }^\mu $$
(81)
where $`\dot{n}=u^\mu n_{,\mu }`$. In terms of $`v(=1/n)`$, the continuity equation can be written as

$$\dot{v}u^\mu v_{,\mu }=vu_{,\mu }^\mu .$$
(82)
If we assume that the bare matter stress tensor satisfies $`\mathrm{\Theta }_{(0),\nu }^{\mu \nu }=0`$ with $`\mathrm{\Theta }_{(0)}^{\mu \nu }=eu^\mu u^\nu ph^{\mu \nu }`$, we find

$$\dot{e}=(e+p)\dot{v}/v,$$
(83)
where we have used (82). Now we are in a position to calculate $`\dot{T}`$ in terms of $`\theta (=u_{,\mu }^\mu )`$, starting with the identity,

$$e(v,T)_{,\mu }=\left(\frac{e}{v}\right)_Tv_{,\mu }+\left(\frac{e}{T}\right)_vT_{,\mu }.$$
(84)
We have

$$\dot{e}=\left(\frac{e}{v}\right)_T\dot{v}+\left(\frac{e}{T}\right)_v\dot{T}.$$
(85)
Hence

$$\dot{T}=\left(\frac{T}{e}\right)_v\left[\dot{e}\left(\frac{e}{v}\right)_T\dot{v}\right].$$
(86)
On replacing $`\dot{e}`$ from (83), $`\left(\frac{e}{v}\right)_T`$ from (80), and $`\dot{v}`$ from (82) in the above relation, we get

$$\dot{T}=u^\mu T_{,\mu }=T\left(\frac{p}{e}\right)_vu_{,\mu }^\mu .$$
(87)
With $`S=aT^4`$, we find that, in the comoving frame, equation (87) reduces to

$$\dot{S}=4S\left(\frac{p}{e}\right)_v\theta .$$
(88)