text:
# A PRIMER ON EULERIAN COMPUTATIONAL FLUID DYNAMICS FOR ASTROPHYSICS

## 1 Introduction

Astrophysical structure formation and the dynamics of astrophysical systems involve nonlinear gas dynamical processes which cannot be modeled analytically but require numerical methods. One would like to address the challenging problem of star formation and how this process produces planetary systems. Observations of the X-ray emission from hot gas in galaxy clusters, the Sunyaev-Zeldovich effect in the CMB spectrum, and the Lyman alpha forest in the spectra of quasars are only meaningful if we understand the gas dynamical processes involved. The evolution of complex systems is best modeled using numerical simulations.

A large class of astrophysical problems involve collisional systems where the mean free path is much smaller than all length scales of interest. Hence, one can appropriately adopt an ideal fluid description of matter where the thermodynamical properties of the fluid obey well known equations of state. Conservation of mass, momentum, and energy allows one to write down the Euler equations which govern fluid mechanics (See Landau & Lifshitz, 1987). This formalism is an ideal basis for simulating astrophysical fluids.

Hydrodynamical simulations are faced with challenging problems, but advancements in the field have made it an important tool for theoretical astrophysics. One of the main challenges in simulating complex fluid flows is the capturing of strong shocks, which frequently occur and play an important role in gas dynamics. The differential Euler equations are ill-defined at shock discontinuities where derivatives are infinite. Much effort has been devoted to solving this problem and a field of work has resulted from it. Computational fluid dynamics (CFD) is a powerful approach to simulating fluid flow with emphasis on high resolution capturing of shocks and prevention of numerical instabilities. Both Eulerian and Lagrangian methods have been developed.

Lagrangian methods based on smoothed particle hydrodynamics (SPH; Gingold & Monaghan, 1977; Lucy, 1977) consider a Monte-Carlo approximation to solving the fluid equations, somewhat analogous to $`N`$-body methods for the Vlasov equation. SPH schemes follow the trajectories of particles of fixed mass which represent fluid elements. The Lagrangian forms of the Euler equations are solved to determine smoothed fluid variables like density, velocity, and temperature. The particle formulation does not naturally capture shocks and artificial viscosity is added to prevent unphysical oscillations. However, the addition of artificial viscosity broadens shocks over several smoothing lengths and degrades the resolution. The Lagrangian approach has a large dynamic range in length but not in mass. It achieves good spatial resolution in high density regions but does poorly in low density regions. SPH schemes must smooth over a large number of neighbouring particles, making it computationally expensive and challenging to implement in parallel.

The standard approach to Eulerian methods is to discretize the problem and solve the integral Euler equations on a Cartesian grid by computing the flux of mass, momentum, and energy across grid cell boundaries. In conservative schemes, the flux taken out of one cell is added to the neighbouring cell and this ensures the correct shock propagation. Flux assignment schemes based on the total variation diminishing condition (Harten, 1983) have been designed for high order accuracy and high resolution capturing of shocks, while preventing unphysical oscillations. The Eulerian approach has a large dynamic range in mass but not in length, opposite to that of Lagrangian schemes. In general, Eulerian algorithms are computationally faster by several orders of magnitude. They are also easy to implement and to parallelize.

The purpose of this paper is to present a pedagogical review of some of the methods employed in Eulerian computational fluid dynamics. In $`\mathrm{\S }2`$ we briefly review the Euler equations and discuss the standard approach to discretizing conservation laws. We describe traditional central differencing methods such as the Lax-Wendroff scheme in $`\mathrm{\S }3`$ and more modern flux assignment methods like the TVD scheme in $`\mathrm{\S }4`$. In $`\mathrm{\S }5`$ we review the relaxing TVD method for systems of conservation laws like the Euler equations, which has been successfully implemented for simulating cosmological astrophysical fluids by Pen (1998). In $`\mathrm{\S }6`$ we apply a self-gravitating hydro code to simulating the formation of blue straggler stars through stellar mergers. A sample 3-D relaxing TVD code is provided in the appendix.

## 2 Eulerian Hydrodynamics

The Euler equations which govern hydrodynamics are a system of conservation laws for mass, momentum, and energy. In differential conservation form, the continuity equation, momentum equation, and energy equation are given as:

$$\frac{\rho }{t}+\frac{}{x_j}(\rho v_j)=0,$$
(1)
$$\frac{(\rho v_i)}{t}+\frac{}{x_j}(\rho v_iv_j+P\delta _{ij})=0,$$
(2)
$$\frac{e}{t}+\frac{}{x_j}[(e+P)v_j]=0.$$
(3)
We have omitted gravitational and other source terms like heating and cooling. The physical state of the fluid is specified by its density $`\rho `$, velocity field $`𝒗`$, and total energy density,

$$e=\frac{1}{2}\rho v^2+\epsilon .$$
(4)
In practice, the thermal energy $`\epsilon `$ is evaluated by subtracting the kinetic energy from the total energy. For an ideal gas, the pressure $`P(\epsilon )`$ is related to the thermal energy by the equation of state,

$$P=(\gamma 1)\epsilon ,$$
(5)
where $`\gamma `$ is the ratio of specific heats. Another thermodynamic variable which is of importance is the sound speed $`c_s`$ which is given by

$$c_s^2\frac{P}{\rho }=\frac{\gamma P}{\rho }.$$
(6)
The thermodynamical properties of an ideal gas obey well known equations of state, which we do not fully list here.

The differential Euler equations require differentiable solutions and therefore, are ill-defined at jump discontinuities where derivatives are infinite. In the literature, nondifferentiable solutions are called weak solutions. The differential form gives a complete description of the flow in smooth regions, but the integral form is needed to properly describe shock discontinuities. In integral conservation form, the rate of change in mass, momentum, and energy is equal to the net flux of those conserved quantities through the surface enclosing a control volume. For simplicity of notation, we will continue to express conservation laws in differential form, as a shorthand for the integral form.

### 2.1 Computational Fluid Dynamics

The standard approach to Eulerian computational fluid dynamics is to discretize time into discrete steps and space into finite volumes or cells, where the conserved quantities are stored. In the simplest case, the integral Euler equations are solved on a Cartesian cubical lattice by computing the flux of mass, momentum, and energy across cell boundaries in discrete time steps. Consider the Euler equations in vector differential conservation form,

$$\frac{𝒖}{t}+\frac{𝑭_i(𝒖)}{x_i}=0,$$
(7)
where $`𝒖=(\rho ,\rho v_x,\rho v_y,\rho v_z,e)`$ contains the conserved physical quantities and $`𝑭(𝒖)`$ represents the flux terms. In practice, the conserved cell-averaged quantities $`𝒖_n𝒖(𝒙_n)`$ and fluxes $`𝑭_n`$ are defined at integer grid cell centres $`𝒙_n`$. The challenge is to use the cell-averaged values to determine the fluxes $`𝑭_{n+1/2}`$ at cell boundaries.

In the following sections, we describe flux assignment methods designed to solve conservation laws like the Euler equations. For ease of illustration, we begin by considering a 1-D scalar conservation law,

$$\frac{u}{t}+\frac{F(u)}{x}=0,$$
(8)
where $`F(u)=vu`$ and $`v`$ is a constant advection velocity. Equation (8) is referred to as a linear advection equation and has the analytical solution,

$$u(x,t)=u(xvt,0).$$
(9)
The linear advection equation describes the transport of the quantity $`u`$ at a constant velocity $`v`$.

In integral flux conservation form, the 1-D scalar conservation law can be written as

$$\frac{}{t}_{x_1}^{x_2}u(x,t)𝑑x+_{x_1}^{x_2}\frac{F(u)}{x}𝑑x=0,$$
(10)
where $`x_1x_{n1/2}`$ and $`x_2x_{n+1/2}`$ for our control cells. Let $`F_{n+1/2}^t`$ denote the flux of $`u`$ through cell boundary $`x_{n+1/2}`$ at time $`t`$. Note then that the second integral is simply equal to $`F_{n+1/2}^tF_{n1/2}^t`$. The rate of change in the cell-integrated quantity $`u𝑑x`$ is equal to the net flux of $`u`$ through the control cell. For a discrete time step, the discretized solution for the cell-averaged quantity $`u_n`$ is given by

$$u_n^{t+\mathrm{\Delta }t}=u_n^t\left(\frac{F_{n+1/2}^tF_{n1/2}^t}{\mathrm{\Delta }x}\right)\mathrm{\Delta }t.$$
(11)
The physical quantity $`u`$ is conserved since the flux taken out of one cell is added to the neighbouring cell which shares the same boundary. Note that Equation (11) has the appearance of being a finite difference scheme for solving the differential form of the 1-D scalar conservation law. This is why the differential form can be used as a shorthand for the integral form.

## 3 Centered Finite-Difference Methods

Central-space finite-difference methods have ease of implementation but at the cost of lower accuracy and stability. For illustrative purposes, we start with a simple first-order centered scheme to solve the linear advection equation. The discretized solution is given by Equation (11) where the fluxes at cell boundaries,

$$F_{n+1/2}^t=\frac{F_{n+1}^t+F_n^t}{2},$$
(12)
are obtained by taking an average of cell-centered fluxes $`F_n^t=vu_n^t`$. The discretized first-order centered scheme can be equivalently written as

$$u_n^{t+\mathrm{\Delta }t}=u_n^t\left(\frac{F_{n+1}^tF_{n1}^t}{2\mathrm{\Delta }x}\right)\mathrm{\Delta }t.$$
(13)
In this form, the discretization has the appearance of using a central difference scheme to approximate spatial derivatives. Hence, centered schemes are often referred to as central difference schemes. In practice when using centered schemes, the discretization is done on the differential conservation equation rather than the integral equation.

This simple scheme is numerically unstable and we can show this using the von Neumann linear stability analysis. Consider writing $`u(x,t)`$ as a discrete Fourier series:

$$u_n^t=\frac{1}{N}\underset{k=N/2}{\overset{N/2}{}}c_k^t\mathrm{exp}\left(\frac{2\pi ikn}{N}\right),$$
(14)
where $`N`$ is the number of cells in our periodic box. In plane-wave solution form, we can write this as

$$u_n^t=\frac{1}{N}\underset{k=N/2}{\overset{N/2}{}}c_k^{}\mathrm{exp}\left[\frac{2\pi i(kn\omega t)}{N}\right],$$
(15)
where $`c_k^{}`$ are the Fourier series coefficients for the initial conditions $`u(x,0)`$. Equivalently, the time evolution of the Fourier series coefficients in Equation (14) can be cast into a plane-wave solution of the form,

$$c_k^t=\mathrm{exp}\left(\frac{2\pi i\omega t}{N}\right)c_k^{},$$
(16)
where the numerical dispersion relation $`\omega (k)`$ is complex in general. The imaginary part of $`\omega `$ represents the growth or decay of the Fourier modes while the real part describes the oscillations. A numerical scheme is linearly stable if $`\text{Im}(\omega )0`$. Otherwise, the Fourier modes will grow exponentially in time and the solution will blow up.

The exact solution to the linear advection equation can be expressed in the form of Equation (9) or by a plane-wave solution where the dispersion relation is given by $`\omega _{}=vk`$. The waves all travel at the same phase velocity $`\omega _{}/k=v`$ in the exact case.

The centrally discretized linear advection equation (Equation 13) is exactly solvable. After $`m`$ times steps, the time evolution of the independent Fourier modes is given by

$$c_k^{m\mathrm{\Delta }t}=\left(1i\lambda \mathrm{sin}\varphi \right)^mc_k^{},$$
(17)
where $`\lambda v\mathrm{\Delta }t/\mathrm{\Delta }x`$ and $`\varphi =2\pi k\mathrm{\Delta }x/N`$. The dispersion relation is given by

$$\omega =\frac{N}{2\pi \mathrm{\Delta }t}\left[\mathrm{tan}^1(\lambda \mathrm{sin}\varphi )+\frac{i}{2}\mathrm{ln}(1+\lambda ^2\mathrm{sin}^2\varphi )\right],$$
(18)
For any time step $`\mathrm{\Delta }t>0`$, the imaginary part of $`\omega `$ will be $`>0`$. The Fourier modes will grow exponentially in time and the solution will blow up. Hence, the first-order centered scheme is numerically unstable.

### 3.1 Lax-Wendroff Scheme

The Lax-Wendroff scheme (Lax & Wendroff, 1960) is second-order accurate in time and space and the idea behind it is to stabilize the unstable first-order scheme from the previous section. Consider a Taylor series expansion for $`u(x,t+\mathrm{\Delta }t)`$:

$$\begin{array}{c}u(x,t+\mathrm{\Delta }t)=u(x,t)\hfill \\ \hfill +\frac{u}{t}\mathrm{\Delta }t+\frac{^2u}{t^2}\frac{\mathrm{\Delta }t^2}{2}+𝒪(\mathrm{\Delta }t^3),\end{array}$$
(19)
and replace the time derivatives with spatial derivatives using the conservation law to obtain

$$\begin{array}{c}u(x,t+\mathrm{\Delta }t)=u(x,t)\hfill \\ \hfill \frac{F}{x}\mathrm{\Delta }t+\frac{}{x}\left(\frac{F}{u}\frac{F}{u}\frac{u}{x}\right)\frac{\mathrm{\Delta }t^2}{2}+𝒪(\mathrm{\Delta }t^3).\end{array}$$
(20)
For the linear advection equation, the eigenvalue of the flux Jacobian is $`F/u=v`$. Discretization using central differences gives

$$\begin{array}{c}u_n^{t+\mathrm{\Delta }t}=u_n^t\frac{F_{n+1}^tF_{n1}^t}{2\mathrm{\Delta }x}\mathrm{\Delta }t\hfill \\ \hfill +\left(\frac{F_{n+1}^tF_n^t}{\mathrm{\Delta }x}\frac{F_n^tF_{n1}^t}{\mathrm{\Delta }x}\right)\frac{v\mathrm{\Delta }t^2}{2\mathrm{\Delta }x}.\end{array}$$
(21)
In conservation form, the solution is given by Equation (11), where the fluxes at cell boundaries are defined as

$$F_{n+1/2}^t=\frac{1}{2}\left(F_{n+1}^t+F_n^t\right)\left(F_{n+1}^tF_n^t\right)\frac{v\mathrm{\Delta }t}{2\mathrm{\Delta }x}.$$
(22)
Compare this with the boundary fluxes for the first-order scheme (Equation 12). The Lax-Wendroff scheme obtains second-order fluxes,

$$F^{(2)}=F^{(1)}\frac{F}{u}\frac{F}{x}\frac{\mathrm{\Delta }t}{2},$$
(23)
by modifying the first-order fluxes $`F^{(1)}`$ with a second-order correction.

The stability of the Lax-Wendroff scheme to solve the linear advection equation can also be determined using the von Neumann analysis. The discretized Lax-Wendroff equation (Equation 21) is exactly solvable and after $`m`$ time steps, the Fourier modes evolve according to

$$c_k^{m\mathrm{\Delta }t}=\left[1\lambda ^2(1\mathrm{cos}\varphi )i\lambda \mathrm{sin}\varphi \right]^mc_k^{},$$
(24)
where $`\lambda v\mathrm{\Delta }t/\mathrm{\Delta }x`$ is called the Courant number and $`\varphi =2\pi k\mathrm{\Delta }x/N`$. The dispersion relation is given by

$$\begin{array}{c}\omega =\frac{N}{2\pi \mathrm{\Delta }t}\mathrm{tan}^1\left[\frac{\lambda \mathrm{sin}\varphi }{1\lambda ^2(1\mathrm{cos}\varphi )}\right]\hfill \\ \hfill +\frac{iN}{4\pi \mathrm{\Delta }t}\mathrm{ln}\left[14\lambda ^2(1\lambda ^2)\mathrm{sin}^4\left(\frac{\varphi }{2}\right)\right].\end{array}$$
(25)
It is important to note three things. First, the Lax-Wendroff scheme is conditionally stable provided that $`\text{Im}(\omega )0`$, which is satisfied if

$$\frac{v\mathrm{\Delta }t}{\mathrm{\Delta }x}1.$$
(26)
This constraint is a particular example of a general stability constraint known as the Courant-Friedrichs-Lewyor (CFL) condition. The Courant number $`\lambda `$ is also referred to as the CFL number. Second, for $`\lambda =1`$ the dispersion relation is exactly identical to that of the exact solution and the numerical advection is exact. This is a special case, however, and it does not test the ability of the Lax-Wendroff scheme to solve general scalar conservation laws. Normally, one chooses $`\lambda <1`$ to satisfy the CFL condition. Lastly, for $`\lambda <1`$ the dispersion relation $`\omega (k)`$ for the Lax-Wendroff solution is different from the exact solution where $`\omega _{}=vk`$. The dispersion relation relative to the exact solution can be parametrized by

$$\mathrm{\Delta }\omega \omega \omega _{}.$$
(27)
The second-order truncation of the Taylor series (Equation 20) results in a phase error $`\text{Re}(\mathrm{\Delta }\omega )`$ which is a function of frequency. In the Lax-Wendroff solution, the waves are damped and travel at different speeds. Hence the scheme is both diffusive and dispersive.

In Figure (1) we plot the phase error $`\text{Re}(\mathrm{\Delta }\omega )`$ and the amplification term $`\text{Im}(\mathrm{\Delta }\omega )`$ for the Lax-Wendroff scheme with parameters $`N=100`$, $`v=1`$, and $`\lambda =0.9`$. A negative value of $`\text{Re}(\mathrm{\Delta }\omega )`$ represents a lagging phase error while a positive value indicates a leading phase error. For the chosen CFL number, the high frequency modes have the largest phase errors but they are highly damped. Some of the modes having lagging phase errors are not highly damped. We will subsequently see how this becomes important.

A rigourous test of the 1-D Lax-Wendroff scheme and other flux assignment schemes we will discuss is the linear advection of a square wave. The challenge is to accurately advect this discontinuous function where the edges mimic Riemann shock fronts. In Figure (2) we show how the Lax-Wendroff scheme does at advecting the square wave once (dashed line) and ten times (dotted line) through a periodic box of 100 grid cells at speed $`v=1`$ and $`\lambda =0.9`$. Note that this scheme produces numerical oscillations. Recall that a square wave can be represented by a sum of Fourier or sine waves. These waves will be damped and disperse when advected using the Lax-Wendroff scheme. Figure (1) shows that the modes having lagging phase errors are not damped away. Hence, the Lax-Wendroff scheme is highly dispersive and the oscillations in Figure (2) are due to dispersion. We leave it as an exercise for the reader to advect a sine wave using the Lax-Wendroff scheme. Since there is only one frequency mode in this case, there will be no spurious oscillations due to dispersion, but a phase error will be present. For a comprehensive discussion on the family of Lax-Wendroff schemes and other centered schemes, see Hirsch (1990) and Laney (1998).

## 4 Upwind Methods

Upwind methods take into account the physical nature of the flow when assigning fluxes for the discrete solution. This class of flux assignment schemes, whose origin dates back to the work of Courant, Isaason, & Reeves (1952), has been shown to be excellent at capturing shocks and also being highly stable.

We start with a simple first-order upwind scheme to solve the linear advection equation. Consider the case where the advection velocity is positive and flow is to the right. The flux of the physical quantity $`u`$ through the cell boundary $`x_{n+1/2}`$ will originate from cell $`n`$. The upwind scheme proposes that, to first-order, the fluxes $`F_{n+1/2}^t`$ at cell boundaries be taken from the cell-centered fluxes $`F_n^t=vu_n^t`$, which is in the upwind direction. If the advection velocity is negative and flow is to the left, the boundary fluxes $`F_{n+1/2}^t`$ are taken from the cell-centered fluxes $`F_{n+1}^t=vu_{n+1}^t`$. The first-order upwind flux assignment scheme can be summarized as follows:

$$F_{n+1/2}^t=\{\begin{array}{cc}F_n^t\hfill & \text{if }v>0\text{ ,}\hfill \\ F_{n+1}^t\hfill & \text{if }v<0\text{ .}\hfill \end{array}$$
(28)
Unlike central difference schemes, upwind schemes are explicitly asymmetric.

The CFL condition for the first-order upwind scheme can be determined from the von Neumann analysis. We consider the case of a positive advection velocity. After $`m`$ time steps, the Fourier modes evolve according to

$$c_k^{m\mathrm{\Delta }t}=\left[1\lambda (1\mathrm{cos}\varphi )i\lambda \mathrm{sin}\varphi \right]^mc_k^{},$$
(29)
where $`\lambda v\mathrm{\Delta }t/\mathrm{\Delta }x`$ and $`\varphi =2\pi k\mathrm{\Delta }x/N`$. The dispersion relation is given by

$$\begin{array}{c}\omega =\frac{N}{2\pi \mathrm{\Delta }t}\mathrm{tan}^1\left[\frac{\lambda \mathrm{sin}\varphi }{1\lambda (1\mathrm{cos}\varphi )}\right]\hfill \\ \hfill +\frac{iN}{4\pi \mathrm{\Delta }t}\mathrm{ln}\left[14\lambda (1\lambda )\mathrm{sin}^2\left(\frac{\varphi }{2}\right)\right].\end{array}$$
(30)
The CFL condition for solving the linear advection equation with this scheme is to have $`\lambda 1`$, identical to that for the Lax-Wendroff scheme. For $`\lambda <1`$ the dispersion relation $`\omega (k)`$ for the first-order upwind scheme is different from the exact solution where $`\omega _{}=vk`$. This scheme is both diffusive and dispersive. Since it is only first-order accurate, the amount of diffusion is large. In Figure (3) we compare the dispersion relation of the upwind scheme to that of the Lax-Wendroff scheme. The Fourier modes in the upwind scheme also have phase errors but they will be damped away. The low frequency modes which contribute to the oscillations in the Lax-Wendroff solution are more damped in the upwind solution. Hence, one does not expect to see oscillations resulting from phase errors.

In Figure (4) we show how the first-order upwind scheme does at advecting the Riemann shock wave. This scheme is well-behaved and produces no spurious oscillations, but since it is only first-order, it is highly diffusive. The first-order upwind scheme has the property of having monotonicity preservation. When applied to the linear advection equation, it does not allow the creation of new extrema in the form of spurious oscillations. The Lax-Wendroff scheme does not have the property of having monotonicity preservation.

The flux assignment schemes that we have discussed so far are all linear schemes. Godunov (1959) showed that all linear schemes are either diffusive or dispersive or a combination of both. The Lax-Wendroff scheme is highly dispersive while the first-order upwind scheme is highly diffusive. Godunov’s theorem also states that linear monotonicity preserving schemes are only first-order accurate. In order to obtain higher order accuracy and prevent spurious oscillations, nonlinear schemes are needed to solve conservation laws.

### 4.1 Total Variation Diminishing Schemes

Harten (1983) proposed the total variation diminishing (TVD) condition which guarantees that a scheme have monotonicity preservation. Applying Godunov’s theorem, we know that all linear TVD schemes are only first-order accurate. In fact, the only linear TVD schemes are the class of first-order upwind schemes. Therefore, higher order accurate TVD schemes must be nonlinear.

The TVD condition is a nonlinear stability condition. The total variation of a discrete solution, defined as

$$TV(u^t)=\underset{i=1}{\overset{N}{}}|u_{i+1}^tu_i^t|,$$
(31)
is a measure of the overall amount oscillations in $`u`$. The direct connection between the total variation and the overall amount of oscillations can be seen in the equivalent definition

$$TV(u^t)=2\left(u_{\text{max}}u_{\text{min}}\right),$$
(32)
where each maxima is counted positively twice and each minima counted negatively twice (See Laney, 1998). The formation of spurious oscillations will contribute new maxima and minima and the total variation will increase. A flux assignment scheme is said to be TVD if

$$TV(u^{t+\mathrm{\Delta }t})TV(u^t),$$
(33)
which signifies that the overall amount of oscillations is bounded. In linear flux-assignment schemes, the von Neumann linear stability condition requires that the Fourier modes remain bounded. In nonlinear schemes, the TVD stability condition requires that the total variation diminishes.

We now describe a nonlinear second-order accurate TVD scheme which builds upon the first-order monotone upwind scheme described in the previous section. The second-order accurate fluxes $`F_{n+1/2}^t`$ at cell boundaries are obtained by taking first-order fluxes $`F_{n+1/2}^{(1),t}`$ from the upwind scheme and modifying it with a second order correction. First consider the case where the advection velocity is positive. The first-order upwind flux $`F_{n+1/2}^{(1),t}`$ comes from the averaged flux $`F_n^t`$ in cell $`n`$. We can define two second-order flux corrections,

$`\mathrm{\Delta }F_{n+1/2}^{L,t}`$ $`={\displaystyle \frac{F_n^tF_{n1}^t}{2}},`$ (34)
$`\mathrm{\Delta }F_{n+1/2}^{R,t}`$ $`={\displaystyle \frac{F_{n+1}^tF_n^t}{2}},`$ (35)
using three local cell-centered fluxes. We use cell $`n`$ and the cells immediately left and right of it. If the advection velocity is negative, the first-order upwind flux comes from the averaged flux $`F_{n+1}^t`$ in cell $`n+1`$. In this case, the second-order flux corrections,

$`\mathrm{\Delta }F_{n+1/2}^{L,t}`$ $`={\displaystyle \frac{F_{n+1}^tF_n^t}{2}},`$ (36)
$`\mathrm{\Delta }F_{n+1/2}^{R,t}`$ $`={\displaystyle \frac{F_{n+2}^tF_{n+1}^t}{2}}.`$ (37)
are based on cell $`n+1`$ and the cells directly adjacent to it. Near extrema where the corrections have opposite signs, we impose no second-order correction and the flux assignment scheme reduces to first-order. A flux limiter $`\varphi `$ is then used to determine the appropriate second-order correction,

$$\mathrm{\Delta }F_{n+1/2}^t=\varphi (\mathrm{\Delta }F_{n+1/2}^{L,t},\mathrm{\Delta }F_{n+1/2}^{R,t}),$$
(38)
which still maintains the TVD condition. The second-order correction is added to the first-order fluxes to get second-order fluxes. The first-order upwind scheme and second-order TVD scheme will be referred to as monotone upwind schemes for conservation laws (MUSCL).

Time integration is performed using a second-order Runge-Kutta scheme. We first do a half time step,

$$u_n^{t+\mathrm{\Delta }t/2}=u_n^t\left(\frac{F_{n+1/2}^tF_{n1/2}^t}{\mathrm{\Delta }x}\right)\frac{\mathrm{\Delta }t}{2},$$
(39)
using the first-order upwind scheme to obtain the half-step values $`u^{t+\mathrm{\Delta }t/2}`$. A full time step is then computed,

$$u_n^{t+\mathrm{\Delta }t}=u_n^t\left(\frac{F_{n+1/2}^{t+\mathrm{\Delta }t/2}F_{n1/2}^{t+\mathrm{\Delta }t/2}}{\mathrm{\Delta }x}\right)\mathrm{\Delta }t,$$
(40)
using the TVD scheme on the half-step fluxes $`F^{t+\mathrm{\Delta }t/2}`$. The reader is encouraged to show that is second-order accurate.

We briefly discuss three TVD limiters. The minmod flux limiter chooses the smallest absolute value between the left and right corrections:

$$minmod(a,b)=\frac{1}{2}[\text{sign}(a)+\text{sign}(b)]\mathrm{min}(|a|,|b|).$$
(41)
The superbee limiter (Roe, 1985) chooses between the larger correction and two times the smaller correction, whichever is smaller in magnitude:

$$superbee(a,b)=\{\begin{array}{cc}minmod(a,2b)\hfill & \text{if }|a||b|,\hfill \\ minmod(2a,b)\hfill & \text{otherwise .}\hfill \end{array}$$
(42)
The Van Leer limiter (Van Leer, 1974) takes the harmonic mean of the left and right corrections:

$$vanleer(a,b)=\frac{2ab}{a+b}.$$
(43)
The minmod limiter is the most moderate of all second-order TVD limiters. In Figure (5) we see that it does not do much better than first-order upwind for the square wave advection test. Superbee chooses the maximum correction allowed under the TVD constraint. It is especially suited for piece-wise linear conditions and is the least diffusive for this particular test, as can be seen in Figure (6). Note that no additional diffusion can be seen by advecting the square wave more than once through the box. It can be shown that the minmod and superbee limiters are extreme cases which bound all other second-order TVD limiters. The Van Leer limiter differs from the previous two in that it is analytic. This symmetrical approach falls somewhere inbetween the other two limiters in terms of moderation and diffusion, as can be seen in Figure (7). It can be shown that the CFL condition for the second-order TVD scheme is to have $`\lambda <1`$. For a comprehensive discussion on TVD limiters, see Hirsch (1990) and Laney (1998).

## 5 Relaxing TVD

We now describe a simple and robust method to solve the Euler equations using the monotone upwind scheme for conservation laws (MUSCL) from the previous section. The relaxing TVD method (Jin & Xin, 1995) provides high resolution capturing of shocks using computationally inexpensive algorithms which are straightforward to implement and to parallelize. It has been successfully implemented for simulating cosmological astrophysical fluids by Pen (1998).

The MUSCL scheme is straightforward to apply to conservation laws like the advection equation since the velocity alone can be used as a marker of the direction of flow. However, applying the MUSCL scheme to solve the Euler equations is made difficult by the fact that the momentum and energy fluxes depend on the pressure. In order to determine the direction upwind of the flow, it becomes necessary to calculate the flux Jacobian eigenvectors using Riemann solvers. This step requires computationally expensive algorithms. The relaxing TVD method offers an attractive alternative.

### 5.1 1-D Scalar Conservation Law

We first present a motivation for the relaxing scheme by again considering the 1-D scalar conservation law. The MUSCL scheme for solving the linear advection equation is explicitly asymmetric in that it depends on the sign of the advection velocity. We now describe a symmetrical approach which applies to a general advection velocity.

The flow can be considered as a sum of a right-moving wave $`u^R`$ and a left-moving wave $`u^L`$. For a positive advection velocity, the amplitude of the left-moving wave is zero and for a negative advection velocity, the amplitude of the right-moving wave is zero. In compact notation, the waves can be defined as:

$`u^R`$ $`=\left({\displaystyle \frac{1+v/c}{2}}\right)u,`$ (44)
$`u^L`$ $`=\left({\displaystyle \frac{1v/c}{2}}\right)u,`$ (45)
where $`c=|v|`$. The two waves are traveling in opposite directions with advection speed $`c`$ and can be described by the advection equations:

$$\frac{u^R}{t}+\frac{}{x}(cu^R)=0,$$
(46)
$$\frac{u^L}{t}\frac{}{x}(cu^L)=0.$$
(47)
The MUSCL scheme is straightforward to apply to solve Equations (46) and (47) since the upwind direction is left for the right-moving wave and right for the left-moving wave. The 1-D relaxing advection equation then becomes

$$\frac{u}{t}+\frac{F^R}{x}\frac{F^L}{x}=0.$$
(48)
where $`F^R=cu^R`$ and $`F^L=cu^L`$. For the discretized solution given by Equation (11), the boundary fluxes $`F_{n+1/2}^t`$ are now a sum of the fluxes $`F_{n+1/2}^{R,t}`$ and $`F_{n+1/2}^{L,t}`$ from the right-moving and left-moving waves, respectively. Note that the relaxing advection equation will correctly reduce to the linear advection equation for any general advection velocity.

Using this symmetrical approach, a general algorithm can be written to solve the linear advection equation for an arbitrary advection velocity. This scheme is indeed inefficient for solving the linear advection equation since one wave will have zero amplitude. However, the Euler equations can have both right-moving and left-moving waves with non-zero amplitudes.

### 5.2 1-D Systems of Conservation Laws

We now discuss the 1-D relaxing TVD scheme and later generalize it to higher spatial dimensions. Consider a 1-D system of conservation laws,

$$\frac{𝒖}{t}+\frac{𝑭(u)}{x}=0,$$
(49)
where for the Euler equations, we have $`𝒖=(\rho ,\rho v,e)`$ and $`𝑭(𝒖)`$ the corresponding flux terms. We now replace the vector conservation law with the relaxation system

$$\frac{𝒖}{t}+\frac{}{x}(c𝒘)=0,$$
(50)
$$\frac{𝒘}{t}+\frac{}{x}(c𝒖)=0,$$
(51)
where $`c(x,t)`$ is a free positive function called the freezing speed. The relaxation system contains two coupled vector linear advection equations. In practice, we set $`𝒘=𝑭(𝒖)/c`$ and use it as an auxiliary vector to calculate fluxes. Equation (50) reduces to our 1-D vector conservation law and Equation (51) is a vector conservation law for $`𝒘`$.

In order to solve the relaxed system, we decouple the equations through a change of variables:

$`𝒖^R`$ $`={\displaystyle \frac{𝒖+𝒘}{2}},`$ (52)
$`𝒖^L`$ $`={\displaystyle \frac{𝒖𝒘}{2}},`$ (53)
which then gives us

$$\frac{𝒖^R}{t}+\frac{}{x}(c𝒖^R)=0,$$
(54)
$$\frac{𝒖^L}{t}\frac{}{x}(c𝒖^L)=0.$$
(55)
Equations (54) and (55) are vector linear advection equations, which can be interpreted as right-moving and left-moving flows with advection speed $`c`$. Note the similarity with their scalar counterparts, Equations (46) and (47). The 1-D vector relaxing conservation law for $`𝒖`$ becomes

$$\frac{𝒖}{t}+\frac{𝑭^R}{x}\frac{𝑭^L}{x}=0,$$
(56)
where $`𝑭^R=c𝒖^R`$ and $`𝑭^L=c𝒖^L`$. The vector relaxing equation can now be solved by applying the MUSCL scheme to Equations (54) and (55). Again, note the similarity between the vector relaxing equation and its scalar counterpart, Equation (48).

The relaxed scheme is TVD under the constraint that the freezing speed $`c`$ be greater than the characteristic speed given by the largest eigenvalue of the flux Jacobian $`𝑭(𝒖)/𝒖`$. For the Euler equations, this is satisfied for

$$c=|v|+c_s.$$
(57)
Jin & Xin (1995) considered the freezing speed to be a positive constant in their relaxing scheme while we generalize it to be a positive function. Time integration is again performed using a second-order Runge-Kutta scheme and the time step is determined by satisfying the CFL condition,

$$\frac{c_{\text{max}}\mathrm{\Delta }t}{\mathrm{\Delta }x}1.$$
(58)
Note that a new freezing speed is computed for each partial step in the Runge-Kutta scheme. The CFL number $`\lambda =c_{max}\mathrm{\Delta }t/\mathrm{\Delta }x`$ should be chosen such that $`c_{\text{max}}`$ will be larger than $`\mathrm{max}(c_n^t)`$ and $`\mathrm{max}(c_n^{t+\mathrm{\Delta }t/2})`$.

We now summarize the steps needed to numerically solve the 1-d Euler equations. At the beginning of each partial step in the Runge-Kutta time integration scheme, we need to calculate the cell-averaged variables defined at grid cell centres. First for the half time step, we calculate the fluxes $`𝑭(u_n^t)`$ and the freezing speed $`c_n^t`$. We then set the auxiliary vector $`𝒘_n^t=𝑭(𝒖_n^t)/c_n^t`$ and construct the right-moving waves $`𝒖_n^{R,t}`$ and left-moving waves $`𝒖_n^{L,t}`$. The half time step is given by

$$𝒖_n^{t+\mathrm{\Delta }t/2}=𝒖_n^t\left(\frac{𝑭_{n+1/2}^t𝑭_{n1/2}^t}{\mathrm{\Delta }x}\right)\frac{\mathrm{\Delta }t}{2},$$
(59)
where

$$𝑭_{n+1/2}^t=𝑭_{n+1/2}^{R,t}𝑭_{n+1/2}^{L,t}.$$
(60)
The first-order upwind scheme is used to compute fluxes at cell boundaries for the right-moving and left-moving waves. For the full time step, we construct the right-moving waves $`𝒖_n^{R,t+\mathrm{\Delta }t/2}`$ and left-moving waves $`𝒖_n^{L,t+\mathrm{\Delta }t/2}`$, using the half-step values of the appropriate variables. The full time step,

$$𝒖_n^{t+\mathrm{\Delta }t}=𝒖_n^t\left(\frac{𝑭_{n+1/2}^{t+\mathrm{\Delta }t/2}𝑭_{n1/2}^{t+\mathrm{\Delta }t/2}}{\mathrm{\Delta }x}\right)\mathrm{\Delta }t,$$
(61)
is computed using the second-order TVD scheme. This completes the updating of $`𝒖^t`$ to $`𝒖^{t+\mathrm{\Delta }t}`$.

We have found that a minor modification to the implementation described above gives more accurate results. Consider writing the flux of the right-moving and left-moving waves as:

$`𝑭^R=c𝑮^R,`$ (62)
$`𝑭^L=c𝑮^L,`$ (63)
where $`𝑮^R`$ is the flux of $`𝝁^R=𝒖^R/c`$ and $`𝑮^L`$ is the flux of $`𝝁^L=𝒖^L/c`$. The linear advection equations for $`𝝁^R`$ and $`𝝁^L`$ are similar to Equations (54) and (55), but where we replace $`𝒖^R`$ with $`𝝁^R`$ and $`𝒖^L`$ with $`𝝁^L`$. For each partial step in the Runge-Kutta scheme, the net fluxes at cell boundaries are then taken to be

$$𝑭_{n+1/2}=c_{n+1/2}(𝑮_{n+1/2}^R𝑮_{n+1/2}^L),$$
(64)
where we use $`c_{n+1/2}=(c_{n+1}+c_n)/2`$. In practice, this modified implementation has been found to resolve shocks with better accuracy in certain cases. Note that the two different implementations of the relaxing TVD scheme are identical when a constant freezing speed is used.

### 5.3 Multi-Dimensional Systems of Conservation Laws

The 1-D relaxing TVD scheme can be generalized to higher dimensions using the dimensional splitting technique by Strang (1968). In three dimensions, the Euler equations can be dimensionally split into three separate 1-D equations which are solved sequentially. Let the operator $`L_i`$ represent the updating of $`u^t`$ to $`u^{t+\mathrm{\Delta }t}`$ by including the flux in the $`i`$ direction. We first complete a forward sweep,

$$u^{t+\mathrm{\Delta }t}=L_zL_yL_xu^t,$$
(65)
and then perform a reverse sweep

$$u^{t+2\mathrm{\Delta }t}=L_xL_yL_zu^{t+\mathrm{\Delta }t},$$
(66)
using the same time step $`\mathrm{\Delta }t`$ to obtain second-order accuracy. We will refer to the combination of the forward and reverse sweeps as a double sweep.

A more symmetrical sweeping pattern can be used by permutating the sweeping order when completing the next double time step. The dimensional splitting or operator splitting technique can be summarized as follows:

$`u^{t_2}`$ $`=u^{t_1+2\mathrm{\Delta }t_1}=L_xL_yL_zL_zL_yL_xu^{t_1},`$ (67)
$`u^{t_3}`$ $`=u^{t_2+2\mathrm{\Delta }t_2}=L_zL_xL_yL_yL_xL_zu^{t_2},`$ (68)
$`u^{t_4}`$ $`=u^{t_3+2\mathrm{\Delta }t_3}=L_yL_zL_xL_xL_zL_yu^{t_3},`$ (69)
where $`\mathrm{\Delta }t_1`$, $`\mathrm{\Delta }t_2`$, and $`\mathrm{\Delta }t_3`$ are newly determined time steps after completing each double sweep.

The CFL condition for the 3-D relaxing TVD scheme is similarly given by Equation (58), but with

$$c_{\mathrm{max}}=\mathrm{max}[(c_x)_{\mathrm{max}},(c_y)_{\mathrm{max}},(c_z)_{\mathrm{max}}].$$
(70)
where $`c_i=|v_i|+c_s`$. Note that since $`\mathrm{max}(|v_i|)`$ is on average a factor of $`\sqrt{3}`$ smaller than $`\mathrm{max}(|𝒗|)`$, a dimensionally split scheme can use a longer time step compared to an un-split scheme.

The dimensional splitting technique also has other advantages. The decomposition into a 1-D problem allows one to write short 1-D algorithms, which are easy to optimize to be cache efficient. A 3-D hydro code is straightforward to implement in parallel. When sweeping in the x direction, for example, one can break up the data into 1-D columns and operate on the independent columns in parallel. A sample 3-D relaxing TVD code, implemented in parallel using OpenMP directives, is provided in the appendix.

## 6 Sedov-Taylor Blast Wave Test for 3-D Hydro

A rigourous and challenging test for any 3-D Eulerian or Lagrangian hydrodynamic code is the Sedov-Taylor blast wave test. We set up the simulation box with a homogeneous medium of density $`\rho _1`$ and negligible pressure and introduce a point-like supply of thermal energy $`E_{}`$ in the centre of the box at time $`t=0`$. The challenge is to accurately capture the strong spherical shock wave which sweeps along material as it propagates out into the ambient medium. The Sedov-Taylor test is used to model nuclear-type explosions. In astrophysics, it is often used as a basic setup to model supernova explosions and the evolution of supernova remnants (See Shu, 1992).

The analytical Sedov solution uses the self-similar nature of the blast wave expansion (See Landau & Lifshitz, 1987). Consider a frame fixed relative to the centre of the explosion. The spherical shock front propagates outward and the distance from the origin is given by

$$r_{sh}(t)=\xi _{}\left(\frac{E_{}t^2}{\rho _1}\right)^{1/5},$$
(71)
where $`\xi _{}=1.15`$ for an ideal gas with $`\gamma =5/3`$. The velocity of the shock $`v_{sh}=r_{sh}/t`$ is given by

$$v_{sh}(t)=\frac{2}{5}\frac{r_{sh}(t)}{t}.$$
(72)
Since the ambient medium has negligible pressure, the shocks will be very strong. The density $`\rho _2`$, velocity $`v_2`$, and pressure $`P_2`$ directly behind the shock front are:

$`\rho _2`$ $`=\left({\displaystyle \frac{\gamma +1}{\gamma 1}}\right)\rho _1,`$ (73)
$`v_2`$ $`=\left({\displaystyle \frac{2}{\gamma +1}}\right)v_{sh},`$ (74)
$`P_2`$ $`=\left({\displaystyle \frac{2}{\gamma +1}}\right)\rho _1v_{sh}^2.`$ (75)
The immediate post-shock gas density is constant in time, while the shocked gas velocity $`v_2`$ and pressure $`P_2`$ decrease as $`t^{3/5}`$ and $`t^{6/5}`$, respectively. The full analytical Sedov-Taylor solutions can be found in Landau & Lifshitz (1987).

The 3-D relaxing TVD code based on the van Leer flux limiter is applied to capturing the Sedov-Taylor blast wave. We set up a box with $`256^3`$ cells and constant initial density $`\rho _1=1`$. At time $`t=0`$, we inject a supply of thermal energy $`E_{}=10^5`$ into one cell at the centre of the box. The simulation is stopped at time $`t=283`$, in which the shock front has propagated out to a distance of $`r_{sh}=110`$ cells from the centre. In Figure (8) and (9) we plot the radial distributions of density, momentum, and pressure, normalized to $`\rho _2`$, $`\rho _2v_2`$, and $`P_2`$ respectively. The data points are taken from a random subset of cells and the solid lines are the analytical Sedov-Taylor solutions. Despite solving a spherically symmetric problem on an explicitly non-rotationally invariant Cartesian grid, the anisotropic scatter in the results is small. The distance of the shock front from the centre of the explosion as a function of time is indeed given by Equation (71), demonstrating that the 3-D relaxing TVD code ensures the correct shock propagation. The resolution of the shock front is roughly two grid cells. The numerical shock jump values of $`\rho _2`$, $`v_2`$, and $`P_2`$ are resolution dependent and come close to the theoretical values for our test with $`256^3`$ cells. We leave it as an exercise for the reader to test the code using the minmod and superbee flux limiters.

## 7 Self-Gravitating Hydro for Astrophysical Applications

For astrophysical applications, both hydrodynamical and gravitational forces are included. The gravitational forces arise from the self-gravity of the fluid and can also come from an external field. The Euler equations with the gravitational source terms included are given as:

$$\frac{\rho }{t}+\frac{}{x_j}(\rho v_j)=0,$$
(76)
$$\frac{(\rho v_i)}{t}+\frac{}{x_j}(\rho v_iv_j+P\delta _{ij})=\rho \frac{\varphi }{x_i},$$
(77)
$$\frac{e}{t}+\frac{}{x_j}[(e+P)v_j]=\rho v_i\frac{\varphi }{x_i}.$$
(78)
where $`\varphi `$ is the gravitational potential. Poisson’s equation,

$$^2\varphi =4\pi G\rho ,$$
(79)
relates the gravitational potential to the density field. The general solution can be written as

$$\varphi (𝒙)=\rho (𝒙^{})w(𝒙𝒙^{})d^3x^{},$$
(80)
where the kernel is given by

$$w(𝒙)=\frac{G}{|𝒙|}.$$
(81)
In the discrete case, the integral in Equation (80) becomes a sum and Poisson’s equation can be solved using fast Fourier transforms (FFT) to do the convolution. The forces are then calculated by finite differencing the potential (See Hockney & Eastwood, 1988).

The addition of gravatitational source terms in the Euler equations is easily handled using the operator splitting technique described in $`\mathrm{\S }`$5.3. Consider the double sweep:

$$u^{t+2\mathrm{\Delta }t}=L_xL_yL_zGGL_zL_yL_xu^t,$$
(82)
where the operator $`L_i`$ represents the updating of $`u`$ by including the flux in the $`i`$ direction and the operator $`G`$ represents the gravitational acceleration of the fluid. During the gravitational step, the flux terms in the Euler equations are ignored. The density distribution does not change and only the fluid momenta and total energy density are updated.

### 7.1 Astrophysical Formation of Blue Stragglers Through Stellar Collisions

The stellar density in the cores of globular and open clusters is high enough for stellar collisions to take place with significant frequency (Hills & Day, 1976). Current observations and simulations suggest that the merger of two main sequence stars produces a blue straggler (Sills et al., 1997; Sandquist, Bolte, & Hernquist, 1997). The blue stragglers are out-lying main sequence stars which lie beyond the main sequence turnoff in the colour-magnitude diagram (CMD) of a star cluster. The blue stragglers are more massive, brighter, and bluer than the turnoff stars. Since more massive stars evolve faster than lower mass stars and are not expected to lie beyond the turnoff, this suggests that blue stragglers must have formed more recently.

In principle the merger of two main sequence stars can produce a young remnant star provided that significant mixing occurs in the process. The mixing produces a higher hydrogen fraction in the core of the remnant than that of the parent stars which have already burnt most of the hydrogen to helium in their cores. Benz & Hills (1987) used low resolution SPH simulations with $`10^3`$ particles to simulate the merging of $`n=3/2`$ polytropes and found that they fully mixed. However, medium resolution SPH simulations with $`10^4`$ particles of $`n=3/2`$ or $`n=3`$ polytropes showed only weak mixing (Lombardi, Rasio, & Shapiro, 1996; Sandquist, Bolte, & Hernquist, 1997). It is worth noting that $`n=3/2`$ polytropes are more representative of low mass main sequence stars with large convective envelopes while $`n=3`$ polytropes resemble main sequence stars near the turnoff which have little mass in their convective envelopes. High resolution SPH simulations involving $`10^510^6`$ particles have now been applied to simulating stellar collisions (Sills et a., 2002).

The merging stars process is mostly subsonic and strong shocks are not expected. In the absence of shocks, SPH particles will follow flow lines of constant entropy due to the Lagrangian nature of the method. As a result, the particles may experience sedimentation. In addition, the mixing can also depend on the adopted smoothing length and the form of artificial viscosity. For a SPH fluid, the Reynolds number is of order $`(N_p/N_s)^{1/3}`$, where $`N_p`$ is the total number of particles and $`N_s`$ is the number of particles over which the smoothing is done. For $`N_p10^5`$ and $`N_s10^2`$, the Reynolds number is $`10`$. However, a fluid with a low Reynolds number will tend to experience laminar flow. Hence, SPH may under mix.

It is a worthwhile exercise to model the merging process using Eulerian hydrodynamical simulations. The differences between Eulerian and Lagrangian approaches may lead to very different results on mixing. As of present, no such work has been reported in the literature.

### 7.2 Numerical Method

We consider the off-axis collision of two main sequence stars with $`M=0.8M_{\mathrm{}}`$ and $`R=0.955R_{\mathrm{}}`$, which are modeled using $`n=3`$ polytropes. A polytrope with polytropic index $`n`$ has equilibrium density and pressure profiles which are related by

$$P\rho ^{1+1/n}.$$
(83)
The density profile is determined by solving the Lane-Emden equation (See Chandrasekhar, 1957). We adopt an ideal gas equation of state with adiabatic index $`\gamma =5/3`$. Note that for an $`n=3`$ polytrope, 90% of the total mass is contained within $`r0.5R`$. We define the dynamical time to be,

$$\tau _{dyn}\frac{1}{\sqrt{G\overline{\rho }}},$$
(84)
where $`\overline{\rho }`$ is the average density. For the chosen parent stars, the dynamical time is approximately one physical hour.

The collision is simulated in a box with $`512\times 512\times 256`$ cells and the orbital plane coincides with the $`xy`$ plane. Initially, each parent star has a radius of 96 grid cells. The stars are set up on zero-energy parabolic orbits with a pericentre separation equal to $`0.25R`$. The initial trajectories are calculated assuming point masses. In an Eulerian simulation, the vacuum cannot have zero density. We set the minimum density of the cells to be $`10^8`$ of the central density of the parent stars. The hydrodynamics is done in a non-periodic box with vacuum boundary conditions.

### 7.3 Numerical Results

A non-trivial test of a self-gravitating Eulerian hydro code is the advection of an object in hydrostatic equilibrium. The challenge is to maintain the equilibrium profile over a large number of time steps. One of the parent stars is placed in a periodic box with $`256^3`$ cells and given some initial momentum. We make the test rigorous by having the polytrope move in all three directions. In Figure 10 we compare the mass and entropy profiles of the initial and advected polytrope. The entropic variable $`AP/\rho ^\gamma `$ is used in place of the specific entropy. The parameter $`A_{}`$ is defined to be the minimum entropy of the parent polytrope. After 1000 timesteps in which the polytrope has moved 256 cells in each direction, the advected polytrope has still retained its equilibrium profile. Shock heating can occur in the outer envelope as the polytrope moves through the false vacuum. However, by setting the density of the false vacuum to be $`10^8`$ of the central density of the polytrope, we can minimize the spurious shock heating.

In Figure 11 we show four snapshots of the merging process taken at time $`t=0`$, 2, 4, and $`8\tau _{dyn}`$. The 2-D density maps are created by averaging over 4 planes taken about the orbital mid-plane. The density contours are spaced logarithmically with 2 per decade and covering 3 decades down from the maximum. The parent stars are initially separated by $`3.75R`$ and placed on zero-energy orbits with a pericentre separation of $`0.25R`$. During the collision process, the outer envelopes of the parent stars are shock heated and material gets ejected. In less than $`10\tau _{dyn}`$, the merger remnant establishes hydrostatic equilibrium. The merger remnant is a rotating oblate with mass approximately 90% of the combined mass of the parent stars. A large fraction of the mass loss is due to the vacuum boundary conditions. Ejected material do not have the opportunity to fall back onto the merger remnant. However, the additional mass loss in the envelope does not present a problem since we are interested in the question of mixing in the interior of the star.

In Figure 12 we plot the thermodynamic profiles of the merger remnant and the parent stars. The central density and pressure in the core of the merger remnant is lower than the corresponding values in the parent stars by approximately half. The entropy floor has risen by a factor of 1.6. Shock heating is expected to be minimal in the core so a change in entropy suggests that some mixing has taken place. However, it is difficult to quantify the amount of mixing from examining the thermodynamic profiles alone.

### 7.4 Future Work

To help address the question of mixing, we are implementing a particle-mesh (PM) scheme where test particles can be used to track passively advected quantities such as chemical composition. Initially, each parent star is assigned a large number of particles with known chemical composition. The test particles are passively advected along velocity field lines. For each time step, the velocity of each particle is interpolated from the grid using a ”cloud-in-cell” (CIC) scheme (Hockney & Eastwood, 1988) and the equations of motions are solved using second-order Runge-Kutta integration. The CIC interpolation scheme is also used to determine the local density, pressure, and entropy associated with each particle. With this setup, we have the benefit of being able to track thermodynamic quantities like in an SPH scheme but avoid the under mixing problem since the fluid equations are solved using the Eulerian scheme.

Future work (Trac, Sills, & Pen 2003) will have higher resolution simulations. Collisions will be simulated in a box with $`1024\times 1024\times 512`$ cells. Each parent star will have a radius of 192 grid cells and be assigned $`256^3`$ test particles. We will also be doing a detailed comparison between Eulerian and SPH simulations of stellar mergers.

The self-gravitating hydro code used for the simulations is very memory friendly. For the $`1024^2\times 512`$ grid, 10 GB is required to store the hydro variables, 2 GB for the potential, and less than 1 GB for the test particles. For every double timestep, approximately 1000 floating point operations per grid cell is needed to carry out the TVD hydro calculations. The potential is computed once for every double step and this requires two FFTs. Since Eulerian codes are very memory friendly, have low floating point counts, are easily parallelized, and scale very well on shared-memory, multiple-processor machines, they can be used to run very high resolution simulations.

## 8 Summary

We have presented several numerical schemes for solving the linear advection equation and given the CFL stability conditions for each scheme. We have implemented the relaxing TVD scheme to solve the Euler system of conservation laws. The second-order accurate TVD scheme provides high resolution capturing of shocks, as can be seen in the Riemann shock test and the Sedov-Taylor blast wave test. The 1-D relaxing TVD scheme can be easily generalized to higher dimensions using the dimensional splitting technique. A dimensionally split scheme can use longer time steps and is straightforward to implement in parallel. We have presented a sample astrophysical application. A 3-D self-gravitating Eulerian hydro code is used to simulate the formation of blue straggler stars through stellar mergers. We hope to have convinced the reader that Eulerian computational fluid dynamics is a powerful approach to simulating complex fluid flows because it is simple, fast, and accurate.

We thank Joachim Stadel and Norm Murray for comments and suggestions on the writing and editing of this paper. We also thank Alison Sills, Phil Arras, and Chris Matzner for discussions on stellar mergers.

## Appendix A 3-D Relaxing TVD code

We provide a sample 3-D relaxing TVD code written in Fortran 90. The code is implemented using OpenMP directives to run in parallel on shared memory machines. The code is fast and memory friendly. The array u(a,i,j,k) stores the five conserved hydro quantities $`𝖺=(\rho ,\rho v_x,\rho v_y,\rho v_z,e)`$ for each cell (i,j,k) in the Cartesian cubical lattice with side length nc. For each sweep, we first call the subroutine timestep to determine the appropriate time step dt which satisfies the CFL condition. The updating of u by including the flux in the x direction is performed by the sweepx subroutine. The data array u is divided into 1-D array sections u1d(a,i) which are operated on by the relaxing subroutine. The independent columns are distributed amongst multiple processors on a shared memory machine by the OpenMP directives.

The relaxing TVD subroutine in this sample code is written for ease of readability and therefore, is not fully optimized. At the beginning of each partial step in the Runge-Kutta time integration scheme, the cell-averaged variables defined at grid cell centres are calculated by the averageflux subroutine. The fluxes at cell boundaries for the right-moving and left-moving waves are stored in fr and fl, respectively. We have implemented the minmod, superbee, and Van Leer flux limiters and the user of the code can easily switch between them.

We have provided some initial conditions for the Sedov-Taylor blast wave test. The reader is encouraged to test the code and compare how the various flux limiters do at resolving strong shocks. This sample code does not implement the modified relaxing TVD scheme described at the end of $`\mathrm{\S }\text{5.2}`$, which has been found work very well with the Van Leer flux limiter but unstable with superbee for the 3-D Sedov Taylor test. We have found that the superbee limiter is often unstable for 3-D fluid simulations. Please contact the authors regarding any questions on the implementation of the relaxing TVD algorithm.

```

program main
  implicit none
  integer, parameter :: nc=64,hc=nc/2
  real, parameter :: gamma=5./3,cfl=0.9

  real, dimension(5,nc,nc,nc) :: u

  integer nsw,stopsim
  real t,tf,dt,E0,rmax

  t=0
  dt=0
  nsw=0
  stopsim=0

  E0=1e5
  rmax=3*hc/4
  tf=sqrt((rmax/1.15)**5/E0)

  call sedovtaylor
  do
     call timestep
     call sweepx
     call sweepy
     call sweepz
     call sweepz
     call sweepy
     call sweepx
     if (stopsim .eq. 1) exit
     call timestep
     call sweepz
     call sweepx
     call sweepy
     call sweepy
     call sweepx
     call sweepz
     if (stopsim .eq. 1) exit
     call timestep
     call sweepy
     call sweepz
     call sweepx
     call sweepx
     call sweepz
     call sweepy
     if (stopsim .eq. 1) exit
  enddo
  call outputresults

contains

  subroutine sedovtaylor
    implicit none
    integer i,j,k

    do k=1,nc
       do j=1,nc
          do i=1,nc
             u(1,i,j,k)=1
             u(2:4,i,j,k)=0
             u(5,i,j,k)=1e-3
          enddo
       enddo
    enddo
    u(5,hc,hc,hc)=E0
    return
  end subroutine sedovtaylor

  subroutine outputresults
    implicit none
    integer i,j,k
    real r,x,y,z

    open(1,file=’sedovtaylor.dat’,recl=200)
    do k=1,nc
       z=k-hc
       do j=1,nc
          y=j-hc
          do i=1,nc
             x=i-hc
             r=sqrt(x**2+y**2+z**2)
             write(1,*) r,u(:,i,j,k)
          enddo
       enddo
    enddo
    close(1)
    return
  end subroutine outputresults

  subroutine timestep
    implicit none
    integer i,j,k
    real P,cs,cmax
    real v(3)

    cmax=1e-5
    !$omp parallel do default(shared) private(i,j,k,v,cs,P) reduction(max:cmax)
    do k=1,nc
       do j=1,nc
          do i=1,nc
             v=abs(u(2:4,i,j,k)/u(1,i,j,k))
             P=max((gamma-1)*(u(5,i,j,k)-u(1,i,j,k)*sum(v**2)/2),0.)
             cs=sqrt(gamma*P/u(1,i,j,k))
             cmax=max(cmax,maxval(v)+cs)
          enddo
       enddo
    enddo
    !$omp end parallel do

    dt=cfl/cmax
    if (t+2*dt .gt. tf) then
       dt=(tf-t)/2
       stopsim=1
    endif
    t=t+2*dt
    nsw=nsw+1
    write(*,”(a7,i3,a8,f7.5,a6,f8.5)”) ’nsw = ’,nsw,’   dt = ’,dt,’   t = ’,t
    return
  end subroutine timestep

  subroutine sweepx
    implicit none
    integer j,k
    real u1d(5,nc)

    !$omp parallel do default(shared) private(j,k,u1d)
    do k=1,nc
       do j=1,nc
          u1d=u(:,:,j,k)
          call relaxing(u1d)
          u(:,:,j,k)=u1d
       enddo
    enddo
    !$omp end parallel do
    return
  end subroutine sweepx

  subroutine sweepy
    implicit none
    integer i,k
    real u1d(5,nc)

    !$omp parallel do default(shared) private(i,k,u1d)
    do k=1,nc
       do i=1,nc
          u1d((/1,3,2,4,5/),:)=u(:,i,:,k)
          call relaxing(u1d)
          u(:,i,:,k)=u1d((/1,3,2,4,5/),:)
       enddo
    enddo
    !$omp end parallel do
    return
  end subroutine sweepy

  subroutine sweepz
    implicit none
    integer i,j
    real u1d(5,nc)

    !$omp parallel do default(shared) private(i,j,u1d)
    do j=1,nc
       do i=1,nc
          u1d((/1,4,3,2,5/),:)=u(:,i,j,:)
          call relaxing(u1d)
          u(:,i,j,:)=u1d((/1,4,3,2,5/),:)
       enddo
    enddo
    !$omp end parallel do
    return
  end subroutine sweepz

  subroutine relaxing(u)
    implicit none
    real, dimension(nc) :: c
    real, dimension(5,nc) :: u,u1,w,fu,fr,fl,dfl,dfr

    !! Do half step using first-order upwind scheme
    call averageflux(u,w,c)
    fr=(u*spread(c,1,5)+w)/2
    fl=cshift(u*spread(c,1,5)-w,1,2)/2
    fu=(fr-fl)
    u1=u-(fu-cshift(fu,-1,2))*dt/2

    !! Do full step using second-order TVD scheme
    call averageflux(u1,w,c)

    !! Right-moving waves
    fr=(u1*spread(c,1,5)+w)/2
    dfl=(fr-cshift(fr,-1,2))/2
    dfr=cshift(dfl,1,2)
    call vanleer(fr,dfl,dfr)
    !call minmod(fr,dfl,dfr)
    !call superbee(fr,dfl,dfr)

    !! Left-moving waves
    fl=cshift(u1*spread(c,1,5)-w,1,2)/2
    dfl=(cshift(fl,-1,2)-fl)/2
    dfr=cshift(dfl,1,2)
    call vanleer(fl,dfl,dfr)
    !call minmod(fl,dfl,dfr)
    !call superbee(fl,dfl,dfr)

    fu=(fr-fl)
    u=u-(fu-cshift(fu,-1,2))*dt
    return
  end subroutine relaxing

  subroutine averageflux(u,w,c)
    implicit none
    integer i
    real P,v
    real u(5,nc),w(5,nc),c(nc)

    !! Calculate cell-centered fluxes and freezing speed
    do i=1,nc
       v=u(2,i)/u(1,i)
       P=max((gamma-1)*(u(5,i)-sum(u(2:4,i)**2)/u(1,i)/2),0.)
       c(i)=abs(v)+max(sqrt(gamma*P/u(1,i)),1e-5)
       w(1,i)=u(1,i)*v
       w(2,i)=(u(2,i)*v+P)
       w(3,i)=u(3,i)*v
       w(4,i)=u(4,i)*v
       w(5,i)=(u(5,i)+P)*v
    enddo
    return
  end subroutine averageflux

  subroutine vanleer(f,a,b)
    implicit none
    real, dimension(5,nc) ::  f,a,b,c

    c=a*b
    where (c .gt. 0)
       f=f+2*c/(a+b)
    endwhere
    return
  end subroutine vanleer

  subroutine minmod(f,a,b)
    implicit none
    real, dimension(nc) :: f,a,b

    f=f+(sign(1.,a)+sign(1.,b))*min(abs(a),abs(b))/2.
    return
  end subroutine minmod

  subroutine superbee(f,a,b)
    implicit none
    real, dimension(5,nc) :: f,a,b

    where (abs(a) .gt. abs(b))
       f=f+(sign(1.,a)+sign(1.,b))*min(abs(a),abs(2*b))/2.
    elsewhere
       f=f+(sign(1.,a)+sign(1.,b))*min(abs(2*a),abs(b))/2.
    endwhere
    return
  end subroutine superbee

end program main

```


format: text

text:
# 1. Introduction

## 1. Introduction

Hess and Philipp  (also see Hess and Philipp ) have constructed what, they claim, is a local hidden-variables model reproducing the quantum mechanical predictions for the singlet state—contrary to the result proved by Bell . The falsity of their claim has been shown by, among others, Myrvold  (also see Gill *et al*  and Mermin ). However, Hess and Philipp  dispute Myrvold’s conclusion. The purpose of this note is to present some additional considerations which, we hope, may help to settle the question.

Hess and Philipp desribe their model in two different places. The two versions are not identical. By some oversight the published version  (version 1) omits a crucial detail. As a result it disagrees with quantum mechanics. It also violates *signal* locality. The unpublished version  of their model (version 2) does not suffer from these deficiencies. However, the same equations which show that version 1 violates signal locality also conclusively demonstrate that version 2 violates parameter independence, as Myrvold has argued.

The concept of parameter independence was analyzed in detail by Jarrett  and Shimony . Consider a pair of spin-$`1/2`$ particles prepared in the singlet state. Particle $`1`$ is sent to station $`S_1`$, which measures its spin in the direction $`𝐚`$, obtaining the outcome $`A=\pm 1`$. Particle $`2`$ is sent to station $`S_2`$, which measures its spin in the direction $`𝐛`$, obtaining the outcome $`B=\pm 1`$. Let $`\lambda `$ denote the complete hidden state of the pair $`1+2`$. Let $`p_1(A|𝐚,𝐛,\lambda )`$ be the probability of obtaining outcome $`A`$ at $`S_1`$ for detector settings $`𝐚`$, $`𝐛`$ and hidden state $`\lambda `$. Similarly, let $`p_2(B|𝐚,𝐛,\lambda )`$ be the probability of obtaining outcome $`B`$ at $`S_2`$ for given $`𝐚`$, $`𝐛`$, $`\lambda `$.

A hidden variables model satisfies the condition of parameter independence (“locality”, in the terminology of Jarrett ) if and only if $`p_1`$ is independent of $`𝐛`$ and $`p_2`$ is independent of $`𝐚`$, so that

$`p_1(A|𝐚,𝐛,\lambda )`$ $`=p_1(A|𝐚,\lambda )`$ (1)
and
$`p(B|𝐚,𝐛,\lambda )`$ $`=p(B|𝐛,\lambda )`$ (2)
for suitable $`p_1(A|𝐚,\lambda )`$ and $`p_2(B|𝐛,\lambda )`$.

Version 2 of the Hess-Philipp model  violates parameter independence. This clearly implies that the model is non-local because it means that if, *per impossibile*, one could prepare an ensemble of pairs all in the same hidden state $`\lambda `$, then experimenters at $`S_1`$ and $`S_2`$ could use these pairs to communicate superluminally.

This point has already been made by Myrvold . However, Myrvold’s paper is largely concerned with a simplified version of the Hess-Philipp model. The extension to the full model is only made briefly, at the end of his paper. In particular, Myrvold does not actually calculate the conditional probabilities appearing in Eqs. (1) and (2). This gives Hess and Philipp  some scope to challenge his conclusion.

In the following we do calculate the conditional probabilities, and explicitly show that they do not satisfy Eqs. (1) and (2)—thereby providing a clear-cut mathematical demonstration that Myrvold is correct. Hess and Philipp are only entitled to disagree with this statement if they can identify a flaw in our *calculations*, and if they can present some alternative, fully explicit calculations leading to expressions which *do* satisfy Eqs. (1) and (2).

In the following we begin by showing that version 1 of the Hess-Philipp model violates signal locality. We then go on to show that, for what is essentially the same reason, version 2 violates the condition of parameter independence.

## 2. Version 1 of the Model Violates Signal Locality

Hess and Philipp construct their model in two stages. The first stage is the theorem proved on pp. 14231–2 of ref. . In Eqs. (22), (23) and (27) they define functions $`A_𝐚(u)`$, $`B_𝐛(v)`$ and a density $`\rho _{\mathrm{𝐚𝐛}}(u,v)`$ with the property

$$A_𝐚(u)B_𝐛(v)\rho _{\mathrm{𝐚𝐛}}(u,v)𝑑u𝑑v=𝐚𝐛$$
(3)
(Eq. (18) in the statement of their theorem). The reader may confirm that one also has111 The reader who does wish to check this statement should bear in mind that there is a missing summation sign in Eq. (26) of ref. , as noted in Hess and Philipp

$`{\displaystyle A_𝐚(u)\rho _{\mathrm{𝐚𝐛}}(u,v)𝑑u𝑑v}`$ $`=𝐚|𝐛|+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{3}{}}}{\displaystyle \underset{r=1}{\overset{n/2}{}}}N_{2r}(|a_k|)\psi _{2r}(|b_k|)`$ (4)
$`{\displaystyle B_𝐛(u)\rho _{\mathrm{𝐚𝐛}}(u,v)𝑑u𝑑v}`$ $`=|𝐚|𝐛`$ (5)
Here $`|𝐚|`$, $`|𝐛|`$ are the vectors with components $`|a_1|,|a_2|,|a_3|`$ and $`|b_1|,|b_2|,|b_3|`$ respectively. $`N_{2r}`$, $`\psi _{2r}`$ and the even integer $`n`$ are the quantities defined in the lemma on p. 14231 of ref. .

The second stage in the construction of the Hess-Philipp model is the complicated combinatoric argument leading to Eq. (35) of ref. . It should be noted that there are two problems with this part of the construction, one major and one minor.

The major problem is the point made by Gill *et al* , that the index $`m`$ featuring in the construction apparently introduces an element of non-locality into the model. In our view Hess and Philipp  do not satisfactorily answer this objection. However, we will not insist on the point here because it is tangential to our argument.

The minor problem is a technical point, concerning the details of the combinatorics. We discuss it in the appendix.

The full probability space for the Hess-Philipp model consists of all pairs $`(\lambda ,\omega )`$, where $`\lambda `$ describes the state of the source particles and $`\omega `$ refers to the detectors (second paragraph on p. 14230 of ref. ). Let $`\nu _\omega `$ be the probability measure on the space of $`\omega `$, and let $`\nu _\lambda `$ be the probability measure on the space of $`\lambda `$. For any random variable $`X`$, let $`E_\lambda (X)=X𝑑\nu _\omega `$ be the conditional expectation value obtained by integrating out $`\omega `$ for a fixed value of $`\lambda `$.

The combinatorics leading to Eq. (35) of ref.  are such that

$$E_\lambda (A_𝐚B_𝐛)=A_𝐚(u)B_𝐛(v)\rho _{\mathrm{𝐚𝐛}}(u,v)𝑑u𝑑v=𝐚𝐛$$
(6)
holds trivially. It is also easily seen that

$$E_\lambda (B_𝐛)=B_𝐛(u)\rho _{\mathrm{𝐚𝐛}}(u,v)𝑑u𝑑v=|𝐚|𝐛$$
(7)
However, the evaluation of $`E_\lambda (A_𝐚)`$ is slightly less straightforward. It is shown in the appendix that

$$E_\lambda (A_𝐚)=𝐚|𝐛|+\frac{1}{2}\left(1|𝐚||𝐛|\right)+\frac{\theta }{16n^2}$$
(8)
with $`0\theta 1`$.

Now let $`E(X)=E_\lambda (X)𝑑\nu _\lambda `$ be the unconditioned expectation value of $`X`$. A striking feature222 It is a curious feature because it means that all the work is being done by the detectors. The source particles might as well not be there.  of Eqs. (68) is that the right-hand sides are independent of $`\lambda `$. Consequently, the unconditioned expectation values are the same as the conditional ones:

$`E(A_𝐚B_𝐛)`$ $`=𝐚𝐛`$ (9)
$`E(B_𝐛)`$ $`=|𝐚|𝐛`$ (10)
$`E(A_𝐚)`$ $`=𝐚|𝐛|+{\displaystyle \frac{1}{2}}\left(1|𝐚||𝐛|\right)+{\displaystyle \frac{\theta }{16n^2}}`$ (11)
whatever the probability measure $`\nu _\lambda `$.

At this stage we notice that the model disagrees with quantum mechanics. For the singlet state quantum mechanics predicts $`E(A_𝐚)=E(B_𝐛)=0`$ for all $`𝐚,𝐛`$. It can be seen that Eqs. (10) and (11) disagree with this prediction.

Not only does the model disagree with quantum mechanics. It also violates *signal* locality333 We are indebted to M. Żukowski for this observation. . To see this, let $`p_2(B=\pm 1|𝐚,𝐛)`$ be the probability of obtaining the measurement outcome $`B=\pm 1`$ at station $`S_2`$, for detector settings $`𝐚,𝐛`$. Then

$$p_2(B=\pm 1|𝐚,𝐛)=\frac{1}{2}\left(1\pm E(B_𝐛)\right)=\frac{1}{2}(1|𝐚|𝐛)$$
(12)
Suppose that Alice at station $`S_1`$ and Bob at station $`S_2`$ have previously agreed that Bob will always measure in the direction $`(1,0,0)`$, while Alice will measure in one of the two alternative directions $`(1,0,0)`$ or $`(0,1,0)`$. Then, if Alice measures in the direction $`(1,0,0)`$ there will be probability 1 of Bob obtaining the result $`+1`$, while if Alice measures in the direction $`(0,1,0)`$ there will only be probability $`1/2`$ of Bob obtaining the result $`+1`$. In this way they can use the arrangement to communicate superluminally, with a probability of error that, with sufficient redundancy, can be made arbitrarily small.

## 3. Version 2 of the Model Violates Parameter Independence

We now turn to version 2 of the Hess-Philipp model, described in the e-print . Unlike version 1, version 2 of the model does reproduce the empirical predictions of quantum mechanics (for the singlet state), and it does not violate signal locality. Nevertheless, it is still non-local. Furthermore, it is non-local for a reason that is closely related to the reason that version 1 violates signal locality—as we now show.

The crucial detail, which is omitted from version 1 of the model, is described at the end of ref. , in Section 5.3. Let $`A_𝐚(\lambda ,\omega )`$, $`B_𝐛(\lambda ,\omega )`$ be the functions describing the measurement outcomes in version 1 of the model. Version 2 is obtained by making the replacements

$`A_𝐚(\lambda ,\omega )`$ $`r(\lambda )A_𝐚(\lambda ,\omega )`$ (13)
$`B_𝐛(\lambda ,\omega )`$ $`r(\lambda )B_𝐛(\lambda ,\omega )`$ (14)
where $`r(\lambda )`$ is a function taking the values $`\pm 1`$, and having the property $`r(\lambda )𝑑\nu _\lambda =0`$ ($`\nu _\lambda `$ being the probability measure on the space of the source variables $`\lambda `$, as before)444 Hess and Philipp stipulate that $`r`$ depends on any “parameter specific to the source (e.g. $`\lambda _1`$, $`\lambda _2`$ or time $`t`$)”. By “time $`t`$” they presumably mean the time at which the particles are emitted by the source. In the definition of parameter independence, as given in Eqs. (1), (2) above, $`\lambda `$ denotes a *complete* specification of the state of the source. On this defintion $`\lambda `$ must be taken to include a specification of the time of emission. .

With these replacements the model does reproduce the correct quantum mechanical predictions for the unconditioned expectation values $`E(A_𝐚B_𝐛)`$, $`E(A_𝐚)`$ , $`E(B_𝐛)`$ (in the singlet state). Consequently, it does not violate signal locality.

Suppose, however, that one considers the *conditional* expectation values $`E_\lambda (A_𝐚)`$, $`E_\lambda (B_𝐛)`$. Then one has, in place of Eqs. (7), (8) above,

$`E_\lambda (A_𝐚)`$ $`=r(\lambda )\left(𝐚|𝐛|+{\displaystyle \frac{1}{2}}\left(1|𝐚||𝐛|\right)+{\displaystyle \frac{\theta }{16n^2}}\right)`$ (15)
$`E_\lambda (B_𝐛)`$ $`=r(\lambda )|𝐚|𝐛`$ (16)
Let $`p_1`$, $`p_2`$ be the probabilities appearing on the left-hand sides of Eqs. (1), (2) above. Then

$`p_1(A=\pm 1|𝐚,𝐛,\lambda )`$ $`={\displaystyle \frac{1}{2}}\left(1\pm E_\lambda (A_𝐚)\right)`$ (17)
$`p_2(B=\pm 1|𝐚,𝐛,\lambda )`$ $`={\displaystyle \frac{1}{2}}\left(1\pm E_\lambda (B_𝐛)\right)`$ (18)
implying

$`p_1(A=\pm 1|𝐚,𝐛,\lambda )`$ $`={\displaystyle \frac{1}{2}}\left(1\pm r(\lambda )\left(𝐚|𝐛|+{\displaystyle \frac{1}{2}}\left(1|𝐚||𝐛|\right)+{\displaystyle \frac{\theta }{16n^2}}\right)\right)`$ (19)
$`p_2(B=\pm 1|𝐚,𝐛,\lambda )`$ $`={\displaystyle \frac{1}{2}}\left(1r(\lambda )|𝐚|𝐛\right)`$ (20)
These expressions clearly fail to satisfy the condition of parameter independence, stated in Eqs. (1), (2) above. It follows that version 2 of the model is non-local.

The violation of parameter independence in version 2 of the model is closely related to the violation of signal locality in version 1. This can be seen by comparing the expression for $`p_2(B=\pm 1|𝐚,𝐛,\lambda )`$ in version 2 (see Eq. (20) above) with the expression for $`p_2(B=\pm 1|𝐚,𝐛)`$ in version 1 (see Eq. (12) above).

A model which violates parameter independence is one which would violate signal locality in an imaginary world, where it was possible to obtain *complete* information regarding the state of the source, including the values of all the quantities which are in fact “hidden”. For version 2 of the Hess-Philipp model complete information regarding of the state of the source would have to include the value of the function $`r(\lambda )`$.

Suppose that we are in such an imaginary world. As in the discussion following Eq. (12) above, consider two experimenters, Alice at $`S_1`$ and Bob at $`S_2`$. Suppose that there is also a third experimenter Xenophon located in the intersection of Alice and Bob’s backward light cones. Xenophon prepares a succession of pairs in the singlet state. For each pair he determines the value of $`r(\lambda )`$ and then sends the particles to Alice and Bob. He also sends the value of $`r(\lambda )`$ to Bob (by telephone, say). Bob always measures in the direction $`(1,0,0)`$. If Alice measures in the direction $`(1,0,0)`$ then Bob obtains the value $`r(\lambda )`$ with probability 1. If, on the other hand, Alice measures in the direction $`(0,1,0)`$ there is only probability $`1/2`$ of Bob obtaining the value $`r(\lambda )`$. The value $`r(\lambda )`$ is known to Bob. It is therefore possible for Alice to send superluminal signals to Bob, just as in the case discussed in the passage following Eq. (12) above.

Finally, let us note that it does not help if one modifies the model again, by making $`r`$ a function of the state of the detectors instead of $`\lambda `$. It is easily verified that the model, thus modified, would still violate Shimony’s  condition of outcome dependence (“completeness” in the terminology of Jarrett )—meaning that the model would still be non-local. Moreover, the model would still allow Alice and Bob to communicate superluminally in an imaginary world were it was possible to obtain complete information regarding the state of any system (although the reading of Bob’s detector is given by $`rB`$, in such a world Bob could still find out the value of $`B`$ by inspecting the hidden state of his detector). This means that, on the level of the hidden variables, there must be superluminal exchanges of information between the two detectors.

## 4. Conclusion

We conclude that the expressions for the probabilities $`p_1`$ and $`p_2`$ calculated in this paper have the inescapable implication that the Hess-Philipp model is non-local. Hess and Philipp are only entitled to challenge this conclusion if they can find an error in our calculations, and if they can provide alternative calculations leading to expressions for $`p_1`$ and $`p_2`$ which do not have the implication of non-locality.

## Acknowledgment

We are grateful to M. Żukowski for some very useful discussions, and for his encouragement.

## Appendix A Remarks on the Combinatoric Part of the Construction

In this appendix we first discuss a minor technical difficulty with the construction leading to Eq. (35) of reference . We then go on to derive the expression for $`E_\lambda (A_𝐚)`$ given in Eq. (8) above.

Hess and Philipp choose in turn each of the $`(n+1)^2`$ squares $`Q_{jk}`$. For each choice of $`Q_{jk}`$ they define $`L`$ measures, where $`L`$ is the binomial coefficient $`{}_{}{}^{9n^2}C_{3n}^{}`$. Each of these $`L`$ measures assigns the weights $`|a_1b_1|`$, $`|a_2b_2|`$, $`|a_3b_3|`$ to $`3`$ of the unit squares $`Q_{jk}`$, and the weights $`(1/2)N_1(|a_1|)\psi _1(|b_1|),\mathrm{},(1/2)N_n(|a_3|)\psi _n(|b_3|)`$ to $`3n`$ of the $`9n^2`$ unit squares not contained in the vertical and horizontal strips defined by $`Q_{jk}`$. It assigns every other unit square weight $`0`$. The assignment is such that each unit square $`Q_{jk}`$ is assigned the weight $`|a_tb_t|`$ by $`L/9`$ measures, and each unit square not in the vertical and horizontal strips defined by $`Q_{jk}`$ is assigned the weight $`(1/2)N_i(|a_t|)\psi _i(|b_t|)`$ by $`L/(9n^2)`$ measures ($`t=1,2,3`$ and $`i=1,\mathrm{},n`$).

The problem with this construction is that it tacitly assumes that the binomial coefficient $`{}_{}{}^{9n^2}C_{3n}^{}`$ is divisible by $`9n^2`$. However, this is typically not the case (the only even integers $`100`$ for which it is true are $`n=10,40,44`$ and $`84`$). Moreover, it is not clear to us that the other requirements can be satisfied even when $`n`$ does satisfy this condition. The difficulty is, however, easily resolved if, instead of taking $`L=^{9n^2}C_{3n}`$, we take it to be the permutation $`{}_{}{}^{9n^2}P_{3n}^{}=(9n^2)!/(9n^23n)!`$.

Now let us turn to the derivation of Eq. (8). Let $`I_{jk}`$ be the set of indices $`m`$ such that $`\mu _m`$ is one of the $`L`$ measures associated with $`Q_{jk}`$. It is easily seen that

$$\underset{mI_{jk}}{}_{Q_{jk}}A_m(u)𝑑\mu _m=L𝐚|𝐛|$$
(21)
Let $`S_{jk}`$ be the subset of $`[3,3n)^2`$ which is obtained by deleting the horizontal and vertical strips defined by $`Q_{jk}`$. For each $`mI_{jk}`$, $`_UA_m(u)𝑑u𝑑v=1`$ for half the unit squares $`US_{jk}`$, and it $`=0`$ for the other half. Consequently

$$\underset{mI_{jk}}{}_{S_{jk}}A_m(u)𝑑\mu _m=\frac{L}{4}\underset{t=3}{\overset{3}{}}\underset{i=1}{\overset{n}{}}N(|a_t|)\psi (|b_t|)=\frac{L}{2}\left(1|𝐚||𝐛|+\frac{\theta }{8n^2}\right)$$
(22)
with $`0\theta 1`$ (where we have used Eq. (21) of ref. ). Hence

$$\underset{mI_{jk}}{}_{[3,3n)^2}A_m(u)𝑑\mu _m=L\left(𝐚|𝐛|+\frac{1}{2}(1|𝐚||𝐛|)+\frac{\theta }{16n^2}\right)$$
(23)
The right-hand side of this equation is independent of $`j,k`$. Consequently, the effect of summing over all $`j`$, $`k`$ is simply to multiply the expression by $`(n+1)^2`$. After dividing by $`N=(n+1)^2L`$ this gives

$$E_\lambda (A_𝐚)=\frac{1}{N}\underset{m=1}{\overset{N}{}}_{[3,3n)^2}A_m(u)𝑑\mu _m=𝐚|𝐛|+\frac{1}{2}(1|𝐚||𝐛|)+\frac{\theta }{16n^2}$$
(24)


text:
# Invariants Associated to Orthogonal ϵ-constants

## 1 Introduction and Background

Brauer group invariants associated to motives have been long studied by mathematicians. In the case of orthogonal motives of even weight, these invariants were first examined by Frohlich in \[F\]. Deligne used these Brauer group invariants and their relationship with certain $`ϵ`$-constants in order to give a proof of the Frohlich-Queyrut Theorem in \[D2\] by showing that certain global orthogonal root numbers were one by interpreting the associated local orthogonal root numbers as Stiefel Whitney classes and then using the local root numbers to define an element of order two in the Brauer group of $``$. Similar results and methods have been done by Saito (in \[tS2\], for example) and others. In this paper, we will define a Brauer group invariant associated to certain motives which are symplectic and have weight one.

In order to construct the relevant motives we first define $`𝒳`$ to be an arithmetic surface of dimension $`2`$ which is flat, regular, and projective over $``$. Throughout this paper we will assume that $`f:𝒳Spec()`$ is the structure morphism. Define $`G`$ to be a finite group which acts tamely on $`𝒳`$. In other words, for each closed point $`xX`$, the order of the inertia group of $`x`$ is relatively prime to the residue characteristic of $`x`$. Let $`𝒴`$ be the quotient scheme $`𝒳/G`$, which we will assume is regular, and that for all finite places $`v`$ the fiber $`𝒴_v=(𝒳_v)/G=𝒴_{}(/p(v))`$ has normal crossings and smooth irreducible components with multiplicities relatively prime to the residue characteristic of $`v`$. Finally, let $`V`$ be a representation of $`G`$ over $`\overline{}`$

We wish to define a class $`\mu (𝒳,G,V)`$ in the global Brauer group $`H^2(,/2)`$ which will depend only on $`𝒳,G,`$ and $`V`$. In particular, $`\mu (𝒳,G,V)`$ will be the element whose local invariant at the place $`v`$ is given by the sign of $`ϵ_v(D^{},V)`$, where $`D^{}`$ is an appropriately chosen horizontal canonical divisor on $`𝒴`$. In other words, $`O_𝒴(D^{}+𝒴_T)`$ is isomorphic to $`\omega _{𝒴/}(𝒴_S^{red})`$ where $`𝒴_T`$ and $`𝒴_S`$ are unions of vertical fibers of $`𝒴`$.

In addition to being a canonical divisor, we must impose an additional condition on our choice of $`D^{}`$ in order for $`\mu (𝒳,G,V)`$ to be well-defined. To do this we recall that Proposition $`3.10`$ of \[G\] showed that if $`F^{}`$ and $`G^{}`$ are components of $`𝒴_v^{red}`$ and $`D^{}`$ is a canonical divisor in the above sense, then there is a canonical isomorphism between $`𝒪_F^{}(D^{}F^{})`$ and $`\omega _F^{}(F^{}G^{})`$. This isomorphism maps the global section $`1\mathrm{\Gamma }(𝒪_F^{}(D^{}F^{}))`$ to an element $`\gamma \mathrm{\Gamma }(\omega _F^{}(F^{}G^{}))`$ such that $`\gamma `$ has a simple pole at all $`xF^{}G^{}`$. Define $`a_x`$ to be the residue of $`\gamma `$ at these points $`x`$. We wish to choose $`D^{}`$ so that all of the residue terms $`a_x`$ are equal to $`1`$ and we will call such a choice of $`D^{}`$ a nice divisor. Nice divisors exist due to a moving lemma which we will make explicit in Section 2. We can now state the main theorem.

###### Theorem 1.1.

In the case where $`V`$ is an orthogonal representation of dimension equal to zero and trivial determinant and where $`D^{}`$ is a nice divisor, the local constant $`ϵ_v(D^{},V)`$ is independent of the choice of $`D^{}`$. In particular, the element $`\mu (𝒳,G,V)`$ in the global Brauer group $`H^2(,/2)`$ whose local invariant at the place $`v`$ is given by the sign of $`ϵ_v(D^{},V)`$ is well-defined.

Section two of this paper will prove theorem 1.1. Section three recalls the definition of the Galois theoretic invariant $`w_2(\pi )`$ defined by Cassou-Nogues, Erez and Taylor in \[CNET\] which comes from a situation similar to the one we are working in and proves a connection between the two invariants. Finally, in section four we prove the moving lemma which we used in section two in order to show that nice divisors always exist.

The author would like to thank Ted Chinburg for his insightful comments and suggestions, in particular on the proof in section four.

## 2 Definition of $`\mu (𝒳,G,V)`$

In \[G\], the author proves the following result about orthogonal $`ϵ`$-constants associated to tame finite group actions on surfaces (We refer the reader to \[D1\], \[CEPT1\], and \[G\] for all relevant definitions.):

###### Theorem 2.1.

If $`V`$ is an orthogonal virtual representation of degree zero and trivial determinant then for all finite places $`v`$ of $``$ we have the following formula:

$`{\displaystyle \frac{ϵ(𝒴_v,V)}{ϵ(D_v^{},V)}}={\displaystyle \underset{i}{}}det(V^{I_i})(\delta _{v,C_i}){\displaystyle \underset{zZ}{}}ϵ_{0,z}(C_{i_1},V^{I_{i_1}})ϵ_{0,z}(C_{i_2},V^{I_{i_2}})ϵ(z,V)`$
where $`Z`$ is the set of crossing points of the components of the fiber $`𝒴_v`$ and $`D^{}`$ is chosen to be a canonical horizontal divisor on $`𝒴`$.

By canonical we mean that $`D^{}`$ is a horizontal divisor on $`𝒴`$ such that $`D^{}+𝒴_T=K_𝒴+𝒴_S^{red}`$, where $`𝒴_S^{red}`$ is the sum of the reductions of the fibers of $`𝒴`$ at the places in the set $`S`$ of bad primes, $`𝒴_T`$ is the sum of the (necessarily reduced) fibers of $`𝒴`$ over the places in $`T`$. We also wish to choose $`D^{}`$ so that it intersects the non-smooth fibers $`𝒴_v`$ of $`𝒴`$ transversally at smooth points on the reduction of $`𝒴_v`$.

Given any such horizontal divisor $`D^{}`$ it is possible to find a nice divisor which is close to it due to the following moving lemma, which is proven in section 4. In particular, this shows that nice divisors $`D^{}`$ always exist and therefore that our class $`\mu (𝒳,G,V)`$ will be well-defined.

###### Lemma 2.2.

There exists a meromorphic function $`h`$ on $`𝒴_v^{red}`$ such that the divisor of $`h`$ intersects the special fibers $`𝒴_v^{red}`$ transversally at smooth points away from $`D_v^{}`$ and such that $`h`$ takes on prescribed values at the singular points of $`𝒴_v^{red}`$. In particular, given a horizontal divisor $`D^{}`$ as in the previous section, the divisor $`D^{}+div(h)`$ will have residue maps equal to one at the crossing points of components of $`𝒴_v^{red}`$.

Proof of theorem 1.1: For any fixed place $`v`$ of $``$ and any component $`C_i`$ of $`𝒴_v^{red}`$, it follows from the definition that $`\delta _{v,C_i}`$ is the class which corresponds to the element $`\delta =(0)(a_x)(_{xC_iZ})(_{xC_iZ}K^{}/U_x^1)`$. This term is independent of the choice of nice divisor and it follows that the right hand side of the equation in Theorem 2.1 is as well. Furthermore, it is clear that $`ϵ(𝒴_v,V)`$ is independent of our choice of $`D^{}`$. Thus, it follows from the theorem that $`ϵ_v(D^{},V)`$ is independent of the choice of $`D^{}`$. Next, we note that $`ϵ_{v,0}(D^{},V)`$ must also be independent of our choice of $`D^{}`$ (by Lemma 3.3 of \[G\], for example). Therefore it must be the case that $`ϵ_v(D^{},V)=ϵ_{v,0}(D^{},V)ϵ(D_v^{},V)`$ is independent of the choice of $`D^{}`$. The product of all of the $`ϵ_v(D^{},V)`$ is equal to $`ϵ(D^{},V)`$, which must be equal to one from the theorem of Fröhlich and Queyrut \[FQ\]. This tells us that we can define an element $`\mu (𝒳,G,V)`$ in $`H^2(,/2)`$ by setting the local component at the prime $`v`$ to be equal to the sign of $`ϵ_v(D^{},V)`$. $`\mathrm{}`$

## 3 The connection to $`w_2(\pi )`$

Let $`\pi :𝒳𝒴`$ be a tamely ramified cover of degree $`n`$, where $`𝒳`$ and $`𝒴`$ are regular schemes and $`𝒴`$ is connected. Furthermore, we must make the technical assumption that the ramification indices are all odd. Cassou-Nogues, Erez, and Taylor use Grothendeick’s equivariant cohomology theory to define an invariant $`w_2(𝒳/𝒴)=w_2(\pi )H^2(𝒴_{et},/2)`$ associated to this situation. Their definition generalizes to define classes $`w_i(\pi )`$ which lie in $`H^i(𝒴_{et},/2)`$ for all positive integers $`i`$, but in this thesis we will only be interested in $`w_2`$. These terms are generalized Stiefel-Whitney classes, and are obtained by pulling back the universal Hasse-Witt classes defined by Jardine using classifying maps related to a quadratic form $`E`$. The precise definition of $`E`$ uses the existence of a locally free sheaf $`𝒟_{𝒳/𝒴}^{1/2}`$ whose square is the inverse different of the covering $`𝒳/𝒴`$. In this section, we will consider the relationship between the class $`\mu (𝒳,G,V)`$ lying in $`H^2(,/2)`$ which we defined in theorem 1.1 and $`w_2(\pi )`$.

In \[CNET\] the following equality in $`H^2(𝒴_{et},/2)`$, which is an analog of a theorem of Serre, is proved:

$$w_2(\pi _{}(𝒟_{𝒳/𝒴}^{1/2},Tr_{𝒳/𝒴}))=w_2(\pi )+(2)(d_{𝒳/𝒴})+\rho (𝒳/𝒴)$$

where $`\rho (𝒳/𝒴)`$ is defined by the ramification of $`𝒳/𝒴`$, $`d_{𝒳/𝒴}`$ is the function field discriminant, and the left hand side of the equation is the Hasse-Witt invariant associated to the square-root of the inverse different bundle. Note that if we look at the one-dimensional version of this formula the middle term on the right hand side becomes trivial. Therefore, in the case of étale covers of curves the formula reduces to

$$w_2(\pi )=w_2(\pi _{}(𝒟_{𝒳/𝒴}^{1/2},Tr_{𝒳/𝒴}))=w_2(E)$$

where $`w_2(E)`$ is the second Hasse-Witt invariant associated to the square root of the inverse different, as described in detail in \[CNET\].

Let $`D^{}`$ be a choice of a canonical divisor on $`𝒴`$ in the sense of the previous sections, and let $`i:D^{}𝒴`$ be the natural inclusion. An étale covering of $`𝒴`$ naturally restricts to give an étale covering of $`D^{}`$. We now have the following natural maps

$$i^{}:H^2(𝒴_{et},/2)H^2(D_{et}^{},/2)$$

$$res:H^2(D_{et}^{},/2)H_{et}^2((D^{}),/2)=H_{gal}^2(\overline{}/(D^{}),/2)$$

$$cor:H_{gal}^2(\overline{}/(D^{}),/2)H^2(,/2)$$

where the latter two maps are restriction and corestriction in the sense of Serre (for details see Chapter VII of \[Se\]). Composing these maps gives a natural map

$$H^2(𝒴_{et},/2)H^2(,/2)$$

We denote the image of the class $`w_2(\pi )H^2(𝒴_{et},/2)`$ under this map by $`\stackrel{~}{w}_2(\pi )H^2(,/2)`$. At first glance it appears as though this element may depend on our choice of canonical divisor $`D^{}`$. However, we want to show that it does not depend on this choice and furthermore that the element $`\stackrel{~}{w}_2(\pi )`$ is connected in a natural way to the element $`\mu (𝒳,G,V)`$. Recall that $`\mu (𝒳,G,V)`$ is defined by letting the local invariant at the place $`v`$ be given by the sign of $`ϵ(D_v^{},V)`$. This also appears to depend on the choice of $`D^{}`$ but turns out to be independent of choice.

The first natural observation is that the class $`\mu (𝒳,G,V)`$ depends on the choice of a representation $`V`$ of $`G`$ while $`\stackrel{~}{w}_2(\pi )`$ does not. The natural representation to consider is the regular representation of the group $`G`$, which we denote by $`R`$. In particular, the nicest possible theorem comparing the invariants would say that $`\mu (𝒳,G,R)=\stackrel{~}{w}_2(\pi )`$. However,we have only shown that $`\mu (𝒳,G,V)`$ is a well-defined class in the case where $`V`$ is of dimension zero and of trivial determinant, neither of which holds for $`R`$. So instead of setting $`V=R`$, we consider the representation $`V=Rdet(R)T^{n1}`$, where $`det(R)`$, the determinant of the regular representation, is a character whose order is either one or two, $`T`$ is the trivial representation and $`n`$ is the degree of the cover $`𝒳/𝒴`$. This choice of $`V`$ is an orthogonal representation, and it has trivial determinant and dimension $`0`$. We can now prove the following theorem:

###### Theorem 3.1.

Assume that we are in the above situation, and in particular that $`V=Rdet(R)T^{n1}`$. Let $`𝒴_1/𝒴`$ be either the trivial cover or the subcover of $`𝒳/𝒴`$ of degree 2, depending on whether $`det(R)`$ is of order $`1`$ or $`2`$ respectively. Then as classes in $`H^2(,/2)`$, we have the equality

$$\mu (𝒳,G,V)=\stackrel{~}{w}_2(𝒳/𝒴)\stackrel{~}{w}_2(𝒴_1/𝒴)(n1)\stackrel{~}{w}_2(𝒴/𝒴)$$

The proof of this theorem relies on the interpretation of each side of the equation as a Stiefel-Whitney class. In particular, Cassou-Nogues, Erez, and Taylor show that the element $`w_2(\pi )`$ is the Hasse-Witt invariant associated to the full covering of surfaces. Thus, when we restrict the class to the one-dimensional divisor $`D^{}`$ we see that the element $`i^{}(w_2(\pi ))H^2(D_{et}^{},/2)`$ is equal to the Stiefel-Whitney class associated to the form $`E^{}=(𝒟_{D/D^{}}^{1/2},Tr_{D/D^{}})`$ on the canonical divisor $`D^{}`$ of $`𝒴`$. This follows as a generalization to étale cohomology of results of Fröhlich in \[F\], which allow us to associate the class $`i^{}(w_2(\pi ))`$ to $`G`$-extensions of the ring of integers of the residue field of the generic point of $`D^{}`$.

Next we make use of the results of Deligne which allow us to interpret local Stiefel-Whitney classes in terms of local root numbers. In particular, the following lemma is shown in \[D2\]:

###### Lemma 3.2.

Let $`d=1`$, so that the fibers $`𝒳_v`$ and $`𝒴_v`$ are all one dimensional schemes. Furthermore, let $`V`$ be an orthogonal virtual representation of dimension zero and trivial determinant. Under these hypotheses, the local root number $`W(V_v)=sign(ϵ_v(𝒴,V))`$ is equal to $`exp(2\pi icl(sw_v))`$, where $`sw_v`$ is the local Stiefel-Whitney class, and $`cl(sw_v)\{0,1/2\}/`$.

In other words, in characteristic not equal to two, the sign of the $`ϵ`$-constants $`ϵ_v(D^{},V)`$ of the representation on the one dimensional horizontal divisor $`D^{}`$ are determined by whether or not the classes $`w_2(\pi )`$ are trivial in the Brauer group, and $`ϵ_v(D^{},V)`$ is automatically positive when $`v=2`$. It turns out that these are exactly the terms that are coming up in the computation of the class of Cassou-Nogues, Erez, and Taylor.

In particular, $`ϵ_v(D^{},V)=ϵ_v(D^{},R)ϵ_v(D^{},det(R))`$ is the same as the local Hasse-Witt invariants. However, we are working with étale covers of curves and so from the results of \[CNET\] discussed above, these Hasse-Witt invariants are simply the images of the appropriate classes $`w_2(\pi )`$. This proves Theorem 3.1. $`\mathrm{}`$

## 4 Proof of Lemma 2.2

The proof of Lemma 2.2 involves a generalized version of Bertini’s Theorem. For now, let us assume that $`X`$ is a smooth curve defined over an infinite field $`k`$ and let us choose a finite set of points $`p_1,\mathrm{},p_m`$ on $`X`$. We define the divisor $`p=_ip_i`$. Furthermore, let us choose constants $`c_i`$ which lie in the residue field $`k(p_i)`$ of the points $`p_i`$. Finally, let us choose $`\mathrm{\Lambda }`$ to be an effective very ample divisor on $`X`$ of large degree which is supported off of $`p`$. We look at the group of global sections $`H^0(X,𝒪_X(\mathrm{\Lambda }))`$ and let $`f_0,\mathrm{},f_t`$ be a basis of this group. This basis gives us a projective embedding from $`X`$ into $`_k^t`$ whose projective coordinates we will write as $`x_0,\mathrm{},x_t`$.

We wish to prove that there exist linear forms $`l_0`$ and $`l_1`$ in the variables $`x_i`$ such that the following properties hold:

1. If $`H_i`$ is the hyperplane defined by $`l_i=0`$ in $`_k^t`$ then $`H_iX`$ is a finite set of closed points which is regular and disjoint from $`\{p_1,\mathrm{},p_m\}`$. Furthermore, we wish to choose the $`l_i`$ so that $`H_1H_2X`$ is empty.
2. It follows from (1) that the function $`l_1/l_0|_X`$ is in $`𝒪_{X,p_i}`$ for each $`i`$. However, we further wish to specify that the image of $`l_1/l_0`$ in $`k(p_i)`$ is the given constant $`c_i`$.

The classical version of Bertini’s Theorem (Theorem $`II\mathrm{.8.18}`$ of \[H\]) tells us that there exist linear forms $`l_0`$ so that $`H_0`$ satisfies condition (1). We now fix one choice of such an $`l_0`$, and we will attempt to construct an $`l_1`$ so that the pair satisfies properties (1) and (2). We begin by looking at the set $`V`$ consisting of all linear forms such that $`\{l_0,l_1\}`$ satisfy condition (2). In other words,

$$V=\{l=a_0x_0+\mathrm{}+a_tx_tj,\frac{l}{l_0}|_X(p_j)=c_jk(p_j)\}$$

This $`V`$ will be an affine space over $`k`$. Furthermore, because we chose the divisor $`\mathrm{\Lambda }`$ to have high degree it follows from a Riemann-Roch argument that $`V`$ is of codimension $`m`$ inside of $`H^0(X,𝒪_X(\mathrm{\Lambda }))`$.

For each point $`xX`$, we now define a set $`V_xV`$ which consists of all linear forms $`lV`$ so that the hyperplane defined by $`l=0`$ has contact order $`>1`$ at $`x`$. In other words, $`V_x`$ will consist of those linear forms who do not intersect $`X`$ nicely at the point $`x`$. We can again use the Riemann-Roch theorem to show that for almost all choices of $`x`$, we get that the dimension of $`V_x`$ is equal to $`dim`$ $`V2`$.

Let $`U=X\{p_1,\mathrm{},p_m\}`$ so that $`U`$ is an affine curve, and define $`TU\times V`$ to be the set of all pairs $`(x,l)`$ such that $`xU`$ and $`lV_x`$. We have seen that the projection map $`\pi :TU`$ is surjective and for almost all $`xU`$ (in particular for those points such that $`k(x)=k`$), we see that the fiber $`\pi ^1(x)`$ is an affine space whose dimension is equal to $`dim`$ $`V2`$. In particular, this shows that $`T`$ is irreducible and that the dimension of $`T`$ is equal to $`dim`$ $`V1`$. But this shows us that the natural projection map $`\gamma :TV`$ must not be surjective.

In particular, we can choose some element $`l_1V`$ which is not in the image of $`\gamma `$. In particular, the hyperplane $`H`$ defined by $`\{l_1=0\}`$ is such that $`HU`$ is regular and, since $`l_1V`$, we know that $`l_1/l_0(p_i)=c_ik(p_i)`$, and thus that $`l_1`$ and $`l_o`$ satisfy conditions (1) and (2) above.

So far, we have only made the argument for the case where $`X`$ is a smooth curve. However, as long as $`X`$ is a reduced curve with smooth irreducible components which have normal crossings, then the same argument will hold as long as we include these crossing points in the set of $`\{p_i\}`$. Instead of using the normal Riemann-Roch theorem we will now use the version for singular curves described on p.298 of \[H\].

In order to prove lemma 2.2 we apply this generalized version of Bertini’s theorem to $`X=𝒴_v^{red}`$. Specifically, choose the set of points $`\{p_1,\mathrm{},p_m\}`$ to include the crossing points of components of $`𝒴_v^{red}`$ as well as the points in $`D^{}𝒴_v^{red}`$. The above argument then allows us to find a meromorphic function $`h`$ where we can specify the values of the function $`h=l_1/l_0`$ at the crossing points of $`𝒴_v^{red}`$ so that the residues that come up when we consider $`D^{\prime \prime }=D^{}+div(h)`$ are all equal to one and $`D^{\prime \prime }`$ intersects $`𝒴_v^{red}`$ in the way we want. $`\mathrm{}`$