text:
# Heat conduction in 1D lattices with on-site potential

## 1 Introduction

Heat conductivity of 1D lattices is well known classical problem related to microscopic foundation of Fourier law. The problem started from the famous work of Fermi, Pasta and Ulam , where abnormal process of heat transfer has been detected in the first time. Non-integrability of the system is necessary condition for normal heat conductivity. As it was demonstrated recently for the FPU lattice , disordered harmonic chain , diatomic 1D gas of colliding particles  and the diatomic Toda lattice , the non-integrability is not sufficient in order to get normal heat conductivity. It leads to linear distribution of the temperature along the chain, but the value of heat flux is proportional to $`1/N^\alpha `$, where $`0<\alpha <1`$ and $`N`$ is the number of particles in the chain. Thus, the coefficient of the heat conductivity diverges in thermodynamical limit $`N\mathrm{}`$. Analytical estimations  have demonstrated that any chain possessing acoustic phonon branch should have infinite heat conductivity in the limit of low temperatures.

From the other side, there are some artificial systems with on-site potential having normal heat conductivity . Heat conductivity of Frenkel-Kontorova chain in the first time was considered in paper . Finite heat conductivity for certain set of parameters was obtained for Frenkel-Kontorova chain , for the chain with sinh-Gordon on-site potential  and for the chain with $`\varphi ^4`$ on-site potential . These models are not invariant with respect to translation and the momentum is not conserved. It was supposed that the on-site potential is extremely significant for normal heat conduction  and that the anharmonicity of the on-site potential is sufficient to ensure the validity of Fourier law . Detailed review of the problem is presented in recent paper .

Peculiarities of the heat conduction of the Frenkel-Kontorova model for complete set of parameters and temperature of the system are not known. The chains with zero average pressure were demonstrated to have normal heat conductivity . In papers  the transition from abnormal to normal heat conductivity has been detected at certain temperature. Paper  contains detailed review of the heat conductivity peculiarities in 1D molecular systems.

There are no detailed investigations similar to mentioned above and concerning the properties of the chains with on-site potential in the whole temperature range. As it was mentioned above, there exists certain incompleteness of the knowledge concerning even the most popular and paradigmatic discrete Frenkel-Kontorova chain. This lattice is of special interest as its continuous counterpart is famous sin-Gordon system (having, of course, divergent heat conductivity). The transition between two regimes with the growth of temperature is expected for discrete system; however it may be dependent also on other parameters of the lattice

The question of special interest is also the mechanism of heat flow scattering which gives rise to finite heat conductivity. For the chain with periodic nearest-neighbor interaction it was demonstrated  that the transition to normal heat conductivity corresponds to abrupt growth of concentration of rotation solitons (rotobreathers), demonstrating certain similarity with phase transition. Namely, the region of the transition temperature corresponds to maximum region of the heat capacity of the lattice. Similarly, it is reasonable to suppose that every lattice with finite heat conductivity has its peculiar mechanism of scattering the heat flow.

The paper is devoted to the detailed simulation of the discrete lattices with on-site nonlinearity and quadratic potential of nearest-neighbor interaction and investigation of their heat conductivity. The lattices are Frenkel-Kontorova, sinh-Gordon, and discrete $`\varphi `$-4. For every case the dependence of the heat conductivity on the temperature and parameters of the lattice will be explored and concrete elementary excitation responsible for the change of regimes will be revealed.

## 2 Description of the model

Let us consider one-dimensional atomic chain arranged along $`x`$ axis. All particles are of equal mass $`M`$, and the nearest-neighbor interaction is described by harmonic potential having rigidity $`K`$. Then the Hamiltonian of the lattice will take a form

$$=\underset{n}{}\{\frac{1}{2}M\dot{x}_n^2+\frac{1}{2}K(x_{n+1}x_n)^2+U(x_n)\},$$
(1)
where the dot denotes the differentiation with respect to time $`t`$, $`x_n`$ is the displacement of the $`n`$-th particle from its equilibrium position and $`U(x)`$ is the on-site potential.

The dimensionless variables are introduced as $`u_n=2\pi x_n/a`$ ($`a`$ is the equilibrium distance between the particles) for the displacement, $`\tau =t\sqrt{K/M}`$ for the time and $`H=4\pi ^2/Ka^2`$ for the energy. Hamiltonian (1) takes the form

$$H=\underset{n}{}\{\frac{1}{2}u_{n}^{}{}_{}{}^{2}+\frac{1}{2}(u_{n+1}u_n)^2+V(u_n)\},$$
(2)
where the stroke denotes the differentiation with respect to the dimensionless time $`\tau `$ and the dimensionless on-site potential is introduced as $`V(u_n)=4\pi ^2U(au_n/2\pi )/Ka^2`$. Natural definition for the dimensionless temperature is $`T=4\pi ^2k_B\mathrm{\Theta }/Ka^2`$, where $`k_B`$ is Boltzmann constant and $`\mathrm{\Theta }`$ is the temperature in common units.

We are going to consider four widely used models for on-site potential: harmonic potential

$$V(u)=\frac{1}{2}\omega _0^2u^2;$$
(3)
sin-Gorgon potential

$$V(u)=ϵ[1+\mathrm{cos}(u)];$$
(4)
$`\varphi `$-4 potential

$$V(u)=2ϵ[(u/\pi )^21]^2$$
(5)
and sinh-Gordon potential

$$V(u)=\omega _0^2[\mathrm{cosh}(u)1].$$
(6)
Parameter $`ϵ>0`$ determines the value of potential barrier between neighboring wells and its inverse $`g=1/ϵ`$ characterizes the cooperativeness of the system. Potentials (4) and (5) have the same distance between neighboring wells equal to 2$`\pi `$ and equal value of the potential barrier 2$`ϵ`$. Parameter $`\omega _0`$ in (3) and (6) corresponds to minimal frequency of harmonic vibrations of the lattice.

## 3 Methods for computation of the heat conduction coefficient

The goal is to simulate the process of heat conduction in finite chain containing N particles. For this sake the left side of the chain ($`n0`$) has to be connected to a thermostat with temperature $`T_+`$, and the right side ($`n>N`$) – to a thermostat with temperature $`T_{}`$ ($`T_+>T_{}`$). For the sake of the simulation we consider the chain of $`N_{}+N+N_+`$ particles, where the first $`N_+`$ particles are attached to the thermostat $`T_+`$, an the last $`N_{}`$ particles – to the thermostat $`T_{}`$ (Fig. 1). The potential of the nearest-neighbor interaction is harmonic, therefore the equilibrium length of the chain does not depend on the temperature. It implies that the boundary conditions at the ends of the lattice has no noticeable effect on the process of the heat conduction and both the conditions of free (Fig.1 (a)) and fixed (Fig. 1 (b)) may be used. Numerical simulations with $`N_\pm =40`$ has demonstrated the absence of the effect depending on concrete choice of the boundary conditions. We’ll use the condition of free ends with $`N_\pm =40`$ for all simulations.

Majority of papers devoted to the topic of the heat conduction  uses deterministic Nose - Hoover thermostat  with $`N_+=N_{}=1`$. However this thermostat has been designed for the description of the equilibrium thermalized system and is not universally suitable for the description of non-equilibrium processes. Therefore our choice is well-known stochastic Langevin thermostat. Detailed comparison of these two thermostats is presented in Appendix 9.3.

Let us consider the chain with free ends ($`1<n<N+N_++N_{}`$) with $`N_\pm `$ particles at both ends attached to Langevin thermostats. The dynamics of the system is described by equations

$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}u_nF(u_n)\gamma u_n^{}+\xi _n^+,`$
$`n`$ $`=`$ $`1,`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n)\gamma u_n^{}+\xi _n^+,`$
$`n`$ $`=`$ $`2,\mathrm{},N_+,`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n),`$ (7)
$`n`$ $`=`$ $`N_++1,\mathrm{},N_++N,`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n)\gamma u_n^{}+\xi _n^{},`$
$`n`$ $`=`$ $`N_++N+1,\mathrm{},N_++N+N_{}1,`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n1}u_nF(u_n)\gamma u_n^{}+\xi _n^{},`$
$`n`$ $`=`$ $`N_++N+N_{},`$
where $`F(u)=dU(u)/du`$, the damping coefficient $`\gamma =1/\tau _r`$, $`\tau _r`$ is the characteristic relaxation time of the particles attached to the thermostat, $`\xi _n^\pm `$ is the random external force corresponding to Gaussian white noise normalized as

$`\xi _n^\pm (\tau )=\xi _n^\pm (\tau _1)\xi _k^{}(\tau _2)=0,`$
$`\xi _n^\pm (\tau _1)\xi _k^\pm (\tau _2)=2\gamma T_\pm \delta _{nk}\delta (\tau _2\tau _1).`$
Details of numeric realization of the Langevin thermostat and random forces are presented in Appendix 9.1.

At every moment the dimensionless temperature of the $`n`$-th particle $`t_n(\tau )=u_{n}^{}{}_{}{}^{2}(\tau )`$. In order to determine the value of the local heat flux $`j_n`$ the energy distribution among the particles of the chain is considered:

$$h_n=\frac{1}{2}\left[\frac{1}{2}(u_{n}^{}{}_{}{}^{2}+u_{n+1}^{}{}_{}{}^{2})+V(u_n)+V(u_{n+1})+(u_{n+1}u_n)^2\right].$$
(8)
By differentiating equation (8) with respect to time $`\tau `$ we get

$$h_n^{}=\frac{1}{2}\left\{u_n^{}[u_n^{\prime \prime }+F(u_n)]+u_{n+1}^{}[u_{n+1}^{\prime \prime }+F(u_{n+1})]\right\}+(u_{n+1}u_n)(u_{n+1}^{}u_n^{}).$$

With account of (7) we obtain

$$h_n^{}=\frac{1}{2}\left[u_{n+1}^{}(u_{n+2}u_n)u_n^{}(u_{n+1}u_{n1})\right].$$
(9)
Taking into account the continuity condition $`h_n^{}=j_nj_{n1}`$ we get the expression for the energy flux: $`j_n=u_n^{}(u_{n+1}u_{n1})/2`$.

System of equations (7) has been integrated numerically. We used the values of $`\gamma =0.1`$, $`N_\pm =40`$, $`N=10`$, 20, 40, 80, 160, 320, 640 and initial conditions corresponding to the ground state of the chain. After the time $`\tau =10^5`$ has elapsed, the end particles achieved thermal equilibrium with the thermostat and stationary heat flux has been formed. Afterwards the dynamics of system (7) has been simulated at the time scale of order $`\tau =10^7`$. The average temperatures of the particles

$$T_n=t_n(\tau )_\tau =\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_0^\tau u_{n}^{}{}_{}{}^{2}(s)𝑑s.$$
(10)
and average values of the heat flux

$$J_n=j_n(\tau )_\tau =\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_0^\tau j_n(s)𝑑s.$$
(11)
were computed for the fragment of the chain between the thermostats.

If the temperature gradient $`\mathrm{\Delta }T=T_+T_{}`$ is small, this method allows avoiding the temperature jumps at the ends of the free fragment of the chain . Characteristic distributions of the heat flux $`J_n`$ and local temperature $`T_n`$ are demonstrated at Fig. 2. At the inner fragment of the chain $`N_+<nN_++N`$ the heat flux is constant $`(J_n=J)`$ and the temperature profile is linear. The coefficient of the heat conductivity is determined using the information concerning the inner fragment of the chain:

$$\kappa (N)=J(N1)/(T_{N_++1}T_{N_++N}).$$
(12)
If the function $`T=\alpha n+b`$ is the best linear approximation for $`T=T_n`$ data at the inner fragment of the chain $`N_+<nN_++N`$, then the value $`\kappa (N)`$ may be calculated with better accuracy by taking $`\kappa (N)=J\alpha `$.

Limit value

$$\kappa =\underset{N\mathrm{}}{lim}\kappa (N)$$
(13)
will correspond to the coefficient of the heat conductivity at temperature $`T=(T_++T_{})/2`$. The question regarding the finiteness of the heat conductivity is reduced to the existence of finite limit (13). I the sequence $`\kappa (N)`$ diverges as $`N\mathrm{}`$ then the chain has infinite heat conductivity at this value of the temperature.

Alternative way to compute the heat conductivity $`\kappa `$ is related to well-known Green-Kubo formula :

$$\kappa =\underset{\tau \mathrm{}}{lim}_0^\tau \underset{N\mathrm{}}{lim}\frac{1}{NT^2}𝐉(s)𝐉(0)ds,$$
(14)
where $`N`$ is already the number of particles in the chain with periodic boundary conditions,

$$𝐉(\tau )=\underset{n=1}{\overset{N}{}}j_n(\tau )$$

is the general heat flux and the averaging $``$ is performed over all thermalized states of the chain. Consequently, the finiteness of the heat conductivity is related to the convergence of the integral

$$_0^{\mathrm{}}C(\tau )𝑑\tau ,$$
(15)
e.g. to the evaluation of the descending rate of the function

$$C(\tau )=\underset{N\mathrm{}}{lim}\frac{1}{NT^2}𝐉(\tau )𝐉(0).$$

Numerically the above autocorrelation function may be found only for finite chain

$$C_N(\tau )=\frac{1}{NT^2}𝐉(s)𝐉(s\tau )_s.$$
(16)
For large enough values of $`N`$ the correlation function $`C_N(\tau )`$ is believed to approximate the function $`C(t)`$ with acceptable accuracy. In order to get stable results the value $`N=4000`$ is usually sufficient. More details concerning the computation of the autocorrelation function are presented in Appendix 9.2.

The methods for computing the heat conductivity coefficient (13), (14) are complementary and allow mutual verification of the results.

## 4 Heat conductivity of the chain with harmonic on-site potential

The chain with harmonic on-site potential (3) is described by linear equations and therefore is completely integrable. The energy transport is performed by non-interacting phonon modes. The heat flux $`J`$ does not depend on the chain length $`N`$, but only on the temperature difference $`\mathrm{\Delta }T`$. Linear thermal profile is not formed. At the inner part of the chain the temperature is nearly constant $`T_n=(T_++T_{})/2`$ (Fig. 3). Therefore according to (12), the heat conductivity coefficient diverges. Correspondingly, the average correlation function $`C(\tau )`$ is constant and integral (15) diverges.

## 5 Heat conduction of the chain with periodic on-site potential

Characteristic features of dynamics of the chain with periodic on-site potential (4) depend on the values of the temperature. As the temperature is small $`Tϵ`$, the on-site potential may be approximated by harmonic single-well potential (3) with $`\omega _0=\sqrt{ϵ}`$. The heat transport is governed by weakly interacting phonons. At the temperature $`Tϵ`$ the chaotic superlattice of topological solitons is formed and the transport properties change drastically. At very high temperatures $`Tϵ`$ the chain is effectively detached from the site and again weakly interacting phonons govern the heat transport. Therefore it is reasonable to investigate the dependence of the heat conductivity on the reduced temperature $`\stackrel{~}{T}=T/ϵ`$.

The behavior of the chain also depends on the cooperativeness parameter $`g=1/ϵ`$. The more the cooperativeness, the less is the density of the soliton superlattice and the phonon scattering effects are less significant. The limit $`g\mathrm{}`$ ($`ϵ0`$) corresponds to completely integrable continuum sin-Gordon equation.

Generally, three limits of discrete Frenkel-Kontorova system correspond to completely integrable systems: at $`\stackrel{~}{T}0`$ the system reduces to the harmonic chain with harmonic on-site potential; at $`\stackrel{~}{T}\mathrm{}`$ – to isolated harmonic chain; at $`g\mathrm{}`$ – to continuum sin-Gordon equation. All these limit systems have diverging heat conductivity. The behavior of the system in the vicinity of these limits is a natural question to be addressed.

Let us start from $`g=1`$ $`(ϵ=1`$, $`\stackrel{~}{T}=T)`$ and investigate the sequence $`\kappa (N)`$ as $`N`$ grows $`(N=10,`$ 20, 40, 80, 160, 320, 640) and different values of $`T`$. As it is may be suggested from Fig. 4 at small ($`T=0.2`$) and large ($`T=200`$) temperatures the heat conductivity coefficient $`\kappa (N)`$ grows as $`N^\alpha `$, at $`T=20`$ – as $`\mathrm{ln}N`$, and at $`T=3`$ converges to finite value $`\kappa =18.5`$. Therefore it may be concluded that at $`T=3`$ the chain has finite heat conductivity. The data related to the other values of the temperature does not allow to draw any conclusions about the behavior of the heat conductivity at larger values of $`N`$. Generally speaking, it may happen that for longer chains $`\kappa (N)`$ will attain certain finite value. Computational tools we use do not allow to investigate higher values of $`N`$. Stil, it is possible to get additional information from the behavior of the autocorrelation function $`C(\tau )`$ at $`\tau \mathrm{}`$.

Numerical simulation demonstrates that for $`T=3`$ the autocorrelation function decreases exponentially (Fig. 5, curve 1). Integral (15) converges and the Green-Kubo formula (14) gives $`\kappa =17.5`$, in good correspondence with $`\kappa =18.5`$ obtained from direct simulation of the heat flux. At $`T=20`$ the autocorrelation function at time scale $`0\tau 800`$ also decreases exponentially (Fig. 5, curve 2). If this trend will persist also for $`\tau >800`$, the Green-Kubo formula will give $`\kappa =77.4`$. It is reasonable to compare this value with the result for $`\kappa (N)`$ presented at Fig. 4 (curve 4). Maximum value of $`\kappa (640)=75.9`$ and no trend towards any finite limit of $`\kappa (N)`$ may be detected. Therefore the likely result is divergence. In order to verify this result the simulation for larger values of $`N`$ (1280, 2560, 5120, 10240) is required, which is beyond our computational possibilities.

The problem for $`T=0.2`$ and $`T=200`$ is even more difficult. The autocorrelation function is presented at Fig. 6. The decrease of the function is very slow and no unambiguous conclusion concerning its character may be drawn out. While extrapolating $`c(\tau )`$ for $`\tau >8000`$ by exponent, the Green-Kubo formula yields $`\kappa =1016`$ for $`T=0.2`$ and $`\kappa =2252`$ for $`T=200`$. In order to get additional information another consuming simulation is required. Still, from the other side, for $`T=200`$ at $`N=640`$ the logarithm of the heat conductivity $`\mathrm{ln}\kappa (N)=7.9>\mathrm{ln}(2252)=7.7`$, and the dependence $`\mathrm{ln}\kappa (N)`$ (Fig. 4, curve 2) does not demonstrate any trend towards convergence. Therefore the most likely result in this case is also the divergence of the heat conductivity.

Let us consider the sequence $`\kappa (N)`$ ($`N=10`$, 20, 40, 80, 160, 320, 640) at other values of the cooperativeness. The results are summarized at Fig. 7. The space of parameters $`(g,\stackrel{~}{T})`$ is divided to two zones denoted by different colors. In the first (gray) zone the sequence $`\kappa (N)`$ converges ($`\kappa (160)\kappa (320)\kappa (640)`$), and in the second (white) zone the sequence grows monotonously. Then, in the first zone Frenkel-Kontorova model has finite heat conductivity, and in the second zone the heat conductivity is either divergent or finite but very high (for $`N640`$ the Fourier law is not valid).

The first zone is limited by certain finite value of $`g`$: for some $`g_0>1`$ and for all $`g>g_0`$ no convergence of $`\kappa (N)`$ was detected. The explanation is that for growing $`g`$ the system becomes closer to continuum integrable sin-Gordon equation. At any fixed $`g<g_0`$ for $`N640`$ the heat conductivity converges only for some finite temperature interval $`0<\stackrel{~}{T}_b<\stackrel{~}{T}<\stackrel{~}{T}_h<\mathrm{}`$. As the cooperativeness decreases ($`g0`$), the upper boundary of this interval tends to infinity ($`\stackrel{~}{T}_h\mathrm{}`$), and the lower boundary tends to zero ($`\stackrel{~}{T}_b0`$) proportionally to $`g`$.

The dependence of $`\kappa `$ on the reduced temperature $`\stackrel{~}{T}`$ is presented at Fig. 8. Within the interval $`[\stackrel{~}{T}_b,\stackrel{~}{T}_h]`$ there exists a critical value $`\stackrel{~}{T}_m`$ corresponding to the minimum of heat conductivity.

In order to reveal the mechanism of the heat conduction it is reasonable to explore the behavior of heat capacity $`c=H/NT`$ ($`H`$ is the average energy of cyclic $`N`$-atomic chain at the temperature $`T`$) on the reduced temperature $`\stackrel{~}{T}`$ (Fig. 9). The heat capacity of classic harmonic chain is always unity, therefore the discrepancy of this value from unity characterizes the significance of nonlinear effects at given temperature. The lattice considered has negative anharmonism and therefore its heat capacity must be more than unity for all temperature diapason. The heat capacity tends to unity as $`\stackrel{~}{T}0`$ and $`\stackrel{~}{T}\mathrm{}`$ and has single maximum at certain temperature $`\stackrel{~}{T}_c`$. This value fairly coincides with the temperature $`\stackrel{~}{T}_m`$, which corresponds to the minimum of the heat conductivity.

Moreover, the increase and decrease of the heat capacity is clearly correlated with the decrease and increase of the heat conductivity respectively. This fact allows suggesting the same physical effects as responsible for both processes. For zero temperature the heat capacity is equal to unity. The increase of the heat capacity at higher temperatures is related to thermal activation of topological solitons (kinks and antikinks) which represent additional degrees of freedom for this system. As a result the dynamical superlattice of solitons appears. The density of this superlattice approaches its maximum at the temperature $`\stackrel{~}{T}_m`$. Further growth of the temperature results in the decrease of the number of degrees of freedom, which is manifested as effective detaching of the chain from the on-site potential. Therefore the heat capacity decreases and tends to unity as the temperature grows.

Correlations between the behavior of the heat capacity and the heat conductivity and especially fair coincidence of $`\stackrel{~}{T}_m`$ and $`\stackrel{~}{T}_c`$ allow us to suppose that the heat transfer is limited by phonon scattering on the soliton superlattice. The effectiveness of such scattering depends on the density of the superlattice as well as on the ability of single kink to scatter phonons. In the strongly cooperative regime $`g>g_0`$ the interaction between solitons and phonons is nearly elastic (close to the case of complete integrability) and therefore the heat conductivity has the trend to diverge. For lower cooperativeness the soliton-phonon interaction is less elastic and finite temperature diapason $`[\stackrel{~}{T}_b,\stackrel{~}{T}_h]`$ of converging heat conductivity appears. For the cases of low $`\stackrel{~}{T}<\stackrel{~}{T}_b`$ and high $`\stackrel{~}{T}>\stackrel{~}{T}_h`$ temperatures the convergence of the heat conductivity cannot be detected in the framework of current experiment. The suggested reason of this effect is that the soliton superlattice effectively disappears.

Let us consider now incommensurate Frenkel - Kontorova chain where the period of the chain is different from the period of on-site potential. The dimensionless on-site potential is periodic function(4) with period $`2\pi `$, and the chain has period $`l=2\pi q`$. Then in system of equations (7) function $`F(u_n)`$ will take the form

$$F(u_n)=\frac{d}{du}U(u_n+nl).$$

For the sake of simulation we choose $`q=l/2\pi =\sqrt{2}`$, corresponding in certain sense to extremely incommensurate case. It is well-known  that such a lattice in its ground state already has soliton superlattice of nonzero density. Therefore the convergence of the heat conductivity is expected to be facilitated as compared to the commensurate case.

Fig. 10 demonstrates the zone in the space of parameters $`(g,\stackrel{~}{T})`$ where the sequence $`\kappa (N)`$, $`N=10`$, 20, 40, 80, 160, 320, 640. converges. For the sake of comparison the boundary for the commensurate case is also presented ($`l=2\pi `$). The result is than no qualitative change of the behavior occurs. The only difference is that the zone with normal heat conductivity moves downwise. This effect is related to presence of superlattice of solitons at any temperature. The transition to normal heat conduction occurs at lower temperature since less solitons should be thermally activated in order to achieve convergence. From the other side, the soliton superlattice facilitates effective detachment of the lattice from the on-site potential (the average coupling energy in the ground state is less) and therefore the upper boundary for the normal heat conduction is also achieved at lower temperatures.

## 6 Heat conductivity of the chain with double-well on-site potential

Let us consider the heat conductivity of the chain with on-site potential $`\varphi `$-4 (5). For this case the analysis of the sequence $`\kappa (N)`$, $`N=10`$, 20, 40, 80, 160, 320, 640, demonstrates that the heat conductivity converges as $`\stackrel{~}{T}>\stackrel{~}{T}_0=3g/2`$ ($`T>1.5`$) - see Fig. 7.

In order to investigate the character of the heat conduction in the temperature range $`\stackrel{~}{T}<\stackrel{~}{T}_0`$ let us consider the temperature behavior of the autocorrelation function $`C(\tau )`$. For $`g=1`$ ($`ϵ=1`$) this behavior is demonstrated at Fig. 11. As $`\tau \mathrm{}`$ the autocorrelation function decreases exponentially. The decrease rate grows as the temperature increases and therefore the conclusion concerning finite heat conductivity at $`\stackrel{~}{T}>\stackrel{~}{T}_0`$ is confirmed. At lower temperatures the decrease rate satisfies a power law $`\tau ^\alpha `$ – see Fig. 12. The degree $`\alpha `$ decreases with the decrease of the temperature. At $`T=1`$ $`\alpha =1.2>1`$, therefore integral (15) converges and the heat conductivity is finite, at $`T=0.5`$ $`\alpha =1.02`$. Within the accuracy available for current numerical possibilities this value corresponds to transition to abnormal heat conduction. It is extremely difficult to obtain reliable data for lower temperatures in order to substantiable this conclusion because of huge computation time required. The reason is that the system is rather close to completely integrable case. New numerical methods based on the latter fact are desirable for investigation of this kind of systems.

The dependence of the heat conductivity $`\kappa `$ on reduced temperature $`\stackrel{~}{T}`$ is presented at Fig. 13. For low cooperativeness ($`g<0.5`$) the heat conductivity approaches local minimum and afterwards local maximum at $`(\stackrel{~}{T}=\stackrel{~}{T}_1)`$ and monotonously decreases to zero as $`\stackrel{~}{T}\mathrm{}`$. The relative value of the maximum decreases as the cooperativeness grows and disappears for certain critical value of $`ϵ`$.

In order to reveal the physical reasons of such behavior of the heat conductivity it is also reasonable to investigate the behavior of the heat capacity $`c`$ (Fig. 9b). As $`\stackrel{~}{T}0`$ the heat capacity $`c1`$. As the temperature grows, the heat capacity grows, achieves its maximum at the temperature $`\stackrel{~}{T}_c`$ and then decreases monotonously to the value less than unity. The value $`\stackrel{~}{T}_c`$ is situated near the maximum point of the heat conductivity $`\stackrel{~}{T}_1`$. Such behavior is related with the peculiarities of $`\varphi `$–4 potential. At low temperatures the main effect is due to negative anharmonism near the ground state (therefore the heat capacity exceeds unity) and for high temperatures ($`\stackrel{~}{T}1`$) the process is governed by positive anharmonism bringing the heat capacity to the value below unity.

Let us now consider the frequency spectrum of vibrations of the chain. The spectrum is computed for $`ϵ=4`$ $`(g=1/4)`$ and three characteristic temperatures $`T=0.4`$, 10, 100. The spectrum of the chain with harmonic on-site potential (3) does not depend on the temperature and has the form

$$E(\omega )=2\omega /\pi \sqrt{(\omega ^2\omega _0^2)(\omega _1^2\omega ^2)},$$
(17)
where maximum frequency $`\omega _1^2=4+\omega _0^2`$. For $`ϵ=4`$, $`\omega _0=4/\pi \sqrt{ϵ}=2.546`$, $`\omega _1=3.238`$. As it is demonstrated at Fig. 14a, for temperature $`T=0.4`$ the spectrum of the chain with on-site $`\varphi `$–4 potential nearly coincides with the spectrum of purely harmonic chain (17). Such spectrum means that at low temperatures only thermalized phonons contribute to the frequency spectrum and other excitations do not play any significant role. For $`T=10>ϵ`$ the distribution crosses the lower boundary of the propagation zone $`\omega _0`$ (Fig. 14b). Such low-frequency component may be associated with intrinsic vibrations of the solitons superlattice. For even higher temperatures $`T=100ϵ`$ the spectrum crosses also the upper boundary of the propagation zone $`\omega _1`$ (Fig. 14c). Such effect may be attributed only to thermalization of high-frequency discrete breathers.

Therefore, for low temperatures $`\stackrel{~}{T}<\stackrel{~}{T}_0=0.5g`$ the dynamics of the system is close to that of harmonic chain. The heat transport is governed by weakly interacting phonons and heat conductivity diverges. For higher temperatures the heat conductivity converges. In the intermediate diapason $`\stackrel{~}{T}_0<\stackrel{~}{T}<\stackrel{~}{T}_1`$ the effective phonon scattering mechanism exists due to the superlattice of topological solitons, and for high temperatures $`\stackrel{~}{T}>\stackrel{~}{T}_1`$ – due to high - frequency discrete breathers. Interplay of two different mechanisms of the phonon scattering explains also the dependence of the heat conductivity on the cooperativeness of the system (Fig. 13). The minimum and maximum of heat conductivity disappear with growth of the cooperativeness since the soliton mechanism of scattering becomes less effective (the soliton-phonon interaction is closer to elastic) and simultaneously the excitation of the discrete breathers becomes easier.

## 7 Heat conductivity of the chain with sinh-Gordon on-site <br>potential

Heat conductivity of this system has been investigated in paper  and here it is reasonable to elucidate the details related to physical mechanisms of the process. The on-site potential (6) is single-well function with positive anharmonism. The analysis of the finite sequence $`\kappa (N)`$, $`N=10`$, 20, 40, 80 160, 320, 640, demonstrates that the heat conductivity converges for high temperatures ($`T>T_0>0`$). This observation is supported by the fact that the autocorrelation function $`C(\tau )`$ at high temperatures for $`\tau \mathrm{}`$ decreases exponentially (Fig. 15), and for low temperatures - by power law (Fig. 16).

The heat conductivity decreases monotonously and for $`T\mathrm{}`$ exponentially tends to zero (Fig. 13, curve 4). Positive anharmonism of the potential leads to monotonous decrease of the heat capacity (Fig. 9, curve 11). The frequency spectrum of vibrations moves towards the upper boundary of the propagation zone with growth of the temperature. These facts allow concluding that the high-frequency discrete breathers provide effective phonon scattering in this model and facilitate the convergence of the heat conductivity. Growing concentration of these breathers with the growth of the temperature leads to monotonous decrease of the heat conductivity coefficient.

Chain with on-site potential

$$V(u)=\beta u^4/4$$
(18)
( positive $`\varphi ^4`$ model) also has finite heat conductivity . Potential (18) as well as sinh-Gordon on - site potential (6) is single - well symmetric function with positive anharmonism. Therefore the mechanism of the phonon scattering is also related to the discrete breathers and $`\kappa (T)0`$ for $`T\mathrm{}`$. For $`\beta =2`$ the heat conductivity $`\kappa (T)T^{1.35}`$ .

## 8 Conclusion

The investigation presented above demonstrates that the anharmonicity of the on-site potential does not constitute sufficient condition for the convergence of the heat conductivity coefficient. The behavior of any concrete model in the above respect depends on its peculiar nonlinear excitations which determine the process of the heat transfer and phonon scattering. Two typical mechanisms of the phonon scattering were revealed in the paper – thermalized soliton superlattice (discrete sin-Gordon and $`\varphi `$–4 models) and discrete high-frequency breathers ($`\varphi `$–4 and sinh-Gordon models). Phonon scattering mechanism may switch with the change of the temperature ($`\varphi `$–4 – model).

For the discrete Frenkel-Kontorova model the zone of the converging heat conductivity for given chain length is limited by low and high temperatures and by high cooperativeness. The numerical possibilities available up to date does not allow establishing unambiguously the character of the heat conductivity outside the zone designated at Fig. 7. Still there is a reason to suggest that infinite chain for certain parameters has diverging heat conductivity, although the zone corresponding to finite heat conductivity will be larger that computed above.

Unlike Frenkel-Kontorova model, for $`\varphi `$–4 model it is possible to demonstrate that for low temperatures the boundary of the transition to abnormal heat conduction may be achieved. It is possible to suppose that there exist a transition from infinite to finite heat conductivity with growth of temperature for any cooperativeness. The probable reason for the divergence to be detectable is the presence of odd-power terms in the on-site potential in the vicinity of the extrema of the potential wells.

The sinh-Gordon model does not allow to detect the divergence of the heat conductivity in current experiment; still, the transition also may be suggested for any cooperativeness.

It is possible to suggest that for any analytic on-site potential for low temperatures the heat conductivity will diverge.

The authors are grateful to Russian foundation of Basic Research (grant 01-03-33122), to RAS Commission for Support of Young Scientists (6th competition, grant no. 123) and to Fund for Support of Young Scientists for financial support.

A.V. Savin is grateful to International Association of Assistance for the promotion of co-operation with scientists from the New Independent States of the former Soviet Union (project INTAS no. 96-158) for financial support.

## 9 Appendix

### 9.1 Numeric realization of the Langevin thermostat

System of equations describing the dynamics of the chain attached to thermostats (7) has been integrated numerically by standard fouth-order Runge-Kutta method with constant step of integration $`\mathrm{\Delta }\tau `$. Numeric realization of delta-function is performed as $`\delta (\tau )=0`$ for $`|\tau |>\mathrm{\Delta }\tau /2`$ and $`\delta (\tau )=1/\mathrm{\Delta }\tau `$ for $`|\tau |\mathrm{\Delta }\tau `$, i.e. the step of integration corresponds to the correlation time of the random forces. That is why in order to get correct description of the Langevin thermostat we must guarantee that the relaxation time $`\tau _r\mathrm{\Delta }\tau `$. In order to fulfill this condition the relaxation time was chosen as $`\tau _r=10`$, and the step of integration for different values of $`N`$, was chosen as $`\mathrm{\Delta }=0.05`$, 0.025, 0.0125.

For every step of integration the random forces $`\xi _n^\pm `$ were taken to be constant. They were computed as independent realizations of the random value $`\xi `$, normally distributed with zero average $`\xi =0`$ and dispersion $`\xi ^2=2T_\pm /\tau _r\mathrm{\Delta }\tau `$. For generating the random value $`\xi `$ program package ZUFALL  was used.

The initial state for the integration of equations (7) was chosen to be equal to ground state of the chain:

$$u_n=u_0,u_{n}^{}{}_{}{}^{}=0,n=1,2,\mathrm{}.,N_++N+N_{},$$
(19)
where $`u_0=0`$ for on-site potentials (3), (6) and $`u_0=\pi `$ for potentials (4), (5). It is convenient to control the accuracy of the simulation through behavior of sequence of average local heat fluxes $`\{J_n\}_{n=N_++1}^{N_++N}`$. If the choice of the integration step $`\mathrm{\Delta }\tau `$ is correct then this sequence should be constant. If the local average heat flux changes from particle to particle then the integration step should be reduced. For growing chain length $`N`$ the step of integration should be also reduced in order to provide sufficient accuracy; the averaging time also grows (see ) and therefore the time of simulation necessary for obtaining reliable results for large $`N`$ turns out to be extremely large.

### 9.2 Computation of the correlation function

In order to compute the autocorrelation function of the heat flux $`C_N(\tau )`$ dynamics of cyclic $`N`$-particle chain was simulated. The thermalized chain with temperature $`T`$ was obtained by integrating Langevin system of equations

$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n)\gamma u_n^{}+\xi _n,`$
$`n=1,2,\mathrm{},N,`$
where $`n+1=1`$ for $`n=N`$ and $`n1=N`$ for $`n=1`$, $`\gamma =0.1`$ (relaxation time $`\tau _r=10`$), $`\xi _n`$ – white Gaussian noise normalized as

$$\xi _n(\tau )=0,\xi _n(\tau _1)\xi _k(\tau _2)=2\gamma T\delta _{nk}\delta (\tau _2\tau _1).$$

System (9.2) has been integrated numerically with initial conditions corresponding to the ground state of the chain. After time $`\tau =10\tau _r`$ the chain approached equilibrium with the thermostat and the coordinates

$$\{u_n(\tau ),u_{n}^{}{}_{}{}^{}(\tau )\}_{n=1}^N$$
(21)
corresponding to the thermalized state at temperature $`T`$.

Afterwards the dynamics of isolated thermalized chain was simulated. For this sake system (9.2) was integrated with zero damping $`\gamma =0`$ and zero external force $`\xi _n0`$). Thermalized state (21) was used as initial condition. The result was the dependence of the general heat flux J on time $`\tau `$. Afterwards with the help of (16) the autocorrelation function $`C_N(\tau )`$ was computed for given thermalized state of the chain. The autocorrelation function depends significantly on concrete realization of the thermalized chain. That is why in order to improve the accuracy this procedure was performed $`10^3÷10^4`$ with independent initial realizations of the thermalized state. Finally the shape of the correlation function was computed as average over all these realizations. It is worth while mentioning that the alternative way of computation (performing of one very large simulation) would not bring about any sufficient gain in the accuracy because of growing integration errors.

In order to verify the independence of the correlation function on the chain length the appropriate calculations were performed for different values of $`N`$. Fig. 16 demonstrates the function $`C_N(\tau )`$ for the chain with sinh-Gordon on-site potential for $`\omega _0=1`$, $`T=2`$ and $`N=500`$, 1000, 2000. It is clear that the autocorrelation function is nearly independent on $`N`$ (the differences are noticeable only for large times and are reduced as the number of realizations used for averaging grows). For given set of parameters $`N=1000`$ provides sufficient accuracy.

### 9.3 Comparison of Langevin and Nose-Hoover thermostats

Unlikely the Langevin thermostat, the Nose-Hoover thermostat (NHT) is not stochastic. Its dynamics is completely determined by the initial conditions. It turned out to be good choice for simulations of FPU system  but its deterministic nature can bring about artifacts in the behavior of the system. We compare this thermostat with the Langevin thermostat (LT) we use for the case of Frenkel - Kontorova model.

Let us consider the chain with fixed ends ($`1<n<N+N_++N_{}`$) with $`N_\pm `$ particles attached to NHT having the temperature $`T_\pm `$. Dynamics of the system is described by equations

$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n)\eta _+u_n^{},`$
$`n`$ $`=`$ $`2,\mathrm{},N_+,`$
$`\eta _{+}^{}{}_{}{}^{}`$ $`=`$ $`{\displaystyle \frac{1}{\tau _r^2}}\left({\displaystyle \frac{1}{(N_+1)T_+}}{\displaystyle \underset{n=2}{\overset{N_+}{}}}u_{n}^{}{}_{}{}^{2}1\right)`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n),`$ (22)
$`n`$ $`=`$ $`N_++1,\mathrm{},N_++N,`$
$`u_n^{\prime \prime }`$ $`=`$ $`u_{n+1}2u_n+u_{n1}F(u_n)\eta _{}u_n^{},`$
$`n`$ $`=`$ $`N_++N+1,\mathrm{},N_++N+N_{}1,`$
$`\eta _{}^{}{}_{}{}^{}`$ $`=`$ $`{\displaystyle \frac{1}{\tau _r^2}}\left({\displaystyle \frac{1}{(N_{}1)T_{}}}{\displaystyle \underset{n=N_++N+1}{\overset{N_++N+N_{}}{}}}u_{n}^{}{}_{}{}^{2}1\right)`$
where $`F(u)=dU(u)/du`$, and $`\tau _r`$ is the relaxation time of the thermostat.

Usually the simulations of the heat conductivity  take $`\tau _r=1`$, and $`N_+=N_{}=2`$ (only end particles are attached to the thermostat $`n=2`$ and $`n=N_++N+N_{}1`$). But, as stated in , such thermostats are not enough random – they cover only a part of the phase space and correspond to strange attractors. In order to reduce this effect we attach to the thermostat $`N_+=N_{}=40`$ particles from every side of the chain.

The dynamics of system (22) is also completely deterministic. It should be mentioned that it is impossible to use the initial condition (19) corresponding to ground state of the system (it is stationary point of system (22)). We take the initial condition

$$u_n=u_0,u_{n}^{}{}_{}{}^{}(0)=4(\xi _n\frac{1}{2})\sqrt{(T_++T_{})/2},$$
(23)
where $`\xi _n`$ – independent realizations of the random variable over the interval .

We choose $`ϵ=1`$ ($`g=1`$), $`T_+=3.05`$, $`T_{}=2.95`$, $`N=80`$ and integrate system (22) numerically with initial condition (23). The distribution of heat fluxes $`J_n`$ and local temperatures $`T_n`$ is presented at Fig. 17 (for the sake of comparison we present also the results received bu using LT - thin lines). Within the left thermostat the heat flux grows linearly and within the other thermostat it decreases linearly with $`n`$. At central part of the chain the value of the heat flux does not depend on $`n`$. Linear temperature profile is formed and the heat conductivity coefficient may be computed according to (12) – $`\kappa (N)=18.4`$. Use of LT gives $`\kappa (N)=18.5`$ (see above), i.e. the value of $`\kappa `$ does not depend on the type of the thermostat.

In addition, it is possible to conclude from Fig. 18 that the frequency distribution of the energy of vibrations also does not depend on the type of thermostat used. It means that for the case of the temperatures close to the value of the potential barrier the choices of NHT or LT bring about equivalent results

The situation is strikingly different if the temperature is lower and the chain is closer to the linear case. The Nose-Hoover thermostat is not effective in this case. In order to illustrate this fact we use the model of harmonic chain. As it is clear from Fig. 3 NHT gives values of the heat flow substantially different from the correct values; at the same times the use of LT secures much better results. That is why in the present paper we used more complicated and consuming LT.

It should be mentioned that sometimes due to its simplicity NHT is used incorporate with LT  (LT is used for the parameters of the model where NHT is not acceptable).


format: text

text:
# On the atmospheric limitations of ground-based submillimetre astronomy using array receivers

## 1 Introduction

Ground-based observations at submillimetre wavelengths are severely hindered by the atmosphere, which absorbs, emits, and refracts the incoming radiation. High, dry sites are needed, such as Mauna Kea in Hawaii. Even then the transmission is generally poor, with only a small number of semi-transparent windows accessible.

To calibrate submillimetre data, the atmospheric opacity must be accurately determined. The main absorber of radiation in this waveband is water vapour, although oxygen and ozone can be significant contributors. Furthermore, the transparency of the atmosphere often changes on short timescales. Thus, frequent measurements of the opacity are crucial, especially for the shorter wavelength windows at 350 and 450 µm. There is a strong dependence of the transmission on wavelength; the shorter windows are more opaque, and deteriorate faster as conditions worsen.

In addition to attenuating the signal, the atmosphere and immediate surroundings of the telescope emit thermal radiation several orders of magnitude larger than the source signal. Spatial and temporal variations in this sky emissivity give rise to ”sky noise”, which can degrade the effective instrument sensitivity by up to an order of magnitude. The thermal DC offset and sky noise variability can be largely removed by the conventional techniques of sky-chopping and telescope nodding. To be effective against sky noise requires the secondary mirror to switch between sky+source and sky faster than the rate at which the sky is varying (typically greater than a few Hz). Nodding the telescope primary to place the source alternately in both chop beams cancels slower varying sky gradients due to chop-beam imbalances and time-dependent telescope spillover signals. It is not practical to nod the primary at the chop rate and so typical frequencies of 0.1 Hz are adopted. However, although these techniques diminish the effects of sky noise, they do not remove the residual signature completely.

Fluctuations in the atmospheric refractive index, as the atmosphere drifts through the telescope beam, cause variations in the path length from source to telescope (Church & Hills 1990; Church et al. 1993). This is sometimes referred to as submillimetre seeing, and causes apparent (short-term) pointing shifts, which ultimately degrade the S/N of an observation. This is particularly problematic for under-sampled arrays such as SCUBA. Until the time that adaptive optics becomes available for the submillimetre, the effects of this seeing must be closely monitored, and the data severely affected removed from a particular observation.

This paper considers the limitations imposed on the quality of submillimetre data by the atmosphere, and presents methods that maximise the accuracy of the data. The layout is as follows: Section 2 provides a technical description of SCUBA and its advantages compared with previous instruments, while Section 3 discusses techniques of determining atmospheric extinction as a function of wavelength. We present our own skydip model and consider how continuously-operating radiometers are crucial for precise measurements of atmospheric opacity. Sky-noise is described in detail in Section 4, particularly how it can degrade instrument sensitivity, and how most of its effects can be removed with state-of-the-art instrumentation such as SCUBA. The effects of refraction noise (submillimetre seeing) are described in Section 5. Using a carefully constructed database of SCUBA observations and atmospheric information, Section 6 analyses how the atmosphere directly effects instrument sensitivity. In Section 7 we conclude with a discussion of future technology that will provide further improvements in the accuracy of submillimetre data.

## 2 The benefits of array receivers: SCUBA

Until recently, observing in the submillimetre was limited to single-pixel devices. UKT14 (Duncan et al. 1990), the forerunner to SCUBA on the JCMT, was detector-noise limited at all wavelengths of operation, except under periods of high sky variability, when the sensitivity at 350 and 450 µm was severely limited by sky-noise.

An array receiver is a major improvement over a single-pixel system. The off-source pixels can be used to measure sky-noise on short timescales. This is crucial for accurate work: Omont et al. (1996) compared the performance of a detector array with a single-pixel system, and found that the poor sky cancellation offered by a single detector could in fact produce fake detections of faint sources.

Within the past few years submillimetre astronomy has witnessed the arrival of the first large format bolometer arrays, together with detector sensitivities which are limited, under stable conditions, by the background photon noise from the sky. The largest and most powerful of this new generation of submm cameras is the Submillimetre Common-User Bolometer Array (Holland et al. 1999), which operates on the 15-m James Clerk Maxwell Telescope, on Mauna Kea. SCUBA is a dual camera system containing 91 pixels in the short-wavelength (SW) array and 37 pixels in the long-wavelength (LW) array. Background-limited performance is achieved by cooling the detectors to $`100`$ mK. Both arrays have approximately the same field-of-view (FOV) on the sky (2.3 arcminutes in diameter) and can be used simultaneously by means of a dichroic beamsplitter. The SW array is optimised for operation at 450 µm, while the LW array is optimised for 850 µm.

The wavelength of operation is selected by a bandpass filter carefully designed to match the transmission window. The filters are multi-layer, metal-mesh interference filters (Hazell 1991) located in a nine-position rotating drum that surrounds the arrays. They have excellent transmission (typically over 80%), and also less than 0.1% out-of-band power leakage. This latter characteristic is particularly important as it ensures that there is minimum contribution to the source signal from extraneous sky emission. The spectral performance of the filters was measured by the University of Lethbridge Fourier Transform Spectrometer (Naylor et al. 1994), and the resultant profiles are overlaid on the Mauna Kea atmospheric transmission curve in Figure 1. For a more detailed look at the submillimetre atmosphere please refer to Naylor et al. (2000).

SCUBA has three basic observing modes (Holland et al. 1999). Photometry, uses a single bolometer to observe point sources. For sources larger than the beam, but smaller than the array FOV, a 64-point jiggle map produces a fully-sampled map for both the LW and the SW arrays, by jiggling the secondary mirror to fill in the undersampled arrays. Finally, sources larger than the array FOV are observed in scan-map mode.

The SCUBA Upgrade Project was designed to improve the sensitivity of the instrument. The first part of the project was completed in October 1999, and included the installation of two wideband filters centred at 450 µm and 850 µm (450W:850W). The wideband filters were designed to be more closely matched to the atmospheric windows (Figure 1), and to be more sensitive than their narrowband predecessors (450N:850N) under all weather conditions. The measured improvement is a few percent at 850 µm and a factor of 2 at 450 µm. The improvement at 850 µm is not due to the width of the filter, but to the blocking filter that was installed to reduce contamination by infrared light. The overall spectral transmission characteristics of the 850-µm waveband are driven by the edge filters and detector feed-horn cut-off. Thus, unfortunately, the measured response of the 850-µm wideband filter is almost identical to that measured for the narrowband filter. For the purposes of this paper, we will concentrate on the 450W:850W filters.

## 3 Atmospheric attenuation

Determining the atmospheric attenuation of a source signal is critical for calibrating submillimetre data. Assuming a plane-parallel atmosphere:

$$I_m=I_{}e^{TauA}$$
(1)
where $`I_m`$ and $`I_{}`$ are the signals incident at the telescope and at the top of the atmosphere respectively, $`A`$ is the airmass (the secant of the zenith distance), and Tau is the zenith sky opacity (Stevens & Robson 1994). Thus, precise measurements of Tau must be taken frequently. This is less crucial at 850 µm; in good weather, Tau850$`<`$0.3, and at low airmass, ($`<`$1.5), a 20% error in Tau850 alters the measured source flux density by 5-10% at most. However, in worse conditions and particularly at 450 µm, an error in Tau can severely affect the measured source flux density. For example, assuming a low airmass $`<`$1.5, a 20$`\%`$ error in Tau450 can alter the measured flux density by $`5080\%`$.

Traditionally, Tau was derived by constructing a secant plot: the signal from the source is measured as a function of airmass, and assuming the sky does not change between measurements, the gradient of the plot gives Tau directly. However, this method requires a meaningful number of measurements of a bright source over a range of airmasses. For any level of accuracy, the sky must remain very stable over a long period of time. Stevens & Robson (1994) demonstrate the difficulties of using this method to derive Tau. We describe here alternative methods that are able to track rapid variations in Tau.

### 3.1 Skydip method

SCUBA estimates the zenith sky opacity at the wavelength and azimuth of observation by performing skydips. Skydips measure the sky brightness temperature as a function of elevation (usually between 80 and 15 degrees), with absolute temperature calibration provided by hot and cold loads. The hot load is ambient temperature Eccosorb, and the cold load is a reflection of the cold optics inside the cryostat, with an effective temperature of $`60`$ K. An aperture plane chopper unit, spinning at 2 Hz, is used to switch rapidly between the sky and the two loads. The temperatures of the hot and cold loads are measured and corrected for the emissivity and reflectivity of the components in the optical path. A model describing both the atmosphere (assuming a plane-parallel form) and the optical system is then fit to the data to calculate the zenith sky opacity:

$$J_{meas}=(1\eta _{tel})J_{tel}+\eta _{tel}J_{atm}bwf\eta _{tel}J_{atm}e^{TauA}$$
(2)
where $`J_{meas}`$ is the measured brightness temperature of the sky, $`\eta _{tel}`$ is the transmission of the telescope, $`J_{tel}`$ is the brightness temperature of a black-body radiating at the temperature of the telescope, $`J_{atm}`$ is the brightness temperature of the atmosphere, $`bwf`$ is the bandwidth factor of the filter being used ($`1bwf`$ is the fraction of the filter bandwidth that is opaque due to atmospheric absorption and, like Tau, is a function of water vapour content), Tau is the zenith sky optical depth and $`A`$ is the airmass of the measurement. Technical details of this model are presented in Appendix A, refer also to Hazell (1991).

Figure 2 shows typical skydips taken with SCUBA at several wavelengths. In practice, the skydip method provides an accurate measurement of Tau. However, it takes $`6`$ minutes to perform a skydip. Given this overhead, it is only practical to perform a skydip every 1.5-2 hours, and quite often the frequency is even less. An on-the-fly skydip mode has recently been commissioned, where data are taken continuously, reducing the time taken to $`23`$ minutes (Coulson 2001).

### 3.2 CSO Tau Monitor

The nearby Caltech Submillimetre Observatory (CSO) operates a 225 GHz (1.25 mm) tipping radiometer, which performs a skydip every 10 minutes, albeit at a fixed azimuth. The precipitable water vapour of the atmosphere is related to the CSO Tau as follows :

$$pwv=20(Tau_{CSO}0.016)$$
(3)
where the $`pwv`$ is in millimetres (Davis et al. 1997).

Comparing skydips taken with SCUBA to those taken at the CSO yields relations between TauCSO & Tau850 and TauCSO & Tau450. If the scatter about these relations is small, the CSO Tau monitor can be used to measure the opacity more frequently than we can with SCUBA, with no additional overhead.

The CSO data show a significant amount of noise, but track long time-scale ($``$2 hour) variations very well. This noise is more than one would expect given the measurement errors quoted in the CSO Tau archive. It is probable that this is due to instrumental noise and does not represent the behaviour of the sky. A polynomial can be fit to the data to track the large-scale variations in CSO Tau as opposed to the (presumably) instrument noise.

In Figure 3, example polynomial fits are presented for two typical nights. Using the relations derived in Section 3.3, we have also plotted the 225 GHz Tau predicted by the 850-µm skydips taken each night. It is striking how well the skydips track the polynomial fits, even when one would imagine that the CSO Tau was moving around too much to be useful.

Producing a composite picture of the night, using the CSO Tau, the SCUBA skydips, and the polynomial fit, provides an additional level of quality control. It can give the observer a much better feel for how the atmosphere was actually behaving on a given night, especially if conditions were apparently unstable. Consider first the lower plot in Figure 3. It is clear that interpolating between the SCUBA skydips would give an erroneous measurement of Tau. In addition, the extreme scatter in CSO Tau indicates that the end of the night was unusable.

The composite picture can provide further information about when data should be treated with care. For example, on some nights the CSO Tau and SCUBA skydip data disagree with each other, on other nights both display unusually high levels of scatter, indicating an inherently unstable night. If, on the other hand, the CSO Tau has a high-level of scatter but the SCUBA skydips follow the polynomial fit, the scatter is unlikely to be representative of the sky itself.

### 3.3 Tau relations

The time-resolution of the polynomial fits to the CSO Tau is impressive, $``$ 2 minutes if the residuals are small, as they typically are. We can capitalise on this by deriving tight relations between the CSO Tau and the opacity measured at SCUBA wavelengths.

These relations have been constructed taking the following points into account:

1. Skydips were discarded where the model failed owing to, for example, the atmosphere losing its plane-parallel nature or a cloud drifting overhead. In these cases the model returns unrealistic values for the fit parameters and/or unusually high fit residuals. This is a more common occurrence at 450 µm, where the shape and height of the atmospheric window are more sensitive to the presence of water vapour (see Figure 1). We discarded $``$20$`\%`$ of the data at 850 µm and $``$50$`\%`$ of the data at 450 µm. Note, we can still calibrate at 450 µm when the skydip fit fails as we can extrapolate from either the CSO or the 850-µm Tau.
2. We ignored data taken when the CSO Tau monitor was off-line.
3. If either the CSO Tau monitor or the skydip indicated a non-physical value of Tau, i.e. negative or zero, the observation was ignored.
4. The datasets were restricted to TauCSO$`<0.2`$ (SCUBA is not used in weather conditions worse than this).
5. A model of the form $`Tau_Y=a(Tau_Xb)`$ was fit to the data to derive the Tau relations.

Figure 4 shows how much the scatter in the relations has been reduced by simply using the polynomial fits to estimate the CSO Tau and ignoring untrustworthy skydips. The final relations are presented in Table 1 and Figure 5. They display relatively little scatter, even in the poorest weather conditions, and can thus provide accurate submillimetre calibration.

Comparing the wideband and narrowband filters, the wideband CSO Tau relations are steeper at 450 µm but are almost identical at 850 µm. At first sight, the difference at 450 µm is perhaps unexpected given the lower central wavelength of the wideband filter (at 850 µm, the narrow and wideband filters have almost identical central wavelengths). However, the 450 µm wideband filter includes an H2O line that is always present, even in very dry weather conditions. Thus, the wideband filter will always have a steeper slope than the narrowband filter.

If the submillimetre opacity is due solely to water vapour, the Tau relations for the different filters are expected to intercept the origin: if there is no water vapour in the atmosphere, TauCSO, Tau850 and Tau450 will all equal zero. However, we have found evidence for non-zero intercepts. Assuming the straight-line model can be extrapolated to the intercept, this could be explained as follows: the relative contribution of ozone vs. water vapour to absorption is greater at 225 GHz than at 850 µm and in turn greater at 850 µm than at 450 µm. The ozone contribution is relatively invariant for long periods of time, and these offsets should be constant.

## 4 Sky-noise

### 4.1 Signatures of sky-noise

The atmosphere emits submillimetre radiation several orders of magnitude larger than the signals we are trying to measure. Variations in this emission result in sky-noise, which manifests itself as a 1/f component in the noise spectrum. This is illustrated in Figure 6, where the 1/f component extends out to 4 Hz and is well above the system noise level. Two-position chopping at higher frequencies, $`8`$ Hz, greatly reduces this sky-noise contamination (Duncan et al. 1995).

As noted earlier, telescope nodding is also essential for the removal of sky-noise. This can be illustrated by the Allan variance, which is simply the variance of the difference of two contiguous measurements (e.g., Allan 1966; Schieder & Kramer 2001). When plotted, a slope of $`1`$ indicates white noise, a slope of 0 denotes 1/f noise, and a slope of +1 is drift noise. Figure 7 shows the Allan variance data as a function of time for a typical observation with and without nodding (in both cases the telescope was chopping at $`8`$ Hz). A clear 1/f component is present if the telescope is not nodded: after around 20 seconds the noise stops integrating down. This component disappears when nodding is employed.

However, it is impractical to nod the telescope at the chop frequency. Thus, even with a chop/nod configuration, short-term temporal variations in the sky emissivity, as well as spatial effects caused by chopping through slightly different atmospheric paths, still exist. This residual sky-noise must be removed. There are several ways to do this; filters can be designed to select the most transparent parts of the atmospheric transmission window. The optical throughput should be single-moded, i.e. the minimum required to couple to a point source with maximum spatial resolution and minimum background (Duncan et al. 1995). Keeping the chop throw as small as possible, and in an azimuthal direction (i.e. chopping through the same atmospheric layer), should also help. In this paper we will concentrate on the benefits of having a bolometer array to aid the removal of sky-noise.

### 4.2 Removal of sky-noise

In the submillimetre, the emissivity variations occur in atmospheric cells that are larger than the array field-of-view (Jenness et al. 1998; Borys et al. 1999). Therefore, the noise should be correlated to a large extent across the array and also between wavelengths.

One of the main advantages of SCUBA is that it contains an array of bolometers. When observing a compact source, the off-source bolometers can be used to measure the sky signal on short timescales and hence remove any residual sky-noise that chopping and nodding have failed to account for.

This is clearly illustrated in Figure 8, where the signal of both the centre on-source 850-µm pixel and the sky (measured by averaging the outer ring of bolometers) are shown. The signals are highly correlated, and subtracting the sky from the source reveals a clear positive signal (in which the jiggle modulation is apparent). Using this sky-subtraction method can improve the S/N significantly. In Figure 8 we also display the 450-µm centre pixel, and, as expected, it is also highly correlated with the 850-µm data.

Failure to remove this residual sky-noise can, in the worst conditions, degrade the sensitivity of the instrument by an order of magnitude. Figure 9 displays the noise integrating down with time for a deep 850-µm photometry observation. Data are shown for the standard chop/nod configuration, and also after the residual sky-noise has been removed using the off-source bolometers. We define the sensitivity, also known as the noise equivalent flux density (NEFD), as the noise reached in 1 sec of integration. It is clear that the NEFD of the instrument is significantly better if the residual sky-noise is removed. Furthermore, at the end of the observation, there was a very obvious, large deterioration in the sky conditions. This is plainly visible in the case where the residual sky-noise was not removed. Removing the residual takes this change of conditions into account, and allows the noise to integrate down $`\sqrt{t}`$.

For scan-map observations, removing the sky-noise is more difficult, as every bolometer could be observing either the source or the sky at a given time. The source can be removed from the datastream by the simple assumption that the source structure is constant over time, while the sky is varying on timescales of a few seconds. The sky-emission noise can then be calculated and removed from the data (Jessop et al. in preparation).

### 4.3 Investigation of sky-noise

For the purposes of this paper, we define sky-noise as the standard deviation of the sky signal, with any errors associated with nodding removed. The dataset used here consists of deep photometry observations only, as it is easy to specify the sky bolometers and the long observations provide good statistics.

The relationship between sky-noise, seeing, and sky opacity is not immediately obvious; the scatter in the data was too large to reveal an obvious trend. To overcome this, we binned the data. This also posed problems as the data were not Gaussian distributed and it was not clear whether the mean (which is skewed by outliers) and associated standard deviation (which is unnaturally large given the scatter in the data) were useful quantities to measure. Instead we measured the mode and full-width half-maximum of the data in each bin. These plots are shown in Figure 10, and although the ‘errors’ are somewhat large, there is evidence for positive correlations between sky-noise and seeing and between sky-noise and Tau. There is, however, no evidence to support the long-standing belief (from single-pixel instruments) that sky-noise often increases in very dry weather conditions.

Considering the 850-µm and 450-µm sky signals for 20 minute photometry observations, Figure 11 presents a histogram of the Spearman-rank correlation coefficients. A value of +1 indicates a pure positive correlation, 0 indicates no correlation, and -1 indicates a pure negative correlation. On the whole, the 850-µm and 450-µm data appear to be highly correlated. However, there are clearly times when the strength of the correlation is low.

Given the known correlation between sky-noise at 850 µm and 450 µm, it has been suggested that the short wavelength array could be used to provide sky cancellation at 850 µm. This would be extremely useful for faint sources where the 450 µm field-of-view is likely to be source-free, but the 850 µm field is expected to contain several sources or extended structure.

The ratio of the 450-µm sky signal to the 850-µm sky signal (corrected for extinction) is displayed in Figure 12. For the narrowband data, the median ratio is $`4.5`$ in agreement with Borys et al. (1999) and Jenness et al. (1998); for wideband data, the median ratio is $`8`$ and the distribution is wider. The ratio shows less dispersion if it is extinction corrected than if not, suggesting that the water vapour creating the sky-noise is higher up in the atmosphere rather than just above the telescope (where zero correction would be more suitable).

In Figure 13 we present the ratio of 450-µm:850-µm sky-noise as a function of the corresponding transmission ratio. Thus, for a given transmission ratio, it is possible to use this plot to correct 850-µm data for sky-noise using 450-µm sky-noise data. It is clear that the sky-noise is larger for the wideband filter set than for the narrowband. This difference is to be expected as the 450-µm wideband filter is double the bandwidth of its narrowband predecessor, and contains an $`\mathrm{H}_2\mathrm{O}`$ absorption line, while the two 850-µm filters have almost identical bandwidths (Duncan 1983 showed analytically that for a sky-limited system $`NEP\mathrm{\Delta }\nu `$). This excess sky-noise is the cost of the improvement in sensitivity offered by increasing the 450-µm bandwidth. However, as we have already demonstrated, this excess is removable with a bolometer array.

A transmission ratio of 1 in Figure 13 corresponds to there being no atmosphere between the telescope and the source of the sky-noise. If the sky-noise ratios in Figure 13 are extrapolated to a transmission ratio of 1, we find values of $`3.5`$ and $`6.5`$ for the narrowband and wideband filters correspondingly. These numbers are in good agreement with the theoretical values derived in Appendix B. The extinction-corrected ratios in Figure 12 are larger however, indicating that we are over-correcting for extinction in calculating these ratios. This in turn suggests that the water vapour responsible for the sky-noise lies somewhere in the middle-upper layers of the atmosphere, not at the very top of the atmosphere itself.

## 5 Submillimetre seeing

### 5.1 Measurement on Mauna Kea

Submillimetre seeing arises from variations in the refractive index of the atmosphere, principally due to the passage of water vapour through the beam. This variation is measured on Mauna Kea by a phase monitor operated by the Smithsonian Astrophysical Observatory (SAO). The SAO device comprises two 1.8 m dishes placed about 100 m apart located near the JCMT. The dishes point low in the eastern sky, and detect a signal at 12 GHz from a geostationary satellite. The difference in path lengths results in a phase difference between the two received signals that changes slowly with the oscillation of the satellite position. Turbulence in the atmosphere adds noise to the phase difference. The phase difference is measured each second and one minute’s worth of data are analysed for the rms scatter to provide the seeing at 12 GHz over a 100 m baseline.

The conversion from the rms of these phase fluctuations to an rms seeing for the 15 m JCMT is suggested by Masson (1991) as being:

$$\mathrm{JCMT}\mathrm{seeing}(\mathrm{arcsec}\mathrm{rms})=0.5\mathrm{SAO}\mathrm{seeing}(\mathrm{degrees}\mathrm{rms})$$

Fifteen minute averages of these seeing data are available at the JCMT as part of the Telescope Management System (Tilanus et al. 1997). By the middle of the night the JCMT seeing typically drops to 0.2 arcsec, and observers consider seeing $`>1`$ arcsec to be ‘poor’.

### 5.2 Impact of seeing: broadening & displacement

If atmospheric turbulence on scales smaller than the telescope diameter is significant then beam-broadening will occur. However, experiments by Church & Hills (1990) indicate that these effects are small at the JCMT. Beam displacement results from large-scale, anomalous refraction, or ‘tip/tilt’ effects, and effective beam-broadening results if observations take longer than the timescale of these fluctuations; the time average of the beam motion being superimposed on the beam profile. Anomalous refraction timescales are of order 1 s (Olmi 2001), which is considerably shorter than standard JCMT observations (minimum 18 s), but comparable with the SAO monitor rate. JCMT jiggle maps (integration times $`>32`$ s) therefore suffer both minimal average pointing shifts and minimal image broadening as the turblence tends to occur on scales larger than the telescope diameter.

Additionally, analysis of JCMT tracking data shows that there is no correlation between seeing and pointing excursions. Tracking data are long jiggle maps of bright point sources, and each integration (32 s worth) is analysed for the location of the image centroid to yield information on pointing offsets as a function of azimuth, for instance. Some tracking datasets cover a sufficiently long time interval that significant changes in seeing occur. Within such data there is no correlation between the change in seeing and excursions in pointing.

Tracking data also tend to validate the conversion factor (0.5) between the SAO phase (rms, degrees) and JCMT-seeing (rms, arcseconds), since the scatter of the pointing offsets about some low-order polynomial fit is usually of the same size as the JCMT-seeing.

### 5.3 Do errors in computed refraction generate pointing errors?

Relatively thin atmospheric layers (thicknesses $`<100`$ m, say) with arbitrary values of temperature, pressure, and humidity have no impact upon the overall refraction: the refraction towards the normal of a ray entering such a layer is nullified by the refraction away from the normal upon exit. Image displacement (a change in pointing) must therefore result from either non-linear atmospheric structures or from uncompensated changes in local humidity.

The latter hypothesis is simple enough to test using archived JCMT/SCUBA pointing data, and no such relationship is found. The formulae in Appendix C.1 show that, to first order and at constant zenith distance, elevation pointing corrections of $`0.0768(h20)`$ arcseconds need to be applied to account for the impact of local humidity, $`h`$(%). In the absence of this correction a plot of the change in elevation pointing against the change in humidity would be expected to have a slope of 0.0768, whereas a plot of the same data perfectly corrected for humidity will have a slope of zero. Analysis of 3000 pairs of consecutive pointing measures yields a slope of $`0.009\pm 0.018`$, commensurate with the latter scenario. The source of residual pointing shifts therefore would seem to be excursions from the idealized atmospheric model - i.e. turbulence - which is probably no great surprise.

### 5.4 Impact of pointing errors upon flux measurements

Whatever their origin, pointing errors of size $`\theta `$ (arcseconds) cause a signal of strength $`S_o`$ to be measured as

$$S=S_oexp\left[4ln(2)\left(\frac{\theta }{FWHM}\right)^2\right]$$
(4)
where FWHM is the full width at half-maximum of the Gaussian beam. These are shown graphically in Figure 14 in the cases of the JCMT operating at 850 µm and 450 µm.

As to the distribution of pointing errors as a result of refraction noise, we have determined that the expected signal, $`<S>`$, is related to the true peak signal, $`S_o`$, by

$$\mathrm{S}_\mathrm{o}/<\mathrm{S}>=1+0.0275\sigma ^2\mathrm{for}\mathrm{a}14.2^{\prime \prime }\mathrm{beam}\mathrm{at}850\mu \mathrm{m}$$

$$\mathrm{S}_\mathrm{o}/<\mathrm{S}>=1+0.1132\sigma ^2\mathrm{for}\mathrm{a}7^{\prime \prime }\mathrm{beam}\mathrm{at}450\mu \mathrm{m}$$

where $`\sigma `$ is the rms refraction noise. The full derivation of these equations can be found in Appendix C.2.

## 6 The effect of the atmosphere on sensitivity

We have recently made a concerted effort to characterise SCUBA’s performance. One aspect of this work was the creation of a database of calibration observations. The key result from studying this data is that the flux conversion factors (FCFs) are quite stable, to an accuracy of 5% at 850 µm and 20% at 450 µm (Jenness et al. 2001). Using this database, we have created a homogeneous dataset of $`4000`$ calibrated photometry observations, to better understand the effect of the atmosphere on sensitivity.

Figure 15 presents the NEFD as a function of normalised sky transmission for a chop/nod configuration, and also using the off-source bolometers to remove the residual sky-noise. The NEFD is significantly lower and is inherently more stable if the off-source bolometers are used for sky-cancellation. Data are shown both for the wideband filters, and the pre-upgrade narrowband filters. It is worth noting that the 450-µm wideband filter (which is better matched to the entire atmospheric window but includes an H2O line) is considerably more sensitive under all conditions. There is little difference between the narrow and wideband 850-µm filters, which is to be expected because, as noted in Section 2, the overall filter transmission profiles are almost identical as measured in SCUBA.

We have also considered how the sensitivity of the instrument is affected by the chop throw. Figure 16 presents the data for two cases: a standard chop/nod method of sky cancellation, and a chop/nod method where the off-source bolometers are used to remove residual sky noise. Although the error bars are somewhat large, the standard chop/nod method seems to support the long-held belief that sensitivity decreases with chop throw (e.g. Church et al. 1993; Duncan et al. 1995). However, when using the array and when off-source bolometers are used for sky cancellation, large chop throws do not degrade the sensitivity. In this case, the NEFD vs. chop throw plot is perfectly flat out to a chop throw of 180 arcsec. The only variation is due to an increase of the FCF with chop throw, of the order of 10% between 45 and 120 arcsec (Jenness et al. 2001).

## 7 Summary and the future

As described in this paper, we have a good grasp of the atmospheric limitations in the submillimetre and how best to overcome them. Our main findings can be summarized as follows:

1. Tau relations between the CSO Tau at 225 GHz and both the 850 µm and the 450 µm SCUBA filters. The relations display relatively little scatter.
2. For accurate sky cancellation, it is essential to use off-source bolometers to monitor and remove the sky-noise.
3. There is evidence for positive correlations between sky-noise and seeing, and sky-noise and sky opacity,but the correlation coefficient varies significantly depending on the dataset Figure 11.
4. 850-µm and 450-µm sky-noise are clearly correlated, but there are times when the correlation is low.
5. The JCMT beam is not significantly broadened by seeing.
6. There is no obvious correlation between seeing and pointing excursions.
7. If off-source bolometers are used for sky cancellation, chopping as far as 180 arcsec (in any direction) does not affect the sensitivity of the instrument.

In the future, as our understanding grows, as technological advances are made, and as the next generation of submillimetre telescopes are destined to higher and drier sites, we will be able to do even better. Some things to look forward to are:

1. Line-of-sight radiometers which give real time estimates of Tau in the direction you are observing. One of these radiometers is now in operation at the JCMT (e.g. Wiedner et al. 2001, Phillips et al. in preparation).
2. Instantaneous estimates of the line-of-sight opacity should be possible using a hot load, cold load, and sky observation at the elevation of the source, given a database of sky temperature versus elevation information (e.g. Smith, Naylor & Feldman 2001 ).
3. High-speed data sampling leading to a more efficient observing mode called ‘DREAM’ (Le Poole & van Someren Greve 1998). In DREAM mode, the functions of chopping and jiggling the secondary mirror are combined to a single step action, eliminating the need to sample empty sky for half the time. Another gain from high-speed sampling is real-time suppression of sky-noise.
4. Adaptive optics for the submillimetre. The Large Millimetre Telescope (LMT) is currently designing a radiometric wave-front sensor which will measure the tilt of the incoming wavefront to compensate for atmosphere-induced pointing errors (For further details refer to the LMT website: http://www-lmt.phast.umass.edu/).
5. DC-coupled fully sampled arrays such as SCUBA 2 (Holland et al. 2000; Robson et al. 2001), which would remove the necessity to chop and nod for sky-noise removal. This is an advantage: there is no chance of chopping onto a nearby source, and there will be a sensitivity gain as all the time is spent looking at the source instead of half the time looking at the sky. In addition, chopping limits the size-scales we see in maps to be no more than a few times the chop throw. Finally, image reconstruction techniques for two-beam chopping tend to propagate noise, and so with only a single-beam on the sky this problem will be minimsed.
6. In the future, submillimetre telescopes will be built at sites where the atmosphere is largely transparent, for example Chajnantor in Chile at 17,000 ft. (ALMA) and the South Pole. The extreme cold at the South Pole results in very small amounts of water vapour in the atmosphere, and thus sky-noise is much lower than other sites (Stark 2001).
7. Interferometers such as ALMA have a unique advantage: the atmospheric signals decorrelate, and therefore sky-noise should not be a problem.

## ACKNOWLEDGMENTS

We wish to thank the Canadian co-op students who have contributed to this work: Ed Chapin, Jeff Wagg, and Karl Kappler. We also wish to thank David Naylor and his group for the SCUBA FTS measurements of the filter profiles. This work has made use of the Tau archive maintained by the Caltech Submillimetre Observatory, and the phase monitor operated by the Smithsonian Astrophysical Observatory. We acknowledge the support software provided by the Starlink Project which is run by CCLRC on behalf of PPARC. The JCMT is operated by the Joint Astronomy Centre, on behalf of the U.K. Particle Physics and Astronomy Research Council, the Netherlands Organisation for Pure Research, and the National Research Council of Canada. Further details and updated information can be found on the SCUBA World-Wide Webpage at URL: http://www.jach.hawaii.edu/JCMT/Continuum$`\mathrm{\_}`$observing/
continuum$`\mathrm{\_}`$observing.html.

## Appendix A Skydip Model

When performing a skydip, SCUBA measures the sky brightness temperature as a function of airmass. We present a multi-layer model of the atmosphere (Hazell 1991) which, when compared with the data, yields the zenith sky opacity. The model at each wavelength takes the form:

$$J_{meas}=(1\eta _{tel})J_{tel}+\eta _{tel}J_{atm}bwf\eta _{tel}J_{atm}e^{TauA}$$
(5)
where $`J_{meas}`$ is the measured brightness temperature of the sky, $`\eta _{tel}`$ is the transmission of the telescope, $`J_{tel}`$ is the brightness temperature of a black-body radiating at the temperature of the telescope, $`J_{atm}`$ is the brightness temperature of the atmosphere, $`bwf`$ is the bandwidth factor of the filter being used ($`1bwf`$ is the fraction of the filter bandwidth that is opaque due to atmospheric absorption and, like Tau, is a function of water vapour content), Tau is the zenith sky optical depth and $`A`$ is the airmass of the measurement.

Of these parameters, $`J_{meas}`$, $`J_{tel}`$ and $`A`$ are known. $`J_{atm}`$ can be estimated from the ambient air temperature at ground level using a model for the behaviour of the observing layer above the telescope, as described below. $`\eta _{tel}`$ may be fitted to the data for every skydip and, because it does not vary with atmospheric conditions, a reliable ‘average’ value can be derived from many observations. Thus, there are two remaining free parameters, $`\tau `$ and $`bwf`$, that must be derived from the fit.

$`J_{atm}`$ is calculated from $`T_{amb}`$, the ambient air temperature, by assuming that the sky emission is dominated by a single absorber/emitter whose density falls exponentially and temperature linearly with height. In this case it can be shown that:

$`J_{atm}=`$ $`J_{amb}{\displaystyle _0^{40}}Akexp({\displaystyle \frac{h}{h_2}})\times `$ (6)
$`exp\left[Akh_2\left(\mathrm{exp}\left({\displaystyle \frac{h}{h_2}}\right)1\right)\right]\left(1{\displaystyle \frac{h}{h_1}}\right)dh`$
where $`h_1`$ is $`J_{amb}/6.5`$ to give a 6.5 K fall in temperature per km height, $`h_2`$ is the scale height of the absorbers (2 km), $`A`$ is the airmass and $`k`$ the extinction per km.

If we approximate the result of the integral by:

$$J_{atm}=J_{amb}X_g\left[1\mathrm{exp}\left(Akh_2\right)\right]$$
(7)
it can be shown that $`X_g`$ has the form:

$$X_g=1+\frac{h_2T_{lapse}}{T_{amb}}\mathrm{exp}\left(\frac{ATau}{X_{gconst}}\right)$$
(8)
where $`T_{lapse}`$ is the temperature drop per kilometre altitude ($`6.5`$K/km) and $`X_{gconst}`$ is a constant determined empirically and has a value of 3.669383.

## Appendix B Theoretical derivation of sky-noise ratio

If we consider two lines of sight separated by an angle theta and identical to each other except for the fact that one line of sight has a cloud of water vapour at a distance to which the line of sight opacity is $`\tau _{cloud}`$, and the cloud itself has optical depth $`\delta \tau `$, then the difference in surface brightness between the two lines of sight (measured in units of brightness temperature) is simply

$$\delta J=exp(\tau _{cloud})\delta \tau (J_{cloud}J_{inc})$$
(9)
where $`J_{cloud}`$ is the physical temperature of the cloud, and $`J_{inc}`$ is the incident flux on the cloud from the upper atmosphere. In all cases it is likely that $`J_{cloud}`$ is greater than $`J_{inc}`$ so that the observed effect of the cloud is to increase the surface brightness along that line of sight.

If the cloud is close to the top of the water vapour column density then $`J_{inc}`$ is negligible compared to $`J_{cloud}`$ and one can simply write:

$$\delta J=exp(\tau _{cloud})\delta \tau J_{cloud}.$$
(10)
If one considers the signal at two wavelengths $`\lambda 1`$, $`\lambda 2`$ then the ratio is

$$\frac{\delta J_{\lambda 1}}{\delta J_{\lambda 2}}=\frac{exp(\tau _{\lambda 1})}{exp(\tau _{\lambda 2})}\frac{\delta \tau _{\lambda 1}}{\delta \tau _{\lambda 2}},$$
(11)
where we have removed the subscripts cloud. One can express this in terms of power by simply invoking the Rayleigh-Jeans law $`B(\lambda )J/(\lambda ^2)`$ so that,

$$\frac{\delta B_{\lambda 1}}{\delta B_{\lambda 2}}=\frac{exp(\tau _{\lambda 1})}{exp(\tau _{\lambda 2})}\frac{\delta \tau _{\lambda 1}}{\delta \tau _{\lambda 2}}\frac{\lambda 1^2}{\lambda 2^2},$$
(12)
and also to measured volts by applying an appropriate FCF

$$\frac{\delta V_{\lambda 1}}{\delta V_{\lambda 2}}=\frac{exp(\tau _{\lambda 1})}{exp(\tau _{\lambda 2})}\frac{\delta \tau _{\lambda 1}}{\delta \tau _{\lambda 2}}\frac{\lambda 1^2}{\lambda 2^2}\frac{FCF_{\lambda 2}}{FCF_{\lambda 1}}.$$
(13)
In the case where the opacities are simply proportional to water vapour column density (as we have shown in this paper) then $`\delta \tau _{\lambda 1}/\delta \tau _{\lambda 2}`$ is a constant, and we expect the ratios of the signal at the two wavelengths to be proportional to the ratios of transmission from the cloud to the telescope at the two wavelengths. One also expects

$$\frac{\delta J_{\lambda 1}}{exp(\tau _{\lambda 1})}\frac{\delta J_{\lambda 2}}{exp(\tau _{\lambda 2})}$$
(14)
where one has to be careful to note that the value of $`\tau `$ is to the cloud, not all the way to the top of the atmosphere.

For the SCUBA filters the ratios of $`\delta \tau _{850n}/\delta \tau _{450n}`$ and $`\delta \tau _{850w}/\delta \tau _{450w}`$ have been measured earlier in this paper. Jenness et al. (2001) have measured the FCFs from point sources using a 60 arcsecond radius aperture, which are also (by symmetry) valid for a uniform extended source of size 60 arcseconds. In the section on sky noise in this paper most of the photometry observations were taken with 60 arcsecond chop throws so that the skynoise measured is a measure of the difference in signal between two largely (but not exactly) overlapping regions of diameter 15 meters, 60 arcseconds apart. If the cloud is at a height, h above the telescope then the diameter of the beam in steradians is 15/h(m) or 25 arcminutes at a height of 2 km. This explains why the sky noise is uniform over the array.

In Fourier space we can consider the beam as a high spatial frequency filter of the form

$$exp\left(\frac{k^2}{15}\right)$$
(15)
i.e. one that filters out spatial scales $`<`$ 15 m and the chop as a low spatial frequency filter of the form

$$sin(2\pi kch)$$
(16)
where $`k`$ is the wavenumber, $`c`$ is the chop throw in steradians, and $`h`$ is the height of the cloud in metres. This can be simplified to $`2\pi kch`$ for the frequencies that are passed by the beam. The effect of having two differing spatial frequency filters is to couple only weakly (at a level $`ch/15`$) to clouds of approximate spatial scale 15 m. Clouds smaller than this are filtered out by the beam, clouds larger than this are filtered out by the chop throw. For typical values ($`h=2`$ km, $`c=60`$ arcseconds), the coupling is only 4 percent - illustrating how effective chopping is as a first stage to removing sky noise.

Although strictly speaking the FCFs calculated by Jenness et al. (2001) are only valid for extended far field objects, one can use the ratio of the values measured as an estimate of the ratio of the FCFs which would be valid for nearfield extended objects (i.e. the cloud in question). The validity of this assumption can then be assessed by comparison with the data.

If one therefore uses the various values discussed in the previous paragraph then one finds:

$$\frac{\delta V_{450W}}{\delta V_{850W}}=6.64\frac{exp(\tau _{450W})}{exp(\tau _{850W})},$$
(17)
$$\frac{\delta V_{450N}}{\delta V_{850N}}=3.32\frac{exp(\tau _{450N})}{exp(\tau _{850N})}$$
(18)
and

$$\frac{\delta B_{450W}}{\delta B_{850W}}=23.2\frac{exp(\tau _{450W})}{exp(\tau _{850W})},$$
(19)
$$\frac{\delta B_{450N}}{\delta B_{850N}}=21.2\frac{exp(\tau _{450N})}{exp(\tau _{850N})}$$
(20)
Given that we have shown we are coupling to sources larger in scale than 60 arcseconds, we could further correct these values by the ratio

$$\frac{\eta _{450(cloud)}}{\eta _{850(cloud)}}\frac{\eta _{850(60\mathrm{"})}}{\eta _{450(60\mathrm{"})}}$$
(21)
where $`\eta `$ is the coupling. However, given uncertainties in the cloud’s angular size (principally due to uncertainties in its height), and the fact that the ratio of the coupling values is likely to vary less significantly than the individual values, we choose not to.

## Appendix C Refraction

### C.1 Refraction Code at JCMT

The primary purpose of the current JCMT refraction model is to provide coarse corrections to pointing based on local atmospheric conditions. It assumes a well-behaved, non-turbulent atmosphere. To first order, refraction displaces images towards the zenith by an amount

$$R=Atan(z)$$
(22)
where $`z`$ is the zenith distance, and $`A`$ is a function of the local atmospheric parameters :

1. $`T`$ \- Temperature (above some mean, chosen as 4 K for Mauna Kea)
2. $`p`$ \- % Pressure change from the MK standard of 624 mb
3. $`h`$ \- % Humidity

Grids of atmospheric models were generated and refraction calculated at several zenith distances by integrating through the atmosphere. Simple functional forms of the results were sought that would enable calculation of refraction from local atmospheric conditions. The optical- and millimetre- $`A`$-terms in the formula above take principally these forms (units are arcseconds) :

$$A=35.8930.00067(h20)0.135(T4)+0.371p$$
(23)
at 0.55 µm and:

$$A=36.800+0.0768(h20)0.0294(T4)+0.371p$$
(24)
at 1 mm.

\- i.e. optical refraction is dependent mostly upon temperature, while submillimetre refraction is dependent mostly upon humidity. The 1 mm formulation used at JCMT contains additional, less significant terms. Its values were compared to the integration results throughout the grid and were found to be accurate to better than 1 arcsec for zenith distances less than 80 degrees under all but the most extreme conditions.

### C.2 Expectation Values of Signal Levels in the Presence of Refraction

This appendix describes the use of the pointing rms values determined from the SAO phase monitor to give an indication of the expectation values of the signal levels, relative to the values for no refraction ’noise’. The assumptions are

1. the refraction noise can be described by a Gaussian process, which in two dimensions gives the Rayleigh distribution:

$$P_\theta =\frac{\theta }{\sigma ^2}exp\left(\frac{\theta ^2}{2\sigma ^2}\right)$$
(25)
where Pθ is the probability of the refraction having a value of $`\theta `$ arcsecs, and $`\sigma `$ is the rms refraction ’noise’.
2. the telescope beam can also be defined as a Gaussian:

$$S(\theta )=S_oexp\left[4ln(2)\frac{\theta ^2}{\theta _{fwhm}^2}\right]$$
(26)
where S is the signal level which would be measured for a movement of the telescope of $`\theta `$ arcsecs, and $`\theta _{fwhm}`$ is the FWHM size of the beam in arcsecs.

The expectation value of the signal level, relative to the no refraction noise case, is given by

$$<S>=\frac{_0^{S_o}SP_S𝑑S}{_0^{S_o}P_S𝑑S}$$
(27)
where PS is the probability distribution for the signal.
Probabilities transform directly so that

$$P_SdS=P_\theta (d\theta )$$
(28)
where use is made of the fact that the probability of S decreasing is a function of the probability of $`\theta `$ increasing.
For simplicity, we rewrite (26) as

$$S(\theta )=S_oexp\left[\frac{\theta ^2}{2\beta ^2}\right]$$
(29)
where $`\beta ^2`$ = $`\frac{\theta _{fwhm}^2}{8ln2}`$ .
Now, from equation (29),

$$\frac{dS}{d\theta }=\frac{\theta }{\beta ^2}S$$
(30)
using this and equations (25),(29),

$$P_S=\frac{\beta ^2}{\sigma ^2}\frac{1}{S}\left[\frac{S}{S_o}\right]^{\frac{\beta ^2}{\sigma ^2}}$$
(31)
It should be noted that, if $`\beta ^2`$ = $`\sigma ^2`$, the probability is uniform.
If we now evaluate (27) using (31) we find that

$$\frac{<S>}{S_o}=\frac{1}{1+\frac{\sigma ^2}{\beta ^2}}=\frac{1}{1+\frac{8ln2\sigma ^2}{\theta _{fwhm}^2}}$$
(32)
In the case of a perfect 15-metre telescope, one can write for $`\theta _{fwhm}`$ in arcsecs :

$$\theta _{fwhm}=[14.025+0.1856T_E]\lambda $$
(33)
where $`T_E`$ is the illumination edge taper in decibels and $`\lambda `$ is the wavelength in mm (Goldsmith 1987). The edge taper determines the level of ground radiation accepted by the feed, and at JCMT $`T_E7.5`$ dB . In the case of a telescope with surface imperfections one must use the measured value of the full width half-maximum beam size in (32).
If we assume a 14.2′′ beam at 850 µm and a 7′′ beam at 450 µm, then (32) becomes

$$\frac{1}{1+0.0275\sigma ^2}\mathrm{for}\mathrm{the}14.2\mathrm{"}\mathrm{beam}$$
(34)
and

$$\frac{1}{1+0.1132\sigma ^2}\mathrm{for}\mathrm{the}7\mathrm{"}\mathrm{beam}$$
(35)
where $`\sigma `$ is the rms refraction noise in arcseconds.