format: text

text:
# Quantum gravity as a theory of quantized area bits fitting together

## Abstract

Non-Abelian Gauss law is interpreted in terms of area bits described in a local frame which fit together into closed surfaces and the Non-Abelian Stokes law in terms of length bits described in a local frame which fit together into closed loops. A new equation relating the area variables and the phase space variables ( or equivalently, angular momentum variables of the lattice Yang-Mills theory and phase space variables of the continuum theory) is obtained. Canonical quantization applied to these variables implies area quantization. A complete orthonormal basis of states satisfying the Gauss constraint is obtained.It has the interpretation of quantized area bits with undefined orientations and edges but fitting together into closed surfaces.

preprint: IMSc./2000/08/47

Quantum gravity has posed many intriguing puzzles which have not been fully unraveled yet. Einstein gravity is a non-renormalizable theory with a dimensionful coupling constant $`\mathrm{}_p=\sqrt{\mathrm{}Gc^3}`$, the Planck length. It is a popular conjecture that this length may play a fundamental role in Nature curing the divergences of quantum gravity and perhaps of matter fields too. Semi-classical analysis of the black holes shows that the area of the horizon has a mysterious significance and it has even been speculated that it may be quantized . In the loop formulation of gravity variables  for area, volume etc. have been proposed. It has been argued that they have non-trivial commutation relations on quantization, leading to quantized values and giving space time a non-commutative geometry.

In this paper we analyze the connection formalism  for gravity (in the real connection formulation ) in a different way. We are able to solve the Gauss law constraint explicitly in the quantum theory and obtain a complete orthonormal set of basis states. This clarifies the precise sense in which classical geometry is relevant to quantum gravity. It also brings out the precise sense in which area is quantized.

Consider Barbero’s formulation  of classical gravity with a real $`SU(2)`$ connection. The 3-dimensional space is taken to be an oriented, analytic manifold $`\mathrm{\Sigma }`$ described by arbitrary local coordinates $`\stackrel{}{x}=\{x^a\}`$.The configuration space consists of all smooth $`SU(2)`$ connections $`\stackrel{}{A}_a(\stackrel{}{x})=\{A_a^i(\stackrel{}{x})\}`$ satisfying appropriate boundary conditions.Here $`i=1,2,3`$ is the internal index (i.e. the label for the adjoint representation of $`SU(2)`$) and $`a=1,2,3`$ is the tangent space index. The momentum space is the space of smooth vector densities $`\stackrel{}{E}_a(\stackrel{}{x})=\{E^{ai}(\stackrel{}{x})\}`$ of weight one on $`\mathrm{\Sigma }`$. The Poisson brackets are

$$[A_a^i(\stackrel{}{x}),E^{bj}(\stackrel{}{y})]_{PB}=\delta _{ij}\delta _{ab}\delta (\stackrel{}{x}\stackrel{}{y})$$
(1)
Other brackets are zero.Only part of this phase space satisfying seven first class constraints is physically accessible and specify the theory. The first three are the non-Abelian Gauss law

$$D_a[\stackrel{}{A}]\stackrel{}{E}^a(\stackrel{}{x})=0.$$
(2)
This embodies the requirement that configurations related by a gauge transformation are equivalent. We will interpret this in a different way below, and it turns out to be of crucial relevance to gravity.

Interpretation of the non-Abelian Gauss law: When $`A_a^i=0`$ we have $`_S𝑑S_aE^{ai}=0`$ for any closed surface $`S`$. $`dS_a`$ is an area element in the local coordinates used. Associate a $`\mathrm{𝑝ℎ𝑦𝑠𝑖𝑐𝑎𝑙}`$ area $`dA^i=dS_aE^{ai}`$ to this area bit.The divergence-free condition means that the vectorial sum of the physical areas is zero, i.e., the surface bits fit together to form a closed surface when described in terms of their physical areas also.

To clarify the interpretation when $`A_a^i0`$, consider a lattice in $`\mathrm{\Sigma }`$ formed by an arbitrary grid of the coordinate axes $`x^a,a=1,2,3`$ close together. Consider a face spanned by the vertices

$$\stackrel{}{x},\stackrel{}{x}+\mathrm{\Delta }x^1\stackrel{}{n}_1,\stackrel{}{x}+\mathrm{\Delta }x^2\stackrel{}{n}_2,\stackrel{}{x}+\mathrm{\Delta }x^1\stackrel{}{n}_1+\mathrm{\Delta }x^2\stackrel{}{n}_2$$
(3)
where $`\stackrel{}{n}_1,\stackrel{}{n}_2`$ are unit vectors along the positive $`1`$\- and $`2`$\- axes and $`\mathrm{\Delta }x^a>0`$. Associate a physical area of magnitude $`\mathrm{\Delta }x^1\mathrm{\Delta }x^2|\stackrel{}{E}^3(\stackrel{}{x})|`$ to this face. However assign different orientations to it in the frames used in the two cells on either side of it. Specifically, assign vectorial areas

$$_\pm ^{3i}=\pm \mathrm{\Delta }x^1\mathrm{\Delta }x^2(exp(\pm \frac{1}{2}\mathrm{\Delta }x^3\stackrel{}{A}_3.\stackrel{}{T}))_j^iE^{3j}$$
(4)
in the cells to its positive and negative side respectively. Here $`(\stackrel{}{T}^i)_{jk}=ϵ_{ijk}`$ are the anti-symmetric matrices denoting the generators in the adjoint representation of $`SU(2)`$. An overall minus sign for $`\stackrel{}{}_{}`$ in (4) is a reflection of the fact that the outward normal (with respect to the appropriate cell ) is used to describe the area element. Thus the same area element is rotated by the orthogonal matrix $`𝒪`$,

$$𝒪=exp(\mathrm{\Delta }x^3\stackrel{}{A}_3.\stackrel{}{T})$$
(5)
in the two adjoining cells. Now demand that the physical areas of the six faces of any cell, when all are described by the frame used for it, vectorially add up to zero. Thus the faces close in the physical space also. This requirement gives the non-Abelian Gauss law, 2 in the limit when the grid is infinitely fine. It is natural to associate $`_{}^{3i}`$ with the dual coordinate $`\stackrel{}{x}^{}`$ pertaining to the cell and $`_+^{3i}`$ with the dual coordinate $`\stackrel{}{x}^{}+\mathrm{\Delta }x^3\stackrel{}{n}_3`$ pertaining to the neighboring cell in the third direction. Similar definitions for other directions are obvious. Now the Gauss law 2 can be rewritten as,

$$_{+x}(\stackrel{}{x}^{})+_x(\stackrel{}{x}^{})+_{+y}(\stackrel{}{x}^{})+_y(\stackrel{}{x}^{})+_{+z}(\stackrel{}{x}^{})+_z(\stackrel{}{x}^{})=0$$
(6)
This form allows an explicit solution of the Gauss law in the quantized case as will be shown later.

Our procedure has an obvious relation to lattice gauge theory . The rotation matrix (5) is associated with the link joining the vertices $`\stackrel{}{x}^{}`$ and $`\stackrel{}{x}^{}+\mathrm{\Delta }x^3\stackrel{}{n}_3`$ of the dual lattice and is the holonomy element associated with this infinitesimal line element. The electric field $`\stackrel{}{E}^3(\stackrel{}{x})`$ is formally associated with the midpoint of the associated dual link $`(\stackrel{}{x}^{},\stackrel{}{x}^{}+\mathrm{\Delta }x^3\stackrel{}{n_3})`$. The new variables $`\stackrel{}{}_\pm `$ are the electric fluxes on the face (3) of our grid, ‘parallel transported’ to the two vertices of the dual link, $`\mathrm{𝑣𝑖𝑎}`$ the connection $`\stackrel{}{A}`$, to the negative and positive side respectively. The factor $`1/2`$ in the exponent in (4) is because the parallel transport is by only half the length of the dual link. Thus all the variables in the Gauss law (6) transform under gauge transformations as vectors at the dual coordinate $`\stackrel{}{x}^{}`$ labeling the cell under consideration. This procedure of splitting the link into two makes the definition of the $`\stackrel{}{}_\pm `$ variables (4) symmetric, but would appear to be inessential to interpret the Gauss law. We could have adopted a convention of associating $`\stackrel{}{E}`$ to one of the associated dual vertices, and defined the other by a parallel transport through the entire dual link to the other dual vertex. However we will see that our definition is significant when we consider the Poisson bracket algebra of the variables 4. It turns out that they can be interpreted as the body- and the face-fixed angular momentum variables of lattice gauge theory .

The interpretation of a lattice regularization is not essential. As emphasized in  we may interpret it as a floating lattice in the continuum theory, used to obtain the later by a limiting procedure. The coordinate grid may be viewed as a special graph convenient to take the limit.

It is to be noted that if $`\stackrel{}{}_\pm `$ are specified for any surface bit, the corresponding $`\stackrel{}{A}`$ on the associated dual link is not completely determined but it will have only two degrees of freedom. This is as it should be, because the Gauss law (2) will not constrain the configuration space completely. It is also important to note that though the cell has the shape of a cuboid in the local coordinates used, it cannot be assigned a definite shape in the physical space.In the frame relevant to one cell, it can be given such a shape by computing the vectors corresponding to the edges from the $`\stackrel{}{}_\pm `$. However these lengths will not match in the neighbouring cubes.(This is discussed further below. ) This discussion is to emphasize that the area bits emerge as the crucial entities, and not other objects. To make the length bits the relevant objects we would need the non-Abelian Stokes law as is discussed below.

Interpretation of the non-Abelian Stokes law: the non-Abelian analogue of the Stokes law is,

$$ϵ_{abc}D_b[\stackrel{}{A}]\stackrel{}{e}_c(\stackrel{}{x})=0.$$
(7)
This is encountered as the condition that the dreibein $`e_a^i(\stackrel{}{x})`$ be torsion-free with respect to the connection 1-form $`A_a^i(\stackrel{}{x})`$ . This too has a geometric meaning. When $`A_a^i=0`$ we have $`_C𝑑x^ae_a^i(\stackrel{}{x})=0`$ for any closed loop $`C`$. $`dx^a`$ is the line element of the loop in the local coordinates used. Associate a $`\mathrm{𝑝ℎ𝑦𝑠𝑖𝑐𝑎𝑙}`$ length vector $`dl^i=dx^ae_a^i`$ to this length bit.The curl free condition means that the vectorial sum of the physical lengths is zero, i.e., the length bits fit together to form a closed loop when described in physical lengths also.

When $`A_a^i0`$, again associate a $`\mathrm{𝑝ℎ𝑦𝑠𝑖𝑐𝑎𝑙}`$ length of magnitude $`|\mathrm{\Delta }l_a|=\mathrm{\Delta }x^a|\stackrel{}{e}_a|`$ (where there is no sum over $`a`$) with each side $`\mathrm{\Delta }x^a`$ of the lattice. First we associate a different frame for each $`\mathrm{𝑓𝑎𝑐𝑒}`$ of the lattice. The above line element is assigned different orientations in each of the frames used for the six faces adjoining it. For the face (3), the edges $`(\stackrel{}{x},\stackrel{}{x}+\mathrm{\Delta }x^a\stackrel{}{n}_a)`$, (where $`a=1,2`$ and there is no sum over $`a`$) attached to the lowest coordinate $`\stackrel{}{x}`$ of this face are assigned directly the physical length vectors $`\mathrm{\Delta }_a\stackrel{}{l}=\mathrm{\Delta }x^a\stackrel{}{e}_a`$. However for the other two edges a parallel transport is required to compare with these vectors in the frame relevant to this face. For ex, the edge $`(\stackrel{}{x}+\mathrm{\Delta }x^1\stackrel{}{n}_1,\stackrel{}{x}+\mathrm{\Delta }x^1\stackrel{}{n}_1+\mathrm{\Delta }x^2\stackrel{}{n}_2)`$ is associated with the length vector $`\mathrm{\Delta }_2l^i=\mathrm{\Delta }x^2(exp(\mathrm{\Delta }x^1\stackrel{}{A}_1.\stackrel{}{T}))_j^ie_2^j`$ obtained by a parallel transport along the 1-axis to the point $`\stackrel{}{x}`$.

In this case it is also possible to choose orientations consistently for each length bit in a frame chosen for each $`\mathrm{𝑐𝑒𝑙𝑙}`$, so that the length bits ( and hence the areas computed from them in the frame) close. The asignment is as follows: for each of the length bits in the three faces attached to the ‘lowest’ coordinate of the cell, the orientations are as chosen above. For the length bits on the remaining three faces of the cell, ‘parallel transport’ all the orientations chosen on the face using the connection on the link connecting the face to the ‘lowest’ coordinate.

This gives the relationship to the Regge calculus. In this case the physical lengths assigned give a consistent simplicial complex. Specification of $`e`$’s in the frame of each face now completely specifies the $`A`$’s. This corresponds to the fact that the torsion-free condition can be used to solve for $`A`$ algebraically as a local expression in $`e`$’s. In this case, the faces are built out of specified length bits in the physical space also.

All this is in contrast to the case of the Gauss law. Consider the problem of constructing the areas out of length bits. Given the configuation of $`\stackrel{}{E}`$s satisfying the Gauss constraint (2) with a specified connection $`\stackrel{}{A}`$, we may construct a dreibein $`\stackrel{}{e}`$ by $`\stackrel{}{E}^a=ϵ_{abc}\stackrel{}{e}_b\times \stackrel{}{e}_c`$. We may also construct the corresponding connection $`\stackrel{}{A}[\stackrel{}{E}]`$ which is torsion-free with respect to which this dreibein. Generically, this will not coincide with the specified connection $`\stackrel{}{A}`$. In this case we get frames in each cell, such that the areas of the faces computed in the frame coincide precisely with the configuation of $`\stackrel{}{E}`$s. But these frames will be related to $`\stackrel{}{A}[\stackrel{}{E}]`$ and not to the specified connection $`\stackrel{}{A}`$. Therefore the area vectors computed will not have the required orientations, though they have the correct magnitude.

The area variables we have defined have a very significant Poisson brackets algebra among each other. With our discretiztion, $`\delta (\stackrel{}{x})(\mathrm{\Delta }x^1\mathrm{\Delta }x^2\mathrm{\Delta }x^3)^1`$. Therefore, $`[\mathrm{\Delta }x^3A_3^i(\stackrel{}{x}),\mathrm{\Delta }x^1\mathrm{\Delta }x^2E^{3j}(\stackrel{}{y})]_{PB}=\delta _{ij}\delta _{\stackrel{}{x},\stackrel{}{y}}`$. Thus the variables in the bracket form canonically conjugate pairs.

Poisson bracket algebra: We now compute the Poisson bracket,

$`[_+^{3a}(\stackrel{}{x}^{}),_\pm ^{3b}(\stackrel{}{y}^{})]_{PB}=\pm \mathrm{\Delta }x^1\mathrm{\Delta }x^2E^{3k}(ϵ_{lik}(exp(\pm {\displaystyle \frac{1}{2}}\mathrm{\Delta }x^3\stackrel{}{A}_3.\stackrel{}{T}))_l^j\pm ϵ_{ljk}(exp({\displaystyle \frac{1}{2}}\mathrm{\Delta }x^3\stackrel{}{A}_3.\stackrel{}{T}))_l^i)`$
Expanding the exponent in powers of $`\mathrm{\Delta }x^3`$ we see that $`(_+^{3a},a=1,2,3)`$ has the algebra of angular momentum among each other for a given cell whereas they commute with $`(_{}^{3b},b=1,2,3)`$ even in the same cell, to order $`O(\mathrm{\Delta }x^3)`$. Thus in the limit of an infinitely fine grid, $`\stackrel{}{}_\pm ^{ai}(\stackrel{}{x}^{})`$ for each $`\pm `$, $`\stackrel{}{x}^{}`$ and $`\stackrel{}{a}`$ have the PB of angular momentum, whereas they have zero PBs with all other variables.This identifies them with the body- and space fixed angular momenta of $`SO(3)`$ lattice gauge theory, which generate an $`SO(3)`$ rotation of the link variable on the left and right respectively . In fact $`\stackrel{}{}_\pm (\stackrel{}{x}^{})`$ respectively generate an $`SO(3)`$ rotation of all Wilson lines starting from the point and having the tangent along the positive and negative 3-direction respectively at this point. Our procedure of splitting the dual link makes this operation precise. If the Wilson line begins with the positive (negative) half of the dual link, it transforms non-trivially under $`\stackrel{}{}_\pm (\stackrel{}{x})`$ respectively. In effect $`\stackrel{}{E}(\stackrel{}{x})`$ in $`\stackrel{}{}_\pm (\stackrel{}{x}^{})`$ (see (4)) removes the orthogonal matrix associated with this half link and it is replaced back by the same orthogonal matrix contained in $`\stackrel{}{}_\pm `$.

Quantization: With this Poisson bracket structure, prescription for quantization is evident. Each $`_+^i`$ or $`_{}^i`$ becomes the three components of independent angular momentum operators. There is however the constraint that $`\stackrel{}{}_+^2(\stackrel{}{x}^{}+\mathrm{\Delta }x^3\stackrel{}{n}_3)=\stackrel{}{}_{}^2(\stackrel{}{x}^{})`$ as the two vectors differ from each other only by an orthogonal rotation. This condition is easy to impose. Ignoring this condition we have the standard basis $`|j_+m_+>|j_{}m_{}>`$ where the total angular momentum and the ‘third’ component are simultaneously diagonalised.Now the constraint is simply solved by choosing a subset of this basis with $`j_+=j_{}(=j`$ say). Thus the label for the basis is $`|j,m_+,m_{}>`$. The reason for this is simple. The link variable may be interpreted  as the configuration of a rigid body described by the rotation $`𝒪`$ from the space- to the body- fixed axes. The corresponding wave function $`\mathrm{\Psi }(𝒪)`$ can be expanded in the basis $`D_{mn}^j`$ of the rotation matrices. $`j`$ labels the total angular momentum which is the same in the two frames. $`m,n`$ label the eigenvalues of the ‘third component’ of the body- and space- fixed angular momenta which can be simultaneously diagonalised.

Gauss law in form (6) is simply the constraint that the six angular momenta associated with the six faces of a cell add up to zero. In quantum theory this can be solved by going to an appropriate coupled basis. This enables us to construct a complete orthonormal basis which satisfies the Gauss law. This procedure is a simple extension of the one used in  for lattice gauge theory in two space dimensions, where the basis was interpreted in terms of triangulations. The main difference is that the $``$s are now interpreted as area bits instead of the length bits.

We choose a convention for the sequence in which the angular momenta are added. Any other convention leads to an equivalent basis, though some of them may not lead to the geometric interpretation given below.For each 3-cell, we add the angular momentum pair-wise in the following sequence repeatedly:

$`|j_{+x}m_{+x}>|j_{+y}m_{+y}>|j_{+z}m_{+z}>|j_xm_x>|j_ym_y>|j_zm_z>`$
$`|j_{+x}j_yj_{+xy}m_{+xy}>|j_{+y}j_zj_{+yz}m_{+yz}>|j_{+z}j_xj_{+zx}m_{+zx}>`$
$`|j_{+x}j_yj_{+y}j_zj_{+xy}j_{+yz}j_{+xy+yz}m_{+xy+yz}>|j_{+z}j_xj_{+zx}m_{+zx}>`$
$`|j_{+x}j_yj_{+y}j_zj_{+z}j_xj_{+xy}j_{+yz}j_{+zx}j_{+xy+yz}j_{+xy+yz+zx}m_{+xy+yz+zx}>`$
In the basis written last it is easy to implement the Gauss constraint. It simply gives $`j_{+xy+yz+zx}=m_{+xy+zx+yz}=0`$. In this case we have necessarily, $`j_{+xy+yz}=j_{+zx}`$. Thus the complete orthonormal basis satisfying the Gauss law (2) is simply labeled by the integer labels, $`|\{j_{+x}(\stackrel{}{x}^{})j_{+y}(\stackrel{}{x}^{})j_{+z}(\stackrel{}{x}^{})j_x(\stackrel{}{x}^{})j_y(\stackrel{}{x}^{})j_z(\stackrel{}{x}^{})j_{+xy}(\stackrel{}{x}^{})j_{+yz}(\stackrel{}{x}^{})j_{+zx}(\stackrel{}{x}^{})\}>`$. It is further understood that $`j_a(\stackrel{}{x}^{})=j_{+a}(\stackrel{}{x}^{}+\mathrm{\Delta }x^a\stackrel{}{n}_a)`$ for each $`a=1,2,3`$.

The interpretation is follows: $`j`$s are the magnitudes of the areas associated with the faces of our cells.They are quantized. Further there are constraints. For example, $`(j_{+x},j_y,j_{+xy})`$ pertaining to one cell satisfy the triangle rule for the addition of angular momenta. This is interpreted to mean that with suitable choice of orientations, the area elements associated with $`j_{+x}`$ and $`j_y`$ add up vectorially to the area element associated with $`j_{+xy}`$. Again $`(j_{+xy},j_{+yz},j_{+zx})`$ satisfy the triangle rule. This is again interpreted to mean that corresponding areas close on appropriate choice of their orientations. The orientations assigned for (say) $`j_{+xy}`$ in these two cases need not match.

It is also possible to interpret the Gauss constraint in the present case in terms of length bits . These form the sides of truncated octohedrons which are stacked end-to end.

Diffeomorphism invariance: We now consider the implications of the next three constraints . They imply that the theory is invariant under diffeomorphisms of 3-dimensional space. Gauss law (2) is form invariant under an arbitrary change of the local coordinates $`\stackrel{}{x}`$ provided $`\stackrel{}{A}`$ and $`\stackrel{}{E}`$ are appropriately transformed. Consider the coordinate grid with the new coordinates. The new $`\stackrel{}{}`$s are related to the earlier ones (in 4) as follows: consider the electric flux on the faces of the new grid instead of the old one. These are related to the earlier ones by a linear transformation. Also, the orthogonal matrix (5) to be used now pertain to the holonomy elements along the new dual links. The appropriate Lie algebra element (i.e. the exponent in (4)) is again a linear combination of the previous ones.The three space diffeomorphism constraints imply that these two configurations are physically equivalent. The way to handle this and its implications will be considered elsewhere.

We now briefly comment on the continuum limit. This corresponds to the limit of an infinite number of cells. It is to be expected that the discreteness of area does not get washed out in this limit at least for certain macroscopic situations such as the black holes.

Our considerations suggest that quantum gravity is a theory of area bits with quantized areas. The magnitudes of these areas can be simultaneously observed, but the orientations are indefinite, due to the quantum mechanical uncertainty. However, the magnitudes are such that it is possible to fit them together into closed surfaces with some choice of orientations. This brings out the specific way in which classical geometry is relevant to quantum gravity. As was noted before, though these area bits were associated with a rectangular bit in the coordinates chosen to describe the manifold, they cannot be associated with specified edges in physical space. As discussed, such a specification would be a consequence of the torsion-free condition instead of the Gauss constraint. This may be again interpreted as a consequence of the quantum mechanical uncertainty which disallows a simultaneous specification of the conjugate variables.

The geometry relevant to the basis of states obtained here will be addressed in more detail elsewhere .

We thank Professor Ramesh Anishetty for discussions on issues and techniques in gravity and gauge theories.


text:
# Quantum authentication using entangled state

## I Introduction

Quantum cryptography is a field which combines quantum theory with information theory. The goal of this field is to use the laws of physics to provide secure information exchange, in contrast to classical methods based on (unproven) complexity assumption. Since the publication of BB84 protocol , quantum key distribution (QKD)  has developed into a well understood application of quantum mechanics to cryptography. In particular, QKD protocols became especially important due to technological advances which allow their implementation in the laboratory . To guarantee the security of the quantum key in practical applications, quantum key authentication is important as well as QKD. Moreover, large quantity of communication tasks in modern society need more reliable authentication systems. However, up to now, the security of practically used authentication systems is based on the computational difficulty, i.e., they rely on limited advancement of computer power, technologies, and mathematical algorithms in the foreseeable future.

Recently, several quantum authentication schemes  have been proposed. The scheme proposed by Dušek et al.  is a combination of classical identification procedure and quantum key distribution. Zeng’s  protocol is using EPR pairs as the first authentication keys, and using classical keys distributed in quantum key distribution procedure after then. Moreover, two very interesting authentication protocols using entanglement and catalysis  were proposed by Barnum  and Jensen et al.  respectively. In their protocol, the two parties in communication have previously shared catalyst (a particular pair of entangled particles). The verifier sends the challenge which is a half part of some entangled state to the identifier, then they can transfer this entangled state to a special state (the state of the catalyst) deterministically by local quantum operation and classical communication (LQCC) with the catalyst. However, this task can not be accomplished by LQCC without the catalyst. So the verifier can authenticate the identifier by measuring the state of the challenge after the identifier sends it back.

On the other hand, entanglement of multiparticle system is a important feature of quantum mechanics. In addition to their central role in discussion of nonlocal quantum correlations, they form the basis of quantum information such as quantum teleportation , quantum key distribution, quantum dense coding  et al. EPR pairs are used as communication channels in protocols mentioned above. In this paper, we present a authentication protocol with EPR state as the key (encoder and decoder). The authentication is accomplished with local controlled-NOT (C-NOT) operations and unitary rotations.

The paper is structured as follows. In Section II, we give the framework of our authentication protocol, and security of this protocol is analyzed in Section III. In Section IV, we consider the nonideal situation and discuss the robustness of this protocol. Section V concludes the paper.

## II Quantum authentication scheme

The general task of authentication is verifying the identification of each other of two parties (Alice and Bob) in communication, using quantum and classical channel. The protocols are such that if Alice and Bob can successfully complete one, Alice is convinced that Bob (or someone who has stolen his identification token) is on the other end of the quantum communication channel. The classical analogue of this can be done by having Bob to reveal, over a classical channel, a secret which Alice and Bob had previously securely shared. The quantum protocol presented here uses shared entangled states as the counterpart of shared secret key. There is protection, via the no-cloning theorem, against copying of the authentication token. And this protocol may provide reusable authentication tokens.

Alice and Bob have previously shared $`2K`$ pairs entangled states in

$$|\mathrm{\Phi }^+=\frac{1}{\sqrt{2}}\left(|00+|11\right).$$
(1)
An authentication round consists of the following steps. First, for example, Alice acts as identifier and Bob acts as verifier.

When the authentication begin, the two parties rotate their particles’ state by $`\theta `$ respectively. The rotation can be described as

$$R(\theta )=\left(\begin{array}{cc}\mathrm{cos}\theta & \mathrm{sin}\theta \\ \mathrm{sin}\theta & \mathrm{cos}\theta \end{array}\right).$$
(2)
The state $`|\mathrm{\Phi }^+`$ does not change under bilateral operation of $`R(\theta )`$. The purpose of this operation is to prevent the eavesdropper’s impersonation. (The detailed interpretation will be given in Section III.) Then Bob prepares $`K^{}`$ ($`K^{}K`$) particles $`\gamma _B^i`$ in arbitrary pure state

$$|\psi _i=a_i|0+b_i|1,\text{ }i=1,3,5,\mathrm{},2K^{}1,$$
(3)
where $`\left|a_i\right|^2+\left|b_i\right|^2=1`$, $`a_i`$ and $`b_i`$ are two complex numbers which are selected randomly. The state $`|\psi _i`$ is only known by Bob himself. Bob sends the challenge particle $`\gamma _B^i`$ to Alice in order of $`i`$ (where $`i`$ is odd). Thus Alice uses the corresponding particle $`\beta _A^i`$ of the entangled pairs and the particle $`\gamma _B^i`$ to do a C-NOT operation ($`\beta _A^i`$ is controller and $`\gamma _B^i`$ is target) and the three particles’ state will be

$$|\mathrm{\Psi }_i=\frac{1}{\sqrt{2}}\left(a_i|000+b_i|001+a_i|111+b_i|110\right),$$
(4)
then she sends back $`\gamma _B^i`$ to Bob. Bob uses his corresponding particles $`\beta _B^i`$ (which entangled with $`\beta _A^i`$) to do a C-NOT operation on $`\gamma _B^i`$ again. Now the state of the key particles and the challenge $`\gamma _B^i`$ is the same as it at the first

$$|\mathrm{\Psi }_i^{}=\frac{1}{\sqrt{2}}\left(|00+|11\right)\left(a_i|0+b_i|1\right)=|\mathrm{\Phi }^+|\psi _i.$$
(5)
Bob measures $`\gamma _B^i`$ in basis $`|\psi _i`$ and $`|\psi _i^{}`$ (state orthogonal to $`|\psi _i`$). If $`\gamma _B^i`$ is in state $`|\psi _i`$, it passes the test; otherwise, it fails and Bob aborts the protocol. Then Alice becomes the verifier, she prepares $`K^{}`$ challenge $`\gamma _A^i`$ and uses her particle $`\beta _A^i`$ ($`i`$ is even) to do the same steps.

The authentication fails if any of the projective measurements in the previous step fails, or if Alice or Bob receive more than $`K^{}`$ requests to send back challenge particles.

If the authentication round succeeds, Alice and Bob retain all $`2K`$ pairs of entangled states and can reuse them in the next time. However, the security of the entangled states used later is a little less than the original one. If the authentication fails, the parties discard all particles used till that point. In this case, Alice and Bob have to start again with new keys (EPR pairs).

## III Security analysis

We now discuss the security of this protocol. First, Eve may impersonate Alice when Alice is not present. When Bob sends out a challenge $`\gamma _B^i`$, Eve intercepts it which she can manipulate using unitary transformation or measurement. Since Eve has not shared the key with Bob, she can not entangle this challenge with Bob’s key particles. Suppose the state Eve send back is

$`\rho _i`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{2}{}}}p_{ik}|\psi _{ik}^{}\psi _{ik}^{}|,`$ ()
$`|\psi _{ik}^{}`$ $`=`$ $`a_{ik}^{}|0+b_{ik}^{}|1.`$ (6)
where $`p_{i1}+p_{i2}=1`$ and $`\left|a_{ik}^{}\right|^2+\left|b_{ik}^{}\right|^2=1`$. After Bob’s C-NOT operation, the fidelity  of the state and the test state $`|\mathrm{\Psi }_i`$ is

$$F=\frac{1}{2}+\underset{k=1}{\overset{2}{}}p_{ik}\left[Re\left(a_i^{}b_ia_{ik}^{}b_{ik}^{}\right)+Re\left(a_i^{}b_ia_{ik}^{}b_{ik}^{}\right)\right],$$
(7)
where $`Re(x)`$ is the real part of complex number $`x`$. Since $`a_i`$ and $`b_i`$ are selected randomly, the average value of $`F`$ is $`\frac{1}{2}`$. That is to say Eve has the probability only $`\frac{1}{2}`$ to pass the test in one time on the average. The probability that Eve is not detected by Bob is $`\left(\frac{1}{2}\right)^K^{}`$ on the average, which can be made arbitrarily small by choosing $`K^{}`$ large enough.

Second, Eve may use the method of denial of service. In this type of attack, Eve deliberately causes the authentication round to fail, and hence causes one party to discard all key particles. Although this protocol is particularly vulnerable to this kind of attack, this is not an essential weakness, since an attacker who controls both quantum and classical communication can always prevent successful authentication between the legitimate parties.

We now look at stronger attacks in which Eve tries to obtain key material which she could then use, e.g., in a later impersonation attack. Eve’s goal is to share pairs of particles in the entangled state $`|\mathrm{\Phi }^+`$ with Alice and/or Bob. For instance, if she succeeds in obtaining a large amount of key material with Bob, she will be able to authenticate herself to Bob without Alice being present. However, if Eve’s presence is detected in a single measurement, all the previously obtained key material she shares with the verifier who performed that measurement will be worthless.

Eve has three choices to attack: (I) She can intercept the challenge Bob sends to Alice, and make her own particle interact with this challenge, then send the challenge or her own particles to Alice. (II) She can pass the challenge to Alice but intercept it when Alice sends it back to Bob, and make her own particles interact with this challenge. (III) She can combine the strategies (I) and (II). We will analyze Eve’s three strategies respectively.

(I) Since the challenge Bob sends to Alice is not entangled with the key particles, Eve can not share key material with Bob or Alice. If she changes the state of the challenge, she will be detected with probability $`1F`$, and the form of $`F`$ is the same as Eq. (6).

(II) In this case, Eve may make his own particles entangle with Alice and Bob’s key pairs. However, the challenge state is selected randomly, Eve can not make her own particle be maximally entangled with Bob and/or Alice’s particle deterministically, and the probability she will be detected is $`1F`$ too.

However, there exist other more powerful strategy.

(III) Eve intercept the challenge Bob send to Alice, and send her own particle in state $`|0`$ or $`|1`$ to Alice. After Alice’s operation, this particle will be entangled with Alice and Bob’s key particles in GHZ state, and Eve can use the key to complete authentication with Alice or Bob as efficient as the key shared by Alice and Bob. We can defend this attack by an additional step at first. Before Bob sends the challenge to Alice, they do a bilateral rotation $`R(\theta _i)=\left(\begin{array}{cc}\mathrm{cos}\theta _i& \mathrm{sin}\theta _i\\ \mathrm{sin}\theta _i& \mathrm{cos}\theta _i\end{array}\right)`$. The state of the maximally entangled two particles will be unchanged when both Alice and Bob rotate the $`i`$th particle by $`\theta _i`$. However if Eve has entangled her particle with Alice and Bob’s particles in state $`|\mathrm{\Phi }_{ABE}=\frac{1}{\sqrt{2}}\left(|000+|111\right)`$. In the second authentication process between Alice and Bob, the state will be changed to

$`|\mathrm{\Phi }_i_{ABE}`$ $`=`$ $`\mathrm{cos}^2\theta _i{\displaystyle \frac{1}{\sqrt{2}}}\left(|000+|111\right)+\mathrm{sin}^2\theta _i{\displaystyle \frac{1}{\sqrt{2}}}\left(|110+|001\right)`$ (7)
$`+\mathrm{sin}\theta _i\mathrm{cos}\theta _i{\displaystyle \frac{1}{\sqrt{2}}}\left(|011|100\right)+\mathrm{sin}\theta _i\mathrm{cos}\theta _i{\displaystyle \frac{1}{\sqrt{2}}}\left(|101|010\right)`$
under the bilateral rotation. When Bob or Alice measures the challenge, the average fidelity will be $`1\frac{1}{2}\mathrm{sin}^2\theta _i`$, and the average probability that they find Eve’s attack will be $`\frac{1}{2}\mathrm{sin}^2\theta _i`$. If $`\theta _i`$ is selected randomly and known only by Alice and Bob (it can be previously shared and reusable), the average fidelity of the challenge when Bob test it is $`\frac{3}{4}`$ and the average probability that Eve will be detected is $`\frac{1}{4}`$. If Eve tries to impersonate Alice or Bob with the entangled state $`|\mathrm{\Phi }_{ABE}`$, it is easy to verify that the average probability for Eve to succeed is less than $`1\frac{1}{2}\mathrm{sin}^2\theta _i`$ (even Eve knows exactly the value of $`\theta _i`$).

Now we consider to use a fixed rotation angle $`\theta `$ for all authentication rounds.

There exist two special angle, that Eve has two corresponding strategies to impersonate or to obtain the authentication key as her particle has been already entangled with Alice and Bob’s key particles in GHZ state. When $`\theta =0`$ (or $`\frac{\pi }{2}`$, $`\pi `$ et al.), Eve can use her particles impersonate Alice without risk that Bob will detect the impersonation. Another special angle is $`\theta =\frac{\pi }{4}`$ (or $`\frac{3\pi }{4}`$, $`\frac{5\pi }{4}`$ et al.). In this case, after Alice and Bob rotate their particle by $`\frac{\pi }{4}`$, Eve can rotate her particle by $`\frac{\pi }{4}`$ too. Then Eve intercepts the challenge Bob sends to Alice and sends her own particle to Alice, after Alice does a C-NOT operation on Eve’s particle, it will be entangled with Bob’s particle in EPR state and Alice’s particle has no quantum correlation with Bob and Eve’s particles. The process can be described as

$`|\mathrm{\Phi }_{ABE}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|000+|111\right]_{ABE}`$ (9)
$`\stackrel{R}{}{\displaystyle \frac{1}{4}}\left[\left(|0|1\right)\left(|0|1\right)\left(|0|1\right)+\left(|0+|1\right)\left(|0+|1\right)\left(|0+|1\right)\right]_{ABE}`$
$`\stackrel{C_{AE}}{}{\displaystyle \frac{1}{2}}\left(|0+|1\right)_A\left(|00+|11\right)_{BE},`$
where $`R=R_A\left(\pi /4\right)R_B\left(\pi /4\right)R_E\left(\pi /4\right)`$ and $`C_{AE}`$ is a C-NOT operation (A controls E). After Alice sent Eve’s particle back, Eve uses it to do a C-NOT operation on Bob’s challenge, then sends the challenge back to Bob for verification. Up to now, Eve has obtained the authentication key without disturbing the authentication process between Alice and Bob.

For an arbitrary $`\theta `$ (we assume that Eve knows this angle), if Eve tries to impersonate, in the case of $`\mathrm{sin}^2\theta <\mathrm{cos}^2\theta `$, the optimal strategy for Eve is to do a C-NOT operation on the challenge directly and send back it. The probability she will succeed is $`1\frac{1}{2}\mathrm{sin}^2\theta `$. In the case of $`\mathrm{sin}^2\theta >\mathrm{cos}^2\theta `$, Eve should do a NOT operation on her particle first, then do a C-NOT operation on the challenge before she sends it back. The probability she will succeed is $`1\frac{1}{2}\mathrm{cos}^2\theta `$. Sum the two cases, the probability Eve will succeed is

$$P_1=max\{1\frac{1}{2}\mathrm{sin}^2\theta ,1\frac{1}{2}\mathrm{cos}^2\theta \}.$$
(10)
If Eve tries to get the key, she does a rotation on her particle, then sends her particle to Alice to do a C-NOT operation on the particle and Eve can do a another proper rotation on her particle when Alice sends it back

$`|\mathrm{\Phi }_{ABE}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|000+|111\right]_{ABE}\stackrel{R_{ABE}}{}|\mathrm{\Phi }(\theta ,\varphi _1)_{ABE}`$ (10)
$`\stackrel{C_{AE}}{}|\mathrm{\Phi }^{}(\theta ,\varphi _1)_{ABE}\stackrel{R_E}{}|\mathrm{\Phi }(\theta ,\varphi _1,\varphi _2)_{ABE},`$
where $`R_{ABE}=R_A\left(\theta \right)R_B\left(\theta \right)R_E\left(\varphi _1\right)`$ and $`R_E`$ is operation $`R_E\left(\varphi _2\right)`$. Then the probability $`P_2`$ that Eve will succeed to obtain the key is

$`P_2`$ $`=`$ $`\underset{\varphi _1,\varphi _2}{\mathrm{max}}\left\{F(|\mathrm{\Phi }^+\mathrm{\Phi }^+|,tr_A\left(|\mathrm{\Phi }(\theta ,\varphi _1,\varphi _2)_{ABE}\mathrm{\Phi }(\theta ,\varphi _1,\varphi _2)|_{ABE}\right))\right\}`$ ()
$`=`$ $`{\displaystyle \frac{1}{2}}\left(\left|\mathrm{cos}\theta \right|+\left|\mathrm{sin}\theta \right|\right)^2.`$ (11)
Now we consider the optimal angle Alice and Bob should select if they use a fixed angle $`\theta `$ for all EPR pairs. Using the principle $`P_1=P_2`$, we can get

$$\left|\mathrm{cos}\theta \right|=\frac{2}{\sqrt{5}}\text{ },\frac{1}{\sqrt{5}};\text{ }P_1=P_2=\frac{9}{10}.$$
(12)
From above we can know, if the authentication process use a fixed rotation angle $`\theta `$, the successful probability of Eve’s eavesdropping will be increased.

## IV Robustness of the protocol

Up to this point, our discussion has assumed that the initial state is ideal maximal entangled state $`|\mathrm{\Phi }^+`$. Suppose, however, that this state is corrupted a little after reused for many times, Alice and Bob have a state described by density matrix

$$\rho =\left(1ϵ\right)|\mathrm{\Phi }^+\mathrm{\Phi }^+|+ϵ\rho _1,$$
(12)
where $`ϵ`$ is a parameter of the deviation of $`\rho `$ from $`|\mathrm{\Phi }^+\mathrm{\Phi }^+|`$ and $`\rho _1`$ is an arbitrary state. Our results are most easily presented using the trace distance, a metric on Hermitian operators defined by $`T(A,B)=Tr\left(\left|AB\right|\right)`$ , where $`\left|X\right|`$ denotes the positive square root of the Hermitian matrix $`X^2`$. From above, we can get that $`T(|\mathrm{\Phi }^+\mathrm{\Phi }^+|,\rho )2\sqrt{ϵ}`$.

Ruskai  has shown that the trace distance contracts under physical processes. If all operations are exact in the next authentication process, since the state $`|\mathrm{\Phi }^+`$ will be unchanged in the process, the density matrix $`\rho `$ will be transformed to

$$\rho ^{}=\left(1ϵ\right)|\mathrm{\Phi }^+\mathrm{\Phi }^+|+ϵ\rho _1^{},$$
(13)
and $`T(|\mathrm{\Phi }^+\mathrm{\Phi }^+|,\rho ^{})2\sqrt{ϵ}`$.

The fidelities $`F(|\mathrm{\Phi }^+\mathrm{\Phi }^+|,\rho )`$ and $`F(|\mathrm{\Phi }^+\mathrm{\Phi }^+|,\rho ^{})`$ are both no less than $`1ϵ`$, so the probability that the authentication fails is no more than $`ϵ`$. In conclusion, we can say that this protocol is robust.

## V Conclusion

A scheme of quantum authentication using entangled state is presented. Two parties share EPR pairs previously as the authentication key which servers as encoder and decoder. The authentication is accomplished with local controlled-NOT operations and unitary rotations. This protocol appears to be secure even in the presence of an eavesdropper who has complete control over both classical and quantum communication at all times. Our protocol is not rely on classical cryptography, and needs not communication of classical information of measurement results or details of operation method except the indices of the particle Bob sends to Alice. Compared to the scheme using catalysis, this protocol uses state only two dimensions, instead of five dimensions. And compare to other authentication schemes, this protocol has the same advantage that the keys are reusable as scheme using catalysis. At the end, we expect that this method can be applied in quantum key distribution (QKD) .

Acknowledgments

We would like to thank Zheng-Wei Zhou for the very helpful discussions. This work was supported by the National Natural Science Foundation of China.


text:
# Model selection for inverse problems: Best choice of basis function and model order selection Presented at the 19th Int. worskhop on Bayesian and Maximum Entropy methods (MaxEnt 1999), Aug. 2-6, 1999, Boise, Idaho, USA

## 1. Introduction

In a very general linear inverse problem, the relation between the data $`𝒚=[y_1,\mathrm{},y_m]^t`$ and the unknown function $`f(.)`$ is

$$y_i=h_i(𝒓)f(𝒓)\text{d}𝒓,i=1,\mathrm{},m,$$
(1)
where $`h_i(𝒓)`$ is the system response for the data $`y_i`$. We assume here that $`h_i(𝒓)`$ are known perfectly. The first step for any numerical processing is the choice of a basis function $`b_j(𝒓)`$ and an order $`k`$, in such a way to be able to write

$$f(𝒓)=\underset{j=1}{\overset{k}{}}x_jb_j(𝒓).$$
(2)
This leads to

$$𝒚=𝑨𝒙+\mathit{ϵ}$$
(3)
with $`𝒚=[y_1,\mathrm{},y_m]^t`$, $`𝒙=[x_1,\mathrm{},x_k]^t`$ and

$$A_{i,j}=h_i(𝒓)b_j(𝒓)\text{d}𝒓,i=1,\mathrm{},m,j=1,\mathrm{},k$$
(4)
where $`\mathit{ϵ}=[ϵ_1,\mathrm{},ϵ_m]^t`$ represents the errors (both the measurement noise and the modeling and the approximation related to the numerical computation of matrix elements $`A_{i,j}`$). Even, when the choice of the basis functions $`b_i(𝒓)`$ and the model order $`k`$ is fixed, obtaining a good estimate for $`𝒙`$ needs other assumptions about the noise $`\mathit{ϵ}`$ and about $`𝒙`$ itself. The Bayesian approach provides a coherent and complete framework to handle the random nature of $`\mathit{ϵ}`$ and the a priori incomplete knowledge of $`𝒙`$.

The first step in a Bayesian approach is to assign the prior probability laws $`p(𝒚|𝒙,\mathit{\varphi },k,l)=p_ϵ(𝒚𝑨𝒙|\mathit{\varphi },k,l)`$, $`p(𝒙|𝝍,k,l)`$, $`p(\mathit{\varphi }|k,l)`$ and $`p(𝝍|k,l)`$, where $`p_ϵ(𝒚𝑨𝒙|\varphi ,k,l)`$ is the probability law of the noise, and $`(\mathit{\varphi },𝝍)`$ the hyperparameters of the problem. Note that $`𝒙`$ represents the unknown parameters, $`k=\text{dim}(𝒙)`$ is the order of the model, $`m=\text{dim}(𝒚)`$ is the number of the data and $`l`$ is an index to a particular choice of basis functions. Note that the elements of the matrix $`𝑨`$ depend on the choice of the basis functions. However, to simplify the notations, we do not write this dependance explicitely. We assume that we have to select one set $`l`$ of basis functions among a finite set (indexed by $`[1:l_{max}]`$) of them. Thus, for a given $`l[1,l_{max}]`$ and a given model order $`k[1,k_{max}]`$, and using the mentionned prior laws, we define the joint probability law

$$p(𝒚,𝒙,\mathit{\varphi },𝝍|k,l)=p(𝒚|𝒙,\mathit{\varphi },k,l)p(𝒙|𝝍,k,l)p(\mathit{\varphi }|k,l)p(𝝍|k,l).$$
(5)
From this probability law, we obtain, either by integration or by summation, any marginal law, and any a posteriori probability law using the Bayes rule.

What we propose in this paper is to consider the following problems:

* Parameter estimation:

$$\widehat{𝒙}=\underset{𝒙}{\text{arg}\mathrm{max}}\left\{p(𝒙|𝒚,\widehat{\mathit{\varphi }},\widehat{𝝍},\widehat{k},\widehat{l})\right\}$$
(6)
where

$$p(𝒙|𝒚,\mathit{\varphi },𝝍,k,l)=p(𝒚,𝒙|\mathit{\varphi },𝝍,k,l)/p(𝒚|\mathit{\varphi },𝝍,k,l),$$
(7)
$$p(𝒚,𝒙|\mathit{\varphi },𝝍,k,l)=p(𝒚|𝒙,\mathit{\varphi },k,l)p(𝒙|𝝍,k,l)$$
(8)
and

$$p(𝒚|\mathit{\varphi },𝝍,k,l)=p(𝒚,𝒙|\mathit{\varphi },𝝍,k,l)\text{d}𝒙.$$
(9)
* Hyperparameter estimation:

$$(\widehat{\mathit{\varphi }},\widehat{𝝍})=\underset{(\mathit{\varphi },𝝍)}{\text{arg}\mathrm{max}}\left\{p(\mathit{\varphi },𝝍|𝒚,\widehat{k},\widehat{l})\right\}$$
(10)
where

$$p(\mathit{\varphi },𝝍|𝒚,k,l)=p(𝒚,\mathit{\varphi },𝝍|k,l)/p(𝒚|k,l)$$
(11)
and

$$p(𝒚|k,l)=p(𝒚,\mathit{\varphi },𝝍|k,l)\text{d}\mathit{\varphi }\text{d}𝝍.$$
(12)
* Model order selection:

$$\widehat{k}=\underset{k}{\text{arg}\mathrm{max}}\left\{p(k|𝒚,\widehat{l})\right\}$$
(13)
where

$$p(k|𝒚,l)=p(𝒚|k,l)p(k)/p(𝒚|l)$$
(14)
and

$$p(𝒚|l)=\underset{k=1}{\overset{k_{max}}{}}p(𝒚|k,l)p(k).$$
(15)
* Basis function selection:

$$\widehat{l}=\underset{l}{\text{arg}\mathrm{max}}\left\{p(l|𝒚)\right\}$$
(16)
where

$$p(l|𝒚)=p(𝒚|l)p(l)/p(𝒚)$$
(17)
and

$$p(𝒚)=\underset{l=1}{\overset{l_{max}}{}}p(𝒚|l)p(l).$$
(18)
* Joint parameter, hyperparameter, model order and basis function estimation:

$$(\widehat{𝒙},\widehat{\mathit{\varphi }},\widehat{𝝍},\widehat{k},\widehat{l})=\underset{(𝒙,\mathit{\varphi },𝝍,k,l)}{\text{arg}\mathrm{max}}\left\{p(𝒚,𝒙,\mathit{\varphi },𝝍|k,l)p(k)p(l)\right\}.$$
(19)

As it can be easily seen, the first problem is, in general, a well posed problem and the solution can be computed, either analytically or numerically. The others (excepted the last) need integrations. These integrals can be done analytically only in the case of Gaussian laws. In other cases, one can either use a numerical integration (either deterministic or stochastic) or to resort to approximations such as the Laplace’s method which allows to obtain a closed-form expression for the criterion to optimise.

Here, we consider these problems for the particular case of Gaussian prior laws:

$`p(𝒚|𝒙,\varphi ,k,l)`$ $`=`$ $`𝒩(𝑨𝒙,{\displaystyle \frac{1}{\varphi }}𝑰)=\left(2\pi /\varphi \right)^{m/2}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\varphi 𝒚𝑨𝒙^2\right]`$ (20)
$`p(𝒙|\psi ,k,l)`$ $`=`$ $`𝒩(\mathrm{𝟎},{\displaystyle \frac{1}{\psi }}𝑰)=\left(2\pi /\psi \right)^{k/2}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\psi 𝒙^2\right]`$ (21)
where $`\frac{1}{\varphi }`$ and $`\frac{1}{\psi }`$ are respectively the variance of the noise and the parameters.

## 2. Parameter estimation

First note that in this special case we have

$$p(𝒚,𝒙|\varphi ,\psi ,k,l)=\left(2\pi /\varphi \right)^{m/2}\left(2\pi /\psi \right)^{k/2}\mathrm{exp}\left[\frac{1}{2}\varphi 𝒚𝑨𝒙^2\frac{1}{2}\psi 𝒙^2\right].$$
(22)
Integration with respect to $`𝒙`$ can be done analytically and we have:

$$p(𝒚|\varphi ,\psi ,k,l)=p(𝒚,𝒙|\varphi ,\psi ,k,l)\text{d}𝒙=𝒩(\mathrm{𝟎},𝑷_y),$$
(23)
with

$$𝑷_y=\frac{1}{\psi }𝑨𝑨^t+\frac{1}{\varphi }𝑰=\frac{1}{\psi }(𝑨𝑨^t+\lambda 𝑰)\text{ and }\lambda =\frac{\psi }{\varphi }.$$
(24)
It is then easy to see that the a posteriori law of $`𝒙`$ is also Gaussian:

$$p(𝒙|𝒚,\varphi ,\psi ,k,l)=𝒩(\widehat{𝒙},\widehat{𝑷})\text{ with }\widehat{𝑷}=\frac{1}{\varphi }(𝑨^t𝑨+\lambda 𝑰)^1\text{ and }\widehat{𝒙}=\varphi \widehat{𝑷}𝑨^t𝒚.$$
(25)
Thus the parameter estimation in this case is straightforward:

$$\widehat{𝒙}=\underset{𝒙}{\text{arg}\mathrm{max}}\left\{p(𝒙|𝒚,\varphi ,\psi ,k,l)\right\}=\underset{𝒙}{\text{arg}\mathrm{min}}\left\{J_1(𝒙)\right\},$$
(26)
with

$$J_1(𝒙)=𝒚𝑨𝒙^2+\lambda 𝒙^2,$$
(27)
which is a quadratic function of $`𝒙`$. The solution is then a linear function of the data $`𝒚`$ and is given by

$$\widehat{𝒙}=𝑲(\lambda )𝒚\text{with}𝑲(\lambda )=(𝑨^t𝑨+\lambda 𝑰)^1𝑨^t.$$
(28)
## 3. Hyperparameter estimation

For the hyperparameter estimation problem we note that:

$`p(\mathit{\varphi },\psi |𝒚,k,l)`$ $`=`$ $`{\displaystyle \frac{p(\varphi |k,l)p(\psi |k,l)}{p(𝒚|k,l)}}p(𝒚|\varphi ,\psi ,k,l)`$ (29)
$`=`$ $`{\displaystyle \frac{p(\varphi |k,l)p(\psi |k,l)}{p(𝒚|k,l)}}(2\pi )^{m/2}\left|𝑷_y\right|^{1/2}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}𝒚^t𝑷_y^1𝒚\right].`$
Thus, the hyperparameter estimation problem becomes:

$$(\widehat{\varphi },\widehat{\psi })=\underset{(\varphi ,\psi )}{\text{arg}\mathrm{max}}\left\{p(\varphi ,\psi |𝒚,\widehat{k},\widehat{l})\right\}=\underset{(\varphi ,\psi )}{\text{arg}\mathrm{min}}\left\{J_2(\varphi ,\psi )\right\}$$
(30)
where

$$J_2(\varphi ,\psi )=\mathrm{ln}p(\varphi |k,l)\mathrm{ln}p(\psi |k,l)+\frac{1}{2}\mathrm{ln}\left|𝑷_y\right|+\frac{1}{2}𝒚^t𝑷_y^1𝒚.$$
(31)
Unfortunately, in general, there is not an analytical expression for the solution, but this optimization can be done numerically. Many works have been investigated to perform this optimization appropriately for particular choices of $`p(\varphi |k,l)`$ and $`p(\psi |k,l)`$. Among the others, we may note the choice of improper prior laws such as Jeffry’s prior $`p(\varphi |k,l)\frac{1}{\varphi }`$ and $`p(𝝍|k,l)\frac{1}{\psi }`$ or proper prior laws of uniform $`p(\varphi |k,l)=\frac{1}{\varphi _{\mathrm{max}}\varphi _{\mathrm{min}}}`$ and $`p(\psi |k,l)=\frac{1}{\psi _{\mathrm{max}}\psi _{\mathrm{min}}}`$ or still the proper Gamma prior laws.

One main issue with improper prior laws is the existence of the solution, because $`p(\mathit{\varphi },\psi |𝒚,k,l)`$ may even not have a maximum or its maximum can be located at the border of the domain of variation of $`(\varphi ,\psi )`$. Here, we propose to use the following proper Gamma priors :

$`p(\varphi )`$ $`=`$ $`𝒢(\alpha _1,\beta _1)\varphi ^{(\alpha _11)}\mathrm{exp}\left[\beta _1\varphi \right]\text{E}\left\{\varphi \right\}=\alpha _1/\beta _1`$ (32)
$`p(\psi )`$ $`=`$ $`𝒢(\alpha _2,\beta )\psi ^{(\alpha _21)}\mathrm{exp}\left[\beta _2\psi \right]\text{E}\left\{\psi \right\}=\alpha _2/\beta _2.`$ (33)
With these priors, we have

$$J_2(\varphi ,\psi )=(1\alpha _1)\mathrm{ln}\varphi +(1\alpha _2)\mathrm{ln}\psi +\beta _1\varphi +\beta _2\psi +\frac{1}{2}\mathrm{ln}\left|𝑷_y\right|+\frac{1}{2}𝒚^t𝑷_y^1𝒚.$$
(34)
The second main issue is the numerical optimization. Many works have been done on this subject. Among the others we can mention those who try to integrate out one of the two parameters directly or after some transformation. For example transforming $`(\varphi ,\psi )(\varphi ,\lambda )`$ and using the identities

$$\left|𝑨𝑨^t+\lambda 𝑰\right|=\lambda ^{mk}\left|𝑨^t𝑨+\lambda 𝑰\right|$$
(35)
and

$$(𝑨𝑨^t+\lambda 𝑰)^1=\frac{1}{\lambda }(𝑰𝑨𝑲(\lambda ))$$
(36)
we have

$$\mathrm{ln}\left|𝑷_y\right|=m\mathrm{ln}\varphi k\mathrm{ln}\lambda +\mathrm{ln}\left|𝑨^t𝑨+\lambda 𝑰\right|$$
(37)
and

$$𝒚^t𝑷_y^1𝒚=\varphi 𝒚^t(𝑰𝑨𝑲(\lambda ))𝒚=\varphi 𝒚^t(𝒚𝑨\widehat{𝒙})=\varphi 𝒚^t(𝒚\widehat{𝒚}).$$
(38)
Then, we obtain

$$\begin{array}{cc}J_2(\varphi ,\psi )=\hfill & (1\alpha _1\frac{mk}{2})\mathrm{ln}\varphi +(1\alpha _2\frac{k}{2})\mathrm{ln}\psi +\beta _1\varphi +\beta _2\psi \hfill \\ & +\frac{1}{2}\mathrm{ln}\left|𝑨^t𝑨+\frac{\psi }{\varphi }𝑰\right|+\frac{\varphi }{2}𝒚^t(𝒚\widehat{𝒚}).\hfill \end{array}$$
(39)
or

$$\begin{array}{cc}J_2(\varphi ,\lambda )=\hfill & (2\alpha _1\alpha _2\frac{m}{2})\mathrm{ln}\varphi +(1\alpha _2\frac{k}{2})\mathrm{ln}\lambda +\beta _1\varphi +\beta _2\varphi \lambda \hfill \\ & +\frac{1}{2}\mathrm{ln}\left|𝑨^t𝑨+\lambda 𝑰\right|+\frac{\varphi }{2}𝒚^t(𝒚\widehat{𝒚}).\hfill \end{array}$$
(40)
For fixed $`\lambda `$, equating to zero the derivative of this expression with respect to $`\varphi `$ has an explicite solution which is

$$\frac{J_2(\varphi ,\lambda )}{\varphi }=0\varphi =(\frac{m}{2}+\alpha _1+\alpha _22)/\left[\beta _1+\lambda \beta _2+\frac{1}{2}𝒚^t(𝒚\widehat{𝒚})\right]$$
(41)
Putting this expression into $`J_2`$ we obtain a criterion depending only on $`\lambda `$ which can be optimized numerically. In addition, it is possible to integrate out $`\varphi `$ to obtain $`p(\lambda |𝒚,k,l)`$, but the expression is too complex to write.

## 4. Joint estimation

One may try to estimate all the unknowns simultaneously by

$$(\widehat{𝒙},\widehat{\varphi },\widehat{\psi },\widehat{k},\widehat{l})=\underset{(𝒙,\varphi ,\psi ,k,l)}{\text{arg}\mathrm{max}}\{p(𝒙,\varphi ,\psi ,k,l|𝒚\}=\underset{(𝒙,\varphi ,\psi ,k,l)}{\text{arg}\mathrm{min}}\left\{J_3(𝒙,\varphi ,\psi ,k,l)\right\}$$
(42)
where

$$\begin{array}{cc}J_3(𝒙,\varphi ,\psi ,k,l)=\hfill & \mathrm{ln}p(k)\mathrm{ln}p(l)(\frac{m}{2}+\alpha _11)\mathrm{ln}\varphi \hfill \\ & (\frac{k}{2}+\alpha _21)\mathrm{ln}\psi +\varphi \left(\beta _1+\frac{1}{2}𝒚𝑨𝒙^2\right)+\psi \left(\beta _2+\frac{1}{2}𝒙^2\right)\hfill \end{array}$$
(43)
The main advantage of this criterion is that we obtain explicit solutions for $`𝒙`$, $`\varphi `$ and $`\psi `$ by equating to zero the derivatives of $`J_3(𝒙,\varphi ,\psi ,k,l)`$ with respect to them:

$$\{\begin{array}{c}\widehat{𝒙}=(𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒚,\text{with}\lambda =\varphi /\psi ;\hfill \\ \widehat{\varphi }=(\frac{m}{2}+\alpha _11)/\left(\beta _1+\frac{1}{2}𝒚𝑨\widehat{𝒙}^2\right);\hfill \\ \widehat{\psi }=(\frac{k}{2}+\alpha _21)/\left(\beta _2+\frac{1}{2}\widehat{𝒙}^2\right).\hfill \end{array}$$
(44)
We can not obtain closed form expressions for $`\widehat{k}`$ and $`\widehat{l}`$ which depend on the particular choice for $`p(k)`$ and $`p(l)`$. These relations suggest an iterative algorithm such as:

Joint MAP estimation algorithm 1
for $`l=1:l_{max}`$
for $`k=1:k_{max}`$
compute the elements of the matrix $`𝑨`$;
initialize $`\widehat{\lambda }=\lambda _0`$;
repeat until convergency:
$$\begin{array}{c}\widehat{𝒙}=(𝑨^t𝑨+\widehat{\lambda }𝑰)^1𝑨^t𝒚;\hfill \\ \{\begin{array}{c}\widehat{\varphi }=(\frac{m}{2}+\alpha _11)/\left(\beta _1+\frac{1}{2}𝒚𝑨\widehat{𝒙}^2\right);\hfill \\ \widehat{\psi }=(\frac{k}{2}+\alpha _21)/\left(\beta _2+\frac{1}{2}\widehat{𝒙}^2\right)\hfill \end{array}\widehat{\lambda }=\widehat{\varphi }/\widehat{\psi }\hfill \end{array}$$
end
compute $`J(k,l)=J_3(\widehat{𝒙},\widehat{\varphi },\widehat{\psi },k,l)`$;
end
end
choose the best model and the best order by
$`(\widehat{l},\widehat{k})=\text{arg}\mathrm{min}_{(k,l)}\left\{J(k,l)\right\}`$

Note however that, for fixed $`𝒙`$, $`\varphi `$ and $`\psi `$, the criteria $`J_3`$ in (43) or $`J_5`$ in (47) are mainly linear functions of $`k`$ if we choose a uniform law for $`p(k)`$. This means that we may not have a minimum for these criteria as a function of $`k`$. The choice of the prior $`p(k)`$ is then important. One possible is the following:

$$p(k)=\{\begin{array}{cc}\frac{2(k_{max}k)}{k_{max}(k_{max}1)}\hfill & 1k<k_{max}\hfill \\ 0\hfill & k>k_{max}\hfill \end{array}$$
(45)
which is a decreasing function of $`k`$ in the range $`k[1,k_{max}]`$ and zero elsewhere. This choice may insure the existence of a minimum if $`k_{max}`$ is choosed appropriately. For $`p(l)`$ we propose to choose a uniform law, because we do not want to give any favor to any model.

Another algorithm can be obtained if we replace the expression of $`\widehat{𝒙}`$ into $`J_3`$ to obtain a criterion depending only on $`(\varphi ,\psi )`$

$$\begin{array}{cc}J_4(\varphi ,\psi ,k,l)=\hfill & \mathrm{ln}p(k)\mathrm{ln}p(l)(\frac{m}{2}+\alpha _11)\mathrm{ln}\varphi (\frac{k}{2}+\alpha _21)\mathrm{ln}\psi \hfill \\ & +\varphi \left(\beta _1+\frac{1}{2}𝒚\widehat{𝒚}(\lambda )^2\right)+\psi \left(\beta _2+\frac{1}{2}\widehat{𝒙}(\lambda )^2\right)\hfill \end{array}$$
(46)
or on $`(\varphi ,\lambda )`$:

$$\begin{array}{cc}J_5(\varphi ,\lambda ,k,l)=\hfill & \mathrm{ln}p(k)\mathrm{ln}p(l)(\frac{m+k}{2}+\alpha _1+\alpha _22)\mathrm{ln}\varphi (\frac{k}{2}+\alpha _21)\mathrm{ln}\lambda \hfill \\ & +\varphi \left(\beta _1+\frac{1}{2}𝒚\widehat{𝒚}(\lambda )^2\right)+(\lambda \varphi )\left(\beta _2+\frac{1}{2}\widehat{𝒙}(\lambda )^2\right)\hfill \end{array}$$
(47)
and then optimize it with respect to them. In the second case, we can again obtain first $`\varphi `$ and put it’s expression

$$\widehat{\varphi }=\left(\frac{m+k}{2}+\alpha _1+\alpha _22\right)/\left[(\beta _1+\frac{1}{2}𝒚\widehat{𝒚}(\lambda )^2)+\lambda (\beta _2+\frac{1}{2}\widehat{𝒙}(\lambda )^2)\right]$$
(48)
in the criterion to obtain another criterion depending only on $`\lambda `$ and optimize it numerically. This gives the following algorithm:

Joint MAP estimation algorithm 2
for $`l=1:l_{max}`$
for $`k=1:k_{max}`$
compute the elements of the matrix $`𝑨`$;
for $`\lambda 10^{[8:1:4]}`$
compute $`\widehat{𝒙}=(𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒚`$ and $`\widehat{𝒚}=𝑨\widehat{𝒙}`$
compute $`\widehat{\varphi }`$ using (eq. 48)
compute $`J(\lambda )=J_5(\widehat{\varphi },\lambda ,k,l)`$ (eq. 47)
end
choose $`\widehat{\lambda }=\text{arg}\mathrm{min}_\lambda \left\{J(\lambda )\right\}`$
compute $`\widehat{𝒙}=(𝑨^t𝑨+\widehat{\lambda }𝑰)^1𝑨^t𝒚`$;
compute $`\widehat{\varphi }`$ using (eq. 48);
compute $`J(k,l)=J_5(\widehat{\varphi },\widehat{\lambda },k,l)`$ (eq. 47)
end
end
choose the best model and the best order by
$`(\widehat{l},\widehat{k})=\text{arg}\mathrm{min}_{(k,l)}\left\{J(k,l)\right\}`$

## 5. Model order selection

The model order selection

$$\widehat{k}=\underset{k}{\text{arg}\mathrm{max}}\left\{p(k|𝒚,l)\right\}=\underset{k}{\text{arg}\mathrm{min}}\left\{J_6(k)\right\}$$
(49)
with

$$J_6(k)=\mathrm{ln}p(k)\mathrm{ln}p(𝒚|k,l)$$
(50)
needs one more integration

$$p(𝒚|k,l)=p(𝒚,\varphi ,\psi |k,l)\text{d}\varphi \text{d}\psi .$$
(51)
or

$$p(𝒚|k,l)=p(𝒚,\varphi ,\lambda |k,l)\text{d}\varphi \text{d}\lambda .$$
(52)
where $`p(𝒚,\varphi ,\lambda |k,l)\mathrm{exp}\left[J_2(\varphi ,\lambda )\right]`$ given by (40). As we mentionned in preceeding section, these integrations can only be down numerically. A good approximation can be obtained using the following:

$$p(𝒚|k,l)=p(𝒚,\varphi ,\psi |k,l)\text{d}\varphi \text{d}\psi \underset{i}{}\underset{j}{}p(𝒚|\varphi _j,\psi _i,k,l)$$
(53)
where $`\{\varphi _j\}`$ and $`\{\psi _i\}`$ are samples generated using the prior laws $`p(\varphi )`$ and $`p(\psi )`$.

## 6. Best basis or model selection

The model selection

$$\widehat{l}=\underset{l}{\text{arg}\mathrm{max}}\left\{p(l|𝒚)\right\}=\underset{l}{\text{arg}\mathrm{min}}\left\{J_7(l)\right\}$$
(54)
with

$$J_7(l)=\mathrm{ln}p(l)\mathrm{ln}p(𝒚|l)$$
(55)
does not need any more integration, but only one summation. Choosing $`p(l)`$ uniform and making the same previous approximations we have

$$J_7(l)=\mathrm{ln}\underset{k=1}{\overset{k_{max}}{}}p(𝒚|k,l)p(k).$$
(56)
## 7. Proposed algorithms

Based on the equations (55), (53), (50), (39) and (40) we propose the following algorithm:
Marginal MAP estimation algorithm 2
Generate a set of samples $`\{\varphi _j\}`$ drawn from $`p(\varphi )`$
Generate a set of samples $`\{\psi _i\}`$ drawn from $`p(\psi )`$
for $`l=1:l_{max}`$
for $`k=1:k_{max}`$
compute the elements of the matrix $`𝑨`$;
for $`\varphi \{\varphi _j\}`$
for $`\psi \{\psi _i\}`$
compute $`\lambda =\varphi /\psi `$, $`\widehat{𝒙}=(𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒚`$ and $`\widehat{𝒚}=𝑨\widehat{𝒙}`$
compute $`p_\psi (i,j,k,l)=\mathrm{exp}\left[J_2(\varphi _j,\psi _i)\right]`$ (eq. 39)
end
normalize $`p_\psi (i,j,k,l)=p_\psi (i,j,k,l)/_ip_\psi (i,j,k,l)`$
compute $`p_\varphi (j,k,l)=_ip_\psi (i,j,k,l)`$
end
normalize $`p_\varphi (j,k,l)=p_\varphi (j,k,l)/_jp_\varphi (j,k,l)`$
compute $`p_k(k,l)=_jp_\varphi (j,k,l)`$
end
normalize $`p_k(k,l)=p_k(k,l)/_kp_k(k,l)`$
compute $`p_l(l)=_kp_k(k,l)`$
end
normalize $`p_l(l)=p_l(l)/_lp(l)`$
choose the best model by $`\widehat{l}=\text{arg}\mathrm{max}_l\left\{p_l(l)\right\}`$
choose the best model order by $`\widehat{k}=\text{arg}\mathrm{max}_k\left\{p_k(k,\widehat{l})\right\}`$
choose the best value for $`\varphi =\varphi _{\widehat{j}}`$ with $`\widehat{j}=\text{arg}\mathrm{max}_j\left\{p_\varphi (j,\widehat{l},\widehat{k})\right\}`$
choose the best value for $`\psi =\psi _{\widehat{i}}`$ with $`\widehat{i}=\text{arg}\mathrm{max}_i\left\{p_\psi (i,\widehat{j},\widehat{l},\widehat{k})\right\}`$
compute $`\widehat{\lambda }=\widehat{\varphi }/\widehat{\psi }`$
compute the elements of the matrix $`𝑨`$ for $`l=\widehat{l}`$ and $`k=\widehat{k}`$
compute $`\widehat{𝒙}=(𝑨^t𝑨+\widehat{\lambda }𝑰)^1𝑨^t𝒚`$.

## 8. Application: Electron scattering data inversion

Elastic electron scattering provides a mean of determining the charge density of a nucleus, $`\rho (r)`$, from the experimentally determined charge form factor, $`F(q)`$. The connection between the charge density and the cross section is well understood and in plane wave Born approximation $`F(q)`$ is just the Fourier transform of $`\rho (r)`$ which for the case of even-even nuclei, which we shall consider, is simply given by

$$F(q)=4\pi _0^{\mathrm{}}r^2J_0(qr)\rho (r)\text{d}r$$
(57)
where $`J_0`$ is the spherical Bessel function of zero order and $`q`$ is the absolute value of the three momentum transfer.

We applied the proposed method with the following usual discretization procedure:

$$\rho (r)=\{\begin{array}{cc}_{j=1}^kx_jb_j(r)\hfill & rR_c\hfill \\ 0\hfill & r>R_c\hfill \end{array}$$
(58)
which results in

$$F(q)=4\pi \underset{j=1}{\overset{k}{}}x_j_0^{R_c}r^2J_0(qr)b_j(r)\text{d}r$$
(59)
and

$$𝒚=𝑨𝒙+\mathit{ϵ}$$
(60)
where $`𝒙`$ is a vector containing the coefficients $`\{x_j,j=1,\mathrm{},k\}`$, $`𝒚`$ is a vector containing the form factor data $`\{F(q_i),i=1,\mathrm{},m\}`$ and $`𝑨`$ an $`(m\times k)`$ matrix containing the coefficients $`A_{i,j}`$ given by

$$A_{i,j}=4\pi _0^{R_c}r^2J_0(q_ir)b_j(r)\text{d}r.$$
(61)
To compute $`A_{i,j}`$ we define a discretization step $`\mathrm{\Delta }r=R_c/N`$, a vector $`𝒓=\{r_n=(n1)\mathrm{\Delta }r,n=1,\mathrm{},N\}`$, a $`(N\times k)`$ matrix $`𝑩`$ with elements $`B_{n,j}=b_j(r_n)`$, a $`(m\times N)`$ matrix $`𝑪`$ with elements $`C_{i,n}=(4\pi \mathrm{\Delta }r)r_n^2J_0(q_ir_n)`$ such that we have $`𝑨=𝑪𝑩`$. Note also that when the vector $`𝒙`$ is determined, we can compute $`𝝆=\{\rho (r_n),n=1,\mathrm{},N\}`$ by $`𝝆=𝑩𝒙`$.

To test the proposed methods, we used the following simulation procedure:

* Select a model type $`l`$ and an order $`k`$ and generate the the matrixes $`𝑩`$, $`𝑪`$ and $`𝑨`$, and for a random set of parameters $`𝒙`$ generate the data $`𝒚=𝑨𝒙`$.
* Add some noise $`\mathit{ϵ}`$ on $`𝒚`$ to obtain $`𝒚=𝑨𝒙+\mathit{ϵ}`$;
* Compute the estimates $`\widehat{l}`$, $`\widehat{k}`$, $`\widehat{𝒙}`$, $`\widehat{𝒚}=𝑨\widehat{𝒙}`$ and $`\widehat{𝝆}=𝑩\widehat{𝒙}`$ and compare them with $`l`$, $`k`$, $`𝒙`$, $`𝒚=𝑨𝒙`$ and $`𝝆=𝑩𝒙`$.

We choosed the following basis functions:

* $`l=1`$ : $`b_j(r)=J_0(q_jr)`$— This is a natural choice due to the integral kernel and the orthogonality property of the bessel functions.
* $`l=2`$ : $`b_j(r)=\text{sinc}(q_jr)`$— This choice is also natural due to the orthogonality and the limited support hypothesis for the function $`\rho (r)`$.
* $`l=3`$ : $`b_j(r)=\mathrm{exp}\left[\frac{1}{2}(q_jr)^2\right]`$— This choice can account for the positivity of the function $`\rho (r)`$ if $`\{x_j\}`$ are constrained to be positive.
* $`l=4`$ : $`b_j(r)=\mathrm{exp}\left[\frac{1}{2}(q_jr)^2\right]J_0(q_jr)`$— This choice combines the first and the third properties.
* $`l=5`$ : $`b_j(r)=1/(\text{cosh}(q_jr))`$— This choice has the same properties of the third one.
* $`l=6`$ : $`b_j(r)=1/(1+(q_jr)^2)`$— This choice has the same properties of the third one.

In all these experiments we choosed $`k=6`$, $`m=20`$, $`N=100`$, $`R_c=8`$ and $`q_i=i\pi /R_c`$. The following figures show typical solutions. Figures 1 and 2 show the details of the procedure for the case $`l=1`$. Figures 3, 4 and 5 show the results for the cases $`l=1`$ to $`l=6`$.

|  | $`b_j(r)`$  $`j=1,\mathrm{},k=6`$ |
| --- | --- |
|  | $`\rho (r)={\displaystyle \underset{j=1}{\overset{k}{}}}x_jb_j(r)`$  with  $`x_1=x_2=\mathrm{}=x_6=1`$ |
|  | $`F(q_i)={\displaystyle \underset{j=1}{\overset{k}{}}}A_{ij}x_j`$  $`A_{ij}={\displaystyle r^2J_0(q_ir)b_j(r)\text{d}r}`$  $`i=1,\mathrm{},m=14`$  $`j=1,\mathrm{},k=6`$ |

Fig. 1: a) basis functions $`b_j(r)`$, b) $`\rho (r)`$, c) data $`F(q_i)`$ in a logarithmic scale.

|  | $`p(k,l)`$  $`p(l)`$ and $`p(k|\widehat{l})`$ |
| --- | --- |
|  | $`\widehat{\rho }(r)={\displaystyle \underset{j=1}{\overset{k}{}}}\widehat{x}_jb_j(r)`$  and  $`\rho (r)={\displaystyle \underset{j=1}{\overset{k}{}}}x_jb_j(r)`$ |
|  | $`\widehat{F}(q_i)={\displaystyle \underset{j=1}{\overset{k}{}}}A_{ij}\widehat{x}_j`$  and  $`F(q_i)={\displaystyle \underset{j=1}{\overset{k}{}}}A_{ij}x_j`$ |
| Fig. 2: a) $`p(k,l|𝒚)`$, $`p(l|𝒚)`$ and $`p(k|𝒚,\widehat{l})`$, b) original $`\rho (r)`$ |  |
| and estimated $`\widehat{\rho }(r)`$, c) original $`F(q_i)`$ and estimated $`\widehat{F}(q_i)`$. |  |

Fig. 3: Left: $`l=1`$ Right: $`l=2`$   a) basis functions $`b_j(r)`$;    b) $`p(k,l|𝒚)`$;    c) $`p(k|𝒚)`$ and $`p(k|𝒚,\widehat{l})`$;    d) $`\rho (r)`$ and $`\widehat{\rho }(r)`$;    e) $`F(q_i)`$ and $`\widehat{F}(q_i)`$.

Fig. 4: Left: $`l=3`$ Right: $`l=4`$   a) basis functions $`b_j(r)`$;    b) $`p(k,l|𝒚)`$;    c) $`p(k|𝒚)`$ and $`p(k|𝒚,\widehat{l})`$;    d) $`\rho (r)`$ and $`\widehat{\rho }(r)`$;    e) $`F(q_i)`$ and $`\widehat{F}(q_i)`$.

Fig. 5: Left: $`l=5`$ Right: $`l=6`$   a) basis functions $`b_j(r)`$;    b) $`p(k,l|𝒚)`$;    c) $`p(k|𝒚)`$ and $`p(k|𝒚,\widehat{l})`$;    d) $`\rho (r)`$ and $`\widehat{\rho }(r)`$;    e) $`F(q_i)`$ and $`\widehat{F}(q_i)`$.

Note that in these tests, we know perfectly the model and generated the data according to our hypothesis. To test the method to a more realistic case, we chose a model for which we can have an exact analytic expression for the integrals. For example, if we choose a symmetric Fermi distribution

$$\rho (r)=\alpha \frac{\mathrm{cosh}(R/d)}{\mathrm{cosh}(R/d)+\mathrm{cosh}(r/d)}$$
(62)
an analytical expression for the corresponding charge form factor can easily be obtained :

$$F(q)=\frac{4\pi ^2\alpha d}{q}\frac{\mathrm{cosh}(R/d)}{\mathrm{sinh}(R/d)}\left[\frac{R\mathrm{cos}(qR)}{\mathrm{sinh}(\pi qd)}\frac{\pi d\mathrm{sin}(qR)\mathrm{cosh}(\pi qd)}{\mathrm{sinh}^2(\pi qd)}\right].$$
(63)
Only two of the parameters $`\alpha `$, $`R`$ and $`d`$ are independent since the charge density must fulfill the normalization condition

$$4\pi r^2\rho (r)\text{d}r=Z.$$
(64)
Figure 6 shows the theoretical charge density $`\rho (r)`$ of <sup>12</sup>C (Z=6) obtained from (62) for $`r[0,8]`$ fm with $`R=1.1`$ $`A`$ and $`d=0.626`$ fm and the theoretical charge form factor $`F(q)`$ obtained by (63) for $`q[0,8]`$ fm<sup>-1</sup> and the 15 simulated data:

$$𝒒=[0.001,.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0]\text{fm}^1$$

which are used as inputs to the inversion method.

Fig. 6: Theoretical charge density $`\rho (r)`$, charge form factor $`\mathrm{log}|F(q)|`$ and the data \[stars\] used for numerical experiments \[right\].

First note that, even with the exact data, there are an infinite number of solutions which fits exactly the data. The following figure shows a few set of these solutions.

Fig. 7: a) $`p(k,l|𝒚)`$;    b) $`p(k|𝒚)`$ and $`p(k|𝒚,\widehat{l})`$;    c) $`\rho (r)`$ and $`\widehat{\rho }(r)`$; d) $`F(q_i)`$ and $`\widehat{F}(q_i)`$.

## 9. Conclusions

We discussed the different steps for a complete resolution of an inverse problem and focused on the choice of a basis function selection and the order of the model. An algorithm based on Bayesian estimation is proposed and tested on simulated data.