text:
# How to measure squeezing and entanglement of Gaussian states without homodyning

## Preliminaries.

Let us begin with introducing the necessary notation and definitions. Let $`r=(x_1,p_1,\mathrm{},x_k,p_k,\mathrm{},x_N,p_N)`$ be the vector of conjugate quadratures of $`N`$ modes which satisfy the canonical commutation relations $`[x_j,p_k]=i\delta _{jk}`$. The Gaussian state is fully described by the vector of mean values $`\xi _j=r_j`$ and the covariance matrix

$$\gamma _{jk}=\mathrm{\Delta }r_j\mathrm{\Delta }r_k+\mathrm{\Delta }r_k\mathrm{\Delta }r_j,$$
(1)
where $`\mathrm{\Delta }r_j=r_j\xi _j`$. The quantum state of the optical field can be fully characterized by a $`s`$-parametrized quasidistribution which provides phase-space representation of the quantum state. For our purposes, it is convenient to utilize the Husimi Q-function. The Q-function of an $`N`$-mode Gaussian state is the Gaussian distribution Perinabook

$$Q(r)=\frac{\pi ^N}{\sqrt{det(\gamma +I)}}\mathrm{exp}\left[(r\xi )^T(\gamma +I)^1(r\xi )\right],$$
(2)
where $`I`$ is the identity matrix. The squeezing properties do not depend on $`\xi `$ and are fully described by $`\gamma `$. The maximal observable squeezing, i.e. the minimal quadrature variance, is called the generalized squeezing variance $`\lambda `$ and can be determined as the minimal eigenvalue of the covariance matrix Simon94 ,

$$\lambda =\mathrm{min}[\mathrm{eig}(\gamma )].$$
(3)
The purity of the mixed state with density matrix $`\rho `$ is defined as $`𝒫=\mathrm{Tr}[\rho ^2]`$. For a Gaussian state with covariance matrix $`\gamma `$, one obtains Kim02

$$𝒫=[det(\gamma )]^{1/2}.$$
(4)
## Single-mode case.

For the sake of simplicity, we first illustrate the procedure on single-mode Gaussian states ($`N=1`$). Consider the setup depicted in Fig. 1(a). The input mode impinges on a beam splitter BS with tunable transmittance $`T`$, and the output mode is measured by a photodetector PD with efficiency $`\eta `$ that is sensitive to single photons (no single-photon resolution is needed). We assume that this realistic detector can be modeled as a beam splitter with transmittance $`\eta `$ followed by an ideal detector that performs a dichotomic measurement described by the POVM elements $`\mathrm{\Pi }_0=|00|`$ and $`\mathrm{\Pi }_1=𝟙\mathbb{\Pi }_\mathrm{𝟘}`$. (In what follows, we assume that the detector is ideal and $`\eta `$ can be taken into account by substituting $`T\eta T`$.) The probability of no-click of an ideal detector PD is given by $`P=\mathrm{Tr}[\rho \mathrm{\Pi }_0]=0|\rho |0=2\pi Q(0)`$, so that inserting $`r=0`$ into Eq. (2) yields

$$P=\frac{2}{\sqrt{det(\gamma ^{}+I)}}\mathrm{exp}\left[\xi ^T(\gamma ^{}+I)^1\xi ^{}\right].$$
(5)
where $`\gamma ^{}`$ and $`\xi ^{}`$ are, respectively, the covariance matrix and the displacement vector of the beam impinging on the photodetector.

Suppose that we set the beam splitter transmittance to the value $`T_j`$. The covariance matrix $`\gamma ^{}`$ of the state after passing the beam splitter reads $`\gamma ^{}=T_j\gamma +(1T_j)I`$. Similarly, the coherent signal is damped to $`\xi ^{}=\sqrt{T_j}\xi `$. On inserting $`\gamma ^{}`$ and $`\xi ^{}`$ into Eq. (5), we obtain

$$P_j=\frac{2}{[det(\stackrel{~}{\gamma }_{T_j})]^{1/2}}\mathrm{exp}\left(T_j\xi ^T\stackrel{~}{\gamma }_{T_j}^1\xi \right).$$
(6)
where $`\stackrel{~}{\gamma }_{T_j}=T_j\gamma +(2T_j)I`$. We thus find that $`P_j`$ depends on $`T_j`$ and four parameters of the state: $`det\gamma `$, $`\mathrm{Tr}(\gamma )`$, $`g_1=\xi _1^2\gamma _{22}+\xi _2^2\gamma _{11}2\gamma _{12}\xi _1\xi _2`$, and $`g_2=\xi _1^2+\xi _2^2`$. This immediately suggests that if we measure $`P_j`$ for at least four different $`T_j`$’s, we might be able to reconstruct the values of these four parameters by solving a system of nonlinear equations. However, numerical simulations reveal that the inversion of these highly nonlinear equations typically leads to extremely large fluctuations of the estimated parameters even for a very large number of measurements for each setting $`T_j`$.

Fortunately, in the important case where the displacement vector is zero ($`\xi =0`$), the scheme provides reliable and well-behaved estimates of these parameters since formula (6) then simplifies to

$$det\stackrel{~}{\gamma }_{T_j}=T^2det(\gamma )+T(2T)\mathrm{Tr}(\gamma )+(2T)^2=4P_j^2.$$
(7)
This results in a system of *linear* equations for $`det(\gamma )`$ and $`\mathrm{Tr}(\gamma )`$. If measurements for two different transmittances $`T_1`$ and $`T_2`$ are performed and the observed probabilities of no-click are $`P_1`$ and $`P_2`$, then the system of Eqs. (7) can easily be solved and yields

$$\mathrm{Tr}(\gamma )=\frac{2}{T_2T_1}\left(\frac{T_2}{T_1P_1^2}\frac{T_1}{T_2P_2^2}\right)+2\frac{2}{T_1}\frac{2}{T_2},$$
(8)
$$det(\gamma )=\frac{2}{T_1T_2}\left(\frac{2T_2}{T_1P_1^2}\frac{2T_1}{T_2P_2^2}\right)+\frac{(2T_1)(2T_2)}{T_1T_2}.$$
(9)
Let us investigate what can be extracted from the knowledge of $`\mathrm{Tr}(\gamma )`$ and $`det(\gamma )`$. As noted above, the squeezing properties of the Gaussian state, namely the generalized squeezing variance $`\lambda `$, can be determined from the eigenvalues of $`\gamma `$, cf. Eq. (3). For a single-mode state, $`\gamma `$ is symmetric $`2\times 2`$ matrix and its eigenvalues can be expressed in terms of $`det(\gamma )`$ and $`\mathrm{Tr}(\gamma )`$, which are both determined by the present method. We find that

$$\lambda =\frac{1}{2}\left[\mathrm{Tr}(\gamma )\sqrt{\mathrm{Tr}^2(\gamma )4det(\gamma )}\right].$$
(10)
Moreover, since our method provides an estimate of $`det\gamma `$, we can also determine the purity from Eq. (4).

If $`\xi 0`$, our scheme is still usable provided that we can perform a collective measurement on two copies of the state, as depicted in Fig. 1(b). The two input modes $`A_I`$ and $`A_{II}`$ prepared in identical Gaussian state interfere on a balanced beam splitter. In the Heisenberg picture, the annihilation operators of the output modes $`A_+`$ and $`A_{}`$ are linear combinations of those of the input modes, $`a_\pm =2^{1/2}(a_I\pm a_{II})`$. The covariance matrix of mode $`A_{}`$ is equal to $`\gamma `$ Browne03  but the coherent signal in $`A_{}`$ vanishes due to the destructive interference, $`\xi _{}=0`$. As shown in Fig. 1(b), the mode $`A_{}`$ is subsequently sent to a direct measurement setup identical to that shown in Fig. 1(a).

## Multimode case.

We now extend this procedure to multimode Gaussian states. A reliable operation again requires two copies $`A_I`$ and $`A_{II}`$. Essentially, we use in parallel $`N`$ setups such as shown in Fig. 1(b). We combine each pair of modes $`A_{I,k}`$ and $`A_{II,k}`$, with $`k=1,\mathrm{},N`$, on a balanced beam splitter. Each “minus” mode $`A_{,k}`$, with $`\xi =0`$, is then sent to a beam splitter of transmittance $`T_j`$ followed by a photodetector. We measure the probability $`P_j=2^N/(det\stackrel{~}{\gamma }_{T_j})^{1/2}`$ that none of the $`N`$ detectors clicks. For an $`N`$-mode state, the determinant of $`\stackrel{~}{\gamma }_T=T\gamma +(2T)I`$ can be expanded as

$$det(\stackrel{~}{\gamma }_T)=\underset{n=1}{\overset{2N}{}}T^n(2T)^{2Nn}f_n(\gamma )+(2T)^{2N},$$
(11)
where $`f_n(\gamma )`$ is an homogeneous polynomial of $`n`$th order in the matrix elements of $`\gamma `$, e.g., $`f_{2N}(\gamma )=det(\gamma )`$ and $`f_1(\gamma )=\mathrm{Tr}(\gamma )`$. The probability $`P_j`$ thus depends on $`T_j`$ and the $`2N`$ parameters $`f_n(\gamma )`$. If we measure $`P_j`$ for $`2N`$ (or more) different transmittances $`T_j`$’s, then we can determine the parameters $`f_n`$ of the Gaussian state by solving a system of linear equations

$$\underset{n=1}{\overset{2N}{}}T_j^n(2T_j)^{2Nn}f_n(\gamma )=2^{2N}P_j^2(2T_j)^{2N}.$$
(12)
Once we know $`f_n(\gamma )`$, we can determine the generalized squeezing variance $`\lambda `$ as the smallest root of the characteristic polynomial $`det(\lambda I\gamma )=0`$. It can be seen from Eq. (11) that the parameters $`f_j(\gamma )`$ are the coefficients of this characteristic polynomial, and we have

$$\lambda ^{2N}+\underset{k=0}{\overset{2N1}{}}\lambda ^k(1)^kf_{2Nk}(\gamma )=0.$$
(13)
We can also determine the purity of the $`N`$-mode Gaussian state from $`f_{2N}(\gamma )`$ with the help of formula (4).

## Entanglement detection.

In the context of quantum information processing with continuous variables, the entanglement properties of Gaussian states deserve particular attention. It has been shown that a two-mode Gaussian state is separable iff it has a positive partial transpose Duan00 ; Simon00 . This property can easily be checked if one knows the covariance matrix

$$\gamma _{AB}=\left(\begin{array}{cc}\gamma _A& \sigma _{AB}\\ \sigma _{AB}^T& \gamma _B\end{array}\right)$$
(14)
of the bipartite state, where $`\gamma _A`$ and $`\gamma _B`$ are the covariance matrices of modes $`A`$ and $`B`$, respectively, while $`\sigma _{AB}`$ captures the intermodal correlations. Moreover, analytical formulas for several entanglement monotones that measure the entanglement of Gaussian states have been given in the literature Vidal02 ; Giedke03 . A particularly simple formula has been obtained for the logarithmic negativity $`E_𝒩`$ of an arbitrary Gaussian state. To calculate $`E_𝒩`$, we must determine the symplectic spectrum of the covariance matrix of the partially transposed state $`\rho _{AB}^{T_A}`$. As shown in Ref. Vidal02 , the symplectic eigenvalues are the two positive roots $`\zeta _1\zeta _2>0`$ of the biquadratic equation

$$\zeta ^4(det\gamma _A+det\gamma _B2det\sigma _{AB})\zeta ^2+det\gamma _{AB}=0.$$
(15)
The solution of Eq. (15) yields

$$\zeta _2^2=\frac{1}{2}\left(D\sqrt{D^24det\gamma _{AB}}\right),$$
(16)
where $`D=det\gamma _A+det\gamma _B2det\sigma _{AB}`$. The two-mode Gaussian state is entangled if and only if $`\zeta _2<1`$. In this case, we have $`E_𝒩=\mathrm{log}(\zeta _2)`$, while $`E_𝒩=0`$ otherwise. The condition $`\zeta _2<1`$ implies the necessary and sufficient entanglement condition $`D>1+det\gamma _{AB}`$ Simon00 ; Giedke01 , which explicitly reads,

$$det\gamma _A+det\gamma _B2det\sigma _{AB}>1+det\gamma _{AB}.$$
(17)
With the use of the method proposed in the present paper we can measure $`det\gamma _{AB}`$, $`det\gamma _A`$, and $`det\gamma _B`$. An upper bound on $`det\sigma _{AB}`$ in terms of these determinants can be derived from the condition that the symplectic eigenvalues $`\stackrel{~}{\zeta }_j`$ of the covariance matrix $`\gamma _{AB}`$ must be greater or equal to one Giedke01 . The lower eigenvalue $`\stackrel{~}{\zeta }_2`$ is given by Eq. (16), where $`D`$ is replaced with $`D^{}=det\gamma _A+det\gamma _B+2det\sigma _{AB}`$. The condition $`\stackrel{~}{\zeta }_2^21`$ yields

$$2det\sigma _{AB}det\gamma _{AB}+1det\gamma _Adet\gamma _B.$$
(18)
This, in turn, implies an upper bound on the lower symplectic eigenvalue $`\zeta _2`$ of the covariance matrix of $`\rho _{AB}^{T_A}`$. On inserting the upper bound on $`2det\sigma _{AB}`$ given by Eq. (18) into Eqs. (16) and (17) we find that $`\zeta _2<1`$, so that the state is entangled when

$$det\gamma _A+det\gamma _B>1+det\gamma _{AB}$$
(19)
holds. Inequality (19) is thus a sufficient condition for entanglement, but it is not necessary as some Gaussian entangled states are not detected by this test. The main advantage of this test is that all determinants appearing in Eq. (19) can be determined by *local* measurements supplemented with classical communication between $`A`$ and $`B`$. Moreover, if we can a priori assume that $`\xi =0`$, then measurements on a single copy of $`\rho _{AB}`$ suffice.

If we now want to exactly determine $`E_𝒩`$, we also need a scheme to measure $`det\sigma _{AB}`$. This can be accomplished provided that we can perform joint measurements on several copies of the state $`\rho _{AB}`$. Unlike the previous one, this scheme requires joint nonlocal measurements on modes $`A`$ and $`B`$, and it involves several steps as schematically illustrated in Fig. 2. Using the scheme of Fig. 2(b), we can measure the determinants of the covariance matrices $`\gamma _+`$ and $`\gamma _{}`$ of modes $`A_+`$ and $`A_{}`$ that are linear combinations of the modes $`A`$ and $`B`$, $`a_\pm =(a\pm b)/\sqrt{2}`$. After a simple algebra, we find that

$$det(\sigma _{AB}+\sigma _{AB}^T)=2det\gamma _++2det\gamma _{}det(\gamma _A+\gamma _B).$$
(20)
In order to determine $`det(\gamma _A+\gamma _B)`$, we have to carry out a joint measurement on two independent copies of the two-mode state, as depicted in Fig. 2(c). By mixing the modes $`A_1`$ and $`B_2`$ on a balanced beam splitter, we prepare an output single-mode state with covariance matrix $`\gamma _{A+B}(\gamma _A+\gamma _B)/2`$. Recall that the modes $`A_1`$ and $`B_2`$ belong to two independent copies of the two-mode state $`\rho _{AB}`$, hence $`A_1`$ and $`B_2`$ are uncorrelated. After the measurement of $`det\gamma _+`$, $`det\gamma _{}`$, and $`det\gamma _{A+B}`$, we calculate $`det(\sigma _{AB}+\sigma _{AB}^T)`$ from Eq. (20). It holds that $`det(\sigma _{AB}+\sigma _{AB}^T)4det\sigma _{AB}`$ and the equality is achieved when $`\sigma _{AB}`$ is symmetric. The matrix $`\sigma _{AB}`$ can be brought to a symmetric form by applying a local phase shift $`\mathrm{exp}(i\varphi b^{}b)`$ to the mode $`B`$ using the phase shifter PS in Fig. 2. This transforms $`\sigma _{AB}`$ to $`\sigma _{AB}U(\varphi )`$ where

$$U(\varphi )=\left(\begin{array}{cc}\mathrm{cos}\varphi & \mathrm{sin}\varphi \\ \mathrm{sin}\varphi & \mathrm{cos}\varphi \end{array}\right).$$

To determine the phase shift $`\varphi `$ that symmetrizes $`\sigma _{AB}`$, we measure $`y_\varphi det[\sigma _{AB}U(\varphi )+U^T(\varphi )\sigma _{AB}^T]`$ for three different phases $`\varphi =0`$, $`\pi /4`$, and $`\pi /2`$. It can be shown that

$$y_\varphi =y_0\mathrm{cos}^2\varphi +y_{\pi /2}\mathrm{sin}^2\varphi +\stackrel{~}{y}_{\pi /4}\mathrm{sin}(2\varphi ),$$

where $`\stackrel{~}{y}_{\pi /4}=y_{\pi /4}(y_0+y_{\pi /2})/2`$. The value of $`det\sigma _{AB}`$ can be found as the maximum of $`y_\varphi `$ over $`\varphi `$ which yields

$$det\sigma _{AB}=\frac{1}{8}\left[y_0+y_{\pi /2}+\sqrt{(y_0y_{\pi /2})^2+4\stackrel{~}{y}_{\pi /4}^2}\right].$$

This finally provides all the information required for the exact calculation of the logarithmic negativity $`E_𝒩`$.

In summary, we have proposed a scheme for the direct measurement of the squeezing, purity, and entanglement of Gaussian states that does not require homodyne detection but only needs beam splitters and photodetectors with single-photon sensitivity. The scheme generally requires joint measurements on two copies of the state, but single-copy measurements suffice if it is a priori known that the mean (or coherent) values of the quadratures vanish, which is, e.g., the case of the squeezed and entangled states generated by spontaneous parametric downconversion. We have shown that, based on Eq. (19), the present method can be used to assess the entanglement of $`1\times 1`$ Gaussian states by means of local measurements, without employing any local oscillator or interferometric schemes. Given the simplicity of the suggested setup, the prospects for an experimental realization in a near future look very good.

Note added: The sufficient condition on entanglement \[Eq. (19)\] has recently and independently been derived by Adesso *et al.* Adesso03 . Besides the lower bound on $`E_𝒩`$ linked to Eq. (19), an upper bound on $`E_𝒩`$ expressed in terms of the determinants of $`\gamma _A`$, $`\gamma _B`$, and $`\gamma _{AB}`$ was also derived in Ref. Adesso03 . Remarkably, these two bounds are typically very close to each other, so the knowledge of the determinants of the covariance matrices provides quite precise quantitative information on the entanglement, making the direct measurement procedure particularly powerful.

We acknowledge financial support from the Communauté Française de Belgique under grant ARC 00/05-251, from the IUAP programme of the Belgian government under grant V-18, from the EU under projects RESQ (IST-2001-37559) and CHIC (IST-2001-32150). JF also acknowledges support from the grant LN00A015 of the Czech Ministry of Education.


format: text

text:
# A Simple Signal of Noncommutative Space

(November 12, 2003)

## Abstract

We propose a simple low-energy classical experiment in which the effects of noncommutativity can be clearly separated from commutative physics. The ensuing bound on the noncommutative scale is remarkable, especially in view of its elementary derivation.

Noncommutative mechanics has recently been under intensive investigation . It uses the following basic commutators for the coordinates $`q_i`$ and momenta $`p_i`$,

$$[q_i,q_j]=i\theta _{ij},[q_i,p_j]=i\delta _{ij},[p_i,p_j]=iF_{ij},$$
(1)
which are a generalization of the Heisenberg commutation relations. The time evolution of an operator follows from its commutator with the prescribed Hamiltonian. In the classical case the commutators in (1) are replaced by Poisson brackets,

$$\{q_i,q_j\}=\theta _{ij},\{q_i,p_j\}=\delta _{ij},\{p_i,p_j\}=F_{ij}.$$
(2)
Above, $`\delta _{ij}`$ is a unit matrix, whereas $`F_{ij}`$ and $`\theta _{ij}`$ are functions of both the coordinates and the momenta, to be restricted only by the Jacobi identities. This note will consider the case in which $`\theta `$ is a nonvanishing constant diagonal matrix, whereas a generic $`F(q_i,p_j)`$ is allowed. A dynamical effect due to the gradient of $`F`$ will be reported. For simplicity in presentation we will work in two dimensions, $`i,j=1,2`$ in (2), although our considerations extend straightforwardly to any space dimensionality. Lower case latin indices will always take two values, e.g. $`i,j,k,s,t=1,2`$, and no distinction will be necessary between up and down indices in the following. In two dimensions the only nontrivial Jacobi identities read, with the notation $`\theta _{12}=\theta _{21}\theta `$,

$$_{p_1}F=\theta _{q_2}F,_{p_2}F=\theta _{q_1}F.$$
(3)
The above relations ensure the invariance of Eqs. (2) under time evolution. Once we disregard the constant $`F`$ solution 111Experimentally $`F`$ is never constant throughout all space, moreover the effect we are after arises when $`F0`$., (3) implies that $`F`$ is a function of the combinations $`\overline{q}_1=q_1+\theta p_2`$ and $`\overline{q}_2=q_2\theta p_1`$, in brief $`F=F(\overline{q}_1,\overline{q}_2)=F(\overline{q}_i)=F(q_i+\theta _{ij}p_j)`$. If $`\theta `$ is then assumed to be small, which is the maximum the real world allows for, one has

$$F(\overline{q}_i)=F(q_i+\theta _{ij}p_j)=F(q_i)+\theta _{st}p_t_{\overline{q}_s}F_{ij}(\overline{q})|_{\overline{q}_i=q_i}.$$
(4)
To see the consequences of such an expression, consider the equations of motion for a free Hamiltonian $`H=\frac{p_i^2}{2m}`$,

$$m\dot{q_i}=m\{q_i,H\}=p_i,m\dot{p_i}=m\{p_i,H\}=F_{ij}(q_s+\theta _{st}p_t)p_j,$$
(5)
where $`\frac{df}{dt}\dot{f}`$. If $`\theta =0`$, $`F_{ij}`$ just provides a magnetic field of intensity $`F_{ij}(q)`$. If $`\theta `$ is small, (4) implies that

$$m\ddot{q}_i=F(q^s+\theta _{st}p_t)_{ij}\dot{q}_jF(q^s)_{ij}\dot{q}_j+m_sF_{ij}(q)\theta _{st}\dot{q}_t\dot{q}_j.$$
(6)
Eq. (6) displays the superposition of a usual electromagnetic background, linear in the velocity, and of a force proportional to the mass and quadratic in the velocities, $`_sF(q)_{mj}\theta _{st}\dot{q}_t\dot{q}_j`$ – echoing gravity. However $`\gamma _{tj}^s_sF(q)_{mj}\theta _{st}`$ simulates a gravitational connection only in a particular reference frame, since $`\gamma _{tj}^s`$ does not behave like a Christoffel symbol under generalized coordinate transformations.

Experimentally speaking, it is important that a nonzero $`\theta `$ requires a nonconstant $`F`$ to be a function of both position and velocity, $`F(q,p)`$. This suggests that noncommutativity could be seen through the effects of the additional term in (6). A simple experiment is actually possible and provides - given its elementary character - a remarkably good bound for $`\theta `$.

Consider the 2-dimensional case of motion on a plane labelled by coordinates $`q_1=x`$, $`q_2=y`$, with velocities $`v_x=\dot{x}=p_1/m`$, $`v_y=\dot{y}=p_2/m`$, and assume a magnetic field of spatial dependence $`F(x,y)=Be(y)`$, with $`e(y)`$ the Heaviside step function, $`e(y)=1`$ if $`y0`$, $`e(y)=0`$ if $`y<0`$, $`\frac{de(y)}{dy}=\delta (y)`$. B is constant. The set-up is sketched in figure 1a.

At the point P a particle with initial velocity $`\stackrel{}{v}_0`$ passes from the lower half-plane (with vanishing magnetic field) into the upper half-plane, which has magnetic field $`B`$. At that same point the particle’s velocity changes to $`\stackrel{}{v}^{}`$ if $`\theta 0`$, as will be seen; subsequently the usual commutative circular motion in a magnetic field $`B`$ follows. While exiting the superior half-plane at point P’, the velocity will change to a final $`\stackrel{}{v}`$. Explicitely, Eq. (6) leads to

$$m\ddot{x}F(y)\dot{y}=m\theta B\delta (y)\dot{x}\dot{y},m\ddot{y}+F(y)\dot{x}=m\theta B\delta (y)\dot{x}^2.$$
(7)
Using $`\delta (y)=\delta (v_yt)=\frac{\delta (t)}{|v_y|}`$, one has

$$m\ddot{x}F(y)v_y=m\theta B\delta (t)\dot{x}\frac{v_y}{|v_y|},m\ddot{y}+F(y)v_x=m\theta B\delta (t)\frac{v_x^2}{|v_y|}.$$
(8)
Integrating from $`t=ϵ`$ to $`t=ϵ`$ (and taking $`ϵ0`$), one obtains

$$v_x^{}=v_x(1\theta B),v_y^{}=v_y(1+\theta B\frac{v_x^2}{v_y^2}).$$
(9)
As long as the velocity has a component orthogonal to the magnetic field gradient, i.e. if $`v_x0`$, an instantaneous velocity change of order $`\theta `$ takes place at the point at which the particle enters the magnetic field. No instantaneous displacement appears. The change in kinetic energy is only of order $`\theta ^2`$, $`(v_x^{})^2+(v_y^{})^2=v_x^2+v_y^2+\theta ^2B^2(v_x^2+v_x^4/v_y^2)`$. The opposite mechanism works while the particle exits the magnetic field region. This would render the correction quadratic in $`\theta `$, hence quite small. To avoid that, figure 1b proposes a slightly modified configuration, which allows the particle to exit the magnetic field along the gradient. In this case no unconventional velocity shift occurs, and the particle keeps the $`O(\theta )`$ correction (9), with $`v_x=v_x^{\prime \prime }=v_x^{}`$, $`v_y=v_y^{\prime \prime }=v_y^{}`$.

It is important that the noncommutative correction (9) depends also on the ratio $`\frac{v_x^2}{v_y^2}`$; this allows us to choose a favorable situation. Assuming that experimental control can be maintained for velocity ratios $`\frac{v_x}{v_y}10^5`$, a field $`B1T`$, and considering a proton of mass $`m_p`$ and charge $`1.6\times 10^{19}C`$, one sees total reflection along the line PP’ if at least

$$\sqrt{\theta }2\times 10^{13}meters.$$
(10)
Although this preliminary bound may not seem a very strong one, compared to the $`10^{18}m`$ presently reached in particle accelerators, it is remarkable that a simple, classical, low-energy - hence cheap - set-up could achieve that. Moreover, the bound can be significantly improved by constructing a system of magnets designed to keep the proton under periodic motion. Individually, each magnet will act on the charged particle as in figure 1b. The velocity shifts will always add, rendering the bound cummulative. Then, $`10^{2N}`$ revolutions will further decrease the bound on $`\sqrt{\theta }`$ by a factor of $`10^N`$. One may also measure corrections of about one per cent in the velocity shift (9), and not of order one as total reflection required. The above considerations lead to a theoretical bound

$$\sqrt{\theta }2\times 10^{14N}meters,$$
(11)
where $`10^{2N}`$ is the number of revolutions performed by the particle inside the system of magnets. According to experience with the much higher energies used at particle accelerators, $`N=5`$ appears to be easily achievable, pushing the noncommutativity scale in the $`TeV`$ regime: $`\mathrm{\Lambda }_{NC}\theta ^{1/2}1TeV`$! This gives a bound similar to the much more complex best estimates performed either via quantum mechanical  of quantum field theoretic  techniques. Stronger claims  are manifestly based on involved considerations of radiative corrections. The bound (11) will of course be affected by experimental errors or difficulties. A realistic set-up will most probably require numerical work in order to follow the trajectories in a real magnetic field, which cannot be as simple as in figures 1a,b.

Although the above low-energy, cummulative, clasical experiment might seem attractive, it is likely that the interference effects characteristic of quantum mechanics, associated with the choice of a convenient set-up, will further improve the bound. This is currently under investigation. The main purpose of this note has been to demonstrate a clear effect of classical noncommutativity – anomalous deflection in magnetic field gradients – in as simple a system as possible.

Acknowledgements

It is a pleasure to thank K.N. Anagnostopoulos, G.G. Athanasiu and T.N. Tomaras for useful discussions. This project was partly supported by a European Community individual Marie Curie fellowship, under contract HPMF-CT-2000-1060. The subsequent hospitality of Prof. Tomaras at Crete University is gratefully acknowledged.


text:
# Ultrahigh energy neutrinos and non-linear QCD dynamics

## I Introduction

The inelastic interaction of neutrinos with nucleons is in general described by the QCD-improved parton model. In ultrahigh energy neutrino interactions one probes a kinematical region which is not accessible in the current collider experiments. It is known so far that these reactions can be sensitive to the domain for the usual Bjorken variable $`xm_{W,Z}^2/E_\nu 10^8`$ at $`E_\nu 10^{12}`$ GeV and electroweak boson virtualities $`Q^2m_{W,Z}^210^3`$ GeV2, where $`m_{W,Z}`$ are the boson masses. In contrast, in the current electron-proton collider DESY HERA the interactions are probed in the typical kinematical range $`x>10^5`$ and $`0Q^2<10^3`$ GeV2 (the variables $`x`$ and photon virtuality $`Q^2`$ are correlated, with small $`x`$ data associated to low virtualities). There, the measurements of the deep inelastic $`ep`$ scattering structure functions have put important constraints on parton distributions, mainly sea-quarks and gluons, in the small $`x`$ region. In perturbative QCD these distributions are expected to grow as $`x`$ decreases, which has been strongly confirmed by HERA data. Therefore, in order to provide accurate predictions for the ultrahigh energy neutrino-nucleon cross sections, a precise extrapolation of the structure functions to the small $`x`$ and large $`Q^2`$ region probed in these interactions is needed.

In the usual QCD-improved parton model, the power growth of the parton distributions functions in very high neutrino energies lead consequently to a power increasing of the neutrino-nucleon total cross section, which in turn implies in a violation of unitarity at high energies. This issue is an outstanding theoretical challenge, which has produced intense work towards a complete understanding of the non-linear (higher twist) QCD corrections in the high parton density regime at fixed impact parameter (for a recent review, see armestoreview  and references therein). At present, one has reliable estimates for the kinematical domain where non-linear (saturation) effects should play an important role. Namely, below a typical transverse momentum scale $`Q_{\mathrm{sat}}^2(x)x^\lambda `$, the so-called saturation scale QSATDEF , the growth of parton (gluon) density is tamed towards the black-disk limit of the target. This slow down of the gluon distribution function prevents the total cross section violating unitarity. An important prediction of the non-linear QCD approaches is the property of geometric scaling SGK  below saturation scale, where the gluon density (and cross sections) scales with the saturation momentum, which grows with $`1/x`$ as a power. It has also been shown that such scaling survives even above the saturation regime IANCUGEO ; STASTOGEO ; MUNIERWALLON : the kinematical region has been estimated to be $`Q_{\mathrm{sat}}^2(x)<Q^2<Q_{\mathrm{sat}}^4(x)/\mathrm{\Lambda }_{\mathrm{QCD}}^2`$ relying on the BFKL equation BFKL  in the presence of saturation IANCUGEO . This feature has important consequences in the underlying QCD dynamics at the nucleus-nucleus accelerators RHIC and LHC and should pose strong constraints to the cross sections beyond current accelerator energies. Concerning high energy neutrino cross sections, there are estimates constrained by HERA data based on DGLAP DGLAP  and effects of BFKL BFKL  linear evolutions DKRS ; GQRS ; RSSSV ; GRVNEUT ; BASU ; KMSNEUT  as well as implementations of non-linear QCD corrections KUTAK . The linear NLO DGLAP and unified BFKL-DGLAP results are found to be consistent with each other presenting somewhat small deviations KMSNEUT . An interesting analytical calculation using an approximate DGLAP solution with initial conditions from a soft non-perturbative model has been recently reported in Ref. FIORE . Saturation effects are shown to be quite small in recent DGLAP analysis BASU , whereas it was found a reduction in the cross section by a factor 2 at $`E_\nu 10^{12}`$ GeV considering the saturation model or screening effects via Balitsky-Kovchegov equation KUTAK . Recently, it has been claimed in Ref. JAMAL  that the geometric scaling property can lead to an enhancement of the neutrino-nucleon total cross section by an order of magnitude compared to the leading twist cross section and would lead it further to its unitarization.

In this work one presents estimates for the ultrahigh energy neutrino-nucleon cross sections, focusing on the saturation models which encode the physics of non-linear QCD evolution in a phenomenological way. The main point is those models provide a suitable dynamical extrapolation to very high energies/virtualities in very contrast with usual DGLAP based approaches, which are not a priori valid in that regime and are dependent on quark and gluon distributions also unknown (no constraint from accelerator data) in that kinematical limit. Based on the color dipole framework, first one presents the result for the saturation model GBW  supplemented by QCD evolution BGBK . The original saturation model GBW  was proposed for either low values of $`Q^2<150`$ GeV2 , whereas high energy neutrino interactions are sensitive to large virtualities $`Q^2m_{W,Z}^2`$. Therefore, it is necessary to perform calculations with complete QCD evolution. Furthermore, one computes estimates of the ultrahigh neutrino-nucleon cross section through an implementation of the dipole cross section considering BFKL formalism in the geometric scaling region IIM .In the last section, the numerical results for both cases are shown and contrasted with recent NLO DGLAP estimates and with unified BFKL-DGLAP formalism without saturation effects.

## II Neutrino-nucleon cross section

Deep inelastic neutrino scattering can proceed via $`W^\pm `$ or $`Z^0`$ exchanges. The first case corresponds to charged current (CC) and the second one to the neutral current (NC) interactions, respectively. The standard kinematical variables which describe the processes above are given by, $`s=2m_NE_\nu `$, $`Q^2=q^2`$, $`x=\frac{Q^2}{2pq}`$ and $`y=\frac{pq}{m_NE_\nu }`$. Here, $`m_N`$ is the nucleon mass, $`E_\nu `$ labels the neutrino energy and $`p`$ and $`q`$ are the four momenta of the nucleon and of the exchanged boson, respectively. The usual Bjorken variable and the momentum transfer are usually labeled as $`x`$ and $`Q^2`$, whereas $`s`$ is the total centre-of-mass energy squared. One assumes an isoscalar target $`N=(p+n)/2`$, which is a good approximation for the present purpose.

At high energies a quite successful framework to describe QCD interactions is provided by the color dipole formalism DIPOLEPIC , which allow an all twist computation (in contrast with the usual leading twist approximation) of the structure functions. The physical picture of the interaction is the deep inelastic scattering at low $`x`$ viewed as the result of the interaction of a color $`q\overline{q}`$ dipole which the gauge bosons fluctuate to with the nucleon target. The interaction is modeled via the dipole-target cross section, whereas the boson fluctuation in a color dipole is given by the corresponding wave function. The DIS structure functions for neutrino-nucleon scattering in the dipole picture read as KUTAK ; JAMAL ,

$`F_{T,L}^{\mathrm{CC},\mathrm{NC}}(x,Q^2)={\displaystyle \frac{Q^2}{4\pi ^2}}{\displaystyle d^2𝒓_0^1𝑑z|\psi _{T,L}^{W^\pm ,Z^0}|^2\sigma _{dip}},`$ (1)
where $`𝒓`$ denotes the transverse size of the color dipole, $`z`$ the longitudinal momentum fraction carried by a quark and $`\psi _{T,L}^{W,Z}`$ are proportional to the wave functions of the (virtual) charged or neutral gauge bosons corresponding to their transverse or longitudinal polarization ($`F_2=F_T+F_L`$). Explicit expressions for the boson ($`W^\pm `$ and $`Z^0`$) wave functions squared are as follows KUTAK ,

$`|\psi _T^{W^\pm }(r,z,Q^2)|^2`$ $`=`$ $`{\displaystyle \frac{6}{\pi ^2}}[z^2+(1z)^2]\epsilon ^2K_1^2(\epsilon 𝒓),`$ (2)
$`|\psi _L^{W^\pm }(r,z,Q^2)|^2`$ $`=`$ $`{\displaystyle \frac{24}{\pi ^2}}z^2(1z)^2Q^2K_0^2(\epsilon 𝒓),`$ (3)
$`|\psi _T^{Z^0}(r,z,Q^2)|^2`$ $`=`$ $`{\displaystyle \frac{3}{2\pi ^2}}K_W[z^2+(1z)^2]\epsilon ^2K_1^2(\epsilon 𝒓),`$ (4)
$`|\psi _L^{Z^0}(r,z,Q^2)|^2`$ $`=`$ $`{\displaystyle \frac{6}{\pi ^2}}K_Wz^2(1z)^2Q^2K_0^2(\epsilon 𝒓),`$ (5)
where one defines the auxiliary variable $`\epsilon =z(1z)Q^2`$ and $`K_{0,1}(x)`$ are the Mc Donald’s functions. Here, $`K_W=(L_u^2+L_d^2+R_u^2+R_d^2)`$ and the the chiral couplings are expressed as functions of the Weinberg angle $`\theta _W`$ as follows,

$`L_u`$ $`=`$ $`1{\displaystyle \frac{4}{3}}\mathrm{sin}^2\theta _W,L_d=1+{\displaystyle \frac{2}{3}}\mathrm{sin}^2\theta _W,`$ (6)
$`R_u`$ $`=`$ $`{\displaystyle \frac{4}{3}}\mathrm{sin}^2\theta _W,R_d={\displaystyle \frac{2}{3}}\mathrm{sin}^2\theta _W.`$ (7)
Following Ref. KUTAK , here one considers only four flavors ($`u,d,s,c`$) and assumes them massless, whereas it has been shown KMSNEUT  that heavy quarks $`(b,t)`$ give relatively small contribution. Moreover, the color dipoles contributing to Cabibbo favored transitions are $`u\overline{d}(d\overline{u})`$, $`c\overline{s}(s\overline{c})`$ for CC interactions and $`u\overline{u},d\overline{d},c\overline{c},s\overline{s}`$ for NC interactions. In our further numerical calculations using the dipole framework, we use the structure function $`F_3^{\mathrm{CC},\mathrm{NC}}`$ given by the usual LO DGLAP expressions NUDGLAP , despite their contributions to be small for the present purpose.

The total CC (NC) neutrino-nucleon cross sections as a function of the neutrino energy are given by the integration over available phase space at the given neutrino energy. They read as,

$`\sigma _{(\nu ,\overline{\nu })}^{\mathrm{CC},\mathrm{NC}}(E_\nu )={\displaystyle _{Q_{\mathrm{min}}^2}^s}𝑑Q^2{\displaystyle _{Q^2/s}^1}𝑑x{\displaystyle \frac{1}{xs}}{\displaystyle \frac{^2\sigma _{(\nu ,\overline{\nu })}^{\mathrm{CC},\mathrm{NC}}}{xy}}`$ (8)
where $`y=Q^2/(xs)`$ and a minimum value $`Q_{\mathrm{min}}^2`$ (of order of a few GeV’s) on $`Q^2`$ is introduced in order to stay in the deep inelastic region.

In the present work we are interested in the ultrahigh energy neutrinos ($`E_\nu 10^7`$ GeV), where the valence quark contribution stays constant and physics is driven by sea quark contributions. In this region, we will use the available phenomenology on non-linear QCD evolutions to study saturation effects in the neutrino-nucleon cross section.

The dipole cross section $`\sigma _{dip}`$, describing the dipole-nucleon interaction, is substantially affected by saturation effects at dipole sizes $`𝒓>1/Q_{\mathrm{sat}}`$. Here, we follow the saturation model GBW , which interpolates between the small and large dipole configurations, providing color transparency behavior, $`\sigma _{dip}𝒓^2`$, as $`𝒓0`$ and constant behavior, $`\sigma _{dip}\sigma _0`$, at large dipoles. The transition is rendered by the saturation phenomenon. The parametrisation for the dipole cross section takes the eikonal-like form,

$`\sigma _{dip}(\stackrel{~}{x},𝒓^2)=\sigma _0\left[\mathrm{\hspace{0.17em}1}\mathrm{exp}\left({\displaystyle \frac{Q_s^2(x)𝒓^2}{4}}\right)\right](1\stackrel{~}{x})^7,`$ (9)
where the saturation scale $`Q_s^2(x)=\left(\frac{x_0}{\stackrel{~}{x}}\right)^\lambda `$ defines the onset of the saturation effects. As known, most of the contribution to ultrahigh energy neutrino cross section comes from very small $`xm_{W,Z}^2/2m_NE_\nu `$ and therefore a sizeable part of the contributions is in the geometric scaling region or even in the saturation region JAMAL . The parameters ($`\sigma _0`$, $`\lambda `$ and $`x_0`$) were obtained from a fit to the HERA data GBW . The variable $`\stackrel{~}{x}=x(\mathrm{\hspace{0.17em}1}+4m_f^2/Q^2)`$ gives a suitable transition to the photoproduction region. The saturation model has been used to compute neutrino-nucleon cross section in Ref. KUTAK . There, it was shown that the saturation model does not include complete DGLAP QCD evolution at high virtualities and tends to be therefore lower than other calculations, probably meaning that it is rather incomplete and inaccurate in the high $`Q^2`$ region.

Recently, a new implementation of the model including QCD evolution BGBK  (labeled BGBK model) has appeared. Now, the dipole cross section depends on the gluon distribution as,

$`\sigma _{dip}=\sigma _0\left[\mathrm{\hspace{0.17em}1}\mathrm{exp}\left({\displaystyle \frac{\pi ^2𝒓^2\alpha _s(\mu ^2)\stackrel{~}{x}G(\stackrel{~}{x},\mu ^2)}{3\sigma _0}}\right)\right],`$ (10)
where the initial condition at $`\mu ^2=1`$ GeV2 is $`xG=A_gx^{\lambda _g}(1x)^{5.6}`$ and $`\mu ^2=C/𝒓^2+\mu _0^2`$. The phenomenological parameters are determined from a fit to small $`x`$ HERA data. The function $`G(x,\mu ^2)`$ is evolved with the leading order DGLAP evolution equation for the gluon density with initial scale $`Q_0^2=1`$ GeV2. The improvement preserves the main features of the low-$`Q^2`$ and transition regions, while providing QCD evolution in the large-$`Q^2`$ domain.

Despite the saturation model to be very successful in describing HERA data, its functional form is only an approximation motivated by the Glauber-Mueller formula, which it does not include impact parameter dependence. On the other hand, an analytical expression for the dipole cross section can be obtained within the BFKL formalism. Currently, intense theoretical studies has been performed towards an understanding of the BFKL approach in the border of the saturation region IANCUGEO ; MUNIERWALLON . In particular, the dipole cross section has been calculated in both LO BFKL  and NLO BFKL NLOBFKL  approach in the geometric scaling region BFKLSCAL . It reads as,

$`\sigma _{dip}=\sigma _0\left[𝒓^2Q_{\mathrm{sat}}^2(x)\right]^{\gamma _{\mathrm{sat}}}\mathrm{exp}\left[{\displaystyle \frac{\mathrm{ln}^2\left(𝒓^2Q_{\mathrm{sat}}^2\right)}{2\beta \overline{\alpha }_sY}}\right],`$ (11)
where $`\sigma _0=2\pi R_p^2`$ ($`R_p`$ is the proton radius) is the overall normalization and the power $`\gamma _{\mathrm{sat}}`$ is the (BFKL) saddle point in the vicinity of the saturation line $`Q^2=Q_{\mathrm{sat}}^2(x)`$ (the anomalous dimension is defined as $`\gamma =1\gamma _{\mathrm{sat}}`$). As usual in the BFKL formalism, $`\overline{\alpha }_s=N_c\alpha _s/\pi `$, $`\beta 28\zeta (3)`$ and the notation $`Y=\mathrm{ln}(1/x)`$. The quadratic diffusion factor in the exponential gives rise to the scaling violations. Eq. (11) has been used in Ref. JAMAL  to show an estimation of the enhancement factor (about a factor 10) in ultrahigh neutrino-nucleon cross section in contrast with the usual leading twist calculations.

The dipole cross section in Eq. (11) does not include an extrapolation from the geometric scaling region to the saturation region. This has been recently implemented in Ref. IIM , where the dipole amplitude $`𝒩(x,𝒓)=\sigma _{dip}/2\pi R_p^2`$ was constructed to smoothly interpole between the limiting behaviors analytically under control: the solution of the BFKL equation for small dipole sizes, $`𝒓1/Q_{\mathrm{sat}}(x)`$, and the Levin-Tuchin law LTLAW  for larger ones, $`𝒓1/Q_{\mathrm{sat}}(x)`$. A fit to the structure function $`F_2(x,Q^2)`$ was performed in the kinematical range of interest, showing that it is not very sensitive to the details of the interpolation. The dipole cross section was parametrized as follows,

$`\sigma _{dip}=\sigma _0\{\begin{array}{cc}𝒩_0\left(\frac{𝒓Q_{\mathrm{sat}}}{2}\right)^{2\left(\gamma _{\mathrm{sat}}+\frac{\mathrm{ln}(2/𝒓Q_{\mathrm{sat}})}{\kappa \lambda Y}\right)},\hfill & (𝒓Q_s2),\hfill \\ 1\mathrm{exp}^{a\mathrm{ln}^2(b𝒓Q_{\mathrm{sat}})},\hfill & (𝒓Q_s>2),\hfill \end{array}`$ (14)
where the expression for $`𝒓Q_{\mathrm{sat}}(x)>2`$ (saturation region) has the correct functional form, as obtained either by solving the Balitsky-Kovchegov (BK) equation BK , or from the theory of the Color Glass Condensate (CGC) CGCREVIEW . Hereafter, we label the model above by CGC. The coefficients $`a`$ and $`b`$ are determined from the continuity conditions of the dipole cross section at $`𝒓Q_{\mathrm{sat}}(x)=2`$. The coefficients $`\gamma _{\mathrm{sat}}=0.63`$ and $`\kappa =9.9`$ are fixed from their LO BFKL values. In our further calculations it will be used the parameters $`R_p=0.641`$ fm, $`\lambda =0.253`$, $`x_0=0.267\times 10^4`$ and $`𝒩_0=0.7`$, which give the best fit result. A large $`x`$ threshold factor $`(1x)^5`$ will be also considered.

## III Results and discussion

In what follows, one computes the ultrahigh energy CC and NC neutrino-nucleon cross sections using the different models containing saturation effects. In Fig. 1 one presents the results for the saturation model (GBW), the BGBK model and also one includes a DGLAP calculation. The result for the saturation model (dot-dashed lines) has been recently computed in Ref. KUTAK , which we reproduce here. The introduction of the DGLAP evolution in the saturation model enhances the cross section by a factor 3 or even more at $`E_\nu 10^{12}`$ GeV. This is shown by the BGBK curves (solid lines). It should be noticed that the power growth for saturation model is milder than the BGBK model and it gives an upper limit for the role played by saturation effects in the cross section. It should be emphasized this assertion concerns to a comparison between phenomenological saturation model and its version including DGLAP evolution. Namely, other non-linear QCD approaches can give a larger amount of saturation in the same kinenatical regime and the phenomenological saturation model probably is inaccurate in this region. For sake of completeness, the DGLAP results (LO and NLO) are also presented. For this purpose, we follow the calculations in Ref. BASU , where a detailed analysis of the NLO corrections to neutrino-nucleon cross sections was performed. The BGBK model has an energy dependence at ultrahigh energy neutrinos similar to the NLO DGLAP (long-dashed lines) for $`E_\nu >10^9`$ GeV.

In order to investigate the BFKL physics in the geometric scaling region, we use the recent fit to the dipole cross section presented in Ref. IIM  (the CGC model), which encodes the BFKL formalism in the geometric scaling regime and an interpolation to the saturation domain, as presented in the last section. It should be noticed that an analysis including DGLAP evolution in the fit would be timely. However, as the CGC model produces an “effective” anomalous dimension, i.e. $`\gamma _{\mathrm{eff}}(x,𝒓Q_{\mathrm{sat}})\gamma _{sat}+(\kappa \lambda Y)^1\mathrm{ln}(4/𝒓^2Q_{\mathrm{sat}}^2)`$ IIM  far from the geometric scaling region $`𝒓^2Q_{\mathrm{sat}}^2(x)1`$, which is closer to the DGLAP one, the estimate presented here would be not strongly spoiled by QCD evolution. The results for the CGC model (solid lines) is presented in Fig. 2, contrasted with the NLO DGLAP calculation (dashed lines) and also with the unified BFKL-DGLAP formalism (dot-dashed lines), which embodies non-leading $`\mathrm{ln}(1/x)`$ contributions KMSNEUT . The latter incorporates both the $`\mathrm{ln}(1/x)`$ BFKL resummation and the complete LO DGLAP evolution, with the dominant non-leading $`\mathrm{ln}(1/x)`$ contributions resummed to all orders. The result found for CGC model is quite similar to both DGLAP and unified BFKL-DGLAP formalisms in the range of energy considered here. It is not observed any enhancement in the neutrino-nucleon cross section in comparison with the leading twist DGLAP calculation, even at very high energies of order $`E_\nu 10^{12}`$ GeV. This fact is in disagreement with Ref. JAMAL , where it has been shown that the geometric scaling lead to an enhancement of the neutrino-nucleon total cross section by an order of magnitude versus the leading twist cross section. We also verify that at $`E_\nu >10^{11}`$ GeV the power on energy starts to slow down and it should be increasingly smoother (logarithmic growth) at higher energies. This expected change of behavior on energy of the cross section is estimated to be reached at $`E_\nu 10^{18}`$ GeV in Ref. JAMAL , suppressing the cross section compared to the leading twist result.

We can qualitatively understand the result above analyzing the ratio between the dipole cross section for BFKL and DGLAP approaches in the geometric scaling regime. As discussed before, the ultrahigh energy neutrino cross section is sensitive to the kinematical region of small $`xm_{W,Z}^2/E_\nu `$ and high $`Q^2m_{W,Z}^2`$. Concerning DGLAP formalism, in such a region we are in the double logarithmic approximation (DLA) limit. Concerning the non-linear QCD approaches, at the vicinity of saturation, the DLA saddle point is no longer equal to one, but it was found to be $`\gamma _{\mathrm{sat}}^{\mathrm{DLA}}1/2`$ IANCUGEO ; STASTOGEO ; GBMS . The ratio of the cross section in the geometric scaling region is given by,

$`R_{\mathrm{scal}}{\displaystyle \frac{\left[𝒓^2Q_{\mathrm{sat}}^2\right]^{\gamma _{\mathrm{sat}}^{\mathrm{BFKL}}}}{\left[𝒓^2Q_{\mathrm{sat}}^2\right]^{\gamma _{\mathrm{sat}}^{\mathrm{DLA}}}}}\left[{\displaystyle \frac{Q_{\mathrm{sat}}^2}{m_{W,Z}^2}}\right]^{\left(\gamma _{\mathrm{sat}}^{\mathrm{BFKL}}\gamma _{\mathrm{sat}}^{\mathrm{DLA}}\right)}𝒪(1),`$
where we have disregarded the diffusion factor in the BFKL expression and consider $`𝒓^21/Q^21/m_{W,Z}^2`$. The result $`R_{\mathrm{scal}}1`$ leads to similar results between the all twist and leading twist calculation in the geometric scaling region. In fact, there is room for a slight suppression since $`(\gamma _{\mathrm{sat}}^{\mathrm{BFKL}}\gamma _{\mathrm{sat}}^{\mathrm{DLA}})0.13`$ and $`Q_{\mathrm{sat}}^2/m_{W,Z}^2<1`$. On the other hand, the argument in Ref. JAMAL  is exactly true in the comparison between CGC and GBW models. In this case, as we will discuss in what follows, $`\gamma _{\mathrm{sat}}^{\mathrm{GBW}}=1`$ in the geometric scaling region and a sizeable enhancement between CGC and GBW is really observed.

The qualitative analysis presented above can also help us to understand the enhancement of the BGBK model (with DGLAP evolution) in relation to GBW saturation model. We can think in an effective “anomalous dimension” for the saturation model in the geometric scaling regime: from Eq. (9), $`\sigma _{dip}^{\mathrm{GBW}}𝒓^2Q_{\mathrm{sat}}^2/4`$ and then one verifies that $`\gamma _{\mathrm{sat}}^{\mathrm{GBW}}=1`$. In this case, the ratio now reads as,

$`R_{\mathrm{scal}}={\displaystyle \frac{\sigma _{dip}(\mathrm{BGBK})}{\sigma _{dip}(\mathrm{GBW})}}\left[{\displaystyle \frac{Q_{\mathrm{sat}}^2}{m_{W,Z}^2}}\right]^{\left(\gamma _{\mathrm{sat}}^{\mathrm{DLA}}\gamma _{\mathrm{sat}}^{\mathrm{GBW}}\right)}1,`$ (15)
whereas at lower energies (beyond the geometric scaling regime) the DGLAP anomalous dimension takes the usual value and the ratio is closer to one, as expected.

As a summary, we have presented the results for the ultrahigh energy neutrino cross section considering the current phenomenology on non-linear QCD dynamics, which play an important role in the correct extrapolation to the regime of very high gluon density. The saturation model with QCD evolution was studied as well as the physics of BFKL dynamics at geometric scaling region. The latter presents suppression in the cross section at $`E_\nu >10^{12}`$ GeV, whereas the BGBK model is consistent with NLO DGLAP calculations. An outstanding advantage in the color dipole framework is the possibility to provide analytical expressions for the extrapolation to ultrahigh energies, towards the black limit of the nucleon target. This cannot be taken into account in a leading twist approach as the QCD-improved parton model. Then these saturation (all twist) approaches play an important role, since a precise knowledge about the neutrino interactions and production rates is essential for estimating background, expected fluxes and detection probabilities KMSNEUT ; STASTO . In particular, reliable estimates for UHE neutrino cross section are strongly needed for the planned experiments KUSENKO  to detect them via nearly horizontal air showers in the Earth atmosphere, as discussed in FIORE .

###### Acknowledgements.

The author is grateful for the warm hospitality and financial support of the High Energy Group in The Abdus Salam International Centre for Theoretical Physics (ICTP) at Trieste and CERN Theory Division at Geneve, where part of this work was performed. The author also would like to thank Alexei Smirnov (ICTP, Trieste) and Krysztof Kutak (II Institute for Theoretical Physics, Hamburg) for helpful discussions. This work was partially supported by CNPq and FAPERGS, Brazil.


text:
# Scaling of dynamics in 2d semi-dilute polymer solutions

## Abstract

We study the viscoelasticity of surface polymer monolayers by measuring the dynamics of thermal concentration fluctuations with surface light scattering. For various systems of proteins and synthetic polymers we find a semi-dilute regime in which both the elastic and viscous components of the dilational modulus show power law dependencies on the concentration. Surprisingly there is a universal relationship between the exponents for the two components: the viscosity scales with a power double that of the elasticity. These results cannot be explained on the basis of theory developed for bulk systems, and a simple explanation for the singular 2d behavior is suggested.

Polymers are flexible long chain molecules, of outstanding importance in diverse fields from processing of materials to biological activity. It is well known that it is possible to effectively confine some polymers to two dimensions, for example by anchoring each monomer to the interface between immiscible fluids Gaines (1960). Polymer dynamics in two dimensions (2d) has remained relatively unexplored compared to bulk solutions, despite there being important examples both in life sciences and in technology where polymer molecules are confined to a plane. We consider the simplest situation where the chains lie flat and are free to move only in a 2d space. Such a surface layer exerts a lateral osmotic pressure $`\mathrm{\Pi }=\gamma _0\gamma `$, which is the amount by which the surface tension $`\gamma `$ is reduced compared to the free interface tension $`\gamma _0`$. This pressure becomes significant above overlap of single chains, where it has a power-law dependence on the concentration that has been known for some time Vilanove and Rondelez (1980). Like its counterpart in 3d bulk solutions, this equilibrium property is well explained within the picture of the polymer semi-dilute regime as an ensemble of independent “blobs” de Gennes (1979). Despite a few decades of experiments, very little is known instead about the physical nature of the divergence of viscosity and elasticity in 2d. These dynamical properties control important processes like foam drainage Stone et al. (2003) and stabilization and flow in emulsions Buzza et al. (1995) when polymers are used as surfactants. The theoretical framework to describe complex non-Newtonian flow and the dynamics of chains is well understood in 3d Doi and Edwards (1986), but cannot be applied straightforwardly in 2d. How does the confinement to the surface, specifically a 2d effect, affect the dynamics of the system?

This Letter presents a study of the dynamics of concentration fluctuations in a wide range of polymer monolayers, performed with surface dynamic light scattering. Polymers are spread onto an interface in dilute conditions, and cannot subsequently submerge and re-surface. Under these conditions polymer chains are unable to cross each other. This is the case for many monolayer systems studied in the literature. The technique used in this work, developed by Langevin and others, see the monograph Langevin (1992), is the only existing probe of thermal concentration fluctuations in monolayers. Despite its use on polymer monolayers by various groups various , some issues concerning the data analysis have been resolved only recently, hence the data which is available for a quantitative comparison is very limited.

Isolated polymer chains are usually modelled as random walks with potential interactions de Gennes (1979). The Flory exponent $`\nu `$ relates the number of monomers $`N`$ to the radius of gyration: $`R_gN^\nu `$. When the concentration in a polymer solution is increased so that individual chains are forced to overlap, the system enters a semi-dilute regime that lasts until the monomer fraction is very high. The semi-dilute region is important because for sufficiently long chains it covers a wide range of concentrations. The equilibrium properties of polymers in this regime are given both in 2d and 3d by scaling laws, with exponents related to $`\nu `$. In particular, fluctuations of the density of monomers are correlated over a length $`\xi `$, defining a region known as a “blob” de Gennes (1979). Inside the blob the chain has the same statistics as an isolated single chain, but different blobs are statistically independent. In 2d $`\xi `$ scales with the concentration $`\mathrm{\Gamma }`$ as:

$$\xi \mathrm{\Gamma }^{\nu /(12\nu )},$$
(1)
decreasing from a value of the order of $`R_g`$ at the overlap concentration to the monomer size. Scaling of the osmotic pressure in a monolayer was first shown experimentally in ref. Vilanove and Rondelez (1980):

$$\mathrm{\Pi }\xi ^2\mathrm{\Gamma }^{y_{eq}},\text{where}y_{eq}=2\nu /(2\nu 1).$$
(2)
Figure 1 shows equilibrium data for some of the systems studied in this work. In 2d there are well known limiting regimes: The “good solvent”, in which the chain behaves as a self-avoiding random walk, hence $`\nu =3/4`$ and $`y_{eq}=3`$, and the “$`\theta `$ conditions” where the effects of excluded volume balance the monomer-monomer preferential attraction. Here the Flory exponent is predicted to be $`\nu =4/7`$ Coniglio et al. (1987) and the power law exponent increases to $`y_{eq}=8`$, making this isotherm much steeper and compact. In contrast to the bulk, where the polymer chains can interpenetrate, in the planar 2d geometry the polymers are segregated and attain disk-like configurations de Gennes (1979). Chains are corralled by their neighbors, as was observed with DNA molecules confined to a plane in ref. Maier and Rädler (1999).

In general, the response to a deformation in an isotropic 2d material is characterized by two elastic moduli: changes in area are controlled by the dilation modulus $`\epsilon `$ and changes in shape by the shear modulus $`G`$ Miller et al. (1996). Polymer monolayers in the semi-dilute regime are fluid-like and the shear modulus is negligible Cicuta et al. (2003). In these conditions it is common practice to determine the equilibrium dilational modulus $`\epsilon _{eq}`$ (the 2d analog of the bulk modulus) from measurements of pressure as a function of area: $`\epsilon _{eq}=d\mathrm{\Pi }/d\mathrm{ln}A`$. By definition, $`\epsilon _{eq}`$ has the same scaling properties as the osmotic pressure, described in Eq. 2. Power laws of $`\epsilon _{eq}`$ as a function of concentration are shown in Fig. 1.

The dilational modulus can also be measured dynamically, in an experiment where the surface area oscillates in time with a frequency $`\omega `$. Then the dynamic complex modulus $`\epsilon ^{}=\epsilon ^{}+i\omega \epsilon ^{\prime \prime }`$ is probed, where $`\epsilon ^{}`$ is the elastic component of the response modulus and $`\epsilon ^{\prime \prime }`$ is the dilational viscosity. $`\epsilon ^{}`$ can be accessed with the surface quasi-elastic light scattering (SQELS) technique, shown in Fig. 2. SQELS measures the time correlation function of light scattered from thermal surface roughness, which acts as a phase grating Langevin (1992). These out of plane fluctuations are underdamped waves with a frequency $`\omega `$, and their motion is affected by the presence of a surface film. To recover surface viscoelaticity, the surface wave dispersion relation $`D(\omega )`$ that relates the wave frequency to the wavelength has to be known:

$`D(\omega )`$ $`=\left[\epsilon ^{}q^2+i\omega \eta \left(q+m\right)\right]\left[\gamma q^2+i\omega \eta \left(q+m\right){\displaystyle \frac{\rho \omega ^2}{q}}\right]+`$
$`\left[i\omega \eta \left(mq\right)\right]^2,\text{ where }m=\sqrt{q^2+i{\displaystyle \frac{\omega \rho }{\eta }}}`$
$`\text{and }\mathrm{Re}(m)>0,`$ (3)
$`\eta `$ is the subphase Newtonian viscosity and $`\rho `$ is the subphase density. Buzza Buzza (2002) recently proved that the model Eq. 3, which is commonly found in the literature Lucassen-Reynders and Lucassen (1969); Langevin (1992), is correct but that even under dynamical conditions the surface tension $`\gamma `$ should be considered as a real quantity, equal to the equilibrium static surface tension 111 The model often used in the literature is over-parametrized. Data analysis with the correct physical parameters is very robust.. As derived in ref. Langevin and Bouchiat (1971), the spectrum of light scattered by thermal roughness is given by:

$`P_q(\omega )={\displaystyle \frac{k_BT}{\pi \omega }}\mathrm{Im}\left[{\displaystyle \frac{i\omega \eta (m+q)+\epsilon q^2}{D(\omega )}}\right].`$ (4)
SQELS data is fitted with Eq. 3 and Eq. 4, with only three physical parameters to be determined. Details of our experimental methods, including calibration, data analysis and limits of the technique, are described in ref. Cicuta and Hopkinson (2003).

In polymer systems, different models are appropriate depending on the lengthscales and timescales that are being observed. The lengthscales $`q^1`$ probed with SQELS satisfy $`q\xi <1`$, so the semidilute solution is expected to be in a “macroscopic” regime where it behaves like a gel de Gennes (1976). The gel’s response is characterized by a rigidity modulus, $`\epsilon ^{}`$, scaling like the number of contacts, hence proportional to the osmotic pressure, and by a viscous dissipation modulus $`\omega \epsilon ^{\prime \prime }`$ describing the friction involving the monomers and the solvent. This is well established in three dimensional solutions, where both a fast and a slow relaxation mode can be probed simultaneously, for example with dynamic light scattering Brown (1993). The fast relaxation can be described within the Rouse model and is related to a cooperative diffusion timescale $`\tau _{coop}`$ de Gennes (1976). This describes fluctuations that become faster as the concentration increases and the correlation length $`\xi `$ decreases. The slow relaxation is related to the self diffusion coefficient for reptation and is determined by the time $`\tau _{rept}`$ required for the chain to diffuse along its length. $`\tau _{rept}`$ is an increasing function of the concentration. The approach of directly measuring the frequency spectrum of light scattered by concentration fluctuations is not possible in monolayers because there is insufficient optical contrast, and dynamics can be probed only by indirect methods such as SQELS.

Figure 3 shows the components of the monolayer viscoelatic moduli, measured with SQELS for different polymer monolayers. Both the elastic and viscous components of $`\epsilon ^{}`$ exhibit a power law dependence on the concentration:

$$\epsilon ^{}\mathrm{\Gamma }^{y_\epsilon ^{}}\text{and}\epsilon ^{\prime \prime }\mathrm{\Gamma }^{y_{\epsilon ^{\prime \prime }}}$$
(5)
A similar scaling behavior was recently reported by Monroy et al. Monroy et al. (1999) but has not been explained 222An attempt in ref. Monroy et al. (1999) to understand the scaling of $`\epsilon ^{\prime \prime }`$ within the Rouse model predictions for $`\tau _{coop}`$ gave a dimensionally wrong result, which clearly has no physical meaning.. Results for all the systems considered in this work are summarized in Table 1. Monolayers of very different compositions, comprising both synthetic polymers and proteins, have been studied, to cover the widest possible range of values of $`\nu `$. The temperature, ionic strength and $`p`$H of the liquid subphase are all controlled, as they affect the chain configuration at the surface.

This data contains a wealth of information. Looking at Fig. 3, it can be clearly seen that for each monolayer the viscosity scales with a higher power than the elastic modulus. On general dimensionality terms, the viscosity can be regarded as the product of the modulus driving the relaxation and a characteristic time de Gennes (1979)$`\epsilon ^{\prime \prime }\epsilon ^{}\tau `$. From this it immediately follows that the timescale $`\tau `$ of these relaxations becomes longer the higher the concentration, hence a slow mode is being probed. The character of this mode can be read from Figure 4, where the power law exponents for $`\epsilon ^{}`$ defined in Eq. 5 are plotted against the equilibrium exponent. It is not surprising to find that $`y_\epsilon ^{}=y_{eq}`$, since the thermal fluctuations probed with SQELS fall within linear response theory. What is of interest is the apparently universally valid relation between the exponents for the elastic and the viscous components of $`\epsilon ^{}`$: $`y_{\epsilon ^{\prime \prime }}=2y_\epsilon ^{}`$. It finally follows that the timescale of these fluctuations, which we shall call $`\tau _{2d}`$, scales like $`\tau _{2d}\xi ^2`$.

We are only able to speculate on the physical origin of this timescale. A simple explanation is that the slowing down of $`\tau _{2d}`$ as the concentration increases is due to friction between the increased number of statistically independent blobs. The number of contacts in between blobs is proportional to $`\xi ^2`$, giving:

$$\epsilon ^{\prime \prime }\xi ^4\mathrm{\Gamma }^{2y_\epsilon ^{}}$$
(6)
which describes the data of Fig. 4 very well. An analogy can be drawn between the diverging viscosity of the close-packed arrangement of blobs in the semi-dilute regime and the case of diverging viscosity of spheres at high packing density Brady (1993), both being determined by the number of contacts.

A question that arises is why this mode is not seen in 3d. In bulk solutions in a good solvent a dynamical mode analogous to that just described for 2d would have a characteristic relaxation time scaling with the concentration $`\varphi `$ like $`\tau \xi ^3\varphi ^{9/4}`$. This is a higher power of the concentration compared to reptation ($`\tau _{rept}\varphi ^{1.5}`$)de Gennes (1979), in agreement with the well known result that in 3d solutions the fluctuations relax via the self diffusion of the polymer chain and not the mechanism outlined above. In 2d relaxation by reptation is hindered by an effective confinement provided by neighboring chains.

In summary, it has been shown that the thermal concentration fluctuations in the semi-dilute regime of polymer monolayers reported in the present paper and in ref. Monroy et al. (1999) have a slow mode of decay with a timescale that had not been previously considered and that describes a specifically 2d process.

###### Acknowledgements.

We thank E.M.Terentjev for very useful comments and discussions.


text:
# The 𝑩_𝒔 and 𝑫_𝒔 decay constants in 3 flavor lattice QCD

## Abstract

Capitalizing on recent advances in lattice QCD, we present a calculation of the leptonic decay constants $`f_{B_s}`$ and $`f_{D_s}`$ that includes effects of one strange sea quark and two light sea quarks. The discretization errors of improved staggered fermion actions are small enough to simulate with 3 dynamical flavors on lattices with spacings around 0.1 fm using present computer resources. By shedding the quenched approximation and the associated lattice scale ambiguity, lattice QCD greatly increases its predictive power. NRQCD is used to simulate heavy quarks with masses between $`1.5m_c`$ and $`m_b`$. We arrive at the following results: $`f_{B_s}=260\pm 7\pm 26\pm 8\pm 5`$ MeV and $`f_{D_s}=290\pm 20\pm 29\pm 29\pm 6`$ MeV. The first quoted error is the statistical uncertainty, and the rest estimate the sizes of higher order terms neglected in this calculation. All of these uncertainties are systematically improvable by including another order in the weak coupling expansion, the nonrelativistic expansion, or the Symanzik improvement program.

We present the first complete calculation of the $`B_s`$ and $`D_s`$ decay constants with 3 flavors of sea quarks with small masses. The $`D_s`$ decay constant, $`f_{D_s}`$, has been measured Chadha:1998zh ; Heister:2002fp , and CLEO-c promises to reduce the experimental errors significantly. The comparison of experimental and lattice results will be a vital test of lattice QCD. The $`B_s`$ and $`B_d`$ decay constants, $`f_{B_s}`$ and $`f_B`$, are necessary in order to constrain $`V_{td}`$ via $`\overline{B^0}B^0`$ mixing. Since experimental measurements of $`f_{B_s}`$ remain elusive, an accurate lattice calculation is critical for improving phenomenological tests of CKM unitarity.

The calculation presented below makes use of lattice Monte Carlo simulations, done by the MILC collaboration Bernard:2001av , which include the proper sea quark content: one dynamical strange quark and two flavors of dynamical quarks with masses as light as $`m_s/4`$. The correct number of flavors is necessary in order for the lattice theory to have the same $`\beta `$-function as real QCD. Only then can we expect, in principle, lattice results to agree with experimental measurements.

Inclusion of light up and down sea quarks is essential for accurate lattice phenomenology. The innovation which allows this on present computers is an improved staggered discretization of the light quark action Bernard:1998gm ; Lepage:1998id ; Bernard:1998mz ; Lepage:1998vj ; Orginos:1998ue ; Toussaint:1998sa ; Orginos:1999cr ; Bernard:1999xx . For each dynamical flavor the fourth root of the fermion determinant is used in the Monte Carlo updating algorithm. Despite some open theoretical issues concerning staggered fermion algorithms, these calculations are free of ambiguities present in quenched simulations, provided the $`u`$ and $`d`$ sea quark masses are light enough. In practice, calculations involving unstable particles or multi-hadronic final states require new or refined techniques before they can be reliably simulated; however, there are several quantities for which numerical calculation is straightforward. Ref. Davies:2003ik  presents lattice results for a variety of these “gold-plated” quantities, including $`f_\pi `$ and $`f_K`$, and several splittings in the $`\mathrm{{\rm Y}}`$ spectrum. The change from $`n_f=0`$ results, which differ substantially from experiment, to $`n_f=3`$ results, which agree with experiment up to the 3% lattice uncertainties, suggests that the improved staggered fermion method correctly simulates QCD. Heavy-light pseudoscalar leptonic decay constants, presented here, fit into the “gold-plated” category Davies:2003ik .

Our calculation of $`f_{B_s}`$ and $`f_{D_s}`$ uses standard lattice QCD methods. Correlation functions were computed using a subset of gauge field configurations generated by the MILC collaboration. The configurations include the effects of 2 dynamical light quarks with equal bare mass, denoted by $`m_{\mathrm{}}^{\mathrm{sea}}`$, and 1 dynamical strange quark, with bare mass $`m_s^{\mathrm{sea}}`$. The lattices we used have spacings of about 1/8 fm and volumes of about $`(2.5\mathrm{fm})^3\times 8.0`$ fm. We focus on the configurations where $`m_{\mathrm{}}^{\mathrm{sea}}/m_s^{\mathrm{sea}}=1/5`$ and 2/5. As we discuss below, the physical strange quark mass, $`m_s`$, obtained from the light hadron spectrum is actually $`4/5`$ of $`m_s^{\mathrm{sea}}`$, so the light sea quark masses are approximately $`m_s/4`$ and $`m_s/2`$.

Details of the Monte Carlo simulations which generated the ensemble of gauge fields are given in Bernard:2001av . In this work we use a level splitting in the $`\mathrm{{\rm Y}}`$ spectrum, e.g. the $`\mathrm{{\rm Y}}(2\mathrm{S}1\mathrm{S})`$ splitting, to determine the lattice spacing, instead of the length scale $`r_1`$, derived from the static quark potential, which was used in Bernard:2001av . The $`\mathrm{{\rm Y}}`$ spectrum will be presented in detail in a future publication Gray:UpsSpect ; however, we note here that as the light sea quark mass is increased, the $`\mathrm{{\rm Y}}(1\mathrm{P}1\mathrm{S})`$ splitting becomes slightly smaller than experiment when the lattice spacing is set by $`\mathrm{{\rm Y}}(2\mathrm{S}1\mathrm{S})`$. With statistical errors between $`12\%`$, the difference between experiment and lattice values for the $`1\mathrm{P}1\mathrm{S}`$ splitting is $`0\sigma `$, $`1\sigma `$, and $`2\sigma `$ for $`m_{\mathrm{}}^{\mathrm{sea}}/m_s`$ 1/4, 1/2, and 3/4, respectively. Except for estimating sea quark mass effects, we use the $`m_{\mathrm{}}^{\mathrm{sea}}/m_s1/4`$ lattice to obtain our results.

The analysis of Bernard:2001av  has also been updated with respect to determining the quark mass corresponding to the physical strange quark mass sector. Rather than using $`\overline{s}s`$ mesons which are either unstable (the $`\varphi `$) or do not exist (the pseudoscalar) we use the results of the partially quenched chiral perturbation theory analysis of $`m_K`$ and $`m_\pi `$ Davies:2003ik ; Aubin:2002ss . On lattices with $`au_0m_{\mathrm{}}^{\mathrm{sea}}=0.007,0.01`$, and $`0.02`$, the physical strange sector is obtained with $`au_0m_s^{\mathrm{val}}=0.040`$. (Here we keep explicit a factor of the mean-link $`u_0`$, which was absorbed into the definition of the light bare quark mass in Bernard:2001av .)

Table 1 lists the bare heavy quark masses used in this work, along with the corresponding NRQCD stabilization parameter $`n`$. For each mass, we estimate the heavy-light meson mass two ways: the “kinetic mass” $`m_{H_s}^{\mathrm{kin}}`$ is extracted from finite momentum correlators using the meson dispersion relation, and the “perturbative mass” $`m_{H_s}^{\mathrm{pert}}`$ is estimated from the zero momentum correlator and the one-loop heavy quark mass renormalization. We find that the bare heavy quark mass $`aM_0=2.8`$ produces the correct experimental mass for the $`\mathrm{{\rm Y}}`$ Gray:2002vk  and the $`B_s`$ within statistical errors.

When we construct correlation functions for heavy-light mesons, we use the equivalence between staggered and naive fermion propagators. The operators we use are the same ones as with Wilson-like discretizations or in the continuum limit. This method for computing the $`B_s`$ mass and decay constant has been presented, along with tests on quenched lattices, in recent work Wingate:2002fh .

The $`B_s`$ decay constant is defined through the axial vector matrix element; $`0|A_\mu |B_s(p_\mu )=f_{B_s}p_\mu `$. Here we use only the temporal component. Up through $`O(1/M_0)`$ three lattice operators contribute to $`A_0`$, the leading-order current, $`J_0`$, and two sub-leading currents $`J_1`$ and $`J_2`$; the matching from the lattice to the continuum is done in perturbation theory Morningstar:1998ep . A one-loop calculation yields

$$A_0=(1+\alpha _s\stackrel{~}{\rho }_0)J_0^{(0)}+(1+\alpha _s\rho _1)J_0^{(1,\mathrm{sub})}+\alpha _s\rho _2J_0^{(2,\mathrm{sub})}$$
(1)
where $`J_0^{(1,\mathrm{sub})}=J_0^{(1)}\alpha _s\zeta _{10}J_0^{(0)}`$, and so too for $`J_0^{(2,\mathrm{sub})}`$. These operators explicitly subtract, through one-loop, the power-law mixing of $`J_0^{(0)}`$ with $`J_0^{(1)}`$ and $`J_0^{(2)}`$ Collins:2000ix . The perturbative calculation which determines $`\stackrel{~}{\rho }_0,\rho _1,\rho _2,\zeta _{10}`$, and $`\zeta _{20}`$ will be presented separately Gulez:HLPT .

From fits to correlation functions, we extract the combination $`\mathrm{\Phi }_{H_s}f_{H_s}\sqrt{m_{H_s}}`$. Let us denote by $`\mathrm{\Phi }^{(i)}`$ the contribution of $`J_0^{(i)}`$ to $`\mathrm{\Phi }_{H_s}`$. For this calculation the matrix elements are computed for mesons at rest, so $`\mathrm{\Phi }^{(1)}=\mathrm{\Phi }^{(2)}`$. Table 2 summarizes fits to the numerical data, converted to physical units using $`1/a=1.59(2)`$ GeV for the $`m_{\mathrm{}}^{\mathrm{sea}}m_s/4`$ lattice and $`1/a=1.61(2)`$ GeV for the $`m_{\mathrm{}}^{\mathrm{sea}}m_s/2`$ lattice. (The quoted uncertainties come from statistical and fitting uncertainties in the $`\mathrm{{\rm Y}}(2\mathrm{S}1\mathrm{S})`$ splitting Gray:UpsSpect .) By comparing $`\mathrm{\Phi }^{(1)}/\mathrm{\Phi }^{(0)}`$ with $`\mathrm{\Phi }^{(1,\mathrm{sub})}/\mathrm{\Phi }^{(0)}`$ one can observe the sizable power law mixing of $`J^{(0)}`$ with $`J^{(1)}`$. The expression (1) absorbs the mixing back into the term proportional to $`J^{(0)}`$, so $`\mathrm{\Phi }^{(1,\mathrm{sub})}/\mathrm{\Phi }^{(0)}`$ represents the physical contribution of $`1/M_0`$ terms to $`\mathrm{\Phi }_{H_s}`$ up to two-loop corrections. The 4% contribution from the $`1/M_0`$ operator we see for the $`B_s`$ ($`aM_0=2.8`$) is the same size seen in quenched studies over a range of lattice spacings Collins:2000ix .

Figure 1 shows the heavy quark mass dependence of $`\mathrm{\Phi }_{H_s}`$ on $`m_{H_s}`$ (plotted as squares). The data are fit well by

$$\mathrm{\Phi }_{H_s}=\mathrm{\Phi }_{H_s}^{\mathrm{stat}}\left(1+\frac{C_1}{m_{H_s}}+\frac{C_2}{m_{H_s}^2}\right)$$
(2)
with a correlated $`\chi ^2`$ per degree-of-freedom of $`0.7/2`$. We find $`\mathrm{\Phi }_{H_s}^{\mathrm{stat}}=0.85(4)`$ GeV3/2, $`C_1=1.82(20)`$ GeV, and $`C_2=1.59(35)\mathrm{GeV}^2`$. It is important to note that most of the mass dependence of $`\mathrm{\Phi }_{H_s}`$ comes through the action and not through the $`1/M_0`$ currents. A fit to the leading order $`\mathrm{\Phi }^{LO}(1+\alpha _s\stackrel{~}{\rho }_0)\mathrm{\Phi }^{(0)}`$ (Table 2) yields similar fit parameters: $`\mathrm{\Phi }^{\mathrm{stat},LO}=0.83(4)`$ GeV3/2, $`C_1^{LO}=1.53(19)`$ GeV, and $`C_2^{LO}=1.32(32)\mathrm{GeV}^2`$. Since the action is responsible for most of the $`m_{H_s}`$ dependence and is accurate through $`O(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_Q^2)`$, performing the operator matching through $`O(1/M_0^2)`$ should yield only small corrections to the values for $`\mathrm{\Phi }_{H_s}^{\mathrm{stat}}`$, $`C_1`$, and $`C_2`$ listed above.

We use the fit above to interpolate $`\mathrm{\Phi }_{H_s}`$ slightly to the physical value of $`1/m_{B_s}`$. The result is $`f_{B_s}=260(7)`$ MeV, where the quoted statistical error combines statistical errors from $`1/a`$ and $`a^{3/2}\mathrm{\Phi }_{B_s}`$. Uncertainties in $`m_{H_s}^{\mathrm{kin}}`$ are small compared to the other uncertainties in the fit. The systematic uncertainties due to the neglect of higher order terms are estimated by assuming coefficients of $`O(1)`$. For example, two-loop terms omitted in (1) are estimated to be 10% effects, taking the coupling constant defined through the plaquette $`\alpha _s=\alpha _P^{n_f=3}(2/a)=0.32`$ Davies:2002mv . This is the largest systematic uncertainty. Leading discretization errors are $`O(\alpha _sa^2\mathrm{\Lambda }_{\mathrm{QCD}}^2)`$, where $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ is the typical scale of nonperturbative dynamics. Taking $`1/a=1.6`$ GeV and $`\mathrm{\Lambda }_{\mathrm{QCD}}=400`$ MeV implies 2% cutoff effects. The operator matching neglects terms $`O(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_b^2)`$ or approximately 1%. We note that the coefficient of the spin-dependent term in the action is the tree level one, giving rise to an $`O(\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/m_Q)`$ error in the action, which is about 3% for bottom quarks. This does not necessarily translate into a comparably large error in decay constants where the spin dependence is not obvious; however, to be conservative, we quote an overall 3% uncertainty due to relativistic corrections.

Extrapolating the fit described above to the physical value for $`1/m_{D_s}`$ requires care since higher order terms in $`1/m_{H_s}`$ become increasingly important. The fit above extrapolates to $`f_{D_s}=290(10)`$ MeV. We use a Bayesian analysis to estimate possible effects due to higher order terms. Allowing terms like $`C_n/m_{H_s}^n`$ with $`n3`$ in the fit (2), with Gaussian priors for $`C_n/\mathrm{GeV}^n=0\pm \delta `$ and $`\delta =14`$, we find higher order terms could lead to a 20 MeV error in $`f_{D_s}`$. This is taken as the statistical and fitting uncertainty. The $`O(\alpha _s^2)`$ perturbative error and the $`O(\alpha _sa^2\mathrm{\Lambda }_{\mathrm{QCD}}^2)`$ discretization error are again estimated to be 10% and 2%, respectively. For charmed mesons $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c=1/4`$, so the $`O(\mathrm{\Lambda }_{\mathrm{QCD}}^2/m_Q^2)`$ corrections to the matching are estimated to be 6%. As mentioned above, the $`O(\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/m_Q)`$ error in the spin dependent term in the action probably does not lead to a proportional error in the decay constant. Nevertheless, we quote 10% as a conservative estimate of possible relativistic corrections.

Figure 2 shows $`\mathrm{\Phi }_{H_s}`$ vs. light sea quark mass for two heavy quark masses. We show the results where the lattice spacing is set from the $`\mathrm{{\rm Y}}(2\mathrm{S}1\mathrm{S})`$ splitting and the same data adjusted as if the spacing were set from the $`\mathrm{{\rm Y}}(1\mathrm{P}1\mathrm{S})`$ splitting. No sea quark mass dependence is observed for $`aM_0=2.8`$, nor is a significant dependence observed for $`aM_0=1.0`$ when the spacing is set with $`\mathrm{{\rm Y}}(2\mathrm{S}1\mathrm{S})`$. Only because there is a 2% ambiguity in setting the scale with $`m_{\mathrm{}}^{\mathrm{sea}}m_s/2`$ does there appear to be some sea quark mass effect in the $`aM_0=1.0`$ data with the $`\mathrm{{\rm Y}}(1\mathrm{P}1\mathrm{S})`$ spacing. No such scale setting ambiguity exists for the data with $`m_{\mathrm{}}^{\mathrm{sea}}m_s/4`$, which is where our main results are obtained.

Experimental measurements of the leptonic decay constants for $`B_s`$ and $`D_s`$ are challenging, so these are rare cases where lattice QCD can lead experiment. The most recent and precise experimental results for $`f_{D_s}`$ are $`280\pm 17\pm 25\pm 34`$ MeV Chadha:1998zh  and $`285\pm 19\pm 40`$ MeV Heister:2002fp , which agree well with the calculation presented here. Shrinking the experimental uncertainties on $`f_{D_s}`$ is a major goal of the CLEO-c program.

Turning to comparison with the existing literature on lattice QCD, we quote recent world averages and comment on another new result. The $`n_f=3`$ results presented here have central values larger than lattice calculations with $`n_f=0`$ or 2. For example, recent averages Ryan:2001ej  of quenched results are $`f_{B_s}^{n_f=0}=200(20)`$ MeV and $`f_{D_s}^{n_f=0}=230(14)`$ MeV, significantly lower than our results. The 2 flavor world averages, $`f_{B_s}^{n_f=2}=230(30)`$ MeV and $`f_{D_s}^{n_f=2}=250(30)`$ MeV, are higher than the quenched and agree within the quoted uncertainties.

Recently the JLQCD collaboration reported a lattice result using two dynamical flavors of improved Wilson fermions with mass between $`0.7m_s`$ and $`2.9m_s`$ Aoki:2003xb . They quote $`f_{B_s}=215(9)\left(\genfrac{}{}{0pt}{}{+0}{2}\right)(13)\left(\genfrac{}{}{0pt}{}{+6}{0}\right)`$ MeV, with the first error statistical, second from chiral extrapolation of the sea quark mass, third from finite lattice spacing combined with truncation of NRQCD and perturbative expansions, and the fourth from ambiguity in setting the strange quark mass. This calculation uses much larger quark masses than ours, uses unstable hadron masses to set the lattice spacing and $`m_s`$, and does not include a dynamical strange quark. Further work will be required to determine which, if any, of these accounts for the differences with our results. For instance, it would be interesting to calculate $`a^1`$ from several $`\mathrm{{\rm Y}}`$ splittings on the configurations of Aoki:2003xb  and compare $`f_{B_s}`$ using those scales to the value quoted in their paper, which uses the scale set by $`m_\rho `$.

Recent sum rule calculations agree within errors for the $`B_s`$ decay constant: e.g. $`f_{B_s}^{\mathrm{s}.\mathrm{r}.}=236(30)`$ MeV Jamin:2001fw ; Narison:2001pu . On the other hand, they do not calculate an increase in the decay constant as the heavy quark mass decreases: e.g. $`f_{D_s}^{\mathrm{s}.\mathrm{r}.}=235(24)`$ MeV Narison:2001pu .

To summarize, we have completed a calculation of the $`B_s`$ and $`D_s`$ decay constants using 3 flavor lattice QCD. The more realistic sea quark content allows a unique lattice spacing to be determined using one of several quantities Davies:2003ik , enabling more reliable prediction of quantities not yet measured experimentally. Our final results are

$`f_{B_s}`$ $`=`$ $`260\pm 7\pm 26\pm 8\pm 5\mathrm{MeV}`$
$`f_{D_s}`$ $`=`$ $`290\pm 20\pm 29\pm 29\pm 6\mathrm{MeV}.`$ (3)
The uncertainties quoted are respectively due to statistics and fitting, perturbation theory, relativistic corrections, and discretization effects. The result for the $`D_s`$ decay constant agrees with experimental determinations, and the result for the $`B_s`$ decay constant is a prediction for future experiments. Improvement of the lattice results requires a two-loop perturbative matching calculation or the use of fully nonperturbative methods. On the other hand, much of the perturbative uncertainty cancels in the ratio $`f_{B_s}/f_{B_d}`$. Work is underway to study $`f_{B_d}`$ using the methods discussed in this paper Wingate:2003ni .

Simulations were performed at NERSC. We thank the MILC collaboration for their gauge field configurations. This work was supported in part by the DOE, NSF, PPARC, and the EU.


text:
# Observation of near-critical reflection of internal waves in a stably stratified fluid

## I Introduction

The reflection of near-critical internal waves over sloping topography plays a crucial role in determining exchanges between the coastal ocean and the adjacent deep waters. Direct measurements of mixing in the ocean, using tracers ledwell12 , have vindicated decades of phenomenological and theoretical inferences. In particular, these measurements have shown that most of the vertical mixing is not taking place inside the ocean, but close to the boundaries and topographic features polzin97 . These results have directed attention to the possible role of internal wave reflection in the boundary-mixing process.

Internal waves have different properties of reflection from a rigid boundary than do sound or light waves phillips . Instead of following the familiar Snell’s law, internal waves reflect off a boundary such that the angle with respect to gravity direction is preserved upon reflection (Fig. 1). This peculiar reflection law leads to a concentration of the reflected energy density into a narrow ray tube upon reflection as displayed in Fig. 1. Theoretical descriptions of this reflection process have been framed largely in terms of linear and stationary wave dynamics phillips ; Gilbert93 . However, when the slope angle, $`\gamma `$, is equal to the incident wave angle, $`\beta `$, these restrictions lead to an unrealistic prediction: The reflected rays lie along the slope with an infinite amplitude and a vanishing group velocity. Theoretical results have recently healed this singularity by taking into account the role of transience and nonlinearity JFM .

Following preliminary oceanographic measurements Sandstom66 ; Ericksen82 , Eriksen eriksen98  has beautifully observed, near the bottom of a steep flank of a tall North Pacific Ocean seamount, an internal wave reflection process leading to a clear departure from a Garett-Munk model garettmunk  for wave frequencies at which ray and bottom slopes match. Several experimental facilities cacchione ; thorpe ; ivey ; imberger97 ; maasommeria ; Sutherland2000  have therefore been dedicated to the understanding of the internal wave reflection and associated instabilities.

However, results of a controlled laboratory experiment close to the critical conditions are still lacking since previous ones with a moving paddle cacchione ; thorpe ; ivey  at one end of a very long tank, or by the vibration of the tank itself maasommeria , does not generate a clear incident wave-beam as needed for a careful study. In addition, a direct comparison with the recent and complete nonlinear theory near the critical reflection would be possible. Finally, the goal is to improve the understanding of the possible mixing mechanism near the sloping topography of ocean as very recently initiated by MacPhee and Kunze McPhee02  by exhibiting the instabilities mechanism leading to mixing. The paper is organized as follows. The experimental setup is described and carefully explained in Sec. II. The main experimental results are presented in Sec. III, and comparisons with the theory is also provided. Finally, Sec. IV contains conclusions and perspectives.

## II Experimental Setup

The experimental setup consists of a 38 cm long Plexiglas tank, 10 cm wide, filled up to 22 cm height with a linearly salt-stratified water obtained by the “double bucket” method oster . The choice of salt, Sodium Nitrate (NaNO3) snow, is motivated due to its highly solubility in water leading to a salty water viscosity close to the fresh one. Moreover, this salt allows to reach a strong stratification: The fluid density ( $`1\rho (z)1.2`$ g/cm3) measured at different altitudes ($`22z>0`$ cm) by a conductimetric probe leads to a linear vertical density profile of slope $`d\rho /dz0.0104`$ g/cm4. A rectangular Plexiglas sheet, $`3`$ mm thick and $`9.6`$ cm wide (to allow exchange of water) is introduced at one end of the tank with an angle $`\gamma =35^{}`$ with respect to the horizontal tank bottom (see Fig. 1) to create the reflective sloping boundary.

Internal waves are generated by a sinusoidal excitation provided by the vertical motion of a horizontal PVC plunging cylinder ($`3.1`$ cm in diameter and $`9.4`$ cm long). The cylinder is located roughly midway between the base tank and the free surface. This wavemaker is driven by an electromagnetic vibration exciter powered by a low frequency power supply. Optical measurements confirm anthony  that the cylinder motion is sinusoidal without distortions for vibrational frequencies $`0.2f0.5`$ Hz and maximal displacement amplitudes up to $`8.5`$ mm (peak to peak).

The well-known and non-intuitive dispersion relation of internal waves in an incompressible, inviscid and linearly stratified fluid reads Lighthill

$$\omega =\pm N\frac{k_{}}{𝐤}=\pm N\mathrm{sin}\beta ,$$
(1)
where $`k_{}`$ is the component of $`𝐤`$ perpendicular to $`z`$-axis, $`N=\sqrt{(g/\rho _0)\rho /z}`$ is the constant buoyancy (or Brunt-Väisälä) frequency, $`\rho (z)`$ the fluid density at altitude $`z`$, $`g=981`$ cm/s2 the acceleration of gravity, and $`\rho _01`$ g/cm3 a reference density. Thus, from Eq. (1), the wave frequency, $`\omega =2\pi f`$, determines the inclination angle $`\beta `$ of the phase surfaces with the vertical, and not the magnitude of the wave vector k. From Eq. (1), $`\beta `$ is also the angle between the group velocity $`𝐜_𝐠=\omega /𝐤`$ and the horizontal, since $`𝐜_𝐠𝐤`$. Thus, for a given stratification and frequency, internal waves of low (higher) frequency propagate at low (steeper) angle.

The outward radiation of energy is thus along four beams oriented at an angle $`\beta `$ with the horizontal, the familiar St Andrews Cross structure Mowbraraity . The beam propagating directly toward the slope has been singled out by adding a grid on the surface of water. This grid strongly damps the three other wave beams that propagate towards the free surface, and thereby prevents reflection of such beams. A planar wave pattern consisting of parallel rays is thus generated from the wavemaker. Although it is well-known that the spatial spectrum of waves generated by an oscillating cylinder is large ApplebyCrighton ; Voisin ; Sutherland2000 , it has been experimentally checked that the dominant wavelength is approximately equal to the cylinder diameter (see below). Moreover, measurements confirm that the low vibrational amplitudes of cylinder do not influence strongly the wavelength generated.

Different visualization methods are used to study the reflection mechanism of such internal gravity waves by a boundary layer. First, the usual shadowgraph technique mer  allows visualization of the qualitative and global 2D evolution of the incident and reflected waves. It involves projecting a point source of light through stratified water onto a screen behind the tank. The optical refractive index variations induced by the fluid density variations, allows to observe isodensity lines (or isopycnals) on the screen, located perpendicularly to the light source and parallel to the longest wall tank. It is therefore possible to measure the group velocity angle, $`\beta `$, and the phase velocity $`v_\phi `$ of the incident wave. For various frequencies of excitation, $`0.2f0.5`$ Hz, the angle $`\beta `$ of the St Andrews Cross is measured on a screen leading to a linear relation between $`\omega `$ and $`\mathrm{sin}\beta `$ as predicted by Eq. (1) with a slope of $`N=3.1\pm 0.1`$ rad/s, for the stratified fluid prepared as above. This value is in good agreement with the above static one obtained from the density profile with a conductimetric probe. The cut-off frequency is then $`f_c=N/(2\pi )0.5`$ Hz. By time of flight measurements between signals delivered by two photodiodes, 1 cm apart, each of 7 mm2 area, located along the propagation direction, a value $`v_\phi =0.6\pm 0.4`$ cm/s was obtained for an excitation frequency $`f=0.25`$ Hz. When these signals are cross-correlated by means of an spectrum analyzer, the averaged dephasing time leads to $`v_\phi =0.7\pm 0.3`$ cm/s, close to the previous value. The wavelength of the incident wave is thus $`\lambda =v_\phi /f3`$ cm, corresponding as expected to the oscillating cylinder’s diameter. However, neither this shadowgraph technique nor one mer  using passive tracers (fluorescein dye) is sufficiently sensitive to observe quantitative and local internal wave properties closed to the reflective boundary layer.

Accordingly, the classical Schlieren method mer ; sutherland  of visualization has been used. Let us just note that behind the tank, the light beam is refocused by a lens a small distance after a slit (instead of the usual knife blade to increase contrast) to filter the rays. The slit is oriented orthogonally to the slope generating straight horizontal fringe lines in the case of no excitation. The image of the observation field (strongly dependent of the $`7.5`$ cm diameter lens) is focused on the screen by a last lens. The internal wave, producing density disturbances, causes lines to distort, this distorting line pattern being recorded by a camera. Note that this experimental technique is sensitive to the index gradient, and therefore to the density gradient, orthogonal to the slit, i.e. parallel to the slope.

## III Experimental results and comparison with theory

The Schlieren technique allows to carefully observe quantitative and local internal wave properties during the reflection process. Critical reflection arises when an incident wave beam with an angle of propagation $`\beta `$ reflects off the slope of angle $`\gamma \beta `$, the reflected wave being then trapped along the plane slope. This corresponds to a critical frequency $`f_c=N\mathrm{sin}\gamma 0.28\pm 0.01`$ Hz, $`N`$ being equal to $`3.1\pm 0.1`$ rad/s for all experiments. It is possible to observe that the isodensity lines (isopycnals), initially horizontal without excitation, are bent for an excitation near $`f_c`$ ($`0.78f/f_c1.14`$), and fold over themselves along the length of the slope. Figures 2 show a time sequence of constant density surfaces, depicting sequential snapshots of the flow throughout one period of its development.

These pictures show strikingly the distortion of the isopycnals in the slightly subcritical case with $`f/f_c=0.78`$. In panel (a), the density disturbance is very small and one can observe essentially the initial horizontal background stratification. This background stratification, usually obtained with dye fluorescein, should be invisible with this schlieren method. However, as previously reported by MacPhee and Kunze McPhee02 , the comparison between shadowgraph and Dye fluorescein visualizations allows to identify these lines with isopycnals. In panel (b), the disturbance generated by the incident wave breaking against the slope begins to ‘fold-up’ the isopycnals. As time progresses (see panel c), wave overturning develops around a front: The buoyancy becomes statically unstable. This overturned region climbs along the slope as time continues, and the folded isopycnals collapse into turbulence that mixes the density field within the breaking region (see panel d). The maximum thickness of this reflected disturbance is of the order of 5 mm, and decreases with time as theoretically predicted JFM . Finally, the flow begins to relaminarize (panel e). One can check that panels (a) and (f) are almost identical, showing that the flow is entirely restratified (panel f) after one period of excitation.

Figures 2 have been analyzed with an image processing software (Scion-Image) to extract the isopycnals from the pictures. A typical experimental result for the distortion of isopycnals is reported in Fig. 3a, and is compared with a theoretical result in Fig. 3b. The analytic solution of the density field of the initial value problem in the critical case reads JFM

$$\rho =\rho _0\left[1\frac{N^2}{g}\left(z\mathrm{cos}\gamma \psi B\sqrt{\frac{𝐤\mathrm{sin}^2(2\beta )}{2\omega _+}}\right)\right]$$
(2)
where

$`B`$ $`=`$ $`\sqrt{{\displaystyle \frac{t}{z}}}J_1\mathrm{sin}\left(\omega _+t𝐤\mathrm{sin}(\beta +\gamma )x\right),`$ (3)
$`J_1`$ $`=`$ $`J_1\left(2\sqrt{2\omega _+\mathrm{cos}^2\beta 𝐤tz}\right),`$ (4)
$`J_1`$ being the Bessel function, $`\omega _+`$ the positive solution of Eq. (1), $`\psi `$ the maximum amplitude of the streamfunction and $`x`$ the horizontal coordinate. Figure 3 show a good qualitative agreement between experimental and theoretical results, the value for the time $`t`$ being arbitrarily chosen. Far from the slope, the density disturbance is very small and one sees essentially the initial background stratification. Closer to the slope, the disturbance folds up the isopycnals, and this produces a region of static instability.

Recording several isopycnals and using image processing, it is also possible to follow the temporal evolution of a single isopycnal during its overturning. As this phenomenon is periodic with a period $`T=1/f`$, it is possible to reconstruct from this temporal evolution the density profile picture at a given time $`t`$. This allows to follow the position, and therefore the propagation velocity of the front along the slope at different times. The front is defined as the inflexion point (represented by the star in Fig. 3b) of the followed isopycnal. Figure 4a shows, during two periods, the isopycnal front position along the slope as a function of time. The periodic evolution of this front position is clearly observed, and the local slope of curves in Fig. 4a allows to roughly measure the front velocity as a function of time, as reported in Fig. 4b.

The front velocity from its creation to its collapsing increases from $`0.5`$ cm/s up to $`3`$ cm/s. The front has thus traveled roughly $`4`$ cm in one period ($`4.5`$ s). This leads to an averaged front velocity in agreement with the phase speed measurement obtained from the shadowgraph visualizations.

The wavemaker frequency is now increased up to $`f=0.32`$ Hz, to have an incident planewave tilted with an angle $`\beta `$ steeper than the slope angle $`\gamma `$. In this slightly supercritical case ($`f/f_c=1.14`$), intrusions are still observed, but the density field doesn’t fold up so abruptly and does not lead anymore to overturning instability (see Fig. 5). Except the value of the frequency $`f`$, all others parameter values have been kept identical to the ones in Fig. 2. The instant of this snapshot has been chosen when the isopycnal distortion is the largest webpage . This distortion is clearly far from leading to an overturning instability as encountered in the subcritical case of Fig. 2c: The reflected wave is not trapped along the slope in the boundary layer, and consequently the isopycnals are not overturned. Moreover the density front velocity is measured roughly constant, $`0.5`$ cm/s, during two periods of vibration. Both differences confirm that the singularity appears only in the critical case JFM .

## IV Conclusions and perspectives

The Schlieren technique allows us to study the spatiotemporal evolution of the internal waves reflection close to the critical reflection, where nonlinear processes occur. The dynamics of isopycnals is then found in agreement with a recent nonlinear theory JFM . Moreover, this experiment confirms the theoretically predicted scenario for the transition to boundary-layer turbulence responsible for boundary mixing: The growth of a density perturbation produces a statically unstable density field which then overturns with small-scale fluctuations inside McPhee02 . Panel (d) of Fig. 2 is characteristic of the onset of turbulence triggered by overturning instability near the slope.

This turbulent mechanism is likely responsible for the formation of highly “stepped” temperature profiles near steep slopes in lakes Caldwell . Moreover, the formation of suspended sediment layers, called nepheloid layers, at continental slopes has been linked to critical angle reflection of internal waves white ; McPhee02 , and suggests that this reflection process plays an important role in the seawards transport of sediments in the fluid. A possible extension of this work, with a smaller slope angle to be closer to real oceanographic situations, would be to study and characterize the long time behavior and diffusion process of such particle layers toward the fluid interior.

###### Acknowledgements.

We warmly thank J.-C. Géminard and J. Sommeria for helpful suggestions. This work has been partially supported by the French Ministère de la Recherche grant ACI jeune chercheur-2001 N 21-31.