J Viquerat

Fitting experimental dispersion data with a simulated annealing method for nano-optics applications

A considerable amount of materials in nanophotonics are dispersive, enabling the propagation of so called surface plasmons at their interfaces with dielectrics. Hence, a reliable fit of frequency-dependent permittivity functions with an appropriate model is a first-order necessity for the accurate design of nano-optics devices with time-domain numerical methods, such as finite-difference time-domain (FDTD) or discontinuous Galerkin time-domain (DGTD). In this paper, we present the necessary ingredients to fit experimental permittivity functions using the simulated annealing (SA) method with a generalized second-order dispersion model implemented in the Diogenes software suite (see https://diogenes. inria.fr/). By scanning through different classes of materials, we come up with effective rules of thumb to make the fitting process fast and accurate.

Introduction

Over the course of its development, the nanophotonics field has seen a tremendous amount of devices designed for applications ranging from biosensors [CLS + 11] to nanostructured photovoltaic cells [START_REF] Atwater | Plasmonics for improved photovoltaic devices[END_REF]. A vast majority of these take advantage of the combination of a specific geometry with selected materials to obtain a desired functionality. Through the years, the design process of these devices has started relying more and more on numerical methods, either in frequency or in time-domain [START_REF] Gallinet | Numerical methods for nanophotonics: standard problems and future challenges[END_REF]. Although it is relatively easy to take dispersive materials into account with frequency-domain methods, with time-domain methods it is necessary to fit experimental data with an adapted frequency-dependent function that can be transformed back in the temporal domain by inverse Fourier transform [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF], and discretized with the chosen numerical method. In section 2, we shortly review existing dispersion models and dispersive materials, and motivate our choice of fitting algorithm. Then, in section 3, the SA technique is outlined. Section 4 represents the core contribution of this paper. We first detail the modifications introduced to improve our fitting results as well as the settings used for the algorithm, before conducting a large panel of fitting experiments on a variety of dispersive materials. Rules of thumb for accurate and efficient material fittings with the SA method are deduced. Eventually, we test the obtained fits in the context of the Diogenes DGTD solver.

Dispersive materials and dispersion models

The development of dispersion models for metals originates in the work of Drude [START_REF] Drude | Zur elektronentheorie der metalle[END_REF], who devised the so-called Drude model based on the kinetic theory of gases:

ε r,Drude (ω) = ε ∞ - ω 2 d ω 2 + iωγ d .
Although this model presents a good adequation with experimental data of noble metals for low visible frequencies, it does not account for interband transitions that arise for certains metals in the high visible frequency range. For this reason, it is common to add one or more Lorentz poles to the Drude model, thus yielding the following form:

ε r,Drude-Lorentz (ω) = ε ∞ - ω 2 d ω 2 + iωγ d - ∆εω 2 l ω 2 -ω 2 l + iωγ l .
This method shows a good agreement with experimental data for noble metals, and has been widely implemented, often with four or more Lorentz poles [START_REF] Hao | Efficient dielectric function for FDTD simulation of the optical properties of silver and gold nanoparticles[END_REF]. However, this model fails to properly fit materials such as transition or post-transition metals, because of the electronic correlation existing in such materials, introducing a retardation effect [START_REF] Wolff | Simple magneto-optic transition metal models for time-domain simulations[END_REF]. During the last decade, advanced dispersion models were designed to tackle dispersive materials more efficiently. Among them, the Critical Points (CP) model [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method[END_REF] and the Complex-Conjugate Pole-Residue Pairs (CCPRP) model [START_REF] Han | Model dispersive media in FDTD method with complexconjugate pole-residue pairs[END_REF] have shown great improvements over the standard Drude-Lorentz model in fitting metals and compounds on broad frequency ranges with a limited amount of poles. As briefly sketched in [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method[END_REF], these two models are mathematically equivalent at the continuous level. While both of these models use coupled first-order poles with complex coefficients, in [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF] we introduced a generalized model exploiting real coefficients only:

ε r (ω) = ε ∞ - np l=1 c l -iωd l ω 2 -e l + iωf l , (1)
where n p is the number of second-order poles. It is easy to show that this model is also equivalent, in the continuous domain, to CP and CCPRP models. In [START_REF] Lanteri | Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics[END_REF], we showed that model (1) is causal if e l ≥ 0 and f l ≥ 0. Additionally, it is easy to obtain a permittivity function with an imaginary part identically equal to zero by imposing d l ≡ and f l ≡ 0.

In figure 1, we regroup a representative sample of dispersive materials regularly used in nano-optics devices, including noble, transition and post-transition metals, but also semiconductors, pure or compounded with metals or other semiconductors. All the experimental data comes from [Pol]. Although these different materials show very different behaviors, we will show that a permittivity function such as (1) allows to fit all of them with sufficient accuracy. In order to fit the free parameters of (1) to experimental data, various techniques can be used. Apart from the well-known least square method, vector fitting techniques (see [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF]) are also well developed for the CCPRP formulation. In [START_REF] Sehmi | Optimizing the Drude-Lorentz model for material permittivity: Method, program, and examples for gold, silver, and copper[END_REF], the authors present a mixed method for fitting CP parameters, using an analytic minimization along with a gradient descent method. For an increasing number of poles, one can be left with a large optimization problem presenting many local minima, likely to trap optimization algorithms away from the desired global minimum. Simulated Annealing (SA) methods have proved to be particularly efficient in finding global minima in these situations, even when the initial guess is far from the optimal solution ([KGV83]). In the following section, we shortly present the SA technique and its possibilities.

The simulated annealing algorithm

The simulated annealing (SA) algorithm was initially proposed by Kirkpatrick et al [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], and is recognized to be a suitable optimization method in the context of a cost function with multiple local minima.

Although it was initially designed to solve combinatorial problems (such as the well-known traveling salesman problem), it can also be exploited to find global minima of continuous functions [START_REF] Locatelli | Simulated annealing algorithms for continuous global optimization[END_REF]. Given a cost function f defined on a continuous domain X, the goal of an optimization algorithm will be to find x * ∈ X such as:

f * ≡ f (x * ) = min x∈X f (x).
The name of the SA algorithm comes from a metallurgic process, in which the controlled cooling of a molten metal allows a final state of lower energy than that obtained if the cooling process is left unconstrained.

As a consequence, one of the main variable of the algorithm is the temperature T , which will control the probability to move from one state to another. In essence, the SA algorithm works as follows: Additionally, we define the following quantities: [Pol] and is selected in the angular frequency range [1, 10] PHz when available.

Algorithm 1 Simulated annealing 1: Set an initial temperature T ← T 0 and choose a final temperature T f < T 0 2: Generate a random initial solution x ← x 0 3: Compute its cost f (x 0 )

4: while T > T f do 5:

for v ← 1, n v do
Loop on a given number of neighbor states 6:

Compute a neighbor solution x n = n(x)

7: Compute its cost f (x n ) 8: if f (x n ) < f (x) then 9:
x ← x n and f (x) ← f (x n ) Accept better solution 10:

else 11:

Compute normal random number r 12:

Compute acceptance a(f (x), f (x n ), T ) T ← c(T ) Decrease temperature 19: end while N M the total number of moves at each temperature step; N A the total number of accepted moves at each temperature step; N AU the total number of accepted uphill moves at each temperature step; N D the total number of denied moves at each temperature step; N CM the total number of cumulated moves.

13: if a(f (x), f (x n ), T ) >
As can be seen, various parameters can be tuned in the SA algorithm, namely:

The initial and final temperature, respectively T 0 and T f ; The number of visited neighbors at each temperature n v ; The initial guess x 0 ; The acceptance function a; The neighbor function n; The cooling function c.

Since the SA algorithm is a well-established method, there is a vast literature on the topic of making wise choices for these parameters and functions. For this reason, we will not expand too long on all the existing possibilities (the reader is referred to the references that follow, and the references therein).

The choice of the initial and final temperatures depend on the considered problem. A typical choice for T 0 is usually 1.0, but we shall see later that in some cases it is interesting to start from higher values. The choice of the final temperature may seem a little arbitrary, and it is possible to replace it with a stopping criterion based on the amount of recently-accepted better solutions at current temperature; Similarly, the number of visited neighbors at each temperature is an arbitrary value. However, it is reasonable to say that the algorithm will perform better if it is allowed to sample a larger amount of states at each temperature, at expense of the computational cost;

The acceptance function, although usually simple, represents the core part of the algorithm: its role is to determine if an uphill solution (i.e. a neighbor solution which cost is higher than current solution) must be accepted or not. The most spread choice of acceptance function is the Metropolis formula [START_REF] Henderson | Theory and practice of simulated annealing[END_REF] :

a(f (x), f (x n ), T ) = min 1, exp f (x) -f (x n ) T . (2)
If the neighbor solution is better, then obviously a = 1. On the other hand, if f (x n ) > f (x), then 0 < a < 1, and this solution will be randomly accepted. Given that T is a decreasing function of the number of iterations of the algorithm, the probability of an uphill solution being accepted will decrease as the algorithm converges toward a minimum (either global if it has succeeded, or local if it failed to escape from it);

The goal of the neighbor function is to build a neighbor state x n from current state x. In discrete problems, this is often simply done by randomly swapping two elements of the solution (in the traveling salesman's problem, two cities are randomly picked and swapped in the ordered list). In continuous problems, however, building neighbors introduce an additional source of freedom for the user, since the neighborhood of a state must be defined. In practice, the loop on the neighbor states can be embedded in a loop over the total number of parameters to optimize. Then, a neighbor state is obtained by randomly moving the current free parameter in a defined range. The size of the neighborhood need not to be fixed throughout the optimization process, several refinements being available in the literature [START_REF] Xinchao | Simulated annealing algorithm with adaptive neighborhood[END_REF]. In the following, we denote by ∆ n the normalized size of the neighborhood of a state. ∆ n = 1 indicates that the whole parameter space can be reached from any state in it, while ∆ n = 0.5 indicates that the neighborhood of a given state is limited to half the size of the parameter space around it;

The choice of the cooling function is the topic of an abundant literature. One of the most common choice might be the exponential schedule [START_REF] Nourani | A comparison of simulated annealing cooling strategies[END_REF]:

T n+1 = αT n ,
where 0 < α < 1 is a constant factor. Typical values range from 0.8 to 0.99, directly impacting the performance and the computational cost of the algorithm.

The variety of choices for each algorithm parameter must be adapted to the considered problem. In the next section, we will focus on the problem of fitting experimental permittivity data to existing models containing up to a few tens of free parameters. 4 Fitting dispersive materials

Cost function

For the following fitting experiments, we define the cost function:

f (x) = Ns j=1 ε r exp (ω j ) -ε r fit (ω j , x) ∆ε r exp + ε i exp (ω j ) -ε i fit (ω j , x) ∆ε i exp , (3)
where N s is the total number of samples in the experimental data (assumed to be much larger than the total number of free parameters of the model), the • r and • i exponents respectively represent the real and imaginary parts, while ∆ε r exp and ∆ε i exp respectively correspond to the range of values of the real and imaginary parts of the experimental permittivity.

By weighting the distance between the fit and the experimental data by 1 ∆ε r exp and 1 ∆ε i exp , we ensure that the algorithm grants a similar importance to the real and the imaginary parts of the permittivity function.

Indeed, in the case where the figure of merit ε i ε r is globally low or high over the frequency range, the absolute distance ε r exp (ω j ) -ε r fit (ω j ) + ε i exp (ω j ) -ε i fit (ω j ) may not be representative of the overall quality of the fit. To illustrate this, we compare in figure 2 the fits obtained with 1 second-order pole for Al2O3, with and without such a weighting. Obviously, the weighting significantly improves the quality of the fit, providing an acceptable solution with a single pole.

Parameters range

A typical range of angular frequencies in plasmonics is ω ∈ [1, 10] PHz. Hence, to reduce the parameter space to be scanned by the algorithm, we work with a normalized angular frequency ω = ω × 10 -15 . Given the form of (1), one would expect c l and e l to be roughly in [0, 100] PHz 2 , while d l and f l would be in [0, 10] PHz. Hence, we rewrite (1) as:

ε r (ω) = ε ∞ - np l=1 cl 2 -iω dl ω2 -ẽl 2 + iω fl , (4)
with the following relations: By doing so, we scale all parameters to the same range of values, thus restricting the parameter space for the SA algorithm to the form [1,

ε max ] × [0, p max ] 4np .
In the remaining of this paper, we set ε max = 10 and p max = 10.

Incorporating constraints to the cost function

In the context of this study, we wish to incorporate parameter constraints to the cost function. In order to respect the parameter range determined above, we use a simple penalization of the form:

f ρ (x) = f (x) + ρ max (0, -x + u, x -v) .
(

) 5
where

u = t (1, 0, ..., 0) , v = t (ε max , p max , ..., p max ) ,
and ρ is a well-chosen penalization parameter. For the current problem, we observed that ρ = 1000 or higher provided satisfying results. In practice, this penalization prevents ε ∞ to be lower than 1, and the parameters of (4) to be negative [START_REF] Gill | Practical Optimization[END_REF].

Sign of ε i fit

As an additional refinement, we modify the cost function in order to discriminate negative imaginary parts.

To do so, at each evaluation of f (x), we evaluate if ε i fit (ω) < 0 for all ω in the range of interest. If so, the cost function is increased in the fashion of previous section so as to rule out current state x. This is particularly important in cases where ε i fit is close to zero on a wide frequency range. As an example, we compare in figure 3 fits of TiO2 with and without this restriction. Although the unrestricted fit provides a better approximation of the experimental data, the negative imaginary part of the permittivity would yield an unstable behavior in any time-domain numerical method. As can be seen, the restricted fit respects this condition by decreasing more slowly to zero.

Non-absorbing materials

In some cases, experimental data can be found where the imaginary part is identically zero. As we shall see, it is enough to set a lower and an upper bound equal to zero for all d l and f l parameters, thus reducing the model to:

ε r (ω) = ε ∞ - np l=1 cl 2 ω2 -ẽl 2 .
Here, one easily recognizes the Sellmeier's law [START_REF] Gooch | Sellmeier Equation[END_REF]. Although it does not perform as well as model (4), it can be very useful in cases where even a very low imaginary part is not acceptable. As an example, figure 4 shows fits obtained with and without this feature for Si3N4 data where the imaginary part is identically zero.

Initial guess and stopping criterion

After performing a vast amount of fittings on various kinds of materials, we observed that in most cases the initial guess x 0 only had a minor influence on the final result. For this reason, in the remaining of this paper, we set x 0 ≡ 0. Regarding the influence of the final temperature T f , we replace it with a stopping criterion: the algorithm exits when no move has been accepted during the last temperature step. Unless stated otherwise, we set the size of the neighborhood to be the full parameter space, and the acceptance function equal to (2). Finally, the number of visited neighbors at each temperature step is chosen equal to 1000 × (4n p + 1), i.e. the free parameters are moved one by one, 1000 trials being made for each of them.

Noble metals

In this section, we assess the influence of the initial temperature T 0 and the decreasing temperature factor α on the obtained f * for noble metals (gold, silver and copper). Results are presented in figure 5. For slow cooling (typically, α = 0.99), one observes that the best cost is remarkably stable with respect to T 0 , with a slight minimum between T 0 = 0.1 and T 0 = 1. As α is decreased, the behavior becomes more erratic, and no clear trend can be deduced. This effect is called quenching, and is known to cause suboptimal results with the SA technique [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. In figure 6, we illustrate the difference in behavior for the algorithm between a slow cooling schedule (α = 0.99), and a quenching situation (α = 0.8) by plotting f * both as a function of T and N CM processed by the algorithm, with T 0 = 1000. As can be seen, with a low value of α the optimization process gets stuck in a local minimum and fails to escape from it, although the stopping criterion is met at a much lower temperature than for a high value of α. This is however less visible for gold and copper than for silver. From these first experiments, we select the couple (α = 0.99, T 0 = 0.1) as a good starting candidate for future fits. Some additional comments on the algorithm insights can be made, especially on the acceptance function and the choice of the neighborhood size. On figure 8(a), we plot N A , N AU and N D as functions of T for gold experimental data. As can be seen, roughly half of the accepted moves are uphill ones, underlining the local-minima-escaping feature of the SA algorithm. Broadly speaking, it is usually recommended to have between 50 and 60% of uphill accepted moves to obtain satisfying results. The best fits obtained with four poles for noble metals are shown in figure 7.

Acceptance function

Regarding the acceptance function, the Metropolis formula (2) is not the only available choice. Another usual option is the Barker function [START_REF] Schuur | Classification of acceptance criteria for the simulated annealing algorithm[END_REF]:

a(f (x), f (x n ), T ) = 1 1 + exp -f (x)-f (xn) T . (6)
With this criterion, some downhill steps might be rejected if they do not substantially improve the solution.

As T gets lower, downhill moves are more likely accepted, while uphill moves get more and more rejected. This is illustrated in figure 8(b), where we can also see that the 50% rule of thumb is respected. To test this acceptance function further, we reproduce the graph of figure 5 for α = 0.99 (computed with Metropolis function) with the Barker function (results are displayed on figure 9). Although we observe some slightly better results for T 0 = 0.01, using the Barker formula as acceptance function brings no major improvement in the final result of the algorithm, and we shall use the Metropolis function in the remaining of this paper.

Neighborhood size

Regarding the size of the neighborhood, no clear strategy seems to emerge from the literature [START_REF] Henderson | Theory and practice of simulated annealing[END_REF].

In figure 10, we assess the impact of the normalized neighborhood size ∆ n on f * for noble metals, with α = 0.99 and T 0 = 0.1. Again, it is hard to draw a clear trend from these results, although we observe no significant variations in the final fit quality over the ∆ n range. For that reason, we will use ∆ n = 1 in the remaining of the paper. 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 2 10 3 0 0.5 Although we can note some differences, the Barker acceptance function yields no clear improvement over the Metropolis one. Please notice that the scale on the y axis is different from that of figure 5.

1 1.5 2 •10 4 T # of moves N A N AU N D (

Transition and post-transition metals

In this section, we perform the same test as before with varying T 0 for transition and post-transition metals (see figure 11). For α = 0.99, the results for transition metals offer more contrast than that of noble metals, with good fits obtained for T 0 = 0.1. However, we see that using a slower cooling schedule with α = 0.999 provides more robust results with regard to the variations of T 0 . It also improves the best fit obtained in the T 0 ∈ [0.1, 1] region. For post-transition metals, results are stable for any value of T 0 above 10 -2 , and significant improvement is observed for higher value of α. The best fits obtained with four poles for transition and post-transition metals are shown in figure 12.

Semiconductors

As previously, we assess the influence of T 0 and α on the fits of semiconductors. As could be expected from figure 1, the permittivity functions of such materials are difficult to fit. Hence, they benefit from slower cooling schedules, as can be seen in figure 13. As for other materials, T 0 = 0.1 is a robust value across all materials. Below, the algorithm provides poor solutions, especially for The best fits obtained with four poles for semiconductors are shown in figure 14.

Computational performances

The SA algorithm requires a large amount of evaluation of the cost function. Each evaluation of the latter implies to loop over all experimental data, and for each data point, evaluate the permittivity provided by the current state x. As a result, the performance of the fitting procedure strongly depends on:

The amount of experimental data points N s provided by the user;

The number of poles n p required by the user;

The number of neighbor states n v .

For reasonable values of N s , the algorithm can be performed on any desktop machine, with running times under a few minutes for standard values of n p . As an example, we show in figure 15 the time required to fit gold (N s = 49), zinc (N s = 148) with different values of n p on a standard desktop machine with a 3 GHz Intel Core i7 processor.

Numerical results

As a test of usability and stability of the obtained fits in time-domain computations, we consider the case of a doubly-periodic slab composed of a given material embedded in air. We successively test the highest stable time-step ∆t mat found with the Diogenes DGTD solver, and compare it with ∆t vac , i.e. the maximal stable time-step for the same configuration without material. In table 1, we report the ratio ∆tmat ∆tvac for the best four-pole fits obtained with all materials presented in this paper in the ω ∈ [1, 10] PHz range. This ratio represents the constraint imposed by the CFL condition inherent to any time-domain numerical method. An unreasonably low ratio would indicate that the set of parameters produced by the fitting procedure is unsuitable for time-domain computations. The restriction to four poles represents a trade-off between the fit accuracy and the memory and CPU overhead induced by each pole in time-domain computations [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method[END_REF]. As can be seen, most fits provide reasonable time-steps. The case of TiO2, for which the obtained time-step is the lowest, can be explained by the large portion of curve for which ε i ≡ 0 (roughly for ω < 4 PHz), which represents a severe constraint on the fitting algorithm.

Conclusions

Our implementation of the simulated annealing technique for the fitting of dispersive material properties has proven to be efficient, providing reliable sets of parameters while never requiring more than a couple of minutes. Among the numerous intrinsic variables of the algorithm, the cooling factor α and the initial temperature T 0 are the prime parameters to tune for materials involved in most nano-optics problems.

The parameter set (α = 0.99, T 0 = 0.1) provided good results for noble and post-transition metals. For transition metals and semi-conductors, we observed that the set (α = 0.999, T 0 = 0.1), although more computationally-intensive, represents a robust choice. Eventually, we provided an evidence for the usability of the fits produced by our fitting method in any timedomain numerical solver by computing the largest stable time-step obtained with our best four-pole fit for a wide range of materials. Even better fits can be obtained by increasing the number of poles, although at the expense of the performance of the fitting algorithm and the time-domain solver. A possible next step in the improvement of the fitting method would be to incorporate the optimization of the time-step value directly in the cost function of the SA method. As of now, our SA fitting algorithm is available in the Diogenes software package (see https://diogenes.inria.fr/).

Figure 1 |

 1 Figure1| Incomplete classification of dispersive materials based on their ε i = f (εr) behavior. All data comes from[Pol] and is selected in the angular frequency range [1, 10] PHz when available.

ω = 10

 10 -15 × ω, cl = 10 -15 × √ c l , dl = 10 -15 × d l , ẽl = 10 -15 × √ e l , fl = 10 -15 × f l .

Figure 3 |Figure 4 |

 34 Figure3| Fit of TiO2 with 3 second-order poles with and without restriction on the sign of ε i fit .

Figure 5

 5 Figure5| of T 0 on f * noble metals with np = 4, x 0 = 0 and variable values of α. For values α close to 1, we observe a remarkable stability of f * with respect to T 0 . When α decreases, the behavior is wilder, and no clear trend can be deduced.

  Figure6| * as a function of the temperature for gold, with α = 0.99 and α = 0.8. On the left figure, we observe the effect of quenching, leading to a poor local minimum despite a very low final temperature. On the right, we see that, because of the fast cooling, very few states are investigated for α = 0.8, leading to the poor solution.

Figure 7 |

 7 Figure7| Best fits obtained with four poles for gold, silver and copper.

Figure 8 |Figure |

 8| Figure 8 | A , N and N D as functions of T for gold, with α = 0.99 and T 0 = 1000. For both acceptance functions, the rule of thumb N AU 0.5 × N A is respected.

Figure 10 |

 10 Figure10| Influence of ∆n on f * for noble metals with np = 4, x 0 = 0, α = 0.99 and T 0 = 0.1. It is difficult to draw a clear trend from these results, although we observe no significant variations.

  Figure11| Influence of T 0 on f * for transition and post-transition metals with np = 4, x 0 = 0 and several values of α. Results for transition metals (left) offer more contrast than for noble ones, with good overall results for T 0 = 0.1 and α = 0.999. For post-transition metals (right), results are remarkably stable for T 0 values above 10 -2 . Notice the different scales of the y axis.

Figure 12 |Figure 13 |Figure 14 |

 121314 Figure12| Best fits obtained with four for nickel, iron, cobalt, aluminium, zinc and lead.

Figure 15 |

 15 Figure 15 | Time required to fit experimental data with variable number of poles for gold and zinc.

Table 1 |

 1 Best ∆t mat ∆tvac ratios obtained with our DGTD solver with four poles for the materials presented in this paper.

	Material	∆t mat ∆tvac
	Vacuum	1
	Dielectric (εr = 4)	0.5
	Ag	0.306
	Au	0.248
	Cu	0.467
	Al	0.421
	Zn	0.791
	Pb	0.307
	Ni	0.467
	Co	0.245
	Fe	0.713
	AlGaAs	0.312
	GaAs	0.276
	Si	0.317
	TiO2	0.0865


Jonathan Viquerat
email: jonathan.viquerat@inria.fr.

Efficient time-domain numerical analysis of waveguides with tailored wideband pulses

Keywords: Discontinuous Galerkin Time-Domain (DGTD) methods, optical waveguides, mode dispersion, wide-band source

A simple procedure for the generation of accurate polychromatic sources in waveguides taking into account mode dispersion is presented. It allows for an efficient use of timedomain solvers in the analysis of guided modes, and can be used wether the mode dispersion and field distribution are known analytically or numerically. This method is implemented in the DIOGENES Discontinuous Galerkin Time-Domain (DGTD) solver (http://diogenes.inria.fr), and is validated on a waveguide mode converter.

I. INTRODUCTION

T HE numerical analysis of optical devices with time- domain methods usually allows the use of polychromatic sources (dipole, plane wave, etc), thus characterizing the structure response in a single simulation. For waveguides, however, the existence of modes brings the task of building wide-band sources to a higher level of complexity [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]. More precisely, two main obstacles must be overcome: (i) for a given mode, the effective index n eff depends on ω, and (ii) the frequency-domain 2D mode field distribution Ê2D also depends on ω.

The bootstrapping method, proposed in [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF], is an effective yet computationally-expensive way to overcome these problems. An erroneous source is injected and propagated in an extended version of the considered waveguide. Along the propagation, the exact mode establishes in the waveguide, while the remaining energy is scattered towards the surrounding PML regions. At the edge of the waveguide, a probe plane is used to record passing fields, providing an accurate source for future computations.

In [START_REF] Yamauchi | Pulse excitation scheme for the analysis of wide-band characteristics of optical waveguides[END_REF], the authors propose another workaround, in which a wide-band source made from the sum of several polychromatic pulses. Each of them is centered around a given frequency and holds an appropriate 2D field distribution associated to this frequency. In this method, the sampling of the 2D field distribution and that of the effective mode index are concurrently made at the selected frequencies, which are most often limited to a handful.

In this paper, we propose a straightforward polychromatic source formulation for waveguide analysis based on an inverse Fourier transform of the frequency-domain solution, incorporating both the variations of the 2D field distribution and that of the mode index. The knowledge of both these quantities can be either analytic or numeric (i.e. computed in a third-party solver or read from a file). Additionally, we describe the way to impose of such a source through a Total-Field/Scattered-Field (TF/SF) interface in the DGTD framework, as well as how to compute figures of merit such as transmission, overlap or mode conversion efficiency. The accuracy of this approach is first tested on a textbook case, before being applied to a waveguide mode converter.

II. DESCRIPTION OF THE SOURCE

A. Source construction

Assuming that Ê2D and n eff are known, we remind that the corresponding monochromatic time-domain source propagating in direction z + at angular frequency ω 0 is obtained as follows:

E(t) = Ê2D (ω 0 )e -jβ(ω0)z e jω0t ,
where β(ω) is the mode propagation constant:

β(ω) = ωn eff (ω) c 0 . (1)
For plane wave analysis, polychromatic signals are most often obtained by modulating a sine or a cosine function at angular frequency ω 0 with the derivative of a gaussian bump centered on ω 0 and with well-suited variance parameter σ. Then, taking the Inverse Fourier Transform (IFT) of the frequency-domain signal yields the desired time-domain pulse. One can follow the same procedure here, taking into account the effective index of the mode:

E(t) = +∞ -∞ Ê2D (ω)σe -(ω-ω 0 ) 2 σ 2 2 e -j(β(ω)z-ωt) (2)
Most often, the latter IFT is not analytic and must be numerically computed using an Inverse Discrete Fourier Transform (IDFT):

φ(t) nω k=1 φ(ω k )e jω k t ∆ω, (3)
where ∆ω is the chosen frequency step:

ω k = ω min + (k -1)∆ω,
and n ω is the number of frequency samples:

n ω = ω max -ω min ∆ω + 1.

B. Sampling and the reciprocal Shannon theorem

The choice of n ω (or equivalently of ∆ω) will have a major impact on the accuracy and efficiency of the procedure. Indeed, if it is chosen too high, the frequency-domain expression under the integral in (3) will be over-sampled, and its evaluation will be uselessly expensive. On the opposite, if it is chosen too low, the obtained time-domain signal will suffer from aliasing, resulting in overlapping incident wave trains. In order to select an adequate value for n ω , we make use of the reciprocal Shannon theorem [START_REF] Khanyan | Sampling theorem in frequency domain for the infinite spectrum[END_REF], stating that the sampling rate ∆ω used to compute (3) should meet the following criterion:

∆ω ≤ 2π T ,
where T is the total duration of the transformed time-domain signal. Although in most cases, T might not be known, in our case it can be well approximated thanks to the Gaussian envelope of the frequency-domain expression (2). To do so, we introduce a parameter , such that fields with amplitudes smaller than are considered to be zero. Then, we seek the width of the gaussian enveloppe for which the signal level is equal to . A few lines of calculus yield:

∆ω = 2π σ 8 ln( 1 )
.

In the following, we choose = 1 × 10 -3 , and therefore:

∆ω 0.845 σ . (5)
C. Total-Field/Scattered-Field injection plane

The discretization of the electric and magnetic fields in the DGTD framework is naturally discontinuous, meaning that the unknowns are duplicated on each face joining two cells of the mesh [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF]. The consistency of the approximation is re-established through the imposition of a numerical flux (E * , H * ) at the interface, which handles the information exchange between cells (the reader is referred to [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF] and [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nanooptics with a high-order discontinuous Galerkin time-domain method[END_REF] for a thorough development of the DGTD method). As a consequence, it is possible to define separate areas in the mesh, in which the numerical unknowns either hold the total fields (E tot , H tot ) or the scattered fields (E sca , H sca ). The interface between these regions is the so-called TF/SF interface. At this interface, the fluxes are modified in order to inject the incident field in the TF region, and to subtract it in the SF region (see figure 1). With this method:

• It is possible to access both the TF and the SF fields in their regions of interest simultaneously; • Incident fields can be imposed virtually anywhere in the mesh; • Accessing both the TF and SF fields at the same place (i.e. on the TF/SF interface) allows for the computation of multiple figures of merit with a single computation; • In the present case, if a mode is properly imposed and propagated through the TF region of an unperturbed waveguide, the fields in the SF regions should be identically zero, making it a reliable method to test the source described in the last section.

Esca

Etot

TF/SF E * + E*,inc E * -E*,inc
Fig. 1. Field injection with the TF/SF method. In the TF zone (in gray), the incident field is added to the numerical flux, whereas in the SF zone, it is substracted.

D. Figures of merit

A brief word must be said about the figures of merit that will be used in this paper, namely the transmission T , the overlap O M and the mode conversion efficiency η M .. These quantities are based on the computation of the averaged Poynting vector π :

π (x, ω) = 1 2 E (x, ω) × H * (x, ω) . (6)
In expression ( 6), the frequency-domain fields are evaluated on control surfaces (here, the TF/SF surfaces) with an onthe-fly Fourier transform during the time-domain computation.

Thanks to the discontinuous nature of the DG field discretization, these quantities can be evaluated either on the total-field or the scattered-field side of the TF/SF faces, allowing to compute the transmission as:

T (ω) = So π tot • n Si π inc • n , (7)
where π tot is computed on the total-field side of the exiting TF/SF surface S o , π inc is computed using the incident field only on the total-field side of the inputting TF/SF surface S i , and n is the unit normal. As such, the transmission computes the normalized transmitted power from the input to the output plane, ranging from 0 to 1. The overlap, also ranging from 0 to 1, measures the similarity of the transmitted fields E and H * to that of a given mode E M and H * M :

O M (ω) =   Π So E, H * M Π So E M , H * Π So E, H *   Π Si E M , H * M , (8)
with: Finally, the mode conversion efficiency, which computes the normalized transmitted power that has been coupled to a given mode M , is easily defined as:

Π S E, H * = S E × H * • n.
η M (ω) = T (ω)O M (ω).

III. NUMERICAL EXPERIMENTS

A. Straight rectangular waveguide

We consider a rectangular waveguide of lateral dimensions a = 10.16 mm and b = 22.86 mm with PEC (Perfect Electric Conductor) walls. The goal of this section is to propagate a TE 01 mode in it, and to assess the efficiency and accuracy of the procedure described above. The operating frequency range is [8, 10] GHz. We exploit a high-order DGTD method [START_REF] Viquerat | Simulation of electromagnetic waves propagation in nanooptics with a high-order discontinuous Galerkin time-domain method[END_REF] on an adequately-refined tetrahedral mesh. The incident field is imposed on a Total-Field/Scattered-Field interface, separating a total field region (in the center of the numerical domain) from a scattered field region (before the input and after the output). The mesh, shown in figure 2, is terminated on both ends with a Perfectly-Matched Layer (PML) region to absorb outgoing field.

As a first experiment, we consider the propagation of an uncorrected polychromatic pulse, where frequency dependence of Ê2D and β is not taken into account. Both the quantities are approximated by their evaluation at the central frequency f c = 9 GHz. To evaluate the quality of the numerical treatment, we simply consider the evaluation of the E y component of the electric field at a chosen probe point in the SF region. Thanks to the use of TF/SF interfaces, for an unperturbed waveguide the field should be identically equal to 0 in the SF regions. In figure 3, we show a 3D snapshot of the field during computation: the presence of field in the right SF region indicates a mismatch between the numerical solution and the expected exact solution. Such a result was to be expected, since most frequencies in the pulse were numerically transported with an erroneous propagation constant.

We now consider the ω dependence of the quantities involved. The effective index of the mode is plotted in figure 4, along with the resulting spectral power of the incident pulse.

In figure 5, we conduct the same experiment as before, but using expression (2) for the omega-dependent quantities. Using expression [START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF] gives n ω = 24 for this computation. A 3D snapshot of the computed electric field is shown in figure 5. For a better comparison, we plot in figure 6 a comparative time-domain plot of E y in the right SF region, •10 -24 f (GHz) Spectral power (a.u.)

Uncorrected approach Corrected approach (b) Spectral power (arbitrary units) of the incident pulse in the regular and corrected cases Fig. 4. TE 01 mode index spectral power of the incident pulse in the uncorrected and corrected approaches.

with and without the IDFT correction. As can be seen, taking into account the ω dependance of the incident field and index greatly improves the mode mismatch. The maximal L 2 errors on the electric field observed in the right SF region, with and without mode correction, are shown in table I, along with the corresponding CPU times. The computed error is 60 times lower with the mode correction, while the CPU times is increased by 34%, thus confirming the interest of this method.

B. TE 01 -TE 02 mode converter

We now consider a TE mode converter, as presented in [START_REF] Zhang | design and analysis for TE20-TE10 rectangular waveguide mode converters[END_REF]. The converter is composed of a bent rectangular PEC waveguide of dimensions a = 22.15 mm and b = 47.55 mm. The bending is tailored for a theoretical power conversion of t (s)

E y (V/m)
Uncorrected approach Corrected approach Fig. 6. Time-domain plots of the Ey field in the right SF region with a corrected and an uncorrected pulse. We observe a very small residual error with the corrected approach, due to the inaccuracies in the Inverse Discrete Fourier Transform (IDFT).

100% at f = 8.5 GHz. The goal of this section is to assess the converter efficiency in the [7.5, 8.5] GHz frequency range in a single run, using the technique presented above. The effective indices of the two modes in the selected frequency range are shown in figure 7. As it shall be useful later, we also plot in the latter figure the effective index of the TE 10 , which is the only other mode which effective index is real in the frequency range of interest. Following the results of the previous section, n ω is set to 24 for this experiment. The mode conversion efficiency from TE 01 to TE 02 is computed following the details of section II-D, and is shown in figure 8 for both pulses. As expected, a peak efficiency of 1 is found for f = 8.5 GHz in both cases, while it decays for lower frequencies. However, the uncorrected approach overestimates the conversion rate from roughly 2.5 % at the minimal frequency. The obtained results are in line with what was found in [START_REF] Zhang | design and analysis for TE20-TE10 rectangular waveguide mode converters[END_REF]. Similarly to the previous section, the CPU time is increased by roughly 30.5% when using the mode-corrected source. Snapshots illustrating the mode conversion are shown in figure 9. The TE 01 mode pulse is incoming from the left, before being converted in a TE 02 that leaves the computational domain from the right (full mesh is not shown here). To complete the study, we compute the conversion efficiencies η TE01 and η TE10 , i.e. the normalized amount of energy coupled from the TE 01 mode at the input to the TE 01 and TE 10 modes at the output. Results, presented in figure 10, show that for lower frequencies, a non-negligible part of the pulse is not converted to the TE 02 mode. On the other hand, the conversion rate to TE 10 mode is identically 0. It should be noted that a very low level of backscattering was observed in the computations, hence it is not plotted here.

IV. CONCLUSION In this paper, we provided a simple expression for the computation of accurate polychromatic sources in waveguides. Details about its practical use were also discussed, such as the minimal frequency-domain sampling or the use of TF/SF faces to impose it in the computational domain. The computation of relevant figures of merit was also sketched. The described procedure was then validated on a microwaveregime rectangular waveguide, on which significantly low error levels were observed. Finally, the method was applied to a reallife mode converter, allowing to assess the device response on a large bandwidth in a single computation.

Fig. 2 .

 2 Fig. 2. Mesh of the rectangular waveguide. Light green faces correspond to PEC boundaries, while yellow ones are for ABC. The pink inner faces correspond to the TF/SF interfaces. The TF cells (not shown) are between the two TF/SF interfaces, while the SF cells (not shown) are on the left and the right of them. The PML cells terminate the mesh on the edges (in dark green).

Fig. 3 .

 3 Fig. 3. Snapshot of the Ey field during computation for an uncorrected polychromatic pulse. The presence of field in the right SF region indicates a mismatch between the numerical solution and the expected exact solution. The field discontinuity at the interface between TF and SF regions is a natural consequence of the TF/SF injection method. Field values are scaled in the range [-1, 1] V/m.

Fig. 5 .

 5 Fig. 5. Snapshot of the Ey field during computation with a corrected polychromatic pulse. The absence of field in the right SF region indicates that the transported numerical solution and the expected exact solution match. Field values are scaled in the range [-1, 1] V/m.

Fig. 9 .

 9 Fig. 9. Successive Ex field snapshots during the conversion of a polychromatic TE 01 mode to a TE 02 mode. Full mesh is not shown here.

Fig. 10 .

 10 Fig. 10. Computed mode conversion efficiency from TE 01 at the input to TE 01 , TE 02 and TE 10 modes at the output.

TABLE I MAXIMAL

 I L 2 ERRORS AND CPU TIMES WITH AND WITHOUT MODE CORRECTION. TE 01 , TE 02 and TE 10 modes effective indices. The selected frequency range is indicated in light blue. Computed mode conversion efficiency from TE 01 to TE 02 in the selected frequency range with the uncorrected and the corrected pulses. As expected, the same conversion efficiency is observed for f = 8.5 GHz. The mode mismatch at the lowest frequency induces a 2.5% discrepancy in the computed conversion rate.

				1
				0.8
			eff (f )	0.6
			n	0.4
				0.2					TE01 TE02
									TE10
				0	6	7	8	9	10	11	12
							f (GHz)
			0.85 0.9 0.95 1 Fig. 7. 7.6 η TE02 (f ) 0.8		7.8	8	8.2 Uncorrected approach 8.4 Corrected approach
							f
			Fig. 8.
		Non-corrected	Corrected
	CPU time (s.)	101	135
	Max. error	1.71 × 10 -1	2.83 × 10 -3


Jean-Alexandre A Long

Vincent Daanen

Alexandre Moreau-Gaudry

Jocelyne Troccaz

Jean-Jacques J Rambeaud

Jean-Luc L Descotes

Prostate biopsies guided by 3-dimension real-time (4D) transrectal ultrasound on a phantom. Comparative study versus 2D transrectal ultrasound

    Prostate biopsies guided by three-dimensional real-time (4-D) transrectal ultrasonography on a phantom: comparative study versus two-dimensional transrectal ultrasound-guided biopsies.

Introduction

The aim of this study is to assess the value of 3 dimensions and real-time (4D mode) ultrasonography in ultrasound-guided prostate biopsies.

We want to determine whether the ability to visualise simultaneously orthogonal planes allows a more accurate and more reproducible sampling, and whether the biopsies are more evenly distributed within the prostate volume.

Our study compares traditional 2D TRUS guided biopsy with 3D real-time (4D) TRUS guided biopsy on a synthetic model (phantom).

The objective of the study is not to assess a biopsy protocol, but rather the effectiveness of each guiding method (2D and 4D) in carrying out a specific protocol.

2D TRUS was considered a promising technique for detecting non-palpable tumor lesions. It turned out that the presence of a hypoechogenic nodule in the peripheral zone did not show sufficient sensitivity and, more importantly, specificity.

In 1989, the team from Stanford was the first to show that sextant biopsies improved a better cancer detection rate than biopsies focusing on ultrasound anomalies [START_REF] Hodge | Random systematic versus directed ultrasound guided transrectal core biopsies of the prostate[END_REF] . Despite the improvement in detection rate, it was shown that cancer was detected in 20% of the cases studied, during a second set following a negative sextant biopsy. This tended to show that the ideal biopsy scheme remained to be found and that it was difficult to determine the features of an ideal protocol [START_REF] Hodge | Random systematic versus directed ultrasound guided transrectal core biopsies of the prostate[END_REF] . This observation was modelled by Chen, who created a statistic model allowing biopsy protocol simulation, using 607 cancer specimens on computer-reconstructed radical prostatectomy specimens. This method showed that 27% of cancer specimens exceeding 0,5 cc would not have been detected using the Stanford protocol [START_REF] Chen | Comparison of prostate biopsy schemes by computer simulation[END_REF] .

Many biopsy schemes have been described to improve the sensitivity of biopsies while limiting the number of samples. The trend is now toward an increase in the number of biopsies and sampling the peripheral zone. Additional morbidity resulting from extensive protocols remains controversial [START_REF] Ravery | Diagnostic value of 10 systematic TRUS-guided prostate biopsies[END_REF] . Three dimensional ultrasonography is commercially available for routine clinical use and has been used on some devices since 1994. Rather than producing an anatomic section, as is the case with traditional 2D ultrasonography, image acquisition is performed as a volume of data with nearly immediate reconstruction and simultaneous display of sectional anatomy in three orthogonal planes--sagittal plane, transverse or coronal plane, or any arbitrary oblique. As a result 3D US allows unrestricted access to an infinite number of viewing planes. The scanning time varies from 3 to 10 seconds, depending on the volume and quality desired. It is also possible to use endocavity probes, in which case the scanning can either be rotational (around the probe's axis) or angular, which is the case with the device tested (General Electric Voluson 730 pro ®).

A 3D-volume acquisition of an organ can also be performed, such that the images can be manipulated and the coordinates of all volume points stored in the memory, with the related grey levels. 3D ultrasound appears to be helpful in predicting extracapsular extension [START_REF] Garg | Staging of prostate cancer using 3-dimensional transrectal ultrasound images: a pilot study[END_REF] .

The latest developments in ultrasound have led to real-time 3D acquisition (4D-mode). It consists of displaying, simultaneously, 3 orthogonal viewing planes.

This guiding mode has proven its effectiveness in terms of accuracy in biopsies on breast tumours [START_REF] Sauer | Ultrasound-guided large-core needle biopsies of breast lesions: analysis of 962 cases to determine the number of samples for reliable tumour classification[END_REF] , liver tumours [START_REF] Won | Value of four-dimensional ultrasonography in ultrasonographically guided biopsy of hepatic masses[END_REF] and pancreatic pseudocysts [START_REF] Polakow | A new diagnostic approach to pancreatic pseudocyst fine-needle puncture: three-dimensional sonography[END_REF] . We wish to study its benefits in prostate biopsy.

Methods

The device used is a VOLUSON® 730 Pro (General Electric) equipped with an RIC5-9 (5.9 MHz) endorectal probe, which can be used in 2D, 3D or 4D mode. The biopsy needle was Tru-Cut 18G, cutting length of 23 mm. It appears important to specify that the benefits of 4D ultrasound reside not in high-accuracy guiding toward a visible zone that is materialised through ultrasound, but rather in better sampling of the gland. Our purpose is not to assess the protocol, but rather the effectiveness in guiding the biopsy toward the zones theoretically set out by the protocol. Testing the device on a synthetic model appeared a necessary preliminary step.

Creation of a synthetic prostate model (phantom)

In a square plastic box (10 cm x 10 cm), a circular orifice spanning 6 cm in diameter was cut, so that a probe could be inserted. Inside the orifice, a cylinder intended to simulate the rectal wall was placed. The ultrasound propagation environment was made of Aquasonic® ultrasound gel (750 ml per phantom). The prostate was made by squeezing candle gel into a silicon mould. The model was the equivalent of a 40-ml prostate. Three landmarks were placed on the convex side so as to more effectively registrate the phantoms with manual methods (figure 1).

After biopsy execution, air was brought in the phantom by the needle. This made it possible to assess the biopsy path using ultrasonographic volume acquisition (3D). The same phantom model was designed for each biopsy (figure 2).

Modelling the biopsy path

The air incorporated into the phantom when the needle entered left a hyperechogenic trail that was easily detectable by ultrasonography (figure 3). A 3D volume acquisition of the entire phantom in which the biopsy was taken was performed.

In three consecutive matching tests, it was possible to superpose the trail left by the air on that of the needle.

Biopsy scheme

The protocol used was a 12-biopsy scheme, in which 2 biopsies (paramedian and lateral) were performed at the apex, in the middle and at the base of each prostate lobe.

Considering the bias that would have been created by the trails from the previous biopsies, it was essential that one phantom be produced per biopsy.

For that reason, 12 phantoms were produced, making it possible to perform a single uninterrupted session. The phantoms were replaced after each session and 336 phantoms were produced.

Each operator performed the first biopsy session using the VOLUSON 730® with 2D mode, as usual, then with 4D mode.

Following each biopsy, an ultrasonographic volume was acquired and recorded.

Processing the ultrasonographic data acquired

Producing a reference prostate.

The first step in processing the biopsy data was to create a segmented reference volume (template), using the image processing software, Analyze®. The reference volume was produced during a 3D-acquisition using standard settings for all of the acquisitions in the trial (17.7 cm-15Hz). Following segmentation, the cloud of points was exported to a software developed specially for our manipulations on Visual C++®, making it possible to visualise the reference prostate volume.

Registration of biopsies with the reference prostate (common reference base)

Volume acquisition was performed on each phantom containing a single biopsy.

The volumes were matched using the manual volume registration function on Analyze®.

Attempts at automatic registration did not yield any results, due to the significant noise and significant variations in interface contrast between the various phantoms.

The biopsy trail in the phantom was manually segmented using Analyze® in the transformed image.

The twelve biopsies carried out in one session were collected and shown in the prostate reference (figures 4 and 5). The methodology is summarised in figure 6.

.

Assessing biopsy results

A quality biopsy assessment protocol was established. Three parameters were studied:

-Entry zone: Based on the anatomical and protocol schemes for the 12 biopsies, a point located on the prostate surface was chosen as the ideal insertion point for the needle in each biopsy. Considering that there was no clear recommendation or obvious target, an entry zone spanning 7 mm in diameter (sphere) around the said puncture point was deemed suitable.

Beyond that zone the error distance was measured. The 7-mm diameter was chosen as it would allow the entry zones to be juxtaposed on the phantom surface without overlapping.

-Target zone: The biopsy tip was also compared to an ideal target zone of 7-mm in diameter. The error was defined as the distance separating the biopsy tip and the sphere (figure 7). -Biopsy distribution in the volume.

Assessing the appropriate biopsy distribution within the prostate volume was not easy.

According to Epstein, cancer is significant above 0.5 ml [START_REF] Epstein | Nonpalpable stage T1c prostate cancer: prediction of insignificant disease using free/total prostate specific antigen levels and needle biopsy findings[END_REF] . We therefore tried to find a method that would make it possible to show the volume "explored" by the biopsies (depending on biopsy length in the prostate and, above all, overlapping between the biopsies). The theoretical hypothesis was that 2 parallel biopsies 10 mm apart could not miss a significant cancer specimen (the diameter of a 0.5-cc ball is 10 mm). We therefore hypothesized that a biopsy "explored" a volume within a 1.2 mm-diameter cylinder (18 gauge needle diameter), in addition to which there was a "security zone" of 10-mm diameter.

Figure 8 shows the oversampling of a zone, in which 2 cylindrical volumes are covered by bringing together 2 biopsies. As a result, the adjacent regions of the prostate are undersampled. In other words, the biopsy-explored prostate volume was equal to the length of the biopsy needle multiplied by the surface of a 5.6 mm-radius disk.

Each voxel (ultrasonographic volume unit) in the biopsy cylinder was allocated a value of 1 or 0; those within the biopsy cylinder were given a value equal to one. The voxels that overlapped with the adjacent biopsy cylinders were counted only once.

We opted for a ratio: the number of voxels explored in the prostate/number of voxels that should have been explored biopsied if the biopsies did not overlap. This method indirectly reflected the distribution, by identifying the ability of each operator to allocate the biopsies so that they would not overlap.

During each session, it was possible to determine the volume identified by the series within the prostate.

Statistical analysis

Descriptive and inferential statistical analysis of the data was performed using STATVIEW 5.0 ® software, in conjunction with the Grenoble Technological Innovative Centre (TIC). The quantitative variables were described by their mean and standard deviation (Mean ± SD) with median and interquartile range. Means comparisons were performed by using suitable parametric tests once the application conditions where checked. When the said conditions were not checked, suitable non-parametric tests were used. As usual, tests results were considered statistically significant for a p-value < 0.05. Accuracy analysis was performed on each 2D and 4D TRUS guided biopsy by analysing distances at the ideal entry zones and ideal theoretical target zones (paired t-test). The differences of distances at the entry zones (and then at the ideal target zones) between the 2D and 4D approaches were analysed taking account of location factors in the sagittal (apex, middle and base) and transversal (paramedian, lateral) sections (One factor ANalysis Of VAriance). Biopsy reproducibility analysis was performed by looking at the deviation from the mean (paired t-test). A non parametric test (Mann-Whitney test) was performed to determine the statistical influence of each biopsy guidance mode on the percentage of the prostate covered by only 1 biopsy.

The standard biopsy protocol was explained to each operator using a summary scheme. The 2D procedure was carried out first. The second session was held within 3 hours to 3 weeks.

Results

In total, 336 biopsies were performed by 14 operators on 12 phantoms per session. The paired series tests involved 158 paired biopsies.

An increase in accuracy in hitting the target zone was identified significantly in 4D as compared to 2D (p=0.0042). In contrast, no significant improvement in biopsy accuracy was found for the entry zone (table 1 and2). Furthermore, the means of differences of distance between 2D and 4D approach, in accuracy in hitting target zone (and then at the entry zone), were not statistically different at each prostate location, as shown by one factor analysis (table 3 and4).

Discussion

We have observed that, on a phantom model, biopsy accuracy on the zones to be sampled improved with 4D guidance. This illustrates the fact that multi-dimensional visualisation improves the operator's ability to hit the virtual target. The improved reproducibility in execution shows the robustness of the guidance. In contrast, no benefits were found when 3D real time reconstruction was used as regards a possible increase in volume sampled by the set and in biopsy dispersion. This would have been a powerful claim for improving the sensitivity of a series of biopsies. Subjectively, it would appear that the 4D guiding technique is more intuitive and more beneficial for less-experienced operators, but no tie has been demonstrated. A number of biases were found throughout the work methodology. The limits inherent in the phantom material do not exactly replicate the tactile conditions of performing a biopsy. More importantly, the registration (matching volumes) introduced a bias in assessing the position of the points segmented in the biopsy. We saw that only biopsy was performed per phantom, to prevent the operator from being able to see his previous biopsies and benefiting from their guidance. For that reason, at each session, 12 registrations had to be performed. As a result, the risk of bad matching was 12 times greater. Lastly, manual registration does not allow for a perfectly reproducible method. However, those biases should cancel each other out in that they are present whatever the guidance method.

The various automatic registrations techniques developed have not, for the time being, yielded significant results. A registration algorithm is still being developed.

It is not easy to find high-quality criteria for a series of prostate biopsies: the purpose of the biopsies is to detect the tumours. The current TRUS guided biopsy technique involves randomising biopsies inside the prostate. In most cases, there is no visible target. We set out two objectives: to be as accurate as possible in producing a selected biopsy protocol and to scatter the biopsies so that they would not overlap. The ability to produce the most accurate protocol possible does not, then, assess the ability of the biopsies to detect cancer, but rather the ability to perform a selected protocol.

A more judicious criterion would be to assess the allocation of samples in the prostate volume, meaning to determine the ability to spread the biopsies across the prostate, without performing a biopsy on an already-sampled area. We assumed that a biopsy "sampled" a volume equal to the inside of a 5.6 mm-radius cylinder. This volume reflects clinical realitythat of a cancer specimen of significant volume, as currently defined by Epstein (volume of 0.5 cc, Gleason score > 6). Our distribution assessment technique explores coverage, but does not actually look at how the biopsies spread across the volume (figure 9). Both A and B biopsies achieved a maximum non-coverage rate, but ignore a cancer specimen of significant size. Consequently, a biopsy session should be judged first on accuracy in placing the needle where planned and secondly the ability to not have the needles enter through the same points. This studied showed that the first point was improved under 4D guidance, but not the second.

Regarding accuracy, the results significantly prove that concurrent visualisation of several viewing planes improves the spatial representation of prostate biopsy. In contrast, the low number of operators (14) makes interpreting the results a thorny task as regards biopsy distribution analysis. No significant difference was shown between 2D and 4D guidance.

The Prospects for 3D Ultrasound in Diagnosing Prostate Cancer

-

Memorising biopsy paths

The same approach to visualising biopsy paths could also apply to living organisms in order to establish a tumour map inside the prostate. The difference in living organisms lies mainly in how to assess the needle path. The trail left by the needle from the haematoma is a haphazard indication. It is therefore mandatory that volume acquisition be performed with the needle in place. Registrating the various biopsy needles would make it possible to display in a common reference base the path of all biopsies that might be segmented manually or more automatically.

The clinical benefit of establishing a biopsy mapping is clear in 2 situations: in a histologically suspicious zone around which additional samples may be taken, or repeated biopsies in response to a high or growing PSA when the first series was negative. In those cases, there would be no use in targeting previously biopsied zones, but rather to turn to other areas.

-Real-time guiding of prostate biopsies:

TIMC Laboratory is developing a biopsy guiding system to a pre-procedure planning. The planning could focus on an optimised protocol [START_REF] Chen | Optimization of prostate biopsy strategy using computer based analysis[END_REF] or on areas previously sampled by a first set of biopsies. A suspicious zone requiring a repeat biopsy could become the target for a guidance system or a special schedule could be applied (biopsies of the transition zone or biopsies of non-sampled zones during a first set, or interspersing the biopsy path between previous biopsies). The aim is to detect cancer specimens that are smaller or less detectable during standard biopsy procedures. The targeted objective schedule differs from recentlyproposed saturation biopsy protocols that involve performing a larger number of biopsies in the gland in an attempt to cover as much volume as possible.

The final step is to be able to carry out a pre-procedure planning focusing on assessment imaging techniques (spectroscopy and injected IRM) and to perform real-time ultrasound guidance after data registration. For the time being, it appears difficult to register imaging procedures without segmentation and real-time matching of the 2 techniques might be a problem, due in large part to processor computation times. The first attempts at computerassisted guidance, developed at TIMC Laboratory, is currently being tested at the Grenoble University Hospital Centre [START_REF] Bucki | PRONAV: a navigation software for prostate biopsies[END_REF] .

Conclusion:

The study is a step in pre-operative assessment of prostate biopsy paths using 3D ultrasound.

It also makes it possible to demonstrate that real-time 3D guidance on phantom biopsies offers greater accuracy in determining the needle target and makes biopsy reproducibility better. No evidence has been shown that biopsy distribution is more even in the prostate volume or that the biopsies were better spread.

Figure 1 :Figure 2 :

 12 Figure 1: Candle-gel prostate phantom

Figure 3 :

 3 Figure 3: Visualisation of 3 biopsy trails in the same phantom

Figure 4 :

 4 Figure 4: Visualisation of a biopsy path in the reference volume. The green ball shows the entry point. The blue ball shows the target hit.

Figure 5 :

 5 Figure 5: Superposition of all biopsies in the template.

Figure 6 :

 6 Figure 6: Summary methodology scheme

Figure 7 :

 7 Figure 7: Creation of ideal entry and target spheres for one biopsy.

Figure 8 :

 8 Figure 8: Overlapping of 2 "biopsy cylinders"

Figure 9 :

 9 Figure 9: Faulting our distribution assessment method. A significant-sized cancer goes undetected between 2 biopsy cylinders that do not overlap in the prostate volume.

Table 1 : Mean distances from entry and target zones in ideal protocol, in 2D and 4D modesTable 2 : Summary table regarding biopsy accuracy in reaching ideal theoretical protocol n Procedure Mean (mm) +/-SD Paired Student t test Distance from ideal target zone

 12

	158	2D	6.79 ± 6.18
		4D	5.1± 4.8	S
	158	2D	5.28 ± 5.71
				NS
	Distance from
	ideal entry zone	4D	5.19 ± 4.83

Table 3 : Interaction curve between differences in accuracy in hitting target zone in 2D and 4D, by location on prostate.

 3

		12
		10
		8
		6
	Moy. cell	0 2 4
		-2
		-4
		-6
		-8
		-10	ApexLateral	ApexMidLateral	BaseLateral	BaseMidLateral	MiddleLateral	MiddleMidLateral

Table 4 : Interaction curve between differences in accuracy at the entry zone in 2D and 4D, by prostate zone.

 4

		10
		8
		6
		4
	Moy. cell	0 2
		-2
		-4
		-6
		-8	ApexLateral	ApexMidLateral	BaseLateral	BaseMidLateral	MiddleLateral	MiddleMidLateral
	A biopsy reproducibility test was performed by looking at the variance in distances at the
	entry and target points. It showed that reproducibility increased significantly in biopsies
	toward the target zone with 4D US mode (p= 0.0385). However, there is no statistical link,
	but a trend toward improvement in reproducibility in needle entry points in the phantom
	(p=0.0623) (table 5).

Table 5 : Biopsy reproducibility study by deviation from mean² (variance) n

 5

			Procedure Deviation from	Standard	Student
				mean (mm)	deviation	t Test
		158	2D	37.95	78.45 0.0385
	Variance in					S
	distance
	compared to		4D	22.9	41.38
	ideal target
	zone
	Variance in	158	2D	32.43	59.66 0.0623
	distance					NS
	compared to
	ideal entry		4D		33.15
	zone
	The biopsy distribution data did not show any significant difference in the prostate volume
	explored under 2D and 4D US guidance (p=0.44), with a mean of 1025 mm³ per biopsy
	(12,3cc) per 4D guided biopsy session), compared to 971 mm³ under 2DUS guidance (11,65
	cc per 2D guided biopsy session). The sampled volume percentage not covered in several
	biopsies (sampled voxels through which only one biopsy runs) was 56.9% under 2D guidance,
	as compared to 56.7% under 4D guidance (p=0.98) (table 6).

Table 6 : Study of biopsy dispersion by measuring volume per biopsy and percentage of the prostate sampled by only one biopsy

 6

		N	Procedure	(median,	mean ± SD	p
				(interquartile
				range)
	Volume		2D	963 (209)	971 mm³
	sampled per	14			±166
	biopsy					NS
			4D	907 (258)	1025
					mm³±184
	Percentage of
	sampled	14	2D		0.569
	prostate

Acknowledgments:

Our thanks to the French Association of Urology (AFU), Sanofi-Aventis for the grant awarded, and the nationwide PHRC Prostate Echo 2003 project, which made this work possible.