Alain Sarlette

Mazyar Mirrahimi

Loss-tolerant parity measurement for distant quantum bits

Keywords:

We propose a scheme to measure the parity of two distant qubits, while ensuring that losses on the quantum channel between them does not destroy coherences within the parity subspaces. This capability enables deterministic preparation of highly entangled qubit states whose fidelity is not limited by the transmission loss.

The key observation is that, for a probe electromagnetic field in a particular quantum state; namely, a superposition of two coherent states of opposite phases, the transmission loss stochastically applies a near-unitary backaction on the probe state. This leads to a parity measurement protocol where the main effect of the transmission losses is a decrease in the measurement strength. By repeating the nondestructive (weak) parity measurement, one achieves a high-fidelity entanglement in spite of a significant transmission loss.

I. INTRODUCTION

The correlation of distant systems through entanglement is a hallmark of quantum physics [1,2] and plays a fundamental role in envisioned quantum technology. Most fundamentally, quantum teleportation [3] shows how entanglement is a resource for effectively transmitting the unknown state of a quantum system, without physically transmitting quantum states. Towards future quantum computers, such teleportation could transport information between few-qubit processing units and memory units which must be well isolated and hence not directly coupled to the rest of the system via tunable physical interactions. This so-called modular architecture for quantum computing provides a viable solution to the major scaling problem for many-qubit quantum information processing [4,5]. In quantum communication, similar ideas would allow quantum repeaters to purify information through local operations only, provided they can consume units of entanglement between the communicating devices [6,7]; such quantum repeaters are a necessary technology for accurate quantum communication, with associated, e.g., cryptographic benefits, over long distances.

A major challenge towards enabling these applications is that generating entangled states between distant systems must rely itself on a quantum channel (see impossibility with local operations and classical communication). Microwave experiments have demonstrated how to deterministically entangle separate quantum subsystems via parity measurements [8], yet with a fidelity directly limited by the quality of the quantum channel: any losses on the probe field imply losses in entanglement. Channel losses can also be made to affect preparation success probability, instead of preparation fidelity. Many theoretical proposals [9][10][11][12][13] and experiments [14][15][16][17][18][19] have considered high-fidelity entanglement generation, despite propagation losses, by heralding the preparation on some rare photodetection events. The success rate of this probabilistic preparation is, however, very low in current experiments, while the preparation fidelity is still limited by imperfections such as dark counts of the photodetector. Recent developments in quantum superconducting circuit setups, like Refs. [20][21][22][23][24][25], show that local operations are becoming extremely reliable and fast, e.g., manipulating "cat states" as efficiently as single photons and having continuous field measurements instead of photon counters. The remaining bottleneck towards entanglement stabilization is then the channel loss.

This paper provides an explicit proposal for solving this problem through the achievement of loss-tolerant degeneracypreserving quantum nondemolition (QND) parity measurements between spatially separated qubits. Our key idea is to transmit over the quantum channel particularly engineered quantum states of light, i.e., "cat states", for which the dominant photon loss errors almost reduce to photon-number parity flips [26]. With this we design the interaction between qubits and probe field such that (i) measuring the probe at the output performs a QND measurement of qubit parity and (ii) photon loss events on the transmitted probe field render the detection less decisive (weak measurement) but, to a very good approximation, do not break the degeneracy of the parity subspaces. This measurement enables an efficient entanglement stabilization scheme, which not only generates highly entangled states but also protects them over arbitrarily long times against various noise channels.

The rather standard abstract setting [Fig. 1(a)] comprises two target qubits |q A , |q B at different locations A, B and possibly embedded in auxiliary quantum machinery, e.g., a cavity in circuit quantum electrodynamics (QED) setups [27]. For each measurement, a source generates a controlled "probe" quantum state |ψ p at A which then interacts with |q A according to a unitary U A , is transmitted over a noisy quantum channel C, interacts with |q B according to U B , and finally hits a detector at B. Those probe states play the role of parity meter. Since the quantum channel is the unequivocal bottleneck for remote entanglement in state-of-the-art technology [28], we focus on this issue and assume in this paper that the local actions (i.e., U A , U B , generating |ψ p , detection at B) are implemented perfectly. As we will show later, imperfections on such operations will have a small effect on the steady state when we stabilize quantum entanglement via feedback. Q but ensures that every eigenspace of Q remains unaffected for all possible detection results. When Q is degenerate, this only ensures that a state inside an eigenspace is sent to a state in the same eigenspace. A degeneracy-preserving quantum nondemolition (DP-QND) measurement is a QND measurement that does not affect any eigenstate of Q. In other words, the DP-QND measurement acts as the identity on each eigenspace [13]. Consider the parity observable

Q = Q + -Q -for two qubits in the canonical basis {|0 ,|1 }, with Q + = |00 00| + |11 11|, Q -= |01 01| + |10 10|.
In a projective measurement of Q, detection result + (-) would project the two qubits onto the even parity manifold span{|00 ,|11 } (the odd parity manifold span{|01 ,|10 }). Less decisive DP-QND measurements can result, for instance, from classical detection uncertainty, e.g., with probability 1ξ an even (odd) parity state gives detection result -(+). Then the probability to detect ± becomes p ± = ξ ψ|Q ± |ψ + (1ξ ) ψ|Q ∓ |ψ with corresponding measurement backaction transforming the initial state |ψ into [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF] K±

(|ψ ψ|) = ξQ ± |ψ ψ|Q ± + (1 -ξ )Q ∓ |ψ ψ|Q ∓ p ± .
Any even-parity (odd-parity) state remains unchanged under such measurement backaction and predominantly gives detection result + (-), i.e., with probability ξ ∈ (1/2,1]. We first sketch the general concept of our proposal for DP-QND measurement of parity, with |ψ p = |q p an abstract probe qubit, similarly to the distillation discussion in e.g., Ref. this can be represented as conditioning each bit flip on an unknown hypothetical state of the environment. For n even, the outcome is as for the perfect channel. For n odd the final probe state is reversed, e.g., input (|11 ± |00 ) A,B |0 p / √ 2 yields output (|11 ± |00 ) A,B |1 p / √ 2, but most importantly, the target qubits' state remains unaffected. This essential property ensures that the expected evolution for n unknown (equivalently, tracing over the unknown states of the environment) remains a DP-QND parity measurement, explicitly described by operators K+ , Kwith ξ = n even Prob(n) < 1. One can easily adapt this ξ to include detection misses and errors. A broad distribution of values of n implies low measurement contrast, pushing ξ close to 1/2, but it does not impede its DP-QND character. Hence when sufficiently many measurements can be repeated within a relevant timescale, a conclusive result is obtained even for ξ close to 1/2 (see details below).

This DP-QND property is, in fact unavoidably, lost under general channel errors (see Appendix A). If, for instance, C includes a phase flip E 1 = |0 0| -|1 1| p , then the initial even-parity Bell state |ψ + = (|11 +|00 ) A,B / √ 2 gets transformed by measurement backaction into |ψ -= (|11 -|00 ) A,B / √ 2; thus this eigenstate of Q is not conserved, breaking the DP-QND character. Not knowing if E 1 was applied or not, the initial maximally entangled state |ψ + gets transformed into a statistical mixture of |ψ + and |ψ -, i.e., entanglement is lost.

III. TOWARDS A PHYSICAL IMPLEMENTATION

These abstract properties indicate a pathway towards losstolerant parity measurement: use a subspace of probe states on which the dominant decoherence channel acts as a unitary root of identity, e.g., a bit flip. Our key observation is that this fits a physical implementation where the probe is an electromagnetic field pulse whose logical states are materialized by so-called "cat states", i.e., mesoscopic superpositions of coherent states. This encoding, proposed earlier for quantum computing [26,[START_REF] Mirrahimi | [END_REF], ensures that photon losses imply to a good approximation coherent logical bit flips. The rest of this paper describes such implementation in detail.

A. Probe's cat states

This physical implementation of loss-tolerant parity measurement is sketched in Fig. 1(c). We denote by

|C ± β = (|β ± |-β )/N ± β , N ± β = 2 ± 2e -2|β| 2 ,
the superpositions between two coherent states |β and |-β , β ∈ R; the normalization constants N ± β rapidly approach √ 2 as the coherent amplitude β increases. The probe field is initially prepared in |C + α p and interacts with two qubitcavity systems in a cascaded manner. Between the two setups it is exposed to losses, modeled as mixing with the vacuum state |0 env of an ancillary mode through a beam-splitter unitary operator U η BS with transmittance √ η and reflectance √ 1η. The unitary operators apply:

U A |0 A |C ± α p = |0 A |C ± α p , U A |1 A |C ± α p = |1 A |C ∓ α p , U B |0 B |C ± √ ηα p = |0 B |C ± √ ηα p , U B |1 B |C ± √ ηα p = |1 B |C ∓ √ ηα p , U η BS |C ± α p |0 env = 1 N ± α (| √ ηα p | 1 -ηα env ± |- √ ηα p |-1 -ηα env ). (1)
Finally, after interaction with the second qubit, a measurement (e.g., of photon-number parity observable) projects the probe's state onto

|C + √ ηα p or |C - √ ηα p . Calling |C + β p and |C - β p
respectively the logical |0 p and |1 p , we recover the scheme where U A , U B implement CNOT gates. Before analyzing in detail the cascade scheme of Fig. 1(c), we want to introduce a corresponding concurrent scheme, depicted in Fig. 1(d). The initial 50/50 beam splitter generates a probe state of two entangled fields, (|C

+ α pA |C - α pB + |C - α pA |C + α pB )/ √ 2, with α = γ / √ 2.
These two output fields are sent along lossy quantum channels towards the respective qubits where they undergo the equivalent of a CNOT. After the final beam-splitter, a similar measurement to the cascaded scheme is applied to the channel A . While the rest of the paper will deal with the analysis of the cascaded scheme, all the derivations can be adapted to this concurrent architecture in a straightforward manner.

B. Performance analysis

We now analyze in detail the performance of the scheme of Fig. 1(c). As we trace on the ancillary field modeling the losses, we can without loss of generality imagine a virtual detector also projecting it to |C + √ 1-ηα env or |C - √ 1-ηα env , but with unread detection result. The measurement outcomes of the two detectors (one real and one virtual) are associated with four Kraus operators M σ,σ and M σ,-σ with σ ∈ {+,-}, modeling the backaction of the measurements on the target qubits: e.g., their state ρ, after measuring even parities for both detectors, is modified to M +,+ ρM † +,+ /Tr(M +,+ ρM † +,+ ). Following the computations of Appendix B, these Kraus operators are

M σ,σ = N σ √ 1-ηα 2 ⎛ ⎝ N σ √ ηα N + α |00 00| + N (-σ ) √ ηα N - α |11 11| ⎞ ⎠ , M σ,-σ = N (-σ ) √ 1-ηα 2 ⎛ ⎝ N σ √ ηα N - α |10 10| + N (-σ ) √ ηα N + α |01 01| ⎞ ⎠ . (2)
Discarding the inaccessible outcome of the virtual detector, the backaction induced by the probe field measurement follows:

K + (ρ) = M +,+ ρM † +,+ + M +,-ρM † +,- Tr(M +,+ ρM † +,+ + M +,-ρM † +,-) , K -(ρ) = M -,+ ρM † -,+ + M -,-ρM † -,- Tr(M -,+ ρM † -,+ + M -,-ρM † -,-) . ( 3
)
In the lossless case (η = 1), M +,-and M -,-vanish as N - 0 = 0; the coefficients in front of |00 00| and |11 11| in M +,+ and in front of |01 01| and |10 10| in M -,+ are identical, equal to 1. This is a projective parity measurement, as described by K+ , Kwith ξ = 1.

The effect of transmission losses (η < 1) is twofold. First, it reduces the measurement strength. Indeed, when the probe field is detected in a given parity (e.g., +), the qubits could be projected to the opposite parity (e.g., by the Kraus operator M +,-). However, each measurement does increase the conditional probability of finding the qubits in the same parity as the one indicated by the detections, because

N + √ 1-ηα min N + √ ηα N + α , N - √ ηα N - α > N - √ 1-ηα max N + √ ηα N - α , N - √ ηα N + α .
A projective parity measurement under perfect transmission (η = 1) is thus replaced by a less decisive measurement for η < 1, where at each shot we gain partial information on the parity. By repeating the measurement the state gets projected onto a well-defined parity subspace. Most importantly, we preserve the QND property: an initial state of definite parity always keeps this parity [e.g., M +,-

(|00 + |11 ) = 0].
However, the second, more harmful effect of transmission loss is a slight perturbation of the DP-QND property, by introducing slow mixing within each parity subspace. This is due to the coherent states |β and |-β not being perfectly orthogonal, so that N + β = N - β . This effectively induces dephasing inside parity subspaces, e.g.,

N + √ ηα /N + α > N - √ ηα /N - α implies that M +,+ drives |B e ± = (|00 ± |11 )/ √ 2 closer to |00 = (|B e + + |B e -)/ √ 2, while M -,-would drive them closer to |11 = (|B e + -|B e -)/ √ 2.
The simulations in Fig. 2 illustrate the competition between parity measurement and undesired dephasing for various |α| 2 (average number of photons in the probe field). Initializing both qubits in the state √ ηα /N + α and N - √ ηα /N - α decreases and this undesired dephasing gets suppressed significantly, at the expense of slower convergence, i.e., weaker parity measurement. As soon as η > 1/2, one can achieve arbitrarily high fidelity in this way.

|+ X = (|0 + |1 )/ √ 2,
In the next section, we give approximate analytical formulas for these two competing rates, providing the maximal achievable fidelity as a function of η, |α| 2 and qubit lifetimes.

C. Measurement and dephasing rates

The rate at which the parity measurement acquires information can be expressed analytically as a function of η and |α| 2 . Define the Lyapunov function

V parity (ρ) = 00|ρ|00 10|ρ|10 + 11|ρ|11 01|ρ|01 = w B0 2 1 -(P B0 ) 2 + w B1 2 1 -(P B1 ) 2 , ( 4
)
where w Bq B = q B |ρ|q B B denotes the population with qubit B in state |q B and P Bq B = Tr(σ z ⊗ σ z ρ|q B q B |)/w Bq B measures the parity, conditioned on qubit B being in state |q B , for q B ∈ {0,1}. One computes that

V parity (ρ k+1 ) = 1 -e -4(1-η)|α| 2
1e -4|α| 2 V parity (ρ k ) , (5) where ρ k is the joint qubits state after the kth measurement and V parity denotes the ensemble average of V parity over measurement realizations. Indeed, we have

E(V parity (ρ k+1 )|ρ k ) = V parity (K + (ρ k ))P(+|ρ k ) + V parity (K -(ρ k ))P(-|ρ k ),
where

P(+|ρ k ) = Tr(M +,+ ρ k M † +,+ + M +,-ρ k M † +,-), P(-|ρ k ) = Tr(M -,+ ρ k M † -,+ + M -,-ρ k M † -,-)
are respectively the conditional probabilities of achieving a positive or negative outcome at kth measurement. From this, we get

E(V parity (ρ k+1 )|ρ k ) = 00|M +,+ ρ k M † +,+ |00 10|M +,-ρ k M † +,-|10 + 11|M +,+ ρ k M † +,+ |11 01|M +,-ρ k M † +,-|01 + 00|M -,-ρ k M † -,-|00 10|M -,+ ρ k M † -,+ |10 + 11|M -,-ρ k M † -,-|11 01|M -,+ ρ k M † -,+ |01 = N + √ 1-ηα N - √ 1-ηα N + α N - α V parity (ρ k ).
Taking the expectation value of both sides with respect to ρ k , we find the result of Eq. ( 5). Thus, V parity exponentially decays to zero at a rate (per measurement iteration)

r parity = 1 2 ln 1 -e -4|α| 2 1 -e -4(1-η)|α| 2 .
In the case of a fully DP-QND measurement, i.e., assuming

N + √ ηα /N + α N - √ ηα /N - α and N + √ ηα /N - α N - √ ηα /N +
α in the Kraus operators, the two terms in the Lyapunov function ( 4) would decay at the same rate and r parity represents precisely the parity measurement strength. Indeed, in this case we would obtain the same result with the alternative Lyapunov function

V parity ideal (ρ) = 1 2 [1 -P (ρ) 2 ] 1/2
, where P (ρ) = Tr(σ z ⊗ σ z ρ). One can also analytically calculate the dephasing rate induced by the transmission loss. To this aim we define the coherence function

C(ρ) = | 00|ρ|11 | + | 01|ρ|10 |.

Similar calculations as above lead to

C(ρ k+1 ) = 1 -e -4η|α| 2
1e -4|α| 2 C(ρ k ) . + and |B o + , in simulations of our measurement scheme, when the two qubits are initialized in the state (|0 + |1 ) A (|0 + |1 ) B /2 and we take an average over 1000 Monte Carlo trajectories. The black (dashed) curve illustrates the average parity converging at the characteristic rate of the Lyapunov function V parity (ρ). The red (dashed) curve, decreasing with the same characteristic rate as the coherence function C(ρ), illustrates the dephasing inside a parity eigenspace. The green (dashed) curve represents the product of both these effects. This theoretical curve fits reasonably well with the simulations (blue solid curve).

Therefore the coherence function C(ρ) exponentially decays to zero at a rate (per measurement iteration)

r dephasing = 1 2 ln 1 -e -4|α| 2 1 -e -4η|α| 2 .
Figure 3 provides a comparison of the simulations to the above analytical results. As can be seen, the performance of the parity measurement protocol can be well explained by using the above two rates. The slight mismatch between the raising rate of the blue curve (fidelity to Bell states) and the theoretical r parity can be explained by the nonlinear relation between the fidelity and the Lyapunov function V parity . This nonlinearity makes it impossible to translate the rate r parity into a precise exponential rate for the raising dynamics of the average fidelity.

The expressions for r parity and r dephasing allow us to analyze the measurement performance as a function of η and |α| 2 . Consider again the evolution depicted in Figs. 2 and3. The fidelity to the closest Bell state is 1/2 times the sum of two terms: dominant parity population, which converges from 1/2 to 1 at roughly a rate r parity , and dominant phase population, which decreases from 1 to 1/2 at a rate r dephasing . The two terms contribute equally to the error after a number of measurements τ = τ meas satisfying e -r parity τ + e -r dephasing τ = 1. The corresponding estimate of Bell-state fidelity is F meas = 1e -r parity τ meas /2. This is a good estimate of the maximum reached by the curves on Figs. 2 and3. Figure 4 illustrates these indicators of our parity measurement's performance. For transmission efficiency as low as 70% we can get F meas as high as 99% in 400 measurement steps, taking |α| 2 = 3.273. Increasing η to 85% one achieves the same F meas with only 17 measurements, taking |α| 2 = 1.63.

In the limit r parity r dephasing (equivalent to e 4(2η-1)|α| 2 1), one can approximately replace e -r dephasing τ by 1r dephasing τ , and therefore the solution to the transcendental equation e -r parity τ + e -r dephasing τ = 1, giving the point of maximal fidelity, is well estimated by

τ meas = 1 r parity W 0 r parity r dephasing ,
where W 0 is the Lambert W function.

IV. FEEDBACK STABILIZATION OF BELL STATES

Our near DP-QND parity measurement can be used to stabilize a particular Bell state with a simple feedback protocol. To stabilize |B e + , (i) apply a π pulse around the X axis on the first qubit whenever the measurements estimate a probability higher than 1/2 to be in the odd-parity subspace, (ii) after that, apply a π/2 pulse on both qubits around the Y axis irrespectively of the detection result. The measurement backaction favors convergence towards the dominant parity, the π pulse correcting the parity whenever this is not even. This pushes the state towards the span of |B e + , |B e -without favoring the target |B e + . The two π/2 pulses then leave |B e + untouched and send the undesired |B e -onto |B o + , such that the next parity measurement stochastically moves the corresponding population as well towards the target. Simulations in Fig. 5 illustrate the performance of this protocol, having fixed η = 0.75 and varying |α| 2 . As expected, a larger |α| 2 allows us to reach higher final fidelity but also implies slower convergence.

The information needed to implement general feedback decisions is typically obtained by a quantum filter, i.e., a computer estimating the state by simulating the evolution (3) associated with the respective detection results. Such filters are known to be stable [31]. In the present case, because the feedback decision requires only to know the respective parity populations, the quantum filter can be simplified significantly; namely, by discarding all off-diagonal components in the Bell-state basis. Indeed, first note that the measurement does not depend on coherences among subspaces of different parity. Moreover, the measurement itself completely destroys, in a single iteration, any coherences between subspaces of different parity. The π/2 pulses "export" coherences possibly present between, e.g., |B e + and |B e -into coherences between |B e + and |B o + , i.e., between parity subspaces, are thus destroyed at the next measurement; while they "import" into, e.g., the even-parity eigenspace, the coherences previously present between |B e + and |B o + , which are none since those two Bell states belong to different-parity eigenspaces and their coherences were thus destroyed by the last measurement. Thus, after one initial weak measurement step at most, no relevant coherences among Bell states will survive. This allows us to update just the populations on the four Bell states, as in a classical filter for a partially observed Markov chain. This filter, updating the populations on the four Bell states, is the only computation required for the feedback: the action itself just requires a binary decision, to switch towards the predominantly populated parity.

To analyze the feedback more explicitly, we note that instead of applying the π/2 pulses to rotate the state in the Schrödinger picture, we can reformulate the dynamics by applying the π/2 pulses to the measurement scheme, in an equivalent Heisenberg picture. The corresponding measurements then alternate between the measurement σ z ⊗ in the z basis, and a parity measurement σ x ⊗ σ x in the x basis. The corresponding state is subject to π pulses in the respective z or x basis. A π pulse in the z basis (x basis) does not change the populations of x-parity subspaces (of z-parity subspaces) and is applied to increase the population in span{|B e + ,|B e -} (span{|B e + ,|B o + }). The stochastic convergence of the measured system towards a definite parity in both x and z coordinates ensures that, in the absence of other effects, the population in span{|B e + ,|B e -} ∩ span{|B e + ,|B o + } = |B e + would increase in expectation until reaching 100%. The dephasing effect limits the actual fidelity for a given α. Taking α larger allows us to increase the fidelity, but slows down the convergence, such that a trade-off value of α must be selected with respect to other decoherence effects acting on the system.

We can quantify the performance of this feedback scheme by using the characteristic convergence rates computed in the previous section. Noting that the π/2 pulses ensure a rapid oscillation between the states |B e -and |B o + , we define ξ (t) to be their combined population at time t. For r parity , r dephasing 1 i.e., a small rate of change per measurement iteration, we can approximate the system dynamics by a continuous-time model for the populations p k of Bell state |k (and ξ = p B e -+ p B o + ):

d dt p B e + = -
r dephasing 2T d p B e + + r parity 2T d + r dephasing 4T d ξ, d dt p B o -= - r parity T d + r dephasing 2T d p B o -+ r dephasing 4T d ξ, d dt ξ = - r parity 2T d + r dephasing 2T d ξ + r dephasing 2T d p B e + + r parity T d + r dephasing 2T d p B o -. (6)
Here t denotes actual time and T d is the time between two consecutive detection results. To obtain these equations, we have used the fact that the measurement destroys any coherences among Bell states, so only populations are relevant. The nonlinear relation between V parity and the fidelities makes it impossible to be more accurate than order of magnitude about r parity , and we have observed in simulations that the true dynamics is more complicated than the linear model (6). However, the model (6) gives valuable indications about the dependence of the whole scheme on α (through r parity and r dephasing ) and allows us to efficiently optimize α for given settings. The steady state of Eq. ( 6) is obtained with 2 , where δ = 2r parity r dephasing + 1, which increases with α as expected.

p B e + = δ 2 (1 + δ)
To practically illustrate the trade-off with slower convergence at larger α, we can add to the picture the qubit relaxation and qubit dephasing:

d dt ρ = k∈{A,B} L k ρL † k - 1 2 (L † k L k ρ + ρL † k L k ) + k∈{A,B} Lk ρ L † k - 1 2 ( L † k Lk ρ + ρ L † k Lk ), with L k = 1 T k 1 |0 1| k , Lk = 1 2T k φ σ k z . ( 7
)
Here |0 1| k is the qubit k's lowering operator, σ k z is just σ z on qubit k, T k 1 is the relaxation time, and T k φ is the pure dephasing time. Note that we apply these decoherences in the laboratory coordinates, where the qubits undergo π/2 rotations after each measurement. One then computes that, in the Bell state basis, when discarding off-diagonal terms in ρ, the effect of this qubit decoherence is equivalent to replacing Eq. ( 6

= r parity T d , 1 = 1 2T A 1 + 1 2T B 1 , 2 = r dephasing T d + 1 2T A 1 + 1 2T B 1 + 1 T A φ + 1 T B φ .
Solving for the steady state gives

p B e + = (1 + δ) 1 / + ( 2 4(1 + δ) 1 / 2 + (1 + δ) 2 , with δ = 2r parity T d 2 + 1. ( 8
)
From this expression, the optimum value of α can easily be computed numerically. We have simulated the actual system dynamics, i.e., adding qubit decoherence (7) to (3), for different values of η and T d /T 1 (we have taken

T A 1 = T B 1 = T 1 and T A,B φ = ∞).
For α we have taken the value that maximizes Eq. ( 8), as well as a few values close to it in order to confirm that we hit close to the actual optimum fidelity. For each set of parameters, we have performed 5000 Monte Carlo simulations in order to estimate the achieved fidelity. Figure 6 shows these simulation results. The value of α obtained by maximizing Eq. ( 8) appears to indeed be close to optimal. While the above simple scheme is certainly not the optimal feedback strategy, the simulation results on Fig. 6 are indicative of the minimal entanglement achievable with our remote parity measurement. They show that if, e.g., 1000 measurements can be performed during a qubit lifetime (comparing 100 μs lifetime to 100 ns duration of the most demanding measurement operation) and transmission fidelity reaches about 85%, then the steady-state entanglement fidelity that is stabilized for all times can be pushed up to 98%.

V. EXPERIMENTAL CONSIDERATIONS

All the required operations for the above proposal have been individually implemented within the framework of quantum superconducting circuits, where the cavity-qubit setting is the standard minimal building block. The strong dispersive coupling of a transmon qubit to a high-Q cavity mode [32] provides universal controllability of the state of the quantum harmonic oscillator modeling the cavity mode [33,34]. This controllability has been experimentally illustrated with circuit QED setups [24,35]. Such a coupling enables us to prepare the probe field in a cat state and to perform the CNOT gates U A and U B of Eq. ( 1) between the qubits and associated intracavity fields. Recent experiments realizing a variable coupling between cavity modes and a transmission line [21][22][23] η T 1 / T d 0. 99 75 0 .9 9 5 0 .9 9 0 . 9 8 0 .9 7 0 .9 5 0 . 9 (7). Here, we have assumed

T A 1 = T B 1 = T 1 and T A,B φ = ∞.
For each value of η (transmission fidelity of the meter quantum channel for the remote parity measurement) and of T 1 /T d (qubit characteristic lifetime expressed in number of measurement iterations), we have estimated the optimal value of α by maximizing Eq. ( 8). We have then run 5000 Monte Carlo simulations with values of α up to two times larger or smaller than this estimated optimum and selected the best result. The latter turned out to be close to the computed α by optimizing Eq. ( 8). To estimate steady-state fidelity, we have averaged the fidelity obtained between t = 3T 1 and t = 6T 1 , when starting at t = 0 from |B e + . We have checked that this initial state does not affect the steady state; but starting close to the expected steady state allows us to minimize the transient effects.

provide the possibility of catching the propagating microwave field, performing the required gate between the qubit and the cavity field, and finally releasing back the cavity photons [36]. Finally, while the measurement of photon-number parity has also been realized in a similar setup [37], one can further simplify the protocol by letting U B map the intracavity field to coherent states |± √ ηα instead of the cat states |C ± √ ηα . This would allow us to replace a photon-number parity measurement by a simple homodyne detection of the released field using a parametric amplifier [20]. A single measurement duration in such an experimental realization depends mainly on the time required to perform the operations U A and U B . As explained in Refs. [33,34] this gate time is roughly the inverse of the dispersive coupling strength. Recent experiments where this coupling strength is more than three orders of magnitude larger than both the qubit and the cavity decay rates [24,35,38] indicate that high entanglement fidelities should be achievable whenever a transmission efficiency of more than 70% is reached.

VI. CONCLUSION

We have shown that it is possible to perform a near DP-QND measurement of the parity of two distant qubits despite an important loss through the transmission channel between them. By preparing the probe field in a quantum superposition of two coherent states with opposite phases, we avoid backaction of the probe losses on the target qubits' state inside any parity subspace. Indeed, such losses mainly decrease the measurement strength but barely affect its DP-QND character. Therefore, even with an inefficient transmission channel, by repeating the measurement many times one efficiently projects the joint qubits state onto a definite-parity subspace. This enables not only to deterministically and hardware-efficiently prepare an entangled state of two distant qubits but, combined with a quantum feedback strategy, it also protects such a state against decoherence channels and systematic errors. The operations required to perform this loss-tolerant parity measurement are within the reach of state-of-the-art experiments. Their implementation will lead to an important step forward for implementing quantum teleportation protocols in a loss-tolerant way, and more particularly towards the modular architecture solution for large-scale quantum information processors [4].

if they have odd parity, while the target qubits remain unaffected. Thus not knowing n makes the final probe state uncertain, as with the bit flip, implying that detection results will not allow us to perfectly discriminate the parity. But the probe ends up in the same (unknown) state for all even-parity states of the target qubit; and it ends up in another same state for all odd-parity states of the target qubit. This ensures preservation of the DP-QND property under channel losses.

APPENDIX B: KRAUS OPERATORS

Here, we summarize the calculations that lead to the Kraus operators of Eq. ( 2). We start with an initial joint state of the two qubits, the probe field, and the ancillary mode modeling the transmission loss, given by

FIG. 1 .

 1 Photon-number parity detector

  [13]. Starting with |q p = |0 p we let U A (U B ) implement a CNOT gate on |q p conditioned by |q A (|q B ); see Fig. 1(b). If the channel C was perfect [E k = Identity for k = 1,2, . . . in Fig. 1(b)], then an initial state (|11 ± |00 ) A,B |0 p / √ 2 would remain unchanged, while an initial state (|10 ± |01 ) A,B |0 p / √ 2 would become (|10 ± |01 ) A,B |1 p / √ 2 just before detection of the probe. Thus the measurement operations correspond to K+ , Kwith ξ = 1. Now let the channel subject the probe to an unknown number n ∈ {0,1,2, . . .} of bit flips E k = |1 0| p + |0 1| p ;

2 ( 1 FIG. 3 .

 213 FIG.3. The performance of the near DP-QND parity measurement is explained by two rates. The measurement strength r parity indicates how fast the state is projected onto a parity eigenspace. The dephasing rate r dephasing indicates how fast the coherence inside a parity eigenspace vanishes. Here, we have fixed η = 0.75 and |α| 2 = 2. The blue (solid) curve illustrates the maximum fidelity to one of the two Bell states |B e + and |B o + , in simulations of our measurement scheme, when the two qubits are initialized in the state (|0 + |1 ) A (|0 + |1 ) B /2 and we take an average over 1000 Monte Carlo trajectories. The black (dashed) curve illustrates the average parity converging at the characteristic rate of the Lyapunov function V parity (ρ). The red (dashed) curve, decreasing with the same characteristic rate as the coherence function C(ρ), illustrates the dephasing inside a parity eigenspace. The green (dashed) curve represents the product of both these effects. This theoretical curve fits reasonably well with the simulations (blue solid curve).

5 FIG. 5 .

 55 FIG. 5. Feedback stabilization of Bell state |B e+ .

|ψ 0 =ηα and β = √ 1 - 1 -

 011 (c 00 |00 A,B + c 11 |11 A,B + c 01 |01 A,B + c 10 |10 A,B )|C + α p |0 env . Following the definition of the unitary operators U A , U B , and U η BS in Eq. (1), just before performing the photon-number parity measurements of the probe and ancillary fields, this joint state has evolved to ηα. Detecting the probe and the ancillary fields both in even parity leads to applying the following projection operator as the measurement backaction:+,+ = Id A,B ⊗ even p ,with Id being the identity map. Since |C + β has even parity and |C - β has odd parity, the projected wave function can be simply read off the above expression, env . In a similar way, we have the following projected wave functions for other measurement outcomes: ηα env . This corresponds to the Kraus operators as defined by Eq. (2).

  FIG. 2. Average evolution, over 1000 Monte Carlo simulations, of the initial state (|0 + |1 ) A (|0 + |1 ) B /2 under repeated approximate DP-QND parity measurements (2), (3), for various probe field intensities |α| 2 and transmission efficiency η = 0.75. The state initially converges towards a definite parity state |B e fidelity to those states by dephasing, due to the slight perturbation of the DP-QND property inside definite parity subspaces. Here and everywhere throughout the paper, the fidelity to a Bell state |B e + is defined as B e + |ρ|B e + .

	+ = (|00 + 2 thanks to the measurement, √ + = (|01 + |10 )/ 2 or |B o √ |11 )/
	then slowly loses

  A (|0 + |1 ) B /2, before dephasing destroys coherence in absence of feedback. Estimate τ meas of the number of measurements after which this highest fidelity is reached (thin red, right axis), giving an indication of measurement strength.

	meas Maximum achievable precision 1-F	1	η=.7 η=.75 η=.8 η=.85 1.5	2	|α | 2	2.5	3	3.5	10 1 10 2 10 3 4 τ meas = #measurements to attain F meas T
	FIG. 4. Estimate F meas of the highest fidelity to the closest
	Bell state (thick black, left axis), obtained by repeated application
	of our near DP-QND parity measurement (2), (3) starting from
	(|0 + |1 )

ACKNOWLEDGMENT

The authors thank the Agence Nationale de la Recherche for financial support under Grant ANR-14-CE26-0018.

APPENDIX A: NECESSARY CONDITION FOR LOSS-TOLERANT DP-QND MEASUREMENTS

We prove that it is impossible to design the measurement setup of Fig. 1(a) in order to obtain a DP-QND parity measurement under arbitrary channel noise, i.e., when E k in Fig. 1(b) can be arbitrary.

(S0) Note that any unconditioned local unitaries on the qubits along the measurement chain can be merged into U A and U B . Indeed, even if physically we apply such ŨA to qubit A after U B , this ŨA commutes with U B since it acts as identity on both |q B and |ψ p . We will thus assume without loss of generality a setting with just two arbitrary unitaries U A and U B , which make the arbitrary probe interact respectively with |q A and |q B .

(S1) To keep invariant the even-parity states |00 A,B and |11 A,B , we have to take U A such that U A |0 A |ψ p = |0 A |φ 0 p and U A |1 A |ψ p = |1 A |φ 1 p for some |φ 0 , |φ 1 . Indeed consider a contrario that, e.g., an initial |0 A was mapped by U A into a state involving |1 A . Then, since no further action is applied on qubit A, at the end of the nominal measurement chain there would unavoidably be a nonzero probability to end up in a state involving this |1 A , which would contradict the DP-QND objective for an initial state |00 A,B .

(S2) If we chose |φ 1 = e iθ |φ some θ , then we could rewrite

with Z θ acting only on qubit A and defined by

This clarifies that the probe state |φ 0 p would carry no information about the state of qubit A, hence we can get no parity information by later measuring the probe system at B.

(S3a) On the other hand, if we have |φ 1 = e iθ |φ 0 , then it is always possible to write a noise channel which keeps |φ 0 in place but, conditionally on the unknown state |q env of some environment variable, moves |φ 1 . In other words, referring to Fig. 1(b),

To have a DP-QND measurement, the remaining action U B has to disentangle (|00 + |11 ) A,B / √ 2 from the probe before detection. In the case where E 1 is not applied (no channel loss), we thus need

for some |φ b . In the case where E 1 in contrast is applied by the lossy channel, we need

for some |φ a . For the experiment we must select a unique U B , without knowing the environment state |q env i.e., without knowing whether E 1 was applied or not. The left parts of Eqs. ( A1 This shows the impossibility for the measurement to be DP-QND with respect to the "conditional E 1 " channel noise. We have taken care to keep our construction fully general, such that we can conclude the following: Whatever our design for |ψ p , U A , and U B , we will obtain a measurement setup either which gains no parity information, or for which some particular noise actions E 1 can destroy the DP-QND character and perturb definite-parity states.

For example, in our proposed construction with CNOT gates:

(1) For S1 we have |φ 0 = |0 , |φ 1 = |1 ; we are obviously in the situation of S3a, not of S2.

(2) Regarding S3a, the phase-flip channel σ z for instance would keep |0 at |0 but move |1 to -|1 . Then to save the situation U B would have to satisfy in particular U B |1 B |1 p = -U B |1 B |1 p , which is not possible.

This example shows that a DP-QND parity measurement is not possible if the transmitted probe state can be subject to both bit flip and phase flip errors. In this paper, we discussed the case, motivated by realistic experimental conditions, where the transmission channel only subjects the probe to an unknown number of bit flips (probe's cat states). Consider indeed a situation where the transmission channel subjects the probe to an unknown number n of unitary operations of the same type U C , where (U C ) N equals identity for some integer N . Now let U A and U B apply conditional V gates on the probe, where V = (U C ) N/2 . Then all the operations U A , U B , U C commute with each other, and an initial probe state |ψ 0 p gets mapped just before detection onto


I Alfaro

J Yvonnet

F Chinesta

E Cueto
email: ecueto@unizar.es

A study on the performance of natural neighbour-based Galerkin methods

Keywords: meshless, natural element method, Sibson interpolation, Laplace interpolation, pseudo-NEM, natural neighbour search, computational cost

In this paper we address the problem of improving natural element simulations in terms of computational cost. Several problems are discussed, that include an efficient natural neighbour search algorithm and a comparison of different natural neighbour-based interpolation algorithms. In particular, we review the so-called pseudo-NEM, a moving least squares-like approximation scheme that employs natural neighbours, and compare it with traditional Sibson and Laplace interpolation schemes in terms of both accuracy and computational cost. Some examples in linear Elasticity and visco-plasticity are analysed in order to test the proposed schemes in engineering problems.

INTRODUCTION

Performance of a numerical method in terms of computational cost is a critical issue. Although computers are growing exponentially in terms of speed of computation and capacity of storage, so does user's needs. The case of meshless methods is particularly noteworthy. It is well-known that nowadays meshless methods greatly alleviate the user time employed in meshing the domain of interest, and significantly reduce the time employed in remeshing, that is nearly completely eliminated in practice, since the distortion of the mesh (or cloud of nodes) does not affect the accuracy of the results.

In a meshless method, instead of constructing a mesh on a given domain (thus placing nodes in appropriate positions so as to obtain well-shaped elements), usually the approach is reversed and the objective is, given a cloud of points, to obtain their connectivity. This connectivity can be sought based on a distance criterion (as in the Element-Free Galerkin Method (EFGM) [START_REF] Belytschko | Element-free Galerkin methods[END_REF] or the Reproducing Kernal Particle Method (RKPM) [START_REF] Liu | Reproducing kernel particle methods[END_REF]) or other type of neighbourhood (natural neighbourhood in the case of the Natural Element Method (NEM), for instance, see [START_REF] Sibson | A vector identity for the Dirichlet tesselation[END_REF]).

In a recent paper devoted to these topics [START_REF] Idelsohn | To mesh or not to mesh: that is the question[END_REF] it has been said that for a meshless method to be competitive, '... the evaluation of the nodal connectivity [should be] bounded in time and linear with the total number of nodes in the domain'. Despite that the constants accompanying the linear law are also important (i.e. the relative amount of time employed in generating the connectivity, despite that it varies linearly or not), it is clear that the generation of the connectivity is a crucial issue in meshless methods.

In the NEM the generation of the connectivity is done by searching for natural neighbours of the nodes. Natural neighbours of a given node are those nodes that share with it an edge of a Delaunay triangle or tetrahedron. Usually, however, natural neighbourhood is sought between an integration point and the nodes whose support cover the mentioned integration point. Thus, efficient natural neighbour search algorithm should be employed in order to make the NEM competitive. This constitutes one of the main objectives of this papers.

In a Galerkin procedure, there is another issue that should be addressed. The computation of the approximation (shape) functions and their derivatives is also important. Since the irruption of the very first meshless methods [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF] it was noted that the computation of the shape functions, and especially their derivatives, can be a very costly process. This was one of the early motivations of the so-called 'diffuse derivatives' in the work of Villon and colleagues before mentioned in the context of moving least squares (MLS) approximations.

In the NEM framework, two different interpolants have been proposed. The original version of the method [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF] employed natural neighbour (Sibson) interpolation [START_REF] Sibson | A vector identity for the Dirichlet tesselation[END_REF]. Sibson interpolation is known to be a costly procedure, although it is also well-known in the approximation community as a high-quality interpolation scheme with many advantages (see, for instance, [START_REF] Hiyoshi | A numerical comparison of the natural neighbor interpolation and the finite element method[END_REF]). Sukumar [START_REF] Sukumar | Natural neighbor Galerkin methods[END_REF] later proposed the use of Laplace (also known as non-Sibsonian) interpolation [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF]. This interpolant is considerably less costly than the original Sibson interpolant, although somewhat less smooth. Some of the advantages of natural neighbour interpolation have been attributed to the particular geometry of the shape functions' support, that remains quite isotropic even for large distortions of the original cloud of points. Based on this, Reddy et al. [START_REF] Schembri | A three-dimensional computational procedure for reproducing meshless methods and the finite element method[END_REF] and Yvonnet and Chinesta [START_REF] Yvonnet | A hybrid element-free Galerkin and natural element meshfree method for direct imposition of boundary conditions and faster three-dimensional computations[END_REF] proposed independently a MLS procedure that employs a weight function whose support is the same as the one of the NEM (we hereafter refer to this approximation as pseudo-NEM). This method runs considerably faster than traditional NEM methods (for equivalent reproducing conditions) and can constitute an interesting alternative to alleviate this high cost of the NEM.

These and other issues are discussed in this paper. Firstly, in Section 2, a review of natural neighbour interpolation is made, with special emphasis on the newly proposed pseudo-NEM approximation before mentioned.

In Section 3, we propose a efficient natural-neighbour search algorithm that runs in linear time and that has proved to be very fast, especially when using a large amount of nodes.

In Section 4, we analyse the computational cost associated to each interpolation procedure in a series of examples in linear elasticity and visco-plasticity. We discuss the relative importance of 2 this cost in linear and non-linear problems, with regular and irregular clouds of different densities. Finally, in Section 5, some conclusions are drawn.

NATURAL NEIGHBOUR INTERPOLATION

Essentially, the NEM is a Galerkin procedure that relies on natural neighbour interpolation to construct the trial and test functions characteristic of this method. As mentioned before, there exist nowadays different interpolation procedures that are based upon natural neighbourhood. In this section we review three of them, although there exist others.

Prior to this definition we introduce some basic geometrical entities that are needed for further developments.

Voronoi/Dirichlet diagrams

Consider a model composed of a cloud of points N = {n 1 , n 2 , . . . , n m } ⊂ R d , for which there is a unique decomposition of the space into regions such that each point within these regions is closer to the node to which the region is associated than to any other in the cloud. This kind of space decomposition is called a Voronoi diagram (also Dirichlet tessellation) of the cloud of points and each Voronoi cell is formally defined as (see Figure 1)

T I = {x ∈ R d : d(x, x I )<d(x, x J ) ∀J = I } (1)
where

d(•, •) is the Euclidean distance function.
The dual structure of the Voronoi diagram is the Delaunay triangulation, ‡ obtained by connecting nodes that share a common (d -1)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there being some so-called degenerate cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in 2D, as depicted in Figure 1 (right)). Another way to define the Delaunay triangulation of a set of nodes is by invoking the empty circumcircle property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called natural neighbours and hence the name of the technique.

Equivalently, the second-order Voronoi diagram of the cloud is defined as

T I J = {x ∈ R d : d(x, x I )<d(x, x J )<d(x, x K ) ∀J = I = K } (2)
Based on these definitions, different natural neighbour interpolation schemes have been proposed. We review some of the most popular.

Thiessen interpolation

The simplest of the natural neighbour-based interpolants is the so-called Thiessen's interpolant [START_REF] Thiessen | Precipitation averages for large areas[END_REF]. Its interpolating functions are defined as

I (x) = 1 if x ∈ T I 0 elsewhere
(3) ‡ Even in three-dimensional spaces, it is common to refer to the Delaunay tetrahedralization with the word triangulation in the vast majority of the literature. The Thiessen interpolant is a piece-wise constant function, defined over each Voronoi cell. It defines a method of interpolation often referred to as nearest neighbour interpolation, since a point is given a value defined by its nearest neighbour. Although it is obviously not valid for the solution of second-order partial differential equations, it can be used to interpolate the pressure in formulations arising from Hellinger-Reissner-like mixed variational principles, as proven in [START_REF] González | Volumetric locking in natural neighbour Galerkin methods[END_REF].

Sibson's interpolation

The most extended natural neighbour interpolation method, however, is the Sibson interpolant [START_REF] Sibson | A vector identity for the Dirichlet tesselation[END_REF][START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF]. Consider the introduction of the point x in the cloud of nodes. Due to this introduction, the Voronoi diagram will be altered, affecting, the Voronoi cells of the natural neighbours of x. Sibson [3] defined the natural neighbour co-ordinates of a point x with respect to one of its neighbours I as the ratio of the cell T I that is transferred to T x when adding x to the initial cloud of points to the total volume of T x . In other words, if (x) and I (x) are the Lebesgue measures of T x and T x I , respectively, the natural neighbour co-ordinates of x with respect to the node I is defined as

I (x) = I (x) (x) (4)
In Figure 2 the shape function associated to node 1 at point x may be expressed as Sibson's interpolation scheme possesses the usual reproducing properties for this class of problems, i.e. verify the partition of unity property (constant consistency), linear consistency (and therefore are suitable for the solution of second-order partial differential equation). Other interesting properties such as the Kronecker delta property [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] and linear interpolation on the boundary [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Yvonnet | A new extension of the natural element method for nonconvex and discontinuous problems: the constrained natural element method[END_REF] are also verified by the NEM.

1 (x) = A ab f e A abcd (5)

Laplace interpolation

Recently, some new interpolation schemes based on the concept of natural neighbours have been proposed [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF][START_REF] Belikov | The non-Sibsonian interpolation: a new method of interpolation of the values of a function on an arbitrary set of points[END_REF]. Its application in the context of the NEM dates back to the work of Sukumar et al. [START_REF] Sukumar | Natural neighbor Galerkin methods[END_REF]. This so-called Laplace or non-Sibsonian interpolation, has received considerable attention, since it involves magnitudes of one order less of the space dimension (i.e. the calculation of areas in three-dimensional problems, for instance, instead of volumes). If we define the cell intersection t I J = {x ∈ T I ∩ T J , J = I } (note that t I J may be an empty set) we can define the value

J (x) = |t xJ | d(x, x J ) (6)
Thus, the shape function related to node 4 at point x in Figure 3 is defined as

ns 4 (x) = 4 (x) n J =1 J (x) = s 4 (x)/ h 4 (x) n J =1 [s J (x)/ h J (x)] (7)
where s J represent the length of the Voronoi segment associated to node J and n represents the number of natural neighbours of the point under consideration, x • n = 4 in this example. Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see [START_REF] Sukumar | Natural neighbor Galerkin methods[END_REF]).

This same interpolation scheme is used in the meshless finite element method (MFEM) context [START_REF] Idelsohn | The meshless finite element method[END_REF], but in this case, the space is divided into polyhedral elements and Laplace interpolation is constructed within each polyhedron. Since the same properties mentioned for Sibson interpolation (i.e. linear completeness, exact interpolation on the nodes and linear interpolation along the boundaries) hold, a conforming method is achieved.

Pseudo-NEM interpolation

As mentioned before, one of the key issues in the high quality of the interpolation achieved by the NEM despite the large distortion of the underlying triangulation is that the support of the shape function remains quite isotropic, since it is composed of the union of circumcircles covering a given evaluation point.

Based on this analysis, Reddy and co-workers [START_REF] Schembri | A three-dimensional computational procedure for reproducing meshless methods and the finite element method[END_REF] and Yvonnet and Chinesta [START_REF] Yvonnet | A hybrid element-free Galerkin and natural element meshfree method for direct imposition of boundary conditions and faster three-dimensional computations[END_REF] proposed independently a MLS technique whose weighting function possesses NEM-like support. But in order to alleviate the burden associated with the computation of natural neighbour shape functions, they substitute it by a simplified one, composed of the union of conical volumes.

Let

u h = p T (x)a(x) ( 8
)
be a polynomial approximation to the essential field of the problem, u, with p T (x) = [1, x, y, z, . . .] a polynomial basis to a given order, and a(x) a vector of unknown coefficients. To determine a(x), the functional

J = 1 2 n i=1 w i (x)[p T (x)a(x) -u i ] 2 (9)
is minimized with respect to a. w i (x) represents a weighting function associated with the node i.

See [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF] for more details. This leads to the linear system

Aa = Bu (10)
where

A jk = n i=1 w i (x) p j (x i ) p k (x i ) ( 11
)
B i j = w i (x) p j (x) (12)
This leads to

u h (x) = p T (x)A -1 Bu (13)
and thus, we can speak of a shape function

(x) = p T (x)A -1 B (14)
The key point in this approach is the choice of the weight function, that in this case is a conical function. Let f be a cone, of unit height, whose basis matches the corresponding Delaunay circumcircle containing the evaluation point x and whose tip projection over the horizontal plane is the node n i (see Figure 4) In order to avoid the overlapping of cone functions, whereby conserving the continuity of the weight function, a cone portion is associated with each of the Delaunay triangles connected to node n i . The cone function is thus non-zero if a point x belongs to the intersection between the Delaunay circumcircle and the portion of the plane such as any point in the basis formed by the origin node n i and the vectors n i n j and n i n k has positive co-ordinates in this basis. n j and n k being the other two vertices of the triangle (see Figure 5). Due to the particular shape of its support, which is defined for any nodal distribution, this weight function guarantees interpolation conditions (w i (x j ) = i j ); as Delaunay circles passes through the nodes. Furthermore, the properties of positiveness and monotonically decrease are verified. Since the cone functions are linear between two nodes, the continuity of the weight function is guaranteed.

f = -→ n i P - -→ n i x -→ n i P (15) with --→ n i P = -2 -→ cn i • -→ n i x -→ n i x • -→ n i x -→ n i x (16)
Sibson, Laplace and pseudo-NEM shape functions for a regular lattice of 5 × 5 nodes are depicted in Figure 6.

Computation of natural neighbour shape functions

The computation of natural neighbour (Sibson or Laplace) shape functions deserves some comments. In three-dimensional settings, it is necessary to compute volumes of Voronoi cells (polyhedra) for the Sibson interpolant and areas of the corresponding Voronoi facets (poligons) for Laplace interpolants. Both schemes share many features, although expressed in different dimensions.

In two dimensions, the algorithm proposed by Watson [20] is especially efficient. Unfortunately, this algorithm seems not to be straightforwardly extended to three-dimensional settings. In our experience, the algorithm proposed by Lasserre [START_REF] Lasserre | An analytical expression and an algorithm for the volume of a convex polyhedron in R n[END_REF] for the computation of the volume of polyhedra seems to be the most appropriate and has been employed in the three-dimensional examples included in this work. For an in-depth discussion on the Watson's algorithm, the interested reader should read [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] or [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF].

We review here Lasserre's algorithm, since we consider that three-dimensional problems are both the most interesting cases and the most time-consuming. Lasserre's algorithm begins by expressing a convex polyhedron's volume (first-or second-order Voronoi cell) in the form of a set of inequalities in R d that may be written as {x|Rx b} [START_REF] Yvonnet | A new extension of the natural element method for nonconvex and discontinuous problems: the constrained natural element method[END_REF] where x represents, as usual, a point in R d , R is a matrix of dimension (m, d) and b is a column vector of dimension m (number of restrictions that define the volume). The volume enclosed by the polytope (d-dimensional counterpart of a polygon) is then given by

V (d, R, b) (18)
the ith face of the polytope given by

{x|(r i • x) = b i , Rx b} (19)
where r i represents the ith row of R. This face's volume, in R d-1 space is denoted by

V i (d -1, R, b) (20)
A traditional way of computing the volume of the polytope is given by

V (d, R, b) = 1 m i=1 d(p, H i ) × V i (d -1, R, b) (21)
d
where p is a fixed point in the space, m the minimum number of restrictions that define the polytope (no redundant restrictions are considered) and d(p, H i ) the distance from point p to the hyperplane H i given by the ith restriction that defines the volume.

This algorithm can be constructed in a recursive way, so the volume computation is performed in the form of a binary tree, beginning by dimension d and leading to the computation of some lengths in R. This volume can then be computed by

V (d, R, b) = 1 d m i=1 b i |r it | V it (d -1, R i,t , b t ) (22)
In this expression R i,t represents the reduced matrix, obtained from R by elimination of the tth variable, by means of the equation r i x = b i ; b t is the reduced vector after this elimination and r it is the tth element of r i . V it represents the volume in dimension d -1 obtained with the reduced matrix R i,t and the reduced vector b t . The formulae used to calculate the r J K and b J in terms of R i,t and b t are

r J K = r jk - r jt r it r ik (23) b J = b j - r jt r it b i (24)
In the work of Braun and Sambridge [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF] the value of t is chosen such that

|r it | = max j |r i j | (25)
Note that in this case the sum is set up to m instead of m since redundant restrictions are automatically eliminated by the algorithm.

In order to calculate the derivatives of the shape functions, Equations ( 22)-( 24) should be differentiated

*V (d, R, b) *x I = * x I V (d, R, b) = 1 d m i=1 b i • * x I V it (d -1, R i,t , b t ) |r it | + V it (d -1, R i,t , b t ) |r it | • * x I b i - b i |r it | * x I |r it | (26)
* x I r J K = * x I r jk -* x I r jt r it r ik - r jt r it * x I r ik (27) * x I b J = * x I b j -* x I r jt r it b i - r jt r it * x I b i ( 28
)
where * x I represents the derivative in the x I direction. 7) has six neighbours. So, it will be necessary to compute six areas corresponding to second-order Voronoi cells, formed by the integration point and its neighbours (highlighted in Figure 7). To this end, we compute the list of neighbours • Neighbours of node 8: 8, 7, 1, 2, 3, 4, 5, 6.

• Neighbours of the integration point: 8, 7, 1, 2, 3, 6.

• Common neighbours: 8, 7, 1, 2, 3, 6.

The equation defining the intersection of the hyperplane between two generic points A and B and the plane R 2 (or, equivalently, z = 0 in R 3 ), that contains these two points, is

(x A -x B ) T x = x A 2 -x B 2 2 (29)
So, computing the different restrictions to form the matrix R, we arrive to

Restriction x A y A x B y B r i x r iy b i 7 -8 0 1 1 1 -1 0 -0.5 1 -8 0 0 1 1 -1 -1 -1 2 -8 1 0 1 1 0 -1 -0.5 3 -8 1.5 0.
R = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x 7 -x 8 y 7 -y 8 x 1 -x 8 y 1 -y 8 x 2 -x 8 y 2 -y 8 x 3 -x 8 y 3 -y 8 x 6 -x 8 y 6 -y 8 x 8 -x p y 8 -y p ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 0 -1 -1 0 -1 0.5 -0.5 -0.4 0.8 0.2 0.2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , b = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x 2 7 + y 2 7 -x 2 8 -y 2 8 2 x 2 1 + y 2 1 -x 2 8 -y 2 8 2 x 2 2 + y 2 2 -x 2 8 -y 2 8 2 x 2 3 + y 2 3 -x 2 8 -y 2 8 2 x 2 6 + y 2 6 -x 2 8 -y 2 8 2 x 2 8 + y 2 8 -x 2 p -y 2 p 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -0.5 -1 -0.5 0.25 0.8 0.36 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (30)
This set of restrictions is shown in a graphic manner in Figure 8. Note that in this case the only row of matrices R and b in Equation ( 30) depending on the position of the evaluation point P is the last one. Thus, by deriving Equation ( 29), defining each row of the restriction set, with respect to the co-ordinates of the evaluation point, we arrive to (i = 1, . . . , m) and thus the derivatives of the restriction matrices are

* *x p [x i -x p ] = -1, * *y p [x i -x p ] = 0, * *x p [y i -y p ] = 0, * *y p [y i -y p ] = -1 (31)
* *x p x 2 i + y 2 i -x 2 p -y 2 p 2 = -x p , * *y p x 2 i + y 2 i -x 2 p -y 2 p 2 = -y p ( 32
)
dR = *R *x p *R *y p = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (33) db = *b *x p *b *y p = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 0 0 0 -0.8 -0.8 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (34)
We start with i = 1 with respect to the sum in Equation [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF]. The first step is to select the pivot and, therefore, the variable to be eliminated in the reduced matrices. In this case, with reference to Equation [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF], |r it | = max j |r 1 j | = |r 11 | = 1 and the new R d-1 space is the y-axis. We compute the x value and introduce it in the rest of equations, giving

R 1,1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ r 1 r 2 r 3 r 4 r 5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ r 22 -
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -1 -1 -0.5 0.8 0.2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , b 1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ b 1 b 2 b 3 b 4 b 5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ b 2 - r 21 r 11 b 1 b 3 - r 31 r 11 b 1 b 4 - r 41 r 11 b 1 b 5 - r 51 r 11 b 1 b 6 - r 61 r 11 b 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -0.5 -0.5 0 1 0.26 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (35) R 1,1 • y b 1 ( 36
)
The reduced derivative matrices can be calculated using Equations ( 27) and ( 28) and the data stored in dR and db matrices, i.e.

* x b 5 = * x b 6 - * x r 61 r 11 b 1 = db 61 - dr 61 r 11 b 1 = -0.8 - -1 -1 (-0.5) = -0.3 (37)
The entire matrices are

dR 1,1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ * x r 1 * y r 1 * x r 2 * y r 2 * x r 3 * y r 3 * x r 4 * y r 4 * x r 5 * y r 5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 0 0 -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , db 1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ * x b 1 * y b 1 * x b 2 * y b 2 * x b 3 * y b 3 * x b 4 * y b 4 * x b 5 * y b 5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 0 -0.3 -0.8 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (38)
y 1.3 (43)
where y 0.5 and y 1.25 are the restrictions that limit the length of the segment, 0.75 (see Figure 9). Let us call y min = b 1 /r 1 = -0.5/-1 = 0.5 and y max = b 4 /r 4 = 1/0.5 = 1.25. Then, the length on the y axis is

L = V it (d -1, R i,t , b t ) = y max -y min = 0.75 (44)
and its derivative is

* x I L = * x I V it (d -1, R i,t , b t ) = * x I y max -* x I y min ( 45
)
The value of * x I y max can be calculated using the chain rule and the dR 1,1 and db 1 matrices, i.e.

* x y max = 1 r 4 (* x b 4 -y max • * x r 4 ) = 0 (46)
In this case * x L and * y L are both zero. In other words, this segment's length (see Figure 9) does not change with a small change in the position of the evaluation point, P.

Introducing this value in the recursive formula ( 22) and ( 26) we obtain

A 8 p = A 8 p + b 1 • L d • |r 11 | = 0 + -0.5 • 0.75 2 • 1 = -0.1875 (47)
* x A 8 p = * x A 8 p + b 1 • * x L d • |r 11 | + L d • |r 11 | * x b 1 - b 1 |r 11 | * x |r 11 | = 0 + -0.5 • 0 2 • 1 + 0.75 2 • 1 0 - -0.5 1 0 = 0 (48) * y A 8 p = * y A 8 p + b 1 • * y L d • |r 11 | + L d • |r 11 | * y b 1 - b 1 |r 11 | * y |r 11 | = 0 + -0.5 • 0 2 • 1 + 0.75 2 • 1 0 - -0.5 1 0 = 0 (49)
This same procedure is then repeated for the m rows of the R matrix. The sum of the m = 6 terms in Equation ( 22) will give the total area of the second-order Voronoi cell corresponding to node 8 and evaluation point P. The shape function 8 (P) can then be obtained as

8 (P) = A 8P A T ( 50
)
where A T = nnn i=1 A i P , being nnn the number of natural neighbours of the considered point. The resulting area A 8P is depicted in Figure 8.

Finally, the derivative of the shape function is obtained by applying the chain rule to Equation (50)

* 8 (P) *x I = 1 A T *A 8P x I -8 (P) *A T *x I with I = 1, 2, . . . , d (51)

AN EFFICIENT NATURAL NEIGHBOUR-SEARCH ALGORITHM

For the computation of the weak form of the problem some kind of numerical integration of the type

f (x) d ≈ nip i=1 i f (x i ) (52)
must be employed, since there is no closed form for the natural neighbour shape functions just described, in general. In Equation (52) i represent the weights associated with the particular quadrature scheme employed, nip the number of integration points and x i the associated quadrature points.

It is well-known that the use of three-points (four in three-dimensional cases) numerical quadrature on each triangle (or, respectively, tetrahedra) leads to a numerical integration error in the NEM and it was judged as the cause of the apparent lack of conformity of the NEM with standard patch-tests. In [START_REF] Gonzalez | Numerical integration in natural neighbour Galerkin methods[END_REF] a study was performed on the influence of this error in the results and how to alleviate it. It was demonstrated that, despite the integration error due to the non-polynomial form of the shape functions, the NEM is frequently more accurate than standard, linear triangular finite element method (FEM).

Whatever method employed to perform the integration, a search for the neighbours of the integration point must be accomplished. This is equivalent to search for all the triangles whose circumcircle contains the point x. In that case the three nodes of the triangle will be natural neighbours of the integration point. The naive approach to this problem is to perform an O(m 2 ) search (being m the number of nodes in the set, as mentioned in the preceding section), by testing if the evaluation point is a natural neighbour of each node in the set. As will be shown later on, this scheme performs extremely badly, and should be employed only when using very coarse clouds of nodes.

In the original work by Braun and Sambridge [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF], the search for natural neighbours was done by employing the walking triangle algorithm, due to Lawson [START_REF] Lawson | Software for C 1 surface interpolation[END_REF]. This method is used to search for the triangle containing the evaluation point.

In our approach this search is not necessary since we employ a method that closely resembles the structure of traditional FE codes. Thus, a 'natural element' will be composed of a triangle, which defines, in general, the integration domain, and a number of nodes whose associated shape function's support covers any of the integration points. A similar definition can be given if stabilized conforming nodal integration [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] is used, for instance, instead of traditional, Gauss-based quadratures.

Recently, a new algorithm has been proposed in [START_REF] Cai | A local search algorithm for natural neighbours in the natural element method[END_REF] for the natural neighbour search. In it, the possible neighbour candidates are limited to a squared region of length 2r around the given integration point. In this way the search is done in an algorithm of order O(n • m) where n<m. Although in principle valid, this algorithm possesses a main drawback derived from the nature of natural neighbourhood itself. Two points can be neighbours even if they are far away from each other, depending on the relative position of the other nodes. This circumstance, that constitutes one of the main differences between NEM and other meshless methods, complicates the practical application of the before mentioned algorithm, so that the number of possible neighbour candidates n to be considered should be high enough. Otherwise, the risk of eliminating natural neighbours of the given point is always present, and these will be omitted, thus altering the approximation result.

We have developed a new searching algorithm that performs the search in an 'expansive' way, starting with the triangle that defines the integration points. We refer in this explanation to two-dimensional settings for simplicity, although the algorithm is extended to three-dimensional settings straightforwardly. Obviously, the Delaunay triangulation of the points must be previously computed, prior to the application of the algorithm. This is done in our case by employing the 'Detri' software [START_REF] Muecke | Shapes and implementations in three-dimensional geometry[END_REF], that has proven to be very fast and efficient. It employs a symbolic perturbation technique to avoid degenerate cases [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]. Nevertheless, this particular choice does not affect the conclusions of the work presented here.

Consider, for instance, a set of nodes whose triangulation gives 11 triangles, as shown in Figure 10. The algorithm starts with the triangle that defines the integration domain, here labelled as 1. Obviously, the three nodes defining triangle 1 are natural neighbours, by definition, of the three integration points.

We then continue (step 2, Figure 11) by moving through the edges of the starting triangle, say, to triangle number 2. The circumcircle of triangle 2 covers some of the integration points of the element, so we add its nodes to the list of neighbours. Repeat this search iteratively by moving through triangle edges, say to triangle number 3. In this case, the circumcircle does not cover any of the integration points and we can stop the search through this branch of the tree and come back to triangle number 2.

All the triangles neighbouring triangle number 2 have been now checked, so we move back to triangle number 1 and perform the search again beginning by other, non-checked, edge (we move, for instance, to triangle number 4). The process finishes when checking triangle number 7 and its associated triangles.

This iterative-search algorithm provides all element's neighbours, while keeping the number of checked triangles proportional to n, the number of nodes (or, equivalently, triangles) on the cloud. This algorithm can be programmed in the form of a tree, whose pseudo-code is summarized in Algorithm 1.

Algorithm 1

Natural neighbour search for a given natural element. triangle = initial triangle for triangles sharing an edge do if any of the integration points is in the circumcircle and not "checked" then add the three nodes to the list of neighbours then label the triangle as "checked" then triangle = next neighbouring triangle else come back to previous triangle end if end for return

The algorithm runs in approximate constant time for each triangle, since the number of checked neighbouring triangles is nearly constant for each of them. Thus, the total time for the complete neighbour search is O(n). For this to be true, it is necessary to perform an storage of triangle neighbouring a given one. This can be done in the triangulation process, thus saving computational time. We refer the reader to [START_REF] Muecke | Shapes and implementations in three-dimensional geometry[END_REF] for further details on the particular data storage algorithm.

Computer savings by employing this algorithm are shown in Figure 12. Savings reach a relative amount of 4500 times for meshes of around 20 000 tetrahedra.

NUMERICAL EXAMPLES

Two-dimensional plate with a hole problem

In order to compare the behaviour of the different approximations previously proposed, we begin by the well-known two-dimensional problem of an infinite plate with a hole under traction. Domain's geometry is shown in Figure 13.

Analytical solution for this problem exists, and it can be found in many classical books like, for instance, in [START_REF] Timoshenko | Teoría de la Elasticidad[END_REF]

u 1 (r, ) = a 8 r a ( + 1) cos + 2 a r ((1 + ) cos + cos 3 ) -2 a 3 r 3 cos 3 (53) u 2 (r, ) = a 8 r a ( -3) sin + 2 a r ((1 -) sin + sin 3 ) -2 a 3 r 3 sin 3 (54)
where

= 3 -4 (55) = 3 - 1 + (56)
is the Kolosov coefficient for plane strain and stress, respectively. In this case, we have assumed E = 1.0 and = 0.25. represents the shear modulus and a is the radius of the circular hole. The traction applied at infinity, 0 , was taken unitary. By imposing appropriate symmetry conditions we model only one quarter of the plate. Analytical tractions are imposed on the boundary of the domain so as to simulate the infinite length of the plate. This problem has been analysed by employing Sibson, Laplace, pseudo-NEM interpolations, and FEM. Comparative results for the convergence rates in L 2 and H 1 norms are shown in Figures 14(a) and(b).

As expected, all methods possess similar convergence rates, indicated by R. If we analyse the computational cost of each simulation, we conclude that FEM simulations are the fastest ones, without significant differences between all the methods, see Figure 15.

Two-dimensional results are not very representative, in our opinion, since there exist competitive algorithms for the computation of natural neighbour shape functions (i.e. the before mentioned Watson's algorithm [START_REF] Watson | An Implementation of Natural Neighbor Interpolation[END_REF]). Pseudo-NEM does not appear to be very competitive, but in threedimensional problems the situation changes drastically, as will be seen in the following section.

Three-dimensional linear elastic compression test

To test the performance of the methods in three-dimensional settings, we have chosen a very simple problem of a cubic block under compression, composed of a linear, elastic material. The schematic representation of this problem is shown in Figure 16. The exact solution of this problem is spanned by each of the four methods, so the analysis of the accuracy will not be performed. Instead, we have focused on the computation time employed in solving the problem, by checking different clouds, both regularly and irregularly distributed over the unit cube. Examples of such clouds are shown in Figure 17. We checked both regular and irregular meshes since the behaviour of the proposed neighbour-search algorithm could differ in its speed of computation. For regular meshes, the search process is somewhat faster.

E,ν

The computational costs generated with this problem are analysed in Figures 18 and19 for regular and irregular meshes, respectively. We analyse both the total time employed in the calculation, as well as the part of this time employed in the computation of the shape function and its derivatives.

As expected, the use of Sibson interpolation becomes prohibitive in three-dimensional cases. But Laplace and pseudo-NEM approaches, despite being more costly than FEM, appear to be competitive. Note that, although employing more computation time, the absence of the need for remeshing is not considered here. Especially noteworthy is the fact that nearly 90% of the computation time is employed in the computation of the shape function when using Sibson interpolation.

Results for irregular meshes are similar. Again, Sibson interpolation becomes prohibitive, while Laplace and pseudo-NEM approaches result to be competitive. While FEM employs a negligible part of the calculation in computing the shape function, Sibson interpolation constitutes up to 90% cent of the computing time.

FEM is a competitive method for this kind of problems as has been proved in the last decades. It has revealed to be among the most used numerical techniques both in industry and academia, for obvious reasons. In our opinion, meshless methods should not be used for problems in which FEM has proven to be a very robust technique. The difficulties associated with FEM begin in problems with large transformations, where remeshing is necessary, especially when transfer of internal variables is necessary from the old meshes to the new ones. This kind of problems is analysed in the following section.

Extrusion of an aluminium profile

Direct aluminium extrusion is one of the most extended forming processes worldwide. Roughly speaking, it is a process used to produce long profiles by pressing a billet of hot aluminium through a hole with a certain shape. A schematic representation of the process is given in Figure 20. Obviously, very large strains are present throughout the process. Numerical simulation of extrusion has been studied largely by the FEM, by employing Arbitrary Lagrangian Eulerian (ALE) or Eulerian approaches. The interested reader is submitted to the early references by Zienkiewicz and co-workers for examples on this topic [START_REF] Zienkiewicz | Flow of solids during forming and extrusion: some aspects of numerical solutions[END_REF][START_REF] Zienkiewicz | Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes[END_REF].

If we neglect elastic deformations when compared to plastic ones, the flow of hot aluminium can be assimilated to that of a non-Newtonian fluid. This approach has received the name of 'flow formulation' [START_REF] Zienkiewicz | Plastic flow in metal forming. (I) Coupled thermal (II) Thin sheet forming. Applications of Numerical Methods to Forming Processes[END_REF]. This is the approach followed here. For an in-depth description of aluminium behaviour and its simulation with NEM, the reader is referred to [START_REF] Alfaro | Three-dimensional simulation of aluminium extrusion by the alpha-shape based natural element method[END_REF].

We consider the balance of momentum equations, without inertia and mass terms

∇ • r = 0 (57)
and the assumed incompressibility of a von Mises-like flow where d represents the strain rate tensor (deviatoric part of the velocity gradient), y represents aluminium's yield stress and d represents the so-called 'effective strain rate' (second invariant of the strain rate tensor). The dependence of the yield stress on the strain rate is given by the well-known Sellars-Tegart equation

∇ • v = 0 (
y (d) = S m arcsinh ⎡ ⎣ d 1 A e Q/RT 1/m ⎤ ⎦ (60)
with d 1 = max{d, d 0 } and d 0 an initial threshold in the strain rate. Parameters governing the aluminium yield stress are summarized in Table II.

With this model we performed the simulation of the extrusion of a cross-shaped profile, whose geometry is shown in Figure 21. Only one-half of the geometry was analysed, taking advantage of the symmetry of the domain. The model was composed of 3021 nodes, 15 843 tetrahedra and 63 372 integration points (four per tetrahedra).

Obviously, large deformations appear during the simulation. For comparison purposes, we ran the simulation using Laplace interpolation. At a given time step, we stopped the simulation and performed the following time step using the four different approximation techniques, namely, finite element, Sibson, Laplace and pseudo-NEM approximations. The geometry of the domain at this time step is shown in Figure 22.

One of the main results of the comparison between the different approximation techniques is the differences found in the results. In Figure 23 a comparison is made of the small area marked in Figure 22. Note the highly distorted mesh appearing in this time step. This produces spurious high levels of strain rate in the FEM solution-in the zone highlighted by a white circle-, far from the extrusion die, due undoubtedly to the distortion of the tetrahedra. At this step, FEM mesh would have necessitated of a remeshing process, with the well-known inherent diffusion of internal variables during the field transfer of internal variables. While Sibson interpolation predicts slightly higher levels strain, Sibson, Laplace and P-NEM interpolations have demonstrated to provide good results when compared to experiments (see [START_REF] Filice | Cross extrusion with asymmetric die: a first preliminary comparison between finite element and meshless formulation predictions[END_REF]). The results seem to be insensible to mesh distortion. Note, however, that near the symmetry planes, the number of natural neighbours for a given tetrahedron is low. In these circumstances, Sibson, Laplace, P-NEM and FEM approximations are very similar (note that Sibson tetrahedron with only four neighbours is a CST tetrahedron, see [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]). This could be the reason for the high level of equivalent strain rate appearing in the four results near the symmetry plane.

The total amount of time employed in the simulations of the before-mentioned time step are summarized in Table III.

These results are even more notorious if we perform the same simulation employing a purely Newtonian behaviour. These results are outlined in Figure 24. In this case the spurious high levels of strain rate present in the FEM simulation are even higher.

It is worth noting that the pseudo-NEM approach seems to present locking in this case, since virtually no strains are predicted in this step. Maybe by employing richer consistency approximations in the MLS adjustment of the approximation would provide locking-free results, but there is a limited number of nodes in the support of the shape functions and this is not possible everywhere with the formulation presented before.

The stiffness matrices in the NEM are banded and sparse, as in the FEM, but the number of neighbours of a given 'element' is usually much higher than in the FEM. Thus, calculation and inversion of the tangent matrix takes around 450 s with the three NEM methods, while in FEM it takes around 80 s. Again, Laplace approximation seems to provide a good compromise between accuracy and speed of computation, with no lack of accuracy detected due to mesh distortion. Is in this type of problems where the authors believe that meshless in general, approaches. Note also that in the before presented results, remeshing time is not considered for FEM results.

CONCLUSIONS

In this paper we have presented some techniques to improve the computational performance of natural neighbour Galerkin methods. In particular, we have developed a new neighbour-search algorithm that proves to run in O(n) time with big efficiency. An in-deep study of the performance of the method has been accomplished, both from the accuracy point of view and from the computational cost one. The P-NEM approach has revealed some very interesting features, especially its efficient computation times, but should be employed, while further analysis is performed, to problems where the LBB condition is not relevant. While Sibson approximation seems to be prohibitive in three-dimensional settings, Laplace interpolation within a Galerkin NEM approach has proven to give an excellent compromise between computational costs (obviously higher than those of the FEM) and efficiency and accuracy, since no need of remeshing is needed.

Thus, the diffusion of internal variables (like those appearing in plasticity, for instance) during field transfers between meshes is avoided. Despite the highly distorted meshes employed in our experiments, no need of remeshing was encountered.

Figure 1 .Figure 2 .

 12 Figure 1. Delaunay triangulation and Voronoi diagram of a cloud of points.

Figure 3 .

 3 Figure 3. Definition of non-Sibsonian co-ordinates.

Figure 4 .

 4 Figure 4. Definition of the conical weighting function f .

Figure 5 .

 5 Figure 5. Zone associated with a cone portion.

Figure 6 .

 6 Figure 6. (a) Sibson; (b) Laplace; and (c) pseudo-NEM shape functions.

Figure 8 .

 8 Figure 8. Graphical representation of the restriction set.

Figure 9 .

 9 Figure 9. Final length of the first segment in the computation.

Figure 10 .

 10 Figure 10. Step 1 in the search algorithm. Integration points are represented by circles.

Figure 11 .

 11 Figure 11. Step 2 in the search algorithm. We move to triangle number 2 and check the neighbourhood.

Figure 12 .Figure 13 .

 1213 Figure 12. Time employed in the neighbour search: (a) for a traditional, O(n 2 ) algorithm and (b) by the proposed O(n) algorithm.

Figure 14 .

 14 Figure 14. Convergence of the results for the different methods in the plate with a hole problem: (a) L 2 -norm error and (b) H 1 -norm error. R denotes the convergence rate found for the different methods.

Figure 15 .

 15 Figure 15. (a) Time employed in the computation of the shape functions and their derivatives and (b) total computation time for the resolution of the plate with a hole problem.

Figure 16 .

 16 Figure 16. Geometry of the compression test.

Figure 17 .

 17 Figure 17. Different clouds of nodes employed in the simulation of the block under compression.

Figure 18 .

 18 Figure 18. (a) Total time employed in the compression test and (b) part of this time employed in the computation of the shape function and its derivatives. Regular clouds.

Figure 19 .

 19 Figure 19. (a) Total time employed in the compression test and (b) part of this time employed in the computation of the shape function and its derivatives. Irregular clouds.

Figure 20 .

 20 Figure 20. Schematic representation of the extrusion process.

Figure 21 .

 21 Figure 21. Geometry of the die for the extrusion of the cross-shaped profile: (a) perspective view and (b) geometry.

Figure 22 .

 22 Figure 22. Geometry of the domain for the extrusion simulation test. The ellipse highlights the zone detailed in Figure 23.

Figure 23 .

 23 Figure 23. Contour plot of the second invariant of the strain rate tensor at the location indicated in Figure 22: (a) FEM; (b) Sibson; (c) Laplace; and (d) pseudo-NEM results.

Figure 24 .

 24 Figure 24. Contour plot of the second invariant of the strain rate tensor at the location indicated in Figure 22. Newtonian behaviour: (a) FEM; (b) Sibson; (c) Laplace; and (d) pseudo-NEM results.

Table I .

 I Nodes for the calculation of Sibson shape function. Configuration of the point set given in TableIand position of the integration point. The areas needed to perform the computation are highlighted.

	Node	x	y
	1	0	0
	2	1	0
	3	1 . 5	0 . 5
	4	2	1
	5	1 . 4	1 . 6
	6	0 . 6	1 . 8
	7	0	1
	8	1	1
	Figure 7.

2.6.1. Example: computation of shape function through Lasserre's algorithm. In order to clarify the previous algorithm, consider the set of points in two dimensions given in Table

I

. We compute Sibson's shape function associated with node 8 at the point (0.8, 0.8).

The evaluation point (labelled P in Figure

Table II .

 II Material parameters for AA6063 Aluminium alloy.

	Parameter	Units	Value
	d 0	(s -1 )	0 .005
	S m	(MPa)	25
	m		2.0
	A	(s -1 )	6 × 10 9
	Q	(J/mol)	1.4 × 10 5
	R	(J/mol K)	8.314

Table III .

 III Total amount of time (seconds) employed for the simulation of one time step in the extrusion problem. NEM in particular, approximations can be competitive with FEM simulations. When large transformations in a Lagrangian framework are present in a problem, possibly with internal variables related to history, NEM seems to provide an attractive alternative to more traditional

		d			d
		1.8		1.8
		1.7		1.7
		1.6		1.6
		1.5		1.5
		1.4		1.4
		1.3		1.3
		1.2		1.2
		1.1		1.1
		1			1
		0.9		0.9
		0.8		0.8
		0.7		0.7
		0.6		0.6
		0.5		0.5
		0.4		0.4
		0.3		0.3
		0.2		0.2
		0.1		0.1
	(a)		(b)
		d			d
		1.8		1.8
		1.7		1.7
		1.6		1.6
		1.5		1.5
		1.4		1.4
		1.3		1.3
		1.2		1.2
		1.1		1.1
		1			1
		0.9		0.9
		0.8		0.8
		0.7		0.7
		0.6		0.6
		0.5		0.5
		0.4		0.4
		0.3		0.3
		0.2		0.2
		0.1		0.1
	(c)		(d)
		Number	Calculated
		of iterations	shape functs.	Total time	Percentage
	FEM	13	0.21	1040.01	0.02019
	NEM-Sibson	15	15 339.81	21 979.48	67.79
	NEM-Laplace	12	70.3	5580.03	1.260
	pseudo-NEM	10	10.99	4674.96	0.2351

and

Contract/grant sponsor: Spanish Ministry of Education and Science; contract/grant number: CICYT-DPI2005-08727-C02-01