∎

11institutetext: SM. Nasiri 22institutetext: 22email<EMAIL_ADDRESS>
_School of ECE, College of Engineering, University of Tehran, Tehran, Iran_
33institutetext: R. Hosseini 44institutetext: Tel.: +98-21-82089799
44email<EMAIL_ADDRESS>
_School of ECE, College of Engineering, University of Tehran, Tehran, Iran_
_School of Computer Science, Institute of Research in Fundamental Sciences
(IPM), Tehran, Iran_
55institutetext: H. Moradi 66institutetext: 66email<EMAIL_ADDRESS>
_School of ECE, College of Engineering, University of Tehran, Tehran, Iran_
_Intelligent Systems Research Institute, SKKU, South Korea_

# Optimal Triangulation Method is Not Really Optimal

Seyed-Mahdi Nasiri  Reshad Hosseini  Hadi Moradi

(Received: date / Accepted: date)

###### Abstract

Triangulation refers to the problem of finding a 3D point from its 2D
projections on multiple camera images. For solving this problem, it is the
common practice to use so-called optimal triangulation method, which we call
the $L_{2}$ method in this paper. But, the method can be optimal only if we
assume no uncertainty in the camera parameters. Through extensive comparison
on synthetic and real data, we observed that the $L_{2}$ method is actually
not the best choice when there is uncertainty in the camera parameters.
Interestingly, it can be observed that the simple mid-point method outperforms
other methods. Apart from its high performance, the mid-point method has a
simple closed formed solution for multiple camera images while the $L_{2}$
method is hard to be used for more than two camera images. Therefore, in
contrast to the common practice, we argue that the simple mid-point method
should be used in structure-from-motion applications where there is
uncertainty in camera parameters.

###### Keywords:

Triangulation Structure-from-Motion Mid-point method

##  1 Introduction

Figure 1: Triangulation methods: a) Mid-point method finds a 3D point which
has the minimum sum of squared distances from the lines of sights. b) $L_{2}$
triangulation method finds a 3D point which has the minimum sum of squared 2D
distances between its projections and the corresponding points on images.

One of the fundamental tasks in 3D vision is reconstructing a point in 3D from
its projections on the camera images in two or multiple views. This task,
which is called triangulation, is used extensively in machine vision and
robotics applications such as stereo vision tippetts2016review , mapping
castle2011wide  and structure-from-motion bartoli2005structure ;
ramalingam2006generic . It is obvious that a 3D point can be simply found by
intersecting the lines of sight of each projection when no noise is present.
But in practice, due to several sources of noise such as uncertainties in
relative camera poses, errors in cameras intrinsic parameters, and subpixel
inaccuracies in the position of matched points, all the lines do not
necessarily intersect at one point or at all. There are many attempts to solve
the triangulation problem in the presence of uncertainties
hartley1997triangulation ; lindstrom2010triangulation ; lee2019closed .

The common method for solving triangulation in two views, called the optimal
method in the literature, is casting it as a nonlinear optimization problem.
In such a problem, a new point is found as close as possible to the measured
point in each view so that the lines of sights for the new points intersect.
In hartley1997triangulation , the authors compared several triangulation
methods on several simulated datasets. The authors observed that the optimal
triangulation method, which we call the $L_{2}$ method in this paper,
outperforms other methods. Apart from good performance, the authors argued
that the $L_{2}$ triangulation method has a nice property to be projective
invariant. In the performed simulations, the authors considered the
uncertainty in the position of corresponding points and not camera parameters.

In the majority of triangulation problems, there are uncertainties in both
cameras’ parameters and measured corresponding points. Thus, a question arises
here, “Is the $L_{2}$ triangulation method still the best performing method?”
To the best of our knowledge, this question has not been investigated in the
literature. In this paper, we evaluate the performance of triangulation
methods in the calibrated structure-from-motion setting. The calibrated
structure-from-motion after the work of nister2004efficient  has become the
natural choice for structure-from-motion applications. In this case, we know
cameras intrinsic parameters, but their extrinsic parameters as well as 3D
points are estimated from the observed points in different views. Knowing
intrinsic calibration, one sees improvement over the accuracy and robustness
of the structure and motion estimates nister2004efficient ;
kukelova2008polynomial .

We will show that when uncertainty exists in camera extrinsic parameters, the
$L_{2}$ triangulation method is no longer the state-of-the-art method.
Interestingly, a simple mid-point method, i.e. the mid-point of lines of
sights in different views, works much better in practice. Fig. 1 depicts the
difference between $L_{2}$ and mid-point triangulation methods. The mid-point
method not only gives better results but it can also be generalized to any
number of views with no difficulties. The $L_{2}$ triangulation method is
normally used for two views where the roots of a polynomial of degree $6$
needs to be computed hartley1997triangulation . For more than two views, the
method becomes computationally expensive and hard. For example, in the case of
three views, the optimal solution is one of the real solutions among the set
of $47$ general roots of a polynomial equation byrod2007improving . The mid-
point method in not projective or affine invariant hartley1997triangulation ,
but this lack of invariance is not important for calibrated reconstruction.

Through extensive simulations, both on synthetic and real datasets, we have
validated the high performance of the mid-point method. The performance is
defined as the accuracy of reconstruction where the optimal similarity
transform is applied to the reconstructed 3D points. We have assessed the
performance when uncertainty in relative cameras poses exists. We have also
assessed the performance when uncertainty is caused by commonly used
structure-from-motion approach, where first the essential matrix111The
essential matrix corresponding to a pair of cameras with relative orientation,
$R$, and translation, $t$, is defined as $E=[t]_{\times}R$. is estimated from
point correspondences and then relative poses are estimated
nister2004efficient  and finally the 3D structure is obtained. In all of these
experiments, we see that the mid-point method outperforms other approaches.
Thus, we suggest that unlike the common practice, the mid-point method should
be used in structure-from-motion applications.

##  2 Related Works

The $L_{2}$ triangulation approach, which is known as the optimal method,
finds the 3D point that minimizes the $L_{2}$ reprojection errors in the image
domain kanatani1996statistical . This leads to finding the optimal, maximum-
likelihood, solution under the assumption of Gaussian noise in the position of
projections. hartley1997triangulation  showed that minimizing the $L_{2}$
reprojection errors, for the case of two images, can be reduced to finding the
stationary points of $6^{th}$ degree polynomial and selecting the best by
evaluating the objective function. A Gröbner basis based algorithm for
minimizing the $L_{2}$ reprojection errors, in the case of three image
observations, is proposed by stewenius2005hard . They showed that the optimal
solution is one of the real solutions among a set of $47$ general roots of a
certain polynomial equation. Since their approach has a significant
computational cost, an alternative method of byrod2007improving  can be used,
which presents techniques that improves the numerical stability of Gröbner
basis solvers and significantly reduces the computational costs.

Because of the non-convexity and complexity of solving the $L_{2}$ norm
hartley2013verifying , other cost functions were considered in the literature.
For instance, a choice which is robust to outliers is to minimize the $L_{1}$
reprojection errors. hartley1997triangulation  find the $L_{1}$ optimal
solution in closed form by solving a polynomial of degree $8$. They also state
that the $L_{1}$ optimization has slightly less 3D error than the $L_{2}$
optimization in real experiments. Another popular approach is to find
$L_{\infty}$ answer which is optimal under the assumption of uniform noise
hartley2007optimal ; hartley2004sub ; li2007practical ; olsson2010outlier ;
sim2006removing . Angular errors were studied in lee2019closed  and closed-
form optimal solutions were derived for $L_{1}$, $L_{2}$, and $L_{\infty}$
angular errors.

In this paper the accuracy of different triangulation approaches in a
calibrated structure-from-motion process toldo2015hierarchical ;
ramalingam2006generic  is investigated. It is shown that the 3D baseline
triangulation approach has less sensitivity to uncertainties in cameras
extrinsic parameters and also has more accuracy in 3D reconstruction in a full
structure-from-motion process where positions and orientations of the cameras
are estimated from observations.

##  3 Preliminaries

###  3.1 Camera Model and Parameters

Let $\mathbf{u}$ be the projection of a point ${\mathbf{p}}$ on a camera’s
image plane. The projection is obtained by $\mathbf{u}=P[{\mathbf{p}};1]$,
where $P$ is the camera matrix. The camera matrix $P$ is given by
$P=K[R\leavevmode\nobreak\ |-R\mathbf{c}]$ , where $K$ is the camera
calibration matrix, and $R$ and $\mathbf{c}$ are the orientation and position
of the camera with respect to a world coordinate system. The line of sight of
the camera image is the line that passes through the camera point $\mathbf{c}$
towards direction $R^{-1}K^{-1}{\mathbf{u}}$.

###  3.2 Optimal Triangulation

2D baseline triangulation, which is known as the optimal triangulation in the
literature, finds a $3D$ point $\mathbf{p}$ so that its projected points on
the cameras, $\hat{\mathbf{u}}_{i},i\in\\{1,\ldots,N_{c}\\}$, have the minimum
sum Euclidean distance from measurements $\mathbf{u}_{i}$s. Hence, it
minimizes the following cost function:

$f(\mathbf{p})=\sum_{i=1}^{N_{c}}d(\hat{\mathbf{u}}_{i},\mathbf{u}_{i})^{2},\quad\hat{\mathbf{u}}_{i}=P[{\mathbf{p}};1],\
i\in\\{1,\ldots,N_{c}\\},$ (1)

in which $d(\hat{\mathbf{u}}_{i},\mathbf{u}_{i})$ is the Euclidean distance
between the projected point and its measurement in the $i^{th}$ image.
Assuming independent Gaussian noise in the image domain and known cameras
positions and orientations, this method provides the maximum-likelihood
estimation of 3D points.

###  3.3 Mid-Point Triangulation

Another simple triangulation method is to find a $3D$ point $\mathbf{p}$ that
minimizes $3D$ distances from the lines of sights. The goal of this method is
to minimize the following cost function:

$f(\mathbf{p})=\sum_{i=1}^{N_{c}}d(\mathbf{p},\mathbf{l}_{i})^{2},$ (2)

in which $\mathbf{l}_{i}$s are the lines of sights and
$d(\mathbf{p},\mathbf{l}_{i})$ is the distance between $\mathbf{p}$ and
$\mathbf{l}_{i}$. For any number of cameras, minimizing (2) is a linear least
squares problem and can be solved in a closed form.

###  3.4 Accuracy of the reconstruction

Point cloud reconstructed by a structure-from-motion procedure is obtained up
to a scaled Euclidean transformation (a more general projective ambiguity
exists in the uncalibrated approach). Suppose that
$\hat{\mathbf{p}}_{i},i\in\\{1,\ldots,N\\}$ are estimated points, and
$\mathbf{p}_{i}$s are the ground truth (with known correspondences). As shown
in Fig. 2, the accuracy of the estimation is obtained by finding a scaled
Euclidean transformation such that the estimated point cloud are aligned to
the ground truth as much as possible. Mathematically speaking

$\min_{R,\mathbf{t},s}\sum_{i=1}^{N}d\bigl{(}sR\hat{\mathbf{p}}_{i}+\mathbf{t},\mathbf{p}_{i}\bigr{)}^{2},$
(3)

where $R$, $\mathbf{t}$, and $s$ are the rotation matrix, translation, and
scale parameters of the scaled Euclidean transformation.

Figure 2: The ground truth point cloud (green circles) and the estimated point
clouds (blue squares and red hexagons). As it is obvious in (a), the blue
points have a smaller sum squared distance to the ground truth versus the red
points. But after finding the best rotation, translation, and scale for both
red and blue points, done for (b), it reveals that the red points are the
better estimation of the ground truth green points.

###  3.5 Reconstruction in a Structure-from-Motion Process

A typical structure-from-motion framework ramalingam2006generic ;
schonberger2016structure ; sweeney2015optimizing  comprises the following
steps:

  1. 1.

Pairwise images registration

     * •

Feature extraction and matching lowe1999object ; bay2008speeded ; lo2009local

     * •

Finding relative rotations and translations between all pairs of images with
the matched features torr2000mlesac ; nister2004efficient ;
kukelova2008polynomial

  2. 2.

Camera pose estimation

     * •

Solving the viewing graph created by the pairwise image registrations to find
camera positions and orientations ozyesil2015robust ; sweeney2015optimizing ;
jiang2013global ; zhu2018very ; hartley2013rotation ; chatterjee2017robust ;
arrigoni2018robust

  3. 3.

Triangulation

     * •

Reconstructing 3D points by triangulating corresponding points.

The obtained camera poses and 3D reconstructed points are usually refined by a
bundle adjustment step. In this paper, we only consider up to the
triangulation adjustment and exclude the bundle adjustment step.

##  4 Experiments

The experiments are two-fold. First, the sensitivity of different
triangulation methods to the error in camera poses are compared. Then, the
accuracy of different triangulation methods are evaluated in a full
reconstruction procedure on synthetic and real datasets.

###  4.1 Sensitivities

In this part, the mid-point method is compared to the $L_{2}$ method
hartley1997triangulation , iteratively reweighted mid-point
yang2019iteratively , minimizing $L_{1}$ reprojection errors
hartley1997triangulation , and minimizing $L_{1}$ and $L_{2}$ angular errors
lee2019closed .The comparison criteria are:

  * •

Position error sensitivity: The sensitivity of the error in the position of a
single triangulated point.

  * •

Distance error sensitivity: The sensitivity of the error in the distance
between two triangulated points.

  * •

Angle error sensitivity: The sensitivity of the error in the angles of a
triangle composed of three triangulated points.

To evaluate the mentioned sensitivities for different triangulation methods,
three configuration are considered for two-views triangulation:

  1. Conf. 1) $\mathbf{c}_{1}=[-5,-1,0]^{T},\mathbf{c}_{2}=[-5,+1,0]^{T}$ and the both cameras point at origin.

  2. Conf. 2) $\mathbf{c}_{1}=[-12,0,0]^{T},\mathbf{c}_{2}=[-2,0,0]^{T}$ and the both cameras point at origin.

  3. Conf. 3) $\mathbf{c}_{1}=[-10,2,-1]^{T},\mathbf{c}_{2}=[-5,-2,1]^{T}$ and the both cameras baselines are aligned with the global coordinate $x$ direction.

The first two configurations are the same as the two configurations in
hartley1997triangulation . In fact, the first configuration simulates a camera
moving forward and looking straight ahead and the second configuration
simulates an aerial imaging procedure. Since these two configurations are
special cases, a more general configuration is added to them.

####  4.1.1 Position Error Sensitivity

In this part the sensitivity of $3D$ error of different triangulation methods
to errors in positions and orientations of cameras are evaluated. For each
configuration, a point $\mathbf{p}$ is placed in an sphere centered at the
origin with the diameter of $0.5$, and the projected points on the cameras are
obtained. The positions of the cameras are perturbed by random Gaussian noise
vectors and the orientation of the cameras are perturbed by random rotations
with random axes and Gaussian random angles. The reconstructed point
$\hat{\mathbf{p}}$, is obtained by different triangulation methods. Euclidean
distance between points is used for computing $3D$ error:

$e=d(\hat{\mathbf{p}},\mathbf{p}).$ (4)

This procedure repeated $100$ times for each configuration and for each noise
level. The standard deviation of the position Gaussian noise is $0.01$, and
the standard deviation of the angle Gaussian noise is $0.1$ degree for the
noise level $1$. The standard deviations are multiplied by the noise levels.
Fig. 3 shows the mean error of different triangulation methods in confs. 1, 2,
and 3.

Figure 3: The position error (4) sensitivity of different methods in - from
top to bottom - configurations 1, 2, 3.

####  4.1.2 Distance Error Sensitivity

To evaluate the distance error sensitivity of different methods, two points
$\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ are randomly placed in the sphere of the
previous part. The projections are computed, cameras positions and
orientations are perturbed by the noise, and different triangulation methods
are applied to find two estimated $3D$ points $\hat{\mathbf{p}}_{1}$ and
$\hat{\mathbf{p}}_{2}$. The error is the absolute value of difference of the
distance between $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$, and the distance
between $\hat{\mathbf{p}}_{1}$ and $\hat{\mathbf{p}}_{2}$:

$e=\left|d(\mathbf{p}_{1},\mathbf{p}_{2})-d(\hat{\mathbf{p}}_{1},\hat{\mathbf{p}}_{2})\right|.$
(5)

This procedure is repeated $100$ times for each configuration and for each
noise level. The mean error of different triangulation methods are shown in
Fig. 4.

Figure 4: The distance error (5) sensitivity of different methods in - from
top to bottom - Configurations 1, 2, 3.

####  4.1.3 Angle Error Sensitivity

To compute the angle error sensitivity of different methods, three points
$\mathbf{p}_{1}$, $\mathbf{p}_{2}$, and $\mathbf{p}_{3}$ are randomly placed
in the aforementioned sphere and the projections are computed. The positions
and orientations of the cameras are perturbed by the noise and three $3D$
points $\hat{\mathbf{p}}_{1}$, $\hat{\mathbf{p}}_{2}$, $\hat{\mathbf{p}}_{3}$
are estimated by different triangulation methods. The error is defined as the
absolute value of difference between the angle between two vectors
$\mathbf{p}_{2}-\mathbf{p}_{1}$ and $\mathbf{p}_{3}-\mathbf{p}_{1}$, and the
angle between two vectors $\hat{\mathbf{p}}_{2}-\hat{\mathbf{p}}_{1}$ and
$\hat{\mathbf{p}}_{3}-\hat{\mathbf{p}}_{1}$:

$e=\left|\measuredangle(\hat{\mathbf{p}}_{2}-\hat{\mathbf{p}}_{1},\hat{\mathbf{p}}_{3}-\hat{\mathbf{p}}_{1})-\measuredangle(\mathbf{p}_{2}-\mathbf{p}_{1},\mathbf{p}_{3}-\mathbf{p}_{1})\right|.$
(6)

This procedure is repeated $100$ times for each configuration and for each
noise level. The mean error of different triangulation methods are shown in
Fig. 5.

Figure 5: The angle error (6) sensitivity of different methods in - from top
to bottom - Configurations 1, 2, 3.

From these experiments, it can be concluded that the mid-point method and its
variant are the best performing methods when there is uncertainty in the
cameras parameters. In reality, there are uncertainty in both cameras
parameters and image points. This is addressed in the following subsection.

###  4.2 Full Reconstruction Procedure on Synthetic Datasets

In the following experiments, the performance of the above methods in a full
structure-from-motion reconstruction procedure is assessed. The datasets have
uncertainties in image points and consequently there is uncertainty in cameras
extrinsic parameters. In this section, first the performance in the case of
two cameras is evaluated, where camera poses are computed through the
essential matrix. Then, the performance for the case of more than two cameras
are evaluated, where an additional viewing graph optimization is needed to be
solved for computing the camera poses.

####  4.2.1 Two Cameras

To evaluate the performance of different triangulation methods for the case of
two cameras, the following steps are done on the simulation setup of Fig. 6:

  1. 1.

$N=20$ points are randomly selected in a box.

  2. 2.

The points are projected on the cameras.

  3. 3.

The projections are displaced by an unbiased Gaussian random noises with the
standard deviation of one pixel.

  4. 4.

The essential matrix between the two cameras is estimated by the method of
kukelova2008polynomial .

  5. 5.

The relative rotation and translation are estimated.

  6. 6.

The cameras poses are calculated from the relative observations.

  7. 7.

The corresponding $3D$ points are reconstructed by triangulation.

  8. 8.

The best rotation, translation, and scale that makes the $3D$ triangulated
points match the selected points in the box are obtained and the errors are
computed.

The procedure is repeated $100$ times.

Figure 6: The configuration of the cameras and the box, and the images of the
box on the cameras. The cameras are placed at
$\mathbf{c}_{1}=[-7,3,0]^{T},\mathbf{c}_{2}=[-10,-3,1]^{T}$. Both cameras
looking at the origin. The dimensions of the box are $3\times 8\times 6$ and,
it is centered at the origin.

The cameras have the same calibration matrix

$K=\begin{bmatrix}300&0&320\\\ 0&300&240\\\ 0&0&1\end{bmatrix},$ (7)

and the images have $640\times 480$ pixels.

Fig. 7 shows the mean and standard deviation of triangulation errors of all
$20$ points for different methods in the first $10$ experiments. Fig. 8 shows
the mean, median, standard deviation, minimum, and maximum of the mean error
of all $100$ experiments.

Figure 7: The mean and standard deviation of 3D triangulation errors for
different methods in the first ten random runs on the synthetic dataset. The
errors are computed for random 20 points selected in a region. The cameras and
the region containing the points are shown in Fig. 6.  Figure 8: The mean,
standard deviation, median, minimum, and maximum of 100 mean 3D triangulation
error for different methods. The errors are computed for random 20 points
selected in a region. The cameras and the region containing the points are
shown in Fig. 6.

####  4.2.2 More Than Two Cameras

If $N_{c}$ cameras ($N_{c}>2$) are involved in the reconstruction process, the
essential matrices and consequently the relative positions and orientations
are computed for every $2$-combinations of $N_{c}$ cameras. The
$\binom{N_{c}}{2}$ relative observations of orientations and directions create
a viewing graph which should be solved to estimate the cameras poses. In the
$4^{th}$ step of the reconstruction procedure, the essential matrices are
estimated for any $2$-combinations of $N_{c}$ cameras, and in $5^{th}$ step, a
viewing graph is created from relative measurements. Obtained viewing graph is
solved in step $6$ to estimate the positions and orientations of the cameras.

To evaluate the performance of different methods on multi-view triangulation,
another camera is added to the mentioned two cameras setup as shown in Fig. 9.
The mean and standard deviation of triangulation error of all $20$ points for
different methods in the first $10$ experiments are shown in Fig. 10. Again
the mean, median, standard deviation, minimum, and maximum of the mean error
of all $100$ experiments are computed and are shown in Fig. 11.

###  4.3 Full Reconstruction Procedure in a Real Dataset

In this part, “Fountain-P11” dataset is used to evaluate the triangulation
methods. The SURF feature correspondences bay2008speeded  are used to find the
essential matrices between all $2$-combinations of cameras222Fountain-P11
dataset consists of $11$ images from different perspectives. We remove the
first and last images, which have a few number of feature correspondences and
find the essential matrices for all $2$-combinations of all other $9$ images..
The test process is the same as the process in the synthetic datasets except
that the projection in step $2$ is replaced by the feature matching and no
noise is added to the points anymore. For each selected pair of cameras, the
process is repeated $10$ times for different random corresponding points.
Figs. 12 and 13 show the results of different triangulation methods in the
mentioned process.

The experiment is repeated for triangulation by three cameras. Again the
process is the same as the process in the synthetic datasets of subsection
4.2.2 with feature matching used in step 2. The process is repeated $10$ time
for any $3$-combinations of cameras. The results are shown in Fig. 14. As it
can be seen in the results, the mid-point method outperforms the other methods
and has less mean 3D reconstruction error in all the experiments.

Figure 9: The configuration of the box and cameras $\mathbf{c}_{1}$ and
$\mathbf{c}_{2}$ are the same as Fig. 6, and the third camera is placed at
$\mathbf{c}_{3}=[-8,0,-2]^{T}$. Figure 10: The mean and standard deviation of
3D triangulation errors for different methods in the first ten random runs on
the synthetic dataset. The errors are computed for random 20 points selected
in a region. The cameras and the region containing the points are shown in
Fig. 9.  Figure 11: The mean, standard deviation, median, minimum, and maximum
of 100 mean 3D triangulation error for different methods. The errors are
computed for random 20 points selected in a region. The cameras and the region
containing the points are shown in Fig. 9.  Figure 12: The mean triangulation
error of different methods for each selected pairs of cameras in all
2-combination of cameras in Fountain-P11 dataset. The results are for 10 runs
of the triangulation procedure, and in each run 20 random corresponding points
using SURF features are used.  Figure 13: The mean, standard deviation,
median, minimum, and maximum of mean 3D triangulation error of all 2-view
experiments on Fountain-P11 dataset.  Figure 14: The mean, standard deviation,
median, minimum, and maximum of mean 3D triangulation error of all 3-view
experiments on Fountain-P11 dataset.

##  5 Conclusion

In this paper, different triangulation methods were evaluated in terms of 3D
reconstruction accuracy in a calibrated structure-from-motion setting. It was
shown that the mid-point triangulation method, which has a closed-form
solution for any number of cameras, is less sensitive to error in the cameras
extrinsic parameters in comparison to the other methods. This results in a
better performance of this triangulation method in structure-from-motion
procedures. The performance of different methods in a structure-from-motion
process were evaluated in synthetic and real datasets through extensive
experiments. It was shown that the mid-point triangulation method outperforms
the commonly used $L_{2}$ triangulation method hartley1997triangulation  in
typical practical applications, where cameras extrinsic parameters are
computed based on image registration and consequently have uncertainties.

## References

  * (1) Arrigoni, F., Rossi, B., Fragneto, P., Fusiello, A.: Robust synchronization in SO(3) and SE(3) via low-rank and sparse matrix decomposition.  Computer Vision and Image Understanding 174, 95–113 (2018)
  * (2) Bartoli, A., Sturm, P.: Structure-from-motion using lines: Representation, triangulation, and bundle adjustment.  Computer Vision and Image Understanding 100(3), 416–441 (2005)
  * (3) Bay, H., Ess, A., Tuytelaars, T., Van Gool, L.: Speeded-up robust features (SURF).  Computer Vision and Image Understanding 110(3), 346–359 (2008)
  * (4) Byröd, M., Josephson, K., Åström, K.: Improving numerical accuracy of Gröbner basis polynomial equation solvers.  In: IEEE International Conference on Computer Vision, pp. 449–456. IEEE (2007)
  * (5) Castle, R.O., Klein, G., Murray, D.W.: Wide-area augmented reality using camera tracking and mapping in multiple regions.  Computer Vision and Image Understanding 115(6), 854–867 (2011)
  * (6) Chatterjee, A., Govindu, V.M.: Robust relative rotation averaging.  IEEE Transactions on Pattern Analysis and Machine Intelligence 40(4), 958–972 (2017)
  * (7) Hartley, R., Kahl, F.: Optimal algorithms in multiview geometry.  In: Asian Conference on Computer Vision, pp. 13–34. Springer (2007)
  * (8) Hartley, R., Kahl, F., Olsson, C., Seo, Y.: Verifying global minima for L2 minimization problems in multiple view geometry.  International Journal of Computer Vision 101(2), 288–304 (2013)
  * (9) Hartley, R., Schaffalitzky, F.: L∞ minimization in geometric reconstruction problems.  In: IEEE Conference on Computer Vision and Pattern Recognition, vol. 1, pp. 504–509. IEEE (2004)
  * (10) Hartley, R., Trumpf, J., Dai, Y., Li, H.: Rotation averaging.  International Journal of Computer Vision 103(3), 267–305 (2013)
  * (11) Hartley, R.I., Sturm, P.: Triangulation.  Computer Vision and Image Understanding 68(2), 146–157 (1997)
  * (12) Jiang, N., Cui, Z., Tan, P.: A global linear method for camera pose registration.  In: IEEE International Conference on Computer Vision, pp. 481–488 (2013)
  * (13) Kanatani, K.: Statistical optimization for geometric computation: theory and practice.  Elsevier, New York, USA (1996)
  * (14) Kukelova, Z., Bujnak, M., Pajdla, T.: Polynomial eigenvalue solutions to the 5-pt and 6-pt relative pose problems.  In: British Machine Vision Conference, vol. 2, pp. 56.1–56.10 (2008)
  * (15) Lee, S.H., Civera, J.: Closed-form optimal two-view triangulation based on angular errors.  In: IEEE International Conference on Computer Vision, pp. 2681–2689 (2019)
  * (16) Li, H.: A practical algorithm for L∞ triangulation with outliers.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–8. IEEE (2007)
  * (17) Lindstrom, P.: Triangulation made easy.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 1554–1561. IEEE (2010)
  * (18) Lo, T.W.R., Siebert, J.P.: Local feature extraction and matching on range images: 2.5D SIFT.  Computer Vision and Image Understanding 113(12), 1235–1250 (2009)
  * (19) Lowe, D.G.: Object recognition from local scale-invariant features.  In: IEEE International Conference on Computer Vision, vol. 2, pp. 1150–1157. IEEE (1999)
  * (20) Nister, D.: An efficient solution to the five-point relative pose problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence 26(6), 756–770 (2004)
  * (21) Olsson, C., Eriksson, A., Hartley, R.: Outlier removal using duality.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 1450–1457. IEEE (2010)
  * (22) Ozyesil, O., Singer, A.: Robust camera location estimation by convex programming.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 2674–2683 (2015)
  * (23) Ramalingam, S., Lodha, S.K., Sturm, P.: A generic structure-from-motion framework.  Computer Vision and Image Understanding 103(3), 218–228 (2006)
  * (24) Schonberger, J.L., Frahm, J.M.: Structure-from-motion revisited.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 4104–4113 (2016)
  * (25) Sim, K., Hartley, R.: Removing outliers using the L∞ norm.  In: IEEE Conference on Computer Vision and Pattern Recognition, vol. 1, pp. 485–494. IEEE (2006)
  * (26) Stewenius, H., Schaffalitzky, F., Nister, D.: How hard is 3-view triangulation really?  In: IEEE International Conference on Computer Vision, vol. 1, pp. 686–693 (2005)
  * (27) Sweeney, C., Sattler, T., Hollerer, T., Turk, M., Pollefeys, M.: Optimizing the viewing graph for structure-from-motion.  In: IEEE International Conference on Computer Vision, pp. 801–809 (2015)
  * (28) Tippetts, B., Lee, D.J., Lillywhite, K., Archibald, J.: Review of stereo vision algorithms and their suitability for resource-limited systems.  Journal of Real-Time Image Processing 11(1), 5–25 (2016)
  * (29) Toldo, R., Gherardi, R., Farenzena, M., Fusiello, A.: Hierarchical structure-and-motion recovery from uncalibrated images.  Computer Vision and Image Understanding 140, 127–143 (2015)
  * (30) Torr, P.H., Zisserman, A.: Mlesac: A new robust estimator with application to estimating image geometry.  Computer Vision and Image Understanding 78(1), 138–156 (2000)
  * (31) Yang, K., Fang, W., Zhao, Y., Deng, N.: Iteratively reweighted midpoint method for fast multiple view triangulation.  IEEE Robotics and Automation Letters 4(2), 708–715 (2019)
  * (32) Zhu, S., Zhang, R., Zhou, L., Shen, T., Fang, T., Tan, P., Quan, L.: Very large-scale global sfm by distributed motion averaging.  In: IEEE Conference on Computer Vision and Pattern Recognition, pp. 4568–4577 (2018)


††institutetext: Physics Department, Moscow State University, Moscow, 119991,
Russia

# Nonperturbative QED vacuum birefringence

V.I.Denisov  E.E.Dolgaya  and V.A.Sokolov<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we represent nonperturbative calculation for one-loop Quantum
Electrodynamics (QED) vacuum birefringence in presence of strong magnetic
field. The dispersion relations for electromagnetic wave propagating in strong
magnetic field point to retention of vacuum birefringence even in case when
the field strength greatly exceeds Sauter-Schwinger limit. This gives a
possibility to extend some predictions of perturbative QED such as
electromagnetic waves delay in pulsars neighbourhood or wave polarization
state changing (tested in PVLAS) to arbitrary magnetic field values. Such
expansion is especially important in astrophysics because magnetic fields of
some pulsars and magnetars greatly exceed quantum magnetic field limit, so the
estimates of perturbative QED effects in this case require clarification.

###### Keywords:

QED vacuum birefringence, Electromagnetic processes in strong field,
Nonperturbative effects

##  1 Introduction

From the very beginning, radiative corrections to the QED Lagrangian coming
from the vacuum fluctuations have been the subject of a great interest. Such
corrections, accounted on background of constant and homogeneous
electromagnetic field, lead to Heisenberg-Euler effective Lagrangian a1 . The
scale parameter for nonlinearities in Heisenberg-Euler electrodynamics
(characteristic quantum electrodynamic induction or Schwinger limit for
electric field) $B_{c}=E_{c}=m^{2}c^{3}/e\hbar=4.41\cdot 10^{13}$G
distinguishes different regimes of the theory.

The perturbative or post-Maxwellian regime takes place for most of the
available in contemporary laboratory electromagnetic fields $E,B\ll
E_{c},B_{c}$. In this case the Heisenberg-Euler Lagrangian can be expanded in
power series of electromagnetic field tensor invariants. The perturbative
regime of QED is well understood and many of its predictions have been
experimentally observed. For instance, electron anomalous magnetic moment
remains a great example of unprecedented correspondence between theoretical
and experimental results a2 ; a3 . Some other predictions such as the Delbrück
scattering a4 , Lamb shift a5  and photon splitting a6 ; a7  are also well-
established. Another manifestation of perturbative QED lies in the phenomenon
of vacuum behavior like a polarizable medium with the cubic nonlinearity in
constitutive relations. Such semiclassical description leads to birefringence
of electromagnetic waves in vacuum in presence of external field. The attempt
of vacuum birefringence experimental observation was made in PVLAS experiment
a8 ; a9 . The measurement of refractive indexes induced by vacuum
birefringence was based on detecting rotation of electromagnetic wave
polarization when it propagated in dipole magnetic field with induction
$B_{ext}=2.5\cdot 10^{4}$G. The results of PVLAS set a new limit on vacuum
magnetic birefringence above the level pointed out by QED a10 ; a9 . According
to the authors, the discrepancy can be explained beyond the Standard Model by
interaction with axiones. The experiment also bounded the coupling constant of
axion-like particles and photons a10 ; a9 . Since the vacuum birefringence is
a very small macroscopic quantum effect it’s detection needs strong enough
magnetic field source, which is difficult to obtain in laboratory. At the same
time, usage of compact astrophysical objects as natural magnetic fields
sources provides wide opportunities for vacuum birefringence investigation due
to the fact that for many pulsars and manetars the magnetic field is close to
or even exceeds $B_{c}$. Vacuum birefringence in this case can be detected by
measuring X- and gamma- ray polarization passing the region of the strong
magnetic field near the pulsar. Due to the difference in wave propagation
velocities induced by vacuum birefringence the time lag between the arrival of
the fast and the slow mode to the detector is proportional to vacuum
refractive index difference for each polarization mode. The calculations in
perturbative QED regime show the detectability of the effect a11  and give the
value for the time lag $\Delta t\sim 10^{-7}$s. The results of calculations
are valid when the pulsar field $B<B_{c}$, however the existence of pulsars
with the overcritical field (for instance B1509-58 with the $B\sim 1.5\cdot
10^{14}$G) provides an attractive possibility to enhance the estimates for the
time lag. Such calculations require of vacuum birefringence consideration
outside the perturbative regime.

Nonperturbative regime of QED arises when we consider external fields which
values are close to or exceed the scale parameter of the Heisenberg-Euler
electrodynamics $E,B\sim E_{c},B_{c}$. This regime shows one of the most
amazing properties – vacuum instability due to electron-positron pair
production. The effect is expected a14  when the electric field exceeds so-
called Schwinger limit $E>E_{c}$ and, although it has never been observed
directly, the advances of high-intensity laser physics and implementation of
such projects as ELI, XFEL and others a12 ; a13  give a great promise in this
area. In general, QED in nonperturbative regime is poorly investigated and
provides new challenges both in theoretical and experimental research.

In this paper we will focus on one-loop QED vacuum birefringence in
nonperturbative regime. In addition to entirely quantum technic proposed in
a15 ; a16 , we represent a semiclassical approach. The maim aim of our
research is to expand vacuum birefringence predictions for different
experiments to values of strong magnetic fields close to character quantum
electrodynamic induction. The paper is organized as follows: in Sec.2 we
derive constitutive relations for QED, Sec.3 is devoted to the weak
electromagnetic wave propagation in strong magnetic field, in Sec.4 we discuss
QED birefringence expansion following from obtained dispersion relations and
possibilities for it’s observation and in conclusive Sec.5 we summarize our
results.

##  2 Constitutive relations for nonperturbative QED

The Lagrangian in QED is represented as a series of corrections to Maxwell
electrodynamics. One-loop QED considers only the first non-vanishing
correction which in non-perturbative regime has the following form b1 :

$L_{1}=-{\alpha B_{c}^{2}\over 8\pi^{2}}\int\limits_{0}^{\infty}{e^{-s}\over
s^{3}}\Big{[}sa\ \cot(sa)\cdot sb\ \coth(sb)-1-{s^{2}\over
3}(b^{2}-a^{2})\Big{]}ds,$ (1)

where $\alpha=e^{2}/\hbar c$ – is a fine structure constant and the parameters
$a$ and $b$ are expressed by the electromagnetic field components:

$a=-{i\over\sqrt{2}B_{c}}\Big{(}\sqrt{F+iG}-\sqrt{F-iG}\Big{)},\quad
b={1\over\sqrt{2}B_{c}}\Big{(}\sqrt{F+iG}+\sqrt{F-iG}\Big{)},$ (2)

where the notations $F=({\bf B}^{2}-{\bf E}^{2})/2$ and $G=({\bf E}{\bf B})$
were used for brevity. Strictly speaking, the Heisenberg-Euler Lagrange
function is valid only in case of constant and homogeneous background fields,
at least at the typical scale of Compton wavelength $\lambda_{c}=h/mc$. Also
it should be noted that the auxiliary parameters $a$ and $b$ of Heisenberg-
Euler Lagrangian (1) are special because $\pm a$ and $\pm b$ are the
eigenvalues of the electromagnetic field tensor $F_{ik}$ when field value is
constant and this makes the problem of QED radiative corrections exactly
solvable arx1 .

In semiclassical approach Heisenberg-Euler theory can be interpreted as
nonlinear electrodynamics of continuous media with special constitutive
relations, in which polarization ${\bf P}$ and magnetization ${\bf M}$ induced
by the external fields ${\bf E}$ and ${\bf B}$, can be expressed from the
Lagrangian:

${\bf P}={\partial L_{1}\over\partial{\bf E}},\qquad{\bf M}={\partial
L_{1}\over\partial{\bf B}}.$ (3)

In order to calculate the explicit values for this vectors, it is useful to
introduce some auxiliary relations:

${\partial a\over\partial{\bf E}}={\partial b\over\partial{\bf B}}={{\bf
E}a+{\bf B}b\over 2\sqrt{F^{2}+G^{2}}},\quad{\partial a\over\partial{\bf
B}}=-{\partial b\over\partial{\bf E}}={{\bf E}b-{\bf B}a\over
2\sqrt{F^{2}+G^{2}}}.$ (4)

substitution of which to (3) finally leads to constitutive relations for
nonberturbative one-loop QED:

${\bf P}={\alpha\over 8\pi^{2}(a^{2}+b^{2})}\Big{[}I_{1}{\bf E}+I_{2}{\bf
B}\Big{]},\quad{\bf M}=-{\alpha\over 8\pi^{2}(a^{2}+b^{2})}\Big{[}I_{1}{\bf
B}-I_{2}{\bf E}\Big{]},$ (5)

where for brevity we have used the notations for the integrals:

$I_{1}=\int\limits_{0}^{\infty}\Big{\\{}{ab\big{[}a\sinh(2sb)-b\sin(2sa)\big{]}\over
2\sinh^{2}(sb)\sin^{2}(sa)}-{2(a^{2}+b^{2})\over 3s}\Big{\\}}e^{-s}ds,$ (6)
$I_{2}=\int\limits_{0}^{\infty}\Big{\\{}{ab\big{[}a\sin(2sa)+b\sinh(2sb)\big{]}\over
2\sinh^{2}(sb)\sin^{2}(sa)}-{(a^{2}+b^{2})\over
s}\cot(sa)\coth(sb)\Big{\\}}e^{-s}ds.$ (7)

It is easy to verify the correspondence to perturbative regime which take
place for relatively weak fields $|{\bf E}|,|{\bf B}|\ll B_{c}$. In this case,
QED correction to Lagrangian (1) and the constitutive relations (5) can be
expanded in a series by the small parameters $a,b\ll 1$. Such expansion leads
to:

$L_{1}={\alpha\over 8\pi B_{c}^{2}}\Big{[}{\eta_{1}({\bf E}^{2}-{\bf
B}^{2})^{2}+2\eta_{2}({\bf E}{\bf B})^{2}}\Big{]},$ (8) ${\bf P}={\xi\over
2\pi}\Big{\\{}\eta_{1}({\bf E}^{2}-{\bf B}^{2}){\bf E}+2\eta_{2}({\bf E}{\bf
B}){\bf B}\Big{\\}},\quad{\bf M}=-{\xi\over 2\pi}\Big{\\{}\eta_{1}({\bf
E}^{2}-{\bf B}^{2}){\bf B}-2\eta_{2}({\bf E}{\bf B}){\bf E}\Big{\\}},$ (9)

where $\xi=1/B_{c}^{2}$, $\eta_{1}=\alpha/45\pi$, and
$\eta_{2}=7\alpha/180\pi$ – so called post-Maxwellian parameters. As the
polarization and magnetization in (9) are cubic on external fields, the
electromagnetic waves and charged particles propagation in perturbative QED
vacuum will possess the properties peculiar to crystal optics with cubic
nonlinearity. This approach leads to the predictions for vacuum birefringence
and dichroism a6 ; a17 ; a18 ; a19 , optical non-reciprocity a20 , light-ray
bending a23 ; a24  and Cherenkov-radiation a21 ; a22  in vacuum at presence of
external electromagnetic field. The expansion of these predictions on
nonperturbative QED regime gives a new insight in understanding of vacuum
nonlinear electrodynamics and helps to enforce the expectations in
experimental manifestations. In this paper we will focus only on vacuum
birefringence expansion. In order to do this, we will derive the dispersion
relations for electromagnetic wave propagating on background of strong
magnetic field.

##  3 Electromagnetic wave propagation in strong magnetic field

Let us consider a weak electromagnetic wave ${\bf e}_{w}$, ${\bf b}_{w}$
propagating in strong external magnetic field ${\bf B}_{0}$. We suppose that
the field intensity in the wave is sufficiently weak, so $|{\bf e}_{w}|,|{\bf
b}_{w}|\ll B_{c},|{\bf B}_{0}|$. As the vacuum nonlinear electrodynamics keeps
the superposition principle, the total field intensities ${\bf B}={\bf
B}_{0}+{\bf b}_{w}$ and ${\bf E}={\bf e}_{w}$ can be used in (2) and in
constitutive relations for nonperturbative QED (5). Here we should take into
account the weakness of the wave and decompose these relations up to the
leading order by ${\bf e}_{w}$ and ${\bf b}_{w}$. Such decomposition gives
linearized constitution relations:

$\displaystyle{\bf D}$ $\displaystyle=$ $\displaystyle{\bf E}+4\pi{\bf P}={\bf
e}_{w}-2\xi\Big{\\{}\eta_{1}Y_{1}{\bf B}_{0}^{2}{\bf
e}_{w}-2\eta_{2}Y_{2}({\bf B}_{0}{\bf e}_{w}){\bf B}_{0}\Big{\\}}$ (10)
$\displaystyle{\bf H}$ $\displaystyle=$ $\displaystyle{\bf B}-4\pi{\bf M}={\bf
B}_{0}+{\bf b}_{w}-2\xi\eta_{1}\Big{\\{}Y_{1}{\bf B}_{0}^{2}{\bf
b}_{w}+2Y_{3}({\bf B}_{0}{\bf b}_{w}){\bf B}_{0}+Y_{1}{\bf B}_{0}^{2}{\bf
B}_{0}\Big{\\}},$

where we use linearized relations for the parameters $a$ and $b$

$a={({\bf B}_{0}{\bf e}_{w})\over B_{c}B_{0}},\quad b={B_{0}\over
B_{c}}\Big{\\{}1+{({\bf B}_{0}{\bf b}_{w})\over B_{0}^{2}}\Big{\\}},$ (11)

introduce the notations for the integrals:

$\displaystyle Y_{1}$ $\displaystyle=$ $\displaystyle-{45\over
4b_{0}^{2}}\int\limits_{0}^{\infty}{e^{-z/b_{0}}\over
z^{2}}\Big{\\{}\coth(z)-{z\over\sinh^{2}(z)}-{2z\over 3}\Big{\\}}dz,$ (12)
$\displaystyle Y_{2}$ $\displaystyle=$ $\displaystyle{45\over
14b_{0}^{2}}\int\limits_{0}^{\infty}{e^{-z/b_{0}}\over
z^{2}}\Big{\\{}{2z^{2}-3\over 3}\coth(z)+{z\over\sinh^{2}(z)}\Big{\\}}dz,$
(13) $\displaystyle Y_{3}$ $\displaystyle=$ $\displaystyle-{45\over
8b_{0}^{2}}\int\limits_{0}^{\infty}{e^{-z/b_{0}}\over
z^{2}}\Big{\\{}{2z^{2}\coth(z)-z\over\sinh^{2}(z)}-\coth(z)\Big{\\}}dz,$ (14)

and the dimensionless parameter ${\bf b}_{0}={\bf B}_{0}/B_{c}$. The
correspondence to the perturbative regime can be obtained when $B_{0}\ll
B_{c}$. This leads to the asymptotic expansion of the integrals:

$Y_{1}=1-{6\over 7}b_{0}^{2}+{16\over 7}b_{0}^{4}+\cdots,\quad
Y_{2}=1-{26\over 49}b_{0}^{2}+{176\over 147}b_{0}^{4}+\cdots,\quad
Y_{3}=1-{12\over 7}b_{0}^{2}+{48\over 7}b_{0}^{4}+\cdots.$ (15)

The graphs for an exact (12) and an approximate (15) $Y$ dependence on
magnetic field strength are represented on Figure 1. Approximate functions are
marked on the graph by the gray line. As it can be seen, even in low field
values $b_{0}\ll 1$ an approximate description can cause significant
inaccuracies whose rectification will require new terms in the expansion (15).

Also it should be noted that, as the unity is the leading term in all of the
listed above expansions, the linearized constitutive relations (10) contain
perturbative relations (9) in low field limit. Besides, the transition from
the perturbative QED to nonperturbative regime is in replacement of "post-
Maxwellian" constants $\eta_{1}$ and $\eta_{2}$ on the functions
$\zeta_{1}=\eta_{1}Y_{1}(b_{0})$, $\zeta_{2}=\eta_{2}Y_{2}(b_{0})$,
$\zeta_{3}=\eta_{1}Y_{3}(b_{0})$ which depend on the external field strength.

Figure 1:  An exact and an approximate $Y$-functions comparison.

To obtain dispersion relations we use the eikonal approximation and represent
the wave field in form: ${\bf e}_{w}={\bf e}\exp{\\{iS({\bf r},t)\\}},$ ${\bf
b}_{w}={\bf b}\exp{\\{iS({\bf r},t)\\}},$ where ${\bf e}$ and ${\bf b}$ are
the amplitudes and $S$ is the eikonal. As usual for the eikonal approximation,
we will suppose that amplitude variations along the wave propagation ray are
small and can be neglected in calculations. Substitution of ${\bf e}_{w}$,
${\bf b}_{w}$ and constitutive relations (10) to the ordinary equations of
continuous media electrodynamics leads to homogeneous equations on wave field
components:

$\Pi_{\alpha\beta}{e}_{w}^{\beta}=0,$ (16)

where indexes enumerate the cartesian components of the electric field vector
$\alpha,\beta=1\ldots 3$ and $\Pi_{\alpha\beta}$ is the polarization tensor:

$\displaystyle\Pi_{\alpha\beta}=\Big{\\{}\Big{(}({\bf\nabla}S)^{2}-({\partial_{0}S})^{2}\Big{)}\delta_{\alpha\beta}-{\partial_{\alpha}S}{\partial_{\beta}S}\Big{\\}}\times\Big{\\{}1-2\zeta_{1}b_{0}^{2}\Big{\\}}$
(17) $\displaystyle-4\zeta_{3}\Big{[}{\bf\nabla}S{\bf
b}_{0}\Big{]}_{\alpha}\Big{[}{\bf\nabla}S{\bf
b}_{0}\Big{]}_{\beta}-4\zeta_{2}({\partial_{0}S})^{2}({\bf
b}_{0})_{\alpha}({\bf b}_{0})_{\beta},$

in which the following notations were used:
$\partial_{0}=\partial/\partial(ct)$ and $\partial_{\alpha}$ denotes spatial
coordinate $x^{\alpha}$ derivative, $\delta_{\alpha\beta}$ – is the Kronecker
symbol, $\nabla$ – is the gradient operator, ${\bf b}_{0}={\bf B}_{0}/B_{c}$
and the square brackets refer to the vector cross-product. In should be noted
that our result for polarization tensor is close to the similar one obtained
in new1 . The existence of nontrivial solutions of (16) requires
$\mbox{det}\|\Pi_{\alpha\beta}\|=0$ which finally leads to dispersion
relations for the electromagnetic wave at the eikonal approximation:

$\displaystyle\Big{\\{}({\bf\nabla}S)^{2}-({\partial_{0}S})^{2}+4\zeta_{3}\Big{[}({\bf
b}_{0}{\bf\nabla}S)^{2}$ $\displaystyle-$
$\displaystyle{b}_{0}^{2}({\bf\nabla}S)^{2}\Big{]}+2\zeta_{1}{b}_{0}^{2}\Big{[}({\partial_{0}S})^{2}-({\bf\nabla}S)^{2}\Big{]}\Big{\\}}$
$\displaystyle\times\Big{\\{}({\bf\nabla}S)^{2}-({\partial_{0}S})^{2}+4\zeta_{2}\Big{[}({\bf
b}_{0}{\bf\nabla}S)^{2}$ $\displaystyle-$
$\displaystyle{b}_{0}^{2}({\partial_{0}S})^{2}\Big{]}+2\zeta_{1}{b}_{0}^{2}\Big{[}({\partial_{0}S})^{2}-({\bf\nabla}S)^{2}\Big{]}\Big{\\}}$
(18) $\displaystyle\times$
$\displaystyle\Big{\\{}2\zeta_{1}{b}_{0}^{2}-1\Big{\\}}=0.$

Multiplicative structure of obtained dispersion relations points to vacuum
birefringence retention even in case of one-loop nonperturbative QED. The
first factor in (3) corresponds to the dispersion law for a normal wave
polarized perpendicular to external magnetic field ($\bot$-mode), whereas the
second multiplier describes the mode polarized along the field ${\bf B}_{0}$
($||$-mode). The last factor in (3) does not depend on the wave parameters and
its equality to zero is expected at huge magnetic field intensities $B_{0}\sim
B_{c}\exp{(1/\alpha)}$. Just for these values QED vacuum instability induced
by magnetic field was predicted a25 . To avoid such a regime in future, we
will consider field intensities much closer to $B_{c}$.

Now let us investigate the properties of normal waves in more detail, and use
the obtained dispersion relations in strengthening some expectations for
vacuum birefringence detection in experiment.

##  4 QED vacuum birefringence extension on nonperturbative regime

There are two traditional approaches to the interpretation of the dispersion
relations for vacuum nonlinear electrodynamics. The first one comes from the
representation of the vacuum as a continuous media. The vacuum birefringence
in this case is explained by the normal modes refraction indexes mismatch
$n_{\bot}\neq n_{||}$. The second approach assumes that the electromagnetic
wave propagates in the space-time with the effective geometry following from
the dispersion relations. This interpretation explains the vacuum
birefringence due to the difference between the effective space-time metric
tensors $G^{ik}_{\bot}\neq G^{ik}_{||}$ for each normal mode. For
completeness, we use each of these approaches in the analysis of obtained
dispersion relations.

For description in terms of refractive indexes one should substitute eikonal
$S(t,{\bf r})=\omega t-({\bf k}{\bf r})$ to dispersion relations (3), and take
into account the relation between the wave vector ${\bf k}$ and frequency
$\omega$ which is ordinary for homogeneous wave in continuous media: ${\bf
k}=\omega n{\bf q}/c$, where $n$ is the refraction index and ${\bf q}$ is the
unity vector in wave propagation direction. Such substitution gives explicit
expressions for the normal modes refraction indexes:

$n_{\bot}^{2}=1+{4\zeta_{3}b_{0}^{2}\sin^{2}{\theta}\over
1-2b_{0}^{2}[\zeta_{1}+2\zeta_{3}\sin^{2}{\theta}]},\qquad
n_{||}^{2}=1+{4\zeta_{2}b_{0}^{2}\sin^{2}{\theta}\over
1-2b_{0}^{2}[\zeta_{1}-2\zeta_{2}\cos^{2}{\theta}]},$ (19)

where $\theta$ is an angle between the wave vector ${\bf k}$ and the magnetic
field ${\bf B}_{0}$. These expression refine the results obtained earlier a26
in which the dependance on angle $\theta$ is more simple and the terms in
denominator are neglected:

$n_{\bot}^{2}\approx 1+{4\zeta_{3}b_{0}^{2}\sin^{2}{\theta}},\qquad
n_{||}^{2}\approx 1+{4\zeta_{2}b_{0}^{2}\sin^{2}{\theta}},$ (20)

which are close to results of a16 ; a26  obtained on the basis of a purely
quantum approach with fixed selection of the gauge. As the angle dependance in
the exact expressions (19) is different, it becomes possible to figure out are
there any conditions under which $n_{\bot}=n_{||}$ and the birefringence is
suppressed. The equality of the refractive indexes leads to the relation which
is valid for any angle $\theta$:

$(\zeta_{3}-\zeta_{2})(1-2\zeta_{1}b_{0}^{2})+4\zeta_{3}\zeta_{2}b_{0}^{2}=0.$
(21)

However, the numerical calculations show that there are no solutions for this
equation at any arbitrary magnetic field close to $B_{c}$ so the vacuum
birefringence remains.

To analyze the dependence of refractive indexes on the magnetic field value,
we will choose $\theta=\pi/2$, which is more valuable for experimental
research because in this case indexes are maximal. As it follows from
numerical calculations, the refractive index $n_{||}$ increases almost
linearly in $1<b_{0}<100$ and shows nonlinear growth in wider field ranges.
Whereas the $n_{\bot}$ tends to saturation at the value
$(n_{\bot})_{sat}-1\approx 4\cdot 10^{-4}$ and ceases to depend on the field
strength. The dependance of $n_{\bot,||}-1$ for $b_{0}<3$ is represented on
the left graph in Figure 2. The right graph shows the exact difference
$n_{||}-n_{\bot}$ coming from (19) and the same difference following from the
perturbative QED: $n_{||}-n_{\bot}\approx 2(\eta_{2}-\eta_{1})b_{0}^{2}$,
which is marked on the graph with a gray line. Despite the proximity of
perturbative description, there is a good accordance to the exact result up to
the field values $b_{0}\approx 2$. This indicates that more then order
discrepancy between perturbative QED birefringence prediction and the
experimental result obtained in PVLAS a8  is not associated with the
inaccuracy of perturbative description and should have a more profound
physical reason.

Figure 2:  The refractive indexes and their difference depending on the
magnetic field value $b_{0}=B_{0}/B_{c}$.

Another traditional approach to the dispersion relations interpretation is
based on effective geometry representation. It assumes that the wave
propagates in curved space-time, which geometry depends on external magnetic
field. The dispersion relations (3), now are interpreted as Hamilton-Jacobi
equation for a massless particle in the effective space-time with the metric
tensor $G^{ik}_{\bot,||}$ correspondent to each normal mode:

$\Big{[}G^{ik}_{\bot}{\partial S\over\partial x^{i}}{\partial S\over\partial
x^{k}}\Big{]}\times\Big{[}G^{mn}_{||}{\partial S\over\partial x^{m}}{\partial
S\over\partial x^{n}}\Big{]}=0.$ (22)

As it follows from (3), the components of the effective metric tensor take the
form:

$\displaystyle G^{00}_{\bot}$ $\displaystyle=$ $\displaystyle g^{00},\quad
G^{\alpha\beta}_{\bot}=g^{\alpha\beta}+{4\zeta_{3}\over
1-2\zeta_{1}b_{0}^{2}}\times(b_{0}^{2}\delta^{\alpha\beta}-b_{0}^{\alpha}b_{0}^{\beta}),$
(23) $\displaystyle G^{00}_{||}$ $\displaystyle=$ $\displaystyle g^{00},\quad
G^{\alpha\beta}_{||}=g^{\alpha\beta}+{4\zeta_{2}\over
1-2\zeta_{1}b_{0}^{2}+4\zeta_{2}b_{0}^{2}}\times(b_{0}^{2}\delta^{\alpha\beta}-b_{0}^{\alpha}b_{0}^{\beta}),$
(24)

where $g^{ik}$ is Minkowski space-time metric tensor and the Greek indexes
enumerate the spatial coordinates, so $\alpha,\beta=1..3$. It is easy to
verify that in low-field limit $b_{0}\ll 1$ the expressions (23) take the form
of perturbative QED effective metric tensor obtained in a11 :

$G^{ik}_{\bot,||}\approx
g^{ik}+4\eta_{1,2}(b_{0}^{2}\delta^{\alpha\beta}-b_{0}^{\alpha}b_{0}^{\beta}).$
(25)

Such correspondence allows us to extend some predictions for vacuum
birefringence manifestations obtained on base of perturbative metric tensor
(25) to non-perturbative region by simple replacement of parameters:

$\eta_{1}\to{\zeta_{3}\over
1-2\zeta_{1}b_{0}^{2}},\quad\eta_{2}\to{\zeta_{2}\over
1-2\zeta_{1}b_{0}^{2}+4\zeta_{2}b_{0}^{2}}.$ (26)

For instance, such extension is especially actual for the pulsars and
magnetars which are one the most attractive objects for QED regime tests. The
magnetic fields of such astrophysical sources can significantly exceed
critical limit $b_{0}=B_{0}/B_{c}\gg 1$. One of observable manifestations of
the vacuum birefringence in pulsar neighbourhood is related to normal mode
delay for hard X-ray and gamma- radiation pulses passing near the pulsar.
Because of the vacuum birefringence, the propagation velocity of ${\bot}$-mode
is greater then ${||}$-mode, so it will reach to the detector earlier. So the
leading part of the pulse coming from the X-ray source to the detector will be
linearly polarized due to the ${\bot}$-mode. This part of the pulse will have
a time duration $\Delta t$. After this time the ${||}$-mode will reach the
detector and the pulse polarization state will change to elliptical. The
maximal estimation for the delay $\Delta t$ in perturbative QED was obtained
in a11 :

$\Delta t={123\pi(\eta_{2}-\eta_{1})b_{0}^{2}R_{s}\over 128c},$ (27)

where $R_{s}$ is the pulsar radius, and $c$ is the speed of light in Maxwell
vacuum. Birefringence expansion (26) allows us to estimate the order of
magnitude for the time delay $\Delta t$ in nonperturbative regime:

$\Delta t\simeq{b_{0}^{2}R_{s}\over c}\Big{[}{\zeta_{3}\over
1-2\zeta_{1}b_{0}^{2}}-{\zeta_{2}\over
1-2\zeta_{1}b_{0}^{2}+4\zeta_{2}b_{0}^{2}}\Big{]}.$ (28)

The dependence of the time delay on magnetic field strength and its estimates
for some pulsars and magnetars are represented on Figure 3. For the estimates
we use the pulsar data from McGill a27  and ATNF a28  catalogs and also
suppose the pulsar radius equal to $R_{s}=10km$.

Figure 3:  Time delay between ${\bot}$ and ${||}$-modes arrival to detector.

The delay value varies in hundredth of microseconds and this is sufficient for
contemporary timing measurements. Furthermore the estimate can be enforced for
the unique object J1808-2024 with the field strength $B_{0}\sim 2.06\cdot
10^{15}G$ for which the delay can reach $\Delta t\sim 0.3\mu s$. Also it
should be noted that there is almost linear dependence between delay and field
strength, instead of quadratic, typical for delay in perturbative regime (27),
therefore the direct interpolation of perturbative description on the region
$b_{0}>1$ will lead to overestimation for the delay.

##  5 Conclusion

In this paper we have investigated expansion of vacuum birefringence on
nonperturbative regime of QED. The obtained constitutive relations (10) for
weak electromagnetic wave propagating on background of a strong magnetic field
indicate that this expansion consists of replacing of perturbative constants
$\eta_{1},\eta_{2}$ by three functions $\zeta_{1},\zeta_{2},\zeta_{3}$, which
depend on the external field strength. For vacuum birefringence description we
have used a semiclassical approach in terms of the wave field strength which
gives gauge independent results unlike most of the other calculations
performed in the fixed gauge selection. In such approach the polarization
tensor (17) and dispersion relations (3) were obtained and both interpreted in
terms of refractive indexes and the effective space-time geometry. The
refractive indexes interpretation confirmed the results obtained earlier by
other authors on the base of quantum field theory methods. The comparison
between the perturbative and nonperturbative refractive indexes, as expected,
reveals an insignificant difference at weak field values $b_{0}\ll 1$, and
this indicates that inaccuracy in QED theoretical description can not be a
cause of discrepancy detected in the PVLAS experiment a9 . Another feature
observed in the refractive indexes analysis is the possible saturation of
$n_{||}$ at the strong external field values $b_{0}\gg 1$. Unfortunately, this
feature can not be verified in conditions of the terrestrial facilities,
however the astrophysical sources such as pulsars and magnetars provide a
wider opportunities in QED features investigations. It is more convenient to
use an effective geometry formalism for vacuum birefringence description in
the neighbourhood of such field sources. The dispersion relations
interpretation in terms of the effective geometry (23) allowed us to extend
the predictions for the normal mode relative delay in the hard radiation
propagating near the pulsar. The estimates for such delay for a different
pulsars, show its detectability with contemporary experimental technique.
Moreover, it was found out that the dependence between the delay and the field
value predicted by nonperturbative QED is close to linear (Figure 3) and can
be roughly expressed as $\Delta t=(5.8b_{0}-7)\times 10^{-3}\mu{s}$ for
$2<b_{0}<100$. This dependence differs from the quadratic one following from
the direct perturbative QED prediction extrapolation to the strong field
region. As the perturbative QED is not valid for such field values, it gives
overestimated result for delay. This feature can be noted in planning of
future astrophysical missions aimed at the QED effects investigation in the
pulsar fields, such as XIPE a30 .

## References

  * (1) W.Heisenberg and H.Euler, _Consequences of Dirac’s Theory of Positrons_ , _Z. Phys._ , 26, (1936) 714.
  * (2) D.Hanneke and S.Fogwell and G.Gabrielse, _New Measurement of the Electron Magnetic Moment and the Fine Structure Constant_ , _Phys. Rev. Lett._ , 100, (2008) 120801.
  * (3) M.Hayakawa, T.Kinoshita and M.Nio, _Revised value of the eighth-order QED contribution to the anomalous magnetic moment of the electron_ , _Phys. Rev. D_ , 77, (2008) 053012.
  * (4) Sh.Zh.Akhmadaliev et al. _Delbruck scattering at energies of 140-450 MeV_ , _Phys.Rev. C_ , 58, (1998) 2844.
  * (5) W.E.Lamb and R.C.Retherford _Fine Structure of the Hydrogen Atom by a Microwave Method_ , _Phys. Rev._ , 72, (1947) 241.
  * (6) S.L.Adler, _Photon splitting and photon dispersion in a strong magnetic field_ , _Ann. Phys._ , 67, (1971) 599.
  * (7) Sh.Zh.Akhmadaliev et al., _Experimental Investigation of High-Energy Photon Splitting in Atomic Fields_ , _Phys. Rev. Lett._ , 89, (2002) 061802.
  * (8) F.Della Valle et al, _The PVLAS experiment: measuring vacuum magnetic birefringence and dichroism with a birefringent Fabry-Perot_ , _Eur. Phys. J. C_ , 76, (2016) 24.
  * (9) F.Della Valle et al., _First results from the new PVLAS apparatus : a new limit on vacuum magnetic birefringence_ , _Phys Rev. D_ , 90, (2014) 092003.
  * (10) S.Villalba-Ch’avez and A.Di Piazza, _Axion-induced birefringence effects in laser driven nonlinear vacuum interaction_ , _Journal of High Energy Physics_ , 11, (2013) 136.
  * (11) V.I.Denisov and S.I.Svertilov, _Nonlinear electromagnetic and gravitational actions of neutron star fields on electromagnetic wave propagation_ , _Phys. Rev. D_ , 71, (2005) 063002.
  * (12) J.Schwinger, _On Gauge Invariance and Vacuum Polarization_ , _Phys. Rev._ , 82, (1951) 664.
  * (13) Sh.Zh.G.V. Dunne, _New strong-field QED effects at extreme light infrastructure_ , _Eur. Phys. J. D_ , 55, (2009) 327.
  * (14) S.Bulanov et al., _On the Schwinger limit attainability with extreme power lasers_ , _Phys. Rev. Lett._ , 105, (2010) 220407.
  * (15) K.Hattori and K.Itakura, _Vacuum birefringence in strong magnetic fields: (I) Photon polarization tensor with all the Landau levels_ , _Ann. of Phys._ , 330, (2013) 23.
  * (16) K.Hattori and K.Itakura, _Vacuum birefringence in strong magnetic fields: (II) Complex refractive index from the lowest Landau level_ , _Ann. of Phys._ , 334, (2013) 58.
  * (17) V.B.Berestetski, E.M.Lifshits, and L.P.Pitaevsky, _Quantum Electrodynamics_ , Pergamon, Oxford (1982).
  * (18) G.V.Dunne, _Heisenberg-Euler Effective Lagrangians: Basics and Extensions_ , arxiv:hep-th/0406216v1.
  * (19) W.Dittrich and H.Gies, _Probing the quantum vacuum. Perturbative effective action approach in quantum electrodynamics and its application_ , _Tracts Mod. Phys._ , 166, (2000).
  * (20) V.Denisov, V.Sokolov, and M.Vasili’ev, _Nonlinear vacuum electrodynamics birefringence effect in a pulsar’s strong magnetic field_ , _Phys. Rev. D_ , 90, (2014) 023011.
  * (21) R.Battesti and C.Rizzo, _Magnetic and electric properties of a quantum vacuum_ _Rept.Prog.Phys._ , 76, (2013) 016401.
  * (22) V.Denisov, _Nonlinear effect of quantum electrodynamics for experiments with a ring laser_ , _J. Opt. A: Pure Appl. Opt._ , 2, (2000) 372.
  * (23) J.Y.Kimand T.Lee, _Nonlinear electrodynamic effect of light bending by charged astronomical objects_ , _Nuovo Cimento. C_ , 1 suppl.1, (1976) 69.
  * (24) V.Denisov and S.Svertilov, _Vacuum nonlinear electrodynamics curvature of photon trajectories in pulsars and magnetars_ , _A &A_, 399, (2003) 39.
  * (25) J.Schinger, Wu-YangTsai and T.Erber, _Classical and quantum theory of synergic Syncrotron-Čherenkov radiation_ , _Ann. Phys._ , 90 , (1976) 303.
  * (26) T.Erber, D.White and Wu-YangTsai, _Experimental aspects of Syncrotron-Čherenkov radiation_ , _Ann. Phys._ , 102, (1976) 405.
  * (27) F.Karbstein and R.Shaisultanov, _Photon propagation in slowly varying inhomogeneous electromagnetic fields_ , _Phys. Rev. D_ , 91, (2015) 085027.
  * (28) A.E.Shabad and V.V.Usov, _Effective Lagrangian in nonlinear electrodynamics and its properties of causality and unitarity_ , _Phys. Rev. D_ , 83, (2011) 105006
  * (29) T.Erber and Wu-YangTsai, _Propogation of photon in homogeneous magnetic fields: Index of refraction_ , _Phys. Rev. D_ , 12, (1975) 1132.
  * (30) S.A.Olausen and V.M. Kaspi, _The McGill magnetar catalog_ , _The Astrophysical Journal Supplement Series_ , 212, (2014) 1.
  * (31) http://www.atnf.csiro.au/research/pulsar/psrcat/
  * (32) P.Soffitta, _XIPE: the X-ray imaging polarimetry explorer_ , _Experimental Astronomy_ , 36, (2013) 523.