text:
# Mesoscopic effects in tunneling between parallel quantum wires

## I Introduction

Tunneling provides a powerful tool to probe electronic properties of matter. Its sensitivity to momentum-resolved spectral features is determined by geometrical details of the tunnel junction. For example, this sensitivity is completely lost when tunneling occurs via a point contact, whereas it is maximal for an extended, clean tunneling barrier. Since the experimental study of tunneling between two separately contacted, parallel, vertically separated two-dimensional (2D) electron systems became possible, electronic structure and interaction effects in low dimensions have been the subject of careful investigation. In the ideal case, conservation of canonical momentum in the plane of the 2D electron systems leads to sharp tunneling resonances; allowing for exploration of electronic subband energies, mapping of the 2D Fermi surface, and life-time measurements of 2D Fermi-liquid quasiparticles. Modification of one of the 2D layers into a superlattice of one-dimensional (1D) quantum wires has been employed to measure vertical tunneling between 1D and 2D electron systems. Constraints on tunneling imposed by the requirement of simultaneous conservation of energy and momentum can be tuned by the transport voltage and external magnetic fields. In certain situations, this makes it possible to observe features of the momentum-resolved single-electron spectral function directly in tunneling transport.

The method of cleaved-edge overgrowth (CEO) makes it possible to create long and clean quantum wires in GaAs/GaAlAs heterostructures. Using the same technique, systems of two parallel quantum wires with a high and extremely clean tunneling barrier between them have been fabricated in double-layer structures. This opens up new possibilities for studying the peculiar dynamics of electrons in interacting 1D systems using 1D–to–1D tunneling. In particular, both the phase-coherence length and the elastic mean free path $`l_{\text{el}}`$ for electrons in these quantum wires usually exceed the wire length. This motivates the present work where we analyze mesoscopic effects in 1D–to–1D transport. In related contexts, phase-coherent transport in double-wire systems was discussed in terms of device applications. For example, a system of two parallel, identical quantum wires coupled within a spatial region of length $`L<l_{\text{el}}`$ via an adjustable tunneling barrier was proposed as a possible realization of a current switch. Wave packets of electrons injected into one of the wires will be coherently transferred to the other one and back with frequency $`2|t|/\mathrm{}`$. (Here we denoted the tunnel-splitting of energy levels in the coupling region by $`2|t|`$.) In steady state, this results in a coherent charge oscillation in real space with wave length $`L_t^{(0)}=\pi \mathrm{}v_\text{F}/|t|`$. Modulation of $`|t|`$ controls the signal at the output of the injecting wire. Ideally, it is maximal (minimal) when the ratio of $`L`$ and $`L_t^{(0)}`$ is (half-)integer. In reality, output characteristics depend sensitively on details of the tunneling barrier. Assuming the feasibility to engineer barrier design, coupled quantum wires were suggested as realizations of quantum logical gates.

In this paper, we consider phase-coherent transport in a system of parallel quantum wires coupled via a finite tunneling barrier. See Fig. 1. Charging effects in the wires caused by applied voltages influence tunneling in an important way because they determine the degree to which 1D subbands are shifted or filled. The basic physics of this interplay is discussed in the following section. Our microscopic model for the double-wire system is introduced in the first part of Sec. III. Apart from capacitance effects, interactions are neglected within our approach, which is therefore valid only for voltages and in-plane magnetic fields probing the 1D electron systems beyond the cut-off for Luttinger-liquid behavior. Results from lowest-order perturbation theory are compared with the exact solution using scattering theory. We calculate linear and differential conductances for 1D–to–1D tunneling transport, and discuss their features in Sec. IV.

## II Effect of an applied voltage

In the typical tunneling experiment, a voltage drop $`V`$ across the barrier drives a current. Microscopically, it is often assumed that the voltage shifts quasiparticle bands in the two subsystems by $`\pm eV/2`$, respectively, as compared to the equilibrium situation where no net current flows. The external circuit is supposed to prevent charging of the subsystems. In general, however, the applied voltage will shift the bands as well as partly fill them. As the I–V curve for 1D–to–1D tunneling depends sensitively on the scenario of band filling vs. band shifting, we discuss this issue here in some detail.

At zero temperature, the free energy (per length) of a 1D system is given by its total energy (per length) $`E_{\text{tot}}`$ which is a functional of particle density $`n`$. In a clean quantum wire, $`n=n_0`$ will be constant. Before applying a voltage, the system is assumed to be charge-neutral, i.e., the uniform electronic charge density $`en_0`$ is compensated by positive background and image charges. It is useful to divide $`E_{\text{tot}}`$ into two parts; $`E_{\text{tot}}=E_{\text{int}}+E_{\text{Coul}}`$. All Coulombic terms (including the Hartree energy of electrons in the wire) are collected in $`E_{\text{Coul}}`$, and $`E_{\text{int}}`$ is the internal energy of the quantum wire comprising kinetic and exchange-correlation contributions. For our purposes, we adopt the simple model with $`E_{\text{Coul}}=(e\mathrm{\Delta }n)^2/2\stackrel{~}{C}`$ where $`\mathrm{\Delta }n`$ is the deviation from the neutralized charge density $`n_0`$, and $`\stackrel{~}{C}`$ denotes the electrostatic capacitance per unit length of the wire. The applied voltage is assumed to lead to a uniform shift $`eV`$ with respect to the equilibrium chemical potential $`\mu _0=E_{\text{tot}}/n|_{n_0}`$ of the wire. The induced change $`\mathrm{\Delta }n`$ in the total density has to be calculated from

$$\mu _0+eV=\frac{e^2}{\stackrel{~}{C}}\mathrm{\Delta }n+\frac{E_{\text{int}}}{n}|_{n=n_0+\mathrm{\Delta }n}.$$
(1)
In the limit of small voltages ($`|eV|\mu _0`$) where linear-response theory is valid, we can use

$$\frac{E_{\text{int}}}{n}|_{n=n_0+\mathrm{\Delta }n}\mu _0\frac{\mathrm{\Delta }n}{D_0}.$$
(2)
Here, $`D_0`$ is the thermodynamic density of states (DOS) defined by $`D_0=n/\mu |_{\mu =\mu _0}=\left(^2E_{\text{int}}/n^2\right)_{n=n_0}^1`$. In the linear-response limit, it is then possible to express $`\mathrm{\Delta }n`$ explicitly in terms of the external voltage,

$$\mathrm{\Delta }n=eV\frac{D_0}{1+\zeta },$$
(3)
where the parameter $`\zeta `$ measures the relative importance of Coulombic and density-of-states effects:

$$\zeta =\frac{e^2}{\stackrel{~}{C}}D_0.$$
(4)
We see that, for $`\zeta 1`$, a voltage will simply fill quasiparticle bands without shifting them (band-filling limit). In particular, this applies in the absence of electron-electron interactions. In the opposite case $`\zeta 1`$, an applied voltage shifts the bands (band-shifting limit). This situation is analogous to that of a bulk metal or a single-electron transistor. To get an idea of the situation realized in our system of interest, we estimate the capacitance of a quantum wire by $`\stackrel{~}{C}=2\pi ϵ/\mathrm{ln}(R/r)`$, with $`R`$ being the distance to surrounding metal gates, and $`r`$ denoting the characteristic transverse dimension of the wire. Typical values are $`10^{10}`$ F/m. With Fermi energies of quantum wires ranging between $`1\mathrm{}10`$ meV, we obtain $`\zeta 1\mathrm{}10`$. Hence, typical quantum wires are in the intermediate regime where both band-filling and band-shifting occurs at the same time. This case is illustrated in Fig. 2. It is important to keep in mind, however, that Eqs. (3) and (4) are only valid when $`\mathrm{\Delta }nn_0`$. In experiment, voltages comparable to and larger than $`\mu _0`$ are applied to probe the full single-particle dispersion relation. Then, for a quantitative comparison between theory and experiment, $`\mathrm{\Delta }n`$ has to be found from Eq. (1). For example, a wire whose density was initially large enough for it to be in the band-filling limit crosses over to the band-shifting limit when it is depleted.

At this point, it is useful to make contact with results obtained for Tomonaga-Luttinger (TL) models of interacting 1D electron systems. Unlike their higher-dimensional counterparts, 1D metals cannot be described within the Fermi-liquid paradigm. Instead, their low-energy properties are represented by effective TL models, and the phenomenology of a Luttinger liquid (LL) applies. Instead of Landau parameters, it is the velocities of certain collective and zero-mode excitations that determine all physical quantities of a LL. In particular, the ratio $`r_\text{N}=v_\text{N}/v_\text{F}`$ of the velocity $`v_\text{N}`$ of the charged zero mode and the bare Fermi velocity enters the expression for the electrostatic capacitance per unit length of a Luttinger liquid: $`\stackrel{~}{C}_{\text{LL}}=e^2D_0/(r_\text{N}1)`$. Here, $`D_0=1/\pi \mathrm{}v_\text{F}`$ is the 1D DOS. Using Eq. (4), we find $`\zeta _{\text{LL}}=r_\text{N}1`$. The non-interacting case where $`r_\text{N}=1`$ corresponds to the band-filling limit, whereas strong Coulomb interactions ($`r_\text{N}\mathrm{}`$) recover the band-shifting limit. We would like to remark that $`r_\text{N}`$ constitutes an independent parameter in the low-energy theory of any given real quasi-1D system. It is unrelated, except in certain special cases, to the famous interaction parameter $`K_\rho `$ that enters power-law expressions for electronic correlation functions. From now on, we consider the model where $`r_\text{N}>1`$ but $`K_\rho =1`$. This approximation is valid to describe current experiments where the wires are probably not long enough for the power-law characteristics of a LL to be observable. Even for infinitely long wires, however, our results apply at energies and wave vectors far enough from the Fermi points where the single-particle spectral function recovers Fermi-liquid-like characteristics.

## III Model and Formalism

We consider two quantum wires of infinite length, labeled U(pper) and L(ower), that are parallel to the $`x`$ direction and located, in the $`yz`$ plane, at $`y=0`$ and $`z=z_{\text{U,L}}`$. The potential barrier between them is assumed to be finite and uniform in the region $`|x|L/2`$ and infinite otherwise. Within the standard notation of second quantization, the Hamiltonian for our system is given by

$`H`$ $`=`$ $`{\displaystyle \underset{\alpha =\text{U,L}}{}}H_\alpha +H_{\text{tun}},`$ (6)
$`H_\alpha `$ $`=`$ $`{\displaystyle \frac{dk}{2\pi }ϵ_\alpha (k)c_{k\alpha }^{}c_{k\alpha }},`$ (7)
$`H_{\text{tun}}`$ $`=`$ $`{\displaystyle _{L/2}^{L/2}}dx\{t\psi _\text{U}^{}(x)\psi _\text{L}(x)+\mathrm{H}.\mathrm{c}.\}`$ (8)
$`=`$ $`{\displaystyle \frac{dk}{2\pi }\frac{dp}{2\pi }\left\{t_{k,p}c_{k\text{U}}^{}c_{p\text{L}}+\text{H.c.}\right\}}.`$ (9)
Here, $`ϵ_{\text{U(L)}}(k)`$ is the electronic dispersion relation in the upper (lower) wire. Modulo an unimportant phase factor, the tunneling matrix element is given by

$$t_{k,p}=2|t|\frac{\mathrm{sin}\left[(pk)L/2\right]}{pk}=|t|\sqrt{2\pi L\delta _L(pk)}.$$
(10)
The second equality in Eq. (10) constitutes the definition of $`\delta _L(pk)`$ which is a finite-size realization of Dirac’s $`\delta `$-function. Tunneling occurs mainly between states with momenta satisfying $`|pk|<2\pi /L`$. Perfect momentum conservation holds only in the limit $`L\mathrm{}`$.

We consider the case where a single (the lowest) 1D subband in each wire is occupied and assume a parabolic subband dispersion. The effect of a magnetic field $`\stackrel{}{B}=B\widehat{y}`$ applied perpendicularly to the plane of the two wires can be included by a shift of kinetic with respect to canonical momentum. Then, dispersion relations read ($`\alpha =\text{U,L}`$)

$$ϵ_\alpha (k)=\frac{\mathrm{}^2}{2m}\left(k\frac{eB}{\mathrm{}}z_\alpha \right)^2+E_{0\alpha }+\nu _\alpha eV_\alpha .$$
(11)
Here, $`m`$ is the effective electron mass in the semiconductor host medium, and $`E_{0\alpha }`$ denotes the energy at the bottom of the respective wire’s lowest 1D subband. The term $`\nu _{\text{U(L)}}eV_{\text{U(L)}}`$ takes into account the shifting of the band in the upper (lower) wire due to an applied voltage. See Fig. 2. For simplicity, we neglect effects due to the mutual capacitance of the two wires, which can be included straightforwardly. In general, the values of $`\nu _\alpha `$ will depend on voltage. Furthermore, except in the limit of weak tunneling where $`|t||eV_\alpha |`$, they have to be determined from a self-consistent treatment of charging effects resulting from tunneling and electrostatics. While this is, in principle, straightforward to implement, we choose to focus here on the weak-tunneling limit which is more relevant for current experiment. In the linear-response regime, we have $`\nu _\alpha =\zeta _\alpha /(1+\zeta _\alpha )`$ with $`\zeta _\alpha `$ defined for each wire in analogy to Eq. (4).

Absolute values of energy and the $`z`$ coordinate are irrelevant; results depend only on the difference of subband energies, $`\mathrm{\Delta }E_0=E_{\text{0L}}E_{\text{0U}}`$, and the wire separation, $`d=z_\text{U}z_\text{L}`$. For simplicity, we choose $`E_{\text{0U}}=0`$ and $`z_\text{U}=0`$ in the following. Also, to avoid cluttering the notation, we have suppressed spin quantum numbers. In typical CEO structures, the effect of Zeeman splitting is negligible for the range of magnetic fields to be considered below. Hence, electron spin leads only to factors of 2 which we include in our final formulae for tunneling current and conductances.

### A Perturbation theory: Lowest order in tunneling

A standard procedure for calculating the tunneling current is to perform perturbation theory in $`H_{\text{tun}}`$. To leading order, the current flowing from the upper to the lower wire is

$`I`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}}}{\displaystyle \frac{dk}{2\pi }\frac{dp}{2\pi }|t_{k,p}|^2_{\mathrm{}}^{\mathrm{}}\frac{dϵ}{2\pi }A_\text{U}(ϵ,k)A_\text{L}(ϵ,p)}`$ (13)
$`\times \left[f_\text{U}(ϵ)f_\text{L}(ϵ)\right],`$
with Fermi functions $`f_\alpha (ϵ)=1/[1+\mathrm{exp}\{(ϵ\mu _0eV_\alpha )/k_\text{B}T\}]`$. The single-particle spectral functions for the wires are given, within the model specified above, by $`A_\alpha (ϵ,k)=2\pi \delta [ϵϵ_\alpha (k)]`$. In the linear-response limit ($`|eV_\alpha |\mu _0`$), we find the tunneling conductance

$$G=\frac{2e^2}{\mathrm{}^3}\frac{|t|^2L}{v_{\text{FU}}v_{\text{FL}}}\underset{\gamma ,\gamma ^{}=\pm 1}{}\delta _L\left(\frac{\pi }{2}[\gamma n_\text{U}\gamma ^{}n_\text{L}]p_\text{B}\right).$$
(14)
Here, $`v_{\text{F}\alpha }`$ and $`n_{\text{F}\alpha }`$ denote the Fermi velocity and electron density of the respective wire at the equilibrium chemical potential $`\mu _0`$. The peak-shape function $`\delta _L`$ has been defined in Eq. (10), and the relative shift of the 1D Fermi seas due to the applied magnetic field is $`p_\text{B}=eBd/\mathrm{}`$.

In analogy to 2D–to–2D tunneling, resonances appear in the tunneling conductance $`G`$ as function of magnetic field whenever parts of the shifted Fermi surfaces of the two wires overlap. As the 1D Fermi surface consists of just two points, the shape of these resonances is that of a smeared delta function of width $`2\pi /L`$. The peak value of the tunneling conductance can be written as

$$G_{\text{pk}}=n^{}\frac{2e^2}{h}\left(\frac{\pi L}{L_t}\right)^2,$$
(15)
where $`L_t=\pi \mathrm{}\sqrt{v_{\text{FU}}v_{\text{FL}}}/|t|`$ is an effective length scale introduced by tunneling, and $`n^{}=1`$ or $`2`$ depending on the number of overlapping Fermi points at peak condition.

In the dirty limit where the length $`L`$ of the tunneling barrier is larger than the mean free path $`l_{\text{el}}`$, Eq. (14) is still valid but the peak width is now given by $`2\pi /l_{\text{el}}`$. Also, the factor $`(\pi L/L_t)^2`$ in Eq. (15) has to be replaced by $`2\pi ^2Ll_{\text{el}}/L_t^2`$. In both the ballistic and diffusive cases, taking the limit of $`L\mathrm{}`$ is unphysical: the conductance through the barrier cannot exceed $`2e^2/h`$ per channel. Hence, the actual small parameter enabling perturbative treatment of tunneling is $`L/L_t`$. Smallness of $`L/L_t`$ means that the time between tunneling events has to be larger than the time it takes electrons to traverse the region where the potential barrier between the wires is finite. Only then it will be possible to neglect higher-order effects due to electrons tunneling coherently back and forth between the wires. Using the exact solution developed in the next subsection, we will find indeed that the perturbative result displayed in Eq. (14) is valid only as long as $`LL_t`$.

### B Exact solution using scattering theory

As the model defined in Eqs. (III) describes two systems of noninteracting fermionic quasiparticles that are coupled via tunneling in a finite region of space, we can use scattering theory for calculating transport. To make this explicit, we rewrite the Hamiltonian of our system in first-quantized notation and real-space representation. It is a $`2\times 2`$ matrix \[because wave functions are two-component spinors $`(\psi _\text{U},\psi _\text{L})^\mathrm{T}`$\]:

$$H=\left(\begin{array}{cc}ϵ_\text{U}(i_x)\hfill & \hfill t(x)\\ t(x)\hfill & \hfill ϵ_\text{L}(i_x)\end{array}\right).$$
(16)
The tunneling matrix element is piecewise constant: $`t(x)=|t|`$ for $`|x|L/2`$ and $`t(x)=0`$ otherwise. Hence, regions with $`|x|>L/2`$ where the wires are independent act as leads where scattering states can be defined. We attach labels 1 through 4 to these leads as shown in Fig. 1. The region $`|x|L/2`$ where tunneling occurs acts as an effective scatterer. The current flowing through the tunnel barrier is then given by

$$I=\frac{2e}{h}\underset{\genfrac{}{}{0pt}{}{m=1,2}{n=3,4}}{}𝑑ϵ\left|T_{m,n}(ϵ)\right|^2\left[f_\text{U}(ϵ)f_\text{L}(ϵ)\right],$$
(17)
where $`T_{m,n}(ϵ)`$ denotes the transmission coefficient for electrons with energy $`ϵ`$ that originate in lead $`m`$ and are scattered into lead $`n`$. We calculate the transmission coefficients by matching scattering states in the leads to the appropriate eigenstates of the Hamiltonian (16) in the region $`|x|L/2`$. As this is a straightforward exercise, and results for the most general case are lengthy, we omit explicit formulae here. Due to the difference of Fermi functions in Eq. (17), only transmission coefficients $`T_{m,n}(ϵ)`$ at energies within the voltage window, i.e., with $`ϵ\mu _0[eV_\text{L},eV_\text{U}]`$, contribute. In the limit of small applied voltage, Eq. (17) yields the linear conductance

$$G=\frac{2e^2}{h}\underset{\genfrac{}{}{0pt}{}{m=1,2}{n=3,4}}{}\left|T_{m,n}(\mu _0)\right|^2.$$
(18)
As expected, $`G`$ obtained from Eq. (18) deviates from the perturbative result \[Eq. (14)\] for long enough barrier length $`L`$, see Fig 3. The oscillatory dependence of $`G`$ on $`L`$ can be tuned by the applied magnetic field, as seen in the inset of Fig. 3. When the effective tunnel splitting in the coupling region is much smaller than the Fermi energy of the quantum wires, the following approximate formula for the linear conductance can be derived (see the Appendix for details):

$$G\frac{2e^2}{h}\underset{\gamma ,\gamma ^{}=\pm 1}{}\frac{\mathrm{sin}^2\left(\pi \sqrt{\left(L/L_t\right)^2+\left(L/L_{\gamma \gamma ^{}}\right)^2}\right)}{1+(L_t/L_{\gamma \gamma ^{}})^2}.$$
(19)
Here, new length scales $`L_{\gamma \gamma ^{}}`$ appear that measure the mismatch of canonical Fermi momentum for pairs of Fermi points from the upper (right-mover $`\gamma =+1`$, left-mover $`\gamma =1`$) and lower ($`\gamma ^{}`$ analogous) wires:

$$L_{\gamma \gamma ^{}}=\frac{2\pi }{\frac{\pi }{2}\left[\gamma n_\text{U}\gamma ^{}n_\text{L}\right]p_\text{B}}.$$
(20)
Exact calculation of $`G`$ in the appropriate limit confirms the validity of Eq. (19); see Fig. 4.

## IV Results and Discussion

Different regimes in the behavior of the linear 1D–to–1D tunneling conductance are distinguished by the interplay of the relevant length scales encountered above. These are the length $`L`$ of the tunnel barrier, $`L_t`$ which is a measure of the strength of tunneling, and lengths $`L_{\gamma \gamma ^{}}`$ which are defined for any pairing of a Fermi point from the upper wire with one from the lower wire. For the following discussion, we consider only the largest $`L_{\gamma \gamma ^{}}`$ of all possible. Comparing $`L_t`$ with the other lengths, a weak-tunneling regime ($`L_t>\mathrm{max}\{L_{\gamma \gamma ^{}},L\}`$) can be distinguished from a strong-tunneling regime ($`L_t<\mathrm{max}\{L_{\gamma \gamma ^{}},L\}`$). Furthermore, we call the system in resonance when the Fermi points of a corresponding $`L_{\gamma \gamma ^{}}`$ are close to each other on the scale of $`2\pi /L`$, i.e., when $`L<L_{\gamma \gamma ^{}}`$. Conversely, the off-resonance limit is reached for $`L>L_{\gamma \gamma ^{}}`$.

In the strong-tunneling regime, the linear conductance oscillates as a function of $`L`$ with wave length $`L_t`$ and maximum amplitude $`2e^2/h`$ ($`4e^2/h`$ for identical wires). Previously, when the feasibility of using the double-wire system as a directional coupler was discussed, the in-resonance limit was considered only. Control of directional-coupler operation is then possible only by varying $`L_t`$, i.e., essentially only by adjusting the barrier height. Here we find that, in the off-resonance limit, the device is tunable, in addition to varying $`t`$, by an applied magnetic field or, equivalently, by adjusting the density mismatch in the two wires. This is seen already in the inset of Fig. 3 and, more clearly, in Fig. 4. It is, therefore, possible to adjust the effective length scale for coherent electron transfer between the wires by applying a magnetic field. When $`L`$ is equal to several times $`L_t`$, the difference between the effective transfer length in a magnetic field and $`L_t`$ leads to an accumulated phase shift over many oscillations that can reach $`\pi /2`$ without concomitant loss in amplitude. In particular for CEO structures, operating the system in off-resonance mode provides a convenient alternative to any (hardly feasible) adjustment of the high tunnel barrier.

The weak-tunneling limit is well-suited for spectroscopic application of 1D–to–1D tunneling. Sharp peaks are exhibited by both the linear and differential conductances for in-resonance conditions. Measured on the scale of resonance peaks, the off-resonance conductance is orders of magnitude smaller. We have calculated the differential conductance as a function of magnetic field and voltage whose resonance condition corresponds to a Fermi point of one wire coinciding with a point on the dispersion curve of the other wire. The exact location of these coincidences in the $`V`$$`B`$ plane depends sensitively on charging effects in the quantum wires. In the following, we focus on the limit of weakly coupled wires where the self-consistent charge profile is not importantly affected by tunneling. Furthermore, we consider the situation with symmetric bias $`V_\text{U}=V_\text{L}=V/2`$.

Figures 5 and 6, respectively, show logarithmic gray-scale plots of the absolute value of the differential tunneling conductance for the ideal band-filling ($`\stackrel{~}{C}=\mathrm{}`$) and band-shifting ($`\stackrel{~}{C}=0`$) cases. Bright lines are formed by points in the $`V`$$`B`$ plane where the above-mentioned resonance condition is fulfilled. Due to the finite length of the tunnel barrier, more maxima appear with peak values being orders of magnitude smaller than at the resonance peaks. In the band-shifting case, resonance lines are direct images of parts of the wires’ electronic dispersion curves. In particular, the extension of the leaf-shaped structure in the positive and negative voltage direction provides a direct measure of the respective wire’s Fermi energy. This is not the case for the band-filling limit, which is characterized by a resonance line running close to the voltage axis when $`E_0\mu _0`$. Its leaf structure is symmetric under voltage reversal, with extension in (positive or negative) voltage direction given by $`2(\mu _0|\mathrm{\Delta }E_0|)`$. In real systems, the capacitance $`\stackrel{~}{C}`$ is finite, and an intermediate picture will be obtained for the differential tunneling conductance. An example is shown in Fig. 7, where $`\stackrel{~}{C}`$ has a value such that $`\zeta _\text{U}=0.5`$ at zero voltage. Depletion of one of the wires for increasing voltage leads to a cross-over to the band-shifting situation. As a result, the ideal band-shifting limit is not easily distinguished from the intermediate case. Quantitative comparison of the measured resonance structure in the differential conductance with results expected from an independent measurement of Fermi energies, electron densities etc. will have to include the effect of the finite $`\stackrel{~}{C}`$. Conversely, for known Fermi-sea parameters, the value of $`\stackrel{~}{C}`$ can be extracted by fitting the measured resonance pattern of the differential conductance for 1D–to–1D tunneling.

## V Conclusions

Motivated by recent experiment, we have investigated linear and differential conductances for 1D–to–1D transport. Our results show that effects due to phase coherence and charging of the wires are important for realistic double-wire structures. Regimes of weak and strong tunneling, as well as in and out of resonance, are distinguished and their key features discussed. We point out new possibilities for device application of 1D–to–1D tunneling and its use for electron-dispersion spectroscopy.

###### Acknowledgements.

Useful discussions with M. Büttiker, G. Schön, and A. Yacoby are gratefully acknowledged. This work was supported by DFG within Sonderforschungsbereich 195, Graduiertenkolleg ”Kollektive Phänomene im Festkörper” (D.B.), and the Emmy-Noether-Programm (A.R.). U.Z. thanks the Braun Submicron Center at the Weizmann Institute, Israel, for hospitality during a visit sponsored by the EU LSF programme.

## Calculation of transmission coefficients

Equations (17) and (18) express tunneling current and linear conductance in terms of transmission coefficients $`T_{m,n}(ϵ)`$. These transmission coefficients can be obtained exactly by matching eigenstates of Hamiltonian $`H`$ \[given in Eq. (16)\] with eigenvalue $`ϵ`$ in the coupling region ($`|x|L/2`$) to appropriate eigenstates in the leads. For example, to calculate $`T_{1,n}`$, we use the ansatz

$`\mathrm{\Psi }_ϵ(x)|_{x<\frac{L}{2}}`$ $`=`$ $`\left(\begin{array}{c}1\\ 0\end{array}\right)e^{ik_\text{U}^{(+)}x}`$ (24)
$`+`$ $`\left(\begin{array}{c}t_{11}\\ 0\end{array}\right)e^{ik_\text{U}^{()}x}+\left(\begin{array}{c}0\\ t_{13}\end{array}\right)e^{ik_\text{L}^{()}x},`$ (29)
$`\mathrm{\Psi }_ϵ(x)|_{x>\frac{L}{2}}`$ $`=`$ $`\left(\begin{array}{c}t_{12}\\ 0\end{array}\right)e^{ik_\text{U}^{(+)}x}+\left(\begin{array}{c}0\\ t_{14}\end{array}\right)e^{ik_\text{L}^{(+)}x},`$ (34)
$`\mathrm{\Psi }_ϵ(x)|_{|x|<\frac{L}{2}}`$ $`=`$ $`{\displaystyle \underset{\alpha =\text{a,b}}{}}\{d_+^{(\alpha )}\left(\begin{array}{c}u_+^{(\alpha )}\\ v_+^{(\alpha )}\end{array}\right)e^{ik_+^{(\alpha )}x}`$ (40)
$`+d_{}^{(\alpha )}\left(\begin{array}{c}u_{}^{(\alpha )}\\ v_{}^{(\alpha )}\end{array}\right)e^{ik_{}^{(\alpha )}x}\}.`$
Here, wave vectors $`k_{\text{U(L)}}^{(\pm )}`$ are solutions of $`ϵϵ_{\text{U(L)}}(k_{\text{U(L)}}^{(\pm )})`$ with positive ($`+`$) and negative ($``$) group velocity $`v_{\text{U(L)}}=_\mathrm{}kϵ_{\text{U(L)}}(k)`$, respectively. With the energy dispersion in the coupling region given by $`ϵ_\pm (k)=\frac{1}{2}\left[ϵ_\text{U}(k)+ϵ_\text{L}(k)\right]\pm |t|\sqrt{1+r^2}`$, we have $`ϵϵ_\pm \left(k_\pm ^{(\text{a,b})}\right)`$. The function $`r=[ϵ_\text{U}(k)ϵ_\text{L}(k)]/2|t|`$ measures the mismatch in the dispersions of the two wires and determines the amplitudes $`u_\pm =\sqrt{(1\pm r/\sqrt{1+r^2})/2}`$, $`v_\pm =\pm \sqrt{(1r/\sqrt{1+r^2})/2}`$. Requiring continuity of the wave function and current conservation at the locations $`x=\pm L/2`$ yields a system of linear equations from which the coefficients $`t_{mn}`$ and $`d_\pm ^{(\alpha )}`$ are found. Transmission coefficients entering Eqs. (17) and (18) are then given by $`T_{m,n}(ϵ)=t_{mn}\sqrt{|v_\text{L}/v_\text{U}|}`$.

The linear conductance is given in terms of transmission coefficients at the equilibrium chemical potential $`\mu _0`$, as expressed in Eq. (18). To derive Eq. (19), we consider, e.g., $`T_{1,4}`$. When the tunnel splitting $`2|t|`$ and Fermi-energy mismatch of the two wires is much smaller than their Fermi energies, only amplitudes for right-moving partial waves in ansatz (Calculation of transmission coefficients) are significantly different from zero. Neglecting left-moving partial waves and the small density mismatch, the matching procedure yields $`t_{14}=i\mathrm{sin}(k_tL)/\sqrt{1+r^2}`$ with $`k_t=t\sqrt{1+r^2}/\mathrm{}v_\text{F}`$. In the regime considered, we have $`rL_t/L_{11}`$ with $`L_{11}`$ defined in Eq. (20). Similar calculations for other transmission coefficients finally yield Eq. (19).


format: text

text:
# Lees-Edwards boundary conditions for Lattice Boltzmann

## 1 Introduction

In 1972 Lees and Edwards published a seminal paper  describing an extension of the periodic boundary conditions for Molecular Dynamics that allows the simulation of a bulk system in a simple shear flow. The uniform shear flow steady state reached is computationally convenient because it reduces finite size effects. If moving solid walls are instead introduced to produce a shear flow, the spatial inhomogeneities induced close to the wall limit the spatial region within the simulated system that can be used to study the bulk behavior of a system under shear. LEbc ensure that the system is spatially homogeneous, so that bulk-like behavior is recovered on smaller simulated systems. For lattice Boltzmann methods there is the additional problem that the maximum shear rate is limited by the existence of a maximum shear rate of order $`1/L_y`$ with $`L_y`$ the transverse system size. This is required to ensure that all fluid velocities remain small compared to unity in lattice units.

The original LEbc were implemented for a system of particles, with positions $`(r_x,r_y,r_z)`$ and velocities $`(v_x,v_y,v_z)`$ within a box of dimension $`(L_x,L_y,L_z)`$. Lees and Edwards showed that by implementing a new kind of boundary condition in a MD simulation a state with a uniform shear velocity profile $`𝐮=\dot{\gamma }y𝐞_x`$ (where $`𝐞_x`$ is the unity vector in x-direction) could be achieved. The method consists of a simple recipe that imposes periodic boundary conditions on particles leaving the simulation box in the directions perpendicular to the velocity gradient $`𝐞_y`$. Particles leaving the box in the $`𝐞_y`$ direction ($`r_y>L_y`$) will be reinserted with their $`x`$-velocity increased by the shear velocity $`\mathrm{\Delta }𝐮=\dot{\gamma }L_y𝐞_x`$ added at a position given by their periodic image displaced by the time-dependent offset $`d_x=t\mathrm{\Delta }u_x`$, $`t`$ being the time elapsed from an appropriate origin of times. If the particles leave the box in the $`𝐞_y`$ direction ($`r_y<0`$) they reappear at the position of their periodic image displaced by $`d_x`$ and have the shear velocity $`\mathrm{\Delta }𝐮`$ subtracted from their $`x`$-velocity.

We can describe the boundary conditions by providing the new positions $`r^{}`$ and velocities $`v^{}`$ of the particles after the particles have been moved by the dynamics of the algorithm

$`r_x^{}`$ $`=`$ $`\{\begin{array}{cc}(r_xmodL_x)+d_x\hfill & \hfill r_zL_z\\ (r_xmodL_x)\hfill & \hfill 0r_z<L_z\\ (r_xmodL_x)d_x\hfill & \hfill r_z<0\end{array}`$ (4)
$`r_y^{}`$ $`=`$ $`(r_ymodL_y)`$ (5)
$`r_z^{}`$ $`=`$ $`(r_zmodL_z)`$
$`v_x^{}`$ $`=`$ $`\{\begin{array}{cc}v_x+s_x\hfill & \hfill r_zL_z\\ v_x\hfill & \hfill 0r_z<L_z\\ v_xs_x\hfill & \hfill r_z<0\end{array}`$ (9)
$`v_y^{}`$ $`=`$ $`v_y`$ (10)
$`v_z^{}`$ $`=`$ $`v_z`$
This is the best known way of implementing boundary conditions that allow the simulation of a bulk material under shear.

This method works well for all particle-based simulation methods. The aim of this paper is to adapt this method for lattice Boltzmann models, where the fundamental quantity is a distribution function, rather than a set of particles, which is defined on a unit lattice. In this case, the generalization of LEbc is not straightforward since the displacement $`d_x`$ will in general not correspond to a full integer multiple of the lattice-spacing and it is not clear how to implement the required velocity shift in a method with a fixed set of velocities. In an earlier paper  one of the authors introduced an interpolation method to deal with the first problem and a partial Galilean transformation which was used to add momentum. This approach required the use of an extended equilibrium distribution function. Hence, it was not applicable to standard lattice-Boltzmann models.

In this paper we develop an alternative approach to impose Lees-Edwards boundary conditions in standard lattice-Boltzmann models. It is based on applying a Galilean-transformation and imposing the appropriate change of momentum across the moving plane. Additionally, we must displace the densities appropriately with an interpolation scheme. The derived method of Galilean transformations is quite general and we we can re-derive the correct boundary conditions for a moving solid wall in contact with the fluid which were originally introduced by A. Ladd .

After describing the lattice-Boltzmann model, we will discuss in detail the implementation of Lees-Edwards boundary conditions, and we will subsequently validate the method by examining implementations for binary fluid mixtures in two dimensions.

## 2 Lattice Boltzmann model

The lattice Boltzmann method can be viewed as a simple discretization of the Boltzmann equation where space is represented by a lattice and the one-particle distribution functions (which we will refer to as ‘densities’) $`f_i(𝐱,t)`$ exist at each node $`𝐱`$ of the underlying lattice, for only a small number of velocities, $`\{𝐯_i\}`$, that (with unit time step) usually connect a small subset of lattice points. These densities evolve following a discrete-time dynamics. In each iteration, the densities are moved along their corresponding velocity vectors during the streaming step to the appropriate neighboring site, while in the subsequent collision step the set of densities at a given node are relaxed toward equilibrium. This relaxation step is performed ensuring mass and momentum conservation at each node. The energy is usually not conserved (although energy conserving schemes also exist ) but it is possible to introduce a nominal temperature, that can be viewed as the effect of an underlying thermostat. This temperature is introduced as a parameter of the equilibrium distribution function (see below). The evolution of the densities in each time step can be written as

$$f_i(𝐱+𝐯_i,t+1)=f_i(𝐱,t)+\frac{1}{\tau }\left(f_i^0(\rho (𝐱,t),𝐮(𝐱,t))f_i(𝐱,t)\right)$$
(11)
where $`\tau `$ is a relaxation parameter, which in the hydrodynamic limit of this equation is related to the viscosity of the fluid. Eq.(11) can be viewed as the discretized version of BGK approximation of the Boltzmann equation. The linearized collision operator can be generalized to allow for several relaxation parameters (and hence the existence of different viscosities for the fluid). In numerical applications the above evolution equation is usually split into two subroutines, a collision and a streaming step. The collision step is local

$$\widehat{f}_i(𝐱,t)=f_i(𝐱,t)+\frac{1}{\tau }\left(f_i^0(\rho (𝐱,t),𝐮(𝐱,t))f_i(𝐱,t)\right)$$
(12)
and the streaming step moves densities along their associated velocity vectors

$$f_i(𝐱+𝐯_i,t+1)=\widehat{f}_i(x,t).$$
(13)
In the above equations, $`f_i^0`$ is the equilibrium distribution function that determines the thermodynamics of the fluid, while $`\rho `$ and $`\rho 𝐮`$ are the local mass and momentum densities. These are the macroscopic quantities of interest, and can be obtained from the densities $`f_i`$ as appropriate moments

$$\rho =\underset{i}{}f_i\rho 𝐮=\underset{i}{}f_i𝐯_i$$
(14)
The choice of the equilibrium densities is one of the key ingredients of the model. They are chosen to recover the equilibrium density and momentum

$$\rho =\underset{i}{}f_i^0\rho 𝐮=\underset{i}{}f_i^0𝐯_i$$
(15)
and the appropriate pressure tensor in equilibrium

$$𝒫_{\alpha \beta }=\underset{i}{}f_i^0(v_{i\alpha }u_\alpha )(v_{i\beta }u_\beta )$$
(16)
where $`𝒫`$ is the sum of the thermodynamic pressure tensor $`P`$ and an additional term that is required for Galilean invariance of the system :

$$𝒫_{\alpha \beta }=P_{\alpha \beta }+\frac{\tau 0.5}{3}(u_\alpha _\beta \rho +u_\beta _\alpha \rho )$$
(17)
Using the thermodynamic pressure tensor $`P`$ is one of the standard ways to impose the desired thermodynamic behavior of the lattice-Boltzmann fluid. In the simplest case, the pressure tensor for the ideal gas is $`P_{\alpha \beta }=\rho \frac{1}{3}\delta _{\alpha \beta }`$. We need an additional constraint for the third moment

$$\underset{i}{}f_i^0v_{i\alpha }v_{i\beta }v_{i\gamma }=\frac{\rho }{3}(u_\alpha \delta _{\beta \gamma }+u_\beta \delta _{\alpha \gamma }+u_\gamma \delta _{\alpha \beta })$$
(18)
to ensure that the anisotropy of the underlying lattice does not affect the effective behavior of the model in the hydrodynamic limit. It is possible to fulfill these constraints by introducing equilibrium densities that depend quadratically on the lattice velocities, $`\{𝐯_i\}`$. Depending on the geometry of the lattice, there is still some freedom in the choice of the equilibrium densities.

A Taylor expansion to second order in the derivatives of the zeroth and first moments in $`𝐯`$ of the lattice Boltzmann equation (11) gives the macroscopic equations which express mass and momentum conservation. From (11) we have

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n!}(_t+𝐯_i)^nf_i=\frac{1}{\tau }(f_i^0f_i)$$
(19)
and we can approximate

$$f_i=f_i^0+\tau (_tf_i^0+v_{i\alpha }_\alpha f_i^0)+O(^2).$$
(20)
The zeroth velocity moment of (11) up to $`O(^2)`$ gives the continuity equation

$$_t\rho +_\alpha (\rho u_\alpha )=0$$
(21)
and the first velocity moment gives the Navier Stokes equations (with $`_\alpha \mathrm{𝐮𝐮𝐮}=O(^3)`$)

$$\rho _tu_\alpha +\rho u_\beta _\beta u_\alpha =_\beta P_{\alpha \beta }+_\beta \left[\rho \nu \left(_\beta u_\alpha +_\alpha u_\beta \frac{2}{D}_\gamma u_\gamma \delta _{\alpha \beta }\right)\right]$$
(22)
where the viscosity is $`\nu =(\tau 0.5)/3`$ and $`D`$ is the number of spatial dimensions. In order to simulate a mixture of two components $`A`$ and $`B`$ we follow Orlandini and define a second lattice Boltzmann equation

$$g_i(𝐱+𝐯_i,t+1)=g_i(x,t)+\frac{1}{\tau _\varphi }\left[g_i^0(\rho (𝐱,t),\varphi (𝐱,t),𝐮(𝐱,t))g_i(𝐱,t)\right]$$
(23)
where the density difference (concentration) $`\varphi =\rho _A\rho _B`$ is defined as

$$\varphi =\underset{i}{}g_i$$
(24)
and is the only additional variable conserved in the dynamics. The total density is now $`\rho =\rho _A+\rho _B`$ and the concentration couples back into the momentum equation, eq.(22), via the pressure tensor $`P_{\alpha \beta }`$ which is now a function of both $`\rho `$ and $`\varphi `$. The relaxation step for the dynamics of $`g_i`$ is dependent on $`\tau _\varphi `$ which is related to the concentration diffusivity in the hydrodynamic limit. The equilibrium densities $`g_i`$ are determined, again by imposing that the local equilibrium concentration is recovered (from eq.(24)). The additional constraints on the $`g_i^0`$ are

$$\varphi u_\alpha =\underset{i}{}g_i^0v_{i\alpha },_{\alpha \beta }=\underset{i}{}g_i^0(v_{i\alpha }u_\alpha )(v_{i\beta }u_\beta )$$
(25)
Analogously to the pressure tensor, $``$ is given by a correction term for Galilean invariance  and the chemical potential $`\mu `$:

$$_{\alpha \beta }=\mathrm{\Gamma }\mu \delta _{\alpha \beta }+\frac{\tau _\varphi 0.5}{3}(u_\alpha _\beta \varphi +u_\beta _\alpha \varphi )$$
(26)
The desired thermodynamics of the fluid is determined by $`\mu `$ which must be consistent with the chosen $`P_{\alpha \beta }`$. These two quantities are normally derived from a given free energy. In this way one ensures their consistency, and can fix beforehand the appropriate thermodynamic behavior of the binary fluid model. It is enough to consider equilibrium densities $`g_i`$ that depend quadratically on the velocities to impose the previous requirements.

Again, we recover the hydrodynamic limit of eq.(23) performing a Taylor expansion of to second order in spatial gradients. In this case we get a convection-diffusion equation for the concentration

$$_t\varphi +_\alpha (\varphi u_\alpha )=_\alpha \left\{(\tau _\varphi 0.5)\left[_\beta \left(\mathrm{\Gamma }\mu (\rho ,\varphi )\delta _{\alpha \beta }\right)+\frac{\varphi }{\rho }_\beta P_{\alpha \beta }\right]\right\}$$
(27)
This model has been used extensively for the simulation of two-phase flows .The following derivation of the LEbc presented in the next section is, however, even more general and can be applied to most lattice Boltzmann methods.

## 3 Lees Edwards’ boundary conditions

The implementation of Lees Edwards boundary conditions for lattice Boltzmann will follow the same steps as in its original implementation: we will apply periodic boundary conditions to the distribution function in the two directions perpendicular to the velocity gradient; and in the third, we will assume that the periodic system is moving as a whole with a velocity $`\mathrm{\Delta }𝐮`$ with respect to the lattice. In a second step, we will have to displace the densities $`f_i`$ (and $`g_i`$ for the binary model), taking into account that we will have to interpolate the densities from the off-lattice positions to which they stream, onto the regular lattice. So this boundary condition will be applied after the collision and before the streaming of the densities.

Let us first focus on the velocity transformation. Since we are dealing now with a distribution function, the way to implement the velocity jump when crossing the system boundary in the velocity gradient direction will be to perform a Galilean transformation, i.e. if the distribution function in the system has momentum $`\rho 𝐮`$, its image node will have momentum $`\rho (𝐮+\mathrm{\Delta }𝐮)`$ if we move in the direction of the gradient, and $`\rho (𝐮\mathrm{\Delta }𝐮)`$ when we leave the system in the direction opposite to the velocity gradient. In this way, we will generate a linear velocity profile with shear rate $`\dot{\gamma }\mathrm{\Delta }𝐮/L_y`$. In order to get the momentum transfer, we need to calculate the difference between a density $`f_i`$ at the system velocity and the one after the Galilean transformation has been applied. All the macroscopic variables and their derivatives will be the same, since the image system is moving as a block at constant velocity. So, using eq.(20), where we retain only the $`u`$-dependence of the $`f`$’s, we obtain

$`\mathrm{\Delta }f_i`$ $`=`$ $`f_i(𝐮+\mathrm{\Delta }𝐮)f_i(𝐮)`$ (28)
$``$ $`f_i^0(𝐮+\mathrm{\Delta }𝐮)f_i^0(𝐮)`$
$`\tau \{_t[f_i^0(𝐮+\mathrm{\Delta }𝐮)f_i^0(𝐮)]+v_{i\alpha }_\alpha [f_i^0(𝐮+\mathrm{\Delta }𝐮)f_i^0(𝐮)]\}`$
where $`\mathrm{\Delta }f_i`$ is the change in the density moving in the $`i`$ direction due to the Galilean transformation. Note that the transformation will only be applied to the densities crossing the LE boundary but not to the rest of the densities at the same node.

The terms linear in $`\tau `$ are of order $`O()`$ and their zeroth and first moments are $`O(^3)`$. These terms are therefore small and enters only into the pressure moments. Therefore we can neglect them. We obtain for the boundary condition

$$f_i^{}=f_i+f_i^0(𝐮+\mathrm{\Delta }𝐮)f_i^0(𝐮).$$
(29)
where $`f_i^{}`$ means the density after the application of the boundary condition. We now need to specify the equilibrium distribution. We will assume the usual quadratic distribution function which can fulfill the conditions (15) and (16):

$$f_i^0=\rho a_0^i+\rho a_1^iv_{i\alpha }u_\alpha +a_2^iv_{i\alpha }v_{i\beta }\mathrm{\Pi }_{\alpha \beta }+a_3^i\text{tr}(\mathrm{\Pi })$$
(30)
where $`\mathrm{\Pi }_{\alpha \beta }=P_{\alpha \beta }^0+\rho u_\alpha u_\beta `$ and we obtain

$`f_i^0(𝐮+\mathrm{\Delta }𝐮)f_i^0(𝐮)`$ $`=`$ $`\rho [a_1^iv_{i\alpha }\mathrm{\Delta }u_\alpha `$ (31)
$`+a_2^iv_{i\alpha }v_{i\beta }(u_\alpha \mathrm{\Delta }u_\beta +u_\beta \mathrm{\Delta }u_\alpha +\mathrm{\Delta }u_\alpha \mathrm{\Delta }u_\beta )`$
$`+a_3^i(2𝐮.\mathrm{\Delta }𝐮+|\mathrm{\Delta }𝐮|^2)]`$
If we have a binary mixture, the same transformation will apply to the densities $`g_i`$ related to the concentration. The equilibrium distribution $`g_i^0`$ can be written in the same form as eq.(30) if we replace $`\rho `$ by $`\varphi `$ and $`P_{\alpha \beta }`$ by $`\mathrm{\Gamma }\mu \delta _{\alpha \beta }`$. In the previous expression $`\mathrm{\Gamma }`$ is a factor which, together with the relaxation parameter $`\tau _\varphi `$, will determine the diffusivity of the concentration; $`\delta _{\alpha \beta }`$ is the Kronecker delta function.

Now that we have defined a suitable Galilean transform, we have to split the densities into an untransformed part and a transformed part that will stream over the boundary, i.e. the densities at the top and bottom boundaries with associated velocities $`v_{iy}>0`$ and $`v_{iy}<0`$ respectively if we take $`𝐞_y`$ as the direction of the shear gradient. It is important to ensure that the two subsets of velocities conserve the node density after the Galilean transformation has been applied. This requires

$`\rho `$ $`=`$ $`{\displaystyle \underset{v_{iy}>0}{}}f_i^{}+{\displaystyle \underset{v_{iy}0}{}}f_i`$
$`=`$ $`{\displaystyle \underset{v_{iy}>0}{}}\left[f_i+f^0(𝐮+𝚫𝐮)f^0(𝐮)\right]+{\displaystyle \underset{v_{iy}0}{}}f_i`$
$`=`$ $`{\displaystyle \underset{v_{iy}>0}{}}\left[f^0(𝐮+𝚫𝐮)f^0(𝐮)\right]+{\displaystyle \underset{i}{}}f_i`$
$`0`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{(v_{iy}^2+v_{iy})}{2}}\left[f^0(𝐮+𝚫𝐮)f^0(𝐮)\right]`$ (32)
$`=`$ $`{\displaystyle \frac{1}{2}}(P_{yy}+\rho u_yu_yP_{yy}\rho u_yu_y)`$
where $`_{v_{iy}>(<)0}`$ is performed over the velocity directions with a positive (negative) component crossing the shear boundary. Eq.(32) is fulfilled if we require that the velocity subset that characterizes the lattice Boltzmann model involves only velocities that do not span beyond the closest site layer, i.e. the factor $`(v_{iy}^2+v_{iy})/2`$ is 1 for velocities with $`v_y>0`$ and $`0`$ otherwise. This is the case for most standard velocity sets and for the velocity set we will consider here (see eq.(37)) in particular. To derive this result we have also used that the velocity-shift is parallel to the boundary, i.e. $`\mathrm{\Delta }u_y=0`$.

It is worth remembering that the total momentum in the direction of the velocity gradient in the system must be zero; otherwise there will be a net mass transfer across the LE boundaries which would result in an increase (or decrease) of momentum in the direction of the velocity gradient. The mean velocity will then steadily increase until it exceeds the maxium velocity and the simulation becomes numerically unstable. This is a general requirement of any Lees-Edwards boundary condition, because the presence of LE boundaries breaks the momentum conservation unless the mass transfer is exactly balanced in both directions.

The idea of the Galilean transformation that we have derived also has more general applications and can be used to derive the boundary conditions for a solid moving wall. At a solid wall, the velocity of the fluid must be that of the wall. This no-slip boundary condition for a stationary wall is implemented through the ‘bounce-back’ method . In it, one replaces the streaming step of eqn.(13) for the densities that would have been streamed through the wall with

$$f_i(𝐱,t+1)=f_i(𝐱,t)$$
(33)
where $`v_i=v_i`$. To generalize this bounce back mechanism for a moving wall we first perform a Galilean transformation of the densities into a reference frame in which the wall is at rest. We then perform the bounce back and then transform the densities back with an inverse Galilean transformation. We then obtain for the bounce back off a wall moving with velocity $`U`$

$$f_i(𝐱,t+1)=f_i(𝐱,t)+2\rho a_1^2v_{i\alpha }U_\alpha $$
(34)
which is a result first derived with a different method by Ladd\[4, eq.(3.3)\]. For this particular case, the change in the distribution function is linear in the velocity jump, rather than quadratic, as in the general case (cf. eq.31) 111 According to the discussion of the previous paragraph, the density at both sides of the boundary is not modified by the bounce in the case of a moving parallel solid wall that slides along one of the boundaries of the system. For more general situations, this is not the case. For example, in the case of a moving rigid object embedded in the fluid, its velocity is not parallel to the boundary. In this case, the density inside the rigid object will be modified as a function of time. However, the presence of a rigid interface separating the inner from the outer motion of the fluid makes it possible to correct for these deviations..

To complete the implementation of Lees-Edwards boundary conditions we also need to move the densities by the displacement of the image system, $`d_x=t\mathrm{\Delta }u_x`$, which is not going to be an integer in general. We must then split the displacement in an integer part $`d_x^I`$ and a real part $`d_x^R`$ with $`0d_x^R<1`$. We then shift the densities that will stream over the boundary with $`v_{iy}>0`$ and $`v_{iy}<0`$ respectively by

$$f^{}(𝐱^\pm ,t)=(1d_x^R)f(𝐱+d_x^I,t)+d_x^Rf(𝐱+d_x^I\pm 1,t)$$
(35)
where $`x^\pm `$ are the positions of the lattice points directly below or above the LE boundary and this shift and interpolation is applied to densities $`f_i`$ with $`v_{iy}=\pm 1`$. That is, we have performed a linear partitioning of the moving density between the two corresponding nearest nodes. If necessary, it might be possible to improve upon this linear interpolation. This simple choice, however, already gives very accurate results.

Lees Edwards boundary conditions for lattice Boltzmann were first derived by Wagner & Yeomans  whose approach also required

$$\mathrm{\Delta }f_i=f_i^0(𝐮)f_i^0(𝐮)$$
(36)
but the authors also insisted on a heuristic transformation rule where the amount of velocity added to the densities that stream over the boundary should correspond to $`\mathrm{\Delta }𝐮`$. In particular they first calculated a separate x-velocity for the particles that crossed the LEbc and then they required that this x-velocity should be independent of the y-velocity in equilibrium. In order to fulfill this additional requirement they needed a more complicated equilibrium distribution (a polynomial fourth order in $`𝐮`$). While this additional requirement is perfectly sensible in the limit of the continuous Boltzmann equation our results show that it is not always helpful to consider the densities to be those of real particles and that it is perfectly feasible to have a Lees Edwards boundary condition for a lattice Boltzmann method with a traditional quadratic equilibrium distribution.

## 4 Numerical Implementation

We implemented a lattice with dimensions $`L_x,L_y`$ and the boundaries of the simulation lattice have periodic boundary conditions. We implemented a D2Q9 model, a two dimensional model with velocities

$$𝐯_i=\{\left(\begin{array}{c}\hfill 0\\ \hfill 0\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 0\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 0\end{array}\right),\left(\begin{array}{c}\hfill 0\\ \hfill 1\end{array}\right),\left(\begin{array}{c}\hfill 0\\ \hfill 1\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 1\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 1\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 1\end{array}\right),\left(\begin{array}{c}\hfill 1\\ \hfill 1\end{array}\right)\}$$
(37)
and with the coefficients given by Table 1.

### 4.1 Multiple LE planes

One of the limitations of Lattice-Boltzmann models to study fluids under shear arises because the maximum velocity cannot be larger than one lattice spacing per time step (and in practice must be $`\stackrel{<}{}0.1`$). This gives an upper bound of order $`1/L_y`$ for the shear rates that can be reached with the method described so far (one cannot simply increase the lattice size while keeping the shear rate constant). This suggests that large systems with Lees-Edwards boundary conditions should be restricted to small shear-rates. We can now overcome this limitation using the present method by introducing a larger number of Lees-Edwards shear planes as one increases the $`L_y`$ dimension of the simulation. We will still require, of course, that the velocity in each of the sub-lattices separated by the LEbc are no larger than about 0.1.

We will allow for $`N^{LE}`$ Lees-Edwards shear planes at almost equidistant positions

$$y=Y^k=\left[(k1/2)L_y/N^{LE}\right]_I\text{ for }k=\{1,\mathrm{},N^{LE}\}$$
(38)
where $`[x]_I`$ is the largest integer smaller than x. The system is then effectively cut into $`N^{LE}1`$ sub-lattices that will move relative to each other with the associated jump-velocities $`U^k`$. The average shear rate for the system is therefore

$$\dot{\gamma }=\frac{_{k=1}^{N^{LE}}U_k}{L_y}$$
(39)
and the stationary solution due to this set of LE shear planes is a linear shear profile of the form

$$𝐮=\left(\begin{array}{c}\dot{\gamma }(yy_0^k)\\ 0\\ 0\end{array}\right)$$
(40)
where $`y_0^k`$ is the position of zero velocity for each sub-lattice.

At this point it may appear that the introduction of the shear-planes does introduce a local velocity profile and biases the flow but, in the absence of errors, the only the global shear rate $`\dot{\gamma }`$, defined in eq.(39) matters to the flow. This is easily seen if one considers a Molecular Dynamics simulation where the introduction of several shear planes is equivalent to describing the different sublattices in different frames of references, with no effect on the real dynamics. Because of this there is, of course, no merit in introducing more than one LE-plane (which establishes $`\dot{\gamma }`$) in MD simulations. The only limit for the use of many LE-planes in lattice Boltzmann lies in the errors induced by the additional internal boundary condition and we discuss them in the following sections.

### 4.2 Transient start of shear and stationary profile

As a first test for the new boundary conditions we will examine the accuracy with which the algorithm simulates the start-up of a shear in a single component fluid. We consider a system with one LEbc and the initial velocity is zero. This corresponds to a system in which there is a periodic stack of slabs of fluid of width $`L_y`$ each of which moves with a different velocity. The velocity of the nth slab is given by $`u_x=n\mathrm{\Delta }u`$. The system will then exhibit viscous friction and the profile will converge toward a linear shear-profile. The density and pressure in this system remain constant and the velocity only has an $`x`$-component and its value only changes in the $`y`$-direction. The Navier-Stokes equation for this system gives

$$_tu_x=_y^2u_x$$
(41)
i.e. a simple diffusion equation for the momentum. The solution for this equation for $`L_y=\mathrm{}`$ is well know \[12, §52\] and given by

$$u_x(y,t)=\text{erf}\left(\frac{y}{2\sqrt{\nu t}}\right)$$
(42)
where

$$\text{erf}(x)=\frac{2}{\sqrt{\pi }}_0^xe^{\xi ^2}𝑑\xi .$$
(43)
Obtaining the solution for finite $`L_y`$ is simply a linear superposition

$$u_x(y,t)=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\frac{k}{|k|}\left[1\text{erf}\left(\frac{y+kL_y}{2\sqrt{\nu t}}\right)\right]$$
(44)
which can be conveniently evaluated numerically. The agreement between the theoretical prediction and the measured velocity is very good as can be seen in Figure 1. The error shown in the inset of Figure 1 is largest for the early time simulations but is always below 1% of the maximum velocity.

### 4.3 Shear in a two component model

The previous simulation was effectively a one-component simulation and it did not test the effectiveness of the interpolation scheme. Unfortunately, it is not easy to find two-component flows under shear for which analytical solutions are known. There is, however, a method to test the errors introduced by the LEbc for a two-component system very accurately. We can choose to have two equidistant LEbc with equal but opposite velocities so that the effective shear-rate is zero. The resulting system is then equivalent to a system at rest. We will consider two periodic stripes of A-rich and B-rich material at rest and simulate it with two LEbc with opposite velocities.

In figure 2 we show the results of a simulation on a 40x60 lattice with a stripe that was 20 lattice spacings wide. There were two LEbc with velocities of 0.005 and -0.005. The figure shows a section of 20x20 after 20000 iterations ($`d_x=0`$) where the stripe interface overlaps the LE boundary. The full system, of course, comprises of two interfaces in the x-direction and two LEbc in the y-direction and there the same error-patterns appear by symmetry.

Without errors this system would show perfectly straight contour lines and velocities of 0.005 in the upper half of figure 2(a) and -0.005 in the lower half. Since the velocities are very nearly that, we have to examine the difference of the real velocity and the predicted velocity shown in figure 2(b). We see that this error corresponds to a deviatoric velocity profile that flows into the interface at the LEbc. We attribute this error to the fact that we use the linear interpolation from eq.(35). The effect of this interpolation is to smear out the interface and make it wider than the equilibrium profile. There will then be a Marangoni flow that tries to relax the profile to its equilibrium shape which we assume is the major contribution to this velocity profile. A close examination of the interface reveals a slight widening around the LEbc. The maximum amplitude of the error in the velocity is 3%. In future work we plan to examine if we can reduce this error by using an interpolation method that will conserve not only the mass but also the derivatives and prevents the additional diffusion introduced by the LEbc. Nonetheless, the accuracy found above is comparable to other sources of error (e.g. anisotropy) present in lattice Boltzmann as currently used.

As a final benchmark of the performance of LEbc’s, we have studied the shape of a drop in a uniform shear flow generated by a single LE plane. In order to check for artifacts induced by the presence of the boundary, we have compared a simulation where the drop does not cross the boundary, and a situation where the LE plane bisects the drop. In Fig.3 we show the comparison for two different capillary numbers $`Ca`$ -defined as the ratio between viscous and interfacial forces, $`Ca=\nu R\dot{\gamma }/\sigma `$, with $`R`$ the radius of the drop and $`\sigma `$ its surface tension. The interface is located in the simulation through an interpolation of the corresponding nearby nodes where the order parameter changes sign, which introduces a small source of error. The drop that is bisected by the boundary will have a larger averaged velocity. Since we are comparing the drop shapes at equal time, we have displaced the drop accordingly to make them overlay. In the figure one can clearly see that the intersection of the boundary does not affect the shape of the drop, within the numerical accuracy already mentioned in the previous paragraph.

## 5 Conclusions

In this paper we have described how to implement Lees-Edwards periodic boundary conditions for lattice-Boltzmann models. Lees Edwards boundary conditions are a useful for the simulation of bulk systems under shear without the complicating influence of walls. In this way, finite size effects are diminished, which has converted the uniform shear flow in an attractive steady configuration to analyze computationally.

The implementation we have described follows the basic steps of the original paper, and is compatible with the standard LB models, that assume an equilibrium distribution for the densities which is quadratic in the velocity. The numerical examples described show that the method is quantitatively accurate for single and even for two component systems.

We have also described how to avoid the limitations on the attainable shear rates that are imposed by the underlying lattice. By adding a number of shearing planes, steady states of large shear rates can then me reached (subject only to the intrinsic stability limitations of any lattice-Boltzmann schemes). Because of the errors introduced by the interpolation scheme, it is not advisable to introduce more LEbc’s than necessary. One could think of a system where the whole lattice is subject to LEbc but this would now lead to a non-isotropic diffusion. We are planing to improve the interpolation scheme in order to reduce the artificial diffusion.

## Acknowledgements

We thank M.E. Cates for his critical reading of this manuscript. The work was funded in part under EPSRC Grant GR/M56234 and EC Access to Research Infrastructure Contract HPRI-1999-CT-00026 TRACS programme at EPCC.


text:
# The ISM Interactions of a Runaway LBV Nebula in the LMC

## 1. Introduction

Small ring nebulae have been detected around luminous blue variables (LBVs) in the Milky Way and the Large Magellanic Cloud (LMC) (Nota et al. 1995). Spectroscopic observations of such nebulae show enrichment of nitrogen and depletion of oxygen, indicating that these nebulae contain processed stellar material ejected by the LBVs (e.g., Smith et al. 1998). A large range of expansion velocities are observed in the dozen or so known LBV nebulae: the $`\eta `$ Car nebula has the largest expansion velocity, $``$600 km s-1; a few nebulae have expansion velocities $``$ 50 km s-1, e.g., P Cyg, HR Car, and AG Car; and the slowest expansion velocities are 20-30 km s-1, e.g., R127 and S119 (Nota et al. 1995). These diverse expansion dynamics indicate that either these LBV nebulae are at different evolutionary stages or they have different formation mechanisms.

Two LBVs projected against the LMC, R127 and S119, are surrounded by ring nebulae with similar angular sizes and expansion velocities (Nota et al. 1995); however, they have very different radial velocities. R127 has a radial velocity consistent with the systemic velocity of the LMC, Vhel = 250–300 km s-1, suggesting it is a member of the LMC (Walborn 1982). S119, on the other hand, has a radial velocity that is blue-shifted $``$100 km s-1 from the LMC systemic velocity. Furthermore, the interstellar Na I absorption in the spectrum of S119 does not show an LMC component (Nota et al. 1994). Thus the membership of S119 in the LMC cannot be automatically assumed.

Whether S119 was formed in the LMC or the Milky Way, its velocity implies that it has been ejected by its host galaxy at a high speed. The evolutionary status of S119 is thus very interesting. We have obtained a high-dispersion, long-slit echelle observation of the nebula around S119 and high-dispersion far ultraviolet spectrum of the star S119. These spectra are used to examine the dynamics of the nebula, study the interaction between the LBV nebula and the ambient interstellar medium (ISM), and determine whether S119 belongs to the Galaxy or the LMC.

## 2. Observations

On 2000 December 9, we obtained a long-slit, high-resolution spectrum of S119 (aka HD 269687, aka Sk$``$69175) using the echelle spectrograph on the CTIO 4-meter telescope. The principal aim of the observing program was to investigate the structure of the Magellanic Cloud ISM near approximately 70 hot stars observed by the Far Ultraviolet Spectroscopic Explorer (FUSE). The data presented here was a 900 second exposure with a 1$`\stackrel{}{\mathrm{.}}`$64$`\times `$4′ slit oriented east-west. A post-slit H$`\alpha `$ filter ($`\lambda _c`$=6563 Å; $`\mathrm{\Delta }\lambda `$=75 Å) restricted our spectral coverage to 6506–6654 Å, a single echelle order which contains H$`\alpha `$ and \[N II\]$`\lambda \lambda `$6548.05, 6583.45 lines. The spectrograph resolution was approximately 12 km s-1$`\times `$2$`\stackrel{}{\mathrm{.}}`$3. The data were reduced using standard techniques within IRAF and are presented in Figure 1.

To supplement our high resolution optical spectrum, we have inspected the FUSE observation of S119. Technical specifications of the FUSE spacecraft are described in Moos et al. (2000) and Sahnow et al. (2000). The single 3780 s FUSE observation (P11748) was taken on 2000 October 2 and covers the wavelength range 905-1187 Å at a resolution comparable to our echelle data ($``$15 km s-1). The FUSE spectral range is unusually rich in ionic, atomic and molecular transitions. We will use the FUSE data to probe gas seen in absorption along the line of sight to S119. Of particular interest are the O I $`\lambda `$1039.23 which probes total H I column and Fe II $`\lambda `$1144.94 which probes ionized material similar to that shown by H$`\alpha `$ and \[N II\]. Additionally, we examine C I $`\lambda `$945.191 and several H2 lines which probe cold gas.

Finally, we consider archival HST imagery of S119 in our analysis. Four WFPC2 exposures from 1998 March observed with the H$`\alpha `$ (F656N) filter were coadded for a total exposure length of 2000 seconds. As the star was centered in the Planetary Camera for this observation, we see the circumstellar nebula at unprecedented resolution. This image is presented in Figure 2.

## 3. Discussion

### 3.1. Expanding Shell

The most obvious feature of Figure 1 is the expanding shell seen in both H$`\alpha `$ and \[N II\]. This shell is easily spatially resolved along the slit and shows a diameter of 7$`\stackrel{}{\mathrm{.}}`$5$`\pm `$0$`\stackrel{}{\mathrm{.}}`$5 or $``$2 pc at the distance of the LMC. In velocity space, the shell shows well-resolved expansion with two sides separated by about 50 km s-1 for an expansion velocity of 25 km s-1. This expansion velocity is on the low side for LBV nebulae but is comparable to that of R127. The red and blue-shifted components of the shell in each line show a width only slightly higher than the resolution of our spectrograph ($``$15 km s-1).

Measuring the east-west diameter of the bubble in Figure 2, we find 7$`\stackrel{}{\mathrm{.}}`$2$`\pm `$0.5″, though the western edge of the nebulosity is somewhat hard to define. In the north-south direction, we get a slightly longer diameter of 7$`\stackrel{}{\mathrm{.}}`$8$`\pm `$0$`\stackrel{}{\mathrm{.}}`$5 (not including the short tail to the southwest). For comparison, the diameter and expansion velocity reported by Nota et al. (1994), measured from ground-based data, are 8″$`\times `$9″ and 25 km s-1.

### 3.2. Where is S119?

By comparing the emission features seen in the optical spectrum with the absorption features in the FUSE spectrum, we can investigate the three-dimensional structure of S119 and its surroundings. By looking at low-ionization gas (O I, Fe II, etc.) we should see in absorption much of the same material seen in emission in H$`\alpha `$.

The upper panel of Figure 3 shows the observed velocity structure in H$`\alpha `$ and \[N II\] extracted from our echelle data. The solid curves show the stellar spectrum in H$`\alpha `$ and \[N II\] and the dashed curve shows the averaged spectrum of an adjacent few arcseconds of ISM in H$`\alpha `$. We see a pair of emission peaks at 140 and 190 km s-1 corresponding to the expanding shell. In the background spectrum, we see a broad band of LMC H$`\alpha `$ emission between 220 and 290 km s-1.

The lower panel of Figure 3 shows the normalized absorption profiles of O I $`\lambda `$1039 and Fe II $`\lambda `$1145 from the FUSE data. At least five distinct velocity components are detected with good agreement between the two species. A strong line at 20 km s-1 arises in local Galactic material. This is consistent with the absorption components at 18 km s-1 seen by Roth & Blades (1997) in Zn II and Cr II, two ions with very similar ionization potentials to Fe II. At 60 km s-1 we see what is most likely an intermediate-velocity cloud (probably located in the Galactic halo).

There are three higher velocity absorption components seen at 135, 215 and 290 km s-1. The strong feature at 135 km s-1 may be associated with the near (blue-shifted) side of the LBV nebula seen in emission at 140 km s-1. However, this feature may also arise in high velocity cloud material from the Galactic halo which is often seen at $``$125 km s-1 toward the LMC (Richter et al. 1999). We do not expect, nor do we detect any absorption features corresponding to the receding side of the bubble at 190 km s-1.

The 215 and 290 km s-1 absorption components represent LMC ISM absorption and are typical of LMC sightlines (Danforth et al. in prep). Roth & Blades see an absorption component at 283 km s-1 which is almost certainly the same as our 290 km s-1 component. A reasonably good correlation is seen between these absorption features and the H$`\alpha `$ emission at 220–290 km s-1; they may trace the same material. Since LMC gas is seen in absorption, S119 must be behind some portion of the LMC ISM.

At first glance, this appears to be in contradiction to the findings of Nota et al. (1994) who detected no absorption in Na I at LMC velocities. Na is far less abundant than either O or Fe in the solar neighborhood (presumably in the metal-poor LMC as well) and thus a non-detection of Na I does not rule out a detectable quantity of O I or Fe II. Furthermore, Na I has an ionization potential of only 5.1 eV and is found in the cold, neutral clouds, while O I is ionized at 13.6 eV and traces the total H I column and can exist in conciderably warmer conditions.

Better tracers of cold gas are C I which has an ionization potential of 11.26 eV and molecular hydrogen which is dissociated at 4.48 eV. The The FUSE data do not show significant H2 or C I $`\lambda `$945.191 absorption at LMC velocities, suggesting that there is a very low column depth of cold neutral gas in our line of sight. The ISM along the sightline toward S119 must be dominated by ionized gas at warm to high temperatures, which can be detected in O I but not Na I or other cold gas tracers.

### 3.3. ISM Interactions

The HST image reveals that the northeastern side of S119’s nebula is brighter and has a sharper outer edge than the southwestern side. The image also shows that S119 is slightly displaced from the nebular center toward the northeast. These characteristics are typically seen in planetary nebulae moving through the ISM (Tweedy & Kwitter 1996). The echelle image reveals that the near (blue-shifted) side of the bubble is brighter than the far (red-shifted) side by a factor of 3.3 in both H$`\alpha `$ and \[N II\] lines. While we may simply be seeing an inhomogeneity in the nebular density or internal extinction in the nebula, we have shown that S119 lies behind a significant amount of LMC ISM and is moving through it at $`>`$100 km s-1. These points seem to indicate that the bubble is interacting with the ISM on the near side in a bow shock. Taking the sharp, bright northeastern shell rim as the leading edge, we can conclude that the star is moving away from the LMC (toward us) and to the northeast at a sizable velocity.

The rapid movement of S119 and its nebula through the LMC ISM is further supported by the detection of the ‘wake’ of the bow shock in the echelle observation. Faintly seen in Figure 1 both in H$`\alpha `$ and the \[N II\] $`\lambda `$6584 lines is a ‘tail’ of emission stretching from the western tip of the expansion feature back toward the bulk of emission representing the LMC gas. This tail is brighter toward the bubble velocity than toward the LMC ISM velocity, indicating that it contains the nebular material that has been decelerated by the LMC ISM.

## 4. Conclusions

S119 has long been known to be an LBV candidate star with an expanding nebular ring. In this letter we have confirmed the expansion velocity of 25 km s-1 and refined the physical scales of the bubble to 7$`\stackrel{}{\mathrm{.}}`$2$`\times `$7$`\stackrel{}{\mathrm{.}}`$8 or 1.8$`\times `$2.0 pc at the distance of the LMC.

Given the anomalous stellar velocity of S119, v$`{}_{}{}^{}`$160 km s-1, membership in the LMC (or in the Galaxy) was in debate. The presence of low-ionization absorption at LMC velocities against the stellar continuum indicates that S119 resides within the LMC. Though Nota et al. (1994) did not see Na I absorption at LMC velocities, a similar dearth of C I and H2 suggest that the S119 sightline contains very little cold neutral gas. The star’s anomalously high velocity might be caused by the supernova explosion of a close companion star (the “sling-shot” effect; Blaauw 1961). It would be extremely interesting to search for a compact companion of S119.

The relative brightening on the near side of the expanding nebular shell as well as the velocity disparity between S119 and the LMC gas leads us to postulate a bowshock interaction. The appearance of a faint ‘tail’ on the western shell edge (Figure 1) and the relative brightening and sharp edge of the northeastern edge of the nebular shell (Figure 2) suggest a general northeastward motion as well.

## 5. Acknowledgements

The authors wish to acknowledge Alex Fullerton and FUSE Hot Star team for use of their proprietary FUSE data. Discussions with Chris Howk, Antonella Nota and Nichole King were particularly helpful as well. This work is based in part on data obtained for the Guaranteed Time Team by the NASA-CNES-CSA FUSE mission operated by the Johns Hopkins University. Generous travel support was provided by CTIO. Financial support has been provided by NASA contract NAS5-32985.


text:
# Scalar gravitational waves and Einstein frame

## Acknowledgments

We wish to thank R. Jackiw for a useful remark on the definition of the energy density in Einstein’s theory of gravity. SB acknowledges support by the “B. Rossi” program of exchange MIT/INFN and thanks warmly MIT-CTP for hospitality.