text:
# 1 Introduction

## 1 Introduction

One of the most intriguing objects of combinatorial geometry is the *$`n^{\text{th}}`$ Birkhoff polytope*

$$_n=\{\left(\begin{array}{ccc}x_{11}& \mathrm{}& x_{1n}\\ \mathrm{}& & \mathrm{}\\ x_{n1}& \mathrm{}& x_{nn}\end{array}\right)^{n^2}:x_{jk}0,\begin{array}{c}_jx_{jk}=1\text{ for all }1kn\hfill \\ _kx_{jk}=1\text{ for all }1jn\hfill \end{array}\},$$

often described as the set of all $`n\times n`$ *doubly stochastic matrices*. $`_n`$ is a convex polytope with integer vertices. It possesses fascinating combinatorial properties  and relates to many mathematical areas . A long-standing open problem is the determination of the relative volume of $`_n`$, which had been known only up to $`n=8`$ . In this paper, we propose a new method of calculating this volume and use it to compute $`\mathrm{vol}_9`$ and $`\mathrm{vol}_{10}`$.

One of the recent attempts to compute $`\mathrm{vol}_n`$ relies on the theory of counting functions for the integer points in polytopes. Ehrhart proved  that for a polytope $`𝒫^d`$ with integral vertices, the number

$$L_𝒫(t):=\mathrm{\#}\left(t𝒫^d\right)$$

is a polynomial in the positive integer variable $`t`$. He showed various other properties for this counting function (in fact, in the more general setting of $`𝒫`$ having rational vertices), of which we mention three here:

* The degree of $`L_𝒫`$ is the dimension of $`𝒫`$.
* The leading term of $`L_𝒫`$ is the relative volume of $`𝒫`$, normalized with respect to the sublattice of $`^d`$ on the affine subspace spanned by $`𝒫`$.
* Since $`L_𝒫`$ is a polynomial, we can evaluate it at nonpositive integers. These evaluations yield

$`L_𝒫(0)=\chi (𝒫),`$ (1)
$`L_𝒫(t)=(1)^{dim𝒫}L_𝒫^{}(t).`$ (2)

Here $`\chi (𝒫)`$ denotes the Euler characteristic, $`𝒫^{}`$ the relative interior of $`𝒫`$. The reciprocity law (2) was in its full generality proved by Macdonald .

The application of this theory to the Birkhoff polytope $`_n`$ incorporates the nice interpretation of the number of integral points in $`t_n`$ as the number of *semi-magic squares*, namely, square matrices whose nonnegative integral coefficients sum up to the same integer $`t`$ along each row and column.

We will denote the Ehrhart polynomial of $`_n`$ by

$$H_n(t):=L__n(t).$$

It is not hard to see that $`dim_n=(n1)^2`$, hence $`H_n`$ is a polynomial in $`t`$ of degree $`(n1)^2`$. The first two of these polynomials are trivial:

$$H_1(t)=1,H_2(t)=t+1,$$

the first nontrivial case was computed by MacMahon  as

$$H_3(t)=3\left(\genfrac{}{}{0pt}{}{t+3}{4}\right)+\left(\genfrac{}{}{0pt}{}{t+2}{2}\right).$$
(3)
The structural properties of $`H_n`$ were first studied in . It is a nice exercise to deduce from (2) that

$$H_n(nt)=(1)^{n1}H_n(t)$$
(4)
and

$$H_n(1)=H_n(2)=\mathrm{}=H_n(n+1)=0.$$

This allows the following strategy of computing $`H_n`$, and therefore, the volume of $`_n`$: compute the first $`\left(\genfrac{}{}{0pt}{}{n1}{2}\right)`$ values of $`H_n`$, use the above symmetry and trivial values of $`H_n`$, and calculate the polynomial $`H_n`$ by interpolation. In fact, as far as we are aware of, the volume of $`_8`$ was computed using essentially this method, combined with some nice computational tricks .

We propose a new, completely different approach of computing $`H_n`$ (and hence $`\mathrm{vol}_n`$). It is based on an analytic method by the first author of computing the Ehrhart polynomial of a polytope . We will introduce the application of this method to the Birkhoff polytope in the following section.

Some recent refreshing approaches—of a more algebraic-geometric/topological flavor—to the problem of computing $`\mathrm{vol}_n`$ and $`H_n`$ can be found in .

## 2 An integral counting integers

A convex polytope $`𝒫^d`$ is an intersection of halfspaces. This allows the compact description

$$𝒫=\{𝐱^d:\mathrm{𝐀𝐱}𝐛\},$$

for some $`(m\times d)`$-matrix $`𝐀`$ and $`m`$-dimensional vector $`𝐛`$. Here the inequality is understood componentwise. In fact, we may convert these inequalities into equalities by introducing ‘slack variables.’ If $`𝒫`$ has rational vertices (those polytopes are called *rational*), we can choose $`𝐀`$ and $`𝐛`$ in such a way that all their entries are integers. In summary, we may assume that a convex rational polytope $`𝒫`$ is given by

$$𝒫=\{𝐱_0^d:\mathrm{𝐀𝐱}=𝐛\},$$
(5)
where $`𝐀M_{m\times d}()`$ and $`𝐛^m`$. (If we are interested in counting the integer points in $`𝒫`$, we may assume that $`𝒫`$ is in the nonnegative orthant, i.e., the points in $`𝒫`$ have nonnegative coordinates, as translation by an integer vector does not change the lattice-point count.) The following straightforward theorem can be found in . We use the standard multivariate notation $`𝐯^𝐰:=v_1^{w_1}\mathrm{}v_n^{w_n}`$.

###### Theorem 1

\[2, Theorem 8\] Suppose the convex rational polytope $`𝒫`$ is given by (5), and denote the columns of $`𝐀`$ by $`𝐜_1,\mathrm{},𝐜_d`$. Then

$$L_𝒫(t)=\frac{1}{(2\pi i)^m}_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_m\right|=ϵ_m}\frac{z_1^{tb_11}\mathrm{}z_m^{tb_m1}}{\left(1𝐳^{𝐜_1}\right)\mathrm{}\left(1𝐳^{𝐜_d}\right)}𝑑𝐳.$$

Here $`0<ϵ_1,\mathrm{},ϵ_m<1`$ are distinct real numbers.

It should be mentioned that $`L_𝒫`$ is in general *not* a polynomial if the vertices of $`𝒫`$ are not integral, but a quasipolynomial, that is, an expression of the form

$$c_d(t)t^d+\mathrm{}+c_1(t)t+c_0(t),$$

where $`c_0,\mathrm{},c_d`$ are periodic functions in $`t`$. (See, for example, \[19, Section 4.4\] for more information about quasipolynomials.) Theorem 1 applies to these slightly more general counting functions; however, in this article we will only deal with polytopes with integer vertices, for which $`L_𝒫`$ is a polynomial.

We also note here that Theorem 1 can be used to quickly compute by hand formulas for certain classes of polytopes (see, for example, ). In this project, we take a slightly different approach and use this theorem to efficiently derive formulas with the help of a computer.

We can view the Birkhoff polytope $`_n`$ as given in the form of (5), where

$$𝐀=\left(\begin{array}{cccccccccccccc}1& \mathrm{}& 1& & & & & & & & & & & \\ & & & 1& \mathrm{}& 1& & & & & & & & \\ & & & & & & \mathrm{}& & & & & & & \\ & & & & & & & 1& \mathrm{}& 1& & & & \\ 1& & & 1& & & & 1& & & & & & \\ & \mathrm{}& & & \mathrm{}& & \mathrm{}& & \mathrm{}& & & & & \\ & & 1& & & 1& & & & 1& & & & \end{array}\right)$$

is a $`(2n\times n^2)`$-matrix and $`𝐛=(1,\mathrm{},1)^{2n}`$. Hence Theorem 1 gives for this special case

$$H_n(t)=\frac{1}{(2\pi i)^{2n}}\mathrm{}\frac{(z_1\mathrm{}z_{2n})^{t1}}{(1z_1z_{n+1})(1z_1z_{n+2})\mathrm{}(1z_nz_{2n})}𝑑𝐳.$$

Here it is understood that each integral is over a circle with radius $`<1`$ centered at 0; all appearing radii should be different. We can separate, say, the last $`n`$ variables and obtain

$$H_n(t)=\frac{1}{(2\pi i)^n}\mathrm{}(z_1\mathrm{}z_n)^{t1}\left(\frac{1}{2\pi i}\frac{z^{t1}}{(1z_1z)\mathrm{}(1z_nz)}𝑑z\right)^n𝑑z_n\mathrm{}𝑑z_1.$$

We may choose the radius of the integration circle of the innermost integral to be smaller than the radii of the other integration paths. Then this innermost integral is easy to compute: It is equal to the residue at 0 of

$$\frac{1}{z^{t+1}(1z_1z)\mathrm{}(1z_nz)}$$

and, by the residue theorem, equal to the negative of the sum of the residues at $`z_1^1,\mathrm{},z_n^1`$. (Note that here we use the fact that $`t>0`$.) The residues at these simple poles are easily computed: the one, say, at $`1/z_1`$ can be calculated as

$`\underset{z1/z_1}{lim}\left(z{\displaystyle \frac{1}{z_1}}\right){\displaystyle \frac{1}{z^{t+1}(1z_1z)\mathrm{}(1z_nz)}}=\underset{z1/z_1}{lim}{\displaystyle \frac{z\frac{1}{z_1}}{1z_1z}}{\displaystyle \frac{1}{z^{t+1}(1z_2z)\mathrm{}(1z_nz)}}`$
$`={\displaystyle \frac{1}{z_1}}{\displaystyle \frac{z_1^{t+1}}{(1z_2/z_1)\mathrm{}(1z_n/z_1)}}={\displaystyle \frac{z_1^{t+n1}}{(z_1z_2)\mathrm{}(z_1z_n)}}.`$
This yields the starting point for our computations.

###### Theorem 2

For any distinct $`0<ϵ_1,\mathrm{},ϵ_n<1`$,

$$H_n(t)=\frac{1}{(2\pi i)^n}_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_n\right|=ϵ_n}(z_1\mathrm{}z_n)^{t1}\left(\underset{k=1}{\overset{n}{}}\frac{z_k^{t+n1}}{_{jk}(z_kz_j)}\right)^n𝑑z_n\mathrm{}𝑑z_1.$$

*Remark*. It can be proved from the form of the integrand that $`H_n`$ is indeed a polynomial in $`t`$: To compute the integral, one has to execute a (huge) number of limit computations, which yield “at worst” powers of $`t`$ (as a consequence of L’Hospital’s Rule). In fact, one can make this property more apparent by noticing that the expression in parenthesis is actually a polynomial; namely

$$H_n(t)=\frac{1}{(2\pi i)^n}_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_n\right|=ϵ_n}(z_1\mathrm{}z_n)^{t1}\left(\underset{m_1+\mathrm{}+m_n=t}{}z_1^{m_1}\mathrm{}z_n^{m_n}\right)^n𝑑z_n\mathrm{}𝑑z_1,$$

where the sum is over all ordered partitions of $`t`$. This formula can also be proved “more directly” combinatorially.333The authors thank Sinai Robins and Frank Sottile for their help in the proof of this equivalence and its combinatorial interpretation.

## 3 Small $`n`$ do not require a computer

We will now illustrate the computation of $`H_n`$ (and hence $`\mathrm{vol}_n`$) by means of Theorem 2 for $`n=3`$ and $`4`$. These calculations “by hand” give an idea what computational tricks one might use in tackling larger $`n`$ with the aid of a computer.

By the theorem,

$$H_3(t)=\frac{1}{(2\pi i)^3}(z_1z_2z_3)^{t1}\left(\frac{z_1^{t+2}}{(z_1z_2)(z_1z_3)}+\frac{z_2^{t+2}}{(z_2z_1)(z_2z_3)}+\frac{z_3^{t+2}}{(z_3z_1)(z_3z_2)}\right)^3𝑑𝐳.$$

We have to order the radii of the integration paths for each variable; we choose $`0<ϵ_3<ϵ_2<ϵ_1<1`$. We heavily use this fact after multiplying out the cubic: integrating, for example, the term

$$\frac{z_1^{t1}z_2^{t1}z_3^{2t+5}}{(z_3z_2)^3(z_3z_1)^3}$$

with respect to $`z_3`$ gives 0, as this function is analytic at the $`z_3`$-origin and $`\left|z_1\right|,\left|z_2\right|>ϵ_3`$. After exploiting this observation for all the terms stemming from the cubic, the only integrals surviving are

$$\frac{1}{(2\pi i)^3}\frac{z_1^{2t+5}z_2^{t1}z_3^{t1}}{(z_1z_2)^3(z_1z_3)^3}𝑑𝐳$$

and

$$\frac{3}{(2\pi i)^3}\frac{z_1^{t+3}z_2z_3^{t1}}{(z_1z_2)^3(z_1z_3)^2(z_2z_3)}𝑑𝐳.$$

The first integral factors and yields, again by residue calculus,

$`{\displaystyle \frac{1}{(2\pi i)^3}}{\displaystyle \frac{z_1^{2t+5}z_2^{t1}z_3^{t1}}{(z_1z_2)^3(z_1z_3)^3}𝑑𝐳}={\displaystyle \frac{1}{(2\pi i)^3}}{\displaystyle z_1^{2t+5}\left(\frac{z^{t1}}{(z_1z)^3}𝑑z\right)^2𝑑z_1}`$
$`={\displaystyle \frac{1}{2\pi i}}{\displaystyle z_1^{2t+5}\left(\frac{1}{2}(t1)(t2)z_1^{t3}\right)^2𝑑z_1}=\left({\displaystyle \genfrac{}{}{0pt}{}{t+2}{2}}\right)^2.`$
For the second integral, it is most efficient to integrate with respect to $`z_2`$ first:

$`{\displaystyle \frac{3}{(2\pi i)^3}}{\displaystyle \frac{z_1^{t+3}z_2z_3^{t1}}{(z_1z_2)^3(z_1z_3)^2(z_2z_3)}𝑑𝐳}={\displaystyle \frac{3}{(2\pi i)^2}}{\displaystyle \frac{z_1^{t+3}z_3^t}{(z_1z_3)^5}𝑑z_3𝑑z_1}`$
$`={\displaystyle \frac{3}{2\pi i}}{\displaystyle z_1^{t+3}\frac{1}{4!}(t)(t1)(t2)(t3)z_1^{t4}𝑑z_1}=3\left({\displaystyle \genfrac{}{}{0pt}{}{t+3}{4}}\right).`$
Adding up the last two lines gives finally

$$H_3(t)=\left(\genfrac{}{}{0pt}{}{t+2}{2}\right)^23\left(\genfrac{}{}{0pt}{}{t+3}{4}\right)=\frac{1}{8}t^4+\frac{3}{4}t^3+\frac{15}{8}t^2+\frac{9}{4}t+1,$$

which is equal to (3). To obtain the volume of $`_3`$, the leading term of $`H_3`$ has to be multiplied by the relative volume of the fundamental domain of the sublattice of $`^9`$ in the affine space spanned by $`_3`$. This volume is 9; hence

$$\mathrm{vol}_3=\frac{9}{8}.$$

In general, it is not hard to prove (see, for example, the appendix of ) that the relative volume of the fundamental domain of the sublattice of $`^{n^2}`$ in the affine space spanned by $`_n`$ is $`n^{n1}`$.

The number of integrals we have to evaluate to compute $`H_4`$ is only slightly higher. By Theorem 2,

$$H_4(t)=\frac{1}{(2\pi i)^4}_{\left|z_1\right|=ϵ_1}_{\left|z_2\right|=ϵ_2}_{\left|z_3\right|=ϵ_3}_{\left|z_4\right|=ϵ_4}(z_1z_2z_3z_4)^{t1}\left(\underset{k=1}{\overset{4}{}}\frac{z_k^{t+3}}{_{jk}(z_kz_j)}\right)^4𝑑𝐳.$$

Again we have a choice of ordering the radii; we use $`0<ϵ_4<ϵ_3<ϵ_2<ϵ_1<1`$. After multiplying out the quartic, we have to calculate five integrals; their evaluation—again straightforward by means of the residue theorem—is as follows. As before, we can ‘save’ computation effort by choosing a particular order with which we integrate and by factoring an integral if possible.

$`{\displaystyle \frac{1}{(2\pi i)^4}}{\displaystyle \frac{z_1^{3t+11}z_2^{t1}z_3^{t1}z_4^{t1}}{(z_1z_2)^4(z_1z_3)^4(z_1z_4)^4}𝑑𝐳}={\displaystyle \frac{1}{(2\pi i)^4}}{\displaystyle z_1^{3t+11}\left(\frac{z^{t1}}{(z_1z)^4}𝑑z\right)^3𝑑z_1}`$
$`=\left({\displaystyle \genfrac{}{}{0pt}{}{t+3}{3}}\right)^3,`$
$`{\displaystyle \frac{4}{(2\pi i)^4}}{\displaystyle \frac{z_1^{2t+8}z_2^2z_3^{t1}z_4^{t1}}{(z_1z_2)^4(z_1z_3)^3(z_1z_4)^3(z_2z_3)(z_2z_4)}𝑑𝐳}`$
$`={\displaystyle \frac{4}{(2\pi i)^4}}{\displaystyle \frac{z_1^{2t+8}z_2^2}{(z_1z_2)^4}\left(\frac{z^{t1}}{(z_1z)^3(z_2z)}𝑑z\right)^2𝑑z_1𝑑z_2}`$
$`={\displaystyle \frac{4}{(2\pi i)^2}}{\displaystyle \frac{z_1^{2t+8}z_2^{t+1}}{(z_1z_2)^7}\left(2\left(\genfrac{}{}{0pt}{}{t+2}{2}\right)\frac{z_1^{t3}}{z_1z_2}+2(t+1)\frac{z_1^{t2}}{(z_1z_2)^2}+2\frac{z_1^{t1}}{(z_1z_2)^3}\frac{z_2^{t1}}{(z_1z_2)^3}\right)𝑑𝐳}`$
$`=8\left({\displaystyle \genfrac{}{}{0pt}{}{t+2}{2}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{7}}\right)+8(t+1)\left({\displaystyle \genfrac{}{}{0pt}{}{t+6}{8}}\right)+8\left({\displaystyle \genfrac{}{}{0pt}{}{t+7}{9}}\right)4\left({\displaystyle \genfrac{}{}{0pt}{}{2t+8}{9}}\right),`$
$`{\displaystyle \frac{4}{(2\pi i)^4}}{\displaystyle \frac{z_1^{2t+8}z_2^{t1}z_3^2z_4^{t1}}{(z_1z_2)^3(z_1z_3)^4(z_1z_4)^3(z_2z_3)(z_3z_4)}𝑑𝐳}`$
$`={\displaystyle \frac{4}{(2\pi i)^3}}{\displaystyle \frac{z_1^{2t+8}z_2^{t+1}z_4^{t1}}{(z_1z_2)^3(z_1z_4)^7(z_2z_4)}𝑑𝐳}`$
$`={\displaystyle \frac{4}{(2\pi i)^2}}{\displaystyle \frac{z_1^{2t+8}z_4^{t1}}{(z_1z_4)^7}\left(\left(\genfrac{}{}{0pt}{}{t+2}{2}\right)\frac{z_1^{t3}}{z_1z_4}+(t+1)\frac{z_1^{t2}}{(z_1z_4)^2}+\frac{z_1^{t1}}{(z_1z_4)^3}\right)𝑑𝐳}`$
$`=4\left(\left({\displaystyle \genfrac{}{}{0pt}{}{t+2}{2}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{7}}\right)+(t+1)\left({\displaystyle \genfrac{}{}{0pt}{}{t+6}{8}}\right)+\left({\displaystyle \genfrac{}{}{0pt}{}{t+7}{9}}\right)\right),`$
$`{\displaystyle \frac{6}{(2\pi i)^4}}{\displaystyle \frac{z_1^{t+5}z_2^{t+5}z_3^{t1}z_4^{t1}}{(z_1z_2)^4(z_1z_3)^2(z_1z_4)^2(z_2z_3)^2(z_2z_4)^2}𝑑𝐳}`$
$`={\displaystyle \frac{6}{(2\pi i)^4}}{\displaystyle \frac{z_1^{t+5}z_2^{t+5}}{(z_1z_2)^4}\left(\frac{z^{t1}}{(z_1z)^2(z_2z)^2}𝑑z\right)^2𝑑z_1𝑑z_2}`$
$`={\displaystyle \frac{6}{(2\pi i)^2}}{\displaystyle \frac{z_1^{t+5}z_2^{t+5}}{(z_1z_2)^4}\left((t+1)^2\frac{z_2^{2t4}}{(z_1z_2)^4}4(t+1)\frac{z_2^{2t3}}{(z_1z_2)^5}+4\frac{z_2^{2t2}}{(z_1z_2)^6}\right)𝑑𝐳}`$
$`=6(t+1)^2\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{7}}\right)24(t+1)\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{8}}\right)+24\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{9}}\right),`$
$`{\displaystyle \frac{12}{(2\pi i)^4}}{\displaystyle \frac{z_1^{t+5}z_2^2z_3^2z_4^{t1}}{(z_1z_2)^3(z_1z_3)^3(z_1z_4)^2(z_2z_3)^2(z_2z_4)(z_3z_4)}𝑑𝐳}`$
$`={\displaystyle \frac{12}{(2\pi i)^3}}{\displaystyle \frac{z_1^{t+5}z_2^2z_4^{t+1}}{(z_1z_2)^3(z_1z_4)^5(z_2z_4)^3}𝑑𝐳}`$
$`={\displaystyle \frac{12}{(2\pi i)^2}}{\displaystyle \frac{z_1^{t+5}z_4^{t+1}}{(z_1z_4)^5}\left(\frac{1}{(z_1z_4)^3}+6\frac{z_4}{(z_1z_4)^4}+6\frac{z_4^2}{(z_1z_4)^5}\right)𝑑𝐳}`$
$`=12\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{7}}\right)72\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{8}}\right)72\left({\displaystyle \genfrac{}{}{0pt}{}{t+5}{9}}\right).`$
Adding them all up gives

$$H_4(t)=\frac{11}{11340}t^9+\frac{11}{630}t^8+\frac{19}{135}t^7+\frac{2}{3}t^6+\frac{1109}{540}t^5+\frac{43}{10}t^4+\frac{35117}{5670}t^3+\frac{379}{63}t^2+\frac{65}{18}t+1$$

and hence

$$\mathrm{vol}_4=4^3\frac{11}{11340}=\frac{176}{2835}.$$

## 4 Larger $`n`$ do

As we have seen in the examples, after multiplying out the integrand of Theorem 2, many of the terms do not contribute to the integral. The following proposition will provide us with a general statement to that effect.

For a rational function $`f`$ in $`n`$ variables $`z_j`$ we use the notation $`d_r(f)`$ for the degree of $`f`$ in the variables $`z_1,\mathrm{},z_r`$.

###### Proposition 3

Suppose $`p_1,\mathrm{},p_n`$ are integers, $`q_{jk}`$ are nonnegative integers ($`1j<kn`$), $`1>ϵ_1>\mathrm{}>ϵ_n>0`$, and

$$f(z_1,\mathrm{},z_n)=\frac{_{1jn}z_j^{p_j}}{_{1j<kn}(z_jz_k)^{q_{jk}}}.$$

If $`1rn`$ and $`d_r(f)<r`$ then

$$_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_n\right|=ϵ_n}f(𝐳)𝑑𝐳=0.$$

*Proof*. We need only show that

$$_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_r\right|=ϵ_r}f(𝐳)𝑑𝐳=0,$$

for then we can integrate over all $`n`$ variables by first integrating over $`z_1\mathrm{}z_r`$.

If $`r=1`$ then $`f`$, considered as a function of $`z_1`$, is analytic outside the circle $`\left|z_1\right|=ϵ_1`$ and has zero residue at infinity since its degree is less than $`1`$, so $`()`$ follows.

We continue by induction, so suppose $`r>1`$ and $`()`$ is true for smaller $`r`$. Suppose that $`d_r(f)<r`$; we may assume $`d_{r1}(f)(r1)`$. Note that $`d_r(f)d_{r1}(f)+d`$ where $`d`$ is the degree of $`f`$ in the single variable $`z_r`$ (the discrepancy is the sum of the exponents $`q_{jr}`$ for $`j<r`$). Hence $`dd_r(f)d_{r1}(f)<1`$. We consider $`f`$ as a function of $`z_r`$ and apply the residue theorem to the region outside the circle $`\left|z_r\right|=ϵ_r`$. As above $`f`$ has zero residue at infinity, so we only need to consider the residues at the poles $`z_j`$ for $`j<r`$. Evaluating these residues converts the integral of $`f`$ into a (possibly huge) linear combination of integrals of functions of the same form as $`f`$, but in the $`n1`$ variables $`z_1,\mathrm{},z_{r1},z_{r+1},\mathrm{},z_n`$. If $`g`$ is any one of these functions then we easily calculate $`d_{r1}(g)=d_r(f)+1<(r1)`$, and, by the induction hypothesis, the integral of $`g`$ is zero. $`\mathrm{}`$

From this proposition we obtain the starting point for our ‘algorithm.’

###### Corollary 4

For $`1>ϵ_1>\mathrm{}>ϵ_n>0`$ and $`t0`$, $`H_n(t)=`$

$$\frac{1}{(2\pi i)^n}_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_n\right|=ϵ_n}(z_1\mathrm{}z_n)^{t1}\underset{m_1+\mathrm{}+m_n=n}{}^{}\left(\genfrac{}{}{0pt}{}{n}{m_1,\mathrm{},m_n}\right)\underset{k=1}{\overset{n}{}}\left(\frac{z_k^{t+n1}}{_{jk}(z_kz_j)}\right)^{m_k}d𝐳,$$

where $`^{}`$ denotes that we only sum over those $`n`$-tuples of non-negative integers satisfying $`m_1+\mathrm{}+m_n=n`$ and $`m_1+\mathrm{}+m_r>r`$ if $`1r<n`$.

*Remark*. The condition on $`m_1,\mathrm{},m_n`$ can be visualized through lattice paths from $`(0,0)`$ to $`(n,n)`$ using the steps $`(1,m_1),(1,m_2),\mathrm{},(1,m_n)`$. The condition means that these paths stay strictly above the diagonal (except at the start and end).444The authors thank Lou Billera for pointing out this lattice-path interpretation.

*Proof*. By Theorem 2,

$`H_n(t)`$ $`={\displaystyle \frac{1}{(2\pi i)^n}}{\displaystyle \mathrm{}(z_1\mathrm{}z_n)^{t1}\left(\underset{k=1}{\overset{n}{}}\frac{z_k^{t+n1}}{_{jk}(z_kz_j)}\right)^n𝑑𝐳}`$
$`={\displaystyle \frac{1}{(2\pi i)^n}}{\displaystyle \mathrm{}(z_1\mathrm{}z_n)^{t1}\underset{m_1+\mathrm{}+m_n=n}{}\left(\genfrac{}{}{0pt}{}{n}{m_1,\mathrm{},m_n}\right)\underset{k=1}{\overset{n}{}}\left(\frac{z_k^{t+n1}}{_{jk}(z_kz_j)}\right)^{m_k}d𝐳}.`$
We select a partition $`m_1+\mathrm{}+m_n=n`$ and rewrite the corresponding integrand in the language of Proposition 3:

$$p_j=(n1)m_j+t(m_j1)1,q_{jk}=m_j+m_k.$$

Now suppose $`1r<n`$. The degree of the denominator in $`z_1,\mathrm{},z_r`$ is

$`{\displaystyle \underset{j=1}{\overset{r}{}}}{\displaystyle \underset{k=j+1}{\overset{n}{}}}(m_j+m_k)=(n1){\displaystyle \underset{j=1}{\overset{r}{}}}m_j+r{\displaystyle \underset{j=r+1}{\overset{n}{}}}m_j=(n1){\displaystyle \underset{j=1}{\overset{r}{}}}m_j+nrr{\displaystyle \underset{j=1}{\overset{r}{}}}m_j`$
$`=(n1){\displaystyle \underset{j=1}{\overset{r}{}}}m_j+nrr{\displaystyle \underset{j=1}{\overset{r}{}}}(m_j1)r^2.`$
We subtract this from $`_{j=1}^rp_j`$ to get

$$d_r(f)=(t+r)\underset{j=1}{\overset{r}{}}(m_j1)r(nr)r.$$

Here $`t`$ is non-negative and $`r(nr)`$ is positive, so Proposition 3 implies that the integral is zero unless $`_{j=1}^r(m_j1)>0`$. $`\mathrm{}`$

Theoretically, Corollary 4 tells us what we have to do to compute $`H_n`$. Thus the integrals we computed in our calculations of $`H_3`$ and $`H_4`$ are exactly the non-zero integrals according to this result. For practical purposes, however, the statement is almost worthless for larger values of $`n`$. The first problem is that the number of terms in the sum equals the $`(n1)^{\text{th}}`$ Catalan number

$$\frac{1}{n}\left(\genfrac{}{}{0pt}{}{2(n1)}{n1}\right)=\frac{(2n2)!}{n!(n1)!},$$

which grows exponentially with $`n`$. Another slippery point is the evaluation of each integral. As we have seen in the examples, and as can be easily seen for the general case, we can compute each integral step by step one variable at a time. However, this means at each step we convert a rational function into a sum of rational functions (of one variable less) by means of the residue theorem. Again, this means that the number of (single-variable) integrals we have to compute grows immensely as $`n`$ increases. In fact, if we just ‘feed’ the statement of Corollary 4 into a computer and tell it to integrate each summand, say, starting with $`z_1`$, then $`z_2`$, and so on, the computation time explodes once one tries $`n=7`$ or $`8`$. We feel that computationally this is as involved as calculating a sufficient number of values of $`H_n`$ and then interpolating this polynomial. However, complex analysis allows us some shortcuts which turn out to speed up the computation by a huge factor and which make Corollary 4 valuable, even from a computational perspective. These ‘tricks’ all showed up already in the examples and include

1. realizing when a function is analytic at the $`z_k`$-origin,
2. trying to choose the most efficient order of variables to integrate (based on estimating how many terms will be generated by the residue calculations, for each available variable), and
3. factoring the integral if some of the variables appear in a symmetric fashion.

Finally, if we are only interested in the volume of $`_n`$, we may also be

1. suppressing a particular integral if it does not contribute to the leading term of $`H_n`$.

It is worth noting that each of these computational ‘speed-ups’ decreases the total computation time substantially. By applying them to Corollary 4, we implemented a C$`++`$ program for the specific functions we have to integrate to compute $`H_n`$. We were able to verify all previously known polynomials ($`n8`$) and to compute $`\mathrm{vol}_9`$, $`\mathrm{vol}_{10}`$, and $`H_9`$. The results of our calculations, including the polynomials $`H_n`$ for $`n9`$, and the source code for our program are available at www.math.binghamton.edu/dennis/Birkhoff. The following table gives the volumes together with the respective computing time (on a 1GHz PC running under Linux). Note that computing the full polynomial $`H_n`$ takes longer, as we cannot make use of shortcut #4. In fact, the computation of $`H_9`$ took about 325 days of computer time, although the elapsed time was only about two weeks since we distributed the calculations among a number of machines (8 – 40, depending on availability).

| $`n`$ | $`\mathrm{vol}_n`$ | time |
| --- | --- | --- |
| 1 | 1 | $`<.01`$ sec |
| 2 | 2 | $`<.01`$ sec |
| 3 | $`\text{}\frac{9}{8}`$ | $`<.01`$ sec |
| 4 | $`\text{}\frac{176}{2835}`$ | $`<.01`$ sec |
| 5 | $`\text{}\frac{23590375}{167382319104}`$ | $`<.01`$ sec |
| 6 | $`\text{}\frac{9700106723}{131928199603210^6}`$ | $`.18`$ sec |
| 7 | $`\text{}\frac{77436678274508929033}{13730296368223523839986892810^8}`$ | $`15`$ sec |
| 8 | $`\text{}\frac{5562533838576105333259507434329}{1258903626009547795008148094269333980330892810^{10}}`$ | $`54`$ min |
| 9 | $`\text{}\frac{559498129702796022246895686372766052475496691}{9269262340995263649896514671280698429605195132920241960610847715334553610^{14}}`$ | $`317`$ hr |

After this article was submitted for publication, we used the idle time on our departmental machines to compute the volume of $`_{10}`$. After a computation time of 6160 days, or almost 17 years (again scaled to a 1GHz processor), we obtained

$`\mathrm{vol}_{10}=`$

$$\frac{727291284016786420977508457990121862548823260052557333386607889}{828160860106766855125676318796872729344622463533089422677980721388055739956270293750883504892820848640000000}.$$

We hope that the reader will take this immense computing time as a challenge to improve our algorithms, and work towards $`n=11`$.

## 5 An outlook towards transportation polytopes

The Birkhoff polytopes are special cases of transportation polytopes, which are defined below. The study of this class of polytopes, which are naturally at least as fascinating as the Birkhoff polytopes, was motivated by problems in linear programming; for combinatorial properties see, for example, . The goal of this section is to show how our methods can be applied in this more general setting.

Fix positive real numbers $`a_1,\mathrm{},a_m,b_1,\mathrm{},b_n`$ such that $`a_1+\mathrm{}+a_m=b_1+\mathrm{}+b_n`$. Let $`𝐚=(a_1,\mathrm{},a_m),𝐛=(b_1,\mathrm{},b_n)`$.

$$𝒯_{𝐚,𝐛}=\{\left(\begin{array}{ccc}x_{11}& \mathrm{}& x_{1n}\\ \mathrm{}& & \mathrm{}\\ x_{m1}& \mathrm{}& x_{mn}\end{array}\right)^{mn}:x_{jk}0,\begin{array}{c}_kx_{jk}=a_j\text{ for all }1jm,\hfill \\ _jx_{jk}=b_k\text{ for all }1kn\hfill \end{array}\}$$

is the set of *solutions to the transportation problem with parameters* $`𝐚,𝐛`$. It is a convex polytope of dimension $`(m1)(n1)`$ in $`^{mn}`$—hence we refer to $`𝒯_{𝐚,𝐛}`$ as a *transportation polytopes*. Note that $`_n=𝒯_{\mathrm{𝟏},\mathrm{𝟏}}`$ where $`\mathrm{𝟏}=(1,\mathrm{},1)^n`$. We will be exclusively interested in the case when $`a_1,\mathrm{},a_m,b_1,\mathrm{},b_n`$ are integers. One reason for this is that the vertices of $`𝒯_{𝐚,𝐛}`$ are then integral. Let

$$T(𝐚,𝐛)=T(a_1,\mathrm{},a_m,b_1,\mathrm{},b_n)=\mathrm{\#}\left(𝒯_{𝐚,𝐛}^{mn}\right)$$

denote the number of integer points in the transportation polytope with integral parameters $`𝐚,𝐛`$. Geometrically, each of these parameters determines the position of a hyperplane bounding the polytope $`𝒯_{𝐚,𝐛}`$. It is well known that $`T(𝐚,𝐛)`$ is a piecewise-defined polynomial in $`a_1,\mathrm{},a_m,b_1,\mathrm{},b_n`$. The regions in which $`T(𝐚,𝐛)`$ is a polynomial depends on the normal cones of the polytopes involved; for details, we refer to .

Moreover, one can derive results for the counting function $`T(𝐚,𝐛)`$ which are ‘higher-dimensional’ extensions of (4), by application of the main theorem in , which in turn is a generalization of the Ehrhart-Macdonald reciprocity theorem (2). To state this theorem, we need to define the following integer-point counting function. Suppose $`𝒫`$ is a rational polytope given in the form (5); denote by $`𝒫(𝐭)=\{𝐱^d:\mathrm{𝐀𝐱}=𝐭\}`$ a polytope which we obtain from $`𝒫=𝒫(𝐛)`$ by translating (some of) its bounding hyperplanes. (Classical Ehrhart dilation is the special case $`𝐭=t𝐛`$.) Let

$$_𝒫(𝐭)=\mathrm{\#}\left(𝒫(𝐭)^d\right).$$

\[3, Theorem 4\] states that $`_𝒫`$ and $`_𝒫^{}`$ are piecewise-defined (multivariable) quasipolynomials satisfying

$$_𝒫(𝐭)=(1)^{dim𝒫}_𝒫^{}(𝐭).$$

As easily as (4) follows from (2), this reciprocity theorem yields a symmetry result for the transportation counting function. Let $`\mathrm{𝟏}_d`$ denote the $`d`$-dimensional vector all of whose entries are one. Then

$$T(𝐚n\mathrm{𝟏}_m,𝐛m\mathrm{𝟏}_n)=(1)^{mn1}T(𝐚,𝐛).$$

We finally turn to the problem of writing $`T(𝐚,𝐛)`$ in form of an integral. As with $`_n`$, we can view $`𝒯_{𝐚,𝐛}`$ as given in the form of (5) and apply the philosophy of Theorem 1 to obtain

$$T(𝐚,𝐛)=\frac{1}{(2\pi i)^{m+n}}\mathrm{}\frac{z_1^{a_11}\mathrm{}z_m^{a_m1}w_1^{b_11}\mathrm{}w_n^{b_n1}}{_{\genfrac{}{}{0pt}{}{1jm}{1kn}}\left(1z_jw_k\right)}𝑑𝐳𝑑𝐰.$$

Again the integral with respect to each (one-dimensional) variable is over a circle with radius $`<1`$ centered at 0, and all appearing radii are different. As in the Birkhoff case, we can separate, say, the $`𝐰`$-variables and obtain

$$T(𝐚,𝐛)=\frac{1}{(2\pi i)^{m+n}}\mathrm{}z_1^{a_11}\mathrm{}z_m^{a_m1}\underset{k=1}{\overset{n}{}}\frac{w_k^{b_k1}}{_{j=1}^m\left(1z_jw_k\right)}𝑑w_k𝑑𝐳.$$

And as with the Birkhoff polytope, the innermost integrals are easy to compute by means of the residue theorem. This yields the following statement which generalizes Theorem 2.

###### Theorem 5

For any distinct $`0<ϵ_1,\mathrm{},ϵ_m<1`$,

$$T(𝐚,𝐛)=\frac{1}{(2\pi i)^m}_{\left|z_1\right|=ϵ_1}\mathrm{}_{\left|z_m\right|=ϵ_m}z_1^{a_11}\mathrm{}z_m^{a_m1}\underset{k=1}{\overset{n}{}}\underset{j=1}{\overset{m}{}}\frac{z_j^{b_k+m1}}{_{lj}(z_jz_l)}d𝐳.$$

*Remark*. As with Theorem 2, it can be proved from the form of the integrand that $`T(𝐚,𝐛)`$ is indeed a piecewise-defined polynomial in $`a_1,\mathrm{},a_m,b_1,\mathrm{},b_n`$.

Theoretically we could now use this theorem to produce formulas for $`T(𝐚,𝐛)`$ just as we did for $`H_n(t)`$. There is one major difference: $`T(𝐚,𝐛)`$ is only a *piecewise-defined* polynomial. In fact, we can see this from the form of the integral in Theorem 5: whether we will get a nonzero contribution at a certain step in the computation depends heavily on the relationship between $`a_1,\mathrm{},a_m,b_1,\mathrm{},b_n`$.

The fact that the counting function $`T(𝐚,𝐛)`$ is of a somewhat more delicate nature naturally has computational consequences. We feel that providing any general results on this function would go beyond the scope of this article and will hopefully be the subject of a future project. On the other hand, we adjusted our algorithm to compute values of $`T(𝐚,𝐛)`$ for three (fixed) pairs $`(𝐚,𝐛)`$, which have been previously computed by Mount  and DeLoera and Sturmfels . The first example is

$`T((3046,5173,6116,10928),(182,778,3635,9558,11110))=`$
$`23196436596128897574829611531938753`$
The authors reported their computation took 20 minutes (Mount) or 10 minutes (DeLoera–Sturmfels) . We computed this number in 0.2 seconds, based on Theorem 5. A similar phenomenon happens with

$`T((338106,574203,678876,1213008),(20202,142746,410755,1007773,1222717))=`$
$`316052820930116909459822049052149787748004963058022997262397`$
and

$`T((30201,59791,70017,41731,58270),(81016,68993,47000,43001,20000))=`$
$`24640538268151981086397018033422264050757251133401758112509495633028`$
which reportedly took 35 minutes/10 days with the DeLoera–Sturmfels algorithm , 0.3/2.9 seconds with ours. These timing comparisons ignore differences in machine speeds and implementation of the algorithms, but suggest that our methods are considerably more efficient.

Acknowledgements. We would like to thank the referees for carefully reading through our paper and software, and for their many helpful comments.

Department of Mathematical Sciences
State University of New York
Binghamton, NY 13902-6000
matthias@math.binghamton.edu
dennis@math.binghamton.edu


format: text

text:
# 1 Introduction

## 1 Introduction

Einstein’s special relativity is based on the two postulates:

1. All physical laws are invariant with respect to inertial reference frames;
2. The speed of light is a universal constant with respect to all inertial reference frames.

The first postulate states that physical phenomena have the same appearance in all inertial reference frames. One obtains the same results for all measurable quantities and clocks are synchronized by means of two-way speed of light measurements. The second postulate states that there is no preferred reference frame and that globally, physical laws appear the same in all inertial frames. For a description of a consistent special relativity theory, it is sufficient to just adopt postulate (1) . However, the speed of light is still a constant but this constant can vary from one reference frame to another. The spacetime coordinates are defined so that the one-way speed of light is constant.

In a Galilean spacetime the physical existence of an absolute time is postulated. The global concept of past, present and future is the same in all inertial frames. Simultaneous events in spacetime can occur in all reference frames. There exists a unique separation of past and future events. The speed of light is only constant in the Newtonian rest frame, which is a preferred frame of reference. Electromagnetic waves are limited not to exceed the speed of light.

The first postulate of Einstein’s special relativity is well established by local experiments, whereas the second one is not. In contrast to Galilean relativity with its absolute time, simultaneity is not an absolute concept, but a relative one, depending on the motion of the observer. However, the first postulate of Einstein’s special relativity may not be so well established at the scale of cosmology or at short distances.

Although the local Lorentz invariance of the laws of physics is a mathematical and physically attractive idea, there have been several reasons to be concerned about its generality at a fundamental level. It should be emphasized that there is no experiment that determines the one-way speed of light, for this would suppose that we can synchronize physical clocks by some way other than by finite speed of light signals. The main reason to be concerned with the generality of Einstein’s special relativity is the conflict between the two cornerstones of modern physics, General relativity (GR) and quantum mechanics.

In GR the concept of time disappears as a physically important quantity, because of the reparameterization invariance of the theory. This has led to considerable research in attempts to understand how to construct a consistent theory of quantum gravity . When we attempt to “quantize” spacetime, we find it difficult to maintain our classical ideas about local causality and special relativity associated with the Minkowski lightcone. Moreover, in quantum mechanics time is an external parameter, whereas in special relativity space and time are on an equal footing.

In cosmology, a comoving coordinate time does appear in the framework of a Friedmann, Robertson and Walker (FRW) spacetime, in which the time $`t`$ appears as a universal time which measures the lifetime of the universe. However, the FRW metric is subject to a diffeomorphism reparameterization transformation to non-comoving coordinates and the notion of a “universal” absolute time disappears.

The idea that a variable speed of light can solve the initial value problems of standard big-bang cosmology has received considerable attention . From our considerations of relativity theory, we learn that it is not possible to have a varying speed of light, without modifying Einstein’s formulation of special relativity and GR.

In our first treatment of a varying speed of light theory , we introduced a formalism for spontaneously breaking the local Lorentz invariance of the vacuum within GR. The symmetry of the homogeneous Lorentz group $`SO(3,1)`$ was spontaneously broken to the Galilean group $`O(3)`$ in the early universe. This introduced the concept of an “absolute” time and a formulation of quantum gravity that did not find itself in conflict with the concept of time in quantum mechanics. As the universe expanded, a phase transition occurred that reinstated the four-dimensional symmetry of the homogenous Lorentz group. However, the speed of light was only allowed to change discontinuously from one constant value near the beginning of the universe to a much smaller value corresponding to the presently measured value of the speed of light $`c_0=\mathrm{299\hspace{0.17em}792\hspace{0.17em}458}m/s`$. A consistent treatment of a varying speed of light within a gravitational theory was lacking. In order to remedy this problem, a bimetric gravity theory (BGT) was constructed, which provided a consistent diffeomorphism invariant formalism for descibing varying speed of light phenomena .

In the following, we shall focus on the scalar-tensor bimetric theory . The theory can be the preliminary stage for a quantum gravity theory that does not conflict with some of the basic properties of quantum mechanics. The physical effects of having two different spacetime metrics may be expected to be observationally detectable only at cosmological scales, when dark matter and dark energy dominate the gravitational effects, or at short distances. The BGT satisfies a matter conservation law that assures that matter test particles travel along geodesics and do not violate the weak equivalence principle, which is observationally well-tested . Recently, the predictions of a primordial fluctuation spectrum were obtained from the scalar-tensor bimetric gravity theory and compared to the CMB data .

The observations of supernovae (SNe Ia) at red shifts $`0.35<z<1.75`$ have shown that they appear to be farther away and therefore fainter than is to be expected from the standard decelerating model of the universe . This has led to the introduction of dark energy models which generate an accelerating expansion of the universe, based on either a cosmological constant or a form of quintessence . The major problem with these models is that they lead to a coincidence problem, and a fine-tuned tiny mass smaller than the Hubble parameter, $`m_Q10^{33}`$ eV, and small couplings to visible matter which must satisfy fifth-force constraints . Moreover, an eternally accelerating universe will produce de Sitter event horizons that cause problems for string theory and quantum gravity theories .

Recently, an alternative explanation of the SNe Ia observations was proposed without the acceleration caused by dark energy, based on the idea of flavor oscillations . An axion and a photon mixing causes an attenuation of the photon flux from distant systems. The model postulates an axion mass $`m_A10^{16}`$ eV and a coupling with a mass scale $`M_S10^{11}`$ GeV. A higher-dimensional model with leakage into the higher dimension has also been proposed as an alternative to the standard cosmological constant or quintessence models .

In the following, we shall consider the problem of the dimming of supernovae within the BGT in the variable speed of light (VSL) frame.

The Friedmann equation obtained from the field equations of BGT in the VSL frame leads to dark matter and dark energy components which give the experimentally observed result $`\mathrm{\Omega }_{\mathrm{Tot}}1`$, without generating an accelerating expansion of the universe. The observed dimming of the SNe Ia is affected by an increase of the red shift by 10%-15%, due to an increase in the speed of light in the past, which increases the luminosity distance to the supernovae.

We consider the consequences for other fundamental “constants” of nature such as Newton’s constant $`G`$ and the fine-structure constant $`\alpha =e^2/\mathrm{}c`$ and discuss the problem of cosmological horizons in BGT.

## 2 Bimetric Scalar-Tensor Gravity

Our scalar-tensor bimetric gravity theory (BGT)  is described by two metrics called the “matter metric” $`\widehat{g}_{\mu \nu }`$ and the “gravitational metric” $`g_{\mu \nu }`$, respectively, connected by the “biscalar” field $`\varphi `$ through the equation

$$\widehat{g}_{\mu \nu }=g_{\mu \nu }+B_\mu \varphi _\nu \varphi ,$$
(1)
where the constant $`B`$ has dimensions of $`[\mathrm{length}]^2`$ and is chosen to be positive. The inverse metrics $`\widehat{g}^{\mu \nu }`$ and $`g^{\mu \nu }`$ satisfy

$$\widehat{g}^{\mu \alpha }\widehat{g}_{\nu \alpha }=\delta _{}^{\mu }{}_{\nu }{}^{},g^{\mu \alpha }g_{\nu \alpha }=\delta _{}^{\mu }{}_{\nu }{}^{}.$$
(2)
We shall restrict ourselves to the physical non-degenerate case: $`\mathrm{Det}(\widehat{g}_{\mu \nu })0`$. The metric $`\widehat{g}_{\mu \nu }`$ is used to construct the matter action and can be said to be the geometry on which matter fields propagate. It is the combination of the gravitational metric and the biscalar field that we consider as being the gravitational fields of our theory. One of the satisfactory features of the theory is that the formalism is generally covariant (diffeomorphism invariant), which guarantees that basic properties of the theoretical structure are consistent.

The model that we introduced in  consisted in a self-gravitating scalar field coupled to matter through the matter metric (1), with the action

$$S=S_{\mathrm{grav}}+S_\varphi +\widehat{S}_\mathrm{M},$$
(3)
where

$$S_{\mathrm{grav}}=\frac{1}{\kappa }𝑑\mu (R[g]+2\mathrm{\Lambda }),$$
(4)
$`\kappa =16\pi G/c_0^4`$, $`\mathrm{\Lambda }`$ is the cosmological constant, and we employ a metric with signature $`(+,,,)`$. We will write, for example, $`d\mu =d^4x\sqrt{g}`$ and $`\mu =\sqrt{g}`$ for the metric density related to the gravitational metric $`g_{\mu \nu }`$, and similar definitions of $`d\widehat{\mu }`$ and $`\widehat{\mu }`$ in terms of the matter metric $`\widehat{g}_{\mu \nu }`$. The minimally-coupled scalar field action is given by

$$S_\varphi =\frac{1}{\kappa }𝑑\mu \left[\frac{1}{2}g^{\mu \nu }_\mu \varphi _\nu \varphi V(\varphi )\right],$$
(5)
where the scalar field $`\varphi `$ has been chosen to be dimensionless. The energy-momentum tensor for the scalar field that we will use is given by

$$T_\varphi ^{\mu \nu }=\frac{1}{\kappa }\left[g^{\mu \alpha }g^{\nu \beta }_\alpha \varphi _\beta \varphi \frac{1}{2}g^{\mu \nu }g^{\alpha \beta }_\alpha \varphi _\beta \varphi +g^{\mu \nu }V(\varphi )\right],$$
(6)
and is the variation of the scalar field action with respect to the gravitational metric:

$$\frac{\delta S_\varphi }{\delta g_{\mu \nu }}=\frac{1}{2}\mu T_\varphi ^{\mu \nu }.$$
(7)
We shall use the metric (1) to construct the matter action $`\widehat{S}_\mathrm{M}`$, resulting in the identification of $`\widehat{g}_{\mu \nu }`$ as the metric that provides the arena on which matter fields interact. The matter action $`\widehat{S}_\mathrm{M}[\psi ^I]=\widehat{S}_\mathrm{M}[\widehat{g},\psi ^I]`$, where $`\psi ^I`$ represents all the matter fields in spacetime. The energy-momentum tensor

$$\frac{\delta S_\mathrm{M}}{\delta \widehat{g}_{\mu \nu }}=\frac{1}{2}\widehat{\mu }\widehat{T}^{\mu \nu },$$
(8)
satisfies the conservation laws

$$\widehat{}_\nu \left[\widehat{\mu }\widehat{T}^{\mu \nu }\right]=0.$$
(9)
This follows as a consequence of the matter field equations only . Here, $`\widehat{}_\mu `$ is the metric compatible covariant derivative determined by the matter metric: $`\widehat{}_\alpha \widehat{g}_{\mu \nu }=0`$.

The gravitational field equations are given by :

$$G^{\mu \nu }=\mathrm{\Lambda }g^{\mu \nu }+\frac{\kappa }{2}T_\varphi ^{\mu \nu }+\frac{\kappa }{2}\frac{\widehat{\mu }}{\mu }\widehat{T}^{\mu \nu }.$$
(10)
The scalar field wave equation is

$$_\mu ^\mu \varphi +V^{}(\varphi )\kappa \frac{\widehat{\mu }}{\mu }\widehat{T}^{\mu \nu }\widehat{}_\mu \widehat{}_\nu \varphi =0,$$
(11)
where $`_\mu `$ is the covariant derivative with respect to the metric $`g_{\mu \nu }`$.

From the definition (1), we obtain the inverses

$$\widehat{g}^{\mu \nu }=g^{\mu \nu }\frac{B}{I}^\mu \varphi ^\nu \varphi ,$$
(12)
and

$$g^{\mu \nu }=\widehat{g}^{\mu \nu }+\frac{B}{K}\widehat{}^\mu \varphi \widehat{}^\nu \varphi ,$$
(13)
where

$$I=1+Bg^{\mu \nu }_\mu \varphi _\nu \varphi ,K=1B\widehat{g}^{\mu \nu }_\mu \varphi _\nu \varphi .$$
(14)
It follows that $`IK=1`$ and we have defined $`^\mu \varphi =g^{\mu \nu }_\nu \varphi `$ and

$$\widehat{}^\mu \varphi =\widehat{g}^{\mu \nu }_\nu \varphi =K^\mu \varphi .$$
(15)
We will assume a perfect fluid form for the matter fields

$$\widehat{T}^{\mu \nu }=\left(\rho +\frac{p}{c_0^2}\right)\widehat{u}^\mu \widehat{u}^\nu p\widehat{g}^{\mu \nu },$$
(16)
where $`\widehat{g}_{\mu \nu }\widehat{u}^\mu \widehat{u}^\nu =1`$ and $`\widehat{u}^\mu =dx^\mu /d\widehat{s}`$.

A fundamental feature of the BGT is that there is a separate frame associated with each metric, and each frame or metric has its own light cone. Because of the biscalar field linkage between the two frames, the two light cones cannot be physically identical. If we choose a small patch of spacetime in the “gravitational” frame with the metric $`g_{\mu \nu }`$, then we can locally make the $`g_{\mu \nu }`$ metric equal the flat, Minkowski metric, i.e. $`g_{\mu \nu }=\eta _{\mu \nu }`$ where $`\eta _{\mu \nu }=\mathrm{diag}(1,1,1,1)`$. This yields

$$ds^2=\eta _{\mu \nu }dx^\mu dx^\nu .$$
(17)
However, from (1) we obtain

$$d\widehat{s}^2=\eta _{\mu \nu }dx^\mu dx^\nu +B_\mu \varphi _\nu \varphi dx^\mu dx^\nu .$$
(18)
If we perform a local Lorentz transformation

$$x^{}_{}{}^{}\mu =\mathrm{\Lambda }_{}^{\alpha }{}_{\nu }{}^{}x^\nu ,$$
(19)
where $`\mathrm{\Lambda }_{}^{\alpha }{}_{\nu }{}^{}`$ are constant tensor coefficients that satisfy

$$\mathrm{\Lambda }_{}^{\alpha }{}_{\mu }{}^{}\mathrm{\Lambda }_{\alpha \nu }=\eta _{\mu \nu },$$
(20)
then (17) and (18) both remain form invariant. This follows if we have $`V_\mu =_\mu \varphi `$, where $`V_\mu ^{}=\mathrm{\Lambda }_{\mu }^{}{}_{}{}^{\alpha }V_\alpha `$ under a Lorentz transformation. If the light cone equation $`ds^2=0`$ is satisfied, then $`d\widehat{s}^20`$ along the same light cone unless $`\varphi =0`$. However, $`d\widehat{s}^2=0`$ along an expanding light cone that encloses the light cone $`ds^2=0`$. We observe that whereas we can transform away the Christoffel symbols $`\mathrm{\Gamma }_{\mu \nu }^\lambda `$ at a spacetime point, we cannot transform away the scalar field $`\varphi `$ at this point, nor the Christoffel symbol $`\widehat{\mathrm{\Gamma }}_{\mu \nu }^\lambda `$.

If we choose the vacuum expectation value $`𝒱_\mu =<V_\mu >_00`$, then the gradient of the scalar field $`\varphi `$ acts as a spontaneous symmetry breaking field, so that it is not possible to maintain the local Lorentz invariance of the vacuum states of both $`\widehat{g}_{\mu \nu }`$ and $`g_{\mu \nu }`$ metrics. Thus, a preferred frame is chosen under the spontaneous symmetry breaking of the vacuum, since $`𝒱_\mu `$ picks out a preferred direction in spacetime. This frame can correspond to a timelike vector $`𝒱_\mu =(𝒱,0,0,0)`$, leading to a spontaneous breaking of the homogeneous Lorentz group $`SO(3,1)O(3)`$ i.e. only the symmetry group of rotations would be preserved. In GR there is only one rigid light cone and one metric and local Lorentz invariance is strictly maintained.

In the BGT developed in earlier papers, we assumed that in the VSL frame all particles move with the the speed of light $`c(t)`$. We can generalize this theory, so that each particle has its own light cone by adopting the metric definition:

$$\widehat{g}_{I\mu \nu }=g_{\mu \nu }+B_I_\mu \varphi _I_\nu \varphi _I,$$
(21)
where the label $`I`$ denotes the specific particle with its associated speed $`c_I(t)`$.

## 3 Bimetric Gravity Cosmology

Let us now consider a cosmological scenario, imposing homogeneity and isotropy on spacetime and writing the “gravitational frame” metric $`g_{\mu \nu }`$ in comoving form

$$ds^2=c_0^2dt^2R^2(t)\left[\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right],$$
(22)
where we employ a dimensionless radial variable $`r`$ and $`k=0,\pm 1`$ for flat, closed and open hyperbolic spatial topologies, respectively. Then, the metric in what we call the VSL frame takes the form

$$d\widehat{s}^2=c^2(t)dt^2R^2(t)\left[\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right],$$
(23)
where

$$c(t)=c_0I^{1/2},$$
(24)
and

$$I=1+\frac{B}{c_0^2}\dot{\varphi }^2.$$
(25)
An overdot indicates a derivative with respect to the time variable $`t`$.

The matter stress-energy tensor is

$$\widehat{T}^{00}=\frac{\rho }{I},\widehat{T}^{0i}=0,\widehat{T}^{ij}=\frac{p}{R^2}\gamma ^{ij},$$
(26)
where $`\gamma ^{ij}(i,j=1,2,3)`$ denotes the spatial metric. The conservation laws become

$$\dot{\rho }+3H\left(\rho +\frac{p}{c_0^2}\right)=0,$$
(27)
where $`H=\dot{R}/R`$ is the Hubble function.

We see that in the VSL frame the speed of light depends on time, whereas in the variable speed of gravitational waves (VSGW) frame, the speed of light is constant and the speed of gravitational waves depends on time. The latter frame is arrived at by writing the matter metric $`\widehat{g}_{\mu \nu }`$ in comoving form with constant $`c_0`$:

$$d\widehat{s}^2=c_0^2dt^2R^2(t)\left[\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right].$$
(28)
We now obtain

$$ds^2=v_g^2dt^2R^2(t)\left[\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right],$$
(29)
where $`v_g`$ is the speed of gravitational waves

$$v_g(t)=c_0K^{1/2},$$
(30)
and

$$K=1\frac{B}{c_0^2}\dot{\varphi }^2.$$
(31)
We see that the speed of gravitational waves $`v_g`$ in the VSGW frame is predicted for $`B>0`$ to be smaller than the presently measured speed of light $`c_0`$.

The laws of physics will be interpreted differently by observers, depending upon whether they perform experiments in the VSL frame or the VSGW frame.

The Friedmann equation in the VSL frame is given by

$$H^2+\frac{kc_0^2}{R^2}=\frac{8\pi G}{3I^{1/2}}\rho +\frac{1}{3}c_0^2\mathrm{\Lambda }+\frac{1}{6}\rho _\varphi ,$$
(32)
where

$$\rho _\varphi =\frac{1}{2}\dot{\varphi }^2+c_0^2V(\varphi ).$$
(33)
The remaining equation is

$$\frac{\ddot{R}}{R}=\frac{4\pi G}{3I^{1/2}}\left(\rho +3I\frac{p}{c_0^2}\right)+\frac{1}{3}c_0^2\mathrm{\Lambda }\frac{1}{12}(\rho _\varphi +3p_\varphi ),$$
(34)
where

$$p_\varphi =\frac{1}{2}\dot{\varphi }^2c_0^2V(\varphi ).$$
(35)
The scalar field equation (11) in the VSL frame is

$$\frac{1}{c_0^2}\left(1\frac{16\pi GB}{c_0^2I^{3/2}}\rho \right)\ddot{\varphi }+\frac{3}{c_0^2}H\dot{\varphi }\left(1+\frac{16\pi GB}{c_0^4I^{1/2}}p\right)+V^{}(\varphi )=0.$$
(36)
The Friedmann equation in the VSGW frame is given by

$$H^2+\frac{c_0^2kK}{R^2}=\frac{8\pi G}{3}K^{3/2}\rho +\frac{1}{3}c_0^2\mathrm{\Lambda }K+\frac{1}{6}\stackrel{~}{\rho }_\varphi .$$
(37)
where we have defined

$$\stackrel{~}{\rho }_\varphi =\frac{1}{2}\dot{\varphi }^2+c_0^2KV(\varphi ),\stackrel{~}{p}_\varphi =\frac{1}{2}\dot{\varphi }^2c_0^2KV(\varphi ).$$
(38)
The remaining equation is

$$\frac{\ddot{R}}{R}=\frac{4\pi G}{3}\left(K^{3/2}\rho +\frac{3}{c_0^2}K^{1/2}p\right)+\frac{1}{3}c_0^2\mathrm{\Lambda }K\frac{1}{12}(\stackrel{~}{\rho }_\varphi +3\stackrel{~}{p}_\varphi )+\frac{1}{2}\frac{\dot{K}}{K}H.$$
(39)
The scalar wave equation in this frame is given by

$$\frac{1}{c_0^2}\left(1\frac{16\pi GB}{c_0^2}K^{3/2}\rho \right)\ddot{\varphi }+\frac{3K}{c_0^2}H\dot{\varphi }\left(1+\frac{16\pi GB}{c_0^4}K^{1/2}p\right)+K^2V^{}(\varphi )=0.$$
(40)
## 4 Calculation of the Red Shift

By adopting the VSL frame and setting $`d\widehat{s}^2=0`$ in (23), we obtain for fixed $`\theta `$ and $`\varphi `$:

$$_{t_1}^{t_0}\frac{dtc(t)}{R(t)}=_0^{r_1}\frac{dr}{\sqrt{1kr^2}}=f(r_1)=\mathrm{const}.,$$
(41)
where $`t_1`$ denotes the time a light wave leaves a supernova or galaxy and $`t_0`$ denotes the time when it reaches us on earth. If the next wave crest leaves $`r_1`$ at time $`t_1+\delta t_1`$, it will arrive on earth at time $`t_0+\delta t_0`$ giving

$$_{t_1+\delta t_1}^{t_0+\delta t_0}\frac{dtc(t)}{R(t)}=f(r_1).$$
(42)
Subtracting (41) from (42) and considering that $`R`$ and $`c`$ vary little during the period of a light signal, we obtain

$$\frac{c(t_0)\delta t_0}{R(t_0)}=\frac{c(t_1)\delta t_1}{R(t_1)}.$$
(43)
The frequency at emission $`\nu _1`$ is then related to the observed frequency $`\nu _0`$ by

$$\frac{\nu _0}{\nu _1}=\frac{\delta t_1}{\delta t_0}=\frac{c(t_0)}{c(t_1)}\frac{R(t_1)}{R(t_0)}.$$
(44)
The red shift

$$z=\frac{\lambda _0\lambda _1}{\lambda _1}$$
(45)
is then given by

$$z=\frac{c(t_1)}{c(t_0)}\frac{R(t_0)}{R(t_1)}1,$$
(46)
where we have

$$\lambda (t)\nu (t)=c_0,$$
(47)
so that $`\lambda _0/\lambda _1=\nu _1/\nu _0`$. This follows from the relation

$$\widehat{k}^\mu \widehat{k}_\mu =\widehat{g}^{\mu \alpha }k_\alpha \widehat{g}_{\mu \beta }k^\beta =k_\beta k^\beta =0,$$
(48)
where we have used (2) and $`k^\mu =(\nu /c_0,1/\stackrel{}{\lambda })`$ and $`\lambda =|\stackrel{}{\lambda }|`$.

We see that if $`c(t_1)>c(t_0)`$, then the observed red shift will appear to be larger due to the speed of light in the past being bigger than the presently observed speed $`c(t_0)=c_0`$ in an expanding universe with $`R(t_0)>R(t_1)`$.

## 5 Apparent Luminosity and Luminosity Distance

The apparent luminosity $`l`$ is the power per unit mirror area in our VSL frame :

$$l\frac{P}{A}=L\left(\frac{c^2(t_0)}{c^2(t_1)}\right)\left(\frac{R^2(t_1)}{R^2(t_0)}\right)\left(\frac{1}{4\pi R^2(t_0)r_1^2}\right),$$
(49)
where $`L`$ is the absolute luminosity, $`A`$ is the proper area of the mirror and $`P`$ is the total power emitted at the source. The luminosity distance of a light source in Euclidean space is

$$d_L=\left(\frac{L}{4\pi l}\right)^{1/2},$$
(50)
which leads to the expression

$$d_L=\left(\frac{c(t_1)}{c(t_0)}\right)\left(\frac{R^2(t_0)r_1}{R(t_1)}\right).$$
(51)
Let us expand $`R(t)`$ as

$$R(t)=R(t_0)[1+H_0(tt_0)\frac{1}{2}q_0H_0^2(tt_0)^2+\mathrm{}],$$
(52)
where

$$H_0=\frac{\dot{R}(t_0)}{R(t_0)},q_0=\frac{\ddot{R}(t_0)R(t_0)}{\dot{R}^2(t_0)}.$$
(53)
We also expand $`c(t)`$ in the power series

$$c(t)=c(t_0)[1+D_0(tt_0)+\frac{1}{2}Q_0D_0^2(tt_0)^2+\mathrm{}],$$
(54)
where

$$D_0=\frac{\dot{c}(t_0)}{c(t_0)},Q_0=\frac{\ddot{c}(t_0)c(t_0)}{\dot{c}^2(t_0)}.$$
(55)
The negative signs of the coefficients $`D_0`$ and $`Q_0`$ are chosen to yield a decreasing speed of light as the universe expands. From (46), we get using (52) and (54) the red shift

$$z=(H_0D_0)(t_0t_1)+[(1+\frac{1}{2}q_0)H_0^2+\frac{1}{2}Q_0D_0^2D_0H_0](t_0t_1)^2+\mathrm{}$$
(56)
Inverting this equation we get

$$t_0t_1=\left(\frac{1}{H_0D_0}\right)z\left(\frac{1}{H_0D_0}\right)^3[(1+\frac{1}{2}q_0)H_0^2+\frac{1}{2}Q_0D_0^2D_0H_0]z^2+\mathrm{}.$$
(57)
From (41), we obtain

$$\left(\frac{c(t_0)}{R(t_0)}\right)[t_0t_1+\frac{1}{2}(H_0D_0)(t_0t_1)^2+..]=r_1+O(r_1^3).$$
(58)
This gives the result for $`r_1`$:

$$r_1=\left(\frac{c(t_0)}{R(t_0)}\right)\{\left(\frac{1}{H_0D_0}\right)z\left(\frac{1}{H_0D_0}\right)^3[(1+\frac{1}{2}q_0)H_0^2+\frac{1}{2}Q_0D_0^2D_0H_0]z^2$$

$$+\frac{1}{2}\left(\frac{1}{H_0D_0}\right)z^2+\mathrm{}\}.$$
(59)
From (46) and (51), we find

$$d_L=c(t_0)r_1(1+z)$$

$$=\left(\frac{c(t_0)}{H_0D_0}\right)\{z\left(\frac{1}{H_0D_0}\right)^2[(1+\frac{1}{2}q_0)H_0^2+\frac{1}{2}Q_0D_0^2D_0H_0$$

$$\frac{3}{2}(H_0D_0)^2]z^2+\mathrm{}\}.$$
(60)
In the limit that $`D_0=Q_0=0`$, we obtain the correct GR limit

$$d_L=\left(\frac{c(t_0)}{H_0}\right)[z+\frac{1}{2}(1q_0)z^2+\mathrm{}].$$
(61)
We now obtain the formula for the apparent luminosity

$$l\frac{L}{4\pi d_L^2}=\left(\frac{L}{4\pi z^2c^2(t_0)}\right)(H_0D_0)^2\{1+\left(\frac{1}{H_0D_0}\right)^2[(1+\frac{1}{2}q_0)H_0^2$$

$$+\frac{1}{2}Q_0D_0^2D_0H_0\frac{3}{2}(H_0D_0)^2]z^2+\mathrm{}\}.$$
(62)
## 6 Dimming of Supernovae

Let us write the Friedmann equation (32) as

$$\mathrm{\Omega }_{0m}+\mathrm{\Omega }_{0k}+\mathrm{\Omega }_{0\mathrm{\Lambda }}+\mathrm{\Omega }_{0\varphi }=1,$$
(63)
where

$$\mathrm{\Omega }_{0m}=\frac{8\pi G\rho _{0m}}{3H_0^2I_0^{1/2}},\mathrm{\Omega }_{0k}=\frac{c_0^2k}{R_0^2H_0^2},\mathrm{\Omega }_{0\mathrm{\Lambda }}=\frac{c_0^2\mathrm{\Lambda }}{3H_0^2},$$

$$\mathrm{\Omega }_{0\varphi }=\frac{\rho _{0\varphi }}{6H_0^2},$$
(64)
where $`\mathrm{\Omega }_0,\rho _0,\varphi _0`$ and $`H_0`$ denote the present values of the corresponding quantities.

Choosing the values

$$\mathrm{\Omega }_{0m}=0.28,\mathrm{\Omega }_{0k}=0,\mathrm{\Omega }_{0\mathrm{\Lambda }}=0,\mathrm{\Omega }_{0\varphi }=0.72,$$
(65)
we obtain a fit to the supernovae data with $`q_0=\frac{1}{2}`$ and $`D_0`$ and $`Q_0`$ chosen for a given red shift $`z`$ to give a 10-15% increase in the luminosity distance (5), corresponding to a 20-30% decrease in the apparent luminosity observed in the supernovae data . For a negligible value of the pressure $`p_m`$, we get from (34) with $`\mathrm{\Lambda }=0`$:

$$\frac{\ddot{R}}{R}=\left[\left(\frac{4\pi G}{3I_0^{1/2}}\right)\rho _{0m}+\frac{1}{12}(\rho _{0\varphi }+3p_{0\varphi })\right].$$
(66)
We see that the scalar field $`\varphi `$ does not produce an acceleration of the universe, when $`\rho _{0\varphi }>0`$ and $`p_{0\varphi }>0`$. The increased speed of light at the observed SNe Ia red shifts $`z`$ has dimmed the supernovae in accordance with the SNe Ia observations. We do not attempt here to explain why $`\mathrm{\Lambda }=0`$.

The Friedmann equation (32) can be rewritten in terms of an effective gravitational constant $`G_{\mathrm{eff}}=G/\sqrt{I(t)}`$, corresponding to a time-varying gravitational constant.

We see that the time variation of $`G`$ depends on the time derivative of the scalar field $`\varphi `$. The cosmic time evolution of $`G`$ should be of the order of the expansion rate of the universe, i.e. $`\dot{G}/GH_0`$, where $`H_0=100h\mathrm{km}s^1\mathrm{Mpc}^1=h\times 10^{10}\mathrm{yr}^1`$ . Current observations of the expansion of the universe yield $`h0.7`$. For $`z1`$ an estimate gives $`\mathrm{\Delta }t4.9\times 10^9\mathrm{yr}`$, so that we obtain $`\mathrm{\Delta }G/G0.3`$. Thus, in the red shift range $`z0.13`$, a 10-15% change in $`G`$ can be accomodated and be consistent with a 10-15% change in $`c(t)=c_0(1+\frac{B}{c_0^2}\dot{\varphi }^2)^{1/2}`$. However, as we approach the time of big-bang nucleosynthesis (BBN) at a red shift $`z10^910^{10}`$, we have $`\mathrm{\Delta }G/G10^910^{10}`$ and the variation of $`G`$ with time becomes much more restrictive. This requires that for red shifts larger than $`z3`$, the scalar field $`\varphi `$ should tend to a constant and $`I(t)1`$ in order not to violate the good agreement of BBN calculations. A more detailed numerical analysis of solutions of the BGT for $`z10^9`$ is needed to determine the behavior of $`I(t)`$. Such an analysis will be considered in a future publication.

If we now view the expansion of the universe from within the gravitational metric or VSGW frame, then the speed of light is constant but gravitational waves will move with a speed different from $`c_0`$. From (37), we obtain the contributions to (63):

$$\mathrm{\Omega }_{0m}=\frac{8\pi GK_0^{3/2}\rho _{0m}}{3H_0^2},\mathrm{\Omega }_{0k}=\frac{c_0^2kK_0}{R_0^2H_0^2},\mathrm{\Omega }_{0\mathrm{\Lambda }}=\frac{c_0^2K_0\mathrm{\Lambda }}{3H_0^2},$$

$$\mathrm{\Omega }_{0\varphi }=\frac{\stackrel{~}{\rho }_{0\varphi }}{6H_0^2}.$$
(67)
We obtain from (37) the deceleration parameter

$$q=\frac{dH^1}{dt}1=\frac{\dot{K}}{2HK}+\frac{1}{2}(1+\mathrm{\Omega }_k)+\frac{\stackrel{~}{p}_\varphi }{4H^2}+\frac{4\pi G}{c_0^2H^2}K^{1/2}p_m\frac{c_0^2K\mathrm{\Lambda }}{2H^2}.$$
(68)
By using the scalar wave equation (40), we find that :

$$\frac{\dot{K}}{2HK}=\frac{3(1K)(1+16\pi GBK^{1/2}p_m)+H^1B\dot{\varphi }KV^{}(\varphi )}{116\pi GBK^{3/2}\rho _m/c_0^2}.$$
(69)
The first term in the numerator and the denominator are positive-definite whenever the metrics have the correct signature, and the second term in the numerator is positive provided that $`H>0`$ and $`\dot{\varphi }KV^{}(\varphi )>0`$.

If we choose $`p_m0`$ and $`k=\mathrm{\Lambda }=0`$, we get a fit to the current data in the VSGW frame with

$$\mathrm{\Omega }_{0m}=0.28,\mathrm{\Omega }_{0\varphi }=0.72,q_0=0.63,$$
(70)
where

$$q_0=\frac{\dot{K}_0}{2H_0K_0}+\frac{1}{2}+\frac{\stackrel{~}{p}_{0\varphi }}{4H_0^2}.$$
(71)
We can achieve an acceleration of the universe, if the first term on the right-hand side of (71) dominates, which means that the cosmic acceleration with $`\ddot{R}>0`$ is caused by the dynamics of the scalar field $`\varphi `$ with $`\mathrm{\Lambda }=0`$ in the VSGW frame. Provided that $`q`$ turns small and positive at an earlier time in the universe’s expansion, then galaxy formation can be achieved.

Thus, in the VSL frame the universe appears to be decelerating and the supernovae are dimmed because of the increase of the speed of light with increasing red shift $`z`$ at least up to $`z34`$, while in the VSGW frame the universe appears to be accelerating and the supernovae appear to be farther away.

We have demonstrated elsewhere, that the early universe standard horizon and flatness problems can be resolved either in the VSL frame with a varying speed of light, or within the VSGW frame with constant speed of light, but with a varying speed of gravitational waves .

## 7 Fine-Structure Constant and Energy Conservation

Let us consider the behavior of the fine-structure constant $`\alpha =e^2/\mathrm{}c`$ and the energy of a system of particles in our BGT. We shall adopt the action for a charged particle moving in an electromagnetic field

$$S_{EM}=𝑑\widehat{s}\left[mc_0\left(\widehat{u}^\mu \widehat{u}^\nu \widehat{g}_{\mu \nu }\right)^{1/2}+\frac{e_0}{c_0}\widehat{u}^\mu \widehat{A}_\mu \right],$$
(72)
where

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu ,$$
(73)
is the electromagnetic field. Moreover, $`\widehat{F}^{\mu \nu }=\widehat{g}^{\mu \alpha }\widehat{g}^{\nu \beta }F_{\alpha \beta }`$ is the electromagnetic field observed in the matter VSL frame, $`\widehat{u}^\mu =\widehat{g}^{\mu \alpha }u_\alpha `$, $`\widehat{A}_\mu =\widehat{g}_{\mu \nu }A^\nu `$, $`m`$ is the particle mass and $`e_0`$ is the constant electron charge. We have in our isotropic and homogeneous cosmology $`A^\mu =(e_0/r,0)`$.

From (2) and (72), we obtain for the Coulomb energy

$$V_{\mathrm{Coul}}=\frac{e_0^2}{c_0r}\frac{dx_0}{ds},$$
(74)
where $`dx_0/ds=\gamma c_0`$ with $`\gamma =1/(1u^2/c_0^2)`$. We obtain in the nonrelativistic limit

$$V_{\mathrm{Coul}}=\frac{e_0^2}{r}.$$
(75)
In particular, the fine-structure constant

$$\alpha \frac{e^2}{\mathrm{}c}=\frac{e_0^2}{\mathrm{}c_0}$$
(76)
is constant in our bimetric gravity theory. This will be true whether or not we observe $`\alpha `$ in the VSL frame or in the VSGW frame.

This does not conform with the observations of absorption line spectra from quasars , which yield

$$\frac{\mathrm{\Delta }\alpha }{\alpha }=0.72\pm 0.18\times 10^5.$$
(77)
However, these data still require to be confirmed by an independent measurement procedure before we can rule out the constancy of $`\alpha `$ predicted by the BGT. If we wish to interpret the dimming of the supernovae by a $`10\%`$ increase in the speed of light in the VSL frame, then we cannot simultaneously demand a similar decrease in the fine-structure constant $`\alpha `$, because this would strongly disagree with the Webb et al. result (77) and lead to serious violations of the weak equivalence principle experimental tests. A recent analysis of the effects of a time varying $`\alpha `$ in the CMB background, gives results which are consistent with no variation in $`\alpha `$ from the epoch of recombination to the present day, and restricts any such variation to be less than about 4% . Forthcoming MAP and Planck experiments will be able to measure variations in $`\alpha `$ to better than a percent accuracy.

Finally, let us consider the conservation of energy in BGT. We have from (2):

$$\widehat{p}^\mu \widehat{p}_\mu =\widehat{g}^{\mu \alpha }p_\alpha \widehat{g}_{\mu \beta }p^\beta =p^\mu p_\mu =(mc_0)^2.$$
(78)
From the identification $`p^\mu =(E/c_0,\stackrel{}{p})`$, we obtain the standard special relativity result

$$E=c_0[(\stackrel{}{p})^2+(mc_0)^2]^{1/2}.$$
(79)
## 8 Cosmological Horizons

According to the models of quintessence , the dark energy of the universe is dominated by the potential $`V(\varphi )`$ of a scalar field $`\varphi `$, which rolls down to its minimum at $`V=0`$. We recall that for an equation of state $`p=w\rho `$ in GR, a cosmological constant corresponds to $`w=1`$, radiation domination to $`w=\frac{1}{3}`$ and matter domination to $`w=0`$. On the other hand, quintessence gives an equation of state with

$$1<w<\frac{1}{3},$$
(80)
while the observational evidence for a cosmological constant is given by the bound:

$$1<w_{\mathrm{observed}}\frac{2}{3}.$$
(81)
In the VSL frame, the proper horizon distance is given by

$$\delta _H(t)=R(t)F,$$
(82)
where

$$F=_{t_0}^{\mathrm{}}\frac{dt^{}c(t^{})}{R(t^{})}.$$
(83)
Whenever $`F`$ diverges there exist no future event horizons in the spacetime geometry. On the other hand, when $`F`$ converges the spacetime geometry exhibits a future horizon, and events whose coordinates at time $`\overline{t}`$ are located beyond $`\delta _H`$ can never communicate with the observer at $`r=0`$.

The variation of the expansion scale factor at large $`R(t)`$, when the curvature becomes negligible, approaches

$$R(t)t^{2/3(1+w)}.$$
(84)
We now have

$$F=c_0_{t_0}^{\mathrm{}}𝑑t^{}t^{[2(n1)/3(1+w)]},$$
(85)
where we have assumed that a solution of the BGT field equations leads to the behavior

$$c(t)c_0R^n,$$
(86)
where $`n`$ is some positive number. In the VSL frame with the metric $`\widehat{g}_{\mu \nu }`$, we can choose $`w>1/3`$, $`\mathrm{\Lambda }=0`$ and a value of $`n`$ such that the integral $`F`$ diverges and there is no future cosmological horizon. On the other hand, in the VSGW frame with the metric $`g_{\mu \nu }`$, the universe is accelerating with $`\ddot{R}>0`$, the velocity of light is constant, $`c(t)=c_0`$, and we expect that the integral $`F`$ will converge and allow a future cosmological horizon.

In an eternally accelerating universe, we are confronted with the difficulty of defining a consistent S-matrix description of quantum field theory and string/M-theories. We have shown that if we restrict ourselves to the VSL frame with a varying speed of light and a decelerating universe, then we can avoid asymptotically all future horizons associated with quintessence models and an accelerating universe. On the other hand, if we perform our experiments and theoretical calculations in the VSGW frame, then we can expect to encounter future cosmological horizons. Thus, it is advantageous for us to restrict ourselves to the VSL frame, since we can construct a viable S-matrix for quantum field theories and string/M-theories with a satisfactory asymptotic null infinity.

## 9 Conclusions

The problem of constructing a self-consistent gravity theory which is diffeomorphism invariant and allows a varying speed of light and a varying speed of gravitational waves can be solved within the BGT formalism. A proper conservation law for the matter tensor can be constructed that leads to particles moving on geodesics in the comoving $`\widehat{g}_{\mu \nu }`$ frame. In this frame, the speed of light is constant and observers can communicate with one another by light signals. An observer in the spacetime determined by the $`g_{\mu \nu }`$ metric, sees the speed of gravitational waves varying with time and the universe accelerate and a future cosmological horizon is expected to exist.

On the other hand, when the metric $`g_{\mu \nu }`$ is chosen to be comoving, so that the speed of gravitational waves is constant, and observers will be able (in the future) to communicate with one another by means of gravitational wave signals, then an observer in the $`\widehat{g}_{\mu \nu }`$ frame will detect a varying speed of light, a dimming of supernovae because of the increase of $`c(t)`$ in the past universe, and a decelerating universe without a cosmological horizon for an increasing speed of light in the future universe.

The observer in the VSL frame will be in an advantageous position when performing physical experiments, because the non-existence of a future cosmological horizon will permit the construction of a meaningful S-matrix for the scattering of particles.

An important feature of the BGT is that the dimensionless fine-structure constant $`\alpha `$ is constant in spacetime, so problems associated with possible violations of the weak equivalence principle tests and massive fine-tuning of vacuum quantum field theory can be avoided . However, future experiments will hopefully decide whether the putative observations of a variation of $`\alpha `$ by analyses of quasar absorption line spectra are valid .

There are many interesting consequences of BGT that can be explored e.g. the behavior of collapsing stars and the nature of dispersion relations at high energy. It is also of interest to investigate the possible role of BGT in quantum gravity. These topics will be the subject of future invetigations.

Acknowledgments

I thank Michael Clayton for helpful and stimulating discussions. This work was supported by the Natural Sciences and Engineering Research Council of Canada.


text:
# Inverse problem for the Landau-Zener effect

## Abstract

We consider the inverse Landau-Zener problem which consists in finding the energy-sweep functions $`W(t)\epsilon _1(t)\epsilon _2(t)`$ resulting in the required time dependences of the level populations for a two-level system crossing the resonance one or more times during the sweep. We find sweep functions of particular forms that let manipulate the system in a required way, including complete switching from the state 1 to the state 2 and preparing the system at the exact ground and excited states at resonance.

The Landau-Zener (LZ) effect , which consists in the quantum-mechanical transition between the two time-dependent levels of a two-level system (states 1 and 2) caused by crossing the resonance, is a well known phenomenon in many areas of physics. Well known applications of the LZ effect are those to molecular collision and dissociation. Recently, LZ effect has been used to measure the tunnel splitting in molecular magnets whereby detecting topological interference effects  which were predicted theoretically . If the system was in state 1 before crossing the resonance, then, for a constant sweep rate $`v`$, the probability to stay in this state after crossing the resonance is given by

$$P=\mathrm{exp}\left(\frac{\pi \mathrm{\Delta }^2}{2\mathrm{}v}\right),$$
(1)
where $`v|\dot{W}(t)|`$, $`W(t)\epsilon _1(t)\epsilon _2(t)`$, $`\epsilon _{1,2}(t)=\psi _{1,2}|\widehat{H}|\psi _{1,2}`$, $`\mathrm{\Delta }2\psi _1|\widehat{H}|\psi _2`$, $`W(t)`$ satisfies $`W(\pm \mathrm{})=\pm \mathrm{}`$ and $`\mathrm{\Delta }`$ is the tunnel level splitting. For fast energy sweep rates $`v`$ the system spends too little time in the vicinity of the resonance so that the tunneling matrix element $`\psi _1|\widehat{H}|\psi _2`$ cannot bring the system into the state 2 (say, onto the other side of a potential barrier), and the system remains in the initial state 1. For slow sweep rates the probability to remain in the initial state becomes exponentially small, and the system travels into the state 2 remaining on the lower exact energy term $`\epsilon _{}(t)`$ of

$$\epsilon _\pm (t)=\frac{1}{2}\left[\epsilon _2(t)+\epsilon _1(t)\pm \sqrt{W^2(t)+\mathrm{\Delta }^2}\right].$$
(2)
The corresponding adiabatic time dependence of the probability to be in state 1 is given by

$$|c_1(t)|^2=|\psi _1|\psi _{}(t)|^2=\frac{1}{2}\left(1\frac{W(t)}{\sqrt{W^2(t)+\mathrm{\Delta }^2}}\right).$$
(3)
In the case where $`\mathrm{}v\mathrm{\Delta }^2`$ and $`v`$ is constant, the solution of the Schrödinger equation for $`|c_1(t)|^2`$ can be found analytically  and it is oscillating at $`t>0`$ slowly approaching $`P`$ given by Eq. (1) at $`t\mathrm{}`$.

The known results above are valid for the linear energy sweep, $`W(t)=vt.`$ One can consider, however, other sweep functions $`W(t)`$ that behave nonlinearly in the vicinity of the resonance, $`W\mathrm{\Delta }`$. If the tunneling resonance is not too narrow and its position is well defined, nonlinear sweeps of a controlled form can be realized experimentally. Although no general analytical solution of the Schrödinger equation is known for the LZ problem with an arbitrary $`W(t),`$ one can look for the functions $`W(t)`$ that result in a desired time dependence of the probabillity $`|c_1(t)|^2`$, thus solving the inverse LZ problem. It is the main motivation of this Letter to study this inverse problem which, for instance, will allow to determine $`W(t)`$ needed to manipulate and to prepare a quantum system in a specific state.

One class of solutions of the time-dependent Schrödinger equation are smooth, non-oscillating functions of time that describe the complete switching of the system from the state 1 to the state 2, e.g., Eq. (3). Another class of solutions are those starting in $`\psi _1`$ at $`t=\mathrm{}`$ off resonance and ending in a given superposition of the exact states $`\psi _\pm `$ at resonance \[$`W(t)=0`$ for $`t0]`$. The most interesting variants of the latter are the system in the ground state $`\psi _{}`$ at resonance and in the excited state $`\psi _+`$ at resonance. These possibilities to manipulate the system in a required way could prove to be useful for the quantum computing (QC). Another interesting possibility is to sweep at a high rate (nonlinearly in time) far from the resonance thus making the whole process of crossing the resonance much shorter. Exploring these new types of solution of the LZ problem is the aim of this Letter.

To this end, we consider the Schrödinger equation for the coefficients in the wave function $`\psi (t)=c_1(t)\psi _1+c_2(t)\psi _2`$ that can be written in the form

$`\mathrm{}\dot{p}(t)`$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{\Delta }q(t),p(\mathrm{})=1`$
$`\mathrm{}\dot{q}(t)`$ $`=`$ $`iW(t)q(t){\displaystyle \frac{i}{2}}\mathrm{\Delta }p(t),q(\mathrm{})=0,W(\mathrm{})=\mathrm{}`$ (4)
with $`p(t)e^{i\epsilon _1t}c_1(t)`$ and $`q(t)e^{i\epsilon _1t}c_2(t)`$. We are looking for sweep functions $`W(t)`$ that result in a required time dependence of the state populations. Before solving this inverse LZ problem analytically, we will illustrate numerically the possibility of time-symmetric solutions of Eq. (Inverse problem for the Landau-Zener effect) going from $`p(\mathrm{})=1`$ to $`p(\mathrm{})=0`$ for strongly nonlinear sweep functions $`W(t)`$. The simplest of these cases is the cubic-parabola sweep with $`W(t)=2.5t+1.83t^3`$ where we have set $`\mathrm{\Delta }=\mathrm{}=1`$ in Eq. (Inverse problem for the Landau-Zener effect). The time dependence of the probability $`|p(t)|^2`$ shown in Fig. 1 can be fitted by $`|p|^2(t)=[1\mathrm{tanh}(\mathrm{sinh}(t))]/2`$ which is, however, not the analytical solution of Eq. (Inverse problem for the Landau-Zener effect).

To find analytical solutions of the inverse LZ problem, one should rewrite Eq. (Inverse problem for the Landau-Zener effect) in terms of the probability $`m|p|^2`$ and some phase variable. We thus define

$$p=\sqrt{m}\mathrm{exp}(i\phi _p),q=\sqrt{1m}\mathrm{exp}(i\phi _q)$$
(5)
and substitute it in Eq. (Inverse problem for the Landau-Zener effect). Evidently the global phase of the system $`(\phi _p+\phi _q)/2`$ is irrelevant and the resulting Schrödinger equation can be formulated in terms of $`m`$ and the phase difference $`\phi =\phi _q\phi _p`$:

$`\mathrm{}\dot{m}`$ $`=`$ $`\mathrm{\Delta }\sqrt{m(1m)}\mathrm{sin}\phi `$
$`\mathrm{}\dot{\phi }`$ $`=`$ $`W+{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \frac{12m}{\sqrt{m(1m)}}}\mathrm{cos}\phi .`$ (6)
This system of equations can be brought into a more elegant form by introducing

$$x2m1,1<x<1,$$
(7)
$$y\mathrm{arcsin}x,\pi /2<y<\pi /2,$$
(8)
and the reduced variables

$$\tau \left(\mathrm{\Delta }/\mathrm{}\right)t,\overline{W}W/\mathrm{\Delta }.$$
(9)
Then Eqs. (Inverse problem for the Landau-Zener effect) can be put into the reduced form

$`_\tau y`$ $`=`$ $`\mathrm{sin}\phi `$
$`_\tau \phi `$ $`=`$ $`\overline{W}\mathrm{tan}y\mathrm{cos}\phi .`$ (10)
The boundary condition for this equation is $`y(\mathrm{})=\pi /2`$ whereas the value of $`\phi (\mathrm{})`$ is irrelevant. In general, far from the resonance $`\phi (\tau )`$ strongly oscillates with time and the derivative $`_\tau y`$ is not small. Below we will consider special kinds of solutions of Eq. (10) which are characterized by a smooth, non-oscillating dependence $`\mathrm{sin}\phi (\tau )`$ and $`_\tau y(\mathrm{})=0`$. We will solve the inverse problem for Eq. (10) which consists in finding $`\overline{W}(\tau )`$ resulting in a required function $`x(\tau )`$. One easily obtains from Eq. (10) a general formula

$$\overline{W}(\tau )=_\tau \phi +\mathrm{tan}y\mathrm{cos}\phi =_\tau ^2y/\mathrm{cos}\phi +\mathrm{tan}y\mathrm{cos}\phi ,\mathrm{cos}\phi =\pm \sqrt{1(_\tau y)^2}.$$
(11)
The sign in front of $`\mathrm{cos}\phi `$ determines the sign of the function $`\overline{W}(\tau )`$. Since $`|c_1(\mathrm{})|^2=1`$ and thus $`x(\mathrm{})=1`$ and $`y(\mathrm{})=\pi /2`$ one has $`\mathrm{tan}y(\mathrm{})=\mathrm{}`$. To comply with the condition $`W(\mathrm{})=\mathrm{}`$ one has to choose, in general, the negative sign:

$$\overline{W}(\tau )=_\tau ^2y/\sqrt{1(_\tau y)^2}\mathrm{tan}y\sqrt{1(_\tau y)^2}.$$
(12)
Before considering particular cases, we will comment on some general properties of the functions $`x(\tau )`$ and $`\overline{W}(\tau )`$. If $`x(\tau )`$ is an odd function of time then $`y(\tau )`$ is also odd and $`_\tau y(\tau )`$ is even, whereas $`_\tau ^2y(\tau )`$ and $`\mathrm{tan}y(\tau )`$ are odd. Then from the first of Eqs. (10) follows that $`\mathrm{sin}\phi (\tau )`$ is even. If the sign of $`\mathrm{cos}\phi `$ does not change, which is the generic situation, $`\mathrm{cos}\phi (\tau )`$ is even. Then from Eq. (11) follows that $`\overline{W}(\tau )`$ is odd. In a similar way one can show that if $`x(\tau )`$ is even (the system returns into the initial state $`\psi _1`$ after crossing the resonance) then $`\overline{W}(\tau )`$ is even, too.

There is an important special case, however, which is realized if the derivative $`_\tau y`$ at some point $`\tau _0`$ attains its maximal value $`|_\tau y|=1`$ which follows from the first of Eqs. (10). At this point $`\mathrm{cos}\phi `$ turns to zero and changes its sign to positive. The latter corresponds to the phase $`\phi `$ going from the interval $`\pi /2<\phi <\pi `$ into the interval $`0<\phi <\pi /2`$. Note that $`_\tau ^2y=0`$ at $`\tau =\tau _0`$ thus Eq. (12) does not diverge. If $`\tau _0=0`$ then $`\mathrm{cos}\phi (\tau )`$ is an odd function: $`\mathrm{cos}\phi <0`$ for $`\tau <0`$ and $`\mathrm{cos}\phi >0`$ for $`\tau >0`$. As a result, the parity of $`\overline{W}(\tau )`$ becomes inverse to the parity of $`x(\tau )`$. If $`\tau _00`$ then for a symmetric function $`x(\tau )`$ the similar happens at $`\tau =\tau _0`$. That is, $`\mathrm{cos}\phi (\tau )`$ changes sign two times and is an even function, thus the parity of $`\overline{W}(\tau )`$ coincides with that of $`x(\tau )`$. It should be noted that parity considerations for the direct LZ problem do not apply: Symmetric forms of $`\overline{W}(\tau )`$ do not necessarily result in symmetric solutions $`x(\tau )`$. If, however, the solution $`x(\tau )`$ satisfies $`x(\pm \mathrm{})=1`$ or $`x(\pm \mathrm{})=\pm 1`$ and $`\overline{W}(\tau )`$ is time symmetric, then $`x(\tau )`$ is also symmetric. This follows from the time reversibility of the Schrödinger equation and the symmetry of the boundary values of $`x(\tau )`$.

Another implication of this analysis is the existence of the maximal rate of the probability change. We have seen that $`|_\tau y|1`$, then from $`_\tau x=\sqrt{1x^2}_\tau y`$ follows $`|_\tau x|1`$. This limitation is physically clear: Since tunneling between $`\psi _1`$ and $`\psi _2`$ is caused by the matrix element $`\mathrm{\Delta }2\psi _1|\widehat{H}|\psi _2`$, the latter sets up the natural scale for the tunneling rate, see Eq. (9). To completely tunnel from $`\psi _1`$ to $`\psi _2`$, the system should spend some minimal time in the vicinity of the resonance, that is, the sweep cannot be too fast. Below we will discuss two most important families of the functions $`x(\tau )`$ and corresponding $`\overline{W}(\tau )`$ which depend on $`\tau `$ in the combination $`u=a\tau `$, the parameter $`a`$ being the measure of the sweep rate. We will see that there is a maximal value $`a=a_{\mathrm{max}}`$ for which $`|_\tau y|=1`$ at some $`\tau =\tau _0`$. We will also see that although the functions $`x(\tau )`$ do formally exist for $`a>a_{\mathrm{max}}`$, one cannot extend $`\overline{W}(\tau )`$ into this region, as a result of the physical limitation on the tunneling rate. Thus solutions $`x(\tau )`$ with $`a>a_{\mathrm{max}}`$ are not realizable.

One of possible solutions is of the form

$$x=(1)^n\frac{u^n}{\sqrt{1+u^{2n}}},ua\tau .$$
(13)
For $`n`$ odd it describes switching of the probability from $`\psi _1`$ to $`\psi _2`$, whereas for $`n`$ even the system returns into the initial state. For $`n=1`$ the form of Eq. (13) corresponds to Eq. (3) by taking into account Eq. (7) with $`m=|c_1|^2`$. However, the reader should note that Eq. (3) is realized with the very slow sweep $`\overline{W}(\tau )`$ of an arbitrary functional form, whereas Eq. (13) is valid for a particular form of $`\overline{W}(\tau )`$ which is not necessarily slow. One obtains from Eqs. (8) and (13)

$`_\tau y`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1x^2}}}_\tau x=(1)^n{\displaystyle \frac{anu^{n1}}{1+u^{2n}}}`$
$`_\tau ^2y`$ $`=`$ $`(1)^n{\displaystyle \frac{a^2nu^{n2}\left[n1(n+1)u^{2n}\right]}{\left(1+u^{2n}\right)^2}},`$ (14)
whereas $`\mathrm{tan}y=x/\sqrt{1x^2}=(1)^nu^n`$. The derivative $`_\tau y`$ has a chance to attain the value 1 at its maximum that is determined by $`_\tau ^2y=0`$, i.e., at

$$u=u_{\mathrm{max}}=\pm \left(\frac{n1}{n+1}\right)^{1/(2n)}.$$
(15)
At this point one has

$$|(_\tau y)_{\mathrm{max}}|=\{\begin{array}{cc}a,\hfill & n=1\hfill \\ a\frac{n+1}{2}\left(\frac{n1}{n+1}\right)^{(n1)/(2n)},\hfill & n>1,\hfill \end{array}$$
(16)
thus the condition

$$aa_{\mathrm{max}}=\{\begin{array}{cc}1,\hfill & n=1\hfill \\ \frac{2}{n+1}\left(\frac{n1}{n+1}\right)^{2n/(n1)},\hfill & n>1\hfill \end{array}$$
(17)
should be fulfilled. For $`a<a_{\mathrm{max}}`$ Eq. (12) yields

$$\overline{W}(\tau )=\frac{(1)^{n1}u^n}{\sqrt{1\left[anu^{n1}/\left(1+u^{2n}\right)\right]^2}}\left[1\frac{a^2n\left[(2n+1)u^{2n2}(n1)u^2\right]}{\left(1+u^{2n}\right)^2}\right],u=a\tau .$$
(18)
This function is plotted for $`n=1`$ and $`n=5`$ and different values of the parameter $`a`$ in Figs. 2$`a`$ and 2$`b`$. For $`aa_{\mathrm{max}}`$ discontinuities in Eq. (18) are formed at $`u=\pm u_{\mathrm{max}}`$. These discontinuities result from the invariably negative choice for $`\mathrm{cos}\phi `$. For $`a=a_{\mathrm{max}}`$ the sign of $`\mathrm{cos}\phi `$ can change at $`u=u_{\mathrm{max}}`$ and then these discontinuities disappear. This form of $`\overline{W}(\tau )`$ includes the factor $`\mathrm{s}gn\left[\left(1+u^{2n}\right)^2\left(anu^{n1}\right)^2\right]`$, in addition to Eq. (18). Note that Eq. (18) cannot be extended into the region $`a>a_{\mathrm{max}}`$ preserving $`\overline{W}(\tau )`$ real, in contrast to Eq. (13). The small- and large-argument forms of $`\overline{W}(\tau )`$ are

$$\overline{W}(\tau )(1)^{n1}\{\begin{array}{cc}u^n,\hfill & |u|1\hfill \\ u^{n2},\hfill & |u|1,n>1\hfill \end{array}.$$
(19)
In the case $`n=1`$, Eq. (18) simplifies to

$$\overline{W}(\tau )=\frac{u}{1+u^2}\frac{13a^2+2u^2+u^4}{\sqrt{1a^2+2u^2+u^4}}.$$
(20)
One can see that for $`u1`$ the function $`\overline{W}(\tau )`$ describes a linear energy sweep with the velocity $`v=a`$. For $`u1`$ the form of $`\overline{W}(\tau )`$ is nonlinear and for $`a<a_{\mathrm{max}}=1`$ the resonance is crossed three times, including $`u=0`$. For $`a=1`$ Eq. (20) yields the even function

$$\overline{W}(\tau )=\frac{22\tau ^2\tau ^4}{(1+\tau ^2)\sqrt{2+\tau ^2}}$$
(21)
This sweep crosses the resonance two times, in the forward and backward directions, which means that the system is in the excited state $`\psi _2`$ at the end of the sweep at $`t=\mathrm{}`$. For $`a<1`$ and far from the resonance Eq. (20) has the form

$$W(t)=\mathrm{\Delta }\overline{W}(\tau )=\mathrm{\Delta }a\tau =\left(\mathrm{\Delta }^2a/\mathrm{}\right)t,$$
(22)
thus the corresponding linear sweep rate is $`v=\mathrm{\Delta }^2a/\mathrm{}`$. The maximal sweep rate that ensures the complete switching to $`\psi _2`$ is realized at $`a=1`$ is $`v_{\mathrm{max}}=\mathrm{\Delta }^2/\mathrm{}`$. For a comparison, in the regular LZ effect with the probability to stay in the initial state given by Eq. (1) the sweep rate $`v\pi \mathrm{\Delta }^2/(2\mathrm{})`$ is required for a complete switching. For $`v=\mathrm{\Delta }^2/\mathrm{}`$ the probability to stay in the initial state is still $`P=0.208`$.

Another solution that describes switching from state 1 to state 2 is

$$x=\mathrm{tanh}u,u=a\tau .$$
(23)
Here the sweep function determined from Eq. (12) reads

$$\overline{W}(\tau )=\frac{\mathrm{sinh}u}{\sqrt{1a^2/\mathrm{cosh}^2u}}\left[1\frac{2a^2}{\mathrm{cosh}^2u}\right].$$
(24)
In the range $`1/\sqrt{2}<a<1`$ the function $`\overline{W}(\tau )`$ changes its sign three times and its derivative at $`\tau =0`$ is negative. With $`a1`$ this derivative increases unlimitedly and a jump of $`\overline{W}(\tau )`$ at $`\tau =0`$ is formed. For $`a=a_{\mathrm{max}}=1`$ the solution with eliminated jump is an even function

$$\overline{W}(\tau )=\frac{2}{\mathrm{cosh}\tau }\mathrm{cosh}\tau .$$
(25)
This sweep function crosses the resonance two times. Again, Eq. (24) cannot be continued into the region $`a>1`$.

Now we return to the solution for the probability $`x`$ described by Eq. (13) and analyze its behavior for large $`n`$. It has a form of three plateaus $`x1`$ for $`u<1`$, $`x0`$ for $`1<u<\mathrm{\hspace{0.17em}1}`$, and $`x1`$ for $`1<u`$ with two narrow transition regions between them. In the central region $`1<u<\mathrm{\hspace{0.17em}1}`$ one has $`|c_2|^2|c_1|^2=(1+x)/21/2`$, i.e., the system is in the mixture of the states $`\psi _1`$ and $`\psi _2`$ with equal weights. The sweep function $`\overline{W}(\tau )`$ of Eq. (18) is fast approaching the resonance at $`u1`$ and remaining in the close vicinity of the latter until $`u1`$. Since the probabilities $`|c_1|^2|c_2|^21/2`$ do not oscillate in the range $`1<u<\mathrm{\hspace{0.17em}1}`$ (unlike the case in which the system is instantaneously brought into resonance), one is lead to the conclusion that the system is in one of the exact quantum states $`\psi _{}`$ or $`\psi _+`$. It is naturally of interest to question in which of these exact states does the system come under the influence of the energy sweep described by Eq. (18). We will see below that both variants can be realized. This opens a possibility to prepare the system in the ground or excited state in resonance by the appropriate choice of $`\overline{W}(\tau )`$. Indeed, for $`n3`$ one has $`\overline{W}(\tau )=_\tau y=_\tau ^2y=0`$ at $`\tau =0,`$ thus from Eq. (10) one can see that the sweep can be stopped at $`\tau =0`$ and the system will remain in this state. To demonstrate this behavior one should rewrite the wave function of the system $`\psi =c_1(t)\psi _1+c_2(t)\psi _2`$ determined above in the basis of the exact states

$$\psi _\pm (t)=\frac{1}{\sqrt{2}}\left[\pm K_\pm (t)\psi _1+K_{}(t)\psi _2\right],K_\pm (t)\sqrt{1\pm \frac{\overline{W}(t)}{\sqrt{1+\overline{W}^2(t)}}}$$
(26)
corresponding to the time-dependent “energy eigenvalues” of Eq. (2) as $`\psi =c_+(t)\psi _+(t)+c_{}(t)\psi _{}(t)`$. The corresponding probabillities are given by

$$|c_\pm (t)|^2=\frac{1}{2}\left\{1\pm \frac{\overline{W}(t)x(t)+\sqrt{1x^2(t)}\mathrm{cos}\phi (t)}{\sqrt{1+\overline{W}^2(t)}}\right\}.$$
(27)
One can see that in resonance at $`t=0`$ ($`\overline{W}=x=0`$, $`\mathrm{cos}\phi =\pm 1`$) the situation depends on the sign of $`\mathrm{cos}\phi `$. For sweep functions given by Eq. (18) with $`a<a_{\mathrm{max}}`$ one has $`\mathrm{cos}\phi =1`$ and the system is in the ground state ($`|c_{}(t)|^2=1`$, $`|c_+(t)|^2=0)`$. For the special sweep functions with $`a=a_{\mathrm{max}}`$ one has $`\mathrm{cos}\phi =1`$ and the system is in the excited state ($`|c_{}(t)|^2=0,`$ $`|c_+(t)|^2=1)`$, the transition from the ground state $`\psi _{}(t)`$ to the excited state $`\psi _+(t)`$ occuring for $`n1`$ in the vicinity of $`u=u_{\mathrm{max}}`$ of Eq. (15). These features are illustrated for the sweep functions of Eqs. (18) and (21) in Figs. 2$`c`$ and 2$`d`$. Note that for $`n=1`$ and $`a=a_{\mathrm{max}}=1`$ the sweep function of Eq. (21) returns to $`\mathrm{}`$ at $`\tau =\mathrm{}`$ and the system ends up in the excited state $`\psi _+=\psi _2`$.

To conclude, we have solved the inverse Landau-Zener problem for a two-level system, i.e., we have established the form of the sweep function $`W(t)\epsilon _1(t)\epsilon _2(t)`$ needed to insure a given probability $`|c_1(t)|^2`$ to be in the state 1 at time $`t`$. We have discussed two examples for which the system performs a complete transition from 1 to 2. In contrast to the direct LZ problem with a linear sweep, this can be accomplished by a finite sweeping rate. Further we have found such sweep functions that prepare the system in its exact ground state or in its exact excited state at resonance. It is not difficult to find $`W(t)`$ that prepare the system in an arbitrary mixture of the exact states at resonance, which can be potentially used in quantum computing based on tunneling units.