text:
# References

1. The aim of this Letter is to propose a quantitative explanation of the tiny anisotropy of the speed of light suggested by the recent, precise Michelson-Morley experiment of Müller et al. . Their result can be conveniently expressed in the form

$$B^{\mathrm{exp}}=(3.1\pm 1.6)10^9$$
(1)
where $`B`$ enters the Robertson-Mansouri-Sexl  parametrization for the two-way speed of light ($`c=2.9979\mathrm{}10^{10}`$ cm/sec)

$$\frac{\overline{c}}{c}=1(A+B\mathrm{sin}^2\theta )\frac{v^2}{c^2}$$
(2)
in a reference frame $`S^{}`$ that moves with speed $`v`$ with respect to a preferred frame $`\mathrm{\Sigma }`$, $`\theta `$ denoting the angle between the direction of $`v`$ and the direction of the light beam. 2. In order to explain the experimental result reported in Eq.(1), as a first step, I’ll adopt the tentative idea that light propagates in a medium with refractive index $`𝒩_{\mathrm{medium}}>1`$ so that there is a small Fresnel’s drag coefficient $`1\frac{1}{𝒩_{\mathrm{medium}}^2}1`$. This provides a general framework to analyze any Michelson-Morley type of experiment (see refs.). In our case, where the medium is the vacuum itself, the physical interpretation of $`𝒩_{\mathrm{vacuum}}`$ will represent a second step and provide a quantitative estimate to be compared with Eq.(1).

In this perspective, I’ll start introducing an isotropical speed of light

$$u=\frac{c}{𝒩_{\mathrm{medium}}}$$
(3)
that refers to the ideal case of a medium that extends to infinity in all spatial directions. Real experiments, however, are performed in finite portions of medium that might even reduce to just fill the arms of an interferometer. In this situation, an observer placed on the earth’s surface has no argument to think that light should propagate with the same velocity $`u`$ in all directions.

However, one may adopt the point of view that any observed anisotropy is due to the earth’s motion with respect to a preferred frame $`\mathrm{\Sigma }`$ where light propagates isotropically. In this case, if $`\mathrm{\Sigma }`$ were identified with the cosmic background radiation, one expects the relevant value of the earth’s velocity to be $`v_{\mathrm{earth}}365`$ km/sec . By adopting this point of view, and recalling that Lorentz transformations are valid both in Special and Lorentzian Relativity, for a frame $`S^{}`$, moving with respect to $`\mathrm{\Sigma }`$ with velocity $`𝐯`$, light will be seen to propagate at a speed ($`\gamma =1/\sqrt{1\frac{𝐯^2}{c^2}}`$)

$$𝐮^{}=\frac{𝐮\gamma 𝐯+𝐯(\gamma 1)\frac{𝐯𝐮}{v^2}}{\gamma (1\frac{𝐯𝐮}{c^2})}$$
(4)
To second order in $`v/u`$, one obtains ($`\theta `$ denotes the angle between $`𝐯`$ and $`𝐮`$)

$$\frac{u^{}(\theta )}{u}=1\alpha \frac{v}{u}\beta \frac{v^2}{u^2}$$
(5)
where

$$\alpha =(1\frac{1}{𝒩_{\mathrm{medium}}^2})\mathrm{cos}\theta +𝒪((𝒩_{\mathrm{medium}}^21)^2)$$
(6)
$$\beta =(1\frac{1}{𝒩_{\mathrm{medium}}^2})P_2(\mathrm{cos}\theta )+𝒪((𝒩_{\mathrm{medium}}^21)^2)$$
(7)
with $`P_2(\mathrm{cos}\theta )=\frac{1}{2}(3\mathrm{cos}^2\theta 1)`$.

Thus, the two-way speed of light is

$$\frac{\overline{u}^{}(\theta )}{u}=\frac{1}{u}\frac{2u^{}(\theta )u^{}(\pi \theta )}{u^{}(\theta )+u^{}(\pi \theta )}=1\frac{v^2}{c^2}(A+B\mathrm{sin}^2\theta )$$
(8)
where

$$A=𝒩_{\mathrm{medium}}^21+𝒪((𝒩_{\mathrm{medium}}^21)^2)$$
(9)
and

$$B=\frac{3}{2}(𝒩_{\mathrm{medium}}^21)+𝒪((𝒩_{\mathrm{medium}}^21)^2)$$
(10)
In this way, using the experimental values $`𝒩_{\mathrm{air}}1.00029`$ or $`𝒩_{\mathrm{helium}}1.000036`$, one can re-analyze  the classical ‘ether-drift’ experiments. For instance, by defining

$$v_{\mathrm{earth}}\sqrt{𝒩_{\mathrm{medium}}^21}=v_{\mathrm{obs}}$$
(11)
(and an in-air-operating optical system) one predicts $`v_{\mathrm{obs}}9`$ km/sec for $`v_{\mathrm{earth}}365`$ km/sec, in good agreement with Miller’s results .

3. To compare with Eq.(1) I’ll now try to provide a quantitative estimate of $`𝒩_{\mathrm{vacuum}}`$, to be used in Eq.(10), starting from the idea of a ‘condensed’ vacuum, as generally accepted in modern elementary particle physics. Indeed, in the physically relevant case of the Standard Model, the situation can be summarized saying  that ”What we experience as empty space is nothing but the configuration of the Higgs field that has the lowest possible energy. If we move from field jargon to particle jargon, this means that empty space is actually filled with Higgs particles. They have Bose condensed.”

In this case, where the condensing quanta are just neutral spinless particles (the ‘phions’ ), the translation from ‘field jargon to particle jargon’, amounts to establish a well defined functional relation (see ref.) $`n=n(\varphi ^2)`$ between the average particle density $`n`$ in the $`𝐤=0`$ mode and the average value of the scalar field $`\mathrm{\Phi }=\varphi `$. Thus, Bose condensation is just a consequence of the minimization condition of the effective potential $`V_{\mathrm{eff}}(\varphi )`$. This has absolute minima at some values $`\varphi =\pm v0`$ for which $`n(v^2)=\overline{n}0`$ .

The symmetric phase, where $`\varphi =0`$ and $`n=0`$, will eventually be re-established at a phase transition temperature $`T=T_c`$. This, in the Standard Model, is so high that one can safely approximate the ordinary vacuum as a zero-temperature system. Thus, the vacuum might be compared to a quantum Bose liquid, a medium where bodies can flow without any apparent friction, as in superfluid 4He, in agreement with the experimental results.

On the other hand, the condensed particle-physics vacuum, while certainly different from the ether of classical physics, is also different from the ‘empty’ space-time of Special Relativity which is assumed at the base of axiomatic quantum field theory. Therefore, following this line of thought, the macroscopic occupation of the same quantum state ($`𝐤=0`$ in a given reference frame) can represent the operative construction of a ‘quantum ether’ whose existence might be detected through a precise ‘ether-drift’ experiment.

On a more formal ground we observe that the coexistence of exact Lorentz covariance and vacuum condensation in effective quantum field theories is not so trivial. In fact, as a consequence of the violations of locality at the energy scale fixed by the ultraviolet cutoff $`\mathrm{\Lambda }`$ , one may be faced with non-Lorentz-covariant infrared effects that depend on the vacuum structure.

To indicate this type of infrared-ultraviolet connection, originating from vacuum condensation in effective quantum field theories, Volovik  has introduced a very appropriate name: reentrant violations of special relativity in the low-energy corner. In the simplest case of spontaneous symmetry breaking in a $`\lambda \mathrm{\Phi }^4`$ theory, where the condensing quanta are just neutral spinless particles, the ‘reentrant’ effects reduce to a small shell of three-momenta, say $`|𝐤|<\delta `$, where the energy spectrum deviates from a Lorentz-covariant form. Namely, by denoting $`M_H`$ as the typical energy scale associated with the Lorentz-covariant part of the energy spectrum, one finds $`\frac{\delta }{M_H}0`$ only when $`\frac{M_H}{\mathrm{\Lambda }}0`$.

The basic ingredient to detect such ‘reentrant’ effects in the broken phase consists in a purely quantum-field-theoretical result: the connected zero-four-momentum propagator $`G^1(k=0)`$ is a two-valued function . In fact, besides the well known solution $`G_a^1(k=0)=M_H^2`$, one also finds $`G_b^1(k=0)=0`$.

The b-type of solution corresponds to processes where assorbing (or emitting) a very small 3-momentum $`𝐤0`$ does not cost any finite energy. This situation is well known in a condensed medium, where a small 3-momentum can be coherently distributed among a large number of elementary constituents, and corresponds to the hydrodynamical regime of density fluctuations whose wavelengths $`2\pi /|𝐤|`$ are larger than $`r_{\mathrm{mfp}}`$, the mean free path for the elementary constituents.

This interpretation  of the gap-less branch, which is very natural on the base of general arguments, is unavoidable in a superfluid medium. In fact, ”Any quantum liquid consisting of particles with integral spin (such as the liquid isotope 4He) must certainly have a spectrum of this type…In a quantum Bose liquid, elementary excitations with small momenta $`𝐤`$ (wavelengths large compared with distances between atoms) correspond to ordinary hydrodynamic sound waves, i.e. are phonons. This means that the energy of such quasi-particles is a linear function of their momentum” . In this sense, a superfluid vacuum provides for $`𝐤0`$ a universal picture. This result does not depend on the details of the short-distance interaction and even on the nature of the elementary constituents. For instance, the same coarse-grained description is found in superfluid fermionic vacua  that, as compared to the Higgs vacuum, bear the same relation of superfluid 3He to superfluid 4He.

Thus there are two possible types of excitations with the same quantum numbers but different energies when the 3-momentum $`𝐤0`$: a single-particle massive one, with $`E_a(𝐤)M_H`$, and a collective gap-less one with $`E_b(𝐤)0`$. In this sense, the situation is very similar to superfluid 4He, where the observed energy spectrum is due to the peculiar transition from the ‘phonon branch’ to the ‘roton branch’ at a momentum scale $`|𝐤_o|`$ where $`E_{\mathrm{phonon}}(𝐤_o)E_{\mathrm{roton}}(𝐤_o)`$. The analog for the scalar condensate amounts to an energy spectrum with the following limiting behaviours :

i) $`E(𝐤)E_b(𝐤)c_s|𝐤|`$ for $`𝐤0`$

ii) $`E(𝐤)E_a(𝐤)M_H+\frac{𝐤^2}{2M_H}`$ for $`|𝐤|\delta `$

where the characteristic momentum scale $`\delta M_H`$, at which $`E_a(\delta )E_b(\delta )`$, marks the transition from collective to single-particle excitations. This occurs for

$$\delta 1/r_{\mathrm{mfp}}$$
(12)
where

$$r_{\mathrm{mfp}}\frac{1}{\overline{n}a^2}$$
(13)
is the phion mean free path, for a given value of the phion density $`n=\overline{n}`$ and a given value of the phion-phion scattering length $`a`$. In terms of the same quantities, one also finds

$$M_H^2\overline{n}a$$
(14)
giving the trend of the dimensionless ratios ($`\mathrm{\Lambda }1/a`$)

$$\frac{\delta }{M_H}\frac{M_H}{\mathrm{\Lambda }}\sqrt{\overline{n}a^3}0$$
(15)
in the continuum limit where $`a0`$ and the mass scale $`\overline{n}a`$ is held fixed.

By taking into account the above results, the physical decomposition of the scalar field in the broken phase can be conveniently expressed as (phys=‘physical’)

$$\mathrm{\Phi }_{\mathrm{phys}}(x)=v_R+h(x)+H(x)$$
(16)
with

$$h(x)=\underset{|𝐤|<\delta }{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{h}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{h}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(17)
and

$$H(x)=\underset{|𝐤|>\delta }{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{H}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{H}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(18)
where $`𝖵`$ is the quantization volume and $`E_k=c_s|𝐤|`$ for $`|𝐤|<\delta `$ while $`E_k=\sqrt{𝐤^2+M_H^2}`$ for $`|𝐤|>\delta `$. Also, $`c_s\delta M_H`$.

Eqs.(16)-(18) replace the more conventional relations

$$\mathrm{\Phi }_{\mathrm{phys}}(x)=v_R+H(x)$$
(19)
where

$$H(x)=\underset{𝐤}{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{H}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{H}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(20)
with $`E_k=\sqrt{𝐤^2+M_H^2}`$. Eqs.(19) and (20) are reobtained in the limit $`\frac{\delta }{M_H}\frac{M_H}{\mathrm{\Lambda }}0`$ where the wavelengths associated to $`h(x)`$ become infinitely large in units of the physical scale set by $`\xi _H=1/M_H`$. In this limit, where for any finite value of $`𝐤`$ the broken phase has only massive excitations, one recovers an exactly Lorentz-covariant theory. 4. In conclusion, for finite values of $`\mathrm{\Lambda }`$ there are long-wavelength density fluctuations of the vacuum and Lorentz-covariance is not exact. Therefore, in the presence of such effects, one can try to detect the existence of the scalar condensate through a precise ‘ether-drift’ experiment. To this end, I observe that a simple physical interpretation of the long-wavelength density fluctuation field

$$\phi (x)\frac{h(x)}{v_R}$$
(21)
has been proposed in refs.. Introducing $`G_F1/v_R^2`$ and choosing the momentum scale $`\delta `$ as

$$\delta =\sqrt{\frac{G_NM_H^2}{G_F}}$$
(22)
($`G_N`$ being the Newton constant) one obtains the identification

$$\phi (x)=U_N(x)+\mathrm{const}.$$
(23)
$`U_N(x)`$ being the Newton potential. Indeed, with the choice in Eq.(22), to first order in $`\phi `$ and in the limits of slow motions, the equations of motion for $`\phi `$ reduce to the Poisson equation for the Newton potential $`U_N`$  so that the deviations from Lorentz covariance are of gravitational strength. If, as in the Standard Model, $`G_F`$ is taken to be the Fermi constant one then finds $`\delta 10^5`$ eV and $`r_{\mathrm{mfp}}1/\delta =𝒪(1)`$ cm. As anticipated, the variation of $`\phi (x)`$ takes place over distances that are larger than $`r_{\mathrm{mfp}}`$ and thus infinitely large on the elementary particle scale. Also, by introducing $`M_{\mathrm{Planck}}=\frac{1}{\sqrt{G_N}}`$, and using Eqs.(15) and (22), one finds $`\mathrm{\Lambda }=q_HM_{\mathrm{Planck}}`$ with $`q_H=\sqrt{G_FM_H^2}=𝒪(1)`$, or $`a1/\mathrm{\Lambda }10^{33}`$ cm.

At the same time, to first order, the observable effects of $`\phi `$ can be re-absorbed  into an effective metric structure

$$ds^2=(1+2\phi )dt^2(12\phi )(dx^2+dy^2+dz^2)$$
(24)
that agrees with the first approximation to the line element of General Relativity . In this perspective, the space-time curvature arises from two sources: i) a re-scaling of the length and time units associated with the modification of any particle mass and ii) a refractive index for light propagation

$$𝒩_{\mathrm{vacuum}}12\phi $$
(25)
needed to preserve the basic particle-wave duality which is intrinsic in the nature of light.

Now, for a centrally symmetric field, and up to a constant, one has $`\phi (R)=\frac{G_NM}{c^2R}`$. Therefore, $`\phi _{\mathrm{earth}}0.710^9`$ (for $`M=M_{\mathrm{earth}}`$ and $`R=R_{\mathrm{earth}}`$) so that, using Eq.(10), I would estimate

$$B^{\mathrm{th}}4.210^9$$
(26)
in good agreement with the experimental result in Eq.(1). 5. Summarizing: the vacuum is not ‘empty’ so that one should check the consistency between exact Lorentz covariance and vacuum condensation in effective quantum field theories. For the specific case of the scalar condensate, the non-locality associated with the presence of the ultraviolet cutoff will also show up at long wavelengths in the form of non-Lorentz-covariant density fluctuations associated with a scalar function $`\phi (x)`$.

If, on the base of refs., these long-wavelength effects are naturally interpreted in terms of the Newton potential $`U_N`$ (with the identification $`\phi =U_N+\mathrm{const}.`$), one obtains the weak-field space-time curvature of General Relativity and a refractive index $`𝒩_{\mathrm{vacuum}}12\phi `$. This value of $`𝒩_{\mathrm{vacuum}}`$, leading to the prediction in Eq.(26), can help to understand the experimental result Eq.(1) obtained by Müller et al. .


format: text

text:
# Monte Carlo errors with less errors.

## 1 Formulation of the problem

In this article we collect some theory and formulas for the practical estimation of statistical and systematic errors in Monte Carlo simulations. Emphasis is put on the estimation of in general nonlinear functions (‘derived quantities’) of the primary expectation values. The strategy focuses on the explicit determination of the relevant autocorrelation functions and times. We shall discuss why this is advantageous compared to the popular binning techniques, which handle autocorrelations only implicitly. A Matlab code is made available, which implements the method described here.

The material given here is only partially new. It has been briefly discussed in the appendix of  and in internal notes of the ALPHA collaboration  and in lecture notes . However, due to the fact that mastering this topic is an indispensable prerequisite for every student or researcher in the popular field of lattice simulation, we felt that there is a gap in the readily citeable literature which we want to fill here.

We assume to have a number of primary observables enumerated by a Greek index with exact statistical mean values $`A_\alpha `$. The object to be estimated is a function of these,

$$Ff(A_1,A_2,\mathrm{})f(A_\alpha ).$$
(1)
A typical example is given by reweighting, where a quotient $`F=A_1/A_2`$ has to be determined. A more involved case are fit-parameters determined by correlation functions over a range of separations.

For the estimation we employ Monte Carlo estimates of the primary observables given by $`a_\alpha ^i`$. For each observable there are $`i=1,2,\mathrm{},N`$ successive estimates separated by executions of some valid update procedure. We assume that the Markov chain has been equilibrated before recording data beginning with $`i=1`$. This has to be secured from case to case, and we in particular recommend visual inspections of initial histories of observables to decide on safely discarding data. A posteriori one should also check that the equilibration time was much longer than autocorrelation times found in the analysis of the supposed equilibrium data.

A key theoretical quantity for error estimation is the autocorrelation  function $`\mathrm{\Gamma }_{\alpha \beta }`$ defined by

$$(a_\alpha ^iA_\alpha )(a_\beta ^{i+n}A_\beta )=\mathrm{\Gamma }_{\alpha \beta }(n).$$
(2)
that correlates the deviation of the $`i`$’th estimate for $`A_\alpha `$ with the deviation for variable $`\beta `$ after performing $`n0`$ updates. Averages $``$ here and throughout this paper refer to an ensemble of identical numerical experiments with independent random numbers and initial states — a concept useful for errors that are themselves statistical in nature. The dependence of $`\mathrm{\Gamma }_{\alpha \beta }`$ on the separation $`n`$ alone is due to being in equilibrium.

We take the simulation to produce configurations $`\varphi _i`$ distributed with the normalized Boltzmannian $`P(\varphi )`$. The update algorithm is specified by transition probabilities $`W(\varphi \varphi ^{})`$ and the measurements by $`a_\alpha ^i=𝒪_\alpha (\varphi _i)`$. Then the autocorrelation function in equilibrium is given as a double sum over configurations111 We here use a discrete notation for sums over configurations. The generalization to continuous fields with integrations is obvious. ,

$$\mathrm{\Gamma }_{\alpha \beta }(n)=\underset{\varphi }{}\underset{\varphi ^{}}{}P(\varphi )W^n(\varphi \varphi ^{})(𝒪_\alpha (\varphi )A_\alpha )(𝒪_\beta (\varphi ^{})A_\beta ).$$
(3)
Here $`W^n`$ characterizes transitions by a succession of $`n`$ update steps. If $`W`$ obeyed detailed balance with respect to $`P`$, we would conclude symmetry in the indices $`\alpha \beta `$. In practice we usually have only stability and no such symmetry seems to be implied in general. It is however natural and useful to define

$$\mathrm{\Gamma }_{\alpha \beta }(n)=\mathrm{\Gamma }_{\beta \alpha }(n)$$
(4)
implying a combined symmetry. Obviously $`\mathrm{\Gamma }_{\alpha \beta }(0)`$ is the symmetric covariance matrix with the ordinary variances on the diagonal. It is a (static) property of the simulated statistical system alone, while $`\mathrm{\Gamma }`$-values for nonzero separation are dynamic and depend on the Monte Carlo algorithm in use.

We end this theoretical setup by a slight generalization of our analysis. We add another index to the estimates of the primary observables, $`a_\alpha ^{i,r}`$. The index $`r=1,2,\mathrm{},R`$ counts a number of statistically independent replica. These are independent simulations which typically arise from repeated runs, or in particular, if parallel computers are used in the ‘trivial parallelization’ mode222The reader is invited to simplify all formulas to $`R=1`$ in a first reading.. For each replicum $`r`$ there are $`i=1,2,\mathrm{},N_r`$ successive equilibrium estimates separated by executions of the same update procedure. Obviously (2) becomes

$$(a_\alpha ^{i,r}A_\alpha )(a_\beta ^{j,s}A_\beta )=\delta _{rs}\mathrm{\Gamma }_{\alpha \beta }(ji).$$
(5)
## 2 Errors and biases

### 2.1 Primary observables

We define the per replicum means

$$\overline{a}_\alpha ^r=\frac{1}{N_r}\underset{i=1}{\overset{N_r}{}}a_\alpha ^{i,r}.$$
(6)
and in terms of them the natural estimator for the primary quantities $`A_\alpha `$

$$\overline{\overline{a}}_\alpha =\frac{1}{N}\underset{r=1}{\overset{R}{}}N_r\overline{a}_\alpha ^r$$
(7)
with the total number of estimates

$$N=\underset{r=1}{\overset{R}{}}N_r.$$
(8)
Introducing the deviations

$$\overline{\delta }_\alpha ^r=\overline{a}_\alpha ^rA_\alpha ,\overline{\overline{\delta }}_\alpha =\overline{\overline{a}}_\alpha A_\alpha .$$
(9)
we state that these estimators are unbiased,

$$\overline{\delta }_\alpha ^r=0=\overline{\overline{\delta }}_\alpha .$$
(10)
We assume a normal distribution333 Note that the $`a_\alpha ^{i,r}`$ themselves need not be Gaussian distributed and that large $`N_r`$ makes the $`\overline{a}_\alpha ^r`$ Gaussian by the central limit theorem. which is then completely defined by the covariance matrix given by

$$\overline{\delta }_\alpha ^r\overline{\delta }_\beta ^s=\frac{1}{N_r^2}\underset{i,j=1}{\overset{N_r}{}}\mathrm{\Gamma }_{\alpha \beta }(ji)\delta _{rs}=\frac{1}{N_r}C_{\alpha \beta }\delta _{rs}\times (1+\text{O}(\tau /N_r))$$
(11)
where

$$C_{\alpha \beta }=\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_{\alpha \beta }(t)$$
(12)
does not depend on the run-length. Here and in the following it is assumed that a finite scale $`\tau `$ characterizes the asymptotic exponential decay of $`\mathrm{\Gamma }_{\alpha \beta }`$ with $`|t|`$. The claimed $`N_r`$-dependence holds for runs with lengths $`N_r\tau `$ and if this condition is violated, an error estimation is hardly possible. For the covariance matrix of $`\overline{\overline{a}}_\alpha `$

$$\overline{\overline{\delta }}_\alpha \overline{\overline{\delta }}_\beta =\frac{1}{N}C_{\alpha \beta }\times (1+\text{O}(R\tau /N))$$
(13)
follows. Roughly speaking, as expected, $`\overline{\overline{a}}_\alpha `$ differs from $`A_\alpha `$ by an error of order $`1/\sqrt{N}`$ and the task of an error analysis amounts to a reasonably accurate estimation of $`C_{\alpha \beta }`$ which includes all autocorrelation effects.

### 2.2 Derived quantities

To determine $`F`$ we consider estimators

$$\overline{\overline{F}}=f(\overline{\overline{a}}_\alpha )$$
(14)
and

$$\overline{F}=\frac{1}{N}\underset{r=1}{\overset{R}{}}N_rf(\overline{a}_\alpha ^r).$$
(15)
We assume the estimates of the primary observables to be accurate enough to justify Taylor expansions of $`f`$ in the fluctuations, for example

$$\overline{\overline{F}}=F+\underset{\alpha }{}f_\alpha \overline{\overline{\delta }}_\alpha +\frac{1}{2}\underset{\alpha \beta }{}f_{\alpha \beta }\overline{\overline{\delta }}_\alpha \overline{\overline{\delta }}_\beta +\mathrm{}$$
(16)
with derivatives

$$f_\alpha =\frac{f}{A_\alpha },f_{\alpha \beta }=\frac{^2f}{A_\alpha A_\beta }$$
(17)
taken at the exact values $`A_1,A_2,\mathrm{}`$. It follows that our estimators for $`F`$ are biased unless $`f`$ is linear,

$$\overline{\overline{F}}F\frac{1}{2N}\underset{\alpha \beta }{}f_{\alpha \beta }C_{\alpha \beta }\frac{1}{R}\overline{F}F,$$
(18)
Usually, this bias is negligible compared to statistical errors due to large enough $`N`$. Replica ($`R2`$) can however be used to control and, if deemed to be appropriate, cancel the leading bias by using

$$\overline{\overline{F}}F\frac{1}{R1}\overline{F}\overline{\overline{F}}$$
(19)
to replace

$$\overline{\overline{F}}\frac{R\overline{\overline{F}}\overline{F}}{R1}.$$
(20)
The error $`\sigma _F`$ is to leading order given by

$$\sigma _F^2=(\overline{\overline{F}}F)^2\frac{1}{N}C_F$$
(21)
with

$$C_F=\underset{\alpha \beta }{}f_\alpha f_\beta C_{\alpha \beta }.$$
(22)
We rewrite this expression as

$$\sigma _F^2=\frac{2\tau _{\mathrm{int},F}}{N}v_F$$
(23)
with the effective ‘naive’ – i.e. disregarding autocorrelations – variance relevant for $`F`$

$$v_F=\underset{\alpha \beta }{}f_\alpha f_\beta \mathrm{\Gamma }_{\alpha \beta }(0)$$
(24)
and the integrated autocorrelation time for $`F`$

$$\tau _{\mathrm{int},F}=\frac{1}{2v_F}\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}\underset{\alpha \beta }{}f_\alpha f_\beta \mathrm{\Gamma }_{\alpha \beta }(t).$$
(25)
This one number encodes the efficiency of the algorithm in use for a determination of $`F`$, if its execution time per update is known. To appreciate this notation for the ratio of the true to the naive error two extremal cases are instructive. In the absence of autocorrelations, $`\mathrm{\Gamma }_{\alpha \beta }(t)\delta _{t,0}`$, we have $`2\tau _{\mathrm{int},F}=1`$ and (23) is the standard error. For a purely exponential behaviour, $`\mathrm{\Gamma }_{\alpha \beta }(t)\mathrm{exp}(|t|/\tau )`$, this scale is recovered, $`\tau _{\mathrm{int},F}=\tau +\text{O}(\tau ^1)`$. In general there are however many contributions and $`\tau _{\mathrm{int},F}`$ is more like a ‘dynamical susceptibility’. Another interpretation of (23) is, that only the reduced number of $`N/(2\tau _{\mathrm{int},F})`$ estimates are effectively independent to bring down statistical errors.

If we have $`R2`$ replica another consistency check is useful. The replica estimates $`f(\overline{a}_\alpha ^r)`$ are assumed to be normally distributed. Their variances are

$$\left(f(\overline{a}_\alpha ^r)F\right)^2=\frac{C_F}{N_r}.$$
(26)
We may imagine to do a $`\chi ^2`$ fit of the replica estimates to a constant $`K`$ by minimizing

$$\chi ^2(K)=\underset{r}{}\frac{(f(\overline{a}_\alpha ^r)K)^2}{C_F/N_r}.$$
(27)
The minimum occurs precisely at $`K=\overline{F}`$ given in (15). The probability to find in such a ‘fit’ a minimal $`\chi ^2x`$ is given by

$$Q=1P(x/2,(R1)/2),$$
(28)
where $`P`$ is the incomplete Gamma function,

$$P(U,n)=\frac{1}{\mathrm{\Gamma }(n)}_0^U𝑑uu^{n1}\mathrm{exp}(u).$$
(29)
Taking the actually observed $`x=\chi ^2(\overline{F})`$ we obtain the so-called goodness of fit. It is very unlikely to see $`Q`$-values much less than 0.1, if the underlying distribution has the required properties, and this is hence a good consistency check . In a histogram plot we may in addition inspect

$$p_r=\frac{f(\overline{a}_\alpha ^r)\overline{F}}{\sigma _F\sqrt{N/N_r1}}$$
(30)
which should have zero mean and unit variance.

In the next section we develop numerical estimators for $`v_F`$ and $`\tau _{\mathrm{int},F}`$, which at present are not yet of practical use, since they so far involve the numerically unknown exact $`A_\alpha `$ and $`\mathrm{\Gamma }_{\alpha \beta }`$.

## 3 Extraction of errors from measured data

### 3.1 Error estimators

We now describe what we would like to call the $`\mathrm{\Gamma }`$-method of error estimation (as opposed to binning, see App.B), where we explicitly estimate the autocorrelation function.

In terms of the given estimates $`a_\alpha ^{i,r}`$ we define an estimator $`\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }`$ of the autocorrelation function

$$\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(t)=\frac{1}{NRt}\underset{r=1}{\overset{R}{}}\underset{i=1}{\overset{N_rt}{}}(a_\alpha ^{i,r}\overline{\overline{a}}_\alpha )(a_\beta ^{i+t,r}\overline{\overline{a}}_\beta )$$
(31)
where $`0tN_r`$ for all $`r`$ is understood. An important principle underlying this choice is to avoid unnecessary noise from terms which vanish (only) on average. Therefore the correlation is first taken within each replicum and only then averaged. Using arguments of the same type as for (11) the leading bias of this estimator is evaluated,

$$\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(t)\mathrm{\Gamma }_{\alpha \beta }(t)\frac{C_{\alpha \beta }}{N}.$$
(32)
It comes from subtracting in (31) the ensemble mean instead of the exact average. For uncorrelated data this is the well known bias in the estimator of the variance cancelled by the $`N1`$ instead of $`N`$ normalization. For the moment we neglect this bias of order $`1/N`$ relative to the $`1/\sqrt{N}`$ statistical error in $`\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }`$ to be worked out in the next subsection, and come back to it later.

For the derived quantity $`F`$ only the projected autocorrelation is required,

$$\overline{\overline{\mathrm{\Gamma }}}_F(t)=\underset{\alpha \beta }{}\overline{\overline{f}}_\alpha \overline{\overline{f}}_\beta \overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(t)$$
(33)
where $`\overline{\overline{f}}_\alpha `$ is the gradient as in (17) but now evaluated at arguments $`\overline{\overline{a}}_1,\overline{\overline{a}}_2,\mathrm{}`$, which inflicts a small relative error of order $`1/N`$ that we do not attempt to cancel.

As practical estimators for $`v_F`$ and $`C_F`$ we take

$$\overline{\overline{v}}_F=\overline{\overline{\mathrm{\Gamma }}}_F(0)$$
(34)
and

$$\overline{\overline{C}}_F(W)=\left[\overline{\overline{\mathrm{\Gamma }}}_F(0)+2\underset{t=1}{\overset{W}{}}\overline{\overline{\mathrm{\Gamma }}}_F(t)\right].$$
(35)
A bias results from the truncation of the autocorrelation sum,

$$\frac{\overline{\overline{C}}_F(W)C_F}{C_F}\mathrm{exp}(W/\tau ).$$
(36)
For a purely exponential $`\mathrm{\Gamma }_F`$ this is an approximate equality, otherwise an order of magnitude statement with the scale $`\tau `$ introduced at the end of section 2.1. The summation window $`W`$ has to be chosen with care. On the one hand it has to be large compared to the decay time $`\tau `$ for a small systematic error, on the other hand with too large $`W`$ one includes terms with negligible signal but excessive noise. In the next subsection we shall define an automatic procedure to choose $`W`$.

In a frequently met case we have many primary observables $`A_\alpha `$ but only want to estimate a few functions of them. Then it is advantageous to actually not compute the full autocorrelation matrices $`\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(t)`$ but to just form $`\overline{\overline{a}}_\alpha `$ and $`\overline{\overline{f}}_\alpha `$ and immediately project the measured data

$$a_\alpha ^{i,r}a_f^{i,r}=\underset{\alpha }{}\overline{\overline{f}}_\alpha a_\alpha ^{i,r}$$
(37)
to base the remaining autocorrelation analysis directly on them.

Another practical point concerns the computation of $`\overline{\overline{f}}_\alpha `$. It may well be inconvenient or even impossible to analytically specify the gradient of $`f`$ and we may prefer to evaluate it numerically based on the supplied routine for $`f`$ itself only. A natural scale $`h_\alpha `$ is extracted from the data,

$$h_\alpha =\sqrt{\frac{\overline{\overline{\mathrm{\Gamma }}}_{\alpha \alpha }(0)}{N}},$$
(38)
and we take

$$\overline{\overline{f}}_\alpha \frac{1}{2h_\alpha }[f(\overline{\overline{a}}_1,\overline{\overline{a}}_2,\mathrm{},\overline{\overline{a}}_\alpha +h_\alpha ,\mathrm{})f(\overline{\overline{a}}_1,\overline{\overline{a}}_2,\mathrm{},\overline{\overline{a}}_\alpha h_\alpha ,\mathrm{})]$$
(39)
as a numerical estimate for the gradient of $`f`$ with negligible errors of $`\text{O}(h_\alpha ^2)1/N`$.

### 3.2 Error of the error

When an error is assessed numerically in a Monte Carlo it is itself affected by statistical error. This is usually not quantified. But to rate algorithms on the basis of the $`\tau _{\mathrm{int},F}`$ they produce for physically interesting $`F`$ it is better to at least have a rough idea. This is even more true, if one wants to falsify a theoretical prediction based on one of these ‘$`2.5\sigma `$ discrepancies’, while maybe the estimate of $`\sigma _F`$ itself is uncertain by a factor 2.

In appendix A the very convenient approximate formula

$$(\overline{\overline{C}}_F(W)C_F)^2\frac{2(2W+1)}{N}C_F^2.$$
(40)
is derived that was to our knowledge first given in . For the quotient

$$\overline{\overline{\tau }}_{\mathrm{int},F}(W)=\frac{\overline{\overline{C}}_F(W)}{2\overline{\overline{v}}_F}$$
(41)
we similarly find444 We here use error propagation and the upper bound in (58). This leads to an upper bound on the error of $`\overline{\overline{\tau }}_{\mathrm{int},F}`$ regarding the terms beyond the one $`W`$ which is both exact and dominant.  with the help of (57) and (58)

$$(\overline{\overline{\tau }}_{\mathrm{int},F}(W)\tau _{\mathrm{int},F})^2\frac{4}{N}(W+1/2\tau _{\mathrm{int},F})\tau _{\mathrm{int},F}^2.$$
(42)
Beside the statistical error of $`\overline{\overline{C}}_F(W)`$ there is the systematic error from the $`W`$-truncation in (36). We define as optimal the $`W`$ value that minimizes the sum of the absolute values of these errors and yields a total relative error

$$\frac{\delta _{\mathrm{tot}}(\overline{\overline{\sigma }}_F)}{\overline{\overline{\sigma }}_F}\frac{1}{2}\underset{W}{\mathrm{min}}\left(\mathrm{exp}(W/\tau )+2\sqrt{W/N}\right)$$
(43)
for the final error estimate

$$\overline{\overline{\sigma }}_F^2=\frac{\overline{\overline{C}}(W)}{N}$$
(44)
evaluated at the optimal $`W`$. It is a function of $`N/\tau `$ and plotted in Fig.1

for relevant arguments. The minimum leads to a simple transcendental equation that we solved numerically, but it can also be well approximated by $`W=\tau \mathrm{ln}(N/\tau )/2`$. The asymptotic behaviour is $`\delta _{\mathrm{tot}}\sqrt{\mathrm{ln}(N)/N}`$. Due to the exponential behaviour (36) the systematic component of the error becomes negligible at large N, albeit only slowly,

$$\frac{\delta _{\mathrm{sys}}(\overline{\overline{\sigma }}_F)}{\delta _{\mathrm{stat}}(\overline{\overline{\sigma }}_F)}=\frac{\tau }{2W}\frac{1}{\mathrm{ln}(N/\tau )}.$$
(45)
A simple and popular alternative to the analysis method just described are binning methods, where data are pre-averaged over sections of the history of length $`B`$. These bin averages are taken as independent estimates. For both ordinary and jackknife binning (see appendix B for more details) a relative systematic error of the error estimate occurs, which decays only proportional to $`\tau /B`$. The relative statistical error of the error due to the finite number of bins is $`\sqrt{2B/N}`$. By balancing these errors similarly as before we find the total error at optimal bin size

$$\frac{\delta _{\mathrm{tot}}(\overline{\sigma }_{F\mathrm{bin}})}{\overline{\sigma }_{F\mathrm{bin}}}\frac{1}{2}\underset{B}{\mathrm{min}}\left(\tau /B+\sqrt{2B/N}\right)=\frac{3}{2}(2N/\tau )^{1/3}$$
(46)
which is reached for

$$B=\tau (2N/\tau )^{1/3}$$
(47)
and is also shown in Fig.1. Here the ratio of systematic to statistical error is constant,

$$\frac{\delta _{\mathrm{sys}}(\overline{\sigma }_{F\mathrm{bin}})}{\delta _{\mathrm{stat}}(\overline{\sigma }_{F\mathrm{bin}})}=\frac{1}{2}.$$
(48)
In comparing the $`\mathrm{\Gamma }`$-method with binning it is clear now that the bin size $`B`$ and the window size $`2W`$ play a very similar role. All advantages of the $`\mathrm{\Gamma }`$-method derive from the exponential in $`W`$ rather than inverse power in $`B`$ behaviour of the systematic error of the error. One could say, that in this way we have an improved estimator of the error whose own error decays with a (slightly) larger power of $`N`$.

We briefly come back to the bias in our error estimate caused by the right hand side of (32). Cancelling this term when estimating $`C_F`$ by (35) amounts to enlarging

$$\overline{\overline{C}}_F(W)\overline{\overline{C}}_F(W)\left(1+\frac{2W+1}{N}\right),$$
(49)
which at the optimal $`W`$ value represents a correction down by a factor $`1/\sqrt{N}`$ (up to logs) relative to both contributions in (43). Since the correction is simple and known in closed form we include it in applications.

### 3.3 Automatic windowing procedure

We present now an algorithm to automatically choose $`W`$ that is similar but not identical to the one proposed in . To determine $`W`$ in an actual data analysis, we start from a hypothesis that $`\tau S\tau _{\mathrm{int},F}`$ with some factor $`S`$. More precisely, we solve

$$2\overline{\overline{\tau }}_{\mathrm{int},F}(W)=\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}(S|t|/\overline{\overline{\tau }}(W))$$
(50)
for $`\overline{\overline{\tau }}`$, to get a formula that is more precise for small $`\overline{\overline{\tau }}`$,

$$\frac{S}{\overline{\overline{\tau }}(W)}=\mathrm{ln}\left(\frac{2\overline{\overline{\tau }}_{\mathrm{int},F}(W)+1}{2\overline{\overline{\tau }}_{\mathrm{int},F}(W)1}\right)=\frac{1}{\overline{\overline{\tau }}_{\mathrm{int},F}}+\frac{1}{12\overline{\overline{\tau }}_{\mathrm{int},F}^3}+\mathrm{}.$$
(51)
If $`\overline{\overline{\tau }}_{\mathrm{int},F}1/2`$, we set $`\overline{\overline{\tau }}(W)`$ to a tiny positive value. Then we compute for $`W=1,2,\mathrm{}`$

$$g(W)=\mathrm{exp}[W/\overline{\overline{\tau }}(W)]\overline{\overline{\tau }}(W)/\sqrt{WN}.$$
(52)
Up to a factor this is the $`W`$-derivative defining the minimum in (43). The first value $`W`$ where $`g(W)`$ is negative (sign change) is taken as our summation window which under our hypothesis is self-consistently close to optimal. For most practical systems, algorithms and observables $`\tau `$ and $`\tau _{\mathrm{int},F}`$ are of the same order, and $`S=1\mathrm{}2`$ is a reasonable choice. A too large $`S`$ leads to (in principle) unnecessarily large statistical errors of the error.

It is important to verify by eye, at least for a representative set of observables, that $`\overline{\overline{\tau }}_{\mathrm{int},F}(W)`$ with its statistical errors (42) as a function of $`W`$ exhibits a plateau around the automatically chosen value. If this is not satisfactory, one has to modify $`S`$ to achieve this. A similar procedure could be set up to semi-objectively choose the bin size $`B`$ although this will be less stable, since the plateau behaviour will be less pronounced.

## 4 Conclusions

We have given a detailed description of the error analysis of Monte Carlo data. We found advantages in explicitly analyzing autocorrelation functions rather than employing binning strategies. They are due to a component in the systematic error of the error estimate from binning which is proportional to $`\tau /B`$. It is a finite bin-size effect in the integrated autocorrelation time from these procedures, which is avoided in the $`\mathrm{\Gamma }`$ method. The discussion includes a careful analysis of the errors of the error estimate, both statistical and systematic, and a minimization of their sum. Their approximate size can be read off from Fig.1 once one has an at least rough estimate of the autocorrelation decay scale $`\tau `$.

A Matlab routine that implements the whole method, including the necessary plots, is described and offered for download in the internet. It is also capable of handling multiple runs (or subdividing one run into several parts) and to judge the compatibility of the multiple estimates relative to the over-all error estimate. The routine adapts itself mostly automatically to the autocorrelation times involved. This works under reasonable assumptions mentioned in the text, but an additional visual inspection of the autocorrelations for at least part of the observables is always recommended.

Acknowledgements. I would like to thank the members of ALPHA for discussions and numerous tests of the software described here and its precursors. Philippe de Forcrand, Andreas Jüttner, Björn Leder, Francesco Knechtli, Thomasz Korzec, Martin Lüscher, Juri Rolf and Rainer Sommer helped with feedback on the manuscript.

## Appendix A The Madras-Sokal formula for the statistical error of the error

To prove (40) we need to investigate

$`\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(s)\overline{\overline{\mathrm{\Gamma }}}_{\gamma \delta }(t)=`$ (53)
$`{\displaystyle \frac{1}{Ns}}{\displaystyle \underset{i=1}{\overset{Ns}{}}}{\displaystyle \frac{1}{Nt}}{\displaystyle \underset{j=1}{\overset{Nt}{}}}(a_\alpha ^i\overline{a}_\alpha )(a_\beta ^{i+s}\overline{a}_\beta )(a_\gamma ^j\overline{a}_\gamma )(a_\delta ^{j+t}\overline{a}_\delta ).`$
This formula is for only one replicum at the moment to avoid even more indices. The essential approximation consists now of approximating this four-point function by its disconnected part,

$`\overline{\overline{\mathrm{\Gamma }}}_{\alpha \beta }(s)\overline{\overline{\mathrm{\Gamma }}}_{\gamma \delta }(t)(\mathrm{\Gamma }_{\alpha \beta }(s)C_{\alpha \beta }/N)(\mathrm{\Gamma }_{\gamma \delta }(t)C_{\gamma \delta }/N)+`$ (54)
$`+{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\left[\mathrm{\Gamma }_{\alpha \gamma }(ji)\mathrm{\Gamma }_{\beta \delta }(ji+ts)+\mathrm{\Gamma }_{\alpha \delta }(ji+t)\mathrm{\Gamma }_{\beta \gamma }(jis)\right]`$
This is plausible for large $`N`$. At least gathering from experience with ordinary statistical systems the connected part makes a smaller contribution since all arguments have to be close to each other. We have also assumed $`|s|,|t|N`$ to extend the sums with relative errors of order $`|s|/N,|t|/N`$. Furthermore we use for a rapidly decaying function $`g`$ the approximation

$$\frac{1}{N^2}\underset{i,j=1}{\overset{N}{}}g(ji)\frac{1}{N}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}g(k)$$

similarly as in (12) to obtain

$$(\overline{\overline{\mathrm{\Gamma }}}_F(s)\mathrm{\Gamma }_F(s))(\overline{\overline{\mathrm{\Gamma }}}_F(t)\mathrm{\Gamma }_F(t))=\frac{1}{N}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_F(k)[\mathrm{\Gamma }_F(k+ts)+\mathrm{\Gamma }_F(kts)].$$
(55)
Summing now $`t,s=W,\mathrm{},W`$ and assuming that each factor $`\mathrm{\Gamma }_F`$ cuts off the sum unless its argument is much smaller in absolute value than $`W`$ we derive

$$(\overline{\overline{C}}_F(W)C_F)^2\frac{2}{N}(2W+1)\underset{k,l=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_F(k)\mathrm{\Gamma }_F(l)=\frac{2}{N}(2W+1)C_F^2$$
(56)
with corrections of $`\text{O}(\tau /W)`$.

By very similar steps one may also derive

$$(\overline{\overline{v}}_Fv_F)(\overline{\overline{C}}_F(W)C_F)\frac{2}{N}C_F^2.$$
(57)
A more complicated case is

$$(\overline{\overline{v}}_Fv_F)^2\frac{2}{N}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}(\mathrm{\Gamma }_F(k))^2\frac{2}{N}v_FC_F,$$
(58)
where the inequality certainly holds if $`0\mathrm{\Gamma }(k)\mathrm{\Gamma }(0)`$ is true. Even if this happens to be violated out in the tail of $`\mathrm{\Gamma }`$, we still expect the inequality to hold for the sum (dominated by small $`|k|`$) in practical applications.

Repeating the above analysis with several replica leads to the same formulas and thus yields (40). In view of the approximations we note that the practical experience is, that the resulting errors of $`\tau _{\mathrm{int},F}`$ are consistent with their scatter when several runs are made.

## Appendix B Binning methods

Here we give the basic formulas for error estimation by binning. We assume data $`a_\alpha ^i`$ where possible replica are sewed together to one history of length $`N=BN_B`$ which we divide into $`N_B`$ sections of $`B`$ consecutive measurements each555 The breaks of autocorrelation at replica boundaries are neglected and a possible remainder is discarded.. We form bins or blocked measurements

$$b_\alpha ^k=\frac{1}{B}\underset{i=1}{\overset{B}{}}a_\alpha ^{(k1)B+i},k=1,\mathrm{},N_B$$
(59)
The mean values are all identical $`\overline{b}=\overline{a}=\overline{\overline{a}}`$. The idea is now to choose $`B`$ large enough such that the bins $`b_\alpha ^k`$ can be regarded as approximately uncorrelated, but small enough such that there remains a large number $`N_B`$ of such ‘events’. As an error estimator one then takes

$$\overline{\sigma }_{F\mathrm{bin}}^2=\frac{1}{N_B(N_B1)}\underset{k=1}{\overset{N_B}{}}(f(b_\alpha ^k)\overline{F})^2.$$
(60)
If we Taylor-expand $`f`$ again, the relevant correlation matrix is

$`(b_\alpha ^k\overline{a}_\alpha )(b_\beta ^k\overline{a}_\beta )={\displaystyle \frac{1}{B^2}}{\displaystyle \underset{u,v}{}}\mathrm{\Gamma }_{\alpha \beta }(uv)`$ (61)
$`{\displaystyle \frac{1}{BN}}{\displaystyle \underset{u,i}{}}(\mathrm{\Gamma }_{\alpha \beta }(ui)+\mathrm{\Gamma }_{\alpha \beta }(iu))+{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j}{}}\mathrm{\Gamma }_{\alpha \beta }(ij)`$
where indices $`u,v`$ run through bin $`k`$ only while $`i,j`$ range from 1 to $`N`$ as before. The result is dominated by the first term

$$\frac{1}{B^2}\underset{u,v}{}\mathrm{\Gamma }_{\alpha \beta }(uv)\frac{1}{B}C_{\alpha \beta }(1\tau /B)$$
(62)
while the other terms approximately sum to $`C_{\alpha \beta }/N`$. The coefficient of $`\tau /B`$ is exact if $`\mathrm{\Gamma }_{\alpha \beta }`$ is simply exponential, but in practice this should be read as O($`\tau /B`$). This term is the analog of (36) for the $`\mathrm{\Gamma }`$-method. We see that the integrated autocorrelation time is effectively constructed by a double sum within one bin of size $`B`$ and the $`\tau /B`$ error is a finite size effect which masks the actual $`\mathrm{exp}(B/\tau )`$ expected from truncation. The many bins are only used to tame statistical fluctuations of the estimate but do not improve the systematic error. The analogous finite size effect for the $`\mathrm{\Gamma }`$-method occurs only at order $`\tau /N`$. A cancellation of the $`\tau /B`$ or $`\mathrm{exp}(W/\tau )`$ effects dominant for the two methods was not found to be practical, as one would have to use a noisy estimate for $`\tau `$ obtained from the data and would not know the precise pre-factor.

In (60) the function $`f`$ has to be evaluated on averages over only $`B`$ events, that is only the $`N_B`$’th fraction of all measurements. This may be affected with stability problems in the case of fits for example. Compared to our earlier discussion, the gradient of $`f`$ is here implicitly formed stochastically with larger average spreads $`\sqrt{N_B}h_\alpha `$. Both problems are cleverly overcome by going to complementary or jackknife bins,

$$c_\alpha ^k=\frac{1}{NB}\left(\underset{i=1}{\overset{N}{}}a_\alpha ^iBb_\alpha ^k\right).$$
(63)
An elementary calculation yields the relation

$$(c_\alpha ^k\overline{a}_\alpha )(c_\beta ^k\overline{a}_\beta )=\frac{1}{(N_B1)^2}(b_\alpha ^k\overline{a}_\alpha )(b_\beta ^k\overline{a}_\beta ).$$
(64)
Hence the resulting jackknife error estimator is

$$\overline{\sigma }_{F\mathrm{jack}}^2=\frac{N_B1}{N_B}\underset{k=1}{\overset{N_B}{}}(f(c_\alpha ^k)\overline{F})^2.$$
(65)
Here all evaluations of $`f`$ practically work on the full statistics, and the spread of its arguments is of the same order $`h_\alpha `$ as taken for the numerical derivative in the $`\mathrm{\Gamma }`$-method. The result of the above analysis of systematic errors is however unchanged for the jackknife estimate.

For completeness we finally mention that also under binning the scatter between $`\overline{F}`$, the function $`f`$ applied to the mean over all data, and the average of $`f(c_\alpha ^k)`$ over all bins can be used to cancel the leading $`\text{O}(1/N)`$ bias in the estimate of the mean of derived observables. In analogy to (20) we have to substitute

$$\overline{F}N_B\overline{F}+(1N_B)\frac{1}{N_B}\underset{k=1}{\overset{N_B}{}}f(c_\alpha ^k).$$
(66)
## Appendix C An implementation of the $`\mathrm{\Gamma }`$-method

In the ALPHA collaboration we have recently performed most data analyses in Matlab, since we found that it combines comfortable programming, a robust library, interactive graphics and acceptable speed suitable for this purpose. An improved version of the well tested core-routine is described in the first subsection. While it has been applied to many realistic data sets in the ALPHA projects, we here develop a simulator, that allows for clinical testing with exactly known errors in a situation that we consider representative for simulations in lattice QCD.

### C.1 Description of the Matlab routine UWerr

The calling sequence of the function is

```
     [value,dvalue,ddvalue,tauint,dtauint,Q] = ...
     UWerr(Data,Stau,Nrep,Name,Quantity,P1,P2,...)

```

with the meaning of the input and output arguments described below. Similar information is obtained from help UWerr under Matlab.
Input:
Data is a matrix filled with the primary estimates $`a_\alpha ^{i,r}`$ from $`R`$ replica with $`N_1,N_2,\mathrm{},N_R`$ measurements each,

$$\mathrm{𝙳𝚊𝚝𝚊}=\left(\begin{array}{cccc}a_1^{1,1}\hfill & a_2^{1,1}\hfill & a_3^{1,1}\hfill & \mathrm{}\\ a_1^{2,1}\hfill & a_2^{2,1}\hfill & a_3^{2,1}\hfill & \mathrm{}\\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\\ a_1^{N_1,1}\hfill & a_2^{N_1,1}\hfill & a_3^{N_1,1}\hfill & \mathrm{}\\ a_1^{1,2}\hfill & a_2^{1,2}\hfill & a_3^{1,2}\hfill & \mathrm{}\\ \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\hfill & \mathrm{}\\ a_1^{N_R,R}\hfill & a_2^{N_R,R}\hfill & a_3^{N_R,R}\hfill & \mathrm{}\end{array}\right)$$

All input arguments that follow can be omitted or set to empty (=\[\]) and are then assigned default values.

If data after reading them into Matlab occur in a different matrix format, e.g. one long vector, the Matlab commands reshape may be employed to transform them very efficiently and safely.
Stau is the estimate $`S=\tau /\tau _{\mathrm{int}}`$ explained in section (3.3). The default is $`S=1.5`$. It is tried to determine an optimal window $`W`$ as described, but it is never taken larger than $`\mathrm{min}\{N_r\}/2`$. If the windowing condition cannot be met, a warning is printed.
Nrep $`=[N_1,N_2,\mathrm{},N_R]`$ is a vector whose elements specify the lengths of the individual replica. As default all data are interpreted as one replicum, Nrep = \[N\]. Even if only one run has been made, one may specify Nrep such as to analyze sections of the history independently and investigate the scatter of sub-averages.
Name can be a string (default: ’NoName’) which is used as a title in the generated plots (see below). If the argument is given but is not a string (an integer for example) no plots are generated.
Quantity can either be an integer from the range of the index $`\alpha `$, in which case the corresponding primary quantity is analyzed. Alternatively it may be a function handle666 At present also a string with the name can be entered, which however is obsolete and announced not to be supported in future versions of Matlab.  @fun  if fun  is the name of the function $`f(A_1,A_2,\mathrm{})`$ which has to be user-defined and has as first argument a row-vector which has as many entries as Data has columns. In this case the corresponding derived quantity is analyzed. Optional further arguments $`P1,P2,\mathrm{}`$ are passed on to the function as second, third, …argument. The default value is Quantity=1.
Output:
value,dvalue,ddvalue,tauint,dtauint,Q have their obvious meaning as best estimates of $`F,\sigma _F`$ as defined in (1),(21), statistical error of the error from (40), $`\tau _{\mathrm{int},F}`$ and $`\delta (\tau _{\mathrm{int},F})`$ given in (25) and extracted from (42). The $`Q`$-value (28) follows if $`R2`$. Note that $`C_F,v_F`$ and the chosen $`W`$ can be reconstructed from the output if necessary.

The leading bias is corrected according to (20) for $`R2`$. If this correction amounts to more than $`\sigma _F/4`$ the user is warned. The routine generates — unless suppressed — plots of our estimate for $`\rho (t)=\mathrm{\Gamma }_F(t)/\mathrm{\Gamma }_F(0)`$ and $`\overline{\overline{\tau }}_{\mathrm{int},F}(W)`$ with errorbars in the relevant range to inspect the required plateau behaviour. Errors of $`\rho (t)`$ have been added in version 6 of the UWerr. They are computed with formula (E.10) in  taking $`\mathrm{\Lambda }=t+W`$. If there are two or more replica their distribution (30) is shown as a histogram with its $`Q`$-value (28) in the title. In addition a histogram of all $`a_\alpha ^{i,r}`$ is displayed in case that a primary observable is analyzed.

The code (m-file) of the Matlab routine UWerr.m described here together with a sample function effmass.m used in the following test can be down-loaded from www-com.physik.hu-berlin.de/ALPHAsoft

### C.2 Test in a simulator

As a test case we think of determining the ‘effective’ mass

$$m=\mathrm{ln}(G(z)/G(z+1))=f(G(z),G(z+1))F$$
(67)
from a correlation $`G(z)=\mathrm{exp}(mz)`$ in euclidean time $`z`$, which we take simply exponential since we are not interested here in systematic errors of $`m`$. To fake noisy (auto)correlated data for $`G`$ we employ normally distributed random numbers $`\eta _i,i=1,2,\mathrm{}`$ with zero mean and unit variance to produce a sequence $`\{\nu _i\}`$ by recursively setting777 I thank Hubert Simma for pointing out an improvement here.

$$\nu _1=\eta _1,\nu _{i+1}=\sqrt{1a^2}\eta _{i+1}+a\nu _i,$$
(68)
with

$$a=\frac{2\tau 1}{2\tau +1},\tau \frac{1}{2}.$$
(69)
A short calculation yields

$$\nu _i=0,\nu _{i+t}\nu _i=a^t,\underset{t=\mathrm{}}{\overset{\mathrm{}}{}}a^{|t|}=2\tau .$$
(70)
We use three independent such sequences $`\nu _i^{(1)},\nu _i^{(2)},\nu _i^{(3)}`$ and choose for them integrated autocorrelation times $`\tau ^{(1)}<\tau ^{(2)}=\tau ^{(3)}`$.

Synthetic data to measure $`m=f(A_1,A_2)`$ are formed,

$$a_1^i=G(0)+q(\nu _i^{(1)}+\nu _i^{(2)}),a_2^i=G(1)+q(\nu _i^{(1)}+\nu _i^{(3)}).$$
(71)
The real parameter $`q`$ sets the strength of the noise, and the presence of $`\nu _i^{(1)}`$ in both observables creates a correlation between them, as in real life.

All relevant quantities of the error analysis can be worked out exactly for this case,

$$v_{A_1}=v_{A_2}=2q^2,\tau _{\mathrm{int},A_1}=\tau _{\mathrm{int},A_2}=(\tau ^{(1)}+\tau ^{(2)})/2,$$
(72)
$$v_F=2q^2(1+\mathrm{exp}(2m)\mathrm{exp}(m)),$$
(73)
$$\tau _{\mathrm{int},F}=\frac{\widehat{m}^2/2}{\widehat{m}^2+1}\tau ^{(1)}+\frac{\widehat{m}^2/2+1}{\widehat{m}^2+1}\tau ^{(2)}$$
(74)
with

$$\widehat{m}=2\mathrm{sinh}(m/2).$$
(75)
Note that in the weighed sum giving $`\tau _{\mathrm{int},F}`$ the time scale $`\tau ^{(1)}`$ contributes with a small weight, due to correlations between $`a_i^1,a_i^2`$ and the cancellation of the corresponding fluctuations in the effective mass.

For our test we now set $`m=0.2,\tau ^{(1)}=4,\tau ^{(2)}=8,q=0.2`$ which leads to $`v_F=0.1016,\tau _{\mathrm{int},F}=7.92`$. We generate $`8\times 1000`$ estimates corresponding to an error in the mass of 0.0142 and get

```
>> Nr=ones(8,1)*(N/8);
>> [value,dvalue,ddvalue,tauint,dtauint,Q]= ...
   UWerr(Data,1,Nr,’mass’,@effmass,1,2);
>> [value,dvalue,ddvalue,tauint,dtauint,Q]
ans =

    0.2128    0.0134    0.0008    7.2909    0.7396   0.2827

```

The generated plots are shown in Fig.2 and in Fig.3. By repeating such an analysis 20000 times we knock down the statistical error of the error far enough to exhibit a systematic error of $`0.5\%`$ which is of the expected order $`\frac{1}{2}\mathrm{exp}(W/\tau )`$ with $`W=37`$ and $`\tau =8`$, while, of course, the actual automatic $`W`$ and estimated $`\tau `$-values fluctuate.

The jackknife analog of the lower part of Fig.2 can be seen in Fig.4 for the same data. It looks consistent with a less pronounced plateau, as expected.

More relevant are tests on real data. The routine passed a number of such tests, but, as usual with software, with probability one some future application will force further evolution.


text:
# New effects in neutrino oscillations in matter and electromagnetic fieldsA short version of invited talk given at the 4th International Bruno Pontecorvo School on ”Neutrino Oscillations, CP and CPT Violation: the Three Windows to Physics Beyond the Standard Model”, Capri, May 26-29, 2003.

Department of Theoretical Physics, Moscow State University, 119992, Moscow, Russia

At present neutrino oscillation phenomenon is the most evidently confirmed window to physics beyond the Standard Model. The thorough studies of neutrino oscillations were launched by the pioneering paper of Bruno Pontecorvo . Indeed, it is surprising that after more than 45 years of researches in this field it is still possible to discover new effects in neutrino oscillations.

It is a pleasure for me to present at the 4th International School dedicated to Bruno Pontecorvo the four new effects in neutrino oscillations that have been recently studied in my research group at the Moscow State University. The whole story began several years ago when we studied neutrino spin oscillations in strong magnetic fields  and to our much surprise realized that at that time in literature there were no even attempts to consider neutrino spin oscillations in any electromagnetic field configurations rather than constant in time and transversal in respect to neutrino motion magnetic fields (see, for example, ). Further more, in all of the studies of neutrino spin and also flavour oscillations in matter, performed before 1995, the matter effect was treated only in the non-relativistic limit.

The first our attempt  to consider neutrino flavour oscillations in matter in the case when matter is moving with relativistic speed was made in 1995. In that our study we tried to apply the Lorentz invariant formalism to describe neutrino flavour oscillations and realized that the value of the matter term in the neutrino effective potential can be significantly changed if matter is moving with relativistic speed. We have continued  our studies on evaluation of the Lorentz invariant formalism in neutrino spin oscillations in 1999 and have developed an approach that enables us to consider neutrino spin oscillations in arbitrary electromagnetic field configurations. We have shown  how the Bargmann-Michel-Telegdi equation , describing a neutral particle spin evolution under the influence of an electromagnetic field, can be generalized for the case of neutrino moving in electromagnetic field and matter. The corresponding Lorentz invariant neutrino spin evolution equation have been derived from the Bargmann-Michel-Telegdi equation by the following substitution: one has to add to the electromagnetic field tensor $`F_{\mu \nu }`$ an anti-symmetric tensor $`G_{\mu \nu }`$ that can be constarcted with use of the neutrino speed, matter speed, and matter polarization four-vectors under the natural assumptions that the equation have to be linear over $`F_{\mu \nu }`$ and over the other mentioned above vectors.

From this new generalized BMT equation for the neutrino spin evolution in an electromagnetic field and matter the corresponding effective Hamiltonian describing neutrino spin oscillations is just straightforward . Thus, the first new effect is the prediction of neutrino spin oscillations in various electromagnetic field configurations. We have derived new resonances in neutrino oscillations for several electromagnetic fields such as the field of an electromagnetic wave and superposition of an electromagnetic wave and constant magnetic field.

We have predicted , using the generalized BMT equation, the second new effect: neutrino spin procession can be stimulated not only by presence of electromagnetic fields (i.e., by electromagnetic interactions of neutrino) but also by weak interactions of neutrino with background matter. Moreover, we have shown  that the neutrino spin precession always occurs in presence of matter (even in the case of non-moving and unpolarized matter) if the initial neutrino state is not longitudinally polarized. We have also considered  the neutrino spin evolution problem within the Lorentz invariant approach in the case of general types of neutrino non-derivative interactions with external fields.

Then we have used the derived new equation for neutrino spin evolution for the case when matter may move with arbitrary (also relativistic) speed and may be polarized. We have predicted  the third effect: the matter motion can drastically change the neutrino oscillation pattern and, in particular, can sufficiently shift the neutrino spin oscillation resonance condition , if compared with the case of non-moving matter. The analogous effect exist  also in neutrino flavour oscillations: the neutrino resonance condition  can be significantly modified if matter is moving with relativistic speed.

The fourth new effect in neutrino oscillations is the prediction  for the new mechanism of electromagnetic radiation by neutrino moving in background matter and/or electromagnetic fields. This new mechanism of electromagnetic radiation originates from the neutrino spin precession that can be produced whether by weak interactions with matter or by electromagnetic interactions with external electromagnetic fields. We have named this radiation as ”spin light of neutrino” ($`SL\nu `$). The total power of the $`SL\nu `$ does not washed out when the emitted photon refractive index is equal to unit and the $`SL\nu `$ can not be considered as the neutrino Cerenkov radiation (see, for example,  and references therein). The discovered important properties of $`SL\nu `$ (such as strong beaming of the radiation along the neutrino momentum, the rapid growth of the total radiation power with the neutrino energy and density of matter, the possibility to emit photons with energies span up to gamma-rays) enables us to predict that this radiation should be important in different astrophysical environments (quasars, gamma-ray bursts etc) and in dense plasma of the early Universe.


text:
# With the relativistic velocity addition law through special relativity

## I Introduction

Each chapter of physics operates with fundamental physical quantities. Kinematics is involved with length and time, dynamics add to them, by tradition, the concept of mass, whereas electromagnetism introduces the concept of electric charge. Combining fundamental physical quantities we obtain new ones, useful in describing physical effects.

All chapters of physics have in common the relativistic postulate, due to Galileo and Einstein, according to which all the true laws of physics are the same for all inertial observers, all of them measuring the value $`c`$ for the velocity of light. When measured from different reference frames some physical quantities have the same magnitude (relativistic invariants). If a physical quantity which characterizes an object is measured by an observer relative to whom the object is in a state of rest, we say that the observer measures its proper value.

An equation which relates physical quantities measured from different reference frames, represents a transformation equation. A transformation equation enable us to derive formulas which account for different relativistic effects. It is considered that the use of the transformation equations reduces the transparency 1,2.

Relativistic equations describing time dilation1,3,4,5, length contraction and addition of relativistic velocities1,6,7 can be derived without using transformation equations, all having as a starting point the time dilation formula. The method works when in one of the involved reference frames a proper physical quantity is measured.

The purpose of our paper is to show that the transformation equations for kinematic, dynamic and electromagnetic physical quantities can be derived using the same strategy, involving the addition law of relativistic velocities as a single relativistic ingredient, which at its turn is a direct consequence of the time dilation effect7.

Three inertial reference frames are involved $`S(XOY)`$, $`S^{}(X^{}O^{}Y^{})`$ and $`S^o(X^oO^0Y^0)`$. $`S^o`$ moves with velocity $`u_x`$ relative to $`S`$, with velocity $`u_x^{}`$ relative to $`S^{}`$, $`S^{}`$ moving at its turn with velocity $`v`$ relative to $`S`$. The axes of the reference frames are parallel to each other, the $`OX`$, $`O^{}X^{}`$ and $`O^oX^o`$ axes are overlapped and all the velocities show in theirs positive directions. The three reference frames have the same space-time origins. In order to keep the problem as simple as possible, we will consider only events which take place on the overlapped axes. The velocities mentioned above are related by

$$u_x=\frac{u_x^{}+v}{1+{\displaystyle \frac{u_x^{}v}{c^2}}}$$
(1)
which solved for $`u_x^{}`$ leads to

$$u_x^{}=\frac{u_xv}{1{\displaystyle \frac{u_xv}{c^2}}}.$$
(2)
Reference frame $`S^o(X^oO^0Y^0)`$ is a special one. In that frame nothing and nobody moves and only events which take place at it origin $`O^o`$ are taken into account. The result is that there, only physical quantities which characterize a stationary body can be detected, like the reading of clock $`C_o^o`$ located at $`O^o(t^o)`$, the mass $`m^o`$ of a particle located at $`O^o`$, the force $`F^o`$ acting on the particle mentioned above and the electric field created by a charged capacitor. Physical quantities like different from zero space coordinates, momentum, mechanical power and magnetic field are meaningless in that frame. All the physical quantities measured in $`S^o`$ are by definition proper physical quantities: proper time reading (date) $`t^o`$, proper mass $`m^o`$, proper force $`F^o`$ and proper electric field $`E^o`$.

## II Relativistic kinematics

In each of the mentioned reference frames, the clocks located at the different points of the overlapped axes are synchronized following a procedure proposed by Einstein. The proposed scenario makes that at a given point of the overlapped axes clock $`C_o^o`$ located at $`O^o`$ meets a clock $`C(x,y=0)`$ of the $`S`$ frame and a clock $`C^{}(x^{},y^{}=0)`$ of the $`S^{}`$ frame. Let $`t^o`$, $`t`$ and $`t^{}`$ be the readings of the three clocks respectively. In accordance with the time dilation formula the three readings are related by

$$t=\frac{t^o}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(3)
and by

$$t^{}=\frac{t^o}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}.$$
(4)
Eliminating $`t^o`$ we obtain that $`t`$ and $`t^{}`$ are related by

$$t=t^{}\frac{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}=t^{}\frac{1+{\displaystyle \frac{u_x^{}v}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}=\frac{t^{}+{\displaystyle \frac{v(u_x^{}t^{})}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}$$
(5)
and by

$$t^{}=t\frac{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}=t\frac{1{\displaystyle \frac{u_xv}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}=\frac{t{\displaystyle \frac{v(u_xt)}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}.$$
(6)
We detect in the right hand side of equations (5) and (6) the presence of the terms ($`u_xt`$) and ($`u_x^{}t^{}`$) respectively. We have to find out names for them. Taking into account that clock $`C_o^o`$ is located at the common origin of time at the common origins of the reference frames, $`x=u_xt`$ can represent the space coordinate of clock $`C_o^o`$ in the $`S`$ frame when the synchronized clocks of that frame read $`t`$ or the distance travelled by $`C_o^o`$ during the time interval ($`t0`$), measured of course in $`S`$. For the same reasons, $`x^{}=u_x^{}t`$ can represent the distance travelled by $`C_o^o`$ in the time interval ($`t^{}0`$), measured in $`S^{}`$ or the space coordinate of clock $`C_o^o`$ in the $`S^{}`$ frame when the synchronized clocks read there $`t^{}`$. With this new notations equations (5) and (6) become

$$t=\gamma _v\left(t^{}+\frac{vx^{}}{c^2}\right);\gamma _v=\frac{1}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}$$
(7)
$$t^{}=\gamma _v\left(t\frac{vx}{c^2}\right).$$
(8)
The space coordinates $`x`$ and $`x^{}`$ transform as

$$x=u_xt=\gamma _vt^{}(u_x^{}+v)=\gamma _v(x^{}+vt^{})$$
(9)
and

$$x^{}=u_x^{}t^{}=\gamma _vt(u_xv)=\gamma _v(xvt).$$
(10)
We have derived so far the Lorenz transformations for the events $`E(x=u_xt,y=0,t)`$ associated with the fact that clock $`C_o^o`$ meets clock $`C(x=u_xt,y=0)`$ and $`E^{}(x^{}=u_x^{}t^{},y^{}=0,t^{})`$ associated with the fact that clock $`C_o^o`$ meets clock $`C^{}(x^{}=u_x^{}t^{},y^{}=0)`$. Because the two events take place at the same point in space, relativists say that they represent the same event.

## III Relativistic dynamics

Consider that a particle is located at rest at the origin $`O^o`$ of $`S^o`$. Observers in that frame measure its rest mass $`m^o`$. Using the addition law of relativistic velocities it was shown13 that as measured from $`S`$ the mass of the particle is

$$m=\frac{m_0}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(11)
whereas as measured from $`S^{}`$ it is

$$m^{}=\frac{m_0}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}},$$
(12)
resulting that $`m`$ and $`m^{}`$ are related by

$$m=m^{}\frac{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}=m^{}\frac{1+{\displaystyle \frac{vu_x^{}}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}=\gamma _v\left(m^{}+\frac{v(m^{}u_x^{})}{c^2}\right)$$
(13)
and by

$$m^{}=m\frac{1{\displaystyle \frac{vu_x}{c^2}}}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}}=\gamma _v\left(m\frac{v(mu_x)}{c^2}\right).$$
(14)
It is considered that $`m`$ and $`m^{}`$ represent the relativistic mass of the same particle as measured from $`S`$ and $`S^{}`$ respectively. We detect in the right hand side of equations (13) and (14) the presence of the terms $`mu_x`$ and $`m^{}u_x^{}`$ which have the physical dimensions of momentum. Physicists baptize physical quantities by a single name and label them by a single notation. So they introduce the notations

$$p_x=mu_x=\frac{m^ou_x}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(15)
$$p_x^{}=m^{}u_x^{}=\frac{m^ou_x^{}}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}},$$
(16)
$`p_x`$ and $`p_x^{}`$ representing the relativistic momentum of the particle with rest mass $`m^o`$ as measured from $`S`$ and $`S^{}`$ respectively.

The relativistic momentum transforms as

$$p_x=mu_x=\gamma _vm^{}(u_x^{}+v)=\gamma _v(p_x^{}+vm^{})$$
(17)
$$p_x^{}=m^{}u_x^{}=\gamma _vm(u_xv)=\gamma _v(p_xvm).$$
(18)
The physical quantity $`E=mc^2`$ ($`E^{}=m^{}c^2`$) has the physical dimensions of energy. Due to the invariance of $`c`$, it transform as mass does

$$E=E^{}\gamma _v\left(1+\frac{vu_x^{}}{c^2}\right)=\gamma _v(E^{}+vp_x^{})$$
(19)
and

$$E^{}=E\gamma _v\left(1\frac{vu_x}{c^2}\right)=\gamma _v(Evp_x).$$
(20)
Expressed as a function of energy, equations (19) and (20) become

$$p_x=\gamma _v\left(p_x^{}+\frac{vE^{}}{c^2}\right)$$
(21)
and

$$p_x^{}=\gamma _v\left(p_x\frac{vE}{c^2}\right).$$
(22)
Equations (19) and (22) put the bane on the concept of relativistic mass as many physicists require8.

A remarkable property of the relativistic energy consists in the fact that a particle at rest is characterized by a rest energy

$$E^o=m^oc^2.$$
(23)
Because the single supplementary energy a free particle moving with constant velocity can possess is its kinetic energy $`E_k`$, we have for it

$$E_k=EE_o=E_o(\gamma _v1)=m^oc^2(\gamma _v1).$$
(24)
It is clear that $`p_x`$ ($`p_x^{}`$), $`E`$ ($`E^{}`$) express the momentum and the energy of the same particle as measured from $`S`$ and $`S^{}`$.

The traditional way to derive the equations presented above is to consider a collision between two particles from two inertial reference frames and imposing that momentum and energy conserve9.

The truth of equation (11) is confirmed by the experiments performed Bucherer10 with electrons moving in a magnetic field, obtaining excellent agreement working with it. This was and continues to be a great proof in favour of the special theory of relativity in all its fields, all being the consequence of Einstein’s postulate of relativity. It is considered10 that to this date no other explanation of Bucherer’s experiment then that supplied by special relativity, is known.

## IV Relativistic electrodynamics

We consider now a capacitor $`C_o`$ the plates of which are parallel to the $`X^oO^oZ^o`$ plane, at rest in $`S^o`$ and charged with electrical charges $`Q`$ and $`Q`$ respectively. It creates an electric field $`E_y^o`$, oriented in the positive direction of the $`O^oY^o`$ axis given by

$$E_y^o=\frac{Q}{\epsilon _0\epsilon _rA^o},$$
(25)
$`A^o`$ representing the surface of one of the plates $`\epsilon _0`$ and $`\epsilon _r`$ representing invariant electrical constants. The electrical neutrality of atoms convinces us that the electric charge is a relativistic invariant. When detected from $`S`$, the sides of the plate parallel with the overlapped axes undergo a Lorentz contraction and so the electric field measured from $`S`$ is

$$E_y=\frac{Q}{\epsilon _0\epsilon _rA^o\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(26)
whereas in $`S^{}`$ it is

$$E_y^{}=\frac{Q}{\epsilon _0\epsilon _rA^o\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}.$$
(27)
The same procedure, as in the previous cases leads to the following relationship between $`E_y`$ and $`E_y^{}`$

$$E_y=E_y^{}\gamma _v\left(1+\frac{vu_x^{}}{c^2}\right)=\gamma _v\left(E_y^{}+\frac{v(u_x^{}E_y^{})}{c^2}\right)$$
(28)
$$E_y^{}=E_y\gamma _v\left(1\frac{vu_x}{c^2}\right)=\gamma _v\left(E_y\frac{v(u_xE_y)}{c^2}\right).$$
(29)
We detect in the right hand side of equations (28) and (29) the presence of the terms ($`u_xE_y`$) and ($`u_x^{}E_y^{}`$) and label them $`B_z`$ and $`B_z^{}`$ respectively, representing the $`OZ`$ ($`OZ`$) components of a vector $`\stackrel{}{B}(0,0,B)`$ in $`S`$ and $`\stackrel{}{B}^{}(0,0,B^{})`$ in $`S^{}`$, baptized magnetic field. Its direction is chosen in such a way that the vector products $`\stackrel{}{u}_x\times \stackrel{}{E}_y`$ and $`\stackrel{}{u}_x^{}\times \stackrel{}{E}_y^{}`$ show in the same direction as $`E_y`$ and $`E_y^{}`$ do, respectively. It transform as

$$B_z=\gamma _vE_y^{}(u_x^{}+v)=\gamma _v(B_z^{}+vE_y^{})$$
(30)
$$B_z^{}=\gamma _vE_y(u_xv)=\gamma _v(B_zvE_y).$$
(31)
Expressed as a function of the magnetic field, equations (28) and (29) become

$$E_y=\gamma _v\left(E_y^{}+\frac{vB_z^{}}{c^2}\right)$$
(32)
$$E_y^{}=\gamma _v\left(E_y\frac{vB_z}{c^2}\right).$$
(33)
With the plates of the capacitor parallel to the $`Y^oO^oZ^o`$, it creates an electric field parallel to the overlapped axes. In that case the sides of the plates are not subjected to the Lorentz contraction and so

$$E_x^0=E_x=E_x^{}.$$
(34)
The transformation equations for the components $`E_z`$ ($`E_z^{}`$) and $`B_y`$ ($`B_y^{}`$) are obtained starting with the capacitor the plates of which are parallel to the $`XOY`$ plane producing an electric field orientated along the $`OZ`$ axis.

The traditional way to derive the equations presented above is to use the transformation properties of force11,12.

## V Conclusions

We have shown that if a proper physical quantity $`F^o`$ (proper time, proper mass, proper electric field) can be defined in a given reference frame $`S^o`$ which moves with velocity $`u_x`$ relative to $`S`$ and with velocity $`u_x^{}`$ relative to $`S^{}`$ and transforms as

$$F=\frac{F^o}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(35)
$$F^{}=\frac{F^o}{\sqrt{1{\displaystyle \frac{u_x^2}{c^2}}}}$$
(36)
respectively, then a transformation equation relating $`F`$ and $`F^{}`$ can be derived using the addition law of velocities as a single relativistic ingredient. The transformation generate a new physical quantity ($`uF`$), ($`u^{}F^{}`$) in $`S`$ and in $`S^{}`$ respectively for which physicists find adequate names and the transformation of which involves the transformation of $`F`$ and $`F^{}`$ and the addition law of velocities.

As compared with the traditional derivations which are very laborious, our derivation is heuristic, being based on facts known from classical physics, uses simple algebra being time saving having a good relationship between time invested and understanding achieved.

We consider that it convinces the reluctants about the truth of special relativity. It does not charge the memory of the teacher or of the learner, who in the traditional approaches, have to learn in each chapter of physics, a new scenario.